;+
;GRAFIT.PRO
; IDL program to
; 1. read in a list of line locations
; 2. figure out best initial parameters for each line
; using spectrum (YCT(XWV)) already read in.
; 3. fit each line to data within the spectrum
; 4. write out the results
;
;USAGE
; read in spectrum and store in variables YCT and XWV
; LWIDTH0=0.05 & LWIDTH_MAX=0.1 & PSROOT='capella'
; .run grafit
;
;expected to be of use in Chandra HETG/LETG analyses
;
;SUBROUTINES
; wc
; setcolor
; nary2dec
; num2str
; fit_levmar
; x3model
; mk_model
; mk_gauss
; mk_lorentz
; adjustie
; levmarq
; lmcoeff
;
;vinay kashyap (Aug99)
;-
; plotting stuff
window,0 & window,1
loadct,3
setcolor,'#ff0000',2
setcolor,'#ffff00',3
; initialize those fit parameters as necessary
if not keyword_set(LWIDTH0) then LWIDTH0=0.01 ;[Ang]
if not keyword_set(LWIDTH_max) then LWIDTH_max=0.05 ;[Ang]
if not keyword_set(functyp) then functyp='lorentz'
if not keyword_set(psroot) then psroot='PLOT'
; is spectrum available?
nxwv=n_elements(xwv) & nyct=n_elements(yct) & ok='ok'
if nxwv eq 0 then ok='XWV not defined' else $
if nyct eq 0 then ok='YCT(XWV) not defined' else $
if nxwv ne nyct then ok='YCT(XWV) and XWV incompatible' else $
if nxwv lt 2 then ok='too few elements to bother!'
if ok ne 'ok' then message,ok else plot,xwv,yct,psym=10,/xs,/ys,$
xtitle='Wavelength [Ang]',ytitle='Counts [ct]'
wset,0
; initialize input line list filename..
if not keyword_set(linefil) then linefil='linelist.txt'
; expect the input file to have the following format:
; wavelength [Ang] (f8.3)
; Element Ion (a10)
; e- configuration (a20)
; level designation (a20)
; counts in line (f8.0)
; of which only the first 2 columns are of relevance
;
;example print statement to make the linelist
;print,format='(f8.3,2x,a10,2x,a50)',wave,z_ion,stuff
; .. and read in the line list
nlin=wc(linefil) & lwvl=fltarr(nlin) & lelem=strarr(nlin)
openr,ulin,linefil,/get_lun
klin=0L
for i=0L,nlin-1L do begin
cc='' & readf,ulin,cc & c1=strmid(cc,0,1)
if c1 ne '#' and c1 ne ';' then begin
aa=0. & bb='' & reads,cc,aa,bb,format='(f8.3,2x,a10)'
lwvl(klin)=aa & lelem(klin)=bb
klin=klin+1L
endif
endfor
if klin lt nlin then begin
lwvl=lwvl(0:klin-1L) & lelem=lelem(0:klin-1L) & nlin=klin
endif
close,ulin & free_lun,ulin
; initialize parameter arrays
lpos=lwvl & lwdt=0.*lwvl+LWIDTH0(0) & lhgt=0.*lwvl+1. ;to hold best-fit values
lpos_e=0.*lpos & lwdt_e=0.*lwdt & lhgt_e=0.*lhgt ;errors
lquality=strarr(n_elements(lpos)) ;quality flag
oldlpos=lpos
xline=intarr(n_elements(lpos)) ;dealt with line(1) or no(0)?
; and array to hold chisq
lchi=0.*lwvl & lrchi=lchi
; run specific initializations
xmin=min(XWV,max=xmax) & nycte=n_elements(YCT_E)
YCT_E=sqrt(abs(YCT)+0.75)+1. ;assume poisson errors
dXWV=XWV(1:*)-XWV & delX=median(dXWV)
if n_elements(dumb) eq 0 then dumb=1 ;see FIT_LEVMAR
; remove lines
bkgval=0.*YCT & bkgerr=bkgval
;bg=linerem(XWV,YCT,sig=YCT_E,nsigma=1,cell=[1,1,1,0,0,0,1,1,1],$
;bg=linerem(XWV,YCT,sig=YCT_E,nsigma=1,cell=[-1.,-1.,0,0,0,0,-1.,-1.],$
; bkgval=bkgval,bkgerr=bkgerr,quiet=0)
; step through and fit
ilin = 0L & dlin=1L
while ilin lt nlin do begin ;{step through linelist
loopbegin: i = ilin
lxrng=lpos(i)+1*LWIDTH0*[-1.,1.] ;range to check for nearby lines
sxrng=lpos(i)+3*LWIDTH0*[-1.,1.] ;range in which to fit the line
pxrng=lpos(i)+7*LWIDTH0*[-1.,1.] ;range in which to display plots
sxrng=((sxrng>xmin) ('+$
strtrim(lpos(ol(j))-2*LWIDTH0,2)+'))'
ties(k+1)=aw+'=(('+aw+' < ('+strtrim(LWIDTH_max,2)+')) > ('+$
strtrim(delX,2)+'))'
ties(k+2)=ah+'=('+ah+' > 0)'
endfor
; fit
fit_levmar,xx,yy,aa,erra=erra,chisq=chisq,funcs='x3model',$
sig=yy_e,/norm,type=functyp,ties=ties,dumb=dumb,$
jumpup=jumpup,jumpdn=jumpdn
lchi(i)=chisq & if mow gt 3 then lrchi(i)=chisq/(mow-3.)
; plot
call_procedure,'x3model',xx,aa,ymod,/norm,type=functyp
plot,XWV,YCT,/xs,psym=10,xr=pxrng,$
title=lelem(i),xtitle='[Ang]',ytitle='[counts]'
for j=0L,mow-1L do oplot,XWV(ow(j))*[1,1],YCT(ow(j))+YCT_E(ow(j))*[-1,1]
xyouts,pxrng(0)+delX,max(YCT(ow))*0.95,$
'!4v!3!u2!dr!n='+strtrim(lrchi(i),2)+$
'; !3C='+strtrim(total(ymod),2)
oplot,xx,ymod+bkgval(ow),col=3,psym=-1
oplot,xx,bkgval(ow),col=2
nextline:
; save the fit parameters
if mol eq 1 then begin ;(one line to worry about
lpos(i)=aa(0) & lpos_e(i)=erra(0)
lwdt(i)=aa(1) & lwdt_e(i)=erra(1)
lhgt(i)=aa(2) & lhgt_e(i)=erra(2)
endif else begin ;)(more than one line
tmp=min(abs(lpos(ol)-lpos(i)),idx)
lpos(i)=aa(3*idx+0) & lpos_e(i)=erra(3*idx+0)
lwdt(i)=aa(3*idx+1) & lwdt_e(i)=erra(3*idx+1)
lhgt(i)=aa(3*idx+2) & lhgt_e(i)=erra(3*idx+2)
;lpos(i)=lpos(i)
;lwdt(i)=LWIDTH0
;lhgt(i)=lhgt(i)*renorm
message,'cannot handle this case yet!',/info
endelse ;MOL)
LWIDTH0=lwdt(i)
; add parameters and errors to plot
xyouts,pxrng(0)+delX,max(YCT(ow))*0.8,$
'POS: '+strtrim(lpos(i),2)+' +- '+strtrim(lpos_e(i),2)
xyouts,pxrng(0)+delX,max(YCT(ow))*0.7,$
'WDT: '+strtrim(lwdt(i),2)+' +- '+strtrim(lwdt_e(i),2)
xyouts,pxrng(0)+delX,max(YCT(ow))*0.6,$
'HGT: '+strtrim(lhgt(i),2)+' +- '+strtrim(lhgt_e(i),2)
endif ;MOW)
ilin0=ilin
if not keyword_set(trust) then begin
print,'hit any key to continue, t/T for thinner/thicker line,'
c='p for previous, q/x to stop, r to repeat, e for errors'
print,c & c=get_kbrd(1)
endif else c='ok'
if c eq 'p' or c eq 'P' then begin
print,'current line is: '+strtrim(ilin,2),lwvl(ilin),' ',lelem(ilin)
ilin = ilin - dlin
xline(ilin)=0
print,'new line is: '+strtrim(ilin,2),lwvl(ilin),' ',lelem(ilin)
goto, loopbegin
endif
if c eq 'q' or c eq 'Q' then c='x'
if c eq 'x' or c eq 'X' then stop,'HALTING. type .CON to continue'
if c eq 't' or c eq 'T' then begin
if c eq 't' then LWIDTH0=LWIDTH0/2. else LWIDTH0=LWIDTH0*1.1
c='x'
endif
if c eq 'e' then begin ;(compute errors
wset,1
erors,xx,yy,aa,erraa,ysig=yy_e,dchi=2.7,verbose=10,$
funcs='x3model',/norm,type=functyp,ties=ties,/dumb,$
jumpup=jumpup,jumpdn=jumpdn
print,'OLD:',erra & erra=erraa & print,'NEW:',erra
wset,0
if mol eq 1 then begin ;(one line to worry about
lpos_e(i)=erra(0)
lwdt_e(i)=erra(1)
lhgt_e(i)=erra(2)
endif else begin ;)(more than one line
tmp=min(abs(lpos(ol)-lpos(i)),idx)
lpos_e(i)=erra(3*idx+0)
lwdt_e(i)=erra(3*idx+1)
lhgt_e(i)=erra(3*idx+2)
message,'cannot handle this case yet!',/info
endelse ;MOL)
endif ;C='E')
if c ne 'r' and c ne 'R' and c ne 'x' and c ne 'X' then ilin=ilin+dlin
if ilin gt ilin0 and mow gt 0 then begin ;(all done, final convulsions
cqual='' & read,prompt='type in a quality comment> ',cqual
lquality(i)=strtrim(cqual,2)
;
xline(ool)=1
;
if keyword_set(psroot) then begin
psfil=strtrim(psroot,2)+'_'+$
num2str(oldlpos(i),len=strlen(strtrim(oldlpos(i),2)))+'.ps'
gifil=strtrim(psroot,2)+'_'+$
num2str(oldlpos(i),len=strlen(strtrim(oldlpos(i),2)))+'.gif'
message,'Writing figure to postscript file: '+psfil,/info
message,'Writing figure to GIF file: '+gifil,/info
tvlct,r,g,b,/get
set_plot,'ps' & device,file=psfil
endif
plot,XWV,YCT,/xs,psym=10,xr=pxrng,thick=2,$
title=lelem(i),xtitle='Wavelength [Ang]',ytitle='counts/pix'
for j=0L,mow-1L do oplot,XWV(ow(j))*[1,1],YCT(ow(j))+YCT_E(ow(j))*[-1,1]
xyouts,pxrng(0)+delX,max(YCT(ow))*0.95,$
'!4v!3!u2!dr!n='+strtrim(lrchi(i),2)+$
'; !3C='+strtrim(total(ymod),2)
oplot,xx,ymod+bkgval(ow),col=3,psym=10 & oplot,xx,bkgval(ow),col=3
; add parameters and errors to plot
xyouts,pxrng(0)+delX,max(YCT(ow))*0.8,$
'POS: '+strtrim(lpos(i),2)+' +- '+strtrim(lpos_e(i),2)
xyouts,pxrng(0)+delX,max(YCT(ow))*0.7,$
'WDT: '+strtrim(lwdt(i),2)+' +- '+strtrim(lwdt_e(i),2)
xyouts,pxrng(0)+delX,max(YCT(ow))*0.6,$
'HGT: '+strtrim(lhgt(i),2)+' +- '+strtrim(lhgt_e(i),2)
if keyword_set(psroot) then begin
device,/close & set_plot,'x'
write_gif,gifil,tvrd(),r,g,b
endif
endif ;ILIN>ILIN0 && MOW>0)
endwhile ;ILIN < NLIN}
; summarize
print,'OUTPUTS in:-' & help,lpos,lwdt,lhgt,lpos_e,lwdt_e,lhgt_e,lchi,lrchi,lelem
print,string("7b)
print,'Fit to spectrum:-' & help,YCT,XWV,YCT_E,functyp,linefil
print,string("7b)
print,'Other variables of interest:-' & help,LWIDTH0,LWIDTH_max,psroot,outfil
if keyword_set(outfil) then begin
message,'saving to file: '+outfil,/info
save,file=outfil
endif
end