;+ 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