; DEF.SPI estimates defocus from power spectra profiles. The program def
;         must be somewhere on the executable path.
;
; Inputs:
;    roo{***} : rotationally averaged power spectra profiles
;    filenums : list of micrographs
;    params   : parameter file
;
; Output:
;    defocus : doc file of defocus estimates. This doc file is not
;              renumbered, but instead is left with the keys
;              corresponding to micrograph numbers.
;

x11 = 60000  ; maximum defocus value 

; ----------- Input files --------------
FR G
[params]../params      ; parameter file	
FR G
[FILENUMS]../filenums     ; file numbers
FR G
[roo]power/roo{****x12}   ; rotational average (doc file)

; ----------- Output files --------------
FR G
[defocus]defocus    ; the output doc file

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

; if the output file already exists, delete it
DE
[defocus]

;  create a 2D transfer function for a range of defoci
;    x axis = spatial frequency
;    y axis = defocus from 1 to max, increment = 100 defocus units per line
;
; NB max defocus value is divided by 100. If max = 40000, the
; ctf file will have 400 elements along the defocus axis.

FR G
[ctf]ctf001

; get values from parameter file
UD  6,x21  ; KV, for computing lambda
[params]  

UD 5,x31   ; pixelsize, for computing max spatial freq.
[params]

UD 7,x33  ; spherical aberration
[params]

x11 = x11 / 100.0
x22 = 12.398 / SQR(x21 * (1022.0 + x21))
x32 = 1.0 / (2.0 * x31)

TF
[ctf]        ; output image
x33        ; spherical aberration coefficient in mm.
x22        ; lambda
1.00       ; lower defocus limit
(x11 * 100)  ; upper defocus limit
(250, x11)     ; (# spat. freq. points, # defocus points)
x32        ; max spatial frequency
0.000      ; source size (set to zero to eliminate envelope)
0.00       ; defocus spread
0.10, 100.    ; (amplitude contrast ratio, Gaussian envelope parameter) 
D          ; Diffractogram
N          ; No frame

; end create ctf -----------------

; labels for columns
SD /     micrograph   defocus
[defocus]
SD E
[defocus]

; get the filenumbers
ud n,x20
[FILENUMS]

DO LB1 x11 = 1,x20
   UD x11,x12
   [FILENUMS]
  ; x12 is now the file number

   VM
   def [ctf].$DATEXT [roo].$DATEXT ]] [defocus].$DATEXT
LB1

DE
[ctf]

EN D