Homework 2 |
time
Figure 4. Traveltime in a medium. | |
---|---|
The program cmp/traveltime.c computes reflection traveltimes in a medium by using different methods.
/* 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) |
Homework 2 |