From f5636b972f823ebc828c31373bebad199e83ab9d Mon Sep 17 00:00:00 2001 From: David Minton Date: Wed, 9 Mar 2022 19:31:57 -0500 Subject: [PATCH] Removed obsolete files --- examples/quasimc-global/craterproduction.py | 216 -------------- examples/quasimc-global/ctem_driver.pro | 233 --------------- examples/quasimc-global/ctem_driver.py | 256 ---------------- examples/quasimc-global/ctem_image_dem.pro | 41 --- .../ctem_image_presentation.pro | 143 --------- .../quasimc-global/ctem_image_regolith.pro | 33 -- .../ctem_image_shaded_relief.pro | 52 ---- .../quasimc-global/ctem_io_read_input.pro | 131 -------- examples/quasimc-global/ctem_io_read_old.pro | 45 --- examples/quasimc-global/ctem_io_readers.py | 146 --------- examples/quasimc-global/ctem_io_writers.py | 282 ------------------ .../quasimc-global/ctem_window_display.pro | 249 ---------------- 12 files changed, 1827 deletions(-) delete mode 100644 examples/quasimc-global/craterproduction.py delete mode 100755 examples/quasimc-global/ctem_driver.pro delete mode 100755 examples/quasimc-global/ctem_driver.py delete mode 100755 examples/quasimc-global/ctem_image_dem.pro delete mode 100755 examples/quasimc-global/ctem_image_presentation.pro delete mode 100755 examples/quasimc-global/ctem_image_regolith.pro delete mode 100755 examples/quasimc-global/ctem_image_shaded_relief.pro delete mode 100755 examples/quasimc-global/ctem_io_read_input.pro delete mode 100755 examples/quasimc-global/ctem_io_read_old.pro delete mode 100644 examples/quasimc-global/ctem_io_readers.py delete mode 100644 examples/quasimc-global/ctem_io_writers.py delete mode 100755 examples/quasimc-global/ctem_window_display.pro diff --git a/examples/quasimc-global/craterproduction.py b/examples/quasimc-global/craterproduction.py deleted file mode 100644 index 4d81b687..00000000 --- a/examples/quasimc-global/craterproduction.py +++ /dev/null @@ -1,216 +0,0 @@ -import numpy as np -from scipy import optimize - -# The Neukum production function for the Moon and Mars -#Lunar PF from: Ivanov, Neukum, and Hartmann (2001) SSR v. 96 pp. 55-86 -#Mars PF from: Ivanov (2001) SSR v. 96 pp. 87-104 -sfd_coef = { - "NPF_Moon" : [-3.0876,-3.557528,+0.781027,+1.021521,-0.156012,-0.444058,+0.019977,+0.086850,-0.005874,-0.006809,+8.25e-4, +5.54e-5], - "NPF_Mars" : [-3.384, -3.197, +1.257, +0.7915, -0.4861, -0.3630, +0.1016, +6.756e-2,-1.181e-2,-4.753e-3,+6.233e-4,+5.805e-5] - } -sfd_range = { - "NPF_Moon" : [0.01,1000], - "NPF_Mars" : [0.015,362] -} -#Exponential time constant (Ga) -tau = 6.93 -Nexp = 5.44e-14 - -time_coef = { - "NPF_Moon" : Nexp, - "NPF_Mars" : Nexp * 10**(sfd_coef.get("NPF_Mars")[0]) / 10**(sfd_coef.get("NPF_Moon")[0]) -} - -def N1(T, pfmodel): - """ - Return the N(1) value as a function of time for a particular production function model - - Parameters - ---------- - T : numpy array - Time in units of Ga - pfmodel : string - The production function model to use - Currently options are 'NPF_Moon' or 'NPF_Mars' - - Returns - ------- - N1 : numpy array - The number of craters per square kilometer greater than 1 km in diameter - """ - T_valid_range = np.where((T <= 4.5) & (T >= 0.0), T, np.nan) - return time_coef.get(pfmodel) * (np.exp(tau * T_valid_range) - 1.0) + 10 ** (sfd_coef.get(pfmodel)[0]) * T_valid_range - -def CSFD(Dkm,planet): - """ - Return the cumulative size-frequency distribution for a particular production function model - - Parameters - ---------- - Dkm : numpy array - Diameters in units of km - pfmodel : string - The production function model to use - Currently options are 'NPF_Moon' or 'NPF_Mars' - - Returns - ------- - CSFD : numpy array - The number of craters per square kilometer greater than Dkm in diameter at T=1 Ga - """ - logCSFD = sum(co * np.log10(Dkm) ** i for i, co in enumerate(sfd_coef.get(planet))) - return 10**(logCSFD) - -def DSFD(Dkm,planet): - """ - Return the differential size-frequency distribution for a particular production function model - - Parameters - ---------- - Dkm : numpy array - Diameters in units of km - pfmodel : string - The production function model to use - Currently options are 'NPF_Moon' or 'NPF_Mars' - - Returns - ------- - DSFD : numpy array - The differential number of craters (dN/dD) per square kilometer greater than Dkm in diameter at T = 1 Ga - """ - dcoef = sfd_coef.get(planet)[1:] - logDSFD = sum(co * np.log10(Dkm) ** i for i, co in enumerate(dcoef)) - return 10**(logDSFD) * CSFD(Dkm,planet) / Dkm - -def Tscale(T,planet): - """ - Return the number density of craters at time T relative to time T = 1 Ga - - Parameters - ---------- - T : numpy array - Time in units of Ga - pfmodel : string - The production function model to use - Currently options are 'NPF_Moon' or 'NPF_Mars' - - Returns - ------- - Tscale : numpy array - N1(T) / CSFD(Dkm = 1.0) - """ - return N1(T,planet) / CSFD(1.0,planet) - -def pf_csfd(T,Dkm,planet): - """ - Return the cumulative size-frequency distribution for a particular production function model as a function of Time - - Parameters - ---------- - T : numpy array - Time in units of Ga - Dkm : numpy array - Diameters in units of km - pfmodel : string - The production function model to use - Currently options are 'NPF_Moon' or 'NPF_Mars' - - Returns - ------- - pf_csfd : numpy array - The cumulative number of craters per square kilometer greater than Dkm in diameter at time T T - """ - D_valid_range = np.where((Dkm >= sfd_range.get(planet)[0]) & (Dkm <= sfd_range.get(planet)[1]),Dkm,np.nan) - return CSFD(D_valid_range,planet) * Tscale(T,planet) - -def pf_dsfd(T,Dkm,planet): - """ - Return the differential size-frequency distribution for a particular production function model as a function of Time - - Parameters - ---------- - T : numpy array - Time in units of Ga - Dkm : numpy array - Diameters in units of km - pfmodel : string - The production function model to use - Currently options are 'NPF_Moon' or 'NPF_Mars' - - Returns - ------- - pf_dsfd : numpy array - The cumulative number of craters per square kilometer greater than Dkm in diameter at time T T - """ - D_valid_range = np.where((Dkm >= sfd_range.get(planet)[0]) & (Dkm <= sfd_range.get(planet)[1]), Dkm, np.nan) - return DSFD(D_valid_range, planet) * Tscale(T, planet) - -def T_from_scale(TS,planet): - """ - Return the time in Ga for the given number density of craters relative to that at 1 Ga. - This is the inverse of Tscale - - Parameters - ---------- - TS : numpy array - number density of craters relative to that at 1 Ga - pfmodel : string - The production function model to use - Currently options are 'NPF_Moon' or 'NPF_Mars' - - Returns - ------- - T_from_scale : numpy array - The time in Ga - """ - def func(S): - return Tscale(S,planet) - TS - return optimize.fsolve(func, np.full_like(TS,4.4),xtol=1e-10) - - -if __name__ == "__main__": - import matplotlib.pyplot as plt - import matplotlib.ticker as ticker - print("Tests go here") - print(f"T = 1 Ga, N(1) = {pf_csfd(1.0,1.00,'NPF_Moon')}") - print(f"T = 4.2 Ga, N(1) = {pf_csfd(4.2,1.00,'NPF_Moon')}") - print("Tscale test: Should return all 1s") - Ttest = np.logspace(-4,np.log10(4.4),num=100) - Tres = T_from_scale(Tscale(Ttest,'NPF_Mars'),'NPF_Mars') - print(Ttest / Tres) - #for i,t in enumerate(Ttest): - # print(t,Tscale(t,'Moon'),Tres[i]) - - CSFDfig = plt.figure(1, figsize=(8, 7)) - ax = {'NPF_Moon': CSFDfig.add_subplot(121), - 'NPF_Mars': CSFDfig.add_subplot(122)} - - tvals = [0.01,1.0,4.0] - x_min = 1e-3 - x_max = 1e3 - y_min = 1e-9 - y_max = 1e3 - Dvals = np.logspace(np.log10(x_min), np.log10(x_max)) - for key in ax: - ax[key].title.set_text(key) - ax[key].set_xscale('log') - ax[key].set_yscale('log') - ax[key].set_ylabel('$\mathregular{N_{>D} (km^{-2})}$') - ax[key].set_xlabel('Diameter (km)') - ax[key].set_xlim(x_min, x_max) - ax[key].set_ylim(y_min, y_max) - ax[key].yaxis.set_major_locator(ticker.LogLocator(base=10.0, numticks=20)) - ax[key].yaxis.set_minor_locator(ticker.LogLocator(base=10.0, subs=np.arange(2,10), numticks=100)) - ax[key].xaxis.set_major_locator(ticker.LogLocator(base=10.0, numticks=20)) - ax[key].xaxis.set_minor_locator(ticker.LogLocator(base=10.0, subs=np.arange(2,10), numticks=100)) - ax[key].grid(True,which="minor",ls="-",lw=0.5,zorder=5) - ax[key].grid(True,which="major",ls="-",lw=1,zorder=10) - for t in tvals: - prod = pf_csfd(t,Dvals,key) - ax[key].plot(Dvals, prod, '-', color='black', linewidth=1.0, zorder=50) - labeli = 15 - ax[key].text(Dvals[labeli],prod[labeli],f"{t:.2f} Ga", ha="left", va="top",rotation=-72) - - plt.tick_params(axis='y', which='minor') - plt.tight_layout() - plt.show() diff --git a/examples/quasimc-global/ctem_driver.pro b/examples/quasimc-global/ctem_driver.pro deleted file mode 100755 index 060d1b8d..00000000 --- a/examples/quasimc-global/ctem_driver.pro +++ /dev/null @@ -1,233 +0,0 @@ -pro ctem_driver - -;------------------------------------------------------------------------- -; Jim Richardson, Arecibo Observatory -; David Minton, Purdue University Dept. of Earth, Atmospheric, & Planetary Sciences -; July 2014 -; Cratered Terrain Evolution Model display module -; -; Inputs are read in from the ctem.in file -; -;------------------------------------------------------------------------- -;------------- Initial Setup ---------------- -;------------------------------------------------------------------------- -Compile_Opt DEFINT32 -!EXCEPT=2 - -; ----------- input file ----------------------- -infilename = 'ctem.in' -DATFILE='ctem.dat' - -; ---------- reading input files ---------- -seedarr = lon64arr(100) -seedn = 1 -totalimpacts = long64(0) -ncount = long64(0) -curyear = 0.d0 -restart = "F" -fracdone = 1.0d0 -masstot = 0.d0 - -ctem_io_read_input,infilename,interval,numintervals,gridsize,pix,seed,numlayers,sfdfile,impfile,maxcrat,ph1,shadedmaxhdefault,shadedminhdefault,shadedminh,shadedmaxh,restart,runtype,popupconsole,saveshaded,saverego,savepres,savetruelist - -seedarr(0) = seed -area = (gridsize * pix)^2 - -;read input data -pnum = file_lines(impfile) -production = dblarr(2,pnum) -productionfunction = dblarr(2,pnum) -openr,LUN,impfile,/GET_LUN -readf,LUN,productionfunction -free_lun,LUN - -;create impactor production population -production(0,*) = productionfunction(0,*) -production(1,*) = productionfunction(1,*)*area*interval - -;write out corrected production population -openw,1,sfdfile -printf,1,production -close,1 -free_lun,1 - -;set up cratering surface grid and display-only grid -surface_dem = dblarr(gridsize,gridsize) -regolith = dblarr(gridsize,gridsize) - -;set up temporary distribution bins -distl = 1 -pdistl = 1 -odist = dblarr(6,distl) -tdist = dblarr(6,distl) -pdist = dblarr(6,pdistl) -pdisttotal = dblarr(6,pdistl) - -datformat = "(I17,1X,I12,1X,E19.12,1X,A1,1X,F9.6,1X,E19.12)" - - -if strmatch(restart,'F',/fold_case) then begin ; Start with a clean slate - print, 'Starting a new run' - curyear = 0.0d0 - totalimpacts = 0 - masstot = 0.d0 - fracdone = 1.0d0 - - if strmatch(runtype,'statistical',/fold_case) then begin - ncount = 1 - openw,LUN,DATFILE,/GET_LUN - printf,LUN,totalimpacts,ncount,curyear,restart,fracdone,masstot,format=datformat - for n=0,seedn-1 do begin - printf,LUN,seedarr(n),format='(I12)' - endfor - free_lun,LUN - endif else begin - ncount = 0 - endelse - - surface_dem(*,*) = 0.0d0 - regolith(*,*) = 0.0d0 - - file_delete, 'tdistribution.dat',/allow_nonexistent - -endif else begin ; continue an old run - print, 'Continuing a previous run' - - ctem_io_read_old,gridsize,surface_dem,regolith,odist,tdist,pdist,mass - - ;read in constants file - openr,LUN,DATFILE,/GET_LUN - readf,LUN,totalimpacts,ncount,curyear,restart,fracdone,masstot,format=datformat - seedn = 0 - while ~ eof(LUN) do begin - readf,LUN,iseed - seedarr(seedn) = long64(iseed) - seedn=seedn+1 - endwhile - -endelse - -openw,FRACDONEFILE,'fracdone.dat',/GET_LUN -openw,REGODEPTHFILE,'regolithdepth.dat',/GET_LUN - -;------------------------------------------------------------------------- -; ---------- begin loops ---------- -;------------------------------------------------------------------------- -print, 'Beginning loops' - -;define number of loop iterations and begin -;numintervals=ceil(endyear/interval)-ncount -while (ncount le numintervals) do begin - - - ; ---------- creating crater population ---------- - if (ncount gt 0) then begin - - fnum = string(ncount,format='(I6.6)') - if (file_test('misc',/DIRECTORY) eq 0) then begin - file_mkdir,'misc' - endif - ; save a copy of the ctem.dat file - fname = 'misc/ctem_' + fnum + '.dat' - file_copy, 'ctem.dat', fname, /OVERWRITE - - print, ncount, ' Calling FORTRAN routine' - ;call fortran program to create & count craters - spawn, './CTEM',/noshell - - ; ---------- reading FORTRAN output ---------- - print, ncount, ' Reading FORTRAN output' - ctem_io_read_old,gridsize,surface_dem,regolith,odist,tdist,pdist,mass - - ;read in constants file - openr,LUN,DATFILE,/GET_LUN - readf,LUN,totalimpacts,ncount,curyear,restart,fracdone,masstot,format=datformat - seedn = 0 - while ~ eof(LUN) do begin - readf,LUN,iseed - seedarr(seedn) = long64(iseed) - seedn = seedn + 1 - endwhile - free_lun,LUN - - curyear = curyear + fracdone * interval - masstot = masstot + mass - printf,FRACDONEFILE,fracdone,curyear - flush,FRACDONEFILE - - printf,REGODEPTHFILE,curyear,mean(regolith),max(regolith),min(regolith) - flush,REGODEPTHFILE - - ;save a copy of the binned observed crater distribution - if (file_test('dist',/DIRECTORY) eq 0) then begin - file_mkdir,'dist' - endif - fname = 'dist/odist_' + fnum + '.dat' - file_copy, 'odistribution.dat', fname, /OVERWRITE - - ; save a copy of the cumulative observed crater distribution - fname = 'dist/ocum_' + fnum + '.dat' - file_copy, 'ocumulative.dat', fname, /OVERWRITE - - ;save a copy of the binned true distribution - fname = 'dist/tdist_' + fnum + '.dat' - file_copy, 'tdistribution.dat', fname, /OVERWRITE - - ;save a copy of the binned idealized production function - fname = 'dist/pdist_' + fnum + '.dat' - file_copy, 'pdistribution.dat', fname, /OVERWRITE - - ; save a copy of the cumulative true crater distribution if the user requests it - if strmatch(savetruelist,'T',/fold_case) then begin - fname = 'dist/tcum_' + fnum + '.dat' - file_copy, 'tcumulative.dat', fname, /OVERWRITE - endif - - ; save a copy of the impacted mass - fname = 'misc/mass_' + fnum + '.dat' - file_copy, 'impactmass.dat', fname, /OVERWRITE - - - endif - - ; Get the accumulated production function - pdisttotal = pdist - pdisttotal(3:5,*) = pdist(3:5,*) * curyear / interval - - ; ---------- displaying results ---------- - print, ncount, ' Displaying results' - - ctem_image_dem,ncount,gridsize,pix,surface_dem,surface_dem_image - if strmatch(saverego,'T',/fold_case) then ctem_image_regolith,ncount,gridsize,pix,regolith,regolith_image - if strmatch(saveshaded,'T',/fold_case) then begin - if (shadedminhdefault eq 1) then shadedminh = min(surface_dem) - if (shadedmaxhdefault eq 1) then shadedmaxh = max(surface_dem) - ctem_image_shaded_relief,ncount,gridsize,pix,surface_dem,surface_dem,shadedminh,shadedmaxh,'shaded',shaded_image - endif - if strmatch(savepres,'T',/fold_case) then ctem_image_presentation,ncount,gridsize,pix,curyear,odist,pdisttotal,tdist,ph1,surface_dem_image - ctem_window_display,ncount,totalimpacts,gridsize,pix,curyear,masstot,odist,pdisttotal,tdist,ph1,surface_dem,regolith,surface_dem_image,popupconsole - - ncount = ncount + 1 - - ;write out the current data file - if (strmatch(runtype,'statistical',/fold_case)) || (ncount eq 1) then begin - restart = 'F' - curyear = 0.0d0 - totalimpacts = 0 - masstot = 0.d0 - file_delete, 'tdistribution.dat',/allow_nonexistent - endif else begin - restart = 'T' - endelse - - openw,LUN,DATFILE,/GET_LUN - printf,LUN,totalimpacts,ncount,curyear,restart,fracdone,masstot,format=datformat - for n=0,seedn-1 do begin - printf,LUN,seedarr(n),format='(I12)' - endfor - free_lun,LUN -endwhile -free_lun,FRACDONEFILE -free_lun,REGODEPTHFILE - -end diff --git a/examples/quasimc-global/ctem_driver.py b/examples/quasimc-global/ctem_driver.py deleted file mode 100755 index 6fb56613..00000000 --- a/examples/quasimc-global/ctem_driver.py +++ /dev/null @@ -1,256 +0,0 @@ -#!/usr/bin/env python3 -# -#Cratered Terrain Evolution Model driver -# -#Original IDL design: Jim Richardson, Arecibo Observatory -#Revised IDL design: David Minton, Purdue University -#Re-engineered design and Python implementation: Matthew Route, Purdue University -#August 2016 - -#Import general purpose modules - -import numpy -import os -import subprocess -import shutil -import pandas -from scipy.interpolate import interp1d - -#Import CTEM modules -import ctem_io_readers -import ctem_io_writers -import craterproduction #craterproduction had to be cp'd to example dir - -#Create and initialize data dictionaries for parameters and options from CTEM.in -notset = '-NOTSET-' -currentdir = os.getcwd() + os.sep - -parameters={'restart': notset, - 'runtype': notset, - 'popupconsole': notset, - 'saveshaded': notset, - 'saverego': notset, - 'savepres': notset, - 'savetruelist': notset, - 'seedn': 1, - 'totalimpacts': 0, - 'ncount': 0, - 'curyear': 0.0, - 'fracdone': 1.0, - 'masstot': 0.0, - 'interval': 0.0, - 'numintervals': 0, - 'pix': -1.0, - 'gridsize': -1, - 'seed': 0, - 'maxcrat': 1.0, - 'shadedminhdefault': 1, - 'shadedmaxhdefault': 1, - 'shadedminh': 0.0, - 'shadedmaxh': 0.0, - 'workingdir': currentdir, - 'ctemfile': 'ctem.in', - 'datfile': 'ctem.dat', - 'impfile': notset, - 'sfdcompare': notset, - 'sfdfile': notset, - 'quasimc': notset, - 'realcraterlist': notset} - -#Read ctem.in to initialize parameter values based on user input -ctem_io_readers.read_ctemin(parameters,notset) - -#Read sfdcompare file -sfdfile = parameters['workingdir'] + parameters['sfdcompare'] -ph1 = ctem_io_readers.read_formatted_ascii(sfdfile, skip_lines = 0) - -#Set up data arrays -seedarr = numpy.zeros(100, dtype = numpy.int) -seedarr[0] = parameters['seed'] -odist = numpy.zeros([1, 6]) -pdist = numpy.zeros([1, 6]) -tdist = numpy.zeros([1, 6]) -surface_dem = numpy.zeros([parameters['gridsize'], parameters['gridsize']], dtype = numpy.float) -regolith = numpy.zeros([parameters['gridsize'], parameters['gridsize']], dtype =numpy.float) - -#Read production function file -impfile = parameters['workingdir'] + parameters['impfile'] -prodfunction = ctem_io_readers.read_formatted_ascii(impfile, skip_lines = 0) - -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, fill_value='extrapolate') - 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()] - - #Export to dat file for Fortran use - ctem_io_writers.write_realcraters(parameters, rclist) - -#Create impactor production population -area = (parameters['gridsize'] * parameters['pix'])**2 -production = numpy.copy(prodfunction) -production[:,1] = production[:,1] * area * parameters['interval'] - -#Write corrected production function to file -ctem_io_writers.write_production(parameters, production) - -#Starting new or old run? -if (parameters['restart'].upper() == 'F'): - print('Starting a new run') - - if (parameters['runtype'].upper() == 'STATISTICAL'): - parameters['ncount'] = 1 - - #Write ctem.dat file - ctem_io_writers.write_ctemdat(parameters, seedarr) - - else: - parameters['ncount'] = 0 - - #Delete tdistribution file, if it exists - tdist_file = parameters['workingdir'] + 'tdistribution.dat' - if os.path.isfile(tdist_file): - os.remove(tdist_file) - -else: - print('Continuing a previous run') - - #Read surface dem(shaded relief) and ejecta data files - dem_file = parameters['workingdir'] + 'surface_dem.dat' - surface_dem = ctem_io_readers.read_unformatted_binary(dem_file, parameters['gridsize']) - ejecta_file = parameters['workingdir'] + 'surface_ejc.dat' - regolith = ctem_io_readers.read_unformatted_binary(ejecta_file, parameters['gridsize']) - - #Read odistribution, tdistribution, and pdistribution files - ofile = parameters['workingdir'] + 'odistribution.dat' - odist = ctem_io_readers.read_formatted_ascii(ofile, skip_lines = 1) - tfile = parameters['workingdir'] + 'tdistribution.dat' - tdist = ctem_io_readers.read_formatted_ascii(tfile, skip_lines = 1) - pfile = parameters['workingdir'] + 'pdistribution.dat' - pdist = ctem_io_readers.read_formatted_ascii(pfile, skip_lines = 1) - - #Read impact mass from file - massfile = parameters['workingdir'] + 'impactmass.dat' - impact_mass = ctem_io_readers.read_impact_mass(massfile) - - #Read ctem.dat file - ctem_io_readers.read_ctemdat(parameters, seedarr) - -#Open fracdonefile and regodepthfile for writing -filename = parameters['workingdir'] + 'fracdone.dat' -fp_frac = open(filename,'w') -filename = parameters['workingdir'] + 'regolithdepth.dat' -fp_reg = open(filename,'w') - -#Begin CTEM processing loops -print('Beginning loops') - -ctem_io_writers.create_dir_structure(parameters) - -while (parameters['ncount'] <= parameters['numintervals']): - - #Create crater population - if (parameters['ncount'] > 0): - - #Move ctem.dat - forig = parameters['workingdir'] + 'ctem.dat' - fdest = parameters['workingdir'] + 'misc' + os.sep + "ctem%06d.dat" % parameters['ncount'] - shutil.copy2(forig, fdest) - - #Create crater population and display CTEM progress on screen - print(parameters['ncount'], ' Calling FORTRAN routine') - with subprocess.Popen([parameters['workingdir']+'CTEM'], stdout=subprocess.PIPE, bufsize=1,universal_newlines=True) as p: - for line in p.stdout: - print(line, end='') - - - #Optional: do not pipe CTEM progress to the screen - #subprocess.check_output([parameters['workingdir']+'CTEM']) - - #Read Fortran output - print(parameters['ncount'], ' Reading Fortran output') - - #Read surface dem(shaded relief) and ejecta data files - dem_file = parameters['workingdir'] + 'surface_dem.dat' - surface_dem = ctem_io_readers.read_unformatted_binary(dem_file, parameters['gridsize']) - ejecta_file = parameters['workingdir'] + 'surface_ejc.dat' - regolith = ctem_io_readers.read_unformatted_binary(ejecta_file, parameters['gridsize']) - - #Read odistribution, tdistribution, and pdistribution files - ofile = parameters['workingdir'] + 'odistribution.dat' - odist = ctem_io_readers.read_formatted_ascii(ofile, skip_lines = 1) - tfile = parameters['workingdir'] + 'tdistribution.dat' - tdist = ctem_io_readers.read_formatted_ascii(tfile, skip_lines = 1) - pfile = parameters['workingdir'] + 'pdistribution.dat' - pdist = ctem_io_readers.read_formatted_ascii(pfile, skip_lines = 1) - - #Read impact mass from file - massfile = parameters['workingdir'] + 'impactmass.dat' - impact_mass = ctem_io_readers.read_impact_mass(massfile) - - #Read ctem.dat file - ctem_io_readers.read_ctemdat(parameters, seedarr) - - #Update parameters: mass, curyear, regolith properties - parameters['masstot'] = parameters['masstot'] + impact_mass - - parameters['curyear'] = parameters['curyear'] + parameters['fracdone'] * parameters['interval'] - template = "%(fracdone)9.6f %(curyear)19.12E\n" - fp_frac.write(template % parameters) - - reg_text = "%19.12E %19.12E %19.12E %19.12E\n" % (parameters['curyear'], - numpy.mean(regolith), numpy.amax(regolith), numpy.amin(regolith)) - fp_reg.write(reg_text) - - #Save copy of crater distribution files - ctem_io_writers.copy_dists(parameters) - - #Display results - print(parameters['ncount'], ' Displaying results') - - #Write surface dem, surface ejecta, shaded relief, and rplot data - ctem_io_writers.image_dem(parameters, surface_dem) - if (parameters['saverego'].upper() == 'T'): - ctem_io_writers.image_regolith(parameters, regolith) - if (parameters['saveshaded'].upper() == 'T'): - ctem_io_writers.image_shaded_relief(parameters, surface_dem) - if (parameters['savepres'].upper() == 'T'): - ctem_io_writers.create_rplot(parameters,odist,pdist,tdist,ph1) - - #Update ncount - parameters['ncount'] = parameters['ncount'] + 1 - - if ((parameters['runtype'].upper() == 'STATISTICAL') or (parameters['ncount'] == 1)): - parameters['restart'] = 'F' - parameters['curyear'] = 0.0 - parameters['totalimpacts'] = 0 - parameters['masstot'] = 0.0 - - #Delete tdistribution file, if it exists - tdist_file = parameters['workingdir'] + 'tdistribution.dat' - if os.path.isfile(tdist_file): - os.remove(tdist_file) - - else: - parameters['restart'] = 'T' - - #Write ctem.dat file - ctem_io_writers.write_ctemdat(parameters, seedarr) - -#Close updateable fracdonefile and regodepthfile files -fp_frac.close() -fp_reg.close() diff --git a/examples/quasimc-global/ctem_image_dem.pro b/examples/quasimc-global/ctem_image_dem.pro deleted file mode 100755 index ad6fbd7e..00000000 --- a/examples/quasimc-global/ctem_image_dem.pro +++ /dev/null @@ -1,41 +0,0 @@ -pro ctem_image_dem,ncount,gridsize,pix,surface_dem,surface_dem_image -; Generates the shaded surface DEM image and saves it as a jpeg output image in the 'surf' directory -; outputs surface_dem_arr, which may be scaled to use as a console image - -; This code is for reading in the x-y positions from the cumulative distribution and separating out into layers -; so that the tallied craters can be drawn as circles -;!PATH = Expand_Path('+/home/campus/daminton/coyote/') + ':' + !PATH - -Compile_Opt DEFINT32 -thisDevice = !D.Name -Set_Plot, 'Z' -Erase -Device, Set_Resolution=[gridsize,gridsize],Set_Pixel_Depth=24, Decomposed=0 -TVLCT, red, green, blue, /GET - -sun = 20.d0 -radsun = sun * !dtor - -surface_dem_arr = dblarr(gridsize,gridsize) - -surface_dem_arr=surface_dem-shift(surface_dem,0,1) -surface_dem_arr = (0.5d0*!dpi) + atan(surface_dem_arr,pix) ; Get average slope -surface_dem_arr = abs(surface_dem_arr - radsun < 0.5d0*!dpi) ; incident angle -surface_dem_arr = 254.0d0 * cos(surface_dem_arr) ; shaded relief surface - -tv, surface_dem_arr, 0, 0, xsize=gridsize, ysize=gridsize, /device -print,'plotting' - -surface_dem_image = TVRD(True=1) -Set_Plot, thisDevice - -; save surface display -if (file_test('surf',/DIRECTORY) eq 0) then begin - file_mkdir,'surf' -endif -fnum = string(ncount,format='(I6.6)') -fname = 'surf/surf' + fnum + '.jpg' - -write_jpeg, fname, surface_dem_image, true=1, quality=100 - -end diff --git a/examples/quasimc-global/ctem_image_presentation.pro b/examples/quasimc-global/ctem_image_presentation.pro deleted file mode 100755 index 05a7d65a..00000000 --- a/examples/quasimc-global/ctem_image_presentation.pro +++ /dev/null @@ -1,143 +0,0 @@ -pro ctem_image_presentation,ncount,gridsize,pix,curyear,odist,pdist,tdist,ph1,map -; Generates the simplified version of the console display for use in presentations -; Plots the R-plot on the left and the image map on the right - -; presentation graphics -presresx = 1024 -presresy = 0.6*presresx -area = (gridsize * pix)^2 - -minx = (pix / 3.0d0) * 1d-3 -maxx = 3 * pix * gridsize * 1d-3 - -;strings for display -dum = string(format='(E10.4)', curyear) -parts = strsplit(dum, 'E', /extract) -fcuryear = strcompress(string(format='(F6.4,"x10!U",i,"!N")', parts[0], parts[1])) -time = 'Time = ' -timeunit = ' yr' - -if (file_test('presentation',/DIRECTORY) eq 0) then begin - file_mkdir,'presentation' -endif - -thisDevice = !D.Name -Set_Plot, 'Z' -Erase -Device, Set_Resolution=[presresx,presresy],Set_Pixel_Depth=24, Decomposed=0 -loadct, 0 -TVLCT, red, green, blue, /GET - -;geometric saturation -geom = dblarr(2,2) -geomem = dblarr(2,2) -geomep = dblarr(2,2) -geomel = dblarr(2,2) -geom(0,0) = minx -geom(0,1) = maxx -geom(1,*) = 3.12636d0 -geomem(0,0) = minx -geomem(0,1) = maxx -geomem(1,*) = 0.156318d0 -geomep(0,0) = minx -geomep(0,1) = maxx -geomep(1,*) = 0.312636d0 -geomel(0,0) = minx -geomel(0,1) = maxx -geomel(1,*) = 0.0312636d0 - -; Remove zeros -nz = where(odist(5,*) ne 0.0,count) -if (count gt 0) then begin - odistnz = dblarr(6,count) - odistnz(*,*)= odist(*,nz) -endif else begin - odistnz = odist -endelse - -nz = where(tdist(5,*) ne 0.0,count) -if (count gt 0) then begin - tdistnz = dblarr(6,count) - tdistnz(*,*)= tdist(*,nz) -endif else begin - tdistnz = tdist -endelse - -nz = where(pdist(5,*) ne 0.0,count) -if (count gt 0) then begin - pdistnz = dblarr(6,count) - pdistnz(*,*)= pdist(*,nz) -endif else begin - pdistnz = pdist -endelse - -; create r-plot array containing exactly 1 crater per bin -area = (gridsize * pix * 1d-3)^2 -plo = 1 -while (sqrt(2.0d0)^plo gt minx) do begin - plo = plo - 1 -endwhile -phi = plo + 1 -while (sqrt(2.0d0)^phi lt maxx) do begin - phi = phi + 1 -endwhile -n = phi - plo -sdist = dblarr(6,n + 1) -p = plo -for i=0,n do begin - sdist(0,i) = sqrt(2.d0)^p - sdist(1,i) = sqrt(2.d0)^(p+1) - sdist(2,i) = sqrt(sdist(0,i) * sdist(1,i)) - sdist(3,i) = 1.0d0 - sdist(5,i) = (sdist(2,i))^3 / (area * (sdist(1,i) - sdist(0,i))) - p = p + 1 -endfor -sdist(4,*) = reverse(total(sdist(3,*),/cumulative),2) - - -plot, odistnz(2,*)*1.0d-3, odistnz(5,*), line=0, color=255, thick=2.0, $ - TITLE='Crater Distribution R-Plot', charsize=1.2, $ - XTITLE='Crater Diameter (km)', $ - XRANGE=[minx,maxx], /XLOG, XSTYLE=1, $ - YTITLE='R Value', $ - YRANGE=[5.0e-4,5.0e0], /YLOG, YSTYLE=1, $ - /device, pos = [0.12,0.12,0.495,0.495]*presresx - -Dfac = sqrt(sqrt(2.0d0)) - -;display observed crater distribution -oplot, tdistnz(2,*)*1.0d-3, tdistnz(5,*), line=0, color=255, thick=1.0 -;oplot, geom(0,*), geom(1,*), line=2, color=255, thick=1.0 -oplot, geomem(0,*), geomem(1,*), line=1, color=255, thick=1.0 -oplot, geomep(0,*), geomep(1,*), line=1, color=255, thick=1.0 -;oplot, geomel(0,*), geomel(1,*), line=1, color=255, thick=1.0 -oplot, pdistnz(2,*)*1.0d-3, pdistnz(5,*), line=3, color=255, thick=1.0 -oplot, odistnz(2,*)*1.0d-3, odistnz(5,*), line=1, color=255, thick=1.0 -oplot, sdist(0,*), sdist(5,*), line=1, color=255, thick=1.0 -oplot, ph1(0,*)*Dfac*1.0d-3, ph1(1,*), psym=1, color=255, thick=1.0 -;oplot, ph2(0,*), ph2(1,*), psym=4, color=255, thick=1.0 -;oplot, ph3(0,*), ph3(1,*), psym=5, color=255, thick=1.0 - -;draw box around the main surface display & fill -displaysize = floor(0.45*presresx) ; size of displayed surface -surfxpos = floor(0.520*presresx) -surfypos = floor(0.12*presresy) -xbox = [surfxpos - 1,surfxpos - 1,surfxpos + displaysize,surfxpos + displaysize,surfxpos - 1] -ybox = [surfypos - 1,surfypos + displaysize,surfypos + displaysize,surfypos - 1,surfypos - 1] -plotS, xbox, ybox, /device, color=255 - -mapscaled = congrid(map,3,displaysize,displaysize,/interp) - -;mapscaled=congrid(map,displaysize,displaysize) ; scale image -tv, mapscaled, surfxpos, surfypos, xsize=displaysize, ysize=displaysize, true=1, /device - -xyouts, 0.310*presresx, 0.04*presresy, time + fcuryear + timeunit, color=255, charsize=2*presresy/720., /device - -snapshot = TVRD(True=1) -Set_Plot, thisDevice -fnum = string(ncount,format='(I6.6)') -fname = 'presentation/presentation' + fnum + '.jpg' -Write_JPEG, fname, snapshot, True=1, Quality=90 - - -end diff --git a/examples/quasimc-global/ctem_image_regolith.pro b/examples/quasimc-global/ctem_image_regolith.pro deleted file mode 100755 index 8d7f07a7..00000000 --- a/examples/quasimc-global/ctem_image_regolith.pro +++ /dev/null @@ -1,33 +0,0 @@ -pro ctem_image_regolith,ncount,gridsize,pix,regolith,regolith_image -; Generates the regolith depth map image and saves it as a jpeg output image in the 'rego' directory -; outputs dregolith, which may be scaled to use as a console image -Compile_Opt DEFINT32 -thisDevice = !D.Name -Set_Plot, 'Z' -Erase -Device, Set_Resolution=[gridsize,gridsize],Set_Pixel_Depth=24, Decomposed=0 -loadct, 39 -TVLCT, red, green, blue, /GET - -minref = pix * 1.0d-4 -regolith_scaled = dblarr(gridsize,gridsize) -maxreg = max(regolith) -minreg = min(regolith) -if minreg lt minref then minreg = minref -if maxreg lt minref then maxreg = minref + 1.0d30 -regolith_scaled = regolith > minreg -regolith_scaled = 254.0d0 * ((alog(regolith_scaled) - alog(minreg)) / (alog(maxreg) - alog(minreg))) - -; save regolith display -tv, regolith_scaled, 0, 0, xsize=gridsize, ysize=gridsize, /device -regolith_image = TVRD(True=1) -Set_Plot, thisDevice - -if (file_test('rego',/DIRECTORY) eq 0) then begin - file_mkdir,'rego' -endif -fnum = string(ncount,format='(I6.6)') -fname = 'rego/rego' + fnum + '.jpg' -write_jpeg, fname, regolith_image, true=1, quality=100 - -end diff --git a/examples/quasimc-global/ctem_image_shaded_relief.pro b/examples/quasimc-global/ctem_image_shaded_relief.pro deleted file mode 100755 index 73986f84..00000000 --- a/examples/quasimc-global/ctem_image_shaded_relief.pro +++ /dev/null @@ -1,52 +0,0 @@ -pro ctem_image_shaded_relief,ncount,gridsize,pix,surface_dem,height_vals,minh,maxh,dirname,shaded_image -; Generates the shaded depth map image and saves it as a jpeg output image in the 'rego' directory -; outputs dshaded, which may be scaled to use as a console image -; Uses the array height_vals for the color -Compile_Opt DEFINT32 -thisDevice = !D.Name -Set_Plot, 'Z' -Erase -Device, Set_Resolution=[gridsize,gridsize],Set_Pixel_Depth=24, Decomposed=0 -loadct, 33 -TVLCT, red, green, blue, /GET - -light=[[1,1,1],[0,0,0],[-1,-1,-1]] - -; convolution of the dem with the 3x3 matrix - -shaded=bytscl(convol(float(surface_dem),float(light))) -if max(shaded) eq 0 then shaded=255 - -; scale the dem to the height -if (minh gt maxh) then begin - tmp = minh - minh = maxh - maxh = tmp -endif - -demscaled = ((height_vals - minh) > 0.0) < (maxh-minh) -if ((maxh - minh) eq 0.0) then begin - demscaled = demscaled * 0.0 -endif else begin - demscaled = demscaled/(maxh-minh)*255.0 -endelse - -tv, demscaled, 0, 0, xsize=gridsize, ysize=gridsize, /device -shaded_image = TVRD(True=1) -shadedscl =float(shaded)/255.0 -shaded_imagearr = dblarr(3,gridsize,gridsize) -shaded_imagearr(0,*,*) = float(shaded_image(0,*,*)) * shadedscl(*,*) -shaded_imagearr(1,*,*) = float(shaded_image(1,*,*)) * shadedscl(*,*) -shaded_imagearr(2,*,*) = float(shaded_image(2,*,*)) * shadedscl(*,*) -shaded_image=round(shaded_imagearr) -Set_Plot, thisDevice - -if (file_test(dirname,/DIRECTORY) eq 0) then begin - file_mkdir,dirname -endif - -fnum = string(ncount,format='(I6.6)') -fname = dirname + '/' + dirname + fnum + '.jpg' -write_jpeg,fname,shaded_image,true=1,quality=90 - -end diff --git a/examples/quasimc-global/ctem_io_read_input.pro b/examples/quasimc-global/ctem_io_read_input.pro deleted file mode 100755 index ea135668..00000000 --- a/examples/quasimc-global/ctem_io_read_input.pro +++ /dev/null @@ -1,131 +0,0 @@ -pro ctem_io_read_input,infilename,interval,numintervals,gridsize,pix,seed,numlayers,sfdfile,impfile,maxcrat,ph1,shadedmaxhdefault,shadedminhdefault,shadedminh,shadedmaxh,restart,runtype,popupconsole,saveshaded,saverego,savepres,savetruelist -Compile_Opt DEFINT32 -print, 'Reading input file' -openr,infile,infilename, /GET_LUN -line="" -comment="!" -interval = 0.d0 -numintervals = 0 -pix=-1.0d0 -gridsize=-1 -seed = 0 -maxcrat = 1.0d0 -shadedmaxhdefault = 1 -shadedminhdefault = 1 -shadedminh = 0.d0 -shademaxnh = 0.d0 - - -; Set required strings to unset value -notset="-----NOTSET----" -sfdfile = notset -impfile = notset -sfdcompare = notset -restart = notset -runtype = notset -popupconsole = notset -saveshaded = notset -saverego = notset -savepres = notset -savetruelist = notset -while (not EOF(infile)) do begin - readf,infile,line - if (~strcmp(line,comment,1)) then begin - substrings = strsplit(line,' ',/extract) - if strmatch(substrings(0),'pix',/fold_case) then reads,substrings(1),pix - if strmatch(substrings(0),'gridsize',/fold_case) then reads,substrings(1),gridsize - if strmatch(substrings(0),'seed',/fold_case) then reads,substrings(1),seed - if strmatch(substrings(0),'sfdfile',/fold_case) then reads,substrings(1),sfdfile - if strmatch(substrings(0),'impfile',/fold_case) then reads,substrings(1),impfile - if strmatch(substrings(0),'maxcrat',/fold_case) then reads,substrings(1),maxcrat - if strmatch(substrings(0),'sfdcompare',/fold_case) then reads,substrings(1),sfdcompare - if strmatch(substrings(0),'interval',/fold_case) then reads,substrings(1),interval - if strmatch(substrings(0),'numintervals',/fold_case) then reads,substrings(1),numintervals - if strmatch(substrings(0),'popupconsole',/fold_case) then reads,substrings(1),popupconsole - if strmatch(substrings(0),'saveshaded',/fold_case) then reads,substrings(1),saveshaded - if strmatch(substrings(0),'saverego',/fold_case) then reads,substrings(1),saverego - if strmatch(substrings(0),'savepres',/fold_case) then reads,substrings(1),savepres - if strmatch(substrings(0),'savetruelist',/fold_case) then reads,substrings(1),savetruelist - if strmatch(substrings(0),'runtype',/fold_case) then reads,substrings(1),runtype - if strmatch(substrings(0),'restart',/fold_case) then reads,substrings(1),restart - if strmatch(substrings(0),'shadedminh',/fold_case) then begin - reads,substrings(1),shadedminh - shadedminhdefault = 0 - endif - if strmatch(substrings(0),'shadedmaxh',/fold_case) then begin - reads,substrings(1),shadedmaxh - shadedmaxhdefault = 0 - endif - end -end -if interval le 0.0d0 then begin - print,'Invalid value for or missing variable INTERVAL in ' + infilename - stop -end -if numintervals le 0 then begin - print,'Invalid value for or missing variable NUMINTERVALS in ' + infilename - stop -end -if pix le 0.0d0 then begin - print,'Invalid value for or missing variable PIX in ' + infilename - stop -end -if gridsize le 0 then begin - print,'Invalid value for or missing variable GRIDSIZE in ' + infilename - stop -end -if seed eq 0 then begin - print,'Invalid value for or missing variable SEED in ' + infilename - stop -end -if strmatch(sfdfile,notset) then begin - print,'Invalid value for or missing variable SFDFILE in ' + infilename - stop -end -if strmatch(impfile,notset) then begin - print,'Invalid value for or missing variable IMPFILE in ' + infilename - stop -end -if strmatch(popupconsole,notset) then begin - print,'Invalid value for or missing variable POPUPCONSOLE in ' + infilename - stop -end -if strmatch(saveshaded,notset) then begin - print,'Invalid value for or missing variable SAVESHADED in ' + infilename - stop -end -if strmatch(saverego,notset) then begin - print,'Invalid value for or missing variable SAVEREGO in ' + infilename - stop -end -if strmatch(savepres,notset) then begin - print,'Invalid value for or missing variable SAVEPRES in ' + infilename - stop -end -if strmatch(savetruelist,notset) then begin - print,'Invalid value for or missing variable SAVETRUELIST in ' + infilename - stop -end -if strmatch(runtype,notset) then begin - print,'Invalid value for or missing variable RUNTYPE in ' + infilename - stop -end -if strmatch(restart,notset) then begin - print,'Invalid value for or missing variable RESTART in ' + infilename - stop -end - -free_lun,infile - -ph1 = dblarr(3,1) -if ~strmatch(sfdcompare,notset) then begin - cnum = file_lines(sfdcompare) - ph1 = dblarr(3,cnum) - openr,COMP,sfdcompare,/GET_LUN - readf,COMP,ph1 - close,COMP - free_lun,COMP -end -free_lun,infile - -end diff --git a/examples/quasimc-global/ctem_io_read_old.pro b/examples/quasimc-global/ctem_io_read_old.pro deleted file mode 100755 index a9013687..00000000 --- a/examples/quasimc-global/ctem_io_read_old.pro +++ /dev/null @@ -1,45 +0,0 @@ -pro ctem_io_read_old,gridsize,surface_dem,regolith,odist,tdist,pdist,mass -Compile_Opt DEFINT32 -openr,LUN,'surface_ejc.dat',/GET_LUN -readu,LUN,regolith -free_lun,LUN - -openr,LUN,'surface_dem.dat',/GET_LUN -readu,LUN,surface_dem -free_lun,LUN - -distl = file_lines('odistribution.dat') - 1 -pdistl = file_lines('pdistribution.dat') - 1 - -;read in observed cumulative distribution -odist = dblarr(6,distl) -line = "temp" -openr,LUN,'odistribution.dat',/GET_LUN -readf,LUN,line -readf,LUN,odist -close,LUN -free_lun,LUN - -;read in true cumulative distribution -tdist = dblarr(6,distl) -openr,LUN,'tdistribution.dat',/GET_LUN -readf,LUN,line -readf,LUN,tdist -close,LUN -free_lun,LUN - -;read in production function -pdist = dblarr(6,pdistl) -openr,LUN,'pdistribution.dat',/GET_LUN -readf,LUN,line -readf,LUN,pdist -close,LUN -free_lun,LUN - -;read in accumulated mass from file -openr,LUN,'impactmass.dat',/GET_LUN -readf,LUN,mass -close,LUN -free_lun,LUN - -end diff --git a/examples/quasimc-global/ctem_io_readers.py b/examples/quasimc-global/ctem_io_readers.py deleted file mode 100644 index 0881ad1a..00000000 --- a/examples/quasimc-global/ctem_io_readers.py +++ /dev/null @@ -1,146 +0,0 @@ -#!/usr/local/bin/python -# -#Cratered Terrain Evolution Model file reading utilities -# -#Original IDL design: Jim Richardson, Arecibo Observatory -#Revised IDL design: David Minton, Purdue University -#Re-engineered design and Python implementation: Matthew Route, Purdue University -#August 2016 -# -#Known issues for operation with CTEM -#1) ctem.in has 1.0d0 value for maxcrat which is not readable by Python - -import numpy - -def read_ctemin(parameters,notset): - #Read and parse ctem.in file - inputfile = parameters['workingdir'] + parameters['ctemfile'] - - #Read ctem.in file - print('Reading input file '+ parameters['ctemfile']) - fp = open(inputfile,'r') - lines = fp.readlines() - fp.close() - - #Process file text - for line in lines: - fields = line.split() - if len(fields) > 0: - if ('pix' == fields[0].lower()): parameters['pix']=float(fields[1]) - if ('gridsize' == fields[0].lower()): parameters['gridsize']=int(fields[1]) - if ('seed' == fields[0].lower()): parameters['seed']=int(fields[1]) - if ('sfdfile' == fields[0].lower()): parameters['sfdfile']=fields[1] - if ('impfile' == fields[0].lower()): parameters['impfile']=fields[1] - if ('maxcrat' == fields[0].lower()): parameters['maxcrat']=float(fields[1]) - if ('sfdcompare' == fields[0].lower()): parameters['sfdcompare']=fields[1] - if ('interval' == fields[0].lower()): parameters['interval']=float(fields[1]) - if ('numintervals' == fields[0].lower()): parameters['numintervals']=int(fields[1]) - if ('popupconsole' == fields[0].lower()): parameters['popupconsole']=fields[1] - if ('saveshaded' == fields[0].lower()): parameters['saveshaded']=fields[1] - if ('saverego' == fields[0].lower()): parameters['saverego']=fields[1] - if ('savepres' == fields[0].lower()): parameters['savepres']=fields[1] - if ('savetruelist' == fields[0].lower()): parameters['savetruelist']=fields[1] - if ('runtype' == fields[0].lower()): parameters['runtype']=fields[1] - if ('restart' == fields[0].lower()): parameters['restart']=fields[1] - if ('quasimc' == fields[0].lower()): parameters['quasimc']=fields[1] - if ('realcraterlist' == fields[0].lower()): parameters['realcraterlist']=fields[1] - if ('shadedminh' == fields[0].lower()): - parameters['shadedminh'] = float(fields[1]) - parameters['shadedminhdefault'] = 0 - if ('shadedmaxh' == fields[0].lower()): - parameters['shadedmaxh'] = float(fields[1]) - parameters['shadedmaxhdefault'] = 0 - - #Test values for further processing - if (parameters['interval'] <= 0.0): - print('Invalid value for or missing variable INTERVAL in '+ inputfile) - if (parameters['numintervals'] <= 0): - print('Invalid value for or missing variable NUMINTERVALS in '+ inputfile) - if (parameters['pix'] <= 0.0): - print('Invalid value for or missing variable PIX in '+ inputfile) - if (parameters['gridsize'] <= 0): - print('Invalid value for or missing variable GRIDSIZE in '+ inputfile) - if (parameters['seed'] == 0): - print('Invalid value for or missing variable SEED in '+ inputfile) - if (parameters['sfdfile'] == notset): - print('Invalid value for or missing variable SFDFILE in '+ inputfile) - if (parameters['impfile'] == notset): - print('Invalid value for or missing variable IMPFILE in '+ inputfile) - if (parameters['popupconsole'] == notset): - print('Invalid value for or missing variable POPUPCONSOLE in '+ inputfile) - if (parameters['saveshaded'] == notset): - print('Invalid value for or missing variable SAVESHADED in '+ inputfile) - if (parameters['saverego'] == notset): - print('Invalid value for or missing variable SAVEREGO in '+ inputfile) - if (parameters['savepres'] == notset): - print('Invalid value for or missing variable SAVEPRES in '+ inputfile) - if (parameters['savetruelist'] == notset): - print('Invalid value for or missing variable SAVETRUELIST in '+ inputfile) - if (parameters['runtype'] == notset): - print('Invalid value for or missing variable RUNTYPE in '+ inputfile) - if (parameters['restart'] == notset): - print('Invalid value for or missing variable RESTART in '+ inputfile) - - return - -def read_formatted_ascii(filename, skip_lines): - #Generalized ascii text reader - #For use with sfdcompare, production, odist, tdist, pdist data files - data = numpy.genfromtxt(filename, skip_header = skip_lines) - return data - -def read_unformatted_binary(filename, gridsize): - #Read unformatted binary files created by Fortran - #For use with surface ejecta and surface dem data files - dt = numpy.float - data = numpy.fromfile(filename, dtype = dt) - data.shape = (gridsize,gridsize) - - return data - -def read_ctemdat(parameters, seedarr): - #Read and parse ctem.dat file - datfile = parameters['workingdir'] + 'ctem.dat' - - #Read ctem.dat file - print('Reading input file '+ parameters['datfile']) - fp = open(datfile,'r') - lines = fp.readlines() - fp.close() - - #Parse file lines and update parameter fields - fields = lines[0].split() - if len(fields) > 0: - parameters['totalimpacts'] = float(fields[0]) - parameters['ncount'] = int(fields[1]) - parameters['curyear'] = float(fields[2]) - parameters['restart'] = fields[3] - parameters['fracdone'] = float(fields[4]) - parameters['masstot'] = float(fields[5]) - - #Parse remainder of file to build seed array - nlines = len(lines) - index = 1 - while (index < nlines): - fields = lines[index].split() - seedarr[index - 1] = float(fields[0]) - index += 1 - - parameters['seedn'] = index - 1 - - return - -def read_impact_mass(filename): - #Read impact mass file - - fp=open(filename,'r') - line=fp.readlines() - fp.close() - - fields = line[0].split() - if (len(fields) > 0): - mass = float(fields[0]) - else: - mass = 0 - - return mass \ No newline at end of file diff --git a/examples/quasimc-global/ctem_io_writers.py b/examples/quasimc-global/ctem_io_writers.py deleted file mode 100644 index a5e72b26..00000000 --- a/examples/quasimc-global/ctem_io_writers.py +++ /dev/null @@ -1,282 +0,0 @@ -#!/usr/local/bin/python -# -#Cratered Terrain Evolution Model file and graphical writing utilities -# -#Original IDL design: Jim Richardson, Arecibo Observatory -#Revised IDL design: David Minton, Purdue University -#Re-engineered design and Python implementation: Matthew Route, Purdue University -#August 2016 -# - -import matplotlib -from matplotlib import pyplot -import numpy -import os -import shutil -import scipy -from scipy import signal - -#Set pixel scaling common for image writing, at 1 pixel/ array element -dpi = 72.0 - -#Write production function to file production.dat -#This file format does not exactly match that generated from IDL. Does it work? -def write_production(parameters, production): - filename = parameters['workingdir'] + parameters['sfdfile'] - numpy.savetxt(filename, production, fmt='%1.8e', delimiter=' ') - - return - -def create_dir_structure(parameters): - #Create directories for various output files if they do not already exist - directories=['dist','misc','rego','rplot','surf','shaded'] - - for directory in directories: - dir_test = parameters['workingdir'] + directory - if not os.path.isdir(dir_test): - os.makedirs(dir_test) - - return - -def write_ctemdat(parameters, seedarr): - #Write various parameters and random number seeds into ctem.dat file - filename = parameters['workingdir'] + parameters['datfile'] - fp = open(filename,'w') - - template = "%(totalimpacts)17d %(ncount)12d %(curyear)19.12E %(restart)s %(fracdone)9.6f %(masstot)19.12E\n" - fp.write(template % parameters) - - #Write random number seeds to the file - for index in range(parameters['seedn']): - fp.write("%12d\n" % seedarr[index]) - - fp.close() - - return - -def copy_dists(parameters): - #Save copies of distribution files - - orig_list = ['odistribution', 'ocumulative', 'pdistribution', 'tdistribution'] - dest_list = ['odist', 'ocum', 'pdist', 'tdist'] - - for index in range(len(orig_list)): - forig = parameters['workingdir'] + orig_list[index] + '.dat' - fdest = parameters['workingdir'] + 'dist' + os.sep + dest_list[index] + "_%06d.dat" % parameters['ncount'] - shutil.copy2(forig, fdest) - - forig = parameters['workingdir'] + 'impactmass.dat' - fdest = parameters['workingdir'] + 'misc' + os.sep + "mass_%06d.dat" % parameters['ncount'] - shutil.copy2(forig, fdest) - - if (parameters['savetruelist'].upper() == 'T'): - forig = parameters['workingdir'] + 'tcumulative.dat' - fdest = parameters['workingdir'] + 'dist' + os.sep + "tcum_%06d.dat" % parameters['ncount'] - shutil.copy2(forig, fdest) - - return - -#Possible references -#http://nbviewer.jupyter.org/github/ThomasLecocq/geophysique.be/blob/master/2014-02-25%20Shaded%20Relief%20Map%20in%20Python.ipynb - -def image_dem(parameters, surface_dem): - - #Create surface dem map - solar_angle = 20.0 - dem_map = numpy.copy(surface_dem) - numpy.roll(surface_dem, 1, 0) - dem_map = (0.5 * numpy.pi) + numpy.arctan2(dem_map, parameters['pix']) - dem_map = dem_map - numpy.radians(solar_angle) * (0.5 * numpy.pi) - numpy.place(dem_map, dem_map > (0.5 * numpy.pi), 0.5 *numpy.pi) - dem_map = numpy.absolute(dem_map) - dem_map = 254.0 * numpy.cos(dem_map) - - #Save image to file - filename = parameters['workingdir'] + 'surf' + os.sep + "surf%06d.png" % parameters['ncount'] - height = parameters['gridsize'] / dpi - width = height - fig = matplotlib.pyplot.figure(figsize = (width, height), dpi = dpi) - fig.figimage(dem_map, cmap = matplotlib.cm.gray, origin = 'lower') - matplotlib.pyplot.savefig(filename) - - return - -def image_regolith(parameters, regolith): - - #Create scaled regolith image - minref = parameters['pix'] * 1.0e-4 - maxreg = numpy.amax(regolith) - minreg = numpy.amin(regolith) - if (minreg < minref): minreg = minref - if (maxreg < minref): maxreg = (minref + 1.0e3) - regolith_scaled = numpy.copy(regolith) - numpy.place(regolith_scaled, regolith_scaled < minref, minref) - regolith_scaled = 254.0 * ((numpy.log(regolith_scaled) - numpy.log(minreg)) / (numpy.log(maxreg) - numpy.log(minreg))) - - #Save image to file - filename = parameters['workingdir'] + 'rego' + os.sep + "rego%06d.png" % parameters['ncount'] - height = parameters['gridsize'] / dpi - width = height - fig = matplotlib.pyplot.figure(figsize = (width, height), dpi = dpi) - fig.figimage(regolith_scaled, cmap = matplotlib.cm.nipy_spectral, origin = 'lower') - matplotlib.pyplot.savefig(filename) - - return - -def image_shaded_relief(parameters, surface_dem): - #The color scale and appearance of this do not quite match the IDL version - - #Create image by convolving DEM with 3x3 illumination matrix - light = numpy.array([[1.0, 1.0, 1.0], [0.0, 0.0, 0.0], [-1.0, -1.0, -1.0]]) - convolved_map = scipy.signal.convolve2d(surface_dem, light, mode = 'same') - - #Adjust output to resemble IDL (north slopes illuminated, south in shadow) - convolved_map = convolved_map * -1.0 - convolved_map[0,:]=0.0 - convolved_map[-1,:] =0.0 - convolved_map[:,0]=0.0 - convolved_map[:,-1]=0.0 - - #If no shadedmin/max parameters are read in from ctem.dat, determine the values from the data - if (parameters['shadedminhdefault'] == 1): shadedminh = numpy.amin(surface_dem) - if (parameters['shadedmaxhdefault'] == 1): shadedmaxh = numpy.amax(surface_dem) - - #If min and max appear to be reversed, then fix them - if (parameters['shadedminh'] > parameters['shadedmaxh']): - temp = parameters['shadedminh'] - shadedminh = parameters['shadedmaxh'] - shadedmaxh = temp - else: - shadedminh = parameters['shadedminh'] - shadedmaxh = parameters['shadedmaxh'] - - #If dynamic range is valid, construct a shaded DEM - dynamic_range = shadedmaxh - shadedminh - if (dynamic_range != 0): - dem_scaled = numpy.copy(surface_dem) - shadedminh - numpy.place(dem_scaled, dem_scaled < 0.0, 0.0) - numpy.place(dem_scaled, dem_scaled > dynamic_range, dynamic_range) - dem_scaled = dem_scaled / dynamic_range - else: - dem_scaled = numpy.copy(surface_dem) * 0.0 - - #Generate shaded depth map with surface_dem color scaling (RGBA) - shaded = scipy.misc.bytescale(convolved_map) - if numpy.amax(shaded) == 0: shaded=255 - shadedscl = shaded / 255.0 - - #shaded_imagearr = matplotlib.cm.jet(dem_scaled) - #print dem_scaled[0:4,0:4] - #shaded_imagearr[:,:,0] = shaded_imagearr[:,:,0] * shadedscl - #shaded_imagearr[:,:,1] = shaded_imagearr[:,:,1] * shadedscl - #shaded_imagearr[:,:,2] = shaded_imagearr[:,:,2] * shadedscl - #shaded_imagearr[:,:,3] = shaded_imagearr[:,:,3] * shadedscl - - #Delivers nearly proper coloring, but no shaded relief - shaded_imagearr = dem_scaled * shadedscl - shaded_imagearr = matplotlib.cm.jet(shaded_imagearr) - shaded_imagearr = numpy.around(shaded_imagearr, decimals = 1) - - #Save image to file - filename = parameters['workingdir'] + 'shaded' + os.sep + "shaded%06d.png" % parameters['ncount'] - height = parameters['gridsize'] / dpi - width = height - fig = matplotlib.pyplot.figure(figsize = (width, height), dpi = dpi) - fig.figimage(shaded_imagearr, cmap = matplotlib.cm.jet, origin = 'lower') - matplotlib.pyplot.savefig(filename) - return - -def create_rplot(parameters,odist,pdist,tdist,ph1): - #Parameters: empirical saturation limit and dfrac - satlimit = 3.12636 - dfrac = 2**(1./4) * 1.0e-3 - - #Calculate geometric saturation - minx = (parameters['pix'] / 3.0) * 1.0e-3 - maxx = 3 * parameters['pix'] * parameters['gridsize'] * 1.0e-3 - geomem = numpy.array([[minx, satlimit / 20.0], [maxx, satlimit / 20.0]]) - geomep = numpy.array([[minx, satlimit / 10.0], [maxx, satlimit / 10.0]]) - - #Create distribution arrays without zeros for plotting on log scale - idx = numpy.nonzero(odist[:,5]) - odistnz = odist[idx] - odistnz = odistnz[:,[2,5]] - - idx = numpy.nonzero(pdist[:,5]) - pdistnz = pdist[idx] - pdistnz = pdistnz[:,[2,5]] - - idx = numpy.nonzero(tdist[:,5]) - tdistnz = tdist[idx] - tdistnz = tdistnz[:,[2,5]] - - #Correct pdist - pdistnz[:,1] = pdistnz[:,1] * parameters['curyear'] / parameters['interval'] - - #Create sdist bin factors, which contain one crater per bin - area = (parameters['gridsize'] * parameters['pix'] * 1.0e-3)**2. - plo = 1 - sq2 = 2**(1./2) - while (sq2**plo > minx): - plo = plo - 1 - phi = plo + 1 - while (sq2**phi < maxx): - phi = phi + 1 - n = phi - plo + 1 - sdist = numpy.zeros([n , 2]) - p = plo - for index in range(n): - sdist[index, 0] = sq2**p - sdist[index, 1] = sq2**(2.0*p + 1.5) / (area * (sq2 - 1)) - p = p + 1 - - #Create time label - tlabel = "%5.4e" % parameters['curyear'] - tlabel = tlabel.split('e') - texp = str(int(tlabel[1])) - timelabel = 'Time = '+ r'${}$ x 10$^{}$'.format(tlabel[0], texp) + ' yrs' - - #Save image to file - filename = parameters['workingdir'] + 'rplot' + os.sep + "rplot%06d.png" % parameters['ncount'] - height = parameters['gridsize'] / dpi - width = height - fig = matplotlib.pyplot.figure(figsize = (width, height), dpi = dpi) - - #Alter background color to be black, and change axis colors accordingly - matplotlib.pyplot.style.use('dark_background') - matplotlib.pyplot.rcParams['axes.prop_cycle'] - - #Plot data - matplotlib.pyplot.plot(odistnz[:,0]*1.0e-3, odistnz[:,1], linewidth=3.0, color = 'blue') - matplotlib.pyplot.plot(pdistnz[:,0]*1.0e-3, pdistnz[:,1], linewidth=2.0, linestyle='dashdot', color = 'white') - matplotlib.pyplot.plot(tdistnz[:,0]*1.0e-3, tdistnz[:,1], linewidth=2.0, color = 'red') - - matplotlib.pyplot.plot(geomem[:,0], geomem[:,1], linewidth=2.0, linestyle =':', color = 'yellow') - matplotlib.pyplot.plot(geomep[:,0], geomep[:,1], linewidth=2.0, linestyle =':', color = 'yellow') - matplotlib.pyplot.plot(sdist[:,0], sdist[:,1], linewidth=2.0, linestyle =':', color = 'yellow') - - matplotlib.pyplot.plot(ph1[:,0] * dfrac, ph1[:,1], 'wo') - - #Create plot labels - matplotlib.pyplot.title('Crater Distribution R-Plot',fontsize=22) - matplotlib.pyplot.xlim([minx, maxx]) - matplotlib.pyplot.xscale('log') - matplotlib.pyplot.xlabel('Crater Diameter (km)',fontsize=18) - matplotlib.pyplot.ylim([5.0e-4, 5.0]) - matplotlib.pyplot.yscale('log') - matplotlib.pyplot.ylabel('R Value', fontsize=18) - matplotlib.pyplot.text(1.0e-2, 1.0, timelabel, fontsize=18) - - matplotlib.pyplot.tick_params(axis='both', which='major', labelsize=14) - matplotlib.pyplot.tick_params(axis='both', which='minor', labelsize=12) - - matplotlib.pyplot.savefig(filename) - - return - -def write_realcraters(parameters, realcraters): - """writes file of real craters for use in quasi-MC runs""" - - filename = parameters['workingdir'] + 'craterlist.dat' - numpy.savetxt(filename, realcraters, fmt='%1.8e', delimiter='\t') - - return \ No newline at end of file diff --git a/examples/quasimc-global/ctem_window_display.pro b/examples/quasimc-global/ctem_window_display.pro deleted file mode 100755 index d6edf5b9..00000000 --- a/examples/quasimc-global/ctem_window_display.pro +++ /dev/null @@ -1,249 +0,0 @@ -pro ctem_window_display,ncount,totalimpacts,gridsize,pix,curyear,masstot,odist,pdist,tdist,ph1,surface_dem,regolith,map,popupconsole -; Produces the console window display. The array 'map' is used to generate the image on the right-hand side -Compile_Opt DEFINT32 - -; console output graphics resolution -minref = pix * 1.0d-4 -conresx = 1280 -;conresy = 0.758*conresx -profsize = 0.00 -conresy = (0.60 + profsize) * conresx -sun = 20.d0 ; sun angle for display - -minx = (pix / 3.0d0) * 1d-3 -maxx = 3 * pix * gridsize * 1d-3 -charfac = 1.6 ; character size multiplier - -if strmatch(popupconsole,'T',/fold_case) then begin - device, decomposed=0, retain=2 - thisDevice = !D.NAME - SET_PLOT, 'Z' - TVLCT, red, green, blue, /GET - SET_PLOT, thisDevice - !P.MULTI = [0, 1, 2] - Window, 0, xsize=conresx, ysize=conresy -endif else begin - SET_PLOT, 'Z', /COPY - device, set_resolution=[conresx,conresy],decomposed=0, Z_Buffer=0 - TVLCT, red, green, blue, /GET - !P.MULTI = [0, 1, 2] - erase -endelse - -;set up first color scale -loadct, 0 - -;strings for display -time = 'Time = ' -timeunit = ' yr' -timp = 'Total impacts (dotdash) = ' -tcrt = 'Total craters (line) = ' -ocrt = 'Countable craters (bold) = ' -epwr = 'Mean regolith depth = ' -mtot = 'Impacted mass = ' -c1lab = 'Min. Elevation' -c2lab = 'Mean Elevation' -c3lab = 'Max. Elevation' - -maxelev = max(surface_dem) -minelev = min(surface_dem) -medelev = mean(surface_dem) -c1 = string(minelev, format = '(F9.1)') + ' m' -c2 = string(medelev, format = '(F9.1)') + ' m' -c3 = string(maxelev, format = '(F9.1)') + ' m' -tslp = mean(regolith) - -;geometric saturation -geom = dblarr(2,2) -geomem = dblarr(2,2) -geomep = dblarr(2,2) -geomel = dblarr(2,2) -geom(0,0) = minx -geom(0,1) = maxx -geom(1,*) = 3.12636d0 -geomem(0,0) = minx -geomem(0,1) = maxx -geomem(1,*) = 0.156318d0 -geomep(0,0) = minx -geomep(0,1) = maxx -geomep(1,*) = 0.312636d0 -geomel(0,0) = minx -geomel(0,1) = maxx -geomel(1,*) = 0.0312636d0 - -; Remove zeros -nz = where(odist(5,*) ne 0.0,count) -if (count gt 0) then begin - odistnz = dblarr(6,count) - odistnz(*,*)= odist(*,nz) -endif else begin - odistnz = odist -endelse - -nz = where(tdist(5,*) ne 0.0,count) -if (count gt 0) then begin - tdistnz = dblarr(6,count) - tdistnz(*,*)= tdist(*,nz) -endif else begin - tdistnz = tdist -endelse - -nz = where(pdist(5,*) ne 0.0,count) -if (count gt 0) then begin - pdistnz = dblarr(6,count) - pdistnz(*,*)= pdist(*,nz) - -endif else begin - pdistnz = pdist -endelse - -; create r-plot array containing exactly 1 crater per bin -area = (gridsize * pix * 1d-3)^2 -plo = 1 -while (sqrt(2.0d0)^plo gt minx) do begin - plo = plo - 1 -endwhile -phi = plo + 1 -while (sqrt(2.0d0)^phi lt maxx) do begin - phi = phi + 1 -endwhile -n = phi - plo -sdist = dblarr(6,n + 1) -p = plo -for i=0,n do begin - sdist(0,i) = sqrt(2.d0)^p - sdist(1,i) = sqrt(2.d0)^(p+1) - sdist(2,i) = sqrt(sdist(0,i) * sdist(1,i)) - sdist(3,i) = 1.0d0 - sdist(5,i) = (sdist(2,i))^3 / (area * (sdist(1,i) - sdist(0,i))) - p = p + 1 -endfor -sdist(4,*) = reverse(total(sdist(3,*),/cumulative),2) - -plot, odistnz(2,*)*1.0d-3, odistnz(5,*), line=0, color=255, thick=2.0, $ - TITLE='Crater Distribution R-Plot', charsize=charfac*conresy/720, $ - XTITLE='Crater Diameter (km)', $ - XRANGE=[minx,maxx], /XLOG, XSTYLE=1, $ - YTITLE='R Value', $ - YRANGE=[5.0e-4,5.0e0], /YLOG, YSTYLE=1, $ - /device, pos = [0.085*conresx,(0.25 + profsize)*conresy,0.44*conresx,0.85*conresy] - -Dfac = sqrt(sqrt(2.0d0)) - -;display observed crater distribution -oplot, tdistnz(2,*)*1.0d-3, tdistnz(5,*), line=0, color=255, thick=1.0 -oplot, geom(0,*), geom(1,*), line=2, color=255, thick=1.0 -oplot, geomem(0,*), geomem(1,*), line=1, color=255, thick=1.0 -oplot, geomep(0,*), geomep(1,*), line=1, color=255, thick=1.0 -oplot, geomel(0,*), geomel(1,*), line=1, color=255, thick=1.0 -oplot, pdistnz(2,*)*1.0d-3, pdistnz(5,*), line=3, color=255, thick=1.0 -oplot, odistnz(2,*)*1.0d-3, odistnz(5,*), line=1, color=255, thick=1.0 -oplot, sdist(0,*), sdist(5,*), line=1, color=255, thick=1.0 -oplot, ph1(0,*)*Dfac*1.0d-3, ph1(1,*), psym=1, color=255, thick=1.0 -;oplot, ph2(0,*), ph2(1,*), psym=4, color=255, thick=1.0 -;oplot, ph3(0,*), ph3(1,*), psym=5, color=255, thick=1.0 - -; set up profile display -surfpos = dblarr(gridsize) -surfpro = dblarr(gridsize) -surfhpro = dblarr(gridsize) -surfrpro = dblarr(gridsize) -for k = 0,gridsize-1 do begin - surfpos(k) = (-0.5d0*(gridsize-1) + double(k)) * pix -endfor -surfpro(*) = surface_dem(*,gridsize/2-1) -surfhpro(*) = regolith(*,gridsize/2-1) -surfrpro(*) = surfpro(*) - surfhpro(*) - -; set up profile -;maxelevp = max(surfpro) -;minelevp = min(surfpro) -;medelevp = mean(surfpro) -;vertrange = 1.5 * abs(maxelevp-minelevp) -;vertadd = 0.5d0 * (vertrange - abs(maxelevp-minelevp)) -;horzrange = vertrange * 5.86666666667d0 -;if (horzrange gt (pix*gridsize)) then horzrange = pix*gridsize -;plot, surfpos(*)/1.0d3, surfpro(*)/1.0d3, line=0, color=255, thick=1.0, $ -; TITLE='Matrix Cross-Section', charsize=1.0, $ -; XTITLE='Horizontal Position (km)', $ -; XRANGE=[(-1.0d0*horzrange/2.0d3),(horzrange/2.0d3)], XSTYLE=1, $ -; YTITLE='Elevation (km)', $ -; YRANGE=[((minelevp - vertadd)/1.0d3),((maxelevp + vertadd)/1.0d3)], YSTYLE=1, $ -; /device, pos = [0.063*conresx,0.049*conresy,0.99*conresx,0.26*conresy] -;oplot, surfpos(*)*1.0d-3, surfrpro(*)*1.0d-3, line=0, color=255, thick=1.0 - -;display text -dum = string(format='(E10.4)', curyear) -parts = strsplit(dum, 'E', /extract) -fcuryear = strcompress(string(format='(F6.4,"x10!U",i,"!N")', parts[0], parts[1])) + ' yr' -ftimp = string(totalimpacts, format = '(I13)') -ftcnt = string(tdist(4,0), format = '(I13)') -focnt = string(odist(4,0), format = '(I13)') -ftslp = string(tslp, format = '(F13.3)') + ' m' - - -dum = string(format='(E10.4)',masstot) -parts = strsplit(dum, 'E', /extract) -fmtot = strcompress(string(format='(F6.4,"x10!U",i,"!N")', parts[0], parts[1])) + ' kg' - -; Display information text below the R-plot -texttop = 0.15 + profsize -textspc = 0.035 -textleft = 0.05*conresx -textbreak = 0.25*conresx -xyouts, textleft, (texttop - 0*textspc)*conresy, time, color=255, charsize=charfac*conresy/720., /device -xyouts, textleft, (texttop - 1*textspc)*conresy, timp, color=255, charsize=charfac*conresy/720., /device -xyouts, textleft, (texttop - 2*textspc)*conresy, tcrt, color=255, charsize=charfac*conresy/720., /device -xyouts, textleft, (texttop - 3*textspc)*conresy, ocrt, color=255, charsize=charfac*conresy/720., /device -xyouts, textleft, (texttop - 4*textspc)*conresy, mtot, color=255, charsize=charfac*conresy/720., /device -xyouts, textleft + textbreak, (texttop - 0*textspc)*conresy, fcuryear, color=255, charsize=charfac*conresy/720., /device -xyouts, textleft + textbreak, (texttop - 1*textspc)*conresy, ftimp, color=255, charsize=charfac*conresy/720., /device -xyouts, textleft + textbreak, (texttop - 2*textspc)*conresy, ftcnt, color=255, charsize=charfac*conresy/720., /device -xyouts, textleft + textbreak, (texttop - 3*textspc)*conresy, focnt, color=255, charsize=charfac*conresy/720., /device -xyouts, textleft + textbreak, (texttop - 4*textspc)*conresy, fmtot, color=255, charsize=charfac*conresy/720., /device - -; Display information text above the R-plot -xyouts, 0.0368*conresx, 0.972*conresy, c1lab, color=255, charsize=charfac*conresy/720., /device -xyouts, 0.195*conresx, 0.972*conresy, c2lab, color=255, charsize=charfac*conresy/720., /device -xyouts, 0.353*conresx, 0.972*conresy, c3lab, color=255, charsize=charfac*conresy/720., /device -xyouts, 0.0368*conresx, 0.944*conresy, c1, color=255, charsize=charfac*conresy/720., /device -xyouts, 0.195*conresx, 0.944*conresy, c2, color=255, charsize=charfac*conresy/720., /device -xyouts, 0.353*conresx, 0.944*conresy, c3, color=255, charsize=charfac*conresy/720., /device -xyouts, 0.100*conresx, 0.917*conresy, epwr, color=255, charsize=charfac*conresy/720., /device -xyouts, 0.278*conresx, 0.917*conresy, ftslp, color=255, charsize=charfac*conresy/720., /device - -;print to screen -print, ncount, ' time = ', curyear -;print, ncount, ' tot impacts = ', tdist -;print, ncount, ' tot craters = ', totalimpacts -;print, ncount, ' obs craters = ', ocrcnt -;print, ncount, ' regolith = ', tslp -;print, ncount, ' prod R = ', mean(prval(1,*)) - -;draw box around the main surface display & fill -displaysize = floor(0.53*conresx) ; size of displayed surface -surfxpos = floor(0.462*conresx) -surfypos = floor((0.05 + profsize)*conresy) -xbox = [surfxpos - 1,surfxpos - 1,surfxpos + displaysize,surfxpos + displaysize,surfxpos - 1] -ybox = [surfypos - 1,surfypos + displaysize,surfypos + displaysize,surfypos - 1,surfypos - 1] -plotS, xbox, ybox, /device, color=255 - -; for display -mapscaled = congrid(map,3,displaysize,displaysize,/interp) - -tv, mapscaled, surfxpos, surfypos, xsize=displaysize, ysize=displaysize, true=1, /device - -; ---------- saving window ---------- - -; save console window -fnum = string(ncount,format='(I6.6)') -; print screen to JPEG file -if (file_test('console',/DIRECTORY) eq 0) then begin - file_mkdir,'console' -endif -print, ncount, ' saving window' -fname = 'console/console' + fnum + '.jpg' -dwin = tvrd(true=1) -write_jpeg, fname, dwin, true=1, quality=90 - -end