Viewing contents of file '../idllib/ghrs/pro/arc.pro'
pro arc,name,gnumber,textout=textout,calspec=calspec,ident=ident,noplot=noplot
;+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
;+
;
;*NAME: ARC
;
;*CLASS:
;
;*CATEGORY: WAVELENGTH CALIBRATION
;
;*PURPOSE:
; A driver for the GHRS WAVECAL procedure to give a more interactive
; representation for the wavelength calibration.
;
;*CALLING SEQUENCE:
; arc,name[,gnumber,
; textout=textout,calspec=calspec,ident=ident,noplot=noplot]
;
;*PARAMETERS:
; name - file name of GHRS WAVECAL observation.
; gnumber - (INT) grating number (range 1-7)
; 1 = G140M
; 2 = G160M
; 3 = G200M
; 4 = G270M
; 5 = G140L
; 6 = Echelle A
; 7 = Echelle B
;
;*Optional Keywords:
; calspec - (logical) - 0 if CALHRS has not been run on WAVECAL obs.
; - 1 if CALHRS has already been run on WAVECAL obs.
; textout - redirects output of log file (i.e. textout=3 creates file
; ARC.PRT). Default is to the screen.
; ident - (string) - identifier for calibrated WAVECAL,(default: _R)
; noplot - (logical) - 0 - perform interactive graphics steps.
; 1 - DO NOT perform interactive graphics steps.
;
;*EXAMPLES:
; 1) Run ARC on archival WAVECAL obs. 3607 (PKS2155-304 WAVCAL/SPYBAL)
; arc,3607
; ...you will be prompted for inputs and plot scaling
; 2) Same as example 1 except we run initial CALHRS prior to calling ARC
; calhrs,3607,'bck=00,errors=0' ; creates h3607_r
; arc,3607,/calspec
; 3) Same as examples 2 except change IDENTifier for calibrated file.
; calhrs,3607,'bck=00,errors=0,ident=_TST ; creates h3607_tst
; arc,3607,/calspec,ident='_tst'
; 4) Same as example 1, except turn off all plotting (non-interactive mode)
; arc,3607,/calspec,/noplot
; 5) Same as example 2 except generate log of session:
; arc,3607,textout=3 ; creates ARC.PRT and ARC_RESULTS.PRT
; 6) Run ARC in full debug mode
; !debug=1 ; diagnostic text and graphics
; !dump=2 ; complete history
; arc,3607,textout=3 ; creates ARC.PRT and ARC_RESULTS.PRT
; 7) Capture TEK mode plots for hardcopy
; plotopen
; arc,3607
; plotprint
;
;*NOTES:
; For grating 5, it uses NBSLIB_G5.TAB.
; For all other gratings, it uses NBSLIB.TAB.
;
; This routine assumes only 1 spectra order for echelle modes.
;
;*MODIFICATION HISTORY:
; R Robinson/CSC - Version 1
; 18-apr-1992 JKF/ACC - moved to GHRS DAF.
; 18-may-1992 JKF/ACC - switched defaults for YESNO.
; 20-jul-1992 JKF/ACC - added NVAR=2 for echelle mode.
;-
;-------------------------------------------------------------------------------
on_error,2
npar=n_params()
if npar eq 0 then begin
message,'Calling Sequence: ARC,name [,gnumber]',/cont
message,'Optional Keywords: ' + $
' textout=textout,calspec=calspec,ident=ident,noplot=noplot'
end
;
; Check for optional keywords.
;
if (not(keyword_set( calspec ))) then calspec = 0 else calspec= 1
if (not(keyword_set( textout ))) then textout = 1
if (not(keyword_set( ident ))) then ident = '_r' else $
if datatype(ident) ne 'STR' then $
message,' Invalid IDENT keyword specified: '+ident
;
; Set up plotting data structure
;
plot_scale= { plot_info, xmin_max: [0.0,0.0], $
ymin_max: [0.0,0.0], $
colors : [0,0,0], $
position: [.2,.2,.9,.9] }
;
; Decipher input naming convention (rootname or database entry) style.
;
name_type = datatype( name )
if name_type eq 'STR' then begin
; OSS file mode...assume ROOTNAME notation (i.e. z0..._r)
cname = strtrim(string(name),2)+ident
end else begin
; Archival data mode...assume DATABASE ENTRY notation (i.e. hnnnn_r)
cname = 'h' + strtrim(string(name),2) + ident
end
;
; Check whether or not the grating number has been specified
; ...grating number is used to set which line library is used.
;
if (npar lt 2) then begin
if name_type eq 'STR' then begin
;
; For user supplied input data, prompt user for Grating Number.
;
print,' '
print,' GHRS Grating Numbers '
print,' 1 = G140M'
print,' 2 = G160M'
print,' 3 = G200M'
print,' 4 = G270M'
print,' 5 = G140L'
print,' 6 = Echelle A'
print,' 7 = Echelle B'
read,' Enter Grating Number(1-7) to be used: ',gnumber
end else begin
;
; For archived data, use the GHRSLOG database.
;
dbopen,'ghrslog'
gnumber= fix(dbval(name,'gnumber'))
dbclose
end
endif
library= getlog('ZAUX') + 'nbslib'
nvar=3
if (gnumber eq 5) then begin
library= library + '_g5'
nvar=2
end
;
;
; Interactive mode. Determine if user has already run CALHRS on the
; WAVECAL observation.
;
if (not (calspec)) then begin
calspec=' '
print,' '
print,' WAVECAL obs.= ', cname
read, $
' Have you already run CALHRS on this WAVECAL obs.? Y/N [Y] ',calspec
yesno,calspec
end
if (not(calspec)) then begin
if !dump gt 0 then $
print,' Running initial CALHRS on WAVECAL observaton'
calhrs,name,'bck=00,errors=0'
endif
;
; Start processing log for WAVECAL
;
if !dump gt 0 then begin
textopen,'arc',textout=textout
if (not(!noprint)) then begin
printf,!textunit,' '
printf,!textunit,' ARC '+ !stime
printf,!textunit,'-----------------------------------------------------
printf,!textunit,' Observation: '+ strtrim(name,2)
printf,!textunit,' Grating Number: ',string(gnumber,'(i2)'), $
' Line Library: ',library
printf,!textunit,' '
endif
endif ; !dump gt 0
;
; Locate and Read calibrated WAVECAL data (output of CALHRS).
;
if name_type eq 'STR' then begin
;
; OSS file mode...assume rootname notation (i.e. z0..._r)
;
status= findfile( cname+'.plh', count=count)
if count le 0 then $
message,' Unable to locate calibrated WAVECAL obs: '+cname
endif else begin
;
; Archival data mode...assume rootname notation (i.e. hnnnn_r)
;
status= findfile( cname+'.plh', count=count)
if count le 0 then $
message,' Unable to locate calibrated WAVECAL obs: '+cname
endelse
hrs_open,cname,io
hrs_read,io,'gro',ih,cflux ; gross
hrs_read,io,'iac',ih,cwave ; lambda's with aperture offset to SSA.
hrs_close,io
csf = max(cflux)
cflux= cflux/csf ; Normalize calibrated flux from WAVECAL.
cflux= cflux(*,0) ; 1st readout only
cwave= cwave(*,0) ; 1st readout only
;
; Read the NBS line list table.
;
if !dump gt 1 then printf,!textunit, $
'Reading WAVELENGTH, INTENSITIES from : '+library
table_ext,library,'wavelength,intensity',lw,li
;
; Lookup the relevant lines from the NBS Line list table which pertain
; to this WAVECAL observation.
;
if !dump gt 1 then printf,!textunit, $
'Extracting relevant WAVELENGTH, INTENSITIES for this WAVECAL'
list = where((lw gt min(cwave)) and (lw lt max(cwave)), found)
if (found le 0) then begin
message, + $
' Unable to locate WAVECAL region in NBS Line List: '+library,/cont
message, + $
' Range of NBS(min,max): '+strtrim(min(lw),2)+ ' , ' + $
strtrim(max(lw),2),/cont
message, + $
' Range of obs.(min,max): '+strtrim(min(cwave),2)+' , ' + $
strtrim(max(cwave),2)
endif
;
; Lookup relevant wavelengths and intensities
;
npts= n_elements(list)
lw = lw(list)
li = li(list)
if !dump gt 0 then printf,!textunit, $
'Number of lines selected from NBS Line List: ',npts
;
; Generate a spectrum containing the lines
;
if !dump gt 1 then printf,!textunit, $
'Creating synthetic spectrum from NBS Line List'
nstep = n_elements(cwave)/500.
wlines= cwave
flines= fltarr(n_elements(wlines))
;
; Select each line and store in synthetic wavelength and flux vectors.
;
for i= 0, npts-1 do begin
wl = lw(i)
tabinv, wlines, wl, ieff
k = fix(ieff + 0.5)
flines(k) = li(i)
endfor
lsf = max(flines)
flines= flines/max(flines) ; normalize the synthetic flux vector
;
; Check for a shift between the two spectra (Cross Correlation)
;
!err=0
search_width = 40
if !dump gt 1 then begin
printf,!textunit, ' '
printf,!textunit, $
'Cross-correlating WAVECAL and synthetic spectrum using'
printf,!textunit, $
' (HRS_OFFSET) with search width of '+ strtrim(search_width,2)
printf,!textunit, ' '
endif
hrs_offset,cflux,flines,offset,0,search_width ; cross-correlate
;
; Display results of cross-correlation in pixels and diodes.
;
doffset = offset/nstep
if !dump gt 1 then printf,!textunit, $
' Cross-correlation offset: ',offset,' PIXELS or ',doffset,' DIODES'
if (not(keyword_set(noplot))) then begin
;
; Display the WAVECAL observation and the synthetic spectrum.
;
print,' '
print,' ARC will display the WAVECAL observation and the synthetic spectrum.'
print,' Rescale the plot (X/Y min. and max.) to examine the fits'
print,' '
wait,1
arc_cplot,cwave,cflux,wlines,flines,/init_plot,plot_info=plot_scale
end
int_lim=0.
lim=0
if (not(keyword_set(noplot))) then begin
; Interactive mode: Find the user specified limits for the useful lines.
; This step allows the user to remove noisy data points from the
; process to reduce confusion.
lim=' '
read,'Do you want to specify LOWER LIMITS for calibration lines Y/N [Y] ',lim
yesno,lim
;
; Let user specify LOWER LIMIT for calibration lines using the cursor.
;
if (lim) then begin
retry:
arc_cplot,cwave,cflux,wlines,flines,init_plot=0,plot_info=plot_scale
print,'Set the cursor at the lower limit for calibration lines'
cursor,xd,yd
int_lim=yd
print,' '
print,'Minimum COUNTS: ',int_lim*csf
print,'Minimum LINE INTENSITY: ',int_lim*lsf
print,' '
ans=' '
read,'Is this OK ? [Y]/N ',ans
yesno,ans
if (not(ans)) then goto,retry
endif
end
;
param = 'linlib='+library
param = param + ',offset=' + strtrim(string(doffset),2)
if lim then begin
param = param + ',MINC=' + strtrim(string(int_lim*csf),2)
param = param + ',MINLIB=' + strtrim(string(int_lim*lsf),2)
end
if gnumber gt 5 then param = param + ',NVAR=2' ; Echelle mode (2 variables)
if !dump gt 0 then begin
printf,!textunit,'WAVECAL parameter history:
printf,!textunit,' '+ param
end
textclose,textout=textout
print,'RUNNING FINAL WAVECAL ANALYSIS'
wavecal,name,param
print,'COMPLETED WAVECAL ANALYSIS'
;
; Display Results of wavecal from dispersion constant table...DC_<id>.
;
for i=0,2 do print,' '
print,'RESULTS OF THE ANALYSIS'
textopen,'arc_results',textout=textout
tabname='DC_'+strtrim(string(name),2)
table_ext,tabname,'gnumber , nlines, nfound, nkept, nfit, rms',$
gnumber,nlines,nfound,nkept,nfit,rms
printf,!textunit,'analysis for grating ',strtrim(gnumber,2)
printf,!textunit,'number of lines in the library = ', $
strtrim(string(nlines),2)
printf,!textunit,'lines found = ', strtrim(string(nfound),2)
printf,!textunit,'lines kept =', strtrim(string(nkept),2)
printf,!textunit,'lines used in the fit = ',strtrim(string(nfit),2)
printf,!textunit,'RMS for the fit = ',strtrim(string(rms),2)
table_ext,tabname,'ave_offset',avoff
printf,!textunit,'Average offset (diodes)= ',strtrim(string(avoff),2)
printf,!textunit,' '
lname='LINE_'+ strtrim(string(name),2)
print,'More information on the lines is in ',lname
print,'More information on the dispersion coef is in ',tabname
print,' '
;
; Display dispersion scatter plots
;
if (not(keyword_set(noplot))) then begin
pltdis=' '
read,'Plot Dispersions [Y]/N ',pltdis
yesno,pltdis
if pltdis then arc_plot_dis,lname,cwave
end
;
textclose,textout=textout
return
end