aboutsummaryrefslogtreecommitdiff
path: root/idl
diff options
context:
space:
mode:
Diffstat (limited to 'idl')
-rw-r--r--idl/cf_obsplot.pro334
-rw-r--r--idl/cf_plot_extract3.pro210
-rw-r--r--idl/cf_plot_rate3.pro227
-rw-r--r--idl/cf_plot_rate_linear.pro175
-rw-r--r--idl/fuse_extract.pro154
-rw-r--r--idl/fuse_sic_only.pro159
6 files changed, 1259 insertions, 0 deletions
diff --git a/idl/cf_obsplot.pro b/idl/cf_obsplot.pro
new file mode 100644
index 0000000..9a4c8c2
--- /dev/null
+++ b/idl/cf_obsplot.pro
@@ -0,0 +1,334 @@
+;+
+; cf_obsplot.pro reads the final, combined *00000all* spectral file
+; for a given observation and produces quick-look plots for
+; each of the four FUSE detectors. A fifth plot shows the
+; best three channels that span the regions 900-1000,
+; 1000-1100, and 1100-1200 A.
+;
+; Author: Ken Sembach
+; Last Update: February 23, 2001
+;
+; V. Dixon, 07/02/2001 - Modified for inclusion in pipeline.
+; V. Dixon, 08/16/2001 - Read and write to current directory.
+; V. Dixon, 08/31/2001 - Add a sixth plot spanning FUSE band.
+; Mark region of worm in LiF 1B.
+; V. Dixon, 11/05/2001 - Compute mean flux with 1-sig errors.
+; V. Dixon, 11/20/2001 - Modify position of mean-flux plots.
+; Add sentence describing error bars.
+; V. Dixon, 01/17/2002 - Deal gracefully with missing data.
+; V. Dixon, 03/12/2002 - Create .gif files if possible.
+; B. Godard, 04/26/2004 - Print units and CalFUSE Ver
+; Don't keep exposures with exptime=0
+; or quality array=0
+; Use quality array to compute min/max
+; V. Dixon, 05/04/2004 - Scale ymax by 1.1 in all plots.
+; Change output file names to match
+; MAST's expectations.
+; V. Dixon, 08/20/2004 - Change A to Angstroms at the bottom
+; of each page.
+; Exclude C III 977 when setting Y
+; scale.
+; V. Dixon, 04/14/2005 - Rewrite the program to use the new
+; *all* files produced by CalFUSE v3.1.
+; V. Dixon, 08/01/2005 - If IDL version > 6.0, generate GIF file.
+; V. Dixon, 03/21/2006 - If channel has no data, set mflux = -1.
+; Use "Exposure" or "Exposures" as
+; appropriate.
+; Consider all four channels for the
+; 1000-1100 A plot.
+; Pass xrange and yscale to TOP_PLOT
+; and BOTTOM_PLOT.
+; V. Dixon, 05/19/2006 - For BKGD targets, omit flux
+; comparison when selecting which
+; detector segments to use.
+; Compute mean flux using same wavelength
+; intervals for each band.
+; V. Dixon, 05/24/2006 - For PC targets, don't include O VI in
+; calculation of mean flux.
+; V. Dixon, 12/12/2006 - If OBSTIME < 10 s for preferred segment
+; and > 100 s for other segment, use other.
+; V. Dixon, 12/19/2006 - Use the same sample regions in all plots.
+; V. Dixon, 04/06/2007 - Comment out calls to set_plot,'X' and
+; CLEANPLOT, as they can cause problems
+; for OPUS.
+; V. Dixon, 04/11/2008 - If obstime < 1000, quote in seconds.
+; V. Dixon, 08/08/2008 - If airglow keyword is set, image files
+; are named *00900*
+;
+; Calling sequence: OBSPLOT, obsname
+;
+; Input: obsname (e.g., P1030604)
+; Output: obsname00000lif1ttagf.jpg, obsname00000lif2ttagf.jpg,
+; obsname00000sic1ttagf.jpg, obsname00000sic2ttagf.jpg,
+; obsname00000specttagf.jpg
+;
+;
+;-
+
+pro TOP_PLOT, data, header, xrange, yscale, aper, sample, mflux, title
+ s=sample
+ !x.range=xrange
+ obstime = SXPAR(header,'OBSTIME')
+ extname = SXPAR(header,'EXTNAME')
+ num_exp = SXPAR(header,'NSPEC')
+ if (num_exp eq 1) then exposures=' Exposure, ' else exposures=' Exposures, '
+ IF (obstime gt 0) THEN BEGIN
+ w = data.wave & f = data.flux
+ loc = WHERE(((w GT s[0]) AND (w LE s[1])) OR ((w GT s[2]) and (w LT s[3])) OR $
+ ((w GT s[4]) and (w LT s[5])))
+ ftmp = SMOOTH(f[loc],6) & nftmp = N_ELEMENTS(ftmp)
+ ymin = MIN(f[loc]) & ymax = MAX(ftmp)
+ ydiv = FIX(1-ALOG10(ymax))
+ ymax = yscale * ymax *10.^ydiv
+ !y.range=[ymin,ymax]
+ mflux = mean(f[loc])
+
+ if (obstime lt 1000) then time_str = STRING(round(obstime),'(I3)')+' sec)' else $
+ time_str = STRING(obstime/1000,'(F4.1)')+' ksec)'
+
+ PLOT,w,f*10.^ydiv,nsum=12,title=title,ytitle='Flux (10!u-'+STRTRIM(ydiv,2)+'!n)',xtitle=' '
+ XYOUTS,(!x.range[1]-!x.range[0])*0.5+!x.range[0],(!y.range[1]-!y.range[0])*0.88+!y.range[0], $
+ extname+' '+aper+' ('+STRTRIM(num_exp,2)+exposures+time_str,charsize=1.3,align=0.5
+ ENDIF ELSE BEGIN
+ !y.range=[0,1]
+ mflux = -1.
+ PLOT,[0,0],[0,0],title=title,ytitle='Flux',xtitle=' '
+ XYOUTS,(!x.range[1]-!x.range[0])*0.5+!x.range[0],0.88, $
+ extname+' '+aper+': No data found',charsize=1.5,align=0.5
+ ENDELSE
+ RETURN
+END
+
+pro BOTTOM_PLOT, data, header, xrange, yscale, aper, sample, mflux, vers
+ s=sample
+ !x.range=xrange
+ obstime = SXPAR(header,'OBSTIME')
+ extname = SXPAR(header,'EXTNAME')
+ num_exp = SXPAR(header,'NSPEC')
+ if (num_exp eq 1) then exposures=' Exposure, ' else exposures=' Exposures, '
+ IF (obstime gt 0) THEN BEGIN
+ w = data.wave & f = data.flux
+ loc = WHERE(((w GT s[0]) AND (w LE s[1])) OR ((w GT s[2]) AND (w LT s[3])) OR $
+ ((w GT s[4]) AND (w LT s[5])))
+ ftmp = SMOOTH(f(loc),6) & nftmp = N_ELEMENTS(ftmp)
+ ymin = MIN(f[loc]) & ymax = MAX(ftmp)
+ ydiv = FIX(1-ALOG10(ymax))
+ ymax = yscale * ymax *10.^ydiv
+ !y.range=[ymin,ymax]
+ mflux = mean(f[loc])
+
+ if (obstime lt 1000) then time_str = STRING(round(obstime),'(I3)')+' sec)' else $
+ time_str = STRING(obstime/1000,'(F4.1)')+' ksec)'
+
+ PLOT,w,f*10.^ydiv,nsum=12,title=' ',ytitle='Flux (10!u-'+STRTRIM(ydiv,2)+'!n)', $
+ xtitle='Wavelength ('+string("305B)+')',/noeras
+ XYOUTS,(!x.range[1]-!x.range[0])*0.5+!x.range[0],(!y.range[1]-!y.range[0])*0.88+!y.range[0], $
+ extname+' '+aper+' ('+STRTRIM(num_exp,2)+exposures+time_str, charsize=1.3,align=0.5
+ IF (STRCMP(extname,'1BLIF',5, /FOLD_CASE) EQ 1) THEN $
+ XYOUTS, 1160., (!y.range[1]-!y.range[0])*0.1+!y.range[0],'---- REGION OF WORM ----',align=0.5
+ ENDIF ELSE BEGIN
+ !y.range=[0,1]
+ mflux = -1.
+ PLOT,[0,0],[0,0],title=' ',ytitle='Flux', xtitle='Wavelength ('+string("305B)+')',/noeras
+ XYOUTS,(!x.range[1]-!x.range[0])*0.5+!x.range[0],0.88, $
+ extname+' '+aper+': No data found',charsize=1.5,align=0.5
+ ENDELSE
+ XYOUTS, 0.05, 0.01, 'Flux units are erg cm!E-2!N s!E-1!N '+string("305B)+'!E-1!N.', alignment=0., /NORMAL
+ XYOUTS, 0.95, 0.01, 'CalFUSE v'+vers, alignment=1., /NORMAL
+ RETURN
+END
+
+PRO CF_OBSPLOT, obsname, airglow=airglow
+
+ ver = float(!version.release)
+
+ if keyword_set(airglow) then fnumber = '00900' else fnumber = '00000'
+ list = FINDFILE(obsname+fnumber+'all*fcal.fit', COUNT=count)
+ lun = fxposit(list[0], 0, /SILENT)
+ foo = MRDFITS(lun,0,h,/SILENT)
+
+ root = SXPAR(h,'PRGRM_ID')+SXPAR(h,'TARG_ID')+SXPAR(h,'SCOBS_ID')
+ root = STRCOMPRESS(root, /REMOVE_ALL)
+ targ = STRCOMPRESS(SXPAR(h,'TARGNAME'), /REMOVE_ALL)
+ objclass = SXPAR(h,'OBJCLASS')
+ sp_type = STRCOMPRESS(SXPAR(h,'SP_TYPE'), /REMOVE_ALL)
+ src_type = STRCOMPRESS(SXPAR(h,'SRC_TYPE'), /REMOVE_ALL)
+ mode = STRLOWCASE(STRCOMPRESS(SXPAR(h,'INSTMODE'), /REMOVE_ALL))
+ aper = STRCOMPRESS(SXPAR(h,'APERTURE'), /REMOVE_ALL)
+ vers = STRCOMPRESS(SXPAR(h,'CF_VERS'), /REMOVE_ALL)
+ ext = '4'
+ IF aper EQ 'MDRS' THEN ext = '2'
+ IF aper EQ 'HIRS' THEN ext = '3'
+ if keyword_set(airglow) then title = ''+root+' [AIRGLOW]' else $
+ title = ''+root+' ['+targ+']'
+ if (src_type eq 'PE') then yscale = 1.0 else yscale = 1.3
+ if (objclass eq 7 or sp_type eq 'BKGD') then bkgd_obs = 1 else bkgd_obs = 0
+
+;
+; Read each channel and its header.
+;
+ L1A = mrdfits(lun,0,HL1A)
+ L1B = mrdfits(lun,0,HL1B)
+ L2B = mrdfits(lun,0,HL2B)
+ L2A = mrdfits(lun,0,HL2A)
+ S1A = mrdfits(lun,0,HS1A)
+ S1B = mrdfits(lun,0,HS1B)
+ S2B = mrdfits(lun,0,HS2B)
+ S2A = mrdfits(lun,0,HS2A)
+
+;
+; Use different sample regions for point and extended sources.
+;
+ S1B_SAMPLE = [912, 935, 955, 970, 980, 985]
+ S2A_SAMPLE = [912, 935, 955, 970, 980, 985]
+
+ if (src_type eq 'PC') then begin
+ L1A_SAMPLE = [0, 0, 0, 0, 1045, 1070]
+ S1A_SAMPLE = [0, 0, 0, 0, 1045, 1070]
+ L2B_SAMPLE = [0, 0, 0, 0, 1045, 1070]
+ S2B_SAMPLE = [0, 0, 0, 0, 1045, 1070]
+ endif else begin
+ L1A_SAMPLE = [0, 0, 1030, 1039, 1045, 1070]
+ S1A_SAMPLE = [0, 0, 1030, 1039, 1045, 1070]
+ L2B_SAMPLE = [0, 0, 1030, 1039, 1045, 1070]
+ S2B_SAMPLE = [0, 0, 1030, 1039, 1045, 1070]
+ endelse
+
+ L1B_SAMPLE = [1095, 1130, 1140, 1165, 1170, 1190]
+ L2A_SAMPLE = [1095, 1130, 1140, 1165, 1170, 1190]
+
+;
+; Plot individual channels. Skip this for airglow files.
+;
+ ;set_plot,'X'
+ ;CLEANPLOT
+ set_plot,'Z'
+ device,set_resolution=[650,580]
+
+ !x.style = 1 & !y.style = 1
+ !x.ticklen=0.045
+ !x.charsize=1.5
+ !y.charsize=1.5
+ !y.minor = 2
+ !p.title = ' '
+ !x.title = ' '
+ !y.title = ' '
+ !p.color = 0
+ !p.background = 255
+
+;
+; LiF 1A and 1B
+;
+ !p.position = [0.13,0.57,0.95,0.95]
+ TOP_PLOT, L1A, HL1A, [985,1085], yscale, aper, L1A_SAMPLE, FL1A, title
+
+ !p.position = [0.13,0.1,0.95,0.48]
+ BOTTOM_PLOT, L1B, HL1B, [1092,1190], yscale, aper, L1B_SAMPLE, FL1B, vers
+
+ if (ver ge 5.4 and ver le 6.0) then write_jpeg,root+fnumber+'lif1'+mode+'f.jpg',TVRD() $
+ else write_gif,root+fnumber+'lif1'+mode+'f.gif',TVRD()
+
+;
+; SiC 1A and 1B
+;
+ !p.position = [0.13,0.57,0.95,0.95]
+ TOP_PLOT, S1A, HS1A, [1000,1094], yscale, aper, S1A_SAMPLE, FS1A, title
+
+ !p.position = [0.13,0.1,0.95,0.48]
+ BOTTOM_PLOT, S1B, HS1B, [910,995], yscale, aper, S1B_SAMPLE, FS1B, vers
+
+ if (ver ge 5.4 and ver le 6.0) then write_jpeg,root+fnumber+'sic1'+mode+'f.jpg',TVRD() $
+ else write_gif,root+fnumber+'sic1'+mode+'f.gif',TVRD()
+
+;
+; LiF 2A and 2B
+;
+ !p.position = [0.13,0.57,0.95,0.95]
+ TOP_PLOT, L2A, HL2A, [1085,1185], yscale, aper, L2A_SAMPLE, FL2A, title
+
+ !p.position = [0.13,0.1,0.95,0.48]
+ BOTTOM_PLOT, L2B, HL2B, [975,1078], yscale, aper, L2B_SAMPLE, FL2B, vers
+
+ if (ver ge 5.4 and ver le 6.0) then write_jpeg,root+fnumber+'lif2'+mode+'f.jpg',TVRD() $
+ else write_gif,root+fnumber+'lif2'+mode+'f.gif',TVRD()
+
+;
+; SiC 2A and 2B
+;
+ !p.position = [0.13,0.57,0.95,0.95]
+ TOP_PLOT, S2A, HS2A, [913,1009], yscale, aper, S2A_SAMPLE, FS2A, title
+
+ !p.position = [0.13,0.1,0.95,0.48]
+ BOTTOM_PLOT, S2B, HS2B, [1013,1106], yscale, aper, S2B_SAMPLE, FS2B, vers
+
+ if (ver ge 5.4 and ver le 6.0) then write_jpeg,root+fnumber+'sic2'+mode+'f.jpg',TVRD() $
+ else write_gif,root+fnumber+'sic2'+mode+'f.gif',TVRD()
+
+;
+; Now plot combined fluxes
+;
+ ;set_plot,'X'
+ ;CLEANPLOT
+ set_plot,'Z'
+ device,set_resolution=[650,870]
+
+ !x.style = 1 & !y.style = 1
+ !x.ticklen=0.045
+ !x.charsize=1.5
+ !y.charsize=1.5
+ !y.minor = 2
+ !p.title = ' '
+ !x.title = ' '
+ !y.title = ' '
+ !p.color = 0
+ !p.background = 255
+
+;
+; 900 - 1000 A
+;
+ !p.position = [0.13,0.70,0.95,0.95]
+ obs2a = SXPAR(HS2A, 'OBSTIME')
+ obs1b = SXPAR(HS1B, 'OBSTIME')
+
+ if ((bkgd_obs eq 1 and FS2A gt -1) or (FS2A gt 0.9*FS1B) or (obs2a gt 100 and obs1b lt 10)) then begin
+ TOP_PLOT, S2A, HS2A, [913,1009], yscale, aper, S2A_SAMPLE, foo, title
+ endif else begin
+ TOP_PLOT, S1B, HS1B, [910,995], yscale, aper, S1B_SAMPLE, foo, title
+ endelse
+
+;
+; 1000 - 1100 A
+;
+ !p.noerase = 1
+ !p.position = [0.13,0.40,0.95,0.65]
+ obl1a = SXPAR(HL1A, 'OBSTIME')
+ obl2b = SXPAR(HL2B, 'OBSTIME')
+
+ if ((bkgd_obs eq 1 and FL1A gt -1) or $
+ ((FL1A gt 0.9*FL2B and FL1A gt 0.7*FS1A and FL1A gt 0.7*FS2B) $
+ and not (obl2b gt 100 and obl1a lt 10))) then $
+ TOP_PLOT, L1A, HL1A, [985,1085], yscale, aper, L1A_SAMPLE, foo, '' else $
+ if ((bkgd_obs eq 1 and FL2B gt -1) or $
+ ((FL2B gt 0.7*FS1A and FL2B gt 0.7*FS2B) and not (obl2b lt 10 and obl1a gt 100))) then $
+ TOP_PLOT, L2B, HL2B, [975,1078], yscale, aper, L2B_SAMPLE, foo, '' else $
+ if ((bkgd_obs eq 1 and FS1A gt -1) or ((FS1A gt FS2B) and not (obl1a lt 10 and obl2b gt 100))) then $
+ TOP_PLOT, S1A, HS1A, [1000,1094], yscale, aper, S1A_SAMPLE, foo, '' else $
+ TOP_PLOT, S2B, HS2B, [1013,1106], yscale, aper, S2B_SAMPLE, foo, ''
+
+;
+; 1100 - 1200 A
+;
+ !p.position = [0.13,0.1,0.95,0.35]
+ obl2a = SXPAR(HL2A, 'OBSTIME')
+ obl1b = SXPAR(HL1B, 'OBSTIME')
+
+ if ((bkgd_obs eq 1 and FL2A gt -1) or (FL2A gt 0.9*FL1B) or (obl2a gt 100 and obl1b lt 10)) then $
+ BOTTOM_PLOT, L2A, HL2A, [1085,1182], yscale, aper, L2A_SAMPLE, foo, vers else $
+ BOTTOM_PLOT, L1B, HL1B, [1092,1190], yscale, aper, L1B_SAMPLE, foo, vers
+
+ if (ver ge 5.4 and ver le 6.0) then write_jpeg,root+fnumber+'spec'+mode+'f.jpg',TVRD() $
+ else write_gif,root+fnumber+'spec'+mode+'f.gif',TVRD()
+
+ RETURN
+END
+
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
diff --git a/idl/cf_plot_rate3.pro b/idl/cf_plot_rate3.pro
new file mode 100644
index 0000000..6a2669b
--- /dev/null
+++ b/idl/cf_plot_rate3.pro
@@ -0,0 +1,227 @@
+;+
+; cf_plot_rate3.pro is a procedure to plot the count rate and screening
+; as a function of time for calfuse 3.0 data.
+;
+; Author: Edward M. Murphy
+; Written: 2000 April 21
+;
+; Calling sequence: cf_plot_rate3, rootname
+;
+; Inputs: rootname will have 'idf.fit' added to it.
+;
+; 02/20/03 JC Hsu Adapt from cf_plot_rate.pro to work on version 3.0 data.
+; 06/23/04 V Dixon Mark low-voltage periods with cross-hatch.
+; 03/16/05 V Dixon Ignore timeline table entries that extend beyond
+; the end of the exposure.
+; 08/01/05 V Dixon If IDL version > 6.0, generate GIF file.
+; 11/03/05 V Dixon Shade OPUS and SAA bad times before drawing rest of plot.
+;-
+
+pro cf_plot_rate3,rootname
+
+ !quiet=1
+
+ ; Open input IDF file
+ filename=rootname+'idf.fit'
+ dummy=mrdfits(filename,0,ah,/SILENT)
+
+ ; Read primary header keywords
+ pdate = fxpar(ah,'DATE')
+ dateobs = fxpar(ah,'DATEOBS')
+ expstart = fxpar(ah,'EXPSTART')
+ expend = fxpar(ah,'EXPEND')
+ detector = fxpar(ah,'DETECTOR')
+
+ exposure = (expend - expstart) * 24. * 3600.
+
+ ; Define plot attributes.
+ !p.multi=[0,1,2]
+ set_plot,'Z'
+ device,set_resolution=[1024,512]
+
+ !x.style=1
+ !y.style=1
+ !x.tickformat='(I6)'
+
+ !x.title="Time (seconds)"
+ !y.title="Count rate"
+ !p.charsize=1.2
+ !p.background=255
+ !p.color=0
+
+ ; Read columns from the timeline extension (ext 3)
+ ftab_ext, filename, 'TIME', time0, exten_no=3
+ ftab_ext, filename, 'STATUS_FLAGS', flags, exten_no=3
+ if (exposure gt 55000) then begin
+ good = where(time0 lt 55000, n)
+ exposure = time0[n-1] - time0[0]
+ endif
+ good = where(time0 lt exposure)
+ time = [0, time0[good], exposure]
+
+ ; Compute the day/night screenings (bit 1)
+ flags1 = ishft(flags, 7)
+ flags1 = ishft(flags1, -7) ; only keep bit 1
+
+ ; Compute the limb angle screening (bit 2)
+ flags2 = ishft(flags, 6) ;
+ flags2 = ishft(flags2, -7) ; only keep bit 2
+
+ ; Compute the SAA screening (bit 3)
+ flags3 = ishft(flags, 5) ;
+ flags3 = ishft(flags3, -7) ; only keep bit 3
+
+ ; Compute the HV screening (bit 4)
+ flags4 = ishft(flags, 4) ;
+ flags4 = ishft(flags4, -7) ; only keep bit 4
+
+ ; Compute burst screening (bit 5)
+ flags5 = ishft(flags, 3) ;
+ flags5 = ishft(flags5, -7) ; only keep bit 5
+
+ ; Compute OPUS screening (bit 6)
+ flags6 = ishft(flags, 2) ;
+ flags6 = ishft(flags6, -7) ; only keep bit 6
+
+ ; Compute jitter screening (bit 7)
+ flags7 = ishft(flags, 1) ;
+ flags7 = ishft(flags7, -7) ; only keep bit 7
+
+ ftab_ext, filename, 'BKGD_CNT_RATE', bkgd, exten_no=3
+ len = size(time0)
+ packet = 1 > len[1]/100 ; determine rebinning size
+ wht = where(time0 eq time0)
+ whbin = where(wht mod packet eq 0)
+
+ for i = 1, 2 do begin
+
+ if (i eq 1) then text = ' LiF' else text = ' SiC'
+ !p.title=filename + text
+
+ ; Draw grid, shade bad times, then over-plot count rate arrays.
+ if (i eq 1) then col = 'LIF' else col = 'SIC'
+ ftab_ext, filename, col+'_CNT_RATE', crate, exten_no=3
+ ymin = 1
+ ymax = max([1000,crate * 2.])
+ crate = crate > ymin
+ !y.range=[ymin, ymax]
+ plot_io, time0, crate, /nodata, xr=[0,exposure]
+
+ ; Shade OPUS
+ flags16 = flags6 * (ymax-ymin) + ymin
+ flags16 = [0, flags16[good], 0]
+ polyfill, time, flags16, COL = 0.75 * !D.N_COLORS
+ ;oplot, time, flags16, psym=10, linestyle=2
+
+ ; Shade SAA violations
+ flags13 = flags3 * (ymax/100-ymin) + ymin
+ flags13 = [0., flags13[good], 0.]
+ polyfill, time, flags13, COL = 0.25 * !D.N_COLORS
+ flags13 = flags3 * (ymax-ymin) + ymin
+ oplot, time, flags13, psym=10, linestyle=2
+
+ ; Shade day/night screenings
+ flags11 = flags1 * (0.2*ymax) + ymax
+ flags11 = [ymax, flags11[good], ymax]
+ oplot, time, flags11, psym=10
+ polyfill, time, flags11
+
+ ; Shade limb angle
+ flags12 = flags2 * (ymax-ymin) + ymin
+ flags12 = [0, flags12[good], 0]
+ oplot, time, flags12, psym=10, linestyle=2
+ polyfill, time, flags12, /line_fill, orientation=0
+
+ ; Shade low-voltage periods.
+ flags14 = flags4 * (ymax-ymin) + ymin
+ flags14 = [0, flags14[good], 0]
+ oplot, time, flags14, psym=10, linestyle=2
+ polyfill, time, flags14, /line_fill, orientation=45
+ polyfill, time, flags14, /line_fill, orientation=135
+
+ ; Shade burst
+ flags15 = flags5 * (ymax-ymin) + ymin
+ flags15 = [0, flags15[good], 0]
+ oplot, time, flags15, psym=10, linestyle=2
+ polyfill, time, flags15, /line_fill, orientation=45
+
+ ; Shade jitter
+ flags17 = flags7 * (ymax-ymin) + ymin
+ flags17 = [0, flags17[good], 0]
+ oplot, time, flags17, psym=10, linestyle=2
+ polyfill, time, flags17, /line_fill, orientation=135
+
+ ; Overplot count-rate arrays.
+ oplot, time0, crate, linestyle=0, psym=10
+ oplot, time0(whbin), bkgd(whbin)+0.001, linestyle=1
+ endfor
+
+ ; Draw legend
+ plot,[0,1],xstyle=7,ystyle=7,/nodata,position=[345-20,5,345+21,25],$
+ /noerase,title=" ",xtitle=" ",ytitle=" ",xrange=[0,1],$
+ yrange=[0,1],/DEVICE
+ polyfill,[0,1,1,0],[0,0,1,1], COL = 0.75 * !D.N_COLORS
+ oplot,[0,1,1,0,0],[0,0,1,1,0],psym=10
+ xyouts,345-17,10,'OPUS',charsize=1,color=0,/DEVICE
+
+ plot,[0,1],xstyle=7,ystyle=7,/nodata,position=[400-20,5,400+21,25],$
+ /noerase,title=" ",xtitle=" ",ytitle=" ",xrange=[0,1],$
+ yrange=[0,1],/DEVICE
+ polyfill,[0,1,1,0],[0,0,0.36,0.36], COL = 0.25 * !D.N_COLORS
+ oplot,[0,1,1,0,0],[0,0,1,1,0],psym=10
+ xyouts,400-12,13,'SAA',charsize=1,color=0,/DEVICE
+
+ plot,[0,1],xstyle=7,ystyle=7,/nodata,position=[450-15,5,450+15,25],$
+ /noerase,title=" ",xtitle=" ",ytitle=" ",xrange=[0,1],$
+ yrange=[0,1],/DEVICE
+ polyfill,[0,1,1,0],[0.8,0.8,1,1]
+ xyouts,450-11,10,'Day',charsize=1,color=0,/DEVICE
+
+ plot,[0,1],xstyle=7,ystyle=7,/nodata,position=[710-15,5,710+15,25],$
+ /noerase,title=" ",xtitle=" ",ytitle=" ",xrange=[0,1],$
+ yrange=[0,1],/DEVICE
+ polyfill,[0,1,1,0],[0,0,1,1],/line_fill,orientation=0
+ oplot,[0,1,1,0,0],[0,0,1,1,0],psym=10
+ xyouts,710+20,10,'Limb',charsize=1,color=0,/DEVICE
+
+ plot,[0,1],xstyle=7,ystyle=7,/nodata,position=[790-15,5,790+15,25],$
+ /noerase,title=" ",xtitle=" ",ytitle=" ",xrange=[0,1],$
+ yrange=[0,1],/DEVICE
+ polyfill,[0,1,1,0],[0,0,1,1],/line_fill,orientation=45
+ oplot,[0,1,1,0,0],[0,0,1,1,0],psym=10
+ xyouts,790+20,10,'Burst',charsize=1,color=0,/DEVICE
+
+ plot,[0,1],xstyle=7,ystyle=7,/nodata,position=[870-15,5,870+15,25],$
+ /noerase,title=" ",xtitle=" ",ytitle=" ",xrange=[0,1],$
+ yrange=[0,1],/DEVICE
+ polyfill,[0,1,1,0],[0,0,1,1],/line_fill,orientation=135
+ oplot,[0,1,1,0,0],[0,0,1,1,0],psym=10
+ xyouts,870+20,10,'Jitter',charsize=1,color=0,/DEVICE
+
+ plot,[0,1],xstyle=7,ystyle=7,/nodata,position=[950-15,5,950+15,25],$
+ /noerase,title=" ",xtitle=" ",ytitle=" ",xrange=[0,1],$
+ yrange=[0,1],/DEVICE
+ polyfill,[0,1,1,0],[0,0,1,1],/line_fill,orientation=45
+ polyfill,[0,1,1,0],[0,0,1,1],/line_fill,orientation=135
+ oplot,[0,1,1,0,0],[0,0,1,1,0],psym=10
+ xyouts,950+20,10,'HV',charsize=1,color=0,/DEVICE
+
+ xyouts,835,27,'Screenings :',charsize=.9,color=0,/DEVICE,alignment=.5
+ xyouts,20,16,'RAW file written at '+pdate,charsize=0.90,color=0,$
+ /DEVICE
+ xyouts,925,495,'DATEOBS: '+dateobs,charsize=.9,color=0,$
+ /DEVICE,alignment=.5
+
+ ; Create a GIF/JPEG file:
+ ver = float(!version.release)
+ if (ver ge 5.4 and ver le 6.0) then begin
+ xyouts,20,6,'JPG file written at '+!stime,charsize=0.90,$
+ color=0,/DEVICE
+ write_jpeg,rootname+'rat.jpg',TVRD()
+ endif else begin
+ xyouts,20,6,'GIF file written at '+!stime,charsize=0.90,$
+ color=0,/DEVICE
+ write_gif,rootname+'rat.gif',TVRD()
+ endelse
+return
+end
diff --git a/idl/cf_plot_rate_linear.pro b/idl/cf_plot_rate_linear.pro
new file mode 100644
index 0000000..e2e3aa7
--- /dev/null
+++ b/idl/cf_plot_rate_linear.pro
@@ -0,0 +1,175 @@
+;+
+; cf_plot_rate_linear.pro is a procedure to plot the count rate and screening
+; as a function of time for calfuse 3.0 data.
+; THIS VERSION PRODUCES A LINEAR PLOT OF COUNT RATE
+;
+; Author: Edward M. Murphy
+; Written: 2000 April 21
+;
+; Calling sequence: cf_plot_rate3, rootname
+;
+; Inputs: rootname will have 'idf.fit' added to it.
+;
+; 02/20/03 JC Hsu Adapt from cf_plot_rate.pro to work on version 3.0 data.
+; 05/17/04 V Dixon Adapt from cf_plot_rate3.pro
+;-
+
+pro cf_plot_rate_linear, rootname
+
+ !quiet=1
+
+ ; Open input IDF file
+ filename=rootname+'idf.fit'
+ dummy=mrdfits(filename,0,ah,/SILENT)
+
+ ; Read primary header keywords
+ pdate = fxpar(ah,'DATE')
+ dateobs = fxpar(ah,'DATEOBS')
+ expstart = fxpar(ah,'EXPSTART')
+ expend = fxpar(ah,'EXPEND')
+ detector = fxpar(ah,'DETECTOR')
+
+ exposure = (expend - expstart) * 24. * 3600.
+
+ ; Define plot attributes.
+ !p.multi=[0,1,2]
+ set_plot,'Z'
+ device,set_resolution=[1024,512]
+
+ !x.style=1
+ !y.style=1
+ !x.tickformat='(I6)'
+
+ !x.title="Time (seconds)"
+ !y.title="Count rate"
+ !p.charsize=1.2
+ !p.background=255
+ !p.color=0
+
+ ; Read columns from the timeline extension (ext 3)
+ ftab_ext, filename, 'TIME', time0, exten_no=3
+ ftab_ext, filename, 'STATUS_FLAGS', flags, exten_no=3
+ time = [0, time0, exposure]
+
+ ; Compute the day/night screenings (bit 1)
+ flags1 = ishft(flags, 7)
+ flags1 = ishft(flags1, -7) ; only keep bit 1
+
+ ; Compute the limb angle screening (bit 2)
+ flags2 = ishft(flags, 6) ;
+ flags2 = ishft(flags2, -7) ; only keep bit 2
+
+ ; Compute the SAA screening (bit 3)
+ flags3 = ishft(flags, 5) ;
+ flags3 = ishft(flags3, -7) ; only keep bit 3
+
+ ; Compute burst screening (bit 5)
+ flags5 = ishft(flags, 3) ;
+ flags5 = ishft(flags5, -7) ; only keep bit 5
+
+ ; Compute jitter screening (bit 7)
+ flags7 = ishft(flags, 1) ;
+ flags7 = ishft(flags7, -7) ; only keep bit 7
+
+ ftab_ext, filename, 'BKGD_CNT_RATE', bkgd, exten_no=3
+ len = size(time0)
+ packet = 1 > len[1]/100 ; determine rebinning size
+ wht = where(time0 eq time0)
+ whbin = where(wht mod packet eq 0)
+
+ for i = 1, 2 do begin
+
+ if (i eq 1) then text = ' LiF' else text = ' SiC'
+ !p.title=filename + text
+
+ ; Plot the count rate
+ if (i eq 1) then col = 'LIF' else col = 'SIC'
+ ftab_ext, filename, col+'_CNT_RATE', crate, exten_no=3
+ ymax = max(crate)*1.2
+ ymin = 0.
+ !y.range=[ymin, ymax]
+ plot, time0, crate, xr=[0,exposure], linestyle=0, psym=10
+ oplot, time0(whbin), bkgd(whbin)*(10./i), linestyle=1
+ if (i eq 1) then text = '5' else text = '10'
+ xyouts,exposure/50,ymax/20,'BKGD scaled by '+text,charsize=1
+
+
+ ; Shade day/night screenings
+ ;;;flags11 = flags1 * (-0.2*ymin) + ymin
+ ;;;flags11 = [0, flags11, 0]
+ flags11 = flags1 * (0.025*ymax) + ymax
+ flags11 = [ymax, flags11, ymax]
+ oplot, time, flags11, psym=10
+ polyfill, time, flags11
+
+ ; Shade limb angle
+ flags12 = flags2 * (ymax-ymin) + ymin
+ flags12 = [ymin, flags12, ymin]
+ oplot, time, flags12, psym=10, linestyle=2
+ polyfill, time, flags12, /line_fill, orientation=0
+
+ ; Shade SAA
+ ;;;flags13 = flags3 * (ymax-ymin) + ymin
+ ;;;flags13 = [ymin, flags13, ymin]
+ ;;;oplot, time, flags13, psym=10, linestyle=2
+
+ ; Shade burst
+ flags15 = flags5 * (ymax-ymin) + ymin
+ flags15 = [ymin, flags15, ymin]
+ oplot, time, flags15, psym=10, linestyle=2
+ polyfill, time, flags15, /line_fill, orientation=45
+
+ ; Shade jitter
+ flags17 = flags7 * (ymax-ymin) + ymin
+ flags17 = [ymin, flags17, ymin]
+ oplot, time, flags17, psym=10, linestyle=2
+ polyfill, time, flags17, /line_fill, orientation=135
+ endfor
+
+ ; Draw legend
+ plot,[0,1],xstyle=7,ystyle=7,/nodata,position=[710-15,5,710+15,25],$
+ /noerase,title=" ",xtitle=" ",ytitle=" ",xrange=[0,1],$
+ yrange=[0,1],/DEVICE
+ polyfill,[0,1,1,0],[0.8,0.8,1,1]
+ xyouts,710+20,10,'Day',charsize=1,color=0,/DEVICE
+
+ plot,[0,1],xstyle=7,ystyle=7,/nodata,position=[790-15,5,790+15,25],$
+ /noerase,title=" ",xtitle=" ",ytitle=" ",xrange=[0,1],$
+ yrange=[0,1],/DEVICE
+ polyfill,[0,1,1,0],[0,0,1,1],/line_fill,orientation=0
+ oplot,[0,1,1,0,0],[0,0,1,1,0],psym=10
+ xyouts,790+20,10,'Limb',charsize=1,color=0,/DEVICE
+
+ plot,[0,1],xstyle=7,ystyle=7,/nodata,position=[870-15,5,870+15,25],$
+ /noerase,title=" ",xtitle=" ",ytitle=" ",xrange=[0,1],$
+ yrange=[0,1],/DEVICE
+ polyfill,[0,1,1,0],[0,0,1,1],/line_fill,orientation=45
+ oplot,[0,1,1,0,0],[0,0,1,1,0],psym=10
+ xyouts,870+20,10,'Burst',charsize=1,color=0,/DEVICE
+
+ plot,[0,1],xstyle=7,ystyle=7,/nodata,position=[950-15,5,950+15,25],$
+ /noerase,title=" ",xtitle=" ",ytitle=" ",xrange=[0,1],$
+ yrange=[0,1],/DEVICE
+ polyfill,[0,1,1,0],[0,0,1,1],/line_fill,orientation=135
+ oplot,[0,1,1,0,0],[0,0,1,1,0],psym=10
+ xyouts,950+20,10,'Jitter',charsize=1,color=0,/DEVICE
+
+ xyouts,835,27,'Screenings :',charsize=.9,color=0,/DEVICE,alignment=.5
+ xyouts,20,16,'RAW file written at '+pdate,charsize=0.90,color=0,$
+ /DEVICE
+ xyouts,925,495,'DATEOBS: '+dateobs,charsize=.9,color=0,$
+ /DEVICE,alignment=.5
+
+ ; Create a GIF/JPEG file:
+ ver = float(!version.release)
+ if (ver ge 5.4) then begin
+ xyouts,20,6,'JPG file written at '+!stime,charsize=0.90,$
+ color=0,/DEVICE
+ write_jpeg,rootname+'lin.jpg',TVRD()
+ endif else begin
+ xyouts,20,6,'GIF file written at '+!stime,charsize=0.90,$
+ color=0,/DEVICE
+ write_gif,rootname+'lin.gif',TVRD()
+ endelse
+return
+end
diff --git a/idl/fuse_extract.pro b/idl/fuse_extract.pro
new file mode 100644
index 0000000..0b784d7
--- /dev/null
+++ b/idl/fuse_extract.pro
@@ -0,0 +1,154 @@
+pro fuse_extract, input
+;+
+; NAME:
+; FUSE_EXTRACT
+;
+;*PURPOSE:
+; To extract full-resolution spectra from the primary LiF and SiC
+; apertures of an intermediate data file.
+;
+;*CATEGORY:
+; INPUT/OUTPUT
+;
+;*CALLING SEQUENCE:
+; FUSE_EXTRACT, input
+;
+;*INPUT:
+; input: intermediate data file
+;
+;*OUTPUT:
+; LiF and SiC spectral files in v2.x format.
+;
+;
+;*PROCEDURES USED:
+; FXADDPAR, FXHMAKE, FXBADDCOL, FXBCREATE, FXBWRITE, FXBFINISH, MRDFITS,
+; SXPAR
+;
+;*HISTORY:
+; Written by: V. Dixon March 2003
+; 03/23/2003 V. Dixon Correct X values for Doppler shift.
+; 03/24/2003 V. Dixon Fix bug in call to where().
+; 06/27/2003 R. Robinson Update to new version of IDF
+; - added fscale keyword to mrdfits call
+; - divide flux by exposure time
+; 06/08/2004 V. Dixon Use EXPTIME rather than RAWTIME.
+; Correct calculation of COUNTS
+; 07/01/2004 V. Dixon Don't read headers of extensions.
+; Employ faster scheme for calculating
+; Doppler shift.
+; 07/02/2004 V. Dixon Don't compute dl per pixel. It suffers
+; from truncation errors. Instead, use
+; value from file header.
+;-
+;-----------------------------------------------------------------------------
+;
+; print calling sequence if no parameters supplied
+;
+ if n_params(0) lt 1 then begin
+ print,'CALLING SEQUENCE: fuse_extract, input'
+ return
+ end
+;
+; Physical constants
+;
+ C = 2.99792458e5 ; km/s
+;
+; Open input file and read keywords from header.
+;
+ dummy = mrdfits(input, 0, hdr, status=status)
+ rootname = SXPAR(hdr, 'ROOTNAME')
+ vhelio = SXPAR(hdr, 'V_HELIO')
+ detector = SXPAR(hdr, 'DETECTOR')
+ aperture = SXPAR(hdr, 'APERTURE')
+ wavecal = SXPAR(hdr, 'WAVE_CAL')
+ exptime = SXPAR(hdr, 'EXPTIME')
+
+ ap_array = ['HIRS', 'MDRS', 'LWRS', 'PINH']
+ ap_code = ['3','2','4','1']
+ for k = 0, 3 do if (strcmp(aperture, ap_array[k], 1)) then ap = k
+ channel = ap+1
+
+ GET_DATE, FITS_date
+;
+; Read photon information in the first extension of the IDF.
+;
+ a = mrdfits(input, 1, status=status, /fscale)
+;
+; Read timeline table from the third extension of the IDF.
+;
+ b = mrdfits(input, 3, status=status, /fscale)
+ nseconds = n_elements(b.time)
+ if (nint(b.time[nseconds-1] - b.time[0]) gt nseconds) then $
+ print, 'WARNING: '+input+' timeline table not continuous.'
+
+ for k = channel, channel+4, 4 do begin
+ print, 'Analyzing extension ', k
+ ;
+ ; Read corresponding wavelength array from WAVE_CAL file.
+ ;
+ wcalfile=strtrim('/data1/fuse/calfuse/v3.0/calfiles/'+wavecal,2)
+ print, 'Reading wavelength data from ', wcalfile
+ w=mrdfits(wcalfile,k,w_hdr)
+ wave = w.wavelength
+ wpc = SXPAR(w_hdr, 'DISPAPIX')
+ ;
+ ; Select photons from aperture of interest.
+ ;
+ g = where (a.channel eq k and (a.timeflgs and not 1B) eq 0 and $
+ (a.loc_flgs and 16B) eq 0, n)
+ print, 'Number of photons in spectrum = ', n
+ x = a.x[g]
+ ;
+ ; Correct X values for Doppler shift.
+ ;
+ print, 'Applying Doppler shift'
+ t = nint(a.time[g] - b.time[0])
+ ovc = (b.ORBITAL_VEL + vhelio)/C
+ for i = 0L, n-1L do x[i] = x[i] + ovc[t[i]] * a.lambda[g[i]]/wpc
+ ;
+ ; Extract spectrum and populate output arrays.
+ ;
+ print, 'Extracting spectrum'
+ x = nint(x)
+ flux = fltarr(16384)
+ error = fltarr(16384)
+ cntserr = fltarr(16384)
+ for i = 0L, n-1L do flux[x[i]] = flux[x[i]] + a.ergcm2[g[i]]
+ flux = flux / exptime / abs(wpc)
+ counts = histogram(x, min=0, max=16383)
+ g = where (counts gt 0)
+ cntserr[g] = sqrt(counts[g])
+ error[g] = flux[g]/cntserr[g]
+ quality = indgen(16384)
+ ;
+ ; Create primary HDU of output file
+ ;
+ if (k lt 5) then ch = 'lif' else ch = 'sic'
+ filename = strcompress(rootname+strlowcase(detector)+ch+ap_code[ap]+'ttagfcal.fit',/remove_all)
+ FXADDPAR, hdr, 'DATE', FITS_date
+ FXADDPAR, hdr, 'FILENAME', filename
+ FXADDPAR, hdr, 'FILETYPE', 'CALIBRATED EXTRACTED SPECTRUM'
+ FXADDPAR, hdr, 'APER_ACT', strupcase(strcompress(ap_array[ap]+'_'+ch,/remove_all))
+ FXHMAKE, hdr, /EXTEND
+ FXWRITE, filename, hdr
+
+ FXBHMAKE, HEADER, 1, 'FUSE 1D Spectrum'
+ FXBADDCOL, iwave, HEADER, wave, 'WAVE', TUNIT='ANGSTROMS'
+ FXBADDCOL, iflux, HEADER, flux, 'FLUX', TUNIT='ERG/CM2/S/A'
+ FXBADDCOL, ierror, HEADER, error, 'ERROR', TUNIT='ERG/CM2/S/A'
+ FXBADDCOL, iquality, HEADER, quality, 'QUALITY', TUNIT='UNITLESS'
+ FXBADDCOL, icounts, HEADER, counts, 'COUNTS', TUNIT='COUNTS'
+ FXBADDCOL, icntserr, HEADER, cntserr, 'CNTSERR', TUNIT='COUNTS'
+
+ FXBCREATE, out, filename, HEADER
+ FXBWRITE, out, wave, iwave, 1
+ FXBWRITE, out, flux, iflux, 1
+ FXBWRITE, out, error, ierror, 1
+ FXBWRITE, out, quality, iquality, 1
+ FXBWRITE, out, counts, icounts, 1
+ FXBWRITE, out, cntserr, icntserr, 1
+ FXBFINISH, out
+ endfor
+ close, /all
+ return
+ end
diff --git a/idl/fuse_sic_only.pro b/idl/fuse_sic_only.pro
new file mode 100644
index 0000000..f0eadc0
--- /dev/null
+++ b/idl/fuse_sic_only.pro
@@ -0,0 +1,159 @@
+pro fuse_sic_only, input
+;+
+; NAME:
+; FUSE_EXTRACT
+;
+;*PURPOSE:
+; To extract full-resolution spectra from the primary LiF and SiC
+; apertures of an intermediate data file.
+;
+; THIS VERSION EXTRACTS ONLY A SIC SPECTRUM.
+;
+;*CATEGORY:
+; INPUT/OUTPUT
+;
+;*CALLING SEQUENCE:
+; FUSE_EXTRACT, input
+;
+;*INPUT:
+; input: intermediate data file
+;
+;*OUTPUT:
+; LiF and SiC spectral files in v2.x format.
+;
+;
+;*PROCEDURES USED:
+; FXADDPAR, FXHMAKE, FXBADDCOL, FXBCREATE, FXBWRITE, FXBFINISH, MRDFITS,
+; SXPAR
+;
+;*HISTORY:
+; Written by: V. Dixon March 2003
+; 03/23/2003 V. Dixon Correct X values for Doppler shift.
+; 03/24/2003 V. Dixon Fix bug in call to where().
+; 06/27/2003 R. Robinson Update to new version of IDF
+; - added fscale keyword to mrdfits call
+; - divide flux by exposure time
+; 06/08/2004 V. Dixon Use EXPTIME rather than RAWTIME.
+; Correct calculation of COUNTS
+; 07/01/2004 V. Dixon Don't read headers of extensions.
+; Employ faster scheme for calculating
+; Doppler shift.
+; 07/02/2004 V. Dixon Don't compute dl per pixel. It suffers
+; from truncation errors. Instead, use
+; value from file header.
+; 07/15/2004 V. Dixon Extract only SiC spectrum.
+;-
+;-----------------------------------------------------------------------------
+;
+; print calling sequence if no parameters supplied
+;
+ if n_params(0) lt 1 then begin
+ print,'CALLING SEQUENCE: fuse_sic_only, input'
+ return
+ end
+;
+; Physical constants
+;
+ C = 2.99792458e5 ; km/s
+;
+; Open input file and read keywords from header.
+;
+ dummy = mrdfits(input, 0, hdr, status=status)
+ rootname = SXPAR(hdr, 'ROOTNAME')
+ vhelio = SXPAR(hdr, 'V_HELIO')
+ detector = SXPAR(hdr, 'DETECTOR')
+ aperture = SXPAR(hdr, 'APERTURE')
+ wavecal = SXPAR(hdr, 'WAVE_CAL')
+ exptime = SXPAR(hdr, 'EXPTIME')
+
+ ap_array = ['HIRS', 'MDRS', 'LWRS', 'PINH']
+ ap_code = ['3','2','4','1']
+ for k = 0, 3 do if (strcmp(aperture, ap_array[k], 1)) then ap = k
+ channel = ap+1
+
+ GET_DATE, FITS_date
+;
+; Read photon information in the first extension of the IDF.
+;
+ a = mrdfits(input, 1, status=status, /fscale)
+;
+; Read timeline table from the third extension of the IDF.
+;
+ b = mrdfits(input, 3, status=status, /fscale)
+ nseconds = n_elements(b.time)
+ if (nint(b.time[nseconds-1] - b.time[0]) gt nseconds) then $
+ print, 'WARNING: '+input+' timeline table not continuous.'
+
+ ; for k = channel, channel+4, 4 do begin
+ k = channel + 4
+ print, 'Analyzing extension ', k
+ ;
+ ; Read corresponding wavelength array from WAVE_CAL file.
+ ;
+ wcalfile=strtrim('/data1/fuse/calfuse/v3.0/calfiles/'+wavecal,2)
+ print, 'Reading wavelength data from ', wcalfile
+ w=mrdfits(wcalfile,k,w_hdr)
+ wave = w.wavelength
+ wpc = SXPAR(w_hdr, 'DISPAPIX')
+ ;
+ ; Select photons from aperture of interest.
+ ;
+ g = where (a.channel eq k and (a.timeflgs and not 1B) eq 0 and $
+ (a.loc_flgs and 16B) eq 0, n)
+ print, 'Number of photons in spectrum = ', n
+ x = a.x[g]
+ ;
+ ; Correct X values for Doppler shift.
+ ;
+ print, 'Applying Doppler shift'
+ t = nint(a.time[g] - b.time[0])
+ ovc = (b.ORBITAL_VEL + vhelio)/C
+ for i = 0L, n-1L do x[i] = x[i] + ovc[t[i]] * a.lambda[g[i]]/wpc
+ ;
+ ; Extract spectrum and populate output arrays.
+ ;
+ print, 'Extracting spectrum'
+ x = nint(x)
+ flux = fltarr(16384)
+ error = fltarr(16384)
+ cntserr = fltarr(16384)
+ for i = 0L, n-1L do flux[x[i]] = flux[x[i]] + a.ergcm2[g[i]]
+ flux = flux / exptime / abs(wpc)
+ counts = histogram(x, min=0, max=16383)
+ g = where (counts gt 0)
+ cntserr[g] = sqrt(counts[g])
+ error[g] = flux[g]/cntserr[g]
+ quality = indgen(16384)
+ ;
+ ; Create primary HDU of output file
+ ;
+ if (k lt 5) then ch = 'lif' else ch = 'sic'
+ filename = strcompress(rootname+strlowcase(detector)+ch+ap_code[ap]+'ttagfcal.fit',/remove_all)
+ FXADDPAR, hdr, 'DATE', FITS_date
+ FXADDPAR, hdr, 'FILENAME', filename
+ FXADDPAR, hdr, 'FILETYPE', 'CALIBRATED EXTRACTED SPECTRUM'
+ FXADDPAR, hdr, 'APER_ACT', strupcase(strcompress(ap_array[ap]+'_'+ch,/remove_all))
+ FXHMAKE, hdr, /EXTEND
+ FXWRITE, filename, hdr
+
+ FXBHMAKE, HEADER, 1, 'FUSE 1D Spectrum'
+ FXBADDCOL, iwave, HEADER, wave, 'WAVE', TUNIT='ANGSTROMS'
+ FXBADDCOL, iflux, HEADER, flux, 'FLUX', TUNIT='ERG/CM2/S/A'
+ FXBADDCOL, ierror, HEADER, error, 'ERROR', TUNIT='ERG/CM2/S/A'
+ FXBADDCOL, iquality, HEADER, quality, 'QUALITY', TUNIT='UNITLESS'
+ FXBADDCOL, icounts, HEADER, counts, 'COUNTS', TUNIT='COUNTS'
+ FXBADDCOL, icntserr, HEADER, cntserr, 'CNTSERR', TUNIT='COUNTS'
+
+ FXBCREATE, out, filename, HEADER
+ FXBWRITE, out, wave, iwave, 1
+ FXBWRITE, out, flux, iflux, 1
+ FXBWRITE, out, error, ierror, 1
+ FXBWRITE, out, quality, iquality, 1
+ FXBWRITE, out, counts, icounts, 1
+ FXBWRITE, out, cntserr, icntserr, 1
+ FXBFINISH, out
+ ;endfor
+
+ close, /all
+ return
+ end
generated by cgit (git 2.34.1) at 2025-09-20 19:51:35 -0400