Merge branch 'master' of https://github.com/ideasrule/pyTPCI

Author: ideasrule

LICENSE

@@ -0,0 +1,30 @@
+pyTPCI is an improved version of The PLUTO-CLOUDY Interface, coupling together hydrodynamics code [PLUTO](https://plutocode.ph.unito.it/) and gas microphysics code [CLOUDY](https://gitlab.nublado.org/cloudy/cloudy).
+
+## Setup
+1) After cloning the pyTPCI directory to your computer, copy the `wrapper` folder for **each instance of pyTPCI you wish to run**. 
+2) Define the environment variable `PLUTO_DIR` and run `$PLUTO_DIR/setup.py` in your folder. Further customization of PLUTO takes place in `pluto_template.ini`.
+3) Run `make` and `make clean` in your folder to produce the PLUTO executable.
+4) Run `make` in `cloudy/source` to produce `cloudy.exe`.
+
+pyTPCI runs PLUTO and then CLOUDY in alternating steps, transferring heating information between them with `write_heating_file`. In order to save computational power and improve speed, 
+by default pyTPCI only calls CLOUDY if there is at least a 10% maximum fractional difference (`max_rel_diff`) in density or pressure between PLUTO files. 
+
+## Running pyTPCI
+1) Create a stellar spectrum at the planet's surface (semimajor axis), with columns of log10(frequency (Hz)) and log10(F_nu), spectral irradiance (ergs/(s cm^2 Hz)). At the end this also contains the band-integrated flux at the planet's surface, and the energy range in Rydbergs. See `spectra_example.ini`.
+2) Edit `tpci.py` to set the system parameters: the stellar spectrum file, the semimajor axis, planet radius and mass, stellar mass, initial temperature, and metallicity. Upon running, this will generate a CLOUDY `params.h` file and an input script ending in `.in`.
+3) Open `screen` and run `python tpci.py` to start at t=0 and use the default timestep dt=0.01. If restarting from a previous pyTPCI run, instead do `python tpci.py global_ind time timestep`, where `global_ind` is the number of the CLOUDY file you wish to start from.
+4) Monitor pyTPCI's progress by plotting PLUTO and CLOUDY files using `plot.py` and `cloudy_plot.py`.
+
+## Plotting information
+- To plot one file, run `python plot.py global_ind` where `global_ind` is the number of the file, such as "54".
+- To plot multiple files, run `python plot.py start stop step`, where `start` and `stop` are the initial and final file numbers to plot, with a given `step`. Plotting files will automatically skip file names that do not exist. 
+
+## Logging and Docs information
+- pyTPCI run information is in `tpci_log.txt`.
+- PLUTO output logs for each PLUTO run are appended in `pluto_log.txt`
+- CLOUDY docs are in `cloudy/docs` and PLUTO docs are in `PLUTO/Doc`
+
+## Troubleshooting and Common Errors
+- Is the path to the `spectra.ini` file correct?
+- Does `definitions.h` contain `#include params.h` at the top?
+- If CLOUDY crashes quickly, try using a smaller dt=0.001

README.md

@@ -1,12 +1,31 @@
 import numpy as np
 import matplotlib.pyplot as plt
 import sys
+import os
 
-unit_length = 2.9 * 6.378136e8
+# replace with planet radius (Earth radii)
+unit_length = 12.54 * 6.378136e8
 
-for filename in sys.argv[1:]:
-    depth, rho, n, mu, Te = np.loadtxt(filename, usecols=range(5), unpack=True)
+
+# start, [stop, step]
+if len(sys.argv) == 2:
+    start = str(sys.argv[1])
+    files = ['cl_data.'+start.zfill(4)+'.over.tab']
+else:
+    start, stop, step = int(sys.argv[1]), int(sys.argv[2]), int(sys.argv[3])
+    file_nums = list(map(str, np.arange(start, stop+1, step)))
+    files = []
+    for i in range(len(file_nums)):
+        fname = 'cl_data.'+file_nums[i].zfill(4)+'.over.tab'
+        if os.path.exists(fname): # skip files that don't exist
+            files.append(fname)
+
+# oldest is purple, proceeds through rainbow to newest (red)
+colors = plt.cm.rainbow(np.linspace(0, 1, len(files)))
+
+for i, file in enumerate(files):
+    depth, rho, n, mu, Te = np.loadtxt(file, usecols=range(5), unpack=True)
     cloudy_radii = 1 + (np.max(depth) - depth) / unit_length
-    plt.plot(cloudy_radii, Te)
+    plt.plot(cloudy_radii, Te, color=colors[i])
 
 plt.show()

hd209458/tpci.py

@@ -1,27 +1,49 @@
 import numpy as np
 import matplotlib.pyplot as plt
 import sys
+import os
 
-for i, filename in enumerate(sys.argv[1:]):
-    #if i % 2 != 1:
-    #    continue
+# PLUTO units
+UNIT_DENSITY = 1e-10
+UNIT_VELOCITY = 1e5
+UNIT_PRESSURE = UNIT_DENSITY * UNIT_VELOCITY**2
+
+# start [stop, step]
+if len(sys.argv) == 2:
+    start = str(sys.argv[1])
+    files = ['data.'+start.zfill(4)+'.tab']
+else:
+    start, stop, step = int(sys.argv[1]), int(sys.argv[2]), int(sys.argv[3])
+    file_nums = list(map(str, np.arange(start, stop+1, step)))
+    files = []
+    for i in range(len(file_nums)):
+        fname = 'data.'+file_nums[i].zfill(4)+'.tab'
+        if os.path.exists(fname): # skip files that don't exist
+            files.append(fname)
+
+# oldest is purple, proceeds through rainbow to newest (red)
+colors = plt.cm.rainbow(np.linspace(0,1, len(files)))
+
+for i, filename in enumerate(files):
 
     r, rho, v, P = np.loadtxt(filename, usecols=(0,2,3,6), unpack=True)
-    #r, rho, v = np.loadtxt(filename, usecols=(0,2,3), unpack=True)
     plt.figure(0)
-    plt.loglog(r, rho)
+    plt.loglog(r, rho*UNIT_DENSITY, color=colors[i])
+    plt.ylabel("Density (g/cm3)")
     plt.title("rho")
 
     plt.figure(1)
-    plt.semilogx(r, v)
+    plt.loglog(r, v, color=colors[i])
+    plt.ylabel("Velocity (km/s)")
     plt.title("v")
 
     plt.figure(2)
-    plt.loglog(r, P)
+    plt.loglog(r, P*UNIT_PRESSURE, color=colors[i])
+    plt.ylabel("Pressure (g/(cm s^2))")
     plt.title("P")
 
     plt.figure(3)
-    plt.loglog(r, P*2.3*1.67e-24/(1.38e-16*rho*1e-10))
-    plt.title("T?")
+    plt.loglog(r, P*2.3*1.67e-24/(1.38e-16*rho*1e-10), color=colors[i])
+    plt.title("Estimated Temperature (K)")
 
 plt.show()