fix psi computation in ecliptic->galactic pointing rotation

[mj-gomes]

Hi,

This PR fixes the error that I found recently and that I briefly mentioned in #542.

I noticed that scanning a map with a non ideal HWP resulted in bizarre output maps (which looked like some weird angle rotation was happening somewhere) compared to when doing the scanning in ecliptic coordinates.

I (+ Eric, + one LLM) managed to find the bug. Currently, when we convert pointings from ecliptic to galactic coordinates, we are only rotating the orientation angle psi considering the rotation of the north pole, and ignoring the fact that the whole reference frame changes (the definition of both north and east directions), and since psi is defined between the north and east directions, the rotation of psi was incorrect.

The reason to why this only showed up when the hwp was non ideal is unclear, but it can come from the fact that systematic effects rotate with psi independently of the hwp rotation angle (rho), while for an ideal hwp, the rotation is always (rho-psi), which perhaps cancels out the effects from an incorrect psi.

Now, for each sample, we first obtain the north and east directions in the ecliptic frame with theta and phi, then we get the orientation vector in this basis, and we rotate both the basis and the vector, and then we obtain psi in galactic coordinates doing the inverse operation from vector to angle in the new reference frame.

I am now getting the same output maps using galactic or ecliptic coordinates.

EDIT: apparently we have to change the expected outputs for some examples in the docs, I'll do that tomorrow.

[mj-gomes]

The error in the map_scanning documentation seems reasonable due to the changes, but the one in beam_convolution (https://github.com/litebird/litebird_sim/actions/runs/29433612982/job/87414564734?pr=543#step:7:115) looks a bit strange.

Namely due to the big difference in the previously expected output TOD and the one we get now:

Expected:
    113.796
    116.186
    118.567
    120.935
    123.300
    125.655
    128.000
    130.335
    132.661
    134.975
Got:
    -113.305
    -101.235
    -87.724
    -73.623
    -59.721
    -46.686
    -35.032
    -25.097
    -17.047
    -10.875

Can someone with more expertise than me in beam convolution (@paganol ?) see if this is expected with the change in psi calculation, or if there is something still wrong ?

[paganol]

Hi @mj-gomes, if I compare the outcome of healpy.Rotator with your fix I got weird results. I'm doing this:

import numpy as np
import healpy as hp
from litebird_sim.coordinates import _rotate_coordinates_and_orientation_e2g_for_one_sample as lbs_rot

rot = hp.Rotator(coord=["E", "G"])

theta_ecl, phi_ecl, psi_ecl = np.random.uniform(0, np.pi), np.random.uniform(0, 2*np.pi), np.random.uniform(0, 2*np.pi)
alpha = rot.angle_ref((theta_ecl, phi_ecl))

psi_gal_hp = (psi_ecl + alpha) % (2 * np.pi)

_, _, psi_gal_lbs = lbs_rot(theta_ecl, phi_ecl, psi_ecl)
psi_gal_lbs %= 2 * np.pi

print(np.isclose(psi_gal_hp, psi_gal_lbs, atol=1e-6))

For master I'm getting perfect match but not for this PR. I'm not saying that a fix is necessary for the HWP case, but this seems to break the agreement with healpy

[mj-gomes]

Hi @mj-gomes, if I compare the outcome of healpy.Rotator with your fix I got weird results. I'm doing this:

import numpy as np
import healpy as hp
from litebird_sim.coordinates import _rotate_coordinates_and_orientation_e2g_for_one_sample as lbs_rot

rot = hp.Rotator(coord=["E", "G"])

theta_ecl, phi_ecl, psi_ecl = np.random.uniform(0, np.pi), np.random.uniform(0, 2*np.pi), np.random.uniform(0, 2*np.pi)
alpha = rot.angle_ref((theta_ecl, phi_ecl))

psi_gal_hp = (psi_ecl + alpha) % (2 * np.pi)

_, _, psi_gal_lbs = lbs_rot(theta_ecl, phi_ecl, psi_ecl)
psi_gal_lbs %= 2 * np.pi

print(np.isclose(psi_gal_hp, psi_gal_lbs, atol=1e-6))

For master I'm getting perfect match but not for this PR. I'm not saying that a fix is necessary for the HWP case, but this seems to break the agreement with healpy

Hum, ok. I will look at this and try to find the problem as soon as I can.

[mj-gomes]

Hi @paganol, you are right! After some headaches, I realized my solution was actually not correct, and was yielding a psi_gal = psi_ecl.

The solution already implemented seems correct, indeed.

Also, I tried reproducing the error again and realized it seems to only exist when using the PTEP IMo in specific (even though I still need to dore more tests to confirm this and understand why). I will keep looking at this to try to figure out if there is indeed a problem or not.

Meanwhile, I will close this PR, since this solution is not only wrong but also not needed. Thank you for pointing it out!