;+

function trace_diffraction_pattern, nxin, nyin, wavelengthin, pixel_size

;NAME:
;     TRACE_DIFFRACTION_PATTERN
;PURPOSE:
;     Calculates the TRACE diffraction pattern.  Useful for then removing
;     the diffraction pattern.
;CATEGORY:
;CALLING SEQUENCE:
;     patt = trace_diffraction_pattern(nx,ny,wavelength,pixel_size)
;INPUTS:
;     nx,ny = size of array
;     wavelength = wavelength in A (string), e.g. '171', '195' or '284'.
;     pixel_size = pixel size in arc sec
;OPTIONAL INPUT PARAMETERS:
;KEYWORD PARAMETERS
;OUTPUTS:
;     patt = array holding the diffraction pattern.  Normalized so that 
;            total(patt) = 1.0.
;COMMON BLOCKS:
;SIDE EFFECTS:
;     For each order, 
;     prints  "order intensityNWSE intensitySWNE distanceNWSE distanceSWNE"
;RESTRICTIONS:
;     This code is still experimental!  It was calibrated for 171 and 195 from 
;     a single flare and assumes that the diffraction pattern does not change
;     with time.  Also, it has not yet been calibrated for 284. 
;PROCEDURE:
;     Interference and diffraction formulae for multiple slits.  The
;     parameters used in the model are the slope of the line connecting the
;     principal maxima, the wire diameter in the filter mesh, and the wire
;     spacing in the filter mesh.  These parameters were set from measurements 
;     of a flare.
;MODIFICATION HISTORY:
;   T. Metcalf  2000-March-04
;-

;      NW     NE
;        \   /
;         \ /
;          x
;         / \
;        /   \
;      SE     SE

   wavelength = float(wavelengthin)

   ; First set the parameters defining the TRACE filter mesh

   ; 70 lines per inch --> 3.6286e6 A wire spacing
   ; 82% transmission --> 3.266e5 A wire diameter

   ; (2.54 cm/in)(10^8 A/cm)/(70 L/in) = 3.6286e6 A/L
   ; Nx = # lines in x direction = X*70, where X is x size in inches
   ; Ny = # lines in y direction = Y*70, where Y is y size in inches
   ; Area of 1 line is L*D, where L is the length and D is the wire diameter
   ; Area of all wires is Nx*Y*D+Ny*X*D = 70*X*Y*D+70*Y*X*D = 140*X*Y*D
   ; Area of mesh is X*Y
   ; Block = 1-transmission = 1-0.82 = 0.18 = 140*X*Y*D/X*Y = 140*D
   ; D = 0.18/140 inches = (0.00129 in)(2.54 cm/in)(10^8 A/cm) = 3.266e5 A


   ; Theoretical values.  Should not be used
   ; since they differ from the  measure values
   ; significantly! 

   ;alpha_NWSE = !pi/4.             ; atan(slope_NWSE) in radians
   ;alpha_SWNE = !pi/4.             ; atan(slope_SWNE) in radians
   ;wire_diameter_NWSE = 3.266e5    ; in A
   ;wire_diameter_SWNE = 3.266e5    ; in A
   ;wire_spacing_NWSE = 3.6286e6    ; in A
   ;wire_spacing_SWNE = 3.6286e6    ; in A

   ; Measurements from the diffraction pattern of a flare

   slope171_NWSE = -1.0346
   slope171_SWNE =  0.9672
   wire_spacing171_NWSE = 3.580e6 ;3.534e6    ; in A
   wire_spacing171_SWNE = 3.580e6 ;3.534e6    ; in A
   wire_diameter171_NWSE = 0.115*wire_spacing171_NWSE
   wire_diameter171_SWNE = 0.115*wire_spacing171_SWNE
   ;wire_diameter171_NWSE = 4.463e5    ; in A
   ;wire_diameter171_SWNE = 4.463e5    ; in A

   slope195_NWSE = -1.1496
   slope195_SWNE =  0.8595
   wire_spacing195_NWSE = 3.591e6    ; in A
   wire_spacing195_SWNE = 3.591e6    ; in A
   wire_diameter195_NWSE = 0.115*wire_spacing195_NWSE
   wire_diameter195_SWNE = 0.115*wire_spacing195_SWNE
   ;wire_diameter195_NWSE = 4.463e5    ; in A
   ;wire_diameter195_SWNE = 4.463e5    ; in A

   ;slope284_NWSE = 
   ;slope284_SWNE =  
   ;wire_spacing284_NWSE =  3.580e6   ; in A
   ;wire_spacing284_SWNE =  3.580e6   ; in A
   ;;wire_diameter284_NWSE = 0.13*wire_spacing195_NWSE
   ;;wire_diameter284_SWNE = 0.13*wire_spacing195_SWNE
   ;wire_diameter284_NWSE = 4.463e5    ; in A
   ;wire_diameter284_SWNE = 4.463e5    ; in A

   if abs(wavelength-171.) LE 6 then {
      slope_NWSE = slope171_NWSE
      slope_SWNE = slope171_SWNE
      wire_diameter_NWSE = wire_diameter171_NWSE
      wire_diameter_SWNE = wire_diameter171_SWNE
      wire_spacing_NWSE = wire_spacing171_NWSE
      wire_spacing_SWNE = wire_spacing171_SWNE
   } else if abs(wavelength-195.) LE 6 then {
      slope_NWSE = slope195_NWSE
      slope_SWNE = slope195_SWNE
      wire_diameter_NWSE = wire_diameter195_NWSE
      wire_diameter_SWNE = wire_diameter195_SWNE
      wire_spacing_NWSE = wire_spacing195_NWSE
      wire_spacing_SWNE = wire_spacing195_SWNE
   } else if abs(wavelength-284.) LE 6 then {
      ty,'Uh oh.  No slope values are available for 284'
      slope_NWSE = slope284_NWSE
      slope_SWNE = slope284_SWNE
      wire_diameter_NWSE = wire_diameter284_NWSE
      wire_diameter_SWNE = wire_diameter284_NWSE
      wire_spacing_NWSE = wire_spacing284_NWSE
      wire_spacing_SWNE = wire_spacing284_NWSE
   }  else ty,'Wavelength is not valid'+string(wavelength)

   ;message,/info,'slope_NWSE: '+string(slope_NWSE)
   ;message,/info,'slope_SWNE: '+string(slope_SWNE)
   ;message,/info,'wire_diameter_NWSE: '+string(wire_diameter_NWSE)
   ;message,/info,'wire_diameter_SWNE: '+string(wire_diameter_SWNE)
   ;message,/info,'wire_spacing_NWSE: '+string(wire_spacing_NWSE)
   ;message,/info,'wire_spacing_SWNE: '+string(wire_spacing_SWNE)

   alpha_NWSE = atan(abs(slope_NWSE)) ;in radians
   alpha_SWNE = atan(abs(slope_SWNE)) ;in radians

   ; Set up some paramters and the ouput array

   I0 = 1.

   s_over_d_NWSE = float(wire_diameter_NWSE)/float(wire_spacing_NWSE)
   s_over_d_SWNE = float(wire_diameter_SWNE)/float(wire_spacing_SWNE)


   ; Pad the PSF array by adding a pixel around the edges.  This is necessary
   ; since a principal max wont't be computed in an edge pixel.  Remove the 
   ; padding later.

   nx = nxin + 2
   ny = nyin + 2

   patt = fltarr(nx,ny)

   ; Find the central pixel

   cx = floor(nx/2) ;(nx-1L)/2L
   cy = floor(ny/2) ;(ny-1L)/2L

   ; Set the central pixel to I0.  We'll renormalize later.

   patt(cx,cy) = I0

   ; Set the principal maxima 

   order = 1

   REPEAT {

      all_orders_out_of_bounds = 1

      phase_NWSE = order*!pi*(1.-s_over_d_NWSE)
      phase_SWNE = order*!pi*(1.-s_over_d_SWNE)
      intensity_NWSE = I0 * (sin(phase_NWSE)/phase_NWSE)^2
      intensity_SWNE = I0 * (sin(phase_SWNE)/phase_SWNE)^2

      sintheta_NWSE = order*wavelength/wire_spacing_NWSE
      sintheta_SWNE = order*wavelength/wire_spacing_SWNE

      if sintheta_NWSE LE 1 and sintheta_SWNE LE 1 then {

         theta_NWSE = asin(sintheta_NWSE)
         theta_SWNE = asin(sintheta_SWNE)

         ; convert theta from radians to arcsec and then to pixels

         theta_NWSE = theta_NWSE*!radeg*3600./pixel_size
         theta_SWNE = theta_SWNE*!radeg*3600./pixel_size

         x_NW = cx - theta_NWSE*cos(alpha_NWSE)
         x_SE = cx + theta_NWSE*cos(alpha_NWSE)
         x_SW = cx - theta_SWNE*cos(alpha_SWNE)
         x_NE = cx + theta_SWNE*cos(alpha_SWNE)

         y_NW = cy + theta_NWSE*sin(alpha_NWSE)
         y_SE = cy - theta_NWSE*sin(alpha_NWSE)
         y_SW = cy - theta_SWNE*sin(alpha_SWNE)
         y_NE = cy + theta_SWNE*sin(alpha_SWNE)

         ; To include the pixelation, embed the diffraction pixel in a much
         ; larger array and then rebin back down to the TRACE pixel size.
         ; This is necessary since the diffraction pattern will fall in 
         ; between pixels.

         tsz = 33
         tsz2 = floor(tsz/2)
         temp = fltarr(tsz*3,tsz*3)

         sample = 0   ; 0: Spread PSF over multiple pixels.  Good.
                      ; 1: Concentrate PSF in nearest pixel.  Bad.

         if x_NW GE 1. and y_NW LT ny-1. then {
            xx = floor(tsz*(1.+x_NW-floor(x_NW))+0.5)
            yy = floor(tsz*(1.+y_NW-floor(y_NW))+0.5)
            zero, temp
            ;temp(xx-tsz2:xx+tsz2,yy-tsz2:yy+tsz2) = intensity_NWSE
            temp(xx:xx+tsz-1,yy:yy+tsz-1) = intensity_NWSE
            patt(floor(x_NW)-1,floor(y_NW)-1) = rebin(temp,3,3,sample=sample)
            all_orders_out_of_bounds = 0
         }
         if x_SE LT nx-1. and y_SE GE 1. then {
            xx = floor(tsz*(1.+x_SE-floor(x_SE))+0.5)
            yy = floor(tsz*(1.+y_SE-floor(y_SE))+0.5)
            zero, temp
            ;temp(xx-tsz2:xx+tsz2,yy-tsz2:yy+tsz2) = intensity_NWSE
            temp(xx:xx+tsz-1,yy:yy+tsz-1) = intensity_NWSE
            patt(floor(x_SE)-1,floor(y_SE)-1) = rebin(temp,3,3,sample=sample)
            all_orders_out_of_bounds = 0
         }
         if x_SW GE 1. and y_SW GE 1. then {
            xx = floor(tsz*(1.+x_SW-floor(x_SW))+0.5)
            yy = floor(tsz*(1.+y_SW-floor(y_SW))+0.5)
            zero, temp
            ;temp(xx-tsz2:xx+tsz2,yy-tsz2:yy+tsz2) = intensity_SWNE
            temp(xx:xx+tsz-1,yy:yy+tsz-1) = intensity_SWNE
            patt(floor(x_SW)-1,floor(y_SW)-1) = rebin(temp,3,3,sample=sample)
            all_orders_out_of_bounds = 0
         }
         if x_NE LT nx-1. and y_NE LT ny-1. then {
            xx = floor(tsz*(1.+x_NE-floor(x_NE))+0.5)
            yy = floor(tsz*(1.+y_NE-floor(y_NE))+0.5)
            zero, temp
            ;temp(xx-tsz2:xx+tsz2,yy-tsz2:yy+tsz2) = intensity_SWNE
            temp(xx:xx+tsz-1,yy:yy+tsz-1) = intensity_SWNE
            patt(floor(x_NE)-1,floor(y_NE)-1) = rebin(temp,3,3,sample=sample)
            all_orders_out_of_bounds = 0
         }

         ;message,/info,strcompress(string(order) + $
          ;                         string(intensity_NWSE) + $
          ;                         string(intensity_SWNE) + $
          ;                         string(theta_NWSE) + $
          ;                        string(theta_SWNE))

         order = order + 1L

      }

   } UNTIL all_orders_out_of_bounds
   
   ty,'Maximum order is ',order-1

   patt = patt(1:nx-2,1:ny-2) ; Remove the padding

   return,patt/total(patt)   ; renormalize to 1 for a deconvolution

end