pro fitsurf,basepath,outfilname ; routine modified from previous versions of fitsurf, 15 Aug 2011 ; This routine fits the variation of the continuum crosstalk in ; the direction along the slit only. Present version fits with ; a 4th-order polynomial (ndeg = 4) ; Smooths the variations of the fit parameters with a running ; boxcar using the IDL smooth proceedure ; INPUT: ; basepath = pathname to level1 data calibrated without ; the residual crosstalk correction ; outfilname = the name for the curent output file, usually named ; after the date of the flat field operation used to generate ; the crosstalk data, for example residxtalk_20120702_114140 ; revised version to use the 2012 eclipse season flat results ; MUST BE RUN UNDER SSW! ; OUTPUT: ; scoef = array of coefficients(3 Stokes,ndeg,nfiles) ; outputs .geny file to be used by SP_PREP ; This routine avoids the locations around the lines where ; the polarization might bias the fit ; Also modified Nov. 2012 to avoid fitting vignetted regions at upper ; and lower ends of slit ; determines surface fit coefficients of Q/I, U/I, V/I ; for the full array ; get the list of input files derf = file_search(basepath + '*.fits') nftot = n_elements(derf) print,'directory: ',basepath,'; number of files = ',nftot pthnam = derf ; output array containing slit positions of data ; The variables below are for checking the results only sltpos = fltarr(nftot) xtalk = fltarr(nftot,3) xtalk100 = xtalk & xtalk500 = xtalk & xtalk900 = xtalk ndeg = 4 ; degree of polynomial fit ; start loop over files for kk = 0,nftot-1 do begin ; read in the data, make sure overflow etc is compensated readl1_sbsp,pthnam(kk),dat,hdr temp = dat(*,*,0) ;nx = sizeof(temp,1) & ny = sizeof(temp,2) sz = size(temp) nx = sz(1) & ny = sz(2) nxm1 = nx - 1 ; get the slit position of observation from the header sltpos(kk) = sxpar(hdr,'SLITPOS') if kk eq 0 then begin ; set up arrays for display bstks = bytarr(nx*8,ny) window,1,xs=nx*8,ys=ny scoef = dblarr(ndeg+1,3,nftot) endif ;; display raw data first bstks(0:nxm1,*) = bytscl(dat(*,*,0)) ; now make variations of Q/I, U/I, V/I along slit for continuum ; pixels soveri = fltarr(ny,4) for jj = 0,ny-2 do begin itemp = temp(*,jj) whr = where(itemp gt median(itemp)) avgi = mean_sbsp(itemp(whr)) for istks = 1,3 do begin sdat = dat(*,jj,istks) sdat = sdat(whr) soveri(jj,istks) = mean_sbsp(sdat)/avgi endfor ; put mean intensity into the first slot soveri(jj,0) = avgi endfor ; loop through Stokes QUV parameters, get fit coefficients along slit yyy = findgen(ny) for istks = 1,3 do begin ; bypass 20 points at each end of the slit because of possible vignetting scoef(*,istks-1,kk) = poly_fit(yyy(20:ny-20),soveri(20:ny-20,istks),ndeg) ; recover the fit smth = fltarr(nx,ny) smth(*,*) = 0. for jj = 0,ny-1 do begin for mj = 0,ndeg do begin smth(*,jj) = smth(*,jj) + $ scoef(mj,istks-1,kk) * yyy(jj)^mj endfor endfor ; use red end of spectrum where there is more continuum xtalk(kk,istks-1) = mean(smth(99:106,10:1000)) xtalk100(kk,istks-1) = mean(smth(99:106,50:150)) xtalk500(kk,istks-1) = mean(smth(99:106,450:550)) xtalk900(kk,istks-1) = mean(smth(99:106,850:950)) ; setup final result display ; insert the raw data in the first 4 columns xstrt = istks*nx mxx = max(dat(*,*,istks)) & mnn = min(dat(*,*,istks)) if istks gt 0 then mxx = 100. & mnn = -100. bstks(xstrt:xstrt+nxm1,*) = bytscl(dat(*,*,istks),max=mxx,min=mnn) xstrt = nx*4 + istks*nx bstks(xstrt:xstrt+nxm1,*) = $ bytscl( (dat(*,*,istks) -smth*dat(*,*,0)),max=mxx,min=mnn) ; end loop over Stokes parameters endfor ;window,1,xs=nx*8,ys=ny tvscl,bstks ; end loop over files endfor ; sort the observations in terms of ascending slit position isrt = sort(sltpos) sltpos = sltpos(isrt) scoef = scoef(*,*,isrt) ; do boxcar averaging in the slit scan direction ; optimal smoothing seems to be about 200 scan steps ; determine smoothing parameter for the present data nsmoothb = fix( 200.*nftot/(max(sltpos)-min(sltpos)) ) print,' boxcar smoothing parameter in measured slit position steps: ',nsmoothb ; smooth the fit coefficients ; bypass smoothing because small changes in high order terms will ruin the fit goto,jump453 for istks = 0,2 do begin for jj = 0,ndeg-1 do begin temp = smooth(scoef(jj,istks,*),nsmoothb,/edge_truncate) scoef(jj,istks,*) = temp endfor endfor jump453: ; output the results save,filename='coef_xtalk_yy.save',nftot,scoef,sltpos,xtalk,ndeg save,filename='residcrossamples.save',xtalk100, xtalk500, xtalk900 ; create the .geny file savegenx,file=outfilname,scoef,ndeg,sltpos,nftot return end