pro iris_master_fit,dir,wavelim=wavelim,minint=minint,extfac=extfac,outdir=outdir,struct=out ;PURPOSE ;this set of routines follows the technique described in ;Schmit et al 2015 to fit a double gaussian+linear function to ;the Mg II h line. The wavelength and intensity of the profile extrema ;are determined from the best fit model. Expect about 20 profile fits ;per second. ; ;INPUTS ;REQUIRED--- ;DIR-directory that contains the raster FITS files. The ;data can be a single FITS (sit/stare) or multiple. The ;code will determine the window of Mg II h. ; ;OPTIONAL--- ;WAVELIM=(ANGSTROM) how big a window do you want to select surrounding ;h-line ; ;MININT=(DN) minimum intensity used in chi^2. This mitigates ;the effects of large negative values that pad data. ; ;EXTFAC=(INTEGER) factor above the spectral platescale for determining ;the extrema position (extfac=5 give spec position to 0.2*platescale ;accuracy) ; ;OUTPUTS ; An IDL SAV file (use RESTORE) will be ; written to OUTDIR (default same as source FITS) ; The SAV file contains a structure (OUT, can also be optional ; output by setting STRUCT), header (HDR), wavelength vectory (WVL), ; section of WVL used for fit (LIM), and the source directory (DIR). ; ;OUT structure makeup ; OUT key arrays match the size of the input FITS: ; (Position along slit by raster step by raster count) ; ;OUT.TYPE= ; 0=Ignored because not enough data ; 1=Double peaked with h1 minima ; 2=Single peaked with h1 minima ; 3=Single peaked without h1 minima ; 4=Double peaked without h1 minima ; 5=Extrema finding did not work ; 6=Minimization did not converge happily ; ;OUT.CHI=reduced Chi^2 of returned model ;OUT.GOODPIX=how many pixels within LIM are >= MININT ;OUT.FITPARM=9 parameters of fit (see Schmit 2015) ; To plot the profile use: ; prof=IRIS_MGFIT(wvl,out.fitparm(x,y,z,*)) ; ;DOUBLE_[H1V,H1R,H2V,H2R,H3]=if the the profile is identified ; as a double peaked profile *AND* the profile has the labeled extrema ; (see TYPE) then this array store the wavelength [0] and intensity [1] of the ; extrema. Profiles that do not qualify have values set to -100 ; ;SINGLE_[H1V,H1R,PEAK]=same as DOUBLE_* but for single peak ; profiles if ~keyword_set(outdir) then outdir=dir f=file_list(dir,'*raster*.fits') a=iris_obj(f(0)) n=a->getnwin() for i=0,n-1 do begin b=a->getlam(i) if(total(abs(b-2803.5) lt .1) gt 0) then begin gd=a->getline_id(i) wvl=b endif endfor read_iris_l2,f,hdr,dat,wave=gd wavesc=hdr(0).cdelt1 if ~keyword_set(wavelim) then wavelim=[2802.3,2804.6] w=where((wvl gt wavelim(0))*(wvl lt wavelim(1))) lim=minmax(w) n=size(dat) if(n(0) eq 3) then rdat=reform(dat,n(1),n(2),1,n(3)) else rdat=dat sdat=rdat(lim(0):lim(1),*,*,*) swvl=wvl(lim(0):lim(1)) n=size(sdat) out={type:intarr(n(2),n(3),n(4)),chi:fltarr(n(2),n(3),n(4))-100.,goodpix:intarr(n(2),n(3),n(4)),fitparm:fltarr(n(2),n(3),n(4),9)-100.,double_h1v:fltarr(n(2),n(3),n(4),2)-100.,double_h1r:fltarr(n(2),n(3),n(4),2)-100.,double_h2v:fltarr(n(2),n(3),n(4),2)-100.,double_h2r:fltarr(n(2),n(3),n(4),2)-100.,double_h3:fltarr(n(2),n(3),n(4),2)-100.,single_h1v:fltarr(n(2),n(3),n(4),2)-100.,single_h1r:fltarr(n(2),n(3),n(4),2)-100.,single_peak:fltarr(n(2),n(3),n(4),2)-100.} if ~keyword_set(extfac) then extfac=5 extwvl=rebin(swvl,n(1)*extfac) extlam=findgen(n(1)*extfac)/extfac laml=indgen(n(1)*extfac) tmp=min(abs(swvl-2803.53),msp) tmp=min(abs(swvl-2803.08),wsp) tmp=min(abs(swvl-2803.3),pt1) f=fltarr(9) f(0)=median(sdat(wsp,*,0,0)) f(1)=(median(sdat(0,*,0,0))-f(0))/wsp f(2)=msp f(4)=msp f(5)=0.029/wavesc^2 f(6)=.9 f(7)=msp f(8)=0.012/wavesc^2 lims={value:0.,limited:[0,0],limits:[0.,0.]} lims=replicate(lims,9) lims.limited=[ [0,1],[0,0], [1,1], [1,0], [1,1], [0,0], [1,1], [1,1], [1,0]] lims.limits=[[0,f(0)*3],[0,0],[msp-6,msp+6],[0.,0],[msp-12,msp+12],[f(5)/4.,f(5)*4.],[0,20],[msp-8,msp+8],[f(8)/3.,f(8)*3.]] mfit1=fltarr(n(2),n(3),n(4),11) lam=findgen(n(1)) if typename(minint) ne 'FLOAT' then minint=0. for i=0,n(2)-1 do begin for j=0,n(3)-1 do begin for k=0,n(4)-1 do begin y=reform(sdat(*,i,j,k)) gdp=total(y gt minint) if(gdp gt n(1)/2.) then begin er=sqrt(y*18.)/18.+1. er(where(y lt 1))=1e5 ; these points are basically ignored f(3)=2.5*y(pt1) lims.value=f p=mpfitfun('iris_mgfit',lam,y,er,parinfo=lims,perror=g,/quiet,maxiter=400,status=st,ftol=1d-9) val=iris_mgfit(lam,p) chi=total((y-val)^2/er^2)/gdp mom=iris_postfit_gen(lam,val) if(mom(10) ne 1) then begin f(5)=0.06/hdr(0).cdelt1^2 f(8)=0.009/hdr(0).cdelt1^2 lims.limits=[[0,f(0)*3],[0,0],[msp-6,msp+6],[0.,0],[msp-12,msp+12],[f(5)/4.,f(5)*4.],[0,20],[msp-8,msp+8],[f(8)/3.,f(8)*3.]] p1=mpfitfun('iris_mgfit',lam,y,er,parinfo=lims,perror=g,/quiet,maxiter=400,status=st1,ftol=1d-9) val1=iris_mgfit(lam,p1) chi1=total((y-val1)^2/er^2)/gdp if(chi1 lt chi) then begin p=p1 chi=chi1 st=st1 endif endif mom=iris_postfit_gen(laml,iris_mgfit(extlam,p)) ; out.fitparm(i,j,k,*)=p out.fitparm(i,j,k,[2,4,7])=interpol(swvl,lam,p([2,4,7])) out.fitparm(i,j,k,[0,3,6])=p([0,3,6]) out.fitparm(i,j,k,[5,8])=p([5,8])*wavesc^2 out.fitparm(i,j,k,1)=p(1)/wavesc out.chi(i,j,k)=chi out.goodpix(i,j,k)=gdp ;print,i,j,k,st,mom(10) if(st lt 6) then mom(10)=6 case mom(10) of 1: begin out.type(i,j,k)=mom(10) out.double_h1v(i,j,k,0:1)=[extwvl(mom(0)),mom(1)] out.double_h1r(i,j,k,0:1)=[extwvl(mom(2)),mom(3)] out.double_h2v(i,j,k,0:1)=[extwvl(mom(4)),mom(5)] out.double_h2r(i,j,k,0:1)=[extwvl(mom(6)),mom(7)] out.double_h3(i,j,k,0:1)=[extwvl(mom(8)),mom(9)] end 2: begin out.type(i,j,k)=mom(10) out.single_h1v(i,j,k,0:1)=[extwvl(mom(0)),mom(1)] out.single_h1r(i,j,k,0:1)=[extwvl(mom(2)),mom(3)] out.single_peak(i,j,k,0:1)=[extwvl(mom(4)),mom(5)] end 3:begin out.type(i,j,k)=mom(10) out.single_peak(i,j,k,0:1)=[extwvl(mom(4)),mom(5)] end 4:begin out.type(i,j,k)=mom(10) out.double_h2v(i,j,k,0:1)=[extwvl(mom(4)),mom(5)] out.double_h2r(i,j,k,0:1)=[extwvl(mom(6)),mom(7)] out.double_h3(i,j,k,0:1)=[extwvl(mom(8)),mom(9)] end 5:out.type(i,j,k)=mom(10) 6:out.type(i,j,k)=mom(10) endcase endif endfor endfor print,i endfor file=outdir+"/"+strmid(hdr(0).startobs,0,10)+"_"+hdr(0).obsid+".sav" save,wvl,out,hdr,lim,dir,file=file return end