diff --git a/examples/examples_fast/Modules/main.py b/examples/examples_fast/Modules/main.py index 0beebf76..06e5fde2 100644 --- a/examples/examples_fast/Modules/main.py +++ b/examples/examples_fast/Modules/main.py @@ -87,7 +87,11 @@ spacetime = xpsi.Spacetime(bounds, - values=dict(frequency = 314.0)) + values=dict(frequency = 314.0), + star_shape="AGM_14") +#Default star shape is an oblate spheroid from AlGendy & Morsink (2014, AGM_14) +#But also a spherical star can be used with: +#spacetime.star_shape = "sphere" # # Hot-spot bounds = dict(super_colatitude = (0.001, math.pi/2 - 0.001), @@ -137,7 +141,6 @@ likelihood.check(None, [-3.1603740790e+04], 1.0e-5, physical_points=[p]) - if __name__ == '__main__': start = time.time() diff --git a/examples/examples_modeling_tutorial/TestRun_NumBeam.py b/examples/examples_modeling_tutorial/TestRun_NumBeam.py index 998e3c2b..db27714c 100644 --- a/examples/examples_modeling_tutorial/TestRun_NumBeam.py +++ b/examples/examples_modeling_tutorial/TestRun_NumBeam.py @@ -109,7 +109,9 @@ def from_response_files(cls, ARF, RMF, channel_edges, max_input, radius = (3.0 * gravradius(1.0), 16.0), # equatorial radius cos_inclination = (0.0, 1.0)) # (Earth) inclination to rotation axis -spacetime = xpsi.Spacetime(bounds=bounds, values=dict(frequency=300.0)) +spacetime = xpsi.Spacetime(bounds=bounds, values=dict(frequency=300.0), star_shape="AGM_14") +#The star shape can also set to be spherical: +#spacetime.star_shape = "sphere" bounds = dict(super_colatitude = (None, None), super_radius = (None, None), diff --git a/examples/examples_modeling_tutorial/modules/CustomPrior_Beaming.py b/examples/examples_modeling_tutorial/modules/CustomPrior_Beaming.py index 964d6aed..3b51345f 100644 --- a/examples/examples_modeling_tutorial/modules/CustomPrior_Beaming.py +++ b/examples/examples_modeling_tutorial/modules/CustomPrior_Beaming.py @@ -67,7 +67,9 @@ def __call__(self, p = None): grav = xpsi.surface_radiation_field.effective_gravity(np.array([1.0, 0.0]), np.array([ref.R] * 2 ), np.array([ref.zeta] * 2), - np.array([ref.epsilon] * 2)) + np.array([ref.epsilon] * 2), + ref.star_shape) + for g in grav: if not 13.7 <= g <= 15.0: return -np.inf diff --git a/xpsi/Elsewhere.py b/xpsi/Elsewhere.py index b3f08160..df47dde3 100644 --- a/xpsi/Elsewhere.py +++ b/xpsi/Elsewhere.py @@ -217,7 +217,8 @@ def _construct_cellMesh(self, st, threads): st.r_s, st.R, st.zeta, - st.epsilon) + st.epsilon, + st.star_shape_ind) def _compute_rays(self, st, threads): """ Trace (integrate) a set of rays. diff --git a/xpsi/Everywhere.py b/xpsi/Everywhere.py index 4e4ad54d..3fe940ec 100644 --- a/xpsi/Everywhere.py +++ b/xpsi/Everywhere.py @@ -405,7 +405,8 @@ def _construct_cellMesh(self, st, threads): st.r_s, st.R, st.zeta, - st.epsilon) + st.epsilon, + st.star_shape_ind) def _calibrate_lag(self, st, photosphere): """ Calibrate lag for cell mesh and normalise by pulse period. """ diff --git a/xpsi/HotRegion.py b/xpsi/HotRegion.py index 86f7dbf4..ce32914d 100644 --- a/xpsi/HotRegion.py +++ b/xpsi/HotRegion.py @@ -823,6 +823,7 @@ def __construct_cellMesh(self, st, fast_total_counts, threads): st.Omega, st.zeta, st.epsilon, + st.star_shape_ind, self['super_radius'], self['cede_radius'], self['super_colatitude'], @@ -853,6 +854,7 @@ def __construct_cellMesh(self, st, fast_total_counts, threads): st.R, st.zeta, st.epsilon, + st.star_shape_ind, self['super_radius'], self['super_colatitude'], self['omit_radius'], @@ -885,6 +887,7 @@ def __construct_cellMesh(self, st, fast_total_counts, threads): st.R, st.zeta, st.epsilon, + st.star_shape_ind, self['cede_radius'], self['cede_colatitude'], self['super_radius'], diff --git a/xpsi/Likelihood.py b/xpsi/Likelihood.py index 0646eddd..b4b70a6a 100644 --- a/xpsi/Likelihood.py +++ b/xpsi/Likelihood.py @@ -303,6 +303,7 @@ def _driver(self, fast_mode=False, synthesise=False, force_update=False, **kwarg signal in self._signals) self._star.update(fast_total_counts, self.threads,force_update=force_update) + except xpsiError as e: if isinstance(e, HotRegion.RayError): print('Warning: HotRegion.RayError raised.') diff --git a/xpsi/Spacetime.py b/xpsi/Spacetime.py index c9d24f8b..e5c004f9 100644 --- a/xpsi/Spacetime.py +++ b/xpsi/Spacetime.py @@ -22,6 +22,11 @@ class Spacetime(ParameterSubspace): no entry for a default parameter, no initial value will be specified, but an exception will be raised if there is also no bound specified. + :param string star_shape: + A string specifying the assumed shape of the star. Options 'AGM_14' + (an oblate spheroid from Algendy & Morsink 2014) or 'sphere' are + currently allowed. + We define a property for parameters and combinations of parameters to shortcut access, given that this subspace is passed to other subspaces for model computation. We would like to access the values with fewer @@ -34,8 +39,9 @@ class Spacetime(ParameterSubspace): 'distance', 'cos_inclination'] - def __init__(self, bounds, values): + def __init__(self, bounds, values, star_shape="AGM_14"): + self.star_shape = star_shape if not isinstance(bounds, dict) or not isinstance(values, dict): raise TypeError("Both bounds and values need to be dictionaries.") @@ -77,6 +83,30 @@ def __init__(self, bounds, values): super(Spacetime, self).__init__(f, M, R, D, cosi) + @property + def star_shape(self): + """ Get the the shape of the star. """ + return self._star_shape + + @star_shape.setter + def star_shape(self,star_shape): + """ Set the shape of the star. """ + allowed_models = ["AGM_14","sphere"] + if star_shape not in allowed_models: + raise TypeError("Invalid star_shape option.") + self._star_shape = star_shape + self._star_shape_ind = allowed_models.index(star_shape) + + @property + def star_shape_ind(self): + """ Get the the index of shape model. """ + return self._star_shape_ind + + @star_shape_ind.setter + def star_shape_ind(self,star_shape_ind): + """ Get the the index of shape model. """ + self._star_shape_ind = star_shape_ind + @property def M(self): """ Get the (rotationally deformed) gravitational mass in SI. """ diff --git a/xpsi/cellmesh/global_mesh.pyx b/xpsi/cellmesh/global_mesh.pyx index 6623bf4a..4338717c 100644 --- a/xpsi/cellmesh/global_mesh.pyx +++ b/xpsi/cellmesh/global_mesh.pyx @@ -25,7 +25,8 @@ def construct_closed_cellMesh(size_t numThreads, double r_s, double R_eq, double zeta, - double epsilon): + double epsilon, + int star_shape_ind): """ Construct a closed photospheric cell mesh. @@ -56,7 +57,7 @@ def construct_closed_cellMesh(size_t numThreads, for k in range(numThreads): w[k] = gsl_integration_cquad_workspace_alloc(100) - cellArea = _2pi * integrateArea(0.0, _hpi, R_eq, epsilon, zeta, 0, w[0]) / ( numCell) + cellArea = _2pi * integrateArea(0.0, _hpi, R_eq, epsilon, zeta, star_shape_ind, 0, w[0]) / ( numCell) cellArea *= 2.0 @@ -74,6 +75,7 @@ def construct_closed_cellMesh(size_t numThreads, R_eq, epsilon, zeta, + star_shape_ind, 0, w[thread]) @@ -97,23 +99,23 @@ def construct_closed_cellMesh(size_t numThreads, thread = threadid() if (i == numCellsCompute - 1): cellColatitudes[i] = integrateArea(parallels[i - 1], _hpi, - R_eq, epsilon, zeta, 1, + R_eq, epsilon, zeta, star_shape_ind, 1, w[thread]) / eta elif (i == 0): cellColatitudes[i] = integrateArea(0.0, parallels[i], - R_eq, epsilon, zeta, 1, + R_eq, epsilon, zeta, star_shape_ind, 1, w[thread]) / eta else: cellColatitudes[i] = integrateArea(parallels[i - 1], parallels[i], - R_eq, epsilon, zeta, 1, + R_eq, epsilon, zeta, star_shape_ind, 1, w[thread]) / eta mu = cos(cellColatitudes[i]) - radius = radiusNormalised(mu, epsilon, zeta) + radius = radiusNormalised(mu, epsilon, zeta, star_shape_ind) cellRadialCoord[i] = radius * R_eq #r_s_over_r = r_s / cellRadialCoord[i] - effGrav[i] = effectiveGravity(mu, R_eq, zeta, epsilon) - f = f_theta(mu, radius, epsilon, zeta) + effGrav[i] = effectiveGravity(mu, R_eq, zeta, epsilon, star_shape_ind) + f = f_theta(mu, radius, epsilon, zeta, star_shape_ind) cos_gamma[i] = 1.0 / sqrt(1.0 + f * f) maxEmissionAngle[i] = _hpi + acos(cos_gamma[i]) diff --git a/xpsi/cellmesh/mesh.pyx b/xpsi/cellmesh/mesh.pyx index 6fb40baa..e4b8410a 100644 --- a/xpsi/cellmesh/mesh.pyx +++ b/xpsi/cellmesh/mesh.pyx @@ -27,6 +27,7 @@ def construct_spot_cellMesh(size_t numThreads, double R_eq, double zeta, double epsilon, + int star_shape_ind, double cedeRadius, double cedeColatitude, double superRadius, @@ -63,6 +64,7 @@ def construct_spot_cellMesh(size_t numThreads, R_eq, epsilon, zeta, + star_shape_ind, 0, w[0]) @@ -101,6 +103,7 @@ def construct_spot_cellMesh(size_t numThreads, R_eq, epsilon, zeta, + star_shape_ind, 0, w[thread]) @@ -207,15 +210,16 @@ def construct_spot_cellMesh(size_t numThreads, R_eq, epsilon, zeta, + star_shape_ind, 1, w[thread]) / eta mu = cos(cellColatitudes[i]) - radius = radiusNormalised(mu, epsilon, zeta) + radius = radiusNormalised(mu, epsilon, zeta, star_shape_ind) cellRadialCoord[i] = radius * R_eq #r_s_over_r = r_s / cellRadialCoord[i] - effGrav[i] = effectiveGravity(mu, R_eq, zeta, epsilon) - f = f_theta(mu, radius, epsilon, zeta) + effGrav[i] = effectiveGravity(mu, R_eq, zeta, epsilon, star_shape_ind) + f = f_theta(mu, radius, epsilon, zeta, star_shape_ind) cos_gamma[i] = 1.0 / sqrt(1.0 + f * f) maxEmissionAngle[i] = _hpi + acos(cos_gamma[i]) @@ -260,6 +264,7 @@ def construct_spot_cellMesh(size_t numThreads, R_eq, epsilon, zeta, + star_shape_ind, cedeColatitude, cedeRadius, superRadius, @@ -298,6 +303,7 @@ def construct_spot_cellMesh(size_t numThreads, R_eq, epsilon, zeta, + star_shape_ind, cedeColatitude, cedeRadius, superRadius, diff --git a/xpsi/cellmesh/mesh_tools.pxd b/xpsi/cellmesh/mesh_tools.pxd index 21d5e859..b3196d63 100644 --- a/xpsi/cellmesh/mesh_tools.pxd +++ b/xpsi/cellmesh/mesh_tools.pxd @@ -27,17 +27,21 @@ cdef double eval_cedeAzi(double theta, double phi, double psi, double THETA) noe cdef double radiusNormalised(double mu, double epsilon, - double zeta) noexcept nogil + double zeta, + int star_shape_ind) noexcept nogil + cdef double f_theta(double mu, double radiusNormed, double epsilon, - double zeta) noexcept nogil + double zeta, + int star_shape_ind) noexcept nogil cdef double integrateArea(double lower_lim, double upper_lim, double R_eq, double epsilon, double zeta, + int star_shape_ind, int av, gsl_integration_cquad_workspace *w) noexcept nogil @@ -48,6 +52,7 @@ cdef double integrateCell(double theta_a, double R_eq, double epsilon, double zeta, + int star_shape_ind, double colat, double cedeRadius, double superRadius, @@ -60,6 +65,7 @@ cdef double integrateSpot(double theta_a, double R_eq, double epsilon, double zeta, + int star_shape_ind, double colat, double cedeRadius, double superRadius, diff --git a/xpsi/cellmesh/mesh_tools.pyx b/xpsi/cellmesh/mesh_tools.pyx index a6284588..d4b259e8 100644 --- a/xpsi/cellmesh/mesh_tools.pyx +++ b/xpsi/cellmesh/mesh_tools.pyx @@ -15,40 +15,51 @@ cdef double _h = 6.62607004e-34 cdef double radiusNormalised(double mu, double epsilon, - double zeta) noexcept nogil: + double zeta, + int star_shape_ind) noexcept nogil: """ Calculate the normalised radius (R(mu)/R_eq) based on the oblateness approximation from AlGendy & - Morsink (2014) (see Eq. 20 and Table 1 for coefficient values). This function models the shape of - a rotating neutron star based on its spin and compactness. It uses the colatitude (`mu`), the - dimensionless spin parameter (`epsilon`), and the compactness (`zeta`) to compute the normalized - radius. + Morsink (2014) (see Eq. 20 and Table 1 for coefficient values), or assuming a spherical star. + In the former case, this function models the shape of a rotating neutron star based on its spin and compactness. + It uses the colatitude (`mu`), the dimensionless spin parameter (`epsilon`), and the compactness + (`zeta`) to compute the normalised radius. - :param double mu: + :param double mu: The cosine of the colatitude angle (theta). This parameter defines the angular position on the star's surface, where mu is between -1 and 1. - :param double epsilon: + :param double epsilon: The dimensionless spin parameter, defined as (omega**2 * R_eq**3) / (G * M) (see Eq. 2 or Eq. 10 from Morsink et al. (2007)). - :param double zeta: + :param double zeta: The compactness parameter, defined as (G * M) / (R_eq * c**2) (see Eq. 1 or Eq. 9 from Morsink et al. (2007)). + :param int star_shape_ind: + An integer flag that corresponds to either an oblate (0) (from Algendy & Morsink 2014) or + to a spherical star (1). + :return: The normalized radius, defined as R(theta)/R_eq, which represents the radial distance from the center of the star at a given colatitude, scaled by the equatorial radius. :rtype: double """ - return 1.0 + epsilon * (-0.788 + 1.030 * zeta) * mu * mu + if(star_shape_ind == 0): # AlGendy & Morsink 2014 oblateness + return 1.0 + epsilon * (-0.788 + 1.030 * zeta) * mu * mu + elif(star_shape_ind == 1): # A spherical star + return 1.0 + else: + raise TypeError("Invalid star_shape option!") cdef double f_theta(double mu, double radiusNormed, double epsilon, - double zeta) noexcept nogil: + double zeta, + int star_shape_ind) noexcept nogil: """ - Calculate f(theta) based on Equation 3 of Morsink et al. (2007). + Calculate f(theta) based on Equation 3 of Morsink et al. (2007), or set f=0 in case of a sphere. :param double mu: The cosine of the colatitude angle (theta). This parameter defines the angular @@ -63,6 +74,10 @@ cdef double f_theta(double mu, :param double zeta: The compactness parameter, defined as (G * M) / (R_eq * c**2) (see Eq. 9). + :param int star_shape_ind: + An integer flag that corresponds to either an oblate (0) (from Algendy & Morsink 2014) or + to a spherical star (1). + :return: The computed value of f_theta, representing the normalized derivative of the radius. :rtype: double @@ -72,7 +87,12 @@ cdef double f_theta(double mu, radiusDerivNormed = -2.0 * epsilon * (-0.788 + 1.030 * zeta) * mu * sqrt(1.0 - mu * mu) - return radiusDerivNormed / (radiusNormed * sqrt(1.0 - 2.0 * zeta / radiusNormed)) + if(star_shape_ind == 0): # AlGendy & Morsink 2014 oblateness + return radiusDerivNormed / (radiusNormed * sqrt(1.0 - 2.0 * zeta / radiusNormed)) + elif(star_shape_ind == 1): # A spherical star + return 0.0 + else: + raise TypeError("Invalid star_shape option!") cdef double integrand(double theta, void *params) noexcept nogil: @@ -85,9 +105,10 @@ cdef double integrand(double theta, void *params) noexcept nogil: double epsilon = ( params)[0] double zeta = ( params)[1] int av = ( params)[2] + int star_shape_ind = ( params)[3] - radius = radiusNormalised(mu, epsilon, zeta) - f = f_theta(mu, radius, epsilon, zeta) + radius = radiusNormalised(mu, epsilon, zeta, star_shape_ind) + f = f_theta(mu, radius, epsilon, zeta, star_shape_ind) if (av == 0): eval_integrand = radius * radius * sqrt(1.0 + f * f) * sin(theta) @@ -101,16 +122,18 @@ cdef double integrateArea(double lower_lim, double R_eq, double epsilon, double zeta, + int star_shape_ind, int av, gsl_integration_cquad_workspace *w) noexcept nogil: - cdef double params[3] + cdef double params[4] cdef double integral, error cdef size_t nevals params[0] = epsilon params[1] = zeta params[2] = av + params[3] = star_shape_ind cdef gsl_function F F.function = integrand @@ -290,12 +313,13 @@ cdef double cell_integrand(double theta, void *params) noexcept nogil: cdef double phi_a = (params)[7][0] cdef double phi_b = (params)[8][0] cdef int verbose = (params)[9][0] + cdef int star_shape_ind = (params)[10][0] cdef double mu = cos(theta) cdef double radius, f - radius = radiusNormalised(mu, epsilon, zeta) - f = f_theta(mu, radius, epsilon, zeta) + radius = radiusNormalised(mu, epsilon, zeta, star_shape_ind) + f = f_theta(mu, radius, epsilon, zeta, star_shape_ind) cdef double aLB, aUB, cLB, cUB, delta_phi @@ -407,6 +431,7 @@ cdef double integrateCell(double theta_a, double R_eq, double epsilon, double zeta, + int star_shape_ind, double colat, double cedeRadius, double superRadius, @@ -414,7 +439,7 @@ cdef double integrateCell(double theta_a, double superCentreColat, gsl_cq_work *w) noexcept nogil: - cdef void *params[10] + cdef void *params[11] cdef double integral, error cdef int verbose = 0 cdef size_t nevals, i @@ -429,6 +454,7 @@ cdef double integrateCell(double theta_a, params[7] = &phi_a params[8] = &phi_b params[9] = &verbose + params[10] = &star_shape_ind cdef gsl_function F F.function = &cell_integrand @@ -712,13 +738,14 @@ cdef double spot_integrand(double theta, void *params) noexcept nogil: cdef int Lorentz = ((params)[7]) cdef double R_eq = (params)[8] cdef double Omega = (params)[9] + cdef int star_shape_ind = ((params)[10]) cdef double r_s_over_r cdef double mu = cos(theta) cdef double radius, f - radius = radiusNormalised(mu, epsilon, zeta) - f = f_theta(mu, radius, epsilon, zeta) + radius = radiusNormalised(mu, epsilon, zeta, star_shape_ind) + f = f_theta(mu, radius, epsilon, zeta, star_shape_ind) cdef double aLB, aUB, cLB, cUB, delta_phi @@ -789,6 +816,7 @@ cdef double integrateSpot(double theta_a, double R_eq, double epsilon, double zeta, + int star_shape_ind, double colat, double cedeRadius, double superRadius, @@ -799,7 +827,7 @@ cdef double integrateSpot(double theta_a, gsl_cq_work *w) noexcept nogil: - cdef double params[10] + cdef double params[11] cdef double integral, error cdef size_t nevals @@ -813,6 +841,7 @@ cdef double integrateSpot(double theta_a, params[7] = Lorentz params[8] = R_eq params[9] = Omega + params[10] = star_shape_ind cdef gsl_function F F.function = &spot_integrand @@ -867,6 +896,7 @@ def allocate_cells(size_t num_cells, double Omega, double zeta, double epsilon, + int star_shape_ind, double superRadius, double cedeRadius, double superColatitude, @@ -927,6 +957,7 @@ def allocate_cells(size_t num_cells, R_eq, epsilon, zeta, + star_shape_ind, 0, w) @@ -935,6 +966,7 @@ def allocate_cells(size_t num_cells, R_eq, epsilon, zeta, + star_shape_ind, superColatitude, superRadius, holeRadius, @@ -1012,6 +1044,7 @@ def allocate_cells(size_t num_cells, R_eq, epsilon, zeta, + star_shape_ind, 0, w) @@ -1025,6 +1058,7 @@ def allocate_cells(size_t num_cells, R_eq, epsilon, zeta, + star_shape_ind, cedeColatitude, cedeRadius, superRadius, diff --git a/xpsi/cellmesh/polar_mesh.pyx b/xpsi/cellmesh/polar_mesh.pyx index 58a5ebde..edebcf27 100644 --- a/xpsi/cellmesh/polar_mesh.pyx +++ b/xpsi/cellmesh/polar_mesh.pyx @@ -27,6 +27,7 @@ def construct_polar_cellMesh(size_t numThreads, double R_eq, double zeta, double epsilon, + int star_shape_ind, double cedeRadius, double cedeColatitude, double superRadius, @@ -67,6 +68,7 @@ def construct_polar_cellMesh(size_t numThreads, R_eq, epsilon, zeta, + star_shape_ind, 0, w[0]) @@ -101,6 +103,7 @@ def construct_polar_cellMesh(size_t numThreads, R_eq, epsilon, zeta, + star_shape_ind, 0, w[thread]) @@ -206,15 +209,16 @@ def construct_polar_cellMesh(size_t numThreads, R_eq, epsilon, zeta, + star_shape_ind, 1, w[thread]) / eta mu = cos(cellColatitudes[i]) - radius = radiusNormalised(mu, epsilon, zeta) + radius = radiusNormalised(mu, epsilon, zeta, star_shape_ind) cellRadialCoord[i] = radius * R_eq #r_s_over_r = r_s / cellRadialCoord[i] - effGrav[i] = effectiveGravity(mu, R_eq, zeta, epsilon) - f = f_theta(mu, radius, epsilon, zeta) + effGrav[i] = effectiveGravity(mu, R_eq, zeta, epsilon, star_shape_ind) + f = f_theta(mu, radius, epsilon, zeta, star_shape_ind) cos_gamma[i] = 1.0 / sqrt(1.0 + f * f) maxEmissionAngle[i] = _hpi + acos(cos_gamma[i]) @@ -270,6 +274,7 @@ def construct_polar_cellMesh(size_t numThreads, R_eq, epsilon, zeta, + star_shape_ind, cedeColatitude, cedeRadius, superRadius, @@ -308,6 +313,7 @@ def construct_polar_cellMesh(size_t numThreads, R_eq, epsilon, zeta, + star_shape_ind, cedeColatitude, cedeRadius, superRadius, diff --git a/xpsi/pixelmesh/geometricConfiguration.pxd b/xpsi/pixelmesh/geometricConfiguration.pxd index 7efaaaea..caf3ba9f 100644 --- a/xpsi/pixelmesh/geometricConfiguration.pxd +++ b/xpsi/pixelmesh/geometricConfiguration.pxd @@ -12,3 +12,4 @@ ctypedef struct _GEOM: double inclination double d double b_max + int star_shape_ind diff --git a/xpsi/pixelmesh/integrator.pyx b/xpsi/pixelmesh/integrator.pyx index db416f9b..dc962403 100644 --- a/xpsi/pixelmesh/integrator.pyx +++ b/xpsi/pixelmesh/integrator.pyx @@ -127,6 +127,7 @@ def integrate(size_t numThreads, GEOM.asq = a * a GEOM.kappa = kappa GEOM.inclination = inclination + GEOM.star_shape_ind=0 #Assuming only AlGendy & Morsink 2014 oblateness for the image plane integrator. if R_eq > 1.5 * r_s: # compare to Schwarzschild photon sphere GEOM.b_max = R_eq / sqrt(1.0 - r_s / R_eq) diff --git a/xpsi/surface_radiation_field/archive/local_variables/PDT_U.pyx b/xpsi/surface_radiation_field/archive/local_variables/PDT_U.pyx index e2ffd8ed..bd19bd6e 100644 --- a/xpsi/surface_radiation_field/archive/local_variables/PDT_U.pyx +++ b/xpsi/surface_radiation_field/archive/local_variables/PDT_U.pyx @@ -175,6 +175,7 @@ cdef int eval_local_variables(double theta, local_vars[1] = effectiveGravity(cos(theta), GEOM.R_eq, GEOM.zeta, - GEOM.epsilon) + GEOM.epsilon, + GEOM.star_shape_ind) return SUCCESS diff --git a/xpsi/surface_radiation_field/archive/local_variables/fieldlines.pyx b/xpsi/surface_radiation_field/archive/local_variables/fieldlines.pyx index a74bde1b..ec2ea0a5 100644 --- a/xpsi/surface_radiation_field/archive/local_variables/fieldlines.pyx +++ b/xpsi/surface_radiation_field/archive/local_variables/fieldlines.pyx @@ -146,6 +146,7 @@ cdef int eval_local_variables(double theta, local_vars[1] = effectiveGravity(cos(theta), GEOM.R_eq, GEOM.zeta, - GEOM.epsilon) + GEOM.epsilon, + GEOM.star_shape_ind) return SUCCESS diff --git a/xpsi/surface_radiation_field/archive/local_variables/r_mode.pyx b/xpsi/surface_radiation_field/archive/local_variables/r_mode.pyx index fea31600..12f914a4 100644 --- a/xpsi/surface_radiation_field/archive/local_variables/r_mode.pyx +++ b/xpsi/surface_radiation_field/archive/local_variables/r_mode.pyx @@ -99,6 +99,7 @@ cdef int eval_local_variables(double theta, local_vars[1] = effectiveGravity(cos(theta), GEOM.R_eq, GEOM.zeta, - GEOM.epsilon) + GEOM.epsilon, + GEOM.star_shape_ind) return SUCCESS diff --git a/xpsi/surface_radiation_field/archive/local_variables/tutorial_spot_and_ring.pyx b/xpsi/surface_radiation_field/archive/local_variables/tutorial_spot_and_ring.pyx index 694690dd..fd694a60 100644 --- a/xpsi/surface_radiation_field/archive/local_variables/tutorial_spot_and_ring.pyx +++ b/xpsi/surface_radiation_field/archive/local_variables/tutorial_spot_and_ring.pyx @@ -91,6 +91,7 @@ cdef int eval_local_variables(double theta, local_vars[1] = effectiveGravity(cos(theta), GEOM.R_eq, GEOM.zeta, - GEOM.epsilon) + GEOM.epsilon, + GEOM.star_shape_ind) return SUCCESS diff --git a/xpsi/surface_radiation_field/core.pyx b/xpsi/surface_radiation_field/core.pyx index 67d54a77..83fda33c 100644 --- a/xpsi/surface_radiation_field/core.pyx +++ b/xpsi/surface_radiation_field/core.pyx @@ -431,6 +431,7 @@ def intensity_from_globals(double[::1] energies, GEOM.R_eq = R_eq GEOM.epsilon = epsilon GEOM.zeta = zeta + GEOM.star_shape_ind = 0 #assuming only oblate stars (Algendy & Morsink 2014) for global variables. cdef fptr_init init_ptr = init_hot cdef fptr_free free_ptr = free_hot @@ -516,8 +517,10 @@ def intensity_from_globals(double[::1] energies, def effective_gravity(double[::1] cos_colatitude, double[::1] R_eq, double[::1] zeta, - double[::1] epsilon): + double[::1] epsilon, + str star_shape = "AGM_14"): """ Approximate local effective gravity using a universal relation. + (or a spherical star if setting star_shape="sphere") See Morsink et al. (2007), AlGendy & Morsink (2014), and Bogdanov et al. (2019, ApJL, 887, L26). @@ -541,6 +544,10 @@ def effective_gravity(double[::1] cos_colatitude, A 1D :class:`numpy.ndarray` of a dimensionless function of stellar properties. See :attr:`~xpsi.Spacetime.Spacetime.epsilon`. + :param double[::1] star_shape: + A string defining the shape model of the star ("AGM_14" or "sphere"). + See :attr:`~xpsi.Spacetime.Spacetime.star_shape`. + :returns: A 1D :class:`numpy.ndarray` of the logarithms (base 10) of the effective gravities in units of cm/s^2. @@ -552,10 +559,14 @@ def effective_gravity(double[::1] cos_colatitude, cdef unsigned int i + allowed_models = ["AGM_14","sphere"] + cdef unsigned int star_shape_ind = allowed_models.index(star_shape) + for i in range(gravity.shape[0]): gravity[i] = effectiveGravity(cos_colatitude[i], R_eq[i], zeta[i], - epsilon[i]) + epsilon[i], + star_shape_ind) return np.asarray(gravity, dtype = np.double, order = 'C') diff --git a/xpsi/surface_radiation_field/effective_gravity_universal.pxd b/xpsi/surface_radiation_field/effective_gravity_universal.pxd index 7cba3be9..08d43736 100644 --- a/xpsi/surface_radiation_field/effective_gravity_universal.pxd +++ b/xpsi/surface_radiation_field/effective_gravity_universal.pxd @@ -1,4 +1,5 @@ cdef double effectiveGravity(double mu, double R_eq, double x, - double epsilon) noexcept nogil + double epsilon, + int star_shape_ind) noexcept nogil diff --git a/xpsi/surface_radiation_field/effective_gravity_universal.pyx b/xpsi/surface_radiation_field/effective_gravity_universal.pyx index b070eaf6..5f1fd0c3 100644 --- a/xpsi/surface_radiation_field/effective_gravity_universal.pyx +++ b/xpsi/surface_radiation_field/effective_gravity_universal.pyx @@ -10,10 +10,12 @@ cdef double c = 2.99792458e8 cdef double effectiveGravity(double mu, double R_eq, double x, - double epsilon) noexcept nogil: + double epsilon, + int star_shape_ind) noexcept nogil: """ Calculate the effective surface gravity log-likelihood value based on the approximation from - AlGendy & Morsink (2014) (see Eq. 21 and Table 5 for coefficient values). + AlGendy & Morsink (2014) (see Eq. 21 and Table 5 for coefficient values), or assuming a spherical + star. This function computes the logarithmic surface gravity (log10(g / g0)) for a rotating neutron star, incorporating the effects of spin and compactness. The result is expressed in units compatible with @@ -34,6 +36,10 @@ cdef double effectiveGravity(double mu, The dimensionless spin parameter, defined as (omega**2 * R_eq**3) / (G * M) (see Eq. 2 or Eq. 10 from Morsink et al. (2007)). + :param int star_shape_ind: + An integer flag that corresponds to either an oblate (0) (from Algendy & Morsink 2014) or + to a spherical star (1). + :return: The effective surface gravities in log10(g*g0). The result is scaled to centimeters (conversion factor 2.0 is applied) to match the format used in tabulated atmosphere models. @@ -57,4 +63,9 @@ cdef double effectiveGravity(double mu, g += (c_p + d_p + f_p - d_60) * esq * mu * mu g += d_60 * esq * fabs(mu) - return log10(g * g_0) + 2.0 + if(star_shape_ind == 0): # AlGendy & Morsink 2014 oblateness + return log10(g * g_0) + 2.0 + elif(star_shape_ind == 1): # A spherical star + return log10(g_0) + 2.0 + else: + raise TypeError("Invalid star_shape option!") diff --git a/xpsi/surface_radiation_field/local_variables.pyx b/xpsi/surface_radiation_field/local_variables.pyx index e2ffd8ed..bd19bd6e 100644 --- a/xpsi/surface_radiation_field/local_variables.pyx +++ b/xpsi/surface_radiation_field/local_variables.pyx @@ -175,6 +175,7 @@ cdef int eval_local_variables(double theta, local_vars[1] = effectiveGravity(cos(theta), GEOM.R_eq, GEOM.zeta, - GEOM.epsilon) + GEOM.epsilon, + GEOM.star_shape_ind) return SUCCESS