Commit 9622e703 authored by Marius Kriegerowski's avatar Marius Kriegerowski
Browse files

fixed indexing

parent 00135715
......@@ -26,20 +26,24 @@ class VelocityShifter(Shifter):
class CakePhaseShifter(Shifter):
earthmodel_1d = Earthmodel1D.T()
tabulated_phases = List.T(gf.meta.TPDef.T())
timing = meta.Timing.T()
def get_vmin(self):
return min(filter(lambda x: x!=0., self.earthmodel_1d.profile('vs')))
def get_table(self, grid, receivers):
phases = [p.split(':')[1] for p in self.timing.phase_defs]
tabulated_phases = [gf.meta.TPDef(id=p, definition=p) for p in phases]
distances = grid.lateral_distances(receivers)
s_depths = grid._coords[2]
x_bounds = num.array(((grid.zmin, grid.zmax),
(min(num.ravel(distances)),
max(num.ravel(distances)))))
x_tolerance = num.array((grid.dz/2., grid.dx/2.))
x_tolerance = num.array((grid.dz/10., grid.dx/10.))
t_tolerance = grid.max_delta()/(self.get_vmin()*20.)
interpolated_tts = {}
for phase_def in self.tabulated_phases:
v_horizontal = phase_def.horizontal_velocities
for phase_def in tabulated_phases:
def evaluate(args):
source_depth, x = args
......@@ -48,13 +52,11 @@ class CakePhaseShifter(Shifter):
phases=phase_def.phases,
distances=[x*cake.m2d],
zstart=source_depth,
zstop=0.) # TODO: receiver depths
zstop=0.)
for ray in rays:
t.append(ray.t)
for v in v_horizontal:
t.append(x/(v*1000.))
if t:
return min(t)
else:
......@@ -63,7 +65,7 @@ class CakePhaseShifter(Shifter):
logger.info('prepare tabulated phases: %s' % phase_def.id)
sptree = spit.SPTree(
f=evaluate,
ftol=self.t_tolerance,
ftol=t_tolerance,
xbounds=x_bounds,
xtols=x_tolerance)
......@@ -71,8 +73,6 @@ class CakePhaseShifter(Shifter):
#sptree.dump(filename='sptree_%s.yaml' % phase_def.id)
n_x = len(distances)
arrivals = num.zeros((len(distances)*len(s_depths), len(receivers)))
def interpolate(phase_id):
......@@ -80,15 +80,15 @@ class CakePhaseShifter(Shifter):
for i_r, r in enumerate(receivers):
distances = grid.lateral_distances([r])
coords = num.array((num.tile(s_depths, len(distances)),
num.repeat(distances, len(s_depths)))).T
a = num.tile(num.ravel(distances), len(s_depths))
b = num.repeat(s_depths, len(distances))
coords = num.vstack((b, a)).T
arr = self.timing.evaluate(interpolate, coords)
arrivals[:, i_r] = arr
arrivals[:, i_r] = num.ravel(arr)
return arrivals
__all__ = [
'Shifter',
'VelocityShifter',
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment