PETra/postprocessing.py at master · tr1010/PETra · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
postprocesses the raw outputs from PETra's main function

Created on Mon Dec 12 15:44:46 2016

@author: tr1010 (Thomas Rees)
"""
import numpy as np
import matplotlib.pyplot as plt
import atmospheres as atmo

def postprocess(t, sol, earth, mass, areas, normals, centroids, I, x0, scLengthScale):
    r2d = 180./np.pi
    ndt = np.size(t)
    # Unpack solution array:

    #Position
    Altitude = (sol[:,0] - earth[1])*1e-3
    Latitude =  sol[:,1]*r2d # Actually Declination
    Longitude = sol[:,2]*r2d

    # Geocentric Velocity components
    u = sol[:,3]
    v = sol[:,4]
    w = sol[:,5]

    # Velocity components to Spherical Coords
    Speed = np.zeros(ndt)
    for i in range(0,ndt):
        Speed[i] = np.linalg.norm(np.array([u[i],v[i],w[i]]))

    FlightPathAngle = np.arcsin(np.divide(w,Speed))*r2d
    Heading = np.arccos(np.divide(u,Speed*np.cos(FlightPathAngle/r2d)))*r2d

    # Quaternion to Euler Angles
    Pitch = np.rad2deg(-np.arcsin(2*(sol[:,7]*sol[:,9] - sol[:,8]*sol[:,6])))
    Yaw = np.rad2deg(np.arctan2(2*(sol[:,8]*sol[:,7]+sol[:,6]*sol[:,9]),sol[:,6]**2+sol[:,7]**2 - sol[:,8]**2-sol[:,9]**2))
    Roll = np.rad2deg(np.arctan2(2*(sol[:,8]*sol[:,9]+sol[:,6]*sol[:,7]),sol[:,6]**2-sol[:,7]**2 - sol[:,8]**2+sol[:,9]**2))

    # Angular Velocities (around body axes)
    wx = sol[:,10]
    wy = sol[:,11]
    wz = sol[:,12]

    # Now calculate interesting outputs --
    # Mach number, Knudsen number, Energy, Force, Dynamic Pressure

    # Get atmospheric quantities at each point (this should be in main function
    # I promise I'll fix it soon)
    rho = np.zeros(ndt)
    P = np.zeros(ndt)
    T = np.zeros(ndt)
    R = np.zeros(ndt)
    mfp = np.zeros(ndt)
    eta = np.zeros(ndt)
    MolW = np.zeros(ndt)
    SoS = np.zeros(ndt)

    for i in range(0,ndt):
        if sol[i,0] <= earth[1]:
            break
        else:
            rho[i], P[i], T[i], R[i], mfp[i], eta[i], MolW[i], SoS[i] = atmo.nrlmsise00(172,0,29000,sol[i,0]-earth[1],sol[i,1],sol[i,2],16,150,150,4)
            #rho[i], P[i], T[i], mfp[i], eta[i], MolW[i], SoS[i]  = atmo.US62_76(sol[i,0],earth[1])

    Mach = np.divide(Speed,SoS)
    Kn = mfp/scLengthScale
    Force = np.diff(Speed)/np.diff(t)/-9.81
    DynamicPressure = 0.5*np.multiply(rho,np.square(Speed))
    Energy = 0.5*mass*np.square(Speed) + mass*9.81*Altitude*1e3

    # Put all results in dictionary
    pp_output = dict([('FlightPathAngle' , FlightPathAngle),
                      ('Heading' , Heading),
                      ('Speed' , Speed),
                      ('Latitude' , Latitude),
                      ('Longitude' , Longitude),
                      ('Altitude' , Altitude),
                      ('Pitch' , Pitch),
                      ('Yaw' , Yaw),
                      ('Roll' , Roll),
                      ('omega_x' , wx),
                      ('omega_y' , wy),
                      ('omega_z' , wz),
                      ('Mach' , Mach),
                      ('Knudsen' , Kn),
                      ('Acceleration' , Force),
                      ('DynamicPressure' , DynamicPressure),
                      ('Energy' , Energy)])

    # Plot Results in Python for user sanity check
    plotfn(0,t,pp_output['FlightPathAngle'], 'time', 'Flight Path Angle', 's', 'deg')
    plotfn(1,t,pp_output['Heading'], 'time', 'Heading', 's', 'deg')
    plotfn(2,t,pp_output['Speed'], 'time', 'Speed', 's', 'm/s')
    plotfn(3,pp_output['Longitude'],pp_output['Latitude'], 'Longitude', 'Latitude', 'deg', 'deg')
    plotfn(4,t,pp_output['Altitude'], 'time', 'Altitude', 's', 'km')
    plotfn(5,t,pp_output['Mach'], 'time', 'Mach Number', 's', '-')
    plotfn(6,t,pp_output['Knudsen'], 'time', 'Knudsen number', 's', '-')
    plotfn(7,t[0:ndt-1],pp_output['Acceleration'], 'time', 'Force', 's', 'g')
    plotfn(8,t,pp_output['DynamicPressure']*1e-3, 'time', 'Dynamic Pressure', 's', 'kPa')
    plotfn(9,t,pp_output['Energy']*1e-9, 'time', 'Energy', 's', 'GJ')

    plt.figure(10)
    fig,ax = plt.subplots()
    plt.title('Euler angles')
    ax.plot(t,pp_output['Pitch'],label='pitch')
    ax.plot(t,pp_output['Yaw'], label='yaw')
    ax.plot(t,pp_output['Roll'],label = 'roll')
    legend = ax.legend()
    plt.xlabel('time [s]')
    plt.ylabel('angle [deg]')
    plt.show()

    return pp_output


def plotfn(figno, x, y, xname, yname, xunits, yunits):

    fig = plt.figure(figno)
    plt.plot(x,y)
    fig.suptitle(xname + ' vs ' + yname)
    plt.xlabel(xname + ' [' + xunits +']')
    plt.ylabel(yname + ' [' + yunits +']')
    plt.savefig('Figures/'+ xname + '-' + yname+'.png', bbox_inches='tight')
    plt.show

    return