diff --git a/Rocket_Simulator.py b/Rocket_Simulator.py index 5ec0fef..f626826 100644 --- a/Rocket_Simulator.py +++ b/Rocket_Simulator.py @@ -4,7 +4,7 @@ #EVERYTHING WORKS EXCEPT FOR THE MAX ACCELERATION. MAX ACC IS ALMOST 10M/S^2 LESS THAN OPEN ROCKET from this import d -from turtle import color +# from turtle import color import matplotlib.pyplot as plt import csv import time as timee @@ -40,22 +40,22 @@ G = 6.67408 * pow(10,-11) # Gravitational constant ME = 5.972 * pow(10,24) # Earth mass in Kg Re = 6371 * pow(10,3) # Radius of earth in m -cg = 0.528 #in cm -cp = 0.728 #in cm +cg = 0.528 #in cm (distance from rocket's nose to center of pressure) +cp = 0.728 #in cm (distance from rocket's nose to center of gravity) mr = 0.618 #Mass without casing or propellant in kg -me = 0.043 -mp = 0.060 -im = 49.61 #in N-s -th = 14.5 #Avg thrust in N(use this) -mth = 25.3 #Max Thrust in N -tb = 3.42 #in s (im/th) -rho = 1.22 #in Kg/m^3, Drops by 10% for each 1000m of altitude for the first 14km or we can use the formula = 1.22 * 0.9^(Altitude/1000) -cd = 0.75 #Avg Cd, but we can calculate it based on our rocket dimensions. -dr = 0.045 #in m +me = 0.043 #Engine case mass +mp = 0.012 #Propellent mass defaullt=0.060 +im = 8 #in N-s default=49.61 +th = 6.7 #Avg thrust in N(use this) default=14.5 +mth = 12 #Max Thrust in N default=25.3 +tb = 3.42 #in s (im/th) +rho = 1.22 #(density of air) in Kg/m^3, Drops by 10% for each 1000m of altitude for the first 14km or we can use the formula = 1.22 * 0.9^(Altitude/1000) +cd = 0.75 #Avg Cd, but we can calculate it based on our rocket dimensions (drag coeff) +dr = 0.045 #in m (rocket diameter) a = 3.14 * pow(dr/2,2) #Area of the rocket dt = 0.01 #in s dm = (mp/tb)*dt#mass flow rate -g = 9.8 #m/s +g = 9.8 #m/s2 tmass = mr + me + mp #Total Mass in kg alt = [] vel = [] @@ -72,7 +72,7 @@ prevHeight = 0 t2apogee = 0 -theta = 20 #in degrees from x axis. +theta = 20 #in degrees from x axis. Given 20 thetaV = 0 # Velocity directions # Currently, we are considering that the angle remains the same throughout the entire flight. @@ -121,11 +121,15 @@ def getThrust(t): def getDrag(v): d = 0.5 * rho * cd * a * pow(v,2) - return d + d_x = d * math.cos(math.radians(theta)) + d_y = d * math.sin(math.radians(theta)) + return [d, d_x, d_y] def getParachuteDrag(v): d = (cd * rho * 3.14 * pow(diameter,2) * pow(v,2))/8 #(diameter /2) ^ 2 ==> (4) and then the entire formula is 1/2 * .... ==> 1/8 - return d + d_x = d * math.cos(math.radians(theta)) + d_y = d * math.sin(math.radians(theta)) + return [d, d_x, d_y] def getAcc(tmass,th,v,velY,altY,t,h,parachuteFlag,ascentFlag): # In the beginning, acceleration is -9.8(Pointing downwards), and then it gradually increases. Until, it reaches a +ve value, @@ -134,41 +138,54 @@ def getAcc(tmass,th,v,velY,altY,t,h,parachuteFlag,ascentFlag): if parachuteFlag == 0 and len(velY) > 0 and velY[-1]<0 and len(altY) > 0 and altY[-1]>0: parachuteFlag = 1 if parachuteFlag == 1: # After parachute is deployed - d = getParachuteDrag(v) + getDrag(v) + d = getParachuteDrag(v)[0] # + getDrag(v) + dx = getParachuteDrag(v)[1] + dy = getParachuteDrag(v)[2] thh = getThrust(t) - acceleration = (d-(tmass*g))/tmass # Sin theta not considered because, once the parachute is deployed, the rocket is going to fall straight down. So the angle is not going to be considered. + thh_x = thh*math.cos(math.radians(theta)) + thh_y = thh*math.sin(math.radians(theta)) + acceleration_x = 0 # wind velocity requires adding here + acceleration_y = ((d - (tmass * g)) / tmass) # Sin theta not considered because, once the parachute is deployed, the rocket is going to fall straight down. So the angle is not going to be considered. + acceleration = acceleration_y # for now, since acc along x is 0 without wind print("Status : Parachute Deployed(Descent)","Time : ",t,"Thrust : ",thh,"Drag : ", d,"Weight : ",tmass*g,"Force : ",(d - (tmass*g) ),"Mass : ",tmass,"Acceleration",acceleration) else: # Before Parachute is deployed # If parachute hasn't deployed, then either the rocket is still on the launch pad, or it is in the ascent phase. thh = getThrust(t) - d = getDrag(v) - acceleration = (thh - ((tmass*g)*math.sin(math.radians(theta))) - d)/tmass # math.sin considers radians, hence theta is first converted to radians, then passed into math.sin + thh_x = thh*math.cos(math.radians(theta)) + thh_y = thh*math.sin(math.radians(theta)) + d = getDrag(v)[0] + dx = getDrag(v)[1] + dy = getDrag(v)[2] + + acceleration_x = ((thh_x - dx)/tmass) # -+ Wind can also be added + acceleration_y = ((thh_y - dy - (tmass * g))/tmass) # ERROR HERE : ACC_Y coming as negative of Acc for some reason + acceleration = ((thh - ((tmass * g) * math.sin(math.radians(theta))) - d) / tmass)# math.sin considers radians, hence theta is first converted to radians, then passed into math.sin # The reason, why we are considering mg.sin(theta) here is because, we are considering that the rocket is going to be launched at a launch angle theta. So, when calculating the net force, at that angle, mg.sin(theta) is used. print("Status : Ascent, Theta : ",theta,"Time : ",t,"Thrust : ",thh,"Drag : ", d,"Weight : ",((tmass*g)*math.sin(math.radians(theta))),"Force : ",(thh - ((tmass*g)*math.sin(math.radians(theta))) - d),"Mass : ",tmass,"Acceleration",acceleration) if ascentFlag == 0 and acceleration >= 0: # If the rocket hasn't begun its ascent phase and the acceleration is +ve, meaning ascent phase is beginning ascentFlag = 1 # Ascent has begun - elif ascentFlag == 1 and acceleration<0: #If the rocket has finished its ascent phase and the acceleration becomes -ve, meaning pointing downwards, means that the ascent phase has ended and the descent phase has begun + elif (ascentFlag == 1) and (acceleration < 0): #If the rocket has finished its ascent phase and the acceleration becomes -ve, meaning pointing downwards, means that the ascent phase has ended and the descent phase has begun ascentFlag = 2 # Ascent has ended and Descent has begun - return [acceleration,ascentFlag] + return [acceleration, ascentFlag, acceleration_x, acceleration_y] -def isCondition(altY,t,parachuteFlag): # Condition to check whether the rocket has landed or not - if len(altY) > 0 and altY[-1]<=0 and t > 0 and parachuteFlag == 1: +def isCondition(altY, t, parachuteFlag): # Condition to check whether the rocket has landed or not + if (len(altY) > 0) and (altY[-1] <= 0) and (t > 0) and (parachuteFlag == 1): return 0 return 1 -while isCondition(altY,t,parachuteFlag1): +while isCondition(altY, t, parachuteFlag1): print("*(***************************",len(velY),len(altY)) if len(velY) != 0 and len(altY) != 0: print("VelY : ", velY[-1], " AltY : ", altY[-1],"Parachute Flag : ", parachuteFlag1) - if parachuteFlag1 == 0 and velY[-1]<0 and altY[-1]>0: # Apogee condition, Checking if the last vertical velocity < 0 , and last vertical altitude is > 0 + if parachuteFlag1 == 0 and velY[-1] < 0 and altY[-1] > 0: # Apogee condition, Checking if the last vertical velocity < 0 , and last vertical altitude is > 0 parachuteFlag1 = 1 #Parachute is deployed t2apogee = t # Time to apogee set at apogee # Gravity modelling(Change in acceleration due to gravity) # g = (G*ME)/(h + Re) try: - theta = math.degrees(math.atan(velY[-1]/velX[-1])) # Theta is the angle at which the rocket is going to be launched. + theta = math.degrees(math.atan(velY[-1] / velX[-1])) # Theta is the angle at which the rocket is going to be launched. print("theta:",theta) except: print("Exception in theta calculations") @@ -177,11 +194,12 @@ def isCondition(altY,t,parachuteFlag): # Condition to check whether the rocket h # Accel, Vel, Height calculations inp = getAcc(tmass,th,v,velY,altY,t,h,parachuteFlag1,ascentFlag) acc = float(inp[0]) + acc_x = float(inp[2]) + acc_y = float(inp[3]) ascentFlag = int(inp[1]) if ascentFlag == 0: v = 0 h = 0 - velY.append(v) accelY.append(acc) else: @@ -190,8 +208,6 @@ def isCondition(altY,t,parachuteFlag): # Condition to check whether the rocket h # Breaking acc into its x and y vectors theta = 90 # Considering that when the parachute is deployed, the rocket will fall straight down. So theta is 90 degrees from x axis and 0 from y # acc_x = acc * math.cos(math.radians(theta)) - acc_x = 0 - acc_y = acc * math.sin(math.radians(theta)) print("COS 90 : " , math.cos(math.radians(90)),"SIN 90 : ", math.sin(math.radians(90)),"ACC_X : ",acc_x,"ACC_Y : ",acc_y) @@ -210,14 +226,14 @@ def isCondition(altY,t,parachuteFlag): # Condition to check whether the rocket h # vy = v_y + acc_y * dt v += acc_y * dt - v_y = v + v_y = v print("++++++++++++++++++v_y : ",vy,"acc_y : ",acc_y ,"dt : ",dt,"acc_y * dt : ",acc_y * dt,"v_y + acc_y * dt : ",v_y + acc_y * dt) # v = math.sqrt(pow(vx,2) + pow(vy,2)) # print(acc,dt,pow(dt,2),acc*pow(dt,2)) - h += v_y*dt + (0.5*acc_y*pow(dt,2)) # Because altitude is basically distance in the y axis + h += v_y * dt + (0.5 * acc_y * pow(dt,2)) # Because altitude is basically distance in the y axis hx += v_x * dt altX.append(hx) @@ -238,8 +254,8 @@ def isCondition(altY,t,parachuteFlag): # Condition to check whether the rocket h else: # Parachute hasn't been deployed yet # Breaking acc into its x and y vectors - acc_x = acc * math.cos(math.radians(theta)) - acc_y = acc * math.sin(math.radians(theta)) + # acc_x = acc * math.cos(math.radians(theta)) + # acc_y = acc * math.sin(math.radians(theta)) # Breaking previous velocity and altitude into its x and y vectors v_x = v * math.cos(math.radians(theta)) # Previous velocity x components @@ -249,8 +265,9 @@ def isCondition(altY,t,parachuteFlag): # Condition to check whether the rocket h vx = v_x + acc_x * dt vy = v_y + acc_y * dt v = math.sqrt(pow(vx,2) + pow(vy,2)) - h += v_y*dt + (0.5*acc_y*pow(dt,2)) # Because altitude is basically distance in the y axis - hx += v_x * dt + h += v_y * dt + (0.5 * acc_y * pow(dt, 2)) # Because altitude is basically distance in the y axis + hx += v_x * dt + (0.5 * acc_x * pow(dt, 2)) + altX.append(hx) altY.append(h) diff --git a/tempCodeRunnerFile.py b/tempCodeRunnerFile.py new file mode 100644 index 0000000..f1cdc78 --- /dev/null +++ b/tempCodeRunnerFile.py @@ -0,0 +1 @@ + , \ No newline at end of file