disp('Running MATLAB script file zarinafig3.m') ;
%  makes figure 2 for Zarina's master paper


%  copied from vire1.m
%

%    Reads already stored Zernike Decomps, and saves as
%    ASCII files, for e.g. SAS analysis
%    Also does first analyses
%


ipat = 0 ;       %  0 - loop through all
                 %  1 - Hickson
                 %  2 - Meadows
                 %  3 - Tesoro 
                 %  4 - Yi     


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


%  Initialize parameters
%
if ipat == 0 ;     %  then loop through all patients
  ips = 1 ;
  ipe = 4 ;
else ;             %  then just do input patients
  ips = ipat ;
  ipe = ipat ;
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 

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

INDATPATH = '\matlab\cornean\indata\' ;
          %  path for input Keratron files
[COMAP, HICURVE] = coacmap ;





  %  First read in data file, created by vire1.m
  %
  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 ;

  [inmat,cnt] = fscanf(fid,forstr) ;
  fclose(fid) ;

  ncol = ninpar + ncoeff ;
  nrow = cnt / ncol ;

  inmat = reshape(inmat,ncol,nrow)' ;
          %   since reshape works column wise

%  inmat(1:10,:)
%  inmat(nrow,:)
          %  used for testing read in



  grandmean = mean(inmat(:,(ninpar+1):ncol)) ;
          %  vector of coefficents of mean of everything


  figure(1) ;
  clf ;



  %  Main loop throught patients
  %
  for ip = ips:ipe ;

    if ip == 1 ;    
      patname = 'Hick' ;
    elseif ip == 2 ;    
      patname = 'Mead' ;
    elseif ip == 3 ;    
      patname = 'Teso' ;
    elseif ip == 4 ;    
      patname = 'Yi' ;
    end ;

    pflag = (inmat(:,2) == ip) ;
          %  2nd column is patient number
    pdat = inmat(pflag,:) ;
          %  matrix of data for this patient
    pmean = mean(pdat(:,(ninpar+1):ncol)) ;
          %  vector of coefficents of mean for this patient


    for ieye = 1:2 ;      %  loop through eyes

      peflag = (pdat(:,5) == ieye) ;
          %  5th column is eye number
      pedat = pdat(peflag,:) ;
          %  matrix of data for this patient, this eye
      pemean = mean(pedat(:,(ninpar+1):ncol)) ;
          %  vector of coefficents of mean for this patient, this eye


disp(['    Patient ' num2str(ip) '  Eye ' num2str(ieye)]) ;



      %  Make overall Mean plot
      %
      plotind = 6 * ieye + ip + 1 ;
      subplot(4,6,plotind) ;
      titstr = [patname ' E_' num2str(ieye) ' Avg.'] ;

      %  Create image with Keratron Color scheme, and put in current axes
      %
      betahat = pemean' ;
      imagout = coazr(betahat) ;
      pcolor(xgrid,ygrid,imagout) ;
        shading flat ;
        axis('square') ;
titstr = ['Avg. Pat. ' num2str(ip) ' Eye ' num2str(ieye)] ;
        title(titstr) ;
          caxis([0 1]) ;
        colormap(COMAP) ;

%      hold on ;
%        tx = -.95 * R0 ;
%        ty = .7 * R0 ;
%        nimag = size(pedat,1) ;
%        iss = nimag * sum((pemean - grandmean).^2) ;
%        tss = sum(sum((pedat(:,(ninpar+1):ncol) - ...
%                                 vec2mat(grandmean,nimag)).^2)) ;
%        ssstr = num2str(100 * iss / tss) ;
%        tstr = [ssstr '% of GMR'] ;
%        text(tx,ty,tstr) ;
%      hold off ;


        orient landscape ;
      print -dpsc \matlab\steve\ps\zarinafig2.ps ;


    end ;



  end ;




