Homework 2

Programming part 1

Instead of an analytical solution for a constant gradient of slowness squared, let us try an analytical solution for a constant gradient of velocity.

1. Change directory

   > cd hw2/eikonal
   

   • Uncomment the part of the SConstruct file that defines a velocity model with the constant velocity gradient.
   • Modify the program analytical.c to implement your equation (16).
   • Recompute the figures and check your results.

/* Analytical first-arrival traveltimes. */
#include <math.h>
#include <rsf.h>

int main (int argc, char* argv[])
{
    char *type;
    int n1, n2, i1, i2;
    float d1,d2, g1,g2, s,v0, x1,x2,gsq,g,s2,z,d;
    float *time;
    sf_file in, out;

    sf_init (argc,argv);
    in = sf_input("in");
    out = sf_output("out");

    /* Get grid dimensions */
    if (!sf_histint(in,"n1",&n1)) sf_error("No n1=");
    if (!sf_histint(in,"n2",&n2)) sf_error("No n1=");
    if (!sf_histfloat(in,"d1",&d1)) sf_error("No d1=");
    if (!sf_histfloat(in,"d2",&d2)) sf_error("No d2=");
    
    if (!sf_getfloat("g1",&g1)) g1 = 0.;    
    /* vertical gradient */
    if (!sf_getfloat("g2",&g2)) g2 = 0.;    
    /* horizontal gradient */
    gsq = g1*g1+g2*g2;
    g = sqrtf(gsq);

    if (!sf_getfloat("v0",&v0)) sf_error("Need v0=");
    /* initial velocity or slowness squared */

    if (!sf_getfloat("s",&s)) s = 0.0;
    /* shot location at the surface */
 
    if (NULL == (type = sf_getstring("case"))) type="c";
    /* case of velocity distribution */

    if (0.0 == g1 && 0.0 == g2) type="const";
    
    time = sf_floatalloc(n1);

    for (i2 = 0; i2 < n2; i2++) {
	x2 = i2*d2;
	for (i1 = 0; i1 < n1; i1++) {
	    x1 = i1*d1;
	    d = x1*x1+(x2-s)*(x2-s);
	    switch (type[0]) {
		case 's':
		    /* slowness squared */
		    s2 = v0+g1*x1+g2*x2;
		    z = 2.0*d/(s2+sqrtf(s2*s2-gsq*d));
		    time[i1] = (s2-gsq*z/6.0)*sqrtf(z);
		    break;
		case 'v': 
		    /* velocity */
		    s2 = 2.0*v0*(v0+g1*x1+g2*x2);
		    /* !!! CHANGE BELOW !!! */
		    time[i1] = hypotf(x2-s,x1)/v0;
		    break;
		case 'c': /* constant velocity */
		default:
		    time[i1] = hypotf(x2-s,x1)/v0;
		    break;
	    }
	}
	sf_floatwrite(time,n1,out);
    }

    exit (0);
}

program Analytical
! Analytical first-arrival traveltimes. */
use rsf

character(len=FSTRLEN) :: type
integer :: n1, n2, i1, i2
real    :: d1,d2, g1,g2, s,v0, x1,x2,gsq,g,s2,z,d
real, dimension (:), allocatable :: time
type (file) :: in, out

call sf_init() ! initialize Madagascar
in = rsf_input("in")
out = rsf_output("out")

! Get grid dimensions 
call from_par(in,"n1",n1)
call from_par(in,"n2",n2)
call from_par(in,"d1",d1)
call from_par(in,"d2",d2)

call from_par("g1",g1,0.) ! vertical gradient 
call from_par("g2",g2,0.) ! horizontal gradient 

gsq = g1*g1+g2*g2
g = sqrt(gsq)

call from_par("v0",v0)
! initial velocity or slowness squared 

call from_par("s",s,0.)
! shot location at the surface 

call from_par("type",type,"constant")
! case of velocity distribution 

if (0.0 == g1 .and. 0.0 == g2) type="const"
    
allocate (time(n1))

do i2 = 1, n2
   x2 = (i2-1)*d2
   do i1 = 1, n1
      x1 = (i1-1)*d1
      d = x1*x1+(x2-s)*(x2-s)
      
      select case (type(1:1))
      case("s") ! slowness squared 
         s2 = v0+g1*x1+g2*x2
         z = 2.0*d/(s2+sqrt(s2*s2-gsq*d))
         time(i1) = (s2-gsq*z/6.0)*sqrt(z)
      case("v") ! velocity 
         s2 = 2.0*v0*(v0+g1*x1+g2*x2)

!!! CHANGE BELOW !!! 
         time(i1) = hypot(x2-s,x1)/v0
      case("c") ! constant velocity 
      case default
         time(i1) = hypot(x2-s,x1)/v0
      end select
   end do
   call rsf_write(out,time)
end do

call exit (0)
end program Analytical

#!/usr/bin/env python

import numpy
from math import sqrt, hypot
import m8r

# initialize parameters
par = m8r.Par()

# input and output
inp=m8r.Input()  
out=m8r.Output()

# get grid dimensions
n1 = inp.int('n1')
n2 = inp.int('n2')
d1 = inp.float('d1')
d2 = inp.float('d2')

g1 = par.float('g1',0.0) # vertical gradient
g2 = par.float('g2',0.0) # horizontal gradient

gsq = g1*g1+g2*g2
g = sqrt(gsq)

v0 = par.float('v0')
# initial velocity or slowness squared

s = par.float('s',0.0)
# shot location at the surface

type = par.string('case','constant')
# case of velocity distribution

if 0.0 == g1 and 0.0 == g2:
    type='const'
    
time = numpy.zeros(n1,'f')

for i2 in range(n2):
   x2 = i2*d2
   for i1 in range(n1):
      x1 = i1*d1
      d = x1*x1+(x2-s)*(x2-s)
      
      if type[0] == 's':
      # slowness squared 
         s2 = v0+g1*x1+g2*x2
         z = 2.0*d/(s2+sqrt(s2*s2-gsq*d))
         time[i1] = (s2-gsq*z/6.0)*sqrt(z)

      elif type[0] == 'v':
      # velocity 
         s2 = 2.0*v0*(v0+g1*x1+g2*x2)

### CHANGE BELOW ### 
         time[i1] = hypot(x2-s,x1)/v0

      else:
         # constant velocity 
         time[i1] = hypot(x2-s,x1)/v0

   out.write(time)

Homework 2