21 Oct 2007

Asymptote using graph3.asy – fig0180

Figure 0018

The spherical harmonics latex2png equation are the angular portion of the solution to Laplace's equation in spherical coordinates where azimuthal symmetry is not present.

The spherical harmonics are defined by:

latex2png equation

where latex2png equation and latex2png equation is the Legendre polynomial.

Source

(Compiled with Asymptote version 2.14svn-r5318)
    
import palette;
import math;
import graph3;

typedef real fct(real);
typedef pair zfct2(real,real);
typedef real fct2(real,real);

real binomial(real n, real k)
{
  return gamma(n+1)/(gamma(n-k+1)*gamma(k+1));
}

real factorial(real n) {
  return gamma(n+1);
}

real[] pdiff(real[] p)
{ // p(x)=p[0]+p[1]x+...p[n]x^n
  // retourne la dérivée de p
  real[] dif;
  for (int i : p.keys) {
    if(i != 0) dif.push(i*p[i]);
  }
  return dif;
}

real[] pdiff(real[] p, int n)
{ // p(x)=p[0]+p[1]x+...p[n]x^n
  // dérivée n-ième de p
  real[] dif={0};
  if(n >= p.length) return dif;
  dif=p;
  for (int i=0; i < n; ++i)
    dif=pdiff(dif);
  return dif;
}

fct operator *(real y, fct f)
{
  return new real(real x){return y*f(x);};
}

zfct2 operator +(zfct2 f, zfct2 g)
{// Défini f+g
  return new pair(real t, real p){return f(t,p)+g(t,p);};
}

zfct2 operator -(zfct2 f, zfct2 g)
{// Défini f-g
  return new pair(real t, real p){return f(t,p)-g(t,p);};
}

zfct2 operator /(zfct2 f, real x)
{// Défini f/x
  return new pair(real t, real p){return f(t,p)/x;};
}

zfct2 operator *(real x,zfct2 f)
{// Défini x*f
  return new pair(real t, real p){return x*f(t,p);};
}

fct fct(real[] p)
{ // convertit le tableau des coefs du poly p en fonction polynôme
  return new real(real x){
    real y=0;
    for (int i : p.keys) {
      y += p[i]*x^i;
    }
    return y;
  };
}

real C(int l, int m)
{
  if(m < 0) return 1/C(l,-m);
  real OC=1;
  int d=l-m, s=l+m;
  for (int i=d+1; i <=s ; ++i) OC *= i;
  return 1/OC;
}

int csphase=-1;
fct P(int l, int m)
{ // Polynôme de Legendre associé
  // http://mathworld.wolfram.com/LegendrePolynomial.html
  if(m < 0) return (-1)^(-m)*C(l,-m)*P(l,-m);
  real[] xl2;
  for (int k=0; k <= l; ++k) {
    xl2.push((-1)^(l-k)*binomial(l,k));
    if(k != l) xl2.push(0);
  }
  fct dxl2=fct(pdiff(xl2,l+m));
  return new real(real x){
    return (csphase)^m/(2^l*factorial(l))*(1-x^2)^(m/2)*dxl2(x);
  };
}

zfct2 Y(int l, int m)
{// http://fr.wikipedia.org/wiki/Harmonique_sph%C3%A9rique#Expression_des_harmoniques_sph.C3.A9riques_normalis.C3.A9es
  return new pair(real theta, real phi) {
    return sqrt((2*l+1)*C(l,m)/(4*pi))*P(l,m)(cos(theta))*expi(m*phi);
  };
}

real xyabs(triple z){return abs(xypart(z));}

size(16cm);
currentprojection=orthographic(0,1,1);

zfct2 Ylm;

triple F(pair z)
{
  //   real r=0.75+dot(0.25*I,Ylm(z.x,z.y));
  //   return r*expi(z.x,z.y);
  real r=abs(Ylm(z.x,z.y))^2;
  return r*expi(z.x,z.y);
}

int nb=4;
for (int l=0; l < nb; ++l) {
  for (int m=0; m <= l; ++m) {
    Ylm=Y(l,m);

    surface s=surface(F,(0,0),(pi,2pi),60);
    s.colors(palette(s.map(xyabs),Rainbow()));

    triple v=(-m,0,-l);
    draw(shift(v)*s);
    label("$Y_"+ string(l) + "^" + string(m) + "$:",shift(X/3)*v);
  }
}

Étiquettes : , , , , ,