bugfixes: formula for gradients in multiple layers for vgrid corrected, fixed issue with integer division that would lead to a wrong spacing

This commit is contained in:
Marcel Paffrath 2015-11-04 10:04:05 +01:00
parent c5fc9ee89c
commit dc4d19ba88

View File

@ -449,9 +449,9 @@ class SeisArray(object):
rtop = Rbt[1] + R rtop = Rbt[1] + R
# need to determine the delta to add two cushion nodes around the min/max values # need to determine the delta to add two cushion nodes around the min/max values
thetaDelta = abs(thetaN - thetaS) / (nTheta - 1) thetaDelta = abs(thetaN - thetaS) / float((nTheta - 1))
phiDelta = abs(phiE - phiW) / (nPhi - 1) phiDelta = abs(phiE - phiW) / float((nPhi - 1))
rDelta = abs(rbot - rtop) / (nR - 1) rDelta = abs(rbot - rtop) / float((nR - 1))
# create a regular grid including +2 cushion nodes in every direction # create a regular grid including +2 cushion nodes in every direction
thetaGrid = np.linspace(thetaS - thetaDelta, thetaN + thetaDelta, num = nTheta + 2) # +2 cushion nodes thetaGrid = np.linspace(thetaS - thetaDelta, thetaN + thetaDelta, num = nTheta + 2) # +2 cushion nodes
@ -487,17 +487,19 @@ class SeisArray(object):
z = -(R + topo - radius) z = -(R + topo - radius)
if z > (topo + 1): if z > (topo + 1):
vel = 0.0 vel = 0.0
elif (topo + 1) >= z > (topo): # cushioning around topography elif (topo + 1) >= z > topo: # cushioning around topography
vel = vtop[0] vel = vtop[0]
else: else:
for index in range(nlayers): for index in range(nlayers):
if (topo + ztop[index]) >= z > (topo + zbot[index]): if (topo + ztop[index]) >= z > (topo + zbot[index]):
vel = (z - topo) / (zbot[index] - topo) * (vbot[index] - vtop[index]) + vtop[index] vel = (z - ztop[index]) / (zbot[index] - ztop[index]) * (vbot[index] - vtop[index]) + vtop[index]
break break
if not (topo + ztop[index]) >= z > (topo + zbot[index]): if not (topo + ztop[index]) >= z > (topo + zbot[index]):
print('ERROR in grid inputfile, could not find velocity for a z-value of %s in the inputfile'%(z - topo)) print('ERROR in grid inputfile, could not find velocity for a z-value of %s in the inputfile'%(z - topo))
return return
count += 1 count += 1
if vel < 0:
print('ERROR, vel <0; z, topo, zbot, vbot, vtop:', z, topo, zbot[index], vbot[index], vtop[index])
outfile.writelines('%10s %10s\n'%(vel, decm)) outfile.writelines('%10s %10s\n'%(vel, decm))
progress = float(count) / float(nTotal) * 100 progress = float(count) / float(nTotal) * 100