pro ota
a=[0.,0.,0.,0.,0.,0.5]
n=256
pm=2.
x1=indgen(n)#replicate(1.,n)
y1=replicate(1.,n)#indgen(n)
x1=(x1-n/2)*(4./n)/pm
y1=(y1-n/2)*(4./n)/pm
r=sqrt(x1^2+y1^2)
t=atan(y1,x1)
tx =1.8992573*r*cos(t)
ty =1.8992573*r*sin(t)
fo =3.8874443*(r^2-0.554450)
ax =2.3137662*r^2*cos(2*t)
ay =2.3137662*r^2*sin(2*t)     
cx =8.3345629*(r^3-0.673796*r)*cos(t)
cy =8.3345629*(r^3-0.673796*r)*sin(t)
clx=2.6701691*r^3*cos(3*t)
cly=2.6701691*r^3*sin(3*t)
sa =16.895979*(r^4-1.108900*r^2+0.241243)                  
i=complex(0.,1.)
wave=a(0)*fo+a(1)*cx+a(2)*cy+a(3)*ax+a(4)*ay+a(5)*sa
pupil=(r lt 1.) and (r gt 0.33) and (x1 lt 0.3)
;;         pupil(where((abs(x1) lt 0.011) or (abs(y1) lt 0.011)))=0.
;         pupil(where((shift(r,0.0,0.889*n*pm/4) lt 0.063) or $
;                     (shift(r,0.826*n*pm/4,-0.340*n*pm/4) lt 0.063) or $
;                     (shift(r,-0.762*n*pm/4,-0.467*n*pm/4) lt 0.063))) =0. 
;
; Orthonormalize
;
print,'ortyho'
fo=fo*pupil &fo=fo/sqrt(total(fo^2))
;
cx=cx*pupil &cx=cx-total(fo*cx)*fo & cx=cx/sqrt(total(cx^2))
;
cy=cy*pupil &cy=cy-total(fo*cy)*fo-total(cx*cy)*cx & cy=cy/sqrt(total(cy^2))
;
ax=ax*pupil &ax=ax-total(fo*ax)*fo-total(cx*ax)*cx-total(cy*ax)*cy 
ax=ax/sqrt(total(ax^2))
;
ay=ay*pupil 
ay=ay-total(fo*ay)*fo-total(cx*ay)*cx-total(cy*ay)*cy-total(ax*ay)*ax 
ay=ay/sqrt(total(ay^2))
;
sa=sa*pupil 
sa=sa-total(fo*sa)*fo-total(cx*sa)*cx-total(cy*sa)*cy-total(ax*sa)*ax-$
      total(ay*sa)*ay
sa=sa/sqrt(total(sa^2))
;
ap=[total(wave*fo),total(wave*cx),total(wave*cy),$
    total(wave*ax),total(wave*ay),total(wave*sa)]
print,'fo,cx,cy,ax,ay,sa',ap
print,total(wave*pupil-ap(0)*fo-ap(1)*cx-ap(2)*cy-ap(3)*ax-ap(4)*ay-ap(5)*sa)
;
;psfc=shift(abs(fft(pupil*exp(2*!pi*i*wav),1))^2,n/2,n/2)
;
stop
end