diff --git a/LoopStructural/interpolators/supports/_2d_structured_grid.py b/LoopStructural/interpolators/supports/_2d_structured_grid.py index 72f641c7..ffc9f6be 100644 --- a/LoopStructural/interpolators/supports/_2d_structured_grid.py +++ b/LoopStructural/interpolators/supports/_2d_structured_grid.py @@ -3,7 +3,6 @@ """ -from ast import Not import logging from typing import Tuple diff --git a/LoopStructural/modelling/__init__.py b/LoopStructural/modelling/__init__.py index 2521ab0f..e720e8dd 100644 --- a/LoopStructural/modelling/__init__.py +++ b/LoopStructural/modelling/__init__.py @@ -6,12 +6,12 @@ __all__ = [ "GeologicalModel", "ProcessInputData", - "Loop3DView", "Map2LoopProcessor", "LoopProjectfileProcessor", ] from ..utils import getLogger from ..utils import LoopImportError +from .core.geological_model import GeologicalModel logger = getLogger(__name__) from ..modelling.input import ( diff --git a/LoopStructural/modelling/features/_geological_feature.py b/LoopStructural/modelling/features/_geological_feature.py index 557134d3..a96179e1 100644 --- a/LoopStructural/modelling/features/_geological_feature.py +++ b/LoopStructural/modelling/features/_geological_feature.py @@ -6,7 +6,7 @@ from ...modelling.features import BaseFeature from ...utils import getLogger from ...modelling.features import FeatureType -from ...interpolators import GeologicalInterpolator, DiscreteInterpolator +from ...interpolators import GeologicalInterpolator import numpy as np from typing import Optional, List, Union from ...datatypes import ValuePoints, VectorPoints diff --git a/LoopStructural/modelling/features/fold/_fold_rotation_angle.py b/LoopStructural/modelling/features/fold/_fold_rotation_angle.py index 518165fc..2b548f9b 100644 --- a/LoopStructural/modelling/features/fold/_fold_rotation_angle.py +++ b/LoopStructural/modelling/features/fold/_fold_rotation_angle.py @@ -2,7 +2,8 @@ from scipy.optimize import curve_fit from ....modelling.features.fold._fold_rotation_angle_feature import ( - fourier_series, + # fourier_series, + trigo_fold_profile, ) from ....modelling.features.fold import SVariogram @@ -53,9 +54,10 @@ def fit_fourier_series( wl = self.svario.find_wavelengths(lags=lags, nlag=nlag, lag=lag) # for now only consider single fold wavelength wl = wl[0] - guess = np.zeros(4) - guess[3] = wl # np.max(limb_wl) - logger.info(f"Guess: {guess[0]} {guess[1]} {guess[2]} {guess[3]}") + guess = np.zeros(3) + guess[1] = wl + guess[2] = np.deg2rad(45) # np.max(limb_wl) + logger.info(f"Guess: {guess[0]} {guess[1]} {guess[2]} ") # mask nans mask = np.logical_or(~np.isnan(self.fold_frame_coordinate), ~np.isnan(self.rotation_angle)) logger.info( @@ -71,7 +73,7 @@ def fit_fourier_series( try: # try fitting using wavelength guess popt, pcov = curve_fit( - fourier_series, + trigo_fold_profile, self.fold_frame_coordinate[mask], np.tan(np.deg2rad(self.rotation_angle[mask])), guess, @@ -81,7 +83,7 @@ def fit_fourier_series( # if fitting failed, try with just 0s logger.info("Running curve fit without initial guess") popt, pcov = curve_fit( - fourier_series, + trigo_fold_profile, self.fold_frame_coordinate[mask], np.tan(np.deg2rad(self.rotation_angle[mask])), ) @@ -90,15 +92,14 @@ def fit_fourier_series( popt = guess logger.error("Could not fit curve to S-Plot, check the wavelength") self.fitted_params = popt - logger.info(f"Fitted: {popt[0]} {popt[1]} {popt[2]} {popt[3]}") + logger.info(f"Fitted: {popt[0]} {popt[1]} {popt[2]} ") self.fold_rotation_function = lambda x: np.rad2deg( np.arctan( - fourier_series( + trigo_fold_profile( x, self.fitted_params[0], self.fitted_params[1], self.fitted_params[2], - self.fitted_params[3], ) ) ) diff --git a/LoopStructural/modelling/features/fold/_fold_rotation_angle_feature.py b/LoopStructural/modelling/features/fold/_fold_rotation_angle_feature.py index cab62603..5abe6ca8 100644 --- a/LoopStructural/modelling/features/fold/_fold_rotation_angle_feature.py +++ b/LoopStructural/modelling/features/fold/_fold_rotation_angle_feature.py @@ -1,6 +1,5 @@ import numpy as np from ....modelling.features import BaseFeature - from ....utils import getLogger logger = getLogger(__name__) @@ -46,6 +45,40 @@ def evaluate_value(self, location): return r +# class BaseFoldProfile(ABCMeta): +# def __init__(self): +# pass + +# @abstractmethod +# def params(self): +# pass + + +# @abstractmethod +# def __call__(self, s): +# pass + + +# class TrigoFoldProfile(BaseFoldProfile): +# def __init__(self, origin, wavelength, inflectionpointangle): +# self.origin = origin +# self.wavelength = wavelength +# self.inflectionpointangle = inflectionpointangle +# self.tan_alpha_delta_half = np.tan(self.inflectionpointangle) + +# self.tan_alpha_shift = 0 + +# def position_in_period(self, s): +# return (s - self.origin) / self.wavelength + +# def __call__(self, s): +# x = self.position_in_period(s) +# return np.rad2deg( +# np.arctan(self.tan_alpha_delta_half * np.sin(2 * np.pi * x) + self.tan_alpha_shift) +# ) +# class FourierSeriesFoldProfile(BaseFoldProfile): + + def fourier_series(x, c0, c1, c2, w): """ @@ -65,3 +98,23 @@ def fourier_series(x, c0, c1, c2, w): # v.fill(c0) v = c0 + c1 * np.cos(2 * np.pi / w * x) + c2 * np.sin(2 * np.pi / w * x) return np.rad2deg(np.arctan(v)) + + +def trigo_fold_profile(s, origin, wavelength, inflectionpointangle): + """ + + Parameters + ---------- + s + origin + wavelength + inflectionpointangle + + Returns + ------- + + """ + tan_alpha_delta_half = np.tan(inflectionpointangle) + tan_alpha_shift = 0 + x = (s - origin) / wavelength + return np.rad2deg(np.arctan(tan_alpha_delta_half * np.sin(2 * np.pi * x) + tan_alpha_shift))