SHTns  3.4.6
SHT_example.py

A Python example using NumPy and SHTns that performs backward and forward Spherical Harmonic Transforms.

1 #
2 # Copyright (c) 2010-2018 Centre National de la Recherche Scientifique.
3 # written by Nathanael Schaeffer (CNRS, ISTerre, Grenoble, France).
4 #
5 # nathanael.schaeffer@univ-grenoble-alpes.fr
6 #
7 # This software is governed by the CeCILL license under French law and
8 # abiding by the rules of distribution of free software. You can use,
9 # modify and/or redistribute the software under the terms of the CeCILL
10 # license as circulated by CEA, CNRS and INRIA at the following URL
11 # "http://www.cecill.info".
12 #
13 # The fact that you are presently reading this means that you have had
14 # knowledge of the CeCILL license and that you accept its terms.
15 #
16 
17 
20 
21 import numpy # numpy for arrays
22 import shtns # shtns module compiled and installed using
23  # ./configure --enable-python && make && make install
24 
25 lmax = 7 # maximum degree of spherical harmonic representation.
26 mmax = 3 # maximum order of spherical harmonic representation.
27 
28 sh = shtns.sht(lmax, mmax) # create sht object with given lmax and mmax (orthonormalized)
29 # sh = shtns.sht(lmax, mmax, mres=2, norm=shtns.sht_schmidt | shtns.SHT_NO_CS_PHASE) # use schmidt semi-normalized harmonics
30 
31 nlat, nphi = sh.set_grid() # build default grid (gauss grid, phi-contiguous)
32 print(sh.nlat, sh.nphi) # displays the latitudinal and longitudinal grid sizes.
33 
34 cost = sh.cos_theta # latitudinal coordinates of the grid as cos(theta)
35 el = sh.l # array of size sh.nlm giving the spherical harmonic degree l for any sh coefficient
36 l2 = el*(el+1) # array l(l+1) that is useful for computing laplacian
37 
38 
44 
45 ylm = sh.spec_array() # a spherical harmonic spectral array, same as numpy.zeros(sh.nlm, dtype=complex)
46 vr = sh.spat_array() # a spatial array, same as numpy.zeros(sh.spat_shape)
47 
48 ylm[sh.idx(1, 0)] = 1.0 # set sh coefficient l=1, m=0 to value 1
49 
50 print(ylm[sh.l == 1]) # print all l=1 coefficients
51 print(ylm[sh.m == 1]) # print all m=1 coefficients
52 
53 ylm = ylm * l2 # multiply by l(l+1)
54 
55 y = sh.synth(ylm) # transform sh description ylm into spatial representation y (scalar transform)
56 
57 print(y) # display spatial field
58 
59 i_theta = sh.nlat//2
60 i_phi = 1
61 print(y[i_theta, i_phi]) # spatial element of coordinate i_theta, i_phi
62 
63 
64 zlm = sh.analys(y) # transform the spatial field back to spectral
65 print(zlm)
66 
67 
68 dy_dt, dy_dp = sh.synth_grad(ylm) # compute gradients of ylm on a spatial grid.
69 
70