disp('Running MATLAB script file PlotZbasis.m') ;
%

%    Makes plot of Zernike Basis Functions



maxsumind = 8 ;    %  maximum diagonal number to consider



%  set up preliminary stuff
%
global R0 NMAX IBASIS IRECENT INDATPATH HICURVE ;
R0 = 1 ;
NMAX = maxsumind ;
IBASIS = 1 ;
IRECENT = 0 ;
INDATPATH = '\matlab\cornean\indata\' ;
[COMAP, HICURVE] = coacmap ;


ngrid = 64 ;     %  Construct surfaces as ngrid x ngrid array
xgrid = linspace(-1, 1, ngrid) ;
ygrid = xgrid ;	        %  Vectors of image coordinate values 

[mxgrid, mygrid] = meshgrid(xgrid,ygrid) ;
              		%  Matrices of image coordinate values 
            
mrgrid = sqrt((mxgrid).^2 + (mygrid).^2) ;
              		%  Matrices of image polar coordinate radii
mthetag = atan2(mxgrid,mygrid) ;
              	%  Matrices of image polar coordinate angles (-pi,pi]
vind = find(mthetag < 0) ;
mthetag(vind) = mthetag(vind) + 2 * pi ;
              	%  Matrices of image polar coordinate angles [0,2pi)
          
mcrop = ones(ngrid,ngrid) ;            
  vindout = find(mrgrid > 1) ;
              	%  Vector of indices where image is outside circle
  vindin = find(mrgrid <= 1) ;
              	%  Vector of indices where image is inside circle
mcrop(vindout) = nan * ones(length(vindout),1) ;
mrgrid(vindout) = 0 * ones(length(vindout), 1) ;
              	%  Cropping Matrix, multiply by this to cut off images
              	%  outside circle of radius R0.    

vrgrid = mrgrid(:) ;
vthetag = mthetag(:) ;





%  get matrix of Zernike basis functions
%  evaluated at those points
%
mzer =  zernike(maxsumind, 1, vrgrid, vthetag) ;





%  get Zernike's for neutral colored sphere
%
sphererad = 7.7 ;
          %  Sphere radius 7.7 gives yellow sphere
spherehts = sphererad * sqrt(1 - (mrgrid(vindin) ./ sphererad).^2) ;
designm = zernike(NMAX,R0,mrgrid(vindin),mthetag(vindin)) ;
          %  construct Zernike design matrix 

spherefeatv = (designm' * designm) \ (designm' * spherehts) ;
          %  do ordinary least squares fit



%  Now loop through Zernike Coefficients
%
for sumind = 0:2:maxsumind ;


  for ifourier = 0:(sumind/2) ;


    ijacob = sumind - ifourier ;
          %  Index of jacobi (in columns)

    termstr = ['J ' num2str(ijacob) ...
               ', F ' num2str(ifourier)] ;

    %  Get Zernike term number
    %
    iz = (ijacob * (ijacob + 1)) / 2 ;
          %  number of Zernikes for all previous Jacobi orders
    if ifourier == 0 ;
      iz = iz + 1 ;
          %  Zernike term number,  since first in this column
    else ;
      izc = iz + ifourier ;
          %  Zernike cos term number
      izs = iz + ifourier + 1;
          %  Zernike sin term number


      iz = izc ;
          %  Just show cosine


      termstr = [termstr ', Angle = 0'] ;


    end ;




    termstr = ['Z ' num2str(iz) ', ' termstr] ;

    disp(['  working on ' termstr]) ;




    vzer = mzer(:,iz) ;
          %  column vector of current Zernike heights
    mvzer = reshape(vzer,ngrid,ngrid) ;

    mvzer = mcrop .* mvzer ;



    %  make an image plot of heights
    %
    subplot(2,2,1) ;
    pcolor(mxgrid,mygrid,mvzer) ;
      shading flat ;
      axis('square') ;
      title(['Height, ' termstr]) ;
      xlabel('x') ;
      ylabel('y') ;
      colormap(COMAP) ;




    %  make a mesh plot of heights
    %
    subplot(2,2,3) ;
    mesh(mxgrid,mygrid,mvzer) ;
      title(['Height, ' termstr]) ;
      xlabel('x') ;
      ylabel('y') ;




    %  make an image plot of curvature
    %
    subplot(2,2,2) ;

    featv = spherefeatv ;
              %  reset to yellow sphere
          cjump = 0.01 * 10^(-ijacob/2) ;
              %  Amount to increase this coefficient by
    featv(iz) = featv(iz) + cjump ;
    mout = coazr(featv) ;

    if iz == 1 ;     %  then are at first Zernike, so remember
                     %  this "constant" curvature
      avgcurv = mean(mout(vindin)) ;
    else ;           %  for all other terms, stretch outwards
                     %  from above average curvature
      sfact = avgcurv / max(abs(mout(vindin) - avgcurv)) ;
      mout = avgcurv + sfact * (mout - avgcurv) ;
    end ;



    pcolor(xgrid,ygrid,mout) ;
      title(['Curvature, ' termstr]) ;
      shading flat ;
      axis('square') ;
        caxis([0 1]) ;
      colormap(COMAP) ;

    
    eval(['print -dpsc \matlab\steve\ps\pzbas' ...
                  num2str(ijacob) num2str(ifourier) '.ps']) ;



%    disp('any key to continue') ;
%    pause ;


  end ;    %  of loop through Fourier orders


end ;    %  of loop through sumind


