%MATLAB Script File: cobas1.m
%
%   For the Cornea Study, this show Surface Plots of Basis Functions
%
clf reset ;
format compact ;
running = 'MATLAB Script File:  cobas1.m'

iplot = 11 ;   %  1 - Simple mesh plot
              %  2 - Mesh plot with contour plot beneath
              %  3 - Contour plot
              %  4 - 3-d Contour plot
              %  5 - Surface plot
              %  6 - Surface plot with contour plot beneath
              %  7 - Pseudocolor plot - faceted
              %  8 - Pseudocolor plot - flat
              %  9 - Pseudocolor plot - interpolated
              %  10 - flat and interpolated on same plot
              %  11- Surface plot - light source
              %  12 - Image plot
              %            (this is not working yet, since seems)
              %            (similar to pcolor, see 2-77 in User's Guide)

ibas = 1 ;    %  1 - Radial Fourier - Legendre
              %  2 - Diameter Fourier - Legendre
              %            (not working yet)
              %  3 - Zernike
              %            (not working yet)


r0 = 4.2 ;    %  Usual radius of analysis

ngrid = 61 ; %  Construct surfaces as ngrid x ngrid array

xgrid = linspace(-r0,r0,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) ;
  vind = find(mthetag < 0) ;
mthetag(vind) = mthetag(vind) + 2 * pi ;
              %  Matrices of image polar coordinate angles [0,2pi)
          
mcrop = ones(ngrid,ngrid) ;            
  vind = find(mrgrid > r0) ;
              %  Vector of indices where image is outside circle
mcrop(vind) = nan * ones(length(vind),1) ;
              %  Cropping Matrix, multiply by this to cut off images
              %  outside circle of radius r_0.       


if ibas == 1 | ibas == 2 ;     %  then use a Fourier - Legendre Basis

  l = 1 ;     %  Legendre Order
  m = 1 ;     %  Fourier Order


  if ibas == 1 ; %  then radial Diameter Fourier - Legendre Basis

  %  Calculate Legendre Part
  polyarg = 2 * (mrgrid ./ r0) - 1 ;
  if l == 0 ;
    mleg = ones(nxgrid) ;
  elseif l == 1 ;
    mleg = polyarg ;
  elseif l == 2 ;
    mleg = (1/2) * (3 * polyarg^2 - 1) ;
  elseif l == 3 ;
    mleg = (1/2) * (5 * polyarg^3 - 3 * polyarg) ;
  elseif l == 4 ;
    mleg = (1/8) * (35 * polyarg^4 - 30 * polyarg^2 + 3) ;
  elseif l == 5 ;
    mleg = (1/8) * (63 * polyarg^5 - 70 * polyarg^3 + 15 * polyarg) ;
  end ;

  mpoly = cos(m * mthetag) .* mleg ;
            %  Combine Legendre Part With Fourier Part
  mpoly = mcrop .* mpoly ;
            %  Cut off part out side circle

  elseif ibas == 2 ; %  then use Diameter Fourier - Legendre Basis
  end ;

elseif ibas == 3 ; %  then use Zernike Basis

end ;

if iplot == 1 ;      %  Then do simple mesh plot
  mesh(xgrid,ygrid,mpoly) ;

elseif iplot == 2 ;  %  Then do mesh with countour beneath
  meshc(xgrid,ygrid,mpoly) ;

elseif iplot == 3 ;  %  Then do contour plot
  contour(xgrid,ygrid,mpoly,20) ;

elseif iplot == 4 ;  %  Then do 3-d contour plot
  contour3(mxgrid,mygrid,mpoly,20) ;

elseif iplot == 5 ;  %  Then do surface plot
  surf(mxgrid,mygrid,mpoly) ;

elseif iplot == 6 ;  %  Then do surface plot with contour beneath
  surfc(mxgrid,mygrid,mpoly) ;

elseif iplot == 7 ;  %  Then do pseudocolor plot - faceted
  pcolor(mxgrid,mygrid,mpoly) ;
  shading faceted ;
  axis('square') ;

elseif iplot == 8 ;  %  Then do pseudocolor plot - flat
  pcolor(mxgrid,mygrid,mpoly) ;
  shading flat ;
  axis('square') ;

elseif iplot == 9 ;  %  Then do pseudocolor plot - interpolated
  pcolor(mxgrid,mygrid,mpoly) ;
  shading interp ;
  axis('square') ;

elseif iplot == 10 ;  %  Then do pseudocolor plot - flat & interpolated
  subplot(1,2,1) ;
    pcolor(mxgrid,mygrid,mpoly) ;
    shading flat ;
    axis('square') ;
            %  Sets all axes to be the same, i.e. correct aspect ratio
  subplot(1,2,2) ;
    pcolor(mxgrid,mygrid,mpoly) ;
    shading interp ;
    axis('square') ;

elseif iplot == 11 ;  %  Then do surface plot
  surfl(xgrid,ygrid,mpoly) ;
  shading flat ;
  colormap(gray) ;
elseif iplot == 12 ;  %  Then do "image" plot
  image(xgrid,ygrid,mpoly) ;

end ;

