pro genresidcross_sbsp,basepath,outfilname ; Routine adapted from fitsurf.pro, the previous version. Name ; of program was no longer appropriate ; 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 ; Modified for smoothing the raw residual crosstalk image for 2013 ; data because amplitude of residual crosstalk is growing rapidly ; 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 = 7 ; degree of polynomial fit ; first cycle through the images to determine the slit scan position ; order for kk = 0,nftot-1 do begin ; read in the data, make sure overflow etc is compensated readl1_sbsp,pthnam(kk),dat,hdr ; define crosstalk image if kk eq 0 then begin temp = reform(dat(*,*,0)) sz = size(temp) nx = sz(1) & ny = sz(2) nxm1 = nx - 1 ; set up arrays for display bstks = bytarr(nx*8,ny) scoef = dblarr(ndeg+1,3,nftot) ; define working arrays xtalkimg = fltarr(nftot,ny,3) xtalkraw = fltarr(nftot,ny,3) xtalkfit = fltarr(nftot,ny,3) soveri = fltarr(ny,4) endif ; get the slit position of observation from the header sltpos(kk) = sxpar(hdr,'SLITPOS') endfor ; sort the observations in terms of ascending slit position isrt = sort(sltpos) ; test to see if sequence of files is in correct order of ascending slitpos iseq = lindgen(nftot) iorder = 0 idiff = isrt-iseq whrdif = where(isrt-iseq ne 0L,countdd) if countdd gt 0 then begin iorder = 1 sltpos = sltpos(isrt) print,'file order not in ascending slit position. Sequence re-ordered' endif ; start processing loop over files for ksl = 0,nftot-1 do begin if iorder ne 0 then kk = isrt(ksl) else kk = iseq(ksl) ; read in the data, make sure overflow etc is compensated readl1_sbsp,pthnam(kk),dat,hdr temp = reform(dat(*,*,0)) ; 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 for jj = 0,ny-2 do begin itemp = reform(temp(*,jj)) whr = where(itemp gt median(itemp)) avgi = mean_sbsp(itemp(whr)) for istks = 1,3 do begin sdat = reform(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 for istks = 1,3 do begin ; load the xtalk images for jj = 0,ny-1 do begin xtalkraw(kk,jj,istks-1) = soveri(jj,istks) endfor endfor ; end loop over files endfor ; smooth the raw results, the fit with 1-D polynomial at each slit position ; use smoothing width of 10 columns in scan direction, 50 rows along slit yyy = findgen(ny) smth = fltarr(nftot,ny,3) xtcorr = fltarr(nx,ny) for istm1 = 0,2 do begin istks = istm1+1 temp = reform(xtalkraw(*,*,istm1)) stemp = smooth(temp,[10,50]) xtalkimg(*,*,istm1) = stemp for kk = 0,nftot-1 do begin ; use red end of spectrum where there is more continuum xtalk100(kk,istm1) = mean(xtalkimg(99:106,50:150,istm1)) xtalk500(kk,istm1) = mean(xtalkimg(99:106,450:550,istm1)) xtalk900(kk,istm1) = mean(xtalkimg(99:106,850:950,istm1)) endfor endfor ; cycle through the data again to review the results of the correction xtalkr = fltarr(nx,ny) window,1,xs=nx*8,ys=ny for ksl = 0,nftot-1 do begin if iorder ne 0 then kk = isrt(ksl) else kk = iseq(ksl) readl1_sbsp,pthnam(kk),dat,hdr ; setup final result display ; insert the raw data in the first 4 panels for istks = 0,3 do begin 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) endfor ; insert the corrected data in the last 4 panels for istks = 0,3 do begin xstrt = nx*4 + istks*nx if istks eq 0 then begin bstks(xstrt:xstrt+nxm1,*) = bstks(0:nxm1,*) goto,jump37 endif else begin for jj = 0,ny-1 do xtalkr(*,jj) = xtalkimg(kk,jj,istks-1) bstks(xstrt:xstrt+nxm1,*) = $ bytscl( (dat(*,*,istks) -xtalkr*dat(*,*,0)),max=mxx,min=mnn) endelse jump37: endfor tvscl,bstks wait,0.4 endfor ; saving original crosstalk image (xtalkraw), the spatially smoothed image ; (xtalksm), and the image reconstructed from polynomial fits along the slit (smth) save,filename='tempjunk.save',xtalkraw,xtalkimg,smth,braw,bsm,bsmth, $ ndeg,sltpos,nftot ; output the results ;save,filename='coef_xtalk_yy.save',nftot,scoef,sltpos,xtalk,ndeg save,filename='residcrossamples.save',xtalk100, xtalk500, xtalk900,xtalkimg, $ xtalkraw,sltpos,nftot ; create the .geny file savegenx,file=outfilname,xtalkimg,sltpos,nftot return end