<jupyter_start><jupyter_text>Math for astropy spherical_offsets_toset up a new spherical coordinate centered at the reference point<jupyter_code>## cartesian coordinates

# pivot point a
x_a = np.cos(dec_a) * np.cos(ra_a)
y_a = np.sin(ra_a) * np.cos(dec_a)
z_a = np.sin(dec_a)

# point b
x_b = np.cos(dec_b) * np.cos(ra_b)
y_b = np.sin(ra_b) * np.cos(dec_b)
z_b = np.sin(dec_b)

# new north pole
x_p = np.cos(dec_p) * np.cos(ra_p)
y_p = np.sin(ra_p) * np.cos(dec_p)
z_p = np.sin(dec_p)


# normal vector of the plane 1 where point a, origin and new north pole reside
A1 = (y_a * z_p) - (z_a * y_p)
B1 = x_a * z_p - z_a * x_p
C1 = y_a * x_p - x_a * y_p

# normal vector of the plane 2 where point b, origin and new north pole reside
A2 = (y_b * z_p) - (z_b * y_p)
B2 = x_b * z_p - z_b * x_p
C2 = y_b * x_p - x_b * y_p

# the intersection angle of plane 1 and plane 2 -> new ra
labc = (Math.sqrt(Math.pow(A1, 2) + Math.pow(B1, 2) + Math.pow(C1, 2))
     * Math.sqrt(Math.pow(A2, 2) + Math.pow(B2, 2) + Math.pow(C2, 2)))
cos_theta = (A1*A2 + B1*B2 + C1*C2)/labc
theta = np.arccos(cos_theta).to(u.deg) * np.sign(ra_b-ra_a)

# the intersection angle of the plane 1 normal and the vector origin-b -> new dec
cos_phi = x_p*x_b + y_p*y_b + z_p*z_b
phi = 90*u.deg - np.arccos(cos_phi).to(u.deg) 

print('my',theta.to('rad'),phi.to('rad'))

a = SkyCoord(ra_a, dec_a , frame='icrs')
b = SkyCoord(ra_b , dec_b , frame='icrs')
dra, ddec = a.spherical_offsets_to(b)
print('astropy',dra.to('rad'),ddec.to('rad'))<jupyter_output>my -0.25069188491905714 rad -0.5099574742224875 rad
astropy -0.250692rad -0.509957rad<jupyter_text>Calculate a precise midpoint and PA on the projected plane centered on the midpoint<jupyter_code>x_a = np.cos(ra_a)*np.cos(dec_a)
x_b = np.cos(ra_b)*np.cos(dec_b)
psedo_x_mid = (x_a + x_b)/2

y_a = np.sin(ra_a)*np.cos(dec_a)
y_b = np.sin(ra_b)*np.cos(dec_b)
psedo_y_mid = (y_a + y_b)/2

z_a = np.sin(dec_a)
z_b = np.sin(dec_b)
psedo_z_mid = (z_a + z_b)/2

scale = psedo_x_mid**2 + psedo_y_mid**2 + psedo_z_mid**2
scale = 1 / np.sqrt(scale)
x_mid = scale * psedo_x_mid
y_mid = scale * psedo_y_mid
z_mid = scale * psedo_z_mid
ra_mid = np.arctan2(y_mid,x_mid).to(u.deg)
dec_mid = np.arctan( z_mid/np.sqrt(x_mid**2 + y_mid**2) ).to(u.deg)
print('a',ra_a.to(u.deg).value, dec_a.value)
print('b',ra_b.to(u.deg).value, dec_b.to(u.deg).value)
print('mid',ra_mid, dec_mid)

p_a = np.array([x_a,y_a,z_a])
p_b = np.array([x_b,y_b,z_b])
p_m = np.array([x_mid,y_mid,z_mid])

# projected plane centered on the midpoint
# (x-x_mid)x_mid+(y-y_mid)y_mid+(z-z_mid)z_mid=0
lm = np.dot(p_m,p_m)
t_a = lm/np.dot(p_a,p_m)
t_b = lm/np.dot(p_b, p_m)

# projected point for a and b and projected vector
cp_a = t_a * p_a
cp_b = t_b * p_b
vec_proj = cp_b - cp_a

# vector on the projected plane towards North
e = np.array([-np.cos(90*u.deg-dec_mid)*np.cos(ra_mid),
                np.cos(90*u.deg-dec_mid)*np.sin(ra_mid),np.sin(90*u.deg-dec_mid)])


ll = np.sqrt(np.dot(vec_proj,vec_proj))
cos_theta = np.dot(e,vec_proj)/ll
theta =  np.arccos(cos_theta)*u.radian
print('PA (intersection angle only; No direction):',theta.to(u.deg))<jupyter_output>a 180.00416666666663 30.0
b 167.49999999999997 0.0
mid 173.3014231030648 deg 15.085223349097275 deg
PA (intersection angle only; No direction): 155.21469941876074 deg<jupyter_text>Math for Keck https://www2.keck.hawaii.edu/inst/common/offset.phpApproximation in a very small FoV<jupyter_code>ra_kmid = (ra_a + ra_b)/2
dec_kmid = (dec_a + dec_b)/2

east_offset = (ra_b-ra_a) * np.cos(dec_kmid)
north_offset = dec_b-dec_a
pa_keck = np.arctan(east_offset/north_offset).to(u.deg)

print('E and N offset:', east_offset.to(u.arcsec),north_offset.to(u.arcsec))
print('PA (E of N):', pa_keck)<jupyter_output>E and N offset: -43481.2arcsec -108000arcsec
PA (E of N): 21.92987608386546 deg
{'fraction_non_alphanumeric': 0.16085017055890843, 'fraction_numerical': 0.05694043558121228, 'mean_word_length': 4.573099415204679, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 12, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 2, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '487046', 'license': 'MIT', 'n_tokens_mistral': 1979, 'n_tokens_neox': 1832, 'n_words': 361, 'path': 'assets/code/Astropy_spherical_offset.ipynb', 'repo_name': 'LittleLin1999/littlelin1999.github.io', 'size': '3811'}