summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
-rw-r--r--manim/harmonics.py67
-rw-r--r--sph_harm.py71
-rw-r--r--sph_harm_sandbox.py32
3 files changed, 170 insertions, 0 deletions
diff --git a/manim/harmonics.py b/manim/harmonics.py
new file mode 100644
index 0000000..a908e25
--- /dev/null
+++ b/manim/harmonics.py
@@ -0,0 +1,67 @@
+from manim import *
+from manim.mobject.opengl.opengl_surface import OpenGLSurface
+from manim.mobject.opengl.opengl_three_dimensions import OpenGLSurfaceMesh
+
+
+import numpy as np
+from scipy.special import sph_harm
+
+# configure style
+config.background_color = '#202020'
+config.tex_template.add_to_preamble(
+ r"\usepackage[p,osf]{scholax}"
+ r"\usepackage{amsmath}"
+ r"\usepackage[scaled=1.075,ncf,vvarbb]{newtxmath}"
+)
+
+class FourierOnS1(ThreeDScene):
+ @staticmethod
+ def harmonic_cos(m, n, theta, phi):
+ # get only the cos() part
+ r = np.real(sph_harm(m, n, theta, phi)) * 5 # Y_n^m
+
+ # coordinates transformation
+ x = r * np.cos(theta) * np.sin(phi)
+ y = r * np.sin(phi) * np.sin(theta)
+ z = r * np.cos(phi)
+
+ return [x, y, z]
+
+ @staticmethod
+ def harmonic_surf(m, n):
+ surf = Surface(lambda u, v: FourierOnS1.harmonic_cos(m, n, u, v),
+ u_range=[0, 2 * PI], v_range=[-PI, PI], resolution=[50, 50],
+ fill_opacity=.8, checkerboard_colors=[BLUE_D, BLUE_E])
+
+ return surf
+
+ def construct(self):
+ axes = ThreeDAxes(x_range=[-1, 1], y_range=[-1, 1])
+ self.add(axes)
+
+ self.set_camera_orientation(theta=30 * DEGREES, phi=75 * DEGREES)
+ self.begin_ambient_camera_rotation(rate=.1)
+
+ last_surf = None
+ last_label = None
+ for n in range(0, 5):
+ for m in range(0, n):
+ surf = FourierOnS1.harmonic_surf(m, n)
+ label = MathTex(r"Y_{%d}^{%d}(\theta, \phi)" % (n, m))
+
+ if last_surf is not None:
+ self.play(ReplacementTransform(last_surf, surf))
+
+ if last_label is not None:
+ self.remove(last_label)
+
+ self.play(Create(surf))
+ self.add_fixed_in_frame_mobjects(label)
+ label.to_corner(UL)
+
+ last_surf = surf
+ last_label = label
+
+ self.wait()
+
+ self.wait(5)
diff --git a/sph_harm.py b/sph_harm.py
new file mode 100644
index 0000000..1d6227b
--- /dev/null
+++ b/sph_harm.py
@@ -0,0 +1,71 @@
+# -*- coding: utf-8 -*-
+"""
+date of creation: Fri Mar 25 16:20:46 2022
+-------------------------------------------------
+@author : Manuel Cattaneo
+@contact: cattaneo.manuel@hotmail.com
+-------------------------------------------------
+
+Description:
+ Program used to play around with the spherical harmonics.
+
+Modules:
+"""
+import numpy as np
+# plot stuff
+import matplotlib.pyplot as plt
+import mpl_toolkits.mplot3d.axes3d as axes3d # to plot 3d
+import matplotlib.gridspec as gs # for grid layout
+# spherical harmonics
+from scipy.special import sph_harm
+"""
+------------------------------------------------------------------------------
+"""
+
+if __name__ == '__main__':
+ # grid spec
+ row = 3
+ col = 2*row-1
+ grid_spec = gs.GridSpec(row,col)
+ fig = plt.figure()
+
+ # creating spherical domain
+ N = 50
+ theta, phi = np.meshgrid(np.linspace(0, 2*np.pi, N), # theta in [0, 2pi]
+ np.linspace(0, np.pi, N)) # phi in [0, pi]
+
+ import matplotlib
+ # m < n
+ for y_grid in range(0,row):
+ for i,x_grid in enumerate(range(col//2-y_grid, col-col//2+y_grid)):
+ # real(...) to obtain only the cos(.) part. This is done in order
+ # to have the same functions as the ones at p.518 chapter 12.2.
+ m = y_grid
+ n = m+i+1
+ r = np.real( sph_harm(m, n, theta, phi) ) # Y_n^m
+ # coordinates transformation
+ x = r*np.cos(theta)*np.sin(phi)
+ y = r*np.sin(phi) *np.sin(theta)
+ z = r *np.cos(phi)
+
+ # adding subplot
+ ax = fig.add_subplot(grid_spec[y_grid,x_grid], projection='3d')
+ ax.set_axis_off()
+ ax.set_title(rf'$Y_{n}^{m}(\theta, \phi)$')
+
+ # colors
+ color_dim = r
+ color_min, color_max = color_dim.min(), color_dim.max()
+ norm = matplotlib.colors.Normalize(color_min,color_max)
+ m = plt.cm.ScalarMappable(norm, cmap='jet')
+ m.set_array([])
+ fc = m.to_rgba(color_dim)
+
+ # plot
+ plot = ax.plot_surface(x, y, z, facecolors = fc,
+ rstride=1, cstride=1,
+ # cmap=plt.get_cmap('jet'),
+ linewidth=0, antialiased=False, alpha=0.5)
+
+
+ plt.show() \ No newline at end of file
diff --git a/sph_harm_sandbox.py b/sph_harm_sandbox.py
new file mode 100644
index 0000000..51828bc
--- /dev/null
+++ b/sph_harm_sandbox.py
@@ -0,0 +1,32 @@
+# -*- coding: utf-8 -*-
+"""
+date of creation: Fri Mar 10 03:20:46 2022
+-------------------------------------------------
+@author : Manuel Cattaneo
+@contact: cattaneo.manuel@hotmail.com
+-------------------------------------------------
+
+Description:
+ Programma per capire come cazzo funziona sta SHExpandDHC(.) di merda.
+
+Modules:
+"""
+from scipy.special import sph_harm
+from pyshtools.expand import SHExpandDHC, SHExpandDH
+
+import numpy as np
+"""
+------------------------------------------------------------------------------
+"""
+
+_theta = np.linspace(0, np.pi , 1001)[1:] # theta in ]0, pi]
+_phi = np.linspace(0, 2*np.pi, 1001)[:-1] # phi in [0, 2pi[
+
+theta, phi = np.meshgrid(_theta, _phi)
+
+# n >! m
+m = 2
+n = 4
+y = sph_harm(m, n, theta=phi, phi=theta) # angles inverted due to convention
+
+cilm = np.abs(SHExpandDHC(y)) \ No newline at end of file