original_docstring: 
      (y, dy_dbeta, dy_dalpha) = spherical_harmonic_all(max_order, alpha, sinbeta, cosbeta)
      	
      Computes a Spherical Harmonic function and it's angular derivatives for
      all (m,n) up to the given maximum order. The algorithm is equivalent to that
      implemented in SphericalHarmonic, but this version avoids repeated calls
      to lpmv, since that is very time consuming.
      	
      Arguments - these should all be scalars:
      r is radius
      alpha is azimuth angle (angle in radians from the positive x axis, with
      rotation around the positive z axis according to the right-hand screw rule)
      beta is polar angle, but it is specified as two arrays of its cos and sin values. 
      max_order is maximum Spherical Harmonic order and should be a non-negative real integer scalar
      	
      Returned data will be vectors of length (max_order+1)^2.
      
      
docstring: (y, dy_dbeta, dy_dalpha) = spherical_harmonic_all(max_order, alpha, sinbeta, cosbeta)
  Computes a Spherical Harmonic function and it's angular derivatives for
  all (m,n) up to the given maximum order. The algorithm is equivalent to that
  implemented in SphericalHarmonic, but this version avoids repeated calls
  to lpmv, since that is very time consuming.
  
  these should all be scalars:
  r is radius
  alpha is azimuth angle (angle in radians from the positive x axis, with
  rotation around the positive z axis according to the right-hand screw rule)
  beta is polar angle, but it is specified as two arrays of its cos and sin values.
  max_order is maximum Spherical Harmonic order and should be a non-negative real integer scalar
  
  Returned data will be vectors of length (max_order+1)^2.
docstring_tokens: ( y dy_dbeta dy_dalpha ) = spherical_harmonic_all ( max_order alpha sinbeta cosbeta ) Computes a Spherical Harmonic function and it ' s angular derivatives for all ( m n ) up to the given maximum order . The algorithm is equivalent to that implemented in SphericalHarmonic but this version avoids repeated calls to lpmv since that is very time consuming . these should all be scalars : r is radius alpha is azimuth angle ( angle in radians from the positive x axis with rotation around the positive z axis according to the right - hand screw rule ) beta is polar angle but it is specified as two arrays of its cos and sin values . max_order is maximum Spherical Harmonic order and should be a non - negative real integer scalar Returned data will be vectors of length ( max_order + 1 ) ^2 .
code: def spherical_harmonic_all (max_order, alpha, sinbeta, cosbeta):
      
      """
      (y, dy_dbeta, dy_dalpha) = spherical_harmonic_all(max_order, alpha, sinbeta, cosbeta)
      	
      Computes a Spherical Harmonic function and it's angular derivatives for
      all (m,n) up to the given maximum order. The algorithm is equivalent to that
      implemented in SphericalHarmonic, but this version avoids repeated calls
      to lpmv, since that is very time consuming.
      	
      Arguments - these should all be scalars:
      r is radius
      alpha is azimuth angle (angle in radians from the positive x axis, with
      rotation around the positive z axis according to the right-hand screw rule)
      beta is polar angle, but it is specified as two arrays of its cos and sin values. 
      max_order is maximum Spherical Harmonic order and should be a non-negative real integer scalar
      	
      Returned data will be vectors of length (max_order+1)^2.
      
      """
      
      # alpha = np.array([[alpha]])
      # cosbeta = np.array([[cosbeta]])
      # sinbeta = np.array([[sinbeta]])
  
      # Preallocate output arrays:
      y =  np.zeros((np.size(alpha),(max_order+1)**2), np.complex128)
      dy_dbeta = np.zeros((np.size(alpha),(max_order+1)**2), np.complex128)
      dy_dalpha = np.zeros((np.size(alpha),(max_order+1)**2), np.complex128)
      
      #% Loop over n and calculate spherical harmonic functions y_nm
      for n in range(max_order+1):
      
          # Compute Legendre function and its derivatives for all m:
          p_n = lpmv(range(0,n+1), n, cosbeta)
          #shape_p_n = np.shape(p_n)
          #p_n = p_n.reshape((shape_p_n[1], shape_p_n[0]))
          
          for m in range(-n, n+1):
      
              # Legendre function its derivatives for |m|:
              p_nm = p_n[:, np.absolute(m)]
              p_nm = p_nm.reshape((np.size(p_nm), ))
  
              if n==0:
                  dPmn_dbeta = 0
              elif m==0:
                  dPmn_dbeta = p_n[:,1]
              elif abs(m)<n:
                  dPmn_dbeta = 0.5*p_n[:,abs(m)+1] - 0.5*(n+abs(m))*(n-abs(m)+1)*p_n[:,abs(m)-1];
                  dPmn_dbeta = dPmn_dbeta.reshape((np.size(dPmn_dbeta), ))
              elif (abs(m)==1) and (n==1):
                  dPmn_dbeta = -cosbeta
                  dPmn_dbeta = dPmn_dbeta.reshape((np.size(dPmn_dbeta), ))
              #elif sinbeta<=np.finfo(float).eps:
                  #dPmn_dbeta = 0
              else:
                  dPmn_dbeta = -abs(m)*cosbeta.reshape((np.size(cosbeta), ))*p_nm/sinbeta.reshape((np.size(sinbeta), )) - (n+abs(m))*(n-abs(m)+1)*p_n[:,abs(m)-1]
                  dPmn_dbeta = dPmn_dbeta.reshape((np.size(dPmn_dbeta), ))
              
              #print (dPmn_dbeta)
              #print (np.shape(dPmn_dbeta))
              # Compute scaling term, including sign factor:
              scaling_term = ((-1)**m) * np.sqrt((2 * n + 1) / (4 * np.pi * np.prod(np.float64(range(n-abs(m)+1, n+abs(m)+1)))))
      
              # Compute exponential term:
              exp_term = np.exp(1j*m*alpha)
              exp_term = exp_term.reshape((np.size(exp_term), ))
              
              # Put it all together:
              i = sub2indSH(m,n)
              #y[:,i] = np.multiply (exp_term, p_nm)
              y[:,i] = scaling_term * exp_term * p_nm
              dy_dbeta[:,i] = scaling_term * exp_term * dPmn_dbeta
              dy_dalpha[:,i] = y[:,i] * 1j * m
              
      return y, dy_dbeta, dy_dalpha
code_tokens: def spherical_harmonic_all ( max_order , alpha , sinbeta , cosbeta ) : y = np . zeros ( ( np . size ( alpha ) , ( max_order + 1 ) ** 2 ) , np . complex128 ) dy_dbeta = np . zeros ( ( np . size ( alpha ) , ( max_order + 1 ) ** 2 ) , np . complex128 ) dy_dalpha = np . zeros ( ( np . size ( alpha ) , ( max_order + 1 ) ** 2 ) , np . complex128 ) for n in range ( max_order + 1 ) : p_n = lpmv ( range ( 0 , n + 1 ) , n , cosbeta ) for m in range ( - n , n + 1 ) : p_nm = p_n [ : , np . absolute ( m ) ] p_nm = p_nm . reshape ( ( np . size ( p_nm ) , ) ) if n == 0 : dPmn_dbeta = 0 elif m == 0 : dPmn_dbeta = p_n [ : , 1 ] elif abs ( m ) < n : dPmn_dbeta = 0.5 * p_n [ : , abs ( m ) + 1 ] - 0.5 * ( n + abs ( m ) ) * ( n - abs ( m ) + 1 ) * p_n [ : , abs ( m ) - 1 ] ; dPmn_dbeta = dPmn_dbeta . reshape ( ( np . size ( dPmn_dbeta ) , ) ) elif ( abs ( m ) == 1 ) and ( n == 1 ) : dPmn_dbeta = - cosbeta dPmn_dbeta = dPmn_dbeta . reshape ( ( np . size ( dPmn_dbeta ) , ) ) else : dPmn_dbeta = - abs ( m ) * cosbeta . reshape ( ( np . size ( cosbeta ) , ) ) * p_nm / sinbeta . reshape ( ( np . size ( sinbeta ) , ) ) - ( n + abs ( m ) ) * ( n - abs ( m ) + 1 ) * p_n [ : , abs ( m ) - 1 ] dPmn_dbeta = dPmn_dbeta . reshape ( ( np . size ( dPmn_dbeta ) , ) ) scaling_term = ( ( - 1 ) ** m ) * np . sqrt ( ( 2 * n + 1 ) / ( 4 * np . pi * np . prod ( np . float64 ( range ( n - abs ( m ) + 1 , n + abs ( m ) + 1 ) ) ) ) ) exp_term = np . exp ( 1j * m * alpha ) exp_term = exp_term . reshape ( ( np . size ( exp_term ) , ) ) i = sub2indSH ( m , n ) y [ : , i ] = scaling_term * exp_term * p_nm dy_dbeta [ : , i ] = scaling_term * exp_term * dPmn_dbeta dy_dalpha [ : , i ] = y [ : , i ] * 1j * m return y , dy_dbeta , dy_dalpha
short_docstring: (y, dy_dbeta, dy_dalpha) = spherical_harmonic_all(max_order, alpha, sinbeta, cosbeta)
  Computes a Spherical Harmonic function and it's angular derivatives for
  all (m,n) up to the given maximum order.
short_docstring_tokens: ( y dy_dbeta dy_dalpha ) = spherical_harmonic_all ( max_order alpha sinbeta cosbeta ) Computes a Spherical Harmonic function and it ' s angular derivatives for all ( m n ) up to the given maximum order .
comment: """
      (y, dy_dbeta, dy_dalpha) = spherical_harmonic_all(max_order, alpha, sinbeta, cosbeta)
      	
      Computes a Spherical Harmonic function and it's angular derivatives for
      all (m,n) up to the given maximum order. The algorithm is equivalent to that
      implemented in SphericalHarmonic, but this version avoids repeated calls
      to lpmv, since that is very time consuming.
      	
      Arguments - these should all be scalars:
      r is radius
      alpha is azimuth angle (angle in radians from the positive x axis, with
      rotation around the positive z axis according to the right-hand screw rule)
      beta is polar angle, but it is specified as two arrays of its cos and sin values. 
      max_order is maximum Spherical Harmonic order and should be a non-negative real integer scalar
      	
      Returned data will be vectors of length (max_order+1)^2.
      
      """$SEP$# alpha = np.array([[alpha]])$SEP$# cosbeta = np.array([[cosbeta]])$SEP$# sinbeta = np.array([[sinbeta]])$SEP$# Preallocate output arrays:$SEP$#% Loop over n and calculate spherical harmonic functions y_nm$SEP$# Compute Legendre function and its derivatives for all m:$SEP$#shape_p_n = np.shape(p_n)$SEP$#p_n = p_n.reshape((shape_p_n[1], shape_p_n[0]))$SEP$# Legendre function its derivatives for |m|:$SEP$#elif sinbeta<=np.finfo(float).eps:$SEP$#dPmn_dbeta = 0$SEP$#print (dPmn_dbeta)$SEP$#print (np.shape(dPmn_dbeta))$SEP$# Compute scaling term, including sign factor:$SEP$# Compute exponential term:$SEP$# Put it all together:$SEP$#y[:,i] = np.multiply (exp_term, p_nm)