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

%    Reads already stored Zernike Decomps, 
%    and does graphics about them individually
%    Improves vire5.m, by showing amplitude and phase
%    for sin-cos terms



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


iimplot = 1 ;    %  1 to add images to Fourier 0 plots
                 %  0 to skip



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


maxsumind = 6 ;    %  maximum of the sum of the indices,
                   %  4 for 9 zernike terms
                   %  6 for 16 zernike terms



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 sumind = 0:2:maxsumind ;


    %  First work with Fourier order 0, and make vire5 type plot 
    %
    ifourier = 0 ;

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


    %  Get Zernike term number
    %
    iz = (ijacob * (ijacob + 1)) / 2 ;
          %  number of Zernikes for all previous Jacobi orders
    iz = iz + 1 ;
          %  Zernike term number,  since first in this column
 


    titstr = ['Zernike, Order ' num2str(ijacob) ',' ...
                     num2str(ifourier)] ;


    disp(['  Working on ' titstr]) ;


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

    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]) ;
      set(gca,'Box','On') ;


    mkde1 = [] ;
    mkde2 = [] ;
    vzpe105 = [] ; 
    vzpe205 = [] ; 
    vzpe1med = [] ; 
    vzpe2med = [] ; 
    vzpe195 = [] ; 
    vzpe295 = [] ; 
    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 + 0.5 * pjump,patname,'color',colstr) ;
        text(tx,ptop - 0.5 * 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


      nzpe1 = length(vzpe1) ;
      nzpe2 = length(vzpe2) ;
      vzpe1o = sort(vzpe1) ;
      vzpe2o = sort(vzpe2) ;
      zpe105 = vzpe1o(round(0.05 * nzpe1)) ;
      zpe205 = vzpe2o(round(0.05 * nzpe2)) ;
      zpe1med = median(vzpe1o) ;
      zpe2med = median(vzpe2o) ;
      zpe195 = vzpe1o(round(0.95 * nzpe1)) ;
      zpe295 = vzpe2o(round(0.95 * nzpe2)) ;

      vzpe105 = [vzpe105 ; zpe105] ; 
      vzpe205 = [vzpe205 ; zpe205] ; 
      vzpe1med = [vzpe1med ; zpe1med] ; 
      vzpe2med = [vzpe2med ; zpe2med] ; 
      vzpe195 = [vzpe195 ; zpe195] ; 
      vzpe295 = [vzpe295 ; zpe295] ; 


    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]) ;
      set(gca,'Box','On') ;
    
    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 ;

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


      hold on ;

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

        plot(vzpe1med(ip),pbot,[colstr 'o']) ;
        plot(vzpe2med(ip),ptop,[colstr 'o']) ;
          vachil = get(gca,'Children') ;
          set(vachil(1),'LineWidth',2) ;
          set(vachil(2),'LineWidth',2) ;
          set(vachil(1),'MarkerSize',8) ;
          set(vachil(2),'MarkerSize',8) ;

        plot([vzpe105(ip); vzpe195(ip)], ...
                              [pbot; pbot],[colstr '-']) ;
        plot([vzpe205(ip); vzpe295(ip)], ...
                              [ptop; ptop],[colstr '-']) ;
          vachil = get(gca,'Children') ;
          set(vachil(1),'LineWidth',3) ;
          set(vachil(2),'LineWidth',3) ;

      hold off ;
      

    end ;
    


    if iimplot == 1 ;    %  then 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 ;

    end ;


      orient landscape ;
    eval(['print -dpsc \matlab\cornean\vire6' num2str(ijacob) ...
                                   num2str(ifourier) '.ps']) ;





    %  Next make plots for other Fourier terms
    %
    for ifourier = 1:(sumind/2) ;

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


      %  Get Zernike term number
      %
      iz = (ijacob * (ijacob + 1)) / 2 ;
          %  number of Zernikes for all previous Jacobi orders
      izc = iz + ifourier ;
          %  Zernike cos term number
      izs = iz + ifourier + 1;
          %  Zernike sin term number



      titstr = ['Zernike, Order ' num2str(ijacob) ',' ...
                       num2str(ifourier)] ;


      disp(['  Working on ' titstr]) ;


      icolc = ninpar + izc ;
          %  index of column for current Zernike, cos
      icols = ninpar + izs ;
          %  index of column for current Zernike, sin

      vzc = mdata(:,icolc) ;    
          %  column vector of Zernike coefficients, cos
      vzs = mdata(:,icols) ;    
          %  column vector of Zernike coefficients, sin


      vzr = sqrt(vzc.^2 + vzs.^2) ;
      vzt = (atan2(vzs,vzc) * 180 / pi) + 180 ;


      zmedian = median(vzr) ;


      vzord = sort(vzr) ;
      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(vzr) ;
      maxz = max(vzr) ;
      left = - 0.05 * maxz ;
      right = 1.05 * maxz ;
      bottom = 0 ;
      top = 1 ;
      center = mean([top; bottom]) ;


      clf ;

      %  make main distribution plot, for radii
      %
      subplot(4,1,1) ;
      axis([left,right,bottom,top]) ;
      set(gca,'Box','On') ;


      mkde1 = [] ;
      mkde2 = [] ;
      vzpe105 = [] ; 
      vzpe205 = [] ; 
      vzpe1med = [] ; 
      vzpe2med = [] ; 
      vzpe195 = [] ; 
      vzpe295 = [] ; 
      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 = vzr(find((vpatnum == ip) & (veyenum == 1))) ;
        vzpe2 = vzr(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 + 0.5 * pjump,patname,'color',colstr) ;
          text(tx,ptop - 0.5 * 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


        nzpe1 = length(vzpe1) ;
        nzpe2 = length(vzpe2) ;
        vzpe1o = sort(vzpe1) ;
        vzpe2o = sort(vzpe2) ;
        zpe105 = vzpe1o(round(0.05 * nzpe1)) ;
        zpe205 = vzpe2o(round(0.05 * nzpe2)) ;
        zpe1med = median(vzpe1o) ;
        zpe2med = median(vzpe2o) ;
        zpe195 = vzpe1o(round(0.95 * nzpe1)) ;
        zpe295 = vzpe2o(round(0.95 * nzpe2)) ;

        vzpe105 = [vzpe105 ; zpe105] ; 
        vzpe205 = [vzpe205 ; zpe205] ; 
        vzpe1med = [vzpe1med ; zpe1med] ; 
        vzpe2med = [vzpe2med ; zpe2med] ; 
        vzpe195 = [vzpe195 ; zpe195] ; 
        vzpe295 = [vzpe295 ; zpe295] ; 


      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 ' radii, ' 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 each = 0.48 for plotting
      kdegrid = linspace(left,right,401)' ;

      subplot(4,1,2) ;
      axis([left,right,bottom,top]) ;
      set(gca,'Box','On') ;

      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 ;

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


        hold on ;

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

          plot(vzpe1med(ip),pbot,[colstr 'o']) ;
          plot(vzpe2med(ip),ptop,[colstr 'o']) ;
            vachil = get(gca,'Children') ;
            set(vachil(1),'LineWidth',2) ;
            set(vachil(2),'LineWidth',2) ;
            set(vachil(1),'MarkerSize',8) ;
            set(vachil(2),'MarkerSize',8) ;

          plot([vzpe105(ip); vzpe195(ip)], ...
                                [pbot; pbot],[colstr '-']) ;
          plot([vzpe205(ip); vzpe295(ip)], ...
                                [ptop; ptop],[colstr '-']) ;
            vachil = get(gca,'Children') ;
            set(vachil(1),'LineWidth',3) ;
            set(vachil(2),'LineWidth',3) ;

        hold off ;
      

      end ;
    




      %  make main distribution plot, for angles
      %
      subplot(4,1,3) ;
      set(gca,'Box','On') ;

        left = -180 ;
        right = 540 ;
        kdeleft = -360 ;
        kderight = 720 ;
      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 = vzt(find((vpatnum == ip) & (veyenum == 1))) ;
        vzpe2 = vzt(find((vpatnum == ip) & (veyenum == 2))) ;

        vzpe1 = [(vzpe1 - 360); vzpe1; (vzpe1 + 360)] ;
        vzpe2 = [(vzpe2 - 360); vzpe2; (vzpe2 + 360)] ;

        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 '-']) ;

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

        hold off ;


        kde1 = gpkde(vzpe1,0,[kdeleft;kderight]) ;
        kde2 = gpkde(vzpe2,0,[kdeleft;kderight]) ;
        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') ;

        plot([0; 0],[bottom; top],'w-') ;
        plot([360; 360],[bottom; top],'w-') ;

      hold off ;

      title(['Distributions for ' titstr ' angles, ' 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 each = 0.48 for plotting
      kdegrid = linspace(kdeleft,kderight,401)' ;

      subplot(4,1,4) ;
      axis([left,right,bottom,top]) ;
      set(gca,'Box','On') ;

      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 ;

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


        hold on ;

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

          plot([0; 0],[bottom; top],'w-') ;
          plot([360; 360],[bottom; top],'w-') ;

        hold off ;
      

      end ;
    




        orient tall ;
      eval(['print -dpsc \matlab\cornean\vire6' num2str(ijacob) ...
                                   num2str(ifourier) '.ps']) ;



    end ;






  end ;








end ;


