From d62b66a6c416876fcc0eab8b01b05dc6e32b3646 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 8 Sep 2023 11:50:21 -0400 Subject: [PATCH] Cleaned up some obsolete files from the morphology test case example. --- .gitignore | 1 + .../complex/crater_profile.py | 176 ------------ .../complex/ctem_draw_surf.py | 1 - .../complex/ctem_driver.pro | 233 ---------------- .../complex/ctem_driver.py | 232 +--------------- .../complex/ctem_image_dem.pro | 41 --- .../complex/ctem_image_presentation.pro | 143 ---------- .../complex/ctem_image_regolith.pro | 33 --- .../complex/ctem_image_shaded_relief.pro | 52 ---- .../complex/ctem_io_read_input.pro | 131 --------- .../complex/ctem_io_read_old.pro | 45 --- .../complex/ctem_io_readers.py | 144 ---------- .../complex/ctem_io_writers.py | 261 ------------------ .../complex/ctem_window_display.pro | 249 ----------------- examples/morphology_test_cases/complex/getjob | 2 - src/Makefile.am | 2 +- 16 files changed, 5 insertions(+), 1741 deletions(-) delete mode 100644 examples/morphology_test_cases/complex/crater_profile.py delete mode 120000 examples/morphology_test_cases/complex/ctem_draw_surf.py delete mode 100755 examples/morphology_test_cases/complex/ctem_driver.pro mode change 100644 => 100755 examples/morphology_test_cases/complex/ctem_driver.py delete mode 100755 examples/morphology_test_cases/complex/ctem_image_dem.pro delete mode 100755 examples/morphology_test_cases/complex/ctem_image_presentation.pro delete mode 100755 examples/morphology_test_cases/complex/ctem_image_regolith.pro delete mode 100755 examples/morphology_test_cases/complex/ctem_image_shaded_relief.pro delete mode 100755 examples/morphology_test_cases/complex/ctem_io_read_input.pro delete mode 100755 examples/morphology_test_cases/complex/ctem_io_read_old.pro delete mode 100644 examples/morphology_test_cases/complex/ctem_io_readers.py delete mode 100644 examples/morphology_test_cases/complex/ctem_io_writers.py delete mode 100755 examples/morphology_test_cases/complex/ctem_window_display.pro delete mode 100755 examples/morphology_test_cases/complex/getjob diff --git a/.gitignore b/.gitignore index 2453fb1c..52cb12cd 100644 --- a/.gitignore +++ b/.gitignore @@ -72,3 +72,4 @@ examples/global-lunar-bombardment/scale.ipynb python/ctem/ctem.egg-info/ .idea/* +*.i90 diff --git a/examples/morphology_test_cases/complex/crater_profile.py b/examples/morphology_test_cases/complex/crater_profile.py deleted file mode 100644 index 12689d84..00000000 --- a/examples/morphology_test_cases/complex/crater_profile.py +++ /dev/null @@ -1,176 +0,0 @@ -import numpy as np -import pandas as pd -import os -from scipy import interpolate -import matplotlib.pyplot as plt -import matplotlib.cm as cm -from matplotlib.colors import LightSource - -def hillshade2(filename,gridsize,pix,array): - dpi = 300.0 #72.0 - ls = LightSource(azdeg=315, altdeg=45) - shadimg = ls.shade(array/pix,cmap=cm.gray) - - height = gridsize / dpi - width = height - fig = plt.figure(figsize=(width, height), dpi=dpi) - fig.figimage(shadimg, origin='lower') - plt.savefig(filename) - print(f'Image saved to {filename}') - return - -def image_dem(filename,gridsize,pix,CTEM_dem): - dpi = 300.0 #72.0 - # Create surface dem map - solar_angle = 20.0 - dem_map = np.copy(CTEM_dem) - np.roll(CTEM_dem, 1, 0) - dem_map = (0.5 * np.pi) + np.arctan2(dem_map, pix) - dem_map = dem_map - np.radians(solar_angle) * (0.5 * np.pi) - np.place(dem_map, dem_map > (0.5 * np.pi), 0.5 * np.pi) - dem_map = np.absolute(dem_map) - dem_map = 254.0 * np.cos(dem_map) - - # Save image to file - height = gridsize / dpi - width = height - fig = plt.figure(figsize=(width, height), dpi=dpi) - fig.figimage(dem_map, cmap=cm.gray, origin='lower') - plt.savefig(filename) - print(f'Image saved to {filename}') - - return - - -def cart2pol(x, y): - rho = np.sqrt(x**2 + y**2) - phi = np.arctan2(y, x) - return(rho, phi) - -def pol2cart(rho, phi): - x = rho * np.cos(phi) - y = rho * np.sin(phi) - return(x, y) - -crater_name = "Copernicus" -path = f'CSV/{crater_name}.csv' - - -unit_correction = 0.2727 #Ratio of Earth to Moon radii -crater = pd.read_csv(path) # * unit_correction * 1e-3 -crater['X'] = crater['X'] * unit_correction * 1e-3 -crater['Y'] = crater['Y'] * unit_correction * 1e-3 -crater['Z'] = crater['Z'] * np.sqrt(unit_correction) * 1e-3 - - -#print('Original input file head') -#print(crater.head()) -# -xshift = (crater['X'].max() + crater['X'].min())/2 -yshift = (crater['Y'].max() + crater['Y'].min())/2 -zshift = -1.0 -crater['X'] = crater['X'] - xshift -crater['Y'] = crater['Y'] - yshift -crater['Z'] = crater['Z'] - zshift -crater = crater.sort_values(by=['Y','X']) - -# -#print('Centered and sorted values head') -#print(crater.head()) -# -gridsize = int(np.sqrt(crater.shape[0])) -# -#print(f'Box size is {gridsize} pixels') -# -xres = 2 * crater['X'].max() / gridsize -yres = 2 * crater['Y'].max() / gridsize -#zres = 0.5 * (xres + yres) -print('X and Y resolutions should be the same') -print(f'X resolution is {xres} km / pix') -print(f'Y resolution is {yres} km / pix') - -xdata = crater['X'].values * xres -ydata = crater['Y'].values * xres -# -pix = xres -xdata = crater['X'].values.reshape(gridsize,gridsize)[0,:] -ydata = crater['Y'].values.reshape(gridsize,gridsize)[:,0] -orig_dem = crater['Z'].values.reshape(gridsize,gridsize) - -f = interpolate.interp2d(xdata, ydata, orig_dem, kind='cubic') - -proffig = plt.figure(1, figsize=(10,2)) -axes = proffig.add_subplot(111) -axes.set_xlabel("Distance from center (km)") -axes.set_ylabel("Elevantion (km)") - -rho = np.arange(0.0,xdata.max(),step=pix) -phi = np.arange(0,2*np.pi,step=30*np.pi/360.0) -zinterp = np.empty([rho.size,phi.size]) - -xprof = np.empty_like(zinterp) -yprof = np.empty_like(zinterp) - -for jdx,p in enumerate(phi): - for idx, r in enumerate(rho): - x,y = pol2cart(r,p) - xprof[idx,jdx], yprof[idx,jdx] = pol2cart(r, p) - zinterp[idx,jdx] = f(x,y) - xprof[idx,jdx] = x - yprof[idx, jdx] = y - axes.plot(rho,zinterp[:,jdx],':') - - -#Compare with CTEM output -gridsize = 2000 -dem_file = 'surface_dem.dat' -CTEM_dem = np.fromfile(dem_file , dtype = np.float64) -CTEM_dem.shape = (gridsize,gridsize) -CTEM_dem *= 1e-3 -pix = 0.200 -gridsize = 2000 -CTEM_dem.shape = (gridsize,gridsize) -xdata, ydata = np.mgrid[-gridsize/2:gridsize/2, -gridsize/2:gridsize/2] * pix - -xdata = np.arange(-gridsize/2,gridsize/2,1) * pix + pix/2 -ydata = np.arange(-gridsize/2,gridsize/2,1) * pix + pix/2 - -f = interpolate.interp2d(xdata, ydata, CTEM_dem, kind='cubic') - -rho = np.arange(0.0,xdata.max(),step=pix) -phi = np.arange(0,2*np.pi,step=30*np.pi/360.0) -zinterp = np.empty([rho.size,phi.size]) - -xprof = np.empty_like(zinterp) -yprof = np.empty_like(zinterp) - -for jdx,p in enumerate(phi): - for idx, r in enumerate(rho): - x,y = pol2cart(r,p) - xprof[idx,jdx], yprof[idx,jdx] = pol2cart(r, p) - zinterp[idx,jdx] = f(x,y) - xprof[idx,jdx] = x - yprof[idx, jdx] = y - axes.plot(rho,zinterp[:,jdx],'-') - - -plt.savefig(f'{crater_name}_profile.png',dpi=300,bbox_inches='tight') -os.system(f'open {crater_name}_profile.png') - -imgfile = f'{crater_name}.png' -image_dem(imgfile,gridsize,pix,orig_dem) -#os.system(f'open {imgfile}') - -#Let's reconstruct the original DEM by interpolating betwen profiles -#xp = xprof.flatten() -#yp = yprof.flatten() -#zp = zinterp.flatten() -# -#grid_x, grid_y = np.mgrid[-gridsize/2:gridsize/2, -gridsize/2:gridsize/2] * pix -#points = (xp,yp) -# -#zrecon = interpolate.griddata(points,zp,(grid_x,grid_y),method='cubic') -# -#imgfile = f'{crater_name}-reconstructed.png' -#image_dem(imgfile,gridsize,pix,zrecon) -#os.system(f'open {imgfile}') - diff --git a/examples/morphology_test_cases/complex/ctem_draw_surf.py b/examples/morphology_test_cases/complex/ctem_draw_surf.py deleted file mode 120000 index a1278a76..00000000 --- a/examples/morphology_test_cases/complex/ctem_draw_surf.py +++ /dev/null @@ -1 +0,0 @@ -../../../python/ctem_draw_surf.py \ No newline at end of file diff --git a/examples/morphology_test_cases/complex/ctem_driver.pro b/examples/morphology_test_cases/complex/ctem_driver.pro deleted file mode 100755 index 071af0f7..00000000 --- a/examples/morphology_test_cases/complex/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/morphology_test_cases/complex/ctem_driver.py b/examples/morphology_test_cases/complex/ctem_driver.py old mode 100644 new mode 100755 index 712a6157..22eb982f --- a/examples/morphology_test_cases/complex/ctem_driver.py +++ b/examples/morphology_test_cases/complex/ctem_driver.py @@ -1,229 +1,3 @@ -#!/usr/local/bin/python -# -#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 CTEM modules -import ctem_io_readers -import ctem_io_writers - -#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} - -#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) - -#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, - 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() +import ctem +sim = ctem.Simulation() +sim.run() diff --git a/examples/morphology_test_cases/complex/ctem_image_dem.pro b/examples/morphology_test_cases/complex/ctem_image_dem.pro deleted file mode 100755 index ad6fbd7e..00000000 --- a/examples/morphology_test_cases/complex/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/morphology_test_cases/complex/ctem_image_presentation.pro b/examples/morphology_test_cases/complex/ctem_image_presentation.pro deleted file mode 100755 index 05a7d65a..00000000 --- a/examples/morphology_test_cases/complex/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/morphology_test_cases/complex/ctem_image_regolith.pro b/examples/morphology_test_cases/complex/ctem_image_regolith.pro deleted file mode 100755 index 8d7f07a7..00000000 --- a/examples/morphology_test_cases/complex/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/morphology_test_cases/complex/ctem_image_shaded_relief.pro b/examples/morphology_test_cases/complex/ctem_image_shaded_relief.pro deleted file mode 100755 index 73986f84..00000000 --- a/examples/morphology_test_cases/complex/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/morphology_test_cases/complex/ctem_io_read_input.pro b/examples/morphology_test_cases/complex/ctem_io_read_input.pro deleted file mode 100755 index ea135668..00000000 --- a/examples/morphology_test_cases/complex/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/morphology_test_cases/complex/ctem_io_read_old.pro b/examples/morphology_test_cases/complex/ctem_io_read_old.pro deleted file mode 100755 index a9013687..00000000 --- a/examples/morphology_test_cases/complex/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/morphology_test_cases/complex/ctem_io_readers.py b/examples/morphology_test_cases/complex/ctem_io_readers.py deleted file mode 100644 index a58ba9ab..00000000 --- a/examples/morphology_test_cases/complex/ctem_io_readers.py +++ /dev/null @@ -1,144 +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 ('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/morphology_test_cases/complex/ctem_io_writers.py b/examples/morphology_test_cases/complex/ctem_io_writers.py deleted file mode 100644 index 673974d0..00000000 --- a/examples/morphology_test_cases/complex/ctem_io_writers.py +++ /dev/null @@ -1,261 +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 -from matplotlib.colors import LightSource -import matplotlib.cm as cm - -#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): - pix = parameters['pix'] - gridsize = parameters['gridsize'] - ve = 1.0 - mode = 'overlay' - azimuth = 300.0 #parameters['azimuth'] - solar_angle = 20.0 #parameters['solar_angle'] - - ls = LightSource(azdeg=azimuth, altdeg=solar_angle) - cmap = cm.cividis - - # If min and max appear to be reversed, then fix them - if (parameters['shadedminh'] > parameters['shadedmaxh']): - temp = parameters['shadedminh'] - parameters['shadedminh'] = parameters['shadedmaxh'] - parameters['shadedmaxh'] = temp - else: - parameters['shadedminh'] = parameters['shadedminh'] - parameters['shadedmaxh'] = parameters['shadedmaxh'] - - # 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) - else: - shadedminh = parameters['shadedminh'] - if (parameters['shadedmaxhdefault'] == 1): - shadedmaxh = numpy.amax(surface_dem) - else: - shadedmaxh = parameters['shadedmaxh'] - - - - dem_img = ls.shade(surface_dem, cmap=cmap,blend_mode=mode, - vert_exag=ve, dx=pix, dy=pix, - vmin=shadedminh, vmax=shadedmaxh) - - #Generate image to put into an array - height = gridsize / dpi - width = gridsize / dpi - fig = matplotlib.pyplot.figure(figsize=(width, height), dpi=dpi) - ax = matplotlib.pyplot.axes([0, 0, 1, 1]) - ax.imshow(dem_img,interpolation="nearest") - matplotlib.pyplot.axis('off') - # Save image to file - filename = parameters['workingdir'] + 'shaded' + os.sep + "shaded%06d.png" % parameters['ncount'] - matplotlib.pyplot.savefig(filename,dpi=dpi,bbox_inches=0) - 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 \ No newline at end of file diff --git a/examples/morphology_test_cases/complex/ctem_window_display.pro b/examples/morphology_test_cases/complex/ctem_window_display.pro deleted file mode 100755 index d6edf5b9..00000000 --- a/examples/morphology_test_cases/complex/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 diff --git a/examples/morphology_test_cases/complex/getjob b/examples/morphology_test_cases/complex/getjob deleted file mode 100755 index 35b613ff..00000000 --- a/examples/morphology_test_cases/complex/getjob +++ /dev/null @@ -1,2 +0,0 @@ -#!/bin/bash -rsync -vah --progress "brown.rcac.purdue.edu:/scratch/brown/daminton/NPF-global-Kd0.0001/" ./ diff --git a/src/Makefile.am b/src/Makefile.am index 69c0278e..ff1e21de 100755 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -6,7 +6,7 @@ HEAPARR = -heap-arrays OPTREPORT = -qopt-report=5 IPRODUCTION = -g -traceback -no-wrap-margin -assume byterecl -O3 -qopt-prefetch=0 -sox $(PAR) $(SIMDVEC) $(HEAPARR) IDEBUG = -O0 -g -traceback -debug all -nogen-interfaces -assume byterecl -m64 -heap-arrays -FR -no-pie -no-ftz -fpe-all=0 -mp1 -fp-model strict -fpe0 -align all -pad -ip -prec-div -prec-sqrt -assume protect-parens -CB -no-wrap-margin -init=snan,arrays -AM_FCFLAGS = $(IPRODUCTION) +AM_FCFLAGS = $(IDEBUG) #ifort debug flags #gfortran optimized flags