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

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



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


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



npc = 4  ;     %  number of PCs 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\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) ;


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


if ipart == 1 ;    %  then make main graphic


  for ipc = 1:npc ;


    titstr = ['PC ' num2str(ipc)] ;


    disp(['  Working on ' titstr]) ;


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

    vpc = mdata(:,icol) ;    
          %  column vector of PCs

    pcmedian = median(vpc) ;

    vpcord = sort(vpc) ;
    pc02 = vpcord(round(0.02 * nrow)) ;
    pc25 = vpcord(round(0.25 * nrow)) ; 
    pc75 = vpcord(round(0.75 * nrow)) ;
    pc98 = vpcord(round(0.98 * nrow)) ;



    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 plot
    %
    subplot(2,1,1) ;
      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 ;

      
      vpcpe1 = vpc(find((vpatnum == ip) & (veyenum == 1))) ;
      vpcpe2 = vpc(find((vpatnum == ip) & (veyenum == 2))) ;

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

      hold on ;
        plot([vpcpe1'; vpcpe1'], ...
                 [pbot; (pbot + pjump)]*ones(1,length(vpcpe1)), ...
                                [colstr '-']) ;
        plot([vpcpe2'; vpcpe2'], ...
                 [ptop; (ptop - pjump)]*ones(1,length(vpcpe2)), ...
                                [colstr '-']) ;
        plot([pc02; pc02],[bottom; top],'w-') ;
        plot([pc25; pc25],[bottom; top],'w-') ;
        plot([pcmedian; pcmedian],[bottom; top],'w-') ;
        plot([pc75; pc75],[bottom; top],'w-') ;
        plot([pc98; pc98],[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(vpcpe1,0,[left;right]) ;
      kde2 = gpkde(vpcpe2,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(2,1,2) ;
      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 ;
    




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


  end ;

end ;


