//
// Spherical harmonics that employ the Condon-Shortley phase
// - cartesian version with l=0,1,2,3,4(s,p,d,f,g)
// - normal version with l=1,2,3
//
// See also my website. There are all renderings.
// http://www.eva.hi-ho.ne.jp/nakajima2001/tadayuki/Scilab/Scilab.html
//
funcprot(0); //nothing special is done when a function is redefined
function [a] = Factorial(n) //n!
    a = 1;
    for i = n:-1:1 
        a = a * i;
    end
endfunction
function [r] = GetPlusValues(r0)
    r = ones(r0)*0;
    for i = 1:size(r0, "r")
        for j = 1:size(r0, "c")
            if real(r0(i,j)) > 0 then r(i,j) = r0(i,j); end
        end
    end
endfunction
function [r] = GetMinusValues(r0)
    r = ones(r0)*0;
    for i = 1:size(r0, "r")
        for j = 1:size(r0, "c")
            if real(r0(i,j)) < 0 then r(i,j) = r0(i,j); end
        end
    end
endfunction
function [r] = SphericalHarmonicY(l, m, t, p)
    r = (-1)^((m+abs(m))/2) ..
    * sqrt((2*l + 1)/(4*%pi) * Factorial(l - abs(m))/Factorial(l + abs(m))) ..
    * legendre(l, abs(m), cos(t), "norm")' * exp(%i*m*p);
endfunction
function [r] = SphericalPlot3D(r0)
    r = GetPlusValues(r0);
    x = r .* (sin(t)' * cos(p));
    y = r .* (sin(t)' * sin(p));
    z = r .* (cos(t)' * ones(p));
    plot3d2(x, y, z, theta = 300, alpha = 60, leg = "x@y@z", flag = [1 6 4]);
    e=gce(); // store the Compound containing the plot
    e.foreground = 5/8*ColormapSize;
    e.hiddencolor = 7/8*ColormapSize;

    r = -GetMinusValues(r0);
    x = r .* (sin(t)' * cos(p));
    y = r .* (sin(t)' * sin(p));
    z = r .* (cos(t)' * ones(p));
    plot3d2(x, y, z, theta = 300, alpha = 60, leg = "x@y@z", flag = [1 6 4]);
    e=gce(); // store the Compound containing the plot
    e.foreground = 0;
    e.hiddencolor = 2/8*ColormapSize;
endfunction


// Real spherical harmonics with l=0 -------------------------------------------


// s
clf(); //To clear the previous plot
xset('colormap', coppercolormap(128));
f = gcf();
f.figure_size = [620, 620];
funcprot(0); //nothing special is done when a function is redefined
function [r] = f(t, p)
	r = SphericalHarmonicY(0, 0, t, p);
endfunction
t = linspace(0.0001, %pi-0.0001, 100);
//p = linspace(0, 2*%pi, 200);
p = linspace(0, %pi, 100);
r = f(t, p);
x = r .* (sin(t)' * cos(p));
y = r .* (sin(t)' * sin(p));
z = r .* (cos(t)' * ones(p));
[xf, yf, zf] = nf3d(x', y', z'); 
plot3d1(xf, yf, zf, theta = 300, alpha = 60, leg = "x@y@z", flag = [-1 6 4]);
xtitle('$s=Y_{0}^{0}$');
a = gca(); //get the current axes
a.title.font_size = 5;
a.margins = [0.1,0.1,0.2,0.1];


// Real spherical harmonics with l=1 -------------------------------------------


// pz
clf(); //To clear the previous plot
xset('colormap', coppercolormap(128));
f = gcf();
f.figure_size = [620, 620];
funcprot(0); //nothing special is done when a function is redefined
function [r] = f(t, p)
    r = SphericalHarmonicY(1, 0, t, p);
endfunction
t = linspace(0.0001, %pi-0.0001, 100);
//p = linspace(0, 2*%pi, 200);
p = linspace(0, %pi, 100);
r = abs(f(t, p));
x = r .* (sin(t)' * cos(p));
y = r .* (sin(t)' * sin(p));
z = r .* (cos(t)' * ones(p));
[xf, yf, zf] = nf3d(x', y', z'); 
plot3d1(xf, yf, zf, theta = 300, alpha = 60, leg = "x@y@z", flag = [-1 6 4]);
xtitle('$\left|p_{z}\right|=\left|Y_{1}^{0}\right|$');
a = gca(); //get the current axes
a.title.font_size = 5;
a.margins = [0.1,0.1,0.2,0.1];


// py
clf(); //To clear the previous plot
xset('colormap', coppercolormap(128));
f = gcf();
f.figure_size = [620, 620];
funcprot(0); //nothing special is done when a function is redefined
function [r] = f(t, p)
    r = (SphericalHarmonicY(1, -1, t, p) + SphericalHarmonicY(1, 1, t, p))*%i/sqrt(2);
endfunction
t = linspace(0.0001, %pi-0.0001, 100);
//p = linspace(0, 2*%pi, 200);
p = linspace(%pi/2, 3*%pi/2, 100);
r = abs(f(t, p));
x = r .* (sin(t)' * cos(p));
y = r .* (sin(t)' * sin(p));
z = r .* (cos(t)' * ones(p));
[xf, yf, zf] = nf3d(x', y', z'); 
plot3d1(xf, yf, zf, theta = 300, alpha = 60, leg = "x@y@z", flag = [-1 5 4], ..
ebox = [-0.5,0.5,-0.5,0.5,-0.25,0.25]);
xtitle('$\left|p_{y}\right|=\left|i\sqrt{1/2}\left(Y_{1}^{-1}+Y_{1}^{1}\right)\right|$');
a = gca(); //get the current axes
a.title.font_size = 5;
a.margins = [0.1,0.1,0.2,0.1];


// px
clf(); //To clear the previous plot
xset('colormap', coppercolormap(128));
f = gcf();
f.figure_size = [620, 620];
funcprot(0); //nothing special is done when a function is redefined
function [r] = f(t, p)
    r = (SphericalHarmonicY(1, -1, t, p) - SphericalHarmonicY(1, 1, t, p))/sqrt(2);
endfunction
t = linspace(0.0001, %pi-0.0001, 100);
//p = linspace(0, 2*%pi, 200);
p = linspace(0, %pi, 100);
r = abs(f(t, p));
x = r .* (sin(t)' * cos(p));
y = r .* (sin(t)' * sin(p));
z = r .* (cos(t)' * ones(p));
[xf, yf, zf] = nf3d(x', y', z'); 
plot3d1(xf, yf, zf, theta = 300, alpha = 60, leg = "x@y@z", flag = [-1 5 4], ..
ebox = [-0.5,0.5,-0.5,0.5,-0.25,0.25]);
xtitle('$\left|p_{x}\right|=\left|\sqrt{1/2}\left(Y_{1}^{-1}-Y_{1}^{1}\right)\right|$');
a = gca(); //get the current axes
a.title.font_size = 5;
a.margins = [0.1,0.1,0.2,0.1];


// Real spherical harmonics with l=2 -------------------------------------------


// dz2
clf(); //To clear the previous plot
xset('colormap', coppercolormap(128));
f = gcf();
f.figure_size = [620, 620];
funcprot(0); //nothing special is done when a function is redefined
function [r] = f(t, p)
    r = SphericalHarmonicY(2, 0, t, p);
endfunction
t = linspace(0.0001, %pi-0.0001, 100);
//p = linspace(0, 2*%pi, 200);
p = linspace(0, %pi, 100);
r = abs(f(t, p));
x = r .* (sin(t)' * cos(p));
y = r .* (sin(t)' * sin(p));
z = r .* (cos(t)' * ones(p));
[xf, yf, zf] = nf3d(x', y', z'); 
plot3d1(xf, yf, zf, theta = 300, alpha = 60, leg = "x@y@z", flag = [-1 6 4]);
xtitle('$\left|d_{z^{2}}\right|=\left|Y_{2}^{0}\right|$');
a = gca(); //get the current axes
a.title.font_size = 5;
a.margins = [0.1,0.1,0.2,0.1];


// dxz
clf(); //To clear the previous plot
xset('colormap', coppercolormap(128));
f = gcf();
f.figure_size = [620, 620];
funcprot(0); //nothing special is done when a function is redefined
function [r] = f(t, p)
    r = (SphericalHarmonicY(2, -1, t, p) - SphericalHarmonicY(2, 1, t, p))/sqrt(2);
endfunction
t = linspace(0.0001, %pi-0.0001, 100);
//p = linspace(0, 2*%pi, 200);
p = linspace(0, %pi, 100);
r = abs(f(t, p));
x = r .* (sin(t)' * cos(p));
y = r .* (sin(t)' * sin(p));
z = r .* (cos(t)' * ones(p));
[xf, yf, zf] = nf3d(x', y', z'); 
plot3d1(xf, yf, zf, theta = 300, alpha = 60, leg = "x@y@z", flag = [-1 5 4], ..
ebox = [-0.3,0.3,-0.3,0.3,-0.3,0.3]);
xtitle('$\left|d_{xz}\right|=\left|\sqrt{1/2}\left(Y_{2}^{-1}-Y_{2}^{1}\right)\right|$');
a = gca(); //get the current axes
a.title.font_size = 5;
a.margins = [0.1,0.1,0.2,0.1];


// dyz
clf(); //To clear the previous plot
xset('colormap', coppercolormap(128));
f = gcf();
f.figure_size = [620, 620];
funcprot(0); //nothing special is done when a function is redefined
function [r] = f(t, p)
    r = (SphericalHarmonicY(2, -1, t, p) + SphericalHarmonicY(2, 1, t, p))*%i/sqrt(2);
endfunction
t = linspace(0.0001, %pi-0.0001, 100);
//p = linspace(0, 2*%pi, 200);
p = linspace(%pi/2, 3*%pi/2, 100);
r = abs(f(t, p));
x = r .* (sin(t)' * cos(p));
y = r .* (sin(t)' * sin(p));
z = r .* (cos(t)' * ones(p));
[xf, yf, zf] = nf3d(x', y', z'); 
plot3d1(xf, yf, zf, theta = 300, alpha = 60, leg = "x@y@z", flag = [-1 5 4], ..
ebox = [-0.3,0.3,-0.3,0.3,-0.3,0.3]);
xtitle('$\left|d_{yz}\right|=\left|i\sqrt{1/2}\left(Y_{2}^{-1}+Y_{2}^{1}\right)\right|$');
a = gca(); //get the current axes
a.title.font_size = 5;
a.margins = [0.1,0.1,0.2,0.1];


// dxy
clf(); //To clear the previous plot
xset('colormap', coppercolormap(128));
f = gcf();
f.figure_size = [620, 620];
funcprot(0); //nothing special is done when a function is redefined
function [r] = f(t, p)
    r = (SphericalHarmonicY(2, -2, t, p) - SphericalHarmonicY(2, 2, t, p))*%i/sqrt(2);
endfunction
//t = linspace(0.0001, %pi-0.0001, 100);
t = linspace(%pi/2, %pi-0.0001, 50);
p = linspace(0, 2*%pi, 200);
r = abs(f(t, p));
x = r .* (sin(t)' * cos(p));
y = r .* (sin(t)' * sin(p));
z = r .* (cos(t)' * ones(p));
[xf, yf, zf] = nf3d(x', y', z'); 
plot3d1(xf, yf, zf, theta = 300, alpha = 60, leg = "x@y@z", flag = [-1 5 4], ..
ebox = [-0.2,0.2,-0.2,0.2,-0.1,0.1]);
xtitle('$\left|d_{xy}\right|=\left|i\sqrt{1/2}\left(Y_{2}^{-2}-Y_{2}^{2}\right)\right|$');
a = gca(); //get the current axes
a.title.font_size = 5;
a.margins = [0.1,0.1,0.2,0.1];


// dx2-y2
clf(); //To clear the previous plot
xset('colormap', coppercolormap(128));
f = gcf();
f.figure_size = [620, 620];
funcprot(0); //nothing special is done when a function is redefined
function [r] = f(t, p)
    r = (SphericalHarmonicY(2, -2, t, p) + SphericalHarmonicY(2, 2, t, p))/sqrt(2);
endfunction
//t = linspace(0.0001, %pi-0.0001, 100);
t = linspace(%pi/2, %pi-0.0001, 50);
p = linspace(0, 2*%pi, 200);
r = abs(f(t, p));
x = r .* (sin(t)' * cos(p));
y = r .* (sin(t)' * sin(p));
z = r .* (cos(t)' * ones(p));
[xf, yf, zf] = nf3d(x', y', z'); 
plot3d1(xf, yf, zf, theta = 300, alpha = 60, leg = "x@y@z", flag = [-1 5 4], ..
ebox = [-0.2,0.2,-0.2,0.2,-0.1,0.1]);
xtitle('$\left|d_{x^{2}-y^{2}}\right|=\left|\sqrt{1/2}\left(Y_{2}^{-2}+Y_{2}^{2}\right)\right|$');
a = gca(); //get the current axes
a.title.font_size = 5;
a.margins = [0.1,0.1,0.2,0.1];


// Real spherical harmonics with l=3 -------------------------------------------


// fz3
clf(); //To clear the previous plot
xset('colormap', coppercolormap(128));
f = gcf();
f.figure_size = [620, 620];
funcprot(0); //nothing special is done when a function is redefined
function [r] = f(t, p)
    r = SphericalHarmonicY(3, 0, t, p);
endfunction
t = linspace(0.0001, %pi-0.0001, 100);
//p = linspace(0, 2*%pi, 200);
p = linspace(0, %pi, 100);
r = abs(f(t, p));
x = r .* (sin(t)' * cos(p));
y = r .* (sin(t)' * sin(p));
z = r .* (cos(t)' * ones(p));
[xf, yf, zf] = nf3d(x', y', z'); 
plot3d1(xf, yf, zf, theta = 300, alpha = 60, leg = "x@y@z", flag = [-1 6 4]);
xtitle('$\left|f_{z^{3}}\right|=\left|Y_{3}^{0}\right|$');
a = gca(); //get the current axes
a.title.font_size = 5;
a.margins = [0.1,0.1,0.2,0.1];


// fxz2
clf(); //To clear the previous plot
xset('colormap', coppercolormap(128));
f = gcf();
f.figure_size = [620, 620];
funcprot(0); //nothing special is done when a function is redefined
function [r] = f(t, p)
    r = (SphericalHarmonicY(3, -1, t, p) - SphericalHarmonicY(3, 1, t, p))/sqrt(2);
endfunction
t = linspace(0.0001, %pi-0.0001, 100);
//p = linspace(0, 2*%pi, 200);
p = linspace(0, %pi, 100);
r = abs(f(t, p));
x = r .* (sin(t)' * cos(p));
y = r .* (sin(t)' * sin(p));
z = r .* (cos(t)' * ones(p));
[xf, yf, zf] = nf3d(x', y', z'); 
plot3d1(xf, yf, zf, theta = 300, alpha = 60, leg = "x@y@z", flag = [-1 5 4], ..
ebox = [-0.25,0.25,-0.25,0.25,-0.3,0.3]);
xtitle('$\left|f_{xz^{2}}\right|=\left|\sqrt{1/2}\left(Y_{3}^{-1}-Y_{3}^{1}\right)\right|$');
a = gca(); //get the current axes
a.title.font_size = 5;
a.margins = [0.1,0.1,0.25,0.1];


// fyz2
clf(); //To clear the previous plot
xset('colormap', coppercolormap(128));
f = gcf();
f.figure_size = [620, 620];
funcprot(0); //nothing special is done when a function is redefined
function [r] = f(t, p)
    r = (SphericalHarmonicY(3, -1, t, p) + SphericalHarmonicY(3, 1, t, p))*%i/sqrt(2);
endfunction
t = linspace(0.0001, %pi-0.0001, 100);
//p = linspace(0, 2*%pi, 200);
p = linspace(%pi/2, 3*%pi/2, 100);
r = abs(f(t, p));
x = r .* (sin(t)' * cos(p));
y = r .* (sin(t)' * sin(p));
z = r .* (cos(t)' * ones(p));
[xf, yf, zf] = nf3d(x', y', z'); 
plot3d1(xf, yf, zf, theta = 300, alpha = 60, leg = "x@y@z", flag = [-1 5 4], ..
ebox = [-0.25,0.25,-0.25,0.25,-0.3,0.3]);
xtitle('$\left|f_{yz^{2}}\right|=\left|i\sqrt{1/2}\left(Y_{3}^{-1}+Y_{3}^{1}\right)\right|$');
a = gca(); //get the current axes
a.title.font_size = 5;
a.margins = [0.1,0.1,0.25,0.1];


// fxyz
clf(); //To clear the previous plot
xset('colormap', coppercolormap(128));
f = gcf();
f.figure_size = [620, 620];
funcprot(0); //nothing special is done when a function is redefined
function [r] = f(t, p)
    r = (SphericalHarmonicY(3, -2, t, p) - SphericalHarmonicY(3, 2, t, p))*%i/sqrt(2);
endfunction
t = linspace(0.0001, %pi-0.0001, 100);
p = linspace(0, 2*%pi, 200);
r = abs(f(t, p))*10;
x = r .* (sin(t)' * cos(p));
y = r .* (sin(t)' * sin(p));
z = r .* (cos(t)' * ones(p));
[xf, yf, zf] = nf3d(x', y', z'); 
plot3d1(xf, yf, zf, theta = 300, alpha = 60, leg = "10x@10y@10z", flag = [-1 5 4], ..
ebox = [-1,1,-1,1,-0.8,0.8]);
xtitle('$\left|f_{xyz}\right|=\left|i\sqrt{1/2}\left(Y_{3}^{-2}-Y_{3}^{2}\right)\right|$');
a = gca(); //get the current axes
a.title.font_size = 5;
a.margins = [0.1,0.1,0.25,0.1];


// fz(x2-y2)
clf(); //To clear the previous plot
xset('colormap', coppercolormap(128));
f = gcf();
f.figure_size = [620, 620];
funcprot(0); //nothing special is done when a function is redefined
function [r] = f(t, p)
    r = (SphericalHarmonicY(3, -2, t, p) + SphericalHarmonicY(3, 2, t, p))/sqrt(2);
endfunction
t = linspace(0.0001, %pi-0.0001, 100);
p = linspace(0, 2*%pi, 200);
//p = linspace(0, %pi, 100);
r = abs(f(t, p))*10;
x = r .* (sin(t)' * cos(p));
y = r .* (sin(t)' * sin(p));
z = r .* (cos(t)' * ones(p));
[xf, yf, zf] = nf3d(x', y', z'); 
plot3d1(xf, yf, zf, theta = 300, alpha = 60, leg = "10x@10y@10z", flag = [-1 5 4], ..
ebox = [-1,1,-1,1,-0.8,0.8]);
xtitle('$\left|f_{z\left(x^{2}-y^{2}\right)}\right|=\left|\sqrt{1/2}\left(Y_{3}^{-2}+Y_{3}^{2}\right)\right|$');
a = gca(); //get the current axes
a.title.font_size = 5;
a.margins = [0.1,0.1,0.25,0.1];


// fz(x2-y2)
clf(); //To clear the previous plot
xset('colormap', coppercolormap(128));
f = gcf();
f.figure_size = [620, 620];
funcprot(0); //nothing special is done when a function is redefined
function [r] = f(t, p)
    r = (SphericalHarmonicY(3, -3, t, p) - SphericalHarmonicY(3, 3, t, p))/sqrt(2);
endfunction
//t = linspace(0.0001, %pi-0.0001, 100);
t = linspace(%pi/2, %pi-0.0001, 50);
p = linspace(0, 2*%pi, 200);
r = abs(f(t, p))*10;
x = r .* (sin(t)' * cos(p));
y = r .* (sin(t)' * sin(p));
z = r .* (cos(t)' * ones(p));
[xf, yf, zf] = nf3d(x', y', z'); 
plot3d1(xf, yf, zf, theta = 300, alpha = 60, leg = "10x@10y@10z", flag = [-1 5 4], ..
ebox = [-0.5,0.5,-0.5,0.5,-0.12,0.12]);
xtitle('$\left|f_{x\left(x^{2}-y^{2}\right)}\right|=\left|\sqrt{1/2}\left(Y_{3}^{-3}-Y_{3}^{3}\right)\right|$');
a = gca(); //get the current axes
a.title.font_size = 5;
a.margins = [0.2,0.1,0.25,0.1];


// fxyz
clf(); //To clear the previous plot
xset('colormap', coppercolormap(128));
f = gcf();
f.figure_size = [620, 620];
funcprot(0); //nothing special is done when a function is redefined
function [r] = f(t, p)
    r = (SphericalHarmonicY(3, -3, t, p) + SphericalHarmonicY(3, 3, t, p))*%i/sqrt(2);
endfunction
//t = linspace(0.0001, %pi-0.0001, 100);
t = linspace(%pi/2, %pi-0.0001, 50);
p = linspace(0, 2*%pi, 200);
r = abs(f(t, p))*10;
x = r .* (sin(t)' * cos(p));
y = r .* (sin(t)' * sin(p));
z = r .* (cos(t)' * ones(p));
[xf, yf, zf] = nf3d(x', y', z'); 
plot3d1(xf, yf, zf, theta = 300, alpha = 60, leg = "10x@10y@10z", flag = [-1 5 4], ..
ebox = [-0.5,0.5,-0.5,0.5,-0.12,0.12]);
xtitle('$\left|f_{y\left(3x^{2}-y^{2}\right)}\right|=\left|i\sqrt{1/2}\left(Y_{3}^{-3}+Y_{3}^{3}\right)\right|$');
a = gca(); //get the current axes
a.title.font_size = 5;
a.margins = [0.2,0.1,0.25,0.1];


// Real spherical harmonics with l=4 -------------------------------------------


// gz4
clf(); //To clear the previous plot
xset('colormap', coppercolormap(128));
f = gcf();
f.figure_size = [620, 620];
funcprot(0); //nothing special is done when a function is redefined
function [r] = f(t, p)
    r = SphericalHarmonicY(4, 0, t, p);
endfunction
t = linspace(0.0001, %pi-0.0001, 100);
//p = linspace(0, 2*%pi, 200);
p = linspace(0, %pi, 100);
r = abs(f(t, p));
x = r .* (sin(t)' * cos(p));
y = r .* (sin(t)' * sin(p));
z = r .* (cos(t)' * ones(p));
[xf, yf, zf] = nf3d(x', y', z'); 
plot3d1(xf, yf, zf, theta = 300, alpha = 60, leg = "x@y@z", flag = [-1 6 4]);
xtitle('$\left|g_{z^{4}}\right|=\left|Y_{4}^{0}\right|$');
a = gca(); //get the current axes
a.title.font_size = 5;
a.margins = [0.1,0.1,0.2,0.1];


// gz3x
clf(); //To clear the previous plot
xset('colormap', coppercolormap(128));
f = gcf();
f.figure_size = [620, 620];
funcprot(0); //nothing special is done when a function is redefined
function [r] = f(t, p)
    r = (SphericalHarmonicY(4, -1, t, p) - SphericalHarmonicY(4, 1, t, p))/sqrt(2);
endfunction
t = linspace(0.0001, %pi-0.0001, 100);
//p = linspace(0, 2*%pi, 200);
p = linspace(0, %pi, 100);
r = abs(f(t, p));
x = r .* (sin(t)' * cos(p));
y = r .* (sin(t)' * sin(p));
z = r .* (cos(t)' * ones(p));
[xf, yf, zf] = nf3d(x', y', z'); 
plot3d1(xf, yf, zf, theta = 300, alpha = 60, leg = "x@y@z", flag = [-1 5 4], ..
ebox = [-0.25,0.25,-0.25,0.25,-0.35,0.35]);
xtitle('$\left|g_{z^{3}x}\right|=\left|\sqrt{1/2}\left(Y_{4}^{-1}-Y_{4}^{1}\right)\right|$');
a = gca(); //get the current axes
a.title.font_size = 5;
a.margins = [0.1,0.1,0.25,0.1];


// gz3y
clf(); //To clear the previous plot
xset('colormap', coppercolormap(128));
f = gcf();
f.figure_size = [620, 620];
funcprot(0); //nothing special is done when a function is redefined
function [r] = f(t, p)
    r = (SphericalHarmonicY(4, -1, t, p) + SphericalHarmonicY(4, 1, t, p))*%i/sqrt(2);
endfunction
t = linspace(0.0001, %pi-0.0001, 100);
//p = linspace(0, 2*%pi, 200);
p = linspace(%pi/2, 3*%pi/2, 100);
r = abs(f(t, p));
x = r .* (sin(t)' * cos(p));
y = r .* (sin(t)' * sin(p));
z = r .* (cos(t)' * ones(p));
[xf, yf, zf] = nf3d(x', y', z'); 
plot3d1(xf, yf, zf, theta = 300, alpha = 60, leg = "x@y@z", flag = [-1 5 4], ..
ebox = [-0.25,0.25,-0.25,0.25,-0.35,0.35]);
xtitle('$\left|g_{z^{3}y}\right|=\left|i\sqrt{1/2}\left(Y_{4}^{-1}+Y_{4}^{1}\right)\right|$');
a = gca(); //get the current axes
a.title.font_size = 5;
a.margins = [0.1,0.1,0.25,0.1];


// gz2xy
clf(); //To clear the previous plot
xset('colormap', coppercolormap(128));
f = gcf();
f.figure_size = [620, 620];
funcprot(0); //nothing special is done when a function is redefined
function [r] = f(t, p)
    r = (SphericalHarmonicY(4, -2, t, p) - SphericalHarmonicY(4, 2, t, p))*%i/sqrt(2);
endfunction
t = linspace(0.0001, %pi-0.0001, 100);
p = linspace(0, 2*%pi, 200);
r = abs(f(t, p))*10;
x = r .* (sin(t)' * cos(p));
y = r .* (sin(t)' * sin(p));
z = r .* (cos(t)' * ones(p));
[xf, yf, zf] = nf3d(x', y', z'); 
plot3d1(xf, yf, zf, theta = 300, alpha = 60, leg = "x@y@z", flag = [-1 5 4], ..
ebox = [-0.6,0.6,-0.6,0.6,-0.6,0.6]);
xtitle('$\left|g_{z^{2}xy}\right|=\left|i\sqrt{1/2}\left(Y_{4}^{-2}-Y_{4}^{2}\right)\right|$');
a = gca(); //get the current axes
a.title.font_size = 5;
a.margins = [0.1,0.1,0.25,0.1];


// gz2(x2-y2)
clf(); //To clear the previous plot
xset('colormap', coppercolormap(128));
f = gcf();
f.figure_size = [620, 620];
funcprot(0); //nothing special is done when a function is redefined
function [r] = f(t, p)
    r = (SphericalHarmonicY(4, -2, t, p) + SphericalHarmonicY(4, 2, t, p))/sqrt(2);
endfunction
t = linspace(0.0001, %pi-0.0001, 100);
p = linspace(0, 2*%pi, 200);
r = abs(f(t, p))*10;
x = r .* (sin(t)' * cos(p));
y = r .* (sin(t)' * sin(p));
z = r .* (cos(t)' * ones(p));
[xf, yf, zf] = nf3d(x', y', z'); 
plot3d1(xf, yf, zf, theta = 300, alpha = 60, leg = "x@y@z", flag = [-1 5 4], ..
ebox = [-0.2,0.2,-0.2,0.2,-0.1,0.1]);
xtitle('$\left|g_{z^{2}\left(x^{2}-y^{2}\right)}\right|=\left|\sqrt{1/2}\left(Y_{4}^{-2}+Y_{4}^{2}\right)\right|$');
a = gca(); //get the current axes
a.title.font_size = 5;
a.margins = [0.1,0.1,0.25,0.1];


// gzx3
clf(); //To clear the previous plot
xset('colormap', coppercolormap(128));
f = gcf();
f.figure_size = [620, 620];
funcprot(0); //nothing special is done when a function is redefined
function [r] = f(t, p)
    r = (SphericalHarmonicY(4, -3, t, p) - SphericalHarmonicY(4, 3, t, p))/sqrt(2);
endfunction
t = linspace(0.0001, %pi-0.0001, 100);
//t = linspace(%pi/2, %pi-0.0001, 50);
p = linspace(0, 2*%pi, 200);
r = abs(f(t, p))*10;
x = r .* (sin(t)' * cos(p));
y = r .* (sin(t)' * sin(p));
z = r .* (cos(t)' * ones(p));
[xf, yf, zf] = nf3d(x', y', z'); 
plot3d1(xf, yf, zf, theta = 300, alpha = 60, leg = "10x@10y@10z", flag = [-1 6 4]);
xtitle('$\left|g_{zx^{3}}\right|=\left|\sqrt{1/2}\left(Y_{4}^{-3}-Y_{4}^{3}\right)\right|$');
a = gca(); //get the current axes
a.title.font_size = 5;
a.margins = [0.1,0.1,0.25,0.1];


// gzy3
clf(); //To clear the previous plot
xset('colormap', coppercolormap(128));
f = gcf();
f.figure_size = [620, 620];
funcprot(0); //nothing special is done when a function is redefined
function [r] = f(t, p)
    r = (SphericalHarmonicY(4, -3, t, p) + SphericalHarmonicY(4, 3, t, p))*%i/sqrt(2);
endfunction
t = linspace(0.0001, %pi-0.0001, 100);
//t = linspace(%pi/2, %pi-0.0001, 50);
p = linspace(0, 2*%pi, 200);
r = abs(f(t, p))*10;
x = r .* (sin(t)' * cos(p));
y = r .* (sin(t)' * sin(p));
z = r .* (cos(t)' * ones(p));
[xf, yf, zf] = nf3d(x', y', z'); 
plot3d1(xf, yf, zf, theta = 300, alpha = 60, leg = "10x@10y@10z", flag = [-1 5 4], ..
ebox = [-0.2,0.2,-0.2,0.2,-0.1,0.1]);
xtitle('$\left|g_{zy^{3}}\right|=\left|i\sqrt{1/2}\left(Y_{4}^{-3}+Y_{4}^{3}\right)\right|$');
a = gca(); //get the current axes
a.title.font_size = 5;
a.margins = [0.1,0.1,0.25,0.1];


// gxy(x2-y2)
clf(); //To clear the previous plot
xset('colormap', coppercolormap(128));
f = gcf();
f.figure_size = [620, 620];
funcprot(0); //nothing special is done when a function is redefined
function [r] = f(t, p)
    r = (SphericalHarmonicY(4, -4, t, p) - SphericalHarmonicY(4, 4, t, p))/sqrt(2);
endfunction
//t = linspace(0.0001, %pi-0.0001, 100);
t = linspace(%pi/2, %pi-0.0001, 50);
p = linspace(0, 2*%pi, 200);
r = abs(f(t, p))*100;
x = r .* (sin(t)' * cos(p));
y = r .* (sin(t)' * sin(p));
z = r .* (cos(t)' * ones(p));
[xf, yf, zf] = nf3d(x', y', z'); 
plot3d1(xf, yf, zf, theta = 300, alpha = 60, leg = "100x@100y@100z", flag = [-1 5 4], ..
ebox = [-0.8,0.8,-0.8,0.8,-0.2,0.2]);
xtitle('$\left|g_{xy\left(x^{2}-y^{2}\right)}\right|=\left|\sqrt{1/2}\left(Y_{4}^{-4}-Y_{4}^{4}\right)\right|$');
a = gca(); //get the current axes
a.title.font_size = 5;
a.margins = [0.2,0.1,0.25,0.1];


// gx4+y4
clf(); //To clear the previous plot
xset('colormap', coppercolormap(128));
f = gcf();
f.figure_size = [620, 620];
funcprot(0); //nothing special is done when a function is redefined
function [r] = f(t, p)
    r = (SphericalHarmonicY(4, -4, t, p) + SphericalHarmonicY(4, 4, t, p))*%i/sqrt(2);
endfunction
//t = linspace(0.0001, %pi-0.0001, 100);
t = linspace(%pi/2, %pi-0.0001, 50);
p = linspace(0, 2*%pi, 200);
r = abs(f(t, p))*100;
x = r .* (sin(t)' * cos(p));
y = r .* (sin(t)' * sin(p));
z = r .* (cos(t)' * ones(p));
[xf, yf, zf] = nf3d(x', y', z'); 
plot3d1(xf, yf, zf, theta = 300, alpha = 60, leg = "100x@100y@100z", flag = [-1 5 4], ..
ebox = [-0.8,0.8,-0.8,0.8,-0.2,0.2]);
xtitle('$\left|g_{x^{4}+y^{4}}\right|=\left|i\sqrt{1/2}\left(Y_{4}^{-4}+Y_{4}^{4}\right)\right|$');
a = gca(); //get the current axes
a.title.font_size = 5;
a.margins = [0.2,0.1,0.25,0.1];


// -----------------------------------------------------------------------------


// Spherical harmonics with l=0 m=0
clf(); //To clear the previous plot
xset('colormap', pinkcolormap(32));
cf = gcf();
cf.figure_size = [1240, 680];
ColormapSize = size(cf.color_map, "r"); //get the colormap size

t = linspace(0.0001, %pi-0.0001, 50);
p = linspace(0, 2*%pi, 100);

subplot(1,3,1)
r = SphericalPlot3D(abs(SphericalHarmonicY(0, 0, t, p)));
xtitle('$\left|Y_{0}^{0}\left(\theta,\varphi\right)\right|$');
a = gca(); //get the current axes
a.title.font_size = 5;

subplot(1,3,2)
r = SphericalPlot3D(real(SphericalHarmonicY(0, 0, t, p)));
xtitle('$Re\, Y_{0}^{0}\left(\theta,\varphi\right)$');
a = gca(); //get the current axes
a.title.font_size = 5;


// Spherical harmonics with l=1 m=0
clf(); //To clear the previous plot
xset('colormap', pinkcolormap(32));
cf = gcf();
cf.figure_size = [1240, 680];
ColormapSize = size(cf.color_map, "r"); //get the colormap size

t = linspace(0.0001, %pi-0.0001, 50);
p = linspace(0, 2*%pi, 100);

subplot(1,3,1)
r = SphericalPlot3D(abs(SphericalHarmonicY(1, 0, t, p)));
xtitle('$\left|Y_{1}^{0}\left(\theta,\varphi\right)\right|$');
a = gca(); //get the current axes
a.title.font_size = 5;

subplot(1,3,2)
r = SphericalPlot3D(real(SphericalHarmonicY(1, 0, t, p)));
xtitle('$Re\, Y_{1}^{0}\left(\theta,\varphi\right)$');
a = gca(); //get the current axes
a.title.font_size = 5;


// Spherical harmonics with l=1 m=1
clf(); //To clear the previous plot
xset('colormap', pinkcolormap(32));
cf = gcf();
cf.figure_size = [1240, 680];
ColormapSize = size(cf.color_map, "r"); //get the colormap size

t = linspace(0.0001, %pi-0.0001, 50);
p = linspace(0, 2*%pi, 100);

subplot(1,3,1)
r = SphericalPlot3D(abs(SphericalHarmonicY(1, 1, t, p)));
xtitle('$\left|Y_{1}^{1}\left(\theta,\varphi\right)\right|$');
a = gca(); //get the current axes
a.data_bounds = [-0.3,-0.3,-0.15;0.3,0.3,0.15];
a.title.font_size = 5;
a.margins = [0.15,0.1,0.2,0.1];

subplot(1,3,2)
r = SphericalPlot3D(real(SphericalHarmonicY(1, 1, t, p)));
xtitle('$Re\, Y_{1}^{1}\left(\theta,\varphi\right)$');
a = gca(); //get the current axes
a.data_bounds = [-0.3,-0.3,-0.15;0.3,0.3,0.15];
a.title.font_size = 5;
a.margins = [0.15,0.1,0.2,0.1];

subplot(1,3,3)
r = SphericalPlot3D(imag(SphericalHarmonicY(1, 1, t, p)));
xtitle('$Im\, Y_{1}^{1}\left(\theta,\varphi\right)$');
a = gca(); //get the current axes
a.data_bounds = [-0.3,-0.3,-0.15;0.3,0.3,0.15];
a.title.font_size = 5;
a.margins = [0.15,0.1,0.2,0.1];



// Spherical harmonics with l=2 m=0
clf(); //To clear the previous plot
xset('colormap', pinkcolormap(32));
cf = gcf();
cf.figure_size = [1240, 680];
ColormapSize = size(cf.color_map, "r"); //get the colormap size

t = linspace(0.0001, %pi-0.0001, 50);
p = linspace(0, 2*%pi, 100);

subplot(1,3,1)
r = SphericalPlot3D(abs(SphericalHarmonicY(2, 0, t, p)));
xtitle('$\left|Y_{2}^{0}\left(\theta,\varphi\right)\right|$');
a = gca(); //get the current axes
a.title.font_size = 5;

subplot(1,3,2)
r = SphericalPlot3D(real(SphericalHarmonicY(2, 0, t, p)));
xtitle('$Re\, Y_{2}^{0}\left(\theta,\varphi\right)$');
a = gca(); //get the current axes
a.title.font_size = 5;


// Spherical harmonics with l=2 m=1
clf(); //To clear the previous plot
xset('colormap', pinkcolormap(32));
cf = gcf();
cf.figure_size = [1240, 680];
ColormapSize = size(cf.color_map, "r"); //get the colormap size

t = linspace(0.0001, %pi-0.0001, 50);
p = linspace(0, 2*%pi, 100);

subplot(1,3,1)
r = SphericalPlot3D(abs(SphericalHarmonicY(2, 1, t, p)));
xtitle('$\left|Y_{2}^{1}\left(\theta,\varphi\right)\right|$');
a = gca(); //get the current axes
a.data_bounds = [-0.2,-0.2,-0.2;0.2,0.2,0.2];
a.title.font_size = 5;
a.margins = [0,0,0.15,0];

subplot(1,3,2)
r = SphericalPlot3D(real(SphericalHarmonicY(2, 1, t, p)));
xtitle('$Re\, Y_{2}^{1}\left(\theta,\varphi\right)$');
a = gca(); //get the current axes
a.data_bounds = [-0.2,-0.2,-0.2;0.2,0.2,0.2];
a.title.font_size = 5;
a.margins = [0,0,0.15,0];

subplot(1,3,3)
r = SphericalPlot3D(imag(SphericalHarmonicY(2, 1, t, p)));
xtitle('$Im\, Y_{2}^{1}\left(\theta,\varphi\right)$');
a = gca(); //get the current axes
a.data_bounds = [-0.2,-0.2,-0.2;0.2,0.2,0.2];
a.title.font_size = 5;
a.margins = [0,0,0.15,0];



// Spherical harmonics with l=2 m=2
clf(); //To clear the previous plot
xset('colormap', pinkcolormap(32));
cf = gcf();
cf.figure_size = [1240, 680];
ColormapSize = size(cf.color_map, "r"); //get the colormap size

t = linspace(0.0001, %pi-0.0001, 50);
p = linspace(0, 2*%pi, 100);

subplot(1,3,1)
r = SphericalPlot3D(abs(SphericalHarmonicY(2, 2, t, p))*10);
xtitle('$\left|Y_{2}^{2}\left(\theta,\varphi\right)\right|$');
a = gca(); //get the current axes
a.data_bounds = [-1.5,-1.5,-0.5;1.5,1.5,0.5];
a.title.font_size = 5;
a.x_label.text = "10x";
a.y_label.text = "10y";
a.z_label.text = "10z";
a.margins = [0.15,0.05,0.2,0.1];

subplot(1,3,2)
r = SphericalPlot3D(real(SphericalHarmonicY(2, 2, t, p))*10);
xtitle('$Re\, Y_{2}^{2}\left(\theta,\varphi\right)$');
a = gca(); //get the current axes
a.data_bounds = [-1.5,-1.5,-0.5;1.5,1.5,0.5];
a.title.font_size = 5;
a.x_label.text = "10x";
a.y_label.text = "10y";
a.z_label.text = "10z";
a.margins = [0.15,0.05,0.2,0.1];

subplot(1,3,3)
r = SphericalPlot3D(imag(SphericalHarmonicY(2, 2, t, p))*10);
xtitle('$Im\, Y_{2}^{2}\left(\theta,\varphi\right)$');
a = gca(); //get the current axes
a.data_bounds = [-1.5,-1.5,-0.5;1.5,1.5,0.5];
a.title.font_size = 5;
a.x_label.text = "10x";
a.y_label.text = "10y";
a.z_label.text = "10z";
a.margins = [0.15,0.05,0.2,0.1];


// Spherical harmonics with l=3 m=0
clf(); //To clear the previous plot
xset('colormap', pinkcolormap(32));
cf = gcf();
cf.figure_size = [1240, 680];
ColormapSize = size(cf.color_map, "r"); //get the colormap size

t = linspace(0.0001, %pi-0.0001, 50);
p = linspace(0, 2*%pi, 100);

subplot(1,3,1)
r = SphericalPlot3D(abs(SphericalHarmonicY(3, 0, t, p)));
xtitle('$\left|Y_{3}^{0}\left(\theta,\varphi\right)\right|$');
a = gca(); //get the current axes
a.title.font_size = 5;

subplot(1,3,2)
r = SphericalPlot3D(real(SphericalHarmonicY(3, 0, t, p)));
xtitle('$Re\, Y_{3}^{0}\left(\theta,\varphi\right)$');
a = gca(); //get the current axes
a.title.font_size = 5;


// Spherical harmonics with l=3 m=1
clf(); //To clear the previous plot
xset('colormap', pinkcolormap(32));
cf = gcf();
cf.figure_size = [1240, 680];
ColormapSize = size(cf.color_map, "r"); //get the colormap size

t = linspace(0.0001, %pi-0.0001, 50);
p = linspace(0, 2*%pi, 100);

subplot(1,3,1)
r = SphericalPlot3D(abs(SphericalHarmonicY(3, 1, t, p)));
xtitle('$\left|Y_{3}^{1}\left(\theta,\varphi\right)\right|$');
a = gca(); //get the current axes
a.data_bounds = [-0.2,-0.2,-0.25;0.2,0.2,0.25];
a.title.font_size = 5;
a.margins = [0,0,0.15,0];

subplot(1,3,2)
r = SphericalPlot3D(real(SphericalHarmonicY(3, 1, t, p)));
xtitle('$Re\, Y_{3}^{1}\left(\theta,\varphi\right)$');
a = gca(); //get the current axes
a.data_bounds = [-0.2,-0.2,-0.25;0.2,0.2,0.25];
a.title.font_size = 5;
a.margins = [0,0,0.15,0];

subplot(1,3,3)
r = SphericalPlot3D(imag(SphericalHarmonicY(3, 1, t, p)));
xtitle('$Im\, Y_{3}^{1}\left(\theta,\varphi\right)$');
a = gca(); //get the current axes
a.data_bounds = [-0.2,-0.2,-0.25;0.2,0.2,0.25];
a.title.font_size = 5;
a.margins = [0,0,0.15,0];



// Spherical harmonics with l=3 m=2
clf(); //To clear the previous plot
xset('colormap', pinkcolormap(32));
cf = gcf();
cf.figure_size = [1240, 680];
ColormapSize = size(cf.color_map, "r"); //get the colormap size

t = linspace(0.0001, %pi-0.0001, 50);
p = linspace(0, 2*%pi, 100);

subplot(1,3,1)
r = SphericalPlot3D(abs(SphericalHarmonicY(3, 2, t, p))*10);
xtitle('$\left|Y_{3}^{2}\left(\theta,\varphi\right)\right|$');
a = gca(); //get the current axes
a.data_bounds = [-0.6,-0.6,-0.5;0.6,0.6,0.5];
a.title.font_size = 5;
a.x_label.text = "10x";
a.y_label.text = "10y";
a.z_label.text = "10z";
a.margins = [0.15,0.05,0.2,0.1];

subplot(1,3,2)
r = SphericalPlot3D(real(SphericalHarmonicY(3, 2, t, p))*10);
xtitle('$Re\, Y_{3}^{2}\left(\theta,\varphi\right)$');
a = gca(); //get the current axes
a.data_bounds = [-0.6,-0.6,-0.5;0.6,0.6,0.5];
a.title.font_size = 5;
a.x_label.text = "10x";
a.y_label.text = "10y";
a.z_label.text = "10z";
a.margins = [0.15,0.05,0.2,0.1];

subplot(1,3,3)
r = SphericalPlot3D(imag(SphericalHarmonicY(3, 2, t, p))*10);
xtitle('$Im\, Y_{3}^{2}\left(\theta,\varphi\right)$');
a = gca(); //get the current axes
a.data_bounds = [-0.6,-0.6,-0.5;0.6,0.6,0.5];
a.title.font_size = 5;
a.x_label.text = "10x";
a.y_label.text = "10y";
a.z_label.text = "10z";
a.margins = [0.15,0.05,0.2,0.1];



// Spherical harmonics with l=3 m=3
clf(); //To clear the previous plot
xset('colormap', pinkcolormap(32));
cf = gcf();
cf.figure_size = [1240, 680];
ColormapSize = size(cf.color_map, "r"); //get the colormap size

t = linspace(0.0001, %pi-0.0001, 50);
p = linspace(0, 2*%pi, 100);

subplot(1,3,1)
r = SphericalPlot3D(abs(SphericalHarmonicY(3, 3, t, p))*10);
xtitle('$\left|Y_{3}^{3}\left(\theta,\varphi\right)\right|$');
a = gca(); //get the current axes
a.title.font_size = 5;
a.x_label.text = "10x";
a.y_label.text = "10y";
a.z_label.text = "10z";
a.margins = [0.15,0.05,0.2,0.1];

subplot(1,3,2)
r = SphericalPlot3D(real(SphericalHarmonicY(3, 3, t, p))*10);
xtitle('$Re\, Y_{3}^{3}\left(\theta,\varphi\right)$');
a = gca(); //get the current axes
a.title.font_size = 5;
a.x_label.text = "10x";
a.y_label.text = "10y";
a.z_label.text = "10z";
a.margins = [0.15,0.05,0.2,0.1];

subplot(1,3,3)
r = SphericalPlot3D(imag(SphericalHarmonicY(3, 3, t, p))*10);
xtitle('$Im\, Y_{3}^{3}\left(\theta,\varphi\right)$');
a = gca(); //get the current axes
a.title.font_size = 5;
a.x_label.text = "10x";
a.y_label.text = "10y";
a.z_label.text = "10z";
a.margins = [0.15,0.05,0.2,0.1];


// Spherical harmonics with l=4 m=0
clf(); //To clear the previous plot
xset('colormap', pinkcolormap(32));
cf = gcf();
cf.figure_size = [1240, 680];
ColormapSize = size(cf.color_map, "r"); //get the colormap size

t = linspace(0.0001, %pi-0.0001, 50);
p = linspace(0, 2*%pi, 100);

subplot(1,3,1)
r = SphericalPlot3D(abs(SphericalHarmonicY(4, 0, t, p)));
xtitle('$\left|Y_{4}^{0}\left(\theta,\varphi\right)\right|$');
a = gca(); //get the current axes
a.title.font_size = 5;

subplot(1,3,2)
r = SphericalPlot3D(real(SphericalHarmonicY(4, 0, t, p)));
xtitle('$Re\, Y_{4}^{0}\left(\theta,\varphi\right)$');
a = gca(); //get the current axes
a.title.font_size = 5;


// Spherical harmonics with l=4 m=1
clf(); //To clear the previous plot
xset('colormap', pinkcolormap(32));
cf = gcf();
cf.figure_size = [1240, 680];
ColormapSize = size(cf.color_map, "r"); //get the colormap size

t = linspace(0.0001, %pi-0.0001, 50);
p = linspace(0, 2*%pi, 100);

subplot(1,3,1)
r = SphericalPlot3D(abs(SphericalHarmonicY(4, 1, t, p))*10);
xtitle('$\left|Y_{4}^{1}\left(\theta,\varphi\right)\right|$');
a = gca(); //get the current axes
a.data_bounds = [-1.5,-1.5,-2.5;1.5,1.5,2.5];
a.title.font_size = 5;
a.x_label.text = "10x";
a.y_label.text = "10y";
a.z_label.text = "10z";
a.margins = [0.1,0.1,0.2,0.1];

subplot(1,3,2)
r = SphericalPlot3D(real(SphericalHarmonicY(4, 1, t, p))*10);
xtitle('$Re\, Y_{4}^{1}\left(\theta,\varphi\right)$');
a = gca(); //get the current axes
a.data_bounds = [-1.5,-1.5,-2.5;1.5,1.5,2.5];
a.title.font_size = 5;
a.x_label.text = "10x";
a.y_label.text = "10y";
a.z_label.text = "10z";
a.margins = [0.1,0.1,0.2,0.1];

subplot(1,3,3)
r = SphericalPlot3D(imag(SphericalHarmonicY(4, 1, t, p))*10);
xtitle('$Im\, Y_{4}^{1}\left(\theta,\varphi\right)$');
a = gca(); //get the current axes
a.data_bounds = [-1.5,-1.5,-2.5;1.5,1.5,2.5];
a.title.font_size = 5;
a.x_label.text = "10x";
a.y_label.text = "10y";
a.z_label.text = "10z";
a.margins = [0.1,0.1,0.2,0.1];


// Spherical harmonics with l=4 m=2
clf(); //To clear the previous plot
xset('colormap', pinkcolormap(32));
cf = gcf();
cf.figure_size = [1240, 680];
ColormapSize = size(cf.color_map, "r"); //get the colormap size

t = linspace(0.0001, %pi-0.0001, 50);
p = linspace(0, 2*%pi, 100);

subplot(1,3,1)
r = SphericalPlot3D(abs(SphericalHarmonicY(4, 2, t, p))*10);
xtitle('$\left|Y_{4}^{2}\left(\theta,\varphi\right)\right|$');
a = gca(); //get the current axes
a.data_bounds = [-0.4,-0.4,-0.4;0.4,0.4,0.4];
a.title.font_size = 5;
a.x_label.text = "10x";
a.y_label.text = "10y";
a.z_label.text = "10z";
a.margins = [0.15,0.05,0.2,0.1];

subplot(1,3,2)
r = SphericalPlot3D(real(SphericalHarmonicY(4, 2, t, p))*10);
xtitle('$Re\, Y_{4}^{2}\left(\theta,\varphi\right)$');
a = gca(); //get the current axes
a.data_bounds = [-0.4,-0.4,-0.4;0.4,0.4,0.4];
a.title.font_size = 5;
a.x_label.text = "10x";
a.y_label.text = "10y";
a.z_label.text = "10z";
a.margins = [0.15,0.05,0.2,0.1];

subplot(1,3,3)
r = SphericalPlot3D(imag(SphericalHarmonicY(4, 2, t, p))*10);
xtitle('$Im\, Y_{4}^{2}\left(\theta,\varphi\right)$');
a = gca(); //get the current axes
a.data_bounds = [-0.4,-0.4,-0.4;0.4,0.4,0.4];
a.title.font_size = 5;
a.x_label.text = "10x";
a.y_label.text = "10y";
a.z_label.text = "10z";
a.margins = [0.15,0.05,0.2,0.1];


// Spherical harmonics with l=4 m=3
clf(); //To clear the previous plot
xset('colormap', pinkcolormap(32));
cf = gcf();
cf.figure_size = [1240, 680];
ColormapSize = size(cf.color_map, "r"); //get the colormap size

t = linspace(0.0001, %pi-0.0001, 50);
p = linspace(0, 2*%pi, 100);

subplot(1,3,1)
r = SphericalPlot3D(abs(SphericalHarmonicY(4, 3, t, p))*100);
xtitle('$\left|Y_{4}^{3}\left(\theta,\varphi\right)\right|$');
a = gca(); //get the current axes
a.data_bounds = [-1.5,-1.5,-0.8;1.5,1.5,0.8];
a.title.font_size = 5;
a.x_label.text = "100x";
a.y_label.text = "100y";
a.z_label.text = "100z";
a.margins = [0.15,0.05,0.2,0.1];

subplot(1,3,2)
r = SphericalPlot3D(real(SphericalHarmonicY(4, 3, t, p))*100);
xtitle('$Re\, Y_{4}^{3}\left(\theta,\varphi\right)$');
a = gca(); //get the current axes
a.data_bounds = [-1.5,-1.5,-0.8;1.5,1.5,0.8];
a.title.font_size = 5;
a.x_label.text = "100x";
a.y_label.text = "100y";
a.z_label.text = "100z";
a.margins = [0.15,0.05,0.2,0.1];

subplot(1,3,3)
r = SphericalPlot3D(imag(SphericalHarmonicY(4, 3, t, p))*100);
xtitle('$Im\, Y_{4}^{3}\left(\theta,\varphi\right)$');
a = gca(); //get the current axes
a.data_bounds = [-1.5,-1.5,-0.8;1.5,1.5,0.8];
a.title.font_size = 5;
a.x_label.text = "100x";
a.y_label.text = "100y";
a.z_label.text = "100z";
a.margins = [0.15,0.05,0.2,0.1];


// Spherical harmonics with l=4 m=4
clf(); //To clear the previous plot
xset('colormap', pinkcolormap(32));
cf = gcf();
cf.figure_size = [1240, 680];
ColormapSize = size(cf.color_map, "r"); //get the colormap size

t = linspace(0.0001, %pi-0.0001, 50);
p = linspace(0, 2*%pi, 100);

subplot(1,3,1)
r = SphericalPlot3D(abs(SphericalHarmonicY(4, 4, t, p))*100);
xtitle('$\left|Y_{4}^{4}\left(\theta,\varphi\right)\right|$');
a = gca(); //get the current axes
a.title.font_size = 5;
a.x_label.text = "100x";
a.y_label.text = "100y";
a.z_label.text = "100z";
a.margins = [0.2,0.05,0.2,0.1];

subplot(1,3,2)
r = SphericalPlot3D(real(SphericalHarmonicY(4, 4, t, p))*100);
xtitle('$Re\, Y_{4}^{4}\left(\theta,\varphi\right)$');
a = gca(); //get the current axes
a.title.font_size = 5;
a.x_label.text = "100x";
a.y_label.text = "100y";
a.z_label.text = "100z";
a.margins = [0.2,0.05,0.2,0.1];

subplot(1,3,3)
r = SphericalPlot3D(imag(SphericalHarmonicY(4, 4, t, p))*100);
xtitle('$Im\, Y_{4}^{4}\left(\theta,\varphi\right)$');
a = gca(); //get the current axes
a.title.font_size = 5;
a.x_label.text = "100x";
a.y_label.text = "100y";
a.z_label.text = "100z";
a.margins = [0.2,0.05,0.2,0.1];