Spherical Harmonics integration

[weberscode]

Dear Community,

is there already a library for the integration of spherical harmonics in Modelica?

I have started something. Could someone proof it? It is important that the function outputs complex-valued later, since I want to work with a linear combination of several bases.

Thanks for the feedback.

Kind regards

Simon Weber, M.Sc.
Research Associate

Institute for Acoustics and Building Physics
Pfaffenwaldring 7
70569 Stuttgart
Germany

Tel: +49 711 685 66301
E-Mail: [email protected]
Website: http://www.iabp.uni-stuttgart.de
https://orcid.org/0000-0002-6268-8547

function sphericalHarmonics
  input Integer M;
  input Integer L;
  input Real theta "Altitude within (0,pi]";
  input Real phi " Azimuth within [0,2pi]";
  output Complex spharm;

protected 
  parameter Real pi = Modelica.Constants.pi;
  Real Nlm = sqrt( (2*L+1)/2 * factorial(L-M)/factorial(L+M));
  Real Plm = if M<0 then (-1)^M * factorial(L-M) / factorial(L+M) * LegendrePolynom(M,L,cos(theta)) else LegendrePolynom(M,L,cos(theta));
  Complex Elm = Modelica.ComplexMath.exp(Complex(0,M*phi));

algorithm 
  spharm :=sqrt(1/2/pi)*Nlm*Plm*Elm;

end sphericalHarmonics;

function LegendrePolynom
  input Integer M;
  input Integer L;
  input Real x "Cos(Altitude) within [1,-1]";
  output Real legPoly;

protected 
  parameter Integer L2Floor = integer(floor(L/2));
  Real dPL = 0;
  Real dx;
algorithm 

  for k in 0:L2Floor loop
    dx := if noEvent(x == 0 and L - 2*k - M == 0) then factorial(L - 2*k)/factorial(L - 2*k - M) elseif noEvent(x == 0) then 0 else factorial(L
       - 2*k)/factorial(L - 2*k - M)*x^(L - 2*k - M) "M-th derivative wrt x of x^(L-2*k); This avoids 0^0 error";
    dPL := dPL + (-1)^k*factorial(2*L - 2*k)/(factorial(k)*factorial(L - k)*factorial(L - 2*k))*dx;
  end for;

  legPoly := if noEvent(1-x^2==0 and M/2==0) then (-1)^M * 1 / 2^L * dPL elseif noEvent(1-x^2==0) then 0 else (-1)^M * (1-x^2)^(M/2) / 2^L * dPL;

end LegendrePolynom;

function factorial "n-th falling factorial"
  input Integer n;
  output Integer f;

protected 
  Integer maxInt = 2147483646 "Max 32-bit integer";

algorithm 
    f := 1;
  for i in 0:(n-1) loop
    assert(f < maxInt, "Integer overflow");
    f := f*(n-i);
  end for;

end factorial;