get_scan_q sometimes fail to estimate timing when planet just pass pointing elevation.

I found get_scan_q sometimes fail to estimate timing when planet just pass pointing elevation.
I put a demonstration code below. Once you modify a context file path, I think this code can work.
In this code, pointing elevation is about 40deg, but we got 47 deg from get_scan_q. This is weird we should get about 40deg if my understanding correctly.

If I had some mistake, please let me know.
Thank you in advance.

import numpy as np
import matplotlib.pyplot as plt
from sotodlib import core
from sotodlib.coords.planets import SlowSource, get_scan_q
from so3g.proj import quat
import so3g

def calc_pointing_planet_radec_azel(aman, planet_str, boresight_offset, s = slice(0, -1, 1)):
    tm = np.mean(aman.timestamps)
    ts = aman.timestamps[s]
    az = aman.boresight.az[s]
    el = aman.boresight.el[s]
    bore = {}
    bore['ts'] = ts
    bore['az'] = az
    bore['el'] = el

    # calc pointing
    q_xieta = quat.rotation_xieta(boresight_offset[0], boresight_offset[1])
    csl_rd = so3g.proj.CelestialSightLine.az_el(ts, az, el, weather="typical")
    csl_ae = so3g.proj.CelestialSightLine.for_horizon(ts, az, el, weather="typical")
    q_bore2rd = csl_rd.Q*q_xieta
    q_bore2ae = csl_ae.Q*q_xieta
    
    th_rd, phi_rd, psi_rd = quat.decompose_iso(q_bore2rd)
    th_ae, phi_ae, psi_ae = quat.decompose_iso(q_bore2ae)

    po = {}
    po['ra'] = phi_rd
    po['dec'] = (np.pi/2 - th_rd)%np.pi
    po['az'] = (2*np.pi - phi_ae)%(2*np.pi)
    po['el'] = (np.pi/2 - th_ae)%np.pi


    # calc planet pointing
    planet = SlowSource.for_named_source(planet_str, tm)
    print(f'Planet positin in RA/Dec [deg]: {np.rad2deg(planet.ra)}/{np.rad2deg(planet.dec)}')
    
    ra_pla, dec_pla = planet.pos(ts)
    q_pla = quat.rotation_lonlat(ra_pla, dec_pla)
    
    q_pla_ae = q_bore2ae*~q_bore2rd * q_pla
    th_pl_ae, phi_pl_ae, psi_pl_ae = quat.decompose_iso(q_pla_ae) # return theta, phi, psi
    el_pla = (np.pi/2 - th_pl_ae)%np.pi
    az_pla = (2*np.pi - phi_pl_ae)%(2*np.pi)

    pl={}
    pl['ra'] = ra_pla
    pl['dec'] = dec_pla
    pl['az'] = az_pla
    pl['el'] = el_pla

    return bore, po, pl

def demo1(ctx):
    iobs_id = 'obs_1709928536_satp3_1111111'
    iws = 'ws2'
    iband = 'f090'
    planet_str = 'jupiter'
    boresight_offset = (np.deg2rad(-12), np.deg2rad(-6))

    print(f'load data: {iobs_id}, {iws}, {iband}, planet = {planet_str}')
    print(f'Boresight_offset, xi/eta [deg] = {np.rad2deg(boresight_offset[0])}/{np.rad2deg(boresight_offset[1])}')

    meta = ctx.get_meta(iobs_id, dets={'wafer.type': 'OPTC', 'wafer_slot': iws, 'wafer.bandpass': iband, })
    meta.restrict('dets', meta.dets.vals[:5])
    aman = ctx.get_obs(meta)
    s = slice(0,-1,100)
    bore, po, pl = calc_pointing_planet_radec_azel(aman, planet_str, boresight_offset, s = s)


    # get index from naive elevation match method
    idx = np.argmin(np.abs(pl['el'] - np.mean(po['el'])))

    # get index from get_scan_q
    X = get_scan_q(aman, planet_str, boresight_offset=boresight_offset, refq=None)
    idx2 = np.argmin(np.abs(bore['ts'] - X['timestamp']))
    # compare planet elevation for both method
    print('')
    print(f'Elevation from focal plane position with boresight_offset = {np.mean(np.rad2deg(po["el"]))} deg')
    print(f'naive elevation match method: index={idx}, elecation of {planet_str} = {np.rad2deg(pl["el"][idx])} deg')
    print(f'get_scan_q: index={idx2}, elecation of {planet_str} = {np.rad2deg(pl["el"][idx2])} deg')
    print(f'Found elevation difference: naive - get_scan_q = {np.rad2deg(pl["el"][idx]) - np.rad2deg(pl["el"][idx2])} deg')

    # plot
    s = slice(0,-1,1)
    fig, ax = plt.subplots(1,2,figsize = (12,6))
    ax[0].plot(np.rad2deg(po['az'][s]), np.rad2deg(po['el'][s]), '.', label = 'Pointing')
    ax[0].plot(np.rad2deg(pl['az'][s]), np.rad2deg(pl['el'][s]), '.', label = 'Jupiter')
    ax[0].set_ylabel('Elevation [deg]')
    ax[0].set_xlabel('Azimuth [deg]')
    ax[0].legend()

    ax[1].plot(np.rad2deg(po['ra'][s]), np.rad2deg(po['dec'][s]), '.', label = 'Pointing')
    ax[1].plot(np.rad2deg(pl['ra'][s]), np.rad2deg(pl['dec'][s]), '*', markersize = 12, label = 'Jupiter')
    ax[1].set_ylabel('Dec [deg]')
    ax[1].set_xlabel('RA [deg]')
    ax[1].legend()

    fig.suptitle(f'{iobs_id}, {iws}, boresight_offset = {np.rad2deg(boresight_offset)} deg', y = 1.02)
    fig.savefig('demo1.png', bbox_inches = 'tight')
    plt.close()

if __name__ == '__main__':
    ctx = core.Context('/scratch/gpfs/SIMONSOBS/so/tracked/metadata/satp3/contexts/use_this_local.yaml')
    demo1(ctx)

Output of this code

Elevation from focal plane position with boresight_offset = 40.72375654178508 deg
naive elevation match method: index=15043, elecation of jupiter = 40.72358800151032 deg
get_scan_q: index=9899, elecation of jupiter = 46.71914332515219 deg
Found elevation difference: naive - get_scan_q = -5.995555323641867 deg