##Design and Analysis of Launch Vehicle and Trajectory

##Importing required directories for the analysis and design

import math as math
import numpy as np
import matplotlib.pyplot as plt

##Computing and analyzing total energy loss due to gravity and drag ascent mission by varying the burn rate of the propellent

m_0     = 80000;                    #kg; Total Mass of the Rocket
m_p     = 60000;                    #kg; Mass of the propellent
m_b     = m_0 - m_p;                #kg; Burnout mass of the Rocket
g_0     = 9.81;                     #m/s^2; Acceleration due to gravity at sea level
I_sp    = 300;                      #sec; Specific Impulse of the rocket engine
ppf     = m_p/m_0;                  #Propellent Fraction
beta    = list(range(300,1601));    #kg/s; Burn rate of the propellent
rho     = 0.225;                    #kg/m^3; Averaged density at the terminal altitude
S_ref   = math.pi;                  #m^2; Net surface area of the bluff body
c_d_0   = 1;                        #Drag coefficient of the body

##Calculating ideal burnout velocity of the rocket (m/s)
v_ideal     = I_sp*g_0*math.log(m_0/m_b);

##Ideal energy generated by the rocket at terminal performance (J)
e_ideal = 0.5*(v_ideal)**2;

##Maximum possible terminal burnout velocity under the effect of gravity (m/s)
v_g     = v_ideal - g_0*(np.divide(ppf*m_0,beta));

##Calculating total drag produced at the terminal end of the trajectory (N)
D       = 0.5*rho*((v_g)**2)*S_ref*c_d_0;

##Deceleration on the rocket body at the terminal altitude (m/s^2)
a_D     = np.divide(D,m_b);

##Average decelaration on the rocket body (m/s^2)
a_avg   = np.divide(a_D,2);

##Maximum possible terminal burnout velocity under the effect of drag and gravity (m/s)
v_D     = v_g - a_avg*(np.divide(ppf*m_0,beta));

##Terminal altitude of the rocket achieved under the gravity (m)
h_g     = ((np.divide(m_0*I_sp*g_0,beta))*((1-ppf)*math.log(1-ppf)+ppf))-(0.5*g_0*(np.divide(ppf*m_0,beta))**2);

##Terminal altitude of the rocket achieved under the drag and gravity (m)
h_D     = h_g - (0.5*a_avg*(np.divide(ppf*m_0,beta))**2);

##Combined loss of energy due to drag and gravity (J)
delta_e = 0.5*(v_ideal)**2 - 0.5*(v_D)**2 - g_0*h_D;

##Percentage change in total energy due to drag and gravity
e_per   = (np.divide(delta_e,e_ideal))*100;

plt.plot(beta,e_per,':', color='b',linewidth=3)
plt.rcParams["figure.figsize"] = (10,6)
plt.title('Total Energy Loss Analysis in Ascent mission due to Gravity and Drag')
plt.grid(color='tab:gray', linestyle=':', linewidth=1)
plt.xlabel('Burn Rate of the Propellent (kg/s)')
plt.ylabel('Percentage change in Total Energy due to Drag and Gravity')


