Viewing contents of file '../idllib/contrib/buie/ccdgain.pro'
;+
; NAME:
; ccdgain
; PURPOSE: (one line)
; Extract and plot CCD gain transfer curve from flat field image data.
; DESCRIPTION:
; CATEGORY:
; CCD data processing
; CALLING SEQUENCE:
; ccdgain,dataset,flatname,biasname,frames
; INPUTS:
; root - Root of the image file name (no path, with period).
; flatname - Filename of master flat field.
; biasname - Filename of bias frame.
; frames - List of frame id's to process.
; OPTIONAL INPUT PARAMETERS:
; KEYWORD PARAMETERS:
; OUTPUTS:
; COMMON BLOCKS:
; SIDE EFFECTS:
; RESTRICTIONS:
; PROCEDURE:
; MODIFICATION HISTORY:
; 93/04/03 - Initial conversion to procedure, Marc W. Buie, Lowell Obs.
; 96/01/06 - MWB, added support for hardcopy under DOS/Windows
; 97/03/15 - MWB, added use of MARKDATA
;-
pro ccdgain,root,start,nframes,calibfile, $
EXCLUDE=exclude,FILKEY=filkey,PATH=path,CALIBPATH=calibpath, $
SECTION=section,EXPKEY=expkey
if badpar(root,7,0,caller='CCDGAIN: (root) ') then return
if badpar(start, [2,3],0,caller='CCDGAIN: (start) ' ) then return
if badpar(nframes,[2,3],0,caller='CCDGAIN: (nframes) ') then return
if badpar(calibfile,[0,7],0,caller='CCDGAIN: (calibfile) ',default='files.cal') then return
if badpar(filkey,[0,7],0,caller='CCDGAIN: (filkey) ',default='FILTER') then return
if badpar(expkey,[0,7],0,caller='CCDGAIN: (expkey) ',default='EXPTIME') then return
if badpar(path,[0,7],0,caller='CCDGAIN: (path) ',default='./') then return
path=addslash(path)
defcalib=path+'calib/'
if badpar(calibpath,[0,7],0,caller='CCDGAIN: (calibpath) ', $
default=defcalib) then return
if badpar(exclude,[0,2,3],[0,1],caller='CCDGAIN: (exclude) ',default=-1) then return
if badpar(section,[0,2,3],1,caller='CCDGAIN: (section) ', $
default=[100,300,100,300],npts=nsect) then return
if nsect ne 4 then $
message,'SECTION must contain four elements'
frames=start+indgen(nframes)
bad=intarr(nframes)
print,'Load frames ',start,' to ',start+nframes-1
for i=0,nframes-1 do begin
z=where(frames[i] eq exclude,count)
if count ne 0 then bad[i]=1
endfor
sel=where(bad eq 0,count)
if count eq 0 then $
message,'Error ** you have excluded all frames, nothing to do.'
; Load calibration information
ldcalib,calibfile,calib,valid,calibpath=calibpath
if not valid then return
sectstr='['+string(section[0])+':'+ $
string(section[1])+','+ $
string(section[2])+':'+ $
string(section[3])+']'
sectstr=strcompress(sectstr,/remove_all)
sectstr='Image section '+sectstr
signal=fltarr(nframes)
variance=fltarr(nframes)
print,' '
print,sectstr
print,' '
print,' filename mean avgdev stddev variance skew kurt nfinal'
for i=0,nframes-1 do begin
if frames[i] ne -1 then begin
filename = root+string(frames[i],format='(i3.3)')
if not exists(path+filename) then begin
print,'Image file ',path+filename,' could not be found. Aborting.'
return
endif
image = readfits(path+filename,hdr,/silent)
filter = sxpar(hdr,filkey,filename)
exptime = sxpar(hdr,expkey,filename)
flatcode = where(filter eq calib.filter)
flatcode=flatcode[0]
if flatcode[0] eq -1 then $
message,'Filter '+filter+' does not have a corresponding flat field'
; Overscan correction.
if calib.xl ge 0 and calib.xr ge 0 and calib.xl lt calib.xr then $
os = mean( image[calib.xl:calib.xr,*] ) $
else $
os = 0.0
if n_elements(calib.dark) eq 1 then $
image = (image[calib.x1:calib.x2,calib.y1:calib.y2] - calib.bias - os) / $
calib.flat[*,*,flatcode] $
else $
image = (image[calib.x1:calib.x2,calib.y1:calib.y2] - calib.bias - $
calib.dark*exptime - os ) / $
calib.flat[*,*,flatcode]
robomean,image[section[0]:section[1],section[2]:section[3]],3.0,0.5, $
avg,avgdev,stddev,var,skew,kurt,nfinal
signal[i] = avg
variance[i] = var
print,filename,exptime,avg,stddev,var, $
format='(a,1x,f6.1,1x,f8.2,1x,f7.2,1x,f9.2)'
endif
endfor
done=0
print=0
while not done do begin
coeff=goodpoly(signal[sel],variance[sel],1,3,yfit,newx,newy)
xlin=minmax(signal[sel])
lin=coeff[0]+coeff[1]*xlin
if coeff[0] ge 0.0 then $
readnoise=sqrt(coeff[0])/coeff[1] $
else $
readnoise=coeff[0]
gain = 1/coeff[1]
if print then begin
devsave=!d.name
pfont=!p.font
set_plot,'ps'
device,/landscape,/Helvetica
!p.font=-1
endif
xr=minmax([0,xlin])
yr=minmax([0,lin,variance[sel]])
plot,xlin,lin,xtit='Signal (DN)',ytit='Variance (DN)!u2', $
title='Gain Transfer Curve for '+root,xtickformat='(i5)',yr=yr,xr=xr
setusym,-1
oplot,signal,variance,psym=8
setusym,1
oplot,newx,newy,psym=8
xyouts,0.15,0.9,'Read-noise '+string(readnoise,format='(f6.1)')+' e-', $
/normal,charsize=1.3
xyouts,0.15,0.85,'Gain '+string(gain,format='(f5.2)')+' e-/ADU', $
/normal,charsize=1.3
xyouts,0.15,0.80,sectstr,/normal,charsize=1.3
for i=0,n_elements(frames)-1 do $
xyouts,signal[i],variance[i],string(frames[i]),align=0.5
if not print then begin
oldbad = bad
markdata,signal,variance,bad, $
xtitle='Variance (DN)',ytitle='Signal (DN)', $
ptitle='CCD Gain transfer curve'
sel = where(bad eq 0,count)
z2bad = where(bad ne oldbad and bad eq 1, countgonebad)
z2good = where(bad ne oldbad and bad eq 0, countgonegood)
if countgonebad ne 0 then print,'Gone bad : ',frames[z2bad]
if countgonegood ne 0 then print,'Gone good: ',frames[z2good]
if countgonebad eq 0 and countgonegood eq 0 then begin
hc='y'
read,prompt='hardcopy? (y/n) ',hc
if hc eq 'y' then $
print=1 $
else $
done=1
endif
endif else begin
done=1
device,/close
if !version.os_family eq 'unix' then begin
cmd='lpr idl.ps'
endif else if !version.os_family eq 'Windows' then begin
cmd = 'copy idl.ps lpt1:'
endif
spawn,cmd
set_plot,devsave
!p.font=pfont
endelse
endwhile
end