GEOPY-2357: Lost control on reference angles for MVI

simpeg_drivers/components/factories/directives_factory.py

@@ -332,9 +332,11 @@ def vector_inversion_directive(self):
         """Directive to update vector model."""
         if self._vector_inversion_directive is None and "vector" in self.factory_type:
             reference_angles = (
-                getattr(self.driver.params, "reference_model", None) is not None,
-                getattr(self.driver.params, "reference_inclination", None) is not None,
-                getattr(self.driver.params, "reference_declination", None) is not None,
+                getattr(self.driver.params.models, "reference_model", None) is not None,
+                getattr(self.driver.params.models, "reference_inclination", None)
+                is not None,
+                getattr(self.driver.params.models, "reference_declination", None)
+                is not None,
             )
 
             self._vector_inversion_directive = directives.VectorInversion(

simpeg_drivers/components/models.py

@@ -220,18 +220,16 @@ def reference_model(self) -> np.ndarray | None:
 
         ref_model = mref.copy()
 
-        if (
-            self.is_sigma
-            and self.driver.params.models.model_type == "Resistivity (Ohm-m)"
-        ):
-            ref_model = 1 / ref_model
+        if self.is_sigma:
+            if self.driver.params.models.model_type == "Resistivity (Ohm-m)":
+                ref_model = 1 / ref_model
 
-        ref_model = np.log(ref_model) if self.is_sigma else ref_model
+            ref_model = np.log(ref_model)
 
         if self.is_vector:
             field_vecs = dip_azimuth2cartesian(
-                self.starting_inclination,
-                self.starting_declination,
+                self.reference_inclination,
+                self.reference_declination,
             )
             ref_model = (field_vecs.T * ref_model).flatten()
 

tests/run_tests/driver_mvi_test.py

@@ -20,6 +20,7 @@
 from geoh5py.objects import Curve
 from geoh5py.workspace import Workspace
 
+from simpeg_drivers.components.factories import DirectivesFactory
 from simpeg_drivers.potential_fields import (
     MVIForwardOptions,
     MVIInversionOptions,
@@ -180,6 +181,58 @@ def test_magnetic_vector_run(
             check_target(output, target_mvi_run)
 
 
+def test_magnetic_vector_reference(
+    tmp_path: Path,
+    n_grid_points=3,
+    refinement=(2,),
+):
+    # Run the forward
+    opts = SyntheticsComponentsOptions(
+        method="magnetic_vector",
+        survey=SurveyOptions(
+            n_stations=n_grid_points, n_lines=n_grid_points, drape=5.0
+        ),
+        mesh=MeshOptions(refinement=refinement),
+        model=ModelOptions(anomaly=0.05),
+    )
+    with get_workspace(tmp_path / "inversion_test.ui.geoh5") as geoh5:
+        components = SyntheticsComponents(geoh5, options=opts)
+
+        tmi = components.survey.add_data(
+            {"tmi": {"values": np.random.randn(components.survey.n_vertices)}}
+        )
+        inducing_field = (50000.0, 90.0, 0.0)
+        params = MVIInversionOptions.build(
+            geoh5=geoh5,
+            mesh=components.mesh,
+            topography_object=components.topography,
+            inducing_field_strength=inducing_field[0],
+            inducing_field_inclination=inducing_field[1],
+            inducing_field_declination=inducing_field[2],
+            tmi_channel=tmi,
+            tmi_uncertainty=5.0,
+            data_object=components.survey,
+            starting_model=components.model,
+            reference_model=1.0,
+            reference_inclination=30,
+            reference_declination=0,
+        )
+    driver = MVIInversionDriver(params)
+
+    directives = DirectivesFactory(driver)
+    assert np.all(directives.vector_inversion_directive.reference_angles)
+    assert np.all(driver.models.reference_inclination == 30)
+    assert np.all(driver.models.reference_declination == 0)
+
+    np.allclose(
+        np.kron(
+            np.r_[0, np.cos(-np.deg2rad(30)), np.sin(-np.deg2rad(30))],
+            np.ones(driver.models.n_active),
+        ),
+        driver.models.reference_model,
+    )
+
+
 if __name__ == "__main__":
     # Full run
     test_magnetic_vector_fwr_run(Path("./"), n_grid_points=20, refinement=(4, 8))