Merge pull request #219 from VirtualPlanetaryLaboratory/dev

Dev

examples/FfdProxCen/1Myr/1Myr.in → examples/FfdProxCen/davenport/star.in

@@ -1,5 +1,5 @@
 #
-sName	                           daven1Myr
+sName	                           star
 saModules	                    flare
 
 #stellar module
@@ -13,4 +13,4 @@ sFlareFFD                           DAVENPORT   # flare frenquency distribution
 dEnergyBin                       4                          # Number of bins between min and max energy
 
 #ouput parameters
-saOutputOrder Time -FlareFreq1 -FlareFreq2 -FlareFreq3 -FlareFreq4 -FlareEnergy1 -FlareEnergy2 -FlareEnergy3 -FlareEnergy4
+saOutputOrder Age -FlareFreq1 -FlareFreq2 -FlareFreq3 -FlareFreq4 -FlareEnergy1 -FlareEnergy2 -FlareEnergy3 -FlareEnergy4

examples/FfdProxCen/1Myr/vpl.in → examples/FfdProxCen/davenport/vpl.in

@@ -1,8 +1,8 @@
 #
-sSystemName           FfdDaven1
+sSystemName           davenport
 iVerbose                      5
 bOverwrite                1
-saBodyFiles              1Myr.in
+saBodyFiles              star.in
 sUnitMass                 solar
 sUnitLength               AU
 sUnitTime                 YEARS
@@ -15,4 +15,3 @@ bVarDt                       1
 dEta                            0.01
 dStopTime                1.1e9
 dOutputTime               1e5
-

examples/FfdProxCen/lacy/lacy.in → examples/FfdProxCen/lacy/star.in

@@ -1,5 +1,5 @@
 #
-sName	                            lacy
+sName	                            star
 saModules	                    flare
 
 #stellar module
@@ -11,9 +11,8 @@ dFlareMaxEnergy           -1.0e37              # Max flare energy
 sFlareFFD                           LACY                 # flare frenquency distribution evolution model
 dEnergyBin                       4                         # Number of bins between min and max energy
 dFlareSlope                      -0.68
-dFlareYInt                        20.9
+dFlareYInt                        -20.9
 #sFlareBandPass              KEPLER #(4000–9000 Å) UV (3000–4300 Å) GOES band (1-8 Å) SXR (1.24 - 1239.85 Å) BOLOMETRIC (Osten and Wolk (2015)) TESSUV (2000-2800 Å)
-#sFlareSlopeUnits            SEC #MIN HOUR DAY #units of the FFD input slopes
 
 #ouput parameters
-saOutputOrder Time -FlareFreq1 -FlareFreq2 -FlareFreq3 -FlareFreq4 -FlareEnergy1 -FlareEnergy2 -FlareEnergy3 -FlareEnergy4
+saOutputOrder Age -FlareFreq1 -FlareFreq2 -FlareFreq3 -FlareFreq4 -FlareEnergy1 -FlareEnergy2 -FlareEnergy3 -FlareEnergy4

examples/FfdProxCen/lacy/vpl.in

@@ -1,8 +1,8 @@
 #
-sSystemName           FfdLacy
+sSystemName           lacy
 iVerbose                      5
 bOverwrite                1
-saBodyFiles              lacy.in
+saBodyFiles              star.in
 sUnitMass                 solar
 sUnitLength               AU
 sUnitTime                 YEARS
@@ -15,4 +15,3 @@ bVarDt                       1
 dEta                            0.01
 dStopTime                1.1e9
 dOutputTime               1e5
-

examples/FfdProxCen/makeplot.py

@@ -41,71 +41,72 @@
 path = pathlib.Path(__file__).parents[0].absolute()
 sys.path.insert(1, str(path.parents[0]))
 
-# mpl.style.use('classic')
 
 run = [
-    "./1Myr/vpl.in",
-    "./10Myr/vpl.in",
-    "./100Myr/vpl.in",
-    "./1000Myr/vpl.in",
+    "./davenport/vpl.in",
     "./lacy/vpl.in",
 ]
 
 # Run the simulations
-for i in range(0, 5):
+for i in range(0, 2):
     vplanet.run(path / run[i], units=False)
 
 # Loading the data
-M = np.loadtxt("./lacy/FfdLacy.lacy.forward")
-MD = np.loadtxt("./1Myr/FfdDaven1.daven1Myr.forward")
-ND = np.loadtxt("./10Myr/FfdDaven10.daven10Myr.forward")
-OD = np.loadtxt("./100Myr/FfdDaven100.daven100Myr.forward")
-PD = np.loadtxt("./1000Myr/FfdDaven1000.daven1000Myr.forward")
 
+data_daven = np.loadtxt("./davenport/davenport.star.forward")
+data_lacy = np.loadtxt("./lacy/lacy.star.forward")
 
-# Organizing the data
-Mdata = []
-MdataD = []
-NdataD = []
-OdataD = []
-PdataD = []
+age_daven = data_daven[:, 0]
 
+lacy = []
 
 for i in range(0, 9):
-    Mdata.append(M[0, i])
-    MdataD.append(MD[0, i])
-    NdataD.append(ND[0, i])
-    OdataD.append(OD[0, i])
-    PdataD.append(PD[0, i])
+    lacy.append(data_lacy[0, i])
+
+Mdata = []
+Ndata = []
+Odata = []
+Pdata = []
+
+for idx, val in enumerate(age_daven):
+    if val == (1e6):
+        Mdata = data_daven[idx, :]
+    if val == (1e7):
+        Ndata = data_daven[idx, :]
+    if val == (1e8):
+        Odata = data_daven[idx, :]
+    if val == (1e9):
+        Pdata = data_daven[idx, :]
+
 
 # Plot
 fig = plt.figure(figsize=(4, 3))
 
 a = 0.5
 
 fp = [1, 1, 1, 1]
-np.interp(Mdata[5:9], Mdata[1:5], fp)
-plt.plot(Mdata[5:9], Mdata[1:5], "-", color="Orange", linewidth=a)
+np.interp(lacy[5:9], lacy[1:5], fp)
+plt.plot(lacy[5:9], lacy[1:5], "-", color="Orange", linewidth=a)
 
-np.interp(MdataD[5:9], MdataD[1:5], fp)
-plt.plot(MdataD[5:9], MdataD[1:5], "-", color="cornflowerblue", linewidth=a)
+np.interp(Mdata[5:9], Mdata[1:5], fp)
+plt.plot(Mdata[5:9], Mdata[1:5], "-", color="cornflowerblue", linewidth=a)
 
-np.interp(NdataD[5:9], NdataD[1:5], fp)
-plt.plot(NdataD[5:9], NdataD[1:5], "-", color="darkviolet", linewidth=a)
+np.interp(Ndata[5:9], Ndata[1:5], fp)
+plt.plot(Ndata[5:9], Ndata[1:5], "-", color="darkviolet", linewidth=a)
 
-np.interp(OdataD[5:9], OdataD[1:5], fp)
-plt.plot(OdataD[5:9], OdataD[1:5], "-", color="darkred", linewidth=a)
+np.interp(Odata[5:9], Odata[1:5], fp)
+plt.plot(Odata[5:9], Odata[1:5], "-", color="darkred", linewidth=a)
 
-np.interp(PdataD[5:9], PdataD[1:5], fp)
-plt.plot(PdataD[5:9], PdataD[1:5], "-", color="red", linewidth=a)
+np.interp(Pdata[5:9], Pdata[1:5], fp)
+plt.plot(Pdata[5:9], Pdata[1:5], "-", color="red", linewidth=a)
 
 
 for i in range(1, 5):
+    plt.plot(lacy[i + 4], lacy[i], label="t")
     plt.plot(Mdata[i + 4], Mdata[i], label="t")
-    plt.plot(MdataD[i + 4], MdataD[i], label="t")
-    plt.plot(NdataD[i + 4], NdataD[i], label="t")
-    plt.plot(OdataD[i + 4], OdataD[i], label="t")
-    plt.plot(PdataD[i + 4], PdataD[i], label="t")
+    plt.plot(Ndata[i + 4], Ndata[i], label="t")
+    plt.plot(Odata[i + 4], Odata[i], label="t")
+    plt.plot(Pdata[i + 4], Pdata[i], label="t")
 
 
 # Limits and scale
@@ -129,15 +130,12 @@
 
 
 # Label
-plt.title("Flare Frequency Distribution")  # , fontsize=37)
-plt.xlabel("log Flare Energy (erg)")  # ,fontsize=40)
-plt.ylabel("Cumulative Flare Freq (#/day)")  # ,fontsize=40)
-# plt.xticks(fontsize=40)
-# plt.yticks(fontsize=40)
-
+plt.title("Flare Frequency Distribution")
+plt.xlabel("log Flare Energy (erg)")
+plt.ylabel("Cumulative Flare Freq (#/day)")
 
 # Saving figure
 if sys.argv[1] == "pdf":
-    fig.savefig("FfdReproduced.pdf", bbox_inches="tight", dpi=600)
+    fig.savefig("FfdProxCen.pdf", bbox_inches="tight", dpi=600)
 if sys.argv[1] == "png":
-    fig.savefig("FfdReproduced.png", bbox_inches="tight", dpi=600)
+    fig.savefig("FfdProxCen.png", bbox_inches="tight", dpi=600)

src/flare.c

@@ -1729,16 +1729,12 @@ double fdLXUVFlare(BODY *body, double dDeltaTime, int iBody) {
   if (body[iBody].iFlareFFD == FLARE_FFD_DAVENPORT) {
     // The coefficient values given here were given by Dr. James Davenport in
     // private comunication
-    dFlareSlope = fdDavenport(-0.07054598,
-                              0.81225239,
-                              -1.07054511,
-                              body[iBody].dAge,
-                              body[iBody].dMass);
-    dFlareYInt  = fdDavenport(2.06012734,
-                              -25.79885288,
-                              34.44115635,
-                              body[iBody].dAge,
-                              body[iBody].dMass);
+    dFlareSlope =
+          fdDavenport(-0.07, 0.79, -1.06, body[iBody].dAge,
+                      body[iBody].dMass); //(-0.07054598,0.81225239,-1.07054511)
+    dFlareYInt = fdDavenport(
+          2.01, -25.15, 33.99, body[iBody].dAge,
+          body[iBody].dMass); //(2.06012734,-25.79885288,34.44115635)
   } else if (body[iBody].iFlareFFD == FLARE_FFD_LACY) {
     dFlareSlope = body[iBody].dFlareSlope;
     dFlareYInt  = body[iBody].dFlareYInt;