#=============================================================================
#
# Marc A. Murison
# Astronomical Applications Dept.
# U.S. Naval Observatory
# 3450 Massachusetts Ave., NW
# Washignton, DC 20392
# murison@riemann.usno.navy.mil
# http://riemann.usno.navy.mil/AESOP/
#
#=============================================================================

if ZSERIES_READ <> 1 then
ZSERIES_READ := 1:
read`d:/optics/AESOP/zerncoef.p`:

ZernikeSeries := proc( OPD::algebraic, N::integer, order_list::list )

    local   n, m, nn, i, k, cA, cB, cAB, nmstr:
    global  S, c:

    c := 'c':
    S := 'S':
    S := 0:
    i := 1:
    for n from 0 to N do
        for m from 0 to n do
            if type(n-m,even) = true then

				if n=0 and m=0 then
					debug_print( procname, `(disk averaged OPD)`, 0 );
				elif n=1 and m=1 then
					debug_print( procname, `(wavefront tilt)`, 0 );
				elif n=2 and m=0 then
					debug_print( procname, `(defocus)`, 0 );
				elif n=2 and m=2 then
					debug_print( procname, `(astigmatism)`, 0 );
				elif n=3 and m=1 then
					debug_print( procname, `(coma)`, 0 );
				elif n=3 and m=3 then
					debug_print( procname, `(trefoil)`, 0 );
				elif n=4 and m=0 then
					debug_print( procname, `(spherical)`, 0 );
				fi;

				nmstr := convert([n,m],string):
				debug_print( procname, `calculating the coefficients `.nmstr.`...`, 0 );
                cAB := ZSeriesCoeff(n,m,OPD):

                cA  := collect( cAB['A'], order_list, factor ):
				if cA <> 0 then
					debug_print( procname, `A`.nmstr.`:`, 1, cA );
				fi;

                cB  := collect( cAB['B'], order_list, factor ):
				if cB <> 0 then
					debug_print( procname, `B`.nmstr.`:`, 1, cB );
				fi;

                if cA <> 0 or cB <> 0 then
                    S := S + cA*ZernikeCos(n,m,rho,phi)
                           + cB*ZernikeSin(n,m,rho,phi):
                    c[i]['index'] := [n,m]:
                    c[i]['A']     := cA:
                    c[i]['B']     := cB:
                    i := i + 1:
                fi:

            fi:
        od:
    od:

    gc():
    RETURN( NULL ):

end:



#=============================================================================
# Given Zernike coeffs c from ZernikeSeries output, extract the [n,m] term.
#=============================================================================
ZernikeTerm := proc( n::integer, m::integer, c::table )
	local k, p;
	p := 0;
	for k from 1 to nops(op(op(c))) do
		if c[k]['index']=[n,m] then
			p := c[k]['A']*ZernikeCos(n,m,rho,phi)
				 + c[k]['B']*ZernikeSin(n,m,rho,phi):
			break:
		fi:
	od:
	RETURN(eval(p));
end:


#=============================================================================
# user help
#=============================================================================
`help/text/ZernikeSeries` := TEXT(
``,
` ZernikeSeries( OPD:algebraic, N:integer, order_list:list, prnflag:integer )`,
``,
` Express a wavefront OPD as a series in Zernike polynomials.  All nonzero`,
` series coefficients, up to order n=N, are determined.  OPD must be `,
` in cylindrical coordinates (rho,phi), where rho has been normalized `,
` to lie in the range [0..1].  Print the coefficients if prnflag=ON.`,
``,
` global variables:`,
`    c      the nonzero Zernike series coefficients`,
`    S      Zernike series, sorted in increasing order of [n,m]`,
``,
` Return: NULL`,
``,
` c is indexed according to `,
`    c[i][index] = [n,m]`,
`    c[i][A]     = A[n,m]`,
`    c[i][B]     = B[n,m]`,
` where i runs from 1 to the total number of nonzero coefficient pairs.`,
``,
` See also: zernike, ZSeriesCoeff`
):

fi;