;--------------------------------------------------------------------------- ; nis_c_back = NIS spectra background utility ; R. Sterner, 2000 Apr 27 ; ; Examples: ; To process a NIS file: nis_c_back,file=file ; Use only once for each NIS file. ; To set a reference spectrum index: nis_c_back,refindex=inr ; Should only need once for each NIS file. ; To set the spectrum to process: nis_c_back,spcindex=ins ; To get the fit to the background: nis_c_back,getfit=fit ; To get the spectrum that was processed: nis_c_back,getspc=spc ; To get the reference spectrum: nis_c_back,getref=ref ; To get the sample index array: nis_c_back,getsmp=smp ; To get the selected background points: nis_c_back,getind=ind ; To get 2-D NIS spectra array: nis_c_back,getdat=dat ; Can set the indices and do all the gets in one call: ; nis_c_back,getfit=fit,getspc=spc,getref=ref,spcindex=93 ; To list number of points and indices: nis_c_back,/list ; To set the smoothing window size: nis_c_back,win=frac ; To plot ref spect and background: nis_c_back,/rplot ; To plot processed spect and background: nis_c_back,/splot ; To check fit: nis_c_back,spcindex=84,/splot ;--------------------------------------------------------------------------- ; A modified version called nis_c_back is called by nis_proc ; ; ABSTRACT ; Goes into sequence and creates a dark background array that is the ; same length as the spectrum arraw (same number of records). program ; averages the space-darks between actual asteroid data, and ; interpolates dark across the data for a single channel. The dark ; array may then be subtracted from the data array for that channel ; for a dark correction. Procedureneeds to be run iteratively for ; each channel, as one run only does a single channel. Ref channel ; is channel 22 ; ; INPUTS ; See Ray's commentary ; ; USES ; A wide variety of Ray's homecooked procedures are used (all replicated ; in the niscal library) ; rayhist - histogram program ; getrun - N'th run of consecutive integers from a given array. ; fixgaps - Linearly interpolate across data gaps in an array ; nruns - Gives # of runs of consecutive integers in given array. ; smooth2 - Do multiple smoothing. Gives near Gaussian smoothing. ; blend - Sets up weighting array for blending arrays together. ; linfill - Fill gaps in array data using linear interpolation. ; ; HISTORY ; Created 4/27/00 by R. Sterner ; Modified 05/08/00 (RS & NRI) Remove Shutter-closed dark data ; Modify for fitting into niscal (v4.1) ; 07/20/01 (NRI) v5.0 Compliant (mainly CNISF 5.0) ;------------------------------------------------------------------------------ ;--------------------------------------------------------------------------- ; nis_c_back_pts.pro = Find NIS background points from reference ;--------------------------------------------------------------------------- pro nis_c_back_pts, r, smp, inb, err npts = n_elements(r) smp = findgen(npts) ; Sample number. ;-------------------------------------------------------------- ; ; Find background points ; ; Work with reference spectrum to find background ; This assumes background for reference spectrum is flat. ; Ideally the min should give the ref background value, ; but using the histogram method allows a bit more noise. ; The flat background will give a spike in the histogram. ; Then the found background value is located in the ref spect. ; These points will occur in segments, shrink the segments ; a bit to avoid getting too close the signal. ;-------------------------------------------------------------- rmx = max(r,min=rmn) ; Find # hist bins needed. mxbins = round(rmx-rmn+5) h = rayhist(r,x,1,maxbins=mxbins) ; Histogram. w = where(h eq max(h)) ; Find histogram peak. ip = w(0) ; Index of peak. xp = x(ip) ; Value of peak. hw = 2. ; Halfwidth around peak. w0 = where(abs(r-xp) le hw) ; Find background. drop = 2 ; Shrink segments this many samples. sm = 2*(drop+0.5) ; Smoothing window to do that. one = r*0. ; Array of 1s. one(w0) = 1. ; Mark backgr pts. one_b = smooth(one,sm) ; Smooth back away end pts. inb = where(one_b eq 1) ; New estimate of backgr pts. ;err=stdev(r(inb)) bkar=r(inb) erar=bkar*0. sz=size(inb) sur=4 ;surrounding window for stdev +/- sur places for erl = 0,sz(1)-1 do begin bot=(erl-sur)*(erl ge sur and (erl lt (sz(1)-sur))) + $ (sz(1)-(2*sur+1))*(erl ge sz(1)-sur) if bot lt 0 then bot = 0 ; Added in 5.0 top=(bot+(2*sur))*(erl lt (sz(1)-sur)) + $ (sz(1)-1) * (erl ge (sz(1)-sur)) if bot lt 0 then bot = 0 if sz(0) eq 1 and top ge sz(1) then top = sz(1)-1 erar(erl)=stdev(bkar(bot:top)) endfor sz2=size(r) err=interpol(erar,inb,findgen(sz2(2))) end ;--------------------------------------------------------------------------- ; nis_c_back_fit.pro = Simple smoothing fit ;--------------------------------------------------------------------------- pro nis_c_back_fit, bx, by, xfit, yf, win=win, help=hlp if (n_params(0) lt 4) or keyword_set(hlp) then begin print,' Do a simple curve smoothing with end effects fixed.' print,' nis_c_back_fit, bx, by, xfit, yfit' print,' bx, by = input samples along curve. in' print,' xfit = array of sample numbers (0 to n). in' print,' yfit = returned smoothed curve. out' print,' Keywords:' print,' WIN=win Smoothing size as fraction of total' print,' curve length (def=0.02).' return endif if n_elements(win) eq 0 then win=0.02 num = n_elements(xfit) ; Size of fit array. swin = fix((win*num)/2)*2+1 ; Force smooth window odd. ed = (swin+1)/2 ; Edge effects size. ;------- Set up connected points as start --------- yf = 0.*xfit ; Array of correct size. yf(bx) = by ; Insert points. yf = fixgaps(yf) ; Connect points. ye=yf ;------- Also set up connected averaged points ------ nr = nruns(bx) mx = fltarr(nr) my = fltarr(nr) yf2 = 0.*xfit ; Array of correct size. for i=0,nr-1 do begin in = getrun(bx,i, loc=inlo) inhi = inlo+n_elements(in)-1 mx(i) = mean(in) my(i) = mean(by(inlo:inhi)) endfor yf2(mx) = my yf2 = fixgaps(yf2) ; Connect points. ;------- Do a simple multi-boxcar smoothing ------- yf = smooth2(yf,swin) ;------- Fix edge effects ------------------------- wt = blend(ed,ed/2,ed-1) ;---- Start --------- yf(0) = yf(0:ed-1)*wt+yf2(0:ed-1)*(1-wt) ;---- End ----------- wt = reverse(wt) ins = num-ed yf(ins) = yf(ins:*)*wt+yf2(ins:*)*(1-wt) end ;--------------------------------------------------------------------------- ; nis_c_back.pro = NIS background utility ;--------------------------------------------------------------------------- pro nis_c_back, seq, refindex=rindx, rplot=rplot, splot=splot, $ spcindex=sindx, list=list, error=err, win=win, $ getfit=gfit, geterr=gerr, getref=gref, getspc=gspc, getsmp=gsmp, $ getind=gind, getdat=gdat, help=hlp, dm=dm ;------------------------------------------------------------------ ; Common to save state. ;------------------------------------------------------------------ common nis_c_back_com, $ nfile, $ ; NIS file with ref spectra. data, $ ; Ref spectra data. smp, $ ; Sample index array. ref, $ ; Reference spectrum. spc, $ ; Spectrum to process. refindex, $ ; Index of ref spectrum. spcindex, $ ; Index of spectrum to process. win0, $ ; Last window size. inb, $ ; Indices of background points. bckx, bcky, $ ; Selected background points (spc). bckfit ; Fitted background (spc). if keyword_set(hlp) then begin print,' Find a smoothed background subtraction for NIS data.' print,' nis_c_back' print,' All args are keywords.' print,' Keywords:' ; print,' FILE=file Name of NIS spectra file.' print,' REFINDEX=rind Reference spectra index (def=73).' print,' SPECINDEX=sind Spectrum to process index (def=105).' print,' WIN=frac Smoothing window size as fraction.' print,' GETFIT=fit Get fitted background.' print,' GETREF=ref Get reference spectrum.' print,' GETSPC=spc Get processed spectrum.' print,' GETSMP=smp Get sample index array.' print,' GETIND=ind Get background pt indices.' print,' GETDAT=dat Get 2-d spectra array.' print,' GETERR=err Get background error (stdev).' print,' dm=dm Dark_mode' print,' /LIST List status.' print,' /PLOT plot spc and fitted background.' print,' Note: Reference spectra is read in and processed when' print,' FILE is given.' return endif ;------------------------------------------------------------------ ; Init some values ;------------------------------------------------------------------ if n_elements(refindex) eq 0 then refindex=-1 if n_elements(spcindex) eq 0 then spcindex=-1 if n_elements(win0) eq 0 then win0=-1 if n_elements(rindx) eq 0 then rindx=-1 if n_elements(sindx) eq 0 then sindx=-1 if n_elements(win) eq 0 then win=-1 if n_elements(dm) eq 0 then dm=5 ;------------------------------------------------------------------ ; Deal with NIS file ;------------------------------------------------------------------ data=seq refindex = -1 ; Clear spectra indices. spcindex = -1 ;endif ;------------------------------------------------------------------ ; Deal with reference spectra ;------------------------------------------------------------------ if n_elements(data) eq 0 then begin print,' Error in nis_c_back: no NIS file given yet.' return endif if rindx eq -1 then rindx=refindex if rindx eq -1 then rindx=153 ; No index given, assume default. if refindex ne rindx then begin ; Processed yet? ;print,' Reading ref spectra ',rindx ref = data(rindx,*) ; Grab spectrum. nis_c_back_pts, ref, smp, inb, err ; Find backgr pts. sh=data(6,*) ; Get shutter info if dm lt 7 and total(where(sh(inb) le 1)) gt -1 then $ inb=inb(where(sh(inb) le 1)) ; Eliminate shutter closed spctra refindex = rindx ; Remember which ref spect. endif ;------------------------------------------------------------------ ; Deal with spectrum to process ;------------------------------------------------------------------ if sindx eq -1 then sindx=spcindex if sindx eq -1 then sindx=185 ; No index given, assume default. if spcindex ne sindx then begin ; Processed yet? ;print,' Reading spectra ',sindx spc = data(sindx,*) ; Grab spectrum. bckx = smp(inb) bcky = spc(inb) spcindex = sindx ; Remember which spect. win0 = -1 ; Force new processing. endif ;------------------------------------------------------------------ ; Deal with smoothing window ;------------------------------------------------------------------ if win eq -1 then win=win0 if win eq -1 then win=0.02 ; No window given, assume default. if win0 ne win then begin ; Processed yet? ;print,' Smoothing with ',win nis_c_back_fit, bckx, bcky, smp, bckfit, win=win win0 = win ; Remember window. endif ;------------------------------------------------------------------ ; List ;------------------------------------------------------------------ if keyword_set(list) then begin print,' ' if n_elements(nfile) eq 0 then txt=' Not given yet' else $ txt=nfile print,' NIS File: '+txt if n_elements(smp) eq 0 then txt=' None' else $ txt=strtrim(n_elements(smp),2) print,' # samples: '+txt if n_elements(refindex) eq 0 then txt=' None' else $ txt=strtrim(refindex,2) print,' Index of Ref: '+txt if n_elements(spcindex) eq 0 then txt=' None' else $ txt=strtrim(spcindex,2) print,' Index of Spec: '+txt if n_elements(inb) eq 0 then txt=' None' else $ txt=strtrim(n_elements(inb),2) print,' # background pts: '+txt endif ;------------------------------------------------------------------ ; Plot reference spectrum ;------------------------------------------------------------------ if keyword_set(rplot) then begin plot,smp,ref,/ynoz,xtitle='Sample',ytitle='Value',title=$ 'Reference Spectrum at index '+strtrim(rindx,2) oplot,smp(inb),ref(inb),psym=3,col=tarclr(255,0,0) endif ;------------------------------------------------------------------ ; Plot processed spectrum ;------------------------------------------------------------------ if keyword_set(splot) then begin plot,smp,spc,/ynoz,xtitle='Sample',ytitle='Value',title=$ 'Processed Spectrum at index '+strtrim(rindx-130,2)+ $ ' Channel '+strtrim(sindx-130,2) oplot,bckx,bcky,psym=2,col=tarclr(255,0,0) oplot,smp,bckfit,col=tarclr(255,255,0) endif ;------------------------------------------------------------------ ; Returned items ;------------------------------------------------------------------ if arg_present(gfit) then gfit=bckfit if arg_present(gref) then gref=ref if arg_present(gspc) then gspc=spc if arg_present(gsmp) then gsmp=smp if arg_present(gind) then gind=inb if arg_present(gdat) then gdat=data if arg_present(gerr) then gerr=err end