PRO wrenorm,date,ncarr
;
;RMS granular contrast calculation of normalized WL images
;restored for the instrumental MTF.
;INPUT:  date = 23 or 25
;	 files XXjun/wnXX.seq, where XX is the date
;OUTPUT: ncarr = 2-column array with mean and rms contrast
;
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/wn'+dat+'.'
im=intarr(1528,1024)

if date eq 23 then begin
	nima=320
	x1=924
	x2=1179
	y1=564
	y2=819 
endif else begin
	nima=359
	x1=231
	x2=486
	y1=1
	y2=256
endelse

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

	stp=0.0829	  ; SCALE (PIXEL SIZE) OF THE IMAGE  (arcsec/pix)
	lambda=8.0d-5     ; WORKING WAVELENGTH (cm)
	d=44.		  ; DIAMETER OF THE TELESCOPE (cm)
	n=256		  ; SIZE OF THE RESTORED FIELD
	cc=25.		  ; WIENER FILTER PARAMETER C

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

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))

; 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, restoration, and statistics -------------------------

ncarr=fltarr(2,nima+1)

for i=0,nima do begin

  openr,1,nami+strtrim(i,2)
  readu,1,im
  close,1
  ima=im(x1:x2,y1:y2)

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

; apodization

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

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

; MEAN AND RMS CONTRAST --------------------------------------------------

  std=stdev(f1(14:240,14:240),me2)
  ncarr(0,i)=me2
  ncarr(1,i)=std/me2
  print,i,me2,ncarr(1,i)

endfor

END