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

%    Reads already stored Zernike Decomps, 
%    and does graphics about them individually
%    Modifies vire6.m, by breaking out day and
%    session terms....



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


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


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



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

        zmedian = median(vz) ;

        vzord = sort(vz) ;
        z10 = vzord(round(0.10 * nrowpe)) ;
        z25 = vzord(round(0.25 * nrowpe)) ; 
        z75 = vzord(round(0.75 * nrowpe)) ;
        z90 = vzord(round(0.90 * nrowpe)) ;


        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 plots
        %
        for day = 1:3 ;

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

          title(['Dist''s for ' titstr ...
                             ', Day ' num2str(day) ', ' date]) ;

          hold on ;
            plot([z10; z10],[bottom; top],'w-') ;
            plot([z25; z25],[bottom; top],'w-') ;
            plot([zmedian; zmedian],[bottom; top],'w-') ;
            plot([z75; z75],[bottom; top],'w-') ;
            plot([z90; z90],[bottom; top],'w-') ;
          hold off ;


          mkdepeds = [] ;
          vzpeds05 = [] ; 
          vzpedsmed = [] ; 
          vzpeds95 = [] ; 

          for sess = 1:3 ;

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

            vzpeds = vz(find((pevdaynum == day) & ...
                                         (pevsessnum == sess))) ;

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


            kdepeds = gpkde(vzpeds,0,[left;right]) ;
            mkdepeds = [mkdepeds, kdepeds] ;


            hold on ;
              plot([vzpeds'; vzpeds'], ...
                   [ptop; (ptop - pjump)]*ones(1,length(vzpeds)), ...
                                  [colstr '-']) ;

                tx = left + .9 * (right - left) ;
              text(tx,ptop - 0.5 * pjump,['sess ' num2str(sess)],'color',colstr) ;

            hold off ;



            nzpeds = length(vzpeds) ;
            vzpedso = sort(vzpeds) ;
            zpeds05 = vzpedso(round(0.05 * nzpeds)) ;
            zpedsmed = median(vzpedso) ;
            zpeds95 = vzpedso(round(0.95 * nzpeds)) ;

            vzpeds05 = [vzpeds05 ; zpeds05] ; 
            vzpedsmed = [vzpedsmed ; zpedsmed] ; 
            vzpeds95 = [vzpeds95 ; zpeds95] ; 


          end ;


          mkdepeds = 0.58 * mkdepeds / max(max(mkdepeds)) ;
            %  make max ht of each = 0.58 for plotting
          kdegrid = linspace(left,right,401)' ;



          for sess = 1:3 ;

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

            hold on ;

              plot(kdegrid,mkdepeds(:,sess),[colstr '-']) ;

              pjump = (top - bottom) * 0.1 ;
              pbot = bottom + pjump * sess ;

              plot(vzpedsmed(sess),pbot,[colstr 'o']) ;
                vachil = get(gca,'Children') ;
                set(vachil(1),'LineWidth',2) ;
                set(vachil(1),'MarkerSize',8) ;

              plot([vzpeds05(sess); vzpeds95(sess)], ...
                              [pbot; pbot],[colstr '-']) ;
                vachil = get(gca,'Children') ;
                set(vachil(1),'LineWidth',3) ;

            hold off ;
      

          end ;

        end ;
    



          orient landscape ;
        eval(['print -dpsc \matlab\cornean\vire7' 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)] ;


          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 ;



          %  make main distribution plot, for radii
          %
          zmedian = median(vzr) ;


          vzord = sort(vzr) ;
          z10 = vzord(round(0.10 * nrowpe)) ;
          z25 = vzord(round(0.25 * nrowpe)) ; 
          z75 = vzord(round(0.75 * nrowpe)) ;
          z90 = vzord(round(0.90 * nrowpe)) ;


          minz = min(vzr) ;
          maxz = max(vzr) ;


          clf ;


          for day = 1:3 ;

            subplot(6,1,day) ;

            left = - 0.05 * maxz ;
            right = 1.05 * maxz ;
            bottom = 0 ;
            top = 1 ;
            center = mean([top; bottom]) ;

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


            title(['Radii for ' titstr ...
                             ', Day ' num2str(day) ', ' date]) ;

            hold on ;
              plot([z10; z10],[bottom; top],'w-') ;
              plot([z25; z25],[bottom; top],'w-') ;
              plot([zmedian; zmedian],[bottom; top],'w-') ;
              plot([z75; z75],[bottom; top],'w-') ;
              plot([z90; z90],[bottom; top],'w-') ;
            hold off ;


            mkdepeds = [] ;
            vzpeds05 = [] ; 
            vzpedsmed = [] ; 
            vzpeds95 = [] ; 

            for sess = 1:3 ;

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


              vzpeds = vzr(find((pevdaynum == day) & ...
                                         (pevsessnum == sess))) ;

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


              kdepeds = gpkde(vzpeds,0,[left;right]) ;
              mkdepeds = [mkdepeds, kdepeds] ;


              hold on ;
                plot([vzpeds'; vzpeds'], ...
                     [ptop; (ptop - pjump)]*ones(1,length(vzpeds)), ...
                                  [colstr '-']) ;

                  tx = left + .9 * (right - left) ;
                text(tx,ptop - 0.5 * pjump,['sess ' num2str(sess)],'color',colstr) ;

              hold off ;


              nzpeds = length(vzpeds) ;
              vzpedso = sort(vzpeds) ;
              zpeds05 = vzpedso(round(0.05 * nzpeds)) ;
              zpedsmed = median(vzpedso) ;
              zpeds95 = vzpedso(round(0.95 * nzpeds)) ;

              vzpeds05 = [vzpeds05 ; zpeds05] ; 
              vzpedsmed = [vzpedsmed ; zpedsmed] ; 
              vzpeds95 = [vzpeds95 ; zpeds95] ; 


            end ;


            mkdepeds = 0.58 * mkdepeds / max(max(mkdepeds)) ;
              %  make max ht of each = 0.58 for plotting
            kdegrid = linspace(left,right,401)' ;


            for sess = 1:3 ;

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

              hold on ;

                plot(kdegrid,mkdepeds(:,sess),[colstr '-']) ;

                pjump = (top - bottom) * 0.1 ;
                pbot = bottom + pjump * sess ;

                plot(vzpedsmed(sess),pbot,[colstr 'o']) ;
                  vachil = get(gca,'Children') ;
                  set(vachil(1),'LineWidth',2) ;
                  set(vachil(1),'MarkerSize',8) ;

                plot([vzpeds05(sess); vzpeds95(sess)], ...
                                [pbot; pbot],[colstr '-']) ;
                  vachil = get(gca,'Children') ;
                  set(vachil(1),'LineWidth',3) ;

              hold off ;
      

            end ;

    


            %  make main distribution plot, for angles
            %
            subplot(6,1,(day + 3)) ;
            set(gca,'Box','On') ;
  
              left = -180 ;
              right = 540 ;
              kdeleft = -360 ;
              kderight = 720 ;
            axis([left,right,bottom,top]) ;


            title(['Angles for ' titstr ...
                             ', Day ' num2str(day) ', ' date]) ;


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


            mkdepeds = [] ;

            for sess = 1:3 ;

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


              vzpeds = vzt(find((pevdaynum == day) & ...
                                         (pevsessnum == sess))) ;
              vzpeds = [(vzpeds - 360); vzpeds; (vzpeds + 360)] ;

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


              kdepeds = gpkde(vzpeds,0,[left;right]) ;
              mkdepeds = [mkdepeds, kdepeds] ;


              hold on ;
                plot([vzpeds'; vzpeds'], ...
                     [ptop; (ptop - pjump)]*ones(1,length(vzpeds)), ...
                                  [colstr '-']) ;

                  tx = left + .9 * (right - left) ;
                text(tx,ptop - 0.5 * pjump,['sess ' num2str(sess)], ...
                                  'color',colstr) ;

              hold off ;


              nzpeds = length(vzpeds) ;
              vzpedso = sort(vzpeds) ;
              zpeds05 = vzpedso(round(0.05 * nzpeds)) ;
              zpedsmed = median(vzpedso) ;
              zpeds95 = vzpedso(round(0.95 * nzpeds)) ;

              vzpeds05 = [vzpeds05 ; zpeds05] ; 
              vzpedsmed = [vzpedsmed ; zpedsmed] ; 
              vzpeds95 = [vzpeds95 ; zpeds95] ; 


            end ;


            mkdepeds = 0.58 * mkdepeds / max(max(mkdepeds)) ;
              %  make max ht of each = 0.58 for plotting
            kdegrid = linspace(left,right,401)' ;


            for sess = 1:3 ;

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

              hold on ;

                plot(kdegrid,mkdepeds(:,sess),[colstr '-']) ;

              hold off ;
      

            end ;

          end ;
    


            orient tall ;
          eval(['print -dpsc \matlab\cornean\vire7' 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 ;  


