PRO restelec_f95,init,nim,x0,y0,xdim,ydim

; RESTORATION of an image (2-D frame) for
; the 2-D theoretical MTF of a diffraction-limited telescope.
; The input and output images are read/written fm/to the disk
; (The working array must have xdim = ydim = even number)
; calling:   RESTELE_F95, 1st_frame_no, No_of_frames,x0,y0,xdim,ydim
;    where x0,y0 represent the lower left corner 
;    and xdim,ydim the x,y-dimensions of the box to be restored.
; Note: program asks for Wiener filter parameters.
; In case of constant S/N, C=15 seems to be adequate.
; In case of variable S/N the approximation is:
; 	1/mtf(f)   for f le fo (in arcsec^-1)
;	mtf(f)/(mtf(f)^2+K*(f-fo))   for f gt fo
; i.e., increasing K we do more filtering.


;on_error,2

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

	stp=0.06218       ; SCALE (PIXEL SIZE) OF THE IMAGE  (arcsec/pix)
	lambda=5.250d-5   ; WORKING WAVELENGTH (cm)
	d=49.5		  ; DIAMETER OF THE TELESCOPE (cm)

; more parameters --------------------------------------------------------

st1='/scratch/sob/lap95/imc/s15gn'		; input name
st2='/scratch/sob/lap95/imc/c_c18Jun93.'	; output name
st3='.lap'

xs =966				; input frame size
ys =756

; ------------------------------------------------------------------------

n=xdim>ydim	; taking larger of the dimensions, so that
		; the basic working frame size will be n,n

a=findgen(n/2+1)/(n*stp)                   ; scale array for plotting
ima1=intarr(xs,ys)			   ; read-in array

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

mtf=dblarr(n)
cc=lambda*2.062648d5/(d*n*stp)   ; 206264.8 arcsec in radian
ss=0

   for i=1,n do begin
	      
      omega=cc*(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

window,2
plot,a,mtf,linestyle=5,xra=[-0.2,6],yra=[-0.2,5],/xst,/yst

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

newc:
print,'Wiener filter with constant S/N (=0) or variable (=1) ?'
if get_kbrd(1) eq '0' then begin
	read,'Filter parameter C =',c
	t1=(c+1)*mtf/(c*mtf*mtf+1)
endif else begin
	read,'Filter parameters K, fo : ',k,fo
	fp=nint(fo*n*stp)	; conversion to pixels
	t1=dblarr(n)		; filter
	uu=findgen(n-fp-1)+1.	; linearly increasing N/S
	uu=uu/n/stp
	t1(0:fp)=1./mtf(0:fp)
	t1(fp+1:n-1)=mtf(fp+1:n-1)/(mtf(fp+1:n-1)^2+k*uu)
endelse

wset,2
plot,a,t1,xra=[-0.2,6],yra=[-0.2,5],/xst,/yst,/noerase
print,'Change filter? (y/n)'
if get_kbrd(1) eq 'y' then goto,newc

print,'Do you want to stop here? (y/n)'
if get_kbrd(1) eq 'y' then begin
	out=0
	goto,fin
endif
print,'I am working...'

t1=float(t1)

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

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
print,'   2-D Wiener filter created'

; creating apodization window apw ----------------------------------------

   perc=2.		  ; percentage of frame size to be damped
   npx=nint(xdim*perc/100.)
   npy=nint(ydim*perc/100.)

   ux=fltarr(xdim)+1.		; x-bell
   ux(0)=0.
   for i=1,npx do ux(i)=(1.-cos(!pi*i/npx))/2.
   ux(xdim-npx:xdim-1)=reverse(ux(1:npx))

   uy=fltarr(ydim)+1.		; y-bell
   uy(0)=0.
   for i=1,npy do uy(i)=(1.-cos(!pi*i/npy))/2.
   uy(ydim-npy:ydim-1)=reverse(uy(1:npy))

   apw=ux#uy		; apodization window (xdim,ydim)
   print,'   Apodization window created'


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

print,nim,'  frames will be processed.'

for k=0,nim-1 do begin    ; loop over frames

   dcn=strtrim(k+init,2)
   if strlen(dcn) eq 3 then seq='0'+ dcn
   if strlen(dcn) eq 2 then seq='00'+dcn
   if strlen(dcn) eq 1 then seq='000'+dcn

   print,'Reading image:  ',st1+seq+st3

   openr,1,st1+seq+st3
   readu,1,ima1
   close,1
   ima=ima1(x0:x0+xdim-1,y0:y0+ydim-1)	; selected box

;window,3,xsize=xdim,ysize=ydim,title='ORIGINAL'
;tvscl,ima
;std=stdev(ima,mea)
;print,'Original image contrast =',std/mea

; apodization

   mea=mean(ima)
   ima=ima-mea
   ima=ima*apw
   ima=ima+mea	 	; apodized frame
  
; creating working array (n,n)

   im=intarr(n,n)+mea		; working array with mean(ima1)
   im(0,0)=ima		; insert image

; 2-D Fourier filtering
   
   f1=fft(im,-1)
   f1=float(fft(f1*t,1))   ; restored frame

   f1=nint(f1(0:xdim-1,0:ydim-1)) ; back to original size and type

;window,0,xsize=xdim,ysize=ydim,title='RESTORED'
;tvscl,f1
;std=stdev(f1,mea)
;print,'Restored image contrast =',std/mea


   print,'Saving image:  ',st2+dcn
   openw,1,st2+dcn
   writeu,1,f1
   close,1

endfor

print,'***** New image size:  ',xdim,ydim

fin:
print,'End of program.'
END