disp('MATLAB Script File: coleg2.m') ;
%
%   For the Cornea Study, this does twiddles with Legendre's
%       including interpretability
%   Copied from coleg2.m, which it modifies by showing
%   several groups on same plot
%
format compact ;

autop = 0 ;   %  1 - for auto print

viimag = [101 110 116] ;  
%  change legend, etc.  when this changes.

              %  1 - Individual 84_1
              %  2 - Individual 85_1
              %  3 - Individual 122_2
              %  4 - Individual 121_2
              %  5 - Individual 55_1
              %  6 - Individual 28_2
              %  7 - Individual 130_1
              %  8 - Individual 139_2
              %  9 - Individual 6_1
              %  10 - Individual 53_1
              %  11 - Individual 43_1

              %  101 - Mean Pre
              %  102 - Mean Pre Right
              %  103 - Mean Pre Left
              %  104 - Mean Diag Ker
              %  105 - Mean Diag Ker Right
              %  106 - Mean Diag Ker Left
              %  107 - Mean Ksuspect
              %  108 - Mean Ksuspect Right
              %  109 - Mean Ksuspect Left
              %  110 - Mean 1 Month
              %  111 - Mean 1 Month Right
              %  112 - Mean 1 Month Left
              %  113 - Mean 6 Month
              %  114 - Mean 6 Month Right
              %  115 - Mean 6 Month Left
              %  116 - Mean 1 Year
              %  117 - Mean 1 Year Right
              %  118 - Mean 1 Year Left
              %  119 - Mean All K
              %  120 - Mean All K Right
              %  121 - Mean All K Left
              %  122 - Mean all P
              %  123 - Mean all P Right
              %  124 - Mean all P Left
              %  125 - Mean Overall
              %  126 - Mean Overall Right
              %  127 - Mean Overall Left


isurf = 1 ;   %  1 - for height map
              %  2 - for height map - sphere  (not working now)
              %  3 - for curvature map


r0 = 4.2 ;    %  Usual radius of analysis

ngrid = 27 ;

rgrid = linspace(0,r0,ngrid)' ;


%  Generate matrix of Legendre's
ell0 = 6 ;
mpolyarg = (2 * rgrid ./ r0) - 1 ;

%  Do ell = 0 Legendre
llleg = ones(size(mpolyarg)) ;
          %  Last, last Legendre for future updating
lllegp = zeros(size(mpolyarg)) ;
          %  Last, last Legendre prime for future updating
lllegpp = zeros(size(mpolyarg)) ;
          %  Last, last Legendre prime, prime for future updating
mleg = llleg ;
mlegp = lllegp ;
mlegpp = lllegpp ;

%  Do ell = 1 Legendre
lleg = mpolyarg ;
          %  Last Legendre for future updating
llegp = ones(size(mpolyarg)) ;
          %  Last Legendre prime for future updating
llegpp = zeros(size(mpolyarg)) ;
          %  Last Legendre prime, prime for future updating
mleg = [mleg lleg] ;
mlegp = [mlegp llegp] ;
mlegpp = [mlegpp llegpp] ;

%  Loop through rest of Legendre's
for ell = 2:ell0 ;
  leg = (1/ell) .* ((2*ell-1) .* mpolyarg .* lleg - (ell-1) .* llleg) ;
              %  Legendre update formula
  legp = (1/ell) .* ((2*ell-1) .* (lleg + mpolyarg .* llegp) ...
                            - (ell-1) .* lllegp) ;
          %  Legendre prime update formula
  legpp = (1/ell) .* ((2*ell-1) .* (2 * llegp + mpolyarg .* llegpp) ...
                            - (ell-1) .* lllegpp) ;
          %  Legendre prime, prime update formula
  mleg = [mleg leg] ;
  mlegp = [mlegp legp] ;
  mlegpp = [mlegpp legpp] ;

  llleg = lleg ;
  lleg = leg ;
  lllegp = llegp ;
  llegp = legp ;
  lllegpp = llegpp ;
  llegpp = legpp ;
          %  Update for next run through loop
end ;

mlegp = mlegp * (2 ./ r0) ;
mlegpp = mlegpp * (2 ./ r0).^2 ;


mout1 = [] ;
mout2 = [] ;
mout3 = [] ;
mout4 = [] ;
mout5 = [] ;
mout6 = [] ;
mout7 = [] ;

%  stuff for loading feature matrices
for ii = 1:length(viimag) ;
  iimag = viimag(ii) ;
  disp(['Working on image no. ' num2str(iimag)]) ;

  if iimag == 1 ;
    imagstr = 'Image 84_1' ;
    fmistr = 'fm420841.mat' ;
    pdist = .123 ;
    pang = 123 ;
    spherad = 7.941 ;  %  Radius of sphere to subtract
    cstep = 6 ;  %  Step size for color scale
  elseif iimag == 2 ;
    imagstr = 'Image 85_1' ;
    fmistr = 'fm420851.mat' ;
    pdist = .254 ;
    pang = 68 ;
    spherad = 7.258 ; %  Radius of sphere to subtract
    cstep = 25 ;  %  Step size for color scale
  elseif iimag == 3 ;
    imagstr = 'Image 122_2' ;
    fmistr = 'fm421222.mat' ;
    pdist = .686 ;
    pang = 89 ;
    spherad = 7.8488 ; %  Radius of sphere to subtract
    cstep = 35 ;  %  Step size for color scale
  elseif iimag == 4 ;
    imagstr = 'Image 121_2' ;
    fmistr = 'fm421212.mat' ;
    pdist = .525 ;
    pang = 108 ;
    spherad = 337.5 / 41 ;  %  Radius of sphere to subtract (Diotpers)
    cstep = 40 ;  %  Step size for color scale
  elseif iimag == 5 ;
    imagstr = 'Image 55_1' ;
    fmistr = 'fm420551.mat' ;
    pdist = .439 ;
    pang = 121 ;
    spherad = 337.5 / 43.5 ;  %  Radius of sphere to subtract (Diotpers)
    cstep = 12 ;  %  Step size for color scale
  elseif iimag == 6 ;
    imagstr = 'Image 28_2' ;
    fmistr = 'fm420282.mat' ;
    pdist = .426 ;
    pang = 115 ;
    spherad = 8.0357 ; %  Radius of sphere to subtract
    cstep = 5 ;  %  Step size for color scale
  elseif iimag == 7 ;
    imagstr = 'Image 130_1' ;
    fmistr = 'fm421301.mat' ;
    pdist = .458 ;
    pang = 112 ;
    spherad = 337.5 / 39.7 ;  %  Radius of sphere to subtract (Diotpers)
    cstep = 10 ;  %  Step size for color scale
  %  spherad = 337.5 / 37.9 ;  %  Radius of sphere to subtract (Diotpers)
  elseif iimag == 8 ;
    imagstr = 'Image 139_2' ;
    fmistr = 'fm421392.mat' ;
    pdist = .417 ;
    pang = 114 ;
    spherad = 8.4375 ; %  Radius of sphere to subtract
    cstep = 10 ;  %  Step size for color scale
  elseif iimag == 9 ;
    imagstr = 'Image 6_1' ;
    fmistr = 'fm420061.mat' ;
    pdist = .454 ;
    pang = 116 ;
    spherad = 7.7586 ; %  Radius of sphere to subtract
    cstep = 6 ;  %  Step size for color scale
  elseif iimag == 10 ;
    imagstr = 'Image 53_1' ;
    fmistr = 'fm420531.mat' ;
    pdist = .375 ;
    pang = 113 ;
    spherad = 8.766 ; %  Radius of sphere to subtract
    cstep = 10 ;  %  Step size for color scale
  elseif iimag == 11 ;
    imagstr = 'Image 43_1' ;
    fmistr = 'fm420431.mat' ;
    pdist = .557 ;
    pang = 91 ;
    spherad = 8.5443 ; %  Radius of sphere to subtract
    cstep = 10 ;  %  Step size for color scale


  elseif iimag == 101 ;
    imagstr = 'Mean Pre' ;
    fmistr = 'mn42pre.mat' ;
    spherad = 8 ;  %  Radius of sphere to subtract
    cstep = 10 ;  %  Step size for color scale
  elseif iimag == 102 ;
    imagstr = 'Mean Pre Right' ;
    fmistr = 'mn42pre.mat' ;
    spherad = 8 ;  %  Radius of sphere to subtract
    cstep = 10 ;  %  Step size for color scale
  elseif iimag == 103 ;
    imagstr = 'Mean Pre Left' ;
    fmistr = 'mn42pre.mat' ;
    spherad = 8 ;  %  Radius of sphere to subtract
    cstep = 10 ;  %  Step size for color scale
  elseif iimag == 104 ;
    imagstr = 'Mean Diag Ker' ;
            %  Diagnosed Keratoconus 
            %          - called ker by NKT
            %          - called kdi by JTZ
    fmistr = 'mn42kdi.mat' ;
    spherad = 8 ;  %  Radius of sphere to subtract
    cstep = 10 ;  %  Step size for color scale
  elseif iimag == 105 ;
    imagstr = 'Mean Diag Ker Right' ;
    fmistr = 'mn42kdi.mat' ;
    spherad = 8 ;  %  Radius of sphere to subtract
    cstep = 10 ;  %  Step size for color scale
  elseif iimag == 106 ;
    imagstr = 'Mean Diag Ker Left' ;
    fmistr = 'mn42kdi.mat' ;
    spherad = 8 ;  %  Radius of sphere to subtract
    cstep = 10 ;  %  Step size for color scale
  elseif iimag == 107 ;
    imagstr = 'Mean Ksuspect' ;
    fmistr = 'mn42ksu.mat' ;
    spherad = 8 ;  %  Radius of sphere to subtract
    cstep = 10 ;  %  Step size for color scale
  elseif iimag == 108 ;
    imagstr = 'Mean Ksuspect Right' ;
    fmistr = 'mn42ksu.mat' ;
    spherad = 8 ;  %  Radius of sphere to subtract
    cstep = 10 ;  %  Step size for color scale
  elseif iimag == 109 ;
    imagstr = 'Mean Ksuspect Left' ;
    fmistr = 'mn42ksu.mat' ;
    spherad = 8 ;  %  Radius of sphere to subtract
    cstep = 10 ;  %  Step size for color scale
  elseif iimag == 110 ;
    imagstr = 'Mean 1 Month' ;
    fmistr = 'mn421mo.mat' ;
    spherad = 8 ;  %  Radius of sphere to subtract
    cstep = 10 ;  %  Step size for color scale
  elseif iimag == 111 ;
    imagstr = 'Mean 1 Month Right' ;
    fmistr = 'mn421mo.mat' ;
    spherad = 8 ;  %  Radius of sphere to subtract
    cstep = 10 ;  %  Step size for color scale
  elseif iimag == 112 ;
    imagstr = 'Mean 1 Month Left' ;
    fmistr = 'mn421mo.mat' ;
    spherad = 8 ;  %  Radius of sphere to subtract
    cstep = 10 ;  %  Step size for color scale
  elseif iimag == 113 ;
    imagstr = 'Mean 6 Month' ;
    fmistr = 'mn426mo.mat' ;
    spherad = 8 ;  %  Radius of sphere to subtract
    cstep = 10 ;  %  Step size for color scale
  elseif iimag == 114 ;
    imagstr = 'Mean 6 Month Right' ;
    fmistr = 'mn426mo.mat' ;
    spherad = 8 ;  %  Radius of sphere to subtract
    cstep = 10 ;  %  Step size for color scale
  elseif iimag == 115 ;
    imagstr = 'Mean 6 Month Left' ;
    fmistr = 'mn426mo.mat' ;
    spherad = 8 ;  %  Radius of sphere to subtract
    cstep = 10 ;  %  Step size for color scale
  elseif iimag == 116 ;
    imagstr = 'Mean 1 Year' ;
    fmistr = 'mn421yr.mat' ;
    spherad = 8 ;  %  Radius of sphere to subtract
    cstep = 10 ;  %  Step size for color scale
  elseif iimag == 117 ;
    imagstr = 'Mean 1 Year Right' ;
    fmistr = 'mn421yr.mat' ;
    spherad = 8 ;  %  Radius of sphere to subtract
    cstep = 10 ;  %  Step size for color scale
  elseif iimag == 118 ;
    imagstr = 'Mean 1 Year Left' ;
    fmistr = 'mn421yr.mat' ;
    spherad = 8 ;  %  Radius of sphere to subtract
    cstep = 10 ;  %  Step size for color scale
  elseif iimag == 119 ;
    imagstr = 'Mean All K Study' ;
            %  All 3 groups in Keratoconus study
            %          - called kstud by NKT
            %          - called ker by JTZ
    fmistr = 'mn42ker.mat' ;
    spherad = 8 ;  %  Radius of sphere to subtract
    cstep = 10 ;  %  Step size for color scale
  elseif iimag == 120 ;
    imagstr = 'Mean All K Study Right' ;
    fmistr = 'mn42ker.mat' ;
    spherad = 8 ;  %  Radius of sphere to subtract
    cstep = 10 ;  %  Step size for color scale
  elseif iimag == 121 ;
    imagstr = 'Mean All K Study Left' ;
    fmistr = 'mn42ker.mat' ;
    spherad = 8 ;  %  Radius of sphere to subtract
    cstep = 10 ;  %  Step size for color scale
  elseif iimag == 122 ;
    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).
    fmistr = 'mn42pos.mat' ;
    spherad = 8 ;  %  Radius of sphere to subtract
    cstep = 10 ;  %  Step size for color scale
  elseif iimag == 123 ;
    imagstr = 'Mean All P Study Right' ;
    fmistr = 'mn42pos.mat' ;
    spherad = 8 ;  %  Radius of sphere to subtract
    cstep = 10 ;  %  Step size for color scale
  elseif iimag == 124 ;
    imagstr = 'Mean All P Study Left' ;
    fmistr = 'mn42pos.mat' ;
    spherad = 8 ;  %  Radius of sphere to subtract
    cstep = 10 ;  %  Step size for color scale
  elseif iimag == 125 ;
    imagstr = 'Overall Mean' ;
            %  Note:  pre group is included twice as "pre" and as "nor"
    fmistr = 'mn42all.mat' ;
    spherad = 8 ;  %  Radius of sphere to subtract
    cstep = 10 ;  %  Step size for color scale
  elseif iimag == 126 ;
    imagstr = 'Overall Mean Right' ;
    fmistr = 'mn42all.mat' ;
    spherad = 8 ;  %  Radius of sphere to subtract
    cstep = 10 ;  %  Step size for color scale
  elseif iimag == 127 ;
    imagstr = 'Overall Mean Left' ;
    fmistr = 'mn42all.mat' ;
    spherad = 8 ;  %  Radius of sphere to subtract
    cstep = 10 ;  %  Step size for color scale
  end ;
  
  ['load \matlab\steve\cornea\fm\' fmistr]
  eval(['load \matlab\steve\cornea\fm\' fmistr]) ;
  
  if iimag < 100 ;
    mfeat = fm' ;    %  Jinting used name fm
  else ;
    if rem(iimag,3) == 2 ;   %  Then this is an "all" image
      mamp = mnamp ;
      mpha = mnpha ;
    elseif rem(iimag,3) == 0 ;   %  Then this is a "right" image
      mamp = mnampr ;
      mpha = mnphar ;
    elseif rem(iimag,3) == 1 ;   %  Then this is a "left" image
      mamp = mnampl ;
      mpha = mnphal ;
    end ;
  
    fm = mamp(1,:) .* cos(mpha(1,:) * pi / 180) ;
            %  Here Jinting used names mnamp & mnpha
    for i = 2:10 ;
      fm = [fm; mamp(i,:) .* cos(mpha(i,:) * pi / 180)] ;
      fm = [fm; mamp(i,:) .* sin(mpha(i,:) * pi / 180)] ;
    end ;
    mfeat = fm' ;
  end ;



  sphht = spherad - sqrt(spherad.^2 - rgrid.^2) ;
  %plot(rgrid,sphht) ;


  legcoef = mfeat(:,1) ;


  %  reconstruct Legendre's 
  mrecon = mleg' .* vec2mat(legcoef,ngrid) ;
  mrecon = cumsum(mrecon)' ;
  mreconp = mlegp' .* vec2mat(legcoef,ngrid) ;
  mreconp = cumsum(mreconp)' ;
  mreconpp = mlegpp' .* vec2mat(legcoef,ngrid) ;
  mreconpp = cumsum(mreconpp)' ;
  
  mcurv = -mreconpp ./ (1 + mreconp.^2).^(3/2) ;
  mcurv = mcurv * 337.5 ;
          %  put on scale of diopters

  if isurf == 1 ;  %  Heights
    surfstr = 'Heights' ;
    mout1 = [mout1 mrecon(:,length(legcoef))] ;
    mout2 = [mout2 mrecon(:,1)] ;
    mout3 = [mout3 mrecon(:,2)-mrecon(:,1)] ;
    mout4 = [mout4 mrecon(:,3)-mrecon(:,2)] ;
    mout5 = [mout5 mrecon(:,4)-mrecon(:,3)] ;
    mout6 = [mout6 mrecon(:,5)-mrecon(:,4)] ;
    mout7 = [mout7 mrecon(:,6)-mrecon(:,5)] ;
  elseif isurf == 2 ;  %  Curvature
    surfstr = 'Curvature' ;
    mout1 = [mout1 mcurv(:,length(legcoef))] ;
    mout2 = [mout2 mcurv(:,3)-mcurv(:,2)] ;
  end ;

end ;

subplot(1,1,1) ;
  plot(rgrid,mout1(:,1),'-',rgrid,mout1(:,2),':',rgrid,mout1(:,3),'-.') ;
    title(['Four 0, Leg 0-6,  ' surfstr]) ;
    legend('Pre','1mo','1yr') ;

if autop ~= 0 ;
  eval(['print -dps le3' num2str(iimag) 'a.ps']) ;
end ;


ommp = mout1(:,2) - mout1(:,1) ;
oymom = mout1(:,3) - mout1(:,2) ;
subplot(1,1,1) ;
  plot(rgrid,ommp,'-',rgrid,oymom,':') ;
    title(['Four 0, Leg 0-6,  ' surfstr]) ;
    legend('1mo - Pre','1yr - 1mo') ;

if autop ~= 0 ;
  eval(['print -dps le3' num2str(iimag) 'b.ps']) ;
end ;


ommp2 = mout2(:,2) - mout2(:,1) ;
oymom2 = mout2(:,3) - mout2(:,2) ;
ommp3 = mout3(:,2) - mout3(:,1) ;
oymom3 = mout3(:,3) - mout3(:,2) ;
ommp4 = mout4(:,2) - mout4(:,1) ;
oymom4 = mout4(:,3) - mout4(:,2) ;
ommp5 = mout5(:,2) - mout5(:,1) ;
oymom5 = mout5(:,3) - mout5(:,2) ;
ommp6 = mout6(:,2) - mout6(:,1) ;
oymom6 = mout6(:,3) - mout6(:,2) ;
ommp7 = mout7(:,2) - mout7(:,1) ;
oymom7 = mout7(:,3) - mout7(:,2) ;

subplot(3,2,1) ;
  plot(rgrid,ommp2,'-',rgrid,oymom2,':') ;
    title(['Four 0, Leg 0,  ' surfstr]) ;
    legend('1mo - Pre','1yr - 1mo') ;

subplot(3,2,2) ;
  plot(rgrid,ommp3,'-',rgrid,oymom3,':') ;
    title(['Four 0, Leg 1,  ' surfstr]) ;
    legend('1mo - Pre','1yr - 1mo') ;
    axis([0 4.2 -.006 .006]) ;

subplot(3,2,3) ;
  plot(rgrid,ommp4,'-',rgrid,oymom4,':') ;
    title(['Four 0, Leg 2,  ' surfstr]) ;
    legend('1mo - Pre','1yr - 1mo') ;
    axis([0 4.2 -.006 .006]) ;

subplot(3,2,4) ;
  plot(rgrid,ommp5,'-',rgrid,oymom5,':') ;
    title(['Four 0, Leg 3,  ' surfstr]) ;
    legend('1mo - Pre','1yr - 1mo') ;
    axis([0 4.2 -.006 .006]) ;

subplot(3,2,5) ;
  plot(rgrid,ommp6,'-',rgrid,oymom6,':') ;
    title(['Four 0, Leg 4,  ' surfstr]) ;
    legend('1mo - Pre','1yr - 1mo') ;
    axis([0 4.2 -.006 .006]) ;

subplot(3,2,6) ;
  plot(rgrid,ommp7,'-',rgrid,oymom7,':') ;
    title(['Four 0, Leg 5,  ' surfstr]) ;
    legend('1mo - Pre','1yr - 1mo') ;
    axis([0 4.2 -.006 .006]) ;

orient tall ;

if autop ~= 0 ;
  eval(['print -dps le3' num2str(iimag) 'c.ps']) ;
end ;
