;+ ; NAME: clean_i ; ; CALLING SEQUENCE: clean_i, i_arr, header, center=ocenter, doplot=doplot, $ ; notelluric=notelluric, maskin=maskin, verbose=verbose ; ; DESCRIPTION: ; Cleans up I spectra by removing pre-filter tilt and ; remove telluric line if it's Fe data. ; ; INPUT PARAMETERS: ; i_arr = I-DAta are assumed to be FLOAT( x, y, w ) ; with the spatial and ; wavelength dimensions given in the header. ; header = header structure, IVM style. use ivm_rd_head('IMGR...') ; outarray = output image array, dimensioned (x, y, nlambda) ; which gives you a nicer x,y,lambda cube ; OPTIONAL INPUT KEYWORDS: ; /notelluric = do not call ivm_telluric at all. Default: ; generally, call it at the end after pre-filter removing, ; and only Fe data, unless limb data. Limb data benefits from ; ivm_telluric first. ; /doplot = plot some diagnostics for debugging purposes. ; /verbose = print progress if >3, etc. ; maskin: byte array of 1s (good to go), 0s (ignore) ; OPTIONAL OUTPUT KEYWORDS: ; center = line center array (crude). ; slope = array of linear fit results for prefilter. ; fringe = array of coefficients from fringe fitting. ; err = Set to 1, or 2 if telluric was done at beginning, if an error was encountered. ; HISTORY: ; Written October 6, 1996 Barry LaBonte ; took ft_gram.pro and hacked it to this - 7/97, KDL ; abbreviated from clean_spectrum.pro 9/97 KDL ; ; TRM 2004-Mar-15 New algorithm that does not requre the ; endpoints to be continuum and does not require the line to be ; centered in the wavelength range. Rather it figures the ; prefilter slope from points symmetric about line center. Also ; allows for a variation of the slope over the FOV by averaging ; by rows and columns rather than the whole FOV. ; ; 2004-Mar-17 TRM added sanity check for line center. If it ; fails, revert to the old algorithm. ; 2004-Oct-25 TRM Added crude line center ouput (in pixels) ; 2004-Oct-27 TRM Compute slope and curvature of prefilter ; 2004-Nov-10 TRM Check status in poly_fit ; 2010-Feb-13 KDL debugged numerous bits: ; (1) the line-center determination was very sensitive to ; how many points were used, and was unstable. Line_cent is ; just inherently very imprecise and susceptible to the existence ; of a roll-over. Now fitting simple Gaussian to 1/line profile, ; appears to be much more stable. ; (2) the variable number of points (e.g., n_elements(wr)) being used to ; fit the blue/red divided=out/normalized spectra (lrb) was leading to ; significant jumps in the quadratic fit. Confirmed less variability ; with linear fits, but linear fits changed the "cleaned" spectra ; much more significantly and not for the better that I could see. ; Modified so that the code now (a) finds line-center for whole FOV, ; then (b) figures minimum extent from line center of any pix (pmax) in FOV, ; and uses that as the limit for all fits, consistently. ; (3) Added a 10% uncertainty evently across lrb for the "measurement errors" ; for poly_fit, and voila! fits almost identical but it never finds a small ; pivot-point or singular matrix (status always returns 0 now). This could ; be played with to decrease the magnitude of the imposed error, but its ; inclusion was a dramatic (and good) result. ; (5) integrated using a cube of very very rebinned down/up spectra rather ; than full-resolution column/row averages. Leads to more realistic ; line-center positions, and avoids serious vertical/horizontal stripes ; in the corrected intensities. ; (6) cleaned up unused variables and double-negative logic ; (7) Added ability to pass in a mask of points to ignore for pukas/sky ; (8) added much-needed comments ; 2010-Feb-24 ELW improved doplot logic, verbose flag ; 2011-Jan-26 ELW got rid of superfluous 'im' array ; 2011-Jan-28 ELW indgen/findgen de-looped, lines plot w/'doplot', formatting ; 2011-June KDL (1) for line-center finding, found that working on binned data and ; binning up the line-center positions resulted in blobs. Now: use smoothed data ; for line-center finding, but do for every point. Retrieves FOV shifts, but ; less sensitive. Added external call to gauss_funct.pro and fortran gauss_funct.f90 ; gauss_sub.f90 fr speed. ; (2) Understand TRM's algorithm, thanks to GB. It does assume linear, not quadratic. ; functional form for prefilter. REMOVED quadratic option, bug fix in ; order of points input to fit (which was making it look like a quadratic was ; needed. Assumes prefilter-function is linear, and symmetric about line-center. ; Assumes solar line is symmetric about line-center. Notes in KDL log book. ; NB: still a problem that line-center was fit for the PF*I, so may be ; off. Ideally could fit for linear part and linecenter at the same time. ; but don't have time to do that now.... ; (3) to use fortran gauss_funct, do: (1) copy Makefile gauss_funct.f90 gauss_funct.pro gauss_sub.f90 (from /home/graham/Gauss/), do 'make', then be sure to .compile IDL's gaussfit.pro, *then* ; compile the local gauss_funct.pro. ; (4) Telluric line removal is now at the end. This was tested significantly, ; and judged best (theoretically & empirically) to do after prefilter slope removal. ; EXCEPT, if data are near the limb, Fe is shallow which screws up the whole ; prefilter-removal algorithm, so in THAT case, it is done at the beginning. ; (5) added the possibility of ivm_fringe (ivm_final_fringe.pro), which is the last half of ; BJL's ivm_prefilter, that removes a residual fringing pattern with ; period of half the wavelength coverage. That code was easiest to ; pull into a different code, now called from here. THIS MAY NOT BE GOOD ; IN NEAR-LIMB DATA, may need more testing. ; 7/2011: more testing => rarely doing good things as shown by bisectors, ; often doing not good things. commented the call out. ; NOTE: there is still smoothly-varying but solar-related structure that is removed. ; bisectors are pretty much straightened out. Little of the IVM data will be ; inverted in ways that variations with height of velocities/etc are useful, ; so taking this to be an OK side effect. ; 2012-Feb-24 ELW fixed flipped notelluric logic to allow normal operation; comments ; 2012-Aug-08 ELW Added err keyword ; 2012 August KDL changes: ; (1) If telluric line(s) are too deep, they negatively influence fitting ; for the prefilter slope. If they are not too deep, removing them first ; is also detrimental. Functions of expected IVM-observed line depth ; (quiet sun) are now installed, and if the telluric line is more than ; half the depth of the Fe line, then ivm_telluric is called first; if not, ; then it is called after prefilter shape removal. ;- function clean_i, i_arr, header, center=ocenter, doplot=doplot, notelluric=notelluric,$ maskin=maskin, verbose=verbose,slope=oslope, fringe=ofringe, err=err if (~keyword_set(verbose)) then verbose=0 if (~keyword_set(doplot)) then doplot=0 do_telluric = ~keyword_set(notelluric) ; ELW20120224 Added ~ err=0 ; generally, gaussfit works well but now that telluric is dealt with later by default, ; it can get confused by telluric line. So, if Fe (not Halpha, Na), do ruff guess: do_ruff=1 if (abs(header.f_wavelength - 630.) gt 1) then begin do_ruff=0 if (do_telluric) then begin if (verbose gt 0) then print,'clean_i: Na? Setting do_telluric to zero.' do_telluric=0 endif endif ; BUT.... if data are near the limb, Fe is shallow which screws up the whole ; prefilter-removal algorithm, so in THAT case, in fact do_telluric *is* usually ; a good bet, but then it can easily get confused, so make it do ruff guess: ; KDL replacing below with function to determine relative line depths ;if (do_telluric AND abs(header.ew) gt 60) then begin ; do_ruff =1 ; do_telluric=2 ; When to do telluric? 2=beginning, 1=end ;endif ;First, get expected line-depth of telluric line(s). ; Function determined empirically using code snippet from ivm_telluric.pro if (do_telluric) then begin stim = FILE2TIME(header.filename, /YOHKOH) EPHEMI, 'Haleakala', ephstr, stim am2o2coef=[0.946207,-0.0377175,0.00201715,-4.82419e-05] o2ld=am2o2coef[0]+am2o2coef[1]*ephstr.airmas + am2o2coef[2]*ephstr.airmas^2+$ am2o2coef[3]*ephstr.airmas^2 ; Next, get expected quiet-sun Fe f(mu) mu = SIN(header.b0)*SIN(header.clat) + COS(header.b0)*COS(header.clat)*COS(header.clng) ; this should be an exponential, but I couldn't get fit to behave and second-order polynomial ; fit the empirical data well enough to about mu=0.35 (70 deg from disk center) ; so this is for figuring out a conservative ratio, not for science... ; NB did convolve w/ 72-mA Gaussian (not Lorentzian, but see above comment) ; mu2fecoef=[0.664296,-0.658001,0.309840] ; feld=mu2fecoef[0]+mu2fecoef[1]*mu+mu2fecoef[2]*mu^2 ; KDL note this predicted much stronger line-depths than seen in IVM data and slightly ; stronger than in disk-center atlas. So, I looked at line-depth as f(mu) for ; *raw* IVM data and found these coefficients. mu2fecoef=[0.699735, -0.0236049,-0.00664622] feld=mu2fecoef[0]+mu2fecoef[1]*mu+mu2fecoef[2]*mu^2 print,feld,o2ld if (o2ld lt 0.88 OR feld gt 0.68) then begin do_telluric=2 ; if O2 is predicted to be deep or fe is predicted to be shallow, do telluric removal first. endif else begin do_telluric=1 ; if O2 is predicted shallow and Fe deep, then do telluric removal last endelse if (verbose ge 1) then print, 'Fe Predicted Line depth: ',feld,$ ' O2 Predicted Line-center intensity: ',o2ld,' Fe Predicted Line-center intensity: ',feld, '. Do telluric: ',do_telluric,' (1=do last, 2=do first)' endif ;if (do_telluric AND abs(header.ew) gt 60) then begin ; do_ruff =1 ; do_telluric=2 ; When to do telluric? 2=beginning, 1=end ;endif if (verbose ge 1) then print, 'CLEAN_I start: ', SYSTIME() ; Handle the input data --------------------- sz = SIZE(i_arr) ncol = sz(1) nrow = sz(2) nwavel = sz(3) ; Constants --------------------------------- ; Nanometers per Fabry-Perot step in IVM (should be spacing, not wavelength) nmfp = IVM_FP_STEP(header.f_wavelength) ; Wavelength scale, nanometers wavel = (header.fp_set(0:nwavel-1) - header.fp_set(0)) * nmfp ; Remove telluric line IF (do_telluric gt 1) THEN BEGIN ; only if doing this ahead.... if (verbose ge 1) then print, ' clean_i: calling ivm_telluric() at beginning' tim = IVM_TELLURIC(i_arr, header); NOTE, scale=1.25 may be good here, don't want to hard-wire it, though endif else begin tim = i_arr endelse ; average wavelength step size, nm; -2 removes the last (often continnum) point from consideration wstep = total((wavel - shift(wavel,1))[1:nwavel-2])/(nwavel-2.) ; set up mask if keyword_set(maskin) then $ mask=maskin $ else $ mask=bytarr(ncol,nrow)+1. mm=where(mask ne 0, mmcount) ; sometimes this format is easier if (mmcount eq 0) then begin ; Check for zero mask - ELW20110621 print, 'clean_i: WARNING! Zero mask was input. Continuing with all-1 mask.' mask=bytarr(ncol,nrow)+1 mm=where(mask ne 0, mmcount) ; sometimes this format is easier endif notmm=where(mask eq 0, notmmcount) ; sometimes this format is easier ; non-masked continuum level ;continuum = mean(((tim[*,*,0])(mm)+(tim[*,*,nwavel-1])(mm))/2.0) continuum = median((tim[*,*,nwavel-1])(mm)) ; KDL 6/2011 peak of 'continuum' wavelength. ; put wavelength scale sorted so that continuum position is in the right location. sswavel = sort(wavel) swavel = wavel[sswavel] ; sorted wavelength ; Algorithm is fairly sensitive to line-center, but do not want to include ; the finest solar structure, only fairly smoooth variation over FOV. ; So, bin I-spectra data small, then smooth & bin up again, to remove ; structure. ; First, make a temporary copy in which all puka (mask) areas are the replicated ; average intensity for that wavelength bgsmtim=fltarr(ncol,nrow,nwavel) for i=0,nwavel-1 do begin tmpim=tim[*,*,i] tmpmean=mean(tmpim(mm)) if (notmmcount ge 1) then $ ; ELW-20100224 tmpim[notmm]=tmpmean ; and smooth it waaaay out bgsmtim[*,*,i]=smoothe(tmpim,ncol/32) ; smooth the result out. endfor ; 6/2011 KDL this next part is commented out; it introduced, still, structure ; that wasn't appropriate. ;; now bin this waaayyy down ;smnx=ncol/32 > 8 ;smny=nrow/32 > 8 ;smtim=rebin(bgsmtim,smnx,smny,nwavel) ;; and then back up to original size. should remove all small-scale structure. ;; FYI, congrid ignores cubic when ndimension > 2. *grrr* ;for i=0,nwavel-1 do begin ; bgsmtim[*,*,i]=congrid(smtim[*,*,i],ncol,nrow,cubic=-0.5,/center) ;endfor ; a few arrays for output: ocenter = fltarr(ncol,nrow) oslope = fltarr(ncol,nrow) ; generally for Fe w/ or w/o O2 still there, need ruff guess. if do_ruff then begin ; if not doing telluric line removal, line-center gets confused. Better w/ rough guess ; ruff ONLY useful if there are strange shapes or secondary lines (like the O2) which destabilize. ; Else, gaussfit does pretty well. ruff=fltarr(6) ; same as terms for gaussfit. oneline=reform(rebin(bgsmtim[*,*,0:nwavel-2],1,1,nwavel-1)) ; take out continuum point ; oneline=oneline/max(oneline) ; oneline=1./oneline ; ruff[1]=5.0 + (where(oneline[5:nwavel-5] eq max(oneline[5:nwavel-5]))) ; very rough line-center ; KDL testing - former is crashing w/ big telluric line... ruff[1]=(lc_find(oneline,5,nwavel-2-5,9))[0] ;lc_find, uses whole line but looks only at ; points 5 away from edges, and uses 9 points to fit a quadratic. simple, but probably more robust... ruff[3] = 1.0 ; fitted line will be normalized so this is best guess at constant term ruff[4] = 0. ; linear term guess ruff[5]= 0. ; quadratic term guess ruff[0] = 0.5 ; 'height' of gaussian, it's 1/ line depth, rarely more than 0.5 ruff[2] = 3. ; width, in steps. endif ; And use this smoothed cube to determine rough line-center position for all pixels. if (verbose ge 4) then $ progress,0.,/reset,label='Clean_I: line center...',maxval=ncol,bigstep=128,smallstep=32 wavelengths = indgen(nwavel) ; ELW-20110128 fwavelengths = findgen(nwavel) ; ELW-20110128 !p.multi=0 ; find line center for i=0,ncol-1 do begin for j=0,nrow-1 do begin line=reform(bgsmtim[i,j,0:nwavel-2]) ; 6/28/2011 take out that damned continuum point... line = line / max(line) ; gfit=gauss_fit(wavelengths, 1./line,aout) ; line_cent is unstable; gauss_fit is slow. ; ocenter[i,j]=aout[4] ; 6/2011 for speed, trying idl's gaussfit. Indicies of returned aout are different for line-center ; see comments re: fortran call_external ; Just as a reminder, line_cent was unstable, and lc_find ALSO gave step-wise answers, and ; that just propagates badly.... if NOT do_ruff then $ gfit=gaussfit(wavelengths[0:nwavel-2], 1./line,aout) $ else $ gfit=gaussfit(wavelengths[0:nwavel-2], 1./line,aout,est=ruff) if (doplot && ((i mod 50) eq 0) && (j mod 50) eq 0) then begin plot,wavelengths[0:nwavel-2],1./line,tit='FINDING LINE CENTER' oplot,wavelengths[0:nwavel-2],gfit,lines=2 wait,0.1 endif ocenter[i,j]=aout[1] if (aout[1] gt nwavel) then ocenter[i,j]=ocenter[i,j-1]; 6/2011: if crazy, set to previous. ;this may still bomb if you get i,j just right. But I've not had it stop here in testing. endfor if (verbose ge 4) then progress,i+1,frequency=5. endfor if (verbose ge 4) then progress,i,/last ; Next, figure out minimum number of wavelength points around line-center ; that can be used, according to how widely the line-center shifts ; over the frame. Try to catch for pathological puka values. ; 6/2011 last part likely not needed if mask (e.g., 'mm' ) is used, but can't hurt. lim1=min([min(fix(ocenter(mm))),min(round(ocenter(mm)))]) ; minimum rounded line-center position over FOV lim2=max([max(fix(ocenter(mm))),max(round(ocenter(mm)))]) ; maximum rounded line-center position over FOV meancent=mean(ocenter(mm)) stdevcent=stdev(ocenter(mm)) ; just checking next that there are no crazy-shifted points: if (lim1 lt meancent-4*stdevcent) then lim1=round(meancent-4*stdevcent) if (lim2 gt meancent+4*stdevcent) then lim2=round(meancent+4*stdevcent) ; how far away from a default line-center do the spectra ever get?: maxexcursion = (abs(nwavel/2. - lim1) > abs(nwavel/2. - lim2)) wavelim=nwavel/2. - maxexcursion ; number of points used for all fitting, below. if wavelim lt 5 then begin ; this is a hard-wired limit! print,"Line too far to the edge ( < 5 steps of edge), aborting" ; If an error was encountered, set err to 1 (or 2 if telluric was done) err = (do_telluric gt 1) ? do_telluric : 1 return,tim endif ; KDL putting an upper limit on #points used appears to minimize influence from ;rollover and other issues *for Fe line only* (Na is wider intrinsically) if (abs(header.f_wavelength - 630.) lt 1) then wavelim=wavelim < 10 fwavelims = findgen(wavelim-1) ; ELW-20110128 !p.multi=[0,2,1] if (verbose ge 4) then progress,0.,/reset,label='clean_I',maxval=ncol,bigstep=128,smallstep=32 for i=0L,ncol-1L do begin for j=0L,nrow-1L do begin if (mask[i,j]) then begin ; get the spectral line line=reform(bgsmtim[i,j,*]) ; normalize line = line / max(line) ; its center, from earlier icenter=ocenter[i,j] if icenter GE nwavel/5. and icenter LE (4*nwavel)/5. then begin ; New algorithm: make sure line is symmetric about line center. ; 6/2011 notes: New was TRM's. ; Assumes prefilter-function is linear, and symmetric about line-center. ; Assumes solar line is symmetric about line-center. Notes in KDL log book. ; NB: still a problem that line-center was fit for the PF*I, so may be ; off. Ideally could fit for linear part and linecenter at the same time. ; but don't have time to do that now.... ; interpolate onto wavelength-based grid center = interpol(swavel,fwavelengths,icenter) p = swavel-center blue = where(p LE 0,nblue) red = where(p GE 0,nred) ; get wavelength grid for red, blue sides of line, up to the wavelength limit set above. wr = center+fwavelims*wstep wb = center-fwavelims*wstep ; get line-intensity array for red & blue interpolated onto wavelength vs. pixels lr = interpol(line,swavel,wr) ; 6/2011 GB suggest maybe use cubic spline? lb = interpol(line,swavel,wb) ; put wavelength grid for blue, red together. w = [rotate(wb[1:*],2),center,wr[1:*]] ; next line is key. Takes ratios of blue side points to counterparts on red side, ; starting at wings, working in, center point is line-center point divided by itself ; so 1.0 by definition, then working from line-center to wing, the ratio red/blue again. ; the latter has no independent information. Leaving it in for the moment. ; indeed, lb[0]=lr[0] ;lrb = [lb[1:*]/lr[1:*],1.0,lr[1:*]/lb[1:*]] ;TRM comment here:"double correction but fixed below" ; KDL 2011 bug above, correct below lrb = [rotate(lb[1:*]/lr[1:*],2),1.0,lr[1:*]/lb[1:*]] ;double correction but fixed below norder = 1 ; 1 or 2. includes curvature of prefilter if norder=2 ; NO, TRM's method assumes symmetry. norder=2 idea actually due to bug in lrb definition. ; KDL 2/2010 I added the measure_errors ; with a fairly small and wavelength-independent error; including this magically ; removes status .ne. 0 problems.... ; the 1% is hard-wired, but my testing shows no difference in outcome if it's 1%, 10%. res = poly_fit(w,lrb,norder,status=status,meas=fltarr(n_elements(w))+0.01*mean(lrb),yfit=yfit) if (doplot && ((i mod 20) eq 0) && ((j mod 20) eq 0)) then begin plot,w,lrb,psym=4,tit='FITTING FOR SLOPE' oplot,w,yfit,lines=2 endif pstatus = 'Linear' if status NE 0 then begin norder = 1 res=[1.0,0.0] ; no correction pstatus = 'No correction' endif ; status ne 0 ; below, 'slope' isn't just the slope, it's the whole fitted function. slope = res(0,norder) ; norder = 1, this is the slope term from the fit ; use non-sorted wavel here ... ; next line is just a computationally efficient way of computing y=res(0) +res(1)*wavel + res(2)*wavel^2. for iorder = norder-1,0,-1 do slope = slope*wavel + res(0,iorder) cslope = interpol(slope[sswavel],swavel,center) ; slope (which is the fitted function) at center ;slope = cslope+(slope-cslope)/2. ; TRM comment here:"Half correction on blue side, half on red side" ; same thing, easier to read: slope = slope/2. + cslope/2. ; is this the fix to 'double correction' above? no idea. KDL 6/2011. endif else begin ; if center NOT w/in 1/5, 4/5 of range, then: ; old algorithm. Assumes there is continuum on both ends of ; the spectrum or that the line is perfectly centered. w1 = MIN( swavel, i1 ) w2 = MAX( swavel, i2 ) res = POLY_FIT( [w1, w2], [line(i1),line(i2)], 1 ) slope = res(0,0) + res(0,1) * wavel ; use non-sorted wavel here ... pstatus = 'Old Linear' endelse ; if/then of line-center w/in range ; slope = slope / MAX(slope) ; 6/2011 KDL testing without the normalization. ; w/ it, the 'continuum' is not changed (probably continuum, one end of the ; spectrm in any case), and everything else is. W/o this renormalization, ; the line-center doesn't change depth, but the rest of the line both sides ; are modified by the 'prefilter' shape. In many ways, this is what the ; code assumes for its fit for the prefilter (line is symmetric about ; line center) , so w/o the re-normalization would in fact be most consistent ; GB convinced me to remove the renormalization for now, since essentially ; no good reason *to* do it. It's a multiplicative factor, rather than additive, ; so it won't change line depth/etc... if (doplot && ((i mod 20) eq 0) && ((j mod 20) eq 0)) then $ plot,swavel,tim[i,j,sswavel],yrange=[0,continuum*1.25],$ tit='-------: orig, -- -- -- --: new' ; Correct the spectra: oslope[i,j]=reform(res(0,norder)); for returning if desired tim[i,j,*] = tim[i,j,*]/slope ; ... since tim is not sorted. if (doplot && ((i mod 20) eq 0) && ((j mod 20) eq 0)) then begin oplot,swavel,tim[i,j,sswavel],line=2,color=128 oplot,swavel,slope[sswavel]*max(tim[i,j,sswavel]),line=1 xyouts,/normal,0.5,0.9,align=0.5,string(i)+' '+string(j)+' '+pstatus wait,0.5 empty endif ; if plot endif ; if mask [i,j] ne 0 endfor ; j loop if (verbose ge 4) then progress,i+1,frequency=5. endfor ; i loop ; remove final fringing. This is from the second part of ; BJL's ivm_prefilter.pro. ; stop ; im not seeing this do good things for sodium.... ;tim=ivm_fringe(tim,header,mask=mask,fringe=ofringe) ; can return fringe coefficients if you want ; 6/2011 and in fact, I'm not seeing anything hugely positive from this, and occasionally ; it's not doing good things, for Na as well as for Fe cubes. removing for now. ; Remove telluric line at the end by default now if (do_telluric eq 1) then begin ; only if doing after, not before PF stuff, and at all... if (verbose ge 1) then print, ' clean_i: calling ivm_telluric()' tim = IVM_TELLURIC(tim, header); NOTE, scale=1.25 may be good here, don't want to hard-wire it, though endif !p.multi=0 if (verbose ge 4) then progress,i,/last if (verbose ge 1) then print, ' CLEAN_I end: ', SYSTIME() return,tim END