diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 9356cb1de9c96876528416199389066580433524..0d9324822c52e8149d001fa245eaadd2c221eb8f 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -1,23 +1,60 @@ stages: - build + - doc - test + +# Change pip's cache directory to be inside the project directory since we can +# only cache local items. +variables: + PIP_CACHE_DIR: "$CI_PROJECT_DIR/.cache/pip" +cache: + paths: + - .cache/pip + - venv/ + +before_script: + - echo "deb http://ftp.debian.org/debian jessie main" > /etc/apt/sources.list + - echo "deb http://ftp.debian.org/debian testing main" >> /etc/apt/sources.list + - echo 'APT::Default-Release "jessie";' >> /etc/apt/apt.conf.d/00local + - apt-get -qq update + - apt-cache policy python3 + - apt-get remove -y binutils + - apt-get -t testing install -y python3 + - apt-get -qq -t testing install -y python3-venv python3-pip python3-tk + - python3 -V + - python3 -m venv venv + - . venv/bin/activate + - pip3 install -r requirements.txt + + +pages: + script: + - echo "Documentation" + - pip3 install sphinx_rtd_theme + - sphinx-apidoc -o docs/source . + - sphinx-build -b html docs/source docs/build + - cp -a docs/build/ public/ + + artifacts: + paths: + - public + + + job_build: - before_script: - - echo "deb http://ftp.debian.org/debian jessie main" > /etc/apt/sources.list - - echo "deb http://ftp.debian.org/debian testing main" >> /etc/apt/sources.list - - echo 'APT::Default-Release "jessie";' >> /etc/apt/apt.conf.d/00local - - apt-get -qq update - - apt-cache policy python3 - - apt-get remove -y binutils - - apt-get -t testing install -y python3 - - apt-get -qq -t testing install -y python3-venv python3-pip python3-tk - - python3 -V - - python3 -m venv venv - - . venv/bin/activate - - pip3 install -r requirements.txt stage: build script: - echo "build" - ./etagement.py + + +job_test: + + stage: test + script: + - echo "test" + - chmod u+rwx Test_Engineering_tools.py + - chmod u+rwx Engineering_tools.py + - ./Test_Engineering_tools.py \ No newline at end of file diff --git a/Data_structure.py b/Data_structure.py index f1741d28e99a8594646a4e1a8e6bae95675b9be8..2c887fd6033ee5e3c27d0a556eeee4aed90fa5ea 100755 --- a/Data_structure.py +++ b/Data_structure.py @@ -1,3 +1,4 @@ +#!/usr/bin/env python # -*- coding: utf-8 -*- """ Created on Sun May 12 19:51:18 2019 @@ -7,79 +8,83 @@ Created on Sun May 12 19:51:18 2019 import Engineering_tools as et +RETROGRADE_ORBIT = 90 + class Rocket(): - """ - """ - - def __init__(self,stage1,stage2,payload,launch_coordinates): - """ - stage1 and stage2 are instances of the class Stage, + + def __init__(self, stage1, stage2, payload, launch_coordinates): + """stage1 and stage2 are instances of the class Stage, payload is an instance of the class Payload, - launch_coordiantes = {latitude(deg),logitude(deg)} + launch_coordinates = {latitude(deg), longitude(deg)} + + :param stage1: [description] + :type stage1: [type] + :param stage2: [description] + :type stage2: [type] + :param payload: [description] + :type payload: Payload + :param launch_coordinates: [description] + :type launch_coordinates: [type] """ self.stage1 = stage1 self.stage2 = stage2 self.payload = payload self.injection_velocity = et.final_target_vel(self.payload.perigee_altitude, - self.payload.target_orbit['a']) + self.payload.target_orbit.a) self.launch_coordinates = launch_coordinates - if self.payload.target_orbit['i'] > 90: # retrograade orbit + if self.payload.target_orbit.i > RETROGRADE_ORBIT: # retrograade orbit self.initial_velocity = et.initial_vel(self.launch_coordinates['latitude'],retrograde=True) else: - self.initial_velocity = et.initial_vel(self.launch_coordinates['latitude']) + self.initial_velocity = et.initial_vel(self.launch_coordinates['longitude']) - def target_Vc(self,losses): - """ - Compute the target critical velocity - losses : velocity losses approximation (m.s-1) + def target_Vc(self, p_losses): + """Compute the target critical velocity losses : velocity losses approximation (m.s-1) + + :param p_losses: [description] + :type p_losses: [type] """ - self.target_Vc = self.injection_velocity-self.initial_velocity+losses + self.target_Vc = self.injection_velocity - self.initial_velocity + p_losses + def __str__(self): + return 'Rocket' class Stage(): - """ - """ - def __init__(self,stage_number,number_of_engines): + def __init__(self, stage_number, number_of_engines = 9, avionics_mass=0): """ stage_number : 1 or 2 """ self.stage_number = stage_number - - - avionics_mass = 0 + self.number_of_engines = number_of_engines + self.avionics_mass = 0 - number_of_engines = 9 - - def tank_mass(self,) + def tank_mass(self): + pass + + def __str__(self): + return 'Stage' class Payload(): - """ - """ - - def __init__(self,mass,target_orbit): - """ - mass in kg - target_orbit = {'a','e','i','capital_omega','small_omega','nu'} - where: - a : semi_major axis (m) - e : eccentricity (adimentional) - i : inclination (deg) - capital_omega : argument of the ascending node (deg) - small_omega : argument of the perigee (deg) - nu : true anomaly (deg) - """ + + def __init__(self, mass, target_orbit): + """ Payload + :param mass: [description] + :type mass: [type] + :param target_orbit: [description] + :type target_orbit: TargetOrbit + """ + self.mass = mass self.target_orbit = target_orbit - self.perigee_altitude = target_orbit['a']*(1-target_orbit['e']) + self.perigee_altitude = target_orbit.a * (1 - target_orbit.e) -# def propagate_orbit(): -# """ -# """ + def __str__(self): + return 'Payload' + class Engine(): """ @@ -87,10 +92,65 @@ class Engine(): cycle_type='' - def __init__(self,oxidiser,fuel,thermo_cycle,exit_section_area): + def __init__(self, oxidiser, fuel, thermo_cycle, exit_section_area): """ """ self.oxidiser = oxidiser self.fuel = fuel self.thermo_cycle = thermo_cycle - self.exit_section_area = exit_section_area \ No newline at end of file + self.exit_section_area = exit_section_area + + def __str__(self): + return 'Engine' + +class TargetOrbit(): + + def __init__(self, a, e, i, capital_omega, small_omega, nu): + """[summary] + + :param a: semi_major axis (m) + :type a: [type] + :param e: eccentricity (adimentional) + :type e: [type] + :param i: inclination (deg) + :type i: [type] + :param capital_omega: argument of the ascending node (deg) + :type capital_omega: [type] + :param small_omega: argument of the perigee (deg) + :type small_omega: [type] + :param nu: [description] + :type nu: [type] + """ + self.a = a + self.e = e + self.i = i + self.capital_omega = capital_omega + self.small_omega = small_omega + self.nu = nu + + def __str__(self): + return ("TargetOrbit object\n" + "Semi major axis : {}\n" + "Exccentricity : {}\n" + "Inclination : {}\n" + "Argument of ascending node : {}\n" + "Argument of periapsis : {}\n" + "nu ? : {}\n").format(self.a,self.e,self.i,self.capital_omega,self.small_omega,self.nu) + + +if __name__ == '__main__': + + v_semiMajorAxis = 400000 # 400 km + v_eccentricity = 15 + v_inclination = 20 + v_ascendingNode = 5 + v_perigee = 5 + v_nu = 0 + + v_targetOrbit = TargetOrbit(v_semiMajorAxis, v_eccentricity, v_inclination, v_ascendingNode, v_perigee, v_nu) + print(v_targetOrbit) + v_payload = Payload(400, v_targetOrbit) + + + Rocket(Stage(1), Stage(2), v_payload, {'latitude':47,'longitude':2}) + \ No newline at end of file diff --git a/Engineering_tools.py b/Engineering_tools.py index 21d3397d92b20ca51def1c26af741ed51814870d..ce63feb554439a85644d839cf6c17bf145e7b820 100755 --- a/Engineering_tools.py +++ b/Engineering_tools.py @@ -1,4 +1,5 @@ -# -*- coding: utf-8 -*- +#!/usr/bin/env python +# - * - coding: utf-8 - * - """ Created on Sun May 12 20:08:52 2019 @@ -7,16 +8,26 @@ Created on Sun May 12 20:08:52 2019 import numpy as np from numpy import linalg as LA - + + + +# Constants +T_SIDE = 86164.10 # sidereal day (s) +MU = 3.9860e14 # earth standard gravitational parameter (=G * M_earth) (m3.s-2) +R_T = 6378e3 # earth radius (m) +G0 = 9.81 # standard gravity (m.s-2) +R = 287 # specific gas constant (J.Kg-1.K-1) + + def initial_vel(latitude,retrograde=False): """ Initial velocity (m.s-1) depending on the launch point latitude (°) - """ - T_side = 86164.10 # sidereal day (s) - R_T = 6378e3 # earth raduis (m) + :param latitude: the launch point latitude in ° + :return: the initial velocity in m.s-1 + """ - initial_velocity = (2*np.pi/T_side)*R_T*np.cos(np.deg2rad(np.abs(latitude))) + initial_velocity = (2 * np.pi / T_SIDE) * R_T * np.cos(np.deg2rad(np.abs(latitude))) if retrograde: return -initial_velocity else: @@ -24,84 +35,77 @@ def initial_vel(latitude,retrograde=False): def final_target_vel(h, a): """ - inputs: - h: injection altitude (m) - a: orbit semi major axis (m) - output: - orbit absolute injection velocity (m.s-1) - """ - - mu = 3.9860e14 # earth standard gravitatinal parameter (=G*M_earth) (m3.s-2) - R_T = 6378e3 # earth raduis (m) - return np.sqrt(2*mu*((1/(R_T+h))-(1/(2*a)))) + .. math:: \\sqrt {2 \\cdot MU \\cdot (\\frac {1}{R_T + h} - \\frac {1}{2 * a})} -def staging_optim(isp1,isp2,eps1,eps2,Vc,Mp): + :param h: injection altitude (m) + :type h: float + :param a: orbit semi major axis (m) + :type a: float + :return: orbit absolute injection velocity (m.s-1) + :rtype: float """ - Optimise mass loading in the case of two stage rocket - inputs: - isp1, isp2 : specific impulse of the 1st and 2nd stage (s) - eps1, eps2 : structural coefficient of the 1st and 2nd stage (adimentional) - Vc : target critical velocity (m.s-1) - Mp : payload mass (Kg) - outputs: - M1 : mass of the 1st stage - M2 : mass of the 2nd stage - - NB: - structure mass - structural coefficient = -------------------------------- - structure mass + propellant mass - - ctritical velocity = final velocity - initial velocity + losses + + v_square = 2 * MU * ((1 / (R_T + h)) - (1 / (2 * a))) + print('{} > {}'.format((R_T + h), (2 * a))) + assert v_square >= 0, 'error: square root only on positiv values' + return np.sqrt(v_square) + +def staging_optim(isp1,isp2,eps1,eps2,Vc,Mp): + + """Optimise mass loading in the case of two stage rocket + + .. math:: structural_coefficient = \\frac {strcture_mass} {structure_mass + propellant_mass} + + :param isp1: specific impulse of the first stage (s) + :type isp1: float + :param isp2: specific impulse of the second stage (s) + :type isp2: float + :param eps1: structural coefficient of the first stage + :type eps1: float + :param eps2: structural coefficient of the second stage + :type eps2: float + :param Vc: target critical velocity (m.s-1) + :type Vc: float + :param Mp: payload mass (Kg) + :type Mp: float + :return: mass of the first and second stage (Kg) + :rtype: float """ -# Notation : -# -# Mt = total mass -# q1 = M1/Mt -# q2 = M2/Mt -# qp = Mp/Mt -# omega = qp + q2 - - g0 = 9.81 # standard gravity (m.s-2) - omega_tab=np.linspace(0.001,0.9,100) # array of omega to test qp_tab=[] for omega in omega_tab: - q2 = omega+eps1*(1-omega) - q2 = np.power(q2,isp1/isp2)*np.exp(Vc/(g0*isp2)) - q2 = 1/q2 - q2 = 1-q2 - q2 *= omega/(1-eps2) - qp = omega-q2 + q2 = omega + eps1 * (1 - omega) + q2 = np.power(q2,isp1 / isp2) * np.exp(Vc / (G0 * isp2)) + q2 = 1 / q2 + q2 = 1 - q2 + q2 *= omega / (1 - eps2) + qp = omega - q2 qp_tab.append(qp) qp=max(qp_tab) omega=omega_tab[qp_tab.index(qp)] - q2=omega-qp - q1=1-q2-qp - M1=(Mp/qp)*q1 - M2=(Mp/qp)*q2 + q2=omega - qp + q1=1 - q2 - qp + M1=(Mp / qp) * q1 + M2=(Mp / qp) * q2 return M1, M2 def standard_atmosphere(altitude): - """ - Compute atmospheric properties according to the International Standaard Atmosphere (ISA) model - input: - altitude (m) - outputs: - temperature (K) - pressure (Pa) - density (Kg.m-3) - """ + """ Compute atmospheric properties according to the International \ + Standard Atmosphere (ISA) model - assert altitude < 84852 , str(altitude)+' is higher than max altitude of the ISA model' + :param altitude: altitude in meter + :type altitude: float + :return: T, P, rho : temperature (K), pressure (Pa), density (Kg.m-3) + :rtype: float, float, float + """ - g0 = 9.81 # standard gravity (m.s-2) - r = 287 # specific gas constant (J.Kg-1.K-1) # ISA model + assert altitude < 84852 , str(altitude) + ' is higher than max altitude of the ISA model' + h0_tab = np.array([0, 11e3, 20e3, 32e3, 47e3, 51e3, 71e3, 84852]) # base altitude (m) a_tab = np.array([-6.5e-3, 0, 1e-3, 2.8e-3, 0, -2.8e-3, -2.0e-3]) # lapse rate (K.m-1) # the base properties of each range are tabulated to reduce computational cost @@ -109,123 +113,129 @@ def standard_atmosphere(altitude): P0_tab = np.array([101325, 22632, 5474.9, 868.02, 110.91, 66.939, 3.9564]) # base pressure (Pa) index = 0 - while altitude > h0_tab[index+1]: - index += 1 + while altitude > h0_tab[index + 1]: + index += 1 a = a_tab[index] T0 = T0_tab[index] P0 = P0_tab[index] - T = T0+a*(altitude-h0_tab[index]) # Toussaint equation + T = T0 + a * (altitude - h0_tab[index]) # Toussaint equation if T == T0 : # isohermal region - P = P0*np.exp(-g0*(altitude-h0_tab[index])/(r*T)) + P = P0 * np.exp(-G0 * (altitude - h0_tab[index]) / (R * T)) else : # gradient region - P = P0*np.power(T/T0,-g0/(a*r)) + P = P0 * np.power(T/T0, -G0/(a*R)) - rho = P/(r*T) + rho = P / (R * T) return T,P,rho - -def laval_nozzle(Pt,Tt,Pamb,Ath,Ae,gamma,r,tolerance=1e-5): - """ - Compute the thrust assuming a started nozzle - inputs: - Pt : total pressure in the combustion chambre (Pa) - Tt : total temperature in the combustion chambre (K) - Pamb : ambient pressure (Pa) - Ath : throat section area (m2) - Ae : exit section area (m2) - gamma : ratio of specific heats (adimentional) - r : specific gas constant (J.Kg-1.K-1) - tolerance : - outputs: - Thrust (N) - Specific impulse (s) - - N.B : Thrust = mass flow rate * exit velocity + (exit pressure - ambiant pressure) * exit area + +def laval_nozzle(Pt,Tt,Pamb,Ath,Ae,gamma,R,tolerance=1e-5): + """Compute the thrust assuming a started nozzle + + .. math:: Thrust = massFlowRate \cdot exitVelocity + (exitPressure - ambiantPressure) \cdot exitArea + + :param Pt: total pressure in the combustion chambre (Pa) + :type Pt: [type] + :param Tt: total temperature in the combustion chambre (K) + :type Tt: [type] + :param Pamb: ambient pressure (Pa) + :type Pamb: [type] + :param Ath: throat section area (m2) + :type Ath: [type] + :param Ae: exit section area (m2) + :type Ae: [type] + :param gamma: ratio of specific heats (adimentional) + :type gamma: [type] + :param R: specific gas constant (J.Kg-1.K-1) + :type R: [type] + :param tolerance: [description], defaults to 1e-5 + :type tolerance: [type], optional + :return: Thrust (N), Specific impulse (s) + :rtype: [type] """ + assert gamma > 1, 'gamma must be > 1' - g0 = 9.81 # standard gravity (m.s-2) - - A_ratio = lambda M,gamma: (1/M)*np.power((2/(gamma+1))*(1+((gamma-1)/2)*M*M),(gamma+1)/(2*gamma-2)) + A_ratio = lambda M,gamma: (1 / M) * np.power((2 / (gamma + 1)) * (1 + ((gamma - 1) / 2) * M * M), (gamma + 1) / (2 * gamma - 2)) + # compute the nozzle exit section area corresponding to the adapted regime - Madapted = np.sqrt((2/(gamma-1))*(np.power(Pamb/Pt,(1-gamma)/gamma)-1)) - Aadapted = Ath*A_ratio(Madapted,gamma) + Madapted = np.sqrt((2 / (gamma - 1)) * (np.power(Pamb / Pt,(1 - gamma) / gamma) - 1)) + Aadapted = Ath * A_ratio(Madapted,gamma) if Aadapted > Ae: # the nozzle is under-expanded # compute the exit Mach number iteratively using dichotomy algorithm - A_ratio_target = Ae/Ath + A_ratio_target = Ae / Ath Me_upper = Madapted Me_lower = 1 - Me = (Me_upper+Me_lower)/2 - error = A_ratio_target-A_ratio(Me,gamma) + Me = (Me_upper + Me_lower) / 2 + error = A_ratio_target - A_ratio(Me,gamma) while np.abs(error) > tolerance: if error > 0: Me_lower = Me else: Me_upper = Me - Me = (Me_upper+Me_lower)/2 - error = A_ratio_target-A_ratio(Me,gamma) - Te = Tt/(1+((gamma-1)/2)*Me*Me) # exit static temperature - Pe = Pt*np.power(Tt/Te,gamma/(1-gamma)) # exit static pressure + Me = (Me_upper + Me_lower) / 2 + error = A_ratio_target - A_ratio(Me,gamma) + Te = Tt / (1 + ((gamma - 1) / 2) * Me * Me) # exit static temperature + Pe = Pt * np.power(Tt / Te,gamma / (1 - gamma)) # exit static pressure elif Aadapted < Ae: # the nozzle is over-expanded Mshock = 1 step = 1 - A_ratio_target = Ae/Ath + A_ratio_target = Ae / Ath Me_upper = 1 Me_lower = 0 - Me = (Me_upper+Me_lower)/2 - error = A_ratio_target-A_ratio(Me,gamma) + Me = (Me_upper + Me_lower) / 2 + error = A_ratio_target - A_ratio(Me,gamma) while np.abs(error) > tolerance: if error < 0: Me_lower = Me else: Me_upper = Me - Me = (Me_upper+Me_lower)/2 - error = A_ratio_target-A_ratio(Me,gamma) - Te = Tt/(1+((gamma-1)/2)*Me*Me) - Pe = Pt*np.power(Tt/Te,gamma/(1-gamma)) - while (Pe-Pamb)*1e-5 > tolerance: + Me = (Me_upper + Me_lower) / 2 + error = A_ratio_target - A_ratio(Me,gamma) + Te = Tt / (1 + ((gamma - 1) / 2) * Me * Me) + Pe = Pt * np.power(Tt / Te,gamma / (1 - gamma)) + while (Pe - Pamb) * 1e-5 > tolerance: Mshock_old = Mshock - Mshock += step + Mshock += step Pe_old = Pe - Pt_ratio = np.power(((gamma+1)*Mshock*Mshock)/((gamma-1)*Mshock*Mshock+2),gamma/(gamma-1))*np.power((gamma+1)/(2*gamma*Mshock*Mshock-gamma+1),1/(gamma-1)) - Pi = Pt*Pt_ratio - Ai = Ath/Pt_ratio - A_ratio_target = Ae/Ai + Pt_ratio = np.power(((gamma + 1) * Mshock * Mshock) / ((gamma - 1) * Mshock * Mshock + 2),gamma / (gamma - 1)) * np.power((gamma + 1) / (2 * gamma * Mshock * Mshock - gamma + 1),1 / (gamma - 1)) + Pi = Pt * Pt_ratio + Ai = Ath / Pt_ratio + A_ratio_target = Ae / Ai Me_upper = 1 Me_lower = 0 - Me = (Me_upper+Me_lower)/2 - error = A_ratio_target-A_ratio(Me,gamma) + Me = (Me_upper + Me_lower) / 2 + error = A_ratio_target - A_ratio(Me,gamma) while np.abs(error) > tolerance: if error < 0: Me_lower = Me else: Me_upper = Me - Me = (Me_upper+Me_lower)/2 - error = A_ratio_target-A_ratio(Me,gamma) - Te = Tt/(1+((gamma-1)/2)*Me*Me) - Pe = Pi*np.power(Tt/Te,gamma/(1-gamma)) - if (Pe > Pamb) and (A_ratio(Mshock,gamma) > Ae/Ath): + Me = (Me_upper + Me_lower) / 2 + error = A_ratio_target - A_ratio(Me,gamma) + Te = Tt / (1 + ((gamma - 1) / 2) * Me * Me) + Pe = Pi * np.power(Tt / Te,gamma / (1 - gamma)) + if (Pe > Pamb) and (A_ratio(Mshock,gamma) > Ae / Ath): # the normal shock is outside the divergent and all the flow inside is supersonic - A_ratio_target = Ae/Ath + A_ratio_target = Ae / Ath Me_upper = Mshock Me_lower = 1 - Me = (Me_upper+Me_lower)/2 - error = A_ratio_target-A_ratio(Me,gamma) + Me = (Me_upper + Me_lower) / 2 + error = A_ratio_target - A_ratio(Me,gamma) while np.abs(error) > tolerance: if error > 0: Me_lower = Me else: Me_upper = Me - Me = (Me_upper+Me_lower)/2 - error = A_ratio_target-A_ratio(Me,gamma) - Te = Tt/(1+((gamma-1)/2)*Me*Me) # exit static temperature - Pe = Pt*np.power(Tt/Te,gamma/(1-gamma)) # exit static pressure + Me = (Me_upper + Me_lower) / 2 + error = A_ratio_target - A_ratio(Me,gamma) + Te = Tt / (1 + ((gamma - 1) / 2) * Me * Me) # exit static temperature + Pe = Pt * np.power(Tt / Te,gamma / (1 - gamma)) # exit static pressure break if Pe < Pamb: Pe = Pe_old @@ -234,26 +244,29 @@ def laval_nozzle(Pt,Tt,Pamb,Ath,Ae,gamma,r,tolerance=1e-5): else: # the nozzle is adapted Me = Madapted - Te = Tt/(1+((gamma-1)/2)*Me*Me) # exit static temperature + Te = Tt / (1 + ((gamma - 1) / 2) * Me * Me) # exit static temperature Pe = Pamb - Ve = Me*np.sqrt(gamma*r*Te) # exit velocity + Ve = Me * np.sqrt(gamma * R * Te) # exit velocity # maximum mass flow rate in the case of a started nozzle - Q = np.power(2/(gamma+1),(gamma+1)/(2*(gamma-1)))*np.sqrt(gamma/(Tt*r))*Pt*Ath - Thrust = Q*Ve+(Pe-Pamb)*Ae - isp = Thrust/(Q*g0) + Q = np.power(2 / (gamma + 1),(gamma + 1) / (2 * (gamma - 1))) * np.sqrt(gamma / (Tt * R)) * Pt * Ath + Thrust = Q * Ve + (Pe - Pamb) * Ae + isp = Thrust / (Q * G0) return Thrust,isp def adiabatic_flame(mech,T,P,X): - """ - Compute adiatabic flame temperature - inputs: - mech : chemical mechanism in cti format - T : initial temperature (K) - P : pressure (Pa) - X : initial chemical composition + """Compute adiabatic flame temperature + + :param mech: chemical mechanism in cti format + :type mech: [type] + :param T: initial temperature (K) + :type T: [type] + :param P: pressure (Pa) + :type P: [type] + :param X: initial chemical composition + :type X: [type] """ # work in progress @@ -292,7 +305,6 @@ def state_to_orbital_elem(r,v): orbel={'a':0,'e':0,'i':0,'capital_omega':0,'small_omega':0,'nu':0} - mu = 3.9860e14 # earth standard gravitatinal parameter (=G*M_earth) (m3.s-2) small = 1e-10 h = np.cross(r,v) # angular momentum @@ -308,19 +320,19 @@ def state_to_orbital_elem(r,v): n = np.array([1,0,0]) nm = 0 else: - n = n/nm + n = n / nm nm = 1 dotrv = np.dot(r,v) - e = ((vm*vm-mu/rm)*r-dotrv*v)/mu + e = ((vm * vm - MU / rm) * r - dotrv * v) / MU em = LA.norm(e) if em < small: e = n em = 0 else: - e = e/em + e = e / em orbel['e'] = em - orbel['a'] = (hm*hm/mu)/(1-em*em) - aux = h[2]/hm + orbel['a'] = (hm * hm / MU) / (1 - em * em) + aux = h[2] / hm if np.abs(aux) > 1: aux = np.sign(aux) # all the angles are first computed in rad orbel['i'] = np.arccos(aux) @@ -328,20 +340,20 @@ def state_to_orbital_elem(r,v): aux = n[0] if np.abs(aux) > 1: aux = np.sign(aux) orbel['capital_omega'] = np.arccos(aux) - if n[1] < 0: orbel['capital_omega'] = 2*np.pi-orbel['capital_omega'] + if n[1] < 0: orbel['capital_omega'] = 2 * np.pi - orbel['capital_omega'] if em > small: aux = np.dot(e,n) if np.abs(aux) > 1: aux = np.sign(aux) orbel['small_omega'] = np.arccos(aux) - if e[2] < 0: orbel['small_omega'] = 2*np.pi-orbel['small_omega'] - aux = np.dot(e,r)/rm + if e[2] < 0: orbel['small_omega'] = 2 * np.pi - orbel['small_omega'] + aux = np.dot(e,r) / rm if np.abs(aux) > 1: aux = np.sign(aux) orbel['nu'] = np.arccos(aux) if em > small: - if dotrv < 0: orbel['nu'] = 2*np.pi-orbel['nu'] + if dotrv < 0: orbel['nu'] = 2 * np.pi - orbel['nu'] else: ha = np.cross(e,r) - if np.dot(h,ha) < 0: orbel['nu'] = 2*np.pi-orbel['nu'] + if np.dot(h,ha) < 0: orbel['nu'] = 2 * np.pi - orbel['nu'] # convert the angles to deg orbel['i'] = np.rad2deg(orbel['i']) orbel['capital_omega'] = np.rad2deg(orbel['capital_omega']) @@ -393,3 +405,11 @@ def ox_rich_staged_combust(): def full_flow_staged_combust(): """ """ + + +if __name__ == "__main__": + + pass + + + \ No newline at end of file diff --git a/Test_Engineering_tools.py b/Test_Engineering_tools.py new file mode 100644 index 0000000000000000000000000000000000000000..dad89bc73afc51e67fe7734732ec41db66b5d980 --- /dev/null +++ b/Test_Engineering_tools.py @@ -0,0 +1,44 @@ +#!/usr/bin/env python + +from Engineering_tools import * + +import unittest + +""" + This is a prototype of test bench for Engineering_tools functions + Each function can be tested using a set of specific parameters. + The results can be evaluated and the success of test is defined by asser**** + function call. +""" +class TestEngineeringTools(unittest.TestCase): + + + + def test_initial_vel(self): + latitude = 30 + v_velocity = initial_vel(latitude) + self.assertTrue(v_velocity is not None) + + + def test_final_target_vel(self): + v_eject_altitude = 100000 + v_orbit_semi_major_axis = 637800000 + + v_f_velocity = final_target_vel(v_eject_altitude,v_orbit_semi_major_axis) + + self.assertTrue(v_f_velocity is not None) + + + def test_staging_optim(self): + v_isp1 = 300 + v_isp2 = 300 + v_eps1 = 0.8 + v_eps2 = 0.7 + v_Vc = 100 + v_Mp = 3000 + v_m1, v_m2 = staging_optim(v_isp1,v_isp2,v_eps1,v_eps2,v_Vc,v_Mp) + + self.assertTrue(v_m1 is not None and v_m2 is not None) + +if __name__ == '__main__': + unittest.main() \ No newline at end of file diff --git a/docs/Makefile b/docs/Makefile new file mode 100644 index 0000000000000000000000000000000000000000..78f01abc4abddd6f2fbedbfe6d76805a0c20a4a4 --- /dev/null +++ b/docs/Makefile @@ -0,0 +1,24 @@ +# Minimal makefile for Sphinx documentation +# + +# You can set these variables from the command line. +SPHINXOPTS = +SPHIXAPIDOC = sphinx-apidoc +SPHINXBUILD = sphinx-build +SOURCEDIR = source +BUILDDIR = build +PYTHONDIR = ../ + +sphinx-apidoc -o docs/source . + +# Put it first so that "make" without argument is like "make help". +help: + @$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) + +.PHONY: help Makefile + +# Catch-all target: route all unknown targets to Sphinx using the new +# "make mode" option. $(O) is meant as a shortcut for $(SPHINXOPTS). +%: Makefile + @$(SPHIXAPIDOC) -o $@ "$(SOURCEDIR)" "$(PYTHONDIR)" + @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) \ No newline at end of file diff --git a/docs/make.bat b/docs/make.bat new file mode 100644 index 0000000000000000000000000000000000000000..cd8aa8379434d9f57f987e8e0733885e5d514e0e --- /dev/null +++ b/docs/make.bat @@ -0,0 +1,39 @@ +@ECHO OFF + +pushd %~dp0 + +REM Command file for Sphinx documentation + +if "%SPHINXBUILD%" == "" ( + set SPHINXBUILD=sphinx-build +) +set SOURCEDIR=source +set BUILDDIR=build +set SPHIXAPIDOC=sphinx-apidoc +set PYTHONDIR=../ + +if "%1" == "" goto help + +%SPHINXBUILD% >NUL 2>NUL +if errorlevel 9009 ( + echo. + echo.The 'sphinx-build' command was not found. Make sure you have Sphinx + echo.installed, then set the SPHINXBUILD environment variable to point + echo.to the full path of the 'sphinx-build' executable. Alternatively you + echo.may add the Sphinx directory to PATH. + echo. + echo.If you don't have Sphinx installed, grab it from + echo.http://sphinx-doc.org/ + exit /b 1 +) + +%SPHIXAPIDOC% -o %SOURCEDIR% %PYTHONDIR% + +%SPHINXBUILD% -M %1 %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% +goto end + +:help +%SPHINXBUILD% -M help %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% + +:end +popd diff --git a/docs/source/conf.py b/docs/source/conf.py new file mode 100644 index 0000000000000000000000000000000000000000..10694a901e8da13372a5494bf090c5e1681ea401 --- /dev/null +++ b/docs/source/conf.py @@ -0,0 +1,59 @@ +# Configuration file for the Sphinx documentation builder. +# +# This file only contains a selection of the most common options. For a full +# list see the documentation: +# http://www.sphinx-doc.org/en/master/config + +# -- Path setup -------------------------------------------------------------- + +# If extensions (or modules to document with autodoc) are in another directory, +# add these directories to sys.path here. If the directory is relative to the +# documentation root, use os.path.abspath to make it absolute, like shown here. +# +import os +import sys +sys.path.insert(0, os.path.abspath('../../')) + + +# -- Project information ----------------------------------------------------- + +project = 'meca_vol' +copyright = '2019, Adastra team' +author = 'Adastra team' + +# The full version, including alpha/beta/rc tags +release = '0.0.0' + + +# -- General configuration --------------------------------------------------- + +# Add any Sphinx extension module names here, as strings. They can be +# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom +# ones. +extensions = [ + 'sphinx.ext.autodoc', + #'sphinx.ext.imgmath' + 'sphinx.ext.mathjax' +] + +# Add any paths that contain templates here, relative to this directory. +templates_path = ['_templates'] + +# List of patterns, relative to source directory, that match files and +# directories to ignore when looking for source files. +# This pattern also affects html_static_path and html_extra_path. +exclude_patterns = [] + +master_doc = 'index' + +# -- Options for HTML output ------------------------------------------------- + +# The theme to use for HTML and HTML Help pages. See the documentation for +# a list of builtin themes. +# +html_theme = 'sphinx_rtd_theme' + +# Add any paths that contain custom static files (such as style sheets) here, +# relative to this directory. They are copied after the builtin static files, +# so a file named "default.css" will overwrite the builtin "default.css". +html_static_path = ['_static'] diff --git a/docs/source/index.rst b/docs/source/index.rst new file mode 100644 index 0000000000000000000000000000000000000000..c5a52dc710ff44d2f7ac5c55ce0dd76ae7410f93 --- /dev/null +++ b/docs/source/index.rst @@ -0,0 +1,21 @@ +.. meca_vol documentation master file, created by + sphinx-quickstart on Sun Jul 7 20:22:14 2019. + You can adapt this file completely to your liking, but it should at least + contain the root `toctree` directive. + +Welcome to meca_vol's documentation! +==================================== + +.. toctree:: + :maxdepth: 2 + :caption: Contents: +name: +Test_Engineering_tools + + +Indices and tables +================== + +* :ref:`genindex` +* :ref:`modindex` +* :ref:`search` diff --git a/requirements.txt b/requirements.txt index dc9d9518f8dfcbc864fe1f2bd853af0ce1fd6e14..e9f817dd1c86d17383f25bd1dac2080eeae13672 100644 --- a/requirements.txt +++ b/requirements.txt @@ -7,3 +7,4 @@ pyparsing==2.4.0 python-dateutil==2.8.0 pytz==2019.1 six==1.12.0 +Sphinx==2.0.1 \ No newline at end of file diff --git a/test_laval_nozzle.py b/test_laval_nozzle.py index a6d925a5801d5644018120f6d7f0222e59b1f4bb..432782d542893cbedfae590b5d805eae2e5d70ce 100755 --- a/test_laval_nozzle.py +++ b/test_laval_nozzle.py @@ -23,14 +23,16 @@ Pamb_tab = np.linspace(1e5,100,100) # array of the ambient pressure (Pa) for Pamb in Pamb_tab: data.loc[len(data)] = (Pamb,)+laval_nozzle(Pt,Tt,Pamb,Ath,Ae,gamma,r) -plt.figure('Laval Nozzle Thrust') -plt.plot(data['Pamb']*1e-5,data['Thrust']*1e-3) -plt.title('Thrust Vs Ambient Pressure') -plt.xlabel('Ambient Pressure (bar)') -plt.ylabel('Thrust (kN)') -plt.figure('Laval Nozzle Isp') -plt.plot(data['Pamb']*1e-5,data['isp']) -plt.title('Isp Vs Ambient Pressure') -plt.xlabel('Ambient Pressure (bar)') -plt.ylabel('Isp (s)') +if __name__ == '__main__': + plt.figure('Laval Nozzle Thrust') + plt.plot(data['Pamb']*1e-5,data['Thrust']*1e-3) + plt.title('Thrust Vs Ambient Pressure') + plt.xlabel('Ambient Pressure (bar)') + plt.ylabel('Thrust (kN)') + + plt.figure('Laval Nozzle Isp') + plt.plot(data['Pamb']*1e-5,data['isp']) + plt.title('Isp Vs Ambient Pressure') + plt.xlabel('Ambient Pressure (bar)') + plt.ylabel('Isp (s)') diff --git a/test_standard_atmosphere.py b/test_standard_atmosphere.py index 2ca63bd160bc91e83378a0c6558b4d68f0b6cfad..3faae58f6a898a11ba1a7e771e3bd13d03aaa32d 100755 --- a/test_standard_atmosphere.py +++ b/test_standard_atmosphere.py @@ -15,23 +15,27 @@ data = pd.DataFrame(columns=['altitude','T','P','rho']) for h in altitude: data.loc[len(data)] = (h,)+standard_atmosphere(h) - -plt.figure('ISA Temperature') -plt.plot(data['T'],data['altitude']*1e-3) -plt.title('Temperature Vs Altitude') -plt.xlabel('Temperature (K)') -plt.ylabel('Altitude (km)') - -plt.figure('ISA Pressure') -plt.plot(data['P'],data['altitude']*1e-3) -plt.xscale('log') -plt.title('Pressure Vs Altitude') -plt.xlabel('Pressure (Pa)') -plt.ylabel('Altitude (km)') - -plt.figure('ISA Density') -plt.plot(data['rho'],data['altitude']*1e-3) -plt.xscale('log') -plt.title('Density Vs Altitude') -plt.xlabel('Density (Kg.m-3)') -plt.ylabel('Altitude (km)') \ No newline at end of file + + +if __name__ == '__main__': + plt.figure('ISA Temperature') + plt.plot(data['T'],data['altitude']*1e-3) + plt.title('Temperature Vs Altitude') + plt.xlabel('Temperature (K)') + plt.ylabel('Altitude (km)') + + plt.figure('ISA Pressure') + plt.plot(data['P'],data['altitude']*1e-3) + plt.xscale('log') + plt.title('Pressure Vs Altitude') + plt.xlabel('Pressure (Pa)') + plt.ylabel('Altitude (km)') + + plt.figure('ISA Density') + plt.plot(data['rho'],data['altitude']*1e-3) + plt.xscale('log') + plt.title('Density Vs Altitude') + plt.xlabel('Density (Kg.m-3)') + plt.ylabel('Altitude (km)') + + plt.show() \ No newline at end of file