From 5574be544836a557594707a7fca059b3e9c1958d Mon Sep 17 00:00:00 2001 From: Jo Bovy Date: Thu, 16 Mar 2023 16:13:23 -0400 Subject: [PATCH 01/18] Prototype of new surface-of-section functionality that includes integrator, SOS function, and plotting function; also working physical outputs and changes for the individual time arrays that these integrations have --- doc/source/reference/orbit.rst | 3 + doc/source/reference/orbitintsos.rst | 4 + doc/source/reference/orbitplotsos.rst | 4 + doc/source/reference/orbitsos.rst | 4 + galpy/orbit/Orbits.py | 508 +++++++++++++++++++++----- galpy/orbit/integrateFullOrbit.py | 107 ++++++ galpy/util/conversion.py | 15 + 7 files changed, 545 insertions(+), 100 deletions(-) create mode 100644 doc/source/reference/orbitintsos.rst create mode 100644 doc/source/reference/orbitplotsos.rst create mode 100644 doc/source/reference/orbitsos.rst diff --git a/doc/source/reference/orbit.rst b/doc/source/reference/orbit.rst index de9072138..5851cc9dd 100644 --- a/doc/source/reference/orbit.rst +++ b/doc/source/reference/orbit.rst @@ -26,6 +26,7 @@ Plotting animate3d plot plot3d + plotSOS In addition to these methods, any calculable attribute listed below can be plotted versus other attributes using ``plotATTR``, where @@ -66,6 +67,7 @@ Methods helioZ integrate integrate_dxdv + integrate_SOS Jacobi jp jr @@ -92,6 +94,7 @@ Methods rguiding rperi SkyCoord + SOS theta time toLinear diff --git a/doc/source/reference/orbitintsos.rst b/doc/source/reference/orbitintsos.rst new file mode 100644 index 000000000..07a86925e --- /dev/null +++ b/doc/source/reference/orbitintsos.rst @@ -0,0 +1,4 @@ +galpy.orbit.Orbit.integrate_SOS +=============================== + +.. automethod:: galpy.orbit.Orbit.integrate_SOS diff --git a/doc/source/reference/orbitplotsos.rst b/doc/source/reference/orbitplotsos.rst new file mode 100644 index 000000000..27d2f3069 --- /dev/null +++ b/doc/source/reference/orbitplotsos.rst @@ -0,0 +1,4 @@ +galpy.orbit.Orbit.plotSOS +========================= + +.. automethod:: galpy.orbit.Orbit.plotSOS diff --git a/doc/source/reference/orbitsos.rst b/doc/source/reference/orbitsos.rst new file mode 100644 index 000000000..f233369ce --- /dev/null +++ b/doc/source/reference/orbitsos.rst @@ -0,0 +1,4 @@ +galpy.orbit.Orbit.SOS +===================== + +.. automethod:: galpy.orbit.Orbit.SOS diff --git a/galpy/orbit/Orbits.py b/galpy/orbit/Orbits.py index 13d7e01cf..0a33249f8 100644 --- a/galpy/orbit/Orbits.py +++ b/galpy/orbit/Orbits.py @@ -34,9 +34,11 @@ from ..util._optional_deps import (_APY3, _APY_COORD_LOADED, _APY_GE_31, _APY_LOADED, _APY_UNITS, _ASTROQUERY_LOADED, _NUMEXPR_LOADED) -from ..util.conversion import physical_compatible, physical_conversion +from ..util.conversion import (physical_compatible, physical_conversion, + physical_conversion_tuple) from ..util.coords import _K -from .integrateFullOrbit import integrateFullOrbit, integrateFullOrbit_c +from .integrateFullOrbit import (integrateFullOrbit, integrateFullOrbit_c, + integrateFullOrbit_sos) from .integrateLinearOrbit import (_ext_loaded, integrateLinearOrbit, integrateLinearOrbit_c) from .integratePlanarOrbit import (integratePlanarOrbit, @@ -64,6 +66,69 @@ except KeyError: import multiprocessing _NUMCORES= multiprocessing.cpu_count() + +# Plot labeling dictionaries +_labeldict_physical= { + 't':r'$t\ (\mathrm{Gyr})$', + 'R':r'$R\ (\mathrm{kpc})$', + 'vR':r'$v_R\ (\mathrm{km\,s}^{-1})$', + 'vT':r'$v_T\ (\mathrm{km\,s}^{-1})$', + 'z':r'$z\ (\mathrm{kpc})$', + 'vz':r'$v_z\ (\mathrm{km\,s}^{-1})$','phi':r'$\phi$', + 'r':r'$r\ (\mathrm{kpc})$', + 'x':r'$x\ (\mathrm{kpc})$','y':r'$y\ (\mathrm{kpc})$', + 'vx':r'$v_x\ (\mathrm{km\,s}^{-1})$', + 'vy':r'$v_y\ (\mathrm{km\,s}^{-1})$', + 'E':r'$E\,(\mathrm{km}^2\,\mathrm{s}^{-2})$', + 'Ez':r'$E_z\,(\mathrm{km}^2\,\mathrm{s}^{-2})$', + 'ER':r'$E_R\,(\mathrm{km}^2\,\mathrm{s}^{-2})$', + 'Enorm':r'$E(t)/E(0.)$', + 'Eznorm':r'$E_z(t)/E_z(0.)$', + 'ERnorm':r'$E_R(t)/E_R(0.)$', + 'Jacobi':r'$E-\Omega_p\,L\,(\mathrm{km}^2\,\mathrm{s}^{-2})$', + 'Jacobinorm':r'$(E-\Omega_p\,L)(t)/(E-\Omega_p\,L)(0)$' +} +_labeldict_internal= { + 't':r'$t$', + 'R':r'$R$', + 'vR':r'$v_R$', + 'vT':r'$v_T$', + 'z':r'$z$', + 'vz':r'$v_z$', + 'phi':r'$\phi$', + 'r':r'$r$', + 'x':r'$x$', + 'y':r'$y$', + 'vx':r'$v_x$', + 'vy':r'$v_y$', + 'E':r'$E$', + 'Enorm':r'$E(t)/E(0.)$', + 'Ez':r'$E_z$', + 'Eznorm':r'$E_z(t)/E_z(0.)$', + 'ER':r'$E_R$', + 'ERnorm':r'$E_R(t)/E_R(0.)$', + 'Jacobi':r'$E-\Omega_p\,L$', + 'Jacobinorm':r'$(E-\Omega_p\,L)(t)/(E-\Omega_p\,L)(0)$' +} +_labeldict_radec= { + 'ra':r'$\alpha\ (\mathrm{deg})$', + 'dec':r'$\delta\ (\mathrm{deg})$', + 'll':r'$l\ (\mathrm{deg})$', + 'bb':r'$b\ (\mathrm{deg})$', + 'dist':r'$d\ (\mathrm{kpc})$', + 'pmra':r'$\mu_\alpha\ (\mathrm{mas\,yr}^{-1})$', + 'pmdec':r'$\mu_\delta\ (\mathrm{mas\,yr}^{-1})$', + 'pmll':r'$\mu_l\ (\mathrm{mas\,yr}^{-1})$', + 'pmbb':r'$\mu_b\ (\mathrm{mas\,yr}^{-1})$', + 'vlos':r'$v_\mathrm{los}\ (\mathrm{km\,s}^{-1})$', + 'helioX':r'$X\ (\mathrm{kpc})$', + 'helioY':r'$Y\ (\mathrm{kpc})$', + 'helioZ':r'$Z\ (\mathrm{kpc})$', + 'U':r'$U\ (\mathrm{km\,s}^{-1})$', + 'V':r'$V\ (\mathrm{km\,s}^{-1})$', + 'W':r'$W\ (\mathrm{km\,s}^{-1})$' +} + # named_objects file def _named_objects_key_formatting(name): # Remove punctuation, spaces, and make lowercase @@ -908,7 +973,11 @@ def __getitem__(self,key): # Also transfer all attributes related to integration if hasattr(self,'orbit'): integrate_kwargs= {} - integrate_kwargs['t']= self.t + # Single vs. individual time arrays + if len(self.t.shape) < len(self.orbit.shape)-1: + integrate_kwargs['t']= self.t + else: + integrate_kwargs['t']= self.t[flat_indx_array] integrate_kwargs['_integrate_t_asQuantity']= \ self._integrate_t_asQuantity integrate_kwargs['orbit']= \ @@ -1041,6 +1110,40 @@ def turn_physical_on(self,ro=None,vo=None): self._vo= vo return None + @staticmethod + def check_integrator(method): + if method.lower() not in ['odeint', 'leapfrog', 'dop853', 'leapfrog_c', + 'symplec4_c', 'symplec6_c', 'rk4_c', 'rk6_c', + 'dopr54_c', 'dop853_c']: + raise ValueError(f'{method:s} is not a valid `method`') + return None + + @staticmethod + def _check_method_c_compatible(method,pot): + if '_c' in method: + if not ext_loaded or not _check_c(pot): + if ('leapfrog' in method or 'symplec' in method): + method= 'leapfrog' + else: + method= 'odeint' + if not ext_loaded: # pragma: no cover + warnings.warn("Cannot use C integration because C extension not loaded (using %s instead)" % (method), galpyWarning) + else: + warnings.warn("Cannot use C integration because some of the potentials are not implemented in C (using %s instead)" % (method), galpyWarning) + return method + + @staticmethod + def _check_method_dissipative_compatible(method,pot): + if _isDissipative(pot) and ('leapfrog' in method + or 'symplec' in method): + if '_c' in method: + method= 'dopr54_c' + else: + method= 'odeint' + warnings.warn("Cannot use symplectic integration because some of the included forces are dissipative (using non-symplectic integrator %s instead)" % (method), galpyWarning) + return method + + def integrate(self,t,pot,method='symplec4_c',progressbar=True, dt=None,numcores=_NUMCORES, force_map=False): @@ -1089,10 +1192,7 @@ def integrate(self,t,pot,method='symplec4_c',progressbar=True, 2018-12-26 - Written to use OpenMP C implementation - Bovy (UofT) """ - if method.lower() not in ['odeint', 'leapfrog', 'dop853', 'leapfrog_c', - 'symplec4_c', 'symplec6_c', 'rk4_c', 'rk6_c', - 'dopr54_c', 'dop853_c']: - raise ValueError(f'{method:s} is not a valid `method`') + self.check_integrator(method) pot= flatten_potential(pot) _check_potential_dim(self,pot) _check_consistent_units(self,pot) @@ -1118,26 +1218,8 @@ def integrate(self,t,pot,method='symplec4_c',progressbar=True, thispot= pot self.t= numpy.array(t) self._pot= thispot - #First check that the potential has C - if '_c' in method: - if not ext_loaded or not _check_c(self._pot): - if ('leapfrog' in method or 'symplec' in method): - method= 'leapfrog' - else: - method= 'odeint' - if not ext_loaded: # pragma: no cover - warnings.warn("Cannot use C integration because C extension not loaded (using %s instead)" % (method), galpyWarning) - else: - warnings.warn("Cannot use C integration because some of the potentials are not implemented in C (using %s instead)" % (method), galpyWarning) - # Now check that we aren't trying to integrate a dissipative force - # with a symplectic integrator - if _isDissipative(self._pot) and ('leapfrog' in method - or 'symplec' in method): - if '_c' in method: - method= 'dopr54_c' - else: - method= 'odeint' - warnings.warn("Cannot use symplectic integration because some of the included forces are dissipative (using non-symplectic integrator %s instead)" % (method), galpyWarning) + method= self._check_method_c_compatible(method,self._pot) + method= self._check_method_dissipative_compatible(method,self._pot) # Implementation with parallel_map in Python if not '_c' in method or not ext_loaded or force_map: if self.dim() == 1: @@ -1230,6 +1312,139 @@ def integrate(self,t,pot,method='symplec4_c',progressbar=True, galpyWarning) return None + def integrate_SOS(self,psi,pot,surface=None,t0=0., + method='symplec4_c',progressbar=True, + dpsi=None,numcores=_NUMCORES, + force_map=False): + """ + NAME: + + integrate_SOS + + PURPOSE: + + integrate this Orbit instance using an independent variable suitable to creating surfaces-of-section + + INPUT: + + psi - increment angles over which to integrate [increments wrt initial angle] (can be Quantity) + + pot - potential instance or list of instances + + surface= (None) surface to punch through (this has no effect in 3D, where the surface is always z=0, but in 2D it can be 'x' or 'y' for x=0 or y=0) + + t0= (0.) initial time (can be Quantity) + + method = 'odeint' for scipy's odeint + 'leapfrog' for a simple leapfrog implementation + 'leapfrog_c' for a simple leapfrog implementation in C + 'symplec4_c' for a 4th order symplectic integrator in C + 'symplec6_c' for a 6th order symplectic integrator in C + 'rk4_c' for a 4th-order Runge-Kutta integrator in C + 'rk6_c' for a 6-th order Runge-Kutta integrator in C + 'dopr54_c' for a 5-4 Dormand-Prince integrator in C + 'dop853' for a 8-5-3 Dormand-Prince integrator in Python + 'dop853_c' for a 8-5-3 Dormand-Prince integrator in C + + progressbar= (True) if True, display a tqdm progress bar when integrating multiple orbits (requires tqdm to be installed!) + + dpsi - if set, force the integrator to use this basic stepsize; must be an integer divisor of output stepsize (only works for the C integrators that use a fixed stepsize) (can be Quantity) + + numcores - number of cores to use for Python-based multiprocessing (pure Python or using force_map=True); default = OMP_NUM_THREADS + + force_map= (False) if True, force use of Python-based multiprocessing (not recommended) + + OUTPUT: + + None (get the actual orbit using getOrbit()) + + HISTORY: + + 2023-03-16 - Written - Bovy (UofT) + + """ + self.check_integrator(method) + pot= flatten_potential(pot) + _check_potential_dim(self,pot) + _check_consistent_units(self,pot) + # Parse psi + if _APY_LOADED and isinstance(psi,units.Quantity): + psi= conversion.parse_angle(psi) + if _APY_LOADED and not dpsi is None and isinstance(dpsi,units.Quantity): + dpsi= conversion.parse_angle(dpsi) + if _APY_LOADED and isinstance(t0,units.Quantity): + t0= conversion.parse_time(t0,ro=self._ro,vo=self._vo) + self._integrate_t_asQuantity= False + from ..potential import MWPotential + if pot == MWPotential: + warnings.warn("Use of MWPotential as a Milky-Way-like potential is deprecated; galpy.potential.MWPotential2014, a potential fit to a large variety of dynamical constraints (see Bovy 2015), is the preferred Milky-Way-like potential in galpy", + galpyWarning) + if not _check_integrate_dt(psi,dpsi): + raise ValueError('dpsi input (integrator stepsize) for Orbit.integrate must be an integer divisor of the output stepsize') + # Delete attributes for interpolation and rperi etc. determination + if hasattr(self,'_orbInterp'): delattr(self,'_orbInterp') + if hasattr(self,'rs'): delattr(self,'rs') + if self.dim() == 2: + thispot= toPlanarPotential(pot) + else: + thispot= pot + self._psi= numpy.array(psi) + self._pot= thispot + method= self._check_method_c_compatible(method,self._pot) + method= self._check_method_dissipative_compatible(method,self._pot) + # Implementation with parallel_map in Python + if not '_c' in method or not ext_loaded or force_map: + if self.dim() == 1: + raise NotImplementedError("SOS integration not implemented for 1D orbits") + out, msg= integrateLinearOrbit(self._pot,self.vxvv,t,method, + progressbar=progressbar, + numcores=numcores,dt=dt) + elif self.dim() == 2: + raise NotImplementedError("SOS integration not implemented for 2D orbits") + out, msg= integratePlanarOrbit(self._pot,self.vxvv,t,method, + progressbar=progressbar, + numcores=numcores,dt=dt) + else: + out, msg= integrateFullOrbit_sos(self._pot,self.vxvv,psi,t0,method, + progressbar=progressbar, + numcores=numcores,dpsi=dpsi) + else: + raise NotImplementedError("SOS integration not implemented for C-based integration") + warnings.warn("Using C implementation to integrate orbits", + galpyWarningVerbose) + if self.dim() == 1: + out, msg= integrateLinearOrbit_c(self._pot, + numpy.copy(self.vxvv), + t,method, + progressbar=progressbar, + dt=dt) + else: + if self.phasedim() == 3 \ + or self.phasedim() == 5: + #We hack this by putting in a dummy phi=0 + vxvvs= numpy.pad(self.vxvv,((0,0),(0,1)), + 'constant',constant_values=0) + else: + vxvvs= numpy.copy(self.vxvv) + if self.dim() == 2: + out, msg= integratePlanarOrbit_c(self._pot,vxvvs, + t,method, + progressbar=progressbar, + dt=dt) + else: + out, msg= integrateFullOrbit_c(self._pot,vxvvs, + t,method, + progressbar=progressbar, + dt=dt) + + if self.phasedim() == 3 \ + or self.phasedim() == 5: + out= out[:,:,:-1] + # Store orbit internally + self.orbit= out[:,:,:-1] + self.t= out[:,:,-1] + return None + def integrate_dxdv(self,dxdv,t,pot,method='dopr54_c', progressbar=True,dt=None, numcores=_NUMCORES,force_map=False, @@ -4490,6 +4705,70 @@ def SkyCoord(self,*args,**kwargs): z_sun=self._zo*units.kpc, galcen_v_sun=v_sun).T + @physical_conversion_tuple(['position','velocity']) + def SOS(self,pot,ncross=500,surface=None,t0=0., + method='dop853_c',progressbar=True, + numcores=_NUMCORES, + force_map=False,**kwargs): + """ + NAME: + + SOS + + PURPOSE: + + calculate the surface of section of the orbit + + INPUT: + + pot - Potential or list of such instances + + ncross= (500) number of times to cross the surface + + surface= (None) surface to punch through (this has no effect in 3D, where the surface is always z=0, but in 2D it can be 'x' or 'y' for x=0 or y=0) + + t0= (0.) time of the initial condition (can be a Quantity) + + integration keyword arguments: + + method = 'odeint' for scipy's odeint + 'leapfrog' for a simple leapfrog implementation + 'leapfrog_c' for a simple leapfrog implementation in C + 'symplec4_c' for a 4th order symplectic integrator in C + 'symplec6_c' for a 6th order symplectic integrator in C + 'rk4_c' for a 4th-order Runge-Kutta integrator in C + 'rk6_c' for a 6-th order Runge-Kutta integrator in C + 'dopr54_c' for a 5-4 Dormand-Prince integrator in C + 'dop853' for a 8-5-3 Dormand-Prince integrator in Python + 'dop853_c' for a 8-5-3 Dormand-Prince integrator in C + + progressbar= (True) if True, display a tqdm progress bar when integrating multiple orbits (requires tqdm to be installed!) + + numcores - number of cores to use for Python-based multiprocessing (pure Python or using force_map=True); default = OMP_NUM_THREADS + + force_map= (False) if True, force use of Python-based multiprocessing (not recommended) + + OUTPUT: + + (R,vR) for 3D orbits, (y,vy) for 2D orbits when surface=='x', (x,vx) for 2D orbits when surface=='y' + + HISTORY: + + 2023-03-16 - Written - Bovy (UofT) + + """ + if self.dim() == 3: + init_psis= numpy.arctan2(self.z(use_physical=False),self.vz(use_physical=False)) + else: + raise NotImplementedError("SOS not implemented for 1D or 2D orbits") + if numpy.any(numpy.fabs(init_psis) > 1e-10): + raise RuntimeError("Use of SOS requires the orbit's initial conditions to be in the plane of the surface of section") + psis= numpy.arange(ncross)*2*numpy.pi + self.integrate_SOS(psis,pot,surface=surface,t0=t0,method=method, + progressbar=progressbar, + numcores=numcores,force_map=force_map) + if self.dim() == 3: + return (self.R(self.t,use_physical=False),self.vR(self.t,use_physical=False)) def __call__(self,*args,**kwargs): """ @@ -4750,8 +5029,10 @@ def _eval(q): # Check those that don't have the exact name of the function if q == 't': # Typically expect this to have same shape as other quantities - return numpy.tile(self.time(self.t,**kwargs), - (len(self.vxvv),1)) + out= self.time(self.t,**kwargs) + if len(self.t.shape) < len(self.orbit.shape)-1: + out= numpy.tile(out,(len(self.vxvv),1)) + return out elif q == 'Enorm': return (self.E(self.t,**kwargs).T/self.E(0.,**kwargs)).T elif q == 'Eznorm': @@ -4827,49 +5108,10 @@ def plot(self,*args,**kwargs): and kwargs.get('ro',self._roSet)) or \ (not 'use_physical' in kwargs \ and kwargs.get('ro',self._roSet)): - labeldict= {'t':r'$t\ (\mathrm{Gyr})$','R':r'$R\ (\mathrm{kpc})$', - 'vR':r'$v_R\ (\mathrm{km\,s}^{-1})$', - 'vT':r'$v_T\ (\mathrm{km\,s}^{-1})$', - 'z':r'$z\ (\mathrm{kpc})$', - 'vz':r'$v_z\ (\mathrm{km\,s}^{-1})$','phi':r'$\phi$', - 'r':r'$r\ (\mathrm{kpc})$', - 'x':r'$x\ (\mathrm{kpc})$','y':r'$y\ (\mathrm{kpc})$', - 'vx':r'$v_x\ (\mathrm{km\,s}^{-1})$', - 'vy':r'$v_y\ (\mathrm{km\,s}^{-1})$', - 'E':r'$E\,(\mathrm{km}^2\,\mathrm{s}^{-2})$', - 'Ez':r'$E_z\,(\mathrm{km}^2\,\mathrm{s}^{-2})$', - 'ER':r'$E_R\,(\mathrm{km}^2\,\mathrm{s}^{-2})$', - 'Enorm':r'$E(t)/E(0.)$', - 'Eznorm':r'$E_z(t)/E_z(0.)$', - 'ERnorm':r'$E_R(t)/E_R(0.)$', - 'Jacobi':r'$E-\Omega_p\,L\,(\mathrm{km}^2\,\mathrm{s}^{-2})$', - 'Jacobinorm':r'$(E-\Omega_p\,L)(t)/(E-\Omega_p\,L)(0)$'} + labeldict= _labeldict_physical.copy() else: - labeldict= {'t':r'$t$','R':r'$R$','vR':r'$v_R$','vT':r'$v_T$', - 'z':r'$z$','vz':r'$v_z$','phi':r'$\phi$', - 'r':r'$r$', - 'x':r'$x$','y':r'$y$','vx':r'$v_x$','vy':r'$v_y$', - 'E':r'$E$','Enorm':r'$E(t)/E(0.)$', - 'Ez':r'$E_z$','Eznorm':r'$E_z(t)/E_z(0.)$', - 'ER':r'$E_R$','ERnorm':r'$E_R(t)/E_R(0.)$', - 'Jacobi':r'$E-\Omega_p\,L$', - 'Jacobinorm':r'$(E-\Omega_p\,L)(t)/(E-\Omega_p\,L)(0)$'} - labeldict.update({'ra':r'$\alpha\ (\mathrm{deg})$', - 'dec':r'$\delta\ (\mathrm{deg})$', - 'll':r'$l\ (\mathrm{deg})$', - 'bb':r'$b\ (\mathrm{deg})$', - 'dist':r'$d\ (\mathrm{kpc})$', - 'pmra':r'$\mu_\alpha\ (\mathrm{mas\,yr}^{-1})$', - 'pmdec':r'$\mu_\delta\ (\mathrm{mas\,yr}^{-1})$', - 'pmll':r'$\mu_l\ (\mathrm{mas\,yr}^{-1})$', - 'pmbb':r'$\mu_b\ (\mathrm{mas\,yr}^{-1})$', - 'vlos':r'$v_\mathrm{los}\ (\mathrm{km\,s}^{-1})$', - 'helioX':r'$X\ (\mathrm{kpc})$', - 'helioY':r'$Y\ (\mathrm{kpc})$', - 'helioZ':r'$Z\ (\mathrm{kpc})$', - 'U':r'$U\ (\mathrm{km\,s}^{-1})$', - 'V':r'$V\ (\mathrm{km\,s}^{-1})$', - 'W':r'$W\ (\mathrm{km\,s}^{-1})$'}) + labeldict= _labeldict_internal.copy() + labeldict.update(_labeldict_radec.copy()) # Cannot be using Quantity output kwargs['quantity']= False #Defaults @@ -4975,36 +5217,10 @@ def plot3d(self,*args,**kwargs): and kwargs.get('ro',self._roSet)) or \ (not 'use_physical' in kwargs \ and kwargs.get('ro',self._roSet)): - labeldict= {'t':r'$t\ (\mathrm{Gyr})$','R':r'$R\ (\mathrm{kpc})$', - 'vR':r'$v_R\ (\mathrm{km\,s}^{-1})$', - 'vT':r'$v_T\ (\mathrm{km\,s}^{-1})$', - 'z':r'$z\ (\mathrm{kpc})$', - 'vz':r'$v_z\ (\mathrm{km\,s}^{-1})$','phi':r'$\phi$', - 'r':r'$r\ (\mathrm{kpc})$', - 'x':r'$x\ (\mathrm{kpc})$','y':r'$y\ (\mathrm{kpc})$', - 'vx':r'$v_x\ (\mathrm{km\,s}^{-1})$', - 'vy':r'$v_y\ (\mathrm{km\,s}^{-1})$'} + labeldict= _labeldict_physical.copy() else: - labeldict= {'t':r'$t$','R':r'$R$','vR':r'$v_R$','vT':r'$v_T$', - 'z':r'$z$','vz':r'$v_z$','phi':r'$\phi$', - 'r':r'$r$','x':r'$x$','y':r'$y$', - 'vx':r'$v_x$','vy':r'$v_y$'} - labeldict.update({'ra':r'$\alpha\ (\mathrm{deg})$', - 'dec':r'$\delta\ (\mathrm{deg})$', - 'll':r'$l\ (\mathrm{deg})$', - 'bb':r'$b\ (\mathrm{deg})$', - 'dist':r'$d\ (\mathrm{kpc})$', - 'pmra':r'$\mu_\alpha\ (\mathrm{mas\,yr}^{-1})$', - 'pmdec':r'$\mu_\delta\ (\mathrm{mas\,yr}^{-1})$', - 'pmll':r'$\mu_l\ (\mathrm{mas\,yr}^{-1})$', - 'pmbb':r'$\mu_b\ (\mathrm{mas\,yr}^{-1})$', - 'vlos':r'$v_\mathrm{los}\ (\mathrm{km\,s}^{-1})$', - 'helioX':r'$X\ (\mathrm{kpc})$', - 'helioY':r'$Y\ (\mathrm{kpc})$', - 'helioZ':r'$Z\ (\mathrm{kpc})$', - 'U':r'$U\ (\mathrm{km\,s}^{-1})$', - 'V':r'$V\ (\mathrm{km\,s}^{-1})$', - 'W':r'$W\ (\mathrm{km\,s}^{-1})$'}) + labeldict= _labeldict_internal.copy() + labeldict.update(_labeldict_radec.copy()) # Cannot be using Quantity output kwargs['quantity']= False #Defaults @@ -5059,6 +5275,98 @@ def plot3d(self,*args,**kwargs): plot._add_ticks() return [line3d] + def plotSOS(self,pot,*args,ncross=500,surface=None,t0=0., + method='dop853_c',progressbar=True, + **kwargs): + """ + NAME: + + plotSOS + + PURPOSE: + + Calculate and plot a surface of section of the orbit + + INPUT: + + pot - Potential or list of such instances + + ncross= (500) number of times to cross the surface + + surface= (None) surface to punch through (this has no effect in 3D, where the surface is always z=0, but in 2D it can be 'x' or 'y' for x=0 or y=0) + + t0= (0.) time of the initial condition (can be a Quantity) + + integration keyword arguments: + + method = 'odeint' for scipy's odeint + 'leapfrog' for a simple leapfrog implementation + 'leapfrog_c' for a simple leapfrog implementation in C + 'symplec4_c' for a 4th order symplectic integrator in C + 'symplec6_c' for a 6th order symplectic integrator in C + 'rk4_c' for a 4th-order Runge-Kutta integrator in C + 'rk6_c' for a 6-th order Runge-Kutta integrator in C + 'dopr54_c' for a 5-4 Dormand-Prince integrator in C + 'dop853' for a 8-5-3 Dormand-Prince integrator in Python + 'dop853_c' for a 8-5-3 Dormand-Prince integrator in C + + progressbar= (True) if True, display a tqdm progress bar when integrating multiple orbits (requires tqdm to be installed!) + + for more control of the integrator, use the SOS method directly and plot its results + + matplotlib.plot inputs+galpy.util.plot.plot inputs + + OUTPUT: + + sends plot to output device + + HISTORY: + + 2023-03-16 - Written - Bovy (UofT) + + """ + if (kwargs.get('use_physical',False) \ + and kwargs.get('ro',self._roSet)) or \ + (not 'use_physical' in kwargs \ + and kwargs.get('ro',self._roSet)): + labeldict= _labeldict_physical.copy() + else: + labeldict= _labeldict_internal.copy() + labeldict.update(_labeldict_radec.copy()) + if self.dim() == 3: + d1= 'R' + d2= 'vR' + x,y= self.SOS(pot,ncross=ncross,surface=surface, + t0=t0,method=method, + progressbar=progressbar,**kwargs) + x= numpy.atleast_2d(x) + y= numpy.atleast_2d(y) + kwargs.pop('ro',None) + kwargs.pop('vo',None) + kwargs.pop('use_physical',None) + kwargs.pop('quantity',None) + auto_scale= not 'xrange' in kwargs and not 'yrange' in kwargs \ + and not kwargs.get('overplot',False) + labels= kwargs.pop('label',[f'Orbit {ii+1}' + for ii in range(self.size)]) + if self.size == 1 and isinstance(labels,str): labels= [labels] + #Plot + if not 'xlabel' in kwargs: + kwargs['xlabel']= labeldict.get(d1,fr'${d1}$') + if not 'ylabel' in kwargs: + kwargs['ylabel']= labeldict.get(d2,fr'${d2}$') + if not 'ls' in kwargs: + kwargs['ls']= 'none' + if not 'marker' in kwargs: + kwargs['marker']= '.' + for ii,(tx,ty) in enumerate(zip(x,y)): + kwargs['label']= labels[ii] + line2d= plot.plot(tx,ty,*args,**kwargs)[0] + kwargs['overplot']= True + if auto_scale: line2d.axes.autoscale(enable=True) + plot._add_ticks() + return [line2d] + def animate(self, **kwargs): #pragma: no cover """ NAME: diff --git a/galpy/orbit/integrateFullOrbit.py b/galpy/orbit/integrateFullOrbit.py index 1d4b13890..9f30e1c76 100644 --- a/galpy/orbit/integrateFullOrbit.py +++ b/galpy/orbit/integrateFullOrbit.py @@ -676,6 +676,88 @@ def integrate_for_map(vxvv): out= out[:,:,:5] return out, numpy.zeros(len(yo)) +def integrateFullOrbit_sos(pot,yo,psi,t0,int_method,rtol=None,atol=None, + numcores=1,progressbar=True,dpsi=None): + """ + NAME: + integrateFullOrbit_sos + PURPOSE: + Integrate an ode for a FullOrbit for integrate_sos + INPUT: + pot - Potential or list of such instances + yo - initial condition [q,p], shape [N,5] or [N,6] + psi - set of increment angles at which one wants the result [increments wrt initial angle] + t0 - initial time + int_method= 'leapfrog', 'odeint', or 'dop853' + rtol, atol= tolerances (not always used...) + numcores= (1) number of cores to use for multi-processing + progressbar= (True) if True, display a tqdm progress bar when integrating multiple orbits (requires tqdm to be installed!) + dpsi= (None) force integrator to use this stepsize (default is to automatically determine one; only for C-based integrators) + OUTPUT: + (y,err) + y : array, shape (N,len(psi),5/6) + Array containing the value of y for each desired angle in psi, \ + with the initial value y0 in the first row. + err: error message, always zero for now + HISTORY: + 2023-03-16 - Written based on integrateFullOrbit - Bovy (UofT) + """ + nophi= False + if len(yo[0]) == 5: + nophi= True + #We hack this by putting in a dummy phi=0 + yo= numpy.pad(yo,((0,0),(0,1)),'constant',constant_values=0) + if not '_c' in int_method: + if rtol is None: rtol= 1e-8 + if int_method.lower() == 'leapfrog': + integrator= symplecticode.leapfrog + extra_kwargs= {'rtol':rtol} + elif int_method.lower() == 'dop853': + integrator= dop853 + extra_kwargs= {} + else: + integrator= integrate.odeint + extra_kwargs= {'rtol':rtol} + def integrate_for_map(vxvv): + #go to the transformed plane: (x,vx,y,vy,A,t) + init_psi= numpy.arctan2(vxvv[3],vxvv[4]) + init= numpy.array([vxvv[0]*numpy.cos(vxvv[5]), + vxvv[1]*numpy.cos(vxvv[5]) + -vxvv[2]*numpy.sin(vxvv[5]), + vxvv[0]*numpy.sin(vxvv[5]), + vxvv[2]*numpy.cos(vxvv[5]) + +vxvv[1]*numpy.sin(vxvv[5]), + numpy.sqrt(vxvv[3]**2.+vxvv[4]**2.), + t0]) + #integrate + intOut= integrator(_SOSEOM,init,t=psi+init_psi,args=(pot,), + **extra_kwargs) + #go back to the cylindrical frame + out= numpy.zeros((len(psi),7)) + out[:,0]= numpy.sqrt(intOut[:,0]**2.+intOut[:,2]**2.) + out[:,5]= numpy.arctan2(intOut[:,2],intOut[:,0]) + out[:,1]= intOut[:,1]*numpy.cos(out[:,5])+intOut[:,3]*numpy.sin(out[:,5]) + out[:,2]= intOut[:,3]*numpy.cos(out[:,5])-intOut[:,1]*numpy.sin(out[:,5]) + out[:,3]= intOut[:,4]*numpy.sin(psi+init_psi) + out[:,4]= intOut[:,4]*numpy.cos(psi+init_psi) + out[:,6]= intOut[:,5] + return out + else: # Assume we are forcing parallel_mapping of a C integrator... + raise NotImplementedError("C-based integrators not implemented for integrateFullOrbit_sos") + def integrate_for_map(vxvv): + return integrateFullOrbit_c(pot,numpy.copy(vxvv), + t,int_method,dt=dt)[0] + if len(yo) == 1: # Can't map a single value... + out= numpy.atleast_3d(integrate_for_map(yo[0]).T).T + else: + out= numpy.array(parallel_map(integrate_for_map,yo,numcores=numcores, + progressbar=progressbar)) + if nophi: + phi_mask= numpy.ones(out.shape[2],dtype='bool') + phi_mask[5]= False + out= out[:,:,phi_mask] + return out, numpy.zeros(len(yo)) + def _RZEOM(y,t,pot,l2): """ NAME: @@ -726,6 +808,31 @@ def _EOM(y,t,pot): _evaluatezforces(pot,y[0],y[4],phi=y[2],t=t, v=[y[1],y[0]*y[3],y[5]])] +def _SOSEOM(y,psi,pot): + """ + NAME: + _SOSEOM + PURPOSE: + implements the EOM, i.e., the right-hand side of the differential + equation, for the SOS integration of a 3D orbit + INPUT: + y - current phase-space position + psi - current angle + pot - (list of) Potential instance(s) + OUTPUT: + dy/dpsi + HISTORY: + 2023-03-16 - Written - Bovy (UofT) + """ + # y = (x,vx,y,vy,A,t) + # Calculate z + sp, cp= numpy.sin(psi), numpy.cos(psi) + z= y[4]*sp + gxyz= _rectForce([y[0],y[2],z],pot,t=y[5]) + psidot= cp**2.-sp/y[4]*gxyz[2] + Adot= y[4]*cp*sp+gxyz[2]*cp + return numpy.array([y[1],gxyz[0],y[3],gxyz[1],Adot,1.])/psidot + def _rectForce(x,pot,t=0.): """ NAME: diff --git a/galpy/util/conversion.py b/galpy/util/conversion.py index 286f844a5..d0954ce56 100644 --- a/galpy/util/conversion.py +++ b/galpy/util/conversion.py @@ -896,6 +896,21 @@ def wrapped(*args,**kwargs): return method(*args,**kwargs) return wrapped return wrapper +def physical_conversion_tuple(quantities,pop=False): + """Decorator to convert to physical coordinates for tuple outputs. + So outputs are a tuple of quantities that each need to be converted, + with possibly different conversions, e.g., (R,vR)""" + def wrapper(method): + @wraps(method) + def wrapped(*args,**kwargs): + rawOut= method(*args,**kwargs) + out= () + for ii in range(len(rawOut)): + # Apply physical conversion by converting a wrapped dummy function that returns the raw output + out= out+(physical_conversion(quantities[ii])(lambda x,**kwargs: rawOut[ii])(args[0],**kwargs),) + return out + return wrapped + return wrapper def potential_physical_input(method): """Decorator to convert inputs to Potential functions from physical to internal coordinates""" From a419f9df54d5eb4ade8cf8d50164d85742c3362b Mon Sep 17 00:00:00 2001 From: Jo Bovy Date: Mon, 20 Mar 2023 10:20:59 -0400 Subject: [PATCH 02/18] Working version of SOS integration in C for 3D orbits --- galpy/orbit/Orbits.py | 70 +++++----- galpy/orbit/integrateFullOrbit.py | 115 +++++++++++++++- galpy/orbit/orbit_c_ext/integrateFullOrbit.c | 134 ++++++++++++++++++- galpy/util/bovy_coords.c | 50 +++++++ galpy/util/bovy_coords.h | 2 + 5 files changed, 329 insertions(+), 42 deletions(-) diff --git a/galpy/orbit/Orbits.py b/galpy/orbit/Orbits.py index 0a33249f8..a51af94f7 100644 --- a/galpy/orbit/Orbits.py +++ b/galpy/orbit/Orbits.py @@ -38,7 +38,8 @@ physical_conversion_tuple) from ..util.coords import _K from .integrateFullOrbit import (integrateFullOrbit, integrateFullOrbit_c, - integrateFullOrbit_sos) + integrateFullOrbit_sos, + integrateFullOrbit_sos_c) from .integrateLinearOrbit import (_ext_loaded, integrateLinearOrbit, integrateLinearOrbit_c) from .integratePlanarOrbit import (integratePlanarOrbit, @@ -1111,10 +1112,19 @@ def turn_physical_on(self,ro=None,vo=None): return None @staticmethod - def check_integrator(method): - if method.lower() not in ['odeint', 'leapfrog', 'dop853', 'leapfrog_c', - 'symplec4_c', 'symplec6_c', 'rk4_c', 'rk6_c', - 'dopr54_c', 'dop853_c']: + def check_integrator(method,no_symplec=False): + valid_methods= [ + 'odeint', 'leapfrog', 'dop853', 'leapfrog_c', + 'symplec4_c', 'symplec6_c', 'rk4_c', 'rk6_c', + 'dopr54_c', 'dop853_c' + ] + if no_symplec: + symplec_methods= [ + 'leapfrog','leapfrog_c', + 'symplec4_c','symplec6_c', + ] + [valid_methods.remove(symplec_method) for symplec_method in symplec_methods] + if method.lower() not in valid_methods: raise ValueError(f'{method:s} is not a valid `method`') return None @@ -1313,7 +1323,7 @@ def integrate(self,t,pot,method='symplec4_c',progressbar=True, return None def integrate_SOS(self,psi,pot,surface=None,t0=0., - method='symplec4_c',progressbar=True, + method='dop853_c',progressbar=True, dpsi=None,numcores=_NUMCORES, force_map=False): """ @@ -1336,10 +1346,6 @@ def integrate_SOS(self,psi,pot,surface=None,t0=0., t0= (0.) initial time (can be Quantity) method = 'odeint' for scipy's odeint - 'leapfrog' for a simple leapfrog implementation - 'leapfrog_c' for a simple leapfrog implementation in C - 'symplec4_c' for a 4th order symplectic integrator in C - 'symplec6_c' for a 6th order symplectic integrator in C 'rk4_c' for a 4th-order Runge-Kutta integrator in C 'rk6_c' for a 6-th order Runge-Kutta integrator in C 'dopr54_c' for a 5-4 Dormand-Prince integrator in C @@ -1363,7 +1369,7 @@ def integrate_SOS(self,psi,pot,surface=None,t0=0., 2023-03-16 - Written - Bovy (UofT) """ - self.check_integrator(method) + self.check_integrator(method,no_symplec=True) pot= flatten_potential(pot) _check_potential_dim(self,pot) _check_consistent_units(self,pot) @@ -1409,10 +1415,10 @@ def integrate_SOS(self,psi,pot,surface=None,t0=0., progressbar=progressbar, numcores=numcores,dpsi=dpsi) else: - raise NotImplementedError("SOS integration not implemented for C-based integration") warnings.warn("Using C implementation to integrate orbits", galpyWarningVerbose) if self.dim() == 1: + raise NotImplementedError("SOS integration not implemented for 1D orbits in C") out, msg= integrateLinearOrbit_c(self._pot, numpy.copy(self.vxvv), t,method, @@ -1427,19 +1433,21 @@ def integrate_SOS(self,psi,pot,surface=None,t0=0., else: vxvvs= numpy.copy(self.vxvv) if self.dim() == 2: + raise NotImplementedError("SOS integration not implemented for 2D orbits") out, msg= integratePlanarOrbit_c(self._pot,vxvvs, t,method, progressbar=progressbar, dt=dt) else: - out, msg= integrateFullOrbit_c(self._pot,vxvvs, - t,method, - progressbar=progressbar, - dt=dt) + out, msg= integrateFullOrbit_sos_c(self._pot,vxvvs,psi,t0, + method,progressbar=progressbar, + dpsi=dpsi) if self.phasedim() == 3 \ or self.phasedim() == 5: - out= out[:,:,:-1] + phi_mask= numpy.ones(out.shape[2],dtype='bool') + phi_mask[5]= False + out= out[:,:,phi_mask] # Store orbit internally self.orbit= out[:,:,:-1] self.t= out[:,:,-1] @@ -4732,15 +4740,11 @@ def SOS(self,pot,ncross=500,surface=None,t0=0., integration keyword arguments: method = 'odeint' for scipy's odeint - 'leapfrog' for a simple leapfrog implementation - 'leapfrog_c' for a simple leapfrog implementation in C - 'symplec4_c' for a 4th order symplectic integrator in C - 'symplec6_c' for a 6th order symplectic integrator in C - 'rk4_c' for a 4th-order Runge-Kutta integrator in C - 'rk6_c' for a 6-th order Runge-Kutta integrator in C - 'dopr54_c' for a 5-4 Dormand-Prince integrator in C - 'dop853' for a 8-5-3 Dormand-Prince integrator in Python - 'dop853_c' for a 8-5-3 Dormand-Prince integrator in C + 'rk4_c' for a 4th-order Runge-Kutta integrator in C + 'rk6_c' for a 6-th order Runge-Kutta integrator in C + 'dopr54_c' for a 5-4 Dormand-Prince integrator in C + 'dop853' for a 8-5-3 Dormand-Prince integrator in Python + 'dop853_c' for a 8-5-3 Dormand-Prince integrator in C progressbar= (True) if True, display a tqdm progress bar when integrating multiple orbits (requires tqdm to be installed!) @@ -5300,15 +5304,11 @@ def plotSOS(self,pot,*args,ncross=500,surface=None,t0=0., integration keyword arguments: method = 'odeint' for scipy's odeint - 'leapfrog' for a simple leapfrog implementation - 'leapfrog_c' for a simple leapfrog implementation in C - 'symplec4_c' for a 4th order symplectic integrator in C - 'symplec6_c' for a 6th order symplectic integrator in C - 'rk4_c' for a 4th-order Runge-Kutta integrator in C - 'rk6_c' for a 6-th order Runge-Kutta integrator in C - 'dopr54_c' for a 5-4 Dormand-Prince integrator in C - 'dop853' for a 8-5-3 Dormand-Prince integrator in Python - 'dop853_c' for a 8-5-3 Dormand-Prince integrator in C + 'rk4_c' for a 4th-order Runge-Kutta integrator in C + 'rk6_c' for a 6-th order Runge-Kutta integrator in C + 'dopr54_c' for a 5-4 Dormand-Prince integrator in C + 'dop853' for a 8-5-3 Dormand-Prince integrator in Python + 'dop853_c' for a 8-5-3 Dormand-Prince integrator in C progressbar= (True) if True, display a tqdm progress bar when integrating multiple orbits (requires tqdm to be installed!) diff --git a/galpy/orbit/integrateFullOrbit.py b/galpy/orbit/integrateFullOrbit.py index 9f30e1c76..09908a598 100644 --- a/galpy/orbit/integrateFullOrbit.py +++ b/galpy/orbit/integrateFullOrbit.py @@ -676,6 +676,119 @@ def integrate_for_map(vxvv): out= out[:,:,:5] return out, numpy.zeros(len(yo)) +def integrateFullOrbit_sos_c(pot,yo,psi,t0,int_method,rtol=None,atol=None, + progressbar=True,dpsi=None): + """ + NAME: + integrateFullOrbit_sos_c + PURPOSE: + Integrate an ode for a FullOrbit for integrate_sos in C + INPUT: + pot - Potential or list of such instances + yo - initial condition [q,p], shape [N,5] or [N,6] + psi - set of increment angles at which one wants the result [increments wrt initial angle] + t0 - initial time + int_method= 'leapfrog', 'odeint', or 'dop853' + rtol, atol= tolerances (not always used...) + progressbar= (True) if True, display a tqdm progress bar when integrating multiple orbits (requires tqdm to be installed!) + dpsi= (None) force integrator to use this stepsize (default is to automatically determine one; only for C-based integrators) + OUTPUT: + (y,err) + y : array, shape (N,len(psi),7) where the last of the last dimension is the time + Array containing the value of y for each desired angle in psi, \ + with the initial value y0 in the first row. + err: error message, always zero for now + HISTORY: + 2023-03-17 - Written based on integrateFullOrbit_c - Bovy (UofT) + """ + if len(yo.shape) == 1: single_obj= True + else: single_obj= False + yo= numpy.atleast_2d(yo) + nobj= len(yo) + rtol, atol= _parse_tol(rtol,atol) + npot, pot_type, pot_args, pot_tfuncs= _parse_pot(pot) + pot_tfuncs= _prep_tfuncs(pot_tfuncs) + int_method_c= _parse_integrator(int_method) + if dpsi is None: + dpsi= -9999.99 + t0= numpy.atleast_1d(t0) + yoo= numpy.empty((nobj,7)) + yoo[:,:6]= yo[:,:6] + if len(t0) == 1: + yoo[:,6]= t0[0] + else: + yoo[:,6]= t0 + + #Set up result array + result= numpy.empty((nobj,len(psi),7)) + err= numpy.zeros(nobj,dtype=numpy.int32) + + #Set up progressbar + progressbar*= _TQDM_LOADED + if nobj > 1 and progressbar: + pbar= tqdm.tqdm(total=nobj,leave=False) + pbar_func_ctype= ctypes.CFUNCTYPE(None) + pbar_c= pbar_func_ctype(pbar.update) + else: # pragma: no cover + pbar_c= None + + #Set up the C code + ndarrayFlags= ('C_CONTIGUOUS','WRITEABLE') + integrationFunc= _lib.integrateFullOrbit_sos + integrationFunc.argtypes= [ctypes.c_int, + ndpointer(dtype=numpy.float64,flags=ndarrayFlags), + ctypes.c_int, + ndpointer(dtype=numpy.float64,flags=ndarrayFlags), + ctypes.c_int, + ndpointer(dtype=numpy.int32,flags=ndarrayFlags), + ndpointer(dtype=numpy.float64,flags=ndarrayFlags), + ctypes.c_void_p, + ctypes.c_double, + ctypes.c_double, + ctypes.c_double, + ndpointer(dtype=numpy.float64,flags=ndarrayFlags), + ndpointer(dtype=numpy.int32,flags=ndarrayFlags), + ctypes.c_int, + ctypes.c_void_p] + + #Array requirements, first store old order + f_cont= [yoo.flags['F_CONTIGUOUS'], + psi.flags['F_CONTIGUOUS']] + yoo= numpy.require(yoo,dtype=numpy.float64,requirements=['C','W']) + psi= numpy.require(psi,dtype=numpy.float64,requirements=['C','W']) + result= numpy.require(result,dtype=numpy.float64,requirements=['C','W']) + err= numpy.require(err,dtype=numpy.int32,requirements=['C','W']) + + #Run the C code + integrationFunc(ctypes.c_int(nobj), + yoo, + ctypes.c_int(len(psi)), + psi, + ctypes.c_int(npot), + pot_type, + pot_args, + pot_tfuncs, + ctypes.c_double(dpsi), + ctypes.c_double(rtol), + ctypes.c_double(atol), + result, + err, + ctypes.c_int(int_method_c), + pbar_c) + + if nobj > 1 and progressbar: + pbar.close() + + if numpy.any(err == -10): #pragma: no cover + raise KeyboardInterrupt("Orbit integration interrupted by CTRL-C (SIGINT)") + + #Reset input arrays + if f_cont[0]: yoo= numpy.asfortranarray(yoo) + if f_cont[1]: psi= numpy.asfortranarray(psi) + + if single_obj: return (result[0],err[0]) + else: return (result,err) + def integrateFullOrbit_sos(pot,yo,psi,t0,int_method,rtol=None,atol=None, numcores=1,progressbar=True,dpsi=None): """ @@ -695,7 +808,7 @@ def integrateFullOrbit_sos(pot,yo,psi,t0,int_method,rtol=None,atol=None, dpsi= (None) force integrator to use this stepsize (default is to automatically determine one; only for C-based integrators) OUTPUT: (y,err) - y : array, shape (N,len(psi),5/6) + y : array, shape (N,len(psi),6/7) where the last of the last dimension is the time Array containing the value of y for each desired angle in psi, \ with the initial value y0 in the first row. err: error message, always zero for now diff --git a/galpy/orbit/orbit_c_ext/integrateFullOrbit.c b/galpy/orbit/orbit_c_ext/integrateFullOrbit.c index b337f73c2..5b0cba239 100644 --- a/galpy/orbit/orbit_c_ext/integrateFullOrbit.c +++ b/galpy/orbit/orbit_c_ext/integrateFullOrbit.c @@ -49,6 +49,8 @@ void evalRectForce(double, double *, double *, int, struct potentialArg *); void evalRectDeriv(double, double *, double *, int, struct potentialArg *); +void evalSOSDeriv(double, double *, double *, + int, struct potentialArg *); void evalRectDeriv_dxdv(double,double *, double *, int, struct potentialArg *); void initMovingObjectSplines(struct potentialArg *, double ** pot_args); @@ -726,6 +728,88 @@ EXPORT void integrateFullOrbit(int nobj, free(potentialArgs); //Done! } +EXPORT void integrateFullOrbit_sos( + int nobj, + double *yo, + int npsi, + double *psi, + int npot, + int * pot_type, + double * pot_args, + tfuncs_type_arr pot_tfuncs, + double dpsi, + double rtol, + double atol, + double *result, + int * err, + int odeint_type, + orbint_callback_type cb){ + //Set up the forces, first count + int ii,jj; + int dim; + int max_threads; + int * thread_pot_type; + double * thread_pot_args; + tfuncs_type_arr thread_pot_tfuncs; + max_threads= ( nobj < omp_get_max_threads() ) ? nobj : omp_get_max_threads(); + // Because potentialArgs may cache, safest to have one / thread + struct potentialArg * potentialArgs= (struct potentialArg *) malloc ( max_threads * npot * sizeof (struct potentialArg) ); +#pragma omp parallel for schedule(static,1) private(ii,thread_pot_type,thread_pot_args,thread_pot_tfuncs) num_threads(max_threads) + for (ii=0; ii < max_threads; ii++) { + thread_pot_type= pot_type; // need to make thread-private pointers, bc + thread_pot_args= pot_args; // these pointers are changed in parse_... + thread_pot_tfuncs= pot_tfuncs; // ... + parse_leapFuncArgs_Full(npot,potentialArgs+ii*npot, + &thread_pot_type,&thread_pot_args,&thread_pot_tfuncs); + } + //Integrate + void (*odeint_func)(void (*func)(double, double *, double *, + int, struct potentialArg *), + int, + double *, + int, double, double *, + int, struct potentialArg *, + double, double, + double *,int *); + void (*odeint_deriv_func)(double, double *, double *, + int,struct potentialArg *); + dim= 7; + odeint_deriv_func= &evalSOSDeriv; + switch ( odeint_type ) { + // case 0: = leapfrog = not supported symplectic method + case 1: //RK4 + odeint_func= &bovy_rk4; + break; + case 2: //RK6 + odeint_func= &bovy_rk6; + break; + // case 3: = symplec4 = not supported symplectic method + // case 4: = symplec6 = not supported symplectic method + case 5: //DOPR54 + odeint_func= &bovy_dopr54; + break; + case 6: //DOP853 + odeint_func= &dop853; + break; + } +#pragma omp parallel for schedule(dynamic,ORBITS_CHUNKSIZE) private(ii,jj) num_threads(max_threads) + for (ii=0; ii < nobj; ii++) { + cyl_to_sos_galpy(yo+dim*ii); + odeint_func(odeint_deriv_func,dim,yo+dim*ii,npsi,dpsi,psi, + npot,potentialArgs+omp_get_thread_num()*npot,rtol,atol, + result+dim*npsi*ii,err+ii); + for (jj=0; jj < npsi; jj++) + sos_to_cyl_galpy(result+dim*jj+dim*npsi*ii); + if ( cb ) // Callback if not void + cb(); + } + //Free allocated memory +#pragma omp parallel for schedule(static,1) private(ii) num_threads(max_threads) + for (ii=0; ii < max_threads; ii++) + free_potentialArgs(npot,potentialArgs+ii*npot); + free(potentialArgs); + //Done! +} // LCOV_EXCL_START void integrateOrbit_dxdv(double *yo, int nt, @@ -801,7 +885,7 @@ void integrateOrbit_dxdv(double *yo, // LCOV_EXCL_STOP void evalRectForce(double t, double *q, double *a, int nargs, struct potentialArg * potentialArgs){ - double sinphi, cosphi, x, y, phi,R,Rforce,phitorque, z, zforce; + double sinphi, cosphi, x, y, phi,R,Rforce,phitorque, z; //q is rectangular so calculate R and phi x= *q; y= *(q+1); @@ -813,15 +897,14 @@ void evalRectForce(double t, double *q, double *a, if ( y < 0. ) phi= 2.*M_PI-phi; //Calculate the forces Rforce= calcRforce(R,z,phi,t,nargs,potentialArgs); - zforce= calczforce(R,z,phi,t,nargs,potentialArgs); phitorque= calcphitorque(R,z,phi,t,nargs,potentialArgs); *a++= cosphi*Rforce-1./R*sinphi*phitorque; *a++= sinphi*Rforce+1./R*cosphi*phitorque; - *a= zforce; + *a= calczforce(R,z,phi,t,nargs,potentialArgs); } void evalRectDeriv(double t, double *q, double *a, int nargs, struct potentialArg * potentialArgs){ - double sinphi, cosphi, x, y, phi,R,Rforce,phitorque,z,zforce,vR,vT; + double sinphi, cosphi, x, y, phi,R,Rforce,phitorque,z,vR,vT; //first three derivatives are just the velocities *a++= *(q+3); *a++= *(q+4); @@ -840,11 +923,50 @@ void evalRectDeriv(double t, double *q, double *a, vT= -*(q+3) * sinphi + *(q+4) * cosphi; //Calculate the forces Rforce= calcRforce(R,z,phi,t,nargs,potentialArgs,vR,vT,*(q+5)); - zforce= calczforce(R,z,phi,t,nargs,potentialArgs,vR,vT,*(q+5)); phitorque= calcphitorque(R,z,phi,t,nargs,potentialArgs,vR,vT,*(q+5)); *a++= cosphi*Rforce-1./R*sinphi*phitorque; *a++= sinphi*Rforce+1./R*cosphi*phitorque; - *a= zforce; + *a= calczforce(R,z,phi,t,nargs,potentialArgs,vR,vT,*(q+5));; +} + +void evalSOSDeriv(double psi, double *q, double *a, + int nargs, struct potentialArg * potentialArgs){ + // q= (x,y,vx,vy,A,t,psi); to save operations, we use the first three + // components of q as (x,y,z) and then re-use a first for the + // rectForce then for the actual RHS + // Note also that we keep track of psi in q+6, not in psi! This is + // such that we can avoid having to convert psi to psi+psi0 + // q+6 starts as psi0 and then just increments as psi (exactly) + double sinpsi,cospsi,psidot,x,y,z,R,phi,sinphi,cosphi,vR,vT,vz,Rforce,phitorque; + sinpsi= sin( *(q+6) ); + cospsi= cos( *(q+6) ); + // Calculate forces, put them in a+3, a+4, a+5 + //q is rectangular so calculate R and phi, vR and vT (for dissipative) + x= *q; + y= *(q+1); + z= *(q+4) * sinpsi; + R= sqrt(x*x+y*y); + phi= atan2( y ,x ); + sinphi= y/R; + cosphi= x/R; + vR= *(q+2) * cosphi + *(q+3) * sinphi; + vT= -*(q+2) * sinphi + *(q+3) * cosphi; + vz= *(q+4) * cospsi; + //Calculate the forces + Rforce= calcRforce(R,z,phi,*(q+5),nargs,potentialArgs,vR,vT,vz); + phitorque= calcphitorque(R,z,phi,*(q+5),nargs,potentialArgs,vR,vT,vz); + *(a+3)= cosphi*Rforce-1./R*sinphi*phitorque; + *(a+4)= sinphi*Rforce+1./R*cosphi*phitorque; + *(a+5)= calczforce(R,z,phi,*(q+5),nargs,potentialArgs,vR,vT,vz); + // Now calculate the RHS of the ODE + psidot= cospsi * cospsi - sinpsi * *(a+5) / ( *(q+4) ); + *(a )= *(q+2) / psidot; + *(a+1)= *(q+3) / psidot; + *(a+2)= *(a+3) / psidot; + *(a+3)= *(a+4) / psidot; + *(a+4)= cospsi * ( *(q+4) * sinpsi + *(a+5) ) / psidot; + *(a+5)= 1./psidot; + *(a+6)= 1.; // dpsi / dpsi to keep track of psi } void initMovingObjectSplines(struct potentialArg * potentialArgs, diff --git a/galpy/util/bovy_coords.c b/galpy/util/bovy_coords.c index 2abeb1cfa..d6d718ab5 100644 --- a/galpy/util/bovy_coords.c +++ b/galpy/util/bovy_coords.c @@ -137,3 +137,53 @@ void rect_to_cyl_galpy(double *vxvv){ *(vxvv+1)= vx * cp + vy * sp; *(vxvv+2)= -vx * sp + vy * cp; } +/* +NAME: cyl_to_sos_galpy +PURPOSE: convert (R,vR,vT,z,vz,phi,t) to (x,y,vx,vy,A,t,psi) coordinates for SOS integration +INPUT: + double * vxvv - (R,vR,vT,z,vz,phi,t) +OUTPUT: + performed in-place +HISTORY: 2023-03-19 - Written - Bovy (UofT) + */ +void cyl_to_sos_galpy(double *vxvv){ + double R,cp,sp,vR,vT; + R = *vxvv; + cp = cos ( *(vxvv+5) ); + sp = sin ( *(vxvv+5) ); + vR = *(vxvv+1); + vT = *(vxvv+2); + *(vxvv+5)= *(vxvv+6); + *(vxvv+6)= atan2( *(vxvv+3) , *(vxvv+4) ); + *(vxvv+4)= sqrt( *(vxvv+3) * *(vxvv+3) + *(vxvv+4) * *(vxvv+4) ); + *vxvv = R * cp; + *(vxvv+1)= R * sp; + *(vxvv+2)= vR * cp - vT * sp; + *(vxvv+3)= vR * sp + vT * cp; +} +/* +NAME: sos_to_cyl_galpy +PURPOSE: convert (x,y,vx,vy,A,t,psi) to (R,vR,vT,z,vz,phi,t) +INPUT: + double * vxvv - (x,y,vx,vy,A,t,psi) +OUTPUT (as arguments): + performed in-place +HISTORY: 2023-03-19 - Written - Bovy (UofT) + */ +void sos_to_cyl_galpy(double *vxvv){ + double x,y,vx,vy,phi,cp,sp; + x = *(vxvv ); + y = *(vxvv+1); + vx= *(vxvv+2); + vy= *(vxvv+3); + phi= atan2( y , x ); + cp = cos ( phi ); + sp = sin ( phi ); + *(vxvv )= sqrt ( x * x + y * y ); + *(vxvv+1)= vx * cp + vy * sp; + *(vxvv+2)= -vx * sp + vy * cp; + *(vxvv+3)= *(vxvv+4) * sin ( *(vxvv+6) ); + *(vxvv+4)= *(vxvv+4) * cos ( *(vxvv+6) ); + *(vxvv+6)= *(vxvv+5); + *(vxvv+5)= phi; +} diff --git a/galpy/util/bovy_coords.h b/galpy/util/bovy_coords.h index 8e4088a06..47b2c11c4 100644 --- a/galpy/util/bovy_coords.h +++ b/galpy/util/bovy_coords.h @@ -47,6 +47,8 @@ void polar_to_rect_galpy(double *); void rect_to_polar_galpy(double *); void cyl_to_rect_galpy(double *); void rect_to_cyl_galpy(double *); +void cyl_to_sos_galpy(double *); +void sos_to_cyl_galpy(double *); #ifdef __cplusplus } #endif From e8c4e168f5a89c9b5faa9364a302f9924e7512d2 Mon Sep 17 00:00:00 2001 From: Jo Bovy Date: Wed, 22 Mar 2023 13:55:37 -0400 Subject: [PATCH 03/18] Allow individual psi integration variables per orbit; don't support dpsi for SOS integration --- galpy/orbit/Orbits.py | 32 +++++++++++--------- galpy/orbit/integrateFullOrbit.py | 31 ++++++++++++------- galpy/orbit/orbit_c_ext/integrateFullOrbit.c | 25 +++++++-------- 3 files changed, 51 insertions(+), 37 deletions(-) diff --git a/galpy/orbit/Orbits.py b/galpy/orbit/Orbits.py index a51af94f7..8ef991bcb 100644 --- a/galpy/orbit/Orbits.py +++ b/galpy/orbit/Orbits.py @@ -1324,7 +1324,7 @@ def integrate(self,t,pot,method='symplec4_c',progressbar=True, def integrate_SOS(self,psi,pot,surface=None,t0=0., method='dop853_c',progressbar=True, - dpsi=None,numcores=_NUMCORES, + numcores=_NUMCORES, force_map=False): """ NAME: @@ -1354,8 +1354,6 @@ def integrate_SOS(self,psi,pot,surface=None,t0=0., progressbar= (True) if True, display a tqdm progress bar when integrating multiple orbits (requires tqdm to be installed!) - dpsi - if set, force the integrator to use this basic stepsize; must be an integer divisor of output stepsize (only works for the C integrators that use a fixed stepsize) (can be Quantity) - numcores - number of cores to use for Python-based multiprocessing (pure Python or using force_map=True); default = OMP_NUM_THREADS force_map= (False) if True, force use of Python-based multiprocessing (not recommended) @@ -1376,8 +1374,6 @@ def integrate_SOS(self,psi,pot,surface=None,t0=0., # Parse psi if _APY_LOADED and isinstance(psi,units.Quantity): psi= conversion.parse_angle(psi) - if _APY_LOADED and not dpsi is None and isinstance(dpsi,units.Quantity): - dpsi= conversion.parse_angle(dpsi) if _APY_LOADED and isinstance(t0,units.Quantity): t0= conversion.parse_time(t0,ro=self._ro,vo=self._vo) self._integrate_t_asQuantity= False @@ -1385,8 +1381,6 @@ def integrate_SOS(self,psi,pot,surface=None,t0=0., if pot == MWPotential: warnings.warn("Use of MWPotential as a Milky-Way-like potential is deprecated; galpy.potential.MWPotential2014, a potential fit to a large variety of dynamical constraints (see Bovy 2015), is the preferred Milky-Way-like potential in galpy", galpyWarning) - if not _check_integrate_dt(psi,dpsi): - raise ValueError('dpsi input (integrator stepsize) for Orbit.integrate must be an integer divisor of the output stepsize') # Delete attributes for interpolation and rperi etc. determination if hasattr(self,'_orbInterp'): delattr(self,'_orbInterp') if hasattr(self,'rs'): delattr(self,'rs') @@ -1411,9 +1405,9 @@ def integrate_SOS(self,psi,pot,surface=None,t0=0., progressbar=progressbar, numcores=numcores,dt=dt) else: - out, msg= integrateFullOrbit_sos(self._pot,self.vxvv,psi,t0,method, + out, msg= integrateFullOrbit_sos(self._pot,self.vxvv,self._psi,t0,method, progressbar=progressbar, - numcores=numcores,dpsi=dpsi) + numcores=numcores) else: warnings.warn("Using C implementation to integrate orbits", galpyWarningVerbose) @@ -1439,9 +1433,8 @@ def integrate_SOS(self,psi,pot,surface=None,t0=0., progressbar=progressbar, dt=dt) else: - out, msg= integrateFullOrbit_sos_c(self._pot,vxvvs,psi,t0, - method,progressbar=progressbar, - dpsi=dpsi) + out, msg= integrateFullOrbit_sos_c(self._pot,vxvvs,self._psi,t0, + method,progressbar=progressbar) if self.phasedim() == 3 \ or self.phasedim() == 5: @@ -4766,13 +4759,24 @@ def SOS(self,pot,ncross=500,surface=None,t0=0., else: raise NotImplementedError("SOS not implemented for 1D or 2D orbits") if numpy.any(numpy.fabs(init_psis) > 1e-10): - raise RuntimeError("Use of SOS requires the orbit's initial conditions to be in the plane of the surface of section") + # Integrate to the next crossing + init_psis= numpy.atleast_1d((init_psis + 2.*numpy.pi) % (2.*numpy.pi)) + psis= numpy.array([numpy.linspace(init_psi,2.*numpy.pi,11) + for init_psi in init_psis]) + self.integrate_SOS(psis,pot,surface=surface,t0=t0,method=method, + progressbar=progressbar, + numcores=numcores,force_map=force_map) + old_vxvv= self.vxvv + self.vxvv= self.orbit[:,-1] psis= numpy.arange(ncross)*2*numpy.pi self.integrate_SOS(psis,pot,surface=surface,t0=t0,method=method, progressbar=progressbar, numcores=numcores,force_map=force_map) if self.dim() == 3: - return (self.R(self.t,use_physical=False),self.vR(self.t,use_physical=False)) + out= (self.R(self.t,use_physical=False),self.vR(self.t,use_physical=False)) + if numpy.any(numpy.fabs(init_psis) > 1e-10): + self.vxvv= old_vxvv + return out def __call__(self,*args,**kwargs): """ diff --git a/galpy/orbit/integrateFullOrbit.py b/galpy/orbit/integrateFullOrbit.py index 09908a598..8ffb3cea3 100644 --- a/galpy/orbit/integrateFullOrbit.py +++ b/galpy/orbit/integrateFullOrbit.py @@ -718,9 +718,10 @@ def integrateFullOrbit_sos_c(pot,yo,psi,t0,int_method,rtol=None,atol=None, yoo[:,6]= t0[0] else: yoo[:,6]= t0 + npsi= len(psi.T) # .T to make npsi always the first dim #Set up result array - result= numpy.empty((nobj,len(psi),7)) + result= numpy.empty((nobj,npsi,7)) err= numpy.zeros(nobj,dtype=numpy.int32) #Set up progressbar @@ -740,6 +741,7 @@ def integrateFullOrbit_sos_c(pot,yo,psi,t0,int_method,rtol=None,atol=None, ctypes.c_int, ndpointer(dtype=numpy.float64,flags=ndarrayFlags), ctypes.c_int, + ctypes.c_int, ndpointer(dtype=numpy.int32,flags=ndarrayFlags), ndpointer(dtype=numpy.float64,flags=ndarrayFlags), ctypes.c_void_p, @@ -759,11 +761,12 @@ def integrateFullOrbit_sos_c(pot,yo,psi,t0,int_method,rtol=None,atol=None, result= numpy.require(result,dtype=numpy.float64,requirements=['C','W']) err= numpy.require(err,dtype=numpy.int32,requirements=['C','W']) - #Run the C code + #Run the C code) integrationFunc(ctypes.c_int(nobj), yoo, - ctypes.c_int(len(psi)), + ctypes.c_int(npsi), psi, + ctypes.c_int(len(psi.shape) > 1), ctypes.c_int(npot), pot_type, pot_args, @@ -831,7 +834,7 @@ def integrateFullOrbit_sos(pot,yo,psi,t0,int_method,rtol=None,atol=None, else: integrator= integrate.odeint extra_kwargs= {'rtol':rtol} - def integrate_for_map(vxvv): + def integrate_for_map(vxvv,psi): #go to the transformed plane: (x,vx,y,vy,A,t) init_psi= numpy.arctan2(vxvv[3],vxvv[4]) init= numpy.array([vxvv[0]*numpy.cos(vxvv[5]), @@ -856,15 +859,21 @@ def integrate_for_map(vxvv): out[:,6]= intOut[:,5] return out else: # Assume we are forcing parallel_mapping of a C integrator... - raise NotImplementedError("C-based integrators not implemented for integrateFullOrbit_sos") - def integrate_for_map(vxvv): - return integrateFullOrbit_c(pot,numpy.copy(vxvv), - t,int_method,dt=dt)[0] + def integrate_for_map(vxvv,psi): + return integrateFullOrbit_sos_c(pot,numpy.copy(vxvv),psi, + t0,int_method,dpsi=dpsi)[0] if len(yo) == 1: # Can't map a single value... - out= numpy.atleast_3d(integrate_for_map(yo[0]).T).T + out= numpy.atleast_3d(integrate_for_map(yo[0],psi).T).T + elif len(psi.shape) > 1: + out= numpy.array(parallel_map( + lambda ii: integrate_for_map(yo[ii],psi[ii]), + range(len(yo)),numcores=numcores,progressbar=progressbar + )) else: - out= numpy.array(parallel_map(integrate_for_map,yo,numcores=numcores, - progressbar=progressbar)) + out= numpy.array(parallel_map( + lambda ii: integrate_for_map(yo[ii],psi), + range(len(yo)),numcores=numcores,progressbar=progressbar + )) if nophi: phi_mask= numpy.ones(out.shape[2],dtype='bool') phi_mask[5]= False diff --git a/galpy/orbit/orbit_c_ext/integrateFullOrbit.c b/galpy/orbit/orbit_c_ext/integrateFullOrbit.c index 5b0cba239..825abff2f 100644 --- a/galpy/orbit/orbit_c_ext/integrateFullOrbit.c +++ b/galpy/orbit/orbit_c_ext/integrateFullOrbit.c @@ -730,19 +730,20 @@ EXPORT void integrateFullOrbit(int nobj, } EXPORT void integrateFullOrbit_sos( int nobj, - double *yo, - int npsi, - double *psi, + double *yo, + int npsi, + double *psi, + int indiv_psi, int npot, - int * pot_type, - double * pot_args, + int * pot_type, + double * pot_args, tfuncs_type_arr pot_tfuncs, - double dpsi, - double rtol, - double atol, - double *result, - int * err, - int odeint_type, + double dpsi, + double rtol, + double atol, + double *result, + int * err, + int odeint_type, orbint_callback_type cb){ //Set up the forces, first count int ii,jj; @@ -795,7 +796,7 @@ EXPORT void integrateFullOrbit_sos( #pragma omp parallel for schedule(dynamic,ORBITS_CHUNKSIZE) private(ii,jj) num_threads(max_threads) for (ii=0; ii < nobj; ii++) { cyl_to_sos_galpy(yo+dim*ii); - odeint_func(odeint_deriv_func,dim,yo+dim*ii,npsi,dpsi,psi, + odeint_func(odeint_deriv_func,dim,yo+dim*ii,npsi,dpsi,psi+npsi*ii*indiv_psi, npot,potentialArgs+omp_get_thread_num()*npot,rtol,atol, result+dim*npsi*ii,err+ii); for (jj=0; jj < npsi; jj++) From 7c2017e5d341a7e09ce3cddce13d15decb21a17b Mon Sep 17 00:00:00 2001 From: Jo Bovy Date: Thu, 23 Mar 2023 14:20:05 -0400 Subject: [PATCH 04/18] Add some initial tests of 3D SOS code --- galpy/orbit/Orbits.py | 17 ++++++++++++---- galpy/orbit/integrateFullOrbit.py | 2 +- tests/test_orbit.py | 34 +++++++++++++++++++++++++++++++ 3 files changed, 48 insertions(+), 5 deletions(-) diff --git a/galpy/orbit/Orbits.py b/galpy/orbit/Orbits.py index 8ef991bcb..279665b8e 100644 --- a/galpy/orbit/Orbits.py +++ b/galpy/orbit/Orbits.py @@ -4761,20 +4761,29 @@ def SOS(self,pot,ncross=500,surface=None,t0=0., if numpy.any(numpy.fabs(init_psis) > 1e-10): # Integrate to the next crossing init_psis= numpy.atleast_1d((init_psis + 2.*numpy.pi) % (2.*numpy.pi)) - psis= numpy.array([numpy.linspace(init_psi,2.*numpy.pi,11) + psis= numpy.array([numpy.linspace(0.,2.*numpy.pi-init_psi,101) for init_psi in init_psis]) self.integrate_SOS(psis,pot,surface=surface,t0=t0,method=method, progressbar=progressbar, numcores=numcores,force_map=force_map) old_vxvv= self.vxvv self.vxvv= self.orbit[:,-1] - psis= numpy.arange(ncross)*2*numpy.pi + if method == 'rk4_c' or method == 'rk6_c': + # Because these are non-adaptive, we need to make sure we + # integrate finely enough + skip= 100 + else: + skip= 1 + psis= numpy.arange(ncross*skip)*2*numpy.pi/skip self.integrate_SOS(psis,pot,surface=surface,t0=t0,method=method, progressbar=progressbar, numcores=numcores,force_map=force_map) + self.t= self.t[:,::skip] + self.orbit= self.orbit[:,::skip] if self.dim() == 3: - out= (self.R(self.t,use_physical=False),self.vR(self.t,use_physical=False)) - if numpy.any(numpy.fabs(init_psis) > 1e-10): + out= (self.R(self.t,use_physical=False), + self.vR(self.t,use_physical=False)) + if numpy.any(numpy.fabs(init_psis) > 1e-7): self.vxvv= old_vxvv return out diff --git a/galpy/orbit/integrateFullOrbit.py b/galpy/orbit/integrateFullOrbit.py index 8ffb3cea3..585c9a777 100644 --- a/galpy/orbit/integrateFullOrbit.py +++ b/galpy/orbit/integrateFullOrbit.py @@ -863,7 +863,7 @@ def integrate_for_map(vxvv,psi): return integrateFullOrbit_sos_c(pot,numpy.copy(vxvv),psi, t0,int_method,dpsi=dpsi)[0] if len(yo) == 1: # Can't map a single value... - out= numpy.atleast_3d(integrate_for_map(yo[0],psi).T).T + out= numpy.atleast_3d(integrate_for_map(yo[0],psi.flatten()).T).T elif len(psi.shape) > 1: out= numpy.array(parallel_map( lambda ii: integrate_for_map(yo[ii],psi[ii]), diff --git a/tests/test_orbit.py b/tests/test_orbit.py index 306fd34a6..9f9f52141 100644 --- a/tests/test_orbit.py +++ b/tests/test_orbit.py @@ -661,6 +661,33 @@ def test_liouville_planar(): or ('Burkert' in p and not ptp.hasC)): break return None +# Test that integrating an orbit in MWPotential2014 using integrate_SOS conserves energy +def test_integrate_SOS_3D(): + pot= potential.MWPotential2014 + o= setup_orbit_energy(pot) + psis= numpy.linspace(0.,20.*numpy.pi,1001) + for method in ['dopr54_c','dop853_c','rk4_c','rk6_c','dop853','odeint']: + o.integrate_SOS(psis,pot,method=method) + Es= o.E(o.t) + assert (numpy.std(Es)/numpy.mean(Es))**2. < 10.**-10, \ + f'Energy is not conserved by integrate_sos for method={method}' + return None + +# Test that the 3D SOS function returns points with z=0, vz > 0 +def test_SOS_3D(): + pot= potential.MWPotential2014 + o= setup_orbit_energy(pot) + psis= numpy.linspace(0.,20.*numpy.pi,1001) + for method in ['dopr54_c','dop853_c','rk4_c','rk6_c','dop853','odeint']: + o.SOS(pot,method=method,ncross=500 if '_c' in method else 20) + zs= o.z(o.t) + vzs= o.vz(o.t) + assert (numpy.fabs(zs) < 10.**-7.).all(), \ + f'z on SOS is not zero for integrate_sos for method={method}' + assert (vzs > 0.).all(), \ + f'vz on SOS is not positive for integrate_sos for method={method}' + return None + # Test that the eccentricity of circular orbits is zero def test_eccentricity(): #return None @@ -4750,6 +4777,13 @@ def test_full_plotting(): else: raise AssertionError("plot3d(d3='vy') applied to RZOrbit did not raise AttributeError") return None +def test_plotSOS(): + pot= potential.MWPotential2014 + o= setup_orbit_energy(pot) + o.plotSOS(pot) + o.plotSOS(pot,use_physical=True) + return None + def test_from_name_values(): from galpy.orbit import Orbit From 85bfb2eb2f28cb45bd5e08e7ca52d3d355c086f4 Mon Sep 17 00:00:00 2001 From: Jo Bovy Date: Thu, 23 Mar 2023 14:45:10 -0400 Subject: [PATCH 05/18] Multiple-orbits tests of 3D SOS --- tests/test_orbit.py | 1 - tests/test_orbits.py | 27 +++++++++++++++++++++++++++ 2 files changed, 27 insertions(+), 1 deletion(-) diff --git a/tests/test_orbit.py b/tests/test_orbit.py index 9f9f52141..d08f98c23 100644 --- a/tests/test_orbit.py +++ b/tests/test_orbit.py @@ -677,7 +677,6 @@ def test_integrate_SOS_3D(): def test_SOS_3D(): pot= potential.MWPotential2014 o= setup_orbit_energy(pot) - psis= numpy.linspace(0.,20.*numpy.pi,1001) for method in ['dopr54_c','dop853_c','rk4_c','rk6_c','dop853','odeint']: o.SOS(pot,method=method,ncross=500 if '_c' in method else 20) zs= o.z(o.t) diff --git a/tests/test_orbits.py b/tests/test_orbits.py index 82fd0cef7..83c00b0ce 100644 --- a/tests/test_orbits.py +++ b/tests/test_orbits.py @@ -409,6 +409,24 @@ def test_integration_dxdv_2d_rectInOut(): assert numpy.amax(numpy.fabs(orbits.getOrbit_dxdv()-numpy.array([o.getOrbit_dxdv() for o in orbits_list]))) < 1e-8, 'Integration of the phase-space volume of multiple orbits as Orbits does not agree with integrating the phase-space volume of multiple orbits' return None +# Test that the 3D SOS function returns points with z=0, vz > 0 +def test_SOS_3D(): + from galpy.orbit import Orbit + times= numpy.linspace(0.,10.,1001) + orbits_list= [Orbit([1.,0.1,1.,0.,0.1,0.]),Orbit([.9,0.3,1.,-0.3,0.4,3.]), + Orbit([1.2,-0.3,0.7,.5,-0.5,6.])] + orbits= Orbit(orbits_list) + pot= potential.MWPotential2014 + for method in ['dopr54_c','dop853_c','rk4_c','rk6_c','dop853','odeint']: + orbits.SOS(pot,method=method,ncross=500 if '_c' in method else 20) + zs= orbits.z(orbits.t) + vzs= orbits.vz(orbits.t) + assert (numpy.fabs(zs) < 10.**-6.).all(), \ + f'z on SOS is not zero for integrate_sos for method={method}' + assert (vzs > 0.).all(), \ + f'vz on SOS is not positive for integrate_sos for method={method}' + return None + # Test slicing of orbits def test_slice_singleobject(): from galpy.orbit import Orbit @@ -2112,6 +2130,15 @@ def test_plotting(): os.plot(d1='t',d2='r*R/vR') return None +def test_plotSOS(): + from galpy.orbit import Orbit + from galpy.potential import LogarithmicHaloPotential + o= Orbit([Orbit([1.,0.1,1.1,0.1,0.2,2.]),Orbit([1.,0.1,1.1,0.1,0.2,2.])]) + pot= potential.MWPotential2014 + o.plotSOS(pot) + o.plotSOS(pot,use_physical=True) + return None + def test_integrate_method_warning(): """ Test Orbits.integrate raises an error if method is invalid """ from galpy.orbit import Orbit From aa7e7920254bb9d681ea188de655b4ea496ca358 Mon Sep 17 00:00:00 2001 From: Jo Bovy Date: Thu, 23 Mar 2023 16:31:35 -0400 Subject: [PATCH 06/18] Add codecov token for hopefully more reliable uploads --- .github/workflows/build.yml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index dec98a968..d18bed9ea 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -298,3 +298,5 @@ jobs: - name: Upload coverage reports to codecov if: ${{ matrix.python-version == env.PYTHON_COVREPORTS_VERSION && matrix.os == 'ubuntu-latest' }} uses: codecov/codecov-action@v3 + with: + token: ${{ secrets.CODECOV_TOKEN }} From 9cf29503f3c029fa6c65e37d311f218a9ef910a5 Mon Sep 17 00:00:00 2001 From: Jo Bovy Date: Thu, 23 Mar 2023 16:41:45 -0400 Subject: [PATCH 07/18] Various SOS testing tweaks to provide fuller coverage, including of quantity support and warnings --- galpy/orbit/integrateFullOrbit.py | 32 +++++++++++++++---------------- tests/test_orbit.py | 20 +++++++++++++++++-- tests/test_orbits.py | 9 +++++++-- tests/test_quantity.py | 28 +++++++++++++++++++++++++++ 4 files changed, 69 insertions(+), 20 deletions(-) diff --git a/galpy/orbit/integrateFullOrbit.py b/galpy/orbit/integrateFullOrbit.py index 585c9a777..c9536cb18 100644 --- a/galpy/orbit/integrateFullOrbit.py +++ b/galpy/orbit/integrateFullOrbit.py @@ -825,16 +825,13 @@ def integrateFullOrbit_sos(pot,yo,psi,t0,int_method,rtol=None,atol=None, yo= numpy.pad(yo,((0,0),(0,1)),'constant',constant_values=0) if not '_c' in int_method: if rtol is None: rtol= 1e-8 - if int_method.lower() == 'leapfrog': - integrator= symplecticode.leapfrog - extra_kwargs= {'rtol':rtol} - elif int_method.lower() == 'dop853': + if int_method.lower() == 'dop853': integrator= dop853 extra_kwargs= {} else: integrator= integrate.odeint extra_kwargs= {'rtol':rtol} - def integrate_for_map(vxvv,psi): + def integrate_for_map(vxvv,psi,t0): #go to the transformed plane: (x,vx,y,vy,A,t) init_psi= numpy.arctan2(vxvv[3],vxvv[4]) init= numpy.array([vxvv[0]*numpy.cos(vxvv[5]), @@ -859,21 +856,24 @@ def integrate_for_map(vxvv,psi): out[:,6]= intOut[:,5] return out else: # Assume we are forcing parallel_mapping of a C integrator... - def integrate_for_map(vxvv,psi): + def integrate_for_map(vxvv,psi,t0): return integrateFullOrbit_sos_c(pot,numpy.copy(vxvv),psi, t0,int_method,dpsi=dpsi)[0] if len(yo) == 1: # Can't map a single value... - out= numpy.atleast_3d(integrate_for_map(yo[0],psi.flatten()).T).T - elif len(psi.shape) > 1: - out= numpy.array(parallel_map( - lambda ii: integrate_for_map(yo[ii],psi[ii]), - range(len(yo)),numcores=numcores,progressbar=progressbar - )) + out= numpy.atleast_3d(integrate_for_map(yo[0],psi.flatten(),t0).T).T else: - out= numpy.array(parallel_map( - lambda ii: integrate_for_map(yo[ii],psi), - range(len(yo)),numcores=numcores,progressbar=progressbar - )) + out= numpy.array( + parallel_map( + lambda ii: integrate_for_map( + yo[ii], + psi[ii] if len(psi.shape) > 1 else psi, + t0[0] if len(t0) == 1 else t0[ii] + ), + range(len(yo)), + numcores=numcores, + progressbar=progressbar + ) + ) if nophi: phi_mask= numpy.ones(out.shape[2],dtype='bool') phi_mask[5]= False diff --git a/tests/test_orbit.py b/tests/test_orbit.py index d08f98c23..7575ee779 100644 --- a/tests/test_orbit.py +++ b/tests/test_orbit.py @@ -664,7 +664,7 @@ def test_liouville_planar(): # Test that integrating an orbit in MWPotential2014 using integrate_SOS conserves energy def test_integrate_SOS_3D(): pot= potential.MWPotential2014 - o= setup_orbit_energy(pot) + o= setup_orbit_energy(pot,axi=True) psis= numpy.linspace(0.,20.*numpy.pi,1001) for method in ['dopr54_c','dop853_c','rk4_c','rk6_c','dop853','odeint']: o.integrate_SOS(psis,pot,method=method) @@ -678,7 +678,11 @@ def test_SOS_3D(): pot= potential.MWPotential2014 o= setup_orbit_energy(pot) for method in ['dopr54_c','dop853_c','rk4_c','rk6_c','dop853','odeint']: - o.SOS(pot,method=method,ncross=500 if '_c' in method else 20) + o.SOS( + pot, + method=method,ncross=500 if '_c' in method else 20, + force_map='rk' in method + ) zs= o.z(o.t) vzs= o.vz(o.t) assert (numpy.fabs(zs) < 10.**-7.).all(), \ @@ -3268,6 +3272,7 @@ def comp_orbfit(of,vxvv,ts,pot,lb=False,radec=False,ro=None,vo=None): def test_MWPotential_warning(): # Test that using MWPotential throws a warning, see #229 ts= numpy.linspace(0.,100.,1001) + psis= numpy.linspace(0.,20.*numpy.pi,1001) o= setup_orbit_energy(potential.MWPotential,axi=False) with pytest.warns(galpyWarning) as record: if PY2: reset_warning_registry('galpy') @@ -3279,6 +3284,17 @@ def test_MWPotential_warning(): # check that the message matches raisedWarning+= (str(rec.message.args[0]) == "Use of MWPotential as a Milky-Way-like potential is deprecated; galpy.potential.MWPotential2014, a potential fit to a large variety of dynamical constraints (see Bovy 2015), is the preferred Milky-Way-like potential in galpy") assert raisedWarning, "Orbit integration with MWPotential should have thrown a warning, but didn't" + # Also test for SOS integration + with pytest.warns(galpyWarning) as record: + if PY2: reset_warning_registry('galpy') + warnings.simplefilter("always",galpyWarning) + o.integrate_SOS(psis,potential.MWPotential) + # Should raise warning bc of MWPotential, might raise others + raisedWarning= False + for rec in record: + # check that the message matches + raisedWarning+= (str(rec.message.args[0]) == "Use of MWPotential as a Milky-Way-like potential is deprecated; galpy.potential.MWPotential2014, a potential fit to a large variety of dynamical constraints (see Bovy 2015), is the preferred Milky-Way-like potential in galpy") + assert raisedWarning, "Orbit integration with MWPotential should have thrown a warning, but didn't" return None # Test the new Orbit.time function diff --git a/tests/test_orbits.py b/tests/test_orbits.py index 83c00b0ce..b15020b9a 100644 --- a/tests/test_orbits.py +++ b/tests/test_orbits.py @@ -418,7 +418,12 @@ def test_SOS_3D(): orbits= Orbit(orbits_list) pot= potential.MWPotential2014 for method in ['dopr54_c','dop853_c','rk4_c','rk6_c','dop853','odeint']: - orbits.SOS(pot,method=method,ncross=500 if '_c' in method else 20) + orbits.SOS( + pot, + method=method,ncross=500 if '_c' in method else 20, + force_map='rk' in method, + t0=numpy.arange(len(orbits)), + ) zs= orbits.z(orbits.t) vzs= orbits.vz(orbits.t) assert (numpy.fabs(zs) < 10.**-6.).all(), \ @@ -2136,7 +2141,7 @@ def test_plotSOS(): o= Orbit([Orbit([1.,0.1,1.1,0.1,0.2,2.]),Orbit([1.,0.1,1.1,0.1,0.2,2.])]) pot= potential.MWPotential2014 o.plotSOS(pot) - o.plotSOS(pot,use_physical=True) + o.plotSOS(pot,use_physical=True,ro=8.,vo=220.) return None def test_integrate_method_warning(): diff --git a/tests/test_quantity.py b/tests/test_quantity.py index f6a40b4e6..4a14626c7 100644 --- a/tests/test_quantity.py +++ b/tests/test_quantity.py @@ -912,6 +912,34 @@ def test_integrate_dxdv_timeAsQuantity_Myr(): assert numpy.all(numpy.fabs(dx-dxc) < 10.**-8.), 'Orbit integrated_dxdv with times specified as Quantity does not agree with Orbit integrated_dxdv with time specified as array' return None +def test_integrate_SOS_psiQuantity(): + import copy + + from galpy.orbit import Orbit + from galpy.potential import MWPotential + from galpy.util import conversion + ro, vo= 8., 200. + o= Orbit([10.*units.kpc,-20.*units.km/units.s,210.*units.km/units.s, + 500.*units.pc,-12.*units.km/units.s,45.*units.deg], + ro=ro,vo=vo) + oc= o() + psis_nounits= numpy.linspace(0.,400.,1001) + psis= units.Quantity(copy.copy(psis_nounits),unit=units.deg) + psis_nounits/= 180./numpy.pi + t0_nounits= 1. + t0= units.Quantity(copy.copy(t0_nounits),unit=units.Gyr) + t0_nounits/= conversion.time_in_Gyr(vo,ro) + # Integrate both with Quantity time and with unitless time + o.integrate_SOS(psis,MWPotential,t0=t0) + oc.integrate_SOS(psis_nounits,MWPotential,t0=t0_nounits) + assert numpy.all(numpy.fabs(o.x(o.t)-oc.x(oc.t)).value < 10.**-8.), 'Orbit SOS integrated with psis specified as Quantity does not agree with Orbit integrated with time specified as array' + assert numpy.all(numpy.fabs(o.y(o.t)-oc.y(oc.t)).value < 10.**-8.), 'Orbit SOS integrated with psis specified as Quantity does not agree with Orbit integrated with time specified as array' + assert numpy.all(numpy.fabs(o.z(o.t)-oc.z(oc.t)).value < 10.**-8.), 'Orbit SOS ntegrated with psis specified as Quantity does not agree with Orbit integrated with time specified as array' + assert numpy.all(numpy.fabs(o.vx(o.t)-oc.vx(oc.t)).value < 10.**-8.), 'Orbit SOS integrated with psis specified as Quantity does not agree with Orbit integrated with time specified as array' + assert numpy.all(numpy.fabs(o.vy(o.t)-oc.vy(oc.t)).value < 10.**-8.), 'Orbit SOS integrated with psis specified as Quantity does not agree with Orbit integrated with time specified as array' + assert numpy.all(numpy.fabs(o.vz(o.t)-oc.vz(oc.t)).value < 10.**-8.), 'Orbit SOS integrated with psis specified as Quantity does not agree with Orbit integrated with time specified as array' + return None + def test_orbit_inconsistentPotentialUnits_error(): from galpy.orbit import Orbit from galpy.potential import IsochronePotential From 13d0d9ddbecb43e4eff005149d2d914fc86a493b Mon Sep 17 00:00:00 2001 From: Jo Bovy Date: Mon, 27 Mar 2023 21:42:24 -0400 Subject: [PATCH 08/18] Add 2D SOS integration, in Python and C --- galpy/orbit/Orbits.py | 61 ++-- galpy/orbit/integrateFullOrbit.py | 2 +- galpy/orbit/integratePlanarOrbit.py | 284 ++++++++++++++++++ galpy/orbit/orbit_c_ext/integrateFullOrbit.c | 3 +- .../orbit/orbit_c_ext/integratePlanarOrbit.c | 159 ++++++++++ galpy/util/bovy_coords.c | 67 +++++ galpy/util/bovy_coords.h | 2 + 7 files changed, 555 insertions(+), 23 deletions(-) diff --git a/galpy/orbit/Orbits.py b/galpy/orbit/Orbits.py index 279665b8e..a275a2891 100644 --- a/galpy/orbit/Orbits.py +++ b/galpy/orbit/Orbits.py @@ -44,7 +44,9 @@ integrateLinearOrbit_c) from .integratePlanarOrbit import (integratePlanarOrbit, integratePlanarOrbit_c, - integratePlanarOrbit_dxdv) + integratePlanarOrbit_dxdv, + integratePlanarOrbit_sos, + integratePlanarOrbit_sos_c) ext_loaded= _ext_loaded if _APY_LOADED: @@ -1400,10 +1402,9 @@ def integrate_SOS(self,psi,pot,surface=None,t0=0., progressbar=progressbar, numcores=numcores,dt=dt) elif self.dim() == 2: - raise NotImplementedError("SOS integration not implemented for 2D orbits") - out, msg= integratePlanarOrbit(self._pot,self.vxvv,t,method, - progressbar=progressbar, - numcores=numcores,dt=dt) + out, msg= integratePlanarOrbit_sos(self._pot,self.vxvv,self._psi,t0,method, + surface=surface,progressbar=progressbar, + numcores=numcores) else: out, msg= integrateFullOrbit_sos(self._pot,self.vxvv,self._psi,t0,method, progressbar=progressbar, @@ -1427,11 +1428,9 @@ def integrate_SOS(self,psi,pot,surface=None,t0=0., else: vxvvs= numpy.copy(self.vxvv) if self.dim() == 2: - raise NotImplementedError("SOS integration not implemented for 2D orbits") - out, msg= integratePlanarOrbit_c(self._pot,vxvvs, - t,method, - progressbar=progressbar, - dt=dt) + out, msg= integratePlanarOrbit_sos_c(self._pot,vxvvs,self._psi,t0, + method,surface=surface, + progressbar=progressbar) else: out, msg= integrateFullOrbit_sos_c(self._pot,vxvvs,self._psi,t0, method,progressbar=progressbar) @@ -1439,7 +1438,7 @@ def integrate_SOS(self,psi,pot,surface=None,t0=0., if self.phasedim() == 3 \ or self.phasedim() == 5: phi_mask= numpy.ones(out.shape[2],dtype='bool') - phi_mask[5]= False + phi_mask[3+2*(self.phasedim()==5)]= False out= out[:,:,phi_mask] # Store orbit internally self.orbit= out[:,:,:-1] @@ -4708,7 +4707,8 @@ def SkyCoord(self,*args,**kwargs): @physical_conversion_tuple(['position','velocity']) def SOS(self,pot,ncross=500,surface=None,t0=0., - method='dop853_c',progressbar=True, + method='dop853_c',skip=100, + progressbar=True, numcores=_NUMCORES, force_map=False,**kwargs): """ @@ -4739,6 +4739,8 @@ def SOS(self,pot,ncross=500,surface=None,t0=0., 'dop853' for a 8-5-3 Dormand-Prince integrator in Python 'dop853_c' for a 8-5-3 Dormand-Prince integrator in C + skip= (100) for non-adaptive integrators, the number of basic steps to take between crossings (these are further refined in the code, but only up to a maximum refinement, so you can use skip to get finer integration in cases where more accuracy is needed) + progressbar= (True) if True, display a tqdm progress bar when integrating multiple orbits (requires tqdm to be installed!) numcores - number of cores to use for Python-based multiprocessing (pure Python or using force_map=True); default = OMP_NUM_THREADS @@ -4756,8 +4758,13 @@ def SOS(self,pot,ncross=500,surface=None,t0=0., """ if self.dim() == 3: init_psis= numpy.arctan2(self.z(use_physical=False),self.vz(use_physical=False)) + elif self.phasedim() == 4: + if not surface is None and surface.lower() == 'y': + init_psis= numpy.arctan2(self.y(use_physical=False),self.vy(use_physical=False)) + else: + init_psis= numpy.arctan2(self.x(use_physical=False),self.vx(use_physical=False)) else: - raise NotImplementedError("SOS not implemented for 1D or 2D orbits") + raise NotImplementedError("SOS not implemented for 1D orbits or 2D orbits without phi") if numpy.any(numpy.fabs(init_psis) > 1e-10): # Integrate to the next crossing init_psis= numpy.atleast_1d((init_psis + 2.*numpy.pi) % (2.*numpy.pi)) @@ -4771,18 +4778,24 @@ def SOS(self,pot,ncross=500,surface=None,t0=0., if method == 'rk4_c' or method == 'rk6_c': # Because these are non-adaptive, we need to make sure we # integrate finely enough - skip= 100 + iskip= skip else: - skip= 1 - psis= numpy.arange(ncross*skip)*2*numpy.pi/skip + iskip= 1 + psis= numpy.arange(ncross*iskip)*2*numpy.pi/iskip self.integrate_SOS(psis,pot,surface=surface,t0=t0,method=method, progressbar=progressbar, numcores=numcores,force_map=force_map) - self.t= self.t[:,::skip] - self.orbit= self.orbit[:,::skip] + self.t= self.t[:,::iskip] + self.orbit= self.orbit[:,::iskip] if self.dim() == 3: out= (self.R(self.t,use_physical=False), self.vR(self.t,use_physical=False)) + elif not surface is None and surface.lower() == 'y': + out= (self.x(self.t,use_physical=False), + self.vx(self.t,use_physical=False)) + else: + out= (self.y(self.t,use_physical=False), + self.vy(self.t,use_physical=False)) if numpy.any(numpy.fabs(init_psis) > 1e-7): self.vxvv= old_vxvv return out @@ -5293,7 +5306,7 @@ def plot3d(self,*args,**kwargs): return [line3d] def plotSOS(self,pot,*args,ncross=500,surface=None,t0=0., - method='dop853_c',progressbar=True, + method='dop853_c',skip=100,progressbar=True, **kwargs): """ NAME: @@ -5323,6 +5336,8 @@ def plotSOS(self,pot,*args,ncross=500,surface=None,t0=0., 'dop853' for a 8-5-3 Dormand-Prince integrator in Python 'dop853_c' for a 8-5-3 Dormand-Prince integrator in C + skip= (100) for non-adaptive integrators, the number of basic steps to take between crossings (these are further refined in the code, but only up to a maximum refinement, so you can use skip to get finer integration in cases where more accuracy is needed) + progressbar= (True) if True, display a tqdm progress bar when integrating multiple orbits (requires tqdm to be installed!) for more control of the integrator, use the SOS method directly and plot its results @@ -5349,8 +5364,14 @@ def plotSOS(self,pot,*args,ncross=500,surface=None,t0=0., if self.dim() == 3: d1= 'R' d2= 'vR' + elif not surface is None and surface.lower() == 'y': + d1= 'x' + d2= 'vx' + else: + d1= 'y' + d2= 'vy' x,y= self.SOS(pot,ncross=ncross,surface=surface, - t0=t0,method=method, + t0=t0,method=method,skip=skip, progressbar=progressbar,**kwargs) x= numpy.atleast_2d(x) y= numpy.atleast_2d(y) diff --git a/galpy/orbit/integrateFullOrbit.py b/galpy/orbit/integrateFullOrbit.py index c9536cb18..d02d4e39e 100644 --- a/galpy/orbit/integrateFullOrbit.py +++ b/galpy/orbit/integrateFullOrbit.py @@ -688,7 +688,7 @@ def integrateFullOrbit_sos_c(pot,yo,psi,t0,int_method,rtol=None,atol=None, yo - initial condition [q,p], shape [N,5] or [N,6] psi - set of increment angles at which one wants the result [increments wrt initial angle] t0 - initial time - int_method= 'leapfrog', 'odeint', or 'dop853' + int_method= 'rk4_c', 'rk6_c', 'dopr54_c', or 'dop853_c' rtol, atol= tolerances (not always used...) progressbar= (True) if True, display a tqdm progress bar when integrating multiple orbits (requires tqdm to be installed!) dpsi= (None) force integrator to use this stepsize (default is to automatically determine one; only for C-based integrators) diff --git a/galpy/orbit/integratePlanarOrbit.py b/galpy/orbit/integratePlanarOrbit.py index 33adb28df..e9c42dd85 100644 --- a/galpy/orbit/integratePlanarOrbit.py +++ b/galpy/orbit/integratePlanarOrbit.py @@ -832,6 +832,242 @@ def integrate_for_map(vxvv): out[...,6]= dvT return out, numpy.zeros(len(yo)) +def integratePlanarOrbit_sos_c(pot,yo,psi,t0,int_method,surface='x', + rtol=None,atol=None, + progressbar=True,dpsi=None): + """ + NAME: + integratePlanarOrbit_sos_c + PURPOSE: + Integrate an ode for a PlanarOrbit for integrate_sos in C + INPUT: + pot - Potential or list of such instances + yo - initial condition [q,p], shape [N,5] or [N,6] + psi - set of increment angles at which one wants the result [increments wrt initial angle] + t0 - initial time + int_method= 'rk4_c', 'rk6_c', 'dopr54_c', or 'dop853_c' + surface= ('x') surface to use ('x' for finding x=0, vx>0; 'y' for finding y=0, vy>0) + rtol, atol= tolerances (not always used...) + progressbar= (True) if True, display a tqdm progress bar when integrating multiple orbits (requires tqdm to be installed!) + dpsi= (None) force integrator to use this stepsize (default is to automatically determine one; only for C-based integrators) + OUTPUT: + (y,err) + y : array, shape (N,len(psi),5) where the last of the last dimension is the time + Array containing the value of y for each desired angle in psi, \ + with the initial value y0 in the first row. + err: error message, always zero for now + HISTORY: + 2023-03-17 - Written based on integrateFullOrbit_sos_c - Bovy (UofT) + """ + if len(yo.shape) == 1: single_obj= True + else: single_obj= False + yo= numpy.atleast_2d(yo) + nobj= len(yo) + rtol, atol= _parse_tol(rtol,atol) + npot, pot_type, pot_args, pot_tfuncs= _parse_pot(pot) + pot_tfuncs= _prep_tfuncs(pot_tfuncs) + int_method_c= _parse_integrator(int_method) + if dpsi is None: + dpsi= -9999.99 + t0= numpy.atleast_1d(t0) + yoo= numpy.empty((nobj,5)) + yoo[:,:4]= yo[:,:4] + if len(t0) == 1: + yoo[:,4]= t0[0] + else: + yoo[:,4]= t0 + npsi= len(psi.T) # .T to make npsi always the first dim + + #Set up result array + result= numpy.empty((nobj,npsi,5)) + err= numpy.zeros(nobj,dtype=numpy.int32) + + #Set up progressbar + progressbar*= _TQDM_LOADED + if nobj > 1 and progressbar: + pbar= tqdm.tqdm(total=nobj,leave=False) + pbar_func_ctype= ctypes.CFUNCTYPE(None) + pbar_c= pbar_func_ctype(pbar.update) + else: # pragma: no cover + pbar_c= None + + #Set up the C code + ndarrayFlags= ('C_CONTIGUOUS','WRITEABLE') + integrationFunc= _lib.integratePlanarOrbit_sos + integrationFunc.argtypes= [ctypes.c_int, + ndpointer(dtype=numpy.float64,flags=ndarrayFlags), + ctypes.c_int, + ndpointer(dtype=numpy.float64,flags=ndarrayFlags), + ctypes.c_int, + ctypes.c_int, + ctypes.c_int, + ndpointer(dtype=numpy.int32,flags=ndarrayFlags), + ndpointer(dtype=numpy.float64,flags=ndarrayFlags), + ctypes.c_void_p, + ctypes.c_double, + ctypes.c_double, + ctypes.c_double, + ndpointer(dtype=numpy.float64,flags=ndarrayFlags), + ndpointer(dtype=numpy.int32,flags=ndarrayFlags), + ctypes.c_int, + ctypes.c_void_p] + + #Array requirements, first store old order + f_cont= [yoo.flags['F_CONTIGUOUS'], + psi.flags['F_CONTIGUOUS']] + yoo= numpy.require(yoo,dtype=numpy.float64,requirements=['C','W']) + psi= numpy.require(psi,dtype=numpy.float64,requirements=['C','W']) + result= numpy.require(result,dtype=numpy.float64,requirements=['C','W']) + err= numpy.require(err,dtype=numpy.int32,requirements=['C','W']) + + #Run the C code) + integrationFunc(ctypes.c_int(nobj), + yoo, + ctypes.c_int(npsi), + psi, + ctypes.c_int(len(psi.shape) > 1), + ctypes.c_int(1 if surface == 'y' else 0), + ctypes.c_int(npot), + pot_type, + pot_args, + pot_tfuncs, + ctypes.c_double(dpsi), + ctypes.c_double(rtol), + ctypes.c_double(atol), + result, + err, + ctypes.c_int(int_method_c), + pbar_c) + + if nobj > 1 and progressbar: + pbar.close() + + if numpy.any(err == -10): #pragma: no cover + raise KeyboardInterrupt("Orbit integration interrupted by CTRL-C (SIGINT)") + + #Reset input arrays + if f_cont[0]: yoo= numpy.asfortranarray(yoo) + if f_cont[1]: psi= numpy.asfortranarray(psi) + + if single_obj: return (result[0],err[0]) + else: return (result,err) + +def integratePlanarOrbit_sos(pot,yo,psi,t0,int_method,surface='x', + rtol=None,atol=None,numcores=1,progressbar=True, + dpsi=None): + """ + NAME: + integratePlanarOrbit_sos + PURPOSE: + Integrate an ode for a PlanarOrbit for integrate_sos + INPUT: + pot - Potential or list of such instances + yo - initial condition [q,p], shape [N,5] or [N,6] + psi - set of increment angles at which one wants the result [increments wrt initial angle] + t0 - initial time + surface= ('x') surface to use ('x' for finding x=0, vx>0; 'y' for finding y=0, vy>0) + int_method= 'leapfrog', 'odeint', or 'dop853' + rtol, atol= tolerances (not always used...) + numcores= (1) number of cores to use for multi-processing + progressbar= (True) if True, display a tqdm progress bar when integrating multiple orbits (requires tqdm to be installed!) + dpsi= (None) force integrator to use this stepsize (default is to automatically determine one; only for C-based integrators) + OUTPUT: + (y,err) + y : array, shape (N,len(psi),4/5) where the last of the last dimension is the time + Array containing the value of y for each desired angle in psi, \ + with the initial value y0 in the first row. + err: error message, always zero for now + HISTORY: + 2023-03-24 - Written based on integrateFullOrbi_sos - Bovy (UofT) + """ + if surface is None: + surface= 'x' + nophi= False + if len(yo[0]) == 3: + nophi= True + #We hack this by putting in a dummy phi=0 + yo= numpy.pad(yo,((0,0),(0,1)),'constant',constant_values=0) + if not '_c' in int_method: + if rtol is None: rtol= 1e-8 + if int_method.lower() == 'dop853': + integrator= dop853 + extra_kwargs= {} + else: + integrator= integrate.odeint + extra_kwargs= {'rtol':rtol} + def integrate_for_map(vxvv,psi,t0): + #go to the transformed plane: (A,t,y,vy) or (x,vx,A,t) + if surface.lower() == 'x': + x= vxvv[0]*numpy.cos(vxvv[3]) + vx= vxvv[1]*numpy.cos(vxvv[3])-vxvv[2]*numpy.sin(vxvv[3]) + init_psi= numpy.arctan2(x,vx) + this_vxvv= numpy.array([ + vxvv[0]*numpy.sin(vxvv[3]), + vxvv[2]*numpy.cos(vxvv[3])+vxvv[1]*numpy.sin(vxvv[3]), + numpy.sqrt(x**2.+vx**2.), + t0, + ]) + #integrate + intOut= integrator(_planarSOSEOMx,this_vxvv, + t=psi+init_psi,args=(pot,), + **extra_kwargs) + #go back to the cylindrical frame + x= intOut[:,2]*numpy.sin(psi+init_psi) + vx= intOut[:,2]*numpy.cos(psi+init_psi) + y= intOut[:,0] + vy= intOut[:,1] + else: + y= vxvv[0]*numpy.sin(vxvv[3]) + vy= vxvv[2]*numpy.cos(vxvv[3])+vxvv[1]*numpy.sin(vxvv[3]) + init_psi= numpy.arctan2(y,vy) + this_vxvv= numpy.array([ + vxvv[0]*numpy.cos(vxvv[3]), + vxvv[1]*numpy.cos(vxvv[3])-vxvv[2]*numpy.sin(vxvv[3]), + numpy.sqrt(y**2.+vy**2.), + t0 + ]) + #integrate + intOut= integrator(_planarSOSEOMy,this_vxvv, + t=psi+init_psi,args=(pot,), + **extra_kwargs) + #go back to the cylindrical frame + x= intOut[:,0] + vx= intOut[:,1] + y= intOut[:,2]*numpy.sin(psi+init_psi) + vy= intOut[:,2]*numpy.cos(psi+init_psi) + out= numpy.zeros((len(psi),5)) + out[:,0]= numpy.sqrt(x**2.+y**2.) + out[:,3]= numpy.arctan2(y,x) + out[:,1]= vx*numpy.cos(out[:,3])+vy*numpy.sin(out[:,3]) + out[:,2]= vy*numpy.cos(out[:,3])-vx*numpy.sin(out[:,3]) + out[:,4]= intOut[:,3] + return out + else: # Assume we are forcing parallel_mapping of a C integrator... + def integrate_for_map(vxvv,psi,t0): + return integratePlanarOrbit_sos_c(pot,numpy.copy(vxvv),psi, + t0,int_method,surface=surface, + dpsi=dpsi)[0] + if len(yo) == 1: # Can't map a single value... + out= numpy.atleast_3d(integrate_for_map(yo[0],psi.flatten(),t0).T).T + else: + out= numpy.array( + parallel_map( + lambda ii: integrate_for_map( + yo[ii], + psi[ii] if len(psi.shape) > 1 else psi, + t0[0] if len(t0) == 1 else t0[ii] + ), + range(len(yo)), + numcores=numcores, + progressbar=progressbar + ) + ) + if nophi: + phi_mask= numpy.ones(out.shape[2],dtype='bool') + phi_mask[3]= False + out= out[:,:,phi_mask] + return out, numpy.zeros(len(yo)) + def _planarREOM(y,t,pot,l2): """ NAME: @@ -932,6 +1168,54 @@ def _planarEOM_dxdv(x,t,pot): dFxdx*x[4]+dFxdy*x[5], dFydx*x[4]+dFydy*x[5]]) +def _planarSOSEOMx(y,psi,pot): + """ + NAME: + _planarSOSEOMx + PURPOSE: + implements the EOM, i.e., the right-hand side of the differential + equation, for integrating a general planar Orbit in the SOS style + INPUT: + y - current phase-space position + psi - current angle + pot - (list of) Potential instance(s) + OUTPUT: + dy/dt + HISTORY: + 2023-03-24 - Written - Bovy (UofT) + """ + # y = (y,vy,A,t) + # Calculate x + sp, cp= numpy.sin(psi), numpy.cos(psi) + gxyz= _planarRectForce([y[2]*sp,y[0]],pot,t=y[3]) + psidot= cp**2.-sp/y[2]*gxyz[0] + Adot= y[2]*cp*sp+gxyz[0]*cp + return numpy.array([y[1],gxyz[1],Adot,1.])/psidot + +def _planarSOSEOMy(y,psi,pot): + """ + NAME: + _planarSOSEOMy + PURPOSE: + implements the EOM, i.e., the right-hand side of the differential + equation, for integrating a general planar Orbit in the SOS style + INPUT: + y - current phase-space position + psi - current angle + pot - (list of) Potential instance(s) + OUTPUT: + dy/dt + HISTORY: + 2023-03-24 - Written - Bovy (UofT) + """ + # y = (x,vx,A,t) + # Calculate y + sp, cp= numpy.sin(psi), numpy.cos(psi) + gxyz= _planarRectForce([y[0],y[2]*sp],pot,t=y[3]) + psidot= cp**2.-sp/y[2]*gxyz[1] + Adot= y[2]*cp*sp+gxyz[1]*cp + return numpy.array([y[1],gxyz[0],Adot,1.])/psidot + def _planarRectForce(x,pot,t=0.): """ NAME: diff --git a/galpy/orbit/orbit_c_ext/integrateFullOrbit.c b/galpy/orbit/orbit_c_ext/integrateFullOrbit.c index 825abff2f..0dc3ce19a 100644 --- a/galpy/orbit/orbit_c_ext/integrateFullOrbit.c +++ b/galpy/orbit/orbit_c_ext/integrateFullOrbit.c @@ -932,8 +932,7 @@ void evalRectDeriv(double t, double *q, double *a, void evalSOSDeriv(double psi, double *q, double *a, int nargs, struct potentialArg * potentialArgs){ - // q= (x,y,vx,vy,A,t,psi); to save operations, we use the first three - // components of q as (x,y,z) and then re-use a first for the + // q= (x,y,vx,vy,A,t,psi); to save operations, we re-use a first for the // rectForce then for the actual RHS // Note also that we keep track of psi in q+6, not in psi! This is // such that we can avoid having to convert psi to psi+psi0 diff --git a/galpy/orbit/orbit_c_ext/integratePlanarOrbit.c b/galpy/orbit/orbit_c_ext/integratePlanarOrbit.c index 4f185dcb3..457270f7d 100644 --- a/galpy/orbit/orbit_c_ext/integratePlanarOrbit.c +++ b/galpy/orbit/orbit_c_ext/integratePlanarOrbit.c @@ -34,6 +34,10 @@ void evalPlanarRectForce(double, double *, double *, int, struct potentialArg *); void evalPlanarRectDeriv(double, double *, double *, int, struct potentialArg *); +void evalPlanarSOSDerivx(double, double *, double *, + int, struct potentialArg *); +void evalPlanarSOSDerivy(double, double *, double *, + int, struct potentialArg *); void evalPlanarRectDeriv_dxdv(double, double *, double *, int, struct potentialArg *); void initPlanarMovingObjectSplines(struct potentialArg *, double ** pot_args); @@ -628,7 +632,98 @@ EXPORT void integratePlanarOrbit(int nobj, free(potentialArgs); //Done! } +EXPORT void integratePlanarOrbit_sos( + int nobj, + double *yo, + int npsi, + double *psi, + int indiv_psi, + int surface, + int npot, + int * pot_type, + double * pot_args, + tfuncs_type_arr pot_tfuncs, + double dpsi, + double rtol, + double atol, + double *result, + int * err, + int odeint_type, + orbint_callback_type cb){ + //Set up the forces, first count + int ii,jj; + int dim; + int max_threads; + int * thread_pot_type; + double * thread_pot_args; + tfuncs_type_arr thread_pot_tfuncs; + max_threads= ( nobj < omp_get_max_threads() ) ? nobj : omp_get_max_threads(); + // Because potentialArgs may cache, safest to have one / thread + struct potentialArg * potentialArgs= (struct potentialArg *) malloc ( max_threads * npot * sizeof (struct potentialArg) ); +#pragma omp parallel for schedule(static,1) private(ii,thread_pot_type,thread_pot_args,thread_pot_tfuncs) num_threads(max_threads) + for (ii=0; ii < max_threads; ii++) { + thread_pot_type= pot_type; // need to make thread-private pointers, bc + thread_pot_args= pot_args; // these pointers are changed in parse_... + thread_pot_tfuncs= pot_tfuncs; // ... + parse_leapFuncArgs(npot,potentialArgs+ii*npot, + &thread_pot_type,&thread_pot_args,&thread_pot_tfuncs); + } + //Integrate + void (*odeint_func)(void (*func)(double, double *, double *, + int, struct potentialArg *), + int, + double *, + int, double, double *, + int, struct potentialArg *, + double, double, + double *,int *); + void (*odeint_deriv_func)(double, double *, double *, + int,struct potentialArg *); + dim= 5; + switch ( odeint_type ) { + // case 0: = leapfrog = not supported symplectic method + case 1: //RK4 + odeint_func= &bovy_rk4; + break; + case 2: //RK6 + odeint_func= &bovy_rk6; + break; + // case 3: = symplec4 = not supported symplectic method + // case 4: = symplec6 = not supported symplectic method + case 5: //DOPR54 + odeint_func= &bovy_dopr54; + break; + case 6: //DOP853 + odeint_func= &dop853; + break; + } + switch ( surface ) { + case 0: // x=0 + odeint_deriv_func= &evalPlanarSOSDerivx; + break; + case 1: // y=0 + odeint_deriv_func= &evalPlanarSOSDerivy; + break; + } +#pragma omp parallel for schedule(dynamic,ORBITS_CHUNKSIZE) private(ii,jj) num_threads(max_threads) + for (ii=0; ii < nobj; ii++) { + polar_to_sos_galpy(yo+dim*ii,surface); + odeint_func(odeint_deriv_func,dim,yo+dim*ii,npsi,dpsi,psi+npsi*ii*indiv_psi, + npot,potentialArgs+omp_get_thread_num()*npot,rtol,atol, + result+dim*npsi*ii,err+ii); + for (jj=0; jj < npsi; jj++) + sos_to_polar_galpy(result+dim*jj+dim*npsi*ii,surface); + if ( cb ) // Callback if not void + cb(); + } + //Free allocated memory +#pragma omp parallel for schedule(static,1) private(ii) num_threads(max_threads) + for (ii=0; ii < max_threads; ii++) + free_potentialArgs(npot,potentialArgs+ii*npot); + free(potentialArgs); + //Done! +} EXPORT void integratePlanarOrbit_dxdv(double *yo, int nt, double *t, @@ -727,6 +822,70 @@ void evalPlanarRectDeriv(double t, double *q, double *a, *a= sinphi*Rforce+1./R*cosphi*phitorque; } +void evalPlanarSOSDerivx(double psi, double *q, double *a, + int nargs, struct potentialArg * potentialArgs){ + // q= (y,vy,A,t,psi); to save operations, we re-use a first for the + // rectForce then for the actual RHS + // Note also that we keep track of psi in q+4, not in psi! This is + // such that we can avoid having to convert psi to psi+psi0 + // q+4 starts as psi0 and then just increments as psi (exactly) + double sinpsi,cospsi,psidot,x,y,R,phi,sinphi,cosphi,Rforce,phitorque; + sinpsi= sin( *(q+4) ); + cospsi= cos( *(q+4) ); + // Calculate forces, put them in a+2, a+3 + //q is rectangular so calculate R and phi + x= *(q+2) * sinpsi;; + y= *(q ); + R= sqrt(x*x+y*y); + phi= atan2( y ,x ); + sinphi= y/R; + cosphi= x/R; + //Calculate the forces + Rforce= calcPlanarRforce(R,phi,*(q+3),nargs,potentialArgs); + phitorque= calcPlanarphitorque(R,phi,*(q+3),nargs,potentialArgs); + *(a+2)= cosphi*Rforce-1./R*sinphi*phitorque; + *(a+3)= sinphi*Rforce+1./R*cosphi*phitorque; + // Now calculate the RHS of the ODE + psidot= cospsi * cospsi - sinpsi * *(a+2) / ( *(q+2) ); + *(a )= *(q+1) / psidot; + *(a+1)= *(a+3) / psidot; + *(a+2)= cospsi * ( *(q+2) * sinpsi + *(a+2) ) / psidot; + *(a+3)= 1./psidot; + *(a+4)= 1.; // dpsi / dpsi to keep track of psi +} + +void evalPlanarSOSDerivy(double psi, double *q, double *a, + int nargs, struct potentialArg * potentialArgs){ + // q= (x,vx,A,t,psi); to save operations, we re-use a first for the + // rectForce then for the actual RHS + // Note also that we keep track of psi in q+4, not in psi! This is + // such that we can avoid having to convert psi to psi+psi0 + // q+4 starts as psi0 and then just increments as psi (exactly) + double sinpsi,cospsi,psidot,x,y,R,phi,sinphi,cosphi,Rforce,phitorque; + sinpsi= sin( *(q+4) ); + cospsi= cos( *(q+4) ); + // Calculate forces, put them in a+2, a+3 + //q is rectangular so calculate R and phi + x= *(q ); + y= *(q+2) * sinpsi; + R= sqrt(x*x+y*y); + phi= atan2( y ,x ); + sinphi= y/R; + cosphi= x/R; + //Calculate the forces + Rforce= calcPlanarRforce(R,phi,*(q+3),nargs,potentialArgs); + phitorque= calcPlanarphitorque(R,phi,*(q+3),nargs,potentialArgs); + *(a+2)= cosphi*Rforce-1./R*sinphi*phitorque; + *(a+3)= sinphi*Rforce+1./R*cosphi*phitorque; + // Now calculate the RHS of the ODE + psidot= cospsi * cospsi - sinpsi * *(a+3) / ( *(q+2) ); + *(a )= *(q+1) / psidot; + *(a+1)= *(a+2) / psidot; + *(a+2)= cospsi * ( *(q+2) * sinpsi + *(a+3) ) / psidot; + *(a+3)= 1./psidot; + *(a+4)= 1.; // dpsi / dpsi to keep track of psi +} + void evalPlanarRectDeriv_dxdv(double t, double *q, double *a, int nargs, struct potentialArg * potentialArgs){ double sinphi, cosphi, x, y, phi,R,Rforce,phitorque; diff --git a/galpy/util/bovy_coords.c b/galpy/util/bovy_coords.c index d6d718ab5..0f559239d 100644 --- a/galpy/util/bovy_coords.c +++ b/galpy/util/bovy_coords.c @@ -187,3 +187,70 @@ void sos_to_cyl_galpy(double *vxvv){ *(vxvv+6)= *(vxvv+5); *(vxvv+5)= phi; } +/* +NAME: polar_to_sos_galpy +PURPOSE: convert (R,vR,vT,phi,t) to (x,vx,A,t,psi) coordinates for SOS integration [or (y,vy,A,t) if surface == 1] +INPUT: + double * vxvv - (R,vR,vT,z,vz,phi,t) + int surface - if 1, convert to (x,vx,A,t) instead of (y,vy,A,t) +OUTPUT: + performed in-place +HISTORY: 2023-03-24 - Written - Bovy (UofT) + */ +void polar_to_sos_galpy(double *vxvv,int surface){ + double R,cp,sp,vR,vT,x,y,vx,vy; + R = *vxvv; + cp = cos ( *(vxvv+3) ); + sp = sin ( *(vxvv+3) ); + vR = *(vxvv+1); + vT = *(vxvv+2); + if ( surface == 1 ) { // surface: y=0, vy > 0 + *(vxvv )= R * cp; + *(vxvv+1)= vR * cp - vT * sp; + y= R * sp; + vy= vR * sp + vT * cp; + *(vxvv+2)= sqrt ( y * y + vy * vy ); + *(vxvv+3)= *(vxvv+4); + *(vxvv+4)= atan2( y , vy ); + } else { // surface: x=0, vx > 0 + *(vxvv )= R * sp; + *(vxvv+1)= vR * sp + vT * cp; + x= R * cp; + vx= vR * cp - vT * sp; + *(vxvv+2)= sqrt ( x * x + vx * vx ); + *(vxvv+3)= *(vxvv+4); + *(vxvv+4)= atan2( x , vx ); + } +} +/* +NAME: sos_to_polar_galpy +PURPOSE: convert (y,vy,A,t,psi) or (x,vx,A,t,psi) to (R,vR,vT,phi,t ) +INPUT: + double * vxvv - (x,y,vx,vy,A,t,psi) + surface - if 1, convert from (y,vy,A,t) instead of (x,vx,A,t) +OUTPUT (as arguments): + performed in-place +HISTORY: 2023-03-24 - Written - Bovy (UofT) + */ +void sos_to_polar_galpy(double *vxvv,int surface){ + double x,y,vx,vy,phi,cp,sp; + if ( surface == 1) { // surface: y=0, vy > 0 + x= *(vxvv ); + vx= *(vxvv+1); + y= *(vxvv+2) * sin ( *(vxvv+4) ); + vy= *(vxvv+2) * cos ( *(vxvv+4) ); + } else { // surface: x=0, vx > 0 + y= *(vxvv ); + vy= *(vxvv+1); + x= *(vxvv+2) * sin ( *(vxvv+4) ); + vx= *(vxvv+2) * cos ( *(vxvv+4) ); + } + phi= atan2( y , x ); + cp = cos ( phi ); + sp = sin ( phi ); + *(vxvv )= sqrt ( x * x + y * y ); + *(vxvv+1)= vx * cp + vy * sp; + *(vxvv+2)= -vx * sp + vy * cp; + *(vxvv+4)= *(vxvv+3); + *(vxvv+3)= phi; +} diff --git a/galpy/util/bovy_coords.h b/galpy/util/bovy_coords.h index 47b2c11c4..86b01c564 100644 --- a/galpy/util/bovy_coords.h +++ b/galpy/util/bovy_coords.h @@ -49,6 +49,8 @@ void cyl_to_rect_galpy(double *); void rect_to_cyl_galpy(double *); void cyl_to_sos_galpy(double *); void sos_to_cyl_galpy(double *); +void polar_to_sos_galpy(double *,int); +void sos_to_polar_galpy(double *,int); #ifdef __cplusplus } #endif From 9cd3015aa8d75489de28a23c1d5660b51f0dbf1a Mon Sep 17 00:00:00 2001 From: Jo Bovy Date: Mon, 27 Mar 2023 21:42:43 -0400 Subject: [PATCH 09/18] Add tests of 2D SOS integration --- tests/test_orbit.py | 57 ++++++++++++++++++++++++++++++++++++++++++++ tests/test_orbits.py | 55 ++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 112 insertions(+) diff --git a/tests/test_orbit.py b/tests/test_orbit.py index 7575ee779..c56de9d50 100644 --- a/tests/test_orbit.py +++ b/tests/test_orbit.py @@ -691,6 +691,57 @@ def test_SOS_3D(): f'vz on SOS is not positive for integrate_sos for method={method}' return None +# Test that integrating an orbit in MWPotential2014 using integrate_SOS conserves energy +def test_integrate_SOS_2D(): + pot= potential.LogarithmicHaloPotential(normalize=1.,q=0.9).toPlanar() + o= setup_orbit_energy(pot,axi=True) + psis= numpy.linspace(0.,20.*numpy.pi,1001) + for method in ['dopr54_c','dop853_c','rk4_c','rk6_c','dop853','odeint']: + for surface in ['x','y']: + o.integrate_SOS(psis,pot,method=method,surface='x') + Es= o.E(o.t) + assert (numpy.std(Es)/numpy.mean(Es))**2. < 10.**-10, \ + f'Energy is not conserved by integrate_sos for method={method} and surface={surface}' + return None + +# Test that the 2D SOS function returns points with x=0, vx > 0 when surface='x' +def test_SOS_2Dx(): + pot= potential.LogarithmicHaloPotential(normalize=1.,q=0.9).toPlanar() + o= setup_orbit_energy(pot) + for method in ['dopr54_c','dop853_c','rk4_c','rk6_c','dop853','odeint']: + o.SOS( + pot, + method=method,ncross=500 if '_c' in method else 20, + force_map='rk' in method, + surface='x' + ) + xs= o.x(o.t) + vxs= o.vx(o.t) + assert (numpy.fabs(xs) < 10.**-6.).all(), \ + f'x on SOS is not zero for integrate_sos for method={method}' + assert (vxs > 0.).all(), \ + f'vx on SOS is not positive for integrate_sos for method={method}' + return None + +# Test that the 2D SOS function returns points with y=0, vy > 0 when surface='y' +def test_SOS_2Dy(): + pot= potential.LogarithmicHaloPotential(normalize=1.,q=0.9).toPlanar() + o= setup_orbit_energy(pot) + for method in ['dopr54_c','dop853_c','rk4_c','rk6_c','dop853','odeint']: + o.SOS( + pot, + method=method,ncross=500 if '_c' in method else 20, + force_map='rk' in method, + surface='y' + ) + ys= o.y(o.t) + vys= o.vy(o.t) + assert (numpy.fabs(ys) < 10.**-7.).all(), \ + f'y on SOS is not zero for integrate_sos for method={method}' + assert (vys > 0.).all(), \ + f'vy on SOS is not positive for integrate_sos for method={method}' + return None + # Test that the eccentricity of circular orbits is zero def test_eccentricity(): #return None @@ -4793,10 +4844,16 @@ def test_full_plotting(): return None def test_plotSOS(): + # 3D pot= potential.MWPotential2014 o= setup_orbit_energy(pot) o.plotSOS(pot) o.plotSOS(pot,use_physical=True) + # 2D + pot= potential.LogarithmicHaloPotential(normalize=1.,q=0.9).toPlanar() + o= setup_orbit_energy(pot) + o.plotSOS(pot) + o.plotSOS(pot,use_physical=True) return None def test_from_name_values(): diff --git a/tests/test_orbits.py b/tests/test_orbits.py index b15020b9a..4c73d9a12 100644 --- a/tests/test_orbits.py +++ b/tests/test_orbits.py @@ -432,6 +432,54 @@ def test_SOS_3D(): f'vz on SOS is not positive for integrate_sos for method={method}' return None +# Test that the 2D SOS function returns points with x=0, vx > 0 +def test_SOS_2Dx(): + from galpy.orbit import Orbit + times= numpy.linspace(0.,10.,1001) + orbits_list= [Orbit([1.,0.1,1.,0.]),Orbit([.9,0.3,1.,3.]), + Orbit([1.2,-0.3,0.7,6.])] + orbits= Orbit(orbits_list) + pot= potential.LogarithmicHaloPotential(normalize=1.,q=0.9).toPlanar() + for method in ['dopr54_c','dop853_c','rk4_c','rk6_c','dop853','odeint']: + orbits.SOS( + pot, + method=method,ncross=500 if '_c' in method else 20, + force_map='rk' in method, + t0=numpy.arange(len(orbits)), + surface='x', + ) + xs= orbits.x(orbits.t) + vxs= orbits.vx(orbits.t) + assert (numpy.fabs(xs) < 10.**-6.).all(), \ + f'x on SOS is not zero for integrate_sos for method={method}' + assert (vxs > 0.).all(), \ + f'vx on SOS is not positive for integrate_sos for method={method}' + return None + +# Test that the 2D SOS function returns points with y=0, vy > 0 +def test_SOS_2Dy(): + from galpy.orbit import Orbit + times= numpy.linspace(0.,10.,1001) + orbits_list= [Orbit([1.,0.1,1.,0.]),Orbit([.9,0.3,1.,3.]), + Orbit([1.2,-0.3,0.7,6.])] + orbits= Orbit(orbits_list) + pot= potential.LogarithmicHaloPotential(normalize=1.,q=0.9).toPlanar() + for method in ['dopr54_c','dop853_c','rk4_c','rk6_c','dop853','odeint']: + orbits.SOS( + pot, + method=method,ncross=500 if '_c' in method else 20, + force_map='rk' in method, + t0=numpy.arange(len(orbits)), + surface='y', + ) + ys= orbits.y(orbits.t) + vys= orbits.vy(orbits.t) + assert (numpy.fabs(ys) < 10.**-6.).all(), \ + f'y on SOS is not zero for integrate_sos for method={method}' + assert (vys > 0.).all(), \ + f'vy on SOS is not positive for integrate_sos for method={method}' + return None + # Test slicing of orbits def test_slice_singleobject(): from galpy.orbit import Orbit @@ -2138,10 +2186,17 @@ def test_plotting(): def test_plotSOS(): from galpy.orbit import Orbit from galpy.potential import LogarithmicHaloPotential + + # 3D o= Orbit([Orbit([1.,0.1,1.1,0.1,0.2,2.]),Orbit([1.,0.1,1.1,0.1,0.2,2.])]) pot= potential.MWPotential2014 o.plotSOS(pot) o.plotSOS(pot,use_physical=True,ro=8.,vo=220.) + # 2D + o= Orbit([Orbit([1.,0.1,1.1,2.]),Orbit([1.,0.1,1.1,2.])]) + pot= LogarithmicHaloPotential(normalize=1.).toPlanar() + o.plotSOS(pot) + o.plotSOS(pot,use_physical=True,ro=8.,vo=220.) return None def test_integrate_method_warning(): From 2b7d71ac091c29e6902f141d1d8f3954c55f704e Mon Sep 17 00:00:00 2001 From: Jo Bovy Date: Tue, 28 Mar 2023 21:29:12 -0400 Subject: [PATCH 10/18] No quantity use when plotting the SOS (similar to plot) --- galpy/orbit/Orbits.py | 1 + 1 file changed, 1 insertion(+) diff --git a/galpy/orbit/Orbits.py b/galpy/orbit/Orbits.py index a275a2891..1ca11dc82 100644 --- a/galpy/orbit/Orbits.py +++ b/galpy/orbit/Orbits.py @@ -5370,6 +5370,7 @@ def plotSOS(self,pot,*args,ncross=500,surface=None,t0=0., else: d1= 'y' d2= 'vy' + kwargs['quantity']= False x,y= self.SOS(pot,ncross=ncross,surface=surface, t0=t0,method=method,skip=skip, progressbar=progressbar,**kwargs) From 65a0429e2fffff089a540c2fccf204dd596e3635 Mon Sep 17 00:00:00 2001 From: Jo Bovy Date: Mon, 27 Mar 2023 23:53:03 -0400 Subject: [PATCH 11/18] SOS test tweaks to increase coverage --- tests/test_orbit.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/tests/test_orbit.py b/tests/test_orbit.py index c56de9d50..0f63047f3 100644 --- a/tests/test_orbit.py +++ b/tests/test_orbit.py @@ -698,7 +698,7 @@ def test_integrate_SOS_2D(): psis= numpy.linspace(0.,20.*numpy.pi,1001) for method in ['dopr54_c','dop853_c','rk4_c','rk6_c','dop853','odeint']: for surface in ['x','y']: - o.integrate_SOS(psis,pot,method=method,surface='x') + o.integrate_SOS(psis,pot,method=method) # default is surface='x' Es= o.E(o.t) assert (numpy.std(Es)/numpy.mean(Es))**2. < 10.**-10, \ f'Energy is not conserved by integrate_sos for method={method} and surface={surface}' @@ -4854,6 +4854,8 @@ def test_plotSOS(): o= setup_orbit_energy(pot) o.plotSOS(pot) o.plotSOS(pot,use_physical=True) + o.plotSOS(pot,surface='y') + o.plotSOS(pot,surface='y',use_physical=True) return None def test_from_name_values(): From 2ec29bae39792ab5e0eb7a8b94448b461d73cd5a Mon Sep 17 00:00:00 2001 From: Jo Bovy Date: Tue, 28 Mar 2023 19:25:33 -0400 Subject: [PATCH 12/18] Check that an orbit leaves the SOS --- galpy/orbit/Orbits.py | 5 +++++ tests/test_orbit.py | 35 +++++++++++++++++++++++++++++++++++ 2 files changed, 40 insertions(+) diff --git a/galpy/orbit/Orbits.py b/galpy/orbit/Orbits.py index 1ca11dc82..cbee43f83 100644 --- a/galpy/orbit/Orbits.py +++ b/galpy/orbit/Orbits.py @@ -4775,6 +4775,11 @@ def SOS(self,pot,ncross=500,surface=None,t0=0., numcores=numcores,force_map=force_map) old_vxvv= self.vxvv self.vxvv= self.orbit[:,-1] + # We are on the SOS now. Let's check that v(x/y/z) > 0 + if ( ( self.dim() == 3 and not self.vz() > 0. ) + or ( self.dim() == 2 and surface.lower() == 'y' and not self.vy() > 0. ) + or ( self.dim() == 2 and surface.lower() == 'x' and not self.vx() > 0. ) ): + raise RuntimeError("Orbit appears to be within the SOS surface. Refusing to perform specialized SOS integration, please use normal integration instead") if method == 'rk4_c' or method == 'rk6_c': # Because these are non-adaptive, we need to make sure we # integrate finely enough diff --git a/tests/test_orbit.py b/tests/test_orbit.py index 0f63047f3..e43100756 100644 --- a/tests/test_orbit.py +++ b/tests/test_orbit.py @@ -742,6 +742,41 @@ def test_SOS_2Dy(): f'vy on SOS is not positive for integrate_sos for method={method}' return None +# Test that the SOS integration returns an error +# when the orbit does not leave the surface +def test_SOS_onsurfaceerror_3D(): + from galpy.orbit import Orbit + o= Orbit([1.,0.1,1.1,0.,0.,0.]) + with pytest.raises(RuntimeError,match="Orbit appears to be within the SOS surface. Refusing to perform specialized SOS integration, please use normal integration instead"): + o.SOS(potential.MWPotential2014,surface='y') + return None + +# Test that the SOS integration returns an error +# when the orbit does not leave the surface +def test_SOS_onsurfaceerror_2D(): + from galpy.orbit import Orbit + + # An orbit considered in the book + def orbit_xvxE(x,vx,E,pot=None): + """Returns Orbit at (x,vx,y=0) with given E""" + return Orbit([x,vx,numpy.sqrt(2.*(E-potential.evaluatePotentials(pot,x,0.,phi=0.) + -vx**2./2)),0.]) + # Need the 2d zvc here (maybe should add to galpy?) + def zvc(x,E,pot=None): + """Returns the maximum v_x at this x and this + energy: the zero-velocity curve""" + return numpy.sqrt(2.*(E-potential.evaluatePotentials(pot,x,0.,phi=0.))) + lp= potential.LogarithmicHaloPotential(normalize=True,b=0.9,core=0.2) + E= -0.87 + x= 0.204 + # This orbit remains in the y=0 plane and psi therefore + # remains zero, thus not increasing + maxvx= zvc(x,E,pot=lp) + o= orbit_xvxE(x,maxvx,E,pot=lp) + with pytest.raises(RuntimeError,match="Orbit appears to be within the SOS surface. Refusing to perform specialized SOS integration, please use normal integration instead"): + o.SOS(lp,surface='y') + return None + # Test that the eccentricity of circular orbits is zero def test_eccentricity(): #return None From a28b8950cb7903748ea410765331739a65e58924 Mon Sep 17 00:00:00 2001 From: Jo Bovy Date: Tue, 28 Mar 2023 19:30:44 -0400 Subject: [PATCH 13/18] Don't support 1D SOS, fix surface parsing issue --- galpy/orbit/Orbits.py | 59 ++++++++++++++++++------------------------- 1 file changed, 24 insertions(+), 35 deletions(-) diff --git a/galpy/orbit/Orbits.py b/galpy/orbit/Orbits.py index cbee43f83..94747d6fe 100644 --- a/galpy/orbit/Orbits.py +++ b/galpy/orbit/Orbits.py @@ -1369,6 +1369,8 @@ def integrate_SOS(self,psi,pot,surface=None,t0=0., 2023-03-16 - Written - Bovy (UofT) """ + if self.dim() == 1: + raise NotImplementedError("SOS integration is not supported for 1D orbits") self.check_integrator(method,no_symplec=True) pot= flatten_potential(pot) _check_potential_dim(self,pot) @@ -1396,12 +1398,7 @@ def integrate_SOS(self,psi,pot,surface=None,t0=0., method= self._check_method_dissipative_compatible(method,self._pot) # Implementation with parallel_map in Python if not '_c' in method or not ext_loaded or force_map: - if self.dim() == 1: - raise NotImplementedError("SOS integration not implemented for 1D orbits") - out, msg= integrateLinearOrbit(self._pot,self.vxvv,t,method, - progressbar=progressbar, - numcores=numcores,dt=dt) - elif self.dim() == 2: + if self.dim() == 2: out, msg= integratePlanarOrbit_sos(self._pot,self.vxvv,self._psi,t0,method, surface=surface,progressbar=progressbar, numcores=numcores) @@ -1412,34 +1409,26 @@ def integrate_SOS(self,psi,pot,surface=None,t0=0., else: warnings.warn("Using C implementation to integrate orbits", galpyWarningVerbose) - if self.dim() == 1: - raise NotImplementedError("SOS integration not implemented for 1D orbits in C") - out, msg= integrateLinearOrbit_c(self._pot, - numpy.copy(self.vxvv), - t,method, - progressbar=progressbar, - dt=dt) + if self.phasedim() == 3 \ + or self.phasedim() == 5: + #We hack this by putting in a dummy phi=0 + vxvvs= numpy.pad(self.vxvv,((0,0),(0,1)), + 'constant',constant_values=0) else: - if self.phasedim() == 3 \ - or self.phasedim() == 5: - #We hack this by putting in a dummy phi=0 - vxvvs= numpy.pad(self.vxvv,((0,0),(0,1)), - 'constant',constant_values=0) - else: - vxvvs= numpy.copy(self.vxvv) - if self.dim() == 2: - out, msg= integratePlanarOrbit_sos_c(self._pot,vxvvs,self._psi,t0, - method,surface=surface, - progressbar=progressbar) - else: - out, msg= integrateFullOrbit_sos_c(self._pot,vxvvs,self._psi,t0, - method,progressbar=progressbar) - - if self.phasedim() == 3 \ - or self.phasedim() == 5: - phi_mask= numpy.ones(out.shape[2],dtype='bool') - phi_mask[3+2*(self.phasedim()==5)]= False - out= out[:,:,phi_mask] + vxvvs= numpy.copy(self.vxvv) + if self.dim() == 2: + out, msg= integratePlanarOrbit_sos_c(self._pot,vxvvs,self._psi,t0, + method,surface=surface, + progressbar=progressbar) + else: + out, msg= integrateFullOrbit_sos_c(self._pot,vxvvs,self._psi,t0, + method,progressbar=progressbar) + + if self.phasedim() == 3 \ + or self.phasedim() == 5: + phi_mask= numpy.ones(out.shape[2],dtype='bool') + phi_mask[3+2*(self.phasedim()==5)]= False + out= out[:,:,phi_mask] # Store orbit internally self.orbit= out[:,:,:-1] self.t= out[:,:,-1] @@ -4777,8 +4766,8 @@ def SOS(self,pot,ncross=500,surface=None,t0=0., self.vxvv= self.orbit[:,-1] # We are on the SOS now. Let's check that v(x/y/z) > 0 if ( ( self.dim() == 3 and not self.vz() > 0. ) - or ( self.dim() == 2 and surface.lower() == 'y' and not self.vy() > 0. ) - or ( self.dim() == 2 and surface.lower() == 'x' and not self.vx() > 0. ) ): + or ( self.dim() == 2 and not surface is None and surface.lower() == 'y' and not self.vy() > 0. ) + or ( self.dim() == 2 and (surface is None or surface.lower() == 'x') and not self.vx() > 0. ) ): raise RuntimeError("Orbit appears to be within the SOS surface. Refusing to perform specialized SOS integration, please use normal integration instead") if method == 'rk4_c' or method == 'rk6_c': # Because these are non-adaptive, we need to make sure we From 426f7ac1086c5342cac8d23e55063016795f64c4 Mon Sep 17 00:00:00 2001 From: Jo Bovy Date: Tue, 28 Mar 2023 19:41:02 -0400 Subject: [PATCH 14/18] Test that slicing of individual time arrays works as expected --- galpy/orbit/Orbits.py | 8 ++++---- tests/test_orbits.py | 19 +++++++++++++++++++ 2 files changed, 23 insertions(+), 4 deletions(-) diff --git a/galpy/orbit/Orbits.py b/galpy/orbit/Orbits.py index 94747d6fe..df863b239 100644 --- a/galpy/orbit/Orbits.py +++ b/galpy/orbit/Orbits.py @@ -4765,10 +4765,10 @@ def SOS(self,pot,ncross=500,surface=None,t0=0., old_vxvv= self.vxvv self.vxvv= self.orbit[:,-1] # We are on the SOS now. Let's check that v(x/y/z) > 0 - if ( ( self.dim() == 3 and not self.vz() > 0. ) - or ( self.dim() == 2 and not surface is None and surface.lower() == 'y' and not self.vy() > 0. ) - or ( self.dim() == 2 and (surface is None or surface.lower() == 'x') and not self.vx() > 0. ) ): - raise RuntimeError("Orbit appears to be within the SOS surface. Refusing to perform specialized SOS integration, please use normal integration instead") + if ( ( self.dim() == 3 and not numpy.all(self.vz() > 0.) ) + or ( self.dim() == 2 and not surface is None and surface.lower() == 'y' and not numpy.all(self.vy() > 0.) ) + or ( self.dim() == 2 and (surface is None or surface.lower() == 'x') and not numpy.all(self.vx() > 0.) ) ): + raise RuntimeError("An orbit appears to be within the SOS surface. Refusing to perform specialized SOS integration, please use normal integration instead") if method == 'rk4_c' or method == 'rk6_c': # Because these are non-adaptive, we need to make sure we # integrate finely enough diff --git a/tests/test_orbits.py b/tests/test_orbits.py index 4c73d9a12..f9a9feabb 100644 --- a/tests/test_orbits.py +++ b/tests/test_orbits.py @@ -685,6 +685,25 @@ def test_slice_physical_issue385(): assert numpy.amax(numpy.fabs((((orbits[ii].phi()-orbits.phi()[ii])+numpy.pi) % (2.*numpy.pi)) - numpy.pi)) < 1e-10, 'Integration of multiple orbits as Orbits does not agree with integrating multiple orbits' return None +# Test that slicing in the case of individual time arrays works as expected +# Currently, the only way individual time arrays occur is through SOS integration +# so we implementing this test using SOS integration +def test_slice_indivtimes(): + from galpy.orbit import Orbit + times= numpy.linspace(0.,10.,1001) + orbits_list= [Orbit([1.,0.1,1.,0.,0.1,0.]),Orbit([.9,0.3,1.,-0.3,0.4,3.]), + Orbit([1.2,-0.3,0.7,.5,-0.5,6.])] + orbits= Orbit(orbits_list) + pot= potential.MWPotential2014 + orbits.SOS(pot,t0=numpy.arange(len(orbits))) + # First check that we actually have individual times + assert len(orbits.t.shape) >= len(orbits.orbit.shape)-1, 'Test should be using individual time arrays, but a single time array was found' + # Now slice single and multiple + assert numpy.all(orbits[0].t == orbits.t[0]), 'Individually sliced orbit with individual time arrays does not produce the correct time array in the slice' + assert numpy.all(orbits[1].t == orbits.t[1]), 'Individually sliced orbit with individual time arrays does not produce the correct time array in the slice' + assert numpy.all(orbits[:2].t == orbits.t[:2]), 'Multiply-sliced orbit with individual time arrays does not produce the correct time array in the slice' + assert numpy.all(orbits[1:4].t == orbits.t[1:4]), 'Multiply-sliced orbit with individual time arrays does not produce the correct time array in the slice' + return None # Test that initializing Orbits with orbits with different phase-space # dimensions raises an error From e377de964a6cdef0f13de874d72e6756fb469325 Mon Sep 17 00:00:00 2001 From: Jo Bovy Date: Tue, 28 Mar 2023 19:56:30 -0400 Subject: [PATCH 15/18] Check that time and psi increase together in SOS integration --- galpy/orbit/Orbits.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/galpy/orbit/Orbits.py b/galpy/orbit/Orbits.py index df863b239..a2087694c 100644 --- a/galpy/orbit/Orbits.py +++ b/galpy/orbit/Orbits.py @@ -1432,6 +1432,11 @@ def integrate_SOS(self,psi,pot,surface=None,t0=0., # Store orbit internally self.orbit= out[:,:,:-1] self.t= out[:,:,-1] + # Quick check that all is well in terms of psi increasing with time + assert numpy.all((numpy.roll(self.t,-1,axis=1)-self.t)[:,:-1] + *(numpy.roll(self._psi.T,-1,axis=0)-self._psi.T)[:-1].T + > 0.), \ + "SOS integration failed (time does not monotonically increase with increasing psi)" return None def integrate_dxdv(self,dxdv,t,pot,method='dopr54_c', From 6d88b85d5fe319311faf93b73a792074a8cb9ec1 Mon Sep 17 00:00:00 2001 From: Jo Bovy Date: Tue, 28 Mar 2023 20:11:05 -0400 Subject: [PATCH 16/18] Move SOS surface check before any integration and test for multiple orbits --- galpy/orbit/Orbits.py | 15 ++++++++++----- tests/test_orbit.py | 6 +++--- tests/test_orbits.py | 9 +++++++++ 3 files changed, 22 insertions(+), 8 deletions(-) diff --git a/galpy/orbit/Orbits.py b/galpy/orbit/Orbits.py index a2087694c..be566598d 100644 --- a/galpy/orbit/Orbits.py +++ b/galpy/orbit/Orbits.py @@ -4759,6 +4759,16 @@ def SOS(self,pot,ncross=500,surface=None,t0=0., init_psis= numpy.arctan2(self.x(use_physical=False),self.vx(use_physical=False)) else: raise NotImplementedError("SOS not implemented for 1D orbits or 2D orbits without phi") + # Let's check that v(x/y/z) != 0 for orbits that are already on the SOS + if ( ( self.dim() == 3 + and not numpy.all((self.vz() != 0.)+(numpy.fabs(init_psis % numpy.pi) > 1e-10)) ) + or ( self.dim() == 2 + and not surface is None and surface.lower() == 'y' + and not numpy.all((self.vy() != 0.)+(numpy.fabs(init_psis % numpy.pi) > 1e-10)) ) + or ( self.dim() == 2 + and (surface is None or surface.lower() == 'x') + and not numpy.all((self.vx() != 0.)+(numpy.fabs(init_psis % numpy.pi) > 1e-10)) ) ): + raise RuntimeError("An orbit appears to be within the SOS surface. Refusing to perform specialized SOS integration, please use normal integration instead") if numpy.any(numpy.fabs(init_psis) > 1e-10): # Integrate to the next crossing init_psis= numpy.atleast_1d((init_psis + 2.*numpy.pi) % (2.*numpy.pi)) @@ -4769,11 +4779,6 @@ def SOS(self,pot,ncross=500,surface=None,t0=0., numcores=numcores,force_map=force_map) old_vxvv= self.vxvv self.vxvv= self.orbit[:,-1] - # We are on the SOS now. Let's check that v(x/y/z) > 0 - if ( ( self.dim() == 3 and not numpy.all(self.vz() > 0.) ) - or ( self.dim() == 2 and not surface is None and surface.lower() == 'y' and not numpy.all(self.vy() > 0.) ) - or ( self.dim() == 2 and (surface is None or surface.lower() == 'x') and not numpy.all(self.vx() > 0.) ) ): - raise RuntimeError("An orbit appears to be within the SOS surface. Refusing to perform specialized SOS integration, please use normal integration instead") if method == 'rk4_c' or method == 'rk6_c': # Because these are non-adaptive, we need to make sure we # integrate finely enough diff --git a/tests/test_orbit.py b/tests/test_orbit.py index e43100756..1c5c47be5 100644 --- a/tests/test_orbit.py +++ b/tests/test_orbit.py @@ -747,8 +747,8 @@ def test_SOS_2Dy(): def test_SOS_onsurfaceerror_3D(): from galpy.orbit import Orbit o= Orbit([1.,0.1,1.1,0.,0.,0.]) - with pytest.raises(RuntimeError,match="Orbit appears to be within the SOS surface. Refusing to perform specialized SOS integration, please use normal integration instead"): - o.SOS(potential.MWPotential2014,surface='y') + with pytest.raises(RuntimeError,match="An orbit appears to be within the SOS surface. Refusing to perform specialized SOS integration, please use normal integration instead"): + o.SOS(potential.MWPotential2014) return None # Test that the SOS integration returns an error @@ -773,7 +773,7 @@ def zvc(x,E,pot=None): # remains zero, thus not increasing maxvx= zvc(x,E,pot=lp) o= orbit_xvxE(x,maxvx,E,pot=lp) - with pytest.raises(RuntimeError,match="Orbit appears to be within the SOS surface. Refusing to perform specialized SOS integration, please use normal integration instead"): + with pytest.raises(RuntimeError,match="An orbit appears to be within the SOS surface. Refusing to perform specialized SOS integration, please use normal integration instead"): o.SOS(lp,surface='y') return None diff --git a/tests/test_orbits.py b/tests/test_orbits.py index f9a9feabb..9407494db 100644 --- a/tests/test_orbits.py +++ b/tests/test_orbits.py @@ -480,6 +480,15 @@ def test_SOS_2Dy(): f'vy on SOS is not positive for integrate_sos for method={method}' return None +# Test that the SOS integration returns an error +# when one orbit does not leave the surface +def test_SOS_onsurfaceerror_3D(): + from galpy.orbit import Orbit + o= Orbit([[1.,0.1,1.1,0.1,0.,0.],[1.,0.1,1.1,0.,0.,0.]]) + with pytest.raises(RuntimeError,match="An orbit appears to be within the SOS surface. Refusing to perform specialized SOS integration, please use normal integration instead"): + o.SOS(potential.MWPotential2014) + return None + # Test slicing of orbits def test_slice_singleobject(): from galpy.orbit import Orbit From 4c10de188ba81664f3ee18dc71c70413408df66a Mon Sep 17 00:00:00 2001 From: Jo Bovy Date: Tue, 28 Mar 2023 20:12:16 -0400 Subject: [PATCH 17/18] Add SOS to HISTORY --- HISTORY.txt | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/HISTORY.txt b/HISTORY.txt index 25c64120c..5e0c31903 100644 --- a/HISTORY.txt +++ b/HISTORY.txt @@ -1,6 +1,10 @@ v1.9.0 (XXXX-XX-XX) =================== +- Add specialized integration method to determine surfaces of sections of orbits and + add Orbit.SOS and Orbit.plotSOS methods to compute and plot the surface of section + of an orbit. + v1.8.3 (2023-03-27) =================== From 9bef72b89695f74caf73415c3a1f44492a349501 Mon Sep 17 00:00:00 2001 From: Jo Bovy Date: Tue, 28 Mar 2023 20:52:15 -0400 Subject: [PATCH 18/18] Add documentation on surfaces of section --- doc/source/images/orbit-sos-2Dbox.png | Bin 0 -> 13373 bytes doc/source/images/orbit-sos-Sun-MWP14.png | Bin 0 -> 12329 bytes doc/source/images/orbit-sos-twoorb-MWP14.png | Bin 0 -> 20082 bytes doc/source/orbit.rst | 101 +++++++++++++++++++ 4 files changed, 101 insertions(+) create mode 100644 doc/source/images/orbit-sos-2Dbox.png create mode 100644 doc/source/images/orbit-sos-Sun-MWP14.png create mode 100644 doc/source/images/orbit-sos-twoorb-MWP14.png diff --git a/doc/source/images/orbit-sos-2Dbox.png b/doc/source/images/orbit-sos-2Dbox.png new file mode 100644 index 0000000000000000000000000000000000000000..0ca62238c3579acde3e4ed86b27ebaac33521ffe GIT binary patch literal 13373 zcmajG1yog0+cmoBjzb8D)S*FADFvwmasWXY9sR!de)sc6iLJ&qb_?5%Q z1)sc8;;;vQBs>-LJhfb`J)fAlTR{)ZJY5}KJRR-KnZ2#tJ?vbZZ}E%r3*Th6_4IW0 zxFaCo^q&*>UEFO1xPH-^gPRbzD(QPb5D_K%7Xx~eP76V|3XureM?UGBGfy3FSY>kV zY-`PkCX`1z-H@kv$?zum_dK)S)XV;@n;%s%@2LF|r0lB1E(v{+8&Wl5@A;12+qznV zYRy_eiSw1(^$5@LA4*)^GD@@^O7*PPCZ>m|h$+0NTORT)5wDQ=5q5VckR0Amq>k=u zR~Pkad)I!NZP+-A%gLQu-b-+Qa&hW=a>2-whYf?l`tz}8AUHgy+&q>I4tH;1Mj{YP z!LsovR2N1F8w?hQg^Ps4;dr~*^z*ezyt+Nr4r7?m`~3Xj@o~N3lYuPhn92QoD8R%#VmzWi8OmoI1@D$*VUgqslDwuRr5fnC!krm zZupwmz{05Sw0#u)SnahkqQ{?RGdNYtjJewmf(W(H7elo1qC%`AZ6)27!oN7Ez!kNZ z(q`a%``w1`lLsWKkxW4%OoM~7nWH@mTcWaT`gyEOfs z8ygHJhAz9Oe!V3d4xTB8TQhf^jt;sG*!=zc*raE`{G^h@S@V8Z++$kL%q7nukRfd0 z3dRn`;9f3ot1yO2v4(P)h`5(Mnr|nfQ1+T1LfC>spc)nK`xr9#ZjbOQ^pz#7JbtG^ zqDtI&e=PFITZKMQz+ikFRYiseCYG)Eu?-p!Jwvlq2#ow-8z-2%Z1O$=Q7up+NCk&QDrt_fX=1nb_!c#T4RxM{}d?1Z@!wp@?=wg1XHl$^j;svj>mF|EnrIt->D zn-ap7&G(f%b=E_bh=n;++?IxG5G-P5W<8E-!J|l(>TZwJikqZl>N1#mpWnkc2UKY; zGnynrnh`{SR}<_v=p>@TN6T=7eXq8RnlO_~SC}y*Cx75V;-YsVR+bOK_DkwGqA6_z ze!e{@s&LlEihiYbEn_MmI`7ia;f(iU0wNq+87>v|d~*A2U*)ZmO8&P|&z2($nuA zOKnEsPTv1C8CiPvde&R_ij0a+jhoqwH2wjkhU89rFt~}~`bhaCa(SE|!sEjb!791- zSVP7tHdm8QcBnLewB8RVaBsk_Z|OS@maHU4?4U}gj1B>YkGqL9-=848n?oH`%w+>o zmpzl#(|+m-?@zmotQZdR)=}+F`BCUi6%vYOfR%>O>5yNR*9rGcL zk3W~(1_2)x{!3Gs7WSq6%t09a?NzqbL1duM7DaT{z`!SOwXw0aCwPgYy*k(rnRLq+ z7u2S%u)JyK6uNnGSFWDCsBEY53oSp{s;{{70 z_l&utTtWl$KSc2!A7*{*RI2VQt+g|eA`;T=Gv&xPKZ=`Mc+)^ThsW|IldZmoKL*C- zv`Zo_rIkp_1fM&zB)lKwKX6FfEPdXV$>!^PH7_2xvyKNv5(N|Wx)}w?V*^xp_vsxS z0H8ppLR`yBNxRJ2rRe}cI1O8Vetw*CUJB1uK3k^QcRDott=~NJr3BR_c@*SF?Fw`kOr|@(a9Ao{Hyd4=T;`nczw%;9oXdy zy<+n{j+xu-?567RJ=RZez40nFd4_ja*9{fkWuv(BalnvJVV*=PX$i#pPSIRY^ycGF zVcUa*6`aU@+I?o1^kd>}b=4`x{Wh1GZuLS0OZm{5NGbkrKTD~(BJ@p_%}o;Dakbxy z?sH)yWU`p8^!1gR&EjG_=`NtSYeJ6!jlM9>>B2#-ASXW>C%#jb^ub#dM-@KXy$o18 z7Qc7w*-*YNS*Cpc&4IO;WFy!&{@M-ajKMocwkC;fKim;NY`|)H|9YmuhF<)Av z_*=<)@2S%lNe~lD>3w?F!oxXtR<;KgK1S&h4bJ@-?NX3HXgj~1EKe^=kMW3q#CGd- zJEi|BA!dqHp!(f8Nrd)*0ljxel^yRYk4h&qoh*9QZ|^KLLiIIP5}HRdtq}5&{L9Jm z@E4prBfIK^@@2kK4+}Lb$-iLIvaxpPHTikfh|7$2VqbPXp6eGQ&>#`^OWGK_eQkw@ zT0tg>q7d8KTKn8~UOvEPakCLmPf$rQSm?;s(NN?ajj9ebJu!+J0CZrhOFwr5A>yp! za(A@EzNKOqd0U{Y6X36@wN4n}ycnmzd`1kXgs$BeX8pxtj~F)`R43KaKF|73Lfv|p z!pnNb-woO%$HOh*$Nk+zu_~Rl)}ST;v)dx`SyH&EPHC0 z_B*~NgZ2`zdt5RknNFzUhDx1(U>X0FJ$(n)`a8Wbj+=M&$a3VpJ4dgxwp+~2`)FKf zLq4__>%j&Jh##c@EzW8V7X+LQ2{5SgP?jB)vSNc`Z%N`IKS)(SGYxYm_&}zqwr*n2 zvvuXnVE`;|-(+HP)JWlN_#UnnqLqOy4GJ*xb*5EDsfRL5LfE~=UOG1an;j6G`&zu4 z(L^4x>CZ>gZaaKApS6Oqd_97GdmEI0o1Bca=jSl}4SQm}PH64U^Qb_3|MSLsj`Nv| zTHF~)c%goWP#{@J$S#g)xNt0RrF5 z_*m=5e`q8kB@fr(aSr>h6`01WTlpwYfJs@69@lS4RKMTnlE2GoSN!nDDh-JggbVRO zZ$5K`j^ms5P%)=Lo+qQGOgFA)wy2t{7xe{~&CK>(FGegC0e11(T{F3TvxxVlHJxu_ zb=|s1pA=HYu;LAc+)KF)ZRcD4x17UAKTen8^}d6cLj{9^YF+pa#xK9Powvs0`95tE zxna(pBlCyd04NCidol4i14Y5c6?><3&#|V)L)dcW!MZ3II$_8d!F=re51ex5AI(cv zwUylNG$;;m>UGCef)g51?Qzku$=?K!6BBd-yQ5pi<#!KkcH&qB6=8MeeY!&l${#7Z(@7ZGfFOkfNfZ@C>!^ zJ2IFC^%54#}Q z{8D5k0urCTJLM0w8C?Kt*cbY z2rU<{YDdVoQDW5CuF}-h{R0CBV0O>kFJ8UsDx#ohQB`YS zh`mOKGKwXxOo9F$ET%h^wfZffqsr_qE8O)UxumV|hBbT62fwN!k^>0i^KQltDbBVgQnWI!Ra}HD`Yed$Lr^4?SugPAtqvTX0C)%=`gpRF7b>Xa z#|7dhdQ34^>!4kO!d6d*#_=5rCC7D>n?(nx)W8)aBYj4wy@6CQ*Ymqi!7zfFa*)&a z&DOl4X5xDHZ&pOoS~7sfEIJ)n!DI<3n2xb-NPJgDESZ80d5QdNgDHb%v!1N})cUo` zL4ZfU{VuSv*kBqFK*whndN9;Z1ni{Ux~1f|H2#YXfF2vf%?=Jdccy@8d^iBuK$AtJ z;#D!9|3S`JRp$ReSw7DeDIssHzkHP!NBVL)Ya%_ep^oBM;0=oFH(D>68Zv?Wc79() z*A2jptQhySEJcefG?kZ#MkH~b(D9pj>6MIJq2nK}ZY@}xInuHJLT}~kuYbH17^WyB zq(C6dtfqOh0v2e&vyqahW9t->;_tc^%|@xxVrYsHn!cOa)$GA&Iyw#`>-Vrs%=QPO zRPt(+>rL|5Eg1?}(wCG{;p9RhWn16}FhW|dSKmXD^ZqJ5Efk&zL=u5kL4zP(H%H$# zwHjt!X?#MabP%JIlvJ4X={4?BGlDztoStKiLjDXjns!VGO&w?HmI@fdY7P#(jXeo{ zFRgaVuHGkQQ$jT#u^C{)@jzEE?(AS&?;9wiGztew)`}kSSW&?s2%m+2oPlN~bBU^I zHdmhrGt+J$T->=UP3Fd$PFe)d3!wb#LxU^|jm?c6_q`Env`xh?Cc3>5R=hv;Bo32-!-+y)MIX*3zY4UnIm12I&33*hy9lTx1V~ zhsCPcDC>LtCHC@*-{a@r4LWbnE*QS?d`L2Ho$1r)$X03;zT8=-|i2 zYBqf0vB3Z$5o~>Pn@Y-w-dlaI*x>G8B5r*PGbnW)RCx^t@`VgFC=sx!-@9lAwW)Rb zbKiQ+$JRXBnIuE`I09vl9NLe6%ivtXVWu3w4liWyXv~!n+>a-Y=eIR6AydHT!Yuv< z_9=Y7F&Ddh!z*Z&9YFZZtfvHsh=Hsa=6q#*=q}^d(nDTYToih*DzA0h&X}k1-{)8Q z#}BNjstxY(h)SDjW%l%OX>0Pd?CWJ;GX1q=5PMa%nxW z^=UowcU(%cyu{~u94}jmnI6Li-@E_N%ZQGRx5c0YeV{PfVM(VJfOWU!NYd@t z&W0)GY#h~r`(X6Ws5_;A-V|pWoyKmv%Z_v<6em`j@cM%r&hD9;jQkSZocT27tw9iO zn2b6wh(MiFBN%jXyV=c~0AvKw?Z-#Ppyy7d#SaYir93V94@Jllg`jjiii zwUF)xB-T|h`PrlJFy-HYCG#^RV4a3%9`Rr-p;Z+i{232id!l6B3G_ZDp}5358kclk z3?nk;^tSB?O_q8_>cd0~!RT>xz2;*ZoPr5jZs&h!^lA&gd^{5ju}uAYm7ceiY!hMp z4QD<|hZ+hPPLaUybFf)c$gE^15^g5!_NMk16i+B4sDc4P{+^l)2M?AhkZu`!Q|rn7 zTi%k>5olT(I{ks3dcrQ=RD!;hM>qOZ=UZy-sc&_b64w&4=+wC}N(i|(97-- zaa)mr&=ZAxpq739bQ*{hbDTf##*e|Ee}1L*LIt&%n+H6{lz}m`rMS6l+tK9b&JMzs zW=oZ>tw^&6$nVcO(+S`GUl6eVT2sVuEM6K{uBelMY3!U&2Gsv#F2L&$d%$bN6nwX> zL<)L;bMV31YT43wB*?f?Tt0!+5(55p+K>BBZJ-wzHyyUVA$8?W zZH&}kY3@fDqXncnt5!TZ_c|LK`x%sROng<~K^=vjYC@I(uH7g)!owH%yygH zZ=NsXvV6@=+-q_lTr(v!AR8{m+4p~G#V`f5`P;W|*~RPuuVn_Qq@4y_mS|Mf8btve zta6VD21_afw;LsmC;$0;#;ufEd~B6t=Em%e*7I32%BnEgB=8-;{-Fs_1tGUz)a5$? z@~_>;VdRhTSuShI%@gw2G8|+z)n{law*G>(C! zB!9~|Dy`Ho8vLa%l<3xg#ISsBwA_MBU~&P)K~6D-B;TBHWB*EUZs%LeQb2+mPMfq= zluiJ7`9%G}PhFT9g7@PCo)UmWE(1{;>O}aEMK!QQsOl;-7B0F$#x~e z%cKm~E1V~RCd+Dif(z|@JQ|qi6{{6|NC7+Pb6)R_@^_nPD>(6Mm;aWCIzhl| z029)sNbHjtB`oga=Zj?7#;5N5)w==!CTOy@jQ37L$>}UUl70k+KCe}Bi!uPrJ%{z) zS0$D3b|q60`&z#w^oIjzV<6UF?#OXHO14rJ=*jWgYBFGaNKU;8L!r|3 z^P9k#JHKt*>bxRINuNH)0!j;hK`Lx9k)0-wprQ*@^A}R$|AJ=NC3d#p)W)TbUIgOa zUXe3z%fcx80{Xu%Ye2L`hu@W&Q=u@8HHrs&+{gD9W#0<>y7_Z@0ZP~jTAmGDQl|75 z&`N4wn#3=oK;Dqoj!mo>X!m_j6+oW1j}PW83GKAKo&s^BFNtI3?byIVqo=2jQ_GK2 zQ{%?Hc5$6PS!|E`W8(*p{oQk-Ufg5c1sSMEBe<7_tC;|vwbQm&G`h0KW769l! z{{6>`(>hFy?0co?#p~ki2moXUR30F8j|H`HwAg47y+*|h-%Ar>B7l6ge}1K#($E|$ z=hS}&#ywv9wCJ-JW^t2C2C;-Cw#usVg}EMa(^Czx^clGF51A+@J(D~4;~+pgeIJg& zGO~oq?C+gsry&qj1Zn=x&qJ>XHq90ax=SHz_t;)Vl3~iIC5@l#!$|7^t95ZX2F}my zM!K{1r*PUWC}!rGE%n{rI}H|hs6$h3lpX&H(hqe4XynrFAL;Kjm)Ac7^tlqyodO@K z@2II&Lm;Sy_gFm&Lmy%TkY76qU~6I4k3(~Ss<#}>pqX@)Jp$FZ@FXx7Tqu{G40tzl z?j#Nux0}$-1Wby5;HZO$MQkzo#@|7={YtFkcsWd$kRUhI9g)%NOhR=Tnhe}Wr%b;4 zHM{~qXS|X`)-sF)MSu-e?&6T;lI0FcMj6Tbaf;z~+p|r534aLx-mNnw=1-AH5nV_w z(aSrfkr4XB`YuVPf@4W0f>3cKLY6n76A*MCg2u4C}2T;-w{Z3a;=q0Hf?T+zvlyo2UEs~!1^8q zzFjQ#5SD4r-ZWJ(UVsZS0x0r|jkp68@1%7$;VMT2>NPrTr}U381E3m8?vqo?fTwwM zdhZU+DDy<MXrKDi52sh(~>(=rGE~7EDVfmtp>fsZ`9t^MiK(x;Tg+T;osJ(L~H`{X55-_kSqv7nIm}k=IAO*e)a)JzrlQK zx&=_uFJ8Rhk$S~C>JEfrdK`U$a`r$-y^;eS65s91miXNB+{DUsVKMoXxdcZ{*YDN1 zzLr4kl3H_s64YDQ7$iKW66AE_2d4x|$B&DC@FaTt)yL$m=_A8frV!T>0{8V+1-V5G z*mJQq10eaC=+o}YP3$tsBq-(sUtyYj9~{l7)baw}@937+kAv(z5UJD#zFN>3jweE;7au z%=~ECEHQzT*C;sP?i`P@#kj9Hz<(@`4+LGu``%NIw9GqhhvvhJg6FLfsMbP{gUFmt zvXaM}CxgGYx5YGN9kpPgiGj&cBs%KqF@J>BNBU^*2%hg4J~3Bu`6M8nhWha4H`9$y zbQY@19tvp&rM5_)uL_(U?uU`naa+Y4oDA_#7vLz$ND8}O1MU~Izb|fTh2y-%Z;+eC z1H}h&&#?ArOmt<*ls?ZQ_RlzQHi`{@1am$Xf03FV2}Isg!!ug02$W~1A}5iZM{tcS zBG!({Ixa1N__>6Tw`jgQ5g#0GY<)B8mX_{A8p?;zSL$cRG1$6hosm*VrW8txF!)YK z#Nbe4W9EJWg31nzZ3B(^YIS~G;L1d1k)0Ulb+jox!S-XrN+W}nSULHmZm z+$`XV6^)tzllE5+dgZ8@F>sM{Xz1*5Z2l$jJ`~sxA_8UM@G7&>?A}3%2&33Px70^^ z-t6z^HPDJ#apWio=Lc{_wg9HUT2lZM*g7bQybGdwcVyS?WxyTdan4Pd_1Pg3A##{6 zJBX-7|5@(b04~hqt&okKCr5MF_0rP!*jUFBh$q3!c+QTa8ozr_W_GIpTMxfMAAx8l z@aPGr4zYy(CVTeUU?PsBQ=h}D{2z~N-YW(MN_hcrHr#M9CtZuM(YPQYfR5a?Cp#T3 zgDe6V3S_ZnYT*{K#48DR_s5iU5d=`(gpG4cSE2}T`wYA~?z8o#r z_&}M98Stv!w<=-@?MezV*bm8fm1=!dUafBq%{n6QVwkft!Xq&qtTON}tEk7sj6xo; z6qQaVvS&(tzNB?IKfHmS1 zILba%sPDQiwmCdNzW64GDa=y{p*^}PH11FQQG$9mrDY{^ zhgwj~_1iM}t>+7$FQ>Qz{+tBOzuOPu!hV3O3@_9Hzzx_60|Es)j8I`2el}pwm{sT1 zG5`bZ(VrXP#t|P304C6$csBg1x4%(q8A}Xv!1y!pgLh}~ap2c;Fd|Qf;KmuHhiv6N zjvZTF!{8Ew>!h_WP(oFHmP%L7yLargw1zZ6zO6Z_~F`r-aen&?kaV+E?|F&HY3-RbdF|Yn2 zVXSWT>GvPvefGfCy1fUSQ8th@4)wcz{qIwl8W(Y=y|gF+(7fGJ_{Uv+3D~Wvob|E) zDAm##(OYUMhYWvo$xhGYY`dlnJ zGF~08;~DKLBi@=)$qV6i*&BD;g9TlBscwO`NI8NR^FNM@iAQrDx;yZ}2{T2q{DIGV zLurP7on_QL)LZTUP_qd9FAw$n6ZdMD-mU?!00Lgf|IbOtK_e6Q+wV_nzkh$q4caRH zS0hNm`(5Im#6=wk1giJ+>Cf#&fY+MM&ms8EMelkHLBJg^91!^rnwTTF&if1>|LCNa z^mvSkE~)^DJ~rLj`jE7C{TaU*w$cpOieJ=<#V06W&Zih}nB|j2kv=*!a<|$+Xh;Ir zmmCT_-4f3Xc-B^mSQ6MD)df3`q<4Mi=eH&i8s=xQxJc+5*dViu39E2_5=BTbS(T1x7Wv~Fqx8#Td_f!5B zJD@QPpNXyU{nf#D`8Li;E<8!XhZIHESLCoNQl2MA$I!%kJJ_|Rqgq-2#WtbVJujhW zou5Y-Criq%B>+9Cmfx>lQgxTVSe9Z4b4K$wM)?=BEsBzok`M0zA}PLoX7Mh;XFGxW z#>J)t1PTV`x>pf}C2>lKg{mOe8*8=q__38BfQ-(L20GP8*Z&8>lEZPb{o~raH}+$- zhWrw&`BLUSRb9T^UR(W;igQ;DD2HK?`vl!~8YPB6b#!r8Y8!1eI=BZ0qt2~nsAeB- zv=!r`nOeB%SOe$3sd|k*lsx83A_qLQod3K42{K_3_m#Lr3J%X1*kx5YK^e zsgQciv@O}CgQ%@!^1qk|I#d9q*(8KU5oG6EJ+%oJgqE)K}JSwvh0AeVS~ z37Ps~4J*e==3vp}5O2{w`LHJKT|ZD;(fD&d=TgUTwQ2Z75iWztQaWSL&Ed)r%%Vj; zO_GgNt~Vxo%`}ysGvz#9J!xP2+5!k5??%z*jXY7_P+1JGc=@X{4GImrV&7QxnNKy= z`^+S2)rRt)T<>o`xVO4fZS&xY5r_GAbkrM=egdBZL;Cu0`3KAnN454C+xVWLVSBWUuuOYhP-|I5;h{SEfcVK))&jSS03eAq1^5^Mj5N48&|nI=c*<@hjtLbe5^Tw=zI{A9ZfWXsi4cBHA&;O zB9mGo{}BWty&nd?`Frfb7l@KeEANoOsA;(w^Em}*q#o1g!i5yR&CR(isQ~vlyM|_K zyc?j^*I5gY)$_j0qO2E+l9aVkmo~N=)6E$-=+kU0R2bm z`K9yu+tD=k&08l2oxtlaTcdVB+|KenD9oWM^FKc$sof7WZDppzQ1r_DSF#3TjfLu& z*dKTA4rRS&A(Z{pf2kuYuZmWFM!6)PooM>wXA>A9mZ2#Q%?9Cmm#b+^Kfk7~b&*t% zD}+@y8cA|6KZ3ahPM!xH3kE+rTCot-Z3W^IkR&LShwG@cFw>bShB10z zR$%V{|Hu`%SmV2{oWOoOMkajdNzpBRGvjxWW$LwKFH{V zQ;8cS-kS0&9Bfm5i-oAr-7sa!U#Y&7)S~1D@L%@#Zd7a>P=+NHNoYL_W(wR;YmVUq zjd8W+HB720<)6ba8> zPrLCI-RfXeZ$P+$9Leb`dYqO#hOPpH$rLC%$=g^|kj*(FzT8rLHY;_YV5Y8b@dq!S zS9b<<9V)#1QHJZ0eY?>!Al)Mgf$@=7$@%y+%OvKLi~Udk?>NNqOI>>7A3?2$hL8+u ziTFsf`z8Nd+=NXumueplhP4x>B`7J-9WolXNg#P5PF_jueDW~JBh9w;Shd+w$u`qx zkZz)=@1xbZ_%+Sm=B)XcTTx5XCiJ+>7<@$3vKtnR+n7v(=`SlSLM_;q*E z1g_)yX)(p)1ZOc`VQ-8T&(puLg?Zz2VML@rFQYA0VFA+hU0n&-M+x6u)kf$?B}*FH zI3EufQc!!mae3K^oof;t(AkzjD{-qoj*d*CA1SO9@1UG2FH}=nv`U`^FHSt_F?M=& zQ)Y`5-N=+Hoqmsjp4XiE?M0fr5*2la|4&fBIsUONeP1n!Xh5toH2rg3LsEcMX6@}2 z{r{34PKJm4Wj%5VS-aiRu1c@%pGyc{{=CWvF`?=V-G61)naOriIL)GejW^p~tHpXz^f-BI!1Wq7K;K5<%C!usFh=@19aacgAj;FZFqL%NMAglQEIj zYW}(@%I3hQgZ7mmX!f4(TU({Wqb%cAFi*spOA`AVe;!|LCRzw+Xh|j-&`=&1d5tZe zzHk_~D2D$=YW|{a(;hX)vF?a1_WCS&DB4ZZ;ZQCxj0x4SG7>F3v_CLY*X%UMT(HqC z&R(88l~bn}@=Kuzwm)GToaIRBl?=CK5PsE3^`X%I%K53G^WJeVSm-n=gZJ8z#hW(ofxw ztnAcls&l;$8rSIi+RPLMh9{?te&O11a#Z!_ZJ1-;{an#pZu8#l5`aR=8w$-1Jg$xK zjJY6h0s#i5&aV7#3LYV7kHb|@7fZ(f;(L3Pf4u7O2V-p-wy=M`E(RM!42)fd=uJ#O zC>$4hO59l=@4l9;tVx=h@}!{j8^%i8e#MJv_}R~c#T1k4kLLB#2pUbD@mt1sJO&5X z-rRXA9@&ZgNPF^fi#Y0OurOIFCf9n{-V^<_LLnta9CT-gyGLKe1Q z=TF%^%ybNSe?&79M~JF7>{95&96XwBR@bS+SQ>zFrO@10h!x|t_%O^j3BL|iZ%|ii6%8pbJpC?nJoC7 zW{<98^!8tx)rV?nF-Bj?#`np4Kr65RTrrib(v7*~_-gY{67h+oOfp!^3XX00D-*Pg znz{{w4$MILws@>05pqMUs$u!8VWJ6v?^;pq>7Sjeg5E9+Sv%gxrtLWC$qq=09~}3%>&Fxwc_5(O~1EgS7#Hga55a3 zfH)bB1B+w|;+_nx>a$B_%p@)HvMS>wSHj?xc0Jj{$^t)a!`J}g$%h+Q?^m1g<{9(n zIYX78aTE?F@Rg8BLit3Aci7ewx9cCOSAK^Y6Hi8!>Dt$5jP!|R3YcZmD7V{$`&6q- zu7@lg@kw>#gygQ(;rPchFfdeBS6AWeDc65Glkj-*nVn}T zzehOx$fmVzy|+q=CkA5&phvLoqBR0v5%$)P&>zYz@~I z@25V8+4s_kY~>pUC-Wl7fo;zE*TksN_Q*{fG+;FJ2WB0DnrPXT7Y5Uy1y!%E#?X}M z|IwXYwn1OdZ1TD<@7#_cJzvgn!G;-sr*o){lTdER+tDZFNo?HNMv zOih9ps8F88vq~KzZz@Qn=(YhBN=Y! zI?MgOSPa=_T3sM36qrZOS&UZ>87LFJx$ z;*QsP$^?B~at{kuqLlxQO%a_DuPM3_ zB-&g5{8(hLWTR`WK&n}-iF)PiX2{Y{zig*YD^Z;JXTbwn4n6h!>DXw-uzT~Vp@g9ois$0$C3bRaV^+f!VdOEgXG`Pkq2rqdylrI!TSlGk z1C;C!<~P2pJiuRj3gGRw^9*xN;5#?X0WUX69p%*5n$7UM^!UlZoEQR@oZG6udc1D# zOKP-^71YZ4X1!Wk*LfPUvZX*eRTAT%E_>!_KKlG5OXu3+sIS~$+s3E)3ON+Neq*E) z?WCXz)ggVx^Sfge>HCOos^dEcBUEyJS+6W(6E$CJ-3!FumrlJcYbc~mT021L$7u4n z6FjdmG>bjtcrw{(Rh9Rq-kqUW@;}QUG2MXwba^7y>Uk61Y5qMHL{4J&QTqU;_Oj6& zSVHMkX$Og?m@U1GB2tfx8l#CaISav+_oyJ42 z>4TUi^WTIh4%)xnN9d=HJ$Y;K%0U`u{u(5puyzU-li~9><`l0zV3kcb&<7~K;LmPjs^ zubq&nM;Zs*28;{K-L6kg(u~ZWlB2RVno*{DK#xJ?7b}N)g;5j6RfQfv4u3a={O5>a72rA$W5Knwr?nW_r>5dh%yg1?71 z3BV_apX|TEKa$>xM&5dEkG=h@JnaBYD{psaH*aSL>pQ-7o?Z@auJ`%H_=RD2o_Kq^ zdr1ihxcp}TzniDM0M{y=H7JD0UB%c707$5=e>gx&8XW)#pjF}j==-N_&6>OGW0N^| zuZAh>HhEx?JJVAJfix}68$YreoEB4q!dp%k0Yob zOPY9CTG91cea8c{+Zh<2nZ9IlTa&JKa- z-cGk-Rro$uW?2>|e_4n1A|WAJU0=tn`Ak&{05dH?jErK=PkNp9nZhS&7bL~brewO2 zNBI7+z?(s=yaD2hSzi=k86Fj@4sJFNY%h1QCx=;mG_`fN6LTg?_V3w z;Wy%k!S1~8wQ_cbZ9cpU&3Y3P6O-I7f18)dJ%k(r$+f7)kt^i&Q>JGBvnFYl(h!~z zVPLa8{Iq4>@1P7XBWWM^0n}u*Q3*~t_T#Ak+VRfPoHY(U{*4gp zKaKH_51-iZ3al$WWSkL)AksU7wguLru;)jjW%hN)0Z0Ua@|BhPs238MD}pS)=NQ6F zp{c%gz6=}Cm@$KlkopO(_uNDin*`2Xx>8fuwHoctbYC>1gU*j$cFP^pkVpjnWFTH= z{ty!l6iHv;@#1l_q0TCG!&_{6_@geO-)QxG%4 z@cKcfUA0Jsbdn))vF{zxKFv0-@aBnAqw%OE>pQtzk)uzHAoPIXva>6myBVT89ovAj z3~1o25(<`L-W@r)i@@p;#Id8GL^N^KEneVDeP-F@ zT|_qDU|L8M!j{M&lPfnm?*mHuKHmZ^vrwCQqkq1~qM_Zq5>thGsrC$?9+nFAz4mZ@BYmM59-rJ+`3cIfY9$jW~pIQ%BdYdG=6*|UkRJ;cG zap~~i7>@RrVq@{@3NZGwWgTd{z@)ZS@sqB&S?dEw_TLZ8S<7fJvv2)szdO~38m}~o zKP3#rxwoVAGAuC(or}?m>XF;{84T>gY0YD>88h)-)aNsR0hyC6jXXej1)@S-^ z*ilcMqm744owWhMd?T`W&15k=zO^!gzwDD6a%)~|_6~7CFM!`bhC&6a7aZZ&m?2M! zOvRw_Ct6zMX3MP}I!M3o&9&yM&8Zy5QCO!Balq;N_olVEg%;$pL|KorNnsK|5AnbhgiDV1cF&nDA1pY&iRH`PVm@ z8AaSB>ntLRq9?D~h2)w}Q~RHqaC(N*kO^^hs`UZ33H*&l zlwwHxgi45y9N5Ftzln0sI2M`?wm&Gn)Cy)Pps8>&mMZDtrs7KX-ZTfqi`(Uy7Ty7Q zX#v4;LYE-3o>ZO?rh5+8@DR0c|O!k~tXr z^s3^g>nbd5{7=_!hTBz*4u{zW#$XE3rQ-OvXYi?{?uWf^sjn8f{VKE2TScSX+vLr- zYzBk(DPLoSSpt$TgouV66=!mKrASZ%PxtKcSr?iO?mV zk@S4qcZ@6WzY?vDJlI}X<4;2R$CTU2+gA!1nb5&H><3>~h2JxMLw%1*LjKQ2P}(yr z0r9%L`1^nE*VDoHNYOuqN>s|~-vDQKxSs_hGI{VdgALw=whw`ln;Z)b=1TMuMu!j! zrI%j@-DUpD4A<*b?2@H@OG_73nQ}HButcj)il!mxJu33n6{4jYc$rhj!&vXaSLGu3 zq>DT}cjYT*=7`&~B}O_Hmag4h_lbyibVqk06`i&T1y~Y9s{4o@&p$o@9U17=NIA9A zPW%u-LvF1j`06z<)7KaAaK%w!9@sXEi?cDfBU9|i3O!S^yNe?4&}z@kTh2Dm@8y!5 zi7iE5uX$3&OUQq#nI85>e`?KL?G^J(3-|64TNG{?Mt^fUnx8f7j}nRM}_pEl)Hf<|LyCYGdEJ&}%vy_%e+$?gAo*j;A%c8x9T*NKn)t{X4Y-_Rju6_Igeh_TE2POEoksrWn*G*R67I zGe4`abFeBuxawF+=h3b5Sav=WG%^`tf$a}_oos(qKiLcpkXKO($$Y-ykY^SNEj8IP zzSZTVQH`9WUvAeq8PFUky_s;n1dWPoZ($)lIz$dCl-YkbgK=3jxs#2%@(HSsJa9J7 zdr{9Qy9F~M0n(YnE{yZcBNzLk(`18IIM8Nm#z*#J5e`P`4B@$k@)3YUT@Y39-sbR6 zW;#eUsk)w6gu?`Nkx1BfD1v^IewA!%0##6b%XR)mK8yJ$sAv(e%}2AlzCzy=yu#6F z#_$Dg_n^X-SvL@omdV{#+gcf!{g_Xq!{jXsjJKM@*fxaA9Ao$1Lb56xXSiL&IqF1{ z#LDR{S9)L1Bww7*&()VsMm(_m$MTJqZGXNRTTnTEK_=fxqj>Wa1RQqW#RpxTujhp z3bBDo7C-z?qhlIy*( z&$u_047q-y9Udg>zw$lyw)mv_OKaYJ=U}v#uxt)LSt|A8PKw6wA#2&o5oxIN>kra` zAIuEeaPWSdb5NK3E`VpH&{)9PxJw6g#L5C@yrw3+Tin!8Os3^sybne;N8rZ&*MEpM z<%P5hqrqcES99IRJ&$Iyh%CVIwepn~?OkJ6Bt;8Z6YMw1d3-uBiU8!7hD_ii_m0#0 zXwAKZC7W|l;PS0bDJ<$`^h`wIIxy2t1yy*7r~4TP(?r)AyBI+)onTPH5t4=GA!nBF zg!K2RLTEF7v>!QvN%#tE#iNRqK;*zE{TIR}1yXjl0M&t3Qsm~e$lAvWs zp$U7@dGDVcFXxWnYaAIXrI{z^hldBCxBy>s<9Cy4PUgl0y*KZBanI_(CzL_3>8ULs zDeFP?n!AKIwL!$$%&{Y|KiAfd1xg&wZCyT(;Z1$9h+j;$!Ny;0_58eRw2~YN+PhcO zO~}ihvQlTtQteU4j=WMkdcwTSp;eVO!g$o;pxc7>KcU&QrR;&<$fS;g3aI!I0vZ;Ye~2|W2e#3G9T_qP=OwH zr-Cu;_`#>YM&SpRT^vRRZ}0`(WW~slehbQv2nTZSW!;+sk8+^TQ%;QAOYj{9tkJ4k z3}*(Xe9Q?uUF*h{+K+BgQdXW0dd8b{P}G?D){8rmi-_PXaCb<&7HU&s@OyfXob9%t zva)ioHBQK3;!&Jtrrqs~&>y|8$QUVn>l@Ye{J!C`pMld<`(j!hsJTM?F3*MLnu7@i z6m3Uy%=qmcKjQ-M6{F3bisVr)zATIU3KW!zVjIj>J5((`-rrd8yovfsGR#qy{er}FZhQmkXEP|3e(!%Le<<#@P zw^2le?xIi!MHu?@Pvi<(@H`!7J2}=)NwjZn#x1SeceAK>F>H*rl1UXrXp9sQlCrIr zWw8ILEDL+yqc{hBVhLO{YG`TUNVCq@>Eno{!!KkCHjDqrVM1aTGcw!`DTgNrdt`)v zMbZPHUK0aRTrrC?&%YukA6r}`6puWFJNSUN@BG|U`+dxtpsnWQ9&Jvma%VP!^}A4$ zuNH(}+yrYz1ifO`xhznq8;Ik2P2THjcjx`De4u!%^))CM)s%eMTRvQ2V`G!IX>OQ% z0Y~u4W=c#t8-y2-rGknJUD`pT7bpSw69Mg?twBp}Q9YBX z+P~ug{5K}$@Tr^+U+;&s+A)kR`0cF3qgw~)GM%hW`Ww0($Vmb)%SdafOj83Raefp9 zB!F&nqrFrhxOB53PW}%pCws3O90UT&3L~Srxj!$JuNBqP-7_v`gXs3$|6TI@Ono`; zt)Nj&XprFnh^!t+Q9`p;89^y6&QFD5{5?Z2VTnph((En)rwRghAk&F&y!6xFi? zt<-}bNmxoNF7MH>m~Q&tzAN!&wx>l6n#GePocKlCyeJmc^ByM}Sj7=sc{b~lF-N+v zqnN|k#uPG)IW6Op9->1+aG_aU$%(3c`jqctgFw&vQC!OL78?573AnnofWq+WPu7yO zfw%huMnEJWx}qI=0}k(>1TWHf*$hmx*_4NiUvrk|!o)<6!2o`%a6L}{(6(BE2*5Oa zz;Nmx;zJ}y&aIRzL3hqMU#~kw130mm-JZjNKsvu8Ey-`m?!z5Sb)8W!O7f6}7*56u z{Qh%i>8hfc^z<;hzl)n7i}~?9!_jo!u4HBVkvKz5Q4z^CL}p=eUE2`GSBH-4Sxb38Nf1-B72 z3q<5yuSCB`WLu@$mI+pd@mSrFjg;GXc!HDR2iWkGormxh(qBT%*`cW6U)35KbMdEoTawzU4Sss$`1Z_T6!% zHceDoB%yV8Cx_@tYywpJ9$`x1D~}-oM)2U9$OnikJOFNM-x~}2Tx#qmd?W^eRCmvX z&ec$J7UX?IuEBnYE;?^h#q&e{h=yFap9;V0vv7EU?-hv(`%$`hV^l;&RT1&rM)n|^ zBaWTs_pq0LDsK!0aka3rRV>^lcvgJ z7F9sgaGzvDrqV;CWn`j%{zIjIQx=zvlX#_$z1xUGrVAcuV~RvZkCR<#8Va>;@A!qEChF(V2+GC2+_`Gdb8d%7 zk3gXCb59}_@tity->9jvnDb+fxFU!>5M8iN4WN1yA(Q!-{lw?RaxfodICaa1PtRg% zV1xiE#xVj@>J?O?cr)=CJcm)Cl_(*W8+gcfWsazcL={|w4g3~mfN>&5KdG+h5(wNh zF5Xm5SN2d$KrhOuL$h|1Ar%UO!zT=y9j7(^injTfmf3|3F$`my`9_8GZW3YO*^!Iw zr^z8-{`dEDV;w&SfF#^cCtrs7)zt#?ewU!!3q~^Za{;Kbj&t*Dr)NFbL)AHZ*&QY^El<{%ATk}uKnQ(V`i#39LI~phL3xz&p^&?qm~6B6D+dkfG@WbNPqz@*wa^ z*Oq-W)zX?HAul0krn@AZ6t{|RdSvcL1wbL-|P<1NSU4~-@p{N@FgFsE0J7xfb5 zJJwUOtq~{`x}Jo>mZGotMan5GJ>(A?CvbVI^W=@PN2pzGFm-R=As|tT5vcl&%!aNC ztAlDw0HMqp7)`bY21tDm%AC%2N{YBb@4p}F%b$5^ts`=1gw-X@us^4Hr9~WzVuuo6 zM*A9lax@(-CyrVWvZ?lwQ!XyCa!91LC(o33TPWn$JMh<47BO-RM_^Ib9JFO(4c zg`b4qrv3l%2^TUFLxS#|u#!=I`t@skkXxl%gg_wde845M^B}2G!zfjp^2;*#A->8=M$p5N? zKLRVm^3{hgmCw)qI-%Sst@L}1zrfTx2fE~!>0vIT3GwK}#-c^KsFyce?8CBAGO>&L zN!uKj(!W_M=3vsIpgQ||!5|uEKWyK{M8(Yl-d;SP>7P5g6Q!B}h2%=pLT z{OdJF`%gDrK%QeRw@}*1TdAEm$!%omQ`o1PMpDoouCvcxBE#R0;{B znOC@0zhkrZWgzS0E)+B^aHhWob{yF=FLjlnuJ_?a)c#Yc`<=^0xc)(l9E5+s_^y-m z^5#qm8F~XaSAuWjN`N;vO->zzaw*TYdQ63S>{TT}r4_vV41+I>IWzcXI%rB2z1dXq zC?BNaIx=T@^81j@EnKuGVxJGhFPwc3Zn$bz!)6@}5kcyR=Q$W-ggGHHYc2aVPGi9- ztK12y5WKTAV46?hb}?gz3=z6UyNaH zFp#3P=9eYMk0cSY2VesrGJpR(ZLinXZgY5|=&0VpCIgo`{bBq5^Y zrH;Kj#5rQcnuvMa@AKcrovRH-UsZ3>uk!WR--yC_NAiFQGTMEBXb>)6?NOP|dtFs6 zR*X8=fp16kTpg&iDoGjD^!^V@)nKE5=x%Mce2WGIg1npxxQ#ySp;xFwIDs_0U0Gls z`r97z#|y*+uf}BRFvUryV3H~{fDM(=hm|8c0|YL^N52f>UlNDM%K z+X4FypdzP#wBBN-=e~U|S~yX@gF>Me@5swWyfbiU^gTq375%`GbHl0oY4Pb1;)=aL zloH8;&Y`-=-}>b+?p@wMj%Oj?$auBW5Oa%ZZ)gJ*Xx$9-Tl=vu$vpYbZWmY_hfiuc=r zuMe#wQQwxqSa6gn$MbX(QuGmWyO0FERY103_itZPObaS#`+RvbDpB$1#*Cq{lZ%n= zsM;W>mRRa&WF&GdcK|fFwmm;}=5gEw>1pwwAeX@0v%?9E!sxxX-(5^KRA)m<3x{Qb zNAiT*2=cYaact+$NHr|nHZB?g9>5FVfe!YX>=~vKSC3Ql!JtZhJXqI{0O`Zc)6Z-t z`B(ezhitA#I=lr(`YFwKHiHElJEoyB0#h^&J&2XH*zabIfQa6dh6`sQ6p`GMZy95R zuMSzGjlU048nSbI&klA;{XT-zSb0_8(}|w_DU{j5)JjNQ3v6Ft+vem8>Z#8W-tfPoNx8 z`-MXXR|rJbQl*$BV7r0UFU-0XqgBt6@s~Tny$s?}rm7(*IOs97$!xfjWxLPNw-3@mHsXpeHjKM&6rEA%&iEhdb_MX;v6G@9?dx)I zYR=%G1Xa4eC3dkOw1YdSj;7~M9}{BDdARc5QRpt`e2g3FLfd_rSepo=A-p*9B^ zc0sr09_&d|)-0~)sL7BHF>R830`kh&b|!PU%5d&`{;&z0@DZSyj8KL1A(fhfeOm#) zv`vbIY0Jrk201mxaJ#;1L~`~pv%7J8C*jaa7)hyne?WTxjsNXltiKDv8O)}qhti)8}G{S3{-iV zXLT6_D5C;;XLDcCAJ`nryE_u;HXzjMOV=qj)veYVlmv}9P^w?2p3^PFKVN4EAf4*r77q>#Z5RI7r1Y&}{EKb* zJvIKPYAQ-m7jHG{#bVoq+h-h%KXK&!MMpOlSjluh`yXl(hg_YDZ%AG>sFI?IuohIq zu(^i19RHQyq*|8sU%g3narM(E-$Dq!(ICq;BCjpPel+-hsy#>n(`0yses2QBrHgpn zK1}b6+fB)~00X#32lw5Q&H%LjA`AD1DC@gnbyQXy!vAwS+D>h>vWFL(_CuFj#ZygX zf1vV0E4(qo9SxjO;?n%O*^MOQLH(2(QXvI-DU*fY7T54CNcv+TDkmJJ7VHTYwYvgguukpD;$Ea1<327UyV}rhp&hGT>Ud^WAB5n zf+uka2(FZ}UjChL(&WTp@~b7%%#7y_TUZ!RU*;?Ns*X+d1F-l8B7=_Y!KyPGKuhxD za-kTK6YtJDNCA)3L{zkh7IQ_pnDf)5FoHyWrZmq7s*id@(npJgB~Hmr1&`NcYoKgV z)mA~}cJygcdRkKV{awpxO<6Y`K!C9O^0e~ZH|5hrW6dW$zC|BaBj*}hj$R?@0$JVy41q$W)TN)J^eVZp*Bv4+h8iP6J zh@p=?0Bh6K@i@T4fnQ=kH^OHopp;SJRvUFlFqp(itsP2{qJJ3p=C>LtwA;%cBZ)0~ z3K+^#g~|?EvyvRQZ07^k?>YScS|S*%mA-E?uXcht#mREnvB&RpMqYzx0WNnp)jeY6 zk++1y;-x^1+}Mg_*ywpGp9)41W`fo#`pt*f@{sj~v?&7Fz;8^$v;+tq{d0(gD*?C*S zK{YRlNk^4eY8SDO10$rVbiXeGLhZkA8U)i6SX4YQo`UDFDzLPHS%*z*!Mf}CeguxQ z9zwb;tE^00+<8K0YGB0^OtiP;pE%YB5;)P3X6)&4=dLGO^#E_86EGXa+Z5+Z&u70K zr+kZ!jkWgnuBu8O#O5OoJz6i87+?ap5{p$5nRvTlm|G1Q$7E=Y4h};;F9{cV94YB< z!!;b~-6T%#O^*r>bM(oM4`Pqq_!Mw9?V=JLDJDgVQOf|S!Q@sFebv%weaFDkQMq>oZe9r4NcQBlKbHW@?{@RE=%y!lItSR)YnZP z43s=lOT0QBQHysm9#aSVBKjI}?ZiWH5VE{bD(KLT8 z$)B%e>{TQh!%af!;$jfFBv*5b3{vyIB< zHFLnmZKI2p>9dBx?E7pw!14Pi<*8&G;pDz39&w9q!$?&tT4G+)*GG~g&oTOv!?Z%*6Ig`A$kG%Y#9@8%K>m%|~y#7_#6<@8;C!&$wV|yQ{0~vm7-PN<29+kyJ=Nhj5U;CcL-~F}YTHYHo0GmZg1rd@VH(R|_i6gvA_?dLr#R6?>$Z(xL!M8@5)g9N-itA0DfmoDRTlf?Vum8b z6(KuaYU7j`IQZkxLJP;J2os;#G~PCEGrJQo0p9txn;{#37@owIA9LiMXnuy5^JYpm z=A3w!q|M6dj4p~+yDUfmceS=g_Bdw9FdgeZak%3OcmV{Tc`*UIVUq^e_*~_?q{$H;>Y4hN+zc-noCTEA)AwB;6vIvjD&mi2^tCuT z&RyvBpl|mlY8i!^evJggl39Au!NSGap#hLJzej`{&-i)%o$%EDIWhT|=jF{ZZ7(|$ zmxyIntZs<|)SiF$AWl$-!LGLSMD>KA4oqoFxr(qs#T&A5DA%~|?W`nBRp(*K@YY@u zR#MEXAB;vxnC*vqu3J)3^}pHV8hx`*UNz3*l-C3)k2CnBv5kxxE|`Uk;rcgk5l-d= z4*VeEm>pw&x;*L;YYg~-siFN9Sa@UkZY7pGFzDESrqm;0P!k*>!-BW?w7(F`#Mtsn z(1@U~{(~8O&$VQ}UT^DYQg)s7k+x4p_Ou7}u;0FmOcT0bpD99w058zlC(TA2O&H!M zO`9Ez#j=%T99KjR`aDT%&4oZDVbEEvS2mQq;aQMICoMees=>Y}h~E?j%=^R;%!pAB z#BB6p`Z@_q=A)jpj^*z&@W$hHn9uY2(sN z`6bF*nlG4&`t)bIZ}i5sNs63}2VI}Rn1sQci#<*nPLx|xz;GWI--c!%^fVmqD5q@y z^k%2Q*PKtYF^ey`zpV$bcBhZ^6nGs=6D?#ZJ-=n$oH^?(8% zd0)pwt1!hCR(mcZiR1IdM0k_mfS_fGt2D_v?Bq?gS_i6ob!LN$z%OgU2Ha0+<)s>@ z0t=3NU6%vf?PCv>pMhBv6*3IG20XS;>rur-+WGab{ep{QVv#)o%QsH zQ;o>QuxzVH=og3d6&+HMve!LAK&vO&Dm{MrG|^k9l2FMSB5LJ^2FvF_XIrfSywCeh zBM|!<)rc8WsFgwOa z_Xlw%#x&}sI9FTS-ufj_YjjYuPwzRh?@igthYa66x5%miO+ZbMv%f+5g=O<38=pEg z$VI`w9vAW8g{qkar96zOQBl0teZb33b_R=l)JETQpY_OtqVx8@!YU~6|7>jk-)6pd ag~K--MtHyM^GA8(9(72ZjkP7X^@a^knWVOck{Ww&+o5y zjKffmn;q9)*IIMUIaip1ocMcWLS!f?sP~c*VoFd@(8Iuw=sN`98%6Bv3E&T2XOM=o zvYnZ;o1vpAl&qn%y|taQwS^J6tEr=tg`F)sGZ!-(Bl&k{XL~0;78aZT?*L{yM{^e1 zeS9O}BuMrWnodwqC|GYl&`@a^_)t)?>ylz$RNOO8I((douYJy5AIeRLrwFCN-LY1h z--7(0$-*h&f|}_XGq`$~krt>O`u>jO;OB2!3Q1kfT%9-~&c6OVJ|6ik z-)CY^Gm@1yJdZXUB3Hgblk@i1#PIL};EOyH5s;E1Hb{{zhoa=|djf?x72uxgLY3^oMN1K=``ypDCJT9 zwGxy&82*d1{GV6;F2TqSdy}UB=;j z@fggQ!Npor>7RpPgIl?$!j}8j_~KO{dQU;6KS!)`bd~){OxhV+cL=Prkhit4%#ukm zDX7&36+mnJ)ta@7T-^0!LHh^G)Xio5dt50yA6Q2vK`|*x|7O-*u=y9Uxp0YJ`0S|G zFmw3E4=1;tVkRxmS#OKIEjwcbnfFLShh=KSS=@YAAS*nnS&F3$T%=)W+D zh|qedX`5tA{?EsYuJ{@{n^*t+f2c%@b`7$>NP2KKU ztG|-+G_YR1WNeXo!-BounYEKuHuL)MgA~7hyj-kk^t#40{%vzpT4z&siI#&nn&!-& z>3pPA;zG>&s~}?5fW>>uI94F_AvGIJq$bGOX=5a4{j77y)wN%zxBB&YZaE4`;J#KGrTt_LcA;glaXpVS#5)I+L z!f>B_hpaE?!#AZ9#2kr|$ZcIKtPr^fjDLVAKT@2c&bw}KoyjH8_hEF43I__OL9+|T zT58^FS{fiJhrbJZeL14;S%Kdgy{i1tF;Qs8eRlZIwCkk(F}CPWShA? z6>KIPF~7Kmq0%6OV;VUWTEG;f>ss!xpy$(vzZJUcNX|_0tf97JSrJzb-OfBEYS-P3 zQj(-fEiHSf1lrG&*UTmogr@kscj?>{?BJhDOMS4$oNiTL0;^^-#zdW3c5e2N`f+vh z8f&JH+i-Y(=1|2lYbcb`PxtiiGH4LD@9}xfY3rY@ZmY9}WZM)2{H8`&*OCp=nCUC~ zb;RSz309kg4^soa&b6vr#FyL29kg+0+QrN?IHsXDJB*;jz5>`8==TxNnIVI-?EBD1It<>WJC;&`%!BKcigX znPHVuA%?)#<}mi%lll}0xhobX6g59hykEyROVyi+rC&fwbm*x%-l^q%GOUtR8q~Ey zR}FmPF9tTo)!>k>r4#+iy#KC-$(+gCgin+t;!`=+@Cn~Fo&SO@IqKMznREoG=ZSc3`Rbh!Z^B<@ikXCzwly(o+88Et?PI%66mIZ|;`MFPG z#(p>-HCTDUala87Fj6zXIrH5WzDALkG_S_%e8E$LQ9y=?IF(@%b!eO~*Q64uhZ94} z>hpAeVffsR{$W(}?y zE^&YElYS1d@IIqMeLfYFtMXdo9mzEwWbHmsrgoA&vi6I7T)Y}rglpdLQ6_j;MPjx} z6>5AoWtz6IvSnO{e;0nEs#dWhCQHU)rAGG=SR%ta&Lvz)pTnvEn3qofkEF{j!l0a$ z&Fc|<@x=su6G{t-R>4t?_4QP}4uLd$F#~Ji#YQ{>B_GsCLrROERS`AiHe7fQ8QG7V zm_})%g7J29<$2it;RS50;0v;@Z_8v)%R@e2OzU9uJ@#mc^w$PMM<2Q47I;&T^g1$# zIq4pMH+DeFEzy|zVSSP$U96e7ZC}~rFA-`>`Sw-pV^_H*es&k{4(eEv>SfU`rrS90 z9sX%=0qAM-(4HJb;dLYZ_{TuhAQT3UGp8zxPa+D)f>11V1J zXY&~nz*VA6EQ!~^u#6r#jKv8owu-;EamVnghIj%JBih&_Q?|GGZ;PLFTEZPGg|ngJ zY5JS3hT;EwkYjnh2t&>9m-};92njH2qQ8$TX>l*=>BA9O@BR>Twh{v^m##C>JJHqu zoy|>}Asz~R2g}auGh(uyVlo*JcLrWcBWfvfxd{prT?+}4F*@8k6=a2MSJK?Za%`dF zU$g#QD`VMSVJ|8mZrki*pMH{!(x>OmqbSjE---s0^%ss4%IJJ_XX=j%`||BuopfHd zsVYsi=NA-ZxlHU~9$jrGW4XN)SsC`Q@2H_+1%lR*l<(byT0;0FT`;`1vP`?ca_}5fG7txhzz7RD} zd`02g`T(@m6dd@CcU?F`}_^~%jqGmyqrs+no) zU=y8c$c^oFPEtVx@^RInOR*V)!_O)#wR-LnaEHhf~*N2){v$QenI`+LdU&`89WR!Y`ObE*?C4kA-g= zml@vXZ;g&mI}jAp@u+RN>8&Lz;*%t2Y<*KTjN}40l;ukXhJM~u9G5w{rQakyFYB6F z8kJDSp|C*NCIR5cjU!q05`LYsn5Q4-KnnP=gVdmf)oh_7CJlSIzVcJ|>-WG}p=;xh zTe46WVAfWLl@gqb!k<9_+8nPydAB#>(4{pLAXi`=76$t&rw(}1#pP&Wy3Ik4>W)T4 zviANT>`LlyX$jg>Prc*(lBAE4PPH`lvcKVk_z(a7 zpArwwBw?D@`cy`}(LO8J+gXun#=0`!)2F*H;wa~J*+LSN$LSVZ=n@09Ka-$LFM}tT zk`}6HzP(>xt3juA7&zYa0*WKGtx(SUxlghtD-{$|NY2nr(96|94#1K`g);U!HZZ%| z_%hgsZKb((zxL8_{NMrazyOT z;}hF;QANrngrmtCHOX7H&XB?H^whSst!*&r@!lGUkW86-)W44v-xjb{A&AW7S? z?CIx$`j*C9oZ-?KNi(TQx{uL2`}5ol%+xw=VM#W$!CRTTwJmnSpYkR3@jA5`<$?6W z(DlSg)^YxON-*J`>aOqb?+^2tyzGrYYqz7%DsCreD(Nx5<9z_}2RqbWd5+i<><~gMq`jiJS}K*|4!TO4 zTeM)6VG(w}h7$I^GPpZkNI*EU(|WDSE=Ll=pkSl0Ft4%2io=!vHE#LX5j5e|l(=U- z+P%D7@psU6x9Z(&cn0n1Paow($BockF{1;MVcXMXvPzA{s*wt6O)PJE$WR}pA!F|~ z+Q#Q3{)R#Py0Po5X&D2;*q1l5Km|3pFt;LN&|W ziWJo6g)tNAO{GdT&2=UnxJCR|OC?>5l4P1)aAcFfdAcp^J8f56vx$tyFFN&|C^T`B zwhdUG9LWh$#HL7Hi*4a~4X)$HFr?z3NH6Teq>#D|HZB$rfz zB|LT6uhhF6j=u6HWQ$lDMJu7mV7D=Ph-dfcX>|3mV{7r>!3DQ!=STC?4c-U7bJ!^3 zf7!>Cuz8m=gS;pu*>1al?rT6@ET2kj3z=A6HvAW4Zi<=J3U8KKg=;D?lR-r=d_0%c zR`CsbK;+qLqpAzf^Tx91F5j$@zCnsPOz70`kAFmA6_rKZJ~PQ5Rr0je!;g-qz;B$0 zk%ZWr&sjVvD%P^cY5X`0caiGS^9{ei(B(VkeaAsmZm+czyKzF`OcFhUkI zK!!u*MqIOUcPaD@a&Q*HoRs>_+0;vZI5fOKwm{7BQGH|4OL?BtcQTFq!)Gu|$*C$R zHKo%(H^K0asqCes({MN^zZ(CoOepoKSX>KAj_u{$dlr{Tf~T=cm=a*~k>F<14Vwd* z#@gE>Hh3>{T-H#QTvGa!b?Yx7&PafQbnO(<9~VeTU^H?;9`>M^!ziEs5slbk!!GB)Z z=l)>ceSh$o?ua8$RqQfayuxLCL?mpZV zaK=A77!qIj_i3f*zFXNAJU-l%4-3Wsl`EeG9bxZ}#i2wKDWcb3BC^13NvnH2pU>xC zzgIBj|3`Ez@T8U-ztg#H=wA6Ghfp?eSG_B*twN4d7*x;6eJP)vGV$qRRMX>tSgr)C z4?F?&NwPOn{&<(Ht1TT^`VIMwlano&gX%{`2C}b`02`4gCXrh111qXCe;uDdNz#Kk zv@FcU5^D_f|M-g#ufyIUtaB|08mCa35^!)hx7mlU2G2E_kII+`Z;womSUi3ia4WmR z%-=G2Z$%guT~`%`0tifdUY zqZ8GkGqle@iYl^7;x(W(4K+RleG|ElU)<1 zkk3I1q<7!1_|o+(K1doV&Gj<}xu0Fx)#^X2fNV0NMlw)Ky(&Yq+R8{4*iiz5G$&E+ zxBu{G4cosQiz8lDu{5D^_uI<^&92k@6iSK!UDDdz>9F9FsD{l6=8@sip23;Z+>t0* zZRdJv0>CDyQCk&<_JJs)6BFtdUVH;)WYnw0>Wb@&Q-XDcu49zPQBCHFT#xhjW?|`$ zaL3XM-xg!%ZL6ISu^8_ZI@#}BeNN+et64{_;9oyT9s>k_a{KS8w*lL?LMt~?MrH~> z9#Z~^BEp~4Sf;5X;{%4<4=a9q?J#+Elr;#tKZB6^KfO$AwWjHIX6o|$^F`2^wqnA_ z;PTw~s-FR+#oZ>B+EGjQO4U_mW4Q$PFh-ek0rH;ew_ylzuc;(+(^eeRsJ|&i?QgJe z1FKoH&R)uZ#&C$=Ob-MRrX}~OR`YDyQvy`NtER^n5@p4_H%(f`X7O(Bd8re5w@N}6 zuMFL}T%QV~q%V`4*6M6d237$@?6DSLk3#fCF`%jC!{O%(X=>I+%2FJEsP2&{GAMa9 zdD42gT1jXcq>Yu*N^w0sbRj486CB=(iLvXKRk~SG+Shg)tb~wF%H}VY1AJ`n>Z^;f zQ7lH=>(M}vsx*mIp!1oV9?n1SC=(erw)THji^(VjKJQDMAfQ(Fwrx)bTPyu>i8JPm zJJl$d6Qd7@z(yZkcg6rQWFc@q6z${v@Yk;>Os8#p#5ED0T;VkfxWbIi=J=m{ae!?t zq=R;7JwQ~sxlnpFtO!%VE~L%CtZ6&GWHnE@x`>w~U2hRiIg9VV?eeW)_3XG30XSc!V3ZSlN3;C*`6BHndxu89 z0OkXBu_=$XVYH#5uEYtRE*|0KB5bw%MUA4HjzTY{(7)I=Ol+a7pIJw#E4ctrYVsiWKL&PVjT+W|31#uKp5+&x!5rZotF{$J+lhn-N(rng700>YJ#d5nw8p( z)*ne@0?nGV5-8g^yAirx&sS%n@mUR~&4qtfssTiA55KNNwIlgy$HOoCT0MpWMzK0= zY01aZ?%GB)68@~(6^~BWd0@bT60Ib@~245Th}KijUr1Y zsIJ1uF%@WL>(D?{V#ic8e(O8-8TVvTqRqs0KCO*8T1uY@2SSVT&a@+%@_h_q+z{-G zB!3&Xs*a{a5756lflP8kb&OxNX(^mw1)BC_HS7f8TP~s_hq3HXS|3vO+a*AXQsMEx zw;E}?rrTIjhCj-ZBt>urk68W4QkPK@V@8;RiW2F1J_~)mS?m@UD}cs_;UYAC;jdG)saU8n{Xn>ro`DTj)?!aE5)sd32xFnZlu_LNPOEY{*`?$BJFOSo>&3Semx!wbZaE*^6@R9A5RiW^M5dcV5X)H z^IU?_ZgDX0N={w9ERWW-zXreSuFsY@a7uC;>G+YZW=+Z)jseMWBe*s>Sgb;>b{gMF zI5KQgo1%6Iw^z2C&{(GH85|SJ$UT6Xy!|_k4<`NQVueWSCYM%f=kL41w#Dv&KJxkc zakOu=9Y`^4WWl=(BYaGhY*gfpIJ~&)_5Cc+&PxIR(^|PC;JHN2ruAFgER{wRLi?n$0f8+)T_ zKkzi;qC8Fu?gGXBS+%5D;V^OcVm)`qr%8BE^YV_If~lKhLy-U(vu4n;%gFYO{xJO6)tHnTHi(XBzH!{l_^h)APBbvF z&*_vvIsDdJ_j-1w=URP5=lI}ZYIRlC&8^967Paf7wx6>9B{rKx(C~*9(mo6>AawLw zY6~+ZqmO0q#&=7Jpgjxt=wfgR0br|qS?H7%*I?g`26-(zk5&pFsN|08gL}o47m)~w zebY10+yYphlA-e=m~o=uWm9ZPcj5}=5z-6li;FKe(GM#C@BJik|e5f6ue1>CN zeo!K(7_A=R&zOqAEr%nmqFRLd7v#yd_AoT=6~44Ofle%~X=`OWQ$PYzq)!4!mA|dh zX#x}BxDMY;}M@%SIf59o+>PIhKTDCkM;(4KlgBJstW z3){9o#$_C>78i>Q#V;qVk^5lBTP-7&6edBep$j3DbSn;}a4&iNTp2At6XtS77oMlW zW-}JGa+g~-kn~>d>_9|f)S9#UD9+#}^V7W+WDPkOVi^Q++g4Ek3B&0O+0%(v@{XjKS zerP4!h*)3M^|l6jk|Y5rC3B{EKme>@w9O_>V%9xO2)I+T5u>NpRJ4ATI=qyA)tWhN zf7Ti`t7QNbPDyR2PZ7IgQkuzqBMU5PH>&5Uul0_?Kxw&rlLLZ3FQAAcD$h&)t(^m; z249bdhC@G7kYZk68hUo7k?hxpUf&*`7BWOVru(1|IKyAh(=&Vy%u&`>0Uq2*>GADh z4S|w*KLy0$yB#OHcA|U7SZ`1*Cf`XyHqaNp2#TOmn!zVW&^lftcy{<;0537GdqjF~ z^1!hzI#Ns59VS<1hukA-gl!QgmWbCtmtQrC8ZPB8 z{8_&ZAepS`-#}3=eOKjkTU5;A(oseFb=^pbC``QbuZS8ToE_~A z?ny*v0SzAV+Nv@CEqh*S1rU)K#a^2b{@eG6$vaCk(8Y2gWG_yVQeTj9hU-Av92oar z&-aMupUr0y}KQrJDX<5D1bo}^DA27W(_$7SsdsjddA=)SzqU0 z+-6?f+E{m~jeozQT*3tw?IY~r`5+%s!y3WIRy}Abf-lqsQp?wXh)JN`wDEj7zp{hh za(wU3HYE0b97T2|j8P6rq6a`ir4|10(H+n$hiLmiHfGS7W6Nw!`D`=+<-99{o9X1x zKX(2wT>#d`eHlBL`CX;mBuA8i_|19A4a}14JQMo{ITOBMMOAH_9R_v&wXyTx&An; z4P_3k6z3L`MLHXA&b#&H*^LNIR{=5?3w6~ z2X-Ss0OX)TF51N(hxt_tPlaG;1-<~j811<{M>lwOqg`L@rSI%L#(f*ittK|b6gdTe zRoEUmiJQ)FvI7?=k6vpe=N$dr#YgsuC4pEB8`Ubpa&lQCkVRNHZoj{rabaK6e**i^ zs!(eCjZv`CQ-tWgMnhS_3Uqldz*%z*0G}1~a-6sy-7!;dW#xUuX{}zwxxdvtP#$Jw z{@WY;f~L9#4I7!a8pH7h9eNH<3SM%#3WyF2lV*{U?D8kzN~(U`UR>iy{6X<;5phMT)PD_8b5EVRb3WS7unc8OX{3O{_sf z8N;SGa$)8_o`cVGL68lk?@4cQK;Tu1Tt2A=SStcTJ(=8dlT5xvXo}$y8fP)scLh&D z5nq(^N6M4z{a3f`$lf@3x%JMKchz!ihukrFv?F5fSlF=A8UQV5=Zy-Xsd$PN=+y@} zfHmAgfCOMrp%kW|Xfo8a631`aB=iS~{H4ndj-RgH1B^Bh2})#PDL)s;+5q(bg7Z^7 zTt6OZ{adsMNQ<;1KthQxmn#+Gj^SQXo$<_JUv%AnnJox4rXLh&C--9+WEMmZD`6Kd zAu)aVqkXApfHWkAO%ExHVimsj5*2%DJr9x%aF(ux-bZ?p@a2F->xgUBAwRTLeft#~ z;*H4LVlbY%29iW0;~K?WG774>ayXK_3u^XN1YUOy;Bj@xj1!V?i;hPvxkpr*; zf*Jc2bcu=4@~2HnOij6yZ+NGhtu0lAcyx&b?TPj@chn`mVh2Xi?Jf63>iZsA1s=+u zB=|A5mhFGRTBI$CTSLZJG>i1m7L)5O6(>fRgLgFZFhbdU&LyMm&K>>l?$SI2e!t8> zM>GEZ_JzVCQWLQpD-J|aOB@Qzt?Ie1_AaN#!_R%SY-v;#*wMfE2^l`CTgB~=+5WP$ zM~mcSCquyHhoPfDiC8pWz+T<(uEv|who4w8@4)?|umK1JGxJf>|5X^7D}5l}Ad}a( zqjX)yz@+}{@sLZ}hLatI-^Yw2eSBr`lfD#~*sr#Nja)j}0$R^ec%dEuEY+T`OyB#j z)>@d%oCBzuk6BtVp_IIya2PYuSLF0`P96pwI`ucP+j77mS1?*QcnWrt`m_4xa zT0CRy!fe(5?^ue|8s;hK)I&x9~9&Y`_H+`Z}RWOnF; zO81D{`0>A0$*)$qAaFd+1REih|GszqmoqPczWZDJRtcXA>)+a>gP!Adu^qD5K#A}n z6^xdmN!ZXKPS-FA5TH{pEy~L)S>YX9iliUtr|$t-*Xtx-85;!jbwM`&C#FSVVYiSc>O-|Rv%|B zJo}?b$dJqi-}p@?1BilKxN4JL-)*%rS~DP0{HoJFS&VC1c8NnbOlsGh-wBa+mmN1raMMpa{|rf_>RQ zt!gaak@oq=7(!W_znelmj>#yJxLGD8zi>V}lJdNVC(if*^H%)cC>)opsTjIqIb^dv zyo&nGSSTe#=}z3zGMsQLds0c5gH^uE5MXcRy@vM-_W|oGg1JfQZ z%2pSXMGI)->b2*d)0=EEJ390Bq!7Lu2|f+dl+RNDl)ed|Mh_%$_5z#q=U>@{k_e=E zI7K2PfI2!J{isw<)q-2*-z^(W8$P7Nh;Aer_lOpp=hei=(*W_xuG8~6Wwvx%Lshq- z2kM^68@U3!73&se>J_)^;GS78{TpF`+UgPO(C$z9)rJ!~WC=7Wd#~M&0%U$PtkKsLHhp+D<9Ld0#>X3WBQ-tYl>j5U49MR?JTJP&>)=VrJ`@) z+}$rGLpb8}XgN?&LPU3*?IQ(1-`2-7bG6M@RExi9*G+6TI;t|f>^`fN@8Kf~rd7II z{tzdZz1c1kWq=-iqSKY!bTu-^C=TE(Zd7wax5Tg+AfupdQ{&&v2DCE%*l}v*B@|G} z2Am4(W?l18eKk5=|MMJHkHc4w2mzYY%wEvu64ZBPXaM9600^1=@PQDuT;7UcRA>d8 zC_V%2Ngsn31}LmVCn=N^mo%W$W$9wTq-gy6x!ckpo+$$=-EV9lW$fnB%H}U>Oc9#K zrQa+el7H*4{h#9SrNXAL?u22j0f21UnF8SYK=N%A(|mu>qmYm=KKjC(_K4ce-y)4T z9oxwJP*yWSd;W8@0t`{SX~f@!7eH|rf!VL`ukOd+B3+)HXr3h?Jt^B@e?^Ko zW%Mi}=$<^2E6|8JYP_|g3k4d%h$d#ozaGVPfXvW;##2Lz<~|}2{akzZj4Pr@ao}17 z6kN`M728M)?g5!W%o>qu-QrQfGVSUCr-uLa*?eq$t&j_;7*cbAWdgC@eXD3UWh|sSU?v+#~5iwiVLKZ znR$>4<*1 z87gu_v8O92^y7oUOL%v@DK0Y$?kk^mQ@ z8jhXDY&=7foQfQ+jjMYf>F70|CVa#mUgPR-oyTdEn@cAUJ2<0n?|C!p-$yuYyn}Z} zp;Tmgu##eE_crR0REsu|9>3S_8#FAkmLpla#&>!tFcFui7pyC2=0ZbPEUVn!nsGaw z&wwYBjTk67t?gYr&~yLYd@7)|U1aLX*N)HPAz($6`$5oI;KK@7cU>oHKP5$oZCA7I zb~!DGuU(ajnGjT7qAj1Kk*`9mp*$~FeZRPyOOh@aKj%eMQ#U+D10gZB7axk)55x|$<_xG0TSc$H%u1cf~t@!Dq~7P18_xs zW{^L3iglkSCVhaFmbtP%cUJe+@ev}&TPd%-+$q-I`j{*dUJeCSbY*;MKM^6^nI#*jUTijAap! zc$n<|rCo#|naes1-eBC|+n_FA8tS<5(*2T=Ld4IUS_2)$6neW?%TnCa!XCx<0e2mC zPnYvl8oe82zoP36*Zy=9OMMPVhi|DpyFk2d?#V*>v#r{Bm4Ewp)6b7tes47vmLTGB&7en01RGfF>g|ELVhBA_rs$^z`$%r;rwsdJ^*zj-%+GaZ5geh&*HJ>y~R zayaQJPRjt013At6$^}R}HTpOWqiHD1u81z!3-mM|b|axbK-2Jq2^PWa?MV zA?N;ojCg)YZyM_Xwcn3+z^#@*(XqPx2xJERffIl%2=56KZ8aA|nEMR1QM@jt4_JB) zZgimlR@`A}8a{cI(Ddt=E&LH3UH695|Cm1DL=GDz%BSPN-go>fT%$u3Nuei*m%gkKX);gejYB4RVM2diICRgg0UEA=?TgGu zUNBqg)K-xY^$L;7mLs<;iG44-KbJld7S8Zzd%S^DAG9B9^sS%g!esflFW!uSfFVsG z8NlD5zmx|&>r6n6b2*XKndak9IZN*IWI*HHkzugQFEM`7hf|0{g|#zy2u3V#bOEhk z%khR{r4!2@zGDK`tG7zCJN$j7L=!It_Xe!9xPgo2#ZpKt}cVvdJ1 zCC+Ych}_V9|F5F;;s6!Rwa$O@0?lOYW(m3?9Z+6q-gbzKuJFZ}(fe{uzrdy&hqDf73{=DYmp1ry2C6h0j!%~ONi4%O5O#q@@$*`j9lh* zPylrxfpHGmRx1mg(_Ats0jQ+s@I2?1mY>$&FwVA+{OWdSW?w#S_@+4UmH}0Sz2CDO ziy^4aVQio;0V+~% zH!wL*%(MYZ`yI={(da_8fIdgttf*wI_G*1qZEfqy-K-*RH{cLZ-j${1zlLT=z%cmPj0Z5KF-40k{e#Wy z-kr+Q!d-dq)4{D32W*sdanIXA0YTU2=4{^SSK2!oaR38ELu4-YTZr_$Cp|H%(oxGA zzg22_d<3RLKm}OR{0AE?~+8r1X%(1@#76 z5VW3doX$7CEG-X)$PaEg(A1T1w4eJ?#X?ad0a7O46}CqmZ{_S;Wqz)qY#@rsQ(=;O z&Tk5o%TX2}cS2woFxPbk=0{vHU5U-^Z14!;;P|3gQyF_LT24)Nf2UGnhWW{rKt(@~ z)cw-H`U>Mt(e`Y7mBRiUqr~@kB*6qKVJ1Fu$S)J5(tzRhP-1MRBSixflq(_PRiNJR zPz_j-hYfw~IvJscP+-!`+8hpQEjwDnyP7SS;s3v3ky^IqJI^~%wz=sA)>q>8Rd>|t zi=8jb-s>tZZS+z>4t)+z>tMO)f&sE4IH5&W8V$U%!1x}%OGK-dd&m9uVrrcVg>P{j z1J3g8v-s0~uF2^QSSe;K!cwho(gC(d_}|G+C635=d437(*9tGh;cy+Mo)%IpaaHzpKgjS%-EY7S7kM~GSO~=O1-yj z_UAet+wHgbcsA(v`b1L*Xo0BW#v1raTEJw&ov=*KFCOTrl2F3IJ5 z5sO}M7}*ig5qqd9ADP8b{*;R3w&ZN&v$16Y@Y|a1YX|=T&#$%aMA_W2Y9wtFJo@4PIO_PdAM;njhmdq3Z?1_BSA%lF#GY`~&dc#X_C|O>l8w@z zD;3>nh*4-*Ec=2n-00k|aEY*og)x=zp?Z1F#~VCyJgjK^i~)Z5n_N%Ony2#Nw9#>` zHd^^GNFI^*@sDhxS=H_cK(2lwh-dr}9wy1x%p@7|`Ujyb8|C-|_=_*Jkjnw_i+t?0 z%G-R?Vy)U+P^bH9)-9Hysba;K^8pET6usAfb?XQw4RiDV97Hwa#GmrNNsD3^{x-oH zR6zI?==&lX1w36)xCvGDcG-DU%$G}(9ObVo0Fg2^w-PEMWf9s~shltsUs3Ssyj-=q zyx+#Hw>x6)Qg2y#2UiM5TwjHf0lJIsH_<+g&qq?o(CpPZ*`(jIyvdl`N@I%i+psqa#OeL1Y)C?_&Zfv94eg=OgEWc z8s4k$^x5!wVBRc3`2-k|QSe`^0GDg0tgX8Y2l}ti0Dc*13@N{q594r7ppuW&n1hnn ztVs=HyGu$l4Cpt%9?A^<1*tAM%&M!B`%E6s(f*x>6AZfZ*8=WJfJ%M@TKb*7dv_X} z^Tl<^Iu4-I9@35JF#+I6s=z04*25fdV;_ADC;O>DFZDU0exj7fPMUWQdYK`CJ{mrk zLbF7t4K1xoqQU~1Zi_Qs!4icZl!{Md?4oQG8=qfPFUcErT>^8ElUn8v96VfKMcDO6 zG@9WKY3FB+B|UeifJVSJ3@8i%=#<-Md{Y2b>Q|Oq^uqwj#{!`;y>=TW5)LQM&z5qO zpW;>z2P@wviNK^NE_ky4<`c$_m?DxV_)`QMzL_5RNujklpX{oV)D9JD@-w2)KllcU zE0H>FK{1N9*_ev#kiL24`Ajz)Z%5IMkP&R{5-bf4JKQ%rgRELso(j|=n<~lIF%6RA z22u!yl}Oud4^NPySA#^+u^P_sPe5~eeZf%5(uTDSI_<9#>&unf5MshnEKSkrW2H{g z=A4O~ozFHqKI=pX`IwB%1RcgFOQI&P?iTy9c_b=mLVe4iaUtc>Lq983lp1a8b_jo# zJ_2;E_L`&rHah4}iQ%~OY09?+afTte6a-WN(2p`NG;cR^>cGU6=YB$FZE^@ivHN}K;1R{6#PtkQ|Tzxmm9d1TMJx3G$vOk*@f z#6)&KAgfiV6)M%-^%EMyuwp<4BKtR_r?fe>r=qr~9>%`hZB2__VAE|aRvgV$JOV8l z@!2CM4Mf#k68y|k4*9@+!rcu){3_^7-Gr^7&y#tj!k9SnNx2+@9HGxBfb~afHN=4! z^Qa#1&ZJ;823wHtLQ z4F5|Tj=uKGw=0CNn)Fz%B*tMLg7QmxzyTPRst^jCZIG4&VbTHOyIKy`K4KApaq9_S z=w8uVhp}cKIA#RKD6#iq4okGVDRWP;N13pcQ21G*i4MnUT#HfT;PDAKmukf6cx?@x z2k|?CU**bKXBbt5>!sxiLL()XU_3$qcU{K{4F!@mDmSzh;PSii5~QeI!DrX$WM}?( zQdWM-A59AbjE)&@qtCsjvxJO$eZ`WdfxQ7dRd1fvf>~9-T;e8FlMcAhE4e5XU*Hs=^Uiu52($QHBLl0rhr5ro$KAEh-QZeQHhUrTRn6fAh2%&Wqm< z0e5@F#JkLMEYB4aE3wa(j2Np4_}dfXf|;qx)MWB?Lh1 zZR%s^%gmUkYTCry&{1+3iog&EvsXNSB&$v3i5)z%!)CIr<*(JlNhj|1yBt%+RxhC> z+LFPj>Co;l&4*PNpo974BtPG58dYaFD3rB$CC$up#R2y5rZhgAVRb8I?Gt;O>|q5) zf(oU(FD}Wl&-(~%>gzvYs_3voV}!xlU|BcLxz@PC03IA;w@wl_IB9`40oLILkw6=Z zA_Zj%q=SHzso`EZd<}2;VL6)QqCZgAtj(toOgCP|k$~l>hhFZr9mFcqlyi0}rbYYs z6X@XEijN|H*ZJ08jTe(8ve*sXGACAANTFJu$5ZDEWeBk0j~(-T0O)U zAw2hFGk4<(9pZ8QaN)*enG+H9GV9as@ z4TBd|2GQ1u6IBDP4^%)lTB@VV&I!49{8OBBmR?eEjA}~MwKe|3rBh~I&^2REVCR!0 zEHLBJVTTS~ZS`uKb$0+BPvr#Z*5U*iQ|ZmIk&rW4l*yGZk7>8Rr1*kR^foWu$eq}0-0k2bKQ`ownO7V9(2b^M^NxT(b9MX$CHWNZ_dpQ>aC zA(9)`7}6)h`XK3)&lmgzZiv9Se;hTD1hhj*<9jLh?}1L}oXO4c@4E`JFkKH*L;rQA zKjO=py+&RYKAO8G9S0}g^%H%0vc&GV5W5YG+?u%sw`C~NJe_s5SJyX0q}n4t9IeB$ zPmBqe_J&E4dQU$r56zH^gmK^`uWB)?P}XhCm(Xw_*}At<@CERgQa*=(N3}!1?70Sf zD(wwQ{wC6u!<^Bf@AC5{4|UL)aFmy}=1o>~sS9U_fiskN_z-BvY&?G)boLpfJ{)4^|n=s_ooE?)}6~?m@J{Z!p~TDssfs)<_c0fH?l4k z4;rhv$E;;t{%PE(-1W?jg){lf5fHU{+SLPSyds$i30^>&v zXT)Q)a)qAH>)LvM`-bl1`Khg4`I`z@EvB{D{_a&)YLxL)(Y=55{1&2P43Uz}VU;8r zvt^1!j!~BjGPWUoe3pDOa@Ln?o>?((UoOfln>@(WEb~(EKtRwP)2vTx*+XB;CH}9? zvNA=*^})dP&b2NsRN(?9I=(Os!vNxX)oTmq!s1w78^taBc5=zV_`>fQcv|{CsIWq= zSJ#N2ZJFu%fcZ9W-}L{@X-J|9T|G+(4AJd?u-pu@h2WYEXk*rLOhb zT9#H5%7HBygoJ@3Nrev>(|?$6lxT0Lw++ruv@LmF-*n4(=Wo|BCL({VS!Ywkcp18`NQaWF@XLMeJN99+ubZ7209ZK2_RQ_KhY(dH~cT@R;s zjz)rj(OJ!_U0S#dh7gDN_P`G9*I5W#n?jQ>eAYAOo$7aH4+Y4E-4d}l^`}%YFv^+` zeWN;;Y*d8p(^WLh@`94#rfOLdW?b`hagYXuP2Hjul8S&O+KQsOS}fHE%t=d0H~x3% zj|*pC;fO*zL7i*PhOYtC`bX8aN)L^(hjjMe4;36#C<1{wbzUC*_T?@ZaYX1E9VfZg zuIzpzh+5Hh`PFpASfJf6JM^B&|6d#D9?e#^$8oA-Zkc+=biAr5tp-Ko)$ypA2?>o> zL?aSav~);bviuk_EO*DiIyQ`sMa@V}(qxJS%Ts(=#EMVEgS$nO``jHslBl3M>R&Pzn z%4*enPN9e22hNidptf1agO+BSmD=0exIH&ky+S`*wl|w>C~9FA+)yQl(;XxSBHB5+ zV1lEbAN>be%fD^DC}+!^PF$z9PF)lkv0LOF=7$z6xm~Zb$qSF7<$bIDb1Pe+Xyjl` zXwZYZ@szIFd_D|ZWgw?y%y^>EJlSyVz^}>-N^!3;B9Bt6I_}-If`;7kHY%=t+%L#* zA7299Mc&P}I-Xt0416(QH~_>yAPDbpZ# z`?A>2*9@?+TjoFRCq`z6@D1}cEj-QoM<;SB9pAn5DZ60c@*fz43f@LvhrR0s=`HN9zje8mjnJg03ni1;{$EH)`OOvfrJ#mnKE=k zKekDyQ^z*QgwLSEJ>`GRg3H3(%%$aLz@q4v_Y{!9xLmXe2I_LqCYxLB*Y_QSs?&h120N<@ZNE~t zpPl7{|8*`dbUMaPQe-&C1r8B(*G!r_SMJZkW#(jx6*KH*pmoz^JKRXc6RrbPE zP~|d z17ls|s@`;qmH;UG-;F<4CfLUDg4?Y=KGr!$i8D7XjE55Gq6Gc#_$P=L7T29*V!YGB zCp01H-}qj>IFlP-Nsoj!=bYI>`|CYj!C8x2Mb4m>JGwPRWG1BTe-UkNo6$j@nt*G$ zENVH$FT7nRpw-HpTxR+}dJ+8ED0gU4@zwNLDQsvf5-LwUi@M@{#SjXslMMLeAOSq6jc9CMP1fInT zpCuc=E?iv4h#UIWpu*Q&#_a0D`On>7KEMCch9_*i?3<0voOwAYJL?DrlYNbaXjC>} z=+kx!F^#%cS!QrW7eAw(W|9YNe$d;W4cVH=1w@@v8U z5)acxrKN8MEUS+ho|%%(5KEpl@R(xTB9bXja*%iOvo;!4Lc^ug+(p0nTTs_h4~b>X zk>=lYQYOXA95vUv^)2)Jn9RO=0KzJ(-~YdfUN*0p<<9lp;br}vNik|zbix$GrTDFy zl4CjXN2*#AYw-qmKb9EIoVHi4K2UPFd4c`AsMY~bE%2#rn95DPC0<@C4FWWIMsQK$ z!x1SUYWI79H8Fk$6Dj+=(bwWeNFG9$TQ^wEc>jT6cG^E4GNAn&Yj=Op-xt}B{hgw&aXUQM-C?LsYMZh=GV_727t z983>d5ebxy?g0q^_A~}rrw!+F%75|-OqHrRR{T7kr;(Kl@NQ5LC44guAs{a9?R&nZ zIq=aZjhkkH9e#DRWHj>lzH2=%;KVZ5Fa?_X2cRr)38>q%v(;w?>nRwVf}c{w7Dn;r z*=Mu0Lbqo$-WlO-)MI+ie;d>q$5nqLLMo}A%XEQr@F3O=c2%wQ-40s_m_Ut-0!)oE z$s}_A8|UT^=^f9!BbVDjr{6zFx{Vf+_r6t8j{d&z_e#`HhGBVlc8Xc@qI!%mR`^=^ zELahH{yMO8MU;}})L1lShARMuPz}xbm=n;U-&>!>>HVTWaxH^t0|kaldFY-`J5B?OQdeZx5CL{5L6Q$1||%pM7uq31Pl9 As{jB1 literal 0 HcmV?d00001 diff --git a/doc/source/orbit.rst b/doc/source/orbit.rst index 8b7529880..c4e1d8353 100644 --- a/doc/source/orbit.rst +++ b/doc/source/orbit.rst @@ -1043,6 +1043,107 @@ same example as above >>> timeit(o.integrate(ts,mp,method='dop853')) # 1.61 s ± 218 ms per loop (mean ± std. dev. of 7 runs, 1 loop each) +**NEW in v1.9** Surfaces of section +----------------------------------- + +``galpy`` can compute surfaces of section for two- and three-dimensional +orbits using a special integration method that exactly determines the +intersection between orbits and a surface of section. The implemented +method follows +`Hunter et al. (1998) `__ +in re-writing the integration in terms of an independent angular variable :math:`\psi` +that is equal to a multiple of :math:`2\pi` at intersections of the orbit with +a surface of section. For example, to compute the surface of section of a +three-dimensional orbit in an axisymmetric potential for the :math:`z=0,\,v_z>0` +surface of section, we re-write the vertical phase-space coordinates as + +.. math:: + + z & = A\,\sin \psi\,,\\ + v_z & = A\,\cos \psi\,. + +When computing the surface of section of a two-dimensional orbit using either the +:math:`x=0,\,v_x>0` or :math:`y=0,\,v_y>0` surfaces, we proceed in a similar way. +As long as the surface of section is a symmetry plane of the potential, the angle +:math:`\psi` increases monotonically with time and we can use :math:`\psi` as the +independent variable in orbit integration rather than time. This allows us to +exactly integrate to integer multiples of :math:`\psi = 2\pi` and compute the +surface of section. ``galpy`` only supports surfaces of section for the three cases +mentioned above, repeated here: + +* 3D orbits: :math:`z=0,\,v_z>0` surface; +* 2D orbits: :math:`x=0,\,v_x>0` or :math:`y=0,\,v_y>0` surfaces. + +While ``galpy`` will happily compute the surface of section for any 3D or 2D orbit, +surfaces of section make most sense in 3D for static, axisymmetric potentials and +in 2D for static, non-axisymmetric potentials, like the examples below. + +As an example, let's consider the orbit of the Sun in ``MWPotential2014``. To plot the +surface of section, simply do:: + + >>> from galpy.orbit import Orbit + >>> from galpy.potential import MWPotential2014 + >>> o= Orbit() + >>> o.plotSOS(MWPotential2014) + +which gives + +.. image:: images/orbit-sos-Sun-MWP14.png + :scale: 100 % + +To get the values plotted here, do:: + + >>> Rs, vRs= o.SOS(MWPotential2014) + +.. WARNING:: + Computing the surface of section leaves the Orbit instance in a state where its + internally-stored integrated orbit is that computed during the surface-of-section + integration (any previously integrated orbit is overwritten). However, the orbit is + only output *at intersections with the surface of section*. Furthermore, with an + Orbit like that of the Sun above whose initial condition is *not* in the surface of + section, the first point along this orbit is not the initial condition, but the first + surface-of-section crossing instead. This is because for orbits like this, internally + ``galpy`` first integrates to the first crossing and then re-starts the integration + to obtain many subsequent crossings to create the surface of section. + +Surfaces of section can also be computed for Orbit instances that contain multiple +orbits, e.g., for the following two orbits defined to have the same energy and +angular momentum:: + + >>> from galpy.orbit import Orbit + >>> from galpy.potential import MWPotential2014, evaluatePotentials + >>> def orbit_RvRELz(R,vR,E,Lz,z_init=0.,pot=None): + """Returns Orbit at (R,vR,phi=0,z=z_init) with given (E,Lz)""" + return Orbit([R,vR,Lz/R,z_init, + numpy.sqrt(2.*(E-evaluatePotentials(pot,R,z_init) + -(Lz/R)**2./2.-vR**2./2)),0.],ro=8.,vo=220.) + >>> twoorb= Orbit([orbit_RvRELz(R,0.,E,Lz,pot=MWPotential2014), + orbit_RvRELz(R,0.1,E,Lz,pot=MWPotential2014,z_init=0.1)]) + >>> twoorb.plotSOS(MWPotential2014) + +which gives + +.. image:: images/orbit-sos-twoorb-MWP14.png + :scale: 100 % + +You can increase the number of points in the surface of section using the ``ncross=`` +keyword. For example, you will see the outer locus fill in with ``ncross=1500``. + +An example in two dimensions is given by the following box orbit in a +non-axisymmetric cored logarithmic potential:: + + >>> from galpy.orbit import Orbit + >>> from galpy.potential import LogarithmicHaloPotential + >>> lp= LogarithmicHaloPotential(normalize=True,b=0.9,core=0.2) + >>> orb= Orbit([0.1,0.,lp.vcirc(0.1,phi=0.),numpy.pi/2.]) + >>> orb.plotSOS(lp,surface='y') # default is surface='x' + +which gives + +.. image:: images/orbit-sos-2Dbox.png + :scale: 100 % + + Integration of the phase-space volume --------------------------------------