disp('Running MATLAB script file vimu3.m') ;
%
%    For Analysis of Vistakon Multifocal Data
%    This does work with no lens images,
%    but centered according to each lnes image
%



ipart = 2 ;      %  1 - Loop through input files, 
                 %        make special centered Zernike decomp.
                 %  2 - Make simple Zernike recon. pic




global INDATPATH NMAX R0 HICURVE IBASIS ;
INDATPATH = '\matlab\cornean\indata\' ;
          %  path for input Keratron files
[COMAP, HICURVE] = coacmap ;
ncomap = size(COMAP,1) ;
          %  number of rows in the color map ;



if ipart == 1 ;      %  read in test images, do special
                     %     Radial decompositions
  ngrid = 64 ;
  nang = 256;    % number of semi-meridians read per radial ring
  R0 = 4.2 ;
  NMAX = 10 ;
  IRECENT = 1 ;    %  work with recentering


  for ifile = 1:10 ;    %  same ordering as in vimu2.m

    disp(['Working on file #' num2str(ifile)]) ;

    if ifile == 1 ;
      titstr = 'Ross,L,304' ;
%      filestr = '000357_2' ;
      filestr = '000352_5' ;
          %  replace filestr with corresponding no lens file
      pupoffamt = 0.119 ;
      pupoffdir = 125 ;
      lencenamt = .99 ;
      lencendir = 50 ;
      lflag = 1 ;
          %  1 for wearing lens
      nlpupoffamt = 0.13 ;
      nlpupoffdir = 101 ;
          %  corresponding no lens pupil center
    elseif ifile == 2 ;
      titstr = 'Ross,L,No' ;
      filestr = '000352_5' ;
      pupoffamt = 0.13 ;
      pupoffdir = 101 ;
      lencenamt = 0 ;
      lencendir = 0 ;
      lflag = 0 ;
          %  1 for wearing lens
    elseif ifile == 3 ;
      titstr = 'Ross,R,No' ;
      filestr = '000352_2' ;
      pupoffamt = 0.248 ;
      pupoffdir = 139 ;
      lencenamt = 0 ;
      lencendir = 0 ;
      lflag = 0 ;
          %  1 for wearing lens
    elseif ifile == 4 ;
      titstr = 'Ross,R,301' ;
%      filestr = '000357_4' ;
      filestr = '000352_2' ;
          %  replace filestr with corresponding no lens file
      pupoffamt = 0.189 ;
      pupoffdir = 144 ;
      lencenamt = 1.3 ;
      lencendir = 118 ;
      lflag = 1 ;
          %  1 for wearing lens
      nlpupoffamt = 0.248 ;
      nlpupoffdir = 139 ;
          %  corresponding no lens pupil center
    elseif ifile == 5 ;
      titstr = 'Duer,L,304' ;
%      filestr = '000358_2' ;
      filestr = '000353_2' ;
          %  replace filestr with corresponding no lens file
      pupoffamt = 0.085 ;
      pupoffdir = 300 ;
      lencenamt = .45 ;
      lencendir = 301 ;
      lflag = 1 ;
          %  1 for wearing lens
      nlpupoffamt = 0.109 ;
      nlpupoffdir = 287 ;
          %  corresponding no lens pupil center
    elseif ifile == 6 ;
      titstr = 'Duer,L,No' ;
      filestr = '000353_2' ;
      pupoffamt = 0.109 ;
      pupoffdir = 287 ;
      lencenamt = 0 ;
      lencendir = 0 ;
      lflag = 0 ;
          %  1 for wearing lens
    elseif ifile == 7 ;
      titstr = 'Duer,R,No' ;
      filestr = '000356_3' ;
      pupoffamt = 0.059 ;
      pupoffdir = 187 ;
      lencenamt = 0 ;
      lencendir = 0 ;
      lflag = 0 ;
          %  1 for wearing lens
    elseif ifile == 8 ;
      titstr = 'Duer,R,301' ;
%      filestr = '000358_1' ;
      filestr = '000356_3' ;
          %  replace filestr with corresponding no lens file
      pupoffamt = 0.099 ;
      pupoffdir = 244 ;
      lencenamt = .89 ;
      lencendir = 145 ;
      lflag = 1 ;
          %  1 for wearing lens
      nlpupoffamt = 0.059 ;
      nlpupoffdir = 187 ;
          %  corresponding no lens pupil center
    elseif ifile == 9 ;
      titstr = 'Duer,L,303' ;
%      filestr = '000360_1' ;
      filestr = '000353_2' ;
          %  replace filestr with corresponding no lens file
      pupoffamt = 0.116 ;
      pupoffdir = 290 ;
      lencenamt = .89 ;
      lencendir = 14 ;
      lflag = 1 ;
          %  1 for wearing lens
      nlpupoffamt = 0.109 ;
      nlpupoffdir = 287 ;
          %  corresponding no lens pupil center
    elseif ifile == 10 ;
      titstr = 'Duer,R,302' ;
%      filestr = '000360_4' ;
      filestr = '000356_3' ;
          %  replace filestr with corresponding no lens file
      pupoffamt = 0.067 ;
      pupoffdir = 261 ;
      lencenamt = .77 ;
      lencendir = 169 ;
      lflag = 1 ;
          %  1 for wearing lens
      nlpupoffamt = 0.059 ;
      nlpupoffdir = 187 ;
          %  corresponding no lens pupil center
    end ;


    if lflag == 1 ;     %  then looking at lens wearing ifile, so proceed


      %  use recenterings from this image, to get shift to lens center
      %
      laa = 2 * pi * lencendir / 360 ;
          %  angle in radians, not degrees
          %  lens center angle, from apex normal
      lra = lencenamt ;
          %  lens center radius, from apex normal
      lxa = lra * cos(laa) ;
      lya = lra * sin(laa) ;
          %  lens center x,y coord's, from apex normal

      paa = 2 * pi * pupoffdir / 360 ;
          %  angle in radians, not degrees
          %  pupil center angle, from apex normal
      pra = pupoffamt ;
          %  pupil center radius, from apex normal
      pxa = pra * cos(paa) ;
      pya = pra * sin(paa) ;
          %  pupil center x,y coord's, from apex normal
    
      lxp = lxa - pxa ;
      lyp = lya - pya ;
          %  lens center x,y coord's, from pupil center


      %  get pupil center of no lens image
      %
      pa = 2 * pi * nlpupoffdir / 360 ;
          %  angle in radians, not degrees
          %  pupil center angle, from apex normal
      pr = nlpupoffamt ;
          %  pupil center radius, from apex normal
      px = pr * cos(pa) ;
      py = pr * sin(pa) ;
          %  pupil center x,y coord's, from apex normal


      %  get point to recenter to:
      %
      xcent = px + lxp ;
      ycent = py + lyp ;



      %  Read in data
      %
      RFILESUF = '.xlb' ;
          %  radii data file suffix 
      HFILESUF = '.zlb' ;
          %  height data file suffix
      CFILESUF = '.tlb' ;
          %  curvature data file suffix
      RISTR = ['\matlab\cornean\indata\vimu\' filestr RFILESUF] ;
      HISTR = ['\matlab\cornean\indata\vimu\' filestr HFILESUF] ;
      CISTR = ['\matlab\cornean\indata\vimu\' filestr CFILESUF] ;


      fidr = fopen(RISTR, 'r', 'l') ;
         %  IEEE format for binary files
      fidh = fopen(HISTR, 'r', 'l') ;
      fidc = fopen(CISTR, 'r', 'l') ;
         %  assigns a file handle, (read only assumed)
      if (fidr == -1) | (fidh == -1) | (fidc == -1) ;
                               %  then a file not found, so reformat
        if (upper(RFILESUF) == RFILESUF);
           RFILESUF = lower(RFILESUF);		% handles case sensitivity
           HFILESUF = lower(HFILESUF);		% .xlb vs. .XLB
           CFILESUF = lower(CFILESUF)
        else				
           RFILESUF = upper(RFILESUF);
           HFILESUF = upper(HFILESUF);
           CFILESUF = upper(CFILESUF);
        end;
        newristr = ['\matlab\cornean\indata\vimu\' filestr RFILESUF];
        newhistr = ['\matlab\cornean\indata\vimu\' filestr HFILESUF];
        newcistr = ['\matlab\cornean\indata\vimu\' filestr CFILESUF];
        fidr = fopen(newristr, 'r', 'l');
        fidh = fopen(newhistr, 'r', 'l');
        fidc = fopen(newcistr, 'r', 'l');
	      %   try binary open with different suffix
      end;

      if (fidr == -1) | (fidc == -1) | (fidc == -1) ;
                                 %  then a file not found, so try ascii
        newrsuf = upper([RFILESUF(1:3) 'a']);
        newhsuf = upper([HFILESUF(1:3) 'a']);
        newcsuf = upper([CFILESUF(1:3) 'a']);
        newristr = ['\matlab\cornean\indata\vimu\' filestr newrsuf];
        newhistr = ['\matlab\cornean\indata\vimu\' filestr newhsuf];
        newcistr = ['\matlab\cornean\indata\vimu\' filestr newcsuf];
        fidr = fopen(newristr);
        fidh = fopen(newhistr);
        fidc = fopen(newcistr);
	        %   new file handles for ascii format

        if (fidr == -1) | (fidc == -1) | (fidc == -1) ;
                                 %  then a file not found, so reformat
          newrsuf = lower(newrsuf);		% handles case sensitivity
          newhsuf = lower(newhsuf);		% .xla vs. .XLA
          newcsuf = lower(newcsuf);
          newristr = ['\matlab\cornean\indata\vimu\' filestr newrsuf];
          newhistr = ['\matlab\cornean\indata\vimu\' filestr newhsuf];
          newcistr = ['\matlab\cornean\indata\vimu\' filestr newcsuf];
          fidr = fopen(newristr);
          fidh = fopen(newhistr);
          fidc = fopen(newcistr);
	       %   try ascii open with different suffix
        end;

        if (fidr == -1) | (fidh == -1) | (fidc == -1) ;
                                  % then file not found, so give error
          disp('   ') ;
          disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
          disp('!!!   Can''t do raw data display, because can''t find file(s):') ;
    
          if (fidr == -1) ;
            disp(['!!!       ' RISTR]) ;
          end ;
 
          if (fidh == -1) ;
            disp(['!!!       ' HISTR]) ;
          end ;

          if (fidc == -1) ;
            disp(['!!!       ' CISTR]) ;
          end ;

          disp('!!!') ;
          disp('!!!   Terminating raw data display') ;
          disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;

        else ;

          %  Continue with ascii data read, do some data preprocessing
          %
          okread = 1; 
          disp('Performing ASCII file read for the input data.');
          disp('');
          rdata = fscanf(fidr,'%g') ;
               %  formatted reads, assumes single column
          cdata = fscanf(fidc,'%g') ;
          hdata = fscanf(fidh,'%g') ;
          fclose(fidr) ;
          fclose(fidh) ;
          fclose(fidc) ;
     
          if rdata(2) == 0 ;		% (no flag on front of data)
            nring = (length(rdata)/nang) + 1;
            if nring - floor(nring) > 0
              disp('!!! File not read in properly, dimensions are off  !!!');
              okread = 0;
            end;

          else
            nring = rdata(2) + 1;
            rdata(1:2) = [];
            cdata(1:2) = [];
            hdata(1:2) = [];
            if nring - floor(nring) > 0
              disp('!!! File not read in properly, dimensions are off  !!!');
              okread = 0;
            end;  
          end;
        end; 

      else ;


        %  Continue with binary data read, do some data preprocessing
        %
        okread = 1;
        disp('Performing binary file read for the input data.');
        disp('');
        rstatus = fseek(fidr, 2, -1);
        rdata = fread(fidr, 'float');
        fclose(fidr);
 
        hstatus = fseek(fidh, 2, -1);
        hdata = fread(fidh, 'float');
        fclose(fidh);

        cstatus = fseek(fidc, 2, -1);
        cdata = fread(fidc, 'float');
        fclose(fidc);

        if rdata(1) == 0 ;		% (no flag on front of data)
          nring = (length(rdata)/nang) + 1;
          if nring - floor(nring) > 0 ; 
            disp('!!! File not read in properly, dimensions are off  !!!');	
            okread = 0;
          end;
        else
          nring = rdata(1) + 1;
          rdata(1) = [];
          hdata(1) = [];
          cdata(1) = [];
          if nring - floor(nring) > 0 ;
            disp('!!! File not read in properly, dimensions are off  !!!');
            okread = 0;
          end;  
        end;
      end;



      rdata = reshape(rdata,nang,nring) ;
      hdata = reshape(hdata,nang,nring) ;
      cdata = reshape(cdata,nang,nring) ;

      adata = linspace(0, (1 - 1/nang) * 2 * pi, nang) ;
      adata = adata' * ones(1,nring) ;
          %  Corresponding matrix of angles.

      rdata = rdata(:,2:nring) ;
      adata = adata(:,2:nring) ;
      hdata = hdata(:,2:nring) ;
      cdata = cdata(:,2:nring) ;
          %  cutoff first columns of all zeros

      %  set values outside R0 to missing
      % 
      vor0ind = (rdata > R0) ;
           %  ones where rdata is outside R0
      nor0ind = sum(sum(vor0ind)) ;
      rdata(vor0ind) = zeros(nor0ind,1) ;
      adata(vor0ind) = zeros(nor0ind,1) ;
      hdata(vor0ind) = zeros(nor0ind,1) ;
      cdata(vor0ind) = zeros(nor0ind,1) ;

      %#############
      % !!! Note:  rinner and cinner draw an inside line for possible missing
      % 	     data on the interior (reinhard).
      % inedger and inedgec have been commented out, and may be used in the
      % future for other types of analyses.

      %  construct vector of data for "outer ring"
      %
      [rinner, rdvmvind] = coamdf(rdata, nring-1, nang) ;
          %  vector of indices before first missing values for rdata
      [cinner, cdvmvind] = coamdf(cdata, nring-1, nang) ;
          %  vector of indices before first missing values for hdata
      vmvind = min([rdvmvind'; cdvmvind'])' ;
          %  vector of indices before first missing values (exterior)
      vinner = max([rinner'; cinner'])' ;
          %  vector of indices after first missing values (interior)
      edgerdata = [] ;
      edgehdata = [] ;
      edgecdata = [] ;
      %  inedger = [];
      %  inedgec = [];
      for iang = 1:nang ;
      %    inedger   = [inedger; rdata(iang, vinner(iang, 1))] ;
      %    inedgec   = [inedgec; cdata(iang, vinner(iang, 1))] ;
        edgerdata = [edgerdata; rdata(iang,vmvind(iang,1))] ;
        edgehdata = [edgehdata; hdata(iang,vmvind(iang,1))] ;
        edgecdata = [edgecdata; cdata(iang,vmvind(iang,1))] ;
      end ;
          %  take values where are inside R0, 
          %  and are before first non-missing
      edgeadata = adata(:,1) ;
          %  each angle

      %##############



      %  remove inappropriate values
      %
      %  first missings in radii files
      vind = find(rdata == 0) ;
           %  vector of indices where have a missing
      rdata(vind) = [] ;
      adata(vind) = [] ;
      hdata(vind) = [] ;
      cdata(vind) = [] ;
           %  Deletes that part of the vectors

      %  now missings in curvature files
      vind = find(cdata == 0) ;
           %  vector of indices where have a missing
      rdata(vind) = [] ;
      adata(vind) = [] ;
      hdata(vind) = [] ;
      cdata(vind) = [] ;
           %  Deletes that part of the vectors
      %disp(['After removing missings, have ' num2str(length(rdata)) ' values']) ;




      %  Now do radial, centered, Zernike decomp., 
      %  using lines from coardec.m
      %
      rdata = rdata' ;
      adata = adata' ;
      hdata = hdata' ;
      cdata = cdata' ;
       %  First turn into column vectors


      %###  Now have original data with all missing values thrown out.
      %###  Need to transform this data, and then get rid of all values
      %###	which are outside of R0.
      %###  also need to transform the boundary region for the plot

      %  #########################################

      if (IRECENT == 1) ;   %  then do recentering

        angc = CANG * pi/180;		%% put on radian scale

        [xdata, ydata] = pol2cart(adata, rdata);
				% converts polar data coords to cartesian
        away = sqrt((xdata - xcent).^2 + (ydata - ycent).^2);
				% distance from each point to center
        [away2, index] = sort(away);	
	         % gives positions of minimum distances
        min1 = [xdata(index(1)), ydata(index(1)), hdata(index(1))];

        j = 2;
        samept = 1;

        while samept ==1;	%% checks that the x/y are not the same (origin)
	       min2 = [xdata(index(j)), ydata(index(j)), hdata(index(j))];
          if(min1(1)==min2(1) & min1(2)==min2(2))
	         j = j + 1;
          else
            samept = 0;
            v1 = min2 - min1;	% first vec to determine plane
          end
        end

        colinear = 1;		% error check for 3 points in plane  
        j = 3;			% start with third smallest distance

        while colinear == 1
          min3 = [xdata(index(j)), ydata(index(j)), hdata(index(j))];
					% try j-th closest to make plane
          v2 = min3 - min1;	% second vector to determine plane
          if rank([v1;v2]) == 2
		      colinear = 0 ;	% found independent vector
          else 
	         j = j + 1;	% found colinear vector
          end
        end

        det1 = [1 0 0]; det2 = [0 1 0]; det3= [0 0 1];
        u1 = det([det1; v1; v2]);
        u2 = det([det2; v1; v2]);
        u3 = det([det3; v1; v2]);
        u = [u1; u2; u3];		
					% normal vector (yields ax+by+cz = ?)
        const = min1 * u;		% constant for ? above
        if u3 ~= 0
          u = -u/u3;		% rescale so that Z-direction is +
          const = -const/u3;	%  "
        end

        %% now normalize the vector

        u = u / norm(u);
        const = const / norm(u);

        %% find approximate height at pupil center
        %% solve for c in the ax+by+cz=d equation found above
        %% take care that min(3)'s third value can be zero!

        hold2 = const - min3(1:2)*u(1:2);
        zp = hold2 / min3(3);

        if isnan(zp)
          hold3 = const - min2(1:2)*u(1:2);
          zp = hold3 / min2(3);		% need another check here?
        end;

        pupil = [xcent ycent zp];		% coordinates of pupil center point.

        %% find pupil center axes
        %% 

        e1 = [1 0 0]'; 
        e2 = [0 1 0]';

        e1hat = e1 - (u'*e1)*u;
        e1hat = e1hat / norm(e1hat);		% normalized "nasal axis"

        e2hat = e2 - (u'*e2)*u - (e1hat'*e2)*e1hat;
        e2hat = e2hat / norm(e2hat);		% normalized "superior axis"

        %% take data and transform the heights H to heights H-twiddle
        %%

        xstar = xcent - xdata;		
        ystar = ycent - ydata;
        zstar = zp - hdata;

        htwiddle = [xstar ystar zstar]*u;	% get transformed heights

        %% get transformed x- and y-axes.
        %% 
        hatmatrix = eye(3) - u*u';
        ztwiddle = hatmatrix * [-xstar -ystar -zstar]';	
					% transformed positions of data

        xtwiddle = ztwiddle' * e1hat;
        ytwiddle = ztwiddle' * e2hat;
					% x, y values in transformed space

        [adata, rdata] = cart2pol(xtwiddle, ytwiddle);
					% expresses values in polar coords
        hdata = htwiddle;
        disp(['Using recentered data for decomposition']);

      end
      %  #########################################
      %  now cutoff the values outside R0
      %
      vind = find(rdata > 1.05*R0) ;
             %  vector of indices where outside
      rdata(vind) = [] ;
      adata(vind) = [] ;
      hdata(vind) = [] ;
           %  Deletes that part of the vectors


      %  Do least square fitting
      %

      rdata = rdata / R0 ;
          %  normalize to Zernike radius 1
      designm = radial(NMAX,rdata,adata) ;
          %  construct radial Zernike design matrix 
      betahat = (designm' * designm) \ (designm' * hdata) ;
          %  do ordinary least squares fit
      IBASIS  = 2 ;
          %  saves a flag that tells us this is a radially symmetric
          %  decomposition

      outfilestr = ['\matlab\cornean\vimu3' num2str(ifile) '.mat'] ;  
      eval(['save ' outfilestr ' betahat IBASIS IRECENT ' ...
                         'R0 NMAX CRAD CANG']) ;


    end ;


  end ;


elseif ipart == 2 ;    %   then make Zernike Recon. pic.

  figure(3) ;

  ngrid = 64 ;
  nang = 256;    % number of semi-meridians read per radial ring
  R0 = 4.2 ;
  NMAX = 10 ;
  IRECENT = 1 ;    %  work with recentering

  xgrid = linspace(-R0, R0, ngrid) ;
  ygrid = xgrid ;	        %  Vectors of image coordinate values 



  for ifile = 1:10 ;

    disp(['Working on file #' num2str(ifile)]) ;


    if ifile == 1 ;
      titstr = 'Ross,L,304' ;
%      filestr = '000357_2' ;
      filestr = '000352_5' ;
          %  replace filestr with corresponding no lens file
      pupoffamt = 0.119 ;
      pupoffdir = 125 ;
      lencenamt = .99 ;
      lencendir = 50 ;
      lflag = 1 ;
          %  1 for wearing lens
      nlpupoffamt = 0.13 ;
      nlpupoffdir = 101 ;
          %  corresponding no lens pupil center
    elseif ifile == 2 ;
      titstr = 'Ross,L,No' ;
      filestr = '000352_5' ;
      pupoffamt = 0.13 ;
      pupoffdir = 101 ;
      lencenamt = 0 ;
      lencendir = 0 ;
      lflag = 0 ;
          %  1 for wearing lens
    elseif ifile == 3 ;
      titstr = 'Ross,R,No' ;
      filestr = '000352_2' ;
      pupoffamt = 0.248 ;
      pupoffdir = 139 ;
      lencenamt = 0 ;
      lencendir = 0 ;
      lflag = 0 ;
          %  1 for wearing lens
    elseif ifile == 4 ;
      titstr = 'Ross,R,301' ;
%      filestr = '000357_4' ;
      filestr = '000352_2' ;
          %  replace filestr with corresponding no lens file
      pupoffamt = 0.189 ;
      pupoffdir = 144 ;
      lencenamt = 1.3 ;
      lencendir = 118 ;
      lflag = 1 ;
          %  1 for wearing lens
      nlpupoffamt = 0.248 ;
      nlpupoffdir = 139 ;
          %  corresponding no lens pupil center
    elseif ifile == 5 ;
      titstr = 'Duer,L,304' ;
%      filestr = '000358_2' ;
      filestr = '000353_2' ;
          %  replace filestr with corresponding no lens file
      pupoffamt = 0.085 ;
      pupoffdir = 300 ;
      lencenamt = .45 ;
      lencendir = 301 ;
      lflag = 1 ;
          %  1 for wearing lens
      nlpupoffamt = 0.109 ;
      nlpupoffdir = 287 ;
          %  corresponding no lens pupil center
    elseif ifile == 6 ;
      titstr = 'Duer,L,No' ;
      filestr = '000353_2' ;
      pupoffamt = 0.109 ;
      pupoffdir = 287 ;
      lencenamt = 0 ;
      lencendir = 0 ;
      lflag = 0 ;
          %  1 for wearing lens
    elseif ifile == 7 ;
      titstr = 'Duer,R,No' ;
      filestr = '000356_3' ;
      pupoffamt = 0.059 ;
      pupoffdir = 187 ;
      lencenamt = 0 ;
      lencendir = 0 ;
      lflag = 0 ;
          %  1 for wearing lens
    elseif ifile == 8 ;
      titstr = 'Duer,R,301' ;
%      filestr = '000358_1' ;
      filestr = '000356_3' ;
          %  replace filestr with corresponding no lens file
      pupoffamt = 0.099 ;
      pupoffdir = 244 ;
      lencenamt = .89 ;
      lencendir = 145 ;
      lflag = 1 ;
          %  1 for wearing lens
      nlpupoffamt = 0.059 ;
      nlpupoffdir = 187 ;
          %  corresponding no lens pupil center
    elseif ifile == 9 ;
      titstr = 'Duer,L,303' ;
%      filestr = '000360_1' ;
      filestr = '000353_2' ;
          %  replace filestr with corresponding no lens file
      pupoffamt = 0.116 ;
      pupoffdir = 290 ;
      lencenamt = .89 ;
      lencendir = 14 ;
      lflag = 1 ;
          %  1 for wearing lens
      nlpupoffamt = 0.109 ;
      nlpupoffdir = 287 ;
          %  corresponding no lens pupil center
    elseif ifile == 10 ;
      titstr = 'Duer,R,302' ;
%      filestr = '000360_4' ;
      filestr = '000356_3' ;
          %  replace filestr with corresponding no lens file
      pupoffamt = 0.067 ;
      pupoffdir = 261 ;
      lencenamt = .77 ;
      lencendir = 169 ;
      lflag = 1 ;
          %  1 for wearing lens
      nlpupoffamt = 0.059 ;
      nlpupoffdir = 187 ;
          %  corresponding no lens pupil center
    end ;


    if lflag == 1 ;     %  then looking at lens wearing ifile, so proceed




      %  Load previously saved data
      %
      infilestr = ['\matlab\cornean\vimu3' num2str(ifile) '.mat'] ;  
      eval(['load ' infilestr]) ;




      %  Now do Recon. using lines from coazrec.m
      %
      mout = coarr(betahat) ;


      if ifile <= 8 ;
        subplot(3,4,ifile) ;
      else ;
        subplot(3,4,ifile+1) ;
      end ;

      %  Generate map
      %
      pcolor(xgrid,ygrid,mout) ;
        shading flat ;
        axis('square') ;
          caxis([0 1]) ;
        colormap(COMAP) ;

      hold on ;
        plot(0,0,'dk') ;
            %  put a diamond, for lens center
          titstr = ['No Lens ---> ' titstr] ;
        title(titstr) ;
      hold off ;

    end ;

  end ;

    orient landscape ;
  print -dpsc '\matlab\cornean\vimu3a.ps' ;



end ;
