;PRO mtelew,a,mtf,wf1 PRO mtelew,a,mtf,wf1,wf2 ; Calculation of 1-D theoretical MTF of a diffraction-limited telescope ; and corresponding Wiener filter (with constant or variable SNR), ; with an optional PS plot. ; calling: MTELEW,a,mtf,wfi ;********************************************************************** STP=0.0405 ; STEP OF RESTORED IMAGES (arcsec/pix) lambda=6.020d-5 ; WORKING WAVELENGTH (cm) ; lambda=4.507d-5 n=1024 ; SIZE (LENGTH) OF ARRAY (pixels) d=97. ; DIAMETER OF THE TELESCOPE (cm) ;********************************************************************** a=findgen(n/2+1)/(n*stp) ; scale array for plotting 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^2) 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,xsiz=512,ysiz=640 plot,a,mtf,xra=[-0.2,12],yra=[-0.2,3.5],/xst,/yst ; creating the 1-D Wiener filter ---------------------------------------- ; constant SNR: read,'Filter parameter C = ',c wf1=(c+1)*mtf/(c*mtf*mtf+1) oplot,a,wf1,linestyle=5 ; oplot,a,wf1*mtf,linestyle=3 ; linear SNR: read,'Filter parameters K, fo : ',k,fo fp=nint(fo*n*stp) ; conversion to pixels wf2=dblarr(n) ; filter uu=findgen(n-fp-1)+1. ; linearly increasing N/S uu=uu/n/stp wf2(0:fp)=1./mtf(0:fp) wf2(fp+1:n-1)=mtf(fp+1:n-1)/(mtf(fp+1:n-1)^2+k*uu) oplot,a,wf2,linestyle=3 ; hardcopy: print,'Do you want a PS figure ? (y/n)' if get_kbrd(1) eq 'y' then begin set_plot,'ps' device,file='wiener3.ps',xsiz=16,ysiz=17,yoffs=9 plot,a,mtf,xra=[-0.2,12],yra=[-0.2,3.5],/xst,/yst, $ xtitle='Frequency (arcsec)!E-1!N', $ ytitle='MTF and Optimum Filter', $ thick=2,xthick=2,ythick=2,charthick=2,charsiz=1.2 oplot,a,wf1,linestyle=5,thick=1 oplot,a,wf2,linestyle=3,thick=2 plots,[-0.2,12],[0,0],linestyle=2 xyouts,0.7,2.7,'!3D = '+strtrim(fix(d),2)+' cm',charthick=2 xyouts,0.7,2.5,'!4k!3 = '+strtrim(fix(lambda*1.e7),2)+' nm',charthick=2 device,/close set_plot,'x' endif end