2.4.10 bugfix release (#698)

* fixes implementation of gravitational redshift.
* fixes an unitialized value problem in gradients resulting in nans for effective temperatures.
* minor updates to passband exporting to support upcoming 2.5 release on the passbands server.
* allows setting SelectParameter to an array or tuple (in addition to a list).

---------

Co-authored-by: Andrej Prsa <aprsa09@gmail.com>

README.md

@@ -85,6 +85,13 @@ To understand how to use PHOEBE, please consult the [tutorials, scripts and manu
 CHANGELOG
 ----------
 
+### 2.4.10
+
+* fixes implementation of gravitational redshift.
+* fixes an unitialized value problem in gradients resulting in nans for effective temperatures.
+* minor updates to passband exporting to support upcoming 2.5 release on the passbands server.
+* allows setting SelectParameter to an array or tuple (in addition to a list).
+
 ### 2.4.9 - asynchronous spots bugfix
 
 * fixes bug introduced in 2.4.8 and ensures that temperatures are recomputed for spots when the star is rotating asynchronously.

phoebe/__init__.py

@@ -17,7 +17,7 @@
 
 """
 
-__version__ = '2.4.9'
+__version__ = '2.4.10'
 
 import os as _os
 import sys as _sys

phoebe/atmospheres/passbands.py

@@ -310,7 +310,7 @@ def save(self, archive, overwrite=True, update_timestamp=True, history_entry='')
         data.append(primary_hdu)
 
         # Tables:
-        ptf_table = Table(self.ptf_table)
+        atms = np.unique([content.split(':')[0] for content in self.content])
         data.append(fits.table_to_hdu(Table(self.ptf_table, meta={'extname': 'PTFTABLE'})))
 
         if 'blackbody:Inorm' in self.content:
@@ -322,7 +322,7 @@ def save(self, archive, overwrite=True, update_timestamp=True, history_entry='')
             data.append(fits.table_to_hdu(Table({'ebv': self._bb_extinct_axes[1]}, meta={'extname': 'BB_EBVS'})))
             data.append(fits.table_to_hdu(Table({'rv': self._bb_extinct_axes[2]}, meta={'extname': 'BB_RVS'})))
 
-        if 'ck2004:Inorm' in self.content:
+        if 'ck2004' in atms:
             ck_teffs, ck_loggs, ck_abuns = self._ck2004_axes
             data.append(fits.table_to_hdu(Table({'teff': ck_teffs}, meta={'extname': 'CK_TEFFS'})))
             data.append(fits.table_to_hdu(Table({'logg': ck_loggs}, meta={'extname': 'CK_LOGGS'})))
@@ -338,7 +338,7 @@ def save(self, archive, overwrite=True, update_timestamp=True, history_entry='')
             data.append(fits.table_to_hdu(Table({'ebv': ck_ebvs}, meta={'extname': 'CK_EBVS'})))
             data.append(fits.table_to_hdu(Table({'rv': ck_rvs}, meta={'extname': 'CK_RVS'})))
 
-        if 'phoenix:Inorm' in self.content:
+        if 'phoenix' in atms:
             ph_teffs, ph_loggs, ph_abuns = self._phoenix_axes
             data.append(fits.table_to_hdu(Table({'teff': ph_teffs}, meta={'extname': 'PH_TEFFS'})))
             data.append(fits.table_to_hdu(Table({'logg': ph_loggs}, meta={'extname': 'PH_LOGGS'})))

phoebe/backend/universe.py

@@ -1734,11 +1734,14 @@ def _populate_rv(self, dataset, **kwargs):
 
         rvs = -1*self.mesh.velocities.for_computations[:,2]
 
-
         # Gravitational redshift
         if self.do_rv_grav:
-            rv_grav = c.G*(self.mass*u.solMass)/(self.instantaneous_rpole*u.solRad)/c.c
-            # rvs are in solRad/d internally
+            # self.mass is in solar masses
+            # self.instantaneous_rpole is in Roche units, i.e. r/a
+            rpole = self.instantaneous_rpole*self.sma*u.solRad
+            rv_grav = c.G*(self.mass*u.solMass)/rpole/c.c
+
+            # rvs are in solRad/d internally, so we need to convert:
             rv_grav = rv_grav.to('solRad/d').value
 
             rvs += rv_grav

phoebe/frontend/default_bundles/default_binary.bundle

@@ -1971,7 +1971,7 @@ null,
 "compute": "phoebe01",
 "kind": "phoebe",
 "context": "compute",
-"description": "Which method to use to handle all irradiation effects (reflection, redistribution)",
+"description": "Which method to use to handle all irradiation effects (reflection)",
 "choices": [
 "none",
 "wilson",

phoebe/frontend/default_bundles/default_contact_binary.bundle

@@ -2162,7 +2162,7 @@ null,
 "compute": "phoebe01",
 "kind": "phoebe",
 "context": "compute",
-"description": "Which method to use to handle all irradiation effects (reflection, redistribution)",
+"description": "Which method to use to handle all irradiation effects (reflection)",
 "choices": [
 "none",
 "wilson",

phoebe/frontend/default_bundles/default_star.bundle

@@ -749,7 +749,7 @@ null,
 "compute": "phoebe01",
 "kind": "phoebe",
 "context": "compute",
-"description": "Which method to use to handle all irradiation effects (reflection, redistribution)",
+"description": "Which method to use to handle all irradiation effects (reflection)",
 "choices": [
 "none",
 "wilson",

phoebe/lib/libphoebe.cpp

@@ -3573,10 +3573,12 @@ static PyObject *rotstar_misaligned_gradOmega_only(PyObject *self, PyObject *arg
 
   auto fname = "rotstar_misaligned_gradOmega_only"_s;
 
+  if (verbosity_level >= 4)
+    report_stream << fname << "::START" << std::endl;
+
   double p[5];
 
   PyObject *o_misalignment;
-
   PyArrayObject *X;
 
   if (!PyArg_ParseTuple(args, "dOO!",
@@ -3588,8 +3590,6 @@ static PyObject *rotstar_misaligned_gradOmega_only(PyObject *self, PyObject *arg
     return NULL;
   }
 
-  Tmisaligned_rot_star<double> b(p);
-
   if (PyFloat_Check(o_misalignment)) {
 
     double s = std::sin(PyFloat_AsDouble(o_misalignment));
@@ -3600,13 +3600,20 @@ static PyObject *rotstar_misaligned_gradOmega_only(PyObject *self, PyObject *arg
 
   } else if (PyArray_Check(o_misalignment)) {
 
-    double *s = (double*) PyArray_DATA((PyArrayObject*)o_misalignment);
-    for (int i = 0; i < 3; ++i) p[i+1] = s[i];
+    double *s = (double *) PyArray_DATA((PyArrayObject *) o_misalignment);
+    for (int i=0; i<3; ++i) p[i+1] = s[i];
   }
+  p[4] = 0;
 
-  double g[3];
+  Tmisaligned_rot_star<double> b(p);
+  double g[3], *r = (double *) PyArray_DATA(X);
 
-  b.grad_only((double*)PyArray_DATA(X), g);
+  b.grad_only(r, g);
+
+  if (verbosity_level>=4) {
+    report_stream << fname + "::g=" << g[0] << '\t' << g[1] << '\t' <<  g[2] << '\n';
+    report_stream << fname << "::END" << std::endl;
+  }
 
   return PyArray_FromVector(3, g);
 }

phoebe/parameters/compute.py

@@ -134,7 +134,7 @@ def phoebe(**kwargs):
 
     # PHYSICS
     # TODO: should either of these be per-dataset... if so: copy_for={'kind': ['rv_dep', 'lc_dep'], 'dataset': '*'}, dataset='_default' and then edit universe.py to pull for the correct dataset (will need to become dataset-dependent dictionary a la ld_func)
-    params += [ChoiceParameter(qualifier='irrad_method', value=kwargs.get('irrad_method', 'horvat'), choices=['none', 'wilson', 'horvat'], description='Which method to use to handle all irradiation effects (reflection, redistribution)')]
+    params += [ChoiceParameter(qualifier='irrad_method', value=kwargs.get('irrad_method', 'horvat'), choices=['none', 'wilson', 'horvat'], description='Which method to use to handle all irradiation effects (reflection)')]
     params += [ChoiceParameter(qualifier='boosting_method', value=kwargs.get('boosting_method', 'none'), choices=['none'], advanced=True, description='Type of boosting method')]
 
     # MESH

phoebe/parameters/parameters.py

@@ -8371,8 +8371,8 @@ def set_value(self, value, run_checks=None, **kwargs):
         if isinstance(value, str):
             value = [value]
 
-        if not isinstance(value, list):
-            raise TypeError("value must be a list of strings, received {}".format(type(value)))
+        if not isinstance(value, (list, tuple, np.ndarray)):
+            raise TypeError("value must be a list, tuple, or array of strings, received {}".format(type(value)))
 
         try:
             value = [str(v) for v in value]

setup.py

@@ -405,7 +405,7 @@ def _env_variable_bool(key, default):
     long_description = "\n".join(long_description_s[long_description_s.index("INTRODUCTION"):])
 
 setup (name = 'phoebe',
-       version = '2.4.9',
+       version = '2.4.10',
        description = 'PHOEBE 2.4',
        long_description=long_description,
        author = 'PHOEBE development team',

tests/nosetests/test_rvgrav/test_rvgrav.py

@@ -0,0 +1,60 @@
+# Test gravitational redshift for the Sun. The theoretically
+# predicted value from GTR is ~633 m/s.
+# 
+# The Sun-Jupiter system
+# expected RV semi-amplitude of the Sun is 12.5 m/s.
+
+# import numpy as np
+
+import phoebe
+from phoebe import u, c
+# import libphoebe
+
+
+def initiate_sun_jupiter_system():
+    b = phoebe.default_binary()
+
+    b['teff@primary'] = 1.*u.solTeff
+    b['requiv@primary'] = 1.*u.solRad
+    b['syncpar@primary'] = 14.61
+
+    b['period@orbit'] = 11.852*u.yr
+    b['q@orbit'] = 1.26687e17/1.3271244e20   # (GM)_J / (GM)_Sun
+    b['sma@orbit'] = 5.2*u.au
+
+    b['teff@secondary'] = (500, 'K')
+    b['requiv@secondary'] = 1.*c.R_jup
+    b['atm@secondary'] = 'blackbody'
+
+    b['ld_mode_bol@secondary'] = 'manual'
+    b['ld_func_bol@secondary'] = 'linear'
+    b['ld_coeffs_bol@secondary'] = [0.5]
+
+    b.add_dataset('rv', times=[0,], dataset='grv', passband='Johnson:V')
+
+    b['rv_grav@secondary@grv'] = True
+    b['ld_mode@secondary@grv'] = 'manual'
+    b['ld_func@secondary@grv'] = 'linear'
+    b['ld_coeffs@secondary@grv'] = [0.5]
+
+    return b
+
+
+def test_rvgrav():
+    b = initiate_sun_jupiter_system()
+
+    # first run without rv_grav:
+    b.run_compute(rv_grav=False, model='without_rvg')
+
+    # second run with rv_grav:
+    b.run_compute(rv_grav=True, model='with_rvg')
+
+    rv_diff = b['value@rvs@primary@with_rvg'][0]-b['value@rvs@primary@without_rvg'][0]
+
+    assert abs(rv_diff-0.6363) < 1e-3
+
+
+if __name__ == '__main__':
+    # logger = phoebe.logger(clevel='INFO')
+
+    test_rvgrav()