Coverage for pygeodyn/shear/templib_temp.py: 38%

1from pygeodyn.shear import blackBoxFormula2temp as bbf 

2import numpy as np 

3 

4 

5def openOrCalculateDBrDtForward( 

6 Lb, 

7 Lsv, 

8 Lq, 

9 tmax, 

10 Nu_real_imag_form, 

11 Br_data_matrix, 

12 u_in_real_imag_form_dict, 

13 gauss_thetas, 

14 gauss_weights, 

15 phis, 

16 u_plus_forward_matrix_complex, 

17 u_minus_forward_matrix_complex, 

18 print_text="Opening or calculating dBr/dt forward matrix...", 

19 debug=True, 

20): 

21 

22 

23 u_plus_data_matrices, u_minus_data_matrices = bbf.calculate_u_data_matrices( 

24 Nu_real_imag_form, 

25 u_in_real_imag_form_dict, 

26 gauss_thetas, 

27 phis, 

28 debug=debug, 

29 u_plus_forward_matrix=u_plus_forward_matrix_complex, 

30 u_minus_forward_matrix=u_minus_forward_matrix_complex, 

31 forward_type="complex", 

32 ) 

33 

34 ( 

35 v_plus_matrices, 

36 v_minus_matrices, 

37 ) = bbf.calculate_v_plus_minus( 

38 Lsv, 

39 Nu_real_imag_form, 

40 Br_data_matrix, 

41 u_plus_data_matrices, 

42 u_minus_data_matrices, 

43 gauss_thetas, 

44 gauss_weights, 

45 tmax, 

46 debug=debug, 

47 ) 

48 

49 ( 

50 SV_in_real_imag_form_dict, 

51 Nsv_real_imag_form, 

52 SV_cos_sin_dict, 

53 Nsv_cossin_form, 

54 ) = bbf.build_real_imag_and_cossin_dict_pos(Lsv) 

55 

56 dBr_dt_matrix_from_u_minus_plus = bbf.calculating_SV_from_v_matrix_ulm( 

57 Lsv, 

58 Nsv_real_imag_form, 

59 SV_in_real_imag_form_dict, 

60 v_plus_matrices, 

61 v_minus_matrices, 

62 ) 

63 

64 return ( 

65 v_plus_matrices, 

66 v_minus_matrices, 

67 dBr_dt_matrix_from_u_minus_plus, 

68 ) 

69 

70 

71def flowToMMatrix(u, Lu): 

72 u_m = np.array(u) 

73 n = 0 

74 for l_ in range(1, Lu + 1): 

75 for m_ in range(0, l_ + 1): 

76 u_m[:, n] = u[:, n] * ((-1) ** m_) 

77 n += 1 

78 if m_ != 0: 

79 u_m[:, n] = u[:, n] * ((-1) ** m_) 

80 n += 1 

81 for l_ in range(1, Lu + 1): 

82 for m_ in range(0, l_ + 1): 

83 u_m[:, n] = u[:, n] * ((-1) ** m_) 

84 n += 1 

85 if m_ != 0: 

86 u_m[:, n] = u[:, n] * ((-1) ** m_) 

87 n += 1 

88 

89 return u_m 

90 

91 

92def flowToM(tlmclm, Lu): 

93 tlmclm_m = np.array(tlmclm.reshape((-1, 1))) 

94 n = 0 

95 for l_ in range(1, Lu + 1): 

96 for m_ in range(0, l_ + 1): 

97 tlmclm_m[n, 0] = tlmclm[n] * ((-1) ** m_) 

98 n += 1 

99 if m_ != 0: 

100 tlmclm_m[n, 0] = tlmclm[n] * ((-1) ** m_) 

101 n += 1 

102 for l_ in range(1, Lu + 1): 

103 for m_ in range(0, l_ + 1): 

104 tlmclm_m[n, 0] = tlmclm[n] * ((-1) ** m_) 

105 n += 1 

106 if m_ != 0: 

107 tlmclm_m[n, 0] = tlmclm[n] * ((-1) ** m_) 

108 n += 1 

109 

110 return tlmclm_m 

111 

112 

113def magnFieldFromRaToRaRcM(glm, Lb): 

114 Br_glm_Rc_2 = np.array(glm.reshape((-1, 1))) 

115 Br_glm_Ra_2 = np.array(glm.reshape((-1, 1))) 

116 Br_glm_Ra = np.array(glm.reshape((-1, 1))) 

117 n = 0 

118 for l_ in range(1, Lb + 1): 

119 for m_ in range(0, l_ + 1): 

120 Br_glm_Ra_2[n, 0] = Br_glm_Ra[n, 0] * ((-1) ** m_) 

121 n += 1 

122 if m_ != 0: 

123 Br_glm_Ra_2[n, 0] = Br_glm_Ra[n, 0] * ((-1) ** m_) 

124 n += 1 

125 n = 0 

126 for l_ in range(1, Lb + 1): 

127 for m_ in range(0, l_ + 1): 

128 Br_glm_Rc_2[n, 0] = ( 

129 Br_glm_Ra[n, 0] * ((bbf.Ra / bbf.Rc) ** (l_ + 2)) * ((-1) ** m_) 

130 ) 

131 n += 1 

132 if m_ != 0: 

133 Br_glm_Rc_2[n, 0] = ( 

134 Br_glm_Ra[n, 0] * ((bbf.Ra / bbf.Rc) ** (l_ + 2)) * ((-1) ** m_) 

135 ) 

136 n += 1 

137 

138 return Br_glm_Rc_2, Br_glm_Ra_2, Br_glm_Ra 

139 

140 

141def magnFieldFromRaToRaRcMMatrix(glm, Lsv): 

142 err_Ra_coeffs = np.array(glm) 

143 err_Rc_coeffs_2 = np.array(glm) 

144 err_Ra_coeffs_2 = np.array(glm) 

145 n = 0 

146 for l_ in range(1, Lsv + 1): 

147 for m_ in range(0, l_ + 1): 

148 err_Ra_coeffs_2[:, n] = err_Ra_coeffs[:, n] * ((-1) ** m_) 

149 n += 1 

150 if m_ != 0: 

151 err_Ra_coeffs_2[:, n] = err_Ra_coeffs[:, n] * ((-1) ** m_) 

152 n += 1 

153 n = 0 

154 for l_ in range(1, Lsv + 1): 

155 for m_ in range(0, l_ + 1): 

156 err_Rc_coeffs_2[:, n] = ( 

157 err_Ra_coeffs[:, n] * ((bbf.Ra / bbf.Rc) ** (l_ + 2)) * ((-1) ** m_) 

158 ) 

159 n += 1 

160 if m_ != 0: 

161 err_Rc_coeffs_2[:, n] = ( 

162 err_Ra_coeffs[:, n] * ((bbf.Ra / bbf.Rc) ** (l_ + 2)) * ((-1) ** m_) 

163 ) 

164 n += 1 

165 

166 return err_Rc_coeffs_2, err_Ra_coeffs_2, err_Ra_coeffs 

167 

168 

169def lvFromVMatrix(v_plus_matrices_2, v_minus_matrices_2, Lsv): 

170 lv_plus_matrix = np.array(v_plus_matrices_2) 

171 lv_minus_matrix = np.array(v_minus_matrices_2) 

172 n = 0 

173 for lv in range(1, Lsv): 

174 for mv in range(-lv, lv + 1): 

175 lv_plus_matrix[n, :] = lv_plus_matrix[n, :] * -lv 

176 lv_minus_matrix[n, :] = lv_minus_matrix[n, :] * -lv 

177 n += 1 

178 lv_plus_matrix[n, :] = lv_plus_matrix[n, :] * -lv 

179 lv_minus_matrix[n, :] = lv_minus_matrix[n, :] * -lv 

180 n += 1 

181 

182 n = 0 

183 for l_ in range(1, Lsv): 

184 for m_ in range(-l_, l_ + 1): 

185 lv_plus_matrix[n, :] = ( 

186 lv_plus_matrix[n, :] * (bbf.Ra ** (l_ + 2)) / (bbf.Rc ** (l_ + 2)) 

187 ) 

188 lv_minus_matrix[n, :] = ( 

189 lv_minus_matrix[n, :] * (bbf.Ra ** (l_ + 2)) / (bbf.Rc ** (l_ + 2)) 

190 ) 

191 n += 1 

192 # if m_ != 0: 

193 lv_plus_matrix[n, :] = ( 

194 lv_plus_matrix[n, :] * (bbf.Ra ** (l_ + 2)) / (bbf.Rc ** (l_ + 2)) 

195 ) 

196 lv_minus_matrix[n, :] = ( 

197 lv_minus_matrix[n, :] * (bbf.Ra ** (l_ + 2)) / (bbf.Rc ** (l_ + 2)) 

198 ) 

199 n += 1 

200 

201 return lv_plus_matrix, lv_minus_matrix 

202 

203 

204def calculate_TlmClm_to_UminusUplus_matrix( 

205 Lu, 

206): 

207 ( 

208 u_in_real_imag_form_dict, 

209 Nu_real_imag_form, 

210 u_cos_sin_dict, 

211 Nu_cossin_form, 

212 ) = bbf.build_real_imag_and_cossin_dict_pos(Lu) 

213 

214 Matrix_multiply_by_i = bbf.building_transformation_matrix_multiply_by_i( 

215 Nu_real_imag_form, u_in_real_imag_form_dict 

216 ) 

217 

218 u_Matrix_from_cossin_to_R_I_form = ( 

219 bbf.building_transformation_matrix_from_cossin_to_R_I_form( 

220 Nu_real_imag_form, Nu_cossin_form, u_in_real_imag_form_dict, u_cos_sin_dict 

221 ) 

222 ) 

223 

224 u_Matrix_from_R_I_from_cossin = ( 

225 bbf.building_transformation_matrix_from_R_I_to_cossin_form( 

226 Nu_real_imag_form, Nu_cossin_form, u_in_real_imag_form_dict, u_cos_sin_dict 

227 ) 

228 ) 

229 

230 ######################################### 

231 

232 TtU = np.zeros((Nu_real_imag_form, Nu_real_imag_form)) 

233 CtU = np.zeros((Nu_real_imag_form, Nu_real_imag_form)) 

234 

235 nu = 0 

236 ntc = 0 

237 for lu in range(1, Lu + 1): 

238 temp = np.sqrt(lu * (lu + 1) / 2) 

239 for mu in range(-lu, lu + 1): 

240 CtU[nu, nu] = temp 

241 TtU[nu, nu] = temp 

242 nu += 1 

243 CtU[nu, nu] = temp 

244 TtU[nu, nu] = temp 

245 nu += 1 

246 

247 TCtUp = np.concatenate((-Matrix_multiply_by_i.dot(TtU), CtU), axis=1) 

248 TCtUm = np.concatenate((Matrix_multiply_by_i.dot(TtU), CtU), axis=1) 

249 

250 tc_from_cossin_to_R_I = np.zeros((Nu_real_imag_form * 2, Nu_cossin_form * 2)) 

251 tc_from_cossin_to_R_I[ 

252 :Nu_real_imag_form, :Nu_cossin_form 

253 ] = u_Matrix_from_cossin_to_R_I_form 

254 tc_from_cossin_to_R_I[ 

255 Nu_real_imag_form:, Nu_cossin_form: 

256 ] = u_Matrix_from_cossin_to_R_I_form 

257 

258 TCtUp_from_cossin = TCtUp.dot(tc_from_cossin_to_R_I) 

259 TCtUm_from_cossin = TCtUm.dot(tc_from_cossin_to_R_I) 

260 

261 TCtUmUp = np.concatenate((TCtUm_from_cossin, TCtUp_from_cossin), axis=0) 

262 

263 return TCtUmUp 

264 

265 

266def calculate_QtT_QtC_Lu_form_Q_Lq(Lq): 

267 ( 

268 Q_in_real_imag_form_dict_Lq, 

269 Nq_real_imag_form_Lq, 

270 Q_cos_sin_dict_Lq, 

271 Nq_cossin_form_Lq, 

272 ) = bbf.build_real_imag_and_cossin_dict_pos(Lq, q=True) 

273 

274 Q_Matrix_from_cossin_to_R_I_form = ( 

275 bbf.building_transformation_matrix_from_cossin_to_R_I_form( 

276 Nq_real_imag_form_Lq, 

277 Nq_cossin_form_Lq, 

278 Q_in_real_imag_form_dict_Lq, 

279 Q_cos_sin_dict_Lq, 

280 ) 

281 ) 

282 

283 Q_Matrix_from_R_I_to_cossin_form = ( 

284 bbf.building_transformation_matrix_from_R_I_to_cossin_form( 

285 Nq_real_imag_form_Lq, 

286 Nq_cossin_form_Lq, 

287 Q_in_real_imag_form_dict_Lq, 

288 Q_cos_sin_dict_Lq, 

289 ) 

290 ) 

291 

292 QtC_matrix_Re_Im_Lq = bbf.build_qlm_to_clm_matrix_Re_Im(Lq) 

293 QtT_matrix_part_Re_Im_Lq = bbf.build_qlm_to_tlm_matrix_Re_Im(Lq) 

294 QtT_matrix_zonal_Re_Im_Lq = bbf.build_qlm_to_tlm_zonal_matrix_Re_Im_v2(Lq) 

295 QtT_matrix_real_imag_form_Lq = np.array(QtT_matrix_zonal_Re_Im_Lq) 

296 QtT_matrix_real_imag_form_Lq[:, 2:] += QtT_matrix_part_Re_Im_Lq 

297 

298 temp = np.zeros(QtT_matrix_real_imag_form_Lq.shape) 

299 temp[: Nq_real_imag_form_Lq - 2, 2:] += QtC_matrix_Re_Im_Lq 

300 QtC_matrix_real_imag_form_Lq = temp 

301 

302 return QtT_matrix_real_imag_form_Lq, QtC_matrix_real_imag_form_Lq 

303 

304 

305def calculate_psi1_to_delta_plus_minus_matrix( 

306 Lu, Lq, QtT_matrix_real_imag_form_Lq, QtC_matrix_real_imag_form_Lq 

307): 

308 

309 ( 

310 Q_in_real_imag_form_dict_Lu, 

311 Nq_real_imag_form_Lu, 

312 Q_cos_sin_dict_Lu, 

313 Nq_cossin_form_Lu, 

314 ) = bbf.build_real_imag_and_cossin_dict_pos(Lu, q=True) 

315 

316 ( 

317 Q_in_real_imag_form_dict_Lq, 

318 Nq_real_imag_form_Lq, 

319 Q_cos_sin_dict_Lq, 

320 Nq_cossin_form_Lq, 

321 ) = bbf.build_real_imag_and_cossin_dict_pos(Lq, q=True) 

322 

323 ( 

324 u_in_real_imag_form_dict, 

325 Nu_real_imag_form, 

326 u_cos_sin_dict, 

327 Nu_cossin_form, 

328 ) = bbf.build_real_imag_and_cossin_dict_pos(Lu) 

329 

330 Matrix_multiply_by_i = bbf.building_transformation_matrix_multiply_by_i( 

331 Nu_real_imag_form, u_in_real_imag_form_dict 

332 ) 

333 

334 ( 

335 ulm_plus_from_qlm_matrix_Lq, 

336 ulm_minus_from_qlm_matrix_Lq, 

337 ) = bbf.calculating_ulm_plus_minus_from_qlm( 

338 Lu, 

339 Nu_real_imag_form, 

340 QtC_matrix_real_imag_form_Lq[:, :], 

341 QtT_matrix_real_imag_form_Lq[:, :], 

342 Matrix_multiply_by_i, 

343 ) 

344 

345 ulm_minus_plus_from_qlm_matrix = np.zeros( 

346 (Nu_real_imag_form * 2, Nq_real_imag_form_Lq) 

347 ) 

348 ulm_minus_plus_from_qlm_matrix[:Nu_real_imag_form, :] = ulm_minus_from_qlm_matrix_Lq 

349 ulm_minus_plus_from_qlm_matrix[Nu_real_imag_form:, :] = ulm_plus_from_qlm_matrix_Lq 

350 

351 psi1_to_psi2_data = bbf.calculate_psi1_to_psi2(Lu) 

352 # psi1_to_psi2_data = calculate_psi1_to_psi2(Lu) 

353 

354 psi1_to_psi2_matrix = bbf.calculate_psi1_to_psi2_matrix(Lu, psi1_to_psi2_data) 

355 

356 psi1_to_psi3_matrix = psi1_to_psi2_matrix + np.eye(Nq_real_imag_form_Lu) * 6 

357 

358 ( 

359 psi3_to_s_tilda_matrix_plus, 

360 psi3_to_s_tilda_matrix_minus, 

361 ) = bbf.calculate_s_tilda_from_qlm_tilda(Lu) 

362 

363 psi1_to_s_tilda_matrix_plus_temp = psi3_to_s_tilda_matrix_plus.dot( 

364 psi1_to_psi3_matrix 

365 ) 

366 psi1_to_s_tilda_matrix_minus_temp = psi3_to_s_tilda_matrix_minus.dot( 

367 psi1_to_psi3_matrix 

368 ) 

369 psi1_to_s_tilda_matrix_plus = psi1_to_s_tilda_matrix_plus_temp 

370 psi1_to_s_tilda_matrix_minus = psi1_to_s_tilda_matrix_minus_temp 

371 

372 ( 

373 psi1_to_s_hat_matrix_plus, 

374 psi1_to_s_hat_matrix_minus, 

375 ) = bbf.calculate_s_hat_from_qlm(Lu) 

376 # psi1_to_s_hat_matrix_plus[: (Lu-1) * (Lu -1 + 2) * 2 + 2, :(Lu-1) * (Lu -1 + 2) * 2 + 2] = psi1_to_s_hat_matrix_plus_temp 

377 # psi1_to_s_hat_matrix_minus[: (Lu-1) * (Lu -1 + 2) * 2 + 2, :(Lu-1) * (Lu -1 + 2) * 2 + 2] = psi1_to_s_hat_matrix_minus_temp 

378 

379 psi1_to_delta_plus_temp = psi1_to_s_tilda_matrix_plus + psi1_to_s_hat_matrix_plus 

380 psi1_to_delta_minus_temp = psi1_to_s_tilda_matrix_minus + psi1_to_s_hat_matrix_minus 

381 psi1_to_delta_plus = np.zeros((Nq_real_imag_form_Lu, Nq_real_imag_form_Lq)) 

382 psi1_to_delta_minus = np.zeros((Nq_real_imag_form_Lu, Nq_real_imag_form_Lq)) 

383 psi1_to_delta_plus = psi1_to_delta_plus_temp[:, :Nq_real_imag_form_Lq] 

384 psi1_to_delta_minus = psi1_to_delta_minus_temp[:, :Nq_real_imag_form_Lq] 

385 

386 return psi1_to_delta_plus, psi1_to_delta_minus 

387 

388 

389def calculate_u_minus_plus_from_u_minus(Lu): 

390 ( 

391 u_in_real_imag_form_dict, 

392 Nu_real_imag_form, 

393 u_cos_sin_dict, 

394 Nu_cossin_form, 

395 ) = bbf.build_real_imag_and_cossin_dict_pos(Lu) 

396 

397 u_plus_from_u_minus = np.zeros((Nu_real_imag_form, Nu_real_imag_form)) 

398 

399 for lu in range(1, Lu + 1): 

400 for mu in range(-lu, lu + 1): 

401 temp = (-1) ** mu 

402 

403 u_pos_reim_plus = list(u_in_real_imag_form_dict.keys())[ 

404 list(u_in_real_imag_form_dict.values()).index((lu, mu, "r")) 

405 ] 

406 u_pos_reim_minus = list(u_in_real_imag_form_dict.keys())[ 

407 list(u_in_real_imag_form_dict.values()).index((lu, -mu, "r")) 

408 ] 

409 

410 u_plus_from_u_minus[u_pos_reim_plus, u_pos_reim_minus] = temp 

411 u_plus_from_u_minus[u_pos_reim_plus + 1, u_pos_reim_minus + 1] = -temp 

412 

413 u_minus_plus_from_u_minus = np.concatenate( 

414 (np.eye(Nu_real_imag_form), u_plus_from_u_minus), axis=0 

415 ) 

416 

417 return u_minus_plus_from_u_minus 

418 

419 

420def calculate_u_minus_from_u_minus_plus(Lu): 

421 u_minus_plus_from_u_minus = calculate_u_minus_plus_from_u_minus(Lu) 

422 u_minus_from_u_minus_plus = np.linalg.inv( 

423 u_minus_plus_from_u_minus.T.dot(u_minus_plus_from_u_minus) 

424 ).dot(u_minus_plus_from_u_minus.T) 

425 return u_minus_from_u_minus_plus 

426 

427 

428def from_Ra_to_Rc(matr, Lsv): 

429 matr_res = np.array(matr) 

430 n = 0 

431 for l_ in range(1, Lsv): 

432 for m_ in range(-l_, l_ + 1): 

433 matr_res[n, :] = matr[n, :] * (bbf.Ra ** (l_ + 2)) / (bbf.Rc ** (l_ + 2)) 

434 n += 1 

435 matr_res[n, :] = matr[n, :] * (bbf.Ra ** (l_ + 2)) / (bbf.Rc ** (l_ + 2)) 

436 n += 1 

437 return matr_res 

438 

439 

440def calculate_corelation(Aij, Bij, thetas, tmax, pmax): 

441 sin_theta_map = np.zeros((tmax, pmax)) 

442 

443 for th in range(tmax): 

444 sin_theta = np.sin(thetas[th]) 

445 for ph in range(pmax): 

446 sin_theta_map[th, ph] = sin_theta 

447 

448 Aij_map = np.array(Aij.real.reshape((tmax, pmax))) 

449 Bij_map = np.array(Bij.real.reshape((tmax, pmax))) 

450 

451 covar_ij = np.zeros((tmax, pmax)) 

452 

453 norm = 0 

454 for th in range(tmax): 

455 for ph in range(pmax): 

456 if th > 0 and th < tmax - 1: 

457 dtheta = np.abs(thetas[th + 1] - thetas[th - 1]) / 2 

458 else: 

459 if th == 0: 

460 dtheta = np.abs(thetas[th + 1] - thetas[th]) 

461 else: 

462 dtheta = np.abs(thetas[th - 1] - thetas[th]) 

463 temp_A2 = dtheta * (Aij_map[th, ph] ** 2) * sin_theta_map[th, ph] 

464 temp_B2 = dtheta * (Bij_map[th, ph] ** 2) * sin_theta_map[th, ph] 

465 

466 norm += np.sqrt(temp_A2 * temp_B2) 

467 

468 norm = norm * 2 * np.pi / pmax 

469 

470 value = 0 

471 for th in range(tmax): 

472 for ph in range(pmax): 

473 if th > 0 and th < tmax - 1: 

474 dtheta = np.abs(thetas[th + 1] - thetas[th - 1]) / 2 

475 else: 

476 if th == 0: 

477 dtheta = np.abs(thetas[th + 1] - thetas[th]) 

478 else: 

479 dtheta = np.abs(thetas[th - 1] - thetas[th]) 

480 temp_AB = ( 

481 dtheta * (Aij_map[th, ph] * Bij_map[th, ph]) * sin_theta_map[th, ph] 

482 ) 

483 

484 value += temp_AB 

485 

486 value = value * 2 * np.pi / pmax 

487 

488 return value / norm 

489 

490 

491def calculate_normalised_misfit(Aij, Bij, thetas, tmax, pmax): 

492 sin_theta_map = np.zeros((tmax, pmax)) 

493 

494 for th in range(tmax): 

495 sin_theta = np.sin(thetas[th]) 

496 for ph in range(pmax): 

497 sin_theta_map[th, ph] = sin_theta 

498 

499 Aij_map = np.array(Aij.real.reshape((tmax, pmax))) 

500 Bij_map = np.array(Bij.real.reshape((tmax, pmax))) 

501 

502 covar_ij = np.zeros((tmax, pmax)) 

503 

504 norm = 0 

505 for th in range(tmax): 

506 for ph in range(pmax): 

507 if th > 0 and th < tmax - 1: 

508 dtheta = np.abs(thetas[th + 1] - thetas[th - 1]) / 2 

509 else: 

510 if th == 0: 

511 dtheta = np.abs(thetas[th + 1] - thetas[th]) 

512 else: 

513 dtheta = np.abs(thetas[th - 1] - thetas[th]) 

514 temp_A2 = dtheta * (Aij_map[th, ph] ** 2) * sin_theta_map[th, ph] 

515 temp_B2 = dtheta * (Bij_map[th, ph] ** 2) * sin_theta_map[th, ph] 

516 

517 norm += np.sqrt(temp_A2 * temp_B2) 

518 # norm += temp_B2 

519 

520 norm = norm * 2 * np.pi / pmax 

521 

522 value = 0 

523 for th in range(tmax): 

524 for ph in range(pmax): 

525 if th > 0 and th < tmax - 1: 

526 dtheta = np.abs(thetas[th + 1] - thetas[th - 1]) / 2 

527 else: 

528 if th == 0: 

529 dtheta = np.abs(thetas[th + 1] - thetas[th]) 

530 else: 

531 dtheta = np.abs(thetas[th - 1] - thetas[th]) 

532 temp_AB = ( 

533 dtheta 

534 * ((Aij_map[th, ph] - Bij_map[th, ph]) ** 2) 

535 * sin_theta_map[th, ph] 

536 ) 

537 

538 value += temp_AB 

539 

540 value = value * 2 * np.pi / pmax 

541 

542 return value / norm 

543 

544 

545def calculate_norm(Aij, thetas, tmax, pmax): 

546 

547 sin_theta_map = np.zeros((tmax, pmax)) 

548 for th in range(tmax): 

549 sin_theta = np.sin(thetas[th]) 

550 for ph in range(pmax): 

551 sin_theta_map[th, ph] = sin_theta 

552 

553 temp = 0 

554 for th in range(tmax): 

555 temp_theta = 0 

556 for ph in range(pmax): 

557 temp_theta += (Aij[th, ph] ** 2) * sin_theta_map[th, ph] 

558 temp += temp_theta 

559 

560 return np.sqrt(temp) 

561 

562 

563def calculate_spectrum_re_im_flow(coef, L): 

564 spectrum = np.zeros((L,)) 

565 iter = 0 

566 for l in range(1, L + 1): 

567 temp = 0 

568 for m in range(-l, l + 1): 

569 c1 = coef[iter] 

570 c2 = coef[iter + 1] 

571 if m != 0: 

572 c1 = c1 / np.sqrt(2) 

573 c2 = c2 / np.sqrt(2) 

574 temp += c1**2 

575 temp += c2**2 

576 iter += 2 

577 spectrum[l - 1] += temp * (l + 1) 

578 for l in range(1, L + 1): 

579 temp = 0 

580 for m in range(-l, l + 1): 

581 c1 = coef[iter] 

582 c2 = coef[iter + 1] 

583 if m != 0: 

584 c1 = c1 / np.sqrt(2) 

585 c2 = c2 / np.sqrt(2) 

586 temp += c1**2 

587 temp += c2**2 

588 iter += 2 

589 spectrum[l - 1] += temp * (l + 1) 

590 

591 return spectrum 

592 

593 

594def calculate_spectrum_re_im_SV2_pos(coef, L): 

595 spectrum = np.zeros((L,)) 

596 iter = 0 

597 for l in range(1, L + 1): 

598 temp = 0 

599 for m in range(0, l + 1): 

600 if m == 0: 

601 temp += coef[iter] ** 2 

602 iter += 1 

603 else: 

604 temp += coef[iter] ** 2 

605 temp += coef[iter + 1] ** 2 

606 iter += 2 

607 spectrum[l - 1] += temp * (l + 1) 

608 

609 return spectrum 

610 

611 

612def calculate_spectrum_re_im_SV2_pos_all_m(coef, L): 

613 spectrum = np.zeros((L,)) 

614 iter = 0 

615 for l in range(1, L + 1): 

616 temp = 0 

617 for m in range(-l, l + 1): 

618 c = coef[iter] 

619 if m != 0: 

620 c = c / np.sqrt(2) 

621 temp += c**2 

622 

623 c = coef[iter + 1] 

624 if m != 0: 

625 c = c / np.sqrt(2) 

626 temp += c**2 

627 iter += 2 

628 

629 spectrum[l - 1] += temp * (l + 1) 

630 

631 return spectrum