//  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 [gr] = CL_sphHarmGrad(pos, a, f, znm, maxnm, inc00)
// Spherical harmonics - gradient
//
// Calling Sequence
// gr = CL_sphHarmGrad(pos, a, f, znm  [, maxnm, inc00])
//
// Description
// <itemizedlist><listitem>
// <p>Computes the gradient of 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. 
// gr : Gradient (3xN) 
//
// Authors
// CNES - DCT/SB (AL)
//
// Bibliography
// 1) On the computation of the spherical harmonic terms ... Leland R. Cunningham 
//
// See also
// CL_sphHarmVal
//
// 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_sphHarmGrad(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 
// Original cunningham algorithm, except for:
// - radius (r) replaced by normalized radius (r/a) to avoid overflow (r^n)
// - complex harmonics znm are Cnm + %i * Snm (and not Cnm - %i * Snm) 
// Also: indices for all arrays are incremented by 1 as indices
// begin at 1 and not 0 (V00 <=> V(1,1)). 
// 
// Vectorized algorithm: 
// - loop over "n"
// - vectorized for pos and "m" 
//
// 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 [gr] = sphHarmGrad(pos, nmax, a, f, znm)

  pos_sph = CL_co_car2sph(pos); 

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

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

  gr = zeros(pos); 

  // Vnm = elementary potentials with normalized 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-2
  // V2(i,j): j:as many as positions; i: tesseral terms
  // V2(1,:) <=> Vn,0
  // V2(2,:) <=> Vn,1
  // V2(3,:) <=> Vn,2

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

  // derivatives: Wx = dV/dx... 
  Wix = zeros(nmax+1, N); 
  Wiy = zeros(nmax+1, N); 
  Wiz = zeros(nmax+1, N); 


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

  // not yet a gradient: will be scaled later (* f/a^2)
  gr(1,:) = real(-znm(1,1)' * real(V1(2,:))); // -Z00 * Re(V11)
  gr(2,:) = real(-znm(1,1)' * imag(V1(2,:))); // -Z00 * Im(V11)
  gr(3,:) = real(-znm(1,1)' * V1(1,:));       // -Z00 * V10


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

  for (n = 2:(nmax+1))
    // 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,:) ;

    // derivatives Wij, computed for i=n-1, j=1:i
    // Wi(j+1,:) = d/dx(Vi,j)
    i = n-1; 
  
    // Wi,0
    Wix(1,:) = -real(V2(2,:));
    Wiy(1,:) = -imag(V2(2,:));
    Wiz(1,:) = -real((i+1) * V2(1,:));

    // Wi,j for j = 1:i
    j = (1:i)'; 
    if (j <> [])
      Wix(j+1,:) = (-V2(j+2,:) + (((i-j+2).*(i-j+1))*ones(1,N)) .* V2(j,:)) / 2;
      Wiy(j+1,:) = ( V2(j+2,:) + (((i-j+2).*(i-j+1))*ones(1,N)) .* V2(j,:)) * %i/2;
      Wiz(j+1,:) = -((i-j+1)*ones(1,N)) .* V2(j+1,:);
    end
  
    // includes effect of Ci,j and Si,j for j=0:i
    // note: znm' : transpose, conjugate
    Z = repmat(znm(i+1,1:i+1)', 1, N); 
    gr(1,:) = gr(1,:) + sum(real(Z .* Wix(1:i+1,:)), "r");  
    gr(2,:) = gr(2,:) + sum(real(Z .* Wiy(1:i+1,:)), "r");  
    gr(3,:) = gr(3,:) + sum(real(Z .* Wiz(1:i+1,:)), "r");  
  
    V0 = V1; 
    V1 = V2; 
  end
  
  // note: division by a^2 : because formulas are computed with normalized radius
  gr = CL_dMult(gr, f ./ a.^2); 

endfunction


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

  pos_sph = CL_co_car2sph(pos); 

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

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

  gr = zeros(pos); 

  // Vnm = elementary potentials with normalized 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-2
  // V(i,j): j:as many as positions; i: tesseral terms
  // V2(1,:) <=> Vn,0

  V0 = zeros(2, N); 
  V1 = zeros(2, N); 
  V2 = zeros(2, N); 

  // V0 = [V00]
  V0(1,:) = 1 ./ s; 

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

  // not yet a gradient: will be scaled later (* f/a^2)
  if (znm(1,1) <> 0)
    gr(1,:) = real(-znm(1,1)' * real(V1(2,:))); // -Z00 * Re(V11)
    gr(2,:) = real(-znm(1,1)' * imag(V1(2,:))); // -Z00 * Im(V11)
    gr(3,:) = real(-znm(1,1)' * V1(1,:));       // -Z00 * V10
  end
  
  // V2 = [V20; V21]
  // then: [V30; V31], etc...

  for (n = 2:(nmax+1))
    V2(1,:) = (((2*n-1) * sinphi ./ s) .* V1(1,:) - ((n-1)./(s.^2)) .* V0(1,:)) / n;
    V2(2,:) = (((2*n-1) * sinphi ./ s) .* V1(2,:) - (n./(s.^2)) .* V0(2,:)) / (n-1);

    if (znm(n,1) <> 0)
      gr(1,:) = gr(1,:) + real(-znm(n,1)' * real(V2(2,:)));
      gr(2,:) = gr(2,:) + real(-znm(n,1)' * imag(V2(2,:)));
      gr(3,:) = gr(3,:) + real(-znm(n,1)' * n * real(V2(1,:))); 
    end
  
    V0 = V1; 
    V1 = V2; 
  end
  
  // note: division by a^2 : because formulas are computed with normalized radius
  gr = CL_dMult(gr, f ./ a.^2); 

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

// replace -1 by dim(1)-1 or dim(2)-1 
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
  gr = sphHarmGrad0(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
  gr = sphHarmGrad(pos, nmax, a, f, Znm); 
end

endfunction