[Debian-astro-commits] [gyoto] 59/221: Implementing specific moving observers for KerrBL/Minkowski.
Thibaut Jean-Claude Paumard
thibaut at moszumanska.debian.org
Fri May 22 20:52:33 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 d0374b75093a82f6dfa419fac880a99ee70217cf
Author: Frederic <frederic at MacFrederic.local>
Date: Wed Oct 15 19:04:09 2014 +0200
Implementing specific moving observers for KerrBL/Minkowski.
Screen: getRayCoord calls Metric::observerTetrad; new keywords in XML
Metric: Metric::observerTetrad checks the normalization
KerrBL and Minkowski: observerTetrad implements Keplerian observer
KerrBL: obserTetrad implements ZAMO observer
---
include/GyotoKerrBL.h | 4 ++
include/GyotoMetric.h | 15 +++++
include/GyotoMinkowski.h | 6 +-
include/GyotoScreen.h | 9 ++-
lib/KerrBL.C | 68 ++++++++++++++++++++
lib/Metric.C | 28 ++++++++-
lib/Minkowski.C | 47 +++++++++++++-
lib/Screen.C | 158 ++++++++++++++++++++++++++++++++---------------
8 files changed, 279 insertions(+), 56 deletions(-)
diff --git a/include/GyotoKerrBL.h b/include/GyotoKerrBL.h
index 5467063..2dfc4b4 100644
--- a/include/GyotoKerrBL.h
+++ b/include/GyotoKerrBL.h
@@ -167,6 +167,10 @@ class Gyoto::Metric::KerrBL : public Metric::Generic {
void setParticleProperties(Worldline* line, const double* coord) const;
virtual int isStopCondition(double const * const coord) const;
+ void observerTetrad(std::string const obskind,
+ double const pos[4], double fourvel[4],
+ double screen1[4], double screen2[4],
+ double screen3[4]) const;
};
#endif
diff --git a/include/GyotoMetric.h b/include/GyotoMetric.h
index d80d257..a0af6b3 100644
--- a/include/GyotoMetric.h
+++ b/include/GyotoMetric.h
@@ -364,6 +364,21 @@ int MyKind::setParameter(std::string name, std::string content, std::string unit
std::string content,
std::string unit);
+ /**
+ * \brief Computes the orthonormal local tetrad of the observer
+ *
+ * \param obskind input: kind of observer (eg: "ZAMO","KeplerianObserver"...)
+ * \param pos input: position,
+ * \param fourvel output: observer 4-velocity (norm -1)
+ * \param screen1 output: first vector in the screen plane
+ * \param screen2 output: second vector in the screen plane
+ * \param screen3 output: vector normal to the screen
+ */
+ virtual void observerTetrad(std::string const obskind,
+ double const pos[4], double fourvel[4],
+ double screen1[4], double screen2[4],
+ double screen3[4]) const ;
+
// Outputs
#ifdef GYOTO_USE_XERCES
diff --git a/include/GyotoMinkowski.h b/include/GyotoMinkowski.h
index b23f5d0..edca0b2 100644
--- a/include/GyotoMinkowski.h
+++ b/include/GyotoMinkowski.h
@@ -65,8 +65,10 @@ class Gyoto::Metric::Minkowski
double gmunu(const double * x, int mu, int nu) const ;
double christoffel(const double coord[8], const int alpha, const int mu,
const int nu) const ;
-
-
+ void observerTetrad(std::string const obskind,
+ double const pos[4], double fourvel[4],
+ double screen1[4], double screen2[4],
+ double screen3[4]) const;
};
#endif
diff --git a/include/GyotoScreen.h b/include/GyotoScreen.h
index a893167..f40590f 100644
--- a/include/GyotoScreen.h
+++ b/include/GyotoScreen.h
@@ -205,6 +205,12 @@ class Gyoto::Screen : protected Gyoto::SmartPointee {
*/
double freq_obs_;
+ /**
+ * \brief What kind of observer are we considering? (At infinity, ZAMO...)
+ *
+ */
+ std::string observerkind_;
+
public:
// Constructors - Destructor
@@ -311,7 +317,8 @@ class Gyoto::Screen : protected Gyoto::SmartPointee {
* system. Content is copied.
*/
void setObserverPos(const double pos[4]);
-
+ void setObserverKind(const std::string kind);
+ std::string getObserverKind();
void setFourVel(const double coord[4]);
///< Sets the observer's 4-velocity
void setScreen1(const double coord[4]);
diff --git a/lib/KerrBL.C b/lib/KerrBL.C
index be13d7c..93a5007 100644
--- a/lib/KerrBL.C
+++ b/lib/KerrBL.C
@@ -1146,6 +1146,74 @@ void KerrBL::setParticleProperties(Worldline * line, const double* coord) const
line -> setCst(cst,5);
}
+void KerrBL::observerTetrad(string const obskind,
+ double const pos[4], double fourvel[4],
+ double screen1[4], double screen2[4],
+ double screen3[4]) const{
+
+ double gtt = gmunu(pos,0,0),
+ grr = gmunu(pos,1,1),
+ gthth = gmunu(pos,2,2),
+ gpp = gmunu(pos,3,3),
+ gtp = gmunu(pos,0,3);
+
+ if (obskind=="ZAMO"){
+ double fact = gtt*gpp - gtp*gtp;
+ if (fact == 0.){
+ throwError("In KerrBL::observerTetrad: "
+ "bad values");
+ }
+ double ut2 = -gpp/fact;
+ if (grr==0. || gthth==0. || gpp==0. || ut2<=0.){
+ throwError("In KerrBL::observerTetrad: "
+ "bad values");
+ }
+ double omega = -gtp/gpp;
+ double ut = sqrt(ut2);
+ double fourv[4]={ut,0.,0.,omega*ut};
+ double e3[4] = {0.,-1./sqrt(grr),0.,0.};
+ double e2[4] = {0.,0.,-1./sqrt(gthth),0.};
+ double e1[4] = {0.,0.,0.,-1/sqrt(gpp)};
+
+ for (int ii=0;ii<4;ii++){
+ fourvel[ii]=fourv[ii];
+ screen1[ii]=e1[ii];
+ screen2[ii]=e2[ii];
+ screen3[ii]=e3[ii];
+ }
+
+ }else if (obskind=="KeplerianObserver"){
+ double omega = 1./(pow(pos[1],1.5)+spin_);
+ double ut2 = -1/(gtt+gpp*omega*omega+2.*gtp*omega);
+ if (ut2 <= 0. || grr<=0. || gthth <=0.) {
+ throwError("In KerrBL::observerTetrad: "
+ "bad values");
+ }
+ double ut = sqrt(ut2);
+ double fourv[4]={ut,0.,0.,omega*ut};
+ double e3[4] = {0.,-1./sqrt(grr),0.,0.};
+ double e2[4] = {0.,0.,-1./sqrt(gthth),0.};
+
+ double fact1 = (gtp+gpp*omega)/(gtt+gtp*omega),
+ fact2 = gtt*fact1*fact1+gpp-2.*gtp*fact1;
+ if (fact2 <= 0.) throwError("In KerrBL::observerTetrad: "
+ "bad values");
+ double a2 = 1./sqrt(fact2), a1 = -a2*fact1;
+ double e1[4] = {-a1,0.,0.,-a2};
+
+ for (int ii=0;ii<4;ii++){
+ fourvel[ii]=fourv[ii];
+ screen1[ii]=e1[ii];
+ screen2[ii]=e2[ii];
+ screen3[ii]=e3[ii];
+ }
+ }else{
+ throwError("In KerrBL::observerTetrad "
+ "unknown observer kind");
+ }
+ Generic::observerTetrad(obskind,pos,fourvel,screen1,screen2,screen3);
+}
+
#ifdef GYOTO_USE_XERCES
void KerrBL::fillElement(Gyoto::FactoryMessenger *fmp) {
fmp -> setParameter("Spin", spin_);
diff --git a/lib/Metric.C b/lib/Metric.C
index 80d5404..7279b33 100644
--- a/lib/Metric.C
+++ b/lib/Metric.C
@@ -490,7 +490,33 @@ void Metric::Generic::setParticleProperties(Worldline*, const double*) const {
# endif
}
-
+void Metric::Generic::observerTetrad(string const obskind,
+ double const coord[4], double fourvel[4],
+ double screen1[4], double screen2[4],
+ double screen3[4]) const{
+ // No general way to define the tetrad, should be defined
+ // in specific metrics. Test below will obviously fail for
+ // a machine-initialized tetrad.
+ double normtol=1e-10;
+ if (fabs(ScalarProd(coord,fourvel,fourvel)+1.)>normtol ||
+ fabs(ScalarProd(coord,screen1,screen1)-1.)>normtol ||
+ fabs(ScalarProd(coord,screen2,screen2)-1.)>normtol ||
+ fabs(ScalarProd(coord,screen3,screen3)-1.)>normtol){
+ cout << "norm= " << ScalarProd(coord,fourvel,fourvel) << " " << ScalarProd(coord,screen1,screen1) << " " << ScalarProd(coord,screen2,screen2) << " " << ScalarProd(coord,screen3,screen3) << endl;
+ throwError("In Metric:observerTetrad: observer's local"
+ " basis is not properly normalized");
+ }
+
+ if (fabs(ScalarProd(coord,fourvel,screen1))>normtol ||
+ fabs(ScalarProd(coord,fourvel,screen2))>normtol ||
+ fabs(ScalarProd(coord,fourvel,screen3))>normtol ||
+ fabs(ScalarProd(coord,screen1,screen2))>normtol ||
+ fabs(ScalarProd(coord,screen1,screen3))>normtol ||
+ fabs(ScalarProd(coord,screen2,screen3))>normtol){
+ throwError("In Metric:observerTetrad: observer's local"
+ " basis is not orthogonal");
+ }
+}
/***************For SmartPointers**************/
diff --git a/lib/Minkowski.C b/lib/Minkowski.C
index e3180e7..c8167fb 100644
--- a/lib/Minkowski.C
+++ b/lib/Minkowski.C
@@ -89,7 +89,7 @@ int Minkowski::christoffel(double dst[4][4][4], const double pos[8]) const {
// second form is given below, as an example.
double Minkowski::gmunu(const double * pos, int mu, int nu) const {
if (mu<0 || nu<0 || mu>3 || nu>3)
- throwError ("KerrKS::gmunu: incorrect value for mu or nu");
+ throwError ("Minkowski::gmunu: incorrect value for mu or nu");
if (mu!=nu) return 0.;
if (mu==0) return -1.;
@@ -158,6 +158,51 @@ double Minkowski::christoffel(const double pos[8], const int alpha, const int mm
}
+void Minkowski::observerTetrad(string const obskind,
+ double const pos[4], double fourvel[4],
+ double screen1[4], double screen2[4],
+ double screen3[4]) const{
+ if (coordKind()!=GYOTO_COORDKIND_SPHERICAL){
+ throwError("In Minkowski::observerTetrad: "
+ "coordinates should be spherical-like");
+ }
+ if (obskind=="KeplerianObserver"){
+ double gtt = gmunu(pos,0,0),
+ grr = gmunu(pos,1,1),
+ gthth = gmunu(pos,2,2),
+ gpp = gmunu(pos,3,3);
+ double omega = 1./(pow(pos[1],1.5));
+ double ut2 = -1/(gtt+gpp*omega*omega);
+ if (ut2 <= 0. || grr<=0. || gthth <=0.) {
+ throwError("In Minkowski::observerTetrad: "
+ "bad values");
+ }
+ double ut = sqrt(ut2);
+ double fourv[4]={ut,0.,0.,omega*ut};
+ double e3[4] = {0.,-1./sqrt(grr),0.,0.};
+ double e2[4] = {0.,0.,-1./sqrt(gthth),0.};
+
+ double fact1 = gpp*omega/gtt,
+ fact2 = gtt*fact1*fact1+gpp;
+ if (fact2 <= 0.) throwError("In Minkowski::observerTetrad: "
+ "bad values");
+ double a2 = 1./sqrt(fact2), a1 = -a2*fact1;
+ double e1[4] = {-a1,0.,0.,-a2};
+
+ for (int ii=0;ii<4;ii++){
+ fourvel[ii]=fourv[ii];
+ screen1[ii]=e1[ii];
+ screen2[ii]=e2[ii];
+ screen3[ii]=e3[ii];
+ }
+ }else{
+ throwError("In Minkowski::observerTetrad "
+ "unknown observer kind");
+ }
+ Generic::observerTetrad(obskind,pos,fourvel,screen1,screen2,screen3);
+}
+
+
// Fillelement is required to be able to export the Metric to an XML file.
#ifdef GYOTO_USE_XERCES
diff --git a/lib/Screen.C b/lib/Screen.C
index de4c181..77d68c3 100644
--- a/lib/Screen.C
+++ b/lib/Screen.C
@@ -21,6 +21,7 @@
#include "GyotoFactoryMessenger.h"
#include "GyotoScreen.h"
#include "GyotoConverters.h"
+#include "GyotoKerrBL.h"
#include <iostream>
#include <cstdlib>
@@ -49,7 +50,8 @@ Screen::Screen() :
distance_(1.), dmax_(GYOTO_SCREEN_DMAX), anglekind_(0),
alpha0_(0.), delta0_(0.),
gg_(NULL), spectro_(NULL),
- freq_obs_(1.)
+ freq_obs_(1.),
+ observerkind_("ObserverAtInfinity")
{
# if GYOTO_DEBUG_ENABLED
GYOTO_DEBUG_EXPR(dmax_);
@@ -72,7 +74,8 @@ Screen::Screen(const Screen& o) :
dmax_(o.dmax_),
anglekind_(o.anglekind_),
alpha0_(o.alpha0_), delta0_(o.delta0_),
- gg_(NULL), spectro_(NULL), freq_obs_(o.freq_obs_)
+ gg_(NULL), spectro_(NULL), freq_obs_(o.freq_obs_),
+ observerkind_(o.observerkind_)
{
if (o.gg_()) gg_=o.gg_->clone();
if (o.spectro_()) spectro_ = o.spectro_ -> clone();
@@ -276,6 +279,13 @@ void Screen::setObserverPos(const double coord[4]) {
computeBaseVectors();
}
+void Screen::setObserverKind(const string kind) {
+ observerkind_=kind;
+}
+string Screen::getObserverKind() {
+ return observerkind_;
+}
+
void Screen::setFourVel(const double coord[4]) {
for (int ii=0;ii<4;ii++)
fourvel_[ii]=coord[ii];
@@ -387,7 +397,7 @@ void Screen::getRayCoord(double alpha, double delta,
{
alpha+=alpha0_; delta+=delta0_; // Screen orientation
-
+ double normtol=1e-10;
int i; // dimension : 0, 1, 2
double pos[4];
# if GYOTO_DEBUG_ENABLED
@@ -489,7 +499,7 @@ void Screen::getRayCoord(double alpha, double delta,
// 4-vector tangent to photon geodesic
- if (fourvel_[0]==0.){
+ if (fourvel_[0]==0. && observerkind_=="ObserverAtInfinity"){
/*
---> Observer local frame not given in XML <---
Assume observer static at infinity ("standard Gyoto")
@@ -548,7 +558,7 @@ void Screen::getRayCoord(double alpha, double delta,
// 0-component of photon tangent 4-vector found by normalizing
gg_ -> nullifyCoord(coord);
- }else{
+ }else{
/*
---> Observer local frame given in XML <---
Express photon tangent 4-vector in the observer basis
@@ -560,52 +570,94 @@ void Screen::getRayCoord(double alpha, double delta,
if (coord[2]==0. || coord[2]==M_PI)
throwError("Please move Screen away from z-axis");
}
-
- //Check orthonormal basis
- double normtol=1e-10;
- if (fabs(gg_->ScalarProd(coord,fourvel_,fourvel_)+1.)>normtol ||
- fabs(gg_->ScalarProd(coord,screen1_,screen1_)-1.)>normtol ||
- fabs(gg_->ScalarProd(coord,screen2_,screen2_)-1.)>normtol ||
- fabs(gg_->ScalarProd(coord,screen3_,screen3_)-1.)>normtol){
- throwError("In Screen:getRayCoord: observer's local"
- " basis is not properly normalized");
+
+ if (fourvel_[0]!=0. && observerkind_!="ObserverAtInfinity"){
+ throwError("In Screen:getRayCoord: "
+ " choose an implemented observer kind OR"
+ " explicitly give the local tetrad in the XML");
+ }
+
+ if (fourvel_[0]==0){
+ // Implemented observer specifid in XML, local tetrad computed by Metric
+ const double fourpos[4]={coord[0],coord[1],coord[2],coord[3]};
+ double fourvel[4], screen1[4], screen2[4], screen3[4];
+ gg_ -> observerTetrad(observerkind_,fourpos,fourvel,screen1,
+ screen2,screen3);
+
+ /*cout << "Vectors in Screen: " << setprecision(17) << endl;
+ for (int ii=0;ii<4;ii++) cout << fourvel[ii] << " ";
+ cout << endl;
+ for (int ii=0;ii<4;ii++) cout << screen1[ii] << " ";
+ cout << endl;
+ for (int ii=0;ii<4;ii++) cout << screen2[ii] << " ";
+ cout << endl;
+ for (int ii=0;ii<4;ii++) cout << screen3[ii] << " ";
+ cout << endl;*/
+
+ /* Photon tagent 4-vector l defined by: l = p + fourvel_
+ where p gives the direction of the photon
+ in the observer's rest space (orthogonal to fourvel).
+ Here we choose the particular tangent 4-vector l
+ that satisfies l.fourvel_=-1
+ Then l = fourvel_ + (orthogonal proj of l onto rest space)
+ = fourvel_ + p
+ and p = vel[0]*screen1_ + vel[1]*screen2_ + vel[2]*screen3_
+ with p.p = 1, thus l.l = 0 as it should
+ */
+
+ coord[4]=vel[0]*screen1[0]
+ +vel[1]*screen2[0]
+ +vel[2]*screen3[0]
+ +fourvel[0];
+ coord[5]=vel[0]*screen1[1]
+ +vel[1]*screen2[1]
+ +vel[2]*screen3[1]
+ +fourvel[1];
+ coord[6]=vel[0]*screen1[2]
+ +vel[1]*screen2[2]
+ +vel[2]*screen3[2]
+ +fourvel_[2];
+ coord[7]=vel[0]*screen1[3]
+ +vel[1]*screen2[3]
+ +vel[2]*screen3[3]
+ +fourvel[3];
+ }else{
+ // Local tetrad given by the user in the XML file. Check it.
+ if (fabs(gg_->ScalarProd(coord,fourvel_,fourvel_)+1.)>normtol ||
+ fabs(gg_->ScalarProd(coord,screen1_,screen1_)-1.)>normtol ||
+ fabs(gg_->ScalarProd(coord,screen2_,screen2_)-1.)>normtol ||
+ fabs(gg_->ScalarProd(coord,screen3_,screen3_)-1.)>normtol){
+ cout << "norm= " << gg_->ScalarProd(coord,fourvel_,fourvel_) << " " << gg_->ScalarProd(coord,screen1_,screen1_) << " " << gg_->ScalarProd(coord,screen2_,screen2_) << " " << gg_->ScalarProd(coord,screen3_,screen3_) << endl;
+ throwError("In Screen:getRayCoord: observer's local"
+ " basis is not properly normalized");
+ }
+
+ if (fabs(gg_->ScalarProd(coord,fourvel_,screen1_))>normtol ||
+ fabs(gg_->ScalarProd(coord,fourvel_,screen2_))>normtol ||
+ fabs(gg_->ScalarProd(coord,fourvel_,screen3_))>normtol ||
+ fabs(gg_->ScalarProd(coord,screen1_,screen2_))>normtol ||
+ fabs(gg_->ScalarProd(coord,screen1_,screen3_))>normtol ||
+ fabs(gg_->ScalarProd(coord,screen2_,screen3_))>normtol)
+ throwError("In Screen:getRayCoord: observer's local"
+ " basis is not orthogonal");
+
+ coord[4]=vel[0]*screen1_[0]
+ +vel[1]*screen2_[0]
+ +vel[2]*screen3_[0]
+ +fourvel_[0];
+ coord[5]=vel[0]*screen1_[1]
+ +vel[1]*screen2_[1]
+ +vel[2]*screen3_[1]
+ +fourvel_[1];
+ coord[6]=vel[0]*screen1_[2]
+ +vel[1]*screen2_[2]
+ +vel[2]*screen3_[2]
+ +fourvel_[2];
+ coord[7]=vel[0]*screen1_[3]
+ +vel[1]*screen2_[3]
+ +vel[2]*screen3_[3]
+ +fourvel_[3];
}
-
- if (fabs(gg_->ScalarProd(coord,fourvel_,screen1_))>normtol ||
- fabs(gg_->ScalarProd(coord,fourvel_,screen2_))>normtol ||
- fabs(gg_->ScalarProd(coord,fourvel_,screen3_))>normtol ||
- fabs(gg_->ScalarProd(coord,screen1_,screen2_))>normtol ||
- fabs(gg_->ScalarProd(coord,screen1_,screen3_))>normtol ||
- fabs(gg_->ScalarProd(coord,screen2_,screen3_))>normtol)
- throwError("In Screen:getRayCoord: observer's local"
- " basis is not orthogonal");
-
- /* Photon tagent 4-vector l defined by: l = p + fourvel_
- where p gives the direction of the photon
- in the observer's rest space (orthogonal to fourvel).
- Here we choose the particular tangent 4-vector l
- that satisfies l.fourvel_=-1
- Then l = fourvel_ + (orthogonal proj of l onto rest space)
- = fourvel_ + p
- and p = vel[0]*screen1_ + vel[1]*screen2_ + vel[2]*screen3_
- with p.p = 1, thus l.l = 0 as it should
- */
- coord[4]=vel[0]*screen1_[0]
- +vel[1]*screen2_[0]
- +vel[2]*screen3_[0]
- +fourvel_[0];
- coord[5]=vel[0]*screen1_[1]
- +vel[1]*screen2_[1]
- +vel[2]*screen3_[1]
- +fourvel_[1];
- coord[6]=vel[0]*screen1_[2]
- +vel[1]*screen2_[2]
- +vel[2]*screen3_[2]
- +fourvel_[2];
- coord[7]=vel[0]*screen1_[3]
- +vel[1]*screen2_[3]
- +vel[2]*screen3_[3]
- +fourvel_[3];
if (fabs(gg_->ScalarProd(coord,coord+4,coord+4))>normtol){
throwError("In Screen::getRayCoord: "
" tangent 4-vector to photon not properly normalized");
@@ -1079,7 +1131,11 @@ SmartPointer<Screen> Screen::Subcontractor(FactoryMessenger* fmp) {
# endif
if (name=="Time") {tobs_tmp = atof(tc); tunit=unit; tobs_found=1;}
else if (name=="Position") {for (int i=0;i<4;++i) pos[i] = strtod(tc, &tc);
- scr -> setObserverPos (pos); }
+ scr -> setObserverPos (pos);
+ }
+ else if (name=="KeplerianObserver" || name=="ZAMO") {
+ scr -> setObserverKind(name);
+ }
else if (name=="FourVelocity") {
fourvel_found=1;
for (int i=0;i<4;++i) pos[i] = strtod(tc, &tc);
--
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