print("Definition of greek characters");//ギリシア文字定義
gth:=Symbol::theta;
gph:=Symbol::phi;

print("Legendre polynomial");//ルジャンドル多項式
Pl:=(x,n)->expand(1/2^n/n!*diff((x^2-1)^n,x$n));

print("Associated Legendre function");//ルジャンドル陪関数
Plm:=(x,l,m)->expand((-1)^m*(1-x^2)^(m/2)*diff(Pl(x,l),x$m));

print("applied absolute value");//絶対値化
Plm2:=(x,l,m)->piecewise([m<0,Plm(x,l,-m)],[Otherwise,Plm(x,l,m)]);
hold(Plm2)(x,3,3)=Plm2(x,3,3);
hold(Plm2)(x,3,-3)=Plm2(x,3,-3);

print("replace x->cos");//置換
Plm3:=(gth,l,m)->expand(subs(Plm2(x,l,m),x^2=1-sin(gth)^2,//ここに入れるのがミソ
                                  x=cos(gth),(sin(gth)^2)^(1/2)=sin(gth),
                                  (sin(gth)^2)^(3/2)=sin(gth)^3,(sin(gth)^2)^(5/2)=sin(gth)^5,
                                  (sin(gth)^2)^(7/2)=sin(gth)^7,(sin(gth)^2)^(9/2)=sin(gth)^9));
hold(Plm2)(x,1,-1)=Plm2(x,1,-1);
hold(Plm3)(gth,1,-1)=Plm3(gth,1,-1);
hold(Plm3)(gth,3,-3)=Plm3(gth,3,-3);
hold(Plm3)(gth,9,-7)=Plm3(gth,9,-7);

print("spherical harmonic function");//球面調和関数
Ylm:=(gth,gph,l,m)->(-1)^((m+abs(m))/2)*sqrt((2*l+1)/(4*PI)*(l-abs(m))!/(l+abs(m))!)*Plm3(gth,l,m)*E^(I*m*gph);

for l from 0 to 5 do
  for m from -l to l do
    print(hold(Ylm)(gth,gph,l,m)=Ylm(gth,gph,l,m));
    //print(hold(Ylm)(gth,gph,l,m)=Ylm(gth,gph,l,m));
  end_for
end_for;

for l from 0 to 5 do
  for m from -l to l do
    S:=Ylm(gth,gph,l,m);
    print(hold(Ylm)(gth,gph,l,m)=S);
    plot(
    plot::Spherical([real(S),gph,gth],gth=-PI..PI,gph=0..PI,
    UMesh=60,VMesh=60,FillColor=[1,0,0],FillColor2=[0,0,1]),
    plot::Spherical([imag(S),gph,gth],gth=-PI..PI,gph=0..PI,
    UMesh=60,VMesh=60,FillColor=[1,1,0],FillColor2=[0,1,1]),
    Scaling=Constrained);
  end_for
end_for;
//plot(plot::Spherical([real(Ylm(gth,gph,3,-1)),gph,gth],
//gth=-PI..PI,gph=0..PI,UMesh=60,VMesh=60,FillColorType=Flat,
//Color=[0.5,0.5,0.5]),Scaling=Constrained);