%MATLAB Script File: corecon5.m
%
%   For the Cornea Study, this reconstructs surfaces from
%                               mean feature matrices
%       Modification of corecon4.m, which it modifes
%           by using correctly averaged feature matrices
%           (after angle problem was fixed)
%
format compact ;
disp('MATLAB Script File:  corecon5.m') ;

autop = 1 ;   %  1 to create a color postscript file

ihires = 0 ;  %  Controls resolution of plot:
              %    0 - for lo res 64 x64  (fastest)
              %    1 - to interpolate to hi res 256 x 256  (slower)
              %    2 - for hi res 256 x 256  (REALLY slow)
              

iimag = 12 ;  
              %  1 - Mean Pre
              %  2 - Mean Pre Right
              %  3 - Mean Pre Left
              %  4 - Mean Diag Ker
              %  5 - Mean Diag Ker Right
              %  6 - Mean Diag Ker Left
              %  7 - Mean Ksuspect
              %  8 - Mean Ksuspect Right
              %  9 - Mean Ksuspect Left
              %  10 - Mean 1 Month
              %  11 - Mean 1 Month Right
              %  12 - Mean 1 Month Left
              %  13 - Mean 6 Month
              %  14 - Mean 6 Month Right
              %  15 - Mean 6 Month Left
              %  16 - Mean 1 Year
              %  17 - Mean 1 Year Right
              %  18 - Mean 1 Year Left
              %  19 - Mean All K
              %  20 - Mean All K Right
              %  21 - Mean All K Left
              %  22 - Mean all P
              %  23 - Mean all P Right
              %  24 - Mean all P Left
              %  25 - Mean Overall
              %  26 - Mean Overall Right
              %  27 - Mean Overall Left

isurf = 3 ;   %  1 - for height map
              %  2 - for height map - sphere
              %  3 - for curvature map


ibas = 1 ;    %  1 - Radial Fourier - Legendre
              %  2 - Diameter Fourier - Legendre
              %            (not working yet)
              %  3 - Zernike
              %            (not working yet)

figure(2) ;
%clg reset ;
%disp('Hit any key to continue after clearing graphics') ;
%pause ;

r0 = 4.2 ;    %  Usual radius of analysis
%spherad = 9.3 ;%  Radius of sphere to subtract
%  now set individually for each image

ngrid = 64 ; %  Construct surfaces as ngrid x ngrid array
if ihires == 2 ;
          %  then change ngrid to larger value
  ngrid = 256 ; %  Construct surfaces as ngrid x ngrid array
end ;
hrngrid = 256 ; %  Interpolate to hrngrid x hrngrid for 
                %  higher resolution part

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

[mxgrid, mygrid] = meshgrid(xgrid,ygrid) ;
              %  Matrices of image coordinate values 
              
mrgrid = sqrt(mxgrid.^2 + mygrid.^2) ;
              %  Matrices of image polar coordinate radii
mthetag = atan2(mxgrid,mygrid) ;
              %  Matrices of image polar coordinate angles (-pi,pi]
  vind = find(mthetag < 0) ;
mthetag(vind) = mthetag(vind) + 2 * pi ;
              %  Matrices of image polar coordinate angles [0,2pi)
          
mcrop = ones(ngrid,ngrid) ;            
  vind = find(mrgrid > r0) ;
              %  Vector of indices where image is outside circle
mcrop(vind) = nan * ones(length(vind),1) ;
              %  Cropping Matrix, multiply by this to cut off images
              %  outside circle of radius r_0.       


nr = 10 ;
nc = 7 ;

if iimag == 1 ;
  imagstr = 'Mean Pre' ;
  load \matlab\steve\cornea\ht42mnpr.mat ;
    fmamp = mnamp ;
    fmpha = mnpha ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 2 ;
  imagstr = 'Mean Pre Right' ;
  load \matlab\steve\cornea\ht42mnpr.mat ;
    fmamp = mnampr ;
    fmpha = mnphar ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 3 ;
  imagstr = 'Mean Pre Left' ;
  load \matlab\steve\cornea\ht42mnpr.mat ;
    fmamp = mnampl ;
    fmpha = mnphal ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 4 ;
  imagstr = 'Mean Diag Ker' ;
          %  Diagnosed Keratoconus 
          %          - called ker by NKT
          %          - called kdi by JTZ
  load \matlab\steve\cornea\ht42mnkd.mat ;
    fmamp = mnamp ;
    fmpha = mnpha ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 5 ;
  imagstr = 'Mean Diag Ker Right' ;
  load \matlab\steve\cornea\ht42mnkd.mat ;
    fmamp = mnampr ;
    fmpha = mnphar ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 6 ;
  imagstr = 'Mean Diag Ker Left' ;
  load \matlab\steve\cornea\ht42mnkd.mat ;
    fmamp = mnampl ;
    fmpha = mnphal ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 7 ;
  imagstr = 'Mean Ksuspect' ;
  load \matlab\steve\cornea\ht42mnku.mat ;
    fmamp = mnamp ;
    fmpha = mnpha ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 8 ;
  imagstr = 'Mean Ksuspect Right' ;
  load \matlab\steve\cornea\ht42mnku.mat ;
    fmamp = mnampr ;
    fmpha = mnphar ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 9 ;
  imagstr = 'Mean Ksuspect Left' ;
  load \matlab\steve\cornea\ht42mnku.mat ;
    fmamp = mnampl ;
    fmpha = mnphal ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 10 ;
  imagstr = 'Mean 1 Month' ;
  load \matlab\steve\cornea\ht42mn1m.mat ;
    fmamp = mnamp ;
    fmpha = mnpha ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 11 ;
  imagstr = 'Mean 1 Month Right' ;
  load \matlab\steve\cornea\ht42mn1m.mat ;
    fmamp = mnampr ;
    fmpha = mnphar ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 12 ;
  imagstr = 'Mean 1 Month Left' ;
  load \matlab\steve\cornea\ht42mn1m.mat ;
    fmamp = mnampl ;
    fmpha = mnphal ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 13 ;
  imagstr = 'Mean 6 Month' ;
  load \matlab\steve\cornea\ht42mn6m.mat ;
    fmamp = mnamp ;
    fmpha = mnpha ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 14 ;
  imagstr = 'Mean 6 Month Right' ;
  load \matlab\steve\cornea\ht42mn6m.mat ;
    fmamp = mnampr ;
    fmpha = mnphar ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 15 ;
  imagstr = 'Mean 6 Month Left' ;
  load \matlab\steve\cornea\ht42mn6m.mat ;
    fmamp = mnampl ;
    fmpha = mnphal ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 16 ;
  imagstr = 'Mean 1 Year' ;
  load \matlab\steve\cornea\ht42mn1y.mat ;
    fmamp = mnamp ;
    fmpha = mnpha ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 17 ;
  imagstr = 'Mean 1 Year Right' ;
  load \matlab\steve\cornea\ht42mn1y.mat ;
    fmamp = mnampr ;
    fmpha = mnphar ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 18 ;
  imagstr = 'Mean 1 Year Left' ;
  load \matlab\steve\cornea\ht42mn1y.mat ;
    fmamp = mnampl ;
    fmpha = mnphal ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 19 ;
  imagstr = 'Mean All K Study' ;
          %  All 3 groups in Keratoconus study
          %          - called kstud by NKT
          %          - called ker by JTZ
  load \matlab\steve\cornea\ht42mnke.mat ;
    fmamp = mnamp ;
    fmpha = mnpha ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 20 ;
  imagstr = 'Mean All K Study Right' ;
  load \matlab\steve\cornea\ht42mnke.mat ;
    fmamp = mnampr ;
    fmpha = mnphar ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 21 ;
  imagstr = 'Mean All K Study Left' ;
  load \matlab\steve\cornea\ht42mnke.mat ;
    fmamp = mnampl ;
    fmpha = mnphal ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 22 ;
  imagstr = 'Mean All P Study' ;
          %  All 4 groups in PRK study
          %          - called pstud by NKT
          %          - called pos by JTZ
          %                   nor by JTZ are not used (same as pre).
  load \matlab\steve\cornea\ht42mnpo.mat ;
    fmamp = mnamp ;
    fmpha = mnpha ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 23 ;
  imagstr = 'Mean All P Study Right' ;
  load \matlab\steve\cornea\ht42mnpo.mat ;
    fmamp = mnampr ;
    fmpha = mnphar ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 24 ;
  imagstr = 'Mean All P Study Left' ;
  load \matlab\steve\cornea\ht42mnpo.mat ;
    fmamp = mnampl ;
    fmpha = mnphal ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 25 ;
  imagstr = 'Overall Mean' ;
          %  Note:  pre group is included twice as "pre" and as "nor"
  load \matlab\steve\cornea\ht42mnal.mat ;
    fmamp = mnamp ;
    fmpha = mnpha ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 26 ;
  imagstr = 'Overall Mean Right' ;
  load \matlab\steve\cornea\ht42mnal.mat ;
    fmamp = mnampr ;
    fmpha = mnphar ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
elseif iimag == 27 ;
  imagstr = 'Overall Mean Left' ;
  load \matlab\steve\cornea\ht42mnal.mat ;
    fmamp = mnampl ;
    fmpha = mnphal ;
  spherad = 8 ;  %  Radius of sphere to subtract
  cstep = 10 ;  %  Step size for color scale
end ;





fmpha = 2 * pi * fmpha / 360 ;
                %  converts degrees to radians

fmamp = fmamp' ;
fmpha = fmpha' ;
           %  puts in same form as above

    %  Convert to c-s form of feature matrix
mfeat = fmamp(:,1) .* cos(fmpha(:,1)) ;
for icol = 2:size(fmamp,2) ;   %  Build up c-s feature matrix
  mfeat = [mfeat (fmamp(:,icol) .* cos(fmpha(:,icol)))] ;
  mfeat = [mfeat (fmamp(:,icol) .* sin(fmpha(:,icol)))] ;
end ;

disp(['feature matrix has size ' num2str(size(mfeat,1)) ' rows ']) ;
disp(['and ' num2str(size(mfeat,2)) ' cols']) ;




disp('Generating surface') ; 
if isurf == 1 | isurf == 2 ;   %  then work with height map
  surfstr = ' Heights' ;

  if ibas == 1 ;     %  then use the radial Fourier - Legendre basis

    mout = corh(mfeat,mrgrid,mthetag,r0) ;
          %  Radial Fourier - Legendre reconstruction of Heights

  elseif ibas == 2 ;     %  then use the diameter Fourier - Legendre basis

  elseif ibas == 3 ;     %  then use the Zernike basis

  end ;

  if isurf == 2 ;   %  Then need to subtract sphere
    surfstr = [surfstr ' - Sphere'] ;

%    minmout = min(min(mout))
%    maxmout = max(max(mout))
            %  use these lines to decide which sphere to subtract

    sphht = spherad - sqrt(spherad.^2 - mrgrid.^2) ;

%    minsph = min(min(sphht))
%    maxsph = max(max(sphht))
%    pause
            %  use these lines to decide which sphere to subtract


    mout = mout + sphht ;
            %  careful about whether this is + or -
  end ;

elseif isurf == 3 ;   %  then work with curvature map
  surfstr = ' Curvatures' ;

  if ibas == 1 ;     %  then use the radial Fourier - Legendre basis

    [mh,mhp,mhpp] = corhpp(mfeat,mrgrid,mthetag,r0) ;

  elseif ibas == 2 ;     %  then use the diameter Fourier - Legendre basis

  elseif ibas == 3 ;     %  then use the Zernike basis
  end ;

  mout = -mhpp ./ (1 + mhp.^2).^(3/2) ;
  mout = mout * 337.5 ;
          %  put on scale of diopters

end ;





if ihires == 1 ;
          %  Then compute higher resolution stuff
  hrxg = linspace(-r0,r0,hrngrid) ;
  hryg = hrxg ;
  [mhrxg,mhryg] = meshgrid(hrxg,hryg) ;
  mhrout = interp2(mxgrid,mygrid,mout,mhrxg,mhryg) ;
  mcrop = ones(hrngrid) ;            
  vind = find(mhrxg.^2 + mhryg.^2 > r0^2) ;
          %  Vector of indices where image is outside circle
  mcrop(vind) = nan * ones(length(vind),1) ;
          %  Cropping Matrix, multiply by this to cut off images
          %  outside circle of radius r_0.       
  mhrout = mhrout .* mcrop ;
end ;






if isurf == 2 ;
          %  Then setup Keratron height color map

  %  These use the "full range" Keratron data, from the email message:
  %             cocolora.msg
  %   which were copied into the ASCII file 
  %             cocolor1.dat
  %  for importing into matlab programs
  %    locurve & hicurve are upper and lower ends of the "range"
  %    for each color
  load cornea\cocolor1.dat ;
  cyan = cocolor1(:,3) ;
  magenta = cocolor1(:,4) ;
  yellow = cocolor1(:,5) ;
  black = cocolor1(:,6) ;
  %  Don't use columns 1 and 2, because those are for curvature

  red   = 1 - cyan + black ;
  green = 1 - magenta + black;
  blue  = 1 - yellow + black ;

  comap = [red green blue] ;
  ncomap = size(comap,1) ;
          %  number of rows in the color map ;

  %  make sure values are positive
  %     (since twiddling "black" can leave this range)
  flag = comap < 0 ;
  fflag = find(flag) ;
  comap(fflag) = zeros(length(fflag),1) ;

  %  make sure values are less than 1
  %     (since twiddling "black" can leave this range)
  flag = comap > 1 ;
  fflag = find(flag) ;
  comap(fflag) = ones(length(fflag),1) ;

  colormap(comap) ;

      cstep = cstep / 1000 ;
          %  cstep is "step size in color scale", different for
          %  each image

      u01 = linspace(0,1,ncomap+1) ;
        lbd = -11.5 * cstep ;
        ubd = lbd + ncomap * cstep ;
      colgrid = linspace(lbd,ubd,ncomap+1) ;
          %  equally spaced grid on desired color scale

    vmout = mout(:) ;
          %  vector version for easier manipulation
          %  NOTE: interp1 seems to need vector argument
      indlo = find(vmout < lbd) ;
    vmout(indlo) = lbd * ones(length(indlo),1) ;
      indup = find(vmout > ubd) ;
    vmout(indup) = lbd * ones(length(indup),1) ;
          %  NOTE: interp1 seems to need values within grid

    colcurv = interp1(colgrid, u01, vmout) ;
          %  use this to determine vector of colors representing
          %  curvature, when plotting.  This transforms curvature 
          %  to get Keratron colors, when using comap as colormap.
  mout = reshape(colcurv,size(mout)) ;  

elseif isurf == 3 ;
          %  Then setup Keratron color map

  %  These use the "full range" Keratron data, from the email message:
  %             cocolora.msg
  %   which were copied into the ASCII file 
  %             cocolor1.dat
  %  for importing into matlab programs
  %    locurve & hicurve are upper and lower ends of the "range"
  %    for each color
  load cornea\cocolor1.dat ;
  locurve = cocolor1(:,1) ;
  hicurve = cocolor1(:,2) ;
  cyan = cocolor1(:,3) ;
  magenta = cocolor1(:,4) ;
  yellow = cocolor1(:,5) ;
  black = cocolor1(:,6) ;

  red   = 1 - cyan + black ;
  green = 1 - magenta + black;
  blue  = 1 - yellow + black ;

  comap = [red green blue] ;
  ncomap = size(comap,1) ;
          %  number of rows in the color map ;

  %  make sure values are positive
  flag = comap < 0 ;
  fflag = find(flag) ;
  comap(fflag) = zeros(length(fflag),1) ;

  %  make sure values are less than 1
  flag = comap > 1 ;
  fflag = find(flag) ;
  comap(fflag) = ones(length(fflag),1) ;

    vmout = mout(:) ;
          %  vector version for easier manipulation
          %  NOTE: interp1 seems to need vector argument
      indlo = find(vmout < 0) ;
    vmout(indlo) = 0 * ones(length(indlo),1) ;
      indup = find(vmout > max(hicurve)) ;
    vmout(indup) = max(hicurve) * ones(length(indup),1) ;
          %  NOTE: interp1 seems to need values within grid

  colormap(comap) ;

    u01 = linspace(0,1,ncomap+1) ;
    colcurv = interp1([0; hicurve], u01, vmout) ;
          %  use this to determine vector of colors representing
          %  curvature, when plotting.  This transforms curvature 
          %  to get Keratron colors, when using comap as colormap.
          %  NOTE: interp1 seems to need vector argument
  mout = reshape(colcurv,size(mout)) ;  
end ;




if ihires == 1 ;
          %  Then compute higher resolution stuff
  hrxg = linspace(-r0,r0,hrngrid) ;
  hryg = hrxg ;
  [mhrxg,mhryg] = meshgrid(hrxg,hryg) ;
  mhrout = interp2(mxgrid,mygrid,mout,mhrxg,mhryg) ;
  mcrop = ones(hrngrid) ;            
  vind = find(mhrxg.^2 + mhryg.^2 > r0^2) ;
          %  Vector of indices where image is outside circle
  mcrop(vind) = nan * ones(length(vind),1) ;
          %  Cropping Matrix, multiply by this to cut off images
          %  outside circle of radius r_0.       
  mhrout = mhrout .* mcrop ;
else ;
  mout = mout .* mcrop ;
          %  Cut off part outside analysis circle
end ;



disp('Making Plot') ;

if ihires == 1 ;
  pcolor(mhrxg,mhryg,mhrout) ;
else ;
  pcolor(mxgrid,mygrid,mout) ;
end ;
  shading flat ;
  axis('square') ;
  title([imagstr surfstr]) ;
  if isurf == 2 | isurf == 3 ;   
          %  if heightmap - sphere or curvature
    caxis([0 1]) ;
  end ;


if autop == 1 ;
  if isurf == 2 ;
    if iimag == 1 ;
      print corec1h.ps -dpsc ;
    elseif iimag == 2 ;
      print corec2h.ps -dpsc ;
    elseif iimag == 3 ;
      print corec3h.ps -dpsc ;
    elseif iimag == 4 ;
      print corec4h.ps -dpsc ;
    elseif iimag == 5 ;
      print corec5h.ps -dpsc ;
    elseif iimag == 6 ;
      print corec6h.ps -dpsc ;
    elseif iimag == 7 ;
      print corec7h.ps -dpsc ;
    elseif iimag == 8 ;
      print corec8h.ps -dpsc ;
    elseif iimag == 9 ;
      print corec9h.ps -dpsc ;
    elseif iimag == 10 ;
      print corec10h.ps -dpsc ;
    elseif iimag == 11 ;
      print corec11h.ps -dpsc ;
    elseif iimag == 12 ;
      print corec12h.ps -dpsc ;
    elseif iimag == 13 ;
      print corec13h.ps -dpsc ;
    elseif iimag == 14 ;
      print corec14h.ps -dpsc ;
    elseif iimag == 15 ;
      print corec15h.ps -dpsc ;
    elseif iimag == 16 ;
      print corec16h.ps -dpsc ;
    elseif iimag == 17 ;
      print corec17h.ps -dpsc ;
    elseif iimag == 18 ;
      print corec18h.ps -dpsc ;
    elseif iimag == 19 ;
      print corec19h.ps -dpsc ;
    elseif iimag == 20 ;
      print corec20h.ps -dpsc ;
    elseif iimag == 21 ;
      print corec21h.ps -dpsc ;
    elseif iimag == 22 ;
      print corec22h.ps -dpsc ;
    elseif iimag == 23 ;
      print corec23h.ps -dpsc ;
    elseif iimag == 24 ;
      print corec24h.ps -dpsc ;
    elseif iimag == 25 ;
      print corec25h.ps -dpsc ;
    elseif iimag == 26 ;
      print corec26h.ps -dpsc ;
    elseif iimag == 27 ;
      print corec27h.ps -dpsc ;
    end ;
  elseif isurf == 3 ;
    if iimag == 1 ;
      print corec1c.ps -dpsc ;
    elseif iimag == 2 ;
      print corec2c.ps -dpsc ;
    elseif iimag == 3 ;
      print corec3c.ps -dpsc ;
    elseif iimag == 4 ;
      print corec4c.ps -dpsc ;
    elseif iimag == 5 ;
      print corec5c.ps -dpsc ;
    elseif iimag == 6 ;
      print corec6c.ps -dpsc ;
    elseif iimag == 7 ;
      print corec7c.ps -dpsc ;
    elseif iimag == 8 ;
      print corec8c.ps -dpsc ;
    elseif iimag == 9 ;
      print corec9c.ps -dpsc ;
    elseif iimag == 10 ;
      print corec10c.ps -dpsc ;
    elseif iimag == 11 ;
      print corec11c.ps -dpsc ;
    elseif iimag == 12 ;
      print corec12c.ps -dpsc ;
    elseif iimag == 13 ;
      print corec13c.ps -dpsc ;
    elseif iimag == 14 ;
      print corec14c.ps -dpsc ;
    elseif iimag == 15 ;
      print corec15c.ps -dpsc ;
    elseif iimag == 16 ;
      print corec16c.ps -dpsc ;
    elseif iimag == 17 ;
      print corec17c.ps -dpsc ;
    elseif iimag == 18 ;
      print corec18c.ps -dpsc ;
    elseif iimag == 19 ;
      print corec19c.ps -dpsc ;
    elseif iimag == 20 ;
      print corec20c.ps -dpsc ;
    elseif iimag == 21 ;
      print corec21c.ps -dpsc ;
    elseif iimag == 22 ;
      print corec22c.ps -dpsc ;
    elseif iimag == 23 ;
      print corec23c.ps -dpsc ;
    elseif iimag == 24 ;
      print corec24c.ps -dpsc ;
    elseif iimag == 25 ;
      print corec25c.ps -dpsc ;
    elseif iimag == 26 ;
      print corec26c.ps -dpsc ;
    elseif iimag == 27 ;
      print corec27c.ps -dpsc ;
    end ;
  end ;
end ;


