code: def real_sph_harm(l, zero_m_only=False, spherical_coordinates=True):
      """
      Computes formula strings of the the real part of the spherical harmonics up to order l (excluded).
      Variables are either cartesian coordinates x,y,z on the unit sphere or spherical coordinates phi and theta.
      """
      if not zero_m_only:
          x = sym.symbols("x")
          y = sym.symbols("y")
          S_m = [x * 0]
          C_m = [1 + 0 * x]
          # S_m = [0]
          # C_m = [1]
          for i in range(1, l):
              x = sym.symbols("x")
              y = sym.symbols("y")
              S_m += [x * S_m[i - 1] + y * C_m[i - 1]]
              C_m += [x * C_m[i - 1] - y * S_m[i - 1]]
  
      P_l_m = associated_legendre_polynomials(l, zero_m_only)
      if spherical_coordinates:
          theta = sym.symbols("theta")
          z = sym.symbols("z")
          for i in range(len(P_l_m)):
              for j in range(len(P_l_m[i])):
                  if type(P_l_m[i][j]) != int:
                      P_l_m[i][j] = P_l_m[i][j].subs(z, sym.cos(theta))
          if not zero_m_only:
              phi = sym.symbols("phi")
              for i in range(len(S_m)):
                  S_m[i] = (
                      S_m[i]
                      .subs(x, sym.sin(theta) * sym.cos(phi))
                      .subs(y, sym.sin(theta) * sym.sin(phi))
                  )
              for i in range(len(C_m)):
                  C_m[i] = (
                      C_m[i]
                      .subs(x, sym.sin(theta) * sym.cos(phi))
                      .subs(y, sym.sin(theta) * sym.sin(phi))
                  )
  
      Y_func_l_m = [["0"] * (2 * j + 1) for j in range(l)]
      for i in range(l):
          Y_func_l_m[i][0] = sym.simplify(
              sph_harm_prefactor(i, 0) * P_l_m[i][0]
          )
  
      if not zero_m_only:
          for i in range(1, l):
              for j in range(1, i + 1):
                  Y_func_l_m[i][j] = sym.simplify(
                      2**0.5
                      * sph_harm_prefactor(i, j)
                      * C_m[j]
                      * P_l_m[i][j]
                  )
          for i in range(1, l):
              for j in range(1, i + 1):
                  Y_func_l_m[i][-j] = sym.simplify(
                      2**0.5
                      * sph_harm_prefactor(i, -j)
                      * S_m[j]
                      * P_l_m[i][j]
                  )
  
      return Y_func_l_m
  
imports: import numpy
  from scipy.optimize import brentq
  from scipy import special
  import torch
  from math import pi
  import sympy
label_desc: Use image node_id 7 for calling a global function with example usage: real_sph_harm(l, zero_m_only, spherical_coordinates) and returns: Y_func_l_m