In [2]:
import pykep as pk
import numpy as np
import json
import pickle as pkl
import requests
import pygmo as pg
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import pyplot as plt
%matplotlib inline
In [3]:
import cascade as csc
from copy import deepcopy
from tqdm.notebook import tqdm
import heyoka as hy
We obtain atmospheric data from the nrlmsise00 model¶
An account is needed at https://amentum.com.au/atmosphere and an API-Key allowing a limited amounts of queries
In [4]:
# This will contain the density in kg/m3
density = []
# This is the grid on the altitude
altitude = np.linspace(78, 1000, 100)
... uncomment and run this cell to retreive the data again, else use the pickled file
In [5]:
#url = "https://atmosphere.amentum.io/nrlmsise00"
#headers = {"API-Key" : "Pmd5s5e36Ny0zpPGhxvDWRtgrmB2rR6K"}
#for alt in altitude:
# params = {
# 'year' : 2021,
# 'month' : 2,
# 'day' : 1,
# 'geodetic_latitude' : 0,
# 'geodetic_longitude' : 42,
# 'altitude' : alt, # km
# 'utc' : 2, # hours
# }
# # handle exceptions
# response = requests.get(url, params=params, headers=headers)
# json_payload = response.json()
# density.append(json_payload['total_mass_density']['value'])
#with open('data/atmosphere.pk', 'wb') as file:
# pkl.dump((altitude, density), file)
In [5]:
with open("data/atmosphere.pk", "rb") as file:
altitude, density = pkl.load(file)
We try and fit the data with some predefined shape¶
It is important that the limit for h that goes to infinity of the density goes to zero. The chosen shape must guarantee this property. We here use:
$$ \rho(h) = \sum_i\alpha_i\exp(-\beta_i(h-\gamma_i)) $$In [6]:
class my_udp:
def __init__(self, X, Y, n=4):
self.X = X
self.Y = Y
self.n = n
def fitness(self,x):
p1 = x[:self.n]
p2 = x[self.n:2*self.n]
p3 = x[2*self.n:]
Y_pred = self._fitting_function(self.X,p1,p2, p3)
return [np.mean(np.abs(np.log10(self.Y)-np.log10(Y_pred)))]
def _fitting_function(self, x, p1,p2, p3):
retval = 0.
for alpha,beta, gamma in zip(p1,p2, p3):
retval += alpha*np.exp(-(x-gamma)*beta)
return retval
def get_bounds(self):
lb = [0]*self.n*3
ub = [1.]*self.n*3
ub[self.n*2:] = [1000]*self.n
return (lb, ub)
In [7]:
udp = my_udp(altitude, density, n=4)
prob = pg.problem(udp)
In [12]:
uda = pg.sade(50, memory=True)
algo = pg.algorithm(uda)
In [13]:
pop =pg.population(prob, 20)
/tmp/ipykernel_3576922/3502744231.py:15: RuntimeWarning: overflow encountered in exp retval += alpha*np.exp(-(x-gamma)*beta)
In [14]:
for i in range(500):
pop = algo.evolve(pop)
print(pop.champion_f[0])
/tmp/ipykernel_3576922/3502744231.py:15: RuntimeWarning: overflow encountered in exp retval += alpha*np.exp(-(x-gamma)*beta) /tmp/ipykernel_3576922/3502744231.py:11: RuntimeWarning: divide by zero encountered in log10 return [np.mean(np.abs(np.log10(self.Y)-np.log10(Y_pred)))]
5.000578643324542 3.026808000514021 2.4888110377559025 2.3741149857341592 2.312471379267102 2.2576780689300335 2.0151634575017283 2.0084685558333772 1.5234155262675437 1.3571149607583286 1.0855638472547917 0.8467652617569728 0.7011640993233682 0.5934220922537995 0.5363040572810933 0.5078961266510804 0.4920599286034672 0.48799479536958273 0.479819505196312 0.4716551608588495 0.46825567427994175 0.4484020342723868 0.445002529220869 0.4377128856788668 0.427695253821798 0.4229629427606548 0.4200896953749397 0.4167008624915533 0.41361431353847544 0.4091325606249492 0.40309691405076675 0.3980713123493298 0.3912029354224048 0.3851995593468093 0.3788508992534813 0.37098361940892416 0.36163285492529007 0.3561494602323464 0.35368312363426146 0.33929542194608464 0.33364284014919504 0.3273158916656956 0.3135478041939112 0.3027096069851567 0.2835032997383042 0.27669080196643236 0.2713205888237881 0.2628590563356096 0.2598657268055497 0.2550407640760979 0.2540268590371742 0.24907796851812145 0.23959344882534878 0.23542166501748768 0.23288154786562096 0.22394633525698684 0.22007610547669385 0.21060876115310465 0.20443083545469265 0.19961557424758344 0.1940756870715757 0.18651570742293846 0.1781016678830613 0.17369337142826374 0.16698619406646173 0.15910149674434648 0.15403364035119416 0.14231691494683466 0.1395002407554455 0.13170647530527865 0.12324990458467226 0.1169945714917992 0.10960988588670306 0.10644681964187655 0.10273929343267044 0.10273929343267044 0.10273921077719367 0.09569923646596111 0.09156711673733456 0.09054846927156188 0.09039879820452842 0.09039879820452842 0.0878004038254625 0.08752095929289397 0.08241040954391866 0.08218034438515787 0.08088739464892686 0.07663325219674079 0.07470603455588634 0.07279300242426047 0.07220101532924912 0.0713367984490804 0.07095171129746897 0.0703295143798569 0.06906259025990472 0.06597661269811318 0.06483987753487981 0.06470492758421352 0.06389135099128224 0.06327347705869955 0.06301510443269336 0.0622108180640035 0.06156054380647122 0.059915027795085854 0.05968036437365327 0.05924887238946277 0.05907321990401259 0.05881136787309061 0.05875917339826256 0.058732871143027755 0.05872052747452336 0.058709197730248734 0.058709097905236235 0.058702553036291086 0.05870104977127192 0.058691466841453724 0.05869059647827932 0.05868644283707711 0.05868415002462084 0.05868300613620628 0.05867944608823571 0.05867891662290886 0.058677703517599766 0.058671307114957244 0.05866864369385569 0.058668374882563164 0.058668374778486125 0.05866701177078011 0.0586665561609856 0.058663610768035925 0.05866212871312979 0.05866212867385527 0.05866207750613922 0.05866142486875221 0.05866142486875221 0.05866053690354258 0.0586589386661164 0.05865696730789123 0.05865696730649865 0.05865535101479122 0.05865393321911934 0.05865393321911934 0.05865350888918775 0.058652849077985285 0.058652849077985285 0.058652800136009375 0.058651379721646905 0.05865094543613008 0.05864985904593143 0.058649859038742536 0.058649583459958084 0.05864875771847456 0.05864875771843256 0.05864838156887919 0.05864666878944143 0.05864666878944143 0.05864666878944143 0.058644611306578234 0.05864461130440735 0.0586444820761525 0.058643661635015464 0.058643453688196706 0.05864338348547432 0.05864280823726556 0.05864218915377721 0.05864128413021084 0.05862422490000797 0.0574803998845354 0.05724386433482603 0.05680018295286984 0.05665319581473068 0.055941904113954166 0.05589392498939354 0.055179323872021 0.05497894924505599 0.05436619807613646 0.05264784067341612 0.0524748080396169 0.052360250644621055 0.05159743519239462 0.05058877763562947 0.05058800831702628 0.05035916889183357 0.050231491941899124 0.049995200893906494 0.049517948744033286 0.047351070521415656 0.04713214475522915 0.04687439174775425 0.046404757523906444 0.04584436900968516 0.04576869111071847 0.04447393771780149 0.04254979269326746 0.04175180061225008 0.04066118368592922 0.039569957391937985 0.03771002904231824 0.035035479829954815 0.03326565059259057 0.03019763774657662 0.028231656943684166 0.027971907804859367 0.027635497613244207 0.02609629843384177 0.026064155692279743 0.025395030573214513 0.025232232498192345 0.025102533258203302 0.024382188721835226 0.022860985466878727 0.02097931975492627 0.020299067881177458 0.019011719016095333 0.018746802137251297 0.017396830053945152 0.017009976425025684 0.017002557939223263 0.016787579636978956 0.016706351432038885 0.016493265311654737 0.016387404819675204 0.01573781403197213 0.015336488002839443 0.015063306051548358 0.014918958536443885 0.014849704688076955 0.014699767901876513 0.014640773670458946 0.014610665032930105 0.014578601029922452 0.014521910016482549 0.0143297258694309 0.014205242717723205 0.014138558934870992 0.014055728536150678 0.013835197883715615 0.0137139283754364 0.01364279396675939 0.013572606381602697 0.013462189532158995 0.013371623946791047 0.013370868318555934 0.013331717072994432 0.013290775730079973 0.01326718850160387 0.013179647050852505 0.013154076132504508 0.013102914616253401 0.013074860533430899 0.012954243194404764 0.012873372920293278 0.012850418758729347 0.01277105517889769 0.012707442643899327 0.012655120207127321 0.012616762227043523 0.012599101378908522 0.01259299729166396 0.012588890105825987 0.01258577805361547 0.01257871155758254 0.012575610148529757 0.012572432689884776 0.012566735821092382 0.012549295319679245 0.012524355958152355 0.01249004048503096 0.012449934035581754 0.012443987538412795 0.012437354180827578 0.012414698424235598 0.012411719018807145 0.012401350643362213 0.012383839195623256 0.012381112754200157 0.012359693838732717 0.012348994714147974 0.01233517340236208 0.012326393056010874 0.012325767048746163 0.012323807642349686 0.012320942722633808 0.012318723449328734 0.01231781931972275 0.012311857400181925 0.012305015664255797 0.012299505573035684 0.012295613793504022 0.012295084424239891 0.012294211214588123 0.012294211214588123 0.01229368229552369 0.012293644254650031 0.012293294115909585 0.012293186961553513 0.012293186961553513 0.012293186606085608 0.012293186606085608 0.012293186606085608 0.012293186606085608 0.012293186606085608 0.012293186606085608 0.012293186606085608 0.012293186606085608 0.012293177702191231 0.012293161044504872 0.012293161044504872 0.012293161044504872 0.012293157424806944 0.012293155368922823 0.012293155368922823 0.012293155368922823 0.012293147982021564 0.012293147982021564 0.012293147982021564 0.012293147982021564 0.012293147982021564 0.012293147982021564 0.012293088416335064 0.012293088416335064 0.012293088416335064 0.012293088416335064 0.012293088416335064 0.012293088416335064 0.012293088416335064 0.012293088416335064 0.012293078896235788 0.012293078896235788 0.012293078896235788 0.012293078896235788 0.012293078896235788 0.012293078896235788 0.012293078896235788 0.012293078896235788 0.012293078896235788 0.012293078896235788 0.012293078896235788 0.012293078896235788 0.012293078896235788 0.012293078896235788 0.012293078896235788 0.012293078896235788 0.012293078896235788 0.012293078896235788 0.012293078896235788 0.012293078896235788 0.012293078896235788 0.012293078896235788 0.012293078896235788 0.012293078896235788 0.012293078896235788 0.012293078896235788 0.012293078896235788 0.012293078896235788 0.012293071314443065 0.012293071314443065 0.012293071314443065 0.012293071314443065 0.012293071314443065 0.012293071314443065 0.012293071314443065 0.012293071314443065 0.012293071314443065 0.012293071314443065 0.012293071314443065 0.012293071314443065 0.012293071314443065 0.012293071314443065 0.012293071314443065 0.012293071314443065 0.012293071314443065 0.012293067948687 0.012293067948687 0.012293067948687 0.012293067230421872 0.012293067230421872 0.012293067230421872 0.012293067230421872 0.012293065250805242 0.012293065250805242 0.012293065250805242 0.012293065250805242 0.012293065250805242 0.012293065250805242 0.012293065250805242 0.012293065250805242 0.012293065064350097 0.012293065064350097 0.012293065064350097 0.012293065064350097 0.012293065064350097 0.012293065064350097 0.012293065064350097 0.012293065064350097 0.012293065064350097 0.012293065064350097 0.012293065064350097 0.012293054890810486 0.012293054890810486 0.012293054890810486 0.012293054890810486 0.012293054890810486 0.012293054890810486 0.012293054890810486 0.012293054890810486 0.012293054815777306 0.012293054815777306 0.012293054815777306 0.012293053805091185 0.012293053805091185 0.012293053805091185 0.01229305270020193 0.01229305270020193 0.012293046694111407 0.012293046694111407 0.012293046694111407 0.012293046694111407 0.012293046694111407 0.012293046694111407 0.012293046694111407 0.012293046694111407 0.012293046694111407 0.012293046694111407 0.012293046694111407 0.012293045729810697 0.012293045729810697 0.012293044977439554 0.012293044977439554 0.012293044977439554 0.012293044977439554 0.012293044977439554 0.012293044977439554 0.012293044977439554 0.012293044977439554 0.012293044977439554 0.012293044977439554 0.012293044977439554 0.012293044350276689 0.012293039029654817 0.012293039029654817 0.012293039029654817 0.012293039029654817 0.012293038375444487 0.012293038375444487 0.012293038375444487 0.012293038375444487 0.012293035745615377 0.012293035745615377 0.012293035745615377 0.012293035745615377 0.012293035745615377 0.012293035745615377 0.012293035745615377 0.012293035745615377 0.012293035745615377 0.012293035745615377 0.012293035745615377 0.012293032524686378 0.012293032524686378 0.012293032524686378 0.012293032524686378 0.012293029298193146 0.012293029298193146 0.012293029298193146 0.012293029298193146 0.012293029298193146 0.012293029298193146 0.012293029298193146 0.012293029298193146 0.012293029298193146 0.012293029298193146 0.012293029298193146 0.012293029298193146 0.012293029298193146 0.012293029001891291 0.012293029001891291 0.012293029001891291 0.012293025642992178 0.012293025642992178 0.012293025642992178 0.012293025642992178 0.012293025642992178 0.012293025642992178 0.012293025642992178 0.012293025305335945 0.012293025305335945 0.012293025305335945 0.012293025305335945 0.012293025305335945 0.012293025305335945 0.012293025305335945 0.012293022270074427 0.012293022270074427 0.012293022270074427 0.012293022270074427 0.012293022270074427 0.012293022270074427 0.012293022270074427 0.012293022270074427
We plot the data and the fit¶
In [15]:
def compute_density(h_in_km, best_x):
"""
returns the atmosheric density in kg.m^3
"""
p1 = np.array(best_x[:4])
p2 = np.array(best_x[4:8])
p3 = np.array(best_x[8:])
retval = 0.
for alpha,beta, gamma in zip(p1,p2, p3):
retval += alpha*np.exp(-(h_in_km-gamma)*beta)
return retval
In [16]:
plt.figure()
altitudes = np.linspace(78, 1500, 100)
plt.semilogy(altitudes, compute_density(altitudes, pop.champion_x), label ='analytical fit')
plt.semilogy(altitude, density, '.', label='data')
plt.legend()
Out[16]:
<matplotlib.legend.Legend at 0x7f5718864040>
In [17]:
with open("data/best_fit_density.pk", "wb") as file:
pkl.dump(pop.champion_x, file)