Skip to content
Snippets Groups Projects
Commit 07172572 authored by Oleg Kalashev's avatar Oleg Kalashev
Browse files

coord_type = long double is used for particle coordinates

parent 425d7b90
Branches
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
src/app/crbeam/CRbeam.cpp
+ 1
− 0
View file @ 07172572
......@@ -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);
......
src/lib/Deflection3D.cpp
+ 17
− 17
View file @ 07172572
......@@ -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);
......
src/lib/Deflection3D.h
+ 4
− 4
View file @ 07172572
......@@ -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,
......
src/lib/Particle.cpp
+ 1
− 1
View file @ 07172572
......@@ -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);
}
......
src/lib/Particle.h
+ 10
− 4
View file @ 07172572
......@@ -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;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment