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

%    Reads already stored Zernike Decomps, 
%    and does graphics about them individually
%



ipart = 1 ;      %  1 - Load data saved by vire1.m, 
                 %               and make graphic


global R0 NMAX IBASIS IRECENT INDATPATH HICURVE ;
NMAX = 10 ;
IBASIS = 1 ;
IRECENT = 0 ;



nz = 9  ;     %  number of Zernikes to treat


iR0 = 2 ;        %  1 - R0 = 4.2
                 %  2 - R0 = 3.0


ips = 1 ;
ipe = 4 ;
          %  loop through all patients


if iR0 == 1 ;
  outsuf = 'a' ;
  R0 = 4.2 ;
elseif iR0 == 2 ;
  outsuf = 'b' ;  
  R0 = 3.0 ;
end ;


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



INDATPATH = '\matlab\cornean\indata\' ;
          %  path for input Keratron files
[COMAP, HICURVE] = coacmap ;







%  Read in full data matrix
%
fid = fopen(['\matlab\cornean\vire1' outsuf '.in'],'rt') ;
          %  'rt' is for "read only" and "text"

forstr = '%f %f %f %f %f' ;
          %  format for reading:
          %      file no.
          %      pat. no.
          %      day no.
          %      session no.
          %      eye no.
ninpar = 5 ;
ncoeff = 66 ;
for iformat = 1:ncoeff ;
  forstr = [forstr ' %f'] ;
end ;

[mdata,cnt] = fscanf(fid,forstr) ;
fclose(fid) ;

ncol = ninpar + ncoeff ;
nrow = cnt / ncol ;

mdata = reshape(mdata,ncol,nrow)' ;
          %   since reshape works column wise

%mdata(1:10,:)
%mdata(nrow,:)
          %  used for testing read in


vfilenum = mdata(:,1) ;
vpatnum = mdata(:,2) ;
vdaynum = mdata(:,3) ;
vsessnum = mdata(:,4) ;
veyenum = mdata(:,5) ;


vmedian = median(mdata(:,ninpar+1:ninpar+ncoeff)) ;


if ipart == 1 ;    %  then make main graphic


  for iz = 1:nz ;


    if iz == 1 ;
      titstr = 'Zernike, Order 0,0' ;
    elseif iz == 2 ;
      titstr = 'Zernike, Order 1,1, cos' ;
    elseif iz == 3 ;
      titstr = 'Zernike, Order 1,1, sin' ;
    elseif iz == 4 ;
      titstr = 'Zernike, Order 2,0' ;
    elseif iz == 5 ;
      titstr = 'Zernike, Order 2,2, cos' ;
    elseif iz == 6 ;
      titstr = 'Zernike, Order 2,2, sin' ;
    elseif iz == 7 ;
      titstr = 'Zernike, Order 3,1, cos' ;
    elseif iz == 8 ;
      titstr = 'Zernike, Order 3,1, sin' ;
    elseif iz == 9 ;
      titstr = 'Zernike, Order 4,0' ;
    end ;

    disp(['  Working on ' titstr]) ;


    icol = ninpar + iz ;
          %  index of column for current Zernikes

    vz = mdata(:,icol) ;    
          %  column vector of Zernike coefficients

    zmedian = median(vz) ;

    vzord = sort(vz) ;
    z02 = vzord(round(0.02 * nrow)) ;
    z25 = vzord(round(0.25 * nrow)) ; 
    z75 = vzord(round(0.75 * nrow)) ;
    z98 = vzord(round(0.98 * nrow)) ;



    minz = min(vz) ;
    maxz = max(vz) ;
    radz = max([(zmedian-minz); (maxz-zmedian)]) ;
    left = zmedian - 1.05 * radz ;
    right = zmedian + 1.05 * radz ;
    bottom = 0 ;
    top = 1 ;
    center = mean([top; bottom]) ;


    clf ;

    %  make main distribution plot
    %
    subplot(3,1,2) ;
      axis([left,right,bottom,top]) ;


    mkde1 = [] ;
    mkde2 = [] ;
    for ip = ips:ipe ;

      if ip == 1 ;    
        patname = 'Hickson' ;
        colstr = 'c' ;
      elseif ip == 2 ;    
        patname = 'Meadows' ;
        colstr = 'm' ;
      elseif ip == 3 ;    
        patname = 'Tesoro' ;
        colstr = 'y' ;
      elseif ip == 4 ;    
        patname = 'Yi' ;
        colstr = 'g' ;
      end ;

      
      vzpe1 = vz(find((vpatnum == ip) & (veyenum == 1))) ;
      vzpe2 = vz(find((vpatnum == ip) & (veyenum == 2))) ;

      pjump = (top - bottom) * 0.1 ;
      pbot = bottom + pjump * (ip - 1) ;
      ptop = top - pjump * (ip - 1) ;

      hold on ;
        plot([vzpe1'; vzpe1'], ...
                 [pbot; (pbot + pjump)]*ones(1,length(vzpe1)), ...
                                [colstr '-']) ;
        plot([vzpe2'; vzpe2'], ...
                 [ptop; (ptop - pjump)]*ones(1,length(vzpe2)), ...
                                [colstr '-']) ;
        plot([z02; z02],[bottom; top],'w-') ;
        plot([z25; z25],[bottom; top],'w-') ;
        plot([zmedian; zmedian],[bottom; top],'w-') ;
        plot([z75; z75],[bottom; top],'w-') ;
        plot([z98; z98],[bottom; top],'w-') ;

          tx = left + .9 * (right - left) ;
        text(tx,pbot + pjump,patname,'color',colstr) ;
        text(tx,ptop - pjump,patname,'color',colstr) ;

      hold off ;


      kde1 = gpkde(vzpe1,0,[left;right]) ;
      kde2 = gpkde(vzpe2,0,[left;right]) ;
      mkde1 = [mkde1, kde1] ;
      mkde2 = [mkde2, kde2] ;
          %  construct kde stuff for use in next lower plot

    end ;


    hold on ;
        tx = left + 0.05 * (right - left) ;
        ty1 = bottom + 0.25 * (top - bottom) ;
        ty2 = top - 0.25 * (top - bottom) ;
      text(tx,ty1,'Eye 1') ;
      text(tx,ty2,'Eye 2') ;
    hold off ;

    title(['Distributions for ' titstr date]) ;



    %  make plot with kde's, etc.
    %
    mkde1 = 0.48 * mkde1 / max(max(mkde1)) ;
    mkde2 = 0.48 * mkde2 / max(max(mkde2)) ;
          %  make max ht of eahc = 0.48 for plotting
    kdegrid = linspace(left,right,401)' ;

    subplot(3,1,3) ;
      axis([left,right,bottom,top]) ;
    
    for ip = ips:ipe ;

      if ip == 1 ;    
        patname = 'Hickson' ;
        colstr = 'c' ;
      elseif ip == 2 ;    
        patname = 'Meadows' ;
        colstr = 'm' ;
      elseif ip == 3 ;    
        patname = 'Tesoro' ;
        colstr = 'y' ;
      elseif ip == 4 ;    
        patname = 'Yi' ;
        colstr = 'g' ;
      end ;

      hold on ;

        plot(kdegrid,bottom + mkde1(:,ip),[colstr '-']) ;
        plot(kdegrid,top - mkde2(:,ip),[colstr '-']) ;

      hold off ;
      

    end ;
    


     


    %  show images expressing variability
    % 
    for iim = 1:5 ;

      vim = vmedian ;

      if iim == 1 ;
        imtitstr = '2% - tile' ;
        vim(iz) = vim(iz) + (z02 - zmedian) ;
      elseif iim == 2 ;
        imtitstr = '25% - tile' ;
        vim(iz) = vim(iz) + (z25 - zmedian) ;
      elseif iim == 3 ;
        imtitstr = '50% - tile' ;
      elseif iim == 4 ;
        imtitstr = '75% - tile' ;
        vim(iz) = vim(iz) + (z75 - zmedian) ;
      elseif iim == 5 ;
        imtitstr = '98% - tile' ;
        vim(iz) = vim(iz) + (z98 - zmedian) ;
      end ;

      %  Create image with Keratron Color scheme, and put in current axes
      %
      subplot(3,5,iim) ;
        betahat = vim' ;
        imagout = coazr(betahat) ;
        pcolor(xgrid,ygrid,imagout) ;
          shading flat ;
          axis('square') ;
          title(imtitstr) ;
            caxis([0 1]) ;
          colormap(COMAP) ;


    end ;



    orient landscape ;
  eval(['print -dpsc \matlab\cornean\vire5' num2str(iz) '.ps']) ;


  end ;

end ;


