[Debian-astro-commits] [gyoto] 214/221: NML: adding the rico() Property read by Polish Doughnut
Thibaut Jean-Claude Paumard
thibaut at moszumanska.debian.org
Fri May 22 20:52:48 UTC 2015
This is an automated email from the git hooks/post-receive script.
thibaut pushed a commit to branch master
in repository gyoto.
commit e8ddcf706e7f80481bcf017ca356940ab86dd5b2
Author: Frederic <frederic at MacFrederic.local>
Date: Thu May 7 15:05:09 2015 +0200
NML: adding the rico() Property read by Polish Doughnut
---
include/GyotoNumericalMetricLorene.h | 3 +++
lib/NumericalMetricLorene.C | 43 +++++++++++++++++++++++++++---------
lib/PolishDoughnut.C | 9 ++++++--
3 files changed, 42 insertions(+), 13 deletions(-)
diff --git a/include/GyotoNumericalMetricLorene.h b/include/GyotoNumericalMetricLorene.h
index 76a17ea..23c1441 100644
--- a/include/GyotoNumericalMetricLorene.h
+++ b/include/GyotoNumericalMetricLorene.h
@@ -73,6 +73,7 @@ class Gyoto::Metric::NumericalMetricLorene
Lorene::Scalar** lorentz_tab_; ///< Lorentz factor at surface (if any)
Lorene::Valeur** hor_tab_; ///< Apparent horizon (if any)
double risco_; ///< ISCO coordinate radius
+ double rico_; ///< Innermost circular orbit coordinate radius
double rmb_; ///< Marginally bound orbit coordinate radius
void free(); ///< deallocate memory
@@ -95,6 +96,8 @@ class Gyoto::Metric::NumericalMetricLorene
void initialTime(double t0);
double horizon() const ;
void horizon(double t0);
+ double rico() const ;
+ void rico(double r0);
bool hasSurface() const;
void hasSurface(bool s);
bool bosonstarcircular() const;
diff --git a/lib/NumericalMetricLorene.C b/lib/NumericalMetricLorene.C
index 511fcbc..8d7fb6e 100644
--- a/lib/NumericalMetricLorene.C
+++ b/lib/NumericalMetricLorene.C
@@ -39,6 +39,7 @@ GYOTO_PROPERTY_BOOL(NumericalMetricLorene,
GYOTO_PROPERTY_BOOL(NumericalMetricLorene, BosonStarCircular, NonBosonStarCircular, bosonstarcircular)
GYOTO_PROPERTY_DOUBLE(NumericalMetricLorene, Horizon, horizon)
GYOTO_PROPERTY_DOUBLE(NumericalMetricLorene, Time, initialTime)
+GYOTO_PROPERTY_DOUBLE(NumericalMetricLorene, Rico, rico)
GYOTO_PROPERTY_VECTOR_DOUBLE(NumericalMetricLorene,
RefineIntegStep, refineIntegStep)
// Keep File last here, so it is processed last in fillElement()
@@ -72,7 +73,8 @@ NumericalMetricLorene::NumericalMetricLorene() :
horizon_(0.),
risco_(0.),
rmb_(0.),
- bosonstarcircular_(false)
+ bosonstarcircular_(false),
+ rico_(0.)
{
GYOTO_DEBUG << endl;
}
@@ -100,7 +102,8 @@ NumericalMetricLorene::NumericalMetricLorene(const NumericalMetricLorene&o) :
horizon_(o.horizon_),
risco_(o.risco_),
rmb_(o.rmb_),
- bosonstarcircular_(o.bosonstarcircular_)
+ bosonstarcircular_(o.bosonstarcircular_),
+ rico_(o.rico_)
{
GYOTO_DEBUG << endl;
if (o.filename_) directory(o.filename_);
@@ -298,23 +301,34 @@ Valeur** NumericalMetricLorene::getHor_tab() const {
return hor_tab_;}
double NumericalMetricLorene::getRms() const {
GYOTO_DEBUG << endl;
- return risco_;}
+ if (rico()!=0.) return rico();
+ else return risco_;
+}
double NumericalMetricLorene::getRmb() const {
GYOTO_DEBUG << endl;
return rmb_;}
double NumericalMetricLorene::getSpecificAngularMomentum(double rr) const {
- // Computes the specific angular momentum \ell = -u_phi / u_t
+ // Computes the Keplerian specific angular momentum \ell = -u_phi / u_t
+ // for circular geodesics,
// for a general axisym metric in the equatorial plane
double pos[4]={0., rr, M_PI/2., 0.};
double gtt_dr =gmunu_up_dr(pos, 0, 0),
gtph_dr =gmunu_up_dr(pos, 0, 3),
gphph_dr=gmunu_up_dr(pos, 3, 3);
-
- return gtph_dr/gphph_dr +
- sqrt(gtph_dr/gphph_dr * gtph_dr/gphph_dr - gtt_dr/gphph_dr) ;
+
+ double lKep = gtph_dr/gphph_dr +
+ sqrt(gtph_dr/gphph_dr * gtph_dr/gphph_dr - gtt_dr/gphph_dr);
+
+ if (lKep!=lKep || lKep==lKep+1.){
+ cerr << "At r= " << rr << endl;
+ throwError("In NML::getSpecificAngMom: lKep not defined here!"
+ " You are probably below the innermost circular orbit.");
+ }
+
+ return lKep;
}
double NumericalMetricLorene::getPotential(double pos[4], double l_cst) const {
@@ -2077,6 +2091,9 @@ void NumericalMetricLorene::specifyMarginalOrbits(bool s) {
}
}
+double NumericalMetricLorene::rico() const {return rico_;}
+void NumericalMetricLorene::rico(double r0) {rico_=r0;}
+
bool NumericalMetricLorene::mapEt() const {return mapet_;}
void NumericalMetricLorene::mapEt(bool s) {
mapet_ = s;
@@ -2156,6 +2173,7 @@ void NumericalMetricLorene::circularVelocity(double const * coor,
double* vel,
double dir,
int indice_time) const {
+ //cout << "IN CIRCULAR" << endl;
if (bosonstarcircular_){
// This expression is related to the ZAMO 3-velocity derived
// in Grandclement+14 boson star paper
@@ -2191,14 +2209,17 @@ void NumericalMetricLorene::circularVelocity(double const * coor,
double ut = 1./(NN*sqrt(1.-Vzamo*Vzamo));
vel[0] = ut; vel[1] = 0.; vel[2] = 0.; vel[3] = Omega*ut;
- //double ell=2.5;
- //double pot = 0.5*log((g_tp*g_tp-g_tt*g_pp)
- ///(g_tt*ell*ell+2.*ell*g_tp+g_pp));
+ double ell=2.5;
+ double pot = 0.5*log((g_tp*g_tp-g_tt*g_pp)
+ /(g_tt*ell*ell+2.*ell*g_tp+g_pp));
//cout << rr << " " << g_tp*g_tp-g_tt*g_pp << " " << g_tt*ell*ell+2.*ell*g_tp+g_pp << " " << pot << endl;
//}
double normtol = 1e-6;
double norm = ScalarProd(coor,vel,vel);
- if (fabs(norm+1.)>normtol) throwError("In NML::circularv: bad norm");
+ if (fabs(norm+1.)>normtol) {
+ cerr << "At rr=" << coor[1] << endl;
+ throwError("In NML::circularv: bad norm");
+ }
return;
}
diff --git a/lib/PolishDoughnut.C b/lib/PolishDoughnut.C
index ef1fc26..5dac89c 100644
--- a/lib/PolishDoughnut.C
+++ b/lib/PolishDoughnut.C
@@ -288,12 +288,17 @@ void PolishDoughnut::angmomrinner(std::vector<double> const &v) {
//cout << "l0,rin= " << l0_ << " " << rintorus_ << endl;
double posin[4]={0.,rintorus_,M_PI/2.,0.};
W_surface_ = gg_->getPotential(posin,l0_);
- double rmin = rintorus_, rmax = 1000.;
+ double rmin;
+ double rico = gg_->getRms();
+ if (rico>0.) rmin=rico;
+ else rmin=rintorus_;
+ double rmax = 1000.;
+ cout << "rmin max= " << rmin << " " << rmax << endl;
r_centre_ = intersection.ridders(rmin, rmax) ;
double posc[4]={0.,r_centre_,M_PI/2.,0.};
W_centre_ = gg_->getPotential(posc,l0_);
DeltaWm1_ = 1./(W_centre_ - W_surface_);
- //cout << "Ws Wc rc= " << W_surface_ << " " << W_centre_ << " " << r_centre_ << endl;
+ cout << "Ws Wc rc= " << W_surface_ << " " << W_centre_ << " " << r_centre_ << endl;
GYOTO_IF_DEBUG;
GYOTO_DEBUG_EXPR(l0_);
GYOTO_DEBUG_EXPR(r_centre_);
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-astro/packages/gyoto.git
More information about the Debian-astro-commits
mailing list