2023-9-28

expansion.py

@@ -1,9 +1,16 @@
+from __future__ import annotations
 import numpy as np
 from dataclasses import dataclass
 from functools import cached_property
 import functions as fn
+import quaternionic
+import spherical
+import time
+import copy
+
 
 Array = np.ndarray
+Quaternion = quaternionic.array
 
 
 class InvalidExpansion(Exception):
@@ -15,14 +22,17 @@ class Expansion:
     """Generic class for storing surface charge expansion coefficients."""
 
     l_array: Array
-    coeffs: Array
+    coefs: Array
+    _starting_coefs: Array = None  # initialized with the __post_init__ method
+    _rotations: Quaternion = Quaternion([1., 0., 0., 0.])
 
     def __post_init__(self):
-        """Validation of the given expansion."""
-        if len(self.coeffs) != np.sum(2 * self.l_array + 1):
+        if self.coefs.shape[-1] != np.sum(2 * self.l_array + 1):
             raise InvalidExpansion('Number of expansion coefficients does not match the provided l_array.')
         if np.all(np.sort(self.l_array) != self.l_array) or np.all(np.unique(self.l_array) != self.l_array):
             raise InvalidExpansion('Array of l values should be unique and sorted.')
+        self.coefs = self.coefs.astype(np.complex128)
+        self._starting_coefs = np.copy(self.coefs)
 
     @property
     def max_l(self) -> int:
@@ -37,17 +47,22 @@ class Expansion:
                 all_m_list.append(-l + i)
         return np.repeat(self.l_array, 2 * self.l_array + 1), np.array(all_m_list)
 
-    def repeat_over_m(self, arr: Array, axis=0):
+    def repeat_over_m(self, arr: Array, axis=0) -> Array:
         if not arr.shape[axis] == len(self.l_array):
             raise ValueError('Array length should be equal to the number of l in the expansion.')
         return np.repeat(arr, 2 * self.l_array + 1, axis=axis)
 
+    def rotate(self, rotations: Quaternion, rotate_existing=False):
+        self._rotations = rotations
+        coefs = self.coefs if rotate_existing else self._starting_coefs
+        self.coefs = expansion_rotation(rotations, coefs, self.l_array)
 
-def rot_sym_expansion(l_array: Array, coeffs: Array) -> Array:
-    """Create full expansion array for rotationally symmetric distributions with only m=0 terms different form 0."""
-    full_coeffs = np.zeros(np.sum(2 * l_array + 1))
-    full_coeffs[np.cumsum(2 * l_array + 1) - l_array - 1] = coeffs
-    return full_coeffs
+    def rotate_euler(self, alpha: Array, beta: Array, gamma: Array, rotate_existing=False):
+        R_euler = quaternionic.array.from_euler_angles(alpha, beta, gamma)
+        self.rotate(R_euler, rotate_existing=rotate_existing)
+
+    def clone(self) -> Expansion:
+        return copy.deepcopy(self)
 
 
 class Expansion24(Expansion):
@@ -66,13 +81,72 @@ class MappedExpansion(Expansion):
                   * np.sqrt(4 * np.pi * (2 * l_array + 1)) * np.power(a_bar, l_array))
         coeffs[0] = sigma0
         coeffs = rot_sym_expansion(l_array, coeffs)
+        # coeffs = np.full((1000, len(coeffs)), coeffs)
         super().__init__(l_array, coeffs)
 
 
+def rot_sym_expansion(l_array: Array, coeffs: Array) -> Array:
+    """Create full expansion array for rotationally symmetric distributions with only m=0 terms different form 0."""
+    full_coeffs = np.zeros(np.sum(2 * l_array + 1))
+    full_coeffs[np.cumsum(2 * l_array + 1) - l_array - 1] = coeffs
+    return full_coeffs
+
+
+def coeffs_fill_missing_l(expansion: Expansion, target_l_array: Array) -> Expansion:
+    missing_l = np.setdiff1d(target_l_array, expansion.l_array, assume_unique=True)
+    fill = np.zeros(np.sum(2 * missing_l + 1))
+    full_l_array1, _ = expansion.lm_arrays
+    # we search for where to place missing coeffs with the help of a boolean array and argmax function
+    comparison_bool = (full_l_array1[:, None] - missing_l[None, :]) > 0
+    indices = np.where(np.any(comparison_bool, axis=0), np.argmax(comparison_bool, axis=0), full_l_array1.shape[0])
+    new_coeffs = np.insert(expansion.coefs, np.repeat(indices, 2 * missing_l + 1), fill, axis=-1)
+    return Expansion(target_l_array, new_coeffs)
+
+
+def expansions_to_common_l(ex1: Expansion, ex2: Expansion) -> (Expansion, Expansion):
+    common_l_array = np.union1d(ex1.l_array, ex2.l_array)
+    return coeffs_fill_missing_l(ex1, common_l_array),  coeffs_fill_missing_l(ex2, common_l_array)
+
+
+def expansion_rotation(rotations: Quaternion, coefs: Array, l_array: Array):
+    """
+    General function for rotations of expansion coefficients using WignerD matrices. Combines all rotations
+    with each expansion given in coefs array.
+    :param rotations: Quaternion array, last dimension is 4
+    :param coefs: array of expansion coefficients
+    :param l_array: array of all ell values of the expansion
+    :return rotated coefficients, output shape is rotations.shape[:-1] + coefs.shape
+    """
+    rot_arrays = rotations.ndarray.reshape((-1, 4))
+    coefs_reshaped = coefs.reshape((-1, coefs.shape[-1]))
+    wigner_matrices = spherical.Wigner(np.max(l_array)).D(rot_arrays)
+    new_coefs = np.zeros((rot_arrays.shape[0],) + coefs_reshaped.shape, dtype=np.complex128)
+    for i, l in enumerate(l_array):
+        Dlmn_slice = np.arange(l * (2 * l - 1) * (2 * l + 1) / 3, (l + 1) * (2 * l + 1) * (2 * l + 3) / 3).astype(int)
+        all_m_indices = np.arange(np.sum(2 * l_array[:i] + 1), np.sum(2 * l_array[:i + 1] + 1))
+        wm = wigner_matrices[:, Dlmn_slice].reshape((-1, 2*l+1, 2*l+1))
+        new_coefs[..., all_m_indices] = np.einsum('rmn, qm -> rqn',
+                                                   wm, coefs_reshaped[:, all_m_indices])
+    return new_coefs.reshape(rotations.ndarray.shape[:-1] + coefs.shape)
+
+
 if __name__ == '__main__':
 
     ex = MappedExpansion(0.44, 3, 1, 10)
-    print(ex.coeffs)
+    print(ex.coefs)
+
+    # new_coeffs = expansion_rotation(Quaternion(np.arange(20).reshape(5, 4)).normalized, ex.coeffs, ex.l_array)
+    # print(new_coeffs.shape)
+    #
+    # newnew_coeffs = expansion_rotation(Quaternion(np.arange(16).reshape(4, 4)).normalized, new_coeffs, ex.l_array)
+    # print(newnew_coeffs.shape)
+
+    rot_angles = np.linspace(0, np.pi, 1000)
+    t0 = time.time()
+    ex.rotate_euler(rot_angles, rot_angles, rot_angles)
+    t1 = time.time()
+    print(ex.coefs.shape)
+    print(t1 - t0)
 
 
 

functions.py

@@ -2,9 +2,6 @@ import numpy as np
 import scipy.special as sps
 
 
-Array = np.ndarray
-
-
 def sph_bessel_i(n, x, **kwargs):
     return sps.spherical_in(n, x, **kwargs)
 
@@ -21,14 +18,9 @@ def interaction_coeff_C(l, p, x):
     return x * sps.iv(l + 1 / 2, x) * sps.iv(p + 1 / 2, x)
 
 
-def coefficient_Cpm(l, x) -> Array:
+def coefficient_Cpm(l, x):
     return x * sps.kv(l + 1 / 2, x) * sps.iv(l + 1 / 2, x)
 
 
 def coeff_C_diff(l, x):
     return 1 / (x * sps.iv(l + 1 / 2, x) * sps.kv(l + 3 / 2, x))
-
-
-if __name__ == '__main__':
-    print(sph_bessel_k(np.array([1, 2, 3])[:, None], np.array([0.1, 0.2, 0.5, 0.6])[None, :]))
-

interactions.py

@@ -1,71 +1,49 @@
 import expansion
 import parameters
 import functions as fn
-import numpy as np
 from py3nj import wigner3j
+import numpy as np
 import time
-import matplotlib.pyplot as plt
+from typing import Literal
+import units_and_constants as uc
 
 
 Array = np.ndarray
 Expansion = expansion.Expansion
 ModelParams = parameters.ModelParams
+EnergyUnit = Literal['kT', 'eV', 'J']
 
 
-# def prefactor(R: float, kappaR: float, c0: float):
-#     return R * kappaR * 1e4 / (12.04 * c0)
-#
-#
-# def c0(R: float, kappaR: float):
-#     return 10 * kappaR ** 2 / (0.329 ** 2 * R ** 2)
-
-
-def expansions_to_common_l(ex1: Expansion, ex2: expansion) -> (Expansion, Expansion):
-    common_l_array = np.union1d(ex1.l_array, ex2.l_array)
-    missing_l1 = np.setdiff1d(common_l_array, ex1.l_array, assume_unique=True)
-    missing_l2 = np.setdiff1d(common_l_array, ex2.l_array, assume_unique=True)
-
-    fill_1 = np.zeros(np.sum(2 * missing_l1 + 1))
-    fill_2 = np.zeros(np.sum(2 * missing_l2 + 1))
-
-    full_l_array1, _ = ex1.lm_arrays
-    full_l_array2, _ = ex2.lm_arrays
-
-    # we search for where to place missing coeffs with the help of a boolean array and argmax function
-    bool1 = (full_l_array1[:, None] - missing_l1[None, :]) > 0
-    bool2 = (full_l_array2[:, None] - missing_l2[None, :]) > 0
-
-    # if all elements of bool are false, i.e. missing_l > all existing_l, additional l values come to the end
-    indices1 = np.where(np.any(bool1, axis=0), np.argmax(bool1, axis=0), full_l_array1.shape[0])
-    indices2 = np.where(np.any(bool2, axis=0), np.argmax(bool2, axis=0), full_l_array2.shape[0])
-
-    new_coeffs1 = np.insert(ex1.coeffs, np.repeat(indices1, 2 * missing_l1 + 1), fill_1)
-    new_coeffs2 = np.insert(ex2.coeffs, np.repeat(indices2, 2 * missing_l2 + 1), fill_2)
+def energy_units(units: EnergyUnit, params: ModelParams) -> float:
+    match units:
+        case 'eV':
+            return 1 / (uc.CONSTANTS.e0 * uc.UNITS.voltage)
+        case 'kT':
+            return 1 / (params.temperature * uc.CONSTANTS.Boltzmann)
+        case 'J':
+            return uc.UNITS.energy
 
-    return Expansion(common_l_array, new_coeffs1), Expansion(common_l_array, new_coeffs2)
 
+def charged_shell_energy(ex1: Expansion, ex2: Expansion, dist: float, params: ModelParams, units: EnergyUnit = 'kT'):
 
-def charged_shell_energy(ex1: Expansion, ex2: Expansion, dist: float, params: ModelParams):
-
-    ex1, ex2 = expansions_to_common_l(ex1, ex2)
+    ex1, ex2 = expansion.expansions_to_common_l(ex1, ex2)
+    dist = dist * params.R
 
     full_l_array, full_m_array = ex1.lm_arrays
 
-    coefficient_C = fn.interaction_coeff_C(ex1.l_array[:, None], ex2.l_array[None, :], params.kappaR)
-    full_coefficient_C = ex1.repeat_over_m(ex2.repeat_over_m(coefficient_C, axis=1), axis=0)
-
-    charge_factor = np.real(ex1.coeffs[:, None] * np.conj(ex2.coeffs[None, :]) +
-                            (-1) ** (full_l_array[:, None] + full_l_array[None, :])
-                            * ex1.coeffs[None, :] * ex1.coeffs[:, None])
-
+    # determine indices of relevant elements in the sum
     indices_l, indices_p = np.nonzero(full_m_array[:, None] == full_m_array[None, :])
     flat_l = full_l_array[indices_l]
     flat_p = full_l_array[indices_p]
     flat_m = full_m_array[indices_l]  # the same as full_m_array[indices_p]
 
+    charge_factor = np.real(ex1.coefs[..., indices_l] * np.conj(ex2.coefs[..., indices_p]) +
+                            (-1) ** (flat_l + flat_p) * ex1.coefs[..., indices_p] * np.conj(ex2.coefs[..., indices_l]))
+
     all_s_array = np.arange(2 * ex1.max_l + 1)
     bessels = fn.sph_bessel_k(all_s_array, params.kappa * dist)
 
+    # additional selection that excludes terms where Wigner 3j symbols are 0 by definition
     s_bool1 = np.abs(flat_l[:, None] - all_s_array[None, :]) <= flat_p[:, None]
     s_bool2 = flat_p[:, None] <= (flat_l[:, None] + all_s_array[None, :])
     indices_lpm, indices_s = np.nonzero(s_bool1 * s_bool2)
@@ -73,27 +51,39 @@ def charged_shell_energy(ex1: Expansion, ex2: Expansion, dist: float, params: Mo
     l_vals = flat_l[indices_lpm]
     p_vals = flat_p[indices_lpm]
     m_vals = flat_m[indices_lpm]
-    C_vals = full_coefficient_C[indices_l, indices_p][indices_lpm]
-    charge_vals = charge_factor[indices_l, indices_p][indices_lpm]
     s_vals = all_s_array[indices_s]
     bessel_vals = bessels[indices_s]
 
+    # While all other arrays are 1D, this one can have extra leading axes corresponding to leading dimensions
+    # of expansion coefficients. The last dimension is the same as other arrays.
+    charge_vals = charge_factor[..., indices_lpm]
+
     lps_terms = (2 * s_vals + 1) * np.sqrt((2 * l_vals + 1) * (2 * p_vals + 1))
 
-    wigner1 = wigner3j(l_vals, s_vals, p_vals,
-                       0, 0, 0)
-    wigner2 = wigner3j(l_vals, s_vals, p_vals,
-                       -m_vals, 0, m_vals)
+    # the same combination of l, s, and p is repeated many times
+    _, unique_indices1, inverse1 = np.unique(np.stack((l_vals, s_vals, p_vals)),
+                                             axis=1, return_inverse=True, return_index=True)
+    wigner1 = wigner3j(2 * l_vals[unique_indices1], 2 * s_vals[unique_indices1], 2 * p_vals[unique_indices1],
+                       0, 0, 0)[inverse1]
+
+    # all the combinations (l, s, p, m) are unique
+    wigner2 = wigner3j(2 * l_vals, 2 * s_vals, 2 * p_vals,
+                       -2 * m_vals, 0, 2 * m_vals)
+
+    constants = params.R ** 2 / (params.kappa * params.epsilon * uc.CONSTANTS.epsilon0) * energy_units(units, params)
+
+    C_vals = fn.interaction_coeff_C(l_vals, p_vals, params.kappaR)
+    lspm_vals = C_vals * (-1) ** (l_vals + m_vals) * lps_terms * bessel_vals * wigner1 * wigner2
+    broadcasted_lspm_vals = np.broadcast_to(lspm_vals, charge_vals.shape)
 
-    return (params.R ** 2 / (params.kappa * params.epsilon * params.epsilon0)
-            * np.sum(C_vals * (-1) ** (l_vals + m_vals) * charge_vals * lps_terms * bessel_vals * wigner1 * wigner2))
+    return 0.5 * constants * np.sum(broadcasted_lspm_vals * charge_vals, axis=-1)
 
 
 if __name__ == '__main__':
 
-    params = ModelParams(1, 3, 1, 1)
-    ex1 = expansion.MappedExpansion(1, params.kappaR, 0.001, max_l=5)
-    ex2 = expansion.MappedExpansion(1, params.kappaR, 0.001, max_l=5)
+    params = ModelParams(R=150, kappaR=3)
+    ex1 = expansion.MappedExpansion(1, params.kappaR, 0.001, max_l=20)
+    ex2 = ex1.clone()
 
     dist = 2.
 

parameters.py

@@ -1,13 +1,41 @@
 from dataclasses import dataclass
+import numpy as np
+import units_and_constants as uc
 
 
 @dataclass
 class ModelParams:
     R: float
-    kappa: float
-    epsilon: float
-    epsilon0: float
+    kappa: float = None
+    kappaR: float = None
+    c0: float = None
+    epsilon: float = 80
+    temperature: float = 293
 
-    @property
-    def kappaR(self):
-        return self.kappa * self.R
+    def __post_init__(self):
+        self.kappa, self.kappaR, self.c0 = screening_calculator(self.R, self.temperature, self.epsilon,
+                                                                self.kappa, self.kappaR, self.c0)
+
+
+def bjerrum(temp: float, epsilon: float) -> float:
+    return uc.CONSTANTS.e0 ** 2 / (4 * np.pi * epsilon * uc.CONSTANTS.epsilon0 * uc.CONSTANTS.Boltzmann * temp)
+
+
+def kappa_from_concentration(c0: float, temp: float, epsilon: float) -> float:
+    return np.sqrt(8 * np.pi * bjerrum(temp, epsilon) * c0)
+
+
+def concentration_from_kappa(kappa: float, temp: float, epsilon: float):
+    return kappa ** 2 / (8 * np.pi * bjerrum(temp, epsilon))
+
+
+def screening_calculator(radius: float, temp: float, epsilon: float,
+                         kappa: float = None, kappaR: float = None, c0: float = None) -> (float, float, float):
+    if kappa is not None:
+        return kappa, kappa * radius, concentration_from_kappa(kappa, temp, epsilon)
+    elif kappaR is not None:
+        return kappaR / radius, kappaR, concentration_from_kappa(kappaR / radius, temp, epsilon)
+    elif c0 is not None:
+        kappa = kappa_from_concentration(c0, temp, epsilon)
+        return kappa, kappa * radius, c0
+    raise ValueError('One of the arguments kappa, kappaR or c0 should be different from None.')

potential_tests.py

@@ -0,0 +1,38 @@
+from interactions import charged_shell_energy
+import expansion
+from parameters import ModelParams
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+def v22_distance_test():
+
+    params = ModelParams(R=10, kappaR=3.29)
+    ex0 = expansion.Expansion24(1, 0, 0)
+    ex1 = ex0.clone()
+
+    ex0.rotate_euler(0, np.array([0, 0, np.pi / 2]), 0)
+    ex1.rotate_euler(0, np.array([0, np.pi / 2, np.pi / 2]), 0)
+
+    dist = np.linspace(2, 3.2, 100)
+    energy_array = np.zeros((dist.shape[0], 3))
+    for i, d in enumerate(dist):
+        energy_array[i, ...] = charged_shell_energy(ex0, ex1, d, params)
+
+    print(charged_shell_energy(ex0, ex1, 2., params))
+
+    plt.plot(dist, energy_array)
+    plt.show()
+
+
+if __name__ == '__main__':
+
+    v22_distance_test()
+
+
+
+
+
+
+
+

potentials.py

@@ -2,6 +2,7 @@ import expansion
 import functions as fn
 import numpy as np
 import parameters
+import units_and_constants as uc
 import matplotlib.pyplot as plt
 
 
@@ -34,22 +35,23 @@ def charged_shell_potential(theta: Array | float,
         raise ValueError('theta and phi arrays should have the same shape.')
     l_array, m_array = ex.lm_arrays
 
+    dist = dist * params.R
     l_factors = (fn.coefficient_Cpm(ex.l_array, params.kappaR) * fn.sph_bessel_k(ex.l_array, params.kappa * dist)
                  / fn.sph_bessel_k(ex.l_array, params.kappaR))
 
-    return (1 / (params.kappa * params.epsilon * params.epsilon0)
-            * np.real(np.sum(ex.repeat_over_m(l_factors)[None, :] * ex.coeffs
+    return (1 / (params.kappa * params.epsilon * uc.CONSTANTS.epsilon0)
+            * np.real(np.sum(ex.repeat_over_m(l_factors)[None, :] * ex.coefs
                              * fn.sph_harm(l_array[None, :], m_array[None, :], theta[:, None], phi[:, None]), axis=1)))
 
 
 if __name__ == '__main__':
 
-    params = ModelParams(1, 3, 1, 1)
-    ex = expansion.MappedExpansion(1, params.kappaR, 0.001, max_l=20)
+    params = ModelParams(R=150, kappaR=3)
+    ex = expansion.MappedExpansion(1, params.kappaR, 0.001, max_l=10)
 
     theta = np.linspace(0, np.pi, 1000)
     phi = 0.
-    dist = 1.
+    dist = 2
 
     potential = charged_shell_potential(theta, phi, dist, ex, params)
     # print(potential)

rotational_path.py

@@ -0,0 +1,66 @@
+import numpy as np
+from dataclasses import dataclass, field
+import quaternionic
+import expansion
+from parameters import ModelParams
+import interactions
+import matplotlib.pyplot as plt
+
+
+Quaternion = quaternionic.array
+Array = np.ndarray
+
+
+@dataclass
+class PairRotationalPath:
+
+    rotations1: list[Quaternion] = field(default_factory=list)
+    rotations2: list[Quaternion] = field(default_factory=list)
+
+    def add(self, rotation1: Quaternion, rotation2: Quaternion):
+        self.rotations1.append(rotation1)
+        self.rotations2.append(rotation2)
+
+    def add_euler(self, alpha1: Array = 0, beta1: Array = 0, gamma1: Array = 0,
+                  alpha2: Array = 0, beta2: Array = 0, gamma2: Array = 0):
+        R1_euler = quaternionic.array.from_euler_angles(alpha1, beta1, gamma1)
+        R2_euler = quaternionic.array.from_euler_angles(alpha2, beta2, gamma2)
+        R1_euler, R2_euler = np.broadcast_arrays(R1_euler, R2_euler)
+        self.add(Quaternion(R1_euler), Quaternion(R2_euler))
+
+    def get_rotations(self) -> (Quaternion, Quaternion):
+        return Quaternion(np.vstack(self.rotations1)), Quaternion(np.vstack(self.rotations2))
+
+
+zero_to_pi_half = np.linspace(0, np.pi/2, 1000)
+
+QuadPath = PairRotationalPath()
+QuadPath.add_euler(beta1=np.pi/2, beta2=zero_to_pi_half)
+QuadPath.add_euler(beta1=zero_to_pi_half[::-1], beta2=zero_to_pi_half[::-1])
+QuadPath.add_euler(beta1=zero_to_pi_half)
+QuadPath.add_euler(beta1=zero_to_pi_half[::-1], beta2=zero_to_pi_half)
+QuadPath.add_euler(beta1=np.pi/2, beta2=zero_to_pi_half, alpha2=np.pi/2)
+QuadPath.add_euler(beta1=np.pi/2, beta2=np.pi/2, alpha2=zero_to_pi_half[::-1])
+
+
+if __name__ == '__main__':
+
+    params = ModelParams(R=150, kappaR=3)
+    ex1 = expansion.MappedExpansion(1, params.kappaR, 0.001, max_l=15)
+    # ex1 = expansion.Expansion24(sigma2=0.001, sigma4=0)
+    ex2 = ex1.clone()
+
+    rotations1, rotations2 = QuadPath.get_rotations()
+    ex1.rotate(rotations1)
+    ex2.rotate(rotations2)
+
+    dist = 3
+
+    energy = interactions.charged_shell_energy(ex1, ex2, dist, params)
+
+    plt.plot(energy)
+    plt.show()
+
+
+
+

units_and_constants.py

@@ -0,0 +1,42 @@
+from dataclasses import dataclass
+
+
+@dataclass
+class BaseUnits:
+    distance: float
+    charge: float
+    voltage: float
+    concentrationM: float
+
+    @property
+    def charge_density(self):
+        return self.charge / (self.distance ** 2)
+
+    @property
+    def energy(self):
+        return self.charge * self.voltage
+
+    @property
+    def concentration(self):
+        return self.concentrationM * 6.02214076 * 1e23 / (10 * self.distance) ** 3
+
+
+@dataclass
+class Constants:
+    Boltzmann: float
+    epsilon0: float
+    e0: float
+
+
+base_units = {'distance': 1e-9,
+              'charge': 1.602176634 * 1e-19,
+              'voltage': 1,
+              'concentrationM': 1e-3}
+
+
+UNITS = BaseUnits(**base_units)
+
+# values of different constants in provided base units
+CONSTANTS = Constants(Boltzmann=1.380649 * 1e-23 / UNITS.energy,
+                      epsilon0=8.8541878128 * 1e-12 * UNITS.distance * UNITS.voltage / UNITS.charge,
+                      e0=1.602176634 * 1e-19 / UNITS.charge)

wigner_test.py

@@ -0,0 +1,11 @@
+import numpy as np
+from py3nj import wigner3j
+
+
+if __name__ == '__main__':
+
+    w3j = wigner3j(0, 38, 38,
+                   0, 38, -38)
+
+    print(w3j)
+