Skip to content

Commit

Permalink
Browse files Browse the repository at this point in the history 
Improved central peak model with Gaussian function that can also do peak rings
  • Loading branch information
@daminton
daminton committed Apr 9, 2020
1 parent c1b16e1 commit fdedbab
Showing 1 changed file with 28 additions and 15 deletions.
43 changes: 28 additions & 15 deletions src/crater/crater_realistic_topography.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -69,12 +69,13 @@ subroutine complex_peak(user,surf,crater,rn,deltaMtot)
real(DP),intent(inout) :: deltaMtot
end subroutine complex_peak

subroutine complex_wall_texture(user,surf,crater,rn,deltaMtot)
subroutine complex_wall_texture(user,surf,crater,domain,rn,deltaMtot)
use module_globals
implicit none
type(usertype),intent(in) :: user
type(surftype),dimension(:,:),intent(inout) :: surf
type(cratertype),intent(inout) :: crater
type(domaintype),intent(in) :: domain
real(DP),dimension(:),intent(in) :: rn
real(DP),intent(inout) :: deltaMtot

Expand All @@ -96,8 +97,7 @@ end subroutine complex_terrace
call random_number(rn)

call complex_terrace(user,surf,crater,rn,deltaMtot)
call complex_wall_texture(user,surf,crater,rn,deltaMtot)
call crater_slope_collapse(user,surf,crater,domain,(0.33_DP * user%pix)**2,deltaMtot)
call complex_wall_texture(user,surf,crater,domain,rn,deltaMtot)
call complex_floor(user,surf,crater,rn,deltaMtot)
call complex_peak(user,surf,crater,rn,deltaMtot)

Expand All @@ -122,20 +122,28 @@ subroutine complex_peak(user,surf,crater,rn,deltaMtot)
integer(I4B) :: i,j,inc,xpi,ypi

! Complex crater peak
real(DP), parameter :: rad = 0.20_DP ! Radial size of central peak relative to frad
integer(I4B), parameter :: num_octaves = 8 ! Number of Perlin noise octaves
integer(I4B), parameter :: num_octaves = 10 ! Number of Perlin noise octaves
integer(I4B), parameter :: offset = 1000 ! Scales the random xy-offset so that each crater's random noise is unique
real(DP), parameter :: xy_noise_fac = 16.0_DP ! Spatial "size" of noise features at the first octave
real(DP), parameter :: xy_noise_fac = 14.0_DP ! Spatial "size" of noise features at the first octave
real(DP), parameter :: freq = 2.0_DP ! Spatial size scale factor multiplier at each octave level
real(DP), parameter :: pers = 0.80_DP ! The relative size scaling at each octave level
real(DP), parameter :: pers = 0.50_DP ! The relative size scaling at each octave level

real(DP) :: FWHM,a,b,c,z ! Gaussian parameters

! Topographic noise parameters
real(DP) :: xynoise, znoise
integer(I4B) :: octave
real(DP) :: noise

!Preliminary empirical fits based on a handful of complex craters
!More work needs to be done to make these models more robust
FWHM = 0.1_DP
a = crater%peakheight
b = 0.003_DP * ((1e-3_DP * crater%fcrat)**(1.75_DP)) / (1e-3_DP * crater%fcrat) ! Make peak rings for sufficiently large craters
c = FWHM / (2 * sqrt(2 * log(2._DP)))
!*********************

inc = max(min(nint(rad * crater%frad / user%pix),PBCLIM*user%gridsize),1) + 1
inc = max(min(nint(0.5_DP * crater%floordiam / user%pix),PBCLIM*user%gridsize),1) + 1
crater%maxinc = max(crater%maxinc,inc)

do j = -inc,inc
Expand All @@ -150,14 +158,16 @@ subroutine complex_peak(user,surf,crater,rn,deltaMtot)
xbar = xpi * user%pix - crater%xl
ybar = ypi * user%pix - crater%yl

areafrac = util_area_intersection(rad * crater%frad,xbar,ybar,user%pix)
areafrac = util_area_intersection(0.5_DP * crater%floordiam,xbar,ybar,user%pix)

r = sqrt(xbar**2 + ybar**2) / crater%frad
! Model central peak / peak ring as Gaussian function
z = a * exp(-(r - b)**2 / (2 * c**2))

noise = 0.0_DP
do octave = 1, num_octaves
xynoise = xy_noise_fac * freq ** (octave - 1) / crater%fcrat
znoise = crater%peakheight * (1.0_DP - r / rad) * pers ** (octave - 1)
znoise = z * pers ** (octave - 1)
noise = noise + util_perlin_noise(xynoise * xbar + offset * rn(1), &
xynoise * ybar + offset * rn(2)) * znoise
end do
Expand Down Expand Up @@ -238,7 +248,7 @@ subroutine complex_floor(user,surf,crater,rn,deltaMtot)
end subroutine complex_floor


subroutine complex_wall_texture(user,surf,crater,rn,deltaMtot)
subroutine complex_wall_texture(user,surf,crater,domain,rn,deltaMtot)
use module_globals
use module_util
use module_crater, EXCEPT_THIS_ONE => crater_realistic_topography
Expand All @@ -248,6 +258,7 @@ subroutine complex_wall_texture(user,surf,crater,rn,deltaMtot)
type(usertype),intent(in) :: user
type(surftype),dimension(:,:),intent(inout) :: surf
type(cratertype),intent(inout) :: crater
type(domaintype),intent(in) :: domain
real(DP),dimension(:),intent(in) :: rn
real(DP),intent(inout) :: deltaMtot

Expand All @@ -265,7 +276,7 @@ subroutine complex_wall_texture(user,surf,crater,rn,deltaMtot)
integer(I4B), parameter :: num_octaves = 10 ! Number of Perlin noise octaves
integer(I4B), parameter :: offset = 4000 ! Scales the random xy-offset so that each crater's random noise is unique
real(DP), parameter :: xy_noise_fac = 4.0_DP ! Spatial "size" of noise features at the first octave
real(DP), parameter :: noise_height = 8.0e-3_DP ! Spatial "size" of noise features at the first octave
real(DP), parameter :: noise_height = 7.0e-3_DP ! Spatial "size" of noise features at the first octave
real(DP), parameter :: freq = 2.0_DP ! Spatial size scale factor multiplier at each octave level
real(DP), parameter :: pers = 1.20_DP ! The relative size scaling at each octave level

Expand Down Expand Up @@ -311,6 +322,8 @@ subroutine complex_wall_texture(user,surf,crater,rn,deltaMtot)
end do
end do

call crater_slope_collapse(user,surf,crater,domain,(0.33_DP * user%pix)**2,deltaMtot)

return
end subroutine complex_wall_texture

Expand Down Expand Up @@ -340,13 +353,13 @@ subroutine complex_terrace(user,surf,crater,rn,deltaMtot)
! Complex crater wall
integer(I4B) :: nscallops ! Approximate number of scallop features on edge of crater
integer(I4B),parameter :: offset = 3000 ! Scales the random xy-offset so that each crater's random noise is unique
real(DP),parameter :: noise_height = 0.15_DP ! Vertical height of noise features as a function of crater radius at the first octave
real(DP),parameter :: noise_height = 0.20_DP ! Vertical height of noise features as a function of crater radius at the first octave
! Higher values make the scallop features broader
real(DP),parameter :: noise_slope = 0.6_DP ! The slope of the power law function that defines the noise shape for scallop features
! Higher values makes the inner sides of the scallop features more rounded, lower values
! make them have sharper points at the inflections
integer(I4B),parameter :: nterraces = 6 ! Number of terraces
real(DP),parameter :: terracefac = 2.0_DP ! Power law scaling for relative terrace sizes
integer(I4B),parameter :: nterraces = 8 ! Number of terraces
real(DP),parameter :: terracefac = 3.0_DP ! Power law scaling for relative terrace sizes
integer(I4B) :: terrace ! The current terrace
real(DP) :: router ! The radius of the terrace outer edge
real(DP) :: rinner ! The radius of the terrace outer edge
Expand Down

0 comments on commit fdedbab

Please sign in to comment.