diff --git a/.github/workflows/on_pr.yml b/.github/workflows/on_pr.yml index 470f3d462..dff012328 100644 --- a/.github/workflows/on_pr.yml +++ b/.github/workflows/on_pr.yml @@ -15,8 +15,8 @@ jobs: strategy: fail-fast: false matrix: - python-version: ["3.7", "3.8", "3.9", "3.10", "3.11"] - os: [ubuntu-latest, macos-latest] # https://docs.github.com/en/actions/using-github-hosted-runners/about-github-hosted-runners + python-version: ["3.8", "3.9", "3.10", "3.11", "3.12"] + os: [ubuntu-22.04, macos-13, macos-14] # https://docs.github.com/en/actions/using-github-hosted-runners/about-github-hosted-runners toolchain: - {compiler: gcc, version: 13} # https://github.com/fortran-lang/setup-fortran @@ -52,24 +52,6 @@ jobs: python -m pip install pipdeptree pipdeptree - - name: Checkout photodynam - uses: actions/checkout@v4 - with: - repository: phoebe-project/photodynam - path: photodynam - - - name: Install photodynam - run: | - cd photodynam - make - sudo cp photodynam /usr/local/bin/ - python setup.py build && python3 setup.py install --user - cd .. - - - name: Test photodynam install - run: | - python -c "import photodynam" - - name: Setup GNU Fortran id: setup-fortran uses: fortran-lang/setup-fortran@v1 diff --git a/README.md b/README.md index 2eeb2b3f9..3c0433016 100644 --- a/README.md +++ b/README.md @@ -85,6 +85,15 @@ To understand how to use PHOEBE, please consult the [tutorials, scripts and manu CHANGELOG ---------- +### 2.4.14 + +* Fix MPI off to not broadcast if never enabled +* Fix warning message in dynesty solver +* Fix multi-compute with enabled/disabled datasets +* Fix error message in compute_ld_coeffs +* Fix segfaults in macos-14 +* Now requires C++14-compatible compiler + ### 2.4.13 * optimization: dynamical RVs avoid unnecessary meshing diff --git a/phoebe/__init__.py b/phoebe/__init__.py index 20e64c77c..f3cb4b25c 100644 --- a/phoebe/__init__.py +++ b/phoebe/__init__.py @@ -17,7 +17,7 @@ """ -__version__ = '2.4.13' +__version__ = '2.4.14' import os as _os import sys as _sys @@ -162,7 +162,7 @@ def on(self, nprocs=None): self.nprocs = nprocs def off(self): - if self.within_mpirun and self.myrank == 0: + if self.within_mpirun and self.enabled and self.myrank == 0: self.comm.bcast({'worker_command': 'release'}, root=0) self._enabled = False diff --git a/phoebe/frontend/bundle.py b/phoebe/frontend/bundle.py index 837f1b9c7..b1dc61494 100644 --- a/phoebe/frontend/bundle.py +++ b/phoebe/frontend/bundle.py @@ -4787,7 +4787,7 @@ def run_checks_solver(self, solver=None, compute=None, solution=None, figure=Non report.add_item(self, "sampling with dataset-scaled can cause unintended issues. Consider using component-coupled and marginalizing over pblum", offending_parameters.to_list()+ - [solver_ps.get_parameter(qualifier='priors' if solver in ['dynesty'] else 'init_from', **_skip_filter_checks)]+ + [solver_ps.get_parameter(qualifier='priors' if solver_kind in ['dynesty'] else 'init_from', **_skip_filter_checks)]+ addl_parameters, False, 'run_solver') @@ -10257,6 +10257,8 @@ def compute_ld_coeffs(self, compute=None, set_value=False, **kwargs): raise TypeError("compute must be a single value (string)") datasets = kwargs.pop('dataset') if 'dataset' in kwargs else self._datasets_where(compute=compute, mesh_needed=True) + if isinstance(datasets, str): + datasets = [datasets] # we'll add 'bol' to the list of default datasets... but only if bolometric is needed for irradiation compute_ps = self.get_compute(compute, **_skip_filter_checks) @@ -12139,7 +12141,7 @@ def _scale_fluxes_cfit(fluxes, scale_factor): continue fluxes = flux_param.get_value(unit=u.W/u.m**2) - fluxes = fluxes/distance**2 + l3s.get(dataset) + fluxes = fluxes/distance**2 + l3s.get(dataset, 0.0) flux_param.set_value(fluxes, ignore_readonly=True) diff --git a/phoebe/lib/Makefile b/phoebe/lib/Makefile index 622ce842b..45da1aa22 100644 --- a/phoebe/lib/Makefile +++ b/phoebe/lib/Makefile @@ -1,14 +1,14 @@ #CXX=icpc CXX=g++ #CXX=clang++ -CXXFLAGS=-O3 -Wall -std=c++11 +CXXFLAGS=-O3 -Wall -Wextra ifdef PYTHON CXXFLAGS+=$(shell $(PYTHON)-config --includes) LDFLAGS+=$(shell $(PYTHON)-config --ldflags) else -CXXFLAGS+=$(shell python-config --includes) -LDFLAGS+=$(shell python-config --ldflags) +CXXFLAGS+=$(shell python3-config --includes) +LDFLAGS+=$(shell python3-config --ldflags) endif LIBP=libphoebe diff --git a/phoebe/lib/bodies.h b/phoebe/lib/bodies.h index ce2a6d3e5..792541505 100644 --- a/phoebe/lib/bodies.h +++ b/phoebe/lib/bodies.h @@ -130,7 +130,7 @@ struct Tsphere { ret[3] = potential-value */ - void grad(T r[3], T ret[4], const bool & precision = false){ + void grad(T r[3], T ret[4], [[maybe_unused]] const bool & precision = false){ for (int i = 0; i < 3; ++i) ret[i] = 2*r[i]; @@ -138,7 +138,7 @@ struct Tsphere { } - void grad_only(T r[3], T ret[3], const bool & precision = false){ + void grad_only(T r[3], T ret[3], [[maybe_unused]] const bool & precision = false){ for (int i = 0; i < 3; ++i) ret[i] = 2*r[i]; @@ -1020,7 +1020,7 @@ struct Tmisaligned_roche { */ Tmisaligned_roche(T *params) { - + q = params[0]; F = params[1]; delta = params[2]; diff --git a/phoebe/lib/cvec.h b/phoebe/lib/cvec.h new file mode 100644 index 000000000..3143f0e7a --- /dev/null +++ b/phoebe/lib/cvec.h @@ -0,0 +1,277 @@ +# pragma once + +/* + Library supporting circular/periodic vectors defined as + by iterators + + it_first -> it_last + + or vector container. + + Author: Martin Horvat, June 2024 +*/ + + +/* + Next interator: it + 1 mod range + + Input: + it : iterator + it_first : iterator, pointer to first element + it_last : iterator, pointer to last element + + Return: + next it +*/ +template +Tit cnext(Tit it, Tit it_first, Tit it_last) { + return it == it_last ? it_first : it + 1; +} + +/* + Next interator: it + 1 mod range of vector + + Input: + it : iterator + P: vector + + Return: + next it +*/ +template +typename V::iterator cnext(typename V::iterator it, V & P) { + return cnext(it, P.begin(), P.end() - 1); +} + + +/* + Previous interator: it - 1 mod range + + Input: + it : iterator + it_first : iterator, pointer to first element + it_last : iterator, pointer to last element + + Return: + previous it +*/ +template +Tit cprev(Tit it, Tit it_first, Tit it_last) { + return it == it_first ? it_last : it - 1; +} + + +/* + Previous interator : it - 1 mod range of vector + + Input: + it : iterator + P : vector + + Return: + iterator: previous it +*/ +template +typename V::iterator cprev(typename V::iterator it, V & P) { + return cprev(it, P.begin(), P.end() - 1); +} + +/* + Copying elements from circular/periodic vector. + + Input: + it0 : iterator, pointer to first element that needs to be copied + it1 : iterator, pointer to last element that needs to be copied + it_first : iterator, pointer to first element of container + it_last : iterator, pointer to last element of container + + Return: + vector: copied elements from t0 to t1 +*/ +template +V ccopy (typename V::iterator it0, typename V::iterator it1, + typename V::iterator it_first, typename V::iterator it_last) { + V r; + + auto it = it0; + + while (true) { + r.push_back(*it); + if (it == it1) break; + it = cnext(it, it_first, it_last); + } + + return r; +} + +template +V ccopy (typename V::iterator it0, typename V::iterator it1, V & P) { + return ccopy(it0, it1, P.begin(), P.end() - 1); +} + +/* + Erasing elements from circular/periodic vector. + + Input: + it0 : iterator, pointer to first element that needs to be erased + it1 : iterator, pointer to last element that needs to be erased + it_first : iterator, pointer to first element of container + it_last : iterator, pointer to last element of container + + Return: + vector of remainder +*/ +template +V cerase (typename V::iterator it0, typename V::iterator it1, + typename V::iterator it_first, typename V::iterator it_last) { + + int size = it_last + 1 - it_first; + + std::vector mask(size, false); + + auto it = it0; + + while (true) { + mask[int(it - it_first)] = true; + if (it == it1) break; + it = cnext(it, it_first, it_last); + } + + V R; + + for (int i = 0; i < size; ++i) { + it = it_first + i; + if (!mask[i]) R.push_back(*it); + } + + return R; +} + +/* + Erasing elements from circular/periodic vector. + + Input: + it0 : iterator, pointer to first element that needs to be erased + it1 : iterator, pointer to last element that needs to be erased + P : vector + + Return: + vector of remainder +*/ +template +V cerase (typename V::iterator it0, typename V::iterator it1, V & P) { + + int size = P.size(); + + std::vector mask(size, false); + + auto it = it0, it_first = P.begin(), it_last = P.end() - 1; + + while (true) { + mask[int(it - it_first)] = true; + if (it == it1) break; + it = cnext(it, it_first, it_last); + } + + V R; + + for (int i = 0; i < size; ++i) { + it = it_first + i; + if (!mask[i]) R.push_back(*it); + } + + return R; +} + + + +/* + Split elements from circular/periodic vector into a pair of vectors + (first, second): + + for it = it0 up to it1 -> first + everything else -> second + + Input: + it0 : iterator, pointer to first element that needs to be copied into first + it1 : iterator, pointer to last element that needs to be copied into first + it_first : iterator, pointer to first element of container vector + it_last : iterator, pointer to last element of container vector + + Return: + pair of vectors (first, second) +*/ +template +std::pair + csplit (typename V::iterator it0, typename V::iterator it1, + typename V::iterator it_first, typename V::iterator it_last) { + + int size = it_last + 1 - it_first; + + std::vector mask(size, false); + + auto it = it0; + + while (true) { + mask[int(it - it_first)] = true; + if (it == it1) break; + it = cnext(it, it_first, it_last); + } + + std::pair p; + + for (int i = 0 ; i < size; ++i) { + it = it_first + i; + if (mask[i]) + p.first.push_back(*it); + else + p.second.push_back(*it); + } + + return p; +} + +/* + Split elements from circular/periodic vector into a pair of vectors + (first, second): + + for it = it0 up to it1 -> first + everything else -> second + + Input: + it0 : iterator, pointer to first element that needs to be copied to first + it1 : iterator, pointer to last element that needs to be copied to first + P : vector + + Return: + pair of vectors (first, second) +*/ +template +std::pair + csplit (typename V::iterator it0, typename V::iterator it1, V & P) { + + int size = P.size(); + + std::vector mask(size, false); + + auto it = it0, it_first = P.begin(), it_last = P.end() - 1 ; + + while (true) { + mask[int(it - it_first)] = true; + if (it == it1) break; + it = cnext(it, it_first, it_last); + } + + std::pair p; + + for (int i = 0 ; i < size; ++i) { + it = it_first + i; + if (mask[i]) + p.first.push_back(*it); + else + p.second.push_back(*it); + } + + return p; +} + diff --git a/phoebe/lib/gen_roche.h b/phoebe/lib/gen_roche.h index 9756021b5..6ae00d18b 100644 --- a/phoebe/lib/gen_roche.h +++ b/phoebe/lib/gen_roche.h @@ -228,7 +228,7 @@ namespace gen_roche { Input: Omega0 - value of potential q - mass ratio M2/M1 - F - synchronicity parameter + F - synchronicity parameter (not used) delta - separation between the two objects Return: @@ -241,7 +241,7 @@ namespace gen_roche { T poleL( const T & Omega0, const T & q, - const T & F = 1, + [[maybe_unused]] const T & F = 1, const T & delta = 1 ) { @@ -576,10 +576,6 @@ namespace gen_roche { const char *fname = "left_lobe_left_xborder"; - const int max_iter = 100; - const T eps = 2*std::numeric_limits::epsilon(); - const T min = 10*std::numeric_limits::min(); - // // Is solution is near to Lagrange point? // @@ -589,16 +585,15 @@ namespace gen_roche { if (q*(1/(1 - l) - l) - 1/l + b*l*l/2 == w) return l; // - // Cases away from Lagrange point + // Two cases away from Lagrange point // - T t; - if (w > 100) { - if (2*q < w){ // w->infty + if (2*q < w){ // assume: w->infty, 2q < w - T q2 = q*q, + T t, + q2 = q*q, s = 1/w, a[8] = {1, q, q2, b/2 + q*(1 + q2), q*(-1 + 2*b + q*(4 + q2)), @@ -608,11 +603,14 @@ namespace gen_roche { }; t = s*(a[0] + s*(a[1] + s*(a[2] + s*(a[3] + s*(a[4] + s*(a[5] + s*(a[6] + s*a[7]))))))); - t = -t; - } else if (q < w) { // w->infty, q ~ w + return polish_xborder(w, q, b, t); + } - T a = b/(1 + q), + if (q < w) { // assume: w->infty, q ~ w, q < w + + T t, + a = b/(1 + q), s = 1/w, f = q*s, f1 = 1 - f, f12 = f1*f1, f13 = f12*f1, @@ -631,20 +629,22 @@ namespace gen_roche { for (int i = 0; i < 8; ++i) C[i] = N[i]/D[i]; t = s/f1*(C[0] + s1*(C[1] + s1*(C[2] + s1*(C[3] + s1*(C[4] + s1*(C[5] + s1*(C[6] + s1*C[7]))))))); - t = -t; + return polish_xborder(w, q, b, t); } - - return polish_xborder(w, q, b, t); } + const int max_iter = 100; + const T eps = 2*std::numeric_limits::epsilon(); + const T min = 10*std::numeric_limits::min(); + const int method = 0; if (method == 0) { // Bisection on [-|l|,0] int it = 0; - long double f, x[2] = {l, 0}; + long double t, f, x[2] = {l, 0}; do { t = (x[0] + x[1])/2; @@ -676,7 +676,6 @@ namespace gen_roche { // grab smallest root positive for (auto && v : roots) if (v > 0) return -v; - } return std::numeric_limits::quiet_NaN(); @@ -705,25 +704,21 @@ namespace gen_roche { const char *fname = "left_lobe_right_xborder"; - const int max_iter = 100; - const T eps = 2*std::numeric_limits::epsilon(); - const T min = 10*std::numeric_limits::min(); - // // Is solution is near to Lagrange point? // - T l = lagrange_point_L1(q, std::sqrt(b/(1 + q)), 1.), t = l; + T l = lagrange_point_L1(q, std::sqrt(b/(1 + q)), 1.); - if (q*(1/(1 - t) - t) + 1/t + b*t*t/2 == w) return t; + if (q*(1/(1 - l) - l) + 1/l + b*l*l/2 == w) return l; // - // Cases away from Lagrange point + // Two cases away from Lagrange point // if (w > 100) { // w->infty - if (2*q < w){ + if (2*q < w){ // assume: w->infty, 2q < w T q2 = q*q, s = 1/w, @@ -734,9 +729,11 @@ namespace gen_roche { q*(1 + b*(3.5 + 21*b/4) + q*(14 + 21*b + q*(42 + q*(35 + 35*b/2 + q*(35 + q2))))) }; - t = s*(a[0] + s*(a[1] + s*(a[2] + s*(a[3] + s*(a[4] + s*(a[5] + s*(a[6] + s*a[7]))))))); + T t = s*(a[0] + s*(a[1] + s*(a[2] + s*(a[3] + s*(a[4] + s*(a[5] + s*(a[6] + s*a[7]))))))); + return polish_xborder(w, q, b, t); + } - } else if (q < w) { + if (q < w) { // assume: w->infty, q ~ w, q < w T a = b/(1 + q), s = 1/w, @@ -756,19 +753,22 @@ namespace gen_roche { for (int i = 0; i < 8; ++i) C[i] = N[i]/D[i]; - t = s/f1*(C[0] + s1*(C[1] + s1*(C[2] + s1*(C[3] + s1*(C[4] + s1*(C[5] + s1*(C[6] + s1*C[7]))))))); + T t = s/f1*(C[0] + s1*(C[1] + s1*(C[2] + s1*(C[3] + s1*(C[4] + s1*(C[5] + s1*(C[6] + s1*C[7]))))))); + return polish_xborder(w, q, b, t); } - - return polish_xborder(w, q, b, t); } + const int max_iter = 100; + const T eps = 2*std::numeric_limits::epsilon(); + const T min = 10*std::numeric_limits::min(); + const int method = 0; if (method == 0) { // Bisection on [0,l] int it = 0; - long double f, x[2] = {0, l}; + long double t, f, x[2] = {0, l}; do { t = (x[0] + x[1])/2; @@ -1699,26 +1699,26 @@ namespace gen_roche { */ template <> const double glq::c2_phi[]={0.999580064408504868392251519574,0.988810441236620370626267163745,0.93792975088501989705468221942,0.814698074016316269208833703789,0.615862835349225523334990208675,0.384137164650774476665009791325,0.185301925983683730791166296211,0.06207024911498010294531778058,0.011189558763379629373732836255,0.0004199355914951316077484804259}; - + template <> const double glq::s2c2_phi[]= {0.0004197592455941272416985899940,0.011064352538060503513192603217,0.058217533289784414692374340457,0.150965122210421127009962838622,0.236575803384838256502140293303,0.236575803384838256502140293303,0.150965122210421127009962838622,0.058217533289784414692374340457,0.011064352538060503513192603217,0.0004197592455941272416985899940}; - + template <> const double glq::sc_phi[]={0.020492330065054378968043365567,0.9997900101563852235169839221736,0.105780710733950116401421690636,0.9943894816602900746042204870849,0.249139015641830179695805173292,0.968467733528081924117038032924,0.430467102092231585658617450885,0.902606267436868770388218111745,0.619787999763446917824362871368,0.784769287975278453643040025427,0.784769287975278453643040025427,0.619787999763446917824362871368,0.902606267436868770388218111745,0.43046710209223158565861745088,0.968467733528081924117038032924,0.24913901564183017969580517329,0.994389481660290074604220487085,0.10578071073395011640142169064,0.9997900101563852235169839221736,0.02049233006505437896804336557}; template <> const double glq::weights[]={0.209454205485130325204687978,0.46951526056054717731928526,0.68828010698191943974759705,0.845926347240509469003431984,0.92841673332168682718764111324,0.92841673332168682718764111324,0.845926347240509469003431984,0.68828010698191943974759705,0.46951526056054717731928526,0.209454205485130325204687978}; - + /* Gauss-Lagrange: long double, N=10 */ template <> const long double glq::c2_phi[]={0.999580064408504868392251519574L,0.988810441236620370626267163745L,0.93792975088501989705468221942L,0.814698074016316269208833703789L,0.615862835349225523334990208675L,0.384137164650774476665009791325L,0.185301925983683730791166296211L,0.06207024911498010294531778058L,0.011189558763379629373732836255L,0.0004199355914951316077484804259L}; - + template <> const long double glq::s2c2_phi[]= {0.0004197592455941272416985899940L,0.011064352538060503513192603217L,0.058217533289784414692374340457L,0.150965122210421127009962838622L,0.236575803384838256502140293303L,0.236575803384838256502140293303L,0.150965122210421127009962838622L,0.058217533289784414692374340457L,0.011064352538060503513192603217L,0.0004197592455941272416985899940L}; - + template <> const long double glq::sc_phi[]={0.020492330065054378968043365567L,0.9997900101563852235169839221736L,0.105780710733950116401421690636L,0.9943894816602900746042204870849L,0.249139015641830179695805173292L,0.968467733528081924117038032924L,0.430467102092231585658617450885L,0.902606267436868770388218111745L,0.619787999763446917824362871368L,0.784769287975278453643040025427L,0.784769287975278453643040025427L,0.619787999763446917824362871368L,0.902606267436868770388218111745L,0.43046710209223158565861745088L,0.968467733528081924117038032924L,0.24913901564183017969580517329L,0.994389481660290074604220487085L,0.10578071073395011640142169064L,0.9997900101563852235169839221736L,0.02049233006505437896804336557L}; template <> @@ -1729,10 +1729,10 @@ namespace gen_roche { */ template <> const double glq::c2_phi[]={0.999911065396171632736007253518,0.997574886997681325232904530099,0.985854212895255771578898079594,0.953879938591329487419660590834,0.890692147886632646280933330332,0.790164039217058712341619458721,0.655400162984313217687403793532,0.5,0.344599837015686782312596206468,0.209835960782941287658380541279,0.109307852113367353719066669668,0.046120061408670512580339409166,0.014145787104744228421101920406,0.0024251130023186747670954699007,0.0000889346038283672639927464819}; - + template <> const double glq::s2c2_phi[] ={0.0000889266944646091553555375073,0.0024192318292446596704491832554,0.013945683811931480320682058875,0.043993001344330973475114779410,0.097359645579729565862303734566,0.165804830345241213606518615509,0.225850789344448888056399859737,0.250000000000000000000000000000,0.225850789344448888056399859737,0.165804830345241213606518615509,0.097359645579729565862303734566,0.043993001344330973475114779410,0.013945683811931480320682058875,0.002419231829244659670449183255,0.0000889266944646091553555375073}; - + template <> const double glq::sc_phi[]={0.009430514504965636496719714287,0.99995553170937138065232597719866,0.049245436360323529649147830490,0.9987867074594461808460040593058,0.118936063095867724375201909234,0.9929019150425966277731892615479,0.214755818102026080080238124559,0.976667772884582020710505864562,0.330617380234868102916887859767,0.943764879557738494192185923899,0.458078553070258123013941167236,0.888911716210928739911629494937,0.587026266035589571311728655155,0.809567886581670936826823880946,0.707106781186547524400844362105,0.707106781186547524400844362105,0.809567886581670936826823880946,0.587026266035589571311728655155,0.888911716210928739911629494937,0.458078553070258123013941167236,0.943764879557738494192185923899,0.33061738023486810291688785977,0.976667772884582020710505864562,0.21475581810202608008023812456,0.992901915042596627773189261548,0.11893606309586772437520190923,0.9987867074594461808460040593058,0.04924543636032352964914783049,0.9999555317093713806523259771987,0.00943051450496563649671971429}; template <> @@ -1745,7 +1745,7 @@ namespace gen_roche { const long double glq::c2_phi[]={0.999911065396171632736007253518L,0.997574886997681325232904530099L,0.985854212895255771578898079594L,0.953879938591329487419660590834L,0.890692147886632646280933330332L,0.790164039217058712341619458721L,0.655400162984313217687403793532L,0.5L,0.344599837015686782312596206468L,0.209835960782941287658380541279L,0.109307852113367353719066669668L,0.046120061408670512580339409166L,0.014145787104744228421101920406L,0.0024251130023186747670954699007L,0.0000889346038283672639927464819L}; template <> const long double glq::s2c2_phi[] = {0.0000889266944646091553555375073L,0.0024192318292446596704491832554L,0.013945683811931480320682058875L,0.043993001344330973475114779410L,0.097359645579729565862303734566L,0.165804830345241213606518615509L,0.225850789344448888056399859737L,0.250000000000000000000000000000L,0.225850789344448888056399859737L,0.165804830345241213606518615509L,0.097359645579729565862303734566L,0.043993001344330973475114779410L,0.013945683811931480320682058875L,0.002419231829244659670449183255L,0.0000889266944646091553555375073L}; - + template <> const long double glq::sc_phi[]={0.009430514504965636496719714287L,0.99995553170937138065232597719866L,0.049245436360323529649147830490L,0.9987867074594461808460040593058L,0.118936063095867724375201909234L,0.9929019150425966277731892615479L,0.214755818102026080080238124559L,0.976667772884582020710505864562L,0.330617380234868102916887859767L,0.943764879557738494192185923899L,0.458078553070258123013941167236L,0.888911716210928739911629494937L,0.587026266035589571311728655155L,0.809567886581670936826823880946L,0.707106781186547524400844362105L,0.707106781186547524400844362105L,0.809567886581670936826823880946L,0.587026266035589571311728655155L,0.888911716210928739911629494937L,0.458078553070258123013941167236L,0.943764879557738494192185923899L,0.33061738023486810291688785977L,0.976667772884582020710505864562L,0.21475581810202608008023812456L,0.992901915042596627773189261548L,0.11893606309586772437520190923L,0.9987867074594461808460040593058L,0.04924543636032352964914783049L,0.9999555317093713806523259771987L,0.00943051450496563649671971429L}; template <> @@ -1823,13 +1823,13 @@ namespace gen_roche { b_dvol = (choice & 4u) == 4u; if (!b_area && !b_vol && !b_dvol) return; - + using real = long double; using G = glq; - + const int dim = G::n + 3; - + real d2 = delta*delta, d3 = d2*delta, d4 = d2*d2, @@ -2003,7 +2003,7 @@ namespace gen_roche { int it; - real + real s1 = t*t, s2 = (t - 1)*(t - 1), s, f, df, ds, f1, f2, g1, g2; @@ -2045,12 +2045,12 @@ namespace gen_roche { if (b_vol) v[1] = d3*y[G::n+1]/2; if (b_dvol) v[2] = d4*y[G::n+2]; } - + /* - A simplified version of directed integrated - + A simplified version of directed integrated + void area_volume_directed_integration() - */ + */ template void area_volume_integration( T v[3], @@ -2081,8 +2081,8 @@ namespace gen_roche { area_volume_directed_integration(v, choice, -1, xrange, Omega0, q, F, d, m, polish); #endif } - - + + /* Parameterize dimensionless potential (omega) of aligned Roche lobe @@ -2165,7 +2165,7 @@ namespace gen_roche { if ((mask & 1U) == 1U) v[0] = (-Wnr-std::sqrt(Wnr*Wnr-Wrr*Wnn))/Wrr; // dr/dnu if ((mask & 2U) == 2U) v[1] = 1/(Wrr*v[0] + Wnr); // d^2(dV/dOmega)/(dnu dphi) } - + /* Computing area of the surface, the volume and derivative of the volume w.r.t. potential of the semi-detached Roche. @@ -2196,8 +2196,8 @@ namespace gen_roche { * https://en.wikipedia.org/wiki/Gaussian_quadrature * https://en.wikipedia.org/wiki/Gauss–Kronrod_quadrature_formula */ - - + + template void critical_area_volume_integration( T v[3], @@ -2225,20 +2225,23 @@ namespace gen_roche { if (b_area) mask |= 4; // + 100b if (b_dvol) mask2 |= 2; // + 010b - + using real = long double; - + using G = glq; - + const int dim = G::n + 3; - - real - y[dim], k[4][dim], w[G::n], W[3], sc_nu[2], sum[3], - r[2], rt, rp, nu, + + real + y[dim], k[4][dim], w[G::n], W[3], sc_nu[2], sum[3], + r[2], rt, rp, nu, q_ = q, d2 = d*d, d3 = d2*d, b = d3*F*F*(1 + q), dnu = utils::m_pi/m; - - + + + // Initialize k to remove warnings "uninitialized in this function" + for (int i = 0; i < 4; ++i) for (int j = 0; j < dim; ++j) k[i][j] = 0; + // // Setup init point // @@ -2250,15 +2253,15 @@ namespace gen_roche { } y[G::n] = y[G::n+1] = y[G::n+2] = 0; } - + // // Integration - // - + // + nu = 0; - + for (int i = 0; i < m; ++i) { - + // 1. step sum[0] = sum[1] = sum[2] = 0; if (i == 0) { // discussing pole = L1 separately, sin(nu) = 0 @@ -2395,38 +2398,38 @@ namespace gen_roche { //#define DEBUG template bool critical_area_volume( - const unsigned &choice, + const unsigned &choice, const T & q, const T & F, const T & delta, T & OmegaC, T av[3]) { - + #if defined(DEBUG) const char *fname = "critical_area_volume"; #endif - + T L1 = lagrange_point_L1(q, F, delta); - + OmegaC = potential_on_x_axis(L1, q, F, delta); - + #if defined(DEBUG) std::cerr.precision(16); - std::cerr << fname << "::OmegaC=" << OmegaC << " L1=" << L1 << '\n'; + std::cerr << fname << "::OmegaC=" << OmegaC << " L1=" << L1 << '\n'; #endif - + // compute volume and d(volume)/dOmega critical_area_volume_integration(av, choice, L1, q, F, delta, 1<<10); - + #if defined(DEBUG) - std::cerr << fname << "::av=" << av[0] << ':' << av[1] << ':' << av[2] << '\n'; - #endif + std::cerr << fname << "::av=" << av[0] << ':' << av[1] << ':' << av[2] << '\n'; + #endif return true; } #if defined(DEBUG) #undef DEBUG #endif - + /* Computing surface area and the volume of the primary Roche lobe in the limit of high w=delta*Omega. It should precise at least up to 5.5 digits for diff --git a/phoebe/lib/ld_models.h b/phoebe/lib/ld_models.h index 4f0e54a98..c0bc5b55c 100644 --- a/phoebe/lib/ld_models.h +++ b/phoebe/lib/ld_models.h @@ -115,7 +115,7 @@ struct TLDuniform: TLDmodel { this->nr_par = 0; } - T D(const T & mu) const { return 1; } + T D([[maybe_unused]] const T & mu) const { return 1; } bool check() const { return true; } bool check_strict() const { return true; } diff --git a/phoebe/lib/libphoebe.cpp b/phoebe/lib/libphoebe.cpp index 3eaf2e479..a1dcdd4b8 100644 --- a/phoebe/lib/libphoebe.cpp +++ b/phoebe/lib/libphoebe.cpp @@ -181,7 +181,7 @@ void raise_exception(const std::string & str){ Input: level */ -static PyObject *setup_verbosity(PyObject *self, PyObject *args, PyObject *keywds) { +static PyObject *setup_verbosity([[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { auto fname = "setup_verbosity"_s; @@ -209,7 +209,7 @@ static PyObject *setup_verbosity(PyObject *self, PyObject *args, PyObject *keywd /* - Insert into dictionary and deferences the inserted object + Insert into dictionary and dereferences the inserted object Ref: https://robinelvin.wordpress.com/2011/03/24/python-c-extension-memory-leak/ */ @@ -359,7 +359,7 @@ void PyArray_To3DPointVector( */ -static PyObject *roche_critical_potential(PyObject *self, PyObject *args, PyObject *keywds) { +static PyObject *roche_critical_potential([[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { // // Reading arguments @@ -488,7 +488,7 @@ static PyObject *roche_critical_potential(PyObject *self, PyObject *args, PyObje theta' - new angle between S and new z-axis ^k' */ -static PyObject *roche_misaligned_transf(PyObject *self, PyObject *args) { +static PyObject *roche_misaligned_transf([[maybe_unused]] PyObject *self, PyObject *args) { auto fname = "roche_misaligned_trans"_s; @@ -558,7 +558,7 @@ static PyObject *roche_misaligned_transf(PyObject *self, PyObject *args) { Omega_critical: float */ -static PyObject *rotstar_critical_potential(PyObject *self, PyObject *args) { +static PyObject *rotstar_critical_potential([[maybe_unused]] PyObject *self, PyObject *args) { // parse input arguments double omega; @@ -576,7 +576,7 @@ static PyObject *rotstar_critical_potential(PyObject *self, PyObject *args) { /* C++ wrapper for Python code: - Calculate critical potential of the rotating star with misaignment + Calculate critical potential of the rotating star with misalignment having potential function Omega(x,y,z; omega, s) = 1/|r| + 1/2 omega^2 |r - s (r.s)|^2 @@ -594,7 +594,7 @@ static PyObject *rotstar_critical_potential(PyObject *self, PyObject *args) { Python: - Omega_crit = rotstar_misaligned_critical_potential(omega, misalignemnt) + Omega_crit = rotstar_misaligned_critical_potential(omega, misalignment) where parameters are @@ -614,7 +614,7 @@ static PyObject *rotstar_critical_potential(PyObject *self, PyObject *args) { Omega_critical: float */ -static PyObject *rotstar_misaligned_critical_potential(PyObject *self, PyObject *args) { +static PyObject *rotstar_misaligned_critical_potential([[maybe_unused]] PyObject *self, PyObject *args) { auto fname = "rotstar_misaligned_critical_potential"_s; @@ -673,7 +673,7 @@ static PyObject *rotstar_misaligned_critical_potential(PyObject *self, PyObject h : height of the lobe's pole */ -static PyObject *roche_pole(PyObject *self, PyObject *args, PyObject *keywds) { +static PyObject *roche_pole([[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { // // Reading arguments @@ -741,7 +741,7 @@ static PyObject *roche_pole(PyObject *self, PyObject *args, PyObject *keywds) { static PyObject *roche_misaligned_pole( - PyObject *self, PyObject *args, PyObject *keywds) { + [[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { auto fname = "roche_misaligned_pole"_s; @@ -809,7 +809,7 @@ static PyObject *roche_misaligned_pole( compact primary Roche lobe: Omega_{min} (x,y,z; q, F, d, misalignment) = - min {Omega(L1(misaligment)), Omega(L2(misalignment)) } + min {Omega(L1(misalignment)), Omega(L2(misalignment)) } Python: @@ -829,7 +829,7 @@ static PyObject *roche_misaligned_pole( Omega0 - value of the Omega at (x,y,z) */ -static PyObject *roche_Omega_min(PyObject *self, PyObject *args, PyObject *keywds) { +static PyObject *roche_Omega_min([[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { auto fname = "roche_Omega_min"_s; @@ -868,7 +868,7 @@ static PyObject *roche_Omega_min(PyObject *self, PyObject *args, PyObject *keywd compact primary Roche lobe: Omega_{min} (q, F, d, misalignment) = - min {Omega(L1(q,F,d,misaligment)), Omega(L2(q,F,d,misalignment)) } + min {Omega(L1(q,F,d,misalignment)), Omega(L2(q,F,d,misalignment)) } Python: @@ -888,7 +888,7 @@ static PyObject *roche_Omega_min(PyObject *self, PyObject *args, PyObject *keywd Omega0 - value of the Omega at (x,y,z) */ -static PyObject *roche_misaligned_Omega_min(PyObject *self, PyObject *args, PyObject *keywds) { +static PyObject *roche_misaligned_Omega_min([[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { auto fname = "roche_misaligned_Omega_min"_s; @@ -972,7 +972,7 @@ static PyObject *roche_misaligned_Omega_min(PyObject *self, PyObject *args, PyOb h : height of the lobe's pole */ -static PyObject *rotstar_pole(PyObject *self, PyObject *args, PyObject *keywds) { +static PyObject *rotstar_pole([[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { // // Reading arguments @@ -1039,7 +1039,7 @@ static PyObject *rotstar_pole(PyObject *self, PyObject *args, PyObject *keywds) h : height of the lobe's pole */ -static PyObject *rotstar_misaligned_pole(PyObject *self, PyObject *args, PyObject *keywds) { +static PyObject *rotstar_misaligned_pole([[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { auto fname = "rotstar_misaligned_pole"_s; @@ -1094,7 +1094,7 @@ static PyObject *rotstar_misaligned_pole(PyObject *self, PyObject *args, PyObjec h : height of the lobe's pole = R */ -static PyObject *sphere_pole(PyObject *self, PyObject *args, PyObject *keywds) { +static PyObject *sphere_pole([[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { auto fname = "sphere_pole"_s; @@ -1149,7 +1149,7 @@ static PyObject *sphere_pole(PyObject *self, PyObject *args, PyObject *keywds) { param_rotstar : 1-rank numpy array = (omega_rotstar, Omega0_rotstar) */ -static PyObject *rotstar_from_roche(PyObject *self, PyObject *args, PyObject *keywds) { +static PyObject *rotstar_from_roche([[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { auto fname = "rotstar_from_roche"_s; // @@ -1242,7 +1242,7 @@ static PyObject *rotstar_from_roche(PyObject *self, PyObject *args, PyObject *ke */ static PyObject *rotstar_misaligned_from_roche_misaligned( - PyObject *self, PyObject *args, PyObject *keywds) { + [[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { auto fname = "rotstar_misaligned_from_roche_misaligned"_s; @@ -1371,7 +1371,7 @@ static PyObject *rotstar_misaligned_from_roche_misaligned( */ -static PyObject *roche_area_volume(PyObject *self, PyObject *args, PyObject *keywds) { +static PyObject *roche_area_volume([[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { auto fname = "roche_area_volume"_s; @@ -1442,7 +1442,7 @@ static PyObject *roche_area_volume(PyObject *self, PyObject *args, PyObject *key double b = (1 + q)*F*F*delta*delta*delta, w0 = 5*(q + std::cbrt(b*b)/4) - 29.554 - 5.26235*std::log(std::min(eps[0], eps[1])), - av[2]; // for results + av[3]; // for results, av[2] not used here if (choice == 0 && delta*Omega0 >= std::max(10., w0)) { @@ -1482,7 +1482,7 @@ static PyObject *roche_area_volume(PyObject *self, PyObject *args, PyObject *key polish = false, adjust = true; - double p[2][2], e, t; + double p[2][3], e, t; // p[*][2] not used here // // one step adjustment of precison for area and volume @@ -1601,7 +1601,7 @@ static PyObject *roche_area_volume(PyObject *self, PyObject *args, PyObject *key float: */ -static PyObject *rotstar_area_volume(PyObject *self, PyObject *args, PyObject *keywds) { +static PyObject *rotstar_area_volume([[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { auto fname = "rotstar_area_volume"_s; @@ -1729,7 +1729,7 @@ static PyObject *rotstar_area_volume(PyObject *self, PyObject *args, PyObject *k float: */ -static PyObject *rotstar_misaligned_area_volume(PyObject *self, PyObject *args, PyObject *keywds) { +static PyObject *rotstar_misaligned_area_volume([[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { auto fname = "rotstar_misaligned_area_volume"_s; @@ -1848,7 +1848,7 @@ static PyObject *rotstar_misaligned_area_volume(PyObject *self, PyObject *args, float: */ -static PyObject *sphere_area_volume(PyObject *self, PyObject *args, PyObject *keywds) { +static PyObject *sphere_area_volume([[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { auto fname = "sphere_area_volume"_s; @@ -1941,7 +1941,7 @@ static PyObject *sphere_area_volume(PyObject *self, PyObject *args, PyObject *ke */ -static PyObject *roche_misaligned_critical_volume(PyObject *self, PyObject *args, PyObject *keywds) { +static PyObject *roche_misaligned_critical_volume([[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { auto fname = "roche_misaligned_critical_volume"_s; @@ -2017,7 +2017,7 @@ static PyObject *roche_misaligned_critical_volume(PyObject *self, PyObject *args /* C++ wrapper for Python code: - Calculate area and volume of the generalied Roche lobes with + Calculate area and volume of the generalized Roche lobes with misaligned spin and orbit velocity vectors. Omega_0 = Omega(x,y,z) @@ -2067,7 +2067,7 @@ static PyObject *roche_misaligned_critical_volume(PyObject *self, PyObject *args */ -static PyObject *roche_misaligned_area_volume(PyObject *self, PyObject *args, PyObject *keywds) { +static PyObject *roche_misaligned_area_volume([[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { auto fname = "roche_misaligned_area_volume"_s; @@ -2209,7 +2209,8 @@ static PyObject *roche_misaligned_area_volume(PyObject *self, PyObject *args, Py bool adjust = true; - double p[2][2], xrange[2], pole, e; + double p[2][2], xrange[2], e, + pole = -1; // set just for control // // Choosing boundaries on x-axis or calculating the pole @@ -2235,7 +2236,7 @@ static PyObject *roche_misaligned_area_volume(PyObject *self, PyObject *args, Py do { for (int i = 0, m = m0; i < 2; ++i, m <<= 1) - if (theta == 0) + if (aligned) gen_roche::area_volume_integration (p[i], res_choice, xrange, Omega0, q, F, d, m); else { @@ -2345,7 +2346,7 @@ static PyObject *roche_misaligned_area_volume(PyObject *self, PyObject *args, Py */ -static PyObject *roche_Omega_at_vol(PyObject *self, PyObject *args, PyObject *keywds) { +static PyObject *roche_Omega_at_vol([[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { auto fname = "roche_Omega_at_vol"_s; @@ -2505,7 +2506,7 @@ static PyObject *roche_Omega_at_vol(PyObject *self, PyObject *args, PyObject *ke return NULL; } // We use the condition on the argument (= Omega) ~ constraining backward error, - // but we could also use condition on the value (= Volume) ~ constraing forward error + // but we could also use condition on the value (= Volume) ~ constraining forward error return PyFloat_FromDouble(Omega); } @@ -2539,7 +2540,7 @@ static PyObject *roche_Omega_at_vol(PyObject *self, PyObject *args, PyObject *ke */ -static PyObject *rotstar_Omega_at_vol(PyObject *self, PyObject *args, PyObject *keywds) { +static PyObject *rotstar_Omega_at_vol([[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { auto fname = "rotstar_Omega_at_vol"_s; @@ -2620,7 +2621,7 @@ static PyObject *rotstar_Omega_at_vol(PyObject *self, PyObject *args, PyObject * */ -static PyObject *rotstar_misaligned_Omega_at_vol(PyObject *self, PyObject *args, PyObject *keywds) { +static PyObject *rotstar_misaligned_Omega_at_vol([[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { auto fname = "rotstar_misaligned_Omega_at_vol"_s; @@ -2708,7 +2709,7 @@ static PyObject *rotstar_misaligned_Omega_at_vol(PyObject *self, PyObject *args, */ -static PyObject *roche_misaligned_Omega_at_vol(PyObject *self, PyObject *args, PyObject *keywds) { +static PyObject *roche_misaligned_Omega_at_vol([[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { auto fname = "roche_misaligned_Omega_at_vol"_s; @@ -3034,7 +3035,7 @@ static PyObject *roche_misaligned_Omega_at_vol(PyObject *self, PyObject *args, P */ -static PyObject *sphere_Omega_at_vol(PyObject *self, PyObject *args, PyObject *keywds) { +static PyObject *sphere_Omega_at_vol([[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { auto fname = "sphere_Omega_at_vol"_s; @@ -3093,7 +3094,7 @@ static PyObject *sphere_Omega_at_vol(PyObject *self, PyObject *args, PyObject *k */ -static PyObject *roche_gradOmega(PyObject *self, PyObject *args) { +static PyObject *roche_gradOmega([[maybe_unused]] PyObject *self, PyObject *args) { double p[4]; @@ -3124,7 +3125,7 @@ static PyObject *roche_gradOmega(PyObject *self, PyObject *args) { /* C++ wrapper for Python code: - Calculate the gradient and the value of the rotenting star potential + Calculate the gradient and the value of the rotating star potential at a given point -grad Omega (x,y,z) @@ -3148,7 +3149,7 @@ static PyObject *roche_gradOmega(PyObject *self, PyObject *args) { */ -static PyObject *rotstar_gradOmega(PyObject *self, PyObject *args) { +static PyObject *rotstar_gradOmega([[maybe_unused]] PyObject *self, PyObject *args) { auto fname = "rotstar_gradOmega"_s; @@ -3225,7 +3226,7 @@ static PyObject *rotstar_gradOmega(PyObject *self, PyObject *args) { */ -static PyObject *rotstar_misaligned_gradOmega(PyObject *self, PyObject *args) { +static PyObject *rotstar_misaligned_gradOmega([[maybe_unused]] PyObject *self, PyObject *args) { auto fname = "rotstar_misaligned_gradOmega"_s; @@ -3295,7 +3296,7 @@ static PyObject *rotstar_misaligned_gradOmega(PyObject *self, PyObject *args) { */ -static PyObject *sphere_gradOmega(PyObject *self, PyObject *args) { +static PyObject *sphere_gradOmega([[maybe_unused]] PyObject *self, PyObject *args) { auto fname = "sphere_gradOmega"_s; @@ -3361,7 +3362,7 @@ static PyObject *sphere_gradOmega(PyObject *self, PyObject *args) { */ -static PyObject *roche_misaligned_gradOmega(PyObject *self, PyObject *args) { +static PyObject *roche_misaligned_gradOmega([[maybe_unused]] PyObject *self, PyObject *args) { auto fname = "roche_misaligned_gradOmega"_s; @@ -3443,7 +3444,7 @@ static PyObject *roche_misaligned_gradOmega(PyObject *self, PyObject *args) { g : 1-rank numpy array = -grad Omega (x,y,z) */ -static PyObject *roche_gradOmega_only(PyObject *self, PyObject *args) { +static PyObject *roche_gradOmega_only([[maybe_unused]] PyObject *self, PyObject *args) { auto fname = "roche_gradOmega_only"_s; @@ -3494,7 +3495,7 @@ static PyObject *roche_gradOmega_only(PyObject *self, PyObject *args) { g : 1-rank numpy array = -grad Omega (x,y,z) */ -static PyObject *rotstar_gradOmega_only(PyObject *self, PyObject *args) { +static PyObject *rotstar_gradOmega_only([[maybe_unused]] PyObject *self, PyObject *args) { auto fname = "rotstar_gradOmega_only"_s; @@ -3569,7 +3570,7 @@ static PyObject *rotstar_gradOmega_only(PyObject *self, PyObject *args) { g : 1-rank numpy array = -grad Omega (x,y,z) */ -static PyObject *rotstar_misaligned_gradOmega_only(PyObject *self, PyObject *args) { +static PyObject *rotstar_misaligned_gradOmega_only([[maybe_unused]] PyObject *self, PyObject *args) { auto fname = "rotstar_misaligned_gradOmega_only"_s; @@ -3639,7 +3640,7 @@ static PyObject *rotstar_misaligned_gradOmega_only(PyObject *self, PyObject *arg g : 1-rank numpy array = -grad Omega (x,y,z) */ -static PyObject *sphere_gradOmega_only(PyObject *self, PyObject *args) { +static PyObject *sphere_gradOmega_only([[maybe_unused]] PyObject *self, PyObject *args) { auto fname = "sphere_gradOmega_only"_s; @@ -3699,7 +3700,7 @@ static PyObject *sphere_gradOmega_only(PyObject *self, PyObject *args) { g : 1-rank numpy array = -grad Omega (x,y,z) */ -static PyObject *roche_misaligned_gradOmega_only(PyObject *self, PyObject *args) { +static PyObject *roche_misaligned_gradOmega_only([[maybe_unused]] PyObject *self, PyObject *args) { auto fname = "roche_misaligned_gradOmega_only"_s; @@ -3781,7 +3782,7 @@ static PyObject *roche_misaligned_gradOmega_only(PyObject *self, PyObject *args) Omega0 - value of the Omega at (x,y,z) */ -static PyObject *roche_Omega(PyObject *self, PyObject *args) { +static PyObject *roche_Omega([[maybe_unused]] PyObject *self, PyObject *args) { double p[4]; @@ -3822,7 +3823,7 @@ static PyObject *roche_Omega(PyObject *self, PyObject *args) { Omega0 - value of the Omega at (x,y,z) */ -static PyObject *rotstar_Omega(PyObject *self, PyObject *args) { +static PyObject *rotstar_Omega([[maybe_unused]] PyObject *self, PyObject *args) { double p[2]; @@ -3876,7 +3877,7 @@ static PyObject *rotstar_Omega(PyObject *self, PyObject *args) { Omega0 - value of the Omega at (x,y,z) */ -static PyObject *rotstar_misaligned_Omega(PyObject *self, PyObject *args) { +static PyObject *rotstar_misaligned_Omega([[maybe_unused]] PyObject *self, PyObject *args) { auto fname = "rotstar_misaligned_Omega"_s; @@ -3941,7 +3942,7 @@ static PyObject *rotstar_misaligned_Omega(PyObject *self, PyObject *args) { Omega0 - value of the Omega at (x,y,z) */ -static PyObject *sphere_Omega(PyObject *self, PyObject *args) { +static PyObject *sphere_Omega([[maybe_unused]] PyObject *self, PyObject *args) { auto fname = "sphere_Omega"_s; @@ -3987,7 +3988,7 @@ static PyObject *sphere_Omega(PyObject *self, PyObject *args) { */ -static PyObject *roche_misaligned_Omega(PyObject *self, PyObject *args) { +static PyObject *roche_misaligned_Omega([[maybe_unused]] PyObject *self, PyObject *args) { auto fname = "roche_misaligned_Omega"_s; @@ -4157,7 +4158,7 @@ static PyObject *roche_misaligned_Omega(PyObject *self, PyObject *args) { */ -static PyObject *roche_marching_mesh(PyObject *self, PyObject *args, PyObject *keywds) { +static PyObject *roche_marching_mesh([[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { auto fname = "roche_marching_mesh"_s; @@ -4519,7 +4520,7 @@ static PyObject *roche_marching_mesh(PyObject *self, PyObject *args, PyObject *k * https://docs.python.org/2/c-api/arg.html#c.PyArg_ParseTupleAndKeywords */ -static PyObject *rotstar_marching_mesh(PyObject *self, PyObject *args, PyObject *keywds) { +static PyObject *rotstar_marching_mesh([[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { auto fname = "rotstar_marching_mesh"_s; @@ -4913,7 +4914,7 @@ static PyObject *rotstar_marching_mesh(PyObject *self, PyObject *args, PyObject * https://docs.python.org/2/c-api/arg.html#c.PyArg_ParseTupleAndKeywords */ -static PyObject *rotstar_misaligned_marching_mesh(PyObject *self, PyObject *args, PyObject *keywds) { +static PyObject *rotstar_misaligned_marching_mesh([[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { auto fname = "rotstar_misaligned_marching_mesh"_s; @@ -5300,7 +5301,7 @@ static PyObject *rotstar_misaligned_marching_mesh(PyObject *self, PyObject *args * https://docs.python.org/2/c-api/arg.html#c.PyArg_ParseTupleAndKeywords */ -static PyObject *sphere_marching_mesh(PyObject *self, PyObject *args, PyObject *keywds) { +static PyObject *sphere_marching_mesh([[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { auto fname = "sphere_marching_mesh"_s; @@ -5683,7 +5684,7 @@ static PyObject *sphere_marching_mesh(PyObject *self, PyObject *args, PyObject * */ -static PyObject *roche_misaligned_marching_mesh(PyObject *self, PyObject *args, PyObject *keywds) { +static PyObject *roche_misaligned_marching_mesh([[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { auto fname = "roche_misaligned_marching_mesh"_s; @@ -5807,9 +5808,11 @@ static PyObject *roche_misaligned_marching_mesh(PyObject *self, PyObject *args, // // Choosing the meshing initial point // - bool rotated, ok, aligned = false; + bool + rotated, ok, aligned = false; - double r[3], g[3], theta, *s = 0; + double r[3], g[3], theta, *s = 0, + eps = 2*std::numeric_limits::epsilon(); if (PyFloat_Check(o_misalignment)) { @@ -5823,7 +5826,7 @@ static PyObject *roche_misaligned_marching_mesh(PyObject *self, PyObject *args, PyArray_TYPE((PyArrayObject *) o_misalignment) == NPY_DOUBLE) { s = (double*) PyArray_DATA((PyArrayObject*)o_misalignment); - aligned = (s[0] == 0 && s[1] == 0); + aligned = (std::abs(s[0]) < eps && std::abs(s[1]) < eps) || (std::abs(1 - s[2]) < eps); // we could work with s[0]==0, calculate aligned case make simple // rotation around x-axis @@ -5857,6 +5860,10 @@ static PyObject *roche_misaligned_marching_mesh(PyObject *self, PyObject *args, return NULL; } + if (verbosity_level>=4) + report_stream << fname + << "::aligned=" << aligned + << " rotated =" << rotated << '\n'; // // Marching triangulation of the Roche lobe and calculate central points // @@ -6061,7 +6068,7 @@ static PyObject *roche_misaligned_marching_mesh(PyObject *self, PyObject *args, * http://folk.uio.no/hpl/scripting/doc/python/NumPy/Numeric/numpy-13.html */ -static PyObject *mesh_visibility(PyObject *self, PyObject *args, PyObject *keywds){ +static PyObject *mesh_visibility([[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds){ auto fname = "mesh_visibility"_s; @@ -6222,7 +6229,7 @@ static PyObject *mesh_visibility(PyObject *self, PyObject *args, PyObject *keywd * http://folk.uio.no/hpl/scripting/doc/python/NumPy/Numeric/numpy-13.html */ -static PyObject *mesh_rough_visibility(PyObject *self, PyObject *args){ +static PyObject *mesh_rough_visibility([[maybe_unused]] PyObject *self, PyObject *args){ // // Storing/Reading arguments @@ -6329,7 +6336,7 @@ static PyObject *mesh_rough_visibility(PyObject *self, PyObject *args){ area: float - area of triangles of mesh */ -static PyObject *mesh_offseting(PyObject *self, PyObject *args, PyObject *keywds){ +static PyObject *mesh_offseting([[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds){ auto fname = "mesh_offseting"_s; @@ -6508,7 +6515,7 @@ static PyObject *mesh_offseting(PyObject *self, PyObject *args, PyObject *keywd volume - volume of body enclosed by triangular mesh */ -static PyObject *mesh_properties(PyObject *self, PyObject *args, PyObject *keywds) { +static PyObject *mesh_properties([[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { // // Reading arguments @@ -6659,7 +6666,7 @@ static PyObject *mesh_properties(PyObject *self, PyObject *args, PyObject *keywd */ -static PyObject *mesh_export_povray(PyObject *self, PyObject *args, PyObject *keywds) { +static PyObject *mesh_export_povray([[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { // // Reading arguments @@ -6938,7 +6945,7 @@ bool LDmodelFromListOfTuples( */ static PyObject *mesh_radiosity_problem( - PyObject *self, PyObject *args, PyObject *keywds) { + [[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { auto fname = "mesh_radiosity_problem"_s; @@ -7162,7 +7169,7 @@ static PyObject *mesh_radiosity_problem( */ static PyObject *mesh_radiosity_problem_nbody_convex( - PyObject *self, PyObject *args, PyObject *keywds) { + [[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { auto fname = "mesh_radiosity_problem_nbody_convex"_s; @@ -7454,7 +7461,7 @@ struct Tmesh_radiosity_redistrib_problem_nbody { } __redistrib_problem_nbody; static PyObject *mesh_radiosity_redistrib_problem_nbody_convex_setup( - PyObject *self, PyObject *args, PyObject *keywds){ + [[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds){ auto fname = "mesh_radiosity_redistrib_problem_nbody_convex_setup"_s; @@ -7491,7 +7498,7 @@ static PyObject *mesh_radiosity_redistrib_problem_nbody_convex_setup( static PyObject *mesh_radiosity_redistrib_problem_nbody_convex( - PyObject *self, PyObject *args, PyObject *keywds) { + [[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { auto fname = "mesh_radiosity_redistrib_problem_nbody_convex"_s; @@ -7825,7 +7832,7 @@ static PyObject *mesh_radiosity_redistrib_problem_nbody_convex( {'update-emittanceB': 2.0410763114593298, 'update-emittanceA': 1.0206982972948087, 'radiosityB': 2.012322893437799, 'radiosityA': 1.014488808106366} */ -static PyObject *radiosity_redistrib_1dmodel(PyObject *self, PyObject *args, PyObject *keywds) { +static PyObject *radiosity_redistrib_1dmodel([[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { auto fname = "radiosity_redistrib_1dmodel"_s; @@ -7958,7 +7965,7 @@ static PyObject *radiosity_redistrib_1dmodel(PyObject *self, PyObject *args, PyO cnormgrads: default false GatC: 1-rank numpy array of norms of gradients at centers */ -static PyObject *roche_central_points(PyObject *self, PyObject *args, PyObject *keywds){ +static PyObject *roche_central_points([[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds){ auto fname = "roche_central_points"_s; @@ -8118,7 +8125,7 @@ static PyObject *roche_central_points(PyObject *self, PyObject *args, PyObject * http://folk.uio.no/hpl/scripting/doc/python/NumPy/Numeric/numpy-13.html */ -static PyObject *roche_reprojecting_vertices(PyObject *self, PyObject *args, PyObject *keywds) { +static PyObject *roche_reprojecting_vertices([[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { // // Reading arguments @@ -8251,7 +8258,7 @@ static PyObject *roche_reprojecting_vertices(PyObject *self, PyObject *args, PyO H: 2-rank numpy array of 3D point on a horizon */ -static PyObject *roche_horizon(PyObject *self, PyObject *args, PyObject *keywds) { +static PyObject *roche_horizon([[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { auto fname = "roche_horizon"_s; @@ -8354,7 +8361,7 @@ static PyObject *roche_horizon(PyObject *self, PyObject *args, PyObject *keywds) H: 2-rank numpy array of floats -- 3D points on a horizon */ -static PyObject *rotstar_horizon(PyObject *self, PyObject *args, PyObject *keywds) { +static PyObject *rotstar_horizon([[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { auto fname = "rotstar_horizon"_s; @@ -8466,7 +8473,7 @@ static PyObject *rotstar_horizon(PyObject *self, PyObject *args, PyObject *keywd H: 2-rank numpy array of floats -- 3D points on a horizon */ -static PyObject *rotstar_misaligned_horizon(PyObject *self, PyObject *args, PyObject *keywds) { +static PyObject *rotstar_misaligned_horizon([[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { auto fname = "rotstar_misaligned_horizon"_s; @@ -8588,7 +8595,7 @@ static PyObject *rotstar_misaligned_horizon(PyObject *self, PyObject *args, PyOb H: 2-rank numpy array of 3D point on a horizon */ -static PyObject *roche_misaligned_horizon(PyObject *self, PyObject *args, PyObject *keywds) { +static PyObject *roche_misaligned_horizon([[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { auto fname = "roche_misaligned_horizon"_s; @@ -8731,7 +8738,7 @@ static PyObject *roche_misaligned_horizon(PyObject *self, PyObject *args, PyObje xrange: 1-rank numpy array of two numbers p */ -static PyObject *roche_xrange(PyObject *self, PyObject *args, PyObject *keywds) { +static PyObject *roche_xrange([[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { auto fname = "roche_xrange"_s; @@ -8840,7 +8847,7 @@ static PyObject *roche_xrange(PyObject *self, PyObject *args, PyObject *keywds) */ -static PyObject *roche_square_grid(PyObject *self, PyObject *args, PyObject *keywds) { +static PyObject *roche_square_grid([[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { auto fname = "roche_square_grid"_s; @@ -9229,7 +9236,7 @@ static PyObject *roche_square_grid(PyObject *self, PyObject *args, PyObject *key value of D(mu) for a given LD model */ -static PyObject *ld_D(PyObject *self, PyObject *args, PyObject *keywds) { +static PyObject *ld_D([[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { auto fname = "ld_D"_s; @@ -9299,7 +9306,7 @@ static PyObject *ld_D(PyObject *self, PyObject *args, PyObject *keywds) { value of integrated D(mu) for a given LD model */ -static PyObject *ld_D0(PyObject *self, PyObject *args, PyObject *keywds) { +static PyObject *ld_D0([[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { auto fname = "ld_D0"_s; @@ -9368,7 +9375,7 @@ static PyObject *ld_D0(PyObject *self, PyObject *args, PyObject *keywds) { 1-rank numpy array of floats: gradient of the function D(mu) w.r.t. parameters */ -static PyObject *ld_gradparD(PyObject *self, PyObject *args, PyObject *keywds) { +static PyObject *ld_gradparD([[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { auto fname = "ld_gradparD"_s; @@ -9449,7 +9456,7 @@ static PyObject *ld_gradparD(PyObject *self, PyObject *args, PyObject *keywds) { int: number of parameters */ -static PyObject *ld_nrpar(PyObject *self, PyObject *args, PyObject *keywds) { +static PyObject *ld_nrpar([[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { auto fname = "ld_nrpar"_s; @@ -9532,7 +9539,7 @@ static PyObject *ld_nrpar(PyObject *self, PyObject *args, PyObject *keywds) { */ -static PyObject *ld_check(PyObject *self, PyObject *args, PyObject *keywds) { +static PyObject *ld_check([[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { auto fname = "ld_check"_s; @@ -9601,7 +9608,7 @@ static PyObject *ld_check(PyObject *self, PyObject *args, PyObject *keywds) { 1-rank numpy array of floats */ -static PyObject *wd_readdata(PyObject *self, PyObject *args, PyObject *keywds) { +static PyObject *wd_readdata([[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { auto fname ="wd_readdata"_s; // // Reading arguments @@ -9651,28 +9658,28 @@ static PyObject *wd_readdata(PyObject *self, PyObject *args, PyObject *keywds) { // // Reading // - + int len[2] = { wd_atm::read_data(PyString_AsString(ofilename_planck), planck_table), wd_atm::read_data(PyString_AsString(ofilename_atm), atm_table) }; - + // // Checks // std::string err_msg; - + if (len[0] < 0) err_msg = "\nProblem opening the planck file:"_s + PyString_AsString(ofilename_planck); else if (len[0] != wd_atm::N_planck) err_msg = "\nWrong size read, len= "_s + std::to_string(len[0]) + " len_expected="_s + std::to_string(wd_atm::N_planck); - + if (len[1] < 0) err_msg += "\nProblem opening the atm file:"_s + PyString_AsString(ofilename_atm); else if (len[1] != wd_atm::N_atm) err_msg += "\nWrong size read, len= "_s + std::to_string(len[1]) + " len_expected="_s + std::to_string(wd_atm::N_atm); - - + + if (err_msg.size() != 0) { raise_exception(fname + "::Problem reading data." + err_msg); delete [] planck_table; @@ -9716,7 +9723,7 @@ static PyObject *wd_readdata(PyObject *self, PyObject *args, PyObject *keywds) { y: float - Planck central intensity */ -static PyObject *wd_planckint(PyObject *self, PyObject *args, PyObject *keywds) { +static PyObject *wd_planckint([[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { // // Reading arguments @@ -9803,7 +9810,7 @@ static PyObject *wd_planckint(PyObject *self, PyObject *args, PyObject *keywds) Note: In the case of errors in calculations ylog/entry in numpy array is NaN. */ -static PyObject *wd_planckint(PyObject *self, PyObject *args, PyObject *keywds) { +static PyObject *wd_planckint([[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { auto fname = "wd_planckint"_s; @@ -9935,7 +9942,7 @@ static PyObject *wd_planckint(PyObject *self, PyObject *args, PyObject *keywds) xint - intensity abunin - the allowed value nearest to the input value. */ -static PyObject *wd_atmint(PyObject *self, PyObject *args, PyObject *keywds) { +static PyObject *wd_atmint([[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { // // Reading arguments @@ -10034,7 +10041,7 @@ static PyObject *wd_atmint(PyObject *self, PyObject *args, PyObject *keywds) { xintlog - log of intensity abunin - allowed value nearest to the input value. */ -static PyObject *wd_atmint(PyObject *self, PyObject *args, PyObject *keywds) { +static PyObject *wd_atmint([[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { auto fname = "wd_atmint"_s; @@ -10274,7 +10281,7 @@ static PyObject *wd_atmint(PyObject *self, PyObject *args, PyObject *keywds) { 2-rank numpy array = MxNv array of interpolated values */ -static PyObject *interp(PyObject *self, PyObject *args, PyObject *keywds) { +static PyObject *interp([[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { char *kwlist[] = { (char*)"req", (char*)"axes", @@ -10414,7 +10421,7 @@ static PyObject *interp(PyObject *self, PyObject *args, PyObject *keywds) { Return: r: 1- rank numpy array */ -static PyObject *scalproj_cosangle(PyObject *self, PyObject *args) { +static PyObject *scalproj_cosangle([[maybe_unused]] PyObject *self, PyObject *args) { auto fname = "vec_proj"_s; @@ -10534,7 +10541,7 @@ static PyObject *scalproj_cosangle(PyObject *self, PyObject *args) { */ -static PyObject *roche_contact_partial_area_volume(PyObject *self, PyObject *args, PyObject *keywds) { +static PyObject *roche_contact_partial_area_volume([[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { auto fname = "roche_contact_partial_area_volume"_s; @@ -10775,7 +10782,7 @@ static PyObject *roche_contact_partial_area_volume(PyObject *self, PyObject *arg */ -static PyObject *roche_contact_neck_min(PyObject *self, PyObject *args, PyObject *keywds) { +static PyObject *roche_contact_neck_min([[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { auto fname = "roche_contact_neck_min"_s; @@ -10856,7 +10863,7 @@ static PyObject *roche_contact_neck_min(PyObject *self, PyObject *args, PyObject value of the Kopal potential for (q,F,d) at which the lobe has the given volume */ -static PyObject *roche_contact_Omega_at_partial_vol(PyObject *self, PyObject *args, PyObject *keywds) { +static PyObject *roche_contact_Omega_at_partial_vol([[maybe_unused]] PyObject *self, PyObject *args, PyObject *keywds) { auto fname = "roche_contact_Omega_at_partial_vol"_s; @@ -11162,7 +11169,7 @@ static PyObject *roche_contact_Omega_at_partial_vol(PyObject *self, PyObject *ar or 2-rank numpy array : if lam and Teff are 1-rank numpy arrays */ -static PyObject *planck_function(PyObject *self, PyObject *args) { +static PyObject *planck_function([[maybe_unused]] PyObject *self, PyObject *args) { const double A = 1.1910429526245747e-16; // = 2 h c^2 [m4 kg / s3]; const double B = 0.014387773538277205; // = hc/k [mK]; @@ -11275,7 +11282,7 @@ static PyObject *planck_function(PyObject *self, PyObject *args) { or 2-rank numpy array: array of two values */ -static PyObject *CCM89_extinction(PyObject *self, PyObject *args) { +static PyObject *CCM89_extinction([[maybe_unused]] PyObject *self, PyObject *args) { const char *fname = "CCM89_extinction"; @@ -11384,7 +11391,7 @@ static PyObject *CCM89_extinction(PyObject *self, PyObject *args) { or 2-rank numpy array: array of two values */ -static PyObject *gordon_extinction(PyObject *self, PyObject *args) { +static PyObject *gordon_extinction([[maybe_unused]] PyObject *self, PyObject *args) { const char *fname = "gordon_extinction"; @@ -11479,10 +11486,13 @@ static PyObject *gordon_extinction(PyObject *self, PyObject *args) { Ref: * https://docs.python.org/2.0/ext/parseTupleAndKeywords.html */ + +#define O2F (PyCFunction)(void(*)(void)) + static PyMethodDef Methods[] = { { "roche_misaligned_transf", - roche_misaligned_transf, + O2F roche_misaligned_transf, METH_VARARGS, "Determine angle parameters of the misaligned Roche lobes from " "the spin angular velocity in the rotating binary system "}, @@ -11491,19 +11501,19 @@ static PyMethodDef Methods[] = { { "roche_critical_potential", - (PyCFunction)roche_critical_potential, + O2F roche_critical_potential, METH_VARARGS|METH_KEYWORDS, "Determine the critical potentials of Kopal potential for given " "values of q, F, and d."}, { "rotstar_critical_potential", - rotstar_critical_potential, + O2F rotstar_critical_potential, METH_VARARGS, "Determine the critical potentials of the rotating star potental " "for given values of omega."}, { "rotstar_misaligned_critical_potential", - rotstar_misaligned_critical_potential, + O2F rotstar_misaligned_critical_potential, METH_VARARGS, "Determine the critical potentials of the rotating star potental " "with misalignment for given values of omega and spin"}, @@ -11511,29 +11521,29 @@ static PyMethodDef Methods[] = { // -------------------------------------------------------------------- { "roche_pole", - (PyCFunction)roche_pole, + O2F roche_pole, METH_VARARGS|METH_KEYWORDS, "Determine the height of the pole of generalized Roche lobes for given " "values of q, F, d and Omega0"}, { "rotstar_pole", - (PyCFunction)rotstar_pole, + O2F rotstar_pole, METH_VARARGS|METH_KEYWORDS, "Determine the height of the pole of rotating star for given a omega."}, { "rotstar_misaligned_pole", - (PyCFunction)rotstar_misaligned_pole, + O2F rotstar_misaligned_pole, METH_VARARGS|METH_KEYWORDS, "Determine the height of the pole of rotating star with misalignment " "for given a omega and spin."}, { "sphere_pole", - (PyCFunction)sphere_pole, + O2F sphere_pole, METH_VARARGS|METH_KEYWORDS, "Determine the height of the pole of sphere for given a R."}, { "roche_misaligned_pole", - (PyCFunction)roche_misaligned_pole, + O2F roche_misaligned_pole, METH_VARARGS|METH_KEYWORDS, "Determine the postion of the pole of generalized Roche lobes with " "misaligned angular spin-orbital angular velocity vectors for given " @@ -11542,14 +11552,14 @@ static PyMethodDef Methods[] = { // -------------------------------------------------------------------- {"roche_Omega_min", - (PyCFunction)roche_Omega_min, + O2F roche_Omega_min, METH_VARARGS|METH_KEYWORDS, "Determine the minimal posible value of the Kopal potential that" "permits existance of the compact Roche lobe for given " "values of q, F and d."}, { "roche_misaligned_Omega_min", - (PyCFunction)roche_misaligned_Omega_min, + O2F roche_misaligned_Omega_min, METH_VARARGS|METH_KEYWORDS, "Determine the minimal posible value of the Kopal potential that" "permits existance of the compact Roche lobe for given " @@ -11557,7 +11567,7 @@ static PyMethodDef Methods[] = { // -------------------------------------------------------------------- { "roche_misaligned_critical_volume", - (PyCFunction)roche_misaligned_critical_volume, + O2F roche_misaligned_critical_volume, METH_VARARGS|METH_KEYWORDS, "Determine the volume of the semi-detached case of the misaligned " "Roche lobe for given values of q, F, F and misalignment (theta or " @@ -11565,14 +11575,14 @@ static PyMethodDef Methods[] = { // -------------------------------------------------------------------- { "rotstar_from_roche", - (PyCFunction)rotstar_from_roche, + O2F rotstar_from_roche, METH_VARARGS|METH_KEYWORDS, "Determine parameters of the rotating stars from parameters Roche " " by matching the poles"}, { "rotstar_misaligned_from_roche_misaligned", - (PyCFunction)rotstar_misaligned_from_roche_misaligned, + O2F rotstar_misaligned_from_roche_misaligned, METH_VARARGS|METH_KEYWORDS, "Determine parameters of the rotating stars with misalignment from " "parameters Roche with misalignment by matching the poles."}, @@ -11580,30 +11590,30 @@ static PyMethodDef Methods[] = { // -------------------------------------------------------------------- { "roche_area_volume", - (PyCFunction)roche_area_volume, + O2F roche_area_volume, METH_VARARGS|METH_KEYWORDS, "Determine the area and volume of the generalized Roche lobes for given " "values of q, F, d and Omega0."}, { "rotstar_area_volume", - (PyCFunction)rotstar_area_volume, + O2F rotstar_area_volume, METH_VARARGS|METH_KEYWORDS, "Determine the area and volume of the rotating star for given a omega " "and Omega0"}, { "rotstar_misaligned_area_volume", - (PyCFunction)rotstar_misaligned_area_volume, + O2F rotstar_misaligned_area_volume, METH_VARARGS|METH_KEYWORDS, "Determine the area and volume of the rotating star with misalignment " "for given a omega and Omega0"}, { "sphere_area_volume", - (PyCFunction)sphere_area_volume, + O2F sphere_area_volume, METH_VARARGS|METH_KEYWORDS, "Determine the area and volume of the sphere for given a R."}, { "roche_misaligned_area_volume", - (PyCFunction)roche_misaligned_area_volume, + O2F roche_misaligned_area_volume, METH_VARARGS|METH_KEYWORDS, "Determine the area and volume of the generalized Roche lobes with " "misaligned spin and orbtal angular velocity vectors for given " @@ -11612,25 +11622,25 @@ static PyMethodDef Methods[] = { // -------------------------------------------------------------------- { "roche_Omega_at_vol", - (PyCFunction)roche_Omega_at_vol, + O2F roche_Omega_at_vol, METH_VARARGS|METH_KEYWORDS, "Determine the value of the generalized Kopal potential at " "values of q, F, d and volume."}, { "rotstar_Omega_at_vol", - (PyCFunction)rotstar_Omega_at_vol, + O2F rotstar_Omega_at_vol, METH_VARARGS|METH_KEYWORDS, "Determine the value of the rotating star potential at " "values of omega and volume."}, { "rotstar_misaligned_Omega_at_vol", - (PyCFunction)rotstar_misaligned_Omega_at_vol, + O2F rotstar_misaligned_Omega_at_vol, METH_VARARGS|METH_KEYWORDS, "Determine the value of the rotating star potential with misalignment at " "values of omega and volume."}, { "roche_misaligned_Omega_at_vol", - (PyCFunction)roche_misaligned_Omega_at_vol, + O2F roche_misaligned_Omega_at_vol, METH_VARARGS|METH_KEYWORDS, "Determine the value of the generalized Kopal potential of " "Roche lobes with with misaligned spin and orbtal angular " @@ -11638,39 +11648,39 @@ static PyMethodDef Methods[] = { "and volume."}, { "sphere_Omega_at_vol", - (PyCFunction)sphere_Omega_at_vol, + O2F sphere_Omega_at_vol, METH_VARARGS|METH_KEYWORDS, "Determine the value of the spherical star potential at given volume."}, // -------------------------------------------------------------------- { "roche_gradOmega", - roche_gradOmega, + O2F roche_gradOmega, METH_VARARGS, "Calculate the gradient and the value of the generalized Kopal potentil" " at given point [x,y,z] for given values of q, F and d."}, { "rotstar_gradOmega", - rotstar_gradOmega, + O2F rotstar_gradOmega, METH_VARARGS, "Calculate the gradient and the value of the rotating star potential" " at given point [x,y,z] for given values of omega."}, { "rotstar_misaligned_gradOmega", - rotstar_misaligned_gradOmega, + O2F rotstar_misaligned_gradOmega, METH_VARARGS, "Calculate the gradient and the value of the rotating star potential " "with misalignment at given point [x,y,z] for given values of omega " "and spin."}, { "sphere_gradOmega", - sphere_gradOmega, + O2F sphere_gradOmega, METH_VARARGS, "Calculate the gradient of the potential of the sphere" " at given point [x,y,z]."}, { "roche_misaligned_gradOmega", - roche_misaligned_gradOmega, + O2F roche_misaligned_gradOmega, METH_VARARGS, "Calculate the gradient of the generalized Kopal potential with " " misaligned angular momenta at given point [x,y,z] for given " @@ -11747,34 +11757,34 @@ static PyMethodDef Methods[] = { // -------------------------------------------------------------------- { "roche_marching_mesh", - (PyCFunction)roche_marching_mesh, + O2F roche_marching_mesh, METH_VARARGS|METH_KEYWORDS, "Determine the triangular meshing of generalized Roche lobes for " "given values of q, F, d and value of the generalized Kopal potential " "Omega0. The edge of triangles used in the mesh are approximately delta."}, { "rotstar_marching_mesh", - (PyCFunction)rotstar_marching_mesh, + O2F rotstar_marching_mesh, METH_VARARGS|METH_KEYWORDS, "Determine the triangular meshing of a rotating star for given " "values of omega and value of the star potential Omega. The edge " "of triangles used in the mesh are approximately delta."}, { "rotstar_misaligned_marching_mesh", - (PyCFunction)rotstar_misaligned_marching_mesh, + O2F rotstar_misaligned_marching_mesh, METH_VARARGS|METH_KEYWORDS, "Determine the triangular meshing of a rotating star with misalignment " "for given values of omega, spin and value of the star potential Omega." "The edge of triangles used in the mesh are approximately delta."}, { "sphere_marching_mesh", - (PyCFunction)sphere_marching_mesh, + O2F sphere_marching_mesh, METH_VARARGS|METH_KEYWORDS, "Determine the triangular meshing of a sphere for given radius R." "The edge of triangles used in the mesh are approximately delta."}, { "roche_misaligned_marching_mesh", - (PyCFunction)roche_misaligned_marching_mesh, + O2F roche_misaligned_marching_mesh, METH_VARARGS|METH_KEYWORDS, "Determine the triangular meshing of generalized Roche lobes with " "misaligned spin and orbital angular velocity vectors for " @@ -11785,7 +11795,7 @@ static PyMethodDef Methods[] = { // -------------------------------------------------------------------- { "mesh_visibility", - (PyCFunction)mesh_visibility, + O2F mesh_visibility, METH_VARARGS|METH_KEYWORDS, "Determine the ratio of triangle surfaces that are visible " "in a triangular mesh."}, @@ -11797,54 +11807,54 @@ static PyMethodDef Methods[] = { "partially hidden and visible"}, { "mesh_offseting", - (PyCFunction)mesh_offseting, + O2F mesh_offseting, METH_VARARGS|METH_KEYWORDS, "Offset the mesh along the normals in vertices to match the " "area with reference area."}, { "mesh_properties", - (PyCFunction)mesh_properties, + O2F mesh_properties, METH_VARARGS|METH_KEYWORDS, "Calculate the properties of the triangular mesh."}, { "mesh_export_povray", - (PyCFunction)mesh_export_povray, + O2F mesh_export_povray, METH_VARARGS|METH_KEYWORDS, "Exporting triangular mesh into a Pov-Ray file."}, { "mesh_radiosity_problem", - (PyCFunction)mesh_radiosity_problem, + O2F mesh_radiosity_problem, METH_VARARGS|METH_KEYWORDS, "Solving the radiosity problem with limb darkening using " "a chosen reflection model."}, { "mesh_radiosity_problem_nbody_convex", - (PyCFunction)mesh_radiosity_problem_nbody_convex, + O2F mesh_radiosity_problem_nbody_convex, METH_VARARGS|METH_KEYWORDS, "Solving the radiosity problem with limb darkening for n separate " "convex bodies using chosen reflection model."}, { "mesh_radiosity_redistrib_problem_nbody_convex", - (PyCFunction)mesh_radiosity_redistrib_problem_nbody_convex, + O2F mesh_radiosity_redistrib_problem_nbody_convex, METH_VARARGS|METH_KEYWORDS, "Solving the radiosity redistribution problem with limb darkening " "for n separate convex bodies using chosen reflection model."}, { "mesh_radiosity_redistrib_problem_nbody_convex_setup", - (PyCFunction)mesh_radiosity_redistrib_problem_nbody_convex_setup, + O2F mesh_radiosity_redistrib_problem_nbody_convex_setup, METH_VARARGS|METH_KEYWORDS, "Background setup of radiosity redistribution problem with limb " "darkening for n separate convex bodies using chosen reflection model."}, { "radiosity_redistrib_1dmodel", - (PyCFunction)radiosity_redistrib_1dmodel, + O2F radiosity_redistrib_1dmodel, METH_VARARGS|METH_KEYWORDS, "Calculating a rough approximate of the surface average updated-exitance " "and radiosity for both bodies of a binary system composed of two spheres."}, // -------------------------------------------------------------------- { "roche_reprojecting_vertices", - (PyCFunction)roche_reprojecting_vertices, + O2F roche_reprojecting_vertices, METH_VARARGS|METH_KEYWORDS, "Reprojecting vertices onto the Roche lobe defined by q,F,d, and the value of" " generalized Kopal potential Omega."}, @@ -11852,33 +11862,33 @@ static PyMethodDef Methods[] = { // -------------------------------------------------------------------- { "roche_central_points", - (PyCFunction)roche_central_points, + O2F roche_central_points, METH_VARARGS|METH_KEYWORDS, "Determining the central points of triangular mesh on the Roche lobe" " defined by q,F,d, and the value of generalized Kopal potential Omega."}, // -------------------------------------------------------------------- { "roche_horizon", - (PyCFunction)roche_horizon, + O2F roche_horizon, METH_VARARGS|METH_KEYWORDS, "Calculating the horizon on the Roche lobe defined by view direction," "q,F,d, and the value of generalized Kopal potential Omega."}, { "rotstar_horizon", - (PyCFunction)rotstar_horizon, + O2F rotstar_horizon, METH_VARARGS|METH_KEYWORDS, "Calculating the horizon on the rotating star defined by view direction," "omega, and the value of the potential"}, { "rotstar_misaligned_horizon", - (PyCFunction)rotstar_misaligned_horizon, + O2F rotstar_misaligned_horizon, METH_VARARGS|METH_KEYWORDS, "Calculating the horizon on the rotating star with misalignment " "defined by view direction, omega, spin and the value of the potential"}, { "roche_misaligned_horizon", - (PyCFunction)roche_misaligned_horizon, + O2F roche_misaligned_horizon, METH_VARARGS|METH_KEYWORDS, "Calculating the horizon on the Roche lobe with misaligned spin and orbital " "angular velocity vectors defined by the view direction," @@ -11886,61 +11896,61 @@ static PyMethodDef Methods[] = { // -------------------------------------------------------------------- { "roche_xrange", - (PyCFunction)roche_xrange, + O2F roche_xrange, METH_VARARGS|METH_KEYWORDS, "Calculating the range of the Roche lobes on x-axis at given" "q, F, d, and the value of generalized Kopal potential Omega."}, // -------------------------------------------------------------------- { "roche_square_grid", - (PyCFunction)roche_square_grid, + O2F roche_square_grid, METH_VARARGS|METH_KEYWORDS, "Calculating the square grid of the interior of the Roche lobes at given" "q, F, d, and the value of generalized Kopal potential Omega."}, // -------------------------------------------------------------------- { "ld_D", - (PyCFunction)ld_D, + O2F ld_D, METH_VARARGS|METH_KEYWORDS, "Calculating the value of the limb darkening function."}, { "ld_D0", - (PyCFunction)ld_D0, + O2F ld_D0, METH_VARARGS|METH_KEYWORDS, "Calculating the integrated limb darkening function."}, { "ld_gradparD", - (PyCFunction)ld_gradparD, + O2F ld_gradparD, METH_VARARGS|METH_KEYWORDS, "Calculating the gradient of the limb darkening function w.r.t. " "parameters."}, { "ld_nrpar", - (PyCFunction)ld_nrpar, + O2F ld_nrpar, METH_VARARGS|METH_KEYWORDS, "Returns the number of required parameters."}, { "ld_check", - (PyCFunction)ld_check, + O2F ld_check, METH_VARARGS|METH_KEYWORDS, "Checking parameters if resulting D(mu) is in the range [0,1] for all mu."}, // -------------------------------------------------------------------- { "wd_readdata", - (PyCFunction)wd_readdata, + O2F wd_readdata, METH_VARARGS|METH_KEYWORDS, "Reading the file with WD coefficients."}, { "wd_planckint", - (PyCFunction)wd_planckint, + O2F wd_planckint, METH_VARARGS|METH_KEYWORDS, "Calculating Planck central intensity at given temperatures," "filter index and array of coefficients"}, { "wd_atmint", - (PyCFunction)wd_atmint, + O2F wd_atmint, METH_VARARGS|METH_KEYWORDS, "Calculating intensity for a given atmospheres at given temperatures," "filter index and array of coefficients"}, @@ -11948,33 +11958,33 @@ static PyMethodDef Methods[] = { // -------------------------------------------------------------------- {"interp", - (PyCFunction)interp, + O2F interp, METH_VARARGS|METH_KEYWORDS, "Multi-dimensional linear interpolation of arrays with gridded data."}, // -------------------------------------------------------------------- {"scalproj_cosangle", - scalproj_cosangle, + O2F scalproj_cosangle, METH_VARARGS, "Calculate normalized projections of vectors."}, // -------------------------------------------------------------------- {"roche_contact_partial_area_volume", - (PyCFunction)roche_contact_partial_area_volume, + O2F roche_contact_partial_area_volume, METH_VARARGS|METH_KEYWORDS, "Determine the area, volume and dvolume/dOmega of the Roche contact lobe" "at given mass ratio q, separation d and and Omega0"}, {"roche_contact_neck_min", - (PyCFunction)roche_contact_neck_min, + O2F roche_contact_neck_min, METH_VARARGS|METH_KEYWORDS, "Determine the minimal distance and position from x axis of the neck " "of the Roche contact lobe at given mass ratio q, separation d and Omega0"}, {"roche_contact_Omega_at_partial_vol", - (PyCFunction)roche_contact_Omega_at_partial_vol, + O2F roche_contact_Omega_at_partial_vol, METH_VARARGS|METH_KEYWORDS, "Determine the value of the potential at a partial volume for the contact " "Roche lobe at given mass ratio q and separation d."}, @@ -11983,18 +11993,18 @@ static PyMethodDef Methods[] = { {"planck_function", - planck_function, + O2F planck_function, METH_VARARGS, "Calculate monochromatic blackbody intensity at a given wavelength and " "temperature."}, {"CCM89_extinction", - CCM89_extinction, + O2F CCM89_extinction, METH_VARARGS, "Calculate CCM89 extinction coefficients for a given wavelength"}, {"gordon_extinction", - gordon_extinction, + O2F gordon_extinction, METH_VARARGS, "Calculate Gordon et al. (2009, UV) and CCM89 (OPT-IR) extinction coefficients for a given wavelength"}, @@ -12002,7 +12012,7 @@ static PyMethodDef Methods[] = { // -------------------------------------------------------------------- {"setup_verbosity", - (PyCFunction)setup_verbosity, + O2F setup_verbosity, METH_VARARGS|METH_KEYWORDS, "Setting the verbosity of libphoebe"}, diff --git a/phoebe/lib/misaligned_roche.h b/phoebe/lib/misaligned_roche.h index 4aa8aad7f..f84008008 100644 --- a/phoebe/lib/misaligned_roche.h +++ b/phoebe/lib/misaligned_roche.h @@ -1,13 +1,13 @@ #pragma once /* - Library dealing with the generalizd Roche potential with misaligned + Library dealing with the generalized Roche potential with misaligned binary system. THERE IS NO SUPPORT FOR (OVER-)CONTACT CASE, AS IT IS NOT PHYSICAL. STILL NOTATION OF TWO LOBES REMAINS. - The surface is defined implicitely with potential + The surface is defined implicitly with potential Omega(x,y,z,params) = 1/r1 + q(1/r2 - x/delta^2) + @@ -636,7 +636,7 @@ template choice : 0 - left 1 - right - 2 - overcontact + 2 - over-contact Omega0 - reference value of the potential q - mass ratio M2/M1 F - synchronicity parameter @@ -671,7 +671,7 @@ template return false; } - // for cases which are not strongy deformed the pole is also + // for cases which are not strongly deformed the pole is also // good point to starts if (!poleL(r, Omega0, q, F, delta, theta)) return false; @@ -899,7 +899,7 @@ template - Using: Integrating surface in spherical coordiantes + Using: Integrating surface in spherical coordinates a. Gauss-Lagrange integration in phi direction b. RK4 in direction of theta @@ -921,7 +921,7 @@ template T v[3], const unsigned & choice, const T & pole, - const T & Omega0, + [[maybe_unused]] const T & Omega0, const T & q, const T & F = 1, const T & d = 1, @@ -942,15 +942,15 @@ template unsigned mask = 3; // = 011b if (b_area) mask |= 4U; // += 100b - + using real = long double; - + using G = glq; - + const int dim = G::n + 3; - - real - W[3], w[G::n], y[dim], k[4][dim], sc_nu[2], sc_th[2], + + real + W[3], w[G::n], y[dim], k[4][dim], sc_nu[2], sc_th[2], q_ = q, d2 = d*d, d3 = d*d2, b = (1 + q)*F*F*d3, dnu = utils::pi()/m, @@ -1170,7 +1170,7 @@ template Wnn = r[0]*(q*(A[0]*(1 - t3) + 3*r[0]*t5*A[1]*A[1]) + b*r[0]*(B00 - B[1]*B[1])), // d^2Omega/dnu^2 Wnr = q*A[1]*(-1 + t5*(1 + r[0]*(A[0] - 2*r[0]))) - 2*b*r[0]*B[0]*B[1]; // d^2Omega/drdnu - + if ((mask & 1U) == 1U) { T det = Wnr*Wnr-Wrr*Wnn; if (det < 0) { @@ -1179,7 +1179,7 @@ template } v[0] = (-Wnr-std::sqrt(det))/Wrr; // dr/dnu } - + if ((mask & 2U) == 2U) v[1] = 1/(Wrr*v[0] + Wnr); // d^2(dV/dOmega)/(dnu dphi) return false; } @@ -1202,7 +1202,7 @@ template spin in plane (x, z) m - number of steps in x - direction - Using: Integrating surface in spherical coordiantes + Using: Integrating surface in spherical coordinates a. Gauss-Lagrange integration in phi direction b. RK4 in direction of theta @@ -1235,16 +1235,16 @@ template const T & th = 0, const int & m = 1 << 14) { - + #if defined(DEBUG) std::cerr.precision(16); const char *fname = "critical_area_volume_integration"; #endif - + // // What is calculated // - + bool b_area = (choice & 1U) == 1U, b_vol = (choice & 2U) == 2U, @@ -1258,14 +1258,14 @@ template if (b_area) mask |= 4; // + 100b if (b_dvol) mask2 |= 2; // + 010b - + using real = long double; using G = glq; - + const int dim = G::n + 3; - - real - w[G::n], y[dim], k[4][dim], sc_nu[2], p[4], W[3], sum[3], + + real + w[G::n], y[dim], k[4][dim], sc_nu[2], p[4], W[3], sum[3], r[2], nu, rt, rp, q_ = q, d2 = d*d, d3 = d2*d, b = (1 + q)*F*F*d3, @@ -1273,45 +1273,45 @@ template #if defined(DEBUG) // - // Checking if the x is a good critical point + // Checking if the x is a good critical point // - + real params[] = {q, F, d, th, 0}; - + Tmisaligned_rotated_roche body(params); - + real B[3][3], g[3], Omega0; - + B[2][0]= x[0]; B[2][1] = 0; B[2][2] = x[1]; - + Omega0 = -body.constrain(B[2]); - + body.grad_only(B[2], g); - + std::cerr.precision(16); - std::cerr - << fname << "::Omega0=" << Omega0 << " grad=" + std::cerr + << fname << "::Omega0=" << Omega0 << " grad=" << g[0] << '\t' << g[1] << '\t' << g[2] << '\n'; - + // // Basis // real t = std::hypot(real(x[0]), real(x[1])); - + B[2][0] /= t; B[2][2] /= t; - + B[0][0] = B[2][2]; B[0][1] = 0; B[0][2] = -B[2][0]; - - B[1][0] = 0; + + B[1][0] = 0; B[1][1] = 1; - B[1][2] = 0; + B[1][2] = 0; #endif - + // // Setting initial values // @@ -1332,9 +1332,9 @@ template w[i] = dnu*G::weights[i]; } y[G::n] = y[G::n + 1] = y[G::n + 2] = real(0); - + #if defined(DEBUG) - for (int i = 0; i < 4; ++i) + for (int i = 0; i < 4; ++i) k[i][G::n] = k[i][G::n + 1] = k[i][G::n + 2] = real(0); #endif } @@ -1349,7 +1349,7 @@ template // 1. step #if defined(DEBUG) W[0]= W[1] = W[2] = 0; - #endif + #endif sum[0] = sum[1] = sum[2] = 0; if (i == 0) { // discussing pole = L1 separately, sin(nu) = 0 @@ -1357,11 +1357,11 @@ template r[0] = y[j], r[1] = r[0]*r[0]; calc_dOmega2_pole(W, mask2, r, G::sc_phi + 2*j, q_, b, p); k[0][j] = dnu*W[0]; - + #if defined(DEBUG) std::cerr << fname << "::k00[" << j << "]=" << k[0][j] << " W=" << W[0] << ':' << W[1] << "\n"; #endif - + if (b_dvol) sum[2] += w[j]*r[1]*W[1]; } @@ -1382,7 +1382,7 @@ template #if defined(DEBUG) std::cerr << fname << "::k0[" << j << "]=" << k[0][j] <<" W=" << W[0] << ':' << W[1] << ':' << W[2] << "\n"; #endif - + if (b_area) { rp = -W[2]/W[0]; // partial_phi r sum[0] += w[j]*r[0]*std::sqrt(rp*rp + sc_nu[0]*sc_nu[0]*(r[1] + rt*rt)); @@ -1412,7 +1412,7 @@ template #if defined(DEBUG) std::cerr << fname << "::k1[" << j << "]=" << k[1][j] << " W=" << W[0] << ':' << W[1] << ':' << W[2] << "\n"; #endif - + if (b_area) { rp = -W[2]/W[0]; // partial_phi r sum[0] += w[j]*r[0]*std::sqrt(rp*rp + sc_nu[0]*sc_nu[0]*(r[1] + rt*rt)); @@ -1427,7 +1427,7 @@ template // 3. step #if defined(DEBUG) W[0]= W[1] = W[2] = 0; - #endif + #endif sum[0] = sum[1] = sum[2] = 0; for (int j = 0; j < G::n; ++j){ r[0] = y[j] + 0.5*k[1][j], r[1] = r[0]*r[0]; @@ -1435,11 +1435,11 @@ template calc_dOmega2(W, mask, r, sc_nu, G::sc_phi + 2*j, q_, b, p); rt = -W[1]/W[0]; // partial_nu r k[2][j] = dnu*rt; - + #if defined(DEBUG) std::cerr << fname << "::k2[" << j << "]=" << k[2][j] << " W=" << W[0] << ':' << W[1] << ':' << W[2] << "\n"; #endif - + if (b_area) { rp = -W[2]/W[0]; // partial_phi r sum[0] += w[j]*r[0]*std::sqrt(rp*rp + sc_nu[0]*sc_nu[0]*(r[1] + rt*rt)); @@ -1454,7 +1454,7 @@ template // 4. step #if defined(DEBUG) W[0]= W[1] = W[2] = 0; - #endif + #endif sum[0] = sum[1] = sum[2] = 0; utils::sincos(nu + dnu, sc_nu, sc_nu+1); for (int j = 0; j < G::n; ++j){ @@ -1483,45 +1483,45 @@ template // final step for (int j = 0; j < dim; ++j) y[j] += (k[0][j] + 2*(k[1][j] + k[2][j]) + k[3][j])/6; - + #if defined(DEBUG) for (int j = 0; j < G::n; ++j) { - + real f, a[3] = {G::sc_phi[2*j+1]*sc_nu[0], G::sc_phi[2*j]*sc_nu[0], sc_nu[1]}, u[3]; - + f = y[j]*d; for (int k = 0; k < 3; ++k) { u[k] = 0; for (int l = 0; l < 3; ++l) u[k]+= f*a[l]*B[l][k]; } - + std::cerr << fname - << "::y[" << j << "]=" << y[j] - << " dOmega=" << + << "::y[" << j << "]=" << y[j] + << " dOmega=" << Omega0 + body.constrain(u) << "\n"; } - - + + for (int j = G::n; j < dim; ++j) std::cerr << << fname << "::y[" << j << "]=" << y[j] << "\n"; #endif - + nu += dnu; } if (b_area) v[0] = 2*d2*y[G::n]; if (b_vol) v[1] = 2*d3*y[G::n + 1]/3; if (b_dvol) v[2] = 2*d3*d*y[G::n + 2]; - + #if defined(DEBUG) if (b_dvol && y[G::n + 2] > 0) { - + std::ofstream f("params.txt",std::ofstream::binary); - + int len = 9; T *buf = new T [len]; - + buf[0]=v[0]; buf[1]=v[1]; buf[2]=v[2]; @@ -1530,14 +1530,14 @@ template buf[5]=q; buf[6]=F; buf[7]=d; - buf[8]=th; - + buf[8]=th; + std::cerr << fname << "::choice=" << choice << " m=" << m << '\n'; - + f.write((char *)buf, sizeof(T)*len); - + f.close(); - + delete [] buf; } #endif @@ -1546,7 +1546,7 @@ template #if defined(DEBUG) #undef DEBUG #endif - + /* Computing volume of the primary generalized Roche lobes with misaligned spin and orbital velocity vector and derivatives of @@ -1777,7 +1777,7 @@ template Return: true - if calculation succeeded, false - otherwise */ - //#define DEBUG + //#define DEBUG template bool lagrange_point( int choice, @@ -1798,11 +1798,11 @@ template case 3: L0 = gen_roche::lagrange_point_L3(q, F, d); break; default: return false; } - + #if defined(DEBUG) std::cerr << "L0=" << L0 << " theta=" << theta << '\n'; #endif - + // // Approximating fixed point using RK4 integration from // position at aligned case @@ -1873,7 +1873,7 @@ template #if defined(DEBUG) #undef DEBUG #endif - + /* Calculate the minimal value of the Kopal potential for which the primary Roche lobe exists. @@ -1884,7 +1884,7 @@ template d - separation between the two objects th - angle between z axis in spin of the object L[2] - Lagrange point limiting the lobe - + Return: minimal permitted Omega, if NaN is returned we have some error */ @@ -1897,12 +1897,12 @@ template const T & th = 0, T *L = 0, T *pth = 0){ - + T W[2], r[2][3], th1 = utils::m_pi*std::abs(std::fmod(th/utils::m_pi + 0.5, 1) - 0.5); - + //std::cerr << "calc_Omega_min::th1=" << th1 << '\n'; - + for (int i = 0; i < 2; ++i) { if (!lagrange_point(i + 1, q, F, d, th1, r[i])) return std::numeric_limits::quiet_NaN(); @@ -1913,16 +1913,16 @@ template W[i] = calc_Omega(r[i], q, F, d, th1); //std::cerr << "calc_Omega_min::W[" << i << "]=" << W[i] << " x=" << r[i][0] << ":" << r[i][2] << '\n'; } - + int ind = (W[0] > W[1] ? 0 : 1); - + if (L) { L[0] = r[ind][0]; L[1] = r[ind][2]; } - + if (pth) *pth = th1; - + return W[ind]; } @@ -1958,43 +1958,43 @@ template const T & th, T & OmegaC, // outputs T av[3]) { - + const char *fname = "critical_area_volume"; - + #if defined(DEBUG) std::cerr.precision(16); - std::cerr + std::cerr << fname << "::START\n" - << fname << "::q=" << q << " F=" << F << " d=" << d << " th=" << th << '\n'; + << fname << "::q=" << q << " F=" << F << " d=" << d << " th=" << th << '\n'; #endif - + if (th == 0) return gen_roche::critical_area_volume(choice, q, F, d, OmegaC, av); - + T x[2], th1; - + OmegaC = calc_Omega_min(q, F, d, th, x, &th1); - + if (std::isnan(OmegaC)) { std::cerr << fname << "::Calculation of Omega_min failed\n"; return false; } - + #if defined(DEBUG) std::cerr << fname << "::x=" << x[0] << ' ' << x[1] << " OmegaC=" << OmegaC << '\n'; #endif - + critical_area_volume_integration(av, choice, x, q, F, d, th1, 1<<10); - + #if defined(DEBUG) - std::cerr + std::cerr << fname << "::av=" << av[0] << ' ' << av[1] << ' ' << av[2] << '\n' << fname << "::critical_volume::END\n"; #endif return true; } - + #if defined(DEBUG) #undef DEBUG #endif diff --git a/phoebe/lib/redistribution.h b/phoebe/lib/redistribution.h index 928ad4e63..f366a1a71 100644 --- a/phoebe/lib/redistribution.h +++ b/phoebe/lib/redistribution.h @@ -243,14 +243,14 @@ class Tredistribution{ /* Calculating weighted redistribution matrix for local redistribution - + Input: - thresh - threshold value + thresh - threshold value weight - weight of the redistribution matrix P - points A - areas TS - redistribution matrix - + Output: TS - redistribution matrix */ @@ -264,52 +264,52 @@ class Tredistribution{ ){ int i, j, N = P.size(); - + std::vector>> D(N); if (thresh == 0) { - + for (i = 0; i < N; ++i) D[i].emplace_back(i, weight); - + } else { - + std::vector S(N, 0); - + T t, t0 = F()(0.0, thresh), *v, a; - + // make an estimate how much space we need j = (N*utils::sqr(thresh))/4; for (i = 0; i < N; ++i) D[i].reserve(j); - + // generate weighted connection matrix and sums of rows for (i = 0; i < N; ++i) { a = A[i]; v = P[i].data; for (j = 0; j < i; ++j) { - + t = F()(utils::__acosf(utils::dot3D(v, P[j].data)), thresh); if (t) { S[i] += t*A[j]; D[i].emplace_back(j, t); - S[j] += t*a; + S[j] += t*a; D[j].emplace_back(i, t); } } - - S[i] += t0*a; + + S[i] += t0*a; D[i].emplace_back(i, t0); } - + // // Calculating redistribution matrix // for (i = 0; i < N; ++i) S[i] = weight*A[i]/S[i]; for (auto && row : D) - for (auto && e : row) e.second *= S[e.first]; + for (auto && e : row) e.second *= S[e.first]; } - + // // Adding matrix D to matrix TS // @@ -317,16 +317,16 @@ class Tredistribution{ TS = D; } else { auto out = TS.begin(); - for (auto && c : D) add_identical(*(out++), c); + for (auto && c : D) add_identical(*(out++), c); } } /* Calculating redistribution matrix for horizontal redistribution - + Input: - thresh - threshold value + thresh - threshold value weight - weight of the redistribution matrix o - direction (unit vector) P - points @@ -335,7 +335,7 @@ class Tredistribution{ Output: TS - redistribution matrix */ - + template void calc_horiz_redistr_matrix( const T & thresh, // input @@ -351,61 +351,61 @@ class Tredistribution{ std::vector>> D(N); if (thresh == 0) { - + for (i = 0; i < N; ++i) D[i].emplace_back(i, weight); - + } else { - - T t, t0 = F()(0.0, thresh), + + T t, t0 = F()(0.0, thresh), c1, s1, c2, s2, a, *p = new T [2*N], *p1, *p2; - + for (i = 0, p1 = p; i < N; ++i){ *(p1++) = t = utils::dot3D(P[i].data, o); *(p1++) = std::sqrt(1 - t*t); } - + std::vector S(N, 0); - + // make an estimate how much space we need j = (2*N*thresh)/utils::pi(); for (i = 0; i < N; ++i) D[i].reserve(j); - + // generate weighted connection matrix and sums of rows for (i = 0, p1 = p; i < N; ++i) { - a = A[i]; + a = A[i]; c1 = *(p1++); s1 = *(p1++); - + for (j = 0, p2 = p; j < i; ++j) { c2 = *(p2++); s2 = *(p2++); - + t = F()(utils::__acosf(s1*s2 + c1*c2), thresh); if (t) { S[i] += t*A[j]; D[i].emplace_back(j, t); - S[j] += t*a; + S[j] += t*a; D[j].emplace_back(i, t); } } - - S[i] += t0*a; + + S[i] += t0*a; D[i].emplace_back(i, t0); } - + delete [] p; - + // // Calculating redistribution matrix // for (i = 0; i < N; ++i) S[i] = weight*A[i]/S[i]; for (auto && row : D) - for (auto && e : row) e.second *= S[e.first]; + for (auto && e : row) e.second *= S[e.first]; } - + // // Adding matrix D to matrix TS // @@ -413,7 +413,7 @@ class Tredistribution{ TS = D; } else { auto out = TS.begin(); - for (auto && c : D) add_identical(*(out++), c); + for (auto && c : D) add_identical(*(out++), c); } } @@ -489,12 +489,12 @@ class Tredistribution{ // add value at duplicated indices std::vector> out; out.reserve(a.size() + b.size()); - + auto it = a.begin(), ite = a.end(); - + int ind = it -> first; T sum = it ->second; - + while (++it != ite) { if (ind == it->first) sum += it-> second; @@ -505,7 +505,7 @@ class Tredistribution{ } } out.emplace_back(ind, sum); - + a = out; } @@ -529,7 +529,7 @@ class Tredistribution{ type - type of surface support V - vectors of vertices Tr - vectors of triangles - N - vectors of normals at triangles/ at vertices + N - vectors of normals at triangles/ at vertices (not used currently) A - vector of areas of triangles Dpars - map of redistribution model parameters which are used to calculate distribution matrices @@ -556,7 +556,7 @@ class Tredistribution{ const Tsupport_type & type, std::vector> & V, std::vector> & Tr, - std::vector> & N, // vertices: NatV, triangle: NatT + [[maybe_unused]] std::vector> & N, // vertices: NatV, triangle: NatT std::vector & A, std::map> & Dpars, std::map & W @@ -578,7 +578,7 @@ class Tredistribution{ std::vector AatV, // areas per vertices (computed if needed) *pAatE = 0; // pointer - + // // number of elements // @@ -657,14 +657,14 @@ class Tredistribution{ std::cerr << fname << "::Projections to sphere failed\n"; return false; } - + calc_horiz_redistr_matrix(h, w.second, o, P, *pAatE, S); - + break; } } } - + trivial_redistr = p.size() == 0 && S.size() == 0; return true; @@ -674,7 +674,7 @@ class Tredistribution{ void mul_add (std::vector & a, std::vector & b){ T sum; - + if (p.size()) { sum = 0; auto ib = b.begin(), eb = b.end(), ip = p.begin(); diff --git a/phoebe/lib/rot_star.h b/phoebe/lib/rot_star.h index 13747e39e..c32f2b9ad 100644 --- a/phoebe/lib/rot_star.h +++ b/phoebe/lib/rot_star.h @@ -214,7 +214,7 @@ namespace rot_star { 2 bit - volume Output: - res[3] = {area, volume} + res[2] = {area, volume} Return: -1 if no results are computed @@ -227,7 +227,7 @@ namespace rot_star { template int area_volume( - double *res, + double res[2], const unsigned &choice, const T & Omega0, const T & omega) { diff --git a/phoebe/lib/triang_marching.h b/phoebe/lib/triang_marching.h index d4b74a922..bb02c4a42 100644 --- a/phoebe/lib/triang_marching.h +++ b/phoebe/lib/triang_marching.h @@ -1,7 +1,7 @@ #pragma once /* - Library for triangulation using maching algorithm specialized for + Library for triangulation using marching algorithm specialized for closed surfaces and connected surfaces. Supports only genus 0 smooth surfaces. @@ -22,13 +22,36 @@ #include "utils.h" #include "triang_mesh.h" +#include "cvec.h" /* - Triangulation of closed surfaces using maching algorithm. + Triangulation of closed surfaces using marching algorithm. */ template struct Tmarching: public Tbody { + + /* + Distance between the two 3D vectors. + + Input: + a,b - vectors + + Return: + |a-b|_2 -- L2 norm of th difference of vectors + */ + + T dist(T *a, T *b){ + // std::hypot(,,) is coming in C++17 + return utils::hypot3(a[0] - b[0], a[1] - b[1], a[2] - b[2]); + } + + T dist2(T *a, T *b){ + T s = 0; + for (int i = 0; i < 3; ++i) s += utils::sqr(a[i] - b[i]); + return s; + } + /* Interval structure for vertex includes the point on the surface and vector base in tangent plane @@ -172,9 +195,9 @@ struct Tmarching: public Tbody { do { n = 0; - + for (int i = 0; i < 3; ++i) r[i] = ri[i]; - + do { // g = (grad F, F) @@ -430,11 +453,11 @@ struct Tmarching: public Tbody { const T min = 10*std::numeric_limits::min(); do { - + nr_iter = 0; - + for (int i = 0; i < 3; ++i) r[i] = ri[i]; - + do { // g = (grad F, F) @@ -475,7 +498,7 @@ struct Tmarching: public Tbody { this->grad_only(r, g, precision); // creating simplified vertex, - // note: std::hypot(,,) is comming in C++17 + // note: std::hypot(,,) is coming in C++17 if (gnorm) fac = 1/(*gnorm = utils::hypot3(g[0], g[1], g[2])); @@ -490,28 +513,6 @@ struct Tmarching: public Tbody { #undef DEBUG #endif - /* - Distance between the two 3D vectors. - - Input: - a,b - vectors - - Return: - |a-b|_2 -- L2 norm of th difference of vectors - */ - - T dist(T *a, T *b){ - // std::hypot(,,) is comming in C++17 - return utils::hypot3(a[0] - b[0], a[1] - b[1], a[2] - b[2]); - } - - T dist2(T *a, T *b){ - T s = 0; - for (int i = 0; i < 3; ++i) s += utils::sqr(a[i] - b[i]); - return s; - } - - int split_angle(Tvertex & v_prev, Tvertex & v, Tvertex & v_next, T *a) { T q[3][2] = {{0.,0.}, {0.,0.}, {0.,0.}}; @@ -849,7 +850,7 @@ struct Tmarching: public Tbody { /* Triangulization using marching method of genus 0 closed and surfaces Has - -- additionals checks + -- additional checks -- support multifronts Input: @@ -865,9 +866,9 @@ struct Tmarching: public Tbody { GatV - norm of the gradient at vertices Return: 0 - no error - 1 - too triangles + 1 - too many triangles 2 - problem with converges - */ + */ int triangulize_full( T init_r[3], T init_g[3], @@ -883,10 +884,10 @@ struct Tmarching: public Tbody { // start with normal precision defined by T precision = false; - - // error + + // error int error = 0; - + V.clear(); Tr.clear(); @@ -897,10 +898,10 @@ struct Tmarching: public Tbody { // Step 0: // typedef std::vector Tfront_polygon; - + // list of frontal polygon, working here as circular list - std::vector lP(1); - + std::vector lP(1); + { Tvertex v, vk; @@ -915,18 +916,18 @@ struct Tmarching: public Tbody { T sa[6], ca[6], qk[3], u[3]; utils::sincos_array(5, utils::m_pi3, sa, ca, delta); - + for (int k = 0; k < 6 && error == 0; ++k){ - - for (int i = 0; i < 3; ++i) + + for (int i = 0; i < 3; ++i) qk[i] = v.r[i] + (u[i] = ca[k]*v.b[0][i] + sa[k]*v.b[1][i]); if (!project_onto_potential(qk, vk, max_iter, v.b[2]) && !slide_over_potential(v.r, v.b[2], u, delta, vk, max_iter)) { std::cerr << "Warning: Projection did not converge for initial frontal polygon.\n"; error = 2; - } - + } + // store points into initial front vk.index = k + 1; // = V.size(); vk.omega_changed = true; @@ -949,7 +950,7 @@ struct Tmarching: public Tbody { // T delta2 = 0.5*delta*delta; // TODO: should be more dynamical - + do { // current front @@ -1178,9 +1179,9 @@ struct Tmarching: public Tbody { T st, ct, qk[3]; Tvertex Pi[6], *vp = Pi; // new front from it_min - + for (int k = 1; k < nt && error == 0; ++k, ++n, ++vp){ - + // rotate in tangent plane ct = c*ca[k] - s*sa[k]; st = c*sa[k] + s*ca[k]; @@ -1214,7 +1215,7 @@ struct Tmarching: public Tbody { << vp->r[0] << ' ' << vp->r[1] << ' ' << vp->r[2] << '\n' << g[0] << ' ' << g[1] << ' ' << g[2] << '\n' << g[3] << '\n'; - + error = 2; } @@ -1245,12 +1246,12 @@ struct Tmarching: public Tbody { } if (Tr.size() >= max_triangles) error = 1; - + } while (error == 0); - - + + } while (lP.size() > 0 && error == 0); - + return error; } @@ -1261,7 +1262,10 @@ struct Tmarching: public Tbody { */ - Tbad_pair check_bad_pairs(Tfront_polygon &P, const T &delta2){ + Tbad_pair check_bad_pairs_old(Tfront_polygon &P, const T &delta2){ + + + if (P.size() <= 3) return Tbad_pair(0, 0); // safeguard int s; @@ -1337,18 +1341,96 @@ struct Tmarching: public Tbody { return Tbad_pair(0, 0); } + Tbad_pair check_bad_pairs_new(Tfront_polygon &P, const T &delta2){ + + + if (P.size() <= 3) return Tbad_pair(0, 0); // safeguard + + int s; + + T a[3]; + + using Tit = typename Tfront_polygon::iterator; + + Tit it_first = P.begin(), it_last = P.end() - 1, buf[6], + *p = buf + 1, *q = buf + 4, p_stop, q_stop; + + p[-1] = it_last; + p[0] = it_first; + p[1] = it_first + 1; + p_stop = it_last - 2; + + while (1) { + + // checking only those it1 - it0 > 1 + q[-1] = p[1]; + q[0] = cnext(p[1], it_first, it_last); + q[1] = cnext(q[0], it_first, it_last); + + q_stop = (p[0] == it_first ? it_last - 1 : it_last); + + while (1) { + + // are on the side the object + if (utils::dot3D(p[0]->b[2], q[0]->b[2]) > 0) { + + // a = q[0] - p[0] + utils::sub3D(q[0]->r, p[0]->r, a); + + // if near enough and looking inside from it and from it1 + if (utils::norm2(a) < delta2) { + + // check if same side of both edges and determine the side + // depending of it_prev -> it -> it_next circle + s = split_angle(*p[-1], *p[0], *p[1], a); + + if (s != 0 && s*split_angle(*q[-1], *q[0], *q[1], a) < 0) { + + // create new last front + #if defined(DEBUG) + std::cerr + << "P.size=" << P.size() + << " i=" << int(it0 - it_begin) + << " j=" << int(it1 - it_begin) + << " len=" << int(it1 + 1 - it) + << std::endl; + #endif + + // std::cerr << "bad_pair1:" + // << p[0] - it_first << ':' << q[0] - it_first << '\n'; + + return Tbad_pair(p[0] - it_first, q[0] - it_first); + } + } + } + + if (q[0] == q_stop) break; + for (int i = 0; i <= 1; ++i) q[i-1] = q[i]; + q[1] = cnext(q[0], it_first, it_last); + } + + if (p[0] == p_stop) break; + for (int i = 0; i <= 1; ++i) p[i-1] = p[i]; + p[1] = cnext(p[0], it_first, it_last); + } + + return Tbad_pair(0, 0); + } + /* Collect all bad points between polygon front P and part of the polygon */ //#define DEBUG - Tbad_pair check_bad_pairs( + Tbad_pair check_bad_pairs_old( Tfront_polygon &P, typename Tfront_polygon::iterator & start, typename Tfront_polygon::iterator & end, const T &delta2 ){ + if (P.size() <= 3) return Tbad_pair(0, 0); // safeguard + int s; T a[3]; @@ -1359,17 +1441,17 @@ struct Tmarching: public Tbody { it_last = P.end() - 1, it0 = start, - it0_next = (it0 == it_last ? it_begin : it0 + 1), - it0_prev = (it0 == it_begin ? it_last : it0 - 1), + it0_next = (it0 == it_last ? it_begin : it0 + 1), // next from it0 + it0_prev = (it0 == it_begin ? it_last : it0 - 1), // prev from it0 it0_prev_start = it0_prev, it0_last = end - 1; while (1) { auto - it1 = (it0_next == it_last ? it_begin : it0_next + 1), - it1_next = (it1 == it_last ? it_begin : it1 + 1), - it1_prev = (it1 == it_begin ? it_last : it1 - 1), + it1 = (it0_next == it_last ? it_begin : it0_next + 1), // next from it0_next + it1_next = (it1 == it_last ? it_begin : it1 + 1), // next from it1 + it1_prev = (it1 == it_begin ? it_last : it1 - 1), // prev from it1 // avoid running twice over the [start,end) it1_last = ( @@ -1444,11 +1526,122 @@ struct Tmarching: public Tbody { #undef DEBUG #endif + /* + Collect all bad points between polygon front P and part + of the polygon + */ + //#define DEBUG + Tbad_pair check_bad_pairs_new( + Tfront_polygon &P, + typename Tfront_polygon::iterator & start, + typename Tfront_polygon::iterator & end, + const T &delta2 ){ + + if (P.size() <= 3) return Tbad_pair(0, 0); // safeguard + + int s; + + T a[3]; + + // checking all combinations of (it0, it1) + auto + it_begin = P.begin(), + it_last = P.end() - 1, + + it0 = start, + it0_next = (it0 == it_last ? it_begin : it0 + 1), // next from it0 + it0_prev = (it0 == it_begin ? it_last : it0 - 1), // prev from it0 + it0_prev_start = it0_prev, + it0_last = end - 1; + + while (1) { + + auto + it1 = (it0_next == it_last ? it_begin : it0_next + 1), // next from it0_next + it1_next = (it1 == it_last ? it_begin : it1 + 1), // next from it1 + it1_prev = (it1 == it_begin ? it_last : it1 - 1), // prev from it1 + + // avoid running twice over the [start,end) + it1_last = ( + it0 == start ? + (it0_prev_start == it_begin ? it_last : it0_prev_start - 1): + it0_prev_start + ); + + while (1) { + + // are on the side the object + if (utils::dot3D(it0->b[2], it1->b[2]) > 0) { + + utils::sub3D(it1->r, it0->r, a); + + // if near enough and looking inside from it and from it1 + if (utils::norm2(a) < delta2) { + + // check if same side of both edges and determine the side + // depending of it_prev -> it -> it_next circle + s = split_angle(*it0_prev, *it0, *it0_next, a); + + if (s != 0 && s*split_angle(*it1_prev, *it1, *it1_next, a) < 0) { + + int ind[2] = { + static_cast(it0 - it_begin), + static_cast(it1 - it_begin) + }; + + // create new last front + #if defined(DEBUG) + std::cerr + << "P.size=" << P.size() + << " i=" << ind[0] + << " j=" << ind[1] + << " len=" << int(it1 + 1 - it0) + << std::endl; + #endif + + //std::cerr << "bad_pair2:" << ind[0] << ':' << ind[0] << '\n'; + + if (ind[0] < ind[1]) return Tbad_pair(ind[0], ind[1]); + return Tbad_pair(ind[1], ind[0]); + } + } + } + + if (it1 == it1_last) break; + + it1_prev = it1; + it1 = it1_next; + + if (it1_next == it_last) + it1_next = it_begin; + else + ++it1_next; + } + + if (it0 == it0_last) break; + + it0_prev = it0; + it0 = it0_next; + + if (it0_next == it_last) + it0_next = it_begin; + else + ++it0_next; + } + + return Tbad_pair(0, 0); + } + #if defined(DEBUG) + #undef DEBUG + #endif + + ALIAS_TEMPLATE_FUNCTION(check_bad_pairs, check_bad_pairs_new) + /* Triangulization using marching method of genus 0 closed and surfaces. Has: - -- additionals checks + -- additional checks -- supports multifronts -- clever detection of bad pairs/points @@ -1464,14 +1657,14 @@ struct Tmarching: public Tbody { NatV - vector of normals at vertices (read N at V) Tr - vector of triangles GatV - norm of the gradient at vertices - + Return: 0 - no error - 1 - too triangles + 1 - too many triangles 2 - problem with converges */ - - int triangulize_full_clever( + + int triangulize_full_clever_old( T init_r[3], T init_g[3], const T & delta, @@ -1480,13 +1673,15 @@ struct Tmarching: public Tbody { std::vector > & NatV, std::vector > & Tr, std::vector * GatV = 0, - const T & init_phi = 0) + const T & init_phi = 0) { + + //std::cerr << "triangulize_full_clever::old\n"; // start with normal precision defined by T precision = false; - - // error + + // error int error = 0; V.clear(); @@ -1524,10 +1719,10 @@ struct Tmarching: public Tbody { T sa[6], ca[6], qk[3], u[3]; utils::sincos_array(5, utils::m_pi3, sa, ca, delta); - + for (int k = 0; k < 6 && error == 0; ++k){ - - for (int i = 0; i < 3; ++i) + + for (int i = 0; i < 3; ++i) qk[i] = v.r[i] + (u[i] = ca[k]*v.b[0][i] + sa[k]*v.b[1][i]); if ( @@ -1558,7 +1753,7 @@ struct Tmarching: public Tbody { // // Triangulization of genus 0 surfaces // - + T delta2 = 0.5*delta*delta; // TODO: should be more dynamical do { @@ -1602,9 +1797,10 @@ struct Tmarching: public Tbody { it0->omega_changed = true; it1->omega_changed = true; - + Tfront_polygon P1(it0, it1 + 1); P.erase(it0 + 1, it1); + B = check_bad_pairs(P, delta2); lP.push_back(P1); @@ -1708,6 +1904,7 @@ struct Tmarching: public Tbody { ) { nt = 1; } + it_prev->omega_changed = true; it_next->omega_changed = true; @@ -1737,9 +1934,9 @@ struct Tmarching: public Tbody { T st, ct, qk[3]; Tvertex Pi[6], *vp = Pi; // new front from it_min - + for (int k = 1; k < nt && error == 0; ++k, ++n, ++vp){ - + // rotate in tangent plane ct = c*ca[k] - s*sa[k]; st = c*sa[k] + s*ca[k]; @@ -1773,7 +1970,7 @@ struct Tmarching: public Tbody { << vp->r[0] << ' ' << vp->r[1] << ' ' << vp->r[2] << '\n' << g[0] << ' ' << g[1] << ' ' << g[2] << '\n' << g[3] << '\n'; - + error = 2; } @@ -1797,11 +1994,10 @@ struct Tmarching: public Tbody { // add vertices to front and replace minimal *(it_min++) = *Pi; - auto - it0 = P.insert(it_min, Pi + 1, Pi + nt - 1), + auto it0 = P.insert(it_min, Pi + 1, Pi + nt - 1), // check if there are any bad pairs - it1 = (--it0) + nt - 1; + it1 = (--it0) + nt - 1; B = check_bad_pairs(P, it0, it1, delta2); @@ -1815,7 +2011,7 @@ struct Tmarching: public Tbody { } if (Tr.size() >= max_triangles) error = 1; - + } while (error == 0); } while (lP.size() > 0 && error == 0); @@ -1823,22 +2019,393 @@ struct Tmarching: public Tbody { return error; } - /* - Calculate the central_points of triangles i.e. barycenters - projected down to surface and normals at that points - - properties: - areas of triangles and a normal - - Input: + int triangulize_full_clever_new( + T init_r[3], + T init_g[3], + const T & delta, + const unsigned & max_triangles, + std::vector > & V, + std::vector > & NatV, + std::vector > & Tr, + std::vector * GatV = 0, + const T & init_phi = 0) + { + + //std::cerr << "triangulize_full_clever::new\n"; + + // start with normal precision defined by T + precision = false; + + // error + int error = 0; + + V.clear(); + Tr.clear(); + + const int max_iter = 100; + + // list of front polygons: front is threated as circular list + std::vector lP(1); + + // list of bad pairs + // pair.first = pair.second means there is no bad pair + std::vector lB; + + // calculate distance between iterators + auto d2 = [&] (auto it0, auto it1) {return dist2(it0->r, it1->r);}; + + // + // Create initial frontal polygon lP[0] and initial bad point lB[0] + // Step 0: + // + + { + Tvertex v, vk; + + Tfront_polygon & P = lP.back(); + + lB.emplace_back(0,0); // no bad pair detected + + // construct the vector base + create_internal_vertex(init_r, init_g, v, init_phi); + + // add vertex to the set, index 0 + V.emplace_back(v.r); // saving only r + if (GatV) GatV->emplace_back(v.norm); // saving g + NatV.emplace_back(v.b[2]); // saving only normal + + T sa[6], ca[6], qk[3], u[3]; + + utils::sincos_array(5, utils::m_pi3, sa, ca, delta); + + for (int k = 0; k < 6 && error == 0; ++k){ + + for (int i = 0; i < 3; ++i) + qk[i] = v.r[i] + (u[i] = ca[k]*v.b[0][i] + sa[k]*v.b[1][i]); + + if ( + !slide_over_potential(v.r, v.b[2], u, delta, vk, max_iter) && + !project_onto_potential(qk, vk, max_iter, v.b[2]) + ) { + std::cerr << "Warning: Projection did not converge for initial frontal polygon!\n"; + error = 2; + } + + // store points into initial front + vk.index = k + 1; // = V.size(); + vk.omega_changed = true; + P.push_back(vk); + + V.emplace_back(vk.r); // saving only r + if (GatV) GatV->emplace_back(vk.norm); // saving norm + NatV.emplace_back(vk.b[2]); // saving only normal + } + + // + // Creating initial hexagon -- triangle faces in Tr + // + for (int k = 0; k < 5; ++k) Tr.emplace_back(0, k + 1, k + 2); + Tr.emplace_back(0, 6, 1); + } + + // + // Triangulization of genus 0 surfaces + // + + T delta2 = 0.5*delta*delta; // TODO: should be more dynamical + + do { + + // current front polygon + Tfront_polygon & P = lP.back(); + Tbad_pair & B = lB.back(); + + do { + + // + // Processing the last three vertices + // + if (P.size() == 3) { + Tr.emplace_back(P[0].index, P[1].index, P[2].index); + + // erasing discussed front + lP.pop_back(); + + // erase discussed possible bad pair + lB.pop_back(); + + break; + } + + // pointers associated to the front + auto it_begin = P.begin(), it_end = P.end(), it_last = it_end - 1; + + // + // If a non-neighboring vertices are too close form new fronts + // Step 2 + // + { + // if bad pair is set do the cut of the front + if (B.first != B.second) { + + // separate fronts P -> P1, P2 + auto + it0 = it_begin + B.first, + it1 = it_begin + B.second; + + it0->omega_changed = true; + it1->omega_changed = true; + + // split front into two parts with common two points it0 and it1 + auto P1 = ccopy(it0, it1, it_begin, it_last); + auto P2 = ccopy(it1, it0, it_begin, it_last); + + // updating referenced P and B + P = P1; + B = check_bad_pairs(P, delta2); + + // add polygonal front and bad_pair to the list + lP.push_back(P2); + lB.push_back(check_bad_pairs(P2, delta2)); + + //std::cerr << "sizes:" << P1.size() << '\t' << P2.size() << '\n'; + + break; + } + } + + // + // Calculate the front angles and choose the point with the smallest + // Step 1 + // + + T omega_min = utils::m_2pi; + + typename Tfront_polygon::iterator it_min; + + { + + T omega, t, tt, c, s, st, ct; + + // set it_prev, it, it_next: as circular list + auto + it = it_begin, + it_next = it + 1, // = cnext(it, it_begin, it_last) + it_prev = it_last; // = cprev(it, it_begin, it_last) + + while (1) { + + if (it -> omega_changed) { // calculate frontal angle if need + + c = s = ct = st = 0; + for (int i = 0; i < 3; ++i) { + t = it_prev->r[i] - it->r[i]; // = dr1[i], dr1 = p_prev - p_cur + c += t*it->b[0][i]; // = dr1[i]*t1[i] + s += t*it->b[1][i]; // = dr1[i]*t2[i] + + tt = it_next->r[i] - it->r[i]; // = dr2[i], dr2 = p_next - p_cur + ct += tt*it->b[0][i]; // = dr2[i]*t1[i] + st += tt*it->b[1][i]; // = dr2[i]*t2[i] + } + + // = arg[ dr1.dr2 + I k.(dr1 x dr2) ] + // omega = atan2(st,ct) - atan2(s,c); + omega = std::atan2(c*st - s*ct, c*ct + s*st); + + // omega = omega mod 2 Pi (offset 0) + if (omega < 0) omega += utils::m_2pi; + + it -> omega = omega; + it -> omega_changed = false; + + } else omega = it -> omega; + + // saving the minimal value of omega + if (omega < omega_min) { + it_min = it; + omega_min = omega; + } + + // cyclic permutation of pointers + it_prev = it; + it = it_next; + + if (it_next == it_begin) break; + + //it_next = cnext(it_next, it_begin, it_last); + if (it_next == it_last) + it_next = it_begin; + else + ++it_next; + } + } + + // + // Discuss the point with the minimal angle + // Step 3 + // + + { + // prepare pointers to vertices in P + // auto it_prev = cprev(it_min, it_begin, it_last); + // auto it_next = cnext(it_min, it_begin, it_last); + + auto + it_prev = it_min, + it_next = it_min; + + if (it_min != it_begin) --it_prev; else it_prev = it_last; + if (it_min != it_last) ++it_next; else it_next = it_begin; + + // number of triangles to be generated + int nt = int(omega_min/utils::m_pi3) + 1; + + T domega = omega_min/nt; + + // correct domega for extreme cases + if (domega < 0.8 && nt > 1) { + domega = omega_min/(--nt); + } else if (nt == 1 && domega > 0.8 && + d2(it_next, it_prev) > 1.4*delta2) { + domega = omega_min/(++nt); + } else if (omega_min < 3 && + (d2(it_min, it_prev) < 0.25*delta2 || + d2(it_min, it_next) < 0.25*delta2) ) { + nt = 1; + } + + it_prev->omega_changed = true; + it_next->omega_changed = true; + + if (nt > 1) { + + // projection of dr = p_next - p_min to tangent space + // c = dr.t1 + // s = dr.t2 + + T c = 0, s = 0, t; + + for (int i = 0; i < 3; ++i){ + t = it_prev->r[i] - it_min->r[i]; // = dr[i] + c += t*it_min->b[0][i]; // = dr[i]*t1[i] + s += t*it_min->b[1][i]; // = dr[i]*t2[i] + } + + // returning fac*(sin(k domega), cos(k domega)) + // where fac = delta/|(c, s)| + + T sa[6], ca[6], u[3]; + + utils::sincos_array(nt - 1, domega, sa, ca, delta/std::hypot(c, s)); + + int n = V.size(); // size of the set of vertices + + T st, ct, qk[3]; + + Tvertex Pi[6], *vp = Pi; // new front from it_min + + for (int k = 1; k < nt && error == 0; ++k, ++n, ++vp){ + + // rotate in tangent plane + ct = c*ca[k] - s*sa[k]; + st = c*sa[k] + s*ca[k]; + + // forming point on tangent plane + for (int i = 0; i < 3; ++i) + qk[i] = it_min->r[i] + (u[i] = it_min->b[0][i]*ct + it_min->b[1][i]*st); + + if (!project_onto_potential(qk, *vp, max_iter, it_min->b[2]) && + !slide_over_potential(it_min->r, it_min->b[2], u, delta, *vp, max_iter)) { + + T g[4]; + + std::cerr << "Warning: Projection did not converge\n"; + + this->grad(qk, g); + + std::cerr.precision(16); + + std::cerr + << "Start\n" + << qk[0] << ' ' << qk[1] << ' ' << qk[2] << '\n' + << g[0] << ' ' << g[1] << ' ' << g[2] << '\n' + << g[3] << '\n'; + + + this->grad(vp->r, g); + + std::cerr + << "End\n" + << vp->r[0] << ' ' << vp->r[1] << ' ' << vp->r[2] << '\n' + << g[0] << ' ' << g[1] << ' ' << g[2] << '\n' + << g[3] << '\n'; + + error = 2; + } + + vp->index = n; // = V.size(); + vp->omega_changed = true; + + // V.emplace_back(vp->r, vp->b[2]); + V.emplace_back(vp->r); // saving only r + if (GatV) GatV->emplace_back(vp->norm); // saving g + NatV.emplace_back(vp->b[2]); // saving only normal + + // add triangle + Tr.emplace_back((k == 1 ? it_prev->index : n - 1), n, it_min->index); + } + + // Note: n = V.size(); + + // add triangle + Tr.emplace_back(n - 1, it_next->index, it_min->index); + + // add vertices to front and replace minimal + *(it_min++) = *Pi; + + auto it0 = P.insert(it_min, Pi + 1, Pi + nt - 1), + + // check if there are any bad pairs + it1 = (--it0) + nt - 1; + + B = check_bad_pairs(P, it0, it1, delta2); + + } else { + // add triangle + Tr.emplace_back(it_prev->index, it_next->index, it_min->index); + + // erase vertex from the front + P.erase(it_min); + } + } + + if (Tr.size() >= max_triangles) error = 1; + + } while (error == 0); + + } while (lP.size() > 0 && error == 0); + + return error; + } + + ALIAS_TEMPLATE_FUNCTION(triangulize_full_clever, triangulize_full_clever_new) + + /* + Calculate the central_points of triangles i.e. barycenters + projected down to surface and normals at that points + + properties: + areas of triangles and a normal + + Input: V - vector of vertices Tr - vector of triangles - connectivity matrix Output: C - central points NatC - normals at central points - GatC - norm of gradient at cetral points - + GatC - norm of gradient at central points + Return: true if ok, false if not ok */ @@ -1872,7 +2439,7 @@ struct Tmarching: public Tbody { const int max_iter = 100; const T eps = 100*std::numeric_limits::epsilon(); - + T *tp, v[3], n[3], q[3], r[3][3]; int i, j; @@ -1899,7 +2466,7 @@ struct Tmarching: public Tbody { if (GatC) GatC->emplace_back(g); } else return false; //std::cerr << "central_points::Warning: Projection did not converge\n"; } - + return true; } diff --git a/phoebe/lib/utils.h b/phoebe/lib/utils.h index b352f33df..51821b493 100644 --- a/phoebe/lib/utils.h +++ b/phoebe/lib/utils.h @@ -14,6 +14,15 @@ #include "sincos.h" +// macro to make aliases of functions +// usage: ALIAS_TEMPLATE_FUNCTION(fun_alias, fun) +#define ALIAS_TEMPLATE_FUNCTION(highLevelF, lowLevelF) \ +template \ +inline auto highLevelF(Args&&... args) -> decltype(lowLevelF(std::forward(args)...)) \ +{ \ + return lowLevelF(std::forward(args)...); \ +} + namespace utils { const double m_pi = 3.14159265358979323846264338327950419984; // pi @@ -27,31 +36,31 @@ namespace utils { template constexpr T pi() {return T(3.14159265358979323846264338327950419984L);}; - + /* Approximation formula for ArcCos - + Input: x in [-1,1] Return: arccos(x) in [0,pi] - - Ref: - p.81 Handbook of Mathematical Functions, by M. Abramowitz and I. Stegun - */ + + Ref: + p.81 Handbook of Mathematical Functions, by M. Abramowitz and I. Stegun + */ float __acosf(const float & x) { - + if (x == 0) return 1.57079632679489; if (x >= 1) return 0; if (x <= -1) return 3.14159265358979; - - float + + float t = std::abs(x), s = std::sqrt(1-t)*(1.5707288 + t*(-0.2121144 + (0.074261 - 0.0187293*t)*t)); - + return (x > 0 ? s : 3.14159265358979 - s); } - + /* Return square of the value. @@ -311,7 +320,6 @@ namespace utils { return x[0]*x[0] + x[1]*x[1] + x[2]*x[2]; } - /* Calculate L2 norm of 3D vector diff --git a/setup.py b/setup.py index 0e7a9b5ad..0ea1fdc9e 100644 --- a/setup.py +++ b/setup.py @@ -6,7 +6,7 @@ 'libphoebe', sources=['phoebe/lib/libphoebe.cpp'], language='c++', - extra_compile_args=["-std=c++11"], + extra_compile_args=["-std=c++14"], include_dirs=[numpy.get_include()] ), Extension( diff --git a/tests/tests/test_dynamics/test_dynamics.py b/tests/tests/test_dynamics/test_dynamics.py deleted file mode 100644 index 01621c2cd..000000000 --- a/tests/tests/test_dynamics/test_dynamics.py +++ /dev/null @@ -1,169 +0,0 @@ -""" -""" - -import phoebe -from phoebe import u -import numpy as np -try: - import matplotlib.pyplot as plt - PLOTTING_ENABLED = True -except ImportError: - PLOTTING_ENABLED = False - - -def _keplerian_v_nbody(b, plot=False): - """ - test a single bundle for the phoebe backend's kepler vs nbody dynamics methods - """ - - # TODO: loop over ltte=True,False (once keplerian dynamics supports the switch) - - # b.add_compute(dynamics_method='bs') - b.set_value('dynamics_method', 'bs') - - times = np.linspace(0, 100, 10000) - nb_ts, nb_us, nb_vs, nb_ws, nb_vus, nb_vvs, nb_vws = phoebe.dynamics.nbody.dynamics_from_bundle(b, times, ltte=False) - k_ts, k_us, k_vs, k_ws, k_vus, k_vvs, k_vws = phoebe.dynamics.keplerian.dynamics_from_bundle(b, times) - - assert np.allclose(nb_ts, k_ts, 1e-8) - for ci in range(len(b.hierarchy.get_stars())): - # TODO: make atol lower (currently 1e-5 solRad which is awfully big, but 1e-6 currently fails!) - if plot: - print("max atol xs:", nb_us[ci] - k_us[ci]) - print("max atol ys:", nb_vs[ci] - k_vs[ci]) - print("max atol zs:", nb_ws[ci] - k_ws[ci]) - print("max atol vxs:", nb_vus[ci] - k_vus[ci]) - print("max atol vys:", nb_vvs[ci] - k_vvs[ci]) - print("max atol vzs:", nb_vws[ci] - k_vws[ci]) - - assert np.allclose(nb_us[ci], k_us[ci], atol=1e-5) - assert np.allclose(nb_vs[ci], k_vs[ci], atol=1e-5) - assert np.allclose(nb_ws[ci], k_ws[ci], atol=1e-5) - assert np.allclose(nb_vus[ci], k_vus[ci], atol=1e-4) - assert np.allclose(nb_vvs[ci], k_vvs[ci], atol=1e-4) - assert np.allclose(nb_vws[ci], k_vws[ci], atol=1e-4) - - -def _phoebe_v_photodynam(b, plot=False): - """ - test a single bundle for phoebe's nbody vs photodynam via the frontend - """ - - times = np.linspace(0, 100, 1000) - - b.add_dataset('orb', times=times, dataset='orb01', component=b.hierarchy.get_stars()) - # photodynam and phoebe should have the same nbody defaults... if for some reason that changes, - # then this will probably fail - b.add_compute('photodynam', compute='pdcompute') - # photodynam backend ONLY works with ltte=True, so we will run the phoebe backend with that as well - # TODO: remove distortion_method='nbody' once that is supported - # NOTE: bs is the exact same as that used in photodynam. Nbody and rebound are slightly different. - b.set_value('dynamics_method', 'bs') - b.set_value('ltte', True) - - b.run_compute('pdcompute', model='pdresults') - b.run_compute('phoebe01', model='phoeberesults') - - for comp in b.hierarchy.get_stars(): - # TODO: check to see how low we can make atol (or change to rtol?) - # TODO: look into justification of flipping x and y for both dynamics (photodynam & phoebe) - # TODO: why the small discrepancy (visible especially in y, still <1e-11) - possibly a difference in time0 or just a precision limit in the photodynam backend since loading from a file?? - - if PLOTTING_ENABLED and plot: - for k in ['us', 'vs', 'ws', 'vus', 'vvs', 'vws']: - plt.cla() - plt.plot(b.get_value('times', model='phoeberesults', component=comp, unit=u.d), b.get_value(k, model='phoeberesults', component=comp), 'r-') - plt.plot(b.get_value('times', model='phoeberesults', component=comp, unit=u.d), b.get_value(k, model='pdresults', component=comp), 'b-') - diff = abs(b.get_value(k, model='phoeberesults', component=comp) - b.get_value(k, model='pdresults', component=comp)) - print("*** max abs ({}): {}".format(k, max(diff))) - plt.xlabel('t') - plt.ylabel(k) - plt.show() - - assert np.allclose(b.get_value('times', dataset='orb01', model='phoeberesults', component=comp, unit=u.d), b.get_value('times', dataset='orb01', model='pdresults', component=comp, unit=u.d), atol=1e-6) - assert np.allclose(b.get_value('us', dataset='orb01', model='phoeberesults', component=comp, unit=u.AU), b.get_value('us', dataset='orb01', model='pdresults', component=comp, unit=u.AU), atol=1e-6) - assert np.allclose(b.get_value('vs', dataset='orb01', model='phoeberesults', component=comp, unit=u.AU), b.get_value('vs', dataset='orb01', model='pdresults', component=comp, unit=u.AU), atol=1e-6) - assert np.allclose(b.get_value('ws', dataset='orb01', model='phoeberesults', component=comp, unit=u.AU), b.get_value('ws', dataset='orb01', model='pdresults', component=comp, unit=u.AU), atol=1e-6) - assert np.allclose(b.get_value('vus', dataset='orb01', model='phoeberesults', component=comp, unit=u.solRad/u.d), b.get_value('vus', dataset='orb01', model='pdresults', component=comp, unit=u.solRad/u.d), atol=1e-6) - assert np.allclose(b.get_value('vvs', dataset='orb01', model='phoeberesults', component=comp, unit=u.solRad/u.d), b.get_value('vvs', dataset='orb01', model='pdresults', component=comp, unit=u.solRad/u.d), atol=1e-6) - assert np.allclose(b.get_value('vws', dataset='orb01', model='phoeberesults', component=comp, unit=u.solRad/u.d), b.get_value('vws', dataset='orb01', model='pdresults', component=comp, unit=u.solRad/u.d), atol=1e-6) - - -def _frontend_v_backend(b, plot=False): - """ - test a single bundle for the frontend vs backend access to both kepler and nbody dynamics - """ - - # TODO: loop over ltte=True,False - - times = np.linspace(0, 100, 21) - b.add_dataset('orb', times=times, dataset='orb01', component=b.hierarchy.get_stars()) - b.rename_compute('phoebe01', 'nbody') - b.set_value('dynamics_method', 'bs') - - b.add_compute('phoebe', dynamics_method='keplerian', compute='keplerian') - - # NBODY - # do backend Nbody - b_ts, b_us, b_vs, b_ws, b_vus, b_vvs, b_vws = phoebe.dynamics.nbody.dynamics_from_bundle(b, times, compute='nbody') - - # do frontend Nbody - b.run_compute('nbody', model='nbodyresults') - - for ci, comp in enumerate(b.hierarchy.get_stars()): - # TODO: can we lower tolerance? - assert np.allclose(b.get_value('times', dataset='orb01', model='nbodyresults', component=comp, unit=u.d), b_ts, atol=1e-6) - assert np.allclose(b.get_value('us', dataset='orb01', model='nbodyresults', component=comp, unit=u.solRad), b_us[ci], atol=1e-5) - assert np.allclose(b.get_value('vs', dataset='orb01', model='nbodyresults', component=comp, unit=u.solRad), b_vs[ci], atol=1e-5) - assert np.allclose(b.get_value('ws', dataset='orb01', model='nbodyresults', component=comp, unit=u.solRad), b_ws[ci], atol=1e-5) - assert np.allclose(b.get_value('vus', dataset='orb01', model='nbodyresults', component=comp, unit=u.solRad/u.d), b_vus[ci], atol=1e-4) - assert np.allclose(b.get_value('vvs', dataset='orb01', model='nbodyresults', component=comp, unit=u.solRad/u.d), b_vvs[ci], atol=1e-4) - assert np.allclose(b.get_value('vws', dataset='orb01', model='nbodyresults', component=comp, unit=u.solRad/u.d), b_vws[ci], atol=1e-4) - - # KEPLERIAN - # do backend keplerian - b_ts, b_xs, b_ys, b_zs, b_vxs, b_vys, b_vzs = phoebe.dynamics.keplerian.dynamics_from_bundle(b, times, compute='keplerian') - - # do frontend keplerian - b.run_compute('keplerian', model='keplerianresults') - - # TODO: loop over components and assert - for ci, comp in enumerate(b.hierarchy.get_stars()): - # TODO: can we lower tolerance? - assert np.allclose(b.get_value('times', dataset='orb01', model='keplerianresults', component=comp, unit=u.d), b_ts, atol=1e-6) - assert np.allclose(b.get_value('us', dataset='orb01', model='keplerianresults', component=comp, unit=u.solRad), b_us[ci], atol=1e-5) - assert np.allclose(b.get_value('vs', dataset='orb01', model='keplerianresults', component=comp, unit=u.solRad), b_vs[ci], atol=1e-5) - assert np.allclose(b.get_value('ws', dataset='orb01', model='keplerianresults', component=comp, unit=u.solRad), b_ws[ci], atol=1e-5) - assert np.allclose(b.get_value('vus', dataset='orb01', model='keplerianresults', component=comp, unit=u.solRad/u.d), b_vus[ci], atol=1e-4) - assert np.allclose(b.get_value('vvs', dataset='orb01', model='keplerianresults', component=comp, unit=u.solRad/u.d), b_vvs[ci], atol=1e-4) - assert np.allclose(b.get_value('vws', dataset='orb01', model='keplerianresults', component=comp, unit=u.solRad/u.d), b_vws[ci], atol=1e-4) - - -def test_binary(plot=False): - """ - """ - phoebe.devel_on() # required for nbody dynamics - - # TODO: grid over orbital parameters - # TODO: once ps.copy is implemented, just send b.copy() to each of these - - b = phoebe.default_binary() - b.get_parameter('dynamics_method')._choices = ['keplerian', 'bs'] - _keplerian_v_nbody(b, plot=plot) - - b = phoebe.default_binary() - b.get_parameter('dynamics_method')._choices = ['keplerian', 'bs'] - _phoebe_v_photodynam(b, plot=plot) - - b = phoebe.default_binary() - b.get_parameter('dynamics_method')._choices = ['keplerian', 'bs'] - _frontend_v_backend(b, plot=plot) - - phoebe.devel_off() # reset for future tests - - -if __name__ == '__main__': - logger = phoebe.logger(clevel='INFO') - test_binary(plot=True) - - # TODO: create tests for both triple configurations (A--B-C, A-B--C) - these should first be default bundles diff --git a/tests/tests/test_dynamics/test_dynamics_grid.py b/tests/tests/test_dynamics/test_dynamics_grid.py deleted file mode 100644 index 0945e029b..000000000 --- a/tests/tests/test_dynamics/test_dynamics_grid.py +++ /dev/null @@ -1,172 +0,0 @@ -""" -""" - -import phoebe -from phoebe import u -import numpy as np -try: - import matplotlib.pyplot as plt - PLOTTING_ENABLED = True -except ImportError: - PLOTTING_ENABLED = False - - -def _keplerian_v_nbody(b, ltte, period, plot=False): - """ - test a single bundle for the phoebe backend's kepler vs nbody dynamics methods - """ - - # TODO: loop over ltte=True,False (once keplerian dynamics supports the switch) - - b.set_value('dynamics_method', 'bs') - - times = np.linspace(0, 5*period, 101) - nb_ts, nb_us, nb_vs, nb_ws, nb_vus, nb_vvs, nb_vws = phoebe.dynamics.nbody.dynamics_from_bundle(b, times, ltte=ltte) - k_ts, k_us, k_vs, k_ws, k_vus, k_vvs, k_vws = phoebe.dynamics.keplerian.dynamics_from_bundle(b, times, ltte=ltte) - - assert np.allclose(nb_ts, k_ts, 1e-8) - for ci in range(len(b.hierarchy.get_stars())): - # TODO: make rtol lower if possible - assert np.allclose(nb_us[ci], k_us[ci], rtol=1e-5, atol=1e-2) - assert np.allclose(nb_vs[ci], k_vs[ci], rtol=1e-5, atol=1e-2) - assert np.allclose(nb_ws[ci], k_ws[ci], rtol=1e-5, atol=1e-2) - - # nbody ltte velocities are wrong so only check velocities if ltte off - if not ltte: - assert np.allclose(nb_vus[ci], k_vus[ci], rtol=1e-5, atol=1e-2) - assert np.allclose(nb_vvs[ci], k_vvs[ci], rtol=1e-5, atol=1e-2) - assert np.allclose(nb_vws[ci], k_vws[ci], rtol=1e-5, atol=1e-2) - - -def _phoebe_v_photodynam(b, period, plot=False): - """ - test a single bundle for phoebe's nbody vs photodynam via the frontend - """ - - times = np.linspace(0, 5*period, 21) - - b.add_dataset('orb', times=times, dataset='orb01', component=b.hierarchy.get_stars()) - # photodynam and phoebe should have the same nbody defaults... if for some reason that changes, - # then this will probably fail - b.add_compute('photodynam', compute='pdcompute') - # photodynam backend ONLY works with ltte=True, so we will run the phoebe backend with that as well - # TODO: remove distortion_method='nbody' once that is supported - b.set_value('dynamics_method', 'bs') - b.set_value('ltte', True) - - b.run_compute('pdcompute', model='pdresults') - b.run_compute('phoebe01', model='phoeberesults') - - for comp in b.hierarchy.get_stars(): - # TODO: check to see how low we can make atol (or change to rtol?) - # TODO: look into justification of flipping x and y for both dynamics (photodynam & phoebe) - # TODO: why the small discrepancy (visible especially in y, still <1e-11) - possibly a difference in time0 or just a precision limit in the photodynam backend since loading from a file?? - - if PLOTTING_ENABLED and plot: - for k in ['us', 'vs', 'ws', 'vus', 'vvs', 'vws']: - plt.cla() - plt.plot(b.get_value('times', model='phoeberesults', component=comp, unit=u.d), b.get_value(k, model='phoeberesults', component=comp), 'r-') - plt.plot(b.get_value('times', model='phoeberesults', component=comp, unit=u.d), b.get_value(k, model='pdresults', component=comp), 'b-') - diff = abs(b.get_value(k, model='phoeberesults', component=comp) - b.get_value(k, model='pdresults', component=comp)) - print("*** max abs: {}".format(max(diff))) - plt.xlabel('t') - plt.ylabel(k) - plt.show() - - assert np.allclose(b.get_value('times', model='phoeberesults', component=comp, unit=u.d), b.get_value('times', model='pdresults', component=comp, unit=u.d), rtol=0, atol=1e-05) - assert np.allclose(b.get_value('us', model='phoeberesults', component=comp, unit=u.AU), b.get_value('us', model='pdresults', component=comp, unit=u.AU), rtol=0, atol=1e-05) - assert np.allclose(b.get_value('vs', model='phoeberesults', component=comp, unit=u.AU), b.get_value('vs', model='pdresults', component=comp, unit=u.AU), rtol=0, atol=1e-05) - assert np.allclose(b.get_value('ws', model='phoeberesults', component=comp, unit=u.AU), b.get_value('ws', model='pdresults', component=comp, unit=u.AU), rtol=0, atol=1e-05) - # assert np.allclose(b.get_value('vxs', model='phoeberesults', component=comp, unit=u.solRad/u.d), b.get_value('vxs', model='pdresults', component=comp, unit=u.solRad/u.d), rtol=0, atol=1e-05) - # assert np.allclose(b.get_value('vys', model='phoeberesults', component=comp, unit=u.solRad/u.d), b.get_value('vys', model='pdresults', component=comp, unit=u.solRad/u.d), rtol=0, atol=1e-05) - # assert np.allclose(b.get_value('vzs', model='phoeberesults', component=comp, unit=u.solRad/u.d), b.get_value('vzs', model='pdresults', component=comp, unit=u.solRad/u.d), rtol=0, atol=1e-05) - - -def _frontend_v_backend(b, ltte, period, plot=False): - """ - test a single bundle for the frontend vs backend access to both kepler and nbody dynamics - """ - - # TODO: loop over ltte=True,False - - times = np.linspace(0, 5*period, 101) - b.add_dataset('orb', times=times, dataset='orb01', component=b.hierarchy.get_stars()) - b.rename_compute('phoebe01', 'nbody') - b.set_value('dynamics_method', 'bs') - b.set_value('ltte', ltte) - - b.add_compute('phoebe', dynamics_method='keplerian', compute='keplerian', ltte=ltte) - - # NBODY - # do backend Nbody - b_ts, b_us, b_vs, b_ws, b_vus, b_vvs, b_vws = phoebe.dynamics.nbody.dynamics_from_bundle(b, times, compute='nbody', ltte=ltte) - - # do frontend Nbody - b.run_compute('nbody', model='nbodyresults') - - for ci, comp in enumerate(b.hierarchy.get_stars()): - # TODO: can we lower tolerance? - assert np.allclose(b.get_value('times', model='nbodyresults', component=comp, unit=u.d), b_ts, rtol=0, atol=1e-6) - assert np.allclose(b.get_value('us', model='nbodyresults', component=comp, unit=u.solRad), b_us[ci], rtol=1e-7, atol=1e-4) - assert np.allclose(b.get_value('vs', model='nbodyresults', component=comp, unit=u.solRad), b_vs[ci], rtol=1e-7, atol=1e-4) - assert np.allclose(b.get_value('ws', model='nbodyresults', component=comp, unit=u.solRad), b_ws[ci], rtol=1e-7, atol=1e-4) - if not ltte: - assert np.allclose(b.get_value('vus', model='nbodyresults', component=comp, unit=u.solRad/u.d), b_vus[ci], rtol=1e-7, atol=1e-4) - assert np.allclose(b.get_value('vvs', model='nbodyresults', component=comp, unit=u.solRad/u.d), b_vvs[ci], rtol=1e-7, atol=1e-4) - assert np.allclose(b.get_value('vws', model='nbodyresults', component=comp, unit=u.solRad/u.d), b_vws[ci], rtol=1e-7, atol=1e-4) - - # KEPLERIAN - # do backend keplerian - b_ts, b_us, b_vs, b_ws, b_vus, b_vvs, b_vws = phoebe.dynamics.keplerian.dynamics_from_bundle(b, times, compute='keplerian', ltte=ltte) - - # do frontend keplerian - b.run_compute('keplerian', model='keplerianresults') - - for ci, comp in enumerate(b.hierarchy.get_stars()): - # TODO: can we lower tolerance? - assert np.allclose(b.get_value('times', model='keplerianresults', component=comp, unit=u.d), b_ts, rtol=0, atol=1e-08) - assert np.allclose(b.get_value('us', model='keplerianresults', component=comp, unit=u.solRad), b_us[ci], rtol=0, atol=1e-08) - assert np.allclose(b.get_value('vs', model='keplerianresults', component=comp, unit=u.solRad), b_vs[ci], rtol=0, atol=1e-08) - assert np.allclose(b.get_value('ws', model='keplerianresults', component=comp, unit=u.solRad), b_ws[ci], rtol=0, atol=1e-08) - assert np.allclose(b.get_value('vus', model='keplerianresults', component=comp, unit=u.solRad/u.d), b_vus[ci], rtol=0, atol=1e-08) - assert np.allclose(b.get_value('vvs', model='keplerianresults', component=comp, unit=u.solRad/u.d), b_vvs[ci], rtol=0, atol=1e-08) - assert np.allclose(b.get_value('vws', model='keplerianresults', component=comp, unit=u.solRad/u.d), b_vws[ci], rtol=0, atol=1e-08) - - -def test_binary(plot=False): - """ - """ - phoebe.devel_on() # required for nbody dynamics - - # TODO: once ps.copy is implemented, just send b.copy() to each of these - - # system = [sma (AU), period (d)] - system1 = [0.05, 2.575] - system2 = [1., 257.5] - system3 = [40., 65000.] - - for system in [system1, system2, system3]: - for q in [0.5, 1.]: - for ltte in [True, False]: - print("test_dynamics_grid: sma={}, period={}, q={}, ltte={}".format(system[0], system[1], q, ltte)) - b = phoebe.default_binary() - b.get_parameter('dynamics_method')._choices = ['keplerian', 'bs'] - - b.set_default_unit_all('sma', u.AU) - b.set_default_unit_all('period', u.d) - - b.set_value('sma@binary', system[0]) - b.set_value('period@binary', system[1]) - b.set_value('q', q) - - _keplerian_v_nbody(b, ltte, system[1], plot=plot) - _frontend_v_backend(b, ltte, system[1], plot=plot) - - phoebe.devel_off() # reset for future tests - - -if __name__ == '__main__': - logger = phoebe.logger(clevel='INFO') - test_binary(plot=True) - - # TODO: create tests for both triple configurations (A--B-C, A-B--C) - these should first be default bundles