Homework 2 |
Instead of an analytical solution for a constant gradient of slowness squared, let us try an analytical solution for a constant gradient of velocity.
> cd hw2/eikonal
/* 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 |