disp('Running MATLAB script file vimu1.m') ;
%
%    For Analysis of Vistakon Multifocal Data
%    This part works with Test Images only
%



ipart = 12 ;      %  1 - Construct the lens, and do Zernike decomp., 
                 %              and make a raw curvature pic.
                 %  2 - Make a raw curvature pic and a Zernike recon.
                 %  3 - Overlay profiles of raw, Zernike recon., 
                 %               and radial recon.
%  Note:  these revealed basic construction was not so good, 
%         so next did:
                 %  11 - Construct the lens, and do Zernike decomp., 
                 %              and make a raw curvature pic.
                 %  12 - Make a raw curvature pic and a Zernike recon.
                 %  13 - Overlay profiles of raw, Zernike recon., 
                 %               and radial recon.



iprofile = 1 ;   %  1 - first test version



global INDATPATH NMAX R0 HICURVE IBASIS ;
INDATPATH = '\matlab\cornean\indata\' ;
          %  path for input Keratron files
[COMAP, HICURVE] = coacmap ;
ncomap = size(COMAP,1) ;
          %  number of rows in the color map ;


if ipart == 1 ;    %  construct test version

  figure(1) ;
  clf ;

  vrad = [1.0; 1.65; 2.125; 2.6; 4.0; 6.5; 7] ;
          %  vector of break points between regions

  if iprofile == 1 ;   
    profstr = 'First Simple Test' ;
    vcurv = [8.4; 8.8; 8.4; 8.8; 8.4; 8.6; 9] ;
          %  vector with radii of curvature for each region
  end ;



  %  Construct curve showing test surface profile
  %
  nsect = length(vrad) ;
  npps = 4 ;
          %  number of pieces per sector
  R0 = 4.2 ;
  NMAX = 10 ;
  ngrid = 64 ;     %  Construct surfaces as ngrid x ngrid array



  centx = 0 ;
  centy = 0 ;
  lastx = vcurv(1) ;
  lasty = 0 ;
  thetalast = 0 ;
  mcurve = [lastx, lasty] ;
          %  first piece of constructed curve
  vcur = vcurv(1) ;
  for isect = 1:nsect ;   %  loop through sectors

    r = vcurv(isect) ;
          %  current radius of curvature
    if isect == 1 ;
      arclength = vrad(1) ;
    else ;
      arclength = vrad(isect) - vrad(isect - 1) ;
    end ;
          %  current arglength
    vcur = [vcur; (r * ones(npps,1))] ;

    theta = arclength / r ;
          %  angle of current sector in radians

    firstx = lastx ;
    firsty = lasty ;
          %  make last corner point into current first corner point

    tgrid = linspace(thetalast,thetalast+theta,npps+1)' ;
          %  grid of thetas for this piece

    vx = centx + r * cos(tgrid) ;
    vy = centy + r * sin(tgrid) ;
          %  generate points on boundary of sector

    hold on ;
      plot([centx; vx; centx],[centy; vy; centy],'-') ;
      plot(centx,centy,'o') ;
    hold off ;
          %  add this sector to plot

    mcurve = [mcurve; [vx(2:npps+1),vy(2:npps+1)]] ;
          %  Build up current curve

    lastx = vx(npps+1) ;
    lasty = vy(npps+1) ;
    thetalast = thetalast + theta ;

    if isect < nsect ;
      centx = centx + (r - vcurv(isect + 1)) * cos(thetalast) ;
      centy = centy + (r - vcurv(isect + 1)) * sin(thetalast) ;
          %  next center point
    end ;

  end ;

  %  Save plot showing result
  %
  axis equal ;
  title(['Profile for ' profstr]) ;

    orient landscape ;
  print -dps '\matlab\cornean\vimu1a.ps' ;



  %  Next get heights and radii of profile, for Zernike decomp.
  %
  vht = mcurve(:,1) ;
          %  old x-coords are heights
  vht = vht - max(vht) ;
          %  record as distance downwards from apex normal
  vrad = mcurve(:,2) ;
          %  radii are old y-coordinates
  vcur = 337.5 ./ vcur ;

  %  make plot on this scale
  %
  figure(2) ;
  clf ;
  plot(vrad,vht) ;
    title(['Height Profile for ' profstr]) ;

    orient landscape ;
  print -dps '\matlab\cornean\vimu1b.ps' ;



  %  Now do Zernike decomposition, using lines from coazdec.m
  %
  nang = 256;    % number of semi-meridians read per radial ring
  nring = length(vrad) ;

  rdata = kron(vrad,ones(nang,1)) ;
  adata = linspace(0, (1 - 1/nang) * 2 * pi, nang)' ;
  adata = kron(ones(nring,1),adata) ;
  hdata = kron(vht,ones(nang,1)) ;
  cdata = kron(vcur,ones(nang,1)) ;
          %  extend to long column vectors

  %  now cutoff the values outside R0
  %
  vind = find(rdata > 1.05*R0) ;
           %  vector of indices where outside
  rdata(vind) = [] ;
  adata(vind) = [] ;
  hdata(vind) = [] ;
  cdata(vind) = [] ;
           %  Deletes that part of the vectors



  %  Construct raw curvature map, using lines from coarawd.m
  %
  [xdata, ydata] = pol2cart(adata, rdata);
  xgrid = linspace(-R0,R0,ngrid) ;
  ygrid = xgrid ;
            %  Vectors of image coordinate values 
  [mxgrid,mygrid,mout] = griddata(xdata,ydata,cdata,xgrid,ygrid') ;


  %  Create the matrix of image pixels
  %
  u01 = linspace(0,1,ncomap+1) ;
          %  integers representing the Keratron colors
  colcurv = interp1([-100; HICURVE], u01, mout(:)) ;
          %  use this to determine vector of colors representing
          %  curvature, when plotting.  This transforms curvature 
          %  to get Keratron colors, when using COMAP as colormap.
          %  NOTE: interp1 seems to need vector argument
  mout = reshape(colcurv,size(mout)) ;  

  %  Crop to image only inside circle
  %
  mcrop = ones(ngrid) ;            
    vind = find(mxgrid.^2 + mygrid.^2 > R0^2) ;
                %  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.       
  moutc = mout .* mcrop ;

  %  Create Image
  %
  figure(3) ;
  clf ;
  subplot(3,4,1) ;
  pcolor(xgrid,ygrid,moutc) ;
    shading flat ;
    axis('square') ;
      titstr = ['Raw C. Lens Only'] ;
    title(titstr) ;
      caxis([0 1]) ;
    colormap(COMAP) ;



    orient landscape ;
  print -dpsc '\matlab\cornean\vimu1c.ps' ;




  %  Do least square fitting, full basis
  %
  rdata = rdata / R0 ;
          %  normalize to Zernike radius 1
  designm = zernike(NMAX,R0,rdata,adata) ;
          %  construct Zernike design matrix 
  betahat = (designm' * designm) \ (designm' * hdata) ;
          %  do ordinary least squares fit
  IBASIS  = 1 ;
	  %  tells us that this is a normal decomposition

  %  Save as a .mat file
  %
  outfilestr = ['\matlab\cornean\vimu1' num2str(iprofile) 'a.mat'] ;  
  eval(['save ' outfilestr ' betahat IBASIS R0 NMAX']) ;



  %  Do least square fitting, radial basis
  %
  designm = radial(NMAX,rdata,adata) ;
          %  construct radial Zernike design matrix 
  betahat = (designm' * designm) \ (designm' * hdata) ;
          %  do ordinary least squares fit
  IBASIS  = 2 ;
	  %  saves a flag that tells us this is a radially symmetric
	  %  decomposition

  %  Save as a .mat file
  %
  outfilestr = ['\matlab\cornean\vimu1' num2str(iprofile) 'b.mat'] ;  
  eval(['save ' outfilestr ' betahat IBASIS R0 NMAX']) ;




elseif ipart == 2 ;


  figure(4) ;

  %  First load full basis file
  %
  infilestr = ['\matlab\cornean\vimu1' num2str(iprofile) 'a.mat'] ;  
  eval(['load ' infilestr]) ;

  
  %  Reconstruct, using lines from coazrec.m
  %
  %  Calculate grid quantities
  %
  ngrid = 64 ;     %  Construct surfaces as ngrid x ngrid array
  xgrid = linspace(-R0, R0, ngrid) ;
  ygrid = xgrid ;	        %  Vectors of image coordinate values 

  %  Create image with Keratron Color scheme
  %
  mout = coazr(betahat) ;
          %  note this uses IBASIS as just read in


  %  Generate map
  %
  subplot(3,4,1) ;
  pcolor(xgrid,ygrid,mout) ;
    shading flat ;
    axis('square') ;
      titstr = ['Full Zernike Rec.'] ;
    title(titstr) ;
    caxis([0 1]) ;
    colormap(COMAP) ;



  %  Next load radial basis file
  %
  infilestr = ['\matlab\cornean\vimu1' num2str(iprofile) 'b.mat'] ;  
  eval(['load ' infilestr]) ;


  
  %  Reconstruct, using lines from coarrec.m
  %
  %  Calculate grid quantities
  %
  ngrid = 64 ;     %  Construct surfaces as ngrid x ngrid array
  xgrid = linspace(-R0, R0, ngrid) ;
  ygrid = xgrid ;	        %  Vectors of image coordinate values 

  %  Create image with Keratron Color scheme
  %
  mout = coarr(betahat) ;
          %  note this uses IBASIS as just read in

  %  Generate map
  %
  subplot(3,4,6) ;
  pcolor(xgrid,ygrid,mout) ;
    shading flat ;
    axis('square') ;
      titstr = ['Radial Zernike Rec.'] ;
    title(titstr) ;
    caxis([0 1]) ;
    colormap(COMAP) ;




    orient landscape ;
  print -dpsc '\matlab\cornean\vimu1d.ps' ;



elseif ipart == 3 ;    %  then compare three types of profiles


  if isempty(mcurve) ;    %  then can't do this, give warning

    disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
    disp('!!!   Warning from vimu1.m, mcurve not found   !!!') ;
    disp('!!!   Need to run part 1 of this first         !!!') ;
    disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;

  else ;    %  then do comparisons


    figure(5) ;

    plot(vrad,vht) ;
      title('Comparison of Height Profiles') ;


    rgrid = linspace(0,R0,421) ;


    %  First load full basis file
    %
    infilestr = ['\matlab\cornean\vimu1' num2str(iprofile) 'a.mat'] ;  
    eval(['load ' infilestr]) ;

    hout = zrecon(betahat, NMAX, rgrid/R0, zeros(1,length(rgrid)), R0) ;
%[ht, dht, ddht] = zrecon(coefs, order, rgrid, angrid, rmax)
% ZRECON, returns a vector of fitted coefficients for a given
%	grid of radii and angles
%
% Inputs -
%	coefs 	- the Zernike coefficients from a particular fit
%	order	- maximal order of Zernike fitting
%	rgrid	- given input grid of radius values
%	angrid	- given grid of angular values
%	rmax	- maximum radius of interest

    hold on ;
      plot(rgrid,hout,'b') ;
    hold off ;



    %  Next load radial basis file
    %
    infilestr = ['\matlab\cornean\vimu1' num2str(iprofile) 'b.mat'] ;  
    eval(['load ' infilestr]) ;

    hout = rrecon(betahat, NMAX, rgrid/R0, zeros(1,length(rgrid)), R0) ;

    hold on ;
      plot(rgrid,hout,'r') ;

      plot(vrad,-2*ones(size(vrad)),'+') ;
          %  a plus where each observation is
        vax = axis ;
      plot([R0,R0],[vax(3),vax(4)],'g') ;
          %  vertical line at R0
    hold off ;

      orient landscape ;
    print -dpsc '\matlab\cornean\vimu1e.ps' ;




    figure(6) ;

    plot(vrad,vcur) ;
      title('Comparison of Curvature Profiles') ;


    rgrid = linspace(0,R0,421) ;


    %  First load full basis file
    %
    infilestr = ['\matlab\cornean\vimu1' num2str(iprofile) 'a.mat'] ;  
    eval(['load ' infilestr]) ;

    [mh,mhp,mhpp] = zrecon(betahat, NMAX, rgrid/R0, zeros(1,length(rgrid)), R0) ;
    mout = -mhpp ./ (1 + mhp.^2).^(3/2) ;
    mout = mout * 337.5 ;

    hold on ;
      plot(rgrid,mout,'b') ;
    hold off ;



    %  Next load radial basis file
    %
    infilestr = ['\matlab\cornean\vimu1' num2str(iprofile) 'b.mat'] ;  
    eval(['load ' infilestr]) ;

    [mh,mhp,mhpp] = rrecon(betahat, NMAX, rgrid/R0, zeros(1,length(rgrid)), R0) ;
    mout = -mhpp ./ (1 + mhp.^2).^(3/2) ;
    mout = mout * 337.5 ;

    hold on ;
      plot(rgrid,mout,'r') ;

      plot(vrad,38*ones(size(vrad)),'+') ;
          %  a plus where each observation is
        vax = axis ;
      plot([R0,R0],[vax(3),vax(4)],'g') ;
          %  vertical line at R0
    hold off ;

      orient landscape ;
    print -dpsc '\matlab\cornean\vimu1f.ps' ;



    axis([vax(1),vax(2),36,44]) ;

      orient landscape ;
    print -dpsc '\matlab\cornean\vimu1g.ps' ;



  end ;



elseif ipart == 11 ;    %  construct test version

  figure(1) ;
  clf ;

  vrad = [1.0; 1.65; 2.125; 2.6; 4.0; 6.5; 7] ;
          %  vector of break points between regions

  if iprofile == 1 ;   
    profstr = 'First Simple Test' ;
    vcurv = [8.4; 8.8; 8.4; 8.8; 8.4; 8.6; 9] ;
          %  vector with radii of curvature for each region
  end ;



  %  Construct curve showing test surface profile
  %
  nsect = length(vrad) ;
  npps = 4 ;
          %  number of pieces per sector
  R0 = 4.2 ;
  NMAX = 10 ;
  ngrid = 64 ;     %  Construct surfaces as ngrid x ngrid array



  centx = 0 ;
  centy = 0 ;
  lastx = vcurv(1) ;
  lasty = 0 ;
  thetalast = 0 ;
  mcurve = [lastx, lasty] ;
          %  first piece of constructed curve
  vcur = vcurv(1) ;
  for isect = 1:nsect ;   %  loop through sectors

    r = vcurv(isect) ;
          %  current radius of curvature
    if isect == 1 ;
      arclength = vrad(1) ;
    else ;
      arclength = vrad(isect) - vrad(isect - 1) ;
    end ;
          %  current arglength
    vcur = [vcur; (r * ones(npps,1))] ;

    theta = arclength / r ;
          %  angle of current sector in radians

    firstx = lastx ;
    firsty = lasty ;
          %  make last corner point into current first corner point

    tgrid = linspace(thetalast,thetalast+theta,npps+1)' ;
          %  grid of thetas for this piece

    vx = centx + r * cos(tgrid) ;
    vy = centy + r * sin(tgrid) ;
          %  generate points on boundary of sector

    %  check for being outside of R0
    %
    if vy(npps+1) <= R0 ;    %  then are still inside, do as above

      hold on ;
        plot([centx; vx; centx],[centy; vy; centy],'-') ;
        plot(centx,centy,'o') ;
      hold off ;
          %  add this sector to plot

      mcurve = [mcurve; [vx(2:npps+1),vy(2:npps+1)]] ;
          %  Build up current curve

      lastx = vx(npps+1) ;
      lasty = vy(npps+1) ;
      thetalast = thetalast + theta ;

      if isect < nsect ;
        centx = centx + (r - vcurv(isect + 1)) * cos(thetalast) ;
        centy = centy + (r - vcurv(isect + 1)) * sin(thetalast) ;
          %  next center point
      end ;

    else ;    %  then this grid goes outside of R0,
              %  need to go to smaller local grid, and end loop


      thetastop = asin((R0 - centy) / r) ;
          %  angle of current sector in radians

      tgrid = linspace(thetalast,thetastop,npps+1)' ;
          %  grid of thetas for this piece

      vx = centx + r * cos(tgrid) ;
      vy = centy + r * sin(tgrid) ;
          %  generate points on boundary of sector


      hold on ;
        plot([centx; vx; centx],[centy; vy; centy],'-') ;
        plot(centx,centy,'o') ;
      hold off ;
          %  add this sector to plot

      mcurve = [mcurve; [vx(2:npps+1),vy(2:npps+1)]] ;
          %  Build up current curve


      break ;
    
    end ;


  end ;

  %  Save plot showing result
  %
  axis equal ;
  title(['Improved Profile for ' profstr]) ;


    orient landscape ;
  print -dps '\matlab\cornean\vimu1h.ps' ;






  %  Next get heights and radii of profile, for Zernike decomp.
  %
  vht = mcurve(:,1) ;
          %  old x-coords are heights
  vht = vht - max(vht) ;
          %  record as distance downwards from apex normal
  vrad = mcurve(:,2) ;
          %  radii are old y-coordinates
  vcur = 337.5 ./ vcur ;


  %  Include sphere, fit to apex normal and last point
  %
  x1 = mcurve(1,1) ;
  x2 = mcurve(size(mcurve,1),1) ;
  y2 = mcurve(size(mcurve,1),2) ;
  dif = x1 - x2 ;
  srad = (dif^2 + y2^2) / (2 * dif) ;
          %  radius of the sphere which goes through endpoints


  stheta = asin(y2 / srad) ;
          %  angle of sector of sphere
  sagrid = linspace(pi/2 - stheta,pi/2,101) ;
  sx = srad * cos(sagrid) ;
  sy = srad * sin(sagrid) - srad ;
          %  sector of sphere for plotting

  %  make plot on this scale
  %
  figure(2) ;
  clf ;
  plot(vrad,vht,'b',sx,sy,'r--') ;
    title(['Improved Height Profile for ' profstr]) ;
    legend('Profile',['Sphere, r = ' num2str(srad)]) ;


    orient landscape ;
  print -dpsc '\matlab\cornean\vimu1i.ps' ;



  %  Now do Zernike decomposition, using lines from coazdec.m
  %
  nang = 256;    % number of semi-meridians read per radial ring
  nring = length(vrad) ;

  rdata = kron(vrad,ones(nang,1)) ;
  adata = linspace(0, (1 - 1/nang) * 2 * pi, nang)' ;
  adata = kron(ones(nring,1),adata) ;
  hdata = kron(vht,ones(nang,1)) ;
  cdata = kron(vcur,ones(nang,1)) ;
          %  extend to long column vectors


  dsx = sx(1:5:length(sx))' ;
  dsy = sy(1:5:length(sy))' ;
          %  decimate to every 5th entry
  srdata = kron(dsx,ones(nang,1)) ;
  sadata = linspace(0, (1 - 1/nang) * 2 * pi, nang)' ;
  sadata = kron(ones(length(dsx),1),sadata) ;
  shdata = kron(dsy,ones(nang,1)) ;
          %  corresponding stuff for sphere

  %  now cutoff the values outside R0
  %
  vind = find(rdata > 1.05*R0) ;
           %  vector of indices where outside
  rdata(vind) = [] ;
  adata(vind) = [] ;
  hdata(vind) = [] ;
  cdata(vind) = [] ;
           %  Deletes that part of the vectors



  %  Construct raw curvature map, using lines from coarawd.m
  %
  [xdata, ydata] = pol2cart(adata, rdata);
  xgrid = linspace(-R0,R0,ngrid) ;
  ygrid = xgrid ;
            %  Vectors of image coordinate values 
  [mxgrid,mygrid,mout] = griddata(xdata,ydata,cdata,xgrid,ygrid') ;


  %  Create the matrix of image pixels
  %
  u01 = linspace(0,1,ncomap+1) ;
          %  integers representing the Keratron colors
  colcurv = interp1([-100; HICURVE], u01, mout(:)) ;
          %  use this to determine vector of colors representing
          %  curvature, when plotting.  This transforms curvature 
          %  to get Keratron colors, when using COMAP as colormap.
          %  NOTE: interp1 seems to need vector argument
  mout = reshape(colcurv,size(mout)) ;  

  %  Crop to image only inside circle
  %
  mcrop = ones(ngrid) ;            
    vind = find(mxgrid.^2 + mygrid.^2 > R0^2) ;
                %  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.       
  moutc = mout .* mcrop ;

  %  Create Image
  %
  figure(3) ;
  clf ;
  subplot(3,4,1) ;
  pcolor(xgrid,ygrid,moutc) ;
    shading flat ;
    axis('square') ;
      titstr = ['Improved Raw C. Lens Only'] ;
    title(titstr) ;
      caxis([0 1]) ;
    colormap(COMAP) ;



    orient landscape ;
  print -dpsc '\matlab\cornean\vimu1j.ps' ;




  %  Do least square fitting, full basis
  %
  rdata = rdata / R0 ;
          %  normalize to Zernike radius 1
  designm = zernike(NMAX,R0,rdata,adata) ;
          %  construct Zernike design matrix 
  betahat = (designm' * designm) \ (designm' * hdata) ;
          %  do ordinary least squares fit
  IBASIS  = 1 ;
	  %  tells us that this is a normal decomposition

  %  also do zernike fit for sphere


  %  Save as a .mat file
  %
  outfilestr = ['\matlab\cornean\vimu1' num2str(iprofile) 'c.mat'] ;  
  eval(['save ' outfilestr ' betahat IBASIS R0 NMAX']) ;



  %  Do least square fitting, radial basis
  %
  designm = radial(NMAX,rdata,adata) ;
          %  construct radial Zernike design matrix 
  betahat = (designm' * designm) \ (designm' * hdata) ;
          %  do ordinary least squares fit
  IBASIS  = 2 ;
	  %  saves a flag that tells us this is a radially symmetric
	  %  decomposition


  %  Do same for sphere data
  %
  srdata = srdata / R0 ;
  designm = radial(NMAX,srdata,sadata) ;
          %  construct radial Zernike design matrix 
  sbetahat = (designm' * designm) \ (designm' * shdata) ;
          %  do ordinary least squares fit
  


  %  Save as a .mat file
  %
  outfilestr = ['\matlab\cornean\vimu1' num2str(iprofile) 'd.mat'] ;  
  eval(['save ' outfilestr ' betahat sbetahat IBASIS R0 NMAX']) ;




elseif ipart == 12 ;


  figure(4) ;

  %  First load full basis file
  %
  infilestr = ['\matlab\cornean\vimu1' num2str(iprofile) 'c.mat'] ;  
  eval(['load ' infilestr]) ;

  
  %  Reconstruct, using lines from coazrec.m
  %
  %  Calculate grid quantities
  %
  ngrid = 64 ;     %  Construct surfaces as ngrid x ngrid array
  xgrid = linspace(-R0, R0, ngrid) ;
  ygrid = xgrid ;	        %  Vectors of image coordinate values 

  %  Create image with Keratron Color scheme
  %
  mout = coazr(betahat) ;
          %  note this uses IBASIS as just read in


  %  Generate map
  %
  subplot(3,4,1) ;
  pcolor(xgrid,ygrid,mout) ;
    shading flat ;
    axis('square') ;
      titstr = ['Improved Full Zernike Rec.'] ;
    title(titstr) ;
    caxis([0 1]) ;
    colormap(COMAP) ;



  %  Next load radial basis file
  %
  infilestr = ['\matlab\cornean\vimu1' num2str(iprofile) 'd.mat'] ;  
  eval(['load ' infilestr]) ;


  
  %  Reconstruct, using lines from coarrec.m
  %
  %  Calculate grid quantities
  %
  ngrid = 64 ;     %  Construct surfaces as ngrid x ngrid array
  xgrid = linspace(-R0, R0, ngrid) ;
  ygrid = xgrid ;	        %  Vectors of image coordinate values 

  %  Create image with Keratron Color scheme
  %
  mout = coarr(betahat) ;
          %  note this uses IBASIS as just read in

  %  Generate map
  %
  subplot(3,4,6) ;
  pcolor(xgrid,ygrid,mout) ;
    shading flat ;
    axis('square') ;
      titstr = ['Improved Radial Zernike Rec.'] ;
    title(titstr) ;
    caxis([0 1]) ;
    colormap(COMAP) ;



  %  Reconstruct Sphere, using similar stuff
  %
  %  Calculate grid quantities
  %
  ngrid = 64 ;     %  Construct surfaces as ngrid x ngrid array
  xgrid = linspace(-R0, R0, ngrid) ;
  ygrid = xgrid ;	        %  Vectors of image coordinate values 

  %  Create image with Keratron Color scheme
  %
  mout = coarr(sbetahat) ;
          %  note this uses IBASIS as just read in

  %  Generate map
  %
  subplot(3,4,11) ;
  pcolor(xgrid,ygrid,mout) ;
    shading flat ;
    axis('square') ;
      titstr = ['Spherical Fit'] ;
    title(titstr) ;
    caxis([0 1]) ;
    colormap(COMAP) ;




    orient landscape ;
  print -dpsc '\matlab\cornean\vimu1k.ps' ;



elseif ipart == 13 ;    %  then compare three types of profiles


  if isempty(mcurve) ;    %  then can't do this, give warning

    disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
    disp('!!!   Warning from vimu1.m, mcurve not found   !!!') ;
    disp('!!!   Need to run part 1 of this first         !!!') ;
    disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;

  else ;    %  then do comparisons


    figure(5) ;

    plot(vrad,vht) ;
      title('Comparison of Improved Height Profiles') ;


    rgrid = linspace(0,R0,421) ;


    %  First load full basis file
    %
    infilestr = ['\matlab\cornean\vimu1' num2str(iprofile) 'c.mat'] ;  
    eval(['load ' infilestr]) ;

    hout = zrecon(betahat, NMAX, rgrid/R0, zeros(1,length(rgrid)), R0) ;
%[ht, dht, ddht] = zrecon(coefs, order, rgrid, angrid, rmax)
% ZRECON, returns a vector of fitted coefficients for a given
%	grid of radii and angles
%
% Inputs -
%	coefs 	- the Zernike coefficients from a particular fit
%	order	- maximal order of Zernike fitting
%	rgrid	- given input grid of radius values
%	angrid	- given grid of angular values
%	rmax	- maximum radius of interest

    hold on ;
      plot(rgrid,hout,'b') ;
    hold off ;



    %  Next load radial basis file
    %
    infilestr = ['\matlab\cornean\vimu1' num2str(iprofile) 'd.mat'] ;  
    eval(['load ' infilestr]) ;

    hout = rrecon(betahat, NMAX, rgrid/R0, zeros(1,length(rgrid)), R0) ;

    hold on ;
      plot(rgrid,hout,'r') ;

      plot(vrad,-.6*ones(size(vrad)),'+') ;
          %  a plus where each observation is
        vax = axis ;
      plot([R0,R0],[vax(3),vax(4)],'g') ;
          %  vertical line at R0
    hold off ;

      orient landscape ;
    print -dpsc '\matlab\cornean\vimu1l.ps' ;




    figure(6) ;

    plot(vrad,vcur) ;
      title('Comparison of Improved Curvature Profiles') ;


    rgrid = linspace(0,R0,421) ;


    %  First load full basis file
    %
    infilestr = ['\matlab\cornean\vimu1' num2str(iprofile) 'c.mat'] ;  
    eval(['load ' infilestr]) ;

    [mh,mhp,mhpp] = zrecon(betahat, NMAX, rgrid/R0, zeros(1,length(rgrid)), R0) ;
    mout = -mhpp ./ (1 + mhp.^2).^(3/2) ;
    mout = mout * 337.5 ;

    hold on ;
      plot(rgrid,mout,'b') ;
    hold off ;



    %  Next load radial basis file
    %
    infilestr = ['\matlab\cornean\vimu1' num2str(iprofile) 'd.mat'] ;  
    eval(['load ' infilestr]) ;

    [mh,mhp,mhpp] = rrecon(betahat, NMAX, rgrid/R0, zeros(1,length(rgrid)), R0) ;
    mout = -mhpp ./ (1 + mhp.^2).^(3/2) ;
    mout = mout * 337.5 ;

    hold on ;
      plot(rgrid,mout,'r') ;

      plot(vrad,38*ones(size(vrad)),'+') ;
          %  a plus where each observation is
        vax = axis ;
      plot([R0,R0],[vax(3),vax(4)],'g') ;
          %  vertical line at R0
    hold off ;

      orient landscape ;
    print -dpsc '\matlab\cornean\vimu1m.ps' ;



    axis([vax(1),vax(2),36,44]) ;

      orient landscape ;
    print -dpsc '\matlab\cornean\vimu1n.ps' ;



  end ;





end ;


