implicit none
integer:: i,j,nx,nz
character*1,allocatable,dimension(:,:) :: icode
real,parameter :: pi=3.1415926
integer :: mm,ang
real :: kk,h,dx,dz
nx=1000 !grid number in x direction
nz=500 !grid number in z direction
dx=1.0 !the horizontal spatial step length (unit: m)
dz=1.0 !the vertical spatial step length (unit: m)
ang=10.0 !Dipping angle of free surface (just for dipping-free-surface model)
h=10.0 !Real thickness of surface layer
allocate(icode(1:nx,1:nz))
mm=(h/(cos(ang*pi/180))/dx)+1
kk=tan(ang*pi/180)
do i=1,(nx-1)/2
do j=1,nz
if(j<kk*(nx-i))then
icode(i,j)='1'
elseif((j>kk*(nx-i)+1).and.(j<kk*(nx-i)+mm))then
icode(i,j)='3'
elseif(j>=kk*(nx-i)+mm)then
icode(i,j)='4'
else
icode(i,j)='2'
endif
enddo
enddo
do i=(nx-1)/2+1,nx
do j=1,nz
if(j<kk*(nx-i))then
icode(i,j)='1'
elseif((j>kk*(nx-i)+1).and.(j<kk*(nx-i)+(mm-1)/2+1))then
icode(i,j)='3'
elseif(j>=kk*(nx-i)+(mm-1)/2+1)then
icode(i,j)='4'
else
icode(i,j)='2'
endif
enddo
enddo
open(2012,file='model.txt',status='unknown')
do j=1,nz
write(2012,1000)(icode(i,j),i=1,nx)
enddo
close(2012)
1000 format(<nx>a1)
1 Like
Hi Mohammed,
Given your previous questions,
- How to make use of cosine Function to make this shape
- How I can draw a surface wave model like the figure below?
- How can I apply this function Y = d+asin(bx-c) in the code shown below? Where a represents the amplitude
- Create a topography, two-layer model I need help with what is the wrong with this code
- How to type this part of the code correctly?
- Create Rectangular Cavity
- I need to model an incline rectangle shape
I’d really encourage you to pursue a more modular and flexible program structure, and also to search for better ways (as in geometrical representations) to describe your models. In most of your threads the convoluted if-elseif blocks and conditions are very difficult to reason about. I have the feeling like you don’t really care about this, since these are probably one-off calculations.
I tried to demonstrate to you a different approach yesterday with the point-in-polygon test. Today I’m sharing with you an implicit line representation I once found in the lectures notes on Analysis of Binary Images by Robyn Owens from the University of Western Australia.
program fault
implicit none
integer, parameter :: nx = 1001, nz = 501
real, parameter :: pi = 4.*atan(1.)
real :: depth, theta_deg, theta, x0, z0, rho
real :: dx, dz
character(1), allocatable :: icode(:,:)
! grid spacing
dx = 1.
dz = 1.
allocate(icode(1:nx,1:nz))
! fault depth
depth = 100.
! fault angle
theta_deg = 10. ! degrees
theta = theta_deg * (pi/180.) ! radian
! fault cuts through this point
x0 = 500.
z0 = -50.
!
! the fault is described by the implicit line
!
! x*sin(phi) - z*cos(phi) + rho = 0
!
! where phi = pi - theta. The meanings of the parameters are as follows
!
! rho ... distance from the origin (0,0) to the closest point
! on the line
! phi ... angle between the line and the x-axis
! theta ... adjacent angle to phi
!
! a full description of the math describing the line is available at
! https://homepages.inf.ed.ac.uk/rbf/CVonline/LOCAL_COPIES/OWENS/LECT2/node3.html#SECTION00031100000000000000
!
! calculate line constant from the known angle, and a point
! lying on the fault
rho = z0*cos(pi - theta) - x0*sin(pi - theta)
! fill icode mask
call two_layers_with_fault( &
ndim = [nx,nz], &
spacing = [dx,dz], &
icode = icode, &
rho = rho, &
theta = theta, &
depth = depth)
call save_txt('model.txt',icode)
call save_pgm('model.pgm',icode,ngray=3)
contains
!> Fill mask array for model with two horizontal layers and a fault
!>
!> The mask values are
!> - 0 or 1 (top layer), depending on the side of the fault line
!> - 2 (bottom layer)
!>
subroutine two_layers_with_fault(ndim,spacing,icode,rho,theta,depth)
integer, intent(in) :: ndim(2) ! (nx, nz)
real, intent(in) :: spacing(2) ! (dx, dz)
character(1), intent(out) :: icode(:,:)
real, intent(in) :: rho, theta, depth
! geometry parameters
real :: xx, zz
integer :: j, i
do j = 1, ndim(2)
do i = 1, ndim(1)
! we assume implicitly that point (i=1,j=1) is the origin (0,0),
! increasing j corresponds to negative z-values in
! physical coordinate,
xx = (i-1) * spacing(1)
zz = 0.0 - (j-1) * spacing(2)
if (abs(zz) < depth) then
! fault layer
icode(i,j) = merge('0','1', &
is_below_fault(xx,zz,rho,theta))
else
! bottom layer
icode(i,j) = '2'
end if
end do
end do
end subroutine
!> Check if a point is below the fault line
!>
!> The line must be prescribed by the implicit form
!>
!> x * sin(pi - theta) - y * cos(pi - theta) + rho = 0
!>
pure function is_below_fault(xx,zz,rho,theta) result(res)
real, intent(in) :: xx, zz ! point of interest
real, intent(in) :: rho ! line constant
real, intent(in) :: theta ! angle in radians
real, parameter :: pi = 4.*atan(1.)
logical :: res
if (xx*sin(pi - theta) - zz*cos(pi - theta) + rho < 0) then
res = .true.
else
res = .false.
end if
end function
!> Save as dumb ASCII txt file (use .txt extension)
subroutine save_txt(filename,icode)
character(len=*), intent(in) :: filename
character(1), intent(in) :: icode(:,:)
integer :: i,j,unit
open(newunit=unit,file=filename,status='unknown')
do j = 1, nz
write(unit,'(*(A1))') (icode(i,j), i = 1, nx)
end do
close(unit)
end subroutine
!> Save as plain Portable GrayMap (use .pgm extension)
subroutine save_pgm(filename,icode,ngray)
character(len=*), intent(in) :: filename
character(1), intent(in) :: icode(:,:)
integer, intent(in) :: ngray
! number of gray values
integer :: i, j, unit
open(newunit=unit,file=filename,status='unknown')
! a description of the format is available at
! https://en.wikipedia.org/wiki/Netpbm
! header part
write(unit,'(A)') 'P2'
write(unit,'(I0,1X,I0)') size(icode,1), size(icode,2)
write(unit,'(I0)') ngray
! pixel data, gray values
do j = 1, size(icode,2)
write(unit,'(*(A1,:,1X))') (icode(i,j), i = 1, size(icode,1))
end do
close(unit)
end subroutine
end program
Picking the right mathematical representation can make the job much easier for you.
Here’s the output I got in the second try, after adjusting a < symbol in the is_below_fault function.
3 Likes
Thank you very much for helping me
1 Like
I study engineering geophysics; I use the surface seismic wave method to analyze the structures near the surface, and this is the first time I have dealt with Fortran language.
1 Like
I hope the snippets shown have been useful.
Besides from the book Computational Seismology, 2016, by Heiner Igel, I’ve also seen a few seismologically-suited methods in A Primer on Radial Basis Functions with Applications to the Geosciences, 2015, by Fornberg & Flyer.
A slightly younger paper, by Martin & Fornberg, 2017, discusses how to derive high-accuracy finite difference approximations for material interfaces (discontinuities).
For a state of the art seismology code available today, you can have a look at SeisSol. This is a C++ code which has received a lot of attention, and is able to run efficiently on large supercomputer. There’s always something to learn by peaking underneath the surface (pun intended).
1 Like