Download Earthquake Catalogs from Global CMT website

In seismology, we always need to download and check the event information about the events. This python script can download the event catalog from the website to the local computer for the given range of time.

Running this program is very simple:
Screen Shot 2017-08-25 at 11.47.09 PM.png

The user just needs to input the time range for the earthquakes e.g., 2000/05-2009/08.

 

Requirements: Python 3

To download the program, please click here.

Solving nth order differential equation for a given boundary condition

We are going to solve the differential equation with the boundary conditions

\psi_{xx} +[100 - \beta ]\psi = 0

\implies \psi_{xx} = [\beta - 100] \psi

Let’s take the simpler boundary condition

\psi(x=\pm 1) =0.

Now, the first thing which we wanna do is to convert this equation in the first order form. After that, we can simply use the 4th order Runge-Kutta to obtain higher order solution.

To transform our 2nd order differential equation into first order, we define two variables. To transform nth order equations we need n such variables. This is the standard process to convert any differential equation into the system of first order differential equations.

Lets take

y_1 = \psi (1)

\implies y_1' = \psi_x

\implies y_1' = y_2

y_2= \psi_x (2) [derivative of \psi]

\implies  y_2' = \psi_{xx}

\implies  y_2' = [\beta - 100] \psi

So, we have transformed the 2nd order differential equation into two first order differential equations.

And our boundary condition will become y_1(\pm 1) =0

Now for the differential equations of the type

\psi_{xx} = [\beta - 100] \psi

\implies \psi_{xx} = a \psi

have the solution \psi = e^{\pm \sqrt{a}x} (for positive a). Now this solution doesn’t satisfy our condition \psi(x=\pm 1) =0. Exponential functions are monotonic and hence cannot be 0 at two points. So, this tells us that a is not positive, i.e., \beta \ngtr 100.

Now, if \beta < 100, then our eqn becomes

\psi_{xx} = - a \psi

The solution to this is of the form of cosines and sines, which has potential to satisfy the boundary condition. So, we have \beta < 100.

Defining Function in MATLAB:

function [ rhs ] = yfunc( x,ic,dummy,beta )
y1=ic(1);
y2=ic(2);

rhs = [y2; (beta -100)*y1];
end

 

Calculating Solutions:

clear; close all; clc
set(0,'defaultlinelinewidth',2);
set(0,'defaultaxeslinewidth',1.5);
set(0,'defaultpatchlinewidth',2);
set(0, 'DefaultAxesFontSize', 14)

tol=10^(-4); % define a tolerance level to be achieved
% by the shooting algorithm

xspan=[-1 1]; %boundary conditions % define the span of the computational domain

A=1; % define the initial slope at x=-1
ic=[0 A]; % initial conditions: x1(-1)=0, x1?(-1)=A

n0=100; % define the parameter n0
dbeta=1;
betasol=[];
beta_start=n0; % beginning value of beta
for jj=1:5 % begin mode loop
    beta=beta_start; % initial value of eigenvalue beta
    dbeta=n0/100; % default step size in beta
    for j=1:1000 % begin convergence loop for beta
        [t,y]=ode45('yfunc',xspan, ic,[],beta); % solve ODEs
       if y(end,1)*((-1)^(jj+1)) > 0 %this checks if the beta needs to be higher or lower for convergence
           beta = beta-dbeta;
       else
           beta = beta + dbeta/2; %uses bisection to converge to the solution
           dbeta=dbeta/2;
       end
       if abs(y(end,1)) < tol % check for convergence
           betasol=[betasol beta]; % write out eigenvalue
           break % get out of convergence loop
       end 
    end
    beta_start=beta-0.1; % after finding eigenvalue, pick
    % new starting value for next mode
    norm=trapz(t,y(:,1).*y(:,1)); % calculate the normalization
    plot(t,y(:,1)/sqrt(norm)); hold on
end
betasol
legend(sprintf('\\beta_1 = %.4f',betasol(1)),sprintf('\\beta_2 = %.4f',betasol(2)),sprintf('\\beta_3 = %.4f',betasol(3)),...
    sprintf('\\beta_4 = %.4f',betasol(4)),sprintf('\\beta_5 = %.4f',betasol(5)));
xlabel('x'); ylabel('\psi (x)','FontWeight','bold')
saveas(gcf,'DiffEqnSol','jpg')

DiffEqnSol

 

Introduction to Genetics Algorithm (GA) (Part 2)

To find a basic introduction of GA, the first part can be found here.

III. Examples using Genetics Algorithm

In these examples, we will use Matlab and its function ga to apply GA for the optimization problem. For the manual of using this function, you can find it at https://www.mathworks.com/help/gads/ga.html or type in Matlab:

help ga

III.1. Optimization of a polynomial function

Problem: find the integer x  value that minimizes function:

y = 0.2x2 + 50/x

First, we define the function:

function y = simple_fitness(x)
 y = 0.2*x.^2+(50./x);

And use Matlab to finish the job:

clear; close all; clc

x=linspace(0,40,100);
y=simple_fitness(x);
plot(x,y);

rng(1,'twister') % for reproducibility
lb=[1];
[x,fval] = ga(@simple_fitness,1,[],[],[],[],lb

xlabel('Size of cup');
ylabel('Objective Function Value');
title(sprintf('Minimum of the function is %.2f at x = %.2f',min(y), x))
grid on

And here is the result:

untitled.jpg

III.2. Find the best location of the earthquakes

Problem: we have a set of 5 seismic stations with coordinates: [-2 3 0; 1 3 0; -2 -1 0; 0 -3 0; 2 -2 0]

The travel time of seismic ray from EQ to stations is calculated from:

t.jpg

To simplify the problem, we choose v = 6 km

clear; close all; clc
stnx=[-2 1 -2 0 2];
stny=[3 3 -1 -3 -2];
stnz=[0 0 0 0 0];
EQ=[2 2 -2];

to=sqrt((EQ(1)-stnx).^2+(EQ(2)-stny).^2+(EQ(3)-stnz).^2)/6
to=to+0.05*randn(size(to))
rms = @(mp) sum(sqrt((1./(length(to)-4)).*((sqrt((mp(1)-stnx).^2+(mp(2)-stny).^2+(mp(3)-stnz).^2)/6)-to).^2));
% rms([2 2 -2])
ub=[10 10 0]; %upper bound
lb=[0 0 -10]; %lower bound
options = gaoptimset('Display','iter','PlotFcn',@gaplotbestf,'Generations',150)
% options = gaoptimset('Display','iter','PlotFcn',@gaplotbestf,'TolCon',1e-3);

mp = ga(rms,3,[],[],[],[],lb,ub,[],options)

figure; hold on; grid on;
scatter3(stnx, stny, stnz, 'b', 'filled', '^');
scatter3(mp(1), mp(2), mp(3), 'r', '*');
scatter3(EQ(1), EQ(2), EQ(3), 'g', '*');
view(3); xlabel('X axis'); ylabel('Y axis'); zlabel('Z axis');
legend('Stations','EQ location from GA','Real EQ location');

untitled3untitled2

Utpal Kumar

IESAS

 

Locating Earthquake using Geiger’s Method

Earthquake location problem is old, however, it is still quite relevant. The problem can be stated as to determine the hypocenter (x0,y0,z0) and origin time (t0) of the rupture of fault on the basis of arrival time data of P and S waves. Here, we have considered hypocenter in the cartesian coordinate system. It can be readily converted to geographical coordinates using appropriate conversion factor. The origin time is the start time of the fault rupture.

Here, our data are the P and S arrival times at the stations located at the surface (xi,yi,0). We have also assumed a point source and a constant velocity medium for the simplicity sake.

The arrival time of the P or S wave at the station is equal to the sum of the origin time of the earthquake and the travel time of the wave to the station.

The arrival time, ti=Ti + t0, where Ti is the travel time of the particular phase and t0 is the origin time of the earthquake.

Ti = sqrt((xi-x0)^2 + (yi-y0)^2 + z0^2)/v; where v is the P or S wave velocity.

Now, we have to determine the hypocenter location and the origin time of the Earthquake. This problem is invariably ill-posed as the data are  contaminated by noise. So, we can pose this problem as an over-determined case. One way to solve this inconsistent problem for the hypocenter and origin time is to minimize the least square error of the observed and predicted arrival times at the stations.

Before that, we need to formulate the problem. As evident from the above equation, the expression for arrival time and eventually the least square error function is nonlinear. This nonlinear inverse problem can be solved using the Newton’s method. We will utilize the derivative of the error function in the vicinity of the trial solution (initializing the problem with the initial wise guess) to devise a better solution. The next generation of the trial solution is obtained by expanding the error function in a Taylor series about the trial solution. This process continues until the error function becomes considerably small. Then we accept that  trial solution as the estimated model parameters.

d = G m, where d is the matrix containing the arrival time data, m is the matrix containing the model parameters and G is the data kernel which maps the data to the model.

There are some instances where the matrix can become underdetermined and the least square would fail such as when the data contains only the P wave arrival times. In such cases the solution can become nonunique. We solve this problem as a damped least square case. Though it can also be using the singular value decomposition which allows the easy identification of the underdetermined case and then partitioning the overdetermined and underdetermined in upper and lower part of the kernel matrix and dealing with both part separately.

m = [G’G + epsilon]G’d

The parameter epsilon is chosen by trial and error to yield a solution that has reasonably small prediction error.

In this problem, we deduce the gradient of the travel time by examining the geometry of the ray as it leaves the source. If the earthquake is moved a small distance s  parallel to the ray in the direction of the receiver, then the travel time is simply decreased by s/v, where v is the velocity of the phase. If the earthquake is moved a small distance perpendicular to the ray, then the change in travel time is negligible since the new path will have nearly same length as the old one.

T = r/v

If the earthquake is moved parallel to the ray path,

T” = (r-s)/v = r/v – s/v

delta T = T” – T = r/v – s/v – r/v  = – s/v

If the earthquake has moved perpendicular to the raypath,

delta T = 0

So, delta T = – s/v, where is a unit vector tangent to the ray at the source and points toward the receiver.

This is the Geiger’s Method, which is used to formulate the data kernel matrix.

Let us take a hypothetical earthquake location problem and solve it using this method.

MATLAB Codes:

clear; close all;
global G epsilon;
epsilon=1.0e-6;

% Velocity parameters
vpvs = 1.78;
vp=6.5;
vs=vp/vpvs;

% defining the region of the problem which contains both the source and the receiver
%We take a 100x100 units^2 area and the depth of 100 units. The surface is 0 units depth.
xmin=-100;
xmax=100;
ymin=-100;
ymax=100;
zmin=-100;
zmax=0;


% stations: x, y, z coordinates (xi, yi, 0)
sxm = [-9, -8, -7, -6, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9]'*ones(1,9)*10;
sym = 10*ones(9,1)*[-9, -8, -7, -6, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9];
sx = sxm(:);
sy = sym(:);
Ns = length(sx); % num of stations
sz = zmax*ones(Ns,1); %zeros

% true earthquakes locations
Ne = 20; % number of earthquakes
M = 4*Ne; % 4 model parameters, x, y, z, and t0, per earthquake
extrue = random('uniform',xmin,xmax,Ne,1); %x0
eytrue = random('uniform',ymin,ymax,Ne,1); %y0
eztrue = random('uniform',zmin,zmax,Ne,1); %z0
t0true = random('uniform',0,0.2,Ne,1); %t0
mtrue = [extrue', eytrue', eztrue', t0true']'; %column matrix containing all the true model parameter values

Nr = Ne*Ns; % number of rays, that is, earthquake-stations pairs
N = 2*Ne*Ns; % total data; Ne*Ns for P and S each

% Generating the true data based on the true earthquake and station locations
dtrue=zeros(N,1); %allocating space for all data
for i = [1:Ns] % loop over stations
for j = [1:Ne] % loop over earthquakes
 dx = mtrue(j)-sx(i); % x component of the displacement obtained by the difference between station and source location x component
 dy = mtrue(Ne+j)-sy(i); % y component of the source- receiver displacement
 dz = mtrue(2*Ne+j)-sz(i); % z component of displacement between source and receiver
 r = sqrt( dx^2 + dy^2 + dz^2 ); % source-receiver distance
 k=(i-1)*Ne+j;
 dtrue(k)=r/vp+mtrue(3*Ne+j); %P arrival time for each station-source pair obtained by summing the travel time with the origin time of the earthquake
 dtrue(Nr+k)=r/vs+mtrue(3*Ne+j); % S arrival time
end
end
 
%  Generating observed data by adding gaussian noise with standard deviation 0.2 
sd = 0.2;
dobs=dtrue+random('normal',0,sd,N,1); % observed data

%% Determining the predicted arrival time using the Geiger's Method
% inital guess of earthquake locations
mest = [random('uniform',xmin,xmax,1,Ne), random('uniform',ymin,ymax,1,Ne), ...
 random('uniform',zmin+2,zmax-2,1,Ne), random('uniform',-0.1,0.1,1,Ne) ]';
% Here, we take a random initial guess


%Formulating the data kernel matrix and estimating the predicted models
for iter=[1:10] %for 10 iterations (termination criteria)
 
 % formulating data kernel
 G = spalloc(N,M,4*N); %N- total num of data,2*Ne*Ns; M is total num of model
 %parameters, 4*Ne
 dpre = zeros(N,1); %allocating space for predicted data matrix
 for i = 1:Ns % loop over stations
 for j = 1:Ne % loop over earthquakes
 dx = mest(j)-sx(i); % x- component of displacement obtained using the initial guess
 dy = mest(Ne+j)-sy(i); % y- component of displacement
 dz = mest(2*Ne+j)-sz(i); % z- component of displacement
 r = sqrt( dx^2 + dy^2 + dz^2 ); % source-receiver distance for each iteration
 k=(i-1)*Ne+j; %index for each ray
 dpre(k)=r/vp+mest(3*Ne+j); %predicted P wave arrival time
 dpre(Nr+k)=r/vs+mest(3*Ne+j); %predicted S wave arrival time
%First half of data kernel matrix correspoding to P wave
 G(k,j) = dx/(r*vp); % first column of data kernel matrix
 G(k,Ne+j) = dy/(r*vp); % second column of data kernel matrix
 G(k,2*Ne+j) = dz/(r*vp); % third column of data kernel matrix
 G(k,3*Ne+j) = 1; % fourth column of data kernel matrix
% Second half of the data kernel matrix corresponding to S wave
 G(Nr+k,j) = dx/(r*vs);
 G(Nr+k,Ne+j) = dy/(r*vs);
 G(Nr+k,2*Ne+j) = dz/(r*vs);
 G(Nr+k,3*Ne+j) = 1;
 end
 end
 
 % solve with dampled least squares
 dd = dobs-dpre;
 dm=bicg(@dlsfun,G'*dd,1e-5,3*M); solving using the biconjugate method
%solving the damped least square equation G'dd = [ G'G + epsilon* I] dm
% We use biconjugate method to reduce the computational cost (see for the dlsfun at the bottom)
 mest = mest+dm; %updated model parameter
 
end

% Generating the final predicted data
dpre=zeros(N,1);
for i = 1:Ns % loop over stations
for j = 1:Ne % loop over earthquakes
 dx = mest(j)-sx(i);
 dy = mest(Ne+j)-sy(i);
 dz = mest(2*Ne+j)-sz(i);
 r = sqrt( dx^2 + dy^2 + dz^2 );
 k=(i-1)*Ne+j;
 dpre(k)=r/vp+mest(3*Ne+j); % S-wave arrival time
 dpre(Nr+k)=r/vs+mest(3*Ne+j); % P- wave arriavl time
end
end

% Calculating the data and model misfit
expre = mest(1:Ne); % x0
eypre = mest(Ne+1:2*Ne); %y0
ezpre = mest(2*Ne+1:3*Ne); %z0
t0pre = mest(3*Ne+1:4*Ne); %t0
dd = dobs-dpre; %residual of observed and predicted arrival time
E = dd'*dd; %error
fprintf('RMS traveltime error: %f\n', sqrt(E/N) );
Emx = (extrue-expre)'*(extrue-expre); %misfit for x0
Emy = (eytrue-eypre)'*(eytrue-eypre); %misfit for y0
Emz = (eztrue-ezpre)'*(eztrue-ezpre); %misfit for z0
Emt = (t0true-t0pre)'*(t0true-t0pre); %misfit for t0
fprintf('RMS model misfit: x %f y %f z %f t0 %f\n', sqrt(Emx/Ne), sqrt(Emy/Ne), sqrt(Emz/Ne), sqrt(Emt/Ne) );
EQloc2.jpg
Earthquake location Example (William Menke, 2012) : The black triangles are the station locations, green circles are true earthquake location and the blue circles are predicted earthquake location.

dlsfun.m:

function y = dlsfun(v,transp_flag)
global G epsilon;
temp = G*v;
y = epsilon * v + G'*temp;
return

Simple 1D velocity model inversion from P arrival time

Refer to Chapter 5, Introduction to Seismology, Shearer 2009.

Problem:

From the P-wave travel time data below (note that the reduction velocity of 8km/s), inverse for the 1D velocity model using T(X) curve fitting (fit the T(X) curve with lines, then invert for the ray parameter p and delay time τ(p), then solve for the velocity and depth).

 

screenshot-from-2016-12-27-16-58-38
Fig 1. P-wave travel time from a seismic experiment.

 

Solution:

We construct the T(X) plot from the data: y = T  – X/8 => T = y + X/8 and then divide the T(X) curve into fitting lines (Fig. 2).

ex5-1fit
Fig 2. T(X) plot with lines fitting to data points.

The ray parameter and delay is calculated for each line by using:

p = dT/dX

τ = T – pX

These parameters p and τ are slope and intercept of a line, respectively.

The τ(p) can be written as:

Screenshot from 2016-12-27 17-14-52.png

as in the matrix form of:

Screenshot from 2016-12-27 17-14-57.png

or, in short,

We can invert for the thickness of each layer using least square inversion:

If we assume ray slowness u(i) = p(i) for each line, we can solve the h and v for the 1D velocity model.

Here is the result of the inversion:

ex5.1v_model.png

The sharp increase in Vp from ~6.6 km/s to 7.7 km/s at ~ 30km depth suggest for the Moho discontinuity.

The codes for the problem is written in python and can be downloaded: source codeτ inversion code, data points.

Nguyen Cong Nghia, IESAS

 

Travel-time curve calculation for the Spherical Earth (PREM Model)

I have used the Preliminary Reference Earth Model of Dziewonski and Anderson (1981) to calculate the travel time curve. This model is very robust in the sense that it is designed to fit a wide variety of data including free-oscillation, surface wave dispersion observations, , travel time data for a number of body wave phases etc. To download the model, click here.

The travel time curve is the plot of first arrival time of a ray from source to receiver vs the epicentral distance (surface to surface) traveled. The sum of the time taken by the ray in each layer from the source to the receiver gives the total travel time. The ray path depends on the ray-parameter (which is constant for a laterally homogeneous, flat layers). Ray with different ray-parameter will propagate along different paths for a given earth model. For a normal velocity model (velocity increases with depth), smaller ray parameters will turn deeper in Earth and hence travel farther. The ray-parameter (p) varies along the travel time curve and is give by its slope at an instant.

I have taken two approach to calculate the travel time; first by modeling the earth using laterally homogeneous layers, and second by modeling using the velocity gradient. For details see Chapter 4, Shearer, 2009.

In the first case, I stack all the layers and the distance traveled by the ray is calculated by summing the travel time in each layer until the ray turns (which is decided by comparing slowness and ray-parameter value). When the slowness (reciprocal of velocity) of the layer is equal to the ray-parameter, the ray turns in that medium and then start traveling towards the surface. There is an inherent uncertainty associated with this method depending on the thickness and large number of the layer.

In different approach, we can parametrize the velocity model at a number of discrete points in depth and evaluate the distance and travel time by assuming an appropriate interpolation function between the model points. Here, for my calculation, I have considered the linear velocity gradient in the layer.

We have used the Earth-flattening transformation to calculate the travel-time for the spherical Earth. The curvature of the travel time curve in a spherical Earth can be simulated in a flat Earth model if a special velocity gradient is introduced in the half-space (Mu¨ller, 1971; Shearer, 2009).

vel_model
The Earth-Flattening Transformation applied to the PREM Model

Fortran Program for layered earth model:

program comp_travel_time
 !Calculates the travel time data for a layered, flat earth for given layer thicknesses and velocities
 implicit none
 real, parameter :: pi = 2*asin(1.0), a=6371
 integer, parameter:: np=401
 real(kind=16), allocatable, dimension(:) :: delztmp, delz, slow, vel, eta, slowS, velS, etaS, utop, ubot, ustop, usbot
 real(kind=16) :: xx0, xx, tt0, tt, delz0, rr, pinit,pend, dp, xx0S, xxS, tt0S, ttS
 integer:: i,j,numl, nlyr
 real(kind=16), dimension(1:np) :: p
 character*11 :: file1 !name of the model file
 CHARACTER(len=32) :: arg1 !number of layers as argument
 CALL get_command_argument(1, arg1)
 read(arg1, '(I2)') numl
 allocate(delztmp(1:numl), delz(1:numl), slow(1:numl), vel(1:numl), eta(1:numl), &
 slowS(1:numl), velS(1:numl), etaS(1:numl),utop(1:numl), ubot(1:numl), ustop(1:numl), usbot(1:numl))

file1='ps_prem.dat' !model file name
 i=1
 open(1,file=file1, status='old')
 do while (i .le. numl)
 read(1,*) delztmp(i), vel(i), velS(i) !layer thickness and vel from file
 call read_data(delztmp(i), vel(i), velS(i))
 i=i+1
 end do
 close(1)
 i=1
 j=1
 do while(i .le. numl) !For flat Earth, calculate the layer thickness, and average velocity
 delz(j)=delztmp(i+1) - delztmp(i)
 utop(j)=vel(i)
 ustop(j)=velS(i)
 ubot(j)=vel(i+1)
 usbot(j)=velS(i+1)
 vel(j)=(utop(j)+ubot(j))/2
 velS(j)=(ustop(j)+usbot(j))/2
 !write(*,'(i2,3f12.2)') j, delz(j), vel(i), vel(i+1)
 nlyr=j
 i=i+1
 j=j+1
 end do
 nlyr=nlyr-1 !total number of layers

pinit=0.0017 !lower bound of ray parameter range
 pend=0.1128 !upper bound of ray parameter range
 dp=(pend-pinit)/np
 !write(*,'(a,1x,a,1x,a,1x,a)') "p in s/km","X(p) in km","X(p) in degrees","Time in mins"
 do j=1,np
 p(j)=pinit + (j-1)*dp !ray paramter increasing by dp 
 xx0=0
 xx0S=0
 tt0=0
 tt0S=0
 !P-wave calculation
 do i=1,nlyr
 slow(i)=1/vel(i) !in s/km
 if (slow(i) .gt. p(j) .and. slow(i) .le. 100000) then
 call travel_time_calc(xx0,tt0,slow(i),p(j),delz(i),eta(i),xx,tt)
 end if
 !S- wave calculation
 slowS(i)=1/velS(i)
 if (slowS(i) .gt. p(j) .and. slowS(i) .le. 100000) then
 call travel_time_calc(xx0S,tt0S,slowS(i),p(j),delz(i),etaS(i),xxS,ttS)
 end if
 end do
 
 open(2, file='p_travel_time.dat')
 write(2,'(f12.4,1X,f12.2,1X,f12.3,1X,f12.2)') p(j), xx, xx*(360/(2*pi*a)), tt/60
 !write(*,'(f12.4,1X,f12.2,1X,f12.3,1X,f12.2)') p(j), xx, xx*(360/(2*pi*a)), tt/60

open(3, file='s_travel_time.dat')
 write(3,'(f12.4,1X,f12.2,1X,f12.3,1X,f12.2)') p(j), xxS, xxS*(360/(2*pi*a)), ttS/60
 !write(*,'(f12.4,1X,f12.2,1X,f12.3,1X,f12.2)') p(j), xxS, xxS*(360/(2*pi*a)), ttS/60
 end do
 
 deallocate(delztmp, delz, slow, vel, eta, slowS, velS, etaS) 
end program comp_travel_time


!!SUBROUTINES


!reading the p, and s wave velocity, and depth from the model and convert it from spherical to flat earth
subroutine read_data(delztmp, vel, velS)
 real, parameter :: pi = 2*asin(1.0), a=6371
 real(kind=16):: delztmp, vel, velS
 real(kind=16) :: rr
 rr=a-delztmp !distance from the center of earth
 delztmp=-a*log(rr/a) !earth flattening transformation
 vel=(a/rr)*vel !earth flattening transformation for velocity
 velS=(a/rr)*velS !earth flattening transformation for shear velocity
 return
end

!Calculation of surface to surface distance of source and receiver and the travel time of ray
subroutine travel_time_calc(xx0,tt0,slow,p,delz,eta,xx,tt)
 real(kind=16) :: eta,slow,p,xx,xx0,delz,tt,tt0
 eta=sqrt(slow**2 - p**2)
 xx=xx0+(2*p*delz*(1/eta))
 tt=tt0+(2*slow**2*delz*(1/eta))
 xx0=xx
 tt0=tt
end

 

To compile and run this program:

datafile=ps_prem.dat
rm -f $datafile *.exe
awk 'NR>2{print $1,$3,$4}' PREM_MODELorg.model > $datafile
if [ ! -f "$datafile" ]; then #When the velocity model is not provided
cat >$datafile<<eof
3 4 3
6 6 5
9 8 7
eof
fi

numl=`wc -l $datafile | awk '{print $1}'`
prog="comp_travel_time.f95"
 
exect=`echo $prog | cut -d"." -f1` #execuatable

gfortran $prog -o ${exect}.exe 
if [ -f "${exect}.exe" ]; then
./${exect}.exe $numl #argument is the number of layers

travel_time.png

For the velocity gradient model, the fortran program is below. For this case, I have made use of the LAYERXT subroutine based on Chris Chapman’s WKBJ Program.

program comp_travel_time_vel_grad
 !Calculates the travel time data for a layered, flat earth for given layer thicknesses and velocities
 implicit none
 real, parameter :: pi = 2*asin(1.0), a=6371
 integer, parameter:: np=501
 real(kind=16), allocatable, dimension(:) :: delztmp, delz, vel, velS, utop, ubot, ustop, usbot
 real(kind=16) :: rr, pinit,pend, dp, delz0, xx, tt, xxS, ttS, xx0,tt0, xx0S, tt0S
 integer:: i,j,numl, nlyr,irtr, irtrS
 real(kind=16), dimension(1:np) :: p
 character*12 :: file1 !name of the model file
 CHARACTER(len=32) :: arg1 !number of layers as argument
 CALL get_command_argument(1, arg1)
 read(arg1, '(I2)') numl
 allocate(delztmp(1:numl), delz(1:numl), vel(1:numl), &
 utop(1:numl), ubot(1:numl), ustop(1:numl), usbot(1:numl), velS(1:numl))


 file1='ps_prem.dat' !model file name
 i=1
 open(1,file=file1, status='old')
 open(5,file='flat_earth.dat')
 do while (i .le. numl)
 read(1,*) delztmp(i), vel(i), velS(i) !layer thickness and vel from file
 rr=a-delztmp(i) !distance from the center of earth
 delztmp(i)=-a*log(rr/a) !earth flattening transformation
 vel(i)=(a/rr)*vel(i) !earth flattening transformation for velocity
 velS(i)=(a/rr)*velS(i) !earth flattening transformation for shear velocity
 write(5,'(i2,3f12.2)') i, delztmp(i), vel(i), velS(i)
 i=i+1
 end do
 close(1)
 i=1
 j=1
 do while(i .le. numl) !For flat Earth
 delz(j)=delztmp(i+1) - delztmp(i)
 utop(j)=1/vel(i)
 ustop(j)=1/velS(i)
 ubot(j)=1/vel(i+1)
 usbot(j)=1/velS(i+1)
 !write(*,'(i2,3f12.2)') j, delz(j), vel(i), vel(i+1)
 nlyr=j
 i=i+1
 j=j+1
 end do
 nlyr=nlyr-1 !because of singularity at earth's center
 write(*,*) "size_delz", nlyr
 pinit=0.0017 !lower bound of ray parameter range
 pend=0.1128 !upper bound of ray parameter range
 dp=(pend-pinit)/(np-1)
 do j=1,np
 p(j)=pinit + (j-1)*dp !ray paramter increasing by dp
 xx0=0
 xx0S=0
 tt0=0
 tt0S=0
 irtr=100
 irtrS=100
 do i=1,nlyr
 if (irtr .ne. 2) then
 call layerxt(p(j),delz(i),utop(i),ubot(i),xx,tt,irtr)
 xx=xx0+xx
 tt=tt0+tt
 xx0=xx
 tt0=tt
 else
 exit
 end if
 end do
 open(2, file='p_travel_time.dat')
 write(2,'(f12.4,1X,f12.2,1X,f12.3,1X,f12.2)') p(j), 2*xx, 2*xx*(360/(2*pi*a)), 2*tt/60
 !write(*,'(f12.4,1X,f12.2,1X,f12.3,1X,f12.2)') p(j), xx, xx*(360/(2*pi*a)), tt/60
 do i=1,nlyr
 if (ustop(i) .le. 100000 .and. usbot(i) .le. 100000) then
 if (irtrS .ne. 2) then
 call layerxt(p(j),delz(i),ustop(i),usbot(i),xxS,ttS,irtrS)
 xxS=xx0S+xxS
 ttS=tt0S+ttS
 xx0S=xxS
 tt0S=ttS
 else
 exit
 end if
 end if
 end do
 open(3, file='s_travel_time.dat')
 write(3,'(f12.4,1X,f12.2,1X,f12.3,1X,f12.2)') p(j), 2*xxS, 2*xxS*(360/(2*pi*a)), 2*ttS/60
 !write(*,'(f12.4,1X,f12.2,1X,f12.3,1X,f12.2)') p(j), xxS, xxS*(360/(2*pi*a)), ttS/60
 end do
 write(*,'(a,1x,a,1x,a,1x,a)') "p in s/km","X(p) in km","X(p) in degrees","Time in mins"
 deallocate(delztmp, delz, vel, &
 utop, ubot, ustop, usbot, velS) 
end program comp_travel_time_vel_grad

!!SUBROUTINES


! LAYERXT calculates dx and dt for a ray in a layer with a linear 
! velocity gradient. This is a highly modified version of a 
! subroutine in Chris Chapman's WKBJ program.
!
! Inputs: p = horizontal slowness
! h = layer thickness
! utop = slowness at top of layer
! ubot = slowness at bottom of layer
! Returns: dx = range offset
! dt = travel time
! irtr = return code
! = -1, zero thickness layer
! = 0, ray turned above layer
! = 1, ray passed through layer
! = 2, ray turned in layer, 1 leg counted in dx,dt
!
subroutine LAYERXT(p,h,utop,ubot,dx,dt,irtr)
 implicit none
 real(kind=16) :: p, h, utop, ubot,dx, dt, u1,u2,v1,v2,b,eta1,x1,tau1, &
 dtau,eta2, x2,tau2
 integer :: irtr
 if (p.ge.utop) then !ray turned above layer
 dx=0.
 dt=0.
 irtr=0
 return 
 else if (h.eq.0.) then !zero thickness layer
 dx=0.
 dt=0.
 irtr=-1
 return 
 end if

u1=utop
 u2=ubot
 v1=1./u1
 v2=1./u2
 b=(v2-v1)/h !slope of velocity gradient

eta1=sqrt(u1**2-p**2)

if (b.eq.0.) then !constant velocity layer
 dx=h*p/eta1
 dt=h*u1**2/eta1
 irtr=1
 return
 end if

x1=eta1/(u1*b*p)
 tau1=(log((u1+eta1)/p)-eta1/u1)/b

if (p.ge.ubot) then !ray turned within layer,
 dx=x1 !no contribution to integral
 dtau=tau1 !from bottom point
 dt=dtau+p*dx
 irtr=2
 return
 end if

irtr=1

eta2=sqrt(u2**2-p**2)
 x2=eta2/(u2*b*p)
 tau2=(log((u2+eta2)/p)-eta2/u2)/b

dx=x1-x2
 dtau=tau1-tau2

dt=dtau+p*dx

return
end

travel_time.png

In the travel-time curve, we can notice the general increase of travel time with delta. We can notice the large shadow zone around 100 degrees. This is because of the P wave transition from the solid mantle to the liquid core; which leads to sharp decrease in velocity. This is the low-velocity zone, where the rays are bent downward. So, no rays originating at the surface can turn within the LVZ. If the seismic rays originate inside the LVZ, it can never get trapped within the zone. If the seismic attenuation is lower in the LVZ layer, then the wave originated in it can propagate a long distance (Optical fiber cables work on the same principle).

xp_curvep
X(p) curve for P-wave

In general, the epicentral distance increases as the p decreases. When the X increases with decreasing p, the branch of the travel time curve is prograde otherwise called retrograde. The retrograde branch can be observed occasionally because of the rapid velocity transition in the Earth. The point of change from prograde to retrograde is called caustics; where X value does not change with p. At this point, the energy is focussed since the rays with different p values (or take off angles) arrive at same range(X). We can also notice in the above X(p) curve that there is sharp change in X around 0.02 and 0.04. Around 0.02, the slope is positive, which could be because of the rapid velocity transition in Earth. Since, for lower p value, the ray travels deeper; this could be the inner core boundary (ICB). This is because at ICB, the P-wave velocity increases drastically.  At point 0.04, the slope is large negative; which could be the core-mantle boundary.

  • Utpal Kumar (IES, Academia Sinica)

Ray tracing through a 1-D velocity model

Refer to Chapter 4 of Shearer, Introduction to Seismology.

screenshot-from-2016-12-21-16-39-44

For a ray piercing through Earth, the ray parameter (or horizontal slowness) p is defined by several expressions:

Screenshot from 2016-12-21 16-38-07.png

where u = 1/v is the slowness, θ is the ray incidence angle, T is the travel time, X is the horizontal range and utp is the slowness at the ray turning point.

Screenshot from 2016-12-21 16-37-05.png

The vertical slowness is defined as:

screenshot-from-2016-12-21-16-40-40

and integral expressions for the surface-to-surface travel time are:

Screenshot from 2016-12-21 16-41-26.png

Screenshot from 2016-12-21 16-41-32.png

With these equations we can calculate the travel time (T) and horizontal distance (X) for a given ray with ray parameter p and velocity model v(z).

Apply these to a problem (Exercise 4.8 Shearer):

(COMPUTER) Consider MARMOD, a velocity-versus-depth model, which is typical of much of the oceanic crust (Table 4.1). Linear velocity gradients are assumed to exist at intermediate depths in the model; for example, the P velocity at 3.75 km is 6.9 km/s. Write a computer program to trace rays through this model and produce a P-wave T(X) curve, using 100 values of the ray parameter p equally spaced between 0.1236 and 0.2217 s/km. You will find it helpful to use subroutine LAYERXT (provided in Fortran in Appendix D and in the supplemental web material as a Matlab script), which gives dx and dt as a function of p for layers with linear velocity gradients. Your program will involve an outer loop over ray parameter and an inner loop over depth in the model. For each ray, set x and t to zero and then, starting with the surface layer and proceeding downward, sum the contributions, dx and dt, from LAYERXT for each layer until the ray turns. This will give x and t for the ray from the surface to the turning point. Multiply by two to obtain the total surface-to-surface values of X(p) and T(p). Now produce plots of: (a) T(X) plotted with a reduction velocity of 8 km/s, (b) X(p), and (c) τ(p). On each plot, label the prograde and retrograde branches. Where might one anticipate that the largest amplitudes will occur?

Screenshot from 2016-12-21 21-32-25.png

Using the LAYERXT subroutine (Appendix D of Shearer), the FORTRAN and MATLAB can be downloaded from here. In this example, we will use python to calculate and plot the result.

ex4-8-1ex4-8-2ex4-8-3

 

 

The main program is as follow:

from ttcal import layerxt
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

v_model = pd.read_csv('v_model.dat')

p1 = 0.1236
p2 = 0.2217
n = 100

result = [['p','X','T']]
for p in np.linspace(p1, p2, n):
    X = 0; T = 0;
    for i in range(0, len(v_model)-1):
            p = p
            h = v_model.h[i+1] - v_model.h[i]
            utop = 1/v_model.vp[i]
            ubot = 1/v_model.vp[i+1]
            dx, dt, irtr = layerxt(p,h,utop,ubot)
            X += dx; T += dt
        #If the ray reflected at layer the dx, dt calculation stops
            if irtr == 2:
                    break
        result.append([p,X,T])

result_df = pd.DataFrame(result, columns=result.pop(0))
result_df = result_df[::-1]

#Multiply by 2 to get total surface-to-surface value ot X(p) and T(p)
result_df['T'] = 2*result_df['T']
result_df['X'] = 2*result_df['X']

plt.figure(1)
result_df['RT'] = result_df['T']-result_df['X']/8
result_df.plot(x='X', y='RT', legend=False)
plt.title('Reduction time')
plt.ylabel('Reduce time T - X/8 (s)')
plt.xlabel('Distance (km)')
plt.text(20,0.80,'Prograde', rotation =40)
plt.text(40,1.6,'Retrograde', rotation =35)
plt.text(80,1.00,'Prograde', rotation =0)
plt.savefig('ex4.8.1.png')

plt.figure(2)
result_df.plot(x='p', y='X', legend=False)
plt.title('X(p)')
plt.ylabel('Distance X (km)')
plt.xlabel('Ray parameter (km/s)')
plt.xlim([p1,p2])
plt.savefig('ex4.8.2.png')

plt.figure(3)
result_df['Taup'] = result_df['T'] - result_df['p']*result_df['X']
result_df.plot(x='p', y='Taup', legend=False)
plt.title('Tau(p)')
plt.ylabel('Tau(p) (s)')
plt.xlabel('Ray parameter (km/s)')
plt.xlim([p1,p2])
plt.savefig('ex4.8.3.png')

plt.figure(4)
v_model.h = v_model.h * (-1)
v_model.plot(x='vp', y='h', legend=False)
plt.ylabel('Depth (km)')
plt.xlabel('Velocity (km/s)')
plt.xlim([4,9])
plt.ylim([-10,0])
plt.savefig('ex4.8.4.png')

#plt.show()

LAYERXT subroutine written in python:

from math import sqrt, log


def layerxt(p,h,utop,ubot):
    ''' LAYERXT calculates dx and dt for ray in layer with linear velocity gradient
 
        Inputs:   p     =  horizontal slowness
             h     =  layer thickness
             utop  =  slowness at top of layer
             ubot  =  slowness at bottom of layer
        Returns:  dx    =  range offset
             dt    =  travel time
             irtr  =  return code
                   = -1,  zero thickness layer
                   =  0,  ray turned above layer
                   =  1,  ray passed through layer
                   =  2,  ray turned within layer, 1 segment counted in dx,dt'''
# ray turned above layer
    if p >= utop:
        dx = 0
        dt = 0
        irtr = 0
        return dx,dt,irtr

#Zero thickness layer
    elif h == 0:
        dx = 0
        dt = 0
        irtr = -1
        return dx,dt,irtr

#Calculate some parameters
    u1 = utop
    u2 = ubot
    v1 = 1.0/u1
    v2 = 1.0/u2
#Slope of velocity gradient
    b = (v2 - v1)/h
    eta1 = sqrt(u1**2 - p**2)

#Constant velocity layer
    if b == 0:
        dx = h*p/eta1
        dt = h*u1**2/eta1
        irtr = 1
        return dx,dt,irtr

    x1 = eta1/(u1*b*p)
    tau1=(log((u1+eta1)/p)-eta1/u1)/b

#Ray turned within layer, no contribution to integral from bottom point
    if p >= ubot:
        dx = x1
        dtau = tau1
        dt = dtau + p*dx
        irtr = 2    
        return dx,dt,irtr

        irtr=1
    eta2=sqrt(u2**2-p**2)
    x2=eta2/(u2*b*p)
    tau2=(log((u2+eta2)/p)-eta2/u2)/b

    dx=x1-x2
    dtau=tau1-tau2

    dt=dtau+p*dx
    return dx,dt,irtr

def flat(zsph,vsph):
    ''' Calculate the flatten transformation from a spherical Earth.
        zf, vf = flat(zsph,vsph)'''
    
# Radius of the Earth
    a = 6371.0

    rsph = a - zsph
    zf = -a*log(rsph/float(a))
    vf = (a/float(rsph))*vsph
    return zf, vf

Input velocity model:

h,vp,vs,den
0.0,4.5,2.4,2
1.5,6.8,3.75,2.8
6.0,7.0,3.5,2.9
6.5,8.0,4.6,3.1
10.0,8.1,4.7,3.1

The codes can be downloaded here: main program, python LAYERXT subroutine, velocity model.

You need to install numpy, pandas and matplotlib package to run this. Normally it can be installed in the terminal by (presume python has been installed):

pip install numpy pandas matplotlib

Nguyen Cong Nghia – IESAS

 

 

 

Normal Mode Summation for calculating the synthetic seismogram for a string

When we pluck a string fixed at both ends, then this will creating a standing wave. We can get some insights on the behavior of the propagating wave by considering normal modes or free oscillations of the string.

A wave propagating along the string follows the 1D scalar wave equation.

Since, the string has fixed ends, it satisfies the boundary conditions (zero displacement at the fixed ends) and this lead to the vibration of the string at only specific frequencies called eigenfrequencies. The eigenfrequencies are separated by pi*v/L, where v is the velocity of the wave and L is the length of the string. The velocity of the wave depends on the material of the string. Since, the eigenfrequencies are inversely proportional to the length of the string, the longer the string, closer the eigenfrequencies become.

Each eigenfrequencies correspond to a solution of space_term*time_term where the spatial term is known as the spatial eigenfunction.

A traveling wave is the weighted sum of the string’s normal modes. Hence, it is the sum of the eigenfunctions times its amplitude (function of the eigenfrequency). Thus, the real propagating wave is the interference of different modes of the string.

The amplitude for each eigenfunction depends on the position of the source that generated the waves and on the behaviour of the source as a function of time.

We can consider a very simple function to model the behavior of the source as a function of time

F(wn)=exp[-(wn*tau)^2/4]    (Stein and Wysession)

The shape of the source depends on the value of tau.

Let us make a program to calculate a synthetic seismogram using the normal mode summation (see Stein and Wysession (2003) pg. 466). We first keep the parameters same as the Stein and Wysession.

We consider a string of length 1m with a wavespeed 1m/s, a source at 0.2 m and the receiver at 0.7 m. Here, we add 200 modes to get the synthetic seismogram. The seismogram length is 1.25 sec and the time steps is 100.

program synt_seis_strings
    implicit none
    real, parameter :: pi = 4*atan(1.0)
    real :: stln,vel,xsrc,xrcvr,tseis,dt,tau,npiL,srctm,rcvrtm,wn,Fwn,spctm,coswt
    integer :: i,n,j
    integer, parameter :: nmode=200    !number of modes (nmode)
    integer, parameter :: nt=100    !number of time steps(nt)
    real(kind=16), dimension(1:nmode) :: Uxt
    real(kind=16), dimension(1:nt) :: Umode
    character*12 :: filename

    
    print *, 'Enter the string length (ex: 1.0):'
    read *, stln

    print *, 'Enter the wave speed: (ex: 1.0):'
    read *, vel
    

    print *, 'Enter the position of the source (in m) (ex: 0.2) :'
    read *, xsrc
    
    print *, 'Enter the position of the receiver (in m) (ex: 0.7) :'
    read *, xrcvr

    print *, 'Enter the duration of the seismogram (in sec) (ex: 1.25) :'
    read *, tseis

    dt=tseis/nt    !time step in sec
    tau=0.02    !source shape term

    
    do 5 i=1,nt
        Uxt(i)=0.0
    5 continue
    do 10 n=1,nmode
        npiL=n*pi/stln
        srctm=sin(npiL*xsrc)    !source term
        rcvrtm=sin(npiL*xrcvr)    !Receiver term
        wn=n*pi*vel/stln    !mode frequency
        Fwn=exp(-((tau*wn)**2)/4)    !weighting factor of different frequencies

        spctm=srctm*rcvrtm*Fwn        !space term
        
        do 15 j=1,nt
            coswt=cos(wn*dt*(j-1))    !time term
            Uxt(j)=Uxt(j) + coswt*spctm    !displacement
        15 continue
    
    10 continue
        open(unit=1,file='seism.dat')
        write(1,'(f6.2,f12.6)') (dt*(j-1),Uxt(j), j=1,nt)    !writing the displacement
        close(1)
        
        
    stop
end program synt_seis_strings

The synthetic seismogram data is saved in the file seism.dat. Now, we run the above program and plot the output.

#!/bin/bash

prog="synt_seis_strings.f" #fortran program name
exect=`echo $prog | cut -d"." -f1` #executable name

rm -f *.dat $exect #remove files

gfortran -ffree-form $prog -o $exect #compiling the fortran program
./$exect #running the executable

datascrpt='modesPlot.gp'
echo "plot [0:1.25][-8:8] 'seism.dat' using 1:2 w lines " >$datascrpt

#######GNUPLOT##################
gnuplot << EOF
fignm='synth_seis.ps'
figtitle='Synthetic Seismogram for a string'
datafile=system("echo $datascrpt")

xlb="Time"
ylb="Displacement"
set terminal postscript color solid #terminal type
set title figtitle #figure title
set xlabel xlb #xlabel
set ylabel ylb #ylabel
set output fignm #output fig name
unset key
load datafile
EOF

test1

In the above figure, we can see three peaks. First peak at 0.5s is the direct arrival from the source to the receiver. The distance between the source and receiver is 0.5. So the time taken from the source to the receiver is 0.5/1=0.5 sec.

The other two inverted peaks are the reflected arrivals. It is the time taken by the wave from the source to reach the end and reflected back to reach the receiver. The first reflected phase we observe is at 0.9s. It it due to the reflection from the left boundary. This can be calculated ((0.2-0)+0.7)/1=0.9 sec. Similarly, the other arrival (reflected phase from the right boundary) could also be calculated.

We have used the exponential source term (F(wn)=exp[-(wn*tau)^2/4] ) which depends on the value of eigenfrequencies and tau. We have modeled its dependency.

fwn_wnfwn_tauwn3-14

In the second figure, we have considered the value of wn=3.14.

Here, also we have considered the source term to be an exponential function. Instead, if we consider the source term to be abs(sin(x)) then we get the following synthetic seismogram.test1.png

Here, we can still see the peaks are 0.5,0.9 and 1.1 seconds but there are some wiggles are it. This is the effect of the source function. We have considered the same number of modes and other parameters.

When we consider the sinc function as the source function then the seismogram wiggles surrounding the peak decreases.

test1.png

Conclusion:

In this study, we have obtained the synthetic seismogram using the normal mode summation method. We now have two ways to think of the displacement of a traveling wave, propagating waves or normal modes. These two ways can give different interesting insights about the properties of the source and the medium. The difference between the two is that in the propagating wave formulation, the travel times are measured to infer the velocity of the medium and in the normal mode formulation, eigenfrequecies are measured to infer the velocity.

The normal mode solution gives us information about the medium in which the waves propagate and the source that generates them. The physical properties of the string control its velocity (or eigenfunctions and eigenfrequencies). The displacement obtained can be used to study the source that excited them.

In normal mode solution, we sum all the modes to obtain the synthetic seismogram. The individual modes, though, are not physically meaningful. Each mode actually start vibrating even before the wave reaches the end of the string. But when all the modes are summed, the resulting wave propagate at the correct velocity.

test1.png

In this study, we have considered a uniform string to calculate the synthetic seismogram. We could generalize this idea to obtain the modes of the non-uniform string. We can consider this case in the future posts. We can even calculate the modes for an sphere.

  • Utpal Kumar (IES, Academia Sinica)

Modeling a wave on a string using finite differences

Based on problem 3.7 chapter 3 of Introduction to Seismology (Shearer)

(COMPUTER) In the case of plane-wave propagation in the x direction within a uniform medium, the homogeneous momentum equation (3.9) for shear waves can be expressed as

screenshot-from-2016-12-14-15-42-21

,where u is the displacement. Write a computer program that uses finite differences to solve this equation for a bar 100 km in length, assuming β = 4 km/s. Use dx = 1 km for the length spacing and dt = 0.1 s for the time spacing. Assume a source-time function  at u (50 km) of the form:

screenshot-from-2016-12-14-15-42-16

Apply a stress-free boundary condition at u(0 km) and a fixed boundary condition at u (100 km).  Approximate the second derivatives using the finite difference scheme:

screenshot-from-2016-12-14-15-42-25

Plot u(x) at 4 s intervals from 1 to 33 s. Verify that the pulses travel at velocities of 4 km/s. What happens to the reflected pulse at each endpoint? What happens when the pulses cross?

Solutions can be downloaded here.

!**********************************************************************
!* *
!* 1-D Plane wave propagation with fixed and free boundaries * 
!* *
!**********************************************************************

!----initialzie t, dx, dt, tlen, beta, and u1, u2, u3 arrays----c

implicit none
 integer:: i, counter, nx, nx1
 integer:: it, ntmax
 real(8):: t, dt, dt2, dx, dx2, beta, beta2, tlen, tmax, rhs
 real*8, allocatable:: u1(:), u2(:), u3(:)
 character(10):: data
 character(4):: nf
!----parameters ----c
! dt must be less than (dx/beta) for numerical stability

nx = 100 ! number of grid
 dt = 0.10d0 ! time interval
 dx = 1.0d0 ! grid interval
 beta = 4.0d0 ! wave velocity
 tlen = 5.0d0 ! time length of wave source
 ntmax = 1000 ! maximum time step (finish at it = ntmax, or)
 tmax = 33.0d0 ! maximum calculation time (finish at t< tmax)

!---- allocate dimension variables ----c

nx1 = nx + 1
 allocate (u1(nx1), u2(nx1), u3(nx1))

!---- initialize t, counter, u1, u2, u3 ----c

t = 0.0d0
 counter = 1
 
 do i = 1, nx1
 u1(i) = 0.0d0
 u2(i) = 0.0d0
 u3(i) = 0.0d0
 end do

!---- calculate c^2, dt^2, dx^2 ----c
 beta2 = beta**2
 dt2 = dt**2
 dx2 = dx**2

! ============= Time marching loop ===========
 
 do it = 1, ntmax

t = t + dt

!---- calculate u3(i) ----c

do i= 2, nx
 rhs=beta2*2*(u2(i+1)-2.0d0*u2(i)+u2(i-1))/dx2
 u3(i)=dt2*rhs+2.0d0*u2(i)-u1(i)
 end do

!---- free boundary (du/dx=0) ----c
 u3(1)=u3(2)
!---- fixed boundary (u=0) ----c
 u3(nx1)=0.0d0
!---- source time function c----c

if (t.le.tlen) then
 u3(51) = sin(3.1415927*t/tlen)**2
 end if

!---- change new and old variables ----c

do i=1,nx1
 u1(i) = u2(i)
 u2(i) = u3(i)
 end do
!---- make data file ----c
 write(nf,'(i3.3)') counter
 data = "data"//nf

open(7,file=data,status='replace')
 do i=1, nx1
 write(7,*) i, u2(i)
 end do
!---- output u2 at desired intervals, stop when t is big enough

write(6,'(a5,1x,i3,1x,f6.3)') 'loop', counter, t
 
 if (t.gt.tmax) exit

counter = counter + 1

end do

stop
 end

The following is used to plot the results:

#!/usr/bin/bash
gnuplot << eof
set term gif animate
set output "animate.gif"
n=9 #n frames
set xrange [0:102]
set yrange [-2:2]
i=0
list = system('ls data*')
do for [file in list] {
 plot file w lines lt 1 lw 1.5 t sprintf("t=%i sec",i/10)
 i = i + 1
}
set output
eof

animate

Nguyen Cong Nghia – IESAS