function ht = corh(fm,r,th,r0) 
% CORH, For the cornea study, this reconstructs the "heights"
%    function H(r,theta), for the Radial Fourier Legendre basis.
% Inputs:
%     fm  - feature matrix (rows ~ Legendre's
%                           cols ~ Fourier's, in cos-sin form
%                                    (so, odd number of these)
%                           for both: "1st entry ~ 0th coeff")
%                           
%     r   - matrix of radii of polar coordinates
%     th  - matrix of angles (radians) of polar coordinates
%                   (same size as r)
%     r0  - radius of analysis region
% Output:
%     ht  - matrix of heights of reconstructed surface
%                   (same size as r)


ht = zeros(size(r)) ;

ell0 = size(fm,1) - 1 ;
          %  Highest Legendre order, 
          %  i.e. number rows in feature matrix - 1
m0 = (size(fm,2) - 1) / 2 ;
          %  Highest Fourier order, 
          %  i.e. (number rows in feature matrix - 1) / 2 

mpolyarg = (2 * r ./ r0) - 1 ;

    %  Do ell = 0 Legendre
llleg = ones(size(mpolyarg)) ;
          %  Last, last Legendre for future updating
trigfact = fm(1,1) ;
for m = 1:m0 ;
  trigfact = trigfact + fm(1,2*m) .* cos(m * th) ;
  trigfact = trigfact + fm(1,2*m+1) .* sin(m * th) ;
end ;
ht = ht + trigfact .* llleg ;

%  Do ell = 1 Legendre
lleg = mpolyarg ;
          %  Last Legendre for future updating
trigfact = fm(2,1) ;
for m = 1:m0 ;
  trigfact = trigfact + fm(2,2*m) .* cos(m * th) ;
  trigfact = trigfact + fm(2,2*m+1) .* sin(m * th) ;
end ;
ht = ht + trigfact .* lleg ;

%  Loop through rest of Legendre's
for ell = 2:ell0 ;
  leg = (1/ell) .* ((2*ell-1) .* mpolyarg .* lleg - (ell-1) .* llleg) ;
          %  Legendre update formula

  trigfact = fm(ell+1,1) ;
  for m = 1:m0 ;
    trigfact = trigfact + fm(ell+1,2*m) .* cos(m * th) ;
    trigfact = trigfact + fm(ell+1,2*m+1) .* sin(m * th) ;
  end ;
  ht = ht + trigfact .* leg ;
      
  llleg = lleg ;
  lleg = leg ;
          %  Update for next run through loop
end ;

