; ctf_compare : compare the power spectrum from a micrograph with
;               a model. Output is half exp, half model power spectrum.

x55 = 6      ; file number
x60 = 35499  ; defocus

x65 =  0.075   ; Gaussian envelope parameter. 
             ; A LARGE number (100) = do not use envelope.
             ; zero = get Gaussian envelope width from parameter file.


; ----- input files -----

FR L
[pow]power/pw_avg{***x55}    ; power spectrum image from micrograph

FR L
[params]../params          ; parameter file

; ----- output files -----

FR L
[model]_11   ; model{***x55}    ; 2D model power spectrum

FR L
[out]split{***x55}

; ----- END BATCH HEADER -----

; Get the size of the experimental power spectrum image (must be square)
FI x11,x12
[pow]
(12,2)

IF (x11.NE.x12) THEN
   VM
   echo "power spectrum file must be a square image"
   EN D
ENDIF

; get the maximum, minimum values of the power spectrum file
FS x21,x22
[pow]

; ----------- make the model power spectrum with TF D ---------------

; get parameters
UD 7,x31   ; spherical aberration
[params]
UD 14,x32   ; lambda
[params]
UD 15,x33   ; max. spatial frequency
[params]
UD 8,x34   ; source size
[params]
UD 9,x35   ; defocus spread
[params]
UD 10,x36   ; astigmatism
[params]
UD 11,x37   ; azimuth
[params]
UD 12,x38   ; amplitude contrast ratio
[params]

x39 = 1000   ; Gaussian envelope halfwidth
IF (x65.EQ.0) THEN
   UD 13,x39
   [params]
ELSE
   x39 = x65
ENDIF

TF D
_1          ; output file
x31         ; Cs
x60,x32     ; defocus, lambda
x11         ; image size
x33         ; maximum spatial frequency
x34,x35     ; source size, defocus spread
x36,x37     ; astigmatism, azimuth
x38,x39     ; amplitude contrast ratio, gaussian envelope parameter
D           ; (D)iffractogram/(E)nvelope/(S)traight

; rescale the model
AR SCA
_1      ; input
_2      ; output
x22,x21 ; new min and max

CP
_2
[model]

@half
[pow]
[model]
[out]

EN D