Merge branch 'realistic' into debugGlassMerge

examples/morphology_test_cases/complex/ctem.in

@@ -21,13 +21,13 @@ numintervals 1                                ! Total number of intervals
 restart      F                                  ! Restart a previous run
 impfile      NPFextrap.dat          ! Impactor SFD rate file (col 1: Dimp (m), col 2: ! impactors > D (m**(-2) y**(-1))
 popupconsole F                                  ! Pop up console window every output interval 
-saveshaded   F                                  ! Output shaded relief images  
+saveshaded   T                                  ! Output shaded relief images  
 saverego     F                                  ! Output regolith map images 
 savepres     F                                  ! Output simplified console display images (presentation-compatible images)
 savetruelist T                                  ! Save the true cumulative crater distribution for each interval (large file size)
 sfdcompare   LOLASethCraterCatalogv8gt20-binned.dat         ! File name for the SFD comparison file used in the console display
-shadedminh   -85.0                              ! Minimum height for shaded relief map (m) (Default - automatic)
-shadedmaxh   85.0                               ! Maximum height for shaded relief map (m) (Default - automatic)
+shadedminh   -5000.0                            ! Minimum height for shaded relief map (m) (Default - automatic)
+shadedmaxh   5000.0                             ! Maximum height for shaded relief map (m) (Default - automatic)
 runtype      single                             ! Run type: options are normal / statistical
                                                 !    single:      craters accumulate in successive intervals
                                                 !    statistical: surface is reset between intervals
@@ -64,7 +64,7 @@ maxcrat         1.00e-98                          ! Fraction of gridsize that ma
 killatmaxcrater F                               ! Stop the run if a crater larger than the maximum is produced - Default F      
 basinimp        35.0e3                          ! Size of impactor to switch to lunar basin scaling law - Default is to ignore
 docollapse      T                               ! Do slope collapse - Default T
-dosoftening     T                               ! Do ejecta softening - Default T
+dosoftening     F                               ! Do ejecta softening - Default T
 doangle         T                               ! Vary the impact angle. Set to F to have only vertical impacts - Default T
 Kd1             0.0001
 psi             2.000

examples/morphology_test_cases/complex/ctem_io_writers.py

@@ -15,6 +15,8 @@
 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
@@ -123,66 +125,51 @@ def image_regolith(parameters, regolith):
     return
 
 def image_shaded_relief(parameters, surface_dem):
-    #The color scale and appearance of this do not quite match the IDL version
+    pix = parameters['pix']
+    gridsize = parameters['gridsize']
+    ve = 1.0
+    mode = 'overlay'
+    azimuth = 300.0 #parameters['azimuth']
+    solar_angle = 20.0  #parameters['solar_angle']
 
-    #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)
+    ls = LightSource(azdeg=azimuth, altdeg=solar_angle)
+    cmap = cm.cividis
 
-    #If min and max appear to be reversed, then fix them
+    # If min and max appear to be reversed, then fix them
     if (parameters['shadedminh'] > parameters['shadedmaxh']):
         temp = parameters['shadedminh']
-        shadedminh = parameters['shadedmaxh']
-        shadedmaxh = temp
+        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']
-        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
+    if (parameters['shadedmaxhdefault'] == 1):
+        shadedmaxh = numpy.amax(surface_dem)
     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
+        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']
-    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)    
+    matplotlib.pyplot.savefig(filename,dpi=dpi,bbox_inches=0)
     return
 
 def create_rplot(parameters,odist,pdist,tdist,ph1):

src/ejecta/ejecta_emplace.f90

@@ -206,6 +206,7 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,age,age_r
 
          lradsq = (crater%xl - xp)**2 + (crater%yl - yp)**2
          lrad = sqrt(lradsq)
+         if (lrad < crater%ejrad) cycle
 
          ! Estimate ejecta pattern distortion due to target surface angle and topography
          ! This must be done iteratively because the ejection distance and ejection angle vary
@@ -222,6 +223,7 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,age,age_r
 
             baseline = ((i * crater%xslp) + (j * crater%yslp)) * user%pix
             craterslope = atan(baseline / lrad)
+            if ((n == 1) .and. abs(craterslope) < epsilon(1._DP)) exit
             if (craterslope > maxslp) maxslp = craterslope
 
             ejheight = erad * sin(craterslope) + crater%melev