pro arte2,psf,mtf,spsf,smtf,wfilt

; a routine to produce arrays for plotting etc.
; psf: input psf (saved in p...), mtf: unsmoothed mtf (output)
; spsf: inverse transform of smtf (output)
; smtf: smoothed mtf by splines (output)
; wfilt: resulting filter (output)
; all arrays are one-dimensional, 512 (513) long.
   
on_error,2

mtf=float(fft(psf,-1))*1024
mtf=mtf(0:512)
 
a=findgen(513)/(0.124*1024)	; arcsec^-1 scale for plotting

no1=[0,5,10,15,20,30,40,50,75,100,150,200]  ; 1st part of nodes
no2=[250,300,400,512]
nel=4
print,'default set of nodes:'
print,[no1,no2]
goto,exesmoo	
 
newsmoo:
print,'old 2nd part of nodes:',no2
read,'no. of elements in the new 2nd part = ',nel
no2=intarr(nel)
read,'new 2nd part of nodes: ',no2

exesmoo:
tension=1.                 ; control of rigidity of splines
nodes=[no1,no2]
smtf=smspline(mtf,tension,12+nel,nodes)

window,0
plot,mtf,psym=3,xra=[-10,520],/xst
oplot,smtf

newdet:
print,'do you want a detail? (y/n)'
if get_kbrd(1) eq 'n' then goto,endsmoo
read,'beginning and end positions: ',i1,i2
i2=i2<512
plot,mtf,psym=1,xra=[i1-5,i2+5],yra=[min(mtf(i1:i2)),max(mtf(i1:i2))],/xst,/yst
oplot,smtf
oplot,intarr(600)
goto,newdet

endsmoo:
print,'smoothing ok? (y/n)
if get_kbrd(1) eq 'n' then goto,newsmoo

smtf=smtf/smtf(0)      ; renormalization

; definition of the cutoff frequency...
mi=min(smtf,pmi)             ; looking for minimum
print,'minimum value of mtf = ',mi,' at position ',pmi

if mi gt 0. then begin
   print,'The minimum is positive. Shall we define it as the cutoff? (y/n)'
   if get_kbrd(1) eq 'n' then stop
   smtf(pmi+1:512)=0.
   cros0=pmi+1
endif else begin
   cros0=min(where(smtf le 0))
   smtf(cros0:512)=0.
endelse

cr0=1/(cros0/(0.124*1024))
print,'cutoff position at',cros0,' pixels, i.e.',cr0,' arcsec'

; smoothed psf generation

spsf=float(fft([smtf,reverse(smtf(1:511))],1))
spsf=spsf(0:512)/1024.
plot,psf(0:20),psym=1
oplot,spsf(0:20)
print,'press any key to continue'
pau=get_kbrd(1)

; --------------- creation of the smoothed Moon limb ----------------

psfs=shift(psf,256)

mol=fltarr(512)   ;"moonlimb"
for i=0,511 do mol(i)=total(psfs(0:i))   ; integration of psf = limb
mol=[mol,reverse(mol)]

fmol=fft(mol,-1)

heavis=fltarr(512)			; Heaviside step function
heavis(256:511)=1.

; ---------------------- filtering ----------------------------------

window,4,xsize=550,ysize=700
plot,a,smtf,linestyle=5,xra=[-0.2,4.2],yra=[-0.2,6],/xst,/yst

newc:
wdelete,0,2,3

read,'Wiener filter parameter C= ',c
wfilt=(c+1)*smtf/(c*smtf^2+1)

wset,4
plot,a,wfilt,xra=[-0.2,4.2],yra=[-0.2,6],/xst,/yst,/noerase

fmols=fmol*[wfilt,reverse(wfilt(1:511))]
rmol=float(fft(fmols,1))

window,6,xsize=600,ysize=700
plot,heavis,linestyle=1,xra=[206,306],yra=[-0.1,1.1],/xst,/yst
oplot,mol,linestyle=1
oplot,rmol

print,'change filter? (y/n)'
if get_kbrd(1) eq 'y' then goto,newc

end