From 85d45cb2a89143fee1f4ac0e188c574b963645c1 Mon Sep 17 00:00:00 2001 From: Nao Pross Date: Sun, 17 Apr 2022 23:55:32 +0200 Subject: Add existing code --- manim/harmonics.py | 67 ++++++++++++++++++++++++++++++++++++++++++++++++++ sph_harm.py | 71 +++++++++++++++++++++++++++++++++++++++++++++++++++++ sph_harm_sandbox.py | 32 ++++++++++++++++++++++++ 3 files changed, 170 insertions(+) create mode 100644 manim/harmonics.py create mode 100644 sph_harm.py create mode 100644 sph_harm_sandbox.py 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 -- cgit v1.2.1