disp('MATLAB Script File: costst1.m') ;
%
%   For the Cornea Study, this does twiddles with Legendre's
%       including interpretability
%   Copied from costst1.m                            
%
format compact ;

ipart = 1 ;   %  1

igroup = 6 ;  %  1 - Pre
              %  2 - 1 month
              %  3 - 6 month
              %  4 - 1 year
              %  5 - Individ 006_1
              %  6 - Individ 028_2

spherad = 9.6 ;


r0 = 4.2 ;    %  Usual radius of analysis

ngrid = 27 ;

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

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

%plot(rgrid,sphht) ;

%  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 ;




if ipart == 1 ;
  if igroup == 1 ;    % Pre
    legcoef = [-4.0504e-01; -6.1574e-01; -2.1762e-01; -7.9496e-03; ...
                   -9.1958e-04; 2.1232e-04; 1.0224e-04] ;
    titstr = 'Pre' ;
  elseif igroup == 2 ;    % 1 month
    legcoef = [-3.7069e-01; -5.7139e-01; -2.1670e-01; -1.9918e-02; ...
                   -2.6055e-03; 2.2705e-03; 9.5838e-04] ;
    titstr = '1 Mo' ;
  elseif igroup == 3 ;    % 6 month
    legcoef = [-3.7237e-01; -5.7312e-01; -2.1570e-01; -1.8871e-02; ...
                   -3.2538e-03; 1.5447e-03; 8.7718e-04] ;
    titstr = '6 Mo' ;
  elseif igroup == 4 ;    % 1 year
    legcoef = [-3.6907e-01; -5.6958e-01; -2.1585e-01; -1.8527e-02; ...
                   -2.2716e-03; 1.7264e-03; 8.1738e-04] ;
    titstr = '1 Yr' ;
  elseif igroup == 5 ;    % Individ 6_1
    legcoef = [-4.0706e-01; -6.1800e-01; -2.1729e-01; -7.2570e-03; ...
                   -8.6380e-04; 2.0735e-04; 1.5811e-04] ;
    titstr = 'Individ 6_1' ;
  elseif igroup == 6 ;    % Individ 28_2
    legcoef = [-3.8558e-01; -5.8558e-01; -2.0641e-01; -7.4671e-03; ...
                   -9.5610e-04; 9.7192e-05; -6.0946e-06] ;
    titstr = 'Individ 28_2' ;
  end ;



  %  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



  subplot(2,2,1) ;
    plot(rgrid,mleg) ;
      title('Legendre poly''s') ;
      legend('l=0','l=1','l=2','l=3','l=4','l=5','l=6',-1) ;

  subplot(2,2,3) ;
    plot(rgrid,mrecon) ;
      title([titstr ', Heights']) ;
      legend('l=0','l=1','l=2','l=3','l=4','l=5','l=6',-1) ;

  subplot(2,2,2) ;
    plot(rgrid,mcurv) ;
      title([titstr ', Curvature']) ;
      legend('l=0','l=1','l=2','l=3','l=4','l=5','l=6',-1) ;




elseif ipart == 2 ;


end ;

