DensityVsAltitude-NRLMSISE00-JB2008

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

import json
import matplotlib.pyplot as plt
import numpy as np
import requests 

"""
Plot variation of density with altitudes as calculated by NRLMSISE00
and JB2008 accessed via the Amentum API
Visit https://atmosphere.amentum.space for API docs
"""

hostname = "https://atmosphere.amentum.space/api/"

# init the plot
fig = plt.figure(figsize=(4.4, 5.6))
ax = fig.add_subplot(111)
ax.set_ylabel("Altitude [km]")
ax.set_xlabel(r"$\rho _{air} \quad [kgm^{-3}]$")

for endpoint in ["nrlmsise00", "jb2008"]:

  # consider altitudes 200-1000 km 
  altitudes = np.arange(200, 900, 10)

  # WARNING: nominal parameters, densities will vary significantly
  # with space weather and position
  payload = {
    'year' : 2018,
    'month' : 12,
    'day' : 31,# NYE party
    'geodetic_latitude' : 33.9207, # above SpaceX HQ
    'geodetic_longitude' : 360.0-118.3278, # long on interval [0,360] deg
    'utc' : 7, # Happy new year! 
  }

  # Only nrlmsise00 requires specification of radio flux, JB2008 fetches indices automatically. 
  if endpoint == "nrlmsise00":
    altitudes = np.arange(0, 900, 10) # valid to ground
    payload.update({
      'f107' : 69, # taken from ftp://ftp.swpc.noaa.gov/pub/warehouse/2018/
      'f107a' : 69, # nominal value only
      'ap' : 7, # from ftp://ftp.gfz-potsdam.de/pub/home/obs/kp-ap/tab/
      })

  densities = []

  for alt in altitudes: 
    # 
    payload.update({
      'altitude' : alt   
    })
    # hit the API
    try:
      response = requests.get(hostname+endpoint, params=payload)
      response.raise_for_status()
    except requests.exceptions.HTTPError as e:  
      print(e)
      raise KeyboardInterrupt

    # extract mass density value from the JSON (response can include composition and temps as well)
    densities.append(response.json()['total_mass_density']['value'])

  # plot altitudes vs densities
  ax.semilogx(
      densities,
      altitudes,
      marker="None",
      linestyle="-",
      label=endpoint.upper(),
  )

plt.ylim(bottom=0)
ax.grid()
ax.legend()
plt.tight_layout()
plt.savefig("air_density_vs_altitude.png")