Skip to content

Update source flags handling #1620

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
mmccrackan wants to merge 1 commit into
base: master
Choose a base branch
from
Draft

Update source flags handling #1620

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
32 changes: 21 additions & 11 deletions sotodlib/coords/planets.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -151,6 +151,7 @@ def get_scan_q(tod, planet, boresight_offset=None, refq=None):
# Get reference elevation...
el = np.median(tod.boresight.el[::10])
az = np.median(tod.boresight.az[::10])
roll = np.median(tod.boresight.roll[::10])
t = (tod.timestamps[0] + tod.timestamps[-1]) / 2
if isinstance(planet, (list, tuple)):
_, ra, dec = planet
Expand All @@ -159,15 +160,15 @@ def get_scan_q(tod, planet, boresight_offset=None, refq=None):
else:
planet = SlowSource.for_named_source(planet, t)

def scan_q_model(t, az, el, planet):
def scan_q_model(t, az, el, roll, planet):
csl = so3g.proj.CelestialSightLine.az_el(
t, az, el, weather='typical', site='so')
t, az, el, roll=roll, weather='typical', site='so')
ra0, dec0 = planet.pos(t)
return csl.Q, ~so3g.proj.quat.rotation_lonlat(ra0, dec0) * csl.Q * q_xieta

def distance(p):
dt, daz = p
q, qnet = scan_q_model(t + dt, az + daz, el, planet)
q, qnet = scan_q_model(t + dt, az + daz, el, roll, planet)
lon, lat, phi = so3g.proj.quat.decompose_lonlat(qnet)
return 90 * coords.DEG - lat

Expand All @@ -177,7 +178,7 @@ def distance(p):
if warnflag != 0:
logger.warning('Source-scan solver failed to converge or otherwise '
f'complained! warnflag={warnflag}')
q, qnet = scan_q_model(t+p[0], az+p[1], el, planet)
q, qnet = scan_q_model(t+p[0], az+p[1], el, roll, planet)
psi = so3g.proj.quat.decompose_xieta(qnet)[2][0]
ra, dec = planet.pos(t+p[0])
rot = ~so3g.proj.quat.rotation_lonlat(ra, dec, psi)
Expand Down Expand Up @@ -591,19 +592,28 @@ def compute_source_flags(tod=None, P=None, mask=None, wrap=None,
the center of the circle, all in degrees.

"""
if P is None:
logger.info('Getting Projection Matrix ...')
P, X = get_scan_P(tod, center_on, res=res, comps='T')
shape, wcs = tuple(P.geom)
if shape[0] * shape[1] > max_pix:
raise ValueError(f'Mask map too large: {shape}')

if isinstance(mask, str):
# Assume it's a filename, and file is simple columns of (x, y,
# radius) in deg. (Deprecated!)
mask = [{'xyr': list(map(float, line.split()))}
for line in open(mask)]

if P is None:
x, y, r = mask['xyr']
shape = 3 * np.ceil(r * coords.DEG / res).astype(int)
# ensure odd dimensions
if shape % 2 == 0:
shape += 1
if shape**2 > max_pix:
raise ValueError(f'Mask map too large: {shape}')
logger.info('Getting Projection Matrix ...')
P, X = get_scan_P(tod, center_on, res=res, comps='T')

geom = enmap.geometry(np.array((0, 0)), shape=(shape, shape),
proj='tan', res=(res, -res))
P.geom = enmap.Geometry(*geom)
shape, wcs = tuple(P.geom)

mask_map = P.zeros()
_add_to_mask(mask, mask_map)
a = P.from_map(mask_map)
Expand Down