# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

import numpy as np


# + [markdown] tags=[]
# # Planetary Orbits
#
# We have the followings relations:
#
# \begin{equation}
# l_2 v_2 = l_1 v_1 \qquad \text{Simplification of Kepler's Law} \\
# v_2^2 - \frac{2GM}{v_1 l_1}v_2 - \left[ v_1^2 - \frac{2GM}{l_1}\right] = 0 \qquad \text{Second order Polynomial}
# \end{equation}
# where $m$ is the mass of the planet, $M = 1.9891 \times 10^{30}\,\text{Kg}$ is the mass of the Sun, and $G = 6.6738 \times 10^{-11}\,\text{m}^3 \text{Kg}^{-1}s^{-2}$ is the universal constant of gravitation.
#
# ---
# The values $v_1$, $l_1$ and $v_2$, $l_2$ need to be the extreme values of the elliptical orbit.
# $l_1$ is the closest point to the sun, it's the perihelion.
# $l_2$ is the farthest point to the sun, it's the aphelion.
#
# To solve this problem we will use Classes to better structure this problem.
# To find the value of $v_2$ we only need to solve a second order polynomial.
# -

class planetSolarSystem:
    # Global Constants, we need to declare them global so we can access
    # them inside the methods of the class.
    global G, M
    G = 6.6738e-11
    M = 1.9891e30

    def __init__(self, v1, l1):
        """Creates the planet values.
        
        :param m: float mass of the planet.
        :param v1: float velocity of the closest point to the sun.
        :param l1: float lowest radius of the orbit.
        
        l1 must be the perihelion and l2 will be the aphelion.
        """
        self.v1 = v1
        self.l1 = l1
        self.v2 = self.calculate_v2()
        self.l2 = self.v1 * self.l1 / self.v2
        self.a = self.semi_major_axis()
        self.b = self.semi_minor_axis()
        self.T = self.orbital_period()
        self.e = self.orbital_eccentricity()
        
    def calculate_v2(self):
        """Calculates v2
        
        V2 is the velocity in the point of the biggest radius of the orbit.
        """
        a = 1
        b = -2 * G * M / (self.v1 * self.l1)
        c = -(self.v1**2 - 2 * G * M / self.l1)
        aux_1 = -b/2
        aux_2 = np.sqrt(b**2/4 - c)
        solution = np.array([aux_1 - aux_2, aux_1 + aux_2])
        return np.amin(solution)
    
    def semi_major_axis(self):
        """Returns the value of the major axis.
        """
        return (self.l1 + self.l2)/2
        
    def semi_minor_axis(self):
        """Returns the value of the minor axis.
        """
        return np.sqrt(self.l1 * self.l2)
    
    def orbital_period(self):
        """Returns the value of the orbital period.
        """
        return 2 * np.pi * self.a * self.b / (self.l1 * self.v1)
        
    def orbital_eccentricity(self):
        """Returns the value of the orbital eccentricity.
        """
        return (self.l2 - self.l1)/(self.l2 + self.l1)
    
    def print_relevant_info(self):
        """Used to show relevant information about the celestial body.
        """
        print("Semi-major axis: {:.3f}".format(self.a))
        print("Semi-minor axis: {:.3f}".format(self.b))
        print("Orbital period: {:.3f}".format(self.T))
        print("Orbital eccentricity: {:.3f}\n".format(self.e))
        pass


# +
v1earth = 3.0287e4
l1earth = 1.4710e11
v1halley = 5.4529e4
l1halley = 8.7830e10
earth = planetSolarSystem(v1earth, l1earth)
halley = planetSolarSystem(v1halley, l1halley)

# print the information asked by the exercise
earth.print_relevant_info()
halley.print_relevant_info()
earthT = earth.T
halleyT = halley.T

# convert to years
earthT = earthT / 60 / 60 / 24 / 365.25
halleyT = halleyT / 60 / 60 / 24 / 365.25
print("{:.1f} years".format(earthT))
print("{:.1f} years".format(halleyT))

StarcoderdataJupyterScriptsDedupFiltered