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

%    Reads already stored Zernike Decomps, 
%    and does graphics about them individually
%    Modifies vire7.m, by an improved graphic
%    that shows more interesting structure



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


ip = 0 ;      %  Patient number
              %  0 - loop through all
              %  1 - Hickson
              %  2 - Meadows
              %  3 - Tesoro
              %  4 - Yi
 
ieye = 0 ;    %   eyenumber,   1  or  2
              %   0 - loop through all


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


tol = 1 ;     %  multiple of standard deviation, to use as
              %  "good fit" tolerance



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



if ip == 0 ;        %  loop through all patients
  ips = 1 ;
  ipe = 4 ;
else ;
  ips = ip ;
  ipe = ip ;
end ;



if ieye == 0 ;        %  loop through all eyes
  ieyes = 1 ;
  ieyee = 2 ;
else ;
  ieyes = ieye ;
  ieyee = ieye ;
end ;



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





%  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) ;






if ipart == 1 ;    %  then make main graphic


  for iip = ips:ipe ;    % loop through patients

    if iip == 1 ;    
      patname = 'Hickson' ;
    elseif iip == 2 ;    
      patname = 'Meadows' ;
    elseif iip == 3 ;    
      patname = 'Tesoro' ;
    elseif iip == 4 ;    
      patname = 'Yi' ;
    end ;


    for iieye = ieyes:ieyee ;    %  loop through eyes



      pemdata = mdata(find(vpatnum == iip & ...
                                     veyenum == iieye),:) ;
          %  data matrix reduced for this patient and eye

      nrowpe = size(pemdata,2) ;

      pevdaynum = pemdata(:,3) ;
      pevsessnum = pemdata(:,4) ;

      pevmedian = median(pemdata(:,ninpar+1:ninpar+ncoeff)) ;



      %  Now loop through Zernike Coefficients
      %
      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 = [patname ', Eye ' num2str(iieye) ...
                      ', Zernike, Order ' num2str(ijacob) ...
                      ',' num2str(ifourier)] ;



        disp(['  Working on ' titstr]) ;


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

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

        zmean = mean(vz) ;
        zstd = std(vz) ;

        tolmin = zmean-tol*zstd ;
        tolmax = zmean+tol*zstd ;

        minz = min(vz) ;
        maxz = max(vz) ;
        radz = max([(zmean-minz); (maxz-zmean)]) ;
        left = zmean - 1.05 * radz ;
        right = zmean + 1.05 * radz ;


        %  make main distribution plot
        %
        [mainkde,mainxgrid] = gpkde(vz,0,[left,right]) ;
        bottom = 0 ;
        top = 1.05 * max(mainkde) ;



        %  Make plot for days
        %
        clf ;
        subplot(2,1,1) ;

        plot(mainxgrid,mainkde,'w-') ;
          set(gca,'Box','On') ;
          title([titstr ', ' date]) ;
          axis([left,right,bottom,top]) ;
          hold on ;
            plot([zmean,zmean],[bottom,top],'-w') ;
            plot([tolmin,tolmin],[bottom,top],'-w') ;
            plot([tolmax,tolmax],[bottom,top],'-w') ;
          hold off ;


        %  Loop through and add stuff for days
        %
        for day = 1:3 ;

          if day == 1 ;    
            colstr = 'c' ;
          elseif day == 2 ;    
            colstr = 'm' ;
          elseif day == 3 ;    
            colstr = 'g' ;
          end ;

          vzped = vz(find(pevdaynum == day)) ;

          pwi = 100 * sum((vzped > tolmin) & ...
                              (vzped < tolmax)) / length(vzped) ;

          pjump = (top - bottom) * 0.2 ;
          ptop = top - pjump * day ;

          kdeped = gpkde(vzped,0,[left;right]) ;
          kdeped = kdeped * length(vzped) / length(vz) ;

          hold on ;

            plot(mainxgrid,kdeped,[colstr '-']) ;
            plot([vzped'; vzped'], ...
                 [ptop; (ptop - pjump/2)]*ones(1,length(vzped)), ...
                                [colstr '-']) ;

              tx = left + .3 * (right - left) ;
              tstr = ['Day ' num2str(day) ', % within = ' num2str(pwi)] ;
            text(tx,ptop - .75 * pjump,tstr,'color',colstr) ;

          hold off ;


        end ;




        %  Make plot for Sessions
        %
        subplot(2,1,2) ;

        plot(mainxgrid,mainkde,'w-') ;
          set(gca,'Box','On') ;
          title([titstr ', ' date]) ;
          axis([left,right,bottom,top]) ;
          hold on ;
            plot([zmean,zmean],[bottom,top],'-w') ;
            plot([tolmin,tolmin],[bottom,top],'-w') ;
            plot([tolmax,tolmax],[bottom,top],'-w') ;
          hold off ;


        %  Loop through and add stuff for sessions
        %
        for sess = 1:3 ;

          if sess == 1 ;    
            colstr = 'c' ;
          elseif sess == 2 ;    
            colstr = 'm' ;
          elseif sess == 3 ;    
            colstr = 'g' ;
          end ;

          vzpes = vz(find(pevsessnum == sess)) ;

          pwi = 100 * sum((vzpes > tolmin) & ...
                              (vzpes < tolmax)) / length(vzpes) ;

          pjump = (top - bottom) * 0.2 ;
          ptop = top - pjump * sess ;

          kdeped = gpkde(vzpes,0,[left;right]) ;
          kdeped = kdeped * length(vzpes) / length(vz) ;

          hold on ;

            plot(mainxgrid,kdeped,[colstr '-']) ;
            plot([vzpes'; vzpes'], ...
                 [ptop; (ptop - pjump/2)]*ones(1,length(vzpes)), ...
                                [colstr '-']) ;

              tx = left + .3 * (right - left) ;
              tstr = ['Time ' num2str(sess) ', % within = ' num2str(pwi)] ;
            text(tx,ptop - .75 * pjump,tstr,'color',colstr) ;

          hold off ;


        end ;



          orient tall ;
          eval(['print -dpsc \matlab\cornean\vire10' ...
                                       num2str(iip) ...
                                       num2str(iieye) ...
                                       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 = [patname ', Eye ' num2str(iieye) ...
                      ', Zernike, Order ' num2str(ijacob) ...
                      ',' num2str(ifourier)] ;



          rtitstr = [patname ', Eye ' num2str(iieye)] ;
          ltitstr = ['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 = pemdata(:,icolc) ;    
            %  column vector of Zernike coefficients, cos
          vzs = pemdata(:,icols) ;    
            %  column vector of Zernike coefficients, sin


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



          rmean = mean(vzr) ;
          rstd = std(vzr) ;
  
          tolmin = rmean-tol*rstd ;
          tolmax = rmean+tol*rstd ;

          minr = min(vzr) ;
          maxr = max(vzr) ;
          radr = max([(rmean-minr); (maxr-rmean)]) ;
          left = rmean - 1.05 * radr ;
          right = rmean + 1.05 * radr ;




          %  make main distribution plot, for radii
          %
          [mainkde,mainxgrid] = gpkde(vzr,0,[left,right]) ;
          bottom = 0 ;
          top = 1.05 * max(mainkde) ;



          %  Make plot for days
          %
          clf ;
          subplot(2,2,1) ;

          plot(mainxgrid,mainkde,'w-') ;
            set(gca,'Box','On') ;
            title([ltitstr ', Magnitudes']) ;
            axis([left,right,bottom,top]) ;
            hold on ;
              plot([rmean,rmean],[bottom,top],'-w') ;
              plot([tolmin,tolmin],[bottom,top],'-w') ;
              plot([tolmax,tolmax],[bottom,top],'-w') ;
            hold off ;


          %  Loop through and add stuff for days
          %
          for day = 1:3 ;

            if day == 1 ;    
              colstr = 'c' ;
            elseif day == 2 ;    
              colstr = 'm' ;
            elseif day == 3 ;    
              colstr = 'g' ;
            end ;

            vrped = vzr(find(pevdaynum == day)) ;

            pwi = 100 * sum((vrped > tolmin) & ...
                                (vrped < tolmax)) / length(vrped) ;

            pjump = (top - bottom) * 0.2 ;
            ptop = top - pjump * day ;

            kdeped = gpkde(vrped,0,[left;right]) ;
            kdeped = kdeped * length(vrped) / length(vzr) ;

            hold on ;

              plot(mainxgrid,kdeped,[colstr '-']) ;
              plot([vrped'; vrped'], ...
                   [ptop; (ptop - pjump/2)]*ones(1,length(vrped)), ...
                                  [colstr '-']) ;

                tx = left + .3 * (right - left) ;
                tstr = ['Day ' num2str(day) ', % within = ' num2str(pwi)] ;
              text(tx,ptop - .75 * pjump,tstr,'color',colstr) ;

            hold off ;


          end ;




          %  Make plot for Sessions
          %
          subplot(2,2,3) ;

          plot(mainxgrid,mainkde,'w-') ;
            set(gca,'Box','On') ;
            title(['Magnitudes, ' date]) ;
            axis([left,right,bottom,top]) ;
            hold on ;
              plot([rmean,rmean],[bottom,top],'-w') ;
              plot([tolmin,tolmin],[bottom,top],'-w') ;
              plot([tolmax,tolmax],[bottom,top],'-w') ;
            hold off ;


          %  Loop through and add stuff for sessions
          %
          for sess = 1:3 ;

            if sess == 1 ;    
              colstr = 'c' ;
            elseif sess == 2 ;    
              colstr = 'm' ;
            elseif sess == 3 ;    
              colstr = 'g' ;
            end ;

            vrpes = vzr(find(pevsessnum == sess)) ;

            pwi = 100 * sum((vrpes > tolmin) & ...
                                (vrpes < tolmax)) / length(vrpes) ;

            pjump = (top - bottom) * 0.2 ;
            ptop = top - pjump * sess ;

            kdeped = gpkde(vrpes,0,[left;right]) ;
            kdeped = kdeped * length(vrpes) / length(vzr) ;

            hold on ;

              plot(mainxgrid,kdeped,[colstr '-']) ;
              plot([vrpes'; vrpes'], ...
                   [ptop; (ptop - pjump/2)]*ones(1,length(vrpes)), ...
                                  [colstr '-']) ;

                tx = left + .3 * (right - left) ;
                tstr = ['Time ' num2str(sess) ', % within = ' num2str(pwi)] ;
              text(tx,ptop - .75 * pjump,tstr,'color',colstr) ;

            hold off ;


          end ;






          tmean = mean(vzt) ;
          tstd = std(vzt) ;
  
          tolmin = tmean-tol*tstd ;
          tolmax = tmean+tol*tstd ;

          mint = min(vzt) ;
          maxt = max(vzt) ;
          radt = max([(tmean-mint); (maxt-tmean)]) ;
          left = tmean - 1.05 * radt ;
          right = tmean + 1.05 * radt ;




          %  make main distribution plot, for angles
          %
          [mainkde,mainxgrid] = gpkde(vzt,0,[left,right]) ;
          bottom = 0 ;
          top = 1.05 * max(mainkde) ;



          %  Make plot for days
          %
          subplot(2,2,2) ;

          plot(mainxgrid,mainkde,'w-') ;
            set(gca,'Box','On') ;
            title([rtitstr ', Angles']) ;
            axis([left,right,bottom,top]) ;
            hold on ;
              plot([tmean,tmean],[bottom,top],'-w') ;
              plot([tolmin,tolmin],[bottom,top],'-w') ;
              plot([tolmax,tolmax],[bottom,top],'-w') ;
            hold off ;


          %  Loop through and add stuff for days
          %
          for day = 1:3 ;

            if day == 1 ;    
              colstr = 'c' ;
            elseif day == 2 ;    
              colstr = 'm' ;
            elseif day == 3 ;    
              colstr = 'g' ;
            end ;

            vtped = vzt(find(pevdaynum == day)) ;

            pwi = 100 * sum((vtped > tolmin) & ...
                                (vtped < tolmax)) / length(vtped) ;

            pjump = (top - bottom) * 0.2 ;
            ptop = top - pjump * day ;

            kdeped = gpkde(vtped,0,[left;right]) ;
            kdeped = kdeped * length(vtped) / length(vzt) ;

            hold on ;

              plot(mainxgrid,kdeped,[colstr '-']) ;
              plot([vtped'; vtped'], ...
                   [ptop; (ptop - pjump/2)]*ones(1,length(vtped)), ...
                                  [colstr '-']) ;

                tx = left + .3 * (right - left) ;
                tstr = ['Day ' num2str(day) ', % within = ' num2str(pwi)] ;
              text(tx,ptop - .75 * pjump,tstr,'color',colstr) ;

            hold off ;


          end ;




          %  Make plot for Sessions
          %
          subplot(2,2,4) ;

          plot(mainxgrid,mainkde,'w-') ;
            set(gca,'Box','On') ;
            title(['Angles, ' date]) ;
            axis([left,right,bottom,top]) ;
            hold on ;
              plot([tmean,tmean],[bottom,top],'-w') ;
              plot([tolmin,tolmin],[bottom,top],'-w') ;
              plot([tolmax,tolmax],[bottom,top],'-w') ;
            hold off ;


          %  Loop through and add stuff for sessions
          %
          for sess = 1:3 ;

            if sess == 1 ;    
              colstr = 'c' ;
            elseif sess == 2 ;    
              colstr = 'm' ;
            elseif sess == 3 ;    
              colstr = 'g' ;
            end ;

            vtpes = vzt(find(pevsessnum == sess)) ;

            pwi = 100 * sum((vtpes > tolmin) & ...
                                (vtpes < tolmax)) / length(vtpes) ;

            pjump = (top - bottom) * 0.2 ;
            ptop = top - pjump * sess ;

            kdeped = gpkde(vtpes,0,[left;right]) ;
            kdeped = kdeped * length(vtpes) / length(vzt) ;

            hold on ;

              plot(mainxgrid,kdeped,[colstr '-']) ;
              plot([vtpes'; vtpes'], ...
                   [ptop; (ptop - pjump/2)]*ones(1,length(vtpes)), ...
                                  [colstr '-']) ;

                tx = left + .3 * (right - left) ;
                tstr = ['Time ' num2str(sess) ', % within = ' num2str(pwi)] ;
              text(tx,ptop - .75 * pjump,tstr,'color',colstr) ;

            hold off ;


          end ;




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





        end ;    %  of loop through Fourier orders


      end ;    %  of loop through sumind


    end ;    %  of loop through eyes



  end ;    % of loop through patients



end ;  


