crarterlist.in now requires crater diameter input (interpolates for <1000 km)

MintonGroup · Apr 13, 2021 · 5fc97bf · 5fc97bf
1 parent 2857ee9
commit 5fc97bf
Show file tree

Hide file tree

Showing 5 changed files with 50 additions and 140 deletions.
diff --git a/.gitignore b/.gitignore
@@ -41,3 +41,5 @@ examples/global-lunar-bombardment/CTEM
 !LOLA*.dat
 !lunar*.dat
 !NPF*.dat
+
+examples/global-lunar-bombardment/scale.ipynb
diff --git a/examples/global-lunar-bombardment/craterlist.in b/examples/global-lunar-bombardment/craterlist.in
@@ -1,4 +1,3 @@
-#  imp   vel   ang   xoff  yoff  t_Ga
-90856.5  15.0e3   90.0  -5.17529e5  1.037286e6  0.02
-18095.0  15.0e3   90.0  -9.38629e5  1.3466084e6 0.03
-32432.1  15.0e3   90.0  -2.19298e5  1.24e6      0.025
+#  Dcrat(m)   vel   ang   xoff  yoff  t_Ga
+1028181  18.3e3   45.0  -5.17529e5  1.037286e6  0.02
+249840  18.3e3   45.0  -9.38629e5  1.3466084e6 0.03
diff --git a/examples/global-lunar-bombardment/ctem.in b/examples/global-lunar-bombardment/ctem.in
@@ -3,8 +3,8 @@
 
 ! Testing input. These are used to perform non-Monte Carlo tests.
 testflag     F                                 ! Set to T to create a single crater with user-defined impactor properties
-testimp      6.444e3 ! 100km crater            ! Diameter of test impactor (m)
-testvel      15.0e3                            ! Velocity of test crater (m/s)
+testimp      90856.5           ! Diameter of test impactor (m)
+testvel      18.3e3                            ! Velocity of test crater (m/s)
 testang      90.0                              ! Impact angle of test crater (deg) - Default 90.0
 testxoffset  0.0e0                             ! x-axis offset of crater center from grid center (m) - Default 0.0
 testyoffset  0.0e0                             ! y-axis offset of crater center from grid center (m) - Default 0.0

diff --git a/examples/global-lunar-bombardment/ctem_driver.py b/examples/global-lunar-bombardment/ctem_driver.py
@@ -13,6 +13,8 @@
 import os
 import subprocess
 import shutil
+import pandas
+from scipy.interpolate import interp1d
 
 #Import CTEM modules
 import ctem_io_readers
@@ -77,9 +79,22 @@
 
 #Read list of real craters and export to dat file for Fortran code
 if (parameters['quasimc'] == 'T'):
+
+    #Read list of real craters
     print("quasi-MC mode is ON")
     craterlistfile = parameters['workingdir'] + parameters['realcraterlist']
     rclist = ctem_io_readers.read_formatted_ascii(craterlistfile, skip_lines = 0)
+
+    #Interpolate craterscale.dat to get impactor sizes from crater sizes given
+    df = pandas.read_csv('craterscale.dat', sep='\s+')
+    df['log(Dc)'] = numpy.log(df['Dcrat(m)'])
+    df['log(Di)'] = numpy.log(df['#Dimp(m)'])
+    xnew = df['log(Dc)'].values
+    ynew = df['log(Di)'].values
+    interp = interp1d(xnew, ynew)
+    rclist[:,0] = numpy.exp(interp(numpy.log(rclist[:,0])))
+
+    #Convert age in Ga to "interval time"
     rclist[:,5] = (parameters['interval'] * parameters['numintervals']) - craterproduction.Tscale(rclist[:,5], 'NPF_Moon')
     rclist = rclist[rclist[:,5].argsort()]
     ctem_io_writers.write_realcraters(parameters, rclist)

diff --git a/examples/global-lunar-bombardment/scale.ipynb b/examples/global-lunar-bombardment/scale.ipynb
@@ -2,7 +2,7 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": 44,
+   "execution_count": 2,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -14,7 +14,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 45,
+   "execution_count": 3,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -23,7 +23,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 46,
+   "execution_count": 4,
    "metadata": {},
    "outputs": [
     {
@@ -129,7 +129,7 @@
        "[1418 rows x 2 columns]"
       ]
      },
-     "execution_count": 46,
+     "execution_count": 4,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -140,16 +140,16 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 47,
+   "execution_count": 5,
    "metadata": {},
    "outputs": [
     {
      "data": {
       "text/plain": [
-       "[<matplotlib.lines.Line2D at 0x7fc208b8a1c0>]"
+       "[<matplotlib.lines.Line2D at 0x7ff09b075e20>]"
       ]
      },
-     "execution_count": 47,
+     "execution_count": 5,
      "metadata": {},
      "output_type": "execute_result"
     },
@@ -172,16 +172,16 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 48,
+   "execution_count": 6,
    "metadata": {},
    "outputs": [
     {
      "data": {
       "text/plain": [
-       "[<matplotlib.lines.Line2D at 0x7fc207a2e610>]"
+       "[<matplotlib.lines.Line2D at 0x7ff09b23ee50>]"
       ]
      },
-     "execution_count": 48,
+     "execution_count": 6,
      "metadata": {},
      "output_type": "execute_result"
     },
@@ -207,164 +207,58 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 49,
+   "execution_count": 7,
    "metadata": {},
    "outputs": [],
    "source": [
-    "xnew = np.linspace(np.min(df['Dcrat(m)']), np.max(df['Dcrat(m)']))\n",
-    "ynew = np.linspace(np.min(df['#Dimp(m)']), np.max(df['#Dimp(m)']))\n",
-    "interp = interp1d(xnew, ynew)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 50,
-   "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "array(41608.16934796)"
-      ]
-     },
-     "execution_count": 50,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
-   "source": [
-    "interp(500000) #crater size as input; impactor size as output"
+    "df['log(Dc)'] = np.log(df['Dcrat(m)'])\n",
+    "df['log(Di)'] = np.log(df['#Dimp(m)'])"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 52,
-   "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "6.4817354382e-06"
-      ]
-     },
-     "execution_count": 52,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
-   "source": [
-    "np.min(df['Dcrat(m)'])"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 53,
-   "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "1028181.1762"
-      ]
-     },
-     "execution_count": 53,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
-   "source": [
-    "np.max(df['Dcrat(m)'])"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    " Generating a test crater\n",
-    " Dimp =    41608.169347960000     \n",
-    " Vimp =    15000.000000000000     \n",
-    " angle =    45.000000000000000     \n",
-    " x offset =    0.0000000000000000     \n",
-    " y offset =    0.0000000000000000     \n",
-    "\n",
-    "  0% (                                                  )  Dcrat =    415695.88033222407"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "^Not 500000 like the interp..."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 54,
+   "execution_count": 8,
    "metadata": {},
    "outputs": [],
    "source": [
-    "intrp2 = interp1d(ynew, xnew)"
+    "xnew = df['log(Dc)'].values\n",
+    "ynew = df['log(Di)'].values\n",
+    "interp = interp1d(xnew, ynew)"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 56,
+   "execution_count": 10,
    "metadata": {},
    "outputs": [
     {
      "data": {
       "text/plain": [
-       "array(499999.99999998)"
+       "35576.54212784875"
       ]
      },
-     "execution_count": 56,
+     "execution_count": 10,
      "metadata": {},
      "output_type": "execute_result"
     }
    ],
    "source": [
-    "intrp2(41608.16934796)"
+    "np.exp(interp(np.log(400000))) #crater size as input; impactor size as output"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
+    " Initializing simulation domain and determining minimum impactor size\n",
     " Generating a test crater\n",
-    " Dimp =    41608.169347960000     \n",
-    " Vimp =    15000.000000000000     \n",
-    " angle =    90.000000000000000     \n",
-    " x offset =    0.0000000000000000     \n",
-    " y offset =    0.0000000000000000     \n",
-    "\n",
-    "  0% (                                                  )  Dcrat =    511635.61214660661     \n",
-    "  \n",
-    "  \n",
-    "  ^Just making sure the angle isn't 90°"
+    " Dimp =    44200.2475301581     \n",
+    " Vimp =    18300.0000000000     \n",
+    " angle =    45.0000000000000     \n",
+    " x offset =   0.000000000000000E+000\n",
+    " y offset =   0.000000000000000E+000\n",
+    "  0% (                                                  )  Dcrat =    500013.835241806    "
    ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Generating a test crater\n",
-    " Dimp =    41608.169347960000     \n",
-    " Vimp =    18300.000000000000     \n",
-    " angle =    45.000000000000000     \n",
-    " x offset =    0.0000000000000000     \n",
-    " y offset =    0.0000000000000000     \n",
-    "\n",
-    "  0% (                                                  )  Dcrat =    468296.05484295724\n",
-    "  \n",
-    "  ^This uses the 'correct' velocity, yet still isn't 500 km like it should be..."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
   }
  ],
  "metadata": {