File: [Development] / JSOC / proj / flatfield / pzt_flat_IDL / flatfield_mask.pro (download) Revision: 1.1, Fri Feb 18 00:21:17 2011 UTC (12 years, 7 months ago) by richard Branch: MAIN CVS Tags: Ver_LATEST, Ver_9-5, Ver_9-41, Ver_9-4, Ver_9-3, Ver_9-2, Ver_9-1, Ver_9-0, Ver_8-8, Ver_8-7, Ver_8-6, Ver_8-5, Ver_8-4, Ver_8-3, Ver_8-2, Ver_8-12, Ver_8-11, Ver_8-10, Ver_8-1, Ver_8-0, Ver_7-1, Ver_7-0, Ver_6-4, Ver_6-3, Ver_6-2, Ver_6-1, Ver_6-0, Ver_5-14, Ver_5-13, HEAD ```IDL package for calculating pzt flatfields 2011.02.17 ```

```pro flatfield_mask, imr, offset, level_spot, level_plage, radius,  masks, flat=flat

nx=(size(imr))[1]
ny=(size(imr))[2]
nl=(size(imr))[3]

lb=fltarr(nx,ny,nl)
apod=lb

sigma=9.0*(nx/1024)
bin=nx/1024

x=rebin(findgen(nx),nx,ny)
y=transpose(rebin(findgen(ny),ny,nx))
;a=[0.340000,      1.37000,     -2.04000,      2.70000,     -1.94000,     0.559000]
;a=[-1.36813,  30.9044, -232.513, 999.353, -2613.42, 4211.68,
;-3921.77, 1512.53, 591.547, -802.067, 226.122]
a=[0.229947, 1.99002,-4.01991, 5.96688, -4.53731, 1.37037]

for k=0, nl-1 do begin
lb[*,*,k]=limb_darkening(nx, radius, offset[k,*], a)  ; limb darkening function: current limb darkening function inaccurate at the limb\
endfor
apodc=apod_circ(nx, 1.09*nx/2.0, 0)
;if not keyword_set(flat) then begin
;c1=readfits('/home/soi/CM/tables/calib/flat/hr/vers_0/1.fits') ;get approximate (large scale) flatfield -- this one is definitely not erfect
;flat=1./c1
;endif

if not keyword_set(flat) then flat=1.0+fltarr(nx,ny)

for k=0L, nl-1 do begin
imsm=smooth(imr[*,*,k]/flat/lb[*,*,k],bin*2)

idx=where(apod[*,*,k] eq 1.0)

print, k, avg(imsm[idx])
threshold=avg(imsm[idx])*level_spot
threshold_plage=avg(imsm[idx])*level_plage

ixx=where((imsm lt threshold or imsm gt threshold_plage) and apod[*,*,k] eq 1.0 and apodc eq 1.0, count)  ; masks out dark regions and bight regions (plages).
mm=replicate(1.0, nx, ny)
nn=n_elements(ixx)
print, nn

if count gt 0 then begin
for i=0L, nn-1 do begin
idd_x=indgen(21)-10+x[ixx[i]]
idd_y=indgen(21)-10+y[ixx[i]]
idd_x=idd_x[where(idd_x ge 0 and idd_x lt nx)]
idd_y=idd_y[where(idd_y ge 0 and idd_y lt ny)]

ibb=transpose(rebin(idd_y,n_elements(idd_y),n_elements(idd_x)))*nx+rebin(idd_x,n_elements(idd_x),n_elements(idd_y))

mm[ibb]=mm[ibb]*(1.0-exp(-((x[ibb]-x[ixx[i]])^2+(y[ibb]-y[ixx[i]])^2)/2.0/sigma))

endfor
endif

; obtain neigborhood for missing pixels
idx=where(mm lt 0.2, ct)
msk=intarr(nx,ny)
;i0x=where(apod[*,*,k] ne 1.0 or apodc ne 1.0)
;msk[i0x]=-1
if ct gt 0 then msk[idx]=1