From 5fc97bf8dbf6a7168f35647aca6a0e4223602ee7 Mon Sep 17 00:00:00 2001 From: Austin Blevins Date: Mon, 12 Apr 2021 21:54:16 -0400 Subject: [PATCH] crarterlist.in now requires crater diameter input (interpolates for <1000 km) --- .gitignore | 2 + .../global-lunar-bombardment/craterlist.in | 7 +- examples/global-lunar-bombardment/ctem.in | 4 +- .../global-lunar-bombardment/ctem_driver.py | 15 ++ examples/global-lunar-bombardment/scale.ipynb | 162 +++--------------- 5 files changed, 50 insertions(+), 140 deletions(-) diff --git a/.gitignore b/.gitignore index f6ef0a31..98a5b65b 100644 --- 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 index 951d95dc..281a248a 100644 --- 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 \ No newline at end of file +# 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 \ No newline at end of file diff --git a/examples/global-lunar-bombardment/ctem.in b/examples/global-lunar-bombardment/ctem.in index 2c8f6c75..7ff9f52d 100755 --- 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 index 3f077046..c291c63f 100644 --- 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 index cd11a2c7..3a10a6e7 100644 --- 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": [ - "[]" + "[]" ] }, - "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": [ - "[]" + "[]" ] }, - "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": {