diff --git a/README.md b/README.md index 11ba1473..e199af1b 100644 --- a/README.md +++ b/README.md @@ -6,16 +6,16 @@ [**Detailed documentation: https://rpoleski.github.io/MulensModel/**](https://rpoleski.github.io/MulensModel/) -[Latest release: 2.16.0](https://github.com/rpoleski/MulensModel/releases/latest) and we're working on further developing the code. +[Latest release: 2.17.0](https://github.com/rpoleski/MulensModel/releases/latest) and we're working on further developing the code. MulensModel can generate a microlensing light curve for a given set of microlensing parameters, fit that light curve to some data, and return a chi2 value. That chi2 can then be input into an arbitrary likelihood function to find the best-fit parameters. If you want to learn more about microlensing, please visit [Microlensing Source website](http://microlensing-source.org/). Currently, MulensModel supports: -* Lens Systems: point lens or binary lens. **New: shear and convergence allowed for both point and binary lenses.** -* Source Stars: single source or binary source. **Xallarap effect is currently developed:** [see this branch](https://github.com/rpoleski/MulensModel/tree/xal). -* Effects: finite source (1-parameter), parallax (satellite & annual), binary lens orbital motion, different parametrizations of microlensing models. +* Lens Systems: point lens or binary lens. Shear and convergence allowed for both point and binary lenses. +* Source Stars: single source or binary source. +* Effects: finite source (1-parameter), parallax (satellite & annual), binary lens orbital motion, **xallarap effect (new)**, different parametrizations of microlensing models. Need more? Open [an issue](https://github.com/rpoleski/MulensModel/issues), start [a discussion](https://github.com/rpoleski/MulensModel/discussions), or send us an e-mail. diff --git a/documents/images/orbit_parameters/NOTES b/documents/images/orbit_parameters/NOTES new file mode 100644 index 00000000..cb9c291d --- /dev/null +++ b/documents/images/orbit_parameters/NOTES @@ -0,0 +1,9 @@ +Image source: +https://en.wikipedia.org/wiki/File:Orbit1.svg + +This file is licensed under the Creative Commons Attribution-Share Alike 3.0 Unported license. +https://creativecommons.org/licenses/by-sa/3.0/deed.en + +Attribution: Lasunncty at the English Wikipedia +https://en.wikipedia.org/wiki/User:Lasunncty + diff --git a/documents/images/orbit_parameters/Orbit1.svg b/documents/images/orbit_parameters/Orbit1.svg new file mode 100644 index 00000000..9f53ccc7 --- /dev/null +++ b/documents/images/orbit_parameters/Orbit1.svg @@ -0,0 +1,44 @@ + + + + + + + + + + + + + + + + + + + + + + + + + Lengte van de klimmende knoopLongitud del node ascendentlongitud del nodo ascendente승교점 경도Longitude of ascending node + Argument van het periapsisArgument del periàpsideargumento de periapsis근일점 편각Argument of periapsis + Ware anomalieAnomalia veritableanomolía verdadera진근점 이각True anomaly + InclinatieInclinacióinclinación경사Inclination + Klimmende knoopNode ascendentnodo ascendente승교점Ascending node + ReferentierichtingDirecció dereferènciadirrección de referencia기준방향Referencedirection + HemellichaamCos celestecuerpo celeste천체Celestial body + ReferentievlakPla de referènciaplano de referencia기준면Plane of reference + BaanÒrbitaórbita궤도Orbit + + + Ω + ω + ν + i + ☊ + ♈︎ + + + \ No newline at end of file diff --git a/documents/images/orbit_parameters/Orbit1_svg.png b/documents/images/orbit_parameters/Orbit1_svg.png new file mode 100644 index 00000000..e39b19fe Binary files /dev/null and b/documents/images/orbit_parameters/Orbit1_svg.png differ diff --git a/documents/logo/logoMM_crop_4_744x520.png b/documents/logo/logoMM_crop_4_744x520.png new file mode 100644 index 00000000..a208d619 Binary files /dev/null and b/documents/logo/logoMM_crop_4_744x520.png differ diff --git a/documents/parameter_names.pdf b/documents/parameter_names.pdf index 929bc924..85c3e26f 100644 Binary files a/documents/parameter_names.pdf and b/documents/parameter_names.pdf differ diff --git a/documents/parameter_names.tex b/documents/parameter_names.tex index 46234ed7..824c74ff 100644 --- a/documents/parameter_names.tex +++ b/documents/parameter_names.tex @@ -2,26 +2,38 @@ \usepackage{geometry} \geometry{textwidth=500pt,top=20mm,bottom=20mm} +\usepackage{pdflscape} % for landscape mode +\usepackage{longtable} +\usepackage{caption} % for captionsetup +\usepackage{graphicx} +\usepackage{hyperref} \newcommand\MM{{\tt MulensModel}} \begin{document} % ######################################## +\begin{center} +%\centering +\includegraphics[width=0.5\textwidth]{logo/logoMM_crop_4_744x520.png} +\end{center} + +\vspace*{3cm} + \begin{center} {\LARGE Microlensing parameters in \MM}\\ \bigskip Radek Poleski\\ -last update: Jan 2022 +last update: Jul 2023 \end{center} \bigskip -Microlensing parameters in \MM\, class {\tt ModelParameters}: +Microlensing parameters in \MM\, class {\tt ModelParameters} are presented on the next page: -\begin{table*}[!h] -\begin{tabular}{l l l p{10cm}} -Parameter & Name in & Unit & Description \\ - & \MM & & \\ +\begin{landscape} +\captionsetup{width=20cm} +\begin{longtable}{l l l p{12cm}} +Parameter & Name in \MM & Unit & Description \\ \hline $t_0$ & {\tt t\_0} & & The time of the closest approach between the source and the lens. \\ $u_0$ & {\tt u\_0} & & The impact parameter between the source and the lens center of mass. \\ @@ -32,26 +44,35 @@ $\pi_{{\rm E}, N}$ & {\tt pi\_E\_N} & & The North component of the microlensing parallax vector. \\ $\pi_{{\rm E}, E}$ & {\tt pi\_E\_E} & & The East component of the microlensing parallax vector. \\ $t_{0,{\rm par}}$ & {\tt t\_0\_par} & & The reference time for parameters in parallax models.$^a$ \\ +$K$ & {\tt convergence\_K} & & External mass sheet convergence. \\ +$G$ & {\tt shear\_G} & & External mass sheet shear; complex valued to represent both the magnitude and angle relative to the binary axis. \\ $s$ & {\tt s} & & The projected separation between the lens primary and its companion as a fraction of the Einstein ring radius. \\ $q$ & {\tt q} & & The mass ratio between the lens companion and the lens primary $q \equiv m_2/m_1$. \\ $\alpha$ & {\tt alpha} & deg. & The angle between the source trajectory and the binary axis. \\ $ds/dt$ & {\tt ds\_dt} & yr$^{-1}$ & The rate of change of the separation. \\ -$K$ (convergence) & {\tt convergence\_K} & & External mass sheet convergence. \\ -$G$ (shear) & {\tt shear\_G} & & External mass sheet shear; complex valued to represent both the magnitude and angle relative to the binary axis. \\ $d\alpha/dt$ & {\tt dalpha\_dt} & deg.~yr$^{-1}$ & The rate of change of $\alpha$. \\ $t_{0,{\rm kep}}$ & {\tt t\_0\_kep} & & The reference time for lens orbital motion calculations.$^a$ \\ $x_{\rm caustic, in}$ & {\tt x\_caustic\_in} & & Curvelinear coordinate of caustic entrance for a binary lens model.$^b$ \\ $x_{\rm caustic, out}$ & {\tt x\_caustic\_out} & & Curvelinear coordinate of caustic exit for a binary lens model.$^b$ \\ $t_{\rm caustic, in}$ & {\tt t\_caustic\_in} & & Epoch of caustic exit for a binary lens model.$^b$ \\ $t_{\rm caustic, out}$ & {\tt t\_caustic\_out} & & Epoch of caustic exit for a binary lens model.$^b$ \\ +$\xi_P$ & \texttt{xi\_period} & d & The orbital period of xallarap.\\ +$\xi_a$ & \texttt{xi\_semimajor\_axis} & & The semimajor axis of a xallarap orbit as a fraction of the Einstein ring.\\ +$\xi_i$ & \texttt{xi\_inclination} & deg & The inclination of a xallarap orbit.$^c$\\ +$\xi_\Omega$ & \texttt{xi\_Omega\_node} & deg & The longitude of the ascending node of a xallarap orbit.$^c$\\ +$\xi_u$ & \texttt{xi\_argument\_of\_latitude\_reference} & deg & The argument of latitude at the reference epoch ($t_{0,\chi}$). The argument of latitude is a sum of true anomaly ($\nu$, changes with time) and the argument of periapsis ($\omega$, orbit parameter, i.e., does not change with time): $u = \nu + \omega$.$^c$\\ +$\xi_e$ & \texttt{xi\_eccentricity} & & The eccentricity of a xallarap orbit.\\ +$\xi_\omega$ & \texttt{xi\_omega\_periapsis} & deg & The argument of periapsis of a xallarap orbit.$^c$\\ +$t_{0,\xi}$ & \texttt{t\_0\_xi} & & The reference epoch for parameters in xallarap models.$^a$\\ \hline -\end{tabular} \caption{Notes: \newline -$^a$ -- $t_{0,{\rm par}}$ and $t_{0,{\rm kep}}$ are reference parameters, hence, do not change these during fitting. \newline -$^b$ -- The four parameters of binary lens in Cassan (2008) parameterization ($x_{\rm caustic, in}$, $x_{\rm caustic, out}$, $t_{\rm caustic, in}$, and $t_{\rm caustic, out}$) are used instead of ($t_0$, $u_0$, $t_{\rm E}$, and $\alpha$). +$^a$ -- $t_{0,{\rm par}}$, $t_{0,{\rm kep}}$, and $t_{0,\chi}$ are reference parameters, hence, do not change these during fitting. \newline +$^b$ -- The four parameters of binary lens in Cassan (2008) parameterization ($x_{\rm caustic, in}$, $x_{\rm caustic, out}$, $t_{\rm caustic, in}$, $t_{\rm caustic, out}$) are used instead of ($t_0$, $u_0$, $t_{\rm E}$, $\alpha$). \newline +$^c$ -- The orbital angles are illustrated in Figure~1. %\label{} } -\end{table*} +\end{longtable} +\end{landscape} Some of the parameters can be defined separately for each of the sources in binary source models. In that case, add {\tt \_1} or {\tt \_2} to parameter name. These are: @@ -72,5 +93,18 @@ \item coordinates of space telescopes -- \texttt{Model.get\_satellite\_coords}. \end{itemize} +\begin{figure} +\centering +\includegraphics[width=0.56\textwidth]{images/orbit_parameters/Orbit1_svg} +\rotatebox[origin=l]{90}{\tiny image by Lasunncty/WikiCommons, license: CC BY-SA 3.0} +\caption{ +Definition of orbital angles. +There are Wikipedia articles that give more details: +\href{https://en.wikipedia.org/wiki/Orbital_elements}{orbital elements} and +\href{https://en.wikipedia.org/wiki/Argument_of_latitude}{argument of latitude}. +For xallarap, the reference direction is the relative lens-source proper motion direction and the reference plane is the plane of the sky. +} +\end{figure} + \end{document} diff --git a/examples/example_16/ulens_model_fit.py b/examples/example_16/ulens_model_fit.py index 19372165..5ea66317 100644 --- a/examples/example_16/ulens_model_fit.py +++ b/examples/example_16/ulens_model_fit.py @@ -38,7 +38,7 @@ except Exception: raise ImportError('\nYou have to install MulensModel first!\n') -__version__ = '0.32.0' +__version__ = '0.33.0' class UlensModelFit(object): @@ -67,7 +67,23 @@ class UlensModelFit(object): starting_parameters: *dict* Starting values of the parameters. It also indicates the EMCEE fitting mode. + There are two possibilities for the information provided. + First, you can provide a name of file with sets of + parameters and names of parameters. For example: + + .. code-block:: python + + { + 'file': 'STARTING_POINTS.txt', + 'parameters': 't_0 u_0 t_E' + } + + In that case the text file has three columns and + at least 2 * 3 = 6 rows with values. + + Second, you can provide distribution that defines distributions + from which values of parameters will be drawn. Keys of this *dict* are microlensing parameters recognized by *MulensModel* and values are *str*. First word indicates the distribution (allowed: ``gauss``, ``uniform``, and @@ -295,7 +311,7 @@ def __init__( ): self._check_MM_version() self._photometry_files = photometry_files - self._starting_parameters = starting_parameters + self._starting_parameters_input = starting_parameters self._prior_limits = prior_limits self._model_parameters = model self._fixed_parameters = fixed_parameters @@ -310,6 +326,7 @@ def __init__( self._set_default_parameters() if self._task == 'fit': self._guess_fitting_method() + self._check_starting_parameters_type() self._set_fit_parameters_unsorted() self._check_imports() @@ -327,7 +344,7 @@ def _which_task(self): Check if input parameters indicate run_fit() or plot_best_model() will be run. """ - if self._starting_parameters is not None: + if self._starting_parameters_input is not None: fit = True elif self._prior_limits is not None: fit = True @@ -373,7 +390,7 @@ def _check_unnecessary_settings_plot(self): Make sure that there arent' too many parameters specified for: self._task == 'plot' """ - keys = ['_starting_parameters', '_min_values', '_max_values', + keys = ['_starting_parameters_input', '_min_values', '_max_values', '_fitting_parameters', '_prior_limits'] for key in keys: if getattr(self, key) is not None: @@ -435,7 +452,7 @@ def _guess_fitting_method(self): guess what is the fitting method based on parameters provided """ method = None - if self._starting_parameters is not None: + if self._starting_parameters_input is not None: method = "EMCEE" if self._prior_limits is not None: if method is not None: @@ -452,19 +469,54 @@ def _guess_fitting_method(self): "starting_parameters or prior_limits, respectively.") self._fit_method = method + def _check_starting_parameters_type(self): + """ + Check if starting parameters are read from file or + will be drawn from distributions specified. + """ + if 'file' in self._starting_parameters_input: + in_type = 'file' + keys_expected = {'file', 'parameters'} + keys = set(self._starting_parameters_input.keys()) + if keys != keys_expected: + error = ('Wrong starting parameters keys. Expected: ' + + str(keys_expected) + '; Provided: ' + str(keys)) + raise KeyError(error) + else: + in_type = 'draw' + + self._starting_parameters_type = in_type + def _set_fit_parameters_unsorted(self): """ Find what are the fitted parameters. It will be sorted later. """ if self._fit_method == "EMCEE": - unsorted_keys = self._starting_parameters.keys() + if self._starting_parameters_type == 'draw': + unsorted_keys = self._starting_parameters_input.keys() + elif self._starting_parameters_type == 'file': + unsorted_keys = self._get_unsorted_starting_parameters() + else: + raise ValueError( + 'unexpected: ' + str(self._starting_parameters_type)) elif self._fit_method == "MultiNest": unsorted_keys = self._prior_limits.keys() else: raise ValueError('unexpected method error') + self._fit_parameters_unsorted = list(unsorted_keys) self._n_fit_parameters = len(self._fit_parameters_unsorted) + def _get_unsorted_starting_parameters(self): + """ + Make sure that a list of parameters is provided + """ + parameters = self._starting_parameters_input['parameters'] + if isinstance(parameters, (str)): + parameters = parameters.split() + + return parameters + def _check_imports(self): """ check if all the required packages are imported @@ -476,6 +528,7 @@ def _check_imports(self): required_packages.add('emcee') elif self._fit_method == "MultiNest": required_packages.add('pymultinest') + if self._plots is not None and 'triangle' in self._plots: required_packages.add('corner') @@ -490,6 +543,7 @@ def _check_imports(self): "\nFor corner package it's enough that you run:\nwget " + "https://raw.githubusercontent.com/dfm/corner.py/" + "v2.0.0/corner/corner.py") + raise ImportError(message) def run_fit(self): @@ -515,10 +569,12 @@ def run_fit(self): self._parse_fit_constraints() if self._fit_method == "EMCEE": self._parse_starting_parameters() + self._check_fixed_parameters() self._make_model_and_event() if self._fit_method == "EMCEE": - self._generate_random_parameters() + self._get_starting_parameters() + self._setup_fit() self._run_fit() self._finish_fit() @@ -1139,7 +1195,13 @@ def _get_n_walkers(self): if 'n_walkers' in self._fitting_parameters: self._n_walkers = self._fitting_parameters['n_walkers'] else: - self._n_walkers = 4 * len(self._starting_parameters) + if self._starting_parameters_type == 'file': + self._n_walkers = None + elif self._starting_parameters_type == 'draw': + self._n_walkers = 4 * self._n_fit_parameters + else: + raise ValueError( + 'Unexpected: ' + self._starting_parameters_type) def _parse_fitting_parameters_MultiNest(self): """ @@ -1443,6 +1505,36 @@ def _read_prior_t_E_data(self): self._prior_t_E_data['function'] = function def _parse_starting_parameters(self): + """ + Read the starting parameters from file or + change the format of provided information. + """ + if self._starting_parameters_type == 'file': + self._read_starting_parameters_from_file() + elif self._starting_parameters_type == 'draw': + self._parse_starting_parameters_to_be_drawn() + else: + raise ValueError( + 'unexpected: ' + str(self._starting_parameters_type)) + + def _read_starting_parameters_from_file(self): + """ + Read starting parameters from a file. + """ + file_name = self._starting_parameters_input['file'] + try: + data = np.loadtxt(file_name, unpack=True) + except Exception: + raise RuntimeError('Error while reading file: ' + file_name) + + if len(data.shape) != 2 or data.shape[0] != self._n_fit_parameters: + raise ValueError( + 'Something wrong in shape of data read from ' + file_name) + + self._starting_parameters_list = data.T.tolist() + self._n_walkers = len(self._starting_parameters_list) + + def _parse_starting_parameters_to_be_drawn(self): """ replace self._starting_parameters with dict that has values [*str*, *float*, ...] @@ -1451,7 +1543,7 @@ def _parse_starting_parameters(self): accepted_types = ['gauss', 'uniform', 'log-uniform'] out = dict() - for (key, value) in self._starting_parameters.items(): + for (key, value) in self._starting_parameters_input.items(): words = value.split() if words[0] not in accepted_types: raise ValueError( @@ -1623,6 +1715,28 @@ def _get_example_parameters_EMCEE(self): """ get sample values of parameters for EMCEE - only to make mm.Model """ + if self._starting_parameters_type == 'file': + return self._extract_example_parameters_EMCEE() + elif self._starting_parameters_type == 'draw': + return self._draw_example_parameters_EMCEE() + else: + raise ValueError( + 'unexpected: ' + str(self._starting_parameters_type)) + + def _extract_example_parameters_EMCEE(self): + """ + Provide example parameters that will be used to initialize + MM.Model and MM.Event. Sets of parameters are already read, + so we just take 0-th set. + """ + keys = self._fit_parameters_unsorted + values = self._starting_parameters_list[0] + return dict(zip(keys, values)) + + def _draw_example_parameters_EMCEE(self): + """ + Randomly draw parameters or EMCEE - only to make mm.Model. + """ parameters = dict() for (key, value) in self._starting_parameters.items(): # We treat Cassan08 case differently so that @@ -1644,6 +1758,21 @@ def _get_example_parameters_EMCEE(self): raise ValueError('internal error: ' + value[0]) return parameters + def _get_starting_parameters(self): + """ + Generate random parameters or read them from file. + Check if there are enough of them and store. + """ + if self._starting_parameters_type == 'file': + starting = self._starting_parameters_list + elif self._starting_parameters_type == 'draw': + starting = self._generate_random_parameters() + else: + raise ValueError( + 'unexpected: ' + str(self._starting_parameters_type)) + + self._check_and_store_generated_random_parameters(starting) + def _generate_random_parameters(self): """ Generate a number of starting parameters values. @@ -1660,9 +1789,7 @@ def _generate_random_parameters(self): max_iteration, settings) starting.append(values) - starting = np.array(starting).T.tolist() - - self._check_generated_random_parameters(starting) + return np.array(starting).T.tolist() def _get_samples_from_distribution(self, n, settings): """ @@ -1681,8 +1808,28 @@ def _get_samples_from_distribution(self, n, settings): values = np.exp(np.random.uniform(beg, end, n)) else: raise ValueError('Unrecognized keyword: ' + settings[0]) + return values + def _check_and_store_generated_random_parameters(self, starting): + """ + Check if the set of points provided has at least self._n_walkers + points inside the prior and save these. + """ + verified = self._check_generated_random_parameters(starting) + + if len(verified) < self._n_walkers: + raise ValueError( + "Couldn't generate required starting points in a prior. " + "Most probably you have to correct at least one of: " + "starting_parameters, min_values, max_values, or " + "fit_constraints.\nGot " + str(len(verified)) + " walkers in " + "the prior, but required " + str(self._n_walkers) + ".\n" + "If you think the code should work with your settings, " + "then please contact Radek Poleski.") + + self._kwargs_EMCEE['initial_state'] = verified + def _check_generated_random_parameters(self, starting): """ Check if the set of points provided has at least self._n_walkers @@ -1698,16 +1845,7 @@ def _check_generated_random_parameters(self, starting): if len(out) == self._n_walkers: break - if len(out) < self._n_walkers: - raise ValueError( - "Couldn't generate required starting points in a prior. " - "Most probably you have to correct at least one of: " - "starting_parameters, min_values, max_values, or " - "fit_constraints.\nGot " + str(len(out)) + " walkers in " - "the prior, but required " + str(self._n_walkers) + ".\n" - "If you think the code should work with your settings, " - "then please contact Radek Poleski.") - self._kwargs_EMCEE['initial_state'] = out + return out def _ln_prob(self, theta): """ diff --git a/source/MulensModel/__init__.py b/source/MulensModel/__init__.py index 75514328..cf306af9 100644 --- a/source/MulensModel/__init__.py +++ b/source/MulensModel/__init__.py @@ -15,6 +15,7 @@ from MulensModel.modelparameters import ModelParameters, which_parameters from MulensModel.mulensdata import MulensData from MulensModel.mulensobjects import * +from MulensModel.orbits import * from MulensModel.pointlens import PointLens, get_pspl_magnification from MulensModel.pointlenswithshear import PointLensWithShear from MulensModel.satelliteskycoord import SatelliteSkyCoord diff --git a/source/MulensModel/fitdata.py b/source/MulensModel/fitdata.py index 76d8c9da..3f19e73a 100644 --- a/source/MulensModel/fitdata.py +++ b/source/MulensModel/fitdata.py @@ -551,6 +551,10 @@ def _check_for_gradient_implementation(self, parameters): raise NotImplementedError('Event.chi2_gradient() is not working ' 'for finite source models yet') + if self.model.parameters.is_xallarap: + raise NotImplementedError('Gradient for xallarap models is not ' + 'implemented yet') + def get_chi2_gradient(self, parameters): """ Fits fluxes and calculates chi^2 gradient (also called Jacobian), i.e., diff --git a/source/MulensModel/model.py b/source/MulensModel/model.py index e9b1d355..56f06345 100644 --- a/source/MulensModel/model.py +++ b/source/MulensModel/model.py @@ -642,27 +642,35 @@ def _plot_single_trajectory(self, times, parameters, satellite_skycoord, self._plt_plot(trajectory.x, trajectory.y, kwargs) if arrow: - if len(times) > 2: - index = int(len(times)/2) - else: - index = 0 - x_0 = trajectory.x[index] - y_0 = trajectory.y[index] - d_x = trajectory.x[index+1] - x_0 - d_y = trajectory.y[index+1] - y_0 - dd = 1e6 * (d_x*d_x + d_y*d_y)**.5 - - xlim = plt.xlim() - ylim = plt.ylim() - width = np.abs(xlim[1]-xlim[0]) * np.abs(ylim[1]-ylim[0]) - width = width**.5 / 100. - - color = kwargs.get('color', 'black') - kwargs_ = {'width': width, 'color': color, 'lw': 0, - 'zorder': -np.inf} - if arrow_kwargs is not None: - kwargs_.update(arrow_kwargs) - plt.arrow(x_0, y_0, d_x/dd, d_y/dd, **kwargs_) + self._plot_arrow(times, trajectory, kwargs, arrow_kwargs) + + def _plot_arrow(self, times, trajectory, kwargs, arrow_kwargs): + """ + Plot arrow for given trajectory. + """ + width_scaling_factor = 0.01 + if len(times) > 2: + index = int(len(times)/2) + else: + index = 0 + + x_0 = trajectory.x[index] + y_0 = trajectory.y[index] + d_x = trajectory.x[index+1] - x_0 + d_y = trajectory.y[index+1] - y_0 + dd = 1e6 * (d_x*d_x + d_y*d_y)**.5 + + xlim = plt.xlim() + ylim = plt.ylim() + width = np.abs(xlim[1]-xlim[0]) * np.abs(ylim[1]-ylim[0]) + width = width**.5 * width_scaling_factor + + color = kwargs.get('color', 'black') + kwargs_ = {'width': width, 'color': color, 'lw': 0, 'zorder': -np.inf} + if arrow_kwargs is not None: + kwargs_.update(arrow_kwargs) + + plt.arrow(x_0, y_0, d_x/dd, d_y/dd, **kwargs_) def plot_source(self, times=None, **kwargs): """ diff --git a/source/MulensModel/modelparameters.py b/source/MulensModel/modelparameters.py index 24f2857a..4fd9b67b 100644 --- a/source/MulensModel/modelparameters.py +++ b/source/MulensModel/modelparameters.py @@ -283,7 +283,7 @@ def _set_type(self, keys): sets self._type property, which indicates what type of a model we have """ types = ['finite source', 'parallax', 'Cassan08', - 'lens 2-parameter orbital motion', 'mass sheet'] + 'lens 2-parameter orbital motion', 'mass sheet', 'xallarap'] out = {type_: False for type_ in types} temp = { @@ -292,7 +292,9 @@ def _set_type(self, keys): 'Cassan08': 'x_caustic_in x_caustic_out t_caustic_in t_caustic_out', 'lens 2-parameter orbital motion': 'dalpha_dt ds_dt', - 'mass sheet': 'convergence_K shear_G'} + 'mass sheet': 'convergence_K shear_G', + 'xallarap': ('xi_period xi_semimajor_axis xi_inclination ' + 'xi_Omega_node xi_argument_of_latitude_reference')} parameter_to_type = dict() for (key, values) in temp.items(): @@ -317,17 +319,7 @@ def _check_types(self, alpha_defined): raise KeyError('Orbital motion of the lens requires two lens ' 'components but only one was provided.') - # Cassan08 prevents some other features: - if self._type['Cassan08']: - if self._type['parallax']: - raise KeyError('Currently we do not allow parallax for ' - 'Cassan (2008) parameterization.') - if self._type['lens 2-parameter orbital motion']: - raise KeyError('Cassan (2008) parameterization and parallax ' - 'are not allowed') - if n_sources > 1: - raise KeyError("Cassan (2008) parameterization doesn't work" - "for multi sources models") + self._check_types_for_Cassan08() if alpha_defined: if self._n_lenses == 1 and not self._type['mass sheet']: @@ -335,6 +327,60 @@ def _check_types(self, alpha_defined): 'You defined alpha for single lens model ' 'without external mass sheet. This is not allowed.') + if n_sources > 1 and self._type['xallarap']: + raise NotImplementedError('We have not yet implemented xallarap ' + 'and multiple luminous sources') + + def _check_valid_combination_1_source(self, keys): + """ + Check that the user hasn't over-defined the ModelParameters. + """ + # Make sure that there are no unwanted keys + allowed_keys = set(( + 't_0 u_0 t_E t_eff rho t_star pi_E pi_E_N pi_E_E t_0_par ' + 's q alpha dalpha_dt ds_dt t_0_kep convergence_K shear_G ' + 't_0_1 t_0_2 u_0_1 u_0_2 rho_1 rho_2 t_star_1 t_star_2 ' + 'x_caustic_in x_caustic_out t_caustic_in t_caustic_out ' + 'xi_period xi_semimajor_axis xi_inclination xi_Omega_node ' + 'xi_argument_of_latitude_reference xi_eccentricity ' + 'xi_omega_periapsis t_0_xi').split()) + difference = set(keys) - allowed_keys + if len(difference) > 0: + derived_1 = ['gamma', 'gamma_perp', 'gamma_parallel'] + if set(keys).intersection(derived_1): + msg = ('You cannot set gamma, gamma_perp, ' + + 'or gamma_parallel. These are derived parameters. ' + + 'You can set ds_dt and dalpha_dt instead.\n') + else: + msg = "" + msg += 'Unrecognized parameters: {:}'.format(difference) + raise KeyError(msg) + + if self._type['Cassan08']: + self._check_valid_combination_1_source_Cassan08(keys) + else: + self._check_valid_combination_1_source_standard(keys) + + def _check_types_for_Cassan08(self): + """ + Check if Cassan08 is used and if so, then make sure that + the trajectory is rectilinear and there is only one source. + """ + if not self._type['Cassan08']: + return + + types = ['parallax', 'xallarap', 'lens 2-parameter orbital motion'] + for type_ in types: + if self._type[type_]: + raise NotImplementedError( + 'Currently we do not allow Cassan (2008) ' + 'parameterization of binary lens and ' + type_) + + if self._n_sources > 1: + raise NotImplementedError( + "Cassan (2008) parameterization doesn't work for " + "multi sources models") + def _divide_parameters(self, parameters): """ Divide an input dict into 2 - each source separately. @@ -396,8 +442,24 @@ def __repr__(self): 't_0_kep': {'width': 13, 'precision': 5, 'unit': 'HJD'}, 'x_caustic_in': {'width': 13, 'precision': 7}, 'x_caustic_out': {'width': 13, 'precision': 7}, - 't_caustic_in': {'width': 19, 'precision': 5, 'unit': 'HJD'}, - 't_caustic_out': {'width': 19, 'precision': 5, 'unit': 'HJD'}, + 't_caustic_in': {'width': 13, 'precision': 5, 'unit': 'HJD'}, + 't_caustic_out': {'width': 13, 'precision': 5, 'unit': 'HJD'}, + 'xi_period': {'width': 10, 'precision': 4, + 'unit': 'd', 'name': 'xallarap period'}, + 'xi_semimajor_axis': {'width': 9, 'precision': 6, + 'name': 'xallarap semimajor axis'}, + 'xi_inclination': {'width': 11, 'precision': 5, 'unit': 'deg', + 'name': 'xallarap inclination'}, + 'xi_Omega_node': {'width': 11, 'precision': 5, 'unit': 'deg', + 'name': 'xallarap Omega node'}, + 'xi_argument_of_latitude_reference': { + 'width': 11, 'precision': 5, 'unit': 'deg', + 'name': 'xallarap argument of latitude reference'}, + 'xi_eccentricity': {'width': 8, 'precision': 6, + 'name': 'xallarap eccentricity'}, + 'xi_omega_periapsis': {'width': 11, 'precision': 5, 'unit': 'deg', + 'name': 'xallarap omega periapsis'}, + 't_0_xi': {'width': 13, 'precision': 5, 'unit': 'HJD'}, } # Add binary source parameters with the same settings. binary_source_keys = ['t_0_1', 't_0_2', 'u_0_1', 'u_0_2', @@ -411,27 +473,54 @@ def __repr__(self): if 'name' in form: raise KeyError('internal issue: {:}'.format(key)) + ordered_keys = [ + 't_0', 't_0_1', 't_0_2', 'u_0', 'u_0_1', 'u_0_2', 't_eff', 't_E', + 'rho', 'rho_1', 'rho_2', 't_star', 't_star_1', 't_star_2', + 'pi_E_N', 'pi_E_E', 't_0_par', 's', 'q', 'alpha', + 'convergence_K', 'shear_G', 'ds_dt', 'dalpha_dt', 't_0_kep', + 'x_caustic_in', 'x_caustic_out', 't_caustic_in', 't_caustic_out', + 'xi_period', 'xi_semimajor_axis', 'xi_inclination', + 'xi_Omega_node', 'xi_argument_of_latitude_reference', + 'xi_eccentricity', 'xi_omega_periapsis', 't_0_xi' + ] + variables = '' values = '' - - for key in formats.keys(): + for key in ordered_keys: if key not in keys: continue - form = formats[key] - fmt_1 = '{:>' + str(form['width']) - fmt_2 = fmt_1 + '.' + str(form['precision']) + 'f} ' - fmt_1 += '} ' - full_name = form.get('name', key) - if 'unit' in form: - full_name += " ({:})".format(form['unit']) + (full_name, value) = self._get_values_for_repr(formats[key], key) + (fmt_1, fmt_2) = self._get_formats_for_repr(formats[key], + full_name) variables += fmt_1.format(full_name) - value = getattr(self, key) - if isinstance(value, u.Quantity): - value = value.value values += fmt_2.format(value) return '{0}\n{1}'.format(variables, values) + def _get_values_for_repr(self, form, key): + """ + Get full name of the parameter and its value (float) + to be used by __rerp__(). + """ + full_name = form.get('name', key) + if 'unit' in form: + full_name += " ({:})".format(form['unit']) + + value = getattr(self, key) + if isinstance(value, u.Quantity): + value = value.value + + return (full_name, value) + + def _get_formats_for_repr(self, form, full_name): + """ + Extract formats to be used by __repr__(). + """ + fmt_1 = '{:>' + str(max([form['width'], len(full_name)])) + fmt_2 = fmt_1 + '.' + str(form['precision']) + 'f} ' + fmt_1 += '} ' + return (fmt_1, fmt_2) + def _check_valid_combination_2_sources(self, keys): """ make sure that there is no conflict between t_0 and t_0_1 etc. @@ -452,6 +541,7 @@ def _check_valid_combination_1_source_standard(self, keys): self._check_valid_combination_1_source_parallax(keys) self._check_valid_combination_1_source_mass_sheet(keys) self._check_valid_combination_1_source_binary_lens(keys) + self._check_valid_combination_1_source_xallarap(keys) def _check_valid_combination_1_source_t_0_u_0(self, keys): """ @@ -459,6 +549,7 @@ def _check_valid_combination_1_source_t_0_u_0(self, keys): """ if 't_0' not in keys: raise KeyError('t_0 must be defined') + if ('u_0' not in keys) and ('t_eff' not in keys): raise KeyError('not enough information to calculate u_0') @@ -474,8 +565,10 @@ def _check_valid_combination_1_source_t_E(self, keys): if (('rho' in keys) and ('t_star' in keys) and ('u_0' in keys) and ('t_eff' in keys)): raise KeyError('You cannot define rho, t_star, u_0, and t_eff') + if ('t_E' in keys) and ('rho' in keys) and ('t_star' in keys): raise KeyError('Only 1 or 2 of (t_E, rho, t_star) may be defined.') + if ('t_E' in keys) and ('u_0' in keys) and ('t_eff' in keys): raise KeyError('Only 1 or 2 of (u_0, t_E, t_eff) may be defined.') @@ -513,6 +606,7 @@ def _check_valid_combination_1_source_mass_sheet(self, keys): if ('shear_G' in keys) and ('alpha' not in keys): raise KeyError( 'A model with external mass sheet shear requires alpha.') + if ('shear_G' not in keys) and ('convergence_K' in keys): if 'alpha' in keys: raise KeyError( @@ -548,6 +642,27 @@ def _check_valid_combination_1_source_binary_lens(self, keys): raise KeyError( 't_0_kep makes sense only when orbital motion is defined.') + def _check_valid_combination_1_source_xallarap(self, keys): + """ + If xallarap parameters are defined, + then make sure there are all required parameters + """ + if not self._type['xallarap']: + return + + required = ('xi_period xi_semimajor_axis xi_inclination ' + 'xi_Omega_node xi_argument_of_latitude_reference').split() + for parameter in required: + if parameter not in keys: + raise KeyError(parameter) + + allowed = set(['xi_eccentricity', 'xi_omega_periapsis']) + n_used = len(set(keys).intersection(allowed)) + if n_used not in [0, len(allowed)]: + raise KeyError( + 'Error in defining xi_eccentricity and xi_omega_periapsis. ' + 'Both of them or neither should be defined.') + def _check_valid_combination_1_source_Cassan08(self, keys): """ Check parameters defined for Cassan 2008 parameterization. @@ -568,33 +683,6 @@ def _check_valid_combination_1_source_Cassan08(self, keys): raise KeyError('Both rho and t_star cannot be defined for ' + 'Cassan 08 parameterization.') - def _check_valid_combination_1_source(self, keys): - """ - Check that the user hasn't over-defined the ModelParameters. - """ - # Make sure that there are no unwanted keys - allowed_keys = set(( - 't_0 u_0 t_E t_eff rho t_star pi_E pi_E_N pi_E_E t_0_par ' - 's q alpha dalpha_dt ds_dt t_0_kep convergence_K shear_G ' - 't_0_1 t_0_2 u_0_1 u_0_2 rho_1 rho_2 t_star_1 t_star_2 ' - 'x_caustic_in x_caustic_out t_caustic_in t_caustic_out').split()) - difference = set(keys) - allowed_keys - if len(difference) > 0: - derived_1 = ['gamma', 'gamma_perp', 'gamma_parallel'] - if set(keys).intersection(derived_1): - msg = ('You cannot set gamma, gamma_perp, ' + - 'or gamma_parallel. These are derived parameters. ' + - 'You can set ds_dt and dalpha_dt instead.\n') - else: - msg = "" - msg += 'Unrecognized parameters: {:}'.format(difference) - raise KeyError(msg) - - if self._type['Cassan08']: - self._check_valid_combination_1_source_Cassan08(keys) - else: - self._check_valid_combination_1_source_standard(keys) - def _check_valid_parameter_values(self, parameters): """ Prevent user from setting negative (unphysical) values for @@ -620,12 +708,24 @@ def _check_valid_parameter_values(self, parameters): msg = "{:} must be a scalar: {:}, {:}" raise TypeError(msg.format(key, value, type(value))) - for name in ['x_caustic_in', 'x_caustic_out', 'q']: + for name in ['x_caustic_in', 'x_caustic_out']: if name in parameters.keys(): if parameters[name] < 0. or parameters[name] > 1.: + msg = "Parameter {:} has to be in [0, 1] range, not {:}" + raise ValueError(msg.format(name, parameters[name])) + + for name in ['q']: + if name in parameters.keys(): + if parameters[name] <= 0. or parameters[name] >= 1.: msg = "Parameter {:} has to be in (0, 1) range, not {:}" raise ValueError(msg.format(name, parameters[name])) + for name in ['xi_eccentricity']: + if name in parameters.keys(): + if parameters[name] < 0. or parameters[name] >= 1.: + msg = "Parameter {:} has to be in [0, 1) range, not {:}" + raise ValueError(msg.format(name, parameters[name])) + if 'shear_G' in parameters.keys(): if not isinstance(parameters['shear_G'], complex): raise TypeError("External shear (shear_G) must be complex") @@ -641,6 +741,14 @@ def _set_parameters(self, parameters): if parameter in self.parameters: self._set_time_quantity(parameter, self.parameters[parameter]) + angle_parameters = [ + 'alpha', 'xi_Omega_node', 'xi_inclination', + 'xi_argument_of_latitude_reference', 'xi_omega_periapsis'] + for parameter in angle_parameters: + if parameter in self.parameters: + self._warn_if_angle_outside_reasonable_range( + self.parameters[parameter], parameter) + def _update_sources(self, parameter, value): """ For multi-source models, update the values for all sources. @@ -651,6 +759,7 @@ def _update_sources(self, parameter, value): if parameter in self._source_1_parameters.parameters: setattr(self._source_1_parameters, parameter, value) + if parameter in self._source_2_parameters.parameters: setattr(self._source_2_parameters, parameter, value) @@ -924,8 +1033,22 @@ def alpha(self, new_alpha): self.parameters['alpha'] = new_alpha else: self.parameters['alpha'] = new_alpha * u.deg + self._warn_if_angle_outside_reasonable_range( + self.parameters['alpha'].to(u.deg).value, 'alpha') self._update_sources('alpha', new_alpha) + def _warn_if_angle_outside_reasonable_range(self, value, name): + """ + Check if value of given angle is in reasonable range and warn if not + """ + min_ = -360. + max_ = 540. + if isinstance(value, u.Quantity): + value = value.to(u.deg).value + if value < min_ or value > max_: + fmt = "Strange value of angle {:}: {:}" + warnings.warn(fmt.format(name, value), RuntimeWarning) + @property def q(self): """ @@ -1273,6 +1396,141 @@ def t_0_kep(self, new): self.parameters['t_0_kep'] = new self._update_sources('t_0_kep', new) + @property + def xi_period(self): + """ + *float* + + Orbital period of the source system (xallarap) in days. + """ + return self.parameters['xi_period'] + + @xi_period.setter + def xi_period(self, new_value): + if new_value < 0.: + raise ValueError('Xallarap period cannot be negative') + self.parameters['xi_period'] = new_value + + @property + def xi_semimajor_axis(self): + """ + *float* + + Semi-major axis of the source orbit (xallarap) in the theta_E units. + """ + return self.parameters['xi_semimajor_axis'] + + @xi_semimajor_axis.setter + def xi_semimajor_axis(self, new_value): + if new_value < 0.: + raise ValueError('Xallarap semimajor axis cannot be negative') + self.parameters['xi_semimajor_axis'] = new_value + + @property + def xi_Omega_node(self): + """ + *float* + + The longitude of the ascending node of the xallarap orbit, i.e., + the angle from relative lens-source proper motion direction + to the ascending node direction. + The units are degrees. + """ + return self.parameters['xi_Omega_node'] + + @xi_Omega_node.setter + def xi_Omega_node(self, new_value): + self._warn_if_angle_outside_reasonable_range(new_value, + 'xi_Omega_node') + self.parameters['xi_Omega_node'] = new_value + + @property + def xi_inclination(self): + """ + *float* + + The inclination of the xallarap orbit, i.e., + the angle between source-orbit plane and the sky plane. + The units are degrees. + """ + return self.parameters['xi_inclination'] + + @xi_inclination.setter + def xi_inclination(self, new_value): + self._warn_if_angle_outside_reasonable_range(new_value, + 'xi_inclination') + self.parameters['xi_inclination'] = new_value + + @property + def xi_argument_of_latitude_reference(self): + """ + *float* + + The argument of latitude for the xallarap orbit at :py:attr:`~t_0_xi`. + The argument of latitude is a sum of the true anomaly and + the argument of periapsis. In standard notation: u = nu + omega. + This parameter is internally used to calculate perihelion passage + (T_0 in standard notation). + The units are degrees. + """ + return self.parameters['xi_argument_of_latitude_reference'] + + @xi_argument_of_latitude_reference.setter + def xi_argument_of_latitude_reference(self, new_value): + self._warn_if_angle_outside_reasonable_range( + new_value, 'xi_argument_of_latitude_reference') + self.parameters['xi_argument_of_latitude_reference'] = new_value + + @property + def xi_eccentricity(self): + """ + *float* + + The eccentricity of the xallarap orbit. Has to be in [0, 1) range. + """ + return self.parameters['xi_eccentricity'] + + @xi_eccentricity.setter + def xi_eccentricity(self, new_value): + if new_value < 0. or new_value >= 1.: + raise ValueError('xallarap eccentricity has to be between 0 and 1') + self.parameters['xi_eccentricity'] = new_value + + @property + def xi_omega_periapsis(self): + """ + *float* + + The argument of periapsis of the xallrap orbit, i.e., the angle + between the ascending node and periapsis measured in + the direction of motion. + The units are degrees. + """ + return self.parameters['xi_omega_periapsis'] + + @xi_omega_periapsis.setter + def xi_omega_periapsis(self, new_value): + self._warn_if_angle_outside_reasonable_range( + new_value, 'xi_omega_periapsis') + self.parameters['xi_omega_periapsis'] = new_value + + @property + def t_0_xi(self): + """ + *float* + + Reference epoch for xallarap orbit. + If not provided, then it defaults to :py:attr:`~t_0`. + """ + if 't_0_xi' not in self.parameters.keys(): + return self.parameters['t_0'] + else: + return self.parameters['t_0_xi'] + + @t_0_xi.setter + def t_0_xi(self, new_value): + self.parameters['t_0_xi'] = new_value + @property def t_0_1(self): """ @@ -1592,7 +1850,7 @@ def n_lenses(self): """ *int* - number of objects in the lens system + Number of objects in the lens system. """ return self._n_lenses @@ -1601,7 +1859,8 @@ def n_sources(self): """ *int* - number of luminous sources; it's possible to be 1 for xallarap model + Number of luminous sources. + It can be be 1 for a xallarap model. """ return self._n_sources @@ -1625,6 +1884,15 @@ def is_external_mass_sheet_with_shear(self): return (('shear_G' in self.parameters.keys()) and (self.parameters['shear_G'] != 0)) + @property + def is_xallarap(self): + """ + *bool* + + Whether the parameters include the xallarap or not. + """ + return self._type['xallarap'] + @property def source_1_parameters(self): """ diff --git a/source/MulensModel/orbits/__init__.py b/source/MulensModel/orbits/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/source/MulensModel/orbits/orbit.py b/source/MulensModel/orbits/orbit.py new file mode 100644 index 00000000..1e62a009 --- /dev/null +++ b/source/MulensModel/orbits/orbit.py @@ -0,0 +1,373 @@ +import numpy as np +import math + +""" +Classes that can be accessed directly: +- Orbit +- OrbitCircular +- OrbitEccentric +""" + + +class Orbit(object): + """ + Class that combines :py:class:`OrbitCircular` and + :py:class:`OrbitEccentric`. + """ + def __new__(self, **kwargs): + for Class_ in [OrbitCircular, OrbitEccentric]: + try: + new = Class_(**kwargs) + except Exception: + pass + else: + return new + raise RuntimeError("Orbit.__new__() failed") + + +class _OrbitAbstract(object): + """ + Abstract class for orbits. + """ + def _check_circular_orbit_parameters(self, semimajor_axis): + """ + Check if period and semimajor axis make physical sense. + """ + if self._period <= 0.: + raise ValueError('Orbital period has to be positive.\n' + 'Provided value: ' + str(self._period)) + if semimajor_axis <= 0.: + raise ValueError('Semimajor axis has to be positive.\n' + 'Provided value: ' + str(semimajor_axis)) + + def _check_for_and_get_periapsis_epoch( + self, periapsis_epoch, argument_of_latitude_reference, + epoch_reference): + """ + Check if arguments properly define epoch of + the periapsis (eccentric orbit) or + the ascending node (circular orbit) passage + """ + if periapsis_epoch is not None: + if argument_of_latitude_reference is not None: + raise RuntimeError( + "periapsis_epoch and argument_of_latitude_reference " + "cannot be both set") + if epoch_reference is not None: + raise RuntimeError( + "periapsis_epoch and epoch_reference cannot be both set") + return periapsis_epoch + + if argument_of_latitude_reference is None or epoch_reference is None: + raise RuntimeError("Not enough arguments to define the epoch of " + "periapsis/ascending node") + + u_reference = argument_of_latitude_reference * np.pi / 180. + return self._get_periapsis_epoch(u_reference, epoch_reference) + + def _set_circular_orbit_parameters(self, period, semimajor_axis, + Omega_node, inclination, + periapsis_epoch): + """ + Set parameters that are used for circular orbits. + """ + self._semimajor_axis = semimajor_axis + Omega = math.pi * Omega_node / 180. + self._rotation_matrix = np.array([[math.cos(Omega), -math.sin(Omega)], + [math.sin(Omega), math.cos(Omega)]]) + self._cos_inclination = math.cos(math.pi * inclination / 180.) + self._periapsis_epoch = periapsis_epoch + + def get_reference_plane_position(self, time): + """ + Calculate positions for given epochs in the reference plane. + + Parameters : + time: *float* or *np.ndarray* + Epochs for which positions are requested. + + Returns : + positions: *np.ndarray* + Positions of the body at given epochs. + """ + projected_position = self._get_projected_position(time) + position = np.matmul(self._rotation_matrix, projected_position) + return position + + def _get_projected_position(self, time): + """ + Get position projected on reference plane, + but not rotated to the reference coordinate system + """ + position = self.get_orbital_plane_position(time) + position[1, ] *= self._cos_inclination + return position + + def _get_eccentric_anomaly(self, time): + """ + Calculate eccentric anomaly (typically indicated by E). + """ + mean_anomaly = self._get_mean_anomaly(time) + anomaly = self._get_normalized_anomaly_minus_pi_pi(mean_anomaly) + eccentric_anomaly = ( + self._get_eccentric_anomaly_from_normalized_mean_anomaly(anomaly)) + return eccentric_anomaly + + def _get_mean_anomaly(self, time): + """ + Calculate mean anomaly, i.e., the one that is linear in time and + typically indicated by M. + """ + anomaly = 2. * math.pi * (time - self._periapsis_epoch) / self._period + return anomaly + + def _get_normalized_anomaly_minus_pi_pi(self, anomaly): + """ + get value normalized to (-pi, pi) range + """ + out = self._get_normalized_anomaly_zero_two_pi(anomaly + np.pi) - np.pi + return out + + def _get_normalized_anomaly_zero_two_pi(self, anomaly): + """ + get value normalized to (0, 2*pi) range + """ + return np.remainder(anomaly, 2 * np.pi) + + def _get_eccentric_anomaly_from_normalized_mean_anomaly(self, + mean_anomaly): + """ + Turn mean anomaly in range (-pi, pi) into eccentric anomaly + """ + anomaly = mean_anomaly + 0. + if self._eccentricity > 0.95: # This limit was manually found. + anomaly = np.pi * np.sign(mean_anomaly) + for _ in range(5): + anomaly += self._get_anomaly_correction(anomaly, mean_anomaly) + return anomaly + + def _get_anomaly_correction(self, anomaly, mean_anomaly): + """ + Calculations needed to solve Kepler's equation. + We use Newton's method. + The input anomaly is current estimate of eccentric anomaly. + """ + numerator = ( + anomaly - self._eccentricity * np.sin(anomaly) - mean_anomaly) + denominator = 1. - self._eccentricity * np.cos(anomaly) + return -numerator / denominator + + +class OrbitCircular(_OrbitAbstract): + """ + Circular orbit. + + Keywords : + period: *float* + Orbital period of binary in days + + semimajor_axis: *float* + Semimajor axis of the orbit. The unit is not specified. + Note that the positions returned by + :py:func:`get_orbital_plane_position()` and + :py:func:`get_reference_plane_position()` + functions will be in the same units. + + Omega_node: *float* + Longitude of the ascending node, i.e., the angle from + the reference direction to the ascending node direction. + + inclination: *float* + Inclination of the orbit relative to plane of the sky. + + ascending_node_epoch: *float* or *None* + Epoch when body is in the ascending node. + It's in days and usually you want to provide full BJD or HJD. + + argument_of_latitude_reference: *float* or *None* + Argument of latitude (i.e., u = omega + nu(t_ref)) for + *epoch_reference*, which together define + *ascending_node_epoch* (omega). + + epoch_reference: *float* or *None* + Reference epoch that together with + *argument_of_latitude_reference* defines + *ascending_node_epoch* (omega). + """ + def __init__(self, period, semimajor_axis, Omega_node, inclination, + ascending_node_epoch=None, + argument_of_latitude_reference=None, epoch_reference=None): + self._period = period + self._check_circular_orbit_parameters(semimajor_axis) + ascending_node_epoch = self._check_for_and_get_periapsis_epoch( + ascending_node_epoch, argument_of_latitude_reference, + epoch_reference) + self._set_circular_orbit_parameters( + period, semimajor_axis, Omega_node, + inclination, ascending_node_epoch) + + def _get_periapsis_epoch(self, u_reference, epoch_reference): + """ + Calculate ascending node epoch + (called periapis in general case of eccentric orbits) based on + the argument_of_latitude (u) at given epoch + """ + time_shift = self._period * u_reference / (2. * np.pi) + return epoch_reference - time_shift + + def get_orbital_plane_position(self, time): + """ + Calculate positions in the orbital plane for given epochs + + Parameters : + time: *float* or *np.ndarray* + Epochs for which positions are requested. + + Returns : + positions: *np.ndarray* + Calculated positions. + """ + anomaly = self._get_mean_anomaly(time) + unit_vector = np.array([np.cos(anomaly), np.sin(anomaly)]) + return self._semimajor_axis * unit_vector + + +class OrbitEccentric(_OrbitAbstract): + """ + Eccentric orbit. + + Keywords : + period: *float* + Orbital period of binary. + + semimajor_axis: *float* + Semimajor axis of the orbit. The unit is not specified. + Note that the positions returned by + :py:func:`get_orbital_plane_position()` and + :py:func:`get_reference_plane_position()` + functions will be in the same units. + + Omega_node: *float* + Longitude of the ascending node, i.e., the angle from + the reference direction to the ascending node direction. + + inclination: *float* + Inclination of the orbit relative to plane of the sky. + + eccentricity: *float* + Eccentricity of the orbit, has to be in (0, 1) range. + + omega_periapsis: *float* + Argument of periapsis in degrees. + + periapsis_epoch: *float* or *None* + Epoch when body is in periapsis. + It's in days and usually you want to provide full BJD or HJD. + + argument_of_latitude_reference: *float* or *None* + Argument of latitude (i.e., u = omega + nu(t_ref)) for + *epoch_reference*, which together define + *periapsis_epoch* (omega). + + epoch_reference: *float* or *None* + Reference epoch that together with + *argument_of_latitude_reference* defines + *periapsis_epoch* (omega). + + """ + def __init__( + self, period, semimajor_axis, Omega_node, inclination, + eccentricity, omega_periapsis, periapsis_epoch=None, + argument_of_latitude_reference=None, epoch_reference=None): + self._period = period + self._omega_periapsis = omega_periapsis * np.pi / 180. + self._eccentricity = eccentricity + self._check_circular_orbit_parameters(semimajor_axis) + periapsis_epoch = self._check_for_and_get_periapsis_epoch( + periapsis_epoch, argument_of_latitude_reference, epoch_reference) + self._set_circular_orbit_parameters( + period, semimajor_axis, Omega_node, inclination, periapsis_epoch) + + def _get_periapsis_epoch(self, u_reference, epoch_reference): + """ + Calculate periapsis epoch (omega) based on + the argument_of_latitude (u) at given reference epoch + """ + true_anomaly = u_reference - self._omega_periapsis + mean_anomaly = self._get_mean_anomaly_from_true_anomaly(true_anomaly) + mean_anomaly = self._get_normalized_anomaly_minus_pi_pi(mean_anomaly) + time_shift = self._period * mean_anomaly / (2. * np.pi) + return epoch_reference - time_shift + + def _get_mean_anomaly_from_true_anomaly(self, true_anomaly): + """ + Calculate mean anomaly (M) based on true anomaly (E) + """ + (sin_E, cos_E) = self._get_sin_cos_eccentric_anomaly(true_anomaly) + E = np.arctan2(sin_E, cos_E) + return E - self._eccentricity * sin_E + + def _get_sin_cos_eccentric_anomaly(self, true_anomaly): + """ + Calculate sin(E) and cos(E) based on true_anomaly (nu) + """ + cos_nu = np.cos(true_anomaly) + denominator = 1 + self._eccentricity * cos_nu + sin_E = (np.sqrt(1 - self._eccentricity**2) * np.sin(true_anomaly) + / denominator) + cos_E = (self._eccentricity + cos_nu) / denominator + return (sin_E, cos_E) + + def get_orbital_plane_position(self, time): + """ + Calculate positions in the orbital plane for given epochs + + Parameters : + time: *float* or *np.ndarray* + Epochs for which positions are requested. + + Returns : + positions: *np.ndarray* + Calculated positions. + """ + eccentric_anomaly = self._get_eccentric_anomaly(time) + out_x = np.cos(eccentric_anomaly) - self._eccentricity + out_y = np.sqrt(1 - self._eccentricity**2) * np.sin(eccentric_anomaly) + return self._semimajor_axis * np.array([out_x, out_y]) + + def get_true_anomaly_deg(self, time): + """ + Calculate true anomaly [deg] for given epochs. + + Parameteres : + time: *float* or *np.ndarray* + Epochs for which positions are requested. + + Returns : + true_anomaly: *float* or *np.ndarray* + Values of true anomaly (nu) for given epochs. + The results are in 0-360 range. + """ + true_anomaly = self._get_true_anomaly(time) + true_anomaly_wrapped = np.remainder(true_anomaly, 2*np.pi) + return true_anomaly_wrapped * (180. / np.pi) + + def _get_true_anomaly(self, time): + """ + Calculate true anomaly for given times. + The result is in (-pi, pi) range. + """ + eccentric_anomaly = self._get_eccentric_anomaly(time) + (sin_nu, cos_nu) = self._get_sin_cos_true_anomaly(eccentric_anomaly) + return np.arctan2(sin_nu, cos_nu) + + def _get_sin_cos_true_anomaly(self, eccentric_anomaly): + """ + Calculate sin(nu) and cos(nu) based on E + """ + cos_E = np.cos(eccentric_anomaly) + denominator = 1. - self._eccentricity * cos_E + cos_nu = (cos_E - self._eccentricity) / denominator + sin_nu = (np.sqrt(1-self._eccentricity**2) * np.sin(eccentric_anomaly) + / denominator) + return (sin_nu, cos_nu) diff --git a/source/MulensModel/tests/test_FitData.py b/source/MulensModel/tests/test_FitData.py index ca681b98..3fcf2769 100644 --- a/source/MulensModel/tests/test_FitData.py +++ b/source/MulensModel/tests/test_FitData.py @@ -621,6 +621,22 @@ def test_photfmt_scaled_2(self): almost(res_errors, self.dataset.err_mag) +class TestGradient(unittest.TestCase): + def test_no_gradient_for_xallarap(self): + """ + Make sure that gradient for xallarap models in not implemented. + """ + data = mm.MulensData(file_name=SAMPLE_FILE_02) + model = mm.Model({ + 't_0': 2456836.22, 'u_0': 0.922, 't_E': 22.87, + 'xi_period': 100., 'xi_semimajor_axis': 0.5, 'xi_Omega_node': 90., + 'xi_inclination': 90., 'xi_argument_of_latitude_reference': 90.}) + fit = mm.FitData(model, data) + + with self.assertRaises(NotImplementedError): + fit.get_chi2_gradient(['t_0', 'u_0', 't_E']) + + # Tests to add: # # test get_chi2_gradient(), chi2_gradient: diff --git a/source/MulensModel/tests/test_Model.py b/source/MulensModel/tests/test_Model.py index 8f20beab..07780e70 100644 --- a/source/MulensModel/tests/test_Model.py +++ b/source/MulensModel/tests/test_Model.py @@ -49,7 +49,7 @@ def test_model_init_1(): class TestModel(unittest.TestCase): def test_negative_t_E(self): with self.assertRaises(ValueError): - my_model = mm.Model({'t_0': 2450000., 'u_0': 0.1, 't_E': -100.}) + mm.Model({'t_0': 2450000., 'u_0': 0.1, 't_E': -100.}) def test_model_parallax_definition(): @@ -373,7 +373,7 @@ def test_model_binary_and_finite_sources(): model_1.set_magnification_methods([t1, finite, t2]) model_2.set_magnification_methods([t3, finite, t4]) - (f_s_1, f_s_2, f_b) = (100., 300., 50.) + (f_s_1, f_s_2) = (100., 300.) time = np.linspace(4900., 5200., 4200) mag_1 = model_1.get_magnification(time) mag_2 = model_2.get_magnification(time) @@ -400,16 +400,16 @@ def test_binary_source_and_fluxes_for_bands(): times_I = np.linspace(4900., 5200., 3000) times_V = np.linspace(4800., 5300., 250) - (f_s_1_I, f_s_2_I, f_b_I) = (10., 20., 3.) - (f_s_1_V, f_s_2_V, f_b_V) = (15., 5., 30.) + (f_s_1_I, f_s_2_I) = (10., 20.) + (f_s_1_V, f_s_2_V) = (15., 5.) q_f_I = f_s_2_I / f_s_1_I q_f_V = f_s_2_V / f_s_1_V (mag_1_I, mag_2_I) = model.get_magnification(times_I, separate=True) (mag_1_V, mag_2_V) = model.get_magnification(times_V, separate=True) effective_mag_I = (mag_1_I + mag_2_I * q_f_I) / (1. + q_f_I) effective_mag_V = (mag_1_V + mag_2_V * q_f_V) / (1. + q_f_V) - flux_I = mag_1_I * f_s_1_I + mag_2_I * f_s_2_I + f_b_I - flux_V = mag_1_V * f_s_1_V + mag_2_V * f_s_2_V + f_b_V + # flux_I = mag_1_I * f_s_1_I + mag_2_I * f_s_2_I + f_b_I + # flux_V = mag_1_V * f_s_1_V + mag_2_V * f_s_2_V + f_b_V # model.set_source_flux_ratio_for_band('I', q_f_I) # model.set_source_flux_ratio_for_band('V', q_f_V) @@ -486,6 +486,150 @@ def test_repr(): assert str(model) == expected +def prepare_xallarap_test( + xi_Omega_node=0., xi_argument_of_latitude_reference=0., + t_0_xi=None, xi_eccentricity=None, xi_omega_periapsis=None): + """ + prepare data for unit tests of xallarap models + """ + t_0 = 2459876.0 + d_time = 4. + u_0 = 0.01357 + xi_a = 0.0123456789 + tau = 0.03838383 + + common = {'t_0': t_0, 't_E': d_time / tau, 'u_0': u_0} + xallarap = { + 'xi_period': 2.*d_time, 'xi_semimajor_axis': xi_a, + 'xi_Omega_node': xi_Omega_node, 'xi_inclination': 0., + 'xi_argument_of_latitude_reference': xi_argument_of_latitude_reference + } + if t_0_xi is not None: + if t_0_xi == "t_0": + t_0_xi = t_0 + elif t_0_xi == "t_0+d_time": + t_0_xi = t_0 + d_time + xallarap['t_0_xi'] = t_0_xi + if xi_eccentricity is not None: + xallarap['xi_eccentricity'] = xi_eccentricity + if xi_omega_periapsis is not None: + xallarap['xi_omega_periapsis'] = xi_omega_periapsis + + model_1 = mm.Model(common) + model_2 = mm.Model({**common, **xallarap}) + return (model_1, model_2, t_0, d_time, tau, u_0, xi_a) + + +def test_xallarap_at_t_0_circular(): + """ + Make sure that xallarap circular and non-xallarap 1L1S models + produce the same magnifications at t_0. + """ + (model_1, model_2, t_0) = prepare_xallarap_test()[:2+1] + + assert model_1.get_magnification(t_0) == model_2.get_magnification(t_0) + + +def test_xallarap_at_t_0_eccentric(): + """ + Make sure that xallarap eccentric and non-xallarap 1L1S models + produce the same magnifications at t_0. + """ + kwargs = {'xi_eccentricity': 0.5, 'xi_omega_periapsis': 12.3456} + (model_1, model_2, t_0) = prepare_xallarap_test(**kwargs)[:2+1] + + assert model_1.get_magnification(t_0) == model_2.get_magnification(t_0) + + +def test_xallarap_at_t_0_plus_half_of_period_1(): + """ + Xallarap - circular orbit, half period after t_0, and Omega+nu_0 = 0. + Expected u is from pen and pencil calculations. + """ + (model, t_0, d_time, tau, u_0, xi_a) = prepare_xallarap_test()[1:] + + u2 = u_0**2 + (tau - 2. * xi_a)**2 + expected = (u2 + 2.) / np.sqrt(u2 * (u2 + 4.)) + assert expected == model.get_magnification(t_0+d_time) + + +def test_xallarap_at_t_0_plus_half_of_period_2(): + """ + Xallarap - circular orbit, half period after t_0, and Omega+nu_0 = 90. + Expected u is from pen and pencil calculations. + """ + (model, t_0, d_time, tau, u_0, xi_a) = prepare_xallarap_test(90., 0.)[1:] + + u2 = (u_0 - 2. * xi_a)**2 + tau**2 + expected = (u2 + 2.) / np.sqrt(u2 * (u2 + 4.)) + assert expected == model.get_magnification(t_0+d_time) + + +def test_xallarap_at_t_0_plus_half_of_period_3(): + """ + Xallarap - circular orbit, half period after t_0, and Omega+nu_0 = 180. + Expected u is from pen and pencil calculations. + """ + (model, t_0, d_time, tau, u_0, xi_a) = prepare_xallarap_test(90., 90.)[1:] + + u2 = u_0**2 + (tau + 2. * xi_a)**2 + expected = (u2 + 2.) / np.sqrt(u2 * (u2 + 4.)) + assert expected == model.get_magnification(t_0+d_time) + + +def test_xallarap_at_t_0_plus_half_of_period_4(): + """ + Xallarap - circular orbit, half period after t_0, and Omega+nu_0 = 180. + The t_0_xi is provided as a parameter and = t_0. + Expected u is from pen and pencil calculations. + """ + args = [90., 90., "t_0"] + (model, t_0, d_time, tau, u_0, xi_a) = prepare_xallarap_test(*args)[1:] + + u2 = u_0**2 + (tau + 2. * xi_a)**2 + expected = (u2 + 2.) / np.sqrt(u2 * (u2 + 4.)) + assert expected == model.get_magnification(t_0+d_time) + + +def test_xallarap_at_t_0_plus_half_of_period_5(): + """ + Xallarap - circular orbit, half period after t_0, and Omega+nu_0 = 180. + The t_0_xi is provided as a parameter and = t_0. + Expected u is from pen and pencil calculations. + """ + args = [90., 90., "t_0+d_time"] + (model, t_0, d_time, tau, u_0, xi_a) = prepare_xallarap_test(*args)[1:] + + u2 = u_0**2 + tau**2 + expected = (u2 + 2.) / np.sqrt(u2 * (u2 + 4.)) + assert expected == model.get_magnification(t_0+d_time) + + +def test_xallarap_at_t_0_plus_half_of_period_6_eccentric(): + """ + Extremely eccentric xallarap orbit checked half period later + """ + kwargs = {'xi_eccentricity': 0.999, 'xi_omega_periapsis': 0.} + (model, t_0, d_time, tau, u_0, xi_a) = prepare_xallarap_test(**kwargs)[1:] + + u2 = u_0**2 + (tau - 2. * xi_a)**2 + expected = (u2 + 2.) / np.sqrt(u2 * (u2 + 4.)) + assert expected == model.get_magnification(t_0+d_time) + + +def test_xallarap_at_t_0_plus_half_of_period_7_eccentric(): + """ + Extremely eccentric xallarap orbit with Omega=270 checked half period later + """ + kwargs = {'xi_eccentricity': 0.999, 'xi_omega_periapsis': 0., + 'xi_Omega_node': 270.} + (model, t_0, d_time, tau, u_0, xi_a) = prepare_xallarap_test(**kwargs)[1:] + + u2 = (u_0 + 2. * xi_a)**2 + tau**2 + expected = (u2 + 2.) / np.sqrt(u2 * (u2 + 4.)) + almost(expected, model.get_magnification(t_0+d_time)) + + # Tests to Add: # # test get_trajectory: diff --git a/source/MulensModel/tests/test_ModelParameters.py b/source/MulensModel/tests/test_ModelParameters.py index ad6e1af4..5eaabb5e 100644 --- a/source/MulensModel/tests/test_ModelParameters.py +++ b/source/MulensModel/tests/test_ModelParameters.py @@ -1,4 +1,6 @@ import unittest +import warnings +import pytest import numpy as np from astropy import units as u @@ -440,3 +442,149 @@ def test_failing_single_lens_with_mass_sheet(self): parameters.convergence_K with self.assertRaises(KeyError): parameters.shear_G = 0. + + +xallarap_parameters = { + 't_0': 0, 't_E': 9., 'u_0': 0.1, 'xi_period': 12.345, + 'xi_semimajor_axis': 0.54321, 'xi_Omega_node': 0.123, + 'xi_inclination': 9.8765, 'xi_argument_of_latitude_reference': 24.68, + 'xi_eccentricity': 0.5, 'xi_omega_periapsis': 12.3456, 't_0_xi': 1.} + + +def setup_xallarap(key): + """ + Setup for xallarap tests. + """ + model = mm.ModelParameters(xallarap_parameters) + return (model, xallarap_parameters[key]) + + +tested_keys_1 = ['xi_period', 'xi_semimajor_axis', 'xi_Omega_node', + 'xi_inclination', 'xi_argument_of_latitude_reference', + 'xi_eccentricity', 'xi_omega_periapsis'] +tested_keys_2 = tested_keys_1 + ['t_0_xi'] + + +@pytest.mark.parametrize("key", tested_keys_2) +def test_xallarap_set_value_1(key): + """ + Check if xallarap settings of xallarap are correctly changed via dictionary + """ + (model, value) = setup_xallarap(key) + new_value = 2.1234 * value + model.parameters[key] = new_value + assert model.parameters[key] == new_value + assert getattr(model, key) == new_value + + +@pytest.mark.parametrize("key", tested_keys_2) +def test_xallarap_set_value_2(key): + """ + Check if xallarap settings are correctly changed via attribute + """ + (model, value) = setup_xallarap(key) + new_value = 1.3579 * value + setattr(model, key, new_value) + assert model.parameters[key] == new_value + assert getattr(model, key) == new_value + + +class TestXallarapErrors(unittest.TestCase): + def test_missing_xallarap_parameters(self): + """Make sure that the required xallarap parameters are all defined""" + for parameter in tested_keys_1: + parameters = {**xallarap_parameters} + parameters.pop(parameter) + with self.assertRaises(KeyError): + mm.ModelParameters(parameters) + print("Test failed (i.e. KeyError was not raised) for ", + parameter) + + def test_negative_period(self): + """ + Make sure that period is positive + """ + key = 'xi_period' + (model, value) = setup_xallarap(key) + new_value = -3.14 * value + with self.assertRaises(ValueError): + setattr(model, key, new_value) + + def test_negative_semimajor_axis(self): + """ + Make sure that period is positive + """ + key = 'xi_semimajor_axis' + (model, value) = setup_xallarap(key) + new_value = -3.14 * value + with self.assertRaises(ValueError): + setattr(model, key, new_value) + + def test_xallarap_and_binary_source(self): + """ + Confirm that binary source and xallrap cannot be defined in + the same model + """ + parameters = { + 't_0_1': 0, 'u_0_1': 1, 't_0_2': 5, 'u_0_2': 0.1, 't_E': 9, + 'xi_period': 12., 'xi_semimajor_axis': 0.5, 'xi_Omega_node': 10., + 'xi_inclination': 50., 'xi_argument_of_latitude_reference': 200.} + with self.assertRaises(NotImplementedError): + mm.ModelParameters(parameters) + + def test_xallarap_and_Cassan08(self): + """ + Confirm that xallrap and binary lens in Cassan 2008 parameterization + cannot be defined jointly + """ + parameters = { + 's': 1, 'q': 0.8, 'x_caustic_in': 0.1, 'x_caustic_out': 0.15, + 't_caustic_in': 1000, 't_caustic_out': 2000., + 'xi_period': 12., 'xi_semimajor_axis': 0.5, 'xi_Omega_node': 10., + 'xi_inclination': 50., 'xi_argument_of_latitude_reference': 200.} + with self.assertRaises(NotImplementedError): + mm.ModelParameters(parameters) + + +@pytest.mark.parametrize( + "parameter", + ['xi_Omega_node', 'xi_inclination', 'xi_argument_of_latitude_reference']) +@pytest.mark.parametrize("value", [-361., 541.]) +def test_warnings_for_xallarap_angles(parameter, value): + """ + Check if xallarap angles in somehow strange range give warning + """ + parameters = { + 't_0': 0, 't_E': 9., 'u_0': 0.1, 'xi_period': 12.345, + 'xi_semimajor_axis': 0.54321, 'xi_Omega_node': 100., + 'xi_inclination': 50., 'xi_argument_of_latitude_reference': 200., + 't_0_xi': 1.} + parameters[parameter] = value + + with warnings.catch_warnings(record=True) as warnings_: + warnings.simplefilter("always") + mm.ModelParameters(parameters) + assert len(warnings_) == 1 + assert issubclass(warnings_[0].category, RuntimeWarning) + + +def test_is_xallarap_1(): + """ + Make sure that .is_xallarap() returns True, when it should. + """ + parameters = { + 't_0': 0, 't_E': 9., 'u_0': 0.1, 'xi_period': 12.345, + 'xi_semimajor_axis': 0.54321, 'xi_Omega_node': 100., + 'xi_inclination': 50., 'xi_argument_of_latitude_reference': 200., + 't_0_xi': 1.} + model_params = mm.ModelParameters(parameters) + assert model_params.is_xallarap + + +def test_is_xallarap_2(): + """ + Check that is_xallarap() returns False for a (static) binary source model. + """ + parameters = {'t_0_1': 0, 'u_0_1': 1, 't_0_2': 5, 'u_0_2': 0.1, 't_E': 9} + model_params = mm.ModelParameters(parameters) + assert not model_params.is_xallarap diff --git a/source/MulensModel/tests/test_Orbit.py b/source/MulensModel/tests/test_Orbit.py new file mode 100644 index 00000000..5c83c502 --- /dev/null +++ b/source/MulensModel/tests/test_Orbit.py @@ -0,0 +1,220 @@ +import numpy as np +from numpy.testing import assert_almost_equal +import unittest + +from MulensModel.orbits.orbit import Orbit, OrbitCircular, OrbitEccentric + + +def test_1_circular(): + """ + circular orbit - simples calculation + """ + orbit = OrbitCircular(period=365., semimajor_axis=1., Omega_node=0., + inclination=0., ascending_node_epoch=0.) + position = orbit.get_orbital_plane_position(time=0.) + assert_almost_equal(position, [1., 0.]) + + +def test_2_circular(): + """ + circular orbit at half period position + """ + orbit = OrbitCircular(200., 1.234, 0., 0., 0.) + position = orbit.get_orbital_plane_position(time=100.) + assert_almost_equal(position, [-1.234, 0.]) + + +def test_3_circular(): + """ + circular orbit, Omega = 180, and half period position + """ + orbit = OrbitCircular(200., 1.234, 180., 0., 0.) + position = orbit.get_reference_plane_position(time=100.) + assert_almost_equal(position, [1.234, 0.]) + + +def test_4_circular(): + """ + circular orbit and non-zero inclination and periapsis epoch + """ + orbit = OrbitCircular(200., 2.345, 0., 60., 123.) + position = orbit.get_reference_plane_position(time=123.+200/4.) + assert_almost_equal(position, [0, 2.345/2.]) + + +def test_5_time_vector_circular(): + """ + circular orbit calculation for multiple periods + """ + n_repeat = 10 + orbit = OrbitCircular(200., 2.345, 0., 60., 123.) + times = (np.arange(n_repeat) + 0.25 - 5.) * 200. + 123. + positions = orbit.get_reference_plane_position(times) + expected = np.repeat([[0.], [2.345/2.]], n_repeat, axis=1) + assert_almost_equal(positions, expected) + + +def test_6_eccentric(): + """ + Eccentric orbit at 1 period + """ + orbit = OrbitEccentric( + period=400., semimajor_axis=5., Omega_node=0., inclination=0., + eccentricity=0.6, omega_periapsis=0., periapsis_epoch=100.) + position = orbit.get_orbital_plane_position(300.) + assert_almost_equal(position, [-8., 0.]) + + +def test_7_time_vector_eccentric(): + """ + Very eccentric orbit at 1 period and half period + """ + orbit = OrbitEccentric(400., 10., 0., 0., 0.999, 0., -100.) + times = np.array([-500, -300]) + positions = orbit.get_orbital_plane_position(times) + expected = np.array([[.01, -19.99], [0., 0.]]) + assert_almost_equal(positions, expected) + + +class TestForWrongValues(unittest.TestCase): + def test_8_negative_period(self): + """ + Check negative period + """ + with self.assertRaises(ValueError): + OrbitCircular(-365., 1., 0., 0., 123.) + + def test_9_negative_semimajor_axis(self): + """ + Test for negative semi-major axis + """ + with self.assertRaises(ValueError): + OrbitCircular(365., -1., 0., 0., 123.) + + +def test_10_true_anomaly_large_eccentricity(): + """ + a few epochs and eccentric orbit + """ + orbit = OrbitEccentric(400., 100., 0., 0., 0.9, 0., 0.) + times = np.array([100., 500., 300., 150., 5.]) + true_anomalies = orbit.get_true_anomaly_deg(times) + value = 167.70030551721663 # Value expected for input of 100 d. + expected = np.array([ + value, value, 360.-value, 174.41306950496173, 101.28599627247006]) + assert_almost_equal(true_anomalies, expected) + + +def test_11_true_anomaly_huge_eccentricity(): + """ + calculate true anomaly for extremely eccentric orbit + """ + orbit = OrbitEccentric(400., 100., 0., 0., 0.999, 0., 0.) + true_anomaly = orbit.get_true_anomaly_deg(5.) + assert_almost_equal(true_anomaly, 173.80472546078212, 3) + + +def test_12_Orbit_class_circular(): + """ + Orbit class and simplest calculation for circular orbit + """ + orbit = Orbit(period=365., semimajor_axis=1., Omega_node=0., + inclination=0., ascending_node_epoch=0.) + position = orbit.get_orbital_plane_position(time=0.) + assert_almost_equal(position, [1., 0.]) + + +def test_13_Orbit_class_eccentric(): + """ + Orbit class and eccentric orbit calculation + """ + orbit = Orbit( + period=400., semimajor_axis=100., Omega_node=0., inclination=0., + eccentricity=0.9, omega_periapsis=0., periapsis_epoch=0.) + true_anomaly = orbit.get_true_anomaly_deg(5.) + assert_almost_equal(true_anomaly, 101.28599627247006) + + +class Test_Orbit_fail(unittest.TestCase): + def test_14_Orbit_not_enough_params(self): + """ + Orbit class with wrong input + """ + with self.assertRaises(RuntimeError): + Orbit(period=10., semimajor_axis=10., Omega_node=0., + inclination=0.) + + +def test_15_OrbitCircular_based_on_argument_of_latitude(): + """ + circular orbit and non-zero argument_of_latitude_reference + """ + orbit = OrbitCircular( + period=365, semimajor_axis=1.5, Omega_node=12.34, inclination=-60., + argument_of_latitude_reference=90, epoch_reference=2456789.01234) + position = orbit.get_orbital_plane_position(time=2456789.01234) + assert_almost_equal(position, [0., 1.5]) + + +class Test_OrbitCircular_fail(unittest.TestCase): + def test_16_added_epoch(self): + """ + missing eccentricity + """ + with self.assertRaises(RuntimeError): + OrbitCircular( + period=123., semimajor_axis=5., Omega_node=90., inclination=0., + ascending_node_epoch=2450000., epoch_reference=2450000.) + + def test_17_added_u(self): + """ + too many arguments for a reference epoch + """ + with self.assertRaises(RuntimeError): + OrbitCircular( + period=123., semimajor_axis=5., Omega_node=90., inclination=0., + ascending_node_epoch=2450000., + argument_of_latitude_reference=90.) + + +def test_18_OrbitEccentric_based_on_argument_of_latitude(): + """ + Eccentric orbit calculation starting from + argument_of_latitude_reference + """ + orbit = OrbitEccentric( + period=365, semimajor_axis=1.5, Omega_node=12.34, inclination=-60., + eccentricity=0.5, omega_periapsis=0., + argument_of_latitude_reference=0, epoch_reference=2456789.01234) + position = orbit.get_orbital_plane_position(2456789.01234) + assert_almost_equal(position, [0.75, 0.]) + + +def test_19_OrbitEccentric_based_on_argument_of_latitude(): + """ + Eccentric orbit calculation starting from non-zero value of + argument_of_latitude_reference + """ + orbit = OrbitEccentric( + period=360., semimajor_axis=1.5, Omega_node=12.34, inclination=-60., + eccentricity=0.5, + omega_periapsis=0., epoch_reference=2456789.01234, + argument_of_latitude_reference=118.81500092699673) + # An independend code says the above value is nu for e=0.5 and t=t_0+P/6. + position = orbit.get_orbital_plane_position(2456789.01234-60.) + assert_almost_equal(position, [0.75, 0.]) + + +def test_20_OrbitEccentric_based_on_argument_of_latitude(): + """ + Eccentric orbit with all argument non-zero and + argument_of_latitude_reference in __init__() + """ + orbit = OrbitEccentric( + period=360., semimajor_axis=1.5, Omega_node=12.34, inclination=-60., + eccentricity=0.5, + omega_periapsis=100., + argument_of_latitude_reference=118.81500092699673+100, + epoch_reference=2456789.01234+60) + position = orbit.get_orbital_plane_position(2456789.01234-180.) + assert_almost_equal(position, [-2.25, 0.]) diff --git a/source/MulensModel/trajectory.py b/source/MulensModel/trajectory.py index cfa629c7..ece50bae 100644 --- a/source/MulensModel/trajectory.py +++ b/source/MulensModel/trajectory.py @@ -7,6 +7,7 @@ from MulensModel import utils from MulensModel.modelparameters import ModelParameters from MulensModel.coordinates import Coordinates +from MulensModel.orbits.orbit import Orbit class Trajectory(object): @@ -151,24 +152,15 @@ def get_xy(self): self.parameters.t_E) vector_u = self.parameters.u_0 * np.ones(self.times.size) - # If parallax is non-zero, apply parallax effects: if self.parameters.pi_E is not None: - if self.coords is None: - raise ValueError("You're trying to calculate trajectory in " + - "a parallax model, but event sky " + - "coordinates were not provided.") - - keys = ['earth_orbital', 'satellite', 'topocentric'] - if set([self.parallax[k] for k in keys]) == set([False]): - raise ValueError( - 'If pi_E value is provided then at least one value ' + - 'of parallax dict has to be True ' + - '(earth_orbital, satellite, or topocentric)') - - self._calculate_delta_N_E() - [delta_tau, delta_u] = self._project_delta() - vector_tau += delta_tau - vector_u += delta_u + shifts = self._get_shifts_parralax() + vector_tau += shifts[0] + vector_u += shifts[1] + + if self.parameters.is_xallarap: + shifts = self._get_shifts_xallarap() + vector_tau += shifts[0] + vector_u += shifts[1] # If 2 lenses, rotate trajectory relative to binary lens axis is_mass_sheet = self.parameters.is_external_mass_sheet_with_shear @@ -196,6 +188,23 @@ def get_xy(self): self._x = vector_x self._y = vector_y + def _get_shifts_parralax(self): + """calculate shifts caused by parallax effect""" + if self.coords is None: + raise ValueError( + "You're trying to calculate trajectory in a parallax model, " + "but event sky coordinates were not provided.") + + keys = ['earth_orbital', 'satellite', 'topocentric'] + if set([self.parallax[k] for k in keys]) == set([False]): + raise ValueError( + 'If pi_E value is provided then at least one value of ' + 'parallax dict has to be True ' + '(earth_orbital, satellite, or topocentric)') + + self._calculate_delta_N_E() + return self._project_delta() + def _calculate_delta_N_E(self): """ Calculate shifts caused by microlensing parallax effect. @@ -303,3 +312,15 @@ def _get_delta_satellite(self): Trajectory._get_delta_satellite_results[index] = delta_satellite return delta_satellite + + def _get_shifts_xallarap(self): + """calculate shifts caused by xallarap effect""" + zip_ = self.parameters.parameters.items() + t_0_xi = self.parameters.t_0_xi + orbit_parameters = { + key[3:]: value for (key, value) in zip_ if key[:3] == "xi_"} + orbit_parameters['epoch_reference'] = t_0_xi + orbit = Orbit(**orbit_parameters) + reference_position = orbit.get_reference_plane_position(t_0_xi) + positions = orbit.get_reference_plane_position(self.times) + return positions - reference_position.reshape((2, 1)) diff --git a/source/MulensModel/version.py b/source/MulensModel/version.py index b15cdcbb..a6b62ff3 100644 --- a/source/MulensModel/version.py +++ b/source/MulensModel/version.py @@ -1 +1 @@ -__version__ = "2.16.7" +__version__ = "2.17.0"