Skip to content
This repository has been archived by the owner on Dec 22, 2021. It is now read-only.

Commit

Permalink
Test jacobians in eqlayer
Browse files Browse the repository at this point in the history
  • Loading branch information
leouieda committed Jun 24, 2016
1 parent b973cf2 commit 3c4ad38
Show file tree
Hide file tree
Showing 2 changed files with 32 additions and 2 deletions.
2 changes: 1 addition & 1 deletion fatiando/gravmag/eqlayer.py
Original file line number Diff line number Diff line change
Expand Up @@ -46,7 +46,7 @@ class EQLBase(Misfit):
"""

def __init__(self, x, y, z, data, grid):
super().__init__(data=data, nparams=grid.size, islinear=True)
super().__init__(data=data, nparams=len(grid), islinear=True)
self.x = x
self.y = y
self.z = z
Expand Down
32 changes: 31 additions & 1 deletion fatiando/tests/gravmag/test_eqlayer.py
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,37 @@
from fatiando import utils, gridder


def test_eql_grav_jacobian():
"EQLGravity produces the right Jacobian matrix for single source"
model = PointGrid([-10, 10, -10, 10], 500, (2, 2))[0]
model.addprop('density', 1)
n = 1000
x, y, z = gridder.scatter([-10, 10, -10, 10], n, z=-100, seed=42)
data = sphere.gz(x, y, z, [model])

eql = EQLGravity(x, y, z, data, [model])
A = eql.jacobian(None)

assert A.shape == (n, 1)
assert_allclose(A[:, 0], data, rtol=0.01)


def test_eql_mag_jacobian():
"EQLTotalField produces the right Jacobian matrix for single source"
inc, dec = -30, 20
model = PointGrid([-10, 10, -10, 10], 500, (2, 2))[0]
model.addprop('magnetization', utils.ang2vec(1, inc, dec))
n = 1000
x, y, z = gridder.scatter([-10, 10, -10, 10], n, z=-100, seed=42)
data = sphere.tf(x, y, z, [model], inc, dec)

eql = EQLTotalField(x, y, z, data, inc, dec, [model])
A = eql.jacobian(None)

assert A.shape == (n, 1)
assert_allclose(A[:, 0], data, rtol=0.01)


def test_eqlgrav_prism_interp():
"EQLGravity can interpolate data from a prism"
model = [Prism(-300, 300, -500, 500, 100, 600, {'density': 400})]
Expand Down Expand Up @@ -56,4 +87,3 @@ def test_eqlayer_polereduce():
calc = sphere.tf(x, y, z, layer, inc=-90, dec=0)

assert_allclose(calc, true, atol=10, rtol=0.05)

0 comments on commit 3c4ad38

Please sign in to comment.