//  Copyright (c) CNES  2008
//
//  This software is part of CelestLab, a CNES toolbox for Scilab
//
//  This software is governed by the CeCILL  license under French law and
//  abiding by the rules of distribution of free software.  You can  use,
//  modify and/ or redistribute the software under the terms of the CeCILL
//  license as circulated by CEA, CNRS and INRIA at the following URL
//  'http://www.cecill.info'.

function [val] = CL_sphHarmVal(pos, a, f, znm, maxnm, inc00)
// Spherical harmonics - value
//
// Calling Sequence
// val = CL_sphHarmVal(pos, a, f, znm  [, maxnm, inc00])
//
// Description
// <itemizedlist><listitem>
// <p>Computes the expression: </p> 
// <imageobject><imagedata fileref="sph_harm.png" /></imageobject> 
// <p></p></listitem>
// <listitem>
// <p>Notes:</p> 
// <p> - <b>maxnm</b> has the form [nmax,mmax]. A value of -1 for nmax or mmax 
// means "all available znm coefficients are used". </p>
// </listitem>
// </itemizedlist>
//
// Parameters
// pos : [x;y;z] Position vector in cartesian coordinates [m] (3xN)
// a : Reference radius (1x1 or 1xN)
// f : Multiplying coefficient (1x1 or 1xN)
// znm : Complex harmonics: znm(n+1, m+1) = Cnm + %i * Snm (PxQ)
// maxnm : (integer, optional) Maximum degree/order to consider: [nmax,mmax] (1x2). Default is [-1,-1].
// inc00 : (boolean, optional) True if harmonic coefficient z00 is included, false otherwise. Default is true. 
// val : Value (1xN) 
//
// Authors
// CNES - DCT/SB (AL)
//
// Bibliography
// 1) On the computation of the spherical harmonic terms ... Leland R. Cunningham 
//
// See also
// CL_sphHarmGrad
//
// Examples
// pos = [[7000.e3; 0; 0], [0; 7000.e3; 0], [0; 0; 7000.e3]];
// er = CL_dataGet("eqRad"); 
// mu = CL_dataGet("mu"); 
// j1jn =  CL_dataGet("j1jn");
// znm = [1; 0; -j1jn(2)];
// CL_sphHarmVal(pos, er, mu, znm)


// Reference: 
// On the computation of the spherical harmonic terms needed during the numerical integration 
// of the orbital motion of an artificial satellite, Leland R. Cunningham, 
// Lockheed Missiles and Space Company, Sunnyvale, Calif., and Astronomy Dept., 
// University of California, Berkeley, Calif., U.S.A.

// --------------------------
// Internal function 1
// --------------------------
// General case 
// NB: size(znm) must be at least: (nmax+1) x (nmax+1)
// znm = Cnm + %i * Snm
// znm(i,j) <>= Cnm i-1,j-1 + %i * Snm i-1,j-1
function [val] = sphHarmVal(pos, nmax, a, f, znm)

  // spherical coordinates
  pos_sph = CL_co_car2sph(pos); 

  lambda = pos_sph(1,:);
  cosphi = cos(pos_sph(2,:)); 
  sinphi = sin(pos_sph(2,:)); 
  r = pos_sph(3,:); 

  // normalized radius (%nan if <= 0 (cannot be < 0))
  s = r ./ a; 
  I = find(s <= 0); 
  s(I) = %nan; 

  val = zeros(r); 

  // Vnm = elementary potentials with normallized radius
  // Vnm = Pnm(sinphi) * exp(%i*m*lambda) / s^(n+1)
  // (n>=0 m=0:n)
  // In algorithm: V2 represents Vn, V1 <=> Vn-1, V0 <=> Vn-1
  // V(i,j): j:as many as positions; i: tesseral terms
  // V2(1,:) <=> Vn,0
  // V2(2,:) <=> Vn,1
  // V2(3,:) <=> Vn,2

  // first index: m
  // nmax+2 : because computes from 0 to nmax+1
  N = size(pos,2); 
  V0 = zeros(nmax+1, N); 
  V1 = zeros(nmax+1, N); 
  V2 = zeros(nmax+1, N); 

  // V0 = [V00]
  V0(1,:) = 1 ./ s; 
  val = real(znm(1,1)' * V0(1,:)); 
  
  // V1 = [V10; V11]
  V1(1,:) = (sinphi./s) .* V0(1,:); 
  V1(2,:) = (cosphi./s) .* V0(1,:) .* exp(%i * lambda); 

  if (nmax >= 1)
  // V10 and V11
  val = val + real(znm(2,1)' * V1(1,:) + znm(2,2)' * V1(2,:));  
  end

  // V2 = [V20; V21; V22]
  // then: [V30; V31; V32; V33], etc...

  for (n = 2:nmax)
    // Vn,n
    V2(n+1,:) = (2*n-1) * (cosphi./s) .* V1(n,:) .* exp(%i * lambda);
  
    // Vn,n-1
    V2(n,:) = (2*n-1) * (sinphi ./ s) .* V1(n,:);
  
    // Vn,m  for m = 0 : n-2
    m = (0:(n-2))';
    V2(m+1,:) = (((2*n-1) * (1 ./ (n-m))) * (sinphi ./ s)) .* V1(m+1,:) - ..
                (((n+m-1) ./ (n-m)) * (1 ./ s.^2))         .* V0(m+1,:) ;


    // note: znm' : transpose, conjugate
    Z = repmat(znm(n+1,1:(nmax+1))', 1, N); 
    val = val + sum(real(Z .* V2(1:(nmax+1),:)), "r");
  
    V0 = V1; 
    V1 = V2; 
  end
  
  // note: division by a : because formulas are computed with normalized radius
  val = (val .* f) ./ a; 

endfunction


// --------------------------
// Internal function 2
// --------------------------
// specialized (simpler) version for mmax = 0 (for efficiency)
// Same inputs / algorithms as above
// znm(1:nmax+1,1) only are used
function [val] = sphHarmVal0(pos, nmax, a, f, znm)

  // spherical coordinates
  pos_sph = CL_co_car2sph(pos); 

  lambda = pos_sph(1,:);
  cosphi = cos(pos_sph(2,:)); 
  sinphi = sin(pos_sph(2,:)); 
  r = pos_sph(3,:); 

  // normalized radius (%nan if <= 0 (cannot be < 0))
  s = r ./ a; 
  I = find(s <= 0); 
  s(I) = %nan; 

  val = zeros(r); 

  // Vnm = elementary potentials with normallized radius
  // Vnm = Pnm(sinphi) * exp(%i*m*lambda) / s^(n+1)
  // (n>=0 m=0)
  // In algorithm: V2 represents Vn, V1 <=> Vn-1, V0 <=> Vn-1
  // V(1,j): j:as many as positions
  // V2(1,:) <=> Vn,0

  N = size(pos,2); 
  V0 = zeros(1, N); 
  V1 = zeros(1, N); 
  V2 = zeros(1, N); 

  // V0 = [V00]
  V0(1,:) = 1 ./ s; 
  val = real(znm(1,1)' * V0(1,:)); 
 
  // V1 = [V10]
  V1(1,:) = (sinphi./s) .* V0(1,:); 
  if (nmax >= 1)
    val = val + real(znm(2,1)' * V1(1,:));  
  end

  // V2 = [V20], then: [V30], etc...
  for (n = 2:nmax) 
    // Vn,0  
    V2(1,:) = (((2*n-1) * (1 ./ (n))) * (sinphi ./ s)) .* V1(1,:) - ..
                (((n-1) ./ n) * (1 ./ s.^2))           .* V0(1,:) ;

    // note: znm' : conjugate
    val = val + real(znm(n+1,1)' * V2(1,:));
  
    V0 = V1; 
    V1 = V2; 
  end
  
  // note: division by a : because formulas are computed with normalized radius
  val = (val .* f) ./ a; 

endfunction


// --------------------------
// MAIN
// --------------------------
if (~exists("maxnm", "local")); maxnm = [-1,-1]; end
if (~exists("inc00", "local")); inc00 = %t; end

// check inputs
if (length(maxnm) <> 2)
  CL__error("Invalid size for ''maxnm''");
end

// NOTE code below allows: maxnm == []
if (maxnm <> [])
  I = find(maxnm <> round(maxnm) | (maxnm < 0 & maxnm <> -1))
  if (I <> [])
    CL__error("Invalid value for ''maxnm''");
  end
end

// check size / resize
[pos, a, f] = CL__checkInputs(pos,3, a,1, f,1); 

// returns empty if inputs are empty
if (znm == [] | pos == [] | a == [] | f == [])
  val = []; 
  return; 
end

dim = size(znm); 

// find maximum values for n and m (nmax and mmax)
if (maxnm == [])
  maxnm = [dim(1)-1, dim(2)-1]; 
end

// maxnm now contains 2 elements 
I = find(maxnm == -1); 
if (I <> []); maxnm(I) = dim(I)-1; end

nmax = maxnm(1);
mmax = maxnm(2); 

if (nmax > dim(1)-1 | mmax > dim(2)-1)
  CL__error("Too few harmonic coefficients"); 
end

// optimize computation if znm contains 0
[i,j] = find(znm <> 0); 
if (i <> [])
  nmax = min(nmax, max(i)-1); 
  mmax = min(mmax, max(j)-1); 
end


if (mmax == 0)
  // z00 = 0 if not included
  if (~inc00); znm(1,1) = 0; end
  val = sphHarmVal0(pos, nmax, a, f, znm); 
else
  // Znm: copy of harmonics used (square matrix, filled with 0) 
  Znm = zeros(nmax+1, nmax+1);
  Znm(1:nmax+1, 1:mmax+1) = znm(1:nmax+1, 1:mmax+1); 
  if (~inc00); Znm(1,1) = 0; end
  val = sphHarmVal(pos, nmax, a, f, Znm); 
end

endfunction