next up previous Homework 2 [pdf]

Next: Completing the assignment Up: Homework 2 Previous: Computational part 2

Programming part 2

time
Figure 4. Traveltime in a $V(z)$ medium. 		time
[pdf] [png] [scons]

The program cmp/traveltime.c computes reflection traveltimes in a $V(z)$ medium by using different methods.

  1. Modify the program to implement approximation (18) using your equation (19).
  2. Modify the program to implement exact traveltime computation by doing shooting iterations with equations (20-21). Using Newton's method, you can find the value of $p$ for a given $h$ by solving the non-linear equation $h(p)=h$ with iterations
    \begin{displaymath} p_{n+1} = p_n - \frac{h(p_n)-h}{h'(p_n)}\;. \end{displaymath} (23)

  3. For the traveltime in Figure 4, find the offset, where the absolute error of approximation (17) exceeds 0.05 s.
  4. For the traveltime in Figure 4, find the offset, where the absolute error of approximation (18) exceeds 0.05 s. 

/* Compute traveltime in a V(z) model. */
#include <rsf.h>

int main(int argc, char* argv[])
{
    char *type;
    int ih, nh, it, nt, ir, nr, *r, iter, niter;
    float h, dh, h0, dt, t0, t2, h2, v2, s, p, hp, tp;
    float *v, *t;
    sf_file vel, tim;

    /* initialize */
    sf_init(argc,argv);

    /* input and output */
    vel = sf_input("in");
    tim = sf_output("out");
    
    /* time axis from input */
    if (!sf_histint(vel,"n1",&nt)) sf_error("No n1=");
    if (!sf_histfloat(vel,"d1",&dt)) sf_error("No d1=");

    /* offset axis from command line */
    if (!sf_getint("nh",&nh)) nh=1;
    /* number of offsets */
    if (!sf_getfloat("dh",&dh)) dh=0.01;
    /* offset sampling */
    if (!sf_getfloat("h0",&h0)) h0=0.0;
    /* first offset */

    /* get reflectors */
    if (!sf_getint("nr",&nr)) nr=1;
    /* number of reflectors */
    r = sf_intalloc(nr);
    if (!sf_getints("r",r,nr)) sf_error("Need r=");

    if (NULL == (type = sf_getstring("type")))
	type = "hyperbolic";
    /* traveltime computation type */

    if (!sf_getint("niter",&niter)) niter=10;
    /* maximum number of shooting iterations */

    /* put dimensions in output */
    sf_putint(tim,"n1",nh);
    sf_putfloat(tim,"d1",dh);
    sf_putfloat(tim,"o1",h0);
    sf_putint(tim,"n2",nr);

    /* read velocity */
    v = sf_floatalloc(nt);
    sf_floatread(v,nt,vel);
    /* convert to velocity squared */
    for (it=0; it < nt; it++) {
	v[it] *= v[it];
    }

    t = sf_floatalloc(nh);

    for (ir=0; ir<nr; ir++) {
	nt = r[ir];
	t0 = nt*dt; /* zero-offset time */
	t2 = t0*t0;

	p = 0.0;

	for (ih=0; ih<nh; ih++) {
	    h = h0+ih*dh; /* offset */
	    h2 = h*h; 

	    switch(type[0]) {
		case 'h': /* hyperbolic approximation */
		    v2 = 0.0;
		    for (it=0; it < nt; it++) {
			v2 += v[it];
		    }
		    v2 /= nt;

		    t[ih] = sqrtf(t2+h2/v2);
		    break;
		case 's': /* shifted hyperbola */

		    /* !!! MODIFY BELOW !!! */

		    s = 0.0;

		    v2 = 0.0;
		    for (it=0; it < nt; it++) {
			v2 += v[it];
		    }
		    v2 /= nt;

		    t[ih] = sqrtf(t2+h2/v2);
		    break;
		case 'e': /* exact */
		    
		    /* !!! MODIFY BELOW !!! */

		    for (iter=0; iter < niter; iter++) {
			hp = 0.0;
			for (it=0; it < nt; it++) {
			    v2 = v[it];
			    hp += v2/sqrtf(1.0-p*p*v2);
			}
			hp *= p*dt;

			/* !!! SOLVE h(p)=h !!! */
		    }

		    tp = 0.0;
		    for (it=0; it < nt; it++) {
			v2 = v[it];
			tp += dt/sqrtf(1.0-p*p*v2);
		    }
		    
		    t[ih] = tp;
		    break;
		default:
		    sf_error("Unknown type");
		    break;
	    }
	}

	sf_floatwrite(t,nh,tim);
    }

    exit(0);
}


program Traveltime
! Compute traveltime in a V(z) model. 
use rsf

character(len=FSTRLEN) :: type
integer :: ih, nh, it, nt, ir, nr, iter, niter
real    :: h, dh, h0, dt, t0, t2, h2, v2, s, p, hp, tp
integer, allocatable, dimension (:) :: r
real,    allocatable, dimension (:) :: v, t
type (file) :: vel, tim

call sf_init() ! initialize Madagascar

! input and output 
vel = rsf_input("in")
tim = rsf_output("out")
    
! time axis from input 
call from_par(vel,"n1",nt)
call from_par(vel,"d1",dt)

! offset axis from command line 

call from_par("nh",nh,1)    ! number of offsets 
call from_par("dh",dh,0.01) ! offset sampling
call from_par("h0",h0,0.0)  ! first offset

! get reflectors 

call from_par("nr",nr,1) ! number of reflectors

allocate (r(nr))

call from_par("r",r)

call from_par("type",type,"hyperbolic")
! traveltime computation type 

call from_par("niter",niter,10)
! maximum number of shooting iterations 

! put dimensions in output 
call to_par(tim,"n1",nh)
call to_par(tim,"d1",dh)
call to_par(tim,"o1",h0)
call to_par(tim,"n2",nr)

! read velocity 
allocate (v(nt))
call rsf_read(vel,v)

! convert to velocity squared 
v = v*v

allocate(t(nh))

do ir=1, nr
   nt = r(ir)
   t0 = nt*dt ! zero-offset time 
   t2 = t0*t0

   p = 0.0;

   do ih=1, nh
      h = h0+(ih-1)*dh ! offset
      h2 = h*h

      select case (type(1:1))
      case ("h") ! hyperbolic approximation 
         v2 = 0.0
         do it=1, nt
            v2 = v2 + v(it)
         end do
         v2 = v2/nt

         t(ih) = sqrt(t2+h2/v2)
      case("s") ! shifted hyperbola 
         
!!! MODIFY BELOW !!! 

         s = 0.0
         v2 = 0.0
         do it=1, nt
            v2 = v2 + v(it)
         end do
         v2 = v2/nt

         t(ih) = sqrt(t2+h2/v2)
      case("e") ! exact 
		    
!!! MODIFY BELOW !!! 

         do iter=1, niter 
            hp = 0.0
            do it=1,nt
               v2 = v(it)
               hp = hp + v2/sqrt(1.0-p*p*v2)
            end do
            hp = hp*p*dt

!!! SOLVE h(p)=h !!! 
	
         end do

         tp = 0.0
         do it=1, nt
            v2 = v(it)
            tp = tp + dt/sqrt(1.0-p*p*v2)
         end do

         t(ih) = tp
      case default
         call sf_error("Unknown type")
      end select
   end do

   call rsf_write(tim,t)
end do

call exit(0)
end program Traveltime


#!/usr/bin/env python

import numpy
from math import sqrt
import m8r

# initialize parameters
par = m8r.Par()

# input and output
vel=m8r.Input()  
tim=m8r.Output() 

# time axis from input
nt = vel.int('n1')
dt = vel.float('d1')

# offset axis from command line
nh = par.int('nh',1)      # number of offsets
dh = par.float('dh',0.01) # offset sampling
h0 = par.float('h0',0.0)  # first offset

# get reflectors
nr = par.int('nr',1) # number of reflectors
r = par.ints('r',nr)

type = par.string('type','hyperbolic')
# traveltime computation type

niter = par.int('niter',10)
# maximum number of shooting iterations

# put dimensions in output
tim.put('n1',nh)
tim.put('d1',dh)
tim.put('o1',h0)
tim.put('n2',nr)

# read velocity
v = numpy.zeros(nt,'f')
vel.read(v)

# convert to velocity squared
v = v*v

t = numpy.zeros(nh,'f')

for ir in range(nr):
    nt = r[ir]
    t0 = nt*dt # zero-offset time 
    t2 = t0*t0

    p = 0.0

    for ih in range(nh):
        h = h0+ih*dh # offset
        h2 = h*h

        if type[0] == 'h':
            # hyperbolic approximation 
            v2 = numpy.sum(v)/nt
            t[ih] = sqrt(t2+h2/v2)

        elif type[0] == 's':
            # shifted hyperbola 
         
            ### MODIFY BELOW ###

            s = 0.0
            v2 = numpy.sum(v)/nt
            t[ih] = sqrt(t2+h2/v2)

        elif type[0] == 'e': 
            # exact 
	    
            ### MODIFY BELOW ###

            for iter in range(niter): 
                hp = 0.0
                for it in range(nt):
                    v2 = v[it]
                    hp += v2/sqrt(1.0-p*p*v2)
                hp *= p*dt

            ### SOLVE h(p)=h ###

            tp = 0.0
            for it in range(nt):
                v2 = v[it]
                tp += dt/sqrt(1.0-p*p*v2)
            t[ih] = tp
        else:
            raise RuntimeError('Unknown type')

    tim.write(t)

next up previous Homework 2 [pdf]

Next: Completing the assignment Up: Homework 2 Previous: Computational part 2

2019-09-26