Coverage for yasfpy/functions/legendre_normalized_trigon.py: 7%

1import numpy as np 

2 

3# https://github.com/disordered-photonics/celes/blob/master/src/mathematics/legendre_normalized_trigon.m 

4 

5 

6def legendre_normalized_trigon(lmax, x: np.ndarray, y: np.ndarray = None): 

7 # if not(isinstance(x, np.ndarray) or np.isscalar(x) or isinstance(x, list)): 

8 # return legendre_normalized_trigon_legacy(x, y, lmax) 

9 

10 # if np.isscalar(x): 

11 # x = np.array([x]) 

12 # elif isinstance(x, list): 

13 # x = np.array(x) 

14 size = x.shape 

15 if y is None: 

16 x = np.ravel(x) 

17 ct = np.cos(x) 

18 st = np.sin(x) 

19 else: 

20 ct = x.ravel() 

21 st = y.ravel() 

22 

23 plm = np.zeros((lmax + 1, lmax + 1, ct.size)) * np.nan 

24 

25 plm[0, 0, :] = np.sqrt(1 / 2) * np.ones_like(ct) 

26 plm[1, 0, :] = np.sqrt(3 / 2) * ct 

27 

28 for l in range(1, lmax): 

29 plm[l + 1, 0, :] = ( 

30 1 / (l + 1) * np.sqrt((2 * l + 1) * (2 * l + 3)) * plm[l, 0, :] * ct 

31 - l / (l + 1) * np.sqrt((2 * l + 3) / (2 * l - 1)) * plm[l - 1, 0, :] 

32 ) 

33 

34 for m in range(1, lmax + 1): 

35 plm[m - 1, m, :] = np.zeros_like(ct) 

36 plm[m, m, :] = ( 

37 np.sqrt((2 * m + 1) / 2 / np.math.factorial(2 * m)) 

38 * np.prod(np.arange(1, 2 * m, 2)) 

39 * np.power(st, m) 

40 ) 

41 for l in range(m, lmax): 

42 plm[l + 1, m, :] = ( 

43 np.sqrt((2 * l + 1) * (2 * l + 3) / (l + 1 - m) / (l + 1 + m)) 

44 * ct 

45 * plm[l, m, :] 

46 - np.sqrt( 

47 (2 * l + 3) 

48 * (l - m) 

49 * (l + m) 

50 / (2 * l - 1) 

51 / (l + 1 - m) 

52 / (l + 1 + m) 

53 ) 

54 * plm[l - 1, m, :] 

55 ) 

56 

57 plm = np.reshape(plm, np.concatenate(([lmax + 1, lmax + 1], size))) 

58 

59 return plm 

60 

61 

62def legendre_normalized_trigon_legacy(x, y=None, lmax=4): 

63 import sympy as sym 

64 

65 plm = sym.zeros(lmax + 1, lmax + 1) 

66 if y == None: 

67 ct = sym.cos(x) 

68 st = sym.sin(x) 

69 elif isinstance(x, sym.core.symbol.Symbol) and isinstance( 

70 y, sym.core.symbol.Symbol 

71 ): 

72 ct = x 

73 st = y 

74 else: 

75 ct = sym.Symbol("ct") 

76 st = sym.Symbol("st") 

77 

78 plm[0, 0] = np.sqrt(1 / 2) 

79 plm[1, 0] = np.sqrt(3 / 2) * ct 

80 

81 for l in range(1, lmax): 

82 plm[l + 1, 0] = ( 

83 1 / (l + 1) * np.sqrt((2 * l + 1) * (2 * l + 3)) * plm[l, 0] * ct 

84 - l / (l + 1) * np.sqrt((2 * l + 3) / (2 * l - 1)) * plm[l - 1, 0] 

85 ) 

86 

87 for m in range(1, lmax + 1): 

88 plm[m - 1, m] = np.zeros_like(ct) 

89 plm[m, m] = ( 

90 np.sqrt((2 * m + 1) / 2 / np.math.factorial(2 * m)) 

91 * np.prod(np.arange(1, 2 * m, 2)) 

92 * st**m 

93 ) 

94 for l in range(m, lmax): 

95 plm[l + 1, m] = ( 

96 np.sqrt((2 * l + 1) * (2 * l + 3) / (l + 1 - m) / (l + 1 + m)) 

97 * ct 

98 * plm[l, m] 

99 - np.sqrt( 

100 (2 * l + 3) 

101 * (l - m) 

102 * (l + m) 

103 / (2 * l - 1) 

104 / (l + 1 - m) 

105 / (l + 1 + m) 

106 ) 

107 * plm[l - 1, m] 

108 ) 

109 

110 return plm