;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, fac)
 ;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 = fac * 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*(1.0-wire_frac_nwse(nwave))
 xq(0) = 1.E-9
 peaks_nwse = ( sin(xq)/xq) ^ 2
 xq = #pi*orders*(1.0-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 =13 

 xim = zero(fltarr(1024,1024)
 ixc = 511.5
 iyc = 511.5
 
 ;loop through the nwse set and insert a gaussian blob at each point
 ;for i=0,norder-1 do {
 ; for j=0,3 do {
 ; ncase j
 ;  { x =  x_nwse(i) + ixc  y =  y_nwse(i) + iyc  fac = peaks_nwse(i) }
 ;  { x = -x_nwse(i) + ixc  y = -y_nwse(i) + iyc  fac = peaks_nwse(i) }
 ;  { x =  x_swne(i) + ixc  y =  y_swne(i) + iyc  fac = peaks_swne(i) }
 ;  { x = -x_swne(i) + ixc  y = -y_swne(i) + iyc  fac = peaks_swne(i) }
 ; endcase
  ;these can be used directly by insert_gauss
  insert_gauss, xim, gwid, nw, x, y, fac
  }
 }
 for i=0,norder-1 do {
 for k=0,norder-1 do {
  for j=0,3 do {
  ncase j
   { x =  x_nwse(i) + ixc + x_swne(k)
     y =  y_nwse(i) + iyc + y_swne(k)
     fac = peaks_nwse(i)*peaks_swne(k) }
   { x = -x_nwse(i) + ixc - x_swne(k)
     y = -y_nwse(i) + iyc - y_swne(k)
     fac = peaks_nwse(i)*peaks_swne(k) }
   { x =  x_swne(i) + ixc - x_nwse(k)
     y =  y_swne(i) + iyc - y_nwse(k)
     fac = peaks_swne(i)*peaks_nwse(k) }
   { x = -x_swne(i) + ixc + x_nwse(k)
     y = -y_swne(i) + iyc + y_nwse(k)
     fac = peaks_swne(i)*peaks_nwse(k) }
  endcase
  ;these can be used directly by insert_gauss
  insert_gauss, xim, gwid, nw, x, y, fac
  }
 }}
 !iorder=0
 tv,alog(xim>1e-6) 
 ;return, xim
 ;endfunc