disp('MATLAB Script File: arvo1c.m') ;
%
%   For the Cornea Study, this makes files for ARVO poster
%         of April 1997
%   This is for Fourier and Legendre plot files
%
format compact ;

ipart = 2     %   1 -  generate Fourier basis functions
              %   2 -  generate Legendre basis functions


if ipart == 1 ;    %  then do Fouriers
  deggrid = linspace(0,360,61)' ;

  mout = ones(length(deggrid),1) ;
  for ifreq = 1:5 ;
    mout = [mout cos(ifreq * 2 * pi * deggrid / 360)] ;
  end ;

  plot(deggrid,mout) ;

  %  WRITE ANSWERS INTO A FILE
    fid = fopen('/matlab/steve/arvo1cf.out','wt') ;
      forstr = '%6.3f, %6.3f, %6.3f, %6.3f, %6.3f, %6.3f, %6.3f\n' ;
    cntbytes = fprintf(fid,forstr,[deggrid mout]') ;
    fclose('all') ;



elseif ipart == 2 ;    %  then do Legendre's



  %  Generate matrix of Legendre's
  ell0 = 6 ;
  mpolyarg = linspace(-1,1,61)' ;

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



  plot(mpolyarg,mleg) ;

  %  WRITE ANSWERS INTO A FILE
    fid = fopen('/matlab/steve/arvo1cl.out','wt') ;
      forstr = '%6.3f, %6.3f, %6.3f, %6.3f, %6.3f, %6.3f, %6.3f, %6.3f\n' ;
    cntbytes = fprintf(fid,forstr,[mpolyarg mleg]') ;
    fclose('all') ;


end ;

