Skip to content

Commit

Permalink
remove redundant variables; add comments
Browse files Browse the repository at this point in the history
  • Loading branch information
@rafaelab
rafaelab committed Nov 8, 2019
1 parent 831c8ec commit e382490
Showing 1 changed file with 14 additions and 28 deletions.
42 changes: 14 additions & 28 deletions grplinst.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,19 +25,17 @@ void PlasmaInstability::initTables() {
}

void PlasmaInstability::setIGMDensity(double rho) {
// Defines the density of the IGM (at z=z?, at z=0?)
// in units of m^-3
// Defines the density of the IGM at z=0, in units of m^-3.
densityIGM = rho;
}

void PlasmaInstability::setBeamLuminosity(double rho) {
// Defines the density of the beam (at z=z?, at z=0?)
// in units of m^-3
// Sets the beam luminosity in units of Watts.
luminosityBeam = rho;
}

void PlasmaInstability::setTemperature(double t) {
// Defines the temperature of the IGM (at z=z?, at z=0?)
// Defines the temperature of the IGM at z=0.
temperatureIGM = t;
}

Expand All @@ -55,6 +53,7 @@ void PlasmaInstability::setModel(char model) {
void PlasmaInstability::process(Candidate *candidate) const {

int id = candidate->current.getId();

// only works for electrons and positrons
if (fabs(id) != 11)
return;
Expand Down Expand Up @@ -103,12 +102,9 @@ double PlasmaInstability::coolingPower(double E, double z) const {

double PlasmaInstability::coolingPowerA(double E, double z) const {

double T = temperatureIGM;
double nIGM = densityIGM;
double L = luminosityBeam;
E /= (1 + z);
double eta = 1.;
double Ethr = 8.7e-6 * pow(1 + z, -13. / 6) * pow(L / 1e38, -1. / 3) * pow(nIGM / 0.1, 1. / 3);
double Ethr = 8.7e-6 * pow(1 + z, -13. / 6) * pow(luminosityBeam / 1e38, -1. / 3) * pow(densityIGM / 0.1, 1. / 3);

double a0, a1, a2, a3, a4;
if (E < Ethr) {
Expand All @@ -124,10 +120,11 @@ double PlasmaInstability::coolingPowerA(double E, double z) const {
a3 = 1. / 3;
a4 = 1. / 6;
}
return a0 * eta * pow(E / TeV, a2) * pow(L / 1e38, a3) * pow(nIGM / 0.1, a4);
return a0 * eta * pow(E / TeV, a2) * pow(luminosityBeam / 1e38, a3) * pow(densityIGM / 0.1, a4);
}

double PlasmaInstability::coolingPowerB(double E, double z) const {

E /= (1 + z);
double d0 = log10(redshift2LightTravelDistance(z) / Mpc); // co-moving?
double w = pow(10, interpolate(d0, _d, _w));
Expand All @@ -136,14 +133,10 @@ double PlasmaInstability::coolingPowerB(double E, double z) const {

double PlasmaInstability::coolingPowerC(double E, double z) const {

double T = temperatureIGM;
double nIGM = densityIGM;
double L = luminosityBeam;

E /= (1 + z);
double eta = 1.;
double Ethr = 7.9e-8 * pow(1 + z, -9. / 4) / sqrt(L / 1e38) * pow(nIGM / 0.1, 0.5) * pow(T / 1e4, 1.);
double F = 1. + 1.25 * log(T / 1e4) - 0.25 * log(nIGM / 0.1) + 0.5 * log(1 + z);
double Ethr = 7.9e-8 * pow(1 + z, -9. / 4) / sqrt(luminosityBeam / 1e38) * pow(densityIGM / 0.1, 0.5) * pow(temperatureIGM / 1e4, 1.);
double F = 1. + 1.25 * log(temperatureIGM / 1e4) - 0.25 * log(densityIGM / 0.1) + 0.5 * log(1 + z);

double a0, a1, a2, a3, a4, b;
if (E < Ethr) {
Expand All @@ -152,7 +145,7 @@ double PlasmaInstability::coolingPowerC(double E, double z) const {
a2 = -1;
a3 = 1. / 3.;
a4 = 5. / 6.;
b = pow(T / 1e4, 2);
b = pow(temperatureIGM / 1e4, 2);
} else {
a0 = 1.4e-23;
a1 = 11. / 3;
Expand All @@ -161,18 +154,14 @@ double PlasmaInstability::coolingPowerC(double E, double z) const {
a4 = 1. / 6.;
b = 1. / F;
}
return a0 * eta * pow(1 + z, a1) * pow(E / TeV, a2) * pow(L / 1e38, a3) * pow(nIGM / 0.1, a4) * b;
return a0 * eta * pow(1 + z, a1) * pow(E / TeV, a2) * pow(luminosityBeam / 1e38, a3) * pow(densityIGM / 0.1, a4) * b;
}

double PlasmaInstability::coolingPowerD(double E, double z) const {

double T = temperatureIGM;
double nIGM = densityIGM;
double L = luminosityBeam;

E /= (1 + z);
double eta = 1.; // for now fixed; should add a function to play with it.
double Ethr = 6.9e-6 * pow(1 + z, -13. / 16) * pow(L / 1e38, -1. / 3) / sqrt(nIGM / 0.1);
double Ethr = 6.9e-6 * pow(1 + z, -13. / 16) * pow(luminosityBeam / 1e38, -1. / 3) / sqrt(densityIGM / 0.1);

double a0, a1, a2, a3, a4, b;
if (E < Ethr) {
Expand All @@ -188,15 +177,12 @@ double PlasmaInstability::coolingPowerD(double E, double z) const {
a3 = 1. / 3;
a4 = 1. / 6;
}
return a0 * eta * pow(1 + z, a1) * pow(E / TeV, a2) * pow(L / 1e38, a3) * pow(nIGM / 0.1, a4);
return a0 * eta * pow(1 + z, a1) * pow(E / TeV, a2) * pow(luminosityBeam / 1e38, a3) * pow(densityIGM / 0.1, a4);
}

double PlasmaInstability::coolingPowerE(double E, double z) const {

double T = temperatureIGM;
double nIGM = densityIGM;
double L = luminosityBeam;
E /= (1 + z);
return 2.7e-20 * pow(0.5 + z / 2., 19. / 6) * E * pow(E / TeV, -1.) * pow(L / 1e38, 1. / 3) * pow(nIGM / 0.1, -1. / 3) * (T / 1e4);
return 2.7e-20 * pow(0.5 + z / 2., 19. / 6) * E * pow(E / TeV, -1.) * pow(luminosityBeam / 1e38, 1. / 3) * pow(densityIGM / 0.1, -1. / 3) * (temperatureIGM / 1e4);

}

0 comments on commit e382490

Please sign in to comment.