<commit_before># celestial_mechanics.py
#   celestial mechanics utilities for exoplanet ephemerides
#
# intellectual property:
#   Copyright 2009 <NAME>.  All rights reserved.	 Watch this space
#   for an open-source license to be inserted later.
#
# comments:
#   - Written for clarity, not speed.  The code is intended to be human-
#     readable.
#
# bugs:
#   - Need to make Fourier expansion functions.
#	- What to do if e is close to 1.0 in eccentric_anomaly
#

from math import pi
import unittest
import sys

import numpy
import numpy as np
from numpy import *
from scipy.special import jn
import matplotlib.pyplot as plt

from astrometry.util.starutil_numpy import *

ihat = array([1.,0.,0.])
jhat = array([0.,1.,0.])
khat = array([0.,0.,1.])
default_tolerance = 1e-15 # (radians) don't set to zero or less
default_maximum_iteration = 1000 # should never hit this limit!
default_order = 32
default_K = 1.0

(Equinox, Solstice, EclipticPole) = ecliptic_basis()

c_au_per_yr = 63239.6717 # google says

def norm(x):
	return np.sqrt(np.dot(x, x))

def deg2rad(x):
	return x * pi/180.
	#return radians(x)

GM_sun = 2.9591310798672560E-04 #AU^3/d^2

def orbital_elements_to_ss_xyz(E, observer, light_travel=True):
	(a,e,i,Omega,pomega,M,GM) = E
	# ugh, it's hard to be units-agnostic.
	# we just assert here so we have to think about this!
	assert(GM == GM_sun)
	# orbital angular velocity  [radians/yr]
	meanfrequency = np.sqrt(GM / a**3)
	# Correct for light-time delay.
	dM = 0.
	lastdM = dM
	for ii in range(100):
		(x,v) = phase_space_coordinates_from_orbital_elements(
			a,e,i,Omega,pomega,M-dM,GM)
		dx = (x - observer)
		if not light_travel:
			break
		r = norm(dx)
		travel = r / c_au_per_yr
		dM = travel * meanfrequency
		if abs(lastdM - dM) < 1e-12:
			break
		lastdM = dM
	#print 'niters', ii
	return x,dx

def orbital_elements_to_xyz(E, observer, light_travel=True):
	(x,dx) = orbital_elements_to_ss_xyz(E, observer, light_travel)
	dx /= norm(dx)
	edx = dx[0] * Equinox + dx[1] * Solstice + dx[2] * EclipticPole
	return edx

# E = (a,e,i,Omega,pomega,M, GM)
# observer = 3-vector
# light_travel: correct for light travel time?
# Returns Ra,Dec in degrees.
def orbital_elements_to_radec(E, observer, light_travel=True):
	xyz = orbital_elements_to_xyz(E, observer, light_travel)
	return xyztoradec(xyz)

# convert orbital elements into vectors in the plane of the orbit.
def orbital_vectors_from_orbital_elements(i, Omega, pomega):
	ascendingnodevector = np.cos(Omega) * ihat + np.sin(Omega) * jhat
	tmpydir= np.cross(khat, ascendingnodevector)
	zhat= np.cos(i) * khat - np.sin(i) * tmpydir
	tmpydir= np.cross(zhat, ascendingnodevector)
	xhat= np.cos(pomega) * ascendingnodevector + np.sin(pomega) * tmpydir
	yhat = np.cross(zhat, xhat)
	return (xhat, yhat, zhat)

def position_from_orbital_vectors(xhat, yhat, a, e, M):
	E = eccentric_anomaly_from_mean_anomaly(M, e)
	cosE = np.cos(E)
	sinE = np.sin(E)
	b = a*np.sqrt(1. - e**2)
	x =  a * (cosE - e)  * xhat + b * sinE        * yhat
	return x

# convert orbital elements to phase-space coordinates
#  a       - semi-major axis (length units)
#  e       - eccentricity
#  i       - inclination (rad)
#  Omega   - longitude of ascending node (rad)
#  pomega  - argument of periapsis (rad)
#  M       - mean anomaly (rad)
#  GM      - Newton's constant times central mass (length units cubed over time units squared)
#  return  - (x,v)
#            position, velocity (length units, length units per time unit)
def phase_space_coordinates_from_orbital_elements(a, e, i, Omega, pomega, M, GM):
	(xhat, yhat, zhat) = orbital_vectors_from_orbital_elements(i, Omega, pomega)
	dMdt = np.sqrt(GM / a**3)
	E = eccentric_anomaly_from_mean_anomaly(M, e)
	cosE = np.cos(E)
	sinE = np.sin(E)
	dEdt = 1.0 / (1.0 - e * cosE) * dMdt
	b = a*np.sqrt(1. - e**2)
	x =  a * (cosE - e)  * xhat + b * sinE        * yhat
	v = -a * sinE * dEdt * xhat + b * cosE * dEdt * yhat
	return (x, v)

class UnboundOrbitError(ValueError):
	pass

def potential_energy_from_position(x, GM):
	return -1. * GM / norm(x)

def energy_from_phase_space_coordinates(x, v, GM):
	return 0.5 * np.dot(v, v) + potential_energy_from_position(x, GM)

# convert phase-space coordinates to orbital elements
#  x       - position (3-vector, length units)
#  v       - velocity (3-vector, length units per time unit)
#  GM      - Newton's constant times central mass (length units cubed over time units squared)
#  return  - (a, e, i, Omega, pomega, M)
#          - see "phase_space_coordinates" for definitions
def orbital_elements_from_phase_space_coordinates(x, v, GM):
	energy = energy_from_phase_space_coordinates(x, v, GM)
	if energy > 0:
		raise UnboundOrbitError('orbital_elements_from_phase_space_coordinates: Unbound orbit')

	angmom = np.cross(x, v)
	zhat = angmom / norm(angmom)
	evec = np.cross(v, angmom) / GM - x / norm(x)
	e = norm(evec)
	if e == 0:
		# by convention:
		xhat = np.cross(jhat, zhat)
		xhat /= norm(xhat)
	else:
		xhat = evec / e
	yhat = np.cross(zhat, xhat)
	a = -0.5 * GM / energy
	dMdt = np.sqrt(GM / a**3)
	i = np.arccos(angmom[2] / norm(angmom))
	if i == 0:
		Omega = 0.0
	else:
		Omega = np.arctan2(angmom[1], angmom[0]) + 0.5 * pi
		if Omega < 0:<commit_msg>make 1-d-ness of norm() explicit<commit_after># celestial_mechanics.py
#   celestial mechanics utilities for exoplanet ephemerides
#
# intellectual property:
#   Copyright 2009 <NAME>.  All rights reserved.	 Watch this space
#   for an open-source license to be inserted later.
#
# comments:
#   - Written for clarity, not speed.  The code is intended to be human-
#     readable.
#
# bugs:
#   - Need to make Fourier expansion functions.
#	- What to do if e is close to 1.0 in eccentric_anomaly
#

from math import pi
import unittest
import sys

import numpy
import numpy as np
from numpy import *
from scipy.special import jn
import matplotlib.pyplot as plt

from astrometry.util.starutil_numpy import *

ihat = array([1.,0.,0.])
jhat = array([0.,1.,0.])
khat = array([0.,0.,1.])
default_tolerance = 1e-15 # (radians) don't set to zero or less
default_maximum_iteration = 1000 # should never hit this limit!
default_order = 32
default_K = 1.0

(Equinox, Solstice, EclipticPole) = ecliptic_basis()

c_au_per_yr = 63239.6717 # google says

def norm1d(x):
	assert(len(x.shape) == 1)
	return np.sqrt(np.dot(x, x))

def deg2rad(x):
	return x * pi/180.
	#return radians(x)

GM_sun = 2.9591310798672560E-04 #AU^3/d^2

def orbital_elements_to_ss_xyz(E, observer, light_travel=True):
	(a,e,i,Omega,pomega,M,GM) = E
	# ugh, it's hard to be units-agnostic.
	# we just assert here so we have to think about this!
	assert(GM == GM_sun)
	# orbital angular velocity  [radians/yr]
	meanfrequency = np.sqrt(GM / a**3)
	# Correct for light-time delay.
	dM = 0.
	lastdM = dM
	for ii in range(100):
		(x,v) = phase_space_coordinates_from_orbital_elements(
			a,e,i,Omega,pomega,M-dM,GM)
		dx = (x - observer)
		if not light_travel:
			break
		r = norm1d(dx)
		travel = r / c_au_per_yr
		dM = travel * meanfrequency
		if abs(lastdM - dM) < 1e-12:
			break
		lastdM = dM
	#print 'niters', ii
	return x,dx

def orbital_elements_to_xyz(E, observer, light_travel=True):
	(x,dx) = orbital_elements_to_ss_xyz(E, observer, light_travel)
	dx /= norm1d(dx)
	edx = dx[0] * Equinox + dx[1] * Solstice + dx[2] * EclipticPole
	return edx

# E = (a,e,i,Omega,pomega,M, GM)
# observer = 3-vector
# light_travel: correct for light travel time?
# Returns Ra,Dec in degrees.
def orbital_elements_to_radec(E, observer, light_travel=True):
	xyz = orbital_elements_to_xyz(E, observer, light_travel)
	return xyztoradec(xyz)

# convert orbital elements into vectors in the plane of the orbit.
def orbital_vectors_from_orbital_elements(i, Omega, pomega):
	ascendingnodevector = np.cos(Omega) * ihat + np.sin(Omega) * jhat
	tmpydir= np.cross(khat, ascendingnodevector)
	zhat= np.cos(i) * khat - np.sin(i) * tmpydir
	tmpydir= np.cross(zhat, ascendingnodevector)
	xhat= np.cos(pomega) * ascendingnodevector + np.sin(pomega) * tmpydir
	yhat = np.cross(zhat, xhat)
	return (xhat, yhat, zhat)

def position_from_orbital_vectors(xhat, yhat, a, e, M):
	E = eccentric_anomaly_from_mean_anomaly(M, e)
	cosE = np.cos(E)
	sinE = np.sin(E)
	b = a*np.sqrt(1. - e**2)
	x =  a * (cosE - e)  * xhat + b * sinE        * yhat
	return x

# convert orbital elements to phase-space coordinates
#  a       - semi-major axis (length units)
#  e       - eccentricity
#  i       - inclination (rad)
#  Omega   - longitude of ascending node (rad)
#  pomega  - argument of periapsis (rad)
#  M       - mean anomaly (rad)
#  GM      - Newton's constant times central mass (length units cubed over time units squared)
#  return  - (x,v)
#            position, velocity (length units, length units per time unit)
def phase_space_coordinates_from_orbital_elements(a, e, i, Omega, pomega, M, GM):
	(xhat, yhat, zhat) = orbital_vectors_from_orbital_elements(i, Omega, pomega)
	dMdt = np.sqrt(GM / a**3)
	E = eccentric_anomaly_from_mean_anomaly(M, e)
	cosE = np.cos(E)
	sinE = np.sin(E)
	dEdt = 1.0 / (1.0 - e * cosE) * dMdt
	b = a*np.sqrt(1. - e**2)
	x =  a * (cosE - e)  * xhat + b * sinE        * yhat
	v = -a * sinE * dEdt * xhat + b * cosE * dEdt * yhat
	return (x, v)

class UnboundOrbitError(ValueError):
	pass

def potential_energy_from_position(x, GM):
	return -1. * GM / norm1d(x)

def energy_from_phase_space_coordinates(x, v, GM):
	return 0.5 * np.dot(v, v) + potential_energy_from_position(x, GM)

# convert phase-space coordinates to orbital elements
#  x       - position (3-vector, length units)
#  v       - velocity (3-vector, length units per time unit)
#  GM      - Newton's constant times central mass (length units cubed over time units squared)
#  return  - (a, e, i, Omega, pomega, M)
#          - see "phase_space_coordinates" for definitions
def orbital_elements_from_phase_space_coordinates(x, v, GM):
	energy = energy_from_phase_space_coordinates(x, v, GM)
	if energy > 0:
		raise UnboundOrbitError('orbital_elements_from_phase_space_coordinates: Unbound orbit')

	angmom = np.cross(x, v)
	zhat = angmom / norm1d(angmom)
	evec = np.cross(v, angmom) / GM - x / norm1d(x)
	e = norm1d(evec)
	if e == 0:
		# by convention:
		xhat = np.cross(jhat, zhat)
		xhat /= norm1d(xhat)
	else:
		xhat = evec / e
	yhat = np.cross(zhat, xhat)
	a = -0.5 * GM / energy
	dMdt = np.sqrt(GM / a**3)
	i = np.arccos(angmom[2] / norm1d(angmom))
	if i == 0:
		Omega = 0.0
	else:
		Omega = np.arctan2(angmom[1], angmom[0]) + 0.5 * pi
		if Omega < 0:
{'fraction_non_alphanumeric': 0.1226145038167939, 'fraction_numerical': 0.020324427480916032, 'mean_word_length': 3.6831992850759607, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3800867', 'n_tokens_mistral': 4395, 'n_tokens_neox': 4098, 'n_words': 1133}