pro imap, spdfile ;### Restore the PIA SPD file and extract headers and data arrays: restore, spdfile buf = phtispd.hdr length = n_elements(buf) time = phtispd.time ; passing the primary header for i=0, length-1 do begin keyword = strmid(buf(i), 0, 8) card = strmid(buf(i), 8, 72) if (keyword eq "ATTRDPTS") then y_step = fix(strmid(card,3,20)) if (keyword eq "ATTRDLNS") then z_step = fix(strmid(card,3,20)) if (keyword eq "ATTRNPTS") then npts = fix(strmid(card,3,20)) if (keyword eq "ATTRNLNS") then nlns = fix(strmid(card,3,20)) if (keyword eq "FILENAME") then tdt = strmid(card,7,8) if (keyword eq "OBSERVER") then obs = strmid(card,2,10) if (keyword eq "OBJECT ") then obj = strmid(card,2,10) endfor n=n_elements(phtispd.filt) fwp=phtispd.filt(n/2) if (fwp eq 1 or fwp eq 14) then filter=0 if (fwp eq 2) then filter=90 if (fwp eq 3) then filter=65 if (fwp eq 4) then filter=60 if (fwp eq 5) then filter=80 if (fwp eq 6) then filter=100 if (fwp eq 7) then filter=105 if (fwp eq 8) then filter=90 if (fwp eq 9) then filter=170 if (fwp eq 10) then filter=200 if (fwp eq 11) then filter=180 if (fwp eq 12) then filter=150 if (fwp eq 13) then filter=120 exactdate=phtispd.unit date=strmid(exactdate,0,11) print, ' ' print, 'Processing ......' print, ' Object: ', obj, ' (P.I. is ', obs,')' print, ' TDT: ', tdt print, ' Detector used: ',phtispd.det+'00' print, ' Acquisition date: ',phtispd.unit if (filter eq 0) then print,' No filter used' else $ print, ' Filter: ',filter,' um' print, ' Map NYxNZ = ', npts,' x ',nlns print, ' Y-axis raster step size = ', y_step,' (arcsec)' print, ' Z-axis raster step size = ', z_step,' (arcsec)' print, ' ' ; determine the y and z coordinates of the detector array center ; at each chopper position using the following: ; phtispd.cpos contains the chopper positions (ranging from -90 to +90 degrees) ; For example, the first 30 elements in phtispd.cpos is the following. ; -89.7902 -74.7187 -59.7765 -44.6666 -29.6613 -14.8003 ; 0.380528 15.8756 31.2478 46.4453 61.5392 76.2771 ; 90.3880 -89.7902 -74.7187 -59.7765 -44.6666 -29.6613 ; -14.8003 0.380528 15.8756 31.2478 46.4453 61.5392 ; 76.2771 90.3880 -89.7902 -74.7187 -59.7765 -44.6666 ; -30.1601 ; ; phtispd.rpid contains the raster point IDs n = n_elements(phtispd.cpos) y0=fltarr(n) z0=fltarr(n) for i=0, n-1 do begin y0(i) = (phtispd.rpid(0, i) - 1) * y_step + phtispd.cpos(i) z0(i) = -(phtispd.rpid(1, i) - 1) * z_step endfor ; determine the y and z coordinates of each pixel: ; see explanation at the end of this loop if phtispd.det eq 'C1' then begin a=43.5/2 xbox=[-a,a,a,-a] ybox=[-a,-a,a,a] pnum =9 ; there are 9 pixels in C100 detector chopnum=13 ; and 13 chopper positions per raster point delt=46.0 ; each pixel has dimensions 43.5" x 43.5" ; this is between the pixel centers yout=fltarr(pnum, n) ; y and z values of each pixel for each ch. pl. zout=fltarr(pnum, n) for k=0, n-1 do begin for i=0, 2 do begin for j=0, 2 do begin p = 3*j + i yout(p, k) = delt*(j - 1) + y0(k) zout(p, k) = delt*(i - 1) + z0(k) endfor endfor if k eq 0 then begin yset=yout(*,0) zset=zout(*,0) endif endfor endif else begin if phtispd.det eq 'C2' then begin a=89.4/2 xbox=[-a,a,a,-a] ybox=[-a,-a,a,a] pnum=4 ; there are 4 pixels in C200 detector chopnum=7 ; and 7 chopper positions per raster point delt = 46.0 ; each pixel has dimensions 89.4" x 89.4" ; this equals half the pixel size + half the gap size yout=fltarr(pnum, n) zout=fltarr(pnum, n) for k=0, n-1 do begin for i=0, 1 do begin for j=0, 1 do begin p = 2 - 2*j + i yout(p, k) = delt*(1 - 2*j) + y0(k) zout(p, k) = delt*(2*i - 1) + z0(k) endfor endfor if k eq 0 then begin yset=yout(*,0) zset=zout(*,0) endif endfor endif endelse ; The center of the detector is noted with X = (y0, z0) ; ; The pixel arrangement for C100 detector. ; ; 3 6 9 +Y to the right ; 2 5=X 8 +Z upwards ; 1 4 7 ; ; j i p=3*j+i yout zout ; --------------------------------------- ; 0 0 0 y0-delt z0-delt ; 1 0 3 y0 z0-delt ; 2 0 6 y0+delt z0-delt ; 0 1 1 y0-delt z0 ; 1 1 4 y0 z0 ; 2 1 7 y0+delt z0 ; 0 2 2 y0-delt z0+delt ; 1 2 5 y0 z0+delt ; 2 2 8 y0+delt z0+delt ; ; ; The pixel arrangement for C200 detector. ; ; 3 4 +Y to the right ; X +Z upwards ; 2 1 ; ; j i p=abs(3*i-j) yout zout ; --------------------------------------- ; 0 0 0 y0+delt z0-delt ; 1 0 1 y0-delt z0-delt ; 0 1 3 y0+delt z0+delt ; 1 1 2 y0-delt z0+delt filt=string(filter) filt=strmid(filt,5,3) head=obj+' '+tdt+' '+date+' '+'Filter: '+filt+'!4l!Xm' window, /free, xsize=500, ysize=600, title='IMAP Results' !P.MULTI=[0,1,1] tvLCT,[255,255,0],[0,255,255],[0,0,0] plot, yout, zout, psym=3, xstyle=2, ystyle=2, $ position=[0.1,0.2,0.95,0.95], background=4, $ xtitle='Y-direction (arcsec)', $ ytitle='Z-direction (arcsec)' usersym,xbox,ybox,color=150,/fill oplot,yset,zset,psym=8,symsize=0.15 oplot, yout, zout, psym=3, color=1 oplot, y0, z0, psym=1, color=2 xyouts,0.1,0.97,/normal,head,charsize=1.5 xyouts,0.1,0.1,/normal,'The dots represent the pixel centers.', $ CharSize=1.5, color=1 xyouts,0.1,0.05,/normal,'The plus signs represent the chopper positions.',$ CharSize=1.5, color=2 print,' ' print,'The gray square boxes represent the individual detector pixels.' print,'These boxes are correctly centered at the pixel centers for the' print,'detector position at the top left most location. Note that these' print,'are NOT drawn to scale, and the aspect ratio of the display is not' print,'correct; This is simply to aid the user with the visualization.' print,' ' ;### Divide the map into cells for masking out deviated data points: print, ' ' print, '### Divide the map into cells for masking out deviated data points ###' print, ' ' cleaning: ; define the cell size print, 'The centers of pixels span a region of ' print, max(yout)-min(yout),' arcsecs by',max(zout)-min(zout),' arcsecs.' print, ' ' input: ; set the default to 50" and 100" for C100 and C200 respectively if phtispd.det eq 'C1' then size_cell = 50 if phtispd.det eq 'C2' then size_cell = 100 print, 'Please input a cell size in arcsec' print, '(e.g. 50 for C100 detector or 100 for C200 detector)' read, size_cell Ks=3.0 print, 'Please input a threshold in units of r.m.s. noise for sigma clipping (e.g. 3.0)' read, Ks ; determine the number of cells: ny_cell = fix(((npts-1)*y_step+180)/size_cell) nz_cell = fix( (nlns-1)*z_step /size_cell) if (ny_cell eq 0 or nz_cell eq 0) then begin print, 'The input cell size is too large, please reenter a smaller size.' goto, input endif ; make the number of the cells odd numbers: if ny_cell eq (ny_cell/2*2) then ny_cell = ny_cell+1 if nz_cell eq (nz_cell/2*2) then nz_cell = nz_cell+1 ; determine the central position of the center cell of the map: y_cell=fltarr(ny_cell, nz_cell) z_cell=fltarr(ny_cell, nz_cell) y_cell(ny_cell/2, nz_cell/2)= 0.5*(npts-1)*y_step z_cell(ny_cell/2, nz_cell/2)=-0.5*(nlns-1)*z_step for i=0, ny_cell-1 do begin for j=0, nz_cell-1 do begin y_cell(i,j) = y_cell(ny_cell/2, nz_cell/2) + (i-ny_cell/2)*size_cell z_cell(i,j) = z_cell(ny_cell/2, nz_cell/2) - (j-nz_cell/2)*size_cell endfor endfor wdelete window, /free, xsize=500, ysize=600, title='IMAP Results' !P.MULTI=[0,1,1] tvLCT,[255,255,0],[0,255,255],[0,0,0] plot, yout, zout, psym=3, xstyle=2, ystyle=2, $ position=[0.1,0.2,0.95,0.95], background=4, $ xtitle='Y-direction (arcsec)', $ ytitle='Z-direction (arcsec)' oplot, yout, zout, psym=3, color=1 oplot, y0, z0, psym=1, color=2 oplot, y_cell, z_cell, psym=2, color=0, SymSize=2.0 xyouts,0.1,0.97,/normal,head,charsize=1.5 xyouts,0.1,0.13,/normal,'The dots represent the pixel centers.', $ CharSize=1.5, color=1 xyouts,0.1,0.08,/normal,'The plus signs represent the chopper positions.',$ CharSize=1.5, color=2 xyouts,0.1,0.03,/normal,'The asterisks represent the cell centers.',$ CharSize=1.5, color=0 ; assign a data point (p, n) to a cell mean =fltarr(ny_cell, nz_cell) sigma1=fltarr(ny_cell, nz_cell) flag=phtispd.flag ; flag indicating status of data point for processing ; 0 = valid, 1 = uncertainties not valid, and ; 2 = not valid mnpw=phtispd.mnpw ;mean in-band power / ch. plateau $ raster pos. [W] print, ' ' print, '... cell statatistics:' print, ' iY, iZ, y_center, z_center, no_of_data_pts, mean, r.m.s., no_of_pts_flagged' print, ' (arcsec) (arcsec) (watts) (watts) ' for i=0, ny_cell-1 do begin for j=0, nz_cell-1 do begin sum = 0. sum2 = 0. goodno = 0 no_reject = 0 for p=0, pnum-1 do begin for n1=0, n-1 do begin if (abs(yout(p, n1)-y_cell(i,j)) lt 0.5*size_cell and $ abs(zout(p, n1)-z_cell(i,j)) lt 0.5*size_cell) then begin if flag(p, n1) eq 0 then begin sum = sum + mnpw(p,n1) sum2= sum2+ mnpw(p,n1)*mnpw(p,n1) goodno = goodno + 1 endif endif endfor endfor if goodno ne 0 then begin mean(i,j) = sum/goodno sigma1(i,j)= sum2 - goodno*mean(i,j)*mean(i,j) sigma1(i,j)= sqrt(sigma1(i,j)/(goodno-1)) endif for p=0, pnum-1 do begin for n1=0, n-1 do begin if (abs(yout(p, n1)-y_cell(i,j)) lt 0.5*size_cell and $ abs(zout(p, n1)-z_cell(i,j)) lt 0.5*size_cell) then begin if (abs(mnpw(p, n1)-mean(i,j)) gt Ks*sigma1(i,j) and $ flag(p, n1) eq 0) then begin flag(p, n1) = 3 no_reject=no_reject+1 endif endif endfor endfor print, i, j, y_cell(i, j), z_cell(i, j), goodno, mean(i, j), sigma1(i,j), no_reject endfor endfor primsg: print, ' ' print, 'Would you like to ignore this result and try using different cell' print, 'size and/or different threshold for the sigma-clipping? [yes/no]' print, 'If you want to save this result and move on, please answer no.' answer = 'string' read, answer if (answer eq 'yes' or answer eq 'y' or $ answer eq 'Yes' or answer eq 'Y') then begin set = where(flag eq 3, count) if (count ne 0) then flag(set) = 0 goto, cleaning endif else begin if (answer eq 'no' or answer eq 'n' or $ answer eq 'No' or answer eq 'N') then begin goto, done_clean endif else begin print, '... Oops, not an acceptable answer. Please answer yes or no.' goto, primsg endelse endelse done_clean: ;### Mask out any bright point sources using a sigma-based threshold and a growing radius print, ' ' print, '### masking out bright point sources if there is any ###' print, ' ' print, 'Please wait while processing.' print, ' ' for i=0, ny_cell-1 do begin for j=0, nz_cell-1 do begin sum = 0. goodno = 0 for p=0, pnum-1 do begin for n1=0, n-1 do begin if (abs(yout(p, n1)-y_cell(i,j)) lt 0.5*size_cell and $ abs(zout(p, n1)-z_cell(i,j)) lt 0.5*size_cell) then begin if flag(p, n1) eq 0 then begin sum = sum + mnpw(p,n1) goodno = goodno + 1 endif endif endfor endfor if goodno ne 0 then begin mean(i,j) = sum/goodno endif endfor endfor ; determine the mean of and the r.m.s. around this mean of the cell means: sum=0. sum2=0. for i=0, ny_cell-1 do begin for j=0, nz_cell-1 do begin sum =sum + mean(i,j) sum2=sum2+ mean(i,j)*mean(i,j) endfor endfor meanmean =sum/(ny_cell*nz_cell) meansigma=sum2 - (ny_cell*nz_cell)*meanmean*meanmean meansigma=sqrt(meansigma/(ny_cell*nz_cell-1)) ; now let's mask out the brightest cells: mask: wdelete window,/free,xsize=800,ysize=800,title='IMAP Results' print,' ' print,'In the current display, the dots represent the centers of each pixel' print,'at each chopper plateau, and the X in the lower right plot (not shown' print,'yet) will show the flagged region.' print,'If there is no X present in the lower right plot, then the two' print,'surface plots are exactly the same except for the scale along' print,'the z-direction.' print,' ' !P.MULTI = [0,2,2] loadct,8 shade_surf, mean, y_cell, z_cell, CharSize=2.0, zrange=[0,max(mean)], $ xtitle='Y-direction (arcsec)',ytitle='Z-direction (arcsec)', $ ztitle='Signal [W]' !P.MULTI = [4,2,2] surface, mean, y_cell, z_cell, CharSize=2.0, zrange=[0,max(mean)], $ /NoErase !P.MULTI = [3,2,2] shade_surf,mean, y_cell, z_cell,xrange=[min(yout),max(yout)],$ xtitle='Y-direction (arcsec)',ytitle='Z-direction (arcsec)',$ yrange=[min(zout),max(zout)],xstyle=1,ystyle=1,zstyle=4,Az=0,Ax=90 !P.MULTI = [3,2,2] plot,yout(where(flag eq 0)),zout(where(flag eq 0)),psym=3,$ xstyle=1,ystyle=1,xrange=[min(yout),max(yout)],$ yrange=[min(zout),max(zout)],/noerase !P.MULTI = [3,2,2] contour,mean, y_cell, z_cell,xrange=[min(yout),max(yout)],$ yrange=[min(zout),max(zout)],xstyle=1,ystyle=1,/noerase xyouts,0.5,0.98,/normal,charsize=1.5, $ alignment=0.5,'Masking out bright point source' set_plot,'ps' device,filename='imap_result.ps',$ xsize=7, ysize=10, /inches, xoffset=0.75, yoffset=0.5 !P.MULTI = [0,2,2] shade_surf, mean, y_cell, z_cell, zrange=[0,max(mean)], $ xtitle='Y-direction (arcsec)',ytitle='Z-direction (arcsec)', $ ztitle='Signal [W]' !P.MULTI = [4,2,2] surface, mean, y_cell, z_cell, zrange=[0,max(mean)], $ /NoErase !P.MULTI = [3,2,2] shade_surf,mean, y_cell, z_cell,xrange=[min(yout),max(yout)],$ xtitle='Y-direction (arcsec)',ytitle='Z-direction (arcsec)',$ yrange=[min(zout),max(zout)],xstyle=1,ystyle=1,zstyle=4,Az=0,Ax=90 !P.MULTI = [3,2,2] plot,yout(where(flag eq 0)),zout(where(flag eq 0)),psym=3,$ xstyle=1,ystyle=1,xrange=[min(yout),max(yout)],$ yrange=[min(zout),max(zout)],/noerase !P.MULTI = [3,2,2] contour,mean, y_cell, z_cell,xrange=[min(yout),max(yout)],$ yrange=[min(zout),max(zout)],xstyle=1,ystyle=1,/noerase xyouts,0.5,0.98,/normal,charsize=1.5, $ alignment=0.5,'Masking out bright point source' set_plot,'x' Kc=2.5 XRc = 1.0 print, ' ' print, 'Now, mask out bright point sources if there is any using a sigma-based' print, 'threshold and a growing radius. First, please input a threshold such' print, 'that any cell whose mean signal deviates by more this number of sigma' print, 'will be masked out (e.g. 2.5).' read, Kc print, 'Now, please input a value for the growth radius around these masked' print, 'out cells (if any) using the units of the cell size. Any points' print, 'found within this region will be masked out (e.g. 1.0). If this' print, 'instruction is not clear, please do try different values and study' print, 'the changes in the display result' read, XRc print, ' ' markcell = fltarr(ny_cell,nz_cell) k1 = 0 k2 = 0 Rc = XRc * size_cell for i=0, ny_cell-1 do begin for j=0, nz_cell-1 do begin if (abs(mean(i,j)-meanmean) gt Kc*meansigma) then begin markcell(i,j) = 0.0 for p=0, pnum-1 do begin for n1=0, n-1 do begin if (abs(yout(p, n1)-y_cell(i,j)) lt Rc and $ abs(zout(p, n1)-z_cell(i,j)) lt Rc and $ flag(p, n1) eq 0) then flag(p, n1) = 4 endfor endfor k1=k1+1 endif else begin markcell(i,j)=mean(i,j) k2=k2+1 endelse endfor endfor ; make a surface plot of the results: temp=n_elements(where(flag eq 4)) if temp gt 1 then begin !P.MULTI = [2,2,2] shade_surf, markcell, y_cell, z_cell, CharSize=2.0, $ zrange=[0,max(mean)], $ xtitle='Y-direction (arcsec)',ytitle='Z-direction (arcsec)', $ ztitle='Signal [W]' !P.MULTI = [2,2,2] surface, markcell, y_cell, z_cell, CharSize=2.0, $ zrange=[0,max(mean)], /NoErase !P.MULTI = [1,2,2] shade_surf,markcell, y_cell, z_cell,xrange=[min(yout),max(yout)],$ xtitle='Y-direction (arcsec)',ytitle='Z-direction (arcsec)',$ yrange=[min(zout),max(zout)],xstyle=1,ystyle=1,zstyle=4,Az=0,Ax=90 !P.MULTI = [1,2,2] plot,yout(where(flag eq 0)),zout(where(flag eq 0)),psym=3,$ xstyle=1,ystyle=1,xrange=[min(yout),max(yout)],$ yrange=[min(zout),max(zout)],/noerase oplot,yout(where(flag eq 4)),zout(where(flag eq 4)),psym=7 !P.MULTI = [1,2,2] contour,markcell, y_cell, z_cell,xrange=[min(yout),max(yout)],$ yrange=[min(zout),max(zout)],/noerase,xstyle=1,ystyle=1 !P.MULTI = [0,1,1] endif else begin !P.MULTI = [2,2,2] shade_surf, mean, y_cell, z_cell, CharSize=2.0, $ zrange=[min(mean),max(mean)], $ xtitle='Y-direction (arcsec)',ytitle='Z-direction (arcsec)', $ ztitle='Signal [W]' !P.MULTI = [2,2,2] surface, mean, y_cell, z_cell, CharSize=2.0, zrange=[min(mean),max(mean)], $ /NoErase !P.MULTI = [1,2,2] shade_surf,mean, y_cell, z_cell,xrange=[min(yout),max(yout)],$ xtitle='Y-direction (arcsec)',ytitle='Z-direction (arcsec)',$ yrange=[min(zout),max(zout)],xstyle=1,ystyle=1,zstyle=4,Az=0,Ax=90 !P.MULTI = [1,2,2] plot,yout(where(flag eq 0)),zout(where(flag eq 0)),psym=3,$ xstyle=1,ystyle=1,xrange=[min(yout),max(yout)],$ yrange=[min(zout),max(zout)],/noerase !P.MULTI = [1,2,2] contour,mean, y_cell, z_cell,xrange=[min(yout),max(yout)],$ yrange=[min(zout),max(zout)],xstyle=1,ystyle=1,/noerase endelse set_plot,'ps' if temp gt 1 then begin !P.MULTI = [2,2,2] shade_surf, markcell, y_cell, z_cell, $ zrange=[0,max(mean)], $ xtitle='Y-direction (arcsec)',ytitle='Z-direction (arcsec)', $ ztitle='Signal [W]' !P.MULTI = [2,2,2] surface, markcell, y_cell, z_cell, $ zrange=[0,max(mean)], /NoErase !P.MULTI = [1,2,2] shade_surf,markcell, y_cell, z_cell,xrange=[min(yout),max(yout)],$ xtitle='Y-direction (arcsec)',ytitle='Z-direction (arcsec)',$ yrange=[min(zout),max(zout)],xstyle=1,ystyle=1,zstyle=4,Az=0,Ax=90 !P.MULTI = [1,2,2] plot,yout(where(flag eq 0)),zout(where(flag eq 0)),psym=3,$ xstyle=1,ystyle=1,xrange=[min(yout),max(yout)],$ yrange=[min(zout),max(zout)],/noerase oplot,yout(where(flag eq 4)),zout(where(flag eq 4)),psym=7 !P.MULTI = [1,2,2] contour,markcell, y_cell, z_cell,xrange=[min(yout),max(yout)],$ yrange=[min(zout),max(zout)],/noerase,xstyle=1,ystyle=1 !P.MULTI = [0,1,1] endif else begin !P.MULTI = [2,2,2] shade_surf, mean, y_cell, z_cell, zrange=[min(mean),max(mean)], $ xtitle='Y-direction (arcsec)',ytitle='Z-direction (arcsec)', $ ztitle='Signal [W]' !P.MULTI = [2,2,2] surface, mean, y_cell, z_cell, zrange=[min(mean),max(mean)], $ /NoErase !P.MULTI = [1,2,2] shade_surf,mean, y_cell, z_cell,xrange=[min(yout),max(yout)],$ xtitle='Y-direction (arcsec)',ytitle='Z-direction (arcsec)',$ yrange=[min(zout),max(zout)],xstyle=1,ystyle=1,zstyle=4,Az=0,Ax=90 !P.MULTI = [1,2,2] plot,yout(where(flag eq 0)),zout(where(flag eq 0)),psym=3,$ xstyle=1,ystyle=1,xrange=[min(yout),max(yout)],$ yrange=[min(zout),max(zout)],/noerase !P.MULTI = [1,2,2] contour,mean, y_cell, z_cell,xrange=[min(yout),max(yout)],$ yrange=[min(zout),max(zout)],xstyle=1,ystyle=1,/noerase endelse ;device,/close_file set_plot,'x' ;spawn,'lpr162 imap_result.ps' print, "... total number of cells: ",ny_cell*nz_cell print, "... number of masked out cells: ",k1 print, ' ' ; see if one more iteration is needed: maskmsg: print, 'Would you like to ignore this result and try using different' print, 'threshold and/or growing radius? [yes/no]' print, 'If you like to save this result and move on, please answer no' answer = 'string' read, answer if (answer eq 'no' or answer eq 'n' or $ answer eq 'No' or answer eq 'N') then goto, done_mask else $ if (answer eq 'yes' or answer eq 'y' or $ answer eq 'Yes' or answer eq 'Y') then begin set = where(flag eq 4, count) if count ne 0 then flag(set) = 0 goto, mask endif else begin print, '... oops, not an acceptable answer. Please answer yes/no.' goto, maskmsg endelse done_mask: ;### Chooper vignetting correction: print, ' ' print, '### statistically determine chopper vignetting ###' cpv = fltarr(pnum, chopnum, n) cpvign = fltarr(pnum, chopnum) m = intarr(pnum, chopnum) for p=0, pnum-1 do begin ; for each pixel (9 for C100) for n1=0, n-1 do begin ; for each chopper plateau if (flag(p, n1) eq 0) then begin ; only for data points with flag=0 i = phtispd.cstp(n1)-1 ; chopper step index to count from 0 cpv(p, i, m(p, i)) = mnpw(p, n1) ; mean in-band power ; per pixel ; per chopper plateau m(p, i) = m(p, i) + 1 endif endfor endfor for p=0, pnum-1 do begin for i=0, chopnum-1 do begin if (m(p, i) ne 0.) then $ cpvign(p, i) = median(cpv(p, i, 0:m(p,i)-1)) ; median of in-band power per pixel per chopper throw endfor endfor ; calculate the relative vignetting correction: for p=0, pnum-1 do begin center=cpvign(p, chopnum/2) for k1=0, chopnum-1 do begin if (center ne 0.) then $ cpvign(p, k1)=cpvign(p, k1)/center endfor endfor if (phtispd.det eq 'C1') then begin print,' ' print,'The pixel #1 : the top left plot' print,'The pixel #2 : the top middle plot' print,'The pixel #3 : the top right plot' print,'The pixel #4 : the middle left plot' print,'The pixel #5 : the center plot' print,'The pixel #6 : the middle right plot' print,'The pixel #7 : the lower left plot' print,'The pixel #8 : the lower middle plot' print,'The pixel #9 : the lower right plot' endif if (phtispd.det eq 'C2') then begin print,' ' print,'The pixel #1 : the top left plot' print,'The pixel #2 : the top right plot' print,'The pixel #3 : the lower left plot' print,'The pixel #4 : the lower right plot' endif print, ' ' ; make polynomial fits and display the results: wdelete window, /free, xsize=800, ysize=800, title='IMAP Results' tvLCT,[255,255,0],[0,255,255],[0,0,0] if (phtispd.det eq 'C1') then begin !P.MULTI = [0,3,3] cs=1.4 xx=[-90,-75,-60,-45,-30,-15,0,15,30,45,60,75,90] endif if (phtispd.det eq 'C2') then begin !P.MULTI = [0,2,2] cs=1.0 xx=[-90,-60,-30,0,30,60,90] endif ;xx = findgen(chopnum) + 1 for p=0, pnum-1 do begin yy = cpvign(p, *) ; perform least-square polynomial fit result = poly_fit(xx, yy, 2, yfit2) ; 2nd degree polynomial fit result = poly_fit(xx, yy, 3, yfit3) ; 3rd degree polynomial fit result = poly_fit(xx, yy, 4, yfit4) ; 4th degree polynomial fit plot, xx, yy, psym=4, xrange=[min(xx)-10, max(xx)+10], $ yrange=[min(yy)-0.05, max(yy)+0.05], $ xstyle=1, ystyle=1, charsize=cs, $ title ='solid line (order=2), dotted (order=3), dashed (order=4)', $ xtitle='chopper position', ytitle='normalized vignetting factor', $ background=3 oplot, xx, yy, psym=4, color=1 oplot, xx, yfit2, line=0, color=1 oplot, xx, yfit3, line=1, color=2 oplot, xx, yfit4, line=2, color=0 ndegree = 2 print, 'For pixel ',p+1,', please select the order of fit (2 to 4 ):' read, ndegree result = poly_fit(xx, yy, ndegree, yfit) cpvign(p, *) = yfit endfor ; finally, make the correction if the user choose so: print, ' ' print, 'Would you like to apply these fittings to the data? [yes/no] ' answer = 'string' read, answer if (answer eq 'yes' or answer eq 'y' or $ answer eq 'Yes' or answer eq 'Y') then begin for p=0, pnum-1 do begin for n1=0, n-1 do begin i = phtispd.cstp(n1)-1 if (cpvign(p, i) ne 0.) then $ mnpw(p, n1) = mnpw(p, n1)/cpvign(p, i) endfor endfor print,' ... done the statistical vignetting correction' endif else begin print,' ... skipped the vignetting correction here' endelse ;### Long-term drift correction using a moving median method: print, ' ' print, '### long-term drift correction using a moving median method ###' mmedian: delta=2*n/3 print, '... the length of the time series = ', n print, '... the suggested width for the moving median = ', delta print, 'Please enter your choice of the width for median (caution: do not make it too small):' read, delta if delta ne (delta/2*2) then delta = delta + 1 median1=fltarr(pnum, n) for p=0, pnum-1 do begin value = 0. set=where(flag(p, 0:delta/4-1) eq 0, count) if(count ne 0) then value = median(mnpw(p,set)) for n1=0, delta/4-1 do begin median1(p, n1)= value endfor xx = fltarr(delta/2) for n1=delta/4, delta/2-1 do begin set=where(flag(p, n1-delta/4:n1-1+delta/4) eq 0, count) if (count ne 0) then begin xx = mnpw(p, n1-delta/4:n1-1+delta/4) median1(p, n1)=median(xx(set)) endif endfor xx = fltarr(delta) for n1=delta/2, n-delta/2 do begin set=where(flag(p, n1-delta/2:n1-1+delta/2) eq 0, count) if (count ne 0) then begin xx = mnpw(p, n1-delta/2:n1-1+delta/2) median1(p, n1)=median(xx(set)) endif endfor xx = fltarr(delta/2) for n1=n-delta/2+1,n-delta/4 do begin set=where(flag(p, n1-delta/4:n1+delta/4-1) eq 0, count) if (count ne 0) then begin xx = mnpw(p, n1-delta/4:n1+delta/4-1) median1(p, n1)=median(xx(set)) endif endfor xx = fltarr(delta/4) value = 0.0 set=where(flag(p, n-delta/4+1:n-1) eq 0, count) if(count ne 0) then begin xx = mnpw(p, n-delta/4+1:n-1) value = median(xx(set)) endif for n1=n-delta/4+1, n-1 do begin median1(p, n1)=value endfor endfor print, ' ' print, '==> plot the result' ; plot the results: wdelete window, /free, xsize=800, ysize=800, title='IMAP Results' tvLCT,[255,255,0],[0,255,255],[0,0,0] if (phtispd.det eq 'C1') then !P.MULTI = [0,3,3] if (phtispd.det eq 'C2') then !P.MULTI = [0,2,2] for p=0, pnum-1 do begin plot, time, mnpw(p,*), psym=3, charsize=1.8, background=3, $ xstyle=1,ystyle=2 set=where(flag(p, 0:n-1) ge 0 and flag(p,0:n-1) le 1, count) if (count ne 0) then oplot, time(set), mnpw(p,set), psym=3, color=2 set=where(flag(p, 0:n-1) eq 3, count) if (count ne 0) then oplot, time(set), mnpw(p,set), psym=7, color=0 set=where(flag(p, 0:n-1) eq 2, count) if (count ne 0) then oplot, time(set), mnpw(p,set), psym=3, color=0 oplot, median1(p,*), line = 0 set=where(flag(p, 0:n-1) eq 4, count) if (count ne 0) then oplot, time(set), mnpw(p,set), psym=1, color=1 endfor set_plot,'ps' ;device,filename='imap_result.ps',xsize=7,ysize=10,$ ; xoffset=0.75,yoffset=0.5,/inches for p=0, pnum-1 do begin plot, time,mnpw(p,*), psym=3, background=3, $ xstyle=1,ystyle=2,$ xtitle='Time [s]', ytitle='Signal [W]' set=where(flag(p, 0:n-1) ge 0 and flag(p,0:n-1) le 1, count) if (count ne 0) then oplot, time(set), mnpw(p,set), psym=3, color=2 set=where(flag(p, 0:n-1) eq 3, count) if (count ne 0) then oplot, time(set), mnpw(p,set), psym=1, color=0 set=where(flag(p, 0:n-1) eq 2, count) if (count ne 0) then oplot, time(set), mnpw(p,set), psym=3, color=0 oplot, median1(p,*), line = 0 set=where(flag(p, 0:n-1) eq 4, count) if (count ne 0) then oplot, time(set), mnpw(p,set), psym=4, color=1 endfor ;device,/close_file set_plot,'x' ;spawn,'lpr162 imap_result.ps' print,' ' print,'----- LEGEND -----' print,'The x-axis = Time [s]' print,'The y-axis = Signal [W]' print,'The data points with flag = 0 and 1 are the green dots.' print,'The masked out bright source is in yellow.' print,'The flagged points using the cells are the red Xs.' print,'The data points with flag = 2 (not valid) are red dots print,'The moving median just calculated is in white.' print,' ' if (phtispd.det eq 'C1') then begin !P.MULTI = [0,3,3] cs=1.4 endif if (phtispd.det eq 'C2') then begin !P.MULTI = [0,2,2] cs=1.0 endif print, 'Do you want to zoom in on the moving median curves? [yes/no]' answer = 'string' read, answer if (answer eq 'yes' or answer eq 'y' or $ answer eq 'Yes' or answer eq 'Y') then begin for p=0, pnum-1 do begin plot, time, median1(p,*), line=0, $ ystyle=1, background=3, CharSize=cs, xstyle=1 oplot, time, mnpw(p,*), psym=3 set=where(flag(p, 0:n-1) ge 0 and flag(p,0:n-1) lt 3, count) if (count ne 0) then oplot, time(set), mnpw(p,set), psym=3, color=2 set=where(flag(p, 0:n-1) eq 3, count) if (count ne 0) then oplot, time(set), mnpw(p,set), psym=7, color=0 set=where(flag(p, 0:n-1) eq 2, count) if (count ne 0) then oplot, time(set), mnpw(p,set), psym=3, color=0 print,' execute this' set=where(flag(p, 0:n-1) eq 4, count) if (count ne 0) then oplot, time(set), mnpw(p,set), psym=1, color=1 endfor endif ; now smooth the medians: print, 'Do you want to heavily smooth these curves? [yes/no]' ; This smooths using a box car average of width, delta/5. ; Each time the user is prompted, the curve has been smoothed 5 times. answer = 'string' read, answer if (answer eq 'yes' or answer eq 'y' or $ answer eq 'Yes' or answer eq 'Y') then begin sdelta =delta/5 smedian= median1 for k=0, 200 do begin for p=0, pnum-1 do begin smedian(p,*) = smooth(smedian(p,*), sdelta) if((k mod 5) eq 0 and k gt 0) then begin plot, time, median1(p,*), line=0, ystyle=1, xstyle=1, $ background=3, CharSize=1.8 set=where(flag(p, 0:n-1) ge 0 and flag(p,0:n-1) lt 3, count) if (count ne 0) then $ oplot, time(set), mnpw(p,set), psym=3, color=2 set=where(flag(p, 0:n-1) eq 3, count) if (count ne 0) then $ oplot, time(set), mnpw(p,set), psym=7, color=0 set=where(flag(p, 0:n-1) eq 2, count) if (count ne 0) then $ oplot, time(set), mnpw(p,set), psym=3, color=0 set=where(flag(p, 0:n-1) eq 4, count) if (count ne 0) then $ oplot, time(set), mnpw(p,set), psym=1, color=1 oplot, time, smedian(p,*), line=2 endif endfor if((k mod 5) eq 0 and k gt 0) then begin print, 'Do you want to continue smoothing? [yes/no] read, answer if (answer eq 'no' or answer eq 'n' or $ answer eq 'No' or answer eq 'N') then goto, smooth_done endif endfor smooth_done: median1 = smedian endif subsigma=1.0 print,'Please input a number of sigma below the smoothed curve such that' print,' the data points below it should be flagged. (The suggest value is' print,' 1.0.) Use a large value if you want to skip this clipping.' read,subsigma for p=0, pnum-1 do begin plot, time, mnpw(p,*), psym=3, charsize=1.8, background=3, $ xstyle=1,ystyle=2 set=where(flag(p, 0:n-1) eq 4, count) if (count ne 0) then oplot, time(set), mnpw(p,set), psym=3, color=2 set=where(flag(p, 0:n-1) eq 2, count) if (count ne 0) then oplot, time(set), mnpw(p,set), psym=3, color=0 set=where(flag(p, 0:n-1) eq 3, count) if (count ne 0) then oplot, time(set), mnpw(p,set), psym=7, color=0 oplot, time, median1(p,*), line = 0 set=where(flag(p, 0:n-1) ge 0 and flag(p,0:n-1) le 1, count) if (count ne 0) then begin oplot, time(set), mnpw(p,set), psym=3, color=2 set1=where(flag(p, 0:n-1) eq 4, count) if (count ne 0) then begin min=min(where(flag(p,0:n-1) eq 4)) max=max(where(flag(p,0:n-1) eq 4)) value=where(set gt min and set lt max and $ mnpw(p,set) lt median1(p,set)-subsigma*meansigma, count) if (count ne 0) then begin set=set(where(set gt min and set lt max and $ mnpw(p,set) lt median1(p,set)-subsigma*meansigma)) oplot, time(set), mnpw(p,set), psym=3, color=0 flag(p,set)=5 endif endif endif endfor set_plot,'ps' ;device,filename='imap_result.ps',/inches,xsize=7,ysize=10,$ ; xoffset=0.75,yoffset=0.5 for p=0, pnum-1 do begin set=where(flag(p, 0:n-1) ge 0 and flag(p,0:n-1) le 1, count) plot, time(set), mnpw(p,set), psym=3, background=3, $ xstyle=1,ystyle=2,$ xtitle='Time [s]', ytitle='Signal [W]' oplot, time, median1(p,*), line = 0, color=1 endfor ;device,/close_file set_plot,'x' ;spawn,'lpr162 imap_result.ps' ; make the correction to the data set if desired by the user answer = 'string' print, 'Please enter one of the following options: ' print, ' enter 1 to accept these fits (i.e., the moving median curves) as the long-term drifts,' print, ' enter 2 to redo the fits using a different width, or ' print, ' enter 0 to skip this long-term drift correction all together:' read, answer if (answer eq '1') then begin for p=0, pnum-1 do begin lastmean=median(median1(p, n-1-delta:n-1)) if (lastmean ne 0.) then begin for n1=0, n-1 do begin median1(p, n1)=median1(p, n1)/lastmean endfor endif endfor for p=0, pnum-1 do begin for n1=0, n-1 do begin mnpw(p, n1)=mnpw(p, n1)/median1(p, n1) endfor endfor endif if (answer eq '2') then goto, mmedian ;### Skyflattening the data if needed: print, ' ' print, '### skyflattening ###' y_min=min(yout) y_max=max(yout) z_min=min(zout) z_max=max(zout) ; get the normalization factor from sky median per pixel: skyflat=fltarr(pnum) for p=0, pnum-1 do begin set=where(yout(p,*) ge y_min+2*delt and yout(p,*) le y_max-2*delt and $ zout(p,*) ge z_min+2*delt and zout(p,*) le z_max-2*delt and $ flag(p, *) eq 0, count) if (count ne 0) then begin skyflat(p)=median(mnpw(p, set)) endif else begin print,' no skyflat defined for pixel ', p endelse endfor meanskyflat=total(skyflat)/pnum for p=0, pnum-1 do begin if (skyflat(p) ne 0.0) then skyflat(p)=meanskyflat/skyflat(p) endfor ; plot the pixel layout and skyflat factors on screen: print, 'Pixel Layout | Multiplying Factor from Skyflat ' if (phtispd.det eq 'C1') then begin print, ' 9 6 3 ', skyflat(8), skyflat(5), skyflat(2) print, ' 8 5 2 ', skyflat(7), skyflat(4), skyflat(1) print, ' 7 4 1 ', skyflat(6), skyflat(3), skyflat(0) endif if (phtispd.det eq 'C2') then begin print, ' 4 3 ', skyflat(3), skyflat(2) print, ' 1 2 ', skyflat(0), skyflat(1) endif ; see if the correction is to be applied to the data: answer = 'string' print, 'Do you want to apply this skyflat to the data? [yes/no]' read, answer if (answer eq 'yes' or answer eq 'y' or $ answer eq 'Yes' or answer eq 'Y') then begin for p=0, pnum-1 do begin for n1=0, n-1 do begin mnpw(p, n1)=mnpw(p, n1)*skyflat(p) endfor endfor endif ; plot the final result print, ' ' print, '==> plot the final result' if (phtispd.det eq 'C1') then begin !P.MULTI = [0,3,3] cs=1.4 endif if (phtispd.det eq 'C2') then begin !P.MULTI = [0,2,2] cs=1.0 endif for p=0, pnum-1 do begin plot, time, mnpw(p,*), psym=3, charsize=cs, background=3, $ xstyle=1,ystyle=2 set=where(flag(p, 0:n-1) ge 0 and flag(p, 0:n-1) le 1, count) if (count ne 0) then oplot, time(set), mnpw(p,set), psym=3, color=1 set=where(flag(p, 0:n-1) eq 3, count) if (count ne 0) then oplot, time(set), mnpw(p,set), psym=3, color=0 set=where(flag(p, 0:n-1) eq 5, count) if (count ne 0) then oplot, time(set), mnpw(p,set), psym=3, color=0 set=where(flag(p, 0:n-1) eq 2, count) if (count ne 0) then oplot, time(set), mnpw(p,set), psym=3, color=0 set=where(flag(p, 0:n-1) eq 4, count) if (count ne 0) then oplot, time(set), mnpw(p,set), psym=3, color=1 endfor set_plot,'ps' ;device,filename='imap_result.ps',xsize=7,ysize=10,/inches,$ ; xoffset=0.75,yoffset=0.5 for p=0, pnum-1 do begin plot, time, mnpw(p,*), psym=3, background=3, $ xstyle=1,ystyle=2,$ xtitle='Time [s]', ytitle='Signal [W]' set=where(flag(p, 0:n-1) ge 0 and flag(p,0:n-1) le 1, count) if (count ne 0) then oplot, time(set), mnpw(p,set), psym=3, color=1 set=where(flag(p, 0:n-1) eq 3, count) if (count ne 0) then oplot, time(set), mnpw(p,set), psym=7, color=0 set=where(flag(p, 0:n-1) eq 5, count) if (count ne 0) then oplot, time(set), mnpw(p,set), psym=7, color=0 set=where(flag(p, 0:n-1) eq 2, count) if (count ne 0) then oplot, time(set), mnpw(p,set), psym=3, color=0 set=where(flag(p, 0:n-1) eq 4, count) if (count ne 0) then oplot, time(set), mnpw(p,set), psym=3, color=1 endfor ;device,/close_file set_plot,'x' ;spawn,'lpr162 imap_result.ps' print, ' ' print, '----- LEGEND -----' print, 'The x-axis = Time [s]' print, 'The y-axis = Signal [W]' print, 'The yellow dots represent the valid data points.' print, 'The red dots respresent the flagged points.' print, ' ' ; Final 3D plot of the result with the given cell size. window,/free,xsize=600,ysize=600,title='IMAP Results' !P.MULTI = [0,2,2] shade_surf, mean, y_cell, z_cell, CharSize=2.0, $ zrange=[0,max(mean)], $ background=3, title='Before...',$ xtitle='Y-direction (arcsec)',ytitle='Z-direction (arcsec)',$ ztitle='Signal [W]' !P.MULTI = [4,2,2] surface, mean, y_cell, z_cell, CharSize=2.0, zrange=[0,max(mean)], $ /NoErase !P.MULTI = [3,2,2] shade_surf,mean, y_cell, z_cell,xrange=[min(y0),max(y0)],$ yrange=[min(z0),max(z0)],xstyle=1,ystyle=1,zstyle=4,Az=0,Ax=90,$ xtitle='Y-direction (arcsec)',ytitle='Z-direction (arcsec)' !P.MULTI = [3,2,2] contour,mean, y_cell, z_cell,xrange=[min(y0),max(y0)],$ yrange=[min(z0),max(z0)],xstyle=1,ystyle=1,/noerase set_plot,'ps' ;device,filename='imap_result.ps',$ ; xsize=7, ysize=9, /inches, xoffset=0.75, yoffset=1 !P.MULTI = [0,2,2] shade_surf, mean, y_cell, z_cell, $ zrange=[0,max(mean)], $ title='Before...',$ xtitle='Y-direction (arcsec)',ytitle='Z-direction (arcsec)',$ ztitle='Signal [W]' !P.MULTI = [4,2,2] surface, mean, y_cell, z_cell, zrange=[0,max(mean)], $ /NoErase !P.MULTI = [3,2,2] shade_surf,mean, y_cell, z_cell,xrange=[min(y0),max(y0)],$ yrange=[min(z0),max(z0)],xstyle=1,ystyle=1,zstyle=4,Az=0,Ax=90,$ xtitle='Y-direction (arcsec)',ytitle='Z-direction (arcsec)' !P.MULTI = [3,2,2] contour,mean, y_cell, z_cell,xrange=[min(y0),max(y0)],$ yrange=[min(z0),max(z0)],xstyle=1,ystyle=1,/noerase set_plot,'x' print,'please wait for the final plots' for i=0, ny_cell-1 do begin for j=0, nz_cell-1 do begin sum = 0. goodno = 0 for p=0, pnum-1 do begin for n1=0, n-1 do begin if (abs(yout(p, n1)-y_cell(i,j)) lt 0.5*size_cell and $ abs(zout(p, n1)-z_cell(i,j)) lt 0.5*size_cell) then begin if (flag(p, n1) eq 0 or flag(p, n1) eq 1 or flag(p,n1) eq 4) $ then begin sum = sum + mnpw(p,n1) goodno = goodno + 1 endif endif endfor endfor if goodno ne 0 then mean(i,j) = sum/goodno endfor endfor !P.MULTI = [2,2,2] shade_surf, mean, y_cell, z_cell, CharSize=2.0, $ zrange=[0,max(mean)], $ background=3, title='After...',$ xtitle='Y-direction (arcsec)',ytitle='Z-direction (arcsec)',$ ztitle='Signal [W]' !P.MULTI = [2,2,2] surface, mean, y_cell, z_cell, CharSize=2.0, zrange=[0,max(mean)], $ /NoErase !P.MULTI = [1,2,2] shade_surf,mean, y_cell, z_cell,xrange=[min(y0),max(y0)],$ yrange=[min(z0),max(z0)],xstyle=1,ystyle=1,zstyle=4,Az=0,Ax=90,$ xtitle='Y-direction (arcsec)',ytitle='Z-direction (arcsec)' !P.MULTI = [1,2,2] contour,mean, y_cell, z_cell,xrange=[min(y0),max(y0)],$ yrange=[min(z0),max(z0)],xstyle=1,ystyle=1,/noerase set_plot,'ps' !P.MULTI = [2,2,2] shade_surf, mean, y_cell, z_cell, $ zrange=[0,max(mean)], $ title='After...',$ xtitle='Y-direction (arcsec)',ytitle='Z-direction (arcsec)',$ ztitle='Signal [W]' !P.MULTI = [2,2,2] surface, mean, y_cell, z_cell, zrange=[0,max(mean)], $ /NoErase !P.MULTI = [1,2,2] shade_surf,mean, y_cell, z_cell,xrange=[min(y0),max(y0)],$ yrange=[min(z0),max(z0)],xstyle=1,ystyle=1,zstyle=4,Az=0,Ax=90,$ xtitle='Y-direction (arcsec)',ytitle='Z-direction (arcsec)' !P.MULTI = [1,2,2] contour,mean, y_cell, z_cell,xrange=[min(y0),max(y0)],$ yrange=[min(z0),max(z0)],xstyle=1,ystyle=1,/noerase device,/close_file set_plot,'x' ;spawn,'lpr162 imap_result.ps' ;### Finishing up: print, ' ' print, '### saving the results to a PIA/SPD file ###' answer = 'string' print, 'Do you want to save the results to a new PIA/SPD file? [yes/no]' read, answer if (answer eq 'y' or answer eq 'yes' or $ answer eq 'Y' or answer eq 'Yes') then begin ; set the flags back to 0's and 2's: for p=0, pnum-1 do begin for n1=0, n-1 do begin if (flag(p, n1) eq 3) then flag(p, n1)=2 if (flag(p, n1) eq 4) then flag(p, n1)=0 if (flag(p, n1) eq 5) then flag(p, n1)=2 endfor endfor ; replace the old data arrays: phtispd.mnpw = mnpw phtispd.flag = flag ; write the updated results to the save file: spdfile=spdfile+'_IMAP' save, phtispd, filename=spdfile print,'... warning: the new SPD file with "_IMAP" ending is created!' endif else begin print,'... warning: the results are NOT saved!' endelse print, ' ' print, 'DONE' print, ' ' END