|
| 1 | +@testmodule LAL begin |
| 2 | + |
| 3 | + """ |
| 4 | + Reproduces the XLALSpinWeightedSphericalHarmonic function from the LALSuite C library: |
| 5 | + https://lscsoft.docs.ligo.org/lalsuite/lal/_spherical_harmonics_8c_source.html#l00042 |
| 6 | + """ |
| 7 | + function LALSpinWeightedSphericalHarmonic( |
| 8 | + theta::Float64, # polar angle (rad) |
| 9 | + phi::Float64, # azimuthal angle (rad) |
| 10 | + s::Int, # spin weight |
| 11 | + l::Int, # mode number l |
| 12 | + m::Int # mode number m |
| 13 | + ) |
| 14 | + # Sanity checks |
| 15 | + if l < abs(s) |
| 16 | + error("Invalid mode s=$s, l=$l, m=$m - require |s| <= l") |
| 17 | + end |
| 18 | + if l < abs(m) |
| 19 | + error("Invalid mode s=$s, l=$l, m=$m - require |m| <= l") |
| 20 | + end |
| 21 | + if s != -2 |
| 22 | + error("Unsupported mode s=$s (only s=-2 implemented)") |
| 23 | + end |
| 24 | + if l < 2 || l > 8 |
| 25 | + error("Unsupported mode l=$l (only l in [2,8] implemented)") |
| 26 | + end |
| 27 | + # Compute real factor |
| 28 | + fac = if l == 2 |
| 29 | + if m == -2 |
| 30 | + sqrt(5.0 / (64.0 * π)) * (1.0 - cos(theta)) * (1.0 - cos(theta)) |
| 31 | + elseif m == -1 |
| 32 | + sqrt(5.0 / (16.0 * π)) * sin(theta) * (1.0 - cos(theta)) |
| 33 | + elseif m == 0 |
| 34 | + sqrt(15.0 / (32.0 * π)) * sin(theta) * sin(theta) |
| 35 | + elseif m == 1 |
| 36 | + sqrt(5.0 / (16.0 * π)) * sin(theta) * (1.0 + cos(theta)) |
| 37 | + elseif m == 2 |
| 38 | + sqrt(5.0 / (64.0 * π)) * (1.0 + cos(theta)) * (1.0 + cos(theta)) |
| 39 | + end |
| 40 | + elseif l == 3 |
| 41 | + if m == -3 |
| 42 | + sqrt(21.0 / (2.0 * π)) * cos(theta/2.0) * sin(theta/2.0)^5 |
| 43 | + elseif m == -2 |
| 44 | + sqrt(7.0 / (4.0 * π)) * (2.0 + 3.0 * cos(theta)) * sin(theta/2.0)^4 |
| 45 | + elseif m == -1 |
| 46 | + sqrt(35.0 / (2.0 * π)) * (sin(theta) + 4.0 * sin(2.0 * theta) - 3.0 * sin(3.0 * theta)) / 32.0 |
| 47 | + elseif m == 0 |
| 48 | + sqrt(105.0 / (2.0 * π)) * cos(theta) * sin(theta)^2 / 4.0 |
| 49 | + elseif m == 1 |
| 50 | + -sqrt(35.0 / (2.0 * π)) * (sin(theta) - 4.0 * sin(2.0 * theta) - 3.0 * sin(3.0 * theta)) / 32.0 |
| 51 | + elseif m == 2 |
| 52 | + sqrt(7.0 / π) * cos(theta/2.0)^4 * (-2.0 + 3.0 * cos(theta)) / 2.0 |
| 53 | + elseif m == 3 |
| 54 | + -sqrt(21.0 / (2.0 * π)) * cos(theta/2.0)^5 * sin(theta/2.0) |
| 55 | + end |
| 56 | + elseif l == 4 |
| 57 | + if m == -4 |
| 58 | + 3.0 * sqrt(7.0 / π) * cos(theta/2.0)^2 * sin(theta/2.0)^6 |
| 59 | + elseif m == -3 |
| 60 | + 3.0 * sqrt(7.0 / (2.0 * π)) * cos(theta/2.0) * (1.0 + 2.0 * cos(theta)) * sin(theta/2.0)^5 |
| 61 | + elseif m == -2 |
| 62 | + 3.0 * (9.0 + 14.0 * cos(theta) + 7.0 * cos(2.0 * theta)) * sin(theta/2.0)^4 / (4.0 * sqrt(π)) |
| 63 | + elseif m == -1 |
| 64 | + 3.0 * (3.0 * sin(theta) + 2.0 * sin(2.0 * theta) + 7.0 * sin(3.0 * theta) - 7.0 * sin(4.0 * theta)) / (32.0 * sqrt(2.0 * π)) |
| 65 | + elseif m == 0 |
| 66 | + 3.0 * sqrt(5.0 / (2.0 * π)) * (5.0 + 7.0 * cos(2.0 * theta)) * sin(theta)^2 / 16.0 |
| 67 | + elseif m == 1 |
| 68 | + 3.0 * (3.0 * sin(theta) - 2.0 * sin(2.0 * theta) + 7.0 * sin(3.0 * theta) + 7.0 * sin(4.0 * theta)) / (32.0 * sqrt(2.0 * π)) |
| 69 | + elseif m == 2 |
| 70 | + 3.0 * cos(theta/2.0)^4 * (9.0 - 14.0 * cos(theta) + 7.0 * cos(2.0 * theta)) / (4.0 * sqrt(π)) |
| 71 | + elseif m == 3 |
| 72 | + -3.0 * sqrt(7.0 / (2.0 * π)) * cos(theta/2.0)^5 * (-1.0 + 2.0 * cos(theta)) * sin(theta/2.0) |
| 73 | + elseif m == 4 |
| 74 | + 3.0 * sqrt(7.0 / π) * cos(theta/2.0)^6 * sin(theta/2.0)^2 |
| 75 | + end |
| 76 | + elseif l == 5 |
| 77 | + if m == -5 |
| 78 | + sqrt(330.0 / π) * cos(theta/2.0)^3 * sin(theta/2.0)^7 |
| 79 | + elseif m == -4 |
| 80 | + sqrt(33.0 / π) * cos(theta/2.0)^2 * (2.0 + 5.0 * cos(theta)) * sin(theta/2.0)^6 |
| 81 | + elseif m == -3 |
| 82 | + sqrt(33.0 / (2.0 * π)) * cos(theta/2.0) * (17.0 + 24.0 * cos(theta) + 15.0 * cos(2.0 * theta)) * sin(theta/2.0)^5 / 4.0 |
| 83 | + elseif m == -2 |
| 84 | + sqrt(11.0 / π) * (32.0 + 57.0 * cos(theta) + 36.0 * cos(2.0 * theta) + 15.0 * cos(3.0 * theta)) * sin(theta/2.0)^4 / 8.0 |
| 85 | + elseif m == -1 |
| 86 | + sqrt(77.0 / π) * (2.0 * sin(theta) + 8.0 * sin(2.0 * theta) + 3.0 * sin(3.0 * theta) + 12.0 * sin(4.0 * theta) - 15.0 * sin(5.0 * theta)) / 256.0 |
| 87 | + elseif m == 0 |
| 88 | + sqrt(1155.0 / (2.0 * π)) * (5.0 * cos(theta) + 3.0 * cos(3.0 * theta)) * sin(theta)^2 / 32.0 |
| 89 | + elseif m == 1 |
| 90 | + sqrt(77.0 / π) * (-2.0 * sin(theta) + 8.0 * sin(2.0 * theta) - 3.0 * sin(3.0 * theta) + 12.0 * sin(4.0 * theta) + 15.0 * sin(5.0 * theta)) / 256.0 |
| 91 | + elseif m == 2 |
| 92 | + sqrt(11.0 / π) * cos(theta/2.0)^4 * (-32.0 + 57.0 * cos(theta) - 36.0 * cos(2.0 * theta) + 15.0 * cos(3.0 * theta)) / 8.0 |
| 93 | + elseif m == 3 |
| 94 | + -sqrt(33.0 / (2.0 * π)) * cos(theta/2.0)^5 * (17.0 - 24.0 * cos(theta) + 15.0 * cos(2.0 * theta)) * sin(theta/2.0) / 4.0 |
| 95 | + elseif m == 4 |
| 96 | + sqrt(33.0 / π) * cos(theta/2.0)^6 * (-2.0 + 5.0 * cos(theta)) * sin(theta/2.0)^2 |
| 97 | + elseif m == 5 |
| 98 | + -sqrt(330.0 / π) * cos(theta/2.0)^7 * sin(theta/2.0)^3 |
| 99 | + end |
| 100 | + elseif l == 6 |
| 101 | + if m == -6 |
| 102 | + (3.0 * sqrt(715.0 / π) * cos(theta/2.0)^4 * sin(theta/2.0)^8) / 2.0 |
| 103 | + elseif m == -5 |
| 104 | + (sqrt(2145.0 / π) * cos(theta/2.0)^3 * (1.0 + 3.0 * cos(theta)) * sin(theta/2.0)^7) / 2.0 |
| 105 | + elseif m == -4 |
| 106 | + (sqrt(195.0 / (2.0 * π)) * cos(theta/2.0)^2 * (35.0 + 44.0 * cos(theta) + 33.0 * cos(2.0 * theta)) * sin(theta/2.0)^6) / 8.0 |
| 107 | + elseif m == -3 |
| 108 | + (3.0 * sqrt(13.0 / π) * cos(theta/2.0) * (98.0 + 185.0 * cos(theta) + 110.0 * cos(2.0 * theta) + 55.0 * cos(3.0 * theta)) * sin(theta/2.0)^5) / 32.0 |
| 109 | + elseif m == -2 |
| 110 | + (sqrt(13.0 / π) * (1709.0 + 3096.0 * cos(theta) + 2340.0 * cos(2.0 * theta) + 1320.0 * cos(3.0 * theta) + 495.0 * cos(4.0 * theta)) * sin(theta/2.0)^4) / 256.0 |
| 111 | + elseif m == -1 |
| 112 | + (sqrt(65.0 / (2.0 * π)) * cos(theta/2.0) * (161.0 + 252.0 * cos(theta) + 252.0 * cos(2.0 * theta) + 132.0 * cos(3.0 * theta) + 99.0 * cos(4.0 * theta)) * sin(theta/2.0)^3) / 64.0 |
| 113 | + elseif m == 0 |
| 114 | + (sqrt(1365.0 / π) * (35.0 + 60.0 * cos(2.0 * theta) + 33.0 * cos(4.0 * theta)) * sin(theta)^2) / 512.0 |
| 115 | + elseif m == 1 |
| 116 | + (sqrt(65.0 / (2.0 * π)) * cos(theta/2.0)^3 * (161.0 - 252.0 * cos(theta) + 252.0 * cos(2.0 * theta) - 132.0 * cos(3.0 * theta) + 99.0 * cos(4.0 * theta)) * sin(theta/2.0)) / 64.0 |
| 117 | + elseif m == 2 |
| 118 | + (sqrt(13.0 / π) * cos(theta/2.0)^4 * (1709.0 - 3096.0 * cos(theta) + 2340.0 * cos(2.0 * theta) - 1320.0 * cos(3.0 * theta) + 495.0 * cos(4.0 * theta))) / 256.0 |
| 119 | + elseif m == 3 |
| 120 | + (-3.0 * sqrt(13.0 / π) * cos(theta/2.0)^5 * (-98.0 + 185.0 * cos(theta) - 110.0 * cos(2.0 * theta) + 55.0 * cos(3.0 * theta)) * sin(theta/2.0)) / 32.0 |
| 121 | + elseif m == 4 |
| 122 | + (sqrt(195.0 / (2.0 * π)) * cos(theta/2.0)^6 * (35.0 - 44.0 * cos(theta) + 33.0 * cos(2.0 * theta)) * sin(theta/2.0)^2) / 8.0 |
| 123 | + elseif m == 5 |
| 124 | + (-sqrt(2145.0 / π) * cos(theta/2.0)^7 * (-1.0 + 3.0 * cos(theta)) * sin(theta/2.0)^3) / 2.0 |
| 125 | + elseif m == 6 |
| 126 | + (3.0 * sqrt(715.0 / π) * cos(theta/2.0)^8 * sin(theta/2.0)^4) / 2.0 |
| 127 | + end |
| 128 | + elseif l == 7 |
| 129 | + if m == -7 |
| 130 | + sqrt(15015.0 / (2.0 * π)) * cos(theta/2.0)^5 * sin(theta/2.0)^9 |
| 131 | + elseif m == -6 |
| 132 | + (sqrt(2145.0 / π) * cos(theta/2.0)^4 * (2.0 + 7.0 * cos(theta)) * sin(theta/2.0)^8) / 2.0 |
| 133 | + elseif m == -5 |
| 134 | + (sqrt(165.0 / (2.0 * π)) * cos(theta/2.0)^3 * (93.0 + 104.0 * cos(theta) + 91.0 * cos(2.0 * theta)) * sin(theta/2.0)^7) / 8.0 |
| 135 | + elseif m == -4 |
| 136 | + (sqrt(165.0 / (2.0 * π)) * cos(theta/2.0)^2 * (140.0 + 285.0 * cos(theta) + 156.0 * cos(2.0 * theta) + 91.0 * cos(3.0 * theta)) * sin(theta/2.0)^6) / 16.0 |
| 137 | + elseif m == -3 |
| 138 | + (sqrt(15.0 / (2.0 * π)) * cos(theta/2.0) * (3115.0 + 5456.0 * cos(theta) + 4268.0 * cos(2.0 * theta) + 2288.0 * cos(3.0 * theta) + 1001.0 * cos(4.0 * theta)) * sin(theta/2.0)^5) / 128.0 |
| 139 | + elseif m == -2 |
| 140 | + (sqrt(15.0 / π) * (5220.0 + 9810.0 * cos(theta) + 7920.0 * cos(2.0 * theta) + 5445.0 * cos(3.0 * theta) + 2860.0 * cos(4.0 * theta) + 1001.0 * cos(5.0 * theta)) * sin(theta/2.0)^4) / 512.0 |
| 141 | + elseif m == -1 |
| 142 | + (3.0 * sqrt(5.0 / (2.0 * π)) * cos(theta/2.0) * (1890.0 + 4130.0 * cos(theta) + 3080.0 * cos(2.0 * theta) + 2805.0 * cos(3.0 * theta) + 1430.0 * cos(4.0 * theta) + 1001.0 * cos(5.0 * theta)) * sin(theta/2.0)^3) / 512.0 |
| 143 | + elseif m == 0 |
| 144 | + (3.0 * sqrt(35.0 / π) * cos(theta) * (109.0 + 132.0 * cos(2.0 * theta) + 143.0 * cos(4.0 * theta)) * sin(theta)^2) / 512.0 |
| 145 | + elseif m == 1 |
| 146 | + (3.0 * sqrt(5.0 / (2.0 * π)) * cos(theta/2.0)^3 * (-1890.0 + 4130.0 * cos(theta) - 3080.0 * cos(2.0 * theta) + 2805.0 * cos(3.0 * theta) - 1430.0 * cos(4.0 * theta) + 1001.0 * cos(5.0 * theta)) * sin(theta/2.0)) / 512.0 |
| 147 | + elseif m == 2 |
| 148 | + (sqrt(15.0 / π) * cos(theta/2.0)^4 * (-5220.0 + 9810.0 * cos(theta) - 7920.0 * cos(2.0 * theta) + 5445.0 * cos(3.0 * theta) - 2860.0 * cos(4.0 * theta) + 1001.0 * cos(5.0 * theta))) / 512.0 |
| 149 | + elseif m == 3 |
| 150 | + -(sqrt(15.0 / (2.0 * π)) * cos(theta/2.0)^5 * (3115.0 - 5456.0 * cos(theta) + 4268.0 * cos(2.0 * theta) - 2288.0 * cos(3.0 * theta) + 1001.0 * cos(4.0 * theta)) * sin(theta/2.0)) / 128.0 |
| 151 | + elseif m == 4 |
| 152 | + (sqrt(165.0 / (2.0 * π)) * cos(theta/2.0)^6 * (-140.0 + 285.0 * cos(theta) - 156.0 * cos(2.0 * theta) + 91.0 * cos(3.0 * theta)) * sin(theta/2.0)^2) / 16.0 |
| 153 | + elseif m == 5 |
| 154 | + -(sqrt(165.0 / (2.0 * π)) * cos(theta/2.0)^7 * (93.0 - 104.0 * cos(theta) + 91.0 * cos(2.0 * theta)) * sin(theta/2.0)^3) / 8.0 |
| 155 | + elseif m == 6 |
| 156 | + (sqrt(2145.0 / π) * cos(theta/2.0)^8 * (-2.0 + 7.0 * cos(theta)) * sin(theta/2.0)^4) / 2.0 |
| 157 | + elseif m == 7 |
| 158 | + -(sqrt(15015.0 / (2.0 * π)) * cos(theta/2.0)^9 * sin(theta/2.0)^5) |
| 159 | + end |
| 160 | + elseif l == 8 |
| 161 | + if m == -8 |
| 162 | + sqrt(34034.0 / π) * cos(theta/2.0)^6 * sin(theta/2.0)^10 |
| 163 | + elseif m == -7 |
| 164 | + sqrt(17017.0 / (2.0 * π)) * cos(theta/2.0)^5 * (1.0 + 4.0 * cos(theta)) * sin(theta/2.0)^9 |
| 165 | + elseif m == -6 |
| 166 | + sqrt(255255.0 / π) * cos(theta/2.0)^4 * (1.0 + 2.0 * cos(theta)) * sin(π/4.0 - theta/2.0) * sin(π/4.0 + theta/2.0) * sin(theta/2.0)^8 |
| 167 | + elseif m == -5 |
| 168 | + (sqrt(12155.0 / (2.0 * π)) * cos(theta/2.0)^3 * (19.0 + 42.0 * cos(theta) + 21.0 * cos(2.0 * theta) + 14.0 * cos(3.0 * theta)) * sin(theta/2.0)^7) / 8.0 |
| 169 | + elseif m == -4 |
| 170 | + (sqrt(935.0 / (2.0 * π)) * cos(theta/2.0)^2 * (265.0 + 442.0 * cos(theta) + 364.0 * cos(2.0 * theta) + 182.0 * cos(3.0 * theta) + 91.0 * cos(4.0 * theta)) * sin(theta/2.0)^6) / 32.0 |
| 171 | + elseif m == -3 |
| 172 | + (sqrt(561.0 / (2.0 * π)) * cos(theta/2.0) * (869.0 + 1660.0 * cos(theta) + 1300.0 * cos(2.0 * theta) + 910.0 * cos(3.0 * theta) + 455.0 * cos(4.0 * theta) + 182.0 * cos(5.0 * theta)) * sin(theta/2.0)^5) / 128.0 |
| 173 | + elseif m == -2 |
| 174 | + (sqrt(17.0 / π) * (7626.0 + 14454.0 * cos(theta) + 12375.0 * cos(2.0 * theta) + 9295.0 * cos(3.0 * theta) + 6006.0 * cos(4.0 * theta) + 3003.0 * cos(5.0 * theta) + 1001.0 * cos(6.0 * theta)) * sin(theta/2.0)^4) / 512.0 |
| 175 | + elseif m == -1 |
| 176 | + (sqrt(595.0 / (2.0 * π)) * cos(theta/2.0) * (798.0 + 1386.0 * cos(theta) + 1386.0 * cos(2.0 * theta) + 1001.0 * cos(3.0 * theta) + 858.0 * cos(4.0 * theta) + 429.0 * cos(5.0 * theta) + 286.0 * cos(6.0 * theta)) * sin(theta/2.0)^3) / 512.0 |
| 177 | + elseif m == 0 |
| 178 | + (3.0 * sqrt(595.0 / π) * (210.0 + 385.0 * cos(2.0 * theta) + 286.0 * cos(4.0 * theta) + 143.0 * cos(6.0 * theta)) * sin(theta)^2) / 4096.0 |
| 179 | + elseif m == 1 |
| 180 | + (sqrt(595.0 / (2.0 * π)) * cos(theta/2.0)^3 * (798.0 - 1386.0 * cos(theta) + 1386.0 * cos(2.0 * theta) - 1001.0 * cos(3.0 * theta) + 858.0 * cos(4.0 * theta) - 429.0 * cos(5.0 * theta) + 286.0 * cos(6.0 * theta)) * sin(theta/2.0)) / 512.0 |
| 181 | + elseif m == 2 |
| 182 | + (sqrt(17.0 / π) * cos(theta/2.0)^4 * (7626.0 - 14454.0 * cos(theta) + 12375.0 * cos(2.0 * theta) - 9295.0 * cos(3.0 * theta) + 6006.0 * cos(4.0 * theta) - 3003.0 * cos(5.0 * theta) + 1001.0 * cos(6.0 * theta))) / 512.0 |
| 183 | + elseif m == 3 |
| 184 | + -(sqrt(561.0 / (2.0 * π)) * cos(theta/2.0)^5 * (-869.0 + 1660.0 * cos(theta) - 1300.0 * cos(2.0 * theta) + 910.0 * cos(3.0 * theta) - 455.0 * cos(4.0 * theta) + 182.0 * cos(5.0 * theta)) * sin(theta/2.0)) / 128.0 |
| 185 | + elseif m == 4 |
| 186 | + (sqrt(935.0 / (2.0 * π)) * cos(theta/2.0)^6 * (265.0 - 442.0 * cos(theta) + 364.0 * cos(2.0 * theta) - 182.0 * cos(3.0 * theta) + 91.0 * cos(4.0 * theta)) * sin(theta/2.0)^2) / 32.0 |
| 187 | + elseif m == 5 |
| 188 | + -(sqrt(12155.0 / (2.0 * π)) * cos(theta/2.0)^7 * (-19.0 + 42.0 * cos(theta) - 21.0 * cos(2.0 * theta) + 14.0 * cos(3.0 * theta)) * sin(theta/2.0)^3) / 8.0 |
| 189 | + elseif m == 6 |
| 190 | + (sqrt(255255.0 / π) * cos(theta/2.0)^8 * (-1.0 + 2.0 * cos(theta)) * sin(theta/2.0)^4) * sin(π/4.0 - theta/2.0) * sin(π/4.0 + theta/2.0); |
| 191 | + elseif m == 7 |
| 192 | + -(sqrt(17017.0 / (2.0 * π)) * cos(theta/2.0)^9 * (-1.0 + 4.0 * cos(theta)) * sin(theta/2.0)^5) |
| 193 | + elseif m == 8 |
| 194 | + sqrt(34034.0 / π) * cos(theta/2.0)^10 * sin(theta/2.0)^6 |
| 195 | + end |
| 196 | + end |
| 197 | + # Include complex phase factor |
| 198 | + if m ≠ 0 |
| 199 | + ans = cis(m*phi) * fac |
| 200 | + else |
| 201 | + ans = fac |
| 202 | + end |
| 203 | + end |
| 204 | + |
| 205 | +end # module LAL |
0 commit comments