the more complete solution

Author: aaronjridley

day_8/Solution/helper_functions.py

@@ -0,0 +1,84 @@
+#!/usr/bin/env python
+
+import numpy as np
+import re
+from datetime import datetime
+import matplotlib.pyplot as plt
+
+#-----------------------------------------------------------------------------
+# Read in the EUV CSV file.  Get out:
+#   - short and long wavelengths
+#   - EUV parameters
+#   - Absorption cross sections for N2, O2, and O
+#-----------------------------------------------------------------------------
+
+def read_euv_csv_file(file):
+
+    fpin = open(file, 'r')
+
+    iColStart = 5
+    iColEnd = -1
+
+    for line in fpin:
+        cols = line.split(',')
+        if (cols[0] == 'Short'):
+            short = np.asarray(cols[iColStart:iColEnd], dtype = float)
+        if (cols[0] == 'Long'):
+            long = np.asarray(cols[iColStart:iColEnd], dtype = float)
+        if (cols[0] == 'F74113'):
+            f74113 = np.asarray(cols[iColStart:iColEnd], dtype = float)
+        if (cols[0] == 'AFAC'):
+            afac = np.asarray(cols[iColStart:iColEnd], dtype = float)
+        if (cols[0] == 'N2' and cols[2] == 'abs'):
+            n2cs = np.asarray(cols[iColStart:iColEnd], dtype = float)
+        if (cols[0] == 'O2' and cols[2] == 'abs'):
+            o2cs = np.asarray(cols[iColStart:iColEnd], dtype = float)
+        if (cols[0] == 'O' and cols[2] == 'abs'):
+            ocs = np.asarray(cols[iColStart:iColEnd], dtype = float)
+    
+    fpin.close()
+
+    data = {'short': short * 1e-10,
+            'long': long * 1e-10,
+            'f74113': f74113 * 1e9 * 1e4, # /cm2 to /m2
+            'afac': afac,
+            'ocross': ocs * 1e-22,
+            'o2cross': o2cs * 1e-22,
+            'n2cross': n2cs * 1e-22}
+
+    return data
+
+#-----------------------------------------------------------------------------
+# Plot out spectrum to a file
+#-----------------------------------------------------------------------------
+
+def plot_spectrum(wavelengths_in_m, intensities, filename):
+
+    fig = plt.figure(figsize = (10,5))
+    ax = fig.add_subplot(111)
+
+    ax.scatter(wavelengths_in_m * 1e10, intensities)
+    ax.set_xlabel('Wavelength (A)')
+    ax.set_ylabel('Intensity (photons/m2/s)')
+
+    fig.savefig(filename)
+    plt.close()
+    
+#-----------------------------------------------------------------------------
+# Plot var as a function of altitutde
+#-----------------------------------------------------------------------------
+
+def plot_value_vs_alt(alts, values, filename, var_name, is_log = False):
+
+    fig = plt.figure(figsize = (5,10))
+    ax = fig.add_subplot(111)
+
+    ax.plot(values, alts)
+    ax.set_xlabel(var_name)
+    ax.set_ylabel('Altitude (km)')
+    if (is_log):
+        ax.set_xscale('log')
+    
+    fig.savefig(filename)
+    plt.close()
+    

day_8/Solution/physics_functions_complete.py

@@ -0,0 +1,183 @@
+#!/usr/bin/env python
+
+import numpy as np
+
+#-----------------------------------------------------------------------------
+# Define some useful constants
+#-----------------------------------------------------------------------------
+
+cKB_ = 1.38064852e-23
+cRE_ = 6371.0
+cG_ = 9.81
+cAMU_ = 1.6726219e-27
+cH_ = 6.6261e-34
+cC_ = 2.9979e8
+cE_ = 1.6022e-19
+cD2R_ = 3.1415927 / 180.0
+
+#-----------------------------------------------------------------------------
+# calculate the scale height, which is KT/mg, where:
+#   k - boltzmann's constant
+#   t = temperature
+#   m = mass
+#   g = gravity
+#  Feed in 
+#-----------------------------------------------------------------------------
+
+def calc_scale_height(mass_in_amu, alt_in_km, temp_in_k):
+
+    gravity = cG_ * (cRE_ / (cRE_ + alt_in_km))
+    mass_in_kg = mass_in_amu * cAMU_
+    h_in_km = cKB_ * temp_in_k / mass_in_kg / gravity / 1000.0;
+
+    return h_in_km
+
+#-----------------------------------------------------------------------------
+# Calculate a hydrostatic solution:
+#  n(i+1) = n(i) * t(i)/t(i+1) * exp(dz / H)
+#  where:
+#  n = density
+#  t = temperature
+#  dz = delta-z (delta altitude)
+#  H = scale-height
+#-----------------------------------------------------------------------------
+
+def calc_hydrostatic(density_bc, scale_height_in_km, temp_in_k, alt_in_km):
+
+    nAlts = len(alt_in_km)
+    density = np.zeros(nAlts)
+
+    # define our boundary condition
+    density[0] = density_bc
+    for iAlt in range(1, nAlts):
+        tRatio = temp_in_k[iAlt - 1] / temp_in_k[iAlt]
+        dz = alt_in_km[iAlt] - alt_in_km[iAlt - 1]
+        h = scale_height_in_km[iAlt]
+        density[iAlt] = density[iAlt-1] * tRatio * np.exp(-dz / h)
+
+    return density
+
+#-----------------------------------------------------------------------------
+# EUVAC - take F107 and F107a and return solar spectrum
+# See Schunk and Nagy, page 259
+#-----------------------------------------------------------------------------
+
+def EUVAC(f74113, afac, f107, f107a):
+    p = (f107 + f107a)/2.0
+    if (p < 80):
+        p = 80.0
+    intensity = f74113 * (1.0 + afac * (p-80))
+    return intensity
+        
+#-----------------------------------------------------------------------------
+# Take a wavelength in m and convert it to energy in Joules:
+# e = plank's constant * speed of light / wavelength
+#-----------------------------------------------------------------------------
+
+def convert_wavelength_to_joules(wavelength):
+    energy = cH_ * cC_ / wavelength
+    return energy
+
+#-----------------------------------------------------------------------------
+# Calculate Tau for a given species, given:
+#  - solar zenity angle (in degrees)
+#  - density in /m3
+#  - scale height in km
+#  - cross sections in /m2
+#-----------------------------------------------------------------------------
+
+def calc_tau_complete(SZA_in_deg, density_in_m3, scale_height_in_km, cross_section):
+
+    # We are only going to do this for a single wavelength for now!
+
+    cs = cross_section[0]
+
+    # convert scale height to m:
+    h = scale_height_in_km * 1000.0
+
+    # convert SZA to radians:
+    sza = SZA_in_deg * cD2R_
+
+    # calculate integrated density (which is density * scale height):
+    integrated_density = density_in_m3 * h
+
+    nWaves = len(cross_section)
+    nAlts = len(density_in_m3)
+    tau = np.zeros((nWaves, nAlts))
+
+    # calculate Tau:
+    for iWave in range(nWaves):
+        tau[iWave][:] = integrated_density * cross_section[iWave]
+
+    return tau
+
+#-----------------------------------------------------------------------------
+# Calculate energy deposition as a function of altitude, given:
+#  - density of a given species in /m3 as a function of altitude
+#  - intensity at infinity as a function of wavelength
+#  - tau as a function of wavelength and altitude
+#  - cross_section of the given species as a function of wavelength
+#  - energies of the wavelengths
+#  - efficiency of the heating (say 30%)
+#-----------------------------------------------------------------------------
+
+def calculate_Qeuv_complete(density_in_m3,
+                            intensity_inf,
+                            tau,
+                            cross_section,
+                            energies,
+                            efficiency):
+
+    nAlts = len(density_in_m3)
+    nWaves = len(intensity_inf)
+
+    Qeuv = np.zeros(nAlts)
+
+    for iWave in range(nWaves):
+        # intensity is a function of altitude (for a given wavelength):
+        intensity = intensity_inf[iWave] * np.exp(-tau[iWave][:])
+        Qeuv = Qeuv + \
+            efficiency * \
+            density_in_m3 * \
+            intensity * \
+            cross_section[iWave] * \
+            energies[iWave]
+
+    return Qeuv
+
+#-----------------------------------------------------------------------------
+# calculate rho given densities of O (and N2 and O2)
+#-----------------------------------------------------------------------------
+
+def calc_rho(density_o, mass_o):
+    rho = density_o * mass_o * cAMU_
+    return rho
+
+#-----------------------------------------------------------------------------
+# calculate cp.  We have hard coded this to be 1500, which is for O, but
+# we could pass densities and vibrational states and get the real Cp.
+#-----------------------------------------------------------------------------
+
+def calculate_cp():
+    cp = 1500.0
+    return cp
+
+#-----------------------------------------------------------------------------
+# calculate dT/dt from Q and rho:
+#-----------------------------------------------------------------------------
+
+def convert_Q_to_dTdt(Qeuv, rho, cp):
+    dTdt = Qeuv / (rho * cp)
+    return dTdt
+
+#-----------------------------------------------------------------------------
+# Calculate the real rho (ignore!!!)
+#-----------------------------------------------------------------------------
+
+def calc_rho_complete(density_o, mass_o,
+                      density_o2, mass_o2,
+                      density_n2, mass_n2):
+    rho = density_o * mass_o * cAMU_
+    rho += density_o2 * mass_o2 * cAMU_
+    rho += density_n2 * mass_n2 * cAMU_
+    return rho