How can I apply this function Y = d+asin(bx-c) in the code shown below? Where a represents the amplitude

1 
program Make_2D_Model
	implicit none
     integer:: i,j,nx,nz 
     integer :: idepth
 character*1,allocatable,dimension(:,:) :: icode
nx=1000  
nz=500
idepth=100
	allocate(icode(1:nx,1:nz))	
	do i=1,nx
		do j=1,nz
			if(j<=idepth)then
				icode(i,j)='1'
			elseif(j>idepth)then
				icode(i,j)='3'		
            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)

    end program Make_2D_Model

Snipaste_2022-06-18_17-41-54

2

If I understand your code:

			if(j<=idepth)then
				icode(i,j)='1'
			elseif(j>idepth)then
				icode(i,j)='3'		
            endif

it defines a flat surface at j=idepth, ‘1’ is red and ‘3’ is blue in your figure. But as its name suggests, idepth depends on i. Just before the if, you should write something like idepth = d + a*sin(b*i-c), with the good values for a, b, c, d.
And if you want also the flat segments, it would be something like idepth = min(d+a*sin(b*i-c), e) where e is the z position of those flat segments.

3

Note also that in modern Fortran you don’t need to define yourself the unit value (2012 in your code). You can write open(newunit=my_file, file=.... where my_file is an integer variable that will automatically receive an available unit number.

4

I’m usin Visual Studio 2013

5

newunit was introduced in Fortran 2008. It depends on the Fortran compiler you use.

6

IMG-20220619-WA0000

1 Like