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

%    Reads already stored Zernike PCA Coeffs, 
%    and does graphics about them individually
%    Modifies vire8.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


npc = 4  ;     %  number of PCs to treat



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 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
%
%  Read in full data matrix
%
fid = fopen(['\matlab\cornean\vire2' 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 = 20 ;
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

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


  pemdata = mdata(find(vpatnum == ip & ...
                                 veyenum == ieye),:) ;
          %  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 Principal Components
  %
  for ipc = 1:npc ;


    titstr = [patname ', Eye ' num2str(ieye) ...
                  ', PC ' num2str(ipc)] ;


    disp(['  Working on ' titstr]) ;


    icol = ninpar + ipc ;
          %  index of column for current PC

    vpc = pemdata(:,icol) ;    
          %  column vector of PC coefficients

    pcmedian = median(vpc) ;

    vpcord = sort(vpc) ;
    pc10 = vpcord(round(0.10 * nrowpe)) ;
    pc25 = vpcord(round(0.25 * nrowpe)) ; 
    pc75 = vpcord(round(0.75 * nrowpe)) ;
    pc90 = vpcord(round(0.90 * nrowpe)) ;


    minpc = min(vpc) ;
    maxpc = max(vpc) ;
    radpc = max([(pcmedian-minpc); (maxpc-pcmedian)]) ;
    left = pcmedian - 1.05 * radpc ;
    right = pcmedian + 1.05 * radpc ;
    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([pc10; pc10],[bottom; top],'w-') ;
        plot([pc25; pc25],[bottom; top],'w-') ;
        plot([pcmedian; pcmedian],[bottom; top],'w-') ;
        plot([pc75; pc75],[bottom; top],'w-') ;
        plot([pc90; pc90],[bottom; top],'w-') ;
      hold off ;


      mkdepeds = [] ;
      vpcpeds05 = [] ; 
      vpcpedsmed = [] ; 
      vpcpeds95 = [] ; 

      for sess = 1:3 ;

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

        vpcpeds = vpc(find((pevdaynum == day) & ...
                                   (pevsessnum == sess))) ;

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


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


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

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

        hold off ;



        npcpeds = length(vpcpeds) ;
        vpcpedso = sort(vpcpeds) ;
        pcpeds05 = vpcpedso(round(0.05 * npcpeds)) ;
        pcpedsmed = median(vpcpedso) ;
        pcpeds95 = vpcpedso(round(0.95 * npcpeds)) ;

        vpcpeds05 = [vpcpeds05 ; pcpeds05] ; 
        vpcpedsmed = [vpcpedsmed ; pcpedsmed] ; 
        vpcpeds95 = [vpcpeds95 ; pcpeds95] ; 


      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(vpcpedsmed(sess),pbot,[colstr 'o']) ;
            vachil = get(gca,'Children') ;
            set(vachil(1),'LineWidth',2) ;
            set(vachil(1),'MarkerSize',8) ;

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

        hold off ;
      

      end ;

    end ;
    


      orient tall ;
    eval(['print -dpsc \matlab\cornean\vire9' num2str(ip) ...
                            num2str(ieye) num2str(ipc) '.ps']) ;



  end ;   %  of loop through PC's


end ;  


