mirror of
https://gitlab.astro-wise.org/micado/atmosphericmodel.git
synced 2026-10-30 05:01:02 +00:00
251 lines
8.7 KiB
Python
251 lines
8.7 KiB
Python
#!/usr/bin/env Python3
|
|
import numpy as np
|
|
import sys
|
|
from AstroAtmosphere import *
|
|
import matplotlib.pyplot as plt
|
|
from netCDF4 import Dataset
|
|
from pdb import set_trace
|
|
from matplotlib import cm
|
|
from scipy.interpolate import interp1d, interp2d
|
|
from matplotlib.colors import LogNorm
|
|
import json
|
|
|
|
# Saves a variable, preferably a list or string to json text file.
|
|
# Useful for long simulations, where you're making a plot
|
|
def saveVariableToTXT(variable, filename):
|
|
variableFile = open(filename, 'w')
|
|
|
|
if type(variable) == np.ndarray:
|
|
json.dump(variable.tolist(), variableFile)
|
|
elif type(variable) in [list, str, float]:
|
|
json.dump(variable, variableFile)
|
|
else:
|
|
print('Not a list, array or string. Check variable type.')
|
|
|
|
variableFile.close()
|
|
return
|
|
|
|
# Load a variable from a json text file.
|
|
def loadVariableFromTXT(filename):
|
|
variableFile = open(filename, 'r')
|
|
variable = json.load(variableFile)
|
|
variableFile.close()
|
|
return variable
|
|
|
|
|
|
def calculateDistanceResolution(lat, lon):
|
|
# Calculate local radius of the earth
|
|
_A = 6378.137 # km
|
|
_B = 6356.752 # km
|
|
_phi = np.deg2rad(lat)
|
|
_rc0 = (_A*_B)**2 / (_A**2 * np.cos(_phi)**2 + _B**2 * np.sin(_phi)**2)**(3/2)
|
|
|
|
# Calculate circumference of the earth at constant latitude.
|
|
latR = 2 * np.pi * np.cos(_phi) * _rc0
|
|
distRes = latR/360 * lon
|
|
return distRes
|
|
|
|
def H_isotherm(T, lat):
|
|
'''
|
|
Calculates the reduced height of the atmosphere, assuming an isotherm distribution.
|
|
'''
|
|
_kb = 1.38064852e-23 # kg m2 s-2 K-1
|
|
_m = 4.809651698e-26 # kg (weight per molecule of air, assuming dry air)
|
|
|
|
_phi = np.deg2rad(lat)
|
|
_c1 = 5.2790414e-3
|
|
_c2 = 2.32718e-5
|
|
_c3 = 1.262e-7
|
|
_c4 = 7e-10
|
|
_g0 = 9.780327 # m/s^2
|
|
_g0_local = _g0 * (1 + _c1 * np.sin(_phi)**2 + _c2 * np.sin(_phi)**4 + _c3 * np.sin(_phi)**6 + _c4 * np.sin(_phi)**8)
|
|
return (_kb*T) / (_m * _g0_local)
|
|
|
|
|
|
def calculatePressure(Th, P0, h, lat):
|
|
redH = H_isotherm(Th, lat)
|
|
Ph = P0 * np.exp(-h/redH)
|
|
return Ph
|
|
|
|
plt.rc('font', family='serif')
|
|
plt.rc('axes', labelsize='large')
|
|
|
|
# Import the real data
|
|
nc = Dataset('air.201806.nc', 'r')
|
|
latitudes = np.asarray(nc.variables['lat'][:])
|
|
longitudes = np.asarray(nc.variables['lon'][:])
|
|
air_temp = nc.variables['air']
|
|
air_press = np.asarray(nc.variables['level'])
|
|
|
|
|
|
|
|
|
|
at = Observatory()
|
|
|
|
# Parameters at Cerro Armazones
|
|
T = air_temp[0,:,0,0][0] # k
|
|
P = air_press[0] # Pa
|
|
H = 0
|
|
xc = 450 # ppm
|
|
lat = 24.44768
|
|
h = 0 # m
|
|
l1 = 1.49 # micron
|
|
l2 = 1.78 # micron
|
|
|
|
rho = at.rho(p=P, T=T, RH=H, xc=xc)
|
|
disp = dispersion(lat, h)
|
|
disp.setReducedHeight(P, rho)
|
|
|
|
air_dens = at.rho(p=air_press, T=air_temp[0,:,0,0], RH=np.asarray([0]*len(air_press)), xc=np.asarray([0]*len(air_press)))
|
|
|
|
heights = disp.calcHeight(air_press, air_temp[0,:,0,0], RH=np.asarray([0]*len(air_press)), xc=np.asarray([0]*len(air_press)), P0=air_press[0]) / 1e3
|
|
distance_lat = calculateDistanceResolution(latitudes[100], longitudes[100,:]) - calculateDistanceResolution(latitudes[100,0], longitudes[100,0])
|
|
temp_selection_ground = air_temp[0,0,100,:]
|
|
temp_selection_upper = air_temp[0,-1,100,:]
|
|
print(distance_lat[1])
|
|
|
|
|
|
#########################################################
|
|
# Extrapolating using the U.S. Standard Atmosphere 1976.
|
|
# The end point of the real data is used to shift the USSA data accordingly (linearly)
|
|
h_ussa1976, _, _, _, air_temp_ussa1976, air_press_ussa1976,_,_,_,_ = np.genfromtxt('../US1976 Standard Atmosphere/si1py_modified.txt', unpack=True, skip_header=6)
|
|
temp_ip = interp1d(h_ussa1976, air_temp_ussa1976)
|
|
press_ip = interp1d(h_ussa1976, air_press_ussa1976)
|
|
|
|
ussa1976_start = temp_ip(heights[-1])
|
|
diff_temp_list = temp_selection_upper - ussa1976_start
|
|
start_i = np.argmax(h_ussa1976 > heights[-1])
|
|
ussa1976_temp_array = np.tile(air_temp_ussa1976[start_i:], (len(temp_selection_upper),1))
|
|
ussa1976_temp_array = np.asarray([ussa1976_temp_array[i] + diff_temp_list[i] for i in range(len(diff_temp_list))]).T
|
|
|
|
# Combining the (masked) arrays
|
|
new_heights = np.append(heights, h_ussa1976[start_i:])
|
|
tempData = np.append(air_temp[0,:,100,:], ussa1976_temp_array, axis=0)
|
|
tempData_mask = np.append(air_temp[0,:,100,:].mask, np.zeros(ussa1976_temp_array.shape), axis=0)
|
|
tempData = np.ma.masked_array(tempData, mask=tempData_mask)
|
|
|
|
# Now the same for the pressure data
|
|
# Combining the (masked) arrays
|
|
start_i = np.argmax(h_ussa1976 > heights[-1])
|
|
new_heights = np.append(heights, h_ussa1976[start_i:])
|
|
pressData = []
|
|
for i in range(len(temp_selection_upper)):
|
|
calc_press_data = calculatePressure(Th=ussa1976_temp_array[:,i], P0=air_press[0], h=h_ussa1976[start_i:]*1000, lat=lat)
|
|
pressData_i = np.append(air_press*100, calc_press_data*100, axis=0)
|
|
pressData.append(pressData_i)
|
|
pressData = np.asarray(pressData).T
|
|
#########################################################
|
|
|
|
# We'd now like to export a subset of the atmospheric data for a distance of ~200 km.
|
|
# Let's use interpolation to create a uniform gridded dataset
|
|
export_heights = np.linspace(0,new_heights.max(),400)
|
|
export_dist = np.linspace(0,distance_lat[136]-distance_lat[130],1000)
|
|
atmosphere_temp_func = interp2d(distance_lat[130:136]-distance_lat[130], new_heights, tempData[:,130:136], kind='linear')
|
|
export_temp_atmosphere = atmosphere_temp_func(export_dist, export_heights)
|
|
|
|
atmosphere_press_func = interp2d(distance_lat[130:136]-distance_lat[130], new_heights, pressData[:,130:136], kind='linear')
|
|
export_press_atmosphere = atmosphere_press_func(export_dist, export_heights)
|
|
|
|
# saveVariableToTXT(export_heights, 'export_heights.txt')
|
|
# saveVariableToTXT(export_dist, 'export_dist.txt')
|
|
# saveVariableToTXT(export_temp_atmosphere, 'export_temp_atmosphere.txt')
|
|
# saveVariableToTXT(export_press_atmosphere, 'export_press_atmosphere.txt')
|
|
|
|
|
|
# Plot the real data + extrapolated standard atmosphere.
|
|
# fig = plt.figure()
|
|
# ax = fig.add_subplot(211)
|
|
# ax2 = fig.add_subplot(212)
|
|
# im1 = ax.pcolormesh(export_dist, export_heights, export_temp_atmosphere, cmap=cm.coolwarm)
|
|
# im2 = ax2.pcolormesh(export_dist, export_heights, export_press_atmosphere, cmap=cm.coolwarm, norm=LogNorm())
|
|
# ax2.set_xlabel('Lateral distance (km)')
|
|
# ax.set_ylabel('Height (km)')
|
|
# ax2.set_ylabel('Height (km)')
|
|
# cbar = fig.colorbar(im1, ax=ax, label='Temperature (K)')
|
|
# lbls = cbar.get_ticks()
|
|
# cbar.ax.invert_yaxis()
|
|
# cbar.set_ticks(lbls[::-1])
|
|
# cbar2 = fig.colorbar(im2, ax=ax2, label='Pressure (Pa)')
|
|
# lbls = cbar2.get_ticks()
|
|
# cbar2.ax.invert_yaxis()
|
|
# cbar2.set_ticks(lbls[::-1])
|
|
# plt.tight_layout()
|
|
# plt.show(block=False)
|
|
# set_trace()
|
|
|
|
|
|
|
|
|
|
# # Plot the real data + extrapolated standard atmosphere.
|
|
# fig = plt.figure()
|
|
# ax = fig.add_subplot(211)
|
|
# ax2 = fig.add_subplot(212)
|
|
# im1 = ax.imshow(tempData[:,130:137], cmap=cm.coolwarm, interpolation= 'bicubic', aspect='auto', origin='lower')
|
|
# # im1 = ax.pcolormesh(distance_lat[130:137], new_heights, tempData[:,130:137], cmap=cm.coolwarm)
|
|
# im2 = ax2.pcolormesh(distance_lat, new_heights, pressData, cmap=cm.coolwarm, norm=LogNorm())
|
|
# ax2.set_xlabel('Lateral distance (km)')
|
|
# ax.set_ylabel('Height (km)')
|
|
# ax2.set_ylabel('Height (km)')
|
|
# cbar = fig.colorbar(im1, ax=ax, label='Temperature (K)')
|
|
# lbls = cbar.get_ticks()
|
|
# cbar.ax.invert_yaxis()
|
|
# cbar.set_ticks(lbls[::-1])
|
|
# cbar2 = fig.colorbar(im2, ax=ax2, label='Pressure (Pa)')
|
|
# lbls = cbar2.get_ticks()
|
|
# cbar2.ax.invert_yaxis()
|
|
# cbar2.set_ticks(lbls[::-1])
|
|
# plt.tight_layout()
|
|
# plt.show(block=False)
|
|
|
|
|
|
# Plot the real data and properties
|
|
fig = plt.figure(figsize=(9,4))
|
|
fig2 = plt.figure()
|
|
ax1 = fig2.add_subplot(221)
|
|
ax2 = fig2.add_subplot(222)
|
|
ax3 = fig2.add_subplot(223)
|
|
ax4 = fig.add_subplot(121)
|
|
ax5 = fig2.add_subplot(224)
|
|
ax6 = fig.add_subplot(122)
|
|
|
|
ax1.grid(alpha=0.3)
|
|
ax2.grid(alpha=0.3)
|
|
ax3.grid(alpha=0.3)
|
|
ax5.grid(alpha=0.3)
|
|
|
|
ax1.plot(air_temp[0,:,0,0], heights)
|
|
ax1.set_xlabel('Temperature (K)')
|
|
ax1.set_ylabel('Height (km)')
|
|
|
|
ax2.plot(air_press, heights)
|
|
ax2.set_xlabel('Pressure (Pa)')
|
|
ax2.set_ylabel('Height (km)')
|
|
|
|
ax3.plot(air_dens, heights)
|
|
ax3.set_xlabel('Density (g cm$^{-3}$)')
|
|
ax3.set_ylabel('Height (km)')
|
|
|
|
im = ax4.imshow(np.rot90(air_temp[0,0,:,:].T), aspect=1, cmap='jet')
|
|
ax4.set_yticks([])
|
|
ax4.set_xticks([])
|
|
# fig.colorbar(im, ax=ax4, label='Temperature (K)')
|
|
|
|
ax5.plot(distance_lat, temp_selection_ground, 'b-', label='ground')
|
|
ax5.plot(distance_lat, temp_selection_upper, 'r-', label='upper atmosphere')
|
|
ax5.legend()
|
|
ax5.set_ylabel('Temperature (K)')
|
|
ax5.set_xlabel('Distance (km)')
|
|
|
|
im2 = ax6.imshow(air_temp[0,:,100,:], aspect='auto', origin='lower', interpolation='bilinear', extent=[0,distance_lat[-1],0,heights[-1] ], cmap='jet')
|
|
ax6.set_xlabel('Distance (km)', fontsize=18)
|
|
ax6.set_ylabel('Height (km)', fontsize=18)
|
|
cbar = fig.colorbar(im2, ax=ax6)#, label='Temperature (K)')
|
|
cbar.ax.tick_params(labelsize=16)
|
|
cbar.set_label('Temperature (K)', size=18)
|
|
lbls = cbar.get_ticks()
|
|
cbar.ax.invert_yaxis()
|
|
cbar.set_ticks(lbls[::-1])
|
|
ax6.tick_params(axis ='both', direction='in', labelsize=16)
|
|
|
|
plt.tight_layout()
|
|
plt.show() |