;S curve  -   P. Bely     June 1990

;set_plot,3
answer=' '
lamref=0.5E-6

read, 'Aperture ? (1=full, 0.666 for 2/3)', aperture

;read, 'Coma in waves (RMS) ? (at 500nm)', coma

D=2.4 
 
range=0.4         ;angle range in arcseconds
rad2sec=180.*3600./!pi
i=complex(0,1)
n1=80                 ;number of points for the wavefront angle range
n2=15                 ; number of points along the mirror radius
n3=9                  ; number of angles around 1/4 the mirror 
theta=(float(indgen(n1)-n1/2)/n1)*range/rad2sec  ; wavefront angle

r=float(indgen(n2))/n2                     ; running mirror radius

alpha=!pi/(4*n3)+float(indgen(n3))*!pi/(2*n3)
esall=fltarr(n1,6)
for coma=0.0,0.5,0.05 do begin
;z_coma=8.3345629*(r^3-0.673796*r)*lamref*coma
z_coma=2.314*(r^2)*lamref*coma
;plot,z_coma/lamref

suma=fltarr(10,n1)
sumb=fltarr(10,n1)
lambda=0.466e-6+indgen(10)*(0.7e-6-0.466e-6)/9.
k=2*!pi/lambda
for ii=0,n1-1 do begin
if((ii/10)*10 eq ii) then print,ii,n1,coma
  sum_a=fltarr(10) & sum_b=sum_a 
  for j=0.33*(n2-1),(n2-1)*aperture  do begin
     for jj=0,n3-1 do begin
       z1=theta(ii)*r(j)*sin(alpha(jj))*d/2.
       zc= z_coma(j)*sin(2*alpha(jj))
       sum_A=sum_A+$
             (abs(exp(i*k*(z1-zc))+exp(i*k*(-z1+zc)-i*!pi/2)))^2*r(j)
       sum_B=sum_B+$
             (abs(exp(i*k*(z1-zc))-exp(i*k*(-z1+zc)-i*!pi/2)))^2*r(j)
     endfor
  endfor
  suma(0,ii)=sum_A
  sumb(0,ii)=sum_b
endfor
es=fltarr(n_elements(suma(0,*)))
for iy=0,n_elements(suma(0,*))-1 do $
        es(iy)=total(suma(*,iy)-sumb(*,iy))/(total(suma(*,iy)+sumb(*,iy)))
esall(0,coma/0.05)=es
; Plotting
;---------

title= 'Aperture' + string(aperture) 
!xtitle= 'Milliarcseconds'
!ytitle='Error signal'
!mtitle=title + 'with/wo astigmatism'
!linetype=coma/0.05
if(coma eq 0.0) then plot,theta*rad2sec*1000.,ES
if(coma ne 0.0) then oplot,theta*rad2sec*1000.,ES
endfor
set_plot,'ps'
for coma=0.,0.25,0.05 do begin
   !linetype=coma/0.05
   if(coma eq 0.0) then device,filename='fgs.ps'
   if(coma eq 0.0) then plot,theta*rad2sec*1000.,ESall(*,coma/0.05)
   if(coma ne 0.0) then oplot,theta*rad2sec*1000.,ESall(*,coma/0.05)
endfor
set_plot,'X'
stop
end