Skip to content

Commit

Permalink
Options to allow spherical star (but still not ignore time delays) (#476
Browse files Browse the repository at this point in the history 
)

* First steps to make a spherical star option.

* Spherical star option finished.

However, not implemented for the image plane integrator and for the module generator.

* Adjusted some minor syntax.

* Forgot to push the syntax change in Spacetime.py

---------

Co-authored-by: Mariska Hoogkamer <mariska.98@hotmail.com>
  • Loading branch information
@thjsal @MHoogkamer
thjsal and MHoogkamer authored Dec 6, 2024
1 parent fac5459 commit cf2d414
Show file tree
Hide file tree
Showing 23 changed files with 182 additions and 55 deletions.
7 changes: 5 additions & 2 deletions examples/examples_fast/Modules/main.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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),
Expand Down Expand Up @@ -137,7 +141,6 @@

likelihood.check(None, [-3.1603740790e+04], 1.0e-5, physical_points=[p])


if __name__ == '__main__':

start = time.time()
Expand Down
4 changes: 3 additions & 1 deletion examples/examples_modeling_tutorial/TestRun_NumBeam.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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),
Expand Down
4 changes: 3 additions & 1 deletion examples/examples_modeling_tutorial/modules/CustomPrior_Beaming.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
3 changes: 2 additions & 1 deletion xpsi/Elsewhere.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down
3 changes: 2 additions & 1 deletion xpsi/Everywhere.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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. """
Expand Down
3 changes: 3 additions & 0 deletions xpsi/HotRegion.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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'],
Expand Down Expand Up @@ -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'],
Expand Down Expand Up @@ -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'],
Expand Down
1 change: 1 addition & 0 deletions xpsi/Likelihood.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.')
Expand Down
32 changes: 31 additions & 1 deletion xpsi/Spacetime.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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.")
Expand Down Expand Up @@ -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. """
Expand Down
18 changes: 10 additions & 8 deletions xpsi/cellmesh/global_mesh.pyx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -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]) / (<double> numCell)
cellArea = _2pi * integrateArea(0.0, _hpi, R_eq, epsilon, zeta, star_shape_ind, 0, w[0]) / (<double> numCell)

cellArea *= 2.0

Expand All @@ -74,6 +75,7 @@ def construct_closed_cellMesh(size_t numThreads,
R_eq,
epsilon,
zeta,
star_shape_ind,
0,
w[thread])

Expand All @@ -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])

Expand Down
12 changes: 9 additions & 3 deletions xpsi/cellmesh/mesh.pyx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -63,6 +64,7 @@ def construct_spot_cellMesh(size_t numThreads,
R_eq,
epsilon,
zeta,
star_shape_ind,
0,
w[0])

Expand Down Expand Up @@ -101,6 +103,7 @@ def construct_spot_cellMesh(size_t numThreads,
R_eq,
epsilon,
zeta,
star_shape_ind,
0,
w[thread])

Expand Down Expand Up @@ -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])

Expand Down Expand Up @@ -260,6 +264,7 @@ def construct_spot_cellMesh(size_t numThreads,
R_eq,
epsilon,
zeta,
star_shape_ind,
cedeColatitude,
cedeRadius,
superRadius,
Expand Down Expand Up @@ -298,6 +303,7 @@ def construct_spot_cellMesh(size_t numThreads,
R_eq,
epsilon,
zeta,
star_shape_ind,
cedeColatitude,
cedeRadius,
superRadius,
Expand Down
10 changes: 8 additions & 2 deletions xpsi/cellmesh/mesh_tools.pxd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -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,
Expand All @@ -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,
Expand Down
Loading

0 comments on commit cf2d414

Please sign in to comment.