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