add scripts

README.md

@@ -1,3 +1,7 @@
+# Modeling Time-Variable Elemental Abundances in Coronal Loop Simulations
+
+A paper by Jeffrey Reep, John Unverferth, Will Barnes, and Sherry Chhabra describing a method to add time-variable elemental abundances due to chromospheric evaporation into coronal loop simulations.
+
 <p align="center">
 <a href="https://github.com/showyourwork/showyourwork">
 <img width = "450" src="https://raw.githubusercontent.com/showyourwork/.github/main/images/showyourwork.png" alt="showyourwork"/>

src/scripts/plot.py

@@ -0,0 +1,90 @@
+"""
+Plots the temperature, density, and abundance factor from an ebtel++ simulation
+
+"""
+import paths
+import numpy as np
+from matplotlib import pyplot as plt
+
+def plot_figures(time, temperature, density, L, H, t, filename=None, xlim=None):
+    """
+    Plot the temperature, density, and abundance factor as a function of time for
+    a given set of simulations
+    
+    Parameters
+    ----------
+    time : 
+        The time-step of a given simulation
+    temperature : 
+        The coronal-averaged electron temperature [K] 
+    density : 
+        The coronal-averaged density [cm**-3]
+    L : `str`
+        The loop length [Mm]
+    H : `str`
+        The maximum heating rate [erg s**-1 cm**-3]
+    t : `str`
+        The total duration of a heating pulse [s]
+    filename : `str`, default: `None`
+        Format for the filename of each figure
+    xlim : 
+        The limits of the x-axis
+    
+    Returns
+    -------
+    None
+
+    """
+    if filename is None:
+        filename = 'L'+L+'_H'+H+'_t'+t+'.png'
+
+    plt.figure(0)
+    plt.plot(time[0][:]/60., temperature[0][:]/1e6, label='Time variable')
+    plt.plot(time[1][:]/60., temperature[1][:]/1e6, label='Photospheric', linestyle='dashed')
+    plt.plot(time[2][:]/60., temperature[2][:]/1e6, label='Coronal', linestyle='dashed')
+    plt.ylabel('Electron Temperature [MK]', fontsize=14)
+    plt.xlabel('Time [min]', fontsize=14)
+    plt.title(r'Length = '+L+' Mm, Heat = '+H+' erg s$^{-1}$ cm$^{-3}$, time = '+t+' s')
+    plt.grid('-')
+    if xlim is not None:
+        plt.xlim(xlim)
+    if H == '0.01':
+        plt.legend(fontsize=18)
+    plt.savefig(paths.figures / ('temperature_' + filename), dpi=300)   
+    plt.close()
+
+    plt.figure(1)
+    plt.plot(time[0][:]/60., density[0][:]/1e8, label='Time variable')
+    plt.plot(time[1][:]/60., density[1][:]/1e8, label='Photospheric', linestyle='dashed')
+    plt.plot(time[2][:]/60., density[2][:]/1e8, label='Coronal', linestyle='dashed')
+    plt.ylabel(r'Density [$10^{8}$ cm$^{-3}$]', fontsize=14)
+    plt.xlabel('Time [min]', fontsize=14)
+    plt.title(r'Length = '+L+' Mm, Heat = '+H+' erg s$^{-1}$ cm$^{-3}$, time = '+t+' s')
+    plt.grid('-')
+    if xlim is not None:
+        plt.xlim(xlim)
+    plt.savefig(paths.figures / ('density_' + filename), dpi=300)
+    plt.close()
+
+    # the abundance factor is not output from ebtel++, so we recalculate it here
+    initial_af = 4.0
+    af = 1.0 + (initial_af - 1.0) * (density[0][0] / density[0][:]);
+    max_n = density[0][0]
+    for i in range(len(af)):
+        if density[0][i] > max_n:
+            max_n = density[0][i]
+        af[i] = 1.0 + (initial_af - 1.0) * (density[0][0] / max_n);
+
+    plt.figure(2)
+    plt.plot(time[0][:]/60., af[:], label='Time variable', )
+    plt.plot(time[1][:]/60., np.full(len(time[1][:]), 1.0), label='Photospheric', linestyle='dashed')
+    plt.plot(time[2][:]/60., np.full(len(time[2][:]), 4.0), label='Coronal', linestyle='dashed')
+    plt.ylabel(r'Abundance Factor', fontsize=14)
+    plt.xlabel('Time [min]', fontsize=14)
+    plt.title(r'Length = '+L+' Mm, Heat = '+H+' erg s$^{-1}$ cm$^{-3}$, time = '+t+' s')
+    plt.ylim(0,5)
+    plt.grid('-')
+    if xlim is not None:
+        plt.xlim(xlim)
+    plt.savefig(paths.figures / ('abundance_' + filename), dpi=300)
+    plt.close()

src/scripts/read.py

@@ -0,0 +1,59 @@
+"""
+Reads the data from the ebtel++ output files
+
+"""
+import numpy as np
+
+def read_ebtel_file(filename='ebtel++_results_file.txt'):
+    """
+    Reads in the output data file from an ebtel++ simulation
+    
+    Parameters
+    ----------
+    filename : `str`, default: `ebtel++_results_file.txt`
+        The name of the data file
+    
+    Returns
+    -------
+    time:
+        An array of the time-steps in the simulation [s]
+    electron_temperature:
+        A time-array of the coronal-averaged electron temperature [K]
+    ion_temperature:
+        A time-array of the coronal-averaged ion temperature [K]
+    density:
+        A time-array of the coronal-averaged density [cm**-3]
+    electron_pressure:
+        A time-array of the coronal-averaged electron pressure [dyn cm**-2]
+    ion_pressure:
+        A time-array of the coronal-averaged ion pressure [dyn cm**-2]
+    velocity:
+        A time-array of the bulk flow velocity [cm s**-1]
+    heating_rate:
+        A time-array of the coronal-averaged heating rate [erg cm**-3 s**-1]
+    
+    """
+    f = open(filename, 'r')
+    time = np.empty(0)
+    electron_temperature = np.empty(0)
+    ion_temperature = np.empty(0)
+    density = np.empty(0)
+    electron_pressure = np.empty(0)
+    ion_pressure = np.empty(0)
+    velocity = np.empty(0)
+    heating_rate = np.empty(0)
+    for line in f:
+        line = line.strip()
+        columns = line.split()
+        time = np.append(time, float(columns[0]))
+        electron_temperature = np.append(electron_temperature, float(columns[1]))
+        ion_temperature = np.append(ion_temperature, float(columns[2]))
+        density = np.append(density, float(columns[3]))
+        electron_pressure = np.append(electron_pressure, float(columns[4]))
+        ion_pressure = np.append(ion_pressure, float(columns[5]))
+        velocity = np.append(velocity, float(columns[6]))
+        heating_rate = np.append(heating_rate, float(columns[7]))     
+
+    f.close()
+
+    return time, electron_temperature, ion_temperature, density, electron_pressure, ion_pressure, velocity, heating_rate

src/scripts/render.py

@@ -0,0 +1,55 @@
+"""
+Script to render the data into figure data
+
+"""
+import numpy as np
+import paths
+from plot import plot_figures
+from read import read_ebtel_file
+
+
+# Let's calculate some cases with single heating events first
+lengths = ['40', '80']
+heat = ['1.0','0.1','0.03','0.01']
+duration = ['20']
+
+for L in lengths:
+    for H in heat:
+        for t in duration:
+            vfile = 'var_L'+L+'_H'+H+'_t'+t+'.txt'
+            pfile = 'pho_L'+L+'_H'+H+'_t'+t+'.txt'
+            cfile = 'cor_L'+L+'_H'+H+'_t'+t+'.txt'
+
+            t_v, T_e_v, T_i_v, n_v, P_e_v, P_i_v, v_v, Q_v = read_ebtel_file(filename=paths.data / vfile)
+            t_p, T_e_p, T_i_p, n_p, P_e_p, P_i_p, v_p, Q_p = read_ebtel_file(filename=paths.data / pfile)
+            t_c, T_e_c, T_i_c, n_c, P_e_c, P_i_c, v_c, Q_c = read_ebtel_file(filename=paths.data / cfile)
+
+            time = [t_v, t_p, t_c]
+            temperature = [T_e_v, T_e_p, T_e_c]
+            density = [n_v, n_p, n_c]
+
+            plot_figures(time, temperature, density, L, H, t)
+
+
+# Let's calculate a couple of nanoflare trains for comparison
+lengths = ['40', '80']
+heat = ['0.01']
+duration = ['20']
+for L in lengths:
+    for H in heat:
+        for t in duration:
+            vfile = 'var_train_L'+L+'_H'+H+'_t'+t+'.txt'
+            pfile = 'pho_train_L'+L+'_H'+H+'_t'+t+'.txt'
+            cfile = 'cor_train_L'+L+'_H'+H+'_t'+t+'.txt'
+
+            t_v, T_e_v, T_i_v, n_v, P_e_v, P_i_v, v_v, Q_v = read_ebtel_file(filename=paths.data / vfile)
+            t_p, T_e_p, T_i_p, n_p, P_e_p, P_i_p, v_p, Q_p = read_ebtel_file(filename=paths.data / pfile)
+            t_c, T_e_c, T_i_c, n_c, P_e_c, P_i_c, v_c, Q_c = read_ebtel_file(filename=paths.data / cfile)
+
+            time = [t_v, t_p, t_c]
+            temperature = [T_e_v, T_e_p, T_e_c]
+            density = [n_v, n_p, n_c]
+
+            plot_figures(time, temperature, density, L, H, t, 
+                            filename = 'train_L'+L+'_H'+H+'_t'+t+'.png',
+                            xlim = [0, 40])