! Change the Zernike coefficients to OSA(rhw) set, to 7th order
! then plot them out
! vec2 holds the coefficients, indexed from 0 up, 
! with 0 being piston, so ignored
!		RHWebb, Spring 2000.  modified, winter 2004

maxterms = 52	
	! need 52 terms from standard to get 36 for OSA to order 7
wave = 1
field = 1
sampling = 1 
	! 32x32
vector = 1  
	!  where the coefficients are put  (vec1(j+8))
zerntype = 1  
	! standard
GETZERNIKE maxterms,wave,field,sampling,vector,zerntype	
	! Got Zernike Standard coefficients, put them in vec1

print "Using OSA coeff's,"
FORMAT 16.6
PRINT "Peak to Valley : ", vec1(1)
PRINT "RMS to zero    : ", vec1(2)
PRINT "RMS to mean    : ", vec1(3)
PRINT "RMS to centroid: ", vec1(4)
PRINT "Variance       : ", vec1(5)
PRINT "Strehl ratio   : ", vec1(6)
PRINT "RMS Fit Error  : ", vec1(7)
PRINT "Max Fit Error  : ", vec1(8)

	! Now shift them in index and put in vec2, in OSA order
vec2(0) = vec1(8+1)
vec4(0) = "piston"
vec2(1) = vec1(8+3)	
vec4(1) = "xtilt  2rsin0"
vec2(2) = vec1(8+2)	
Vec4(2) = "ytilt  2rcos0"
vec2(3) = vec1(8+5)	
vec4(3) = "H astig sqrt(6)*r^2sin20"
vec2(4) = vec1(8+4)	
vec4(4) = "defocus sqrt(3)*(2r^2 -1)"
vec2(5) = vec1(8+6)	
vec4(5) = "45 astig sqrt(6)*r^2cos20"

vec2(6) = vec1(8+9)	
        !r^3sin30
vec2(7) = vec1(8+7)	
       !(3r^3-2r)sin30
vec2(8) = vec1(8+8)	
       !(3r^3-2r)cos30
vec2(9) = vec1(8+10)	
       !r^3cos30

vec2(10) = vec1(8+15)	
       !r^4sin40
vec2(11) = vec1(8+13)	
       !(4r^4-3r^2)sin20
vec2(12) = vec1(8+11)	
       !6r^4-6r^2+1
vec2(13) = vec1(8+12)	
       !(4r^4-3r^2)cos20
vec2(14) = vec1(8+14)	
       !r^4cos40

vec2(15) = vec1(8+21)	
       !r^5sin50
vec2(16) = vec1(8+19)	
       !(5r^5-4r^3)sin30
vec2(17) = vec1(8+17)		
       !(10r^5-12r^3+3r)sin0
vec2(18) = vec1(8+16)		
       !(10r^5-12r^3+3r)cos0
vec2(19) = vec1(8+18)	
       !(5r^5-4r^3)cos30
vec2(20) = vec1(8+20)	
       !r^5cos50

vec2(21) = vec1(8+27)		
       !r^6sin60
vec2(22) = vec1(8+25)	
       !(6r^6-5r^4)sin40
vec2(23) = vec1(8+23)	
       !(15r^6-20r^4+6r^2)sin20
vec2(24) = vec1(8+22)	
       !20r^6-30r^4+12r^2-1
vec2(25) = vec1(8+24)	
       !(15r^6-20r^4+6r^2)cos20
vec2(26) = vec1(8+26)	
       !(6r^6-5r^4)cos40
vec2(27) = vec1(8+28)		
       !r^6cos60

vec2(28) = vec1(8+35)		
       !r^7sin70
vec2(29) = vec1(8+33)		
       !r^7...sin50
vec2(30) = vec1(8+31)	
       !r^7...sin30
vec2(31) = vec1(8+47)	
       !r^7...sin0
vec2(32) = vec1(8+46)	
       !r^7...cos0
vec2(33) = vec1(8+48)	
       !r^7...cos30
vec2(34) = vec1(8+50)
       !r^7...cos50
vec2(35) = vec1(8+52)		
       !r^7cos70

cmax = 0
cmin = 0

k = 0
	n = inte(0.5*(sqrt(8*k+1)-1))
	m = 2*k-n*(n+2)
	FORMAT 2.0
	PRINT "C(",k,"[",n,",",m,"]) = ",
		if (vec2(k) > cmax) then cmax = vec2(k)
		if (vec2(k) < cmin) then cmin = vec2(k)
	FORMAT 10.6
	PRINT vec2(k),
	print "   piston"
k=1 
	n = inte(0.5*(sqrt(8*k+1)-1))
	m = 2*k-n*(n+2)
	FORMAT 2.0
	PRINT "C(",k,"[",n,",",m,"]) = ",
		if (vec2(k) > cmax) then cmax = vec2(k)
		if (vec2(k) < cmin) then cmin = vec2(k)
	FORMAT 10.6
	PRINT vec2(k),
	print "   ytilt  2rsin0"
k=2
	n = inte(0.5*(sqrt(8*k+1)-1))
	m = 2*k-n*(n+2)
	FORMAT 2.0
	PRINT "C(",k,"[",n,",",m,"]) = ",
		if (vec2(k) > cmax) then cmax = vec2(k)
		if (vec2(k) < cmin) then cmin = vec2(k)
	FORMAT 10.6
	PRINT vec2(k),
	print "   xtilt  2rcos0"
k=3

	n = inte(0.5*(sqrt(8*k+1)-1))
	m = 2*k-n*(n+2)
	FORMAT 2.0
	PRINT "C(",k,"[",n,",",m,"]) = ",
		if (vec2(k) > cmax) then cmax = vec2(k)
		if (vec2(k) < cmin) then cmin = vec2(k)
	FORMAT 10.6
	PRINT vec2(k),
	print "   oblique astig sqrt(6)*r^2sin20"
k=4 
	n = inte(0.5*(sqrt(8*k+1)-1))
	m = 2*k-n*(n+2)
	FORMAT 2.0
	PRINT "C(",k,"[",n,",",m,"]) = ",
		if (vec2(k) > cmax) then cmax = vec2(k)
		if (vec2(k) < cmin) then cmin = vec2(k)
	FORMAT 10.6
	PRINT vec2(k),
	print "   defocus sqrt(3)*(2r^2 -1)"
k=5
	n = inte(0.5*(sqrt(8*k+1)-1))
	m = 2*k-n*(n+2)
	FORMAT 2.0
	PRINT "C(",k,"[",n,",",m,"]) = ",
		if (vec2(k) > cmax) then cmax = vec2(k)
		if (vec2(k) < cmin) then cmin = vec2(k)
	FORMAT 10.6
	PRINT vec2(k),
	print "   orthog astig sqrt(6)*r^2cos20"
k = 6
	n = inte(0.5*(sqrt(8*k+1)-1))
	m = 2*k-n*(n+2)
	FORMAT 2.0
	PRINT "C(",k,"[",n,",",m,"]) = ",
		if (vec2(k) > cmax) then cmax = vec2(k)
		if (vec2(k) < cmin) then cmin = vec2(k)
	FORMAT 10.6
	PRINT vec2(k),
	print "   trefoil1 sqrt(8)*sin30"
k=7 
	n = inte(0.5*(sqrt(8*k+1)-1))
	m = 2*k-n*(n+2)
	FORMAT 2.0
	PRINT "C(",k,"[",n,",",m,"]) = ",
		if (vec2(k) > cmax) then cmax = vec2(k)
		if (vec2(k) < cmin) then cmin = vec2(k)
	FORMAT 10.6
	PRINT vec2(k),
	print "   coma1 sqrt(8)* (3r^3 -2r)sin0"
k=8
	n = inte(0.5*(sqrt(8*k+1)-1))
	m = 2*k-n*(n+2)
	FORMAT 2.0
	PRINT "C(",k,"[",n,",",m,"]) = ",
		if (vec2(k) > cmax) then cmax = vec2(k)
		if (vec2(k) < cmin) then cmin = vec2(k)
	FORMAT 10.6
	PRINT vec2(k),
	print "   coma2 sqrt(8)* (3r^3 -2r)cos0"
k=9

	n = inte(0.5*(sqrt(8*k+1)-1))
	m = 2*k-n*(n+2)
	FORMAT 2.0
	PRINT "C(",k,"[",n,",",m,"]) = ",
		if (vec2(k) > cmax) then cmax = vec2(k)
		if (vec2(k) < cmin) then cmin = vec2(k)
	FORMAT 10.6
	PRINT vec2(k),
	print "   trefoil2 sqrt(8)*cos30"
k=10 
	n = inte(0.5*(sqrt(8*k+1)-1))
	m = 2*k-n*(n+2)
	FORMAT 2.0
	PRINT "C(",k,"[",n,",",m,"]) = ",
		if (vec2(k) > cmax) then cmax = vec2(k)
		if (vec2(k) < cmin) then cmin = vec2(k)
	FORMAT 10.6
	PRINT vec2(k),
	print "   sqrt(10)*r^4sin40"
k=11
	n = inte(0.5*(sqrt(8*k+1)-1))
	m = 2*k-n*(n+2)
	FORMAT 2.0
	PRINT "C(",k,"[",n,",",m,"]) = ",
		if (vec2(k) > cmax) then cmax = vec2(k)
		if (vec2(k) < cmin) then cmin = vec2(k)
	FORMAT 10.6
	PRINT vec2(k),
	print "   sqrt(10)*(4r^4-3r^2)sin20"
k=12

	n = inte(0.5*(sqrt(8*k+1)-1))
	m = 2*k-n*(n+2)
	FORMAT 2.0
	PRINT "C(",k,"[",n,",",m,"]) = ",
		if (vec2(k) > cmax) then cmax = vec2(k)
		if (vec2(k) < cmin) then cmin = vec2(k)
	FORMAT 10.6
	PRINT vec2(k),
	print "   spherical  sqrt(5)*(6r^4-6r^2+1)"
k=13 
	n = inte(0.5*(sqrt(8*k+1)-1))
	m = 2*k-n*(n+2)
	FORMAT 2.0
	PRINT "C(",k,"[",n,",",m,"]) = ",
		if (vec2(k) > cmax) then cmax = vec2(k)
		if (vec2(k) < cmin) then cmin = vec2(k)
	FORMAT 10.6
	PRINT vec2(k),
	print "   sqrt(10)*(4r^4-3r^2)cos20"
k=14
	n = inte(0.5*(sqrt(8*k+1)-1))
	m = 2*k-n*(n+2)
	FORMAT 2.0
	PRINT "C(",k,"[",n,",",m,"]) = ",
		if (vec2(k) > cmax) then cmax = vec2(k)
		if (vec2(k) < cmin) then cmin = vec2(k)
	FORMAT 10.6
	PRINT vec2(k),
	print "   sqrt(10)*r^4cos40"


for k=15,35,1
	! OSA uses the scaling of the Standard coefficients
	n = inte(0.5*(sqrt(8*k+1)-1))
	m = 2*k-n*(n+2)
	FORMAT 2.0
	PRINT "C(",k,"[",n,",",m,"]) = ",
		if (vec2(k) > cmax) then cmax = vec2(k)
		if (vec2(k) < cmin) then cmin = vec2(k)
	FORMAT 10.6
	PRINT vec2(k)
!	PRINT vec4(k)
next
crange = cmax-cmin
cscale = cmax
if (cmax < -cmin) then cscale = -cmin
       ! cmax is the biggest positive coefficient
       ! cmin is the biggest negative coefficient

FORMAT 16.6
LC = logT(cscale)

IC = inte(LC)
C = expT(IC)
CC = cscale/C
IC = inte(CC)
C = C*IC
Ctick = C/cscale
       ! now scale for graphing
for k=0,35,1
	FORMAT 2.0
	vec3(k) = vec2(k)/cscale
next

!	************************************************
graphics screen
xmx = xmax()
xmn = xmin()
ymx = ymax()
ymn = ymin()
xwidth = xmx-xmn
ywidth = ymx-ymn
xleft = xmn + (.05*xwidth)
xrigh = xmn + (.95*xwidth)
ytopp = ymn + (.05*ywidth)
ybott = ymn + (.75*ywidth)
!  so frame goes ybott to ytopp,  which is .6*ywidth
line xleft, ytopp-.01*ywidth, xrigh, ytopp-.01*ywidth
line xrigh, ytopp-.01*ywidth, xrigh, ybott+.01*ywidth
line xrigh, ybott+.01*ywidth, xleft, ybott+.01*ywidth
line xleft, ybott+.01*ywidth, xleft, ytopp-.01*ywidth

!gtitle "The Rain in Spain falls mainly on the Plain"
glensname
gdate

ymid = ymn +  (.35*ywidth)
! half way up is ymn + .3*ywidth: position of x axis
dx = xmx/50
line  xleft+dx, ymid, xleft+36*dx, ymid

!Yticks
!  this is the scale line, with label
yt = 0.3*ywidth*Ctick
FORMAT .3
line xleft,ymid-yt,xleft+0.7*xwidth,ymid-yt
gtext xleft+0.72*xwidth,ymid-yt,0,C

! Xticks
FORMAT 2.0
dy = 0.01*ywidth
for k = 1,35,1
	line xleft+k*dx+dx/2,ymid-dy,xleft+k*dx+dx/2,ymid+dy	
! extra long at k = 4
	if k==4
		line xleft+k*dx+dx/2,ymid-2*dy,xleft+k*dx+dx/2,ymid+2*dy			
		if vec3(k)>0 
			gtext xleft+dx/2+k*dx,ybott-2*dy,0,k
		else 
			gtext xleft+dx/2+k*dx,ytopp+2*dy,0,k
		endif
	endif
! extra long at k = 12
	if k==12
		line xleft+k*dx+dx/2,ymid-2*dy,xleft+k*dx+dx/2,ymid+2*dy			
		if vec3(k)>0 
			gtext xleft+dx/2+k*dx,ybott-2*dy,0,k
		else 
			gtext xleft+dx/2+k*dx,ytopp+2*dy,0,k
		endif
	endif
! extra long at k = 24
	if k==24
		line xleft+k*dx+dx/2,ymid-2*dy,xleft+k*dx+dx/2,ymid+2*dy			
		if vec3(k)>0 
			gtext xleft+dx/2+k*dx,ybott-2*dy,0,k
		else 
			gtext xleft+dx/2+k*dx,ytopp+2*dy,0,k
		endif
	endif

! plot
!  * for unused
!	if (vec3(k) == 999)
!		line xleft+k*dx, ymid, xleft+k*dx, ymid
!		line xleft+k*dx, ymid, xleft+k*dx+dx, ymid
	!	line xleft+k*dx+dx, ymid, xleft+k*dx+dx, ymid	
	!	gtext xleft+dx/2+k*dx,ymid,0,"*"	
!	else
! this is the coefficient bar
		v = 0.3*ywidth*vec3(k)
		line xleft+k*dx, ymid, xleft+k*dx, ymid - v
		line xleft+k*dx, ymid-v, xleft+k*dx+dx, ymid - v
		line xleft+k*dx+dx, ymid-v, xleft+k*dx+dx, ymid
		
		dd = dx/10
		for j = 1,9,1
			line xleft+k*dx+j*dd, ymid, xleft+k*dx+j*dd, ymid - v
		next
!	endif

next
!
graphics off