;adapted from Tom Metcalf's TRACE_DIFFRACTION_PATTERN.pro
;simplified here for eventual porting to C or Java

;Tom defines the two directions of the grid empirically and allows the wire
;thickness and spacing to be different for each, two lines of diffraction
;points are defined. The finite thickness of the wire causes the intensity
;of the orders to vary. The order spacing depends on the angular scale
;on the detector. This was also done empirically but gives a wire spacing
;for a known wavelength. Apparently these are quite different from the
;theoretical values based on the physical grid spacing. Why?

subr insert_gauss(big,gwid,nw,x,y)
 ;get the ranges
 i1 = rfix(x - nw/2)
 j1 = rfix(y - nw/2)
 ; the center has to be at (x,y)
 gaussx = gauss( (indgen(fltarr(nw)) + i1 - x)/gwid )
 gaussy = gauss( (indgen(fltarr(nw)) + j1 - y)/gwid )
 redim, gaussx, nw, 1
 redim, gaussy, 1, nw
 spot = gaussx * gaussy
 insert, big, spot, i1, j1
 endsubr
;func trace_diffract(wave)
wave = 171
 ;Tom's values are packed into arrays here. The order is 171, 195

 slope_NWSE = [ -1.0346, -1.1496]
 slope_SWNE = [ 0.9672, 0.8595]

 wire_frac_NWSE = [ 0.115, 0.115]
 wire_frac_SWNE = [ 0.115, 0.115]

 wire_spacing_NWSE = [3.580e6, 3.591e6]; in A
 wire_spacing_SWNE = [3.580e6, 3.591e6]; in A

 ;the spacing of the orders on the detector depends on the pixel size, adopt
 ;Tom's value here for consistency. The spacing is about 20 pixels so we need
 ;about 25 orders to cover the detector

 pixel_size = 0.4993

 if abs(wave-171) le 10 then nwave = 0 else
  if abs(wave-195) le 10 then nwave = 1 else
    { ty,'input wavelength not legal: ', wave return, 0 }

 ;estimate # of orders needed to fill this 1024x1024 image, assuming equal
 ;spacing
 xq = (wave/wire_spacing_NWSE(nwave)) * #r_d *3600./pixel_size
 norder = rfix(1.5*512/xq)
 ty,'norder =', norder

 ;convert spacing to pixels in the 2 directions, assume equally spaced on
 ;the detector, the angles are very small. Tom took the asin of equal spacing
 ;but I'm not sure if that is really any more accurate, a full treatment
 ;is much more complicated. Anyhow, the largest angle is about 0.0023 radians
 ;so we should be very OK with equal spacing.
 xq = #r_d*3600./pixel_size
 s_nwse = (indgen(fltarr(norder))*wave/wire_spacing_NWSE(nwave))* xq
 s_swne = (indgen(fltarr(norder))*wave/wire_spacing_SWNE(nwave))* xq
 xq = atan(slope_NWSE(nwave))
 x_nwse = s_nwse * cos(xq)
 y_nwse = s_nwse * sin(xq)
 xq = atan(slope_SWNE(nwave))
 x_swne = s_swne * cos(xq)
 y_swne = s_swne * sin(xq)

 ;compute the intensities of the two sets of peaks
 orders = indgen(fltarr(norder))
 xq = #pi*orders*wire_frac_nwse(nwave)
 xq(0) = 1.E-9
 $peaks_nwse = ( sin(xq)/xq) ^ 2
 xq = #pi*orders*wire_frac_swne(nwave)
 xq(0) = 1.E-9
 $peaks_swne = ( sin(xq)/xq) ^ 2

 ;each of these will be a gaussian of width gwid (in pixels)
 gwid = 1.0
 ;full size is nw
 nw =21
 gaussx = gauss( (indgen(fltarr(nw))- (nw-1)/2.)/gwid )
 gaussy = gaussx
 redim, gaussx, nw, 1
 redim, gaussy, 1, nw
 spot = gaussx * gaussy
 
 xim = zero(fltarr(1024,1024)
 
 ;loop through the nwse set and insert a gaussian blob at each point
 for i=0,norder-1 do {
  x = x_nwse(i)
  y = y_nwse(i)
  ;get the corner position for the gaussian square
  x = x -.5*(nw-1)
  dx = x - rfix(x)
  ty,'x, dx =', x, dx
  
 }

 ;return, x
 ;endfunc