diff --git a/src/app/crbeam/CRbeam.cpp b/src/app/crbeam/CRbeam.cpp
index f475bfb0032bd14cb483176f4332cc9c1f1d275b..697428352f4e4eef789d6556bd1bf70082828057 100644
--- a/src/app/crbeam/CRbeam.cpp
+++ b/src/app/crbeam/CRbeam.cpp
@@ -234,6 +234,7 @@ CRbeam::CRbeam(int argc, char** argv):
fTauPrint(false),
cmd(argc,argv,commands)
{
+ std::cout << "sizeof(coord_type) = " << sizeof(coord_type) << std::endl;
if( cmd.has_param("-h") )
{
cmd.printHelp(std::cerr);
diff --git a/src/lib/Deflection3D.cpp b/src/lib/Deflection3D.cpp
index 3febb8daf09a8917140a10cbc05051cf6b986b38..d24de774f2a5dd84cb48529ec518c8e30bacc086 100644
--- a/src/lib/Deflection3D.cpp
+++ b/src/lib/Deflection3D.cpp
@@ -52,7 +52,7 @@ namespace Interactions {
{
//dn/dt = qeB/E, where q is charge in units of electron charge; e=sqrt(alpha)-electron charge
double mult = aParticle.ElectricCharge()*units.Gauss*units.e/aParticle.Energy;
- const double* N = aParticle.Pdir;
+ const coord_type* N = aParticle.Pdir;
ASSERT(outRate.size()==3);
outRate[0] = (N[1]* aBgauss[2]- aBgauss[1]*N[2])*mult;
outRate[1] = (N[2]* aBgauss[0]- aBgauss[2]*N[0])*mult;
@@ -88,28 +88,28 @@ namespace Interactions {
dx = deltaT/nSteps;
std::vector<double> deflRate(3);
std::vector<double> curB(3);
- double newPdir[3];
+ coord_type newPdir[3];
int totSteps=0;//debug info
for(u_int stepB = 0; stepB<nSteps; stepB++){
pB->GetValueGauss(aParticle.X, aParticle.Time, curB);
GetRotationRate(aParticle, curB, deflRate);
- double absRate = sqrt(deflRate[0]*deflRate[0]+deflRate[1]*deflRate[1]+deflRate[2]*deflRate[2]);
+ cosmo_time absRate = sqrt(deflRate[0]*deflRate[0]+deflRate[1]*deflRate[1]+deflRate[2]*deflRate[2]);
cosmo_time maxStep = 0.05*M_PI/absRate/fAccuracy;//max step corresponds to 9/fAccuracy degrees turn
u_int nSubSteps = (u_int)(dx/maxStep+1.);
cosmo_time dXsubStep = dx/nSubSteps;
for(u_int substep=0; substep<nSubSteps; substep++){
if(substep)
GetRotationRate(aParticle, curB, deflRate);
- double absPdir = 0.;
+ cosmo_time absPdir = 0.;
for(int i=0; i<3; i++){
- double val = aParticle.Pdir[i] + deflRate[i]*dXsubStep;
+ coord_type val = aParticle.Pdir[i] + deflRate[i]*dXsubStep;
newPdir[i] = val;
absPdir += val*val;
}
absPdir = sqrt(absPdir);//norm of new vector
for(int i=0; i<3; i++){
- double val = newPdir[i]/absPdir;//make sure rotated vector has unit norm
+ coord_type val = newPdir[i]/absPdir;//make sure rotated vector has unit norm
aParticle.X[i] += 0.5*(aParticle.Pdir[i] + val)*beta*dXsubStep;
aParticle.Pdir[i] = val;
}
@@ -175,7 +175,7 @@ namespace Interactions {
fAmplitude[2] = std::complex<double>(-cosAlpha* fSinTheta, 0)* fAbsBgauss;
}
- void MonochromaticMF::GetValueGauss(const double *x, const CosmoTime &aTime,
+ void MonochromaticMF::GetValueGauss(const coord_type *x, const CosmoTime &aTime,
std::vector<double> &outValue) const
{
ASSERT(outValue.size()==3);
@@ -194,12 +194,12 @@ namespace Interactions {
int MagneticField::Print(double lambdaMpc, std::ostream& aOutput) {
std::vector<double> B(3);
CosmoTime t(0);
- double step=lambdaMpc*units.Mpc/20.;
- double L=lambdaMpc*units.Mpc*3;
- for(double x=0.; x<=L; x+=step) {
- for (double y = 0.; y <= L; y+=step) {
- for (double z = 0.; z <= L; z += step) {
- const double r[] = {x, y, z};
+ coord_type step=lambdaMpc*units.Mpc/20.;
+ coord_type L=lambdaMpc*units.Mpc*3;
+ for(coord_type x=0.; x<=L; x+=step) {
+ for (coord_type y = 0.; y <= L; y+=step) {
+ for (coord_type z = 0.; z <= L; z += step) {
+ const coord_type r[] = {x, y, z};
GetValueGauss(r, t, B);
aOutput << x / units.Mpc << "\t" << y / units.Mpc << "\t" << z / units.Mpc << "\t"
<< B[0] << "\t" << B[1] << "\t" << B[2] << "\n";
@@ -212,7 +212,7 @@ namespace Interactions {
return 0;
}
- double MagneticField::CorrelationLength(std::vector<double> x, std::vector<double> aDirection, const CosmoTime &aTime,
+ double MagneticField::CorrelationLength(std::vector<coord_type> x, std::vector<coord_type> aDirection, const CosmoTime &aTime,
double aMaxValMpc) {
double step = 0.1*MinVariabilityScale(aTime);
int max_steps = (int)(aMaxValMpc*units.Mpc/step + 0.5);
@@ -280,7 +280,7 @@ namespace Interactions {
return 0;
}
- void TurbulentMF::GetValueGauss(const double *x, const CosmoTime &aTime, std::vector<double> &outValue) const {
+ void TurbulentMF::GetValueGauss(const coord_type *x, const CosmoTime &aTime, std::vector<double> &outValue) const {
std::vector<double> wave(3, 0.);
outValue.assign(3, 0.);
for(std::vector<SafePtr<MonochromaticMF> >::const_iterator pw = fWaves.begin(); pw!=fWaves.end(); pw++){
@@ -323,8 +323,8 @@ namespace Interactions {
double TurbulentMF::MeanCorLength(Randomizer &aRandomizer, double aLmin_Mpc, double aLmax_Mpc, double aVariance_Gauss,
int aNmodes, const CosmoTime &aTime, int aNsamples, double aMaxValMpc){
double l_cor_sum = 0.;
- auto x = std::vector<double>({0.,0.,0.});
- auto dx = std::vector<double>({0.,0.,1.});
+ auto x = std::vector<coord_type>({0.,0.,0.});
+ auto dx = std::vector<coord_type>({0.,0.,1.});
for(int i=0; i<aNsamples; i++){
TurbulentMF B(aRandomizer, aLmin_Mpc, aLmax_Mpc, aVariance_Gauss, aNmodes);
l_cor_sum += B.CorrelationLength(x, dx, aTime, aMaxValMpc);
diff --git a/src/lib/Deflection3D.h b/src/lib/Deflection3D.h
index dd1dcfa35c8cf0c826088870fc85f547b5f0df54..28f71b03f6074c24c236ef6dec27b0996daad2ac 100644
--- a/src/lib/Deflection3D.h
+++ b/src/lib/Deflection3D.h
@@ -39,13 +39,13 @@ namespace Interactions {
class MagneticField : public SmartReferencedObj{
public:
- virtual void GetValueGauss(const double *x, const CosmoTime &aTime, std::vector<double>& outValue) const = 0;
+ virtual void GetValueGauss(const coord_type *x, const CosmoTime &aTime, std::vector<double>& outValue) const = 0;
virtual double MinVariabilityScale(const CosmoTime &aTime) const = 0;
int Print(double lambdaMpc, std::ostream& aOutput=std::cout);
- double CorrelationLength(std::vector<double> x, std::vector<double> aDirection, const CosmoTime &aTime,
+ double CorrelationLength(std::vector<coord_type> x, std::vector<coord_type> aDirection, const CosmoTime &aTime,
double aMaxValMpc=10000);
};
@@ -84,7 +84,7 @@ namespace Interactions {
MonochromaticMF(double aLamdbaMpc, double aAmplitudeGauss, double aBetaK,
double aAlphaK, double aTheta, double aPhi);
- virtual void GetValueGauss(const double *x, const CosmoTime &aTime,
+ virtual void GetValueGauss(const coord_type *x, const CosmoTime &aTime,
std::vector<double> &outValue) const;
virtual double MinVariabilityScale(const CosmoTime &aTime) const;
@@ -111,7 +111,7 @@ namespace Interactions {
public:
TurbulentMF(Randomizer &aRandomizer, double aLmin_Mpc, double aLmax_Mpc, double aVariance_Gauss, int aNmodes);
- virtual void GetValueGauss(const double *x, const CosmoTime &aTime, std::vector<double> &outValue) const;
+ virtual void GetValueGauss(const coord_type *x, const CosmoTime &aTime, std::vector<double> &outValue) const;
virtual double MinVariabilityScale(const CosmoTime &aTime) const;
static int UnitTest();
static double MeanCorLength(Randomizer &aRandomizer, double aLmin_Mpc, double aLmax_Mpc, double aVariance_Gauss,
diff --git a/src/lib/Particle.cpp b/src/lib/Particle.cpp
index 62f28a05d294206c78bfeed813dbab871b62a409..1a1e1a9b8bbf22cc4d8d1ad0b2233ee12f5a2350 100644
--- a/src/lib/Particle.cpp
+++ b/src/lib/Particle.cpp
@@ -144,7 +144,7 @@ const char* Particle::ParticleNames[ParticleTypeEOF] =
void Particle::PropagateFreely(cosmo_time dt)
{
Time += dt;
- double b=beta();
+ coord_type b=beta();
for(int i=0; i<3; i++)
X[i]+=(b*Pdir[i]*dt);
}
diff --git a/src/lib/Particle.h b/src/lib/Particle.h
index 54a5700e7aa8472c3e302ff73029f694819e9aa9..773adceadf4bb9cad3d18b8db87cd51408957cbe 100644
--- a/src/lib/Particle.h
+++ b/src/lib/Particle.h
@@ -108,6 +108,9 @@ namespace mcray
EndNuclei,
ParticleTypeEOF = EndNuclei //must be the last
};
+
+ typedef cosmo_time coord_type;
+
#define ParticleCustomDataSize 16
class Particle
{
@@ -132,8 +135,11 @@ namespace mcray
inline Particle& operator=(const Particle& aParticle) {
memcpy(this, &aParticle, sizeof(Particle)); return *this;
}
- inline double beta() const {
- double m=Mass(); return sqrt(1.-m*m/Energy/Energy);
+ inline coord_type beta() const {
+ coord_type m=Mass();
+ coord_type gamma_1 = m/Energy;
+ coord_type beta2 = 1. - gamma_1*gamma_1;
+ return sqrt(beta2);
}
inline void SetParent(const Particle& aParticle){
@@ -152,8 +158,8 @@ namespace mcray
CosmoTime Time;
double Energy;
long Ninteractions;//number of interactions in the interaction chain
- double X[3];//comoving coordinates, source location: (0,0,0)
- double Pdir[3];//momentum direction, initial value: (0,0,1)
+ coord_type X[3];//comoving coordinates, source location: (0,0,0)
+ coord_type Pdir[3];//momentum direction, initial value: (0,0,1)
unsigned long long id;
unsigned long long fPrevId;