([filter-limit],[maxspfreq],[next-iter])
;
; SOURCE: enhance.pam       
;
; PURPOSE: Calculate 1D PW, then create scattering curve in the same format,
;          Search for close frequences -- intensity can be normalized / 
;          multiplied or divided to get it closer to PW values
;
; INPUT REGISTERS:
;    [filter-limit]          Filter limit (in Fourier pixels), i.e., resolution cutoff
;    [maxspfreq]             Max. spatial frequency (0.5/pixel_size)
;    [next-iter]             Iteration counter
;
; INPUT FILES:
;    final/val{**[next-iter]}          [vol]               Volume to be corrected
;    input/scattering          [scattering_doc]    X-ray scattering power spectrum
;
; OUTPUT FILES:
;    work/enhance_doc_{**[next-iter]}  [enhance_doc]       Enhancement doc file
;    work/vol_ps_{**[next-iter]}       [vol_power_spec]    Temp file (deleted)
;
;.......................................................................

FR                                   ; Input file
?volume to be corrected?[vol]

FR L
[vol_power_spec]rojo{**[next-iter]}          ; Temp output file


x80 = 0.5 / [maxspfreq]  ; WAS 2.82          ; Pixel size (in A)

; Create doc file with 1D power spectrum of input volume

PW                ; Power spectrum of input volume
[vol]
_1

SQ                ; Square the values in vol
_1
_2

RO                ; Rotational average
_2
_3

DE                ; Delete existing file
[vol_power_spec]

LI D              ; List volume to doc file
_3
[vol_power_spec]  ; Output doc file
R                 ; Real format
1                 ; Register number

UD N [size]          ; Find size of doc file
[vol_power_spec]

; Create output enhancement curve
[x55] = 0.0
[x56] = 100

DE
[enhance_doc]

DO LB1 [x21]=2,[filter-limit]      ; Curve goes out to filter limit

   UD IC [x21], x54
   [vol_power_spec]

   [x55] = ([x21]-1)/(2*([size]-1))
   [x56] = x80/[x55]     ; Spatial frequency in Angstroms,
   ; Output [x77]: corresponding intensity (from 2nd col in scatter doc file)
   ; The filter is scaled to have 1 at 100A


   UD N [npts]           ; Number of points in scattering file
   [scattering_doc]

   [x52] = INT([npts]/2)

   DO LB11 [ipt]=1,[npts]

      UD IC [x52], [x31],[x22]
      [scattering_doc]

      [x52] = [x52] + 1

      IF ([x52] .GT. 1500) GOTO LB12

      UD IC [x52], [x32],x23
      [scattering_doc]

      IF ([x56] .LE. [x31]) THEN
         IF ([x56] .GE.[x32])  THEN
            x77 = ([x31]-[x56])**2
            [x78] = ([x32]-[x56])**2
            IF ([x78] .LT. x77) [x22] = x23
            GOTO LB12
         ENDIF
      ENDIF

      IF ([x56] .LT.[x32]) THEN
         [ipt]  = [x52]
         [x52]  = INT(([npts]-[x52])/2+[x52])
      ELSE
         [npts] = [x52]
         [x52]  = INT(([x52]-[ipt])/2+[ipt])
      ENDIF

   LB11

   LB12
   [x22] = [x22] / 37891.

   ; Xray/EM
   [x22] = SQR([x22]/x54)
   [x78] = LON([x22])

   SD [x21], [x22],[x55],[x56],[x78]
   [enhance_doc]

LB1

UD ICE
[vol_power_spec]
UD ICE
[scattering_doc]

DE               ; No longer needed
[vol_power_spec] ; Temp output file

RE