PRO resteli,date
;
;Restoration of stabilized WL images for the instrumental MTF.
;INPUT:  date = 23 or 25
;	 files XXjun/isXX.seq, where XX is the date
;OUTPUT: files XXjun/irXX.seq (224x224), where XX is the date
;
on_error,1

; IMAGE PARAMETERS ----------------------------------------------------

if (date ne 23) and (date ne 25) then  $
   message,'Sorry, I do not know how to handle this dataset'
dat=strtrim(date,2)
nami=dat+'jun/is'+dat+'.'
namo=dat+'jun/ir'+dat+'.'
im=intarr(224,224)

if date eq 23 then nima=320 else nima=323

; PARAMETERS OF RESTORATION ---------------------------------------------

	stp=0.3616	  ; SCALE OF THE IMAGE  (arcsec/pix)
	lambda=1.5542d-4  ; WORKING WAVELENGTH (cm)
	d=44.		  ; DIAMETER OF THE TELESCOPE (cm)
	n=224		  ; SIZE OF THE RESTORED FIELD
	cc=20.		  ; WIENER FILTER PARAMETER C

; 1-D calculation of MTF ------------------------------------------------

a=findgen(n/2+1)/(n*stp)                   ; scale array for plotting
mtf=dblarr(n)
c=lambda*2.062648d5/(d*n*stp)   ; 206264.8 arcsec in radian
ss=0

   for i=1,n do begin
	      
      omega=c*(i-1)
      if omega lt 1.d0 then begin
        tt=acos(omega)-omega*sqrt(1.-omega*omega)
      endif else begin
        ss=ss+1
	if ss eq 1 then print,'MTF cutoff: ',i-1
	tt=0.d0
      endelse
      mtf(i-1)=2.*tt/3.141592654d0

   endfor

; creating the 1-D Wiener filter ----------------------------------------

t1=float((cc+1)*mtf/(cc*mtf*mtf+1))

window,0,xsiz=512,ysiz=384
plot,a,mtf,linestyle=5,xra=[-0.2,3],yra=[-0.2,3],/xst,/yst
oplot,a,t1

; transition to 2-D Wiener filter ---------------------------------------

t=fltarr(n,n)
f1=fltarr(n/2+1,n/2+1)

cros0=min(where(t1 le 0))		; zero crossing point of t1,
					; i.e. we can compute only
for j=0,cros0 do begin			; the non-zero part. The rest
	for i=0,cros0 do begin		; are zeros by definition.
		ii=float(i)
		jj=float(j)
		r=sqrt(ii*ii+jj*jj)
		r1=fix(r)
		if r1 ge cros0 then begin
		  f1(i,j)=0.
		endif else begin
		  f1(i,j)=t1(r1)+(t1(r1+1)-t1(r1))*(r-r1)
		endelse
	endfor
endfor

f2=reverse(f1(1:n/2-1,*),1)

t(0:n/2,0:n/2)=f1			; t - 2-D Wiener filter
t(n/2+1:n-1,0:n/2)=f2
t(n/2+1:n-1,n/2+1:n-1)=reverse(f2(*,1:n/2-1),2)
t(0:n/2,n/2+1:n-1)=reverse(f1(*,1:n/2-1),2)
f1=0
f2=0

;apodization window

   perc=5.		  ; percentage of frame size to be damped
   np=nint(n*perc/100.)
   u=fltarr(n)+1.
   u(0)=0.
   for i=1,np do u(i)=(1.-cos(!pi*i/np))/2.
   u(n-np:n-1)=reverse(u(1:np))
   u=u#u		  	; 2-D apodization window

;loop over images, rebinning, and restoration -------------------------

for i=0,nima do begin

  openr,1,nami+strtrim(i,2)
  readu,1,im
  close,1
  ima=float(im)
  print,i

; RESTORATION ------------------------------------------------------------

; apodization

   mea=mean(ima)
   ima=ima-mea  
   ima=ima*u
   ima=ima+mea		  ; apodized frame

; 2-D Fourier restoration/filtering
   
   f1=fft(ima,-1)
   f1=float(fft(f1*t,1))

; OUTPUT --------------------------------------------------

   openw,1,namo+strtrim(i,2)
   writeu,1,nint(ima)
   close,1

endfor

wdelete,0
END