Initial commit

1 files changed, 210 insertions, 0 deletions

diff --git a/idl/cf_plot_extract3.pro b/idl/cf_plot_extract3.pro
new file mode 100644
index 0000000..c201754
--- /dev/null
+++ b/idl/cf_plot_extract3.pro
@@ -0,0 +1,210 @@
+;+
+;  cf_plot_extract3.pro is a procedure to plot the extraction windows over a 
+;                   greyscale image for a given IDF to a standard image file. 
+;
+;  Author:  Edward M. Murphy
+;  Written: 1999 November 24
+;
+;       Calling sequence: cf_plot_extract3, rootname
+;
+;       Inputs:    rootname will have 'idf.fit' added to it.
+;
+; 02/12/03 jch: Adapt from cf_plot_extract.pro to work for CalFUSE 3.0.
+; 03/10/03 wvd: Silence calls to mrdfits.  Simplify image construction.
+; 04/08/03 wvd: Shift apertures in X to match FPA positions.
+; 06/02/03 wvd: Don't plot non-target apertures in HIST mode.
+; 10/10/03 wvd: Exclude photons from bad times and those outside PHA limits.
+; 12/09/03 wvd: Use extended apertures when appropriate.
+; 02/26/04 wvd: Use strcmp to compare strings.
+; 05/03/04 wvd: In HIST mode, don't plot non-target apertures.
+;		Read expected spectral centroid from CHID_CAL file header.
+; 05/14/04 wvd: Don't use strcmp to compare strings.
+;		Modify boundaries between LiF and SiC channels.
+; 06/14/04 wvd: Don't shift apertures in X to match FPA positions.
+; 07/12/05 wvd: Clean up call to HIST_2D.
+; 08/01/05 wvd: If IDL version > 6.0, generate GIF file.
+; 08/30/05 wvd: If there are no good photons, plot all photons.
+; 11/01/05 wvd: If DAYNIGHT = NIGHT, exclude daytime photons.
+; 03/14/08 wvd: Require LOC_FLAG < 16B for good data.
+;-
+
+pro cf_plot_extract3,rootname
+
+	!quiet=1
+	filename=rootname+'idf.fit'
+	a=mrdfits(filename,0,ah,/SILENT)
+
+	; Read input file's primary header keywords
+	wave_cal = strcompress(fxpar(ah,'WAVE_CAL'),/rem)
+	chid_cal = strcompress(fxpar(ah,'CHID_CAL'),/rem)
+	filedate = strcompress(fxpar(ah,'DATE'),/rem)
+	obsdate = strcompress(fxpar(ah,'DATEOBS'),/rem)
+	obstime = strcompress(fxpar(ah,'TIMEOBS'),/rem)
+	cfvers = strcompress(fxpar(ah,'CF_VERS'),/rem)
+	exptime = fxpar(ah,'EXPTIME')
+	instmode = strcompress(fxpar(ah,'INSTMODE'),/rem)
+	detector = strcompress(fxpar(ah,'DETECTOR'),/rem)
+	aperture = strcompress(fxpar(ah,'APERTURE'),/rem)
+	srctype = strmid(fxpar(ah,'SRC_TYPE'),0,1)
+	daynight = strcompress(fxpar(ah,'DAYNIGHT'),/rem)
+
+        ; Interpret aperture and srctype keywords.
+        if (aperture EQ 'MDRS') then begin
+		aplif = 2
+		apsic = 6
+        endif else if (aperture EQ 'HIRS') then begin
+		aplif = 1
+		apsic = 5
+        endif else begin	; assume LWRS
+		aplif = 3
+		apsic = 7
+        endelse
+	if (srctype EQ 'P') then extend = 0 else extend = 8
+
+	cal_path = getenv("CF_CALDIR")
+	slitname=["HIRS","MDRS","LWRS","PINH","HIRS","MDRS","LWRS","PINH"]
+	specname=["LiF","LiF","LiF","LiF","SiC","SiC","SiC","SiC"]    
+
+	; Set the plot resolution.
+	set_plot,'Z'
+	device,set_resolution=[1124,612]
+
+	; Define plot attributes
+	!x.style=1
+	!y.style=1
+	!x.tickformat='(I5)'
+	!p.title=filename
+	!x.title="X pixel"
+	!y.title="Y pixel"
+	!p.charsize=1.2
+	cs=0.8
+
+        ; Read the X and Y columns from the first extension (photon list),
+        ; exclude bad events, and generate an image from the data.
+        ftab_ext, filename, 'X,Y,TIMEFLGS,LOC_FLGS', xarr, yarr, timeflag, loc_flag
+	if (daynight EQ 'NIGHT') then begin
+            good = where (timeflag eq 0 and loc_flag lt 16, n) 
+	endif else begin
+            good = where ((timeflag and not 1B) eq 0 and loc_flag lt 16, n) 
+	endelse
+        if (n gt 0) then begin
+            xarr = xarr[good]
+            yarr = yarr[good]
+        endif
+        a = HIST_2D(xarr,yarr,bin1=16,bin2=2,min1=0,max1=16383,min2=0,max2=1023)
+        b = HIST_EQUAL(a)
+
+	!p.charsize=1.3
+	x0=85
+	y0=65
+	x1=1024+x0
+	y1=512+y0
+	dummy_x = indgen(2)*16383
+	dummy_y = indgen(2)*1023
+
+	; Plot the pixel axes and the photon image.
+	plot,dummy_x,dummy_y,/nodata,position=[x0,y0,x1,y1],/DEVICE,$
+		background=255,color=0,xticks=8
+	tv,255-b,x0,y0
+	plot,dummy_x,dummy_y,/nodata,/noerase,position=[x0,y0,x1,y1],/DEVICE,$
+		background=255,color=0,xticks=8,xminor=4
+
+	miny=10000
+	maxy=0.0
+
+	; Plot the extraction windows.
+	for i=1,7 do begin
+	    ; Skip the pinhole and any non-target apertures for HIST data.
+	    if ((i ne 4) and ((instmode EQ 'TTAG') $
+                or (i eq aplif) or (i eq apsic))) then begin
+
+		; Read aperture limits from calibration file.
+		j = 8; use extended apertures for non-target apertures
+		if ((i eq aplif) or (i eq apsic)) then j = extend
+		xw=mrdfits(cal_path+'/'+chid_cal,i+j,chidhdr,/SILENT)
+		yhigh = xw.yhigh
+		ylow  = xw.ylow 
+		default_y_centroid = fxpar(chidhdr,'CENTROID')
+
+		; Get the YCENT keywords and adjust the windows accordingly.
+		key = 'YCENT' + strcompress(string(i),/rem)
+		ycent = fxpar(ah, key)
+		shift = nint(default_y_centroid - ycent)
+		yhigh = yhigh - shift
+		ylow  = ylow - shift
+
+		; Plot the extraction windows.
+		oplot,yhigh,color=0
+		oplot,ylow,color=0
+		miny=min([miny,ylow])
+		maxy=max([maxy,yhigh])
+		xyouts,15600,ylow[16350]+5,specname[i-1],$
+			color=0,charsize=0.8
+		xyouts,15600,ylow[16350]-20,slitname[i-1],$
+			color=0,charsize=0.8
+	    endif
+	endfor
+
+	; Write history information.
+	xyouts,50,20,'Observation date='+obsdate+'T'+obstime,$
+		charsize=0.90,color=0,/DEVICE
+	xyouts,50,5,'IDF file date= '+filedate,charsize=0.90,color=0,/DEVICE
+	xyouts,1124/2,5,'CALFUSE version '+cfvers,$
+		charsize=0.90,color=0,/DEVICE,align=0.5
+	xyouts,1124-50,20,'CHID calibration file='+chid_cal,$
+		charsize=0.90,color=0,/DEVICE,align=1.0
+	xyouts,1124-50,5,'Exposure time='+$
+		strcompress(string(exptime,format="(F10.2)"),/rem),$
+		charsize=0.90,color=0,/DEVICE,align=1.0
+
+	;  Read in and plot the wavelength scale
+	wl_lif=mrdfits(cal_path+'/'+wave_cal,4,ah,/SILENT)
+	wl_sic=mrdfits(cal_path+'/'+wave_cal,8,ah,/SILENT)
+
+	wlsc_lif=[wl_lif[0].wavelength,wl_lif[16383].wavelength]
+	wlsc_sic=[wl_sic[0].wavelength,wl_sic[16383].wavelength]
+
+	detector=strtrim(detector,2)
+
+	; Determine locations of the wavelength axes.
+	if (detector EQ '1A') then begin
+		hip=maxy+25
+		midp=425 
+		lowp=miny-15
+	endif else if (detector EQ '1B') then begin 
+		hip=maxy+25
+		midp=425
+		lowp=miny-10
+	endif else if (detector EQ '2A') then begin
+		hip=maxy+100
+		midp=500
+		lowp=miny-100
+	endif else if (detector EQ '2B') then begin
+		hip=maxy+100
+		midp=500
+		lowp=miny-100
+	endif else begin
+		hip=maxy+25
+		midp=512
+		lowp=miny-15
+	endelse
+
+	; Draw wavelength axes.
+	axis,0,hip,xaxis=1,xrange=wlsc_lif,xstyle=1,charsize=cs,$
+		color=0,xtitle="Wavelength"
+	axis,0,midp,xaxis=0,xrange=wlsc_lif,xstyle=1,charsize=0.001,$
+		color=0,xtitle="Wavelength"
+
+	axis,0,midp,xaxis=1,xrange=wlsc_sic,xs=1,charsize=0.001,$
+		color=0,xtitle=""
+	axis,0,lowp,xaxis=0,xrange=wlsc_sic,xs=1,charsize=cs,$
+		color=0,xtitle="Wavelength"
+
+	; Create a GIF/JPEG file: 
+	ver = float(!version.release)
+	if (ver ge 5.4 and ver le 6.0) then write_jpeg,rootname+'ext.jpg',TVRD() $
+	else write_gif,rootname+'ext.gif',TVRD()
+
+EXIT:
+return
+end