Spherical harmonics

One way to represent the spherical harmonics is to consider their real and imaginary parts and assign a tone corresponding to their values at each point on the sphere. The first figure in each table shows the real part of \(Y_\ell^m\) with the jet colormap that assigns reddish colors to higher values and bluish colors to lower values. The second figure shows the result with lighting and a single color in three dimensions, with the radius proportional to the absolute value of the real part. Thanks to the identity \(\sin\varphi=\cos(\pi/2-\varphi)\), taking imaginary parts results in similar figures with rotation and symmetry. The formula for \(Y_\ell^m\) in terms of the spherical coordinates \(\theta\) and \(\varphi\) is shown to the right of the figures.

Due to the symmetry \(Y_\ell^{-m}=(-1)^m \overline{Y_\ell^{m}}\) it is only considered the case \(m\ge 0\). Adjusting the parameters in the code below one can change the ranges of \(\ell\) and \(m\). There are actually two programs: one in Octave for generating the images and another in Sagemath for automatically building this page from an HTML document template.

The code




for ll = [0:5]

  for mm = [0:ll]

    figure(1)

    % real part (surface map)

    [Xc, Yc, Zc, Z, Xc2, Yc2, Zc2, Z2] = plot_sph_harm(ll,mm);

    surf(Xc, Yc, Zc, Z, 'EdgeColor','none');

    colormap(jet);

    axis equal off;

    saveas(gcf, sprintf("./images/spha%d%d.jpg", ll, mm))

    figure(2)

    % absolute values of the real part as radius

    [Xc, Yc, Zc, Z, Xc2, Yc2, Zc2, Z2] = plot_sph_harm(ll,mm);

    surf(Xc2, Yc2, Zc2, 0*Z2, 'EdgeColor','none');

    colormap(jet);

    aa = 0.53;

    axis([-aa aa -aa aa -aa aa]);    axis equal off;

    light;

    lighting gouraud; material dull;

    saveas(gcf, sprintf("./images/spha%d%d_2.jpg", ll, mm))

  end

end



%colormaps: viridis, jet, hsv, hot, cool, and gray


% plot_sph_harm.m



function [Xc, Yc, Zc, Z, Xc2, Yc2, Zc2, Z2] = plot_sph_harm(l, m)



  [PHI, THETA] = meshgrid(linspace(0, 2*pi, 640), linspace(0, pi, 480));



  Y = sph_harm_octave(l, m, THETA, PHI);



  Z = real(Y);

  Z2 = abs(Z);

  R = 0.5+0.0*Z2;

  if l*l+m*m>0 % exlude the case l=m=0

    MM = max(Z2(:));

    mm = min(Z2(:));

    R = 0.5 * (Z2 - mm) / (MM-mm); % automatic scale

  end



  % convert to Cartesian

  Xc = sin(THETA) .* cos(PHI);

  Yc= sin(THETA) .* sin(PHI);

  Zc = cos(THETA);

  Xc2 = R .* Xc;

  Yc2 = R .* Yc;

  Zc2 = R .* Zc;

end





% This part was done by an AI

function Y = sph_harm_octave(l, m, theta, phi)

  % Computes complex spherical harmonic Y_l^m(θ,φ)

  % Uses Condon-Shortley phase (includes (-1)^m) and Schmidt normalization.

  % theta,phi may be matrices (same size)

  if m < 0

    % relation for negative m: Y_l^{-m} = (-1)^m * conj(Y_l^{m})

    Ypos = sph_harm_octave(l, -m, theta, phi);

    Y = (-1)^m * conj(Ypos);

    return;

  end



  % compute associated Legendre P_l^m(cos theta)

  x = cos(theta);

  Plm = assoc_legendre(l, m, x);   % same size as theta



  % normalization factor: sqrt((2l+1)/(4π) * (l-m)!/(l+m)!)

  norm = sqrt((2*l+1)/(4*pi) * factorial(l-m)/factorial(l+m));



  Y = norm * Plm .* exp(1i * m * phi);

end



function P = assoc_legendre(l, m, x)

  % Compute associated Legendre polynomials P_l^m(x) for scalar l,m and array x

  % using recurrence. Handles 0 <= m <= l.

  x = double(x);

  sz = size(x);

  x = x(:);



  % handle m==0 separately

  if m == 0

    P0 = ones(size(x));

    if l == 0

      P = reshape(P0, sz); return;

    end

    P1 = x; % P_1^0(x) = x

    if l == 1

      P = reshape(P1, sz); return;

    end

    Pn_2 = P0; Pn_1 = P1;

    for n = 2:l

      Pn = ((2*n-1).*x.*Pn_1 - (n-1).*Pn_2)./n;

      Pn_2 = Pn_1; Pn_1 = Pn;

    end

    P = reshape(Pn, sz); return;

  end



  % m > 0: use standard initial Pmm

  % P_m^m(x) = (-1)^m * (2m-1)!! * (1 - x^2)^(m/2)

  % compute double factorial (2m-1)!! safely

  df = 1;

  for k = 1:m

    df = df * (2*k - 1);

  end

  Pmm = (-1)^m * df .* (1 - x.^2).^(m/2);



  if m == l

    P = reshape(Pmm, sz); return;

  end



  % P_{m+1}^m = x*(2m+1)*Pmm

  Pm1 = x .* (2*m + 1) .* Pmm;

  if (m+1) == l

    P = reshape(Pm1, sz); return;

  end



  % upward recurrence

  Pn_2 = Pmm; Pn_1 = Pm1;

  for n = m+2 : l

    Pn = ( (2*n - 1) .* x .* Pn_1 - (n + m - 1) .* Pn_2 ) ./ (n - m);

    Pn_2 = Pn_1; Pn_1 = Pn;

  end

  P = reshape(Pn, sz);

end


# Make sphhart.html



theta = var('theta')



N = 5



# Generate Legendre polynomials



Le = [1+0*x,x]



for k in srange(1,N):

  Le.append(  ( (2*k+1)/(k+1)*x*Le[k]-k/(k+1)*Le[k-1] ).expand()  )





# HTML TEMPLATE FOR TABLE

###

### LLL MMM IMGRNAME IMGINAME LATXDF

templ_html = r'<table style="border-style:solid;border-width:1px;margin-left: 30px; margin-right: auto;"><thead>'

templ_html += '\n' + r'<tr><th colspan="3" style="border-style:solid;border-width:1px;">'

templ_html += '\n' + r'<span class="math inline">'

templ_html += '\n' + r'\(\ell=LLL,\qquad m=MMM\)'

templ_html += '\n' + r'</span></th></tr></thead><tbody><tr>'

templ_html += '\n' + r'<td>'

templ_html += '\n' + r'<img style="clip-path: circle(27%);width: 550px; margin: -70px -110px -75px -130px;" alt="IMGRNAME" src="./IMGRNAME.jpg">'

#templ_html += '\n' + r'<div style="overflow: hidden;width: 230px;height: 200px;">'

#templ_html += '\n' + r'<img style="margin: -120px 0px -300px -270px;" alt="IMGRNAME" src="./IMGRNAME.jpg">'

#templ_html += '\n' + r'</div>'

templ_html += '\n' + r'</td><td>'

templ_html += '\n' + r'<img style="clip-path: circle(27%);width: 550px; margin: -70px -110px -75px -130px;" alt="IMGINAME" src="./IMGINAME.jpg">'

#templ_html += '\n' + r'<img style="width: 130px;" alt="IMGINAME" src="./IMGINAME.jpg">'

templ_html += '\n' + r'</td> <td>'

templ_html += '\n' + r' <span class="math display">'

templ_html += '\n' + r'LATXDF'

templ_html += '\n' + r'</span></td></tr></tbody></table>'

templ_html += '\n' + r'<br>'





# MAIN LOOP



# k is l

table = ''

for k in srange(len(Le)):

  for m in srange(k+1):

    mpol = Le[k]

    if m!=0: mpol = diff( Le[k], x, m)

      

    co = (ZZ['x'](mpol*2^k)).content()

    mpol = 2^k/co*mpol

    c = co/2^k



      

      

    res = latex(mpol(x=cos(theta)))



    # azimuthal

    res = r'\big(' + res.strip() + r'\big)\,e^{' +str(m) + r'i\varphi}'



    # constant

    c *= (-1)^m*sqrt( (2*k+1)*factorial(k-m)/4/pi/factorial(k+m) )

    c = c.simplify()

    c = c.simplify_full()



    # sine

    if m!=0: c = c* (sin(theta))^m

    

    

    res = latex(c) + res

    

    #clean

    res = res.replace(r'\cos\left(\theta\right)^',r'\left(\cos\theta\right)^')

    res = res.replace(r'\cos\left(\theta\right)',r'\cos\theta')

    res = res.replace(r'\, \left(\cos',r'\left(\cos')

    res = res.replace(r'\big(1\big)',r'\,')

    res = res.replace(r'\sin\left(\theta\right)^',r'\left(\sin\theta\right)^')

    res = res.replace(r'\sin\left(\theta\right)',r'\sin\theta')  

    res = res.replace(r'{1i\varphi}',r'{i\varphi}')

    res = res.replace(r'e^{0i\varphi}',r'')  

    res = res.replace(r'\big(\cos\theta\big)',r'(\cos\theta)')  

    

    



    

    res = r'\[' + '\n' + r' ' + res + '\n' + r'\]'



    it_table = templ_html

    it_table = it_table.replace('LLL', str(k))

    it_table = it_table.replace('MMM', str(m))

    it_table = it_table.replace('IMGRNAME', r'images/spha' + str(k) +str(m))

    it_table = it_table.replace('IMGINAME', r'images/spha' + str(k) +str(m)+'_2')

    it_table = it_table.replace('LATXDF', res)

    table += it_table +'\n'

  

#print(table)







def f_save( text ):

  try:

    with open( 'sphhar.html' , 'w') as sali:

      sali.write( text )

  except IOError as e:

    sys.exit('Error writing on sphhar.html')





try:

  with open('sphhart.html', encoding='UTF-8') as entr:

#  with open('sphhart.html', encoding='iso-8859-15') as entr:

    text_file = entr.read()

    L = text_file.split('<!--REPLACE-->',1)

    text_file = L[0] + table + L[1]

    f_save( text_file )



except IOError as e:

  sys.exit('sphhart.html not found')