commit/enzo-dev: jwise77: Merged in aemerick/enzo-dev-constants (pull request #417)
8 views
Skip to first unread message
commits...@bitbucket.org
Feb 21, 2019, 4:43:33 AM2/21/19
to enzo...@googlegroups.com
1 new commit in enzo-dev:
https://bitbucket.org/enzo/enzo-dev/commits/7d8fe75451e0/
Changeset: 7d8fe75451e0
Branch: week-of-code
User: jwise77
Date: 2019-02-21 09:43:17+00:00
Summary: Merged in aemerick/enzo-dev-constants (pull request #417)
Moving all local constants to those in header file
Approved-by: Nathan Goldbaum <ngol...@illinois.edu>
Approved-by: Brian O'Shea <os...@msu.edu>
Approved-by: John Regan <johnanth...@gmail.com>
Affected #: 194 files
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/AMRH5writer.C
--- a/src/enzo/AMRH5writer.C
+++ b/src/enzo/AMRH5writer.C
@@ -24,6 +24,7 @@
#include "typedefs.h"
#include "AMRH5writer.h"
#include "StarParticleData.h"
+#include "phys_constants.h"
AMRHDF5Writer::AMRHDF5Writer() :
h5DataType(H5T_NATIVE_FLOAT),
@@ -151,14 +152,14 @@
error = true;
}
- // The 3.14159 header is to test the endian since this is a binary
+ // The pi header is to test the endian since this is a binary
// file. The root cell spacing is necessary to calculate the dx for
// all other grids.
- float myPi = 3.14159;
char nFields8bit = (char) NumberOfAllFields;
char stag8bit = (char) stag;
int field;
- fwrite(&myPi, sizeof(float), 1, index_file);
+ double test_pi = pi;
+ fwrite(&test_pi, sizeof(float), 1, index_file);
fwrite(&root_dx, sizeof(float), 1, index_file);
fwrite(&nFields8bit, sizeof(char), 1, index_file);
fwrite(&stag8bit, sizeof(char), 1, index_file);
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/AdiabaticExpansionInitialize.C
--- a/src/enzo/AdiabaticExpansionInitialize.C
+++ b/src/enzo/AdiabaticExpansionInitialize.C
@@ -29,6 +29,7 @@
#include "Hierarchy.h"
#include "TopGridData.h"
#include "CosmologyParameters.h"
+#include "phys_constants.h"
/* Set the mean molecular mass as in Grid_ComputeTemperatureField.C */
@@ -135,7 +136,7 @@
ENZO_FAIL("Error in GetUnits.");
}
PressureUnits = DensityUnits * (LengthUnits/TimeUnits)*(LengthUnits/TimeUnits);
- MagneticUnits = sqrt(PressureUnits*4.0*M_PI);
+ MagneticUnits = sqrt(PressureUnits*4.0*pi);
CRUnits = DensityUnits * (LengthUnits/TimeUnits)*(LengthUnits/TimeUnits);
for (int dim = 0; dim < MAX_DIMENSION; dim++)
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/AdjustMustRefineParticlesRefineToLevel.C
--- a/src/enzo/AdjustMustRefineParticlesRefineToLevel.C
+++ b/src/enzo/AdjustMustRefineParticlesRefineToLevel.C
@@ -41,6 +41,7 @@
#include "TopGridData.h"
#include "LevelHierarchy.h"
#include "Star.h"
+#include "phys_constants.h"
int GetUnits(float *DensityUnits, float *LengthUnits,
float *TemperatureUnits, float *TimeUnits,
@@ -61,7 +62,7 @@
ENZO_FAIL("Your setup of assymmetrical topgrid is never tested.\n");
}
- const double pc = 3.086e18, ln2 = 0.69314718;
+ const double ln2 = 0.69314718;
int MustRefineParticlesRefineToLevel_prev = MustRefineParticlesRefineToLevel;
/* Compute Units. */
@@ -80,7 +81,7 @@
MustRefineParticlesRefineToLevel =
(int) floor( log( LengthUnits / MetaData->TopGridDims[0]
- / MustRefineParticlesRefineToLevelAutoAdjust / pc) / ln2);
+ / MustRefineParticlesRefineToLevelAutoAdjust / pc_cm) / ln2);
if (MyProcessorNumber == ROOT_PROCESSOR) {
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/AnalysisBaseClass.C
--- a/src/enzo/AnalysisBaseClass.C
+++ b/src/enzo/AnalysisBaseClass.C
@@ -15,6 +15,7 @@
#include "TopGridData.h"
#include "CosmologyParameters.h"
#include "AnalysisBaseClass.h"
+#include "phys_constants.h"
std::map<HierarchyEntry *, int> OriginalGridID;
#ifdef USE_HDF5_GROUPS
@@ -31,16 +32,8 @@
HierarchyEntry *topgrid,
int level,
FLOAT *Left,
- FLOAT *Right):
-#ifdef USE_BAD_VALUES
- G(6.67e-8),
- MPC_CM(3.0856e+24),
- MSOLAR_G(1.989e+33)
-#else
- G(6.6742e-8),
- MPC_CM(3.0856776e+24),
- MSOLAR_G(1.989e+33)
-#endif
+ FLOAT *Right)
+
{
int i;
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/AnalysisBaseClass.h
--- a/src/enzo/AnalysisBaseClass.h
+++ b/src/enzo/AnalysisBaseClass.h
@@ -33,10 +33,6 @@
protected:
- const float G;
- const float MPC_CM;
- const float MSOLAR_G;
-
TopGridData *MetaData;
HierarchyEntry *TopGrid;
int MaximumAnalysisLevel;
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/CalculateSubtractionParameters.C
--- a/src/enzo/CalculateSubtractionParameters.C
+++ b/src/enzo/CalculateSubtractionParameters.C
@@ -29,6 +29,7 @@
#include "TopGridData.h"
#include "LevelHierarchy.h"
#include "CommunicationUtilities.h"
+#include "phys_constants.h"
int CommunicationBroadcastValue(int *Value, int BroadcastProcessor);
@@ -43,8 +44,7 @@
float &Radius, double &Subtraction)
{
- const double pc = 3.086e18, Msun = 1.989e33, Grav = 6.673e-8, yr = 3.1557e7, Myr = 3.1557e13,
- k_b = 1.38e-16, m_h = 1.673e-24, c = 3.0e10, sigma_T = 6.65e-25, h=0.70;
+ const double h=0.70;
float mdot, AccretedMass, SafetyFactor;
float MassEnclosed = 0, Metallicity = 0, ColdGasMass = 0, OneOverRSquaredSum, AvgVelocity[MAX_DIMENSION];
@@ -129,7 +129,7 @@
density = star_CurrentGrid->ReturnBaryonField(DensNum)[index];
if (MultiSpecies == 0) {
- number_density = density * DensityUnits / (Mu * m_h);
+ number_density = density * DensityUnits / (Mu * mh);
mu = Mu;
} else {
number_density =
@@ -157,7 +157,7 @@
because only integer can be broadcasted at the moment;
anyway, this shouldn't make a huge difference overall. */
- c_s = (int)(sqrt(Gamma * k_b * temperature[index] / (mu * m_h)));
+ c_s = (int)(sqrt(Gamma * kboltz * temperature[index] / (mu * mh)));
fprintf(stdout, "index = %d, temp = %g, mu = %g, density = %g, number_density = %g, c_s = %d\n",
index, temperature[index], mu, density, number_density, c_s);
@@ -182,7 +182,7 @@
/* When calculating Bondi radius, let's not bother to get accurate temperature and mu
just use the MBHAccretionFixedTemperature and default mu */
- c_s = (int)(sqrt(Gamma * k_b * MBHAccretionFixedTemperature / (Mu * m_h)));
+ c_s = (int)(sqrt(Gamma * kboltz * MBHAccretionFixedTemperature / (Mu * mh)));
#endif
@@ -192,14 +192,14 @@
old_mass = (float)(star_mass);
mdot = isnan(star_last_accretion_rate) ? 0.0 : star_last_accretion_rate;
- AccretedMass = mdot * dtForThisStar * TimeUnits; //in Msun
+ AccretedMass = mdot * dtForThisStar * TimeUnits; //in SolarMass
if (MBHAccretionRadius > 0)
- Radius = MBHAccretionRadius * pc / LengthUnits;
+ Radius = MBHAccretionRadius * pc_cm / LengthUnits;
else {
SafetyFactor = -MBHAccretionRadius;
- Radius = SafetyFactor * 2.0 * Grav * old_mass * Msun / (c_s * c_s) / LengthUnits;
+ Radius = SafetyFactor * 2.0 * GravConst * old_mass * SolarMass / (c_s * c_s) / LengthUnits;
}
Radius = min(max(Radius, 4*StarLevelCellWidth), 100*StarLevelCellWidth);
@@ -284,8 +284,8 @@
here Subtraction is the density that should be subtracted.
if one wants to use this, change Grid_SAMFS.C accordingly. */
- Subtraction = -AccretedMass * Msun /
- (4*M_PI/3.0 * pow(Radius*LengthUnits, 3)) / DensityUnits;
+ Subtraction = -AccretedMass * SolarMass /
+ (4*pi/3.0 * pow(Radius*LengthUnits, 3)) / DensityUnits;
#endif
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/ClusterSMBHSumGasMass.C
--- a/src/enzo/ClusterSMBHSumGasMass.C
+++ b/src/enzo/ClusterSMBHSumGasMass.C
@@ -84,7 +84,7 @@
float ColdGasMassMsun=ClusterSMBHColdGasMass*MassUnits/SolarMass;
if (ClusterSMBHCalculateGasMass >1 ){ //2-calculate&remove,3-calculate&remove&re-orient,4-Bondi
if (ClusterSMBHCalculateGasMass == 3){
- ClusterSMBHJetDim = floor(Time*TimeUnits/(1.0e6*3.1557e7*ClusterSMBHJetPrecessionPeriod)); //ClusterSMBHJetPrecessionPeriod is now the Dim changing period.
+ ClusterSMBHJetDim = floor(Time*TimeUnits/(1.0e6*yr_s*ClusterSMBHJetPrecessionPeriod)); //ClusterSMBHJetPrecessionPeriod is now the Dim changing period.
}
if (ColdGasMassMsun > 0.000001) {
ClusterSMBHFeedbackSwitch = TRUE;
@@ -94,7 +94,7 @@
}
else {
ClusterSMBHJetMdot = (ColdGasMassMsun/(ClusterSMBHAccretionTime*1e6))/2.0; // AccretionTime from Myr to yr; reset Mdot, still in Msun/yr. Devide it by 2 because Mdot is for only one jet.
- ClusterSMBHJetEdot = (ClusterSMBHAccretionEpsilon*ClusterSMBHJetMdot * SolarMass/3.1557e7) * POW(clight,2)/1.0e44; //for one jet
+ ClusterSMBHJetEdot = (ClusterSMBHAccretionEpsilon*ClusterSMBHJetMdot * SolarMass/yr_s) * POW(clight,2)/1.0e44; //for one jet
}
}
else
@@ -108,7 +108,7 @@
ClusterSMBHFeedbackSwitch = (ColdGasMassMsun <= ClusterSMBHEnoughColdGas && LastClusterSMBHFeedbackSwitch == FALSE) ? FALSE : TRUE;
if (LastClusterSMBHFeedbackSwitch == FALSE && ClusterSMBHFeedbackSwitch == TRUE) {
- ClusterSMBHStartTime = Time + ClusterSMBHTramp*1.0e6*3.1557e7/TimeUnits;
+ ClusterSMBHStartTime = Time + ClusterSMBHTramp*1.0e6*yr_s/TimeUnits;
if (ClusterSMBHJetPrecessionPeriod < 0.00001 & ClusterSMBHJetPrecessionPeriod > -0.00001) //if precession off (set to 0), change the angle of the jets
ClusterSMBHJetAnglePhi += 0.5;
if (ClusterSMBHJetPrecessionPeriod < -0.00001) //if precession negative, change the jet dimension
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/CommunicationMergeStarParticle.C
--- a/src/enzo/CommunicationMergeStarParticle.C
+++ b/src/enzo/CommunicationMergeStarParticle.C
@@ -33,6 +33,7 @@
#include "TopGridData.h"
#include "Hierarchy.h"
#include "LevelHierarchy.h"
+#include "phys_constants.h"
/* function prototypes */
@@ -148,8 +149,7 @@
GetUnits(&DensityUnits, &LengthUnits, &TemperatureUnits, &TimeUnits,
&VelocityUnits, &MassUnits, 1.0);
- double Msun = 1.989e33;
- if(UsePhysicalUnit) MassUnits = DensityUnits*pow(LengthUnits,3)/Msun;
+ if(UsePhysicalUnit) MassUnits = DensityUnits*pow(LengthUnits,3)/SolarMass;
/*
if (debug) {
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/ConductionBubbleInitialize.C
--- a/src/enzo/ConductionBubbleInitialize.C
+++ b/src/enzo/ConductionBubbleInitialize.C
@@ -29,6 +29,7 @@
#include "Grid.h"
#include "Hierarchy.h"
#include "TopGridData.h"
+#include "phys_constants.h"
int GetUnits(float *DensityUnits, float *LengthUnits,
float *TemperatureUnits, float *TimeUnits,
@@ -108,7 +109,7 @@
total energy (that's not done automatically, since we don't know for sure
what the user is going to do). */
if (HydroMethod == MHD_RK){
- float MagneticUnits = sqrt(DensityUnits*4.0*M_PI)*VelocityUnits;
+ float MagneticUnits = sqrt(DensityUnits*4.0*pi)*VelocityUnits;
ConductionBubbleInitialUniformBField[0] /= MagneticUnits;
ConductionBubbleInitialUniformBField[1] /= MagneticUnits;
ConductionBubbleInitialUniformBField[2] /= MagneticUnits;
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/ConductionCloudInitialize.C
--- a/src/enzo/ConductionCloudInitialize.C
+++ b/src/enzo/ConductionCloudInitialize.C
@@ -25,6 +25,7 @@
#include "Grid.h"
#include "Hierarchy.h"
#include "TopGridData.h"
+#include "phys_constants.h"
int GetUnits(float *DensityUnits, float *LengthUnits,
float *TemperatureUnits, float *TimeUnits,
@@ -123,11 +124,11 @@
return FAIL;
}
- float Boltzmann = 1.38e-16, mu = 0.6, mh=1.67e-24;
+ float mu = 0.6;
if(ConductionCloudTemperature > 1.0){
- ConductionCloudTotalEnergy = (Boltzmann*ConductionCloudTemperature)/((Gamma - 1.0)*mu*mh);
+ ConductionCloudTotalEnergy = (kboltz*ConductionCloudTemperature)/((Gamma - 1.0)*mu*mh);
ConductionCloudTotalEnergy /= (VelocityUnits*VelocityUnits);
printf("ConductionCloudTotalEnergy is %e and ConductionCloudTemperature is %e\n\n",ConductionCloudTotalEnergy, ConductionCloudTemperature);
fflush(stdout);
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/ConductionTestInitialize.C
--- a/src/enzo/ConductionTestInitialize.C
+++ b/src/enzo/ConductionTestInitialize.C
@@ -27,6 +27,7 @@
#include "Grid.h"
#include "Hierarchy.h"
#include "TopGridData.h"
+#include "phys_constants.h"
int GetUnits(float *DensityUnits, float *LengthUnits,
float *TemperatureUnits, float *TimeUnits,
@@ -105,12 +106,12 @@
return FAIL;
}
- float Boltzmann = 1.38e-16, mu = 0.6, mh=1.67e-24;
+ float mu = 0.6;
// User can set
if(ConductionTestTemperature > 1.0){
- ConductionTestTotalEnergy = (Boltzmann*ConductionTestTemperature)/((Gamma - 1.0)*mu*mh);
+ ConductionTestTotalEnergy = (kboltz*ConductionTestTemperature)/((Gamma - 1.0)*mu*mh);
ConductionTestTotalEnergy /= (VelocityUnits*VelocityUnits);
ConductionTestGasEnergy = ConductionTestTotalEnergy;
printf("ConductionTestTotalEnergy is %e and ConductionTestTemperature is %e\n\n",ConductionTestTotalEnergy, ConductionTestTemperature);
@@ -118,7 +119,7 @@
}
if (HydroMethod == MHD_RK){
- float MagneticUnits = sqrt(DensityUnits*4.0*M_PI)*VelocityUnits;
+ float MagneticUnits = sqrt(DensityUnits*4.0*pi)*VelocityUnits;
ConductionTestInitialUniformBField[0] /= MagneticUnits;
ConductionTestInitialUniformBField[1] /= MagneticUnits;
ConductionTestInitialUniformBField[2] /= MagneticUnits;
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/CosmoIonizationInitialize.C
--- a/src/enzo/CosmoIonizationInitialize.C
+++ b/src/enzo/CosmoIonizationInitialize.C
@@ -29,6 +29,7 @@
#include "Hierarchy.h"
#include "LevelHierarchy.h"
#include "TopGridData.h"
+#include "phys_constants.h"
// default constants
#define DEFAULT_MU 0.6 // mean molecular mass
@@ -134,8 +135,7 @@
// convert input temperature to internal energy
RadHydroTemperature = max(RadHydroTemperature,MIN_TEMP); // enforce minimum
- float mp = 1.67262171e-24; // proton mass [g]
- float kb = 1.3806504e-16; // boltzmann constant [erg/K]
+
float nH, HI, HII, nHe, HeI, HeII, HeIII, ne, num_dens, mu;
if (RadHydroChemistry == 1) {
HI = 1.0 - RadHydroInitialFractionHII;
@@ -157,7 +157,7 @@
mu = 1.0/num_dens;
}
// compute the internal energy
- float RadHydroIEnergy = kb*RadHydroTemperature/mu/mp/(Gamma-1.0);
+ float RadHydroIEnergy = kboltz*RadHydroTemperature/mu/mh/(Gamma-1.0);
// set up the grid(s) on this level
HierarchyEntry *Temp = &TopGrid;
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/CosmologyGetUnits.C
--- a/src/enzo/CosmologyGetUnits.C
+++ b/src/enzo/CosmologyGetUnits.C
@@ -39,6 +39,7 @@
#include "ErrorExceptions.h"
#include "macros_and_parameters.h"
#include "CosmologyParameters.h"
+#include "phys_constants.h"
/* function prototypes */
@@ -66,7 +67,7 @@
*DensityUnits = 1.8788e-29*OmegaMatterNow*POW(HubbleConstantNow,2)*
POW(1 + CurrentRedshift,3);
- *LengthUnits = 3.085678e24*ComovingBoxSize/HubbleConstantNow/
+ *LengthUnits = Mpc_cm*ComovingBoxSize/HubbleConstantNow/
(1 + CurrentRedshift);
*TemperatureUnits = 1.81723e6*POW(ComovingBoxSize,2)*OmegaMatterNow*
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/CosmologySimulationInitialize.C
--- a/src/enzo/CosmologySimulationInitialize.C
+++ b/src/enzo/CosmologySimulationInitialize.C
@@ -39,6 +39,7 @@
#include "CosmologyParameters.h"
#include "fortran.def"
#include "CommunicationUtilities.h"
+#include "phys_constants.h"
// Function prototypes
@@ -419,7 +420,7 @@
ENZO_FAIL("Error in GetUnits.");
}
PressureUnits = DensityUnits * VelocityUnits*VelocityUnits;
- MagneticUnits = sqrt(PressureUnits*4.0*M_PI);
+ MagneticUnits = sqrt(PressureUnits*4.0*pi);
for (int dim = 0; dim < MAX_DIMENSION; dim++) {
if (CosmologySimulationInitialUniformBField[dim] != 0.0 && HydroMethod != 4 && HydroMethod != 6)
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/ExternalBoundary_SetWengenCollidingFlowBoundary.C
--- a/src/enzo/ExternalBoundary_SetWengenCollidingFlowBoundary.C
+++ b/src/enzo/ExternalBoundary_SetWengenCollidingFlowBoundary.C
@@ -27,6 +27,7 @@
#include "Grid.h"
#include "MHD2DTestGlobalData.h"
#include "hydro_rk/EOS.h"
+#include "phys_constants.h"
/* function prototypes */
@@ -117,7 +118,7 @@
EOSSoundSpeed*EOSSoundSpeed*rho ;
EOS(pres, rho, eintl, h, cs, dpdrho, dpde, 0, 1); // compute eintl
ramp = 1./(1.+exp(-2/RampWidth*(y-0.5)));
- f = cos(2.*M_PI*x*10.)*exp(-fabs(y-0.5)*10.)*cos(2.*M_PI*x*3);
+ f = cos(2.*pi*x*10.)*exp(-fabs(y-0.5)*10.)*cos(2.*pi*x*3);
vx = f * (vxl + ramp*(vxu-vxl)) ;
vy = vyl + ramp*(vyu - vyl);
bx = (Bxl+ ramp*(Bxu-Bxl)) ;
@@ -152,7 +153,7 @@
EOSSoundSpeed*EOSSoundSpeed*rho ;
EOS(pres, rho, eintl, h, cs, dpdrho, dpde, 0, 1); // compute eintl
ramp = 1./(1.+exp(-2/RampWidth*(y-0.5)));
- f = cos(2.*M_PI*x*10.)*exp(-fabs(y-0.5)*10.)*cos(2.*M_PI*x*3);
+ f = cos(2.*pi*x*10.)*exp(-fabs(y-0.5)*10.)*cos(2.*pi*x*3);
vx = f * (vxl + ramp*(vxu-vxl)) ;
vy = vyl + ramp*(vyu - vyl);
bx = (Bxl+ ramp*(Bxu-Bxl)) ;
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/FOF.C
--- a/src/enzo/FOF.C
+++ b/src/enzo/FOF.C
@@ -26,6 +26,7 @@
#include "macros_and_parameters.h"
#include "typedefs.h"
#include "global_data.h"
+#include "phys_constants.h"
#include "Fluxes.h"
#include "GridList.h"
@@ -148,9 +149,9 @@
void set_units(FOFData &AllVars) {
- AllVars.UnitLength_in_cm = 3.085678e21;
- AllVars.UnitMass_in_g = 1.989e43;
- AllVars.UnitVelocity_in_cm_per_s = 1.0e5;
+ AllVars.UnitLength_in_cm = kpc_cm;
+ AllVars.UnitMass_in_g = SolarMass * 1.0E10;
+ AllVars.UnitVelocity_in_cm_per_s = km_cm;
AllVars.Theta = 0.8; /* opening angle for potential computation */
AllVars.DesDensityNgb = 32;
@@ -1395,7 +1396,7 @@
FILE *in, *out;
int i, Ngroups, *Npart_g;
float *CM, *Mass;
- double SolarMass = UnitMass_in_g / 1.989e33;
+ double Msun = UnitMass_in_g / SolarMass;
float RootBoxSize[3];
if (MyProcessorNumber == 0) {
@@ -1427,7 +1428,7 @@
"Halo", "Npart", "Mass [Msun]", "CM(x)", "CM(y)", "CM(z)");
for (i = 0; i < Ngroups; i++)
fprintf(out, "%8d %12d %12.6g %12.6f %12.6f %12.6f\n",
- i, Npart_g[i], Mass[i]*SolarMass,
+ i, Npart_g[i], Mass[i]*Msun,
CM[3*i+0] / RootBoxSize[0] + leftEdge[0],
CM[3*i+1] / RootBoxSize[1] + leftEdge[1],
CM[3*i+2] / RootBoxSize[2] + leftEdge[2]);
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/FOF_Initialize.C
--- a/src/enzo/FOF_Initialize.C
+++ b/src/enzo/FOF_Initialize.C
@@ -33,6 +33,7 @@
#include "Hierarchy.h"
#include "LevelHierarchy.h"
#include "CosmologyParameters.h"
+#include "phys_constants.h"
#include "FOF_allvars.h"
#include "FOF_nrutil.h"
@@ -82,13 +83,13 @@
D.Time = a / (1 + InitialRedshift);
}
else {
- D.BoxSize = LengthUnits / 3.086e21;
+ D.BoxSize = LengthUnits / kpc_cm;
D.Time = 1.0;
}
// Critical density in units of Msun / kpc^3
D.RhoCritical0 = 1.4775867e31 *
- ((3 * pow(100 * HubbleConstantNow / 3.086e19, 2)) / (8 * M_PI * GRAVITY));
+ ((3.0 * POW(100 * HubbleConstantNow / 3.086e19, 2)) / (8.0 * pi * GRAVITY));
//D.RhoCritical /= pow(D.Time, 3);
// Sometimes MassUnits is infinite (in cgs) when using single
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/FOF_ngbtree.C
--- a/src/enzo/FOF_ngbtree.C
+++ b/src/enzo/FOF_ngbtree.C
@@ -12,7 +12,7 @@
#include "macros_and_parameters.h"
#include "typedefs.h"
#include "global_data.h"
-
+#include "phys_constants.h"
#include "FOF_allvars.h"
#include "FOF_ngbtree.h"
@@ -54,7 +54,7 @@
break;
}
- sr = th->len * pow((3.0/(4*M_PI)*SECFACTOR) * desngb /
+ sr = th->len * pow((3.0/(4*pi)*SECFACTOR) * desngb /
((float)(th->count)), 1.0/3);
} // ENDELSE
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/FOF_properties.C
--- a/src/enzo/FOF_properties.C
+++ b/src/enzo/FOF_properties.C
@@ -6,7 +6,7 @@
#include "macros_and_parameters.h"
#include "typedefs.h"
#include "global_data.h"
-
+#include "phys_constants.h"
#include "FOF_allvars.h"
#include "FOF_proto.h"
@@ -134,7 +134,7 @@
// kpc^3, as is rho_crit.
e2 = D.OmegaLambda + D.Omega * pow(D.Time, -3.0);
// 1/a^3 factor converts comoving radius to proper.
- factor = 1e10 / (4*M_PI/3.0) * pow(D.Time, -3.0);
+ factor = 1e10 / (4*pi/3.0) * pow(D.Time, -3.0);
rho178 = 178.0 * D.RhoCritical0 * e2;
len4 = len/4;
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/FSProb_ComputeOpacityLW.C
--- a/src/enzo/FSProb_ComputeOpacityLW.C
+++ b/src/enzo/FSProb_ComputeOpacityLW.C
@@ -21,12 +21,12 @@
************************************************************************/
#ifdef TRANSFER
#include "FSProb.h"
+#include "phys_constants.h"
int FSProb::ComputeOpacityLW(float *H2Density)
{
const float H2ISigma = 3.71e-18; // cm^-2
- const float mh = 1.673e-24;
double factor = H2ISigma * (DenUnits/mh);
int i, j, k, dim, size;
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/FreeExpansionInitialize.C
--- a/src/enzo/FreeExpansionInitialize.C
+++ b/src/enzo/FreeExpansionInitialize.C
@@ -33,6 +33,7 @@
#include "Grid.h"
#include "Hierarchy.h"
#include "TopGridData.h"
+#include "phys_constants.h"
int GetUnits(float *DensityUnits, float *LengthUnits,
float *TemperatureUnits, float *TimeUnits,
@@ -128,7 +129,7 @@
&VelocityUnits, 0.0);
PressureUnits = DensityUnits * (LengthUnits/TimeUnits)*(LengthUnits/TimeUnits);
- MagneticUnits = sqrt(PressureUnits*4.0*M_PI);
+ MagneticUnits = sqrt(PressureUnits*4.0*pi);
for (dim = 0; dim < MAX_DIMENSION; dim++)
FreeExpansionBField[dim] /= MagneticUnits;
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/FreezeRateData.C
--- a/src/enzo/FreezeRateData.C
+++ b/src/enzo/FreezeRateData.C
@@ -30,6 +30,7 @@
#include "ExternalBoundary.h"
#include "Grid.h"
#include "Hierarchy.h"
+#include "phys_constants.h"
/* function prototypes */
int GetUnits(float *DensityUnits, float *LengthUnits,
@@ -120,9 +121,7 @@
// put the temperature together
float mu = rho/num_density;
- float mp = 1.67262171e-24; // proton mass [g]
- float kb = 1.3806504e-16; // boltzmann constant [erg/K]
- float Temp = max((Gamma-1.0)*mu*mp*eint/kb, 1.0);
+ float Temp = max((Gamma-1.0)*mu*mh*eint/kboltz, 1.0);
// find temperature bin
float lamT = 3.15614e5/Temp;
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/Grid_AddExternalAcceleration.C
--- a/src/enzo/Grid_AddExternalAcceleration.C
+++ b/src/enzo/Grid_AddExternalAcceleration.C
@@ -25,7 +25,8 @@
#include "ExternalBoundary.h"
#include "Grid.h"
#include "SphericalInfall.h"
-
+#include "phys_constants.h"
+
/* function prototypes */
int CosmologyComputeExpansionFactor(FLOAT time, FLOAT *a, FLOAT *dadt);
@@ -39,7 +40,7 @@
if (ProblemType == 24) {
- float Pi = 3.14159, accel;
+ float accel;
FLOAT a = 1, dadt, rcubed, xpos, ypos = 0, zpos = 0;
int i, j, k, n = 0;
@@ -82,7 +83,7 @@
(i.e. 1/a^2 * a = 1/a). */
accel = GravitationalConstant*SphericalInfallFixedMass/
- (4.0*Pi*rcubed*a);
+ (4.0*pi*rcubed*a);
/* Apply force. */
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/Grid_AddFeedbackSphere.C
--- a/src/enzo/Grid_AddFeedbackSphere.C
+++ b/src/enzo/Grid_AddFeedbackSphere.C
@@ -28,6 +28,7 @@
#include "Grid.h"
#include "Hierarchy.h"
#include "CosmologyParameters.h"
+#include "phys_constants.h"
#ifdef CONFIG_BFLOAT_4
#define TOLERANCE 1e-06
@@ -51,8 +52,6 @@
{
const float WhalenMaxVelocity = 35; // km/s
- const double Msun = 1.989e33;
- const double c = 3.0e10;
int dim, i, j, k, index;
int sx, sy, sz;
@@ -117,7 +116,7 @@
if (MetalNum > 0 && SNColourNum > 0 && cstar->type == PopIII)
MetalNum = SNColourNum;
- float BubbleVolume = (4.0 * M_PI / 3.0) * radius * radius * radius;
+ float BubbleVolume = (4.0 * pi / 3.0) * radius * radius * radius;
/***********************************************************************
SUPERNOVAE
@@ -144,7 +143,7 @@
// Correct if the volume with 27 cells is larger than the energy bubble volume
#ifdef UNUSED
float BoxVolume = 27 * CellWidth[0][0] * CellWidth[0][0] * CellWidth[0][0];
- float BubbleVolume = (4.0 * M_PI / 3.0) * radius * radius * radius;
+ float BubbleVolume = (4.0 * pi / 3.0) * radius * radius * radius;
//printf("BoxVolume = %lg, BubbleVolume = %lg\n", BoxVolume, BubbleVolume);
if (BoxVolume > BubbleVolume) {
//printf("Reducing ejecta density by %g\n", BubbleVolume / BoxVolume);
@@ -174,7 +173,7 @@
// printf("grid::AFS: before: cstar->Mass = %lf\n", cstar->Mass);
if (cstar->FeedbackFlag != SUPERNOVA) {
float old_mass = (float)(cstar->Mass);
- cstar->Mass -= EjectaDensity * DensityUnits * BubbleVolume * pow(LengthUnits,3.0) / Msun;
+ cstar->Mass -= EjectaDensity * DensityUnits * BubbleVolume * pow(LengthUnits,3.0) / SolarMass;
float frac = old_mass / cstar->Mass;
cstar->vel[0] *= frac;
cstar->vel[1] *= frac;
@@ -466,7 +465,7 @@
metallicity_inside = 0.0, colour_inside = 0.0, rho_inside, rho_metal_inside, rho_colour_inside;
float m_cell_edge = 0.0, metal_cell_edge = 0.0, colour_cell_edge = 0.0,
metallicity_edge = 0.0, colour_edge = 0.0, rho_jet, rho_metal_jet, rho_colour_jet;
- float L_x, L_y, L_z, L_s, nx_L = 0.0, ny_L = 0.0, nz_L = 0.0, costheta = cos(3.1415926/3.9);
+ float L_x, L_y, L_z, L_s, nx_L = 0.0, ny_L = 0.0, nz_L = 0.0, costheta = cos(pi/3.9);
float EjectaMass, EjectaMetalMass, MBHJetsVelocity;
if (cstar->FeedbackFlag == MBH_JETS) {
@@ -492,7 +491,7 @@
/* Calculate mass that has accumulated thus far; since EjectaDensity is calculated
for isotropic MBH_THERMAL, we use it only to compute EjectaMass, not apply it directly. */
- cstar->NotEjectedMass += EjectaDensity * DensityUnits * BubbleVolume * pow(LengthUnits,3.0) / Msun;
+ cstar->NotEjectedMass += EjectaDensity * DensityUnits * BubbleVolume * pow(LengthUnits,3.0) / SolarMass;
// fprintf(stdout, "d1, d2, d3, i, j, k = %d %d %d / %d %d %d\n",
// GridDimension[0], GridDimension[1], GridDimension[2], i,j,k);
@@ -526,7 +525,7 @@
/* Find ejecta mass */
- EjectaMass = cstar->NotEjectedMass * Msun / DensityUnits / pow(LengthUnits,3.0);
+ EjectaMass = cstar->NotEjectedMass * SolarMass / DensityUnits / pow(LengthUnits,3.0);
EjectaMetalMass = EjectaMass * MBHFeedbackMetalYield;
OutputWhenJetsHaveNotEjected = FALSE;
@@ -561,8 +560,8 @@
//add random noise to theta and phi
srand(time(NULL));
- theta += MaximumNoiseAngle * M_PI / 180.0 * ((2.0*(float)rand()/((float)(RAND_MAX)+(float)(1))) - 1.0);
- phi += MaximumNoiseAngle * M_PI / 180.0 * ((2.0*(float)rand()/((float)(RAND_MAX)+(float)(1))) - 1.0);
+ theta += MaximumNoiseAngle * pi / 180.0 * ((2.0*(float)rand()/((float)(RAND_MAX)+(float)(1))) - 1.0);
+ phi += MaximumNoiseAngle * pi / 180.0 * ((2.0*(float)rand()/((float)(RAND_MAX)+(float)(1))) - 1.0);
//get back to Cartesian coordinate; some tricks needed to preserve the signs of nx_L and ny_L
nx_L = sign(nx_L) * fabs(cos(theta))*sin(phi);
@@ -620,7 +619,7 @@
} // END jj-direction
} // END kk-direction
-// printf("EjectaM in Msun = %g, EjectaM = %g, EjectaMetalM = %g, m_cell_edge = %g, n_cell_edge = %d\n",
+// printf("EjectaM in SolarMass = %g, EjectaM = %g, EjectaMetalM = %g, m_cell_edge = %g, n_cell_edge = %d\n",
// cstar->NotEjectedMass, EjectaMass, EjectaMetalMass, m_cell_edge, n_cell_edge);
/* Calculate the jet density */
@@ -654,13 +653,13 @@
Note that EjectaThermalEnergy is never used; MBHFeedbackEnergyCoupling
should now be calculated considering gravitational redshift (Kim et al. 2010) */
- MBHJetsVelocity = c * sqrt( 2 * MBHFeedbackEnergyCoupling * MBHFeedbackRadiativeEfficiency
+ MBHJetsVelocity = clight * sqrt( 2 * MBHFeedbackEnergyCoupling * MBHFeedbackRadiativeEfficiency
/ MBHFeedbackMassEjectionFraction ) / VelocityUnits;
- if (MBHJetsVelocity * VelocityUnits > 0.99*c) {
+ if (MBHJetsVelocity * VelocityUnits > 0.99*clight) {
ENZO_VFAIL("grid::AddFS: MBHJetsVelocity is ultra-relativistic! (%g/ %g/ %g/ %g c)\n",
MBHFeedbackEnergyCoupling, MBHFeedbackRadiativeEfficiency,
- MBHFeedbackMassEjectionFraction, MBHJetsVelocity * VelocityUnits / c);
+ MBHFeedbackMassEjectionFraction, MBHJetsVelocity * VelocityUnits / clight);
}
/* Finally, add the jet feedback at the edges (outer part of the supercell) */
@@ -757,7 +756,7 @@
// old_mass, cstar->Mass, old_vel1, cstar->vel[1]);
fprintf(stdout, "grid::AddFS: jets injected (EjectaM = %g Ms, JetsVelocity = %g, grid v = %g, rho_jet = %g) along n_L = (%g, %g, %g)\n",
- EjectaMass * DensityUnits * pow(LengthUnits,3.0) / Msun,
+ EjectaMass * DensityUnits * pow(LengthUnits,3.0) / SolarMass,
MBHJetsVelocity, BaryonField[Vel3Num][n_cell_edge-1], rho_jet, nx_L, ny_L, nz_L);
} // END MBH_JETS
@@ -854,7 +853,7 @@
float kph, kheat;
sigma_HI *= (double) TimeUnits / ((double)LengthUnits * (double)LengthUnits)
- / (4.0 * M_PI);
+ / (4.0 * pi);
kph = (float) (Q_HI * sigma_HI);
kheat = kph * deltaE;
@@ -988,7 +987,7 @@
/* Now it's done, unmark. */
//cstar->FeedbackFlag = NO_FEEDBACK;
-
+
return SUCCESS;
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/Grid_AddH2DissociationFromSources.C
--- a/src/enzo/Grid_AddH2DissociationFromSources.C
+++ b/src/enzo/Grid_AddH2DissociationFromSources.C
@@ -29,12 +29,12 @@
#include "Hierarchy.h"
#include "CosmologyParameters.h"
#include "Star.h"
+#include "phys_constants.h"
#ifdef FUTUREAP
#include "ActiveParticle.h"
#include "ActiveParticle_RadiationParticle.h"
#include "ActiveParticle_SmartStar.h"
#endif
-#include "phys_constants.h"
#define THRESHOLD_DENSITY_DB36 1e14
#define THRESHOLD_DENSITY_DB37 5e14
@@ -59,6 +59,7 @@
int DensNum, GENum, TENum, Vel1Num, Vel2Num, Vel3Num, B1Num, B2Num, B3Num;
int DeNum, HINum, HIINum, HeINum, HeIINum, HeIIINum, HMNum, H2INum, H2IINum,
DINum, DIINum, HDINum;
+
double l_char = 0.0, N_H2 = 0.0;
double H2ISigma = 3.71e-18;
@@ -130,7 +131,7 @@
// where the gas is fully molecular. If the solver still breaks,
// then this parameter can be safely increased to ~100 AU or turning
// on H2 self-shielding with RadiationShield = 2.
- float dilutionRadius = 4.848e-6 * pc / (double) LengthUnits; // 1 AU
+ float dilutionRadius = 4.848e-6 * pc_cm / (double) LengthUnits; // 1 AU
float dilRadius2 = dilutionRadius * dilutionRadius;
float LightTravelDist = TimeUnits * clight / LengthUnits;
@@ -196,9 +197,9 @@
int TemperatureField = 0;
/* Pre-compute some quantities to speed things up */
TemperatureField = this->GetTemperatureFieldNumberForH2Shield();
- kdiss_r2 = (float) (LWLuminosity * H2ISigma / (4.0 * M_PI));
- kph_hm = (float) (IRLuminosity * HMSigma / (4.0 * M_PI));
- kdiss_H2II = (float) (H2IILuminosity * H2IISigma / (4.0 * M_PI));
+ kdiss_r2 = (float) (LWLuminosity * H2ISigma / (4.0 * pi));
+ kph_hm = (float) (IRLuminosity * HMSigma / (4.0 * pi));
+ kdiss_H2II = (float) (H2IILuminosity * H2IISigma / (4.0 * pi));
for (k = 0; k < ActiveDims[2]; k++) {
for (j = 0; j < ActiveDims[1]; j++) {
radius2_yz = ddr2[1][j] + ddr2[2][k];
@@ -291,7 +292,8 @@
double radius2, radius2_yz;
- kdiss_r2 = (float) (LWLuminosity * H2ISigma / (4.0 * M_PI));
+ kdiss_r2 = (float) (LWLuminosity * H2ISigma / (4.0 * pi));
+
for (k = 0; k < ActiveDims[2]; k++) {
for (j = 0; j < ActiveDims[1]; j++) {
radius2_yz = ddr2[1][j] + ddr2[2][k];
@@ -365,9 +367,9 @@
int TemperatureField = 0;
/* Pre-compute some quantities to speed things up */
TemperatureField = this->GetTemperatureFieldNumberForH2Shield();
- kdiss_r2 = (float) (LWLuminosity * H2ISigma / (4.0 * M_PI));
- kph_hm = (float) (IRLuminosity * HMSigma / (4.0 * M_PI));
- kdiss_H2II = (float) (H2IILuminosity * H2IISigma / (4.0 * M_PI));
+ kdiss_r2 = (float) (LWLuminosity * H2ISigma / (4.0 * pi));
+ kph_hm = (float) (IRLuminosity * HMSigma / (4.0 * pi));
+ kdiss_H2II = (float) (H2IILuminosity * H2IISigma / (4.0 * pi));
for (k = 0; k < ActiveDims[2]; k++) {
for (j = 0; j < ActiveDims[1]; j++) {
radius2_yz = ddr2[1][j] + ddr2[2][k];
@@ -455,7 +457,7 @@
{
float jeans_length = 0.0;
- jeans_length = 15*kboltz*T/(4.0*M_PI*GravConst*mh*dens*density_units);
+ jeans_length = 15*kboltz*T/(4.0*pi*GravConst*mh*dens*density_units);
return sqrt(jeans_length);
}
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/Grid_AddH2DissociationFromTree.C
--- a/src/enzo/Grid_AddH2DissociationFromTree.C
+++ b/src/enzo/Grid_AddH2DissociationFromTree.C
@@ -23,6 +23,7 @@
#include "Hierarchy.h"
#include "CosmologyParameters.h"
#include "Star.h"
+#include "phys_constants.h"
#define MIN_OPENING_ANGLE 0.2 // 0.2 = arctan(11.3 deg)
@@ -38,8 +39,6 @@
int grid::AddH2DissociationFromTree(void)
{
- const double pc = 3.086e18, clight = 3e10;
-
int i, j, k, index, dim, ci;
FLOAT pos[MAX_DIMENSION];
FLOAT radius2;
@@ -77,7 +76,7 @@
H2ISigma *= (double)TimeUnits / ((double)LengthUnits * (double)LengthUnits);
// Dilution factor (prevent breaking in the rate solver near the star)
- float dilutionRadius = 10.0 * pc / (double) LengthUnits;
+ float dilutionRadius = 10.0 * pc_cm / (double) LengthUnits;
float dilRadius2 = dilutionRadius * dilutionRadius;
// Convert from #/s to RT units
@@ -87,7 +86,7 @@
/* Find sources in the tree that contribute to the cells */
SuperSourceEntry *Leaf;
- float factor = LConv_inv * H2ISigma / (4.0 * M_PI);
+ float factor = LConv_inv * H2ISigma / (4.0 * pi);
float angle;
Leaf = SourceClusteringTree;
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/Grid_AddMagneticSupernovaeToList.C
--- a/src/enzo/Grid_AddMagneticSupernovaeToList.C
+++ b/src/enzo/Grid_AddMagneticSupernovaeToList.C
@@ -56,22 +56,22 @@
// below based on the resolution of the grid at the highest refinement level
if (UseMagneticSupernovaFeedback > 1){
sn_duration = 5.0 * dtFixed;
- MagneticSupernovaDuration = sn_duration * TimeUnits / SecondsPerYear;
+ MagneticSupernovaDuration = sn_duration * TimeUnits / yr_s;
sn_radius = 1.5 * this->CellWidth[0][0];
- MagneticSupernovaRadius = sn_radius * LengthUnits / pc;
+ MagneticSupernovaRadius = sn_radius * LengthUnits / pc_cm;
}
else {
// Converting time from years to seconds, then internal units
- sn_duration = MagneticSupernovaDuration * SecondsPerYear / TimeUnits;
+ sn_duration = MagneticSupernovaDuration * yr_s / TimeUnits;
if(sn_duration < 5.0 * dtFixed) {
printf("WARNING: Magnetic supernova feedback duration is less than the recommended minimum of 5 x dtFixed\n");
- printf("Current dtFixed = %e years \n", dtFixed * TimeUnits / SecondsPerYear);
+ printf("Current dtFixed = %e years \n", dtFixed * TimeUnits / yr_s);
}
// Converting radius from parsecs to cm, then internal units
- sn_radius = MagneticSupernovaRadius * pc / LengthUnits;
+ sn_radius = MagneticSupernovaRadius * pc_cm / LengthUnits;
if(sn_radius < 1.5 * this->CellWidth[0][0]){
printf("WARNING: Magnetic supernova feedback radius is less than the recommended minimum of 1.5 x CellWidth\n");
- printf("Current CellWidth = %e pc \n", this->CellWidth[0][0] * LengthUnits / pc);
+ printf("Current CellWidth = %e pc \n", this->CellWidth[0][0] * LengthUnits / pc_cm);
}
}
@@ -93,7 +93,7 @@
random_u = (float)(mt_random()%32768)/32768.0; // random variable from 0 to 1
random_v = (float)(mt_random()%32768)/32768.0;
- random_phi = 2*M_PI*random_u; // 0 to 2pi
+ random_phi = 2*pi*random_u; // 0 to 2pi
random_theta = acos(2*random_v-1); // 0 to pi
// Setting up randomly oriented magnetic feedback of supernova
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/Grid_AddRadiationImpulse.C
--- a/src/enzo/Grid_AddRadiationImpulse.C
+++ b/src/enzo/Grid_AddRadiationImpulse.C
@@ -23,6 +23,7 @@
#include "Hierarchy.h"
#include "CosmologyParameters.h"
#include "Star.h"
+#include "phys_constants.h"
int FindField(int f, int farray[], int n);
int GetUnits(float *DensityUnits, float *LengthUnits,
@@ -33,8 +34,6 @@
FLOAT BirthTime, FLOAT* pos)
{
- const double pc = 3.086e19, clight = 3e10;
-
int i, j, k, index, dim, FieldNum;
float kdiss_r2;
FLOAT delx, dely, delz, DomainWidth[MAX_DIMENSION];
@@ -62,7 +61,7 @@
DomainWidth[dim] = DomainRightEdge[dim] - DomainLeftEdge[dim];
// Dilution factor (prevent breaking in the rate solver near the star)
- float dilutionRadius = 10.0 * pc / (double) LengthUnits;
+ float dilutionRadius = 10.0 * pc_cm / (double) LengthUnits;
float dilRadius2 = dilutionRadius * dilutionRadius;
float LightTravelDist = TimeUnits * clight / LengthUnits;
@@ -73,7 +72,7 @@
/* Loop over cells */
index = 0;
- kdiss_r2 = (float) (Luminosity * sigma / (4.0 * M_PI));
+ kdiss_r2 = (float) (Luminosity * sigma / (4.0 * pi));
for (k = 0; k < GridDimension[2]; k++) {
delz = fabs(CellLeftEdge[2][k] + 0.5*CellWidth[2][k] - pos[2]);
delz = min(delz, DomainWidth[2]-delz);
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/Grid_CRShockTubesInitializeGrid.C
--- a/src/enzo/Grid_CRShockTubesInitializeGrid.C
+++ b/src/enzo/Grid_CRShockTubesInitializeGrid.C
@@ -10,6 +10,7 @@
#include "ExternalBoundary.h"
#include "Grid.h"
#include "CosmologyParameters.h"
+#include "phys_constants.h"
int GetUnits(float *DensityUnits, float *LengthUnits,
float *TemperatureUnits, float *TimeUnits,
@@ -216,7 +217,7 @@
for (i = 0; i < GridDimension[0]; i++) {
x = CellLeftEdge[0][i] + 0.5*CellWidth[0][i];
x = x - x0; // translate gaussian to center
- BaryonField[iCRD][i] = 1.0/sqrt(4.0*3.14159*CRkappa*t0)
+ BaryonField[iCRD][i] = 1.0/sqrt(4.0*pi*CRkappa*t0)
* PEXP( -x*x/(4.0*CRkappa*t0));
} // end i for
} // end if
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/Grid_ChangeParticleTypeBeforeSN.C
--- a/src/enzo/Grid_ChangeParticleTypeBeforeSN.C
+++ b/src/enzo/Grid_ChangeParticleTypeBeforeSN.C
@@ -23,6 +23,7 @@
#include "GridList.h"
#include "ExternalBoundary.h"
#include "Grid.h"
+#include "phys_constants.h"
int GetUnits(float *DensityUnits, float *LengthUnits,
float *TemperatureUnits, float *TimeUnits,
@@ -40,7 +41,6 @@
if (Stars == NULL)
return SUCCESS;
- const float pc = 3.086e18, mh = 1.673e-24, Msun = 1.989e33;
const float PISNLowerMass = 140.0, PISNUpperMass = 260.0;
const float StartRefineAtTime = 0.99; // Percentage of stellar lifetime
const float EndRefineAtTime = 1.0;
@@ -98,7 +98,7 @@
if (factor < 1) {
Diameter = 2 * BufferZone * PopIIISupernovaRadius *
- (pc/LengthUnits);
+ (pc_cm/LengthUnits);
}
@@ -166,9 +166,6 @@
/* ASSUMES POP III PAIR-INSTABILITY SUPERNOVA FOR ENERGY */
- const float pc = 3.086e18, mh = 1.673e-24, Msun = 1.989e33;
- const float kb = 1.38e-16;
-
FLOAT StartTime, SoundCrossingTime;
float HeliumCoreMass, TransitionTime, ShockVelocity, STradius;
float BlastWaveRadius, SoundSpeed, SNTemperature;
@@ -183,14 +180,14 @@
// With a constant shock velocity, calculate time it takes for
// blastwave to travel from r=0 to r=PopIIISupernovaRadius
- ShockVelocity = 1.165 * sqrt(2.0 * SNenergy / (Mass * Msun));
- StartTime = PopIIISupernovaRadius * pc / ShockVelocity;
+ ShockVelocity = 1.165 * sqrt(2.0 * SNenergy / (Mass * SolarMass));
+ StartTime = PopIIISupernovaRadius * pc_cm / ShockVelocity;
// Because we inject thermal energy, the blastwave is delayed by a
// sound crossing time.
- SNTemperature = min(double(Mass*Msun) / double(mh) / kb, 1e8);
- SoundSpeed = sqrt(kb * SNTemperature / (0.6*mh));
- SoundCrossingTime = PopIIISupernovaRadius * pc / SoundSpeed;
+ SNTemperature = min(double(Mass*SolarMass) / double(mh) / kboltz, 1e8);
+ SoundSpeed = sqrt(kboltz * SNTemperature / (0.6*mh));
+ SoundCrossingTime = PopIIISupernovaRadius * pc_cm / SoundSpeed;
// units in cm
STradius = 3.62e19 * powf(Mass/100., 1./3.) * powf(n0, -1./3.);
@@ -212,7 +209,7 @@
// it slows by a factor of 2 near the transition to the ST phase)
if (Time < TransitionTime) {
- BlastWaveRadius = PopIIISupernovaRadius * pc + ShockVelocity * Time;
+ BlastWaveRadius = PopIIISupernovaRadius * pc_cm + ShockVelocity * Time;
printf("Free expansion\n");
} // ENDIF Free expansion phase
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/Grid_ClusterInitializeGrid.C
--- a/src/enzo/Grid_ClusterInitializeGrid.C
+++ b/src/enzo/Grid_ClusterInitializeGrid.C
@@ -21,6 +21,7 @@
#include "ExternalBoundary.h"
#include "Grid.h"
#include "CosmologyParameters.h"
+#include "phys_constants.h"
int GetUnits(float *DensityUnits, float *LengthUnits,
float *TemperatureUnits, float *TimeUnits,
@@ -95,8 +96,8 @@
/* Set various units. */
- const double Mpc = 3.0856e24, SolarMass = 1.989e33, GravConst = 6.67e-8,
- pi = 3.14159, mh = 1.67e-24, kboltz = 1.381e-16, keV=1.1604e7;
+ const double keV=1.1604e7;
+
float DensityUnits, LengthUnits, TemperatureUnits = 1, TimeUnits,
VelocityUnits, CriticalDensity = 1, BoxLength = 1, mu = 0.6;
double MassUnits = 1;
@@ -110,7 +111,7 @@
BoxLength = ComovingBoxSize*ExpansionFactor/HubbleConstantNow; // in Mpc
} else {
CriticalDensity = 2.78e11*POW(0.74,2); // in Msolar/Mpc^3 for h=0.74
- BoxLength = LengthUnits / 3.086e24;
+ BoxLength = LengthUnits / Mpc_cm;
HubbleConstantNow = 1.0;
OmegaMatterNow = 1.0;
}
@@ -161,32 +162,32 @@
x1 = NFWRadius[i]/SphereCoreRadius[sphere];
NFWDensity[i] = SphereDensity[sphere]/(x1*(1.0+x1)*(1.0+x1)); // DM Density
if (SphereType[sphere]>=6 && SphereType[sphere] <= 8) { //aka 6, 7, 8: Perseus Cluster
- rkpc=NFWRadius[i]*LengthUnits/(1.0e-3*Mpc);
+ rkpc=NFWRadius[i]*LengthUnits/(1.0e-3*Mpc_cm);
/* Initial Temperature */
if (rkpc > 300.0){
NFWTemp[i]=7.0594*1.3*POW((1.0+1.5*NFWRadius[i]/SphereRadius[sphere]),-1.6)*keV; // in K
} else{
- NFWTemp[i]=8.12e7*(1.0+POW(NFWRadius[i]*LengthUnits/(1.0e-3*Mpc*71),3))/(2.3 + pow(NFWRadius[i]*LengthUnits/(1.0e-3*Mpc*71),3));
+ NFWTemp[i]=8.12e7*(1.0+POW(NFWRadius[i]*LengthUnits/(1.0e-3*Mpc_cm*71),3))/(2.3 + pow(NFWRadius[i]*LengthUnits/(1.0e-3*Mpc_cm*71),3));
if (SphereType[sphere]==8)
NFWTemp[i]=8.12e7;
}
/* Set Gravity. NFW Dark Matter, BCG+BH */
Allg[i]=GravConst*PointSourceGravityConstant*SolarMass*
((log(1.0+x1)-x1/(1.0+x1)) /(log(1.0+1.0)-1.0/(1.0+1.0)))/POW(NFWRadius[i]*LengthUnits, 2.0) +
- ClusterSMBHBCG*POW((POW(POW(NFWRadius[i]*LengthUnits/(1.0e-3*Mpc),0.5975)/3.206e-7,0.9) +
- POW(POW(NFWRadius[i]*LengthUnits/(1.0e-3*Mpc), 1.849)/1.861e-6, 0.9)), -1.0/0.9) +
+ ClusterSMBHBCG*POW((POW(POW(NFWRadius[i]*LengthUnits/(1.0e-3*Mpc_cm),0.5975)/3.206e-7,0.9) +
+ POW(POW(NFWRadius[i]*LengthUnits/(1.0e-3*Mpc_cm), 1.849)/1.861e-6, 0.9)), -1.0/0.9) +
GravConst*SolarMass* ClusterSMBHMass / POW(NFWRadius[i]*LengthUnits, 2) ;
}//end Perseus
else{ //Elliptical galaxies
if (SphereType[sphere]==1) { //NGC 4472
- NFWTemp[i]=1.16059e7*(1.17-(1.17-0.6)*exp(-NFWRadius[i]*LengthUnits/(2.0*5*1.0e-3*Mpc)));
+ NFWTemp[i]=1.16059e7*(1.17-(1.17-0.6)*exp(-NFWRadius[i]*LengthUnits/(2.0*5*1.0e-3*Mpc_cm)));
}
if (SphereType[sphere]==2) { //NGC 6482
- NFWTemp[i]=1.16059e7*(0.4+0.4*exp(-NFWRadius[i]*LengthUnits/(2.0*8.0*1.0e-3*Mpc)));
+ NFWTemp[i]=1.16059e7*(0.4+0.4*exp(-NFWRadius[i]*LengthUnits/(2.0*8.0*1.0e-3*Mpc_cm)));
}
Allg[i]=GravConst*PointSourceGravityConstant*SolarMass*
((log(1.0+x1)-x1/(1.0+x1)) /(log(1.0+1.0)-1.0/(1.0+1.0)))/POW(NFWRadius[i]*LengthUnits, 2.0) +
- GravConst*(ClusterSMBHBCG*SolarMass*1.0e11)/POW(NFWRadius[i]*LengthUnits+EllipticalGalaxyRe*1.0e-3*Mpc/1.8153, 2) + //ClusterSMBHBCG is M_* here
+ GravConst*(ClusterSMBHBCG*SolarMass*1.0e11)/POW(NFWRadius[i]*LengthUnits+EllipticalGalaxyRe*1.0e-3*Mpc_cm/1.8153, 2) + //ClusterSMBHBCG is M_* here
GravConst*SolarMass* ClusterSMBHMass / POW(NFWRadius[i]*LengthUnits, 2) ;
}
@@ -197,7 +198,7 @@
}
dpdr_old=dpdr;
if (i>0) {
- NFWPressure[i]= (NFWPressure[i-1] + (NFWRadius[i]-NFWRadius[i-1])*16.0*Mpc*0.5*dpdr_old)/(1.0+(NFWRadius[i]-NFWRadius[i-1])*16.0*Mpc*0.5*Allg[i]/(kboltz*NFWTemp[i]/(mu * mh)));
+ NFWPressure[i]= (NFWPressure[i-1] + (NFWRadius[i]-NFWRadius[i-1])*16.0*Mpc_cm*0.5*dpdr_old)/(1.0+(NFWRadius[i]-NFWRadius[i-1])*16.0*Mpc_cm*0.5*Allg[i]/(kboltz*NFWTemp[i]/(mu * mh)));
GasDensity[i]=NFWPressure[i]/(kboltz * NFWTemp[i]/(mu * mh));
dpdr = -Allg[i]* GasDensity[i];
}
@@ -275,8 +276,8 @@
temp1 = NFWTemp[m] + (NFWTemp[m-1] - NFWTemp[m])*
(r - NFWRadius[m])/(NFWRadius[m-1] - NFWRadius[m]); // in K
if (SphereType[sphere] == 6) {
- gas_dens1 =mh*(0.0192/(1.0+POW(r*LengthUnits/(18.0e-3*Mpc),3.0))+0.046/pow((1.0+pow(r*LengthUnits/(57.0e-3*Mpc), 2.0)), 1.8)+
- 0.0048/POW((1.0+pow(r*LengthUnits/(200.0e-3*Mpc), 2.0)), 1.1))/DensityUnits/0.88;
+ gas_dens1 =mh*(0.0192/(1.0+POW(r*LengthUnits/(18.0e-3*Mpc_cm),3.0))+0.046/pow((1.0+pow(r*LengthUnits/(57.0e-3*Mpc_cm), 2.0)), 1.8)+
+ 0.0048/POW((1.0+pow(r*LengthUnits/(200.0e-3*Mpc_cm), 2.0)), 1.1))/DensityUnits/0.88;
}
else{
gas_dens1 = (GasDensity[m] + (GasDensity[m-1] - GasDensity[m])*
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/Grid_ClusterSMBHEachGridGasMass.C
--- a/src/enzo/Grid_ClusterSMBHEachGridGasMass.C
+++ b/src/enzo/Grid_ClusterSMBHEachGridGasMass.C
@@ -65,7 +65,7 @@
int dim = 0;
float DiskRadius; //ClusterSMBHDiskRadius = 0.5; ClusterSMBHDiskRadiu now a parameter.
- DiskRadius = ClusterSMBHDiskRadius*kpc/LengthUnits; //from kpc to codeunits
+ DiskRadius = ClusterSMBHDiskRadius*kpc_cm/LengthUnits; //from kpc to codeunits
for (dim = 0; dim < GridRank; dim++) {
DiskCenter[dim] = PointSourceGravityPosition[dim];
DiskLeftCorner[dim] = PointSourceGravityPosition[dim]- DiskRadius;
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/Grid_ClusterSMBHFeedback.C
--- a/src/enzo/Grid_ClusterSMBHFeedback.C
+++ b/src/enzo/Grid_ClusterSMBHFeedback.C
@@ -85,7 +85,7 @@
float JetVelocity, FastJetVelocity; // Jet Velocity in km/s (should make parameter)-- gets value from parameter ClusterSMBHJetVelocity
JetScaleRadius = ClusterSMBHJetRadius/2.0; //JetScaleRadius is half the radius of the jet launch region in cellwidths
float DiskRadius; //ClusterSMBHDiskRadius = 0.5; //ClusterSMBHDiskRadiu now a parameter
- DiskRadius = ClusterSMBHDiskRadius*kpc/LengthUnits; //from kpc to codeunits
+ DiskRadius = ClusterSMBHDiskRadius*kpc_cm/LengthUnits; //from kpc to codeunits
for (dim = 0; dim < GridRank; dim++) {
JetCenter[dim] = PointSourceGravityPosition[dim];
@@ -100,8 +100,8 @@
DiskRightCorner[dim] = PointSourceGravityPosition[dim] + DiskRadius;
}
- JetLeftCorner[jet_dim] -= ClusterSMBHJetLaunchOffset*kpc/LengthUnits;
- JetRightCorner[jet_dim] += ClusterSMBHJetLaunchOffset*kpc/LengthUnits; //from kpc to codeunits
+ JetLeftCorner[jet_dim] -= ClusterSMBHJetLaunchOffset*kpc_cm/LengthUnits;
+ JetRightCorner[jet_dim] += ClusterSMBHJetLaunchOffset*kpc_cm/LengthUnits; //from kpc to codeunits
/* Compute indices of jet launch region. */
@@ -143,7 +143,7 @@
} // end: loop over dim
/* Compute mass and momentum to be put into cells in code units if Jet is on this grid. */
- JetMdot = (ClusterSMBHJetMdot*SolarMass/3.1557e7)/(MassUnits/TimeUnits); // from M_sun/yr to code units
+ JetMdot = (ClusterSMBHJetMdot*SolarMass/yr_s)/(MassUnits/TimeUnits); // from M_sun/yr to code units
if (JetOnGrid == true){
float JetNormalization = 0.0, density_normalization, radius, Tramp;
if (jet_dim == 2){
@@ -179,9 +179,9 @@
/* Convert to code units. */
density_normalization = (JetMdot/JetNormalization)*dtFixed/POW(CellWidth[0][0], 3);
- Tramp = ClusterSMBHTramp*1.0e6*3.1557e7/TimeUnits; // from Myr to code units
+ Tramp = ClusterSMBHTramp*1.0e6*yr_s/TimeUnits; // from Myr to code units
- JetVelocity = sqrt((ClusterSMBHJetEdot*1.0e44*ClusterSMBHKineticFraction*2)/(ClusterSMBHJetMdot*SolarMass/3.1557e7))/VelocityUnits;
+ JetVelocity = sqrt((ClusterSMBHJetEdot*1.0e44*ClusterSMBHKineticFraction*2)/(ClusterSMBHJetMdot*SolarMass/yr_s))/VelocityUnits;
JetVelocity *= min((Time-ClusterSMBHStartTime)/Tramp, 1.0); //linear ramp
/* Clip edge of jet launching disk so we don't set cell off the edge of the grid. */
@@ -212,7 +212,7 @@
JetVelocity_x = JetVelocity; // mutiplied by sincos later
JetVelocity_y = JetVelocity; // mutiplied by sincos later
if (ClusterSMBHJetPrecessionPeriod > 0.00001 || ClusterSMBHJetPrecessionPeriod < -0.00001)
- ClusterSMBHJetAnglePhi = Time*2.0/(ClusterSMBHJetPrecessionPeriod*1.0e6*3.1557e7/TimeUnits); // ClusterSMBHJetPrecessionPeriod from Myr to codeunit; *2.0 instead of 2*pi because pi is used later
+ ClusterSMBHJetAnglePhi = Time*2.0/(ClusterSMBHJetPrecessionPeriod*1.0e6*yr_s/TimeUnits); // ClusterSMBHJetPrecessionPeriod from Myr to codeunit; *2.0 instead of 2*pi because pi is used later
}
else {
JetVelocity_z = JetVelocity * ClusterSMBHJetOpenAngleRadius / sqrt(POW(ClusterSMBHJetOpenAngleRadius, 2) + POW(radius, 2));
@@ -273,7 +273,7 @@
JetVelocity_z = JetVelocity; // mutiplied by sincos later
JetVelocity_y = JetVelocity; // mutiplied by sincos later
if (ClusterSMBHJetPrecessionPeriod > 0.00001)
- ClusterSMBHJetAnglePhi = Time*2.0/(ClusterSMBHJetPrecessionPeriod*1.0e6*3.1557e7/TimeUnits); // ClusterSMBHJetPrecessionPeriod from Myr to codeunit; *2.0 instead of 2*pi because pi is used later
+ ClusterSMBHJetAnglePhi = Time*2.0/(ClusterSMBHJetPrecessionPeriod*1.0e6*yr_s/TimeUnits); // ClusterSMBHJetPrecessionPeriod from Myr to codeunit; *2.0 instead of 2*pi because pi is used later
}
else {
JetVelocity_x = JetVelocity * ClusterSMBHJetOpenAngleRadius / sqrt(POW(ClusterSMBHJetOpenAngleRadius, 2) + POW(radius, 2));
@@ -335,7 +335,7 @@
JetVelocity_z = JetVelocity; // mutiplied by sincos later
JetVelocity_x = JetVelocity; // mutiplied by sincos later
if (ClusterSMBHJetPrecessionPeriod > 0.00001)
- ClusterSMBHJetAnglePhi = Time*2.0/(ClusterSMBHJetPrecessionPeriod*1.0e6*3.1557e7/TimeUnits); // ClusterSMBHJetPrecessionPeriod from Myr to codeunit; *2.0 instead of 2*pi because pi is used later
+ ClusterSMBHJetAnglePhi = Time*2.0/(ClusterSMBHJetPrecessionPeriod*1.0e6*yr_s/TimeUnits); // ClusterSMBHJetPrecessionPeriod from Myr to codeunit; *2.0 instead of 2*pi because pi is used later
}
else {
JetVelocity_y = JetVelocity * ClusterSMBHJetOpenAngleRadius / sqrt(POW(ClusterSMBHJetOpenAngleRadius, 2) + POW(radius, 2));
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/Grid_CollapseTestInitializeGrid.C
--- a/src/enzo/Grid_CollapseTestInitializeGrid.C
+++ b/src/enzo/Grid_CollapseTestInitializeGrid.C
@@ -24,6 +24,7 @@
#include "ExternalBoundary.h"
#include "Grid.h"
#include "CosmologyParameters.h"
+#include "phys_constants.h"
#define NTHETA 1000
#define NR 1000
@@ -157,8 +158,6 @@
/* Set various units. */
- const double Mpc = 3.0856e24, SolarMass = 1.989e33, GravConst = 6.67e-8,
- pi = 3.14159, mh = 1.67e-24, kboltz = 1.381e-16, LightSpeed = 2.9979e10;
float DensityUnits, LengthUnits, TemperatureUnits, TimeUnits,
VelocityUnits, CriticalDensity = 1, BoxLength = 1, mu = 0.6;
@@ -168,11 +167,11 @@
if (ComovingCoordinates) {
CosmologyComputeExpansionFactor(Time, &a, &dadt);
ExpansionFactor = a/(1.0+InitialRedshift);
- CriticalDensity = 2.78e11*pow(HubbleConstantNow, 2); // in Msolar/Mpc^3
+ CriticalDensity = 2.78e11*POW(HubbleConstantNow, 2); // in Msolar/Mpc^3
BoxLength = ComovingBoxSize*ExpansionFactor/HubbleConstantNow; // in Mpc
} else {
- CriticalDensity = 2.78e11*pow(0.74,2); // in Msolar/Mpc^3 for h=0.74
- BoxLength = LengthUnits / 3.086e24;
+ CriticalDensity = 2.78e11*POW(0.74,2); // in Msolar/Mpc^3 for h=0.74
+ BoxLength = LengthUnits / Mpc_cm;
HubbleConstantNow = 1.0;
OmegaMatterNow = 1.0;
}
@@ -193,20 +192,20 @@
NFWPressure[0] = 1.0 * kboltz * InitialTemperature / (mu * mh);
FILE *fptr = fopen("NFWProfile.out", "w");
for (i = 0; i < NFW_POINTS; i++) {
- NFWRadius[i] = SphereRadius[sphere]*pow(10, -3*(float(i)/NFW_POINTS));
+ NFWRadius[i] = SphereRadius[sphere]*POW(10, -3*(float(i)/NFW_POINTS));
x1 = NFWRadius[i]/SphereCoreRadius[sphere];
NFWDensity[i] = SphereDensity[sphere]/(x1*(1.0+x1)*(1.0+x1));
NFWMass[i] = 4.0*pi*SphereDensity[sphere]*
- (CriticalDensity/pow(ExpansionFactor, 3)) *
- pow(SphereCoreRadius[sphere]*BoxLength, 3) *
+ (CriticalDensity/POW(ExpansionFactor, 3)) *
+ POW(SphereCoreRadius[sphere]*BoxLength, 3) *
(log(1.0+x1) - x1/(x1+1.0)); // in Msolar
dpdr_old = dpdr;
dpdr = GravConst * NFWMass[i] * SolarMass *
NFWDensity[i] /
- pow(NFWRadius[i]*BoxLength*Mpc, 2);
+ POW(NFWRadius[i]*BoxLength*Mpc_cm, 2);
if (i > 0)
NFWPressure[i] = NFWPressure[i-1] -
- 0.5*(dpdr+dpdr_old)*(NFWRadius[i]-NFWRadius[i-1])*BoxLength*Mpc;
+ 0.5*(dpdr+dpdr_old)*(NFWRadius[i]-NFWRadius[i-1])*BoxLength*Mpc_cm;
NFWTemp[i] = NFWPressure[i]*mu*mh/(kboltz*NFWDensity[i]); // in K
NFWSigma[i] = sqrt(kboltz * NFWTemp[i] / (mu * mh)); // in cm/s
float mean_overdensity = 3.0*SphereDensity[sphere] / (x1*x1*x1) *
@@ -318,7 +317,7 @@
switch (SphereType[sphere]) {
case 1:
- SphereMass = (4*pi/3)*pow((SphereRadius[sphere]*LengthUnits), 3) *
+ SphereMass = (4*pi/3.0)*POW((SphereRadius[sphere]*LengthUnits), 3) *
(SphereDensity[sphere]*DensityUnits);
printf("mass = %"GSYM", lunit = %"GSYM", dunit = %"GSYM", rho = %"GSYM", r = %"GSYM"\n",
SphereMass, LengthUnits, DensityUnits, SphereDensity[sphere],
@@ -332,32 +331,32 @@
case 4:
// Integral of a Gaussian
term1 = 0.25 * sqrt(pi) * 0.3536 * // 2^(-3/2) = 0.3536
- pow((SphereCoreRadius[sphere]), 3) *
+ POW((SphereCoreRadius[sphere]), 3) *
ERF(SphereRadius[sphere]/SphereCoreRadius[sphere]);
term2 = 0.25 * SphereRadius[sphere] *
- pow(SphereCoreRadius[sphere], 2) *
- PEXP(-0.5 * pow((SphereRadius[sphere]/SphereCoreRadius[sphere]), 2));
+ POW(SphereCoreRadius[sphere], 2) *
+ PEXP(-0.5 * POW((SphereRadius[sphere]/SphereCoreRadius[sphere]), 2));
SphereMass = (4*pi*SphereDensity[sphere]*DensityUnits) *
- pow(LengthUnits, 3) * (term1 - term2);
+ POW(LengthUnits, 3) * (term1 - term2);
break;
case 5:
SphereCoreDens = (SphereDensity[sphere]*DensityUnits) *
- pow(SphereCoreRadius[sphere] / SphereRadius[sphere], -2);
- SphereCoreMass = (4*pi/3) *
- pow((SphereCoreRadius[sphere]*LengthUnits),3) * (SphereCoreDens);
+ POW(SphereCoreRadius[sphere] / SphereRadius[sphere], -2);
+ SphereCoreMass = (4*pi/3.0) *
+ POW((SphereCoreRadius[sphere]*LengthUnits),3) * (SphereCoreDens);
SphereMass = (4*pi*SphereDensity[sphere]*DensityUnits *
- pow(SphereRadius[sphere]*LengthUnits, 2) *
+ POW(SphereRadius[sphere]*LengthUnits, 2) *
((SphereRadius[sphere] - SphereCoreRadius[sphere]) *
LengthUnits)) + SphereCoreMass;
break;
case 6:
- SphereMass = pow(VelocitySound[sphere],3) /
- (sqrt(4*pi*pow((GravConst)/(4*pi),3) *
+ SphereMass = POW(VelocitySound[sphere],3) /
+ (sqrt(4*pi*POW((GravConst)/(4*pi),3) *
(SphereDensity[sphere]))) *
q(SphereCutOff[sphere]);
- SphereMass = SphereMass*DensityUnits*pow(LengthUnits,3);
+ SphereMass = SphereMass*DensityUnits*POW(LengthUnits,3);
break;
case 7:
@@ -376,7 +375,7 @@
break;
case 8:
- SphereMass = (4*pi/3)*pow((SphereRadius[sphere]*LengthUnits), 3) *
+ SphereMass = (4*pi/3.0)*POW((SphereRadius[sphere]*LengthUnits), 3) *
(SphereDensity[sphere]*DensityUnits);
ComputeRadialVelocity(SphereDensity[sphere], SphereMass,
SphereRadius[sphere], VelocityUnits,
@@ -402,15 +401,15 @@
BHMass = SphereMass / SolarMass; // in solar masses
SchwarzschildRadius = 2.0 * GravConst * SphereMass /
- (LightSpeed*LightSpeed);
- CavityRadius = 117.0 * SchwarzschildRadius * pow((BHMass/1e6), 0.08);
- InnerDensity = 4.31e-10 * pow((BHMass/1e8), -0.8) *
- pow((CavityRadius/SchwarzschildRadius/1e3), -0.6) / DensityUnits;
- InnerTemperature = 1.7e6 * pow((BHMass/1e6), -0.28);
- InnerScaleHeight = 0.46 * CavityRadius * pow(BHMass/1e6, -0.12) /
+ (clight*clight);
+ CavityRadius = 117.0 * SchwarzschildRadius * POW((BHMass/1e6), 0.08);
+ InnerDensity = 4.31e-10 * POW((BHMass/1e8), -0.8) *
+ POW((CavityRadius/SchwarzschildRadius/1e3), -0.6) / DensityUnits;
+ InnerTemperature = 1.7e6 * POW((BHMass/1e6), -0.28);
+ InnerScaleHeight = 0.46 * CavityRadius * POW(BHMass/1e6, -0.12) /
LengthUnits;
ThickenTransitionRadius = 1.9e3 * SchwarzschildRadius *
- pow((BHMass/1e6), 2.0/21) / LengthUnits;
+ POW((BHMass/1e6), 2.0/21) / LengthUnits;
CavityRadius /= LengthUnits;
printf("cgs: %"GSYM" %"GSYM" %"GSYM" %"GSYM" %"GSYM" %"GSYM"\n", SchwarzschildRadius,
CavityRadius*LengthUnits, InnerDensity*DensityUnits,
@@ -508,9 +507,9 @@
else if (xpos > 0 && ypos < 0)
theta = 2*pi + atan(ypos/xpos);
} else if (xpos == 0 && ypos > 0)
- theta = 3.14159 / 2.0;
+ theta = pi / 2.0;
else if (xpos == 0 && ypos < 0)
- theta = (3*3.14159) / 2.0;
+ theta = (3*pi) / 2.0;
else
theta = 0.0;
@@ -529,7 +528,7 @@
RotVelocity[2] = 2*pi*(xpos*sin(a)) /
SphereRotationalPeriod[sphere];
if (SphereType[sphere] == 7 || SphereType[sphere] == 9)
- GasRotVelocityCorrection = pow(rcyl/SphereRadius[sphere],-1.5);
+ GasRotVelocityCorrection = POW(rcyl/SphereRadius[sphere],-1.5);
for (dim = 0; dim < MAX_DIMENSION; dim++) {
Velocity[dim] = GasRotVelocityCorrection * RotVelocity[dim];
DMVelocity[dim] = DMRotVelocityCorrection * RotVelocity[dim];
@@ -566,10 +565,10 @@
if (SphereType[sphere] == 1)
dens1 = SphereDensity[sphere];
- /* 2) r^-2 power law */
+ /* 2) r^-2 POWer law */
if (SphereType[sphere] == 2)
- dens1 = SphereDensity[sphere]*pow(r/SphereRadius[sphere], -2);
+ dens1 = SphereDensity[sphere]*POW(r/SphereRadius[sphere], -2);
/* 3) NFW profile (use look-up table for temperature and
velocity dispersion)*/
@@ -591,17 +590,17 @@
if (SphereType[sphere] == 4) {
dens1 = SphereDensity[sphere]*
- PEXP(-0.5*pow(r/SphereCoreRadius[sphere], 2));
+ PEXP(-0.5*POW(r/SphereCoreRadius[sphere], 2));
}
- /* 5) r^-2 power law with core radius */
+ /* 5) r^-2 POWer law with core radius */
if (SphereType[sphere] == 5) {
if (r < SphereCoreRadius[sphere]) {
- dens1 = SphereDensity[sphere]*pow(SphereCoreRadius[sphere]/
+ dens1 = SphereDensity[sphere]*POW(SphereCoreRadius[sphere]/
SphereRadius[sphere], -2);
} else {
- dens1 = SphereDensity[sphere]*pow(r/SphereRadius[sphere], -2);
+ dens1 = SphereDensity[sphere]*POW(r/SphereRadius[sphere], -2);
}
}
@@ -638,14 +637,14 @@
/* 9) Circumbinary BH accretion disk (Lippai et al. 2008) */
if (SphereType[sphere] == 9) {
- MidplaneDensity = InnerDensity * pow(rcyl/CavityRadius, -0.6);
+ MidplaneDensity = InnerDensity * POW(rcyl/CavityRadius, -0.6);
MidplaneTemperature =
- InnerTemperature * pow(rcyl/CavityRadius, -0.9);
+ InnerTemperature * POW(rcyl/CavityRadius, -0.9);
if (rcyl < ThickenTransitionRadius)
ScaleHeight = InnerScaleHeight;
else
ScaleHeight = InnerScaleHeight *
- pow(rcyl/ThickenTransitionRadius, 1.05);
+ POW(rcyl/ThickenTransitionRadius, 1.05);
ScaleHeight = max(ScaleHeight, 2*CellWidth[0][0]);
// printf("r=%"FSYM", z=%"FSYM", h=%"GSYM", rho=%"GSYM", T=%"GSYM"\n",
@@ -708,7 +707,7 @@
if (dim == 0)
dens1 = SphereDensity[sphere]*PEXP(-drad/ScaleHeightR)/
- pow(cosh(zheight/max(ScaleHeightz, CellWidth[0][0])), 2);
+ POW(cosh(zheight/max(ScaleHeightz, CellWidth[0][0])), 2);
// if (dens1 < density)
// break;
@@ -723,11 +722,11 @@
float accel = PointSourceGravityConstant;
if (PointSourceGravity == 1)
accel = PointSourceGravityConstant/
- (pow(drad,3) + pow(PointSourceGravityCoreRadius, 3));
+ (POW(drad,3) + POW(PointSourceGravityCoreRadius, 3));
if (PointSourceGravity == 2) {
x1 = drad/PointSourceGravityCoreRadius;
accel = PointSourceGravityConstant*(log(1+x1)-x1/(1+x1))/
- pow(drad, 3);
+ POW(drad, 3);
}
float vel = sqrt(accel);
@@ -748,10 +747,10 @@
} // end: disk
- /* 11) stellar wind, with r^-2 power law density*/
+ /* 11) stellar wind, with r^-2 POWer law density*/
if (SphereType[sphere] == 11) {
radial_velocity = StellarWindSpeed/VelocityUnits;
- dens1 = StellarWindDensity*pow(r/StellarWindRadius, -2);
+ dens1 = StellarWindDensity*POW(r/StellarWindRadius, -2);
Velocity[0] += radial_velocity * xpos / r;
Velocity[1] += radial_velocity * ypos / r;
Velocity[2] += radial_velocity * zpos / r;
@@ -902,13 +901,13 @@
if (MultiSpecies > 1) {
BaryonField[HMNum][n] = CollapseTestInitialFractionHM *
- BaryonField[HIINum][n]* pow(temperature,float(0.88));
+ BaryonField[HIINum][n]* POW(temperature,float(0.88));
BaryonField[H2IINum][n] = CollapseTestInitialFractionH2II *
- 2.0*BaryonField[HIINum][n]* pow(temperature,float(1.8));
+ 2.0*BaryonField[HIINum][n]* POW(temperature,float(1.8));
if (ComovingCoordinates)
BaryonField[H2INum][n] = H2I_Fraction *
- BaryonField[0][n]*CoolData.HydrogenFractionByMass*pow(301.0,5.1)*
- pow(OmegaMatterNow, float(1.5))/
+ BaryonField[0][n]*CoolData.HydrogenFractionByMass*POW(301.0,5.1)*
+ POW(OmegaMatterNow, float(1.5))/
(OmegaMatterNow*BaryonMeanDensity)/
// CosmologySimulationOmegaBaryonNow/
HubbleConstantNow*2.0;
@@ -968,7 +967,7 @@
if (HydroMethod != Zeus_Hydro)
for (dim = 0; dim < GridRank; dim++)
- BaryonField[1][n] += 0.5*pow(BaryonField[ivel+dim][n], 2);
+ BaryonField[1][n] += 0.5*POW(BaryonField[ivel+dim][n], 2);
/* Set particles if being used (generate a number of particle
proportional to density). */
@@ -977,11 +976,11 @@
if (i >= GridStartIndex[0] && i <= GridEndIndex[0] &&
j >= GridStartIndex[1] && j <= GridEndIndex[1] &&
k >= GridStartIndex[2] && k <= GridEndIndex[2] ) {
- ParticleCount += density/pow(float(RefineBy), GridRank*level);
+ ParticleCount += density/POW(float(RefineBy), GridRank*level);
while (ParticleCount > 1) {
if (SetupLoopCount > 0) {
ParticleMass[npart] = ParticleMeanDensity*
- pow(float(RefineBy), GridRank*level);
+ POW(float(RefineBy), GridRank*level);
ParticleNumber[npart] = CollapseTestParticleCount++;
ParticleType[npart] = PARTICLE_TYPE_DARK_MATTER;
@@ -1065,8 +1064,8 @@
double BE(double r)
{
double factor;
- factor = 4.8089e-04*pow(r,5) - 1.0173e-02*pow(r,4) + 7.7899e-02*pow(r,3) -
- 2.3299e-01*pow(r,2) + 1.4721e-02*r + 1.0008e+00;
+ factor = 4.8089e-04*POW(r,5) - 1.0173e-02*POW(r,4) + 7.7899e-02*POW(r,3) -
+ 2.3299e-01*POW(r,2) + 1.4721e-02*r + 1.0008e+00;
return factor;
}
@@ -1075,8 +1074,8 @@
double q(double r)
{
double factor;
- factor = 0.0015970*pow(r,5) - 0.0229113*pow(r,4) + 0.0386709*pow(r,3) +
- 0.7350457*pow(r,2) - 0.5490283*r + 0.0872061;
+ factor = 0.0015970*POW(r,5) - 0.0229113*POW(r,4) + 0.0386709*POW(r,3) +
+ 0.7350457*POW(r,2) - 0.5490283*r + 0.0872061;
return factor;
}
@@ -1092,13 +1091,10 @@
double Maxwellian(double c_tilda, double vel_unit, double mu, double gamma)
{
// Set constants
- double mh = 1.67e-24;
- double kboltz = 1.38e-16;
- double pi = 3.14159;
// Compute temperature in cgs units
double c = c_tilda*vel_unit;
- double T = (pow(c,2)*mu*mh)/(kboltz*gamma);
+ double T = (POW(c,2)*mu*mh)/(kboltz*gamma);
// Compute random Maxwellian velocity
double mean = 0;
@@ -1114,8 +1110,8 @@
double ERF(double x)
{
- return (2.0 / sqrt(M_PI)) *
- (x - pow(x,3)/3.0 + pow(x,5)/10.0 - pow(x,7)/42.0 + pow(x,9)/216.0);
+ return (2.0 / sqrt(pi)) *
+ (x - POW(x,3)/3.0 + POW(x,5)/10.0 - POW(x,7)/42.0 + POW(x,9)/216.0);
}
int ComputeRadialVelocity(float density, double mass, float r_init,
@@ -1124,8 +1120,7 @@
double exterior_rho[], int Npts)
{
- const float theta0 = 0.5*M_PI, theta1 = 1.9*M_PI;
- const float Grav = 6.673e-8, Mpc = 3.086e24, yr = 3.1557e7, Msun = 1.989e33;
+ const float theta0 = 0.5*pi, theta1 = 1.9*pi;
const float delta_i = 0.5;
float dtheta, delta_t[NTHETA], Theta[NTHETA], Lambda[NTHETA], Chi;
float beta[NTHETA], little_d[NTHETA];
@@ -1138,44 +1133,44 @@
dtheta = (theta1 - theta0) / (NTHETA-1.0);
for (i = 0; i < NTHETA; i++) {
Theta[i] = theta0 + i*dtheta;
- Lambda[i] = pow(sin(Theta[i]/2), 2) *
- pow((Theta[i] - sin(Theta[i])) / M_PI, -8.0/9);
- delta_t[i] = 4.5 * pow((Theta[i] - sin(Theta[i])), 2) /
- pow((1.0 - cos(Theta[i])), 3) - 1.0;
+ Lambda[i] = POW(sin(Theta[i]/2), 2) *
+ POW((Theta[i] - sin(Theta[i])) / pi, -8.0/9);
+ delta_t[i] = 4.5 * POW((Theta[i] - sin(Theta[i])), 2) /
+ POW((1.0 - cos(Theta[i])), 3) - 1.0;
beta[i] = sin(0.5*Theta[i]) * sin(0.5*Theta[i]);
little_d[i] = 0.75 * (Theta[i] - sin(Theta[i]));
if (delta_t[i] < density)
iinit = i;
- if (Theta[i] < M_PI) // Turnaround
+ if (Theta[i] < pi) // Turnaround
iturn = i;
}
- t_init = 5.38e8 * yr * pow((1+InitialRedshift) / 10.0, -1.5);
- z_vir = pow((Theta[iinit] - sin(Theta[iinit])) / (2*M_PI), 2.0/3) *
+ t_init = 5.38e8 * yr_s * POW((1+InitialRedshift) / 10.0, -1.5);
+ z_vir = POW((Theta[iinit] - sin(Theta[iinit])) / (2*pi), 2.0/3.0) *
(1 + InitialRedshift) - 1.0;
- t_vir = 5.38e8 * yr * pow((1+z_vir) / 10.0, -1.5);
+ t_vir = 5.38e8 * yr_s * POW((1+z_vir) / 10.0, -1.5);
t_ta = 0.5 * t_vir;
- t_i = t_ta * pow(delta_i, 1.5) / (3*M_PI/4);
- rho_ci = 1.0 / (6 * M_PI * Grav * t_i*t_i);
- r_i = pow(mass / (4*M_PI/3 * rho_ci), 1.0/3);
+ t_i = t_ta * POW(delta_i, 1.5) / (3*pi/4);
+ rho_ci = 1.0 / (6 * pi * GravConst * t_i*t_i);
+ r_i = POW(mass / (4*pi/3.0 * rho_ci), 1.0/3.0);
- dm_rotvel_corr = 2 * pow(2*r_init/r_i, 2);
+ dm_rotvel_corr = 2 * POW(2*r_init/r_i, 2);
// dm_rotvel_corr = 1;
- r_ta_now = pow((3*M_PI/4), -8.0/9) * pow(delta_i, 1.0/3) * r_i
- * pow(t_init / t_i, 8.0/9);
+ r_ta_now = POW((3*pi/4), -8.0/9) * POW(delta_i, 1.0/3.0) * r_i
+ * POW(t_init / t_i, 8.0/9);
for (i = 0; i < NTHETA; i++) {
local_Vr[NTHETA-i-1] = (Lambda[i] * r_ta_now / t_init) *
sin(Theta[i]) * (Theta[i] - sin(Theta[i])) /
- pow(1.0 - cos(Theta[i]), 2);
+ POW(1.0 - cos(Theta[i]), 2);
Vl = local_Vr[NTHETA-i-1] / (r_ta_now / t_init);
// Subtract Hubble flow and convert to code units
- local_Vr[NTHETA-i-1] = (Vl - (2.0/3)*Lambda[i]) * r_ta_now / t_init;
+ local_Vr[NTHETA-i-1] = (Vl - (2.0/3.0)*Lambda[i]) * r_ta_now / t_init;
local_Vr[NTHETA-i-1] /= VelocityUnits;
Chi = 1.0 - 1.5 * Vl / Lambda[i];
- local_rho[NTHETA-i-1] = pow(little_d[i],2) / pow(beta[i],3) /
+ local_rho[NTHETA-i-1] = POW(little_d[i],2) / POW(beta[i],3) /
(1.0 + 3*Chi);
local_radius[NTHETA-i-1] = Lambda[i] * r_ta_now / LengthUnits;
if (local_radius[NTHETA-i-1] > r_init)
@@ -1211,7 +1206,7 @@
exterior_rho[i] = density;
}
// printf("%"ISYM": r = %"GSYM" pc, v_r = %"GSYM" km/s, rho = %"GSYM"\n",
-// i, radius_vr[i]*LengthUnits/3.086e18, Vr[i]*VelocityUnits/1e5,
+// i, radius_vr[i]*LengthUnits/pc, Vr[i]*VelocityUnits/1e5,
// exterior_rho[i]);
}
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/Grid_ComputeAccelerationFieldExternal.C
--- a/src/enzo/Grid_ComputeAccelerationFieldExternal.C
+++ b/src/enzo/Grid_ComputeAccelerationFieldExternal.C
@@ -185,8 +185,8 @@
((log(1.0+x )-x /(1.0+x )) /
(log(1.0+1.0)-1.0/(1.0+1.0))) /
POW(radius*LengthUnits, 2.0) / AccelUnits +
- ClusterSMBHBCG*POW((POW(POW(radius*LengthUnits/(1.0e-3*Mpc),0.5975)/3.206e-7,0.9) +
- POW(POW(radius*LengthUnits/(1.0e-3*Mpc), 1.849)/1.861e-6, 0.9)), -1.0/0.9) / AccelUnits +
+ ClusterSMBHBCG*POW((POW(POW(radius*LengthUnits/(1.0e-3*Mpc_cm),0.5975)/3.206e-7,0.9) +
+ POW(POW(radius*LengthUnits/(1.0e-3*Mpc_cm), 1.849)/1.861e-6, 0.9)), -1.0/0.9) / AccelUnits +
GravConst*SolarMass*ClusterSMBHMass / POW(radius*LengthUnits, 2) / AccelUnits; // + BCG + BH mass
/*Bondi*/
if(ClusterSMBHCalculateGasMass == 4){
@@ -194,8 +194,8 @@
((log(1.0+x )-x /(1.0+x )) /
(log(1.0+1.0)-1.0/(1.0+1.0))) /
POW(radius*LengthUnits, 2.0) / AccelUnits +
- ClusterSMBHBCG*POW((POW(POW(radius*LengthUnits/(1.0e-3*Mpc),0.5975)/3.206e-7,0.9) +
- POW(POW(radius*LengthUnits/(1.0e-3*Mpc), 1.849)/1.861e-6, 0.9)), -1.0/0.9) / AccelUnits +
+ ClusterSMBHBCG*POW((POW(POW(radius*LengthUnits/(1.0e-3*Mpc_cm),0.5975)/3.206e-7,0.9) +
+ POW(POW(radius*LengthUnits/(1.0e-3*Mpc_cm), 1.849)/1.861e-6, 0.9)), -1.0/0.9) / AccelUnits +
GravConst*SolarMass*ClusterSMBHMass/POW(radius*LengthUnits - 2.0*GravConst*SolarMass/POW(clight,2), 2)/ AccelUnits;
}
/*Elliptical Galaxy Fixed Gravity*/
@@ -204,7 +204,7 @@
((log(1.0+x )-x /(1.0+x )) /
(log(1.0+1.0)-1.0/(1.0+1.0))) /
POW(radius*LengthUnits, 2.0) / AccelUnits +
- GravConst*(ClusterSMBHBCG*SolarMass*1.0e11)/POW(radius*LengthUnits+EllipticalGalaxyRe*1.0e-3*Mpc/1.8153, 2)/AccelUnits +
+ GravConst*(ClusterSMBHBCG*SolarMass*1.0e11)/POW(radius*LengthUnits+EllipticalGalaxyRe*1.0e-3*Mpc_cm/1.8153, 2)/AccelUnits +
GravConst*SolarMass*ClusterSMBHMass/POW(radius*LengthUnits - 2.0*GravConst*SolarMass/POW(clight,2), 2)/ AccelUnits;
}
}
@@ -304,8 +304,8 @@
((log(1.0+x )-x /(1.0+x )) /
(log(1.0+1.0)-1.0/(1.0+1.0))) /
POW(radius*LengthUnits, 2.0) / AccelUnits +
- ClusterSMBHBCG*POW((POW(POW(radius*LengthUnits/(1.0e-3*Mpc),0.5975)/3.206e-7,0.9) +
- POW(POW(radius*LengthUnits/(1.0e-3*Mpc), 1.849)/1.861e-6, 0.9)), -1.0/0.9) / AccelUnits +
+ ClusterSMBHBCG*POW((POW(POW(radius*LengthUnits/(1.0e-3*Mpc_cm),0.5975)/3.206e-7,0.9) +
+ POW(POW(radius*LengthUnits/(1.0e-3*Mpc_cm), 1.849)/1.861e-6, 0.9)), -1.0/0.9) / AccelUnits +
GravConst*SolarMass*ClusterSMBHMass / POW(radius*LengthUnits, 2) / AccelUnits; // + BCG + BH mass
/*Bondi*/
if(ClusterSMBHCalculateGasMass == 4){
@@ -313,8 +313,8 @@
((log(1.0+x )-x /(1.0+x )) /
(log(1.0+1.0)-1.0/(1.0+1.0))) /
POW(radius*LengthUnits, 2.0) / AccelUnits +
- ClusterSMBHBCG*POW((POW(POW(radius*LengthUnits/(1.0e-3*Mpc),0.5975)/3.206e-7,0.9) +
- POW(POW(radius*LengthUnits/(1.0e-3*Mpc), 1.849)/1.861e-6, 0.9)), -1.0/0.9) / AccelUnits +
+ ClusterSMBHBCG*POW((POW(POW(radius*LengthUnits/(1.0e-3*Mpc_cm),0.5975)/3.206e-7,0.9) +
+ POW(POW(radius*LengthUnits/(1.0e-3*Mpc_cm), 1.849)/1.861e-6, 0.9)), -1.0/0.9) / AccelUnits +
GravConst*SolarMass*ClusterSMBHMass/POW(radius*LengthUnits - 2.0*GravConst*SolarMass/POW(clight,2), 2)/ AccelUnits;
}
/*Elliptical Galaxy Fixed Gravity*/
@@ -323,7 +323,7 @@
((log(1.0+x )-x /(1.0+x )) /
(log(1.0+1.0)-1.0/(1.0+1.0))) /
POW(radius*LengthUnits, 2.0) / AccelUnits +
- GravConst*(ClusterSMBHBCG*SolarMass*1.0e11)/POW(radius*LengthUnits+EllipticalGalaxyRe*1.0e-3*Mpc/1.8153, 2)/AccelUnits +
+ GravConst*(ClusterSMBHBCG*SolarMass*1.0e11)/POW(radius*LengthUnits+EllipticalGalaxyRe*1.0e-3*Mpc_cm/1.8153, 2)/AccelUnits +
GravConst*SolarMass*ClusterSMBHMass/POW(radius*LengthUnits - 2.0*GravConst*SolarMass/POW(clight,2), 2)/ AccelUnits;
}
}
@@ -427,21 +427,21 @@
rcyl = sqrt(xpos1*xpos1 + ypos1*ypos1 + zpos1*zpos1);
radius = radius*LengthUnits;
rcyl = rcyl*LengthUnits;
- accelsph = (GravConst)*MBulge*SolarMass/POW(radius+rBulge*Mpc,2)
- + pi*GravConst*densDMConst*POW(rDMConst*Mpc,3)/POW(radius,2)
- *(-2.0*atan(radius/rDMConst/Mpc)
- +2.0*log(1.0+radius/rDMConst/Mpc)
- +log(1.0+POW(radius/rDMConst/Mpc,2))
+ accelsph = (GravConst)*MBulge*SolarMass/POW(radius+rBulge*Mpc_cm,2)
+ + pi*GravConst*densDMConst*POW(rDMConst*Mpc_cm,3)/POW(radius,2)
+ *(-2.0*atan(radius/rDMConst/Mpc_cm)
+ +2.0*log(1.0+radius/rDMConst/Mpc_cm)
+ +log(1.0+POW(radius/rDMConst/Mpc_cm,2))
);
accelcylR = GravConst*MSDisk*SolarMass*rcyl/sqrt(POW(POW(rcyl,2)
- + POW(SDiskScaleHeightR*Mpc+sqrt(POW(zheight*LengthUnits,2)
- + POW(SDiskScaleHeightz*Mpc,2)),2),3));
+ + POW(SDiskScaleHeightR*Mpc_cm+sqrt(POW(zheight*LengthUnits,2)
+ + POW(SDiskScaleHeightz*Mpc_cm,2)),2),3));
accelcylz = GravConst*MSDisk*SolarMass/sqrt(POW(zheight*LengthUnits,2)
- + POW(SDiskScaleHeightz*Mpc,2))*zheight*LengthUnits/sqrt(POW(POW(rcyl,2)
- + POW(SDiskScaleHeightR*Mpc+sqrt(POW(zheight*LengthUnits,2)
- + POW(SDiskScaleHeightz*Mpc,2)),2),3))
- *( SDiskScaleHeightR*Mpc+sqrt(POW(zheight*LengthUnits,2)
- + POW(SDiskScaleHeightz*Mpc,2))
+ + POW(SDiskScaleHeightz*Mpc_cm,2))*zheight*LengthUnits/sqrt(POW(POW(rcyl,2)
+ + POW(SDiskScaleHeightR*Mpc_cm+sqrt(POW(zheight*LengthUnits,2)
+ + POW(SDiskScaleHeightz*Mpc_cm,2)),2),3))
+ *( SDiskScaleHeightR*Mpc_cm+sqrt(POW(zheight*LengthUnits,2)
+ + POW(SDiskScaleHeightz*Mpc_cm,2))
)/AccelUnits;
accelsph = (radius ==0.0?0.0:fabs(accelsph )/(radius/LengthUnits)/AccelUnits);
@@ -500,21 +500,21 @@
radius = radius*LengthUnits;
rcyl = rcyl*LengthUnits;
- accelsph = (GravConst)*MBulge*SolarMass/POW(radius+rBulge*Mpc,2)
- + pi*GravConst*densDMConst*POW(rDMConst*Mpc,3)/POW(radius,2)
- * (-2.0*atan(radius/rDMConst/Mpc)
- +2.0*log(1.0+radius/rDMConst/Mpc)
- +log(1.0+POW(radius/rDMConst/Mpc,2))
+ accelsph = (GravConst)*MBulge*SolarMass/POW(radius+rBulge*Mpc_cm,2)
+ + pi*GravConst*densDMConst*POW(rDMConst*Mpc_cm,3)/POW(radius,2)
+ * (-2.0*atan(radius/rDMConst/Mpc_cm)
+ +2.0*log(1.0+radius/rDMConst/Mpc_cm)
+ +log(1.0+POW(radius/rDMConst/Mpc_cm,2))
);
accelcylR = GravConst*MSDisk*SolarMass*rcyl/sqrt(POW(POW(rcyl,2)
- + POW(SDiskScaleHeightR*Mpc+sqrt(POW(zheight*LengthUnits,2)
- + POW(SDiskScaleHeightz*Mpc,2)),2),3));
+ + POW(SDiskScaleHeightR*Mpc_cm+sqrt(POW(zheight*LengthUnits,2)
+ + POW(SDiskScaleHeightz*Mpc_cm,2)),2),3));
accelcylz = GravConst*MSDisk*SolarMass/sqrt(POW(zheight*LengthUnits,2)
- + POW(SDiskScaleHeightz*Mpc,2))*zheight*LengthUnits/sqrt(POW(POW(rcyl,2)
- + POW(SDiskScaleHeightR*Mpc+sqrt(POW(zheight*LengthUnits,2)
- + POW(SDiskScaleHeightz*Mpc,2)),2),3))
- *( SDiskScaleHeightR*Mpc+sqrt(POW(zheight*LengthUnits,2)
- + POW(SDiskScaleHeightz*Mpc,2))
+ + POW(SDiskScaleHeightz*Mpc_cm,2))*zheight*LengthUnits/sqrt(POW(POW(rcyl,2)
+ + POW(SDiskScaleHeightR*Mpc_cm+sqrt(POW(zheight*LengthUnits,2)
+ + POW(SDiskScaleHeightz*Mpc_cm,2)),2),3))
+ *( SDiskScaleHeightR*Mpc_cm+sqrt(POW(zheight*LengthUnits,2)
+ + POW(SDiskScaleHeightz*Mpc_cm,2))
)/AccelUnits;
accelsph = (radius ==0.0?0.0:fabs(accelsph )/(radius/LengthUnits)/AccelUnits);
@@ -561,16 +561,15 @@
FLOAT xc = 0.5, yc = 0.5, zc = 0.5;
double rs = rvir / c;
- double Mvir = 4.0*M_PI*rhoc*POW(rs,3)*(log(1.0+c)-c/(1.0+c));
+ double Mvir = 4.0*pi*rhoc*POW(rs,3)*(log(1.0+c)-c/(1.0+c));
float DensityUnits = 1.0, LengthUnits = 1.0, TemperatureUnits = 1,
TimeUnits = 1.0, VelocityUnits = 1.0;
GetUnits(&DensityUnits, &LengthUnits, &TemperatureUnits,
&TimeUnits, &VelocityUnits, &MassUnits, Time);
double AccelerationUnits = LengthUnits / POW(TimeUnits,2);
- double CGSGravConst = 6.672e-8;
- printf("rhoc=%g, rvir=%g, Mvir=%g\n", rhoc, rvir, Mvir/1.989e33);
+ printf("rhoc=%g, rvir=%g, Mvir=%g\n", rhoc, rvir, Mvir/SolarMass);
FLOAT x, y, z, xpos, ypos, zpos, r;
int n = 0;
@@ -603,12 +602,12 @@
if (r < rvir/LengthUnits) {
x1 = r*LengthUnits/rs;
- M = 4.0*M_PI*rhoc*POW(rs,3)*(log(1.0+x1)-x1/(1.0+x1));
+ M = 4.0*pi*rhoc*POW(rs,3)*(log(1.0+x1)-x1/(1.0+x1));
}
else {
M = Mvir;
}
- g = CGSGravConst*M/POW(r*LengthUnits,2);
+ g = GravConst*M/POW(r*LengthUnits,2);
g /= AccelerationUnits;
if (dim == 0) {
AccelerationField[0][n] += -g*xpos/r;
@@ -636,12 +635,12 @@
if (r < rvir/LengthUnits) {
x1 = r*LengthUnits/rs;
- M = 4.0*M_PI*rhoc*POW(rs,3)*(log(1.0+x1)-x1/(1.0+x1));
+ M = 4.0*pi*rhoc*POW(rs,3)*(log(1.0+x1)-x1/(1.0+x1));
}
else {
M = Mvir;
}
- g = CGSGravConst*M/POW(r*LengthUnits,2);
+ g = GravConst*M/POW(r*LengthUnits,2);
g /= AccelerationUnits;
ParticleAcceleration[0][i] += -g*xpos/r;
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/Grid_ComputeCoolingTime.C
--- a/src/enzo/Grid_ComputeCoolingTime.C
+++ b/src/enzo/Grid_ComputeCoolingTime.C
@@ -25,6 +25,7 @@
#include "Grid.h"
#include "fortran.def"
#include "CosmologyParameters.h"
+#include "phys_constants.h"
/* This parameter controls whether the cooling function recomputes
the metal cooling rates. It is reset by RadiationFieldUpdate. */
@@ -300,8 +301,7 @@
#ifdef TRANSFER
/* unit conversion from Enzo RT units to CGS */
- const float ev2erg = 1.60217653E-12;
- float rtunits = ev2erg / TimeUnits;
+ float rtunits = erg_eV / TimeUnits;
if ( RadiativeTransfer ){
my_fields.RT_HI_ionization_rate = BaryonField[kphHINum];
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/Grid_ComputeElementalDensity.C
--- a/src/enzo/Grid_ComputeElementalDensity.C
+++ b/src/enzo/Grid_ComputeElementalDensity.C
@@ -37,6 +37,7 @@
#include "ExternalBoundary.h"
#include "Grid.h"
#include "CosmologyParameters.h"
+#include "phys_constants.h"
/* function prototypes */
@@ -316,9 +317,9 @@
}
/* Loop over grid and compute emissivity
- (fh is hydrogen fraction by mass; 1.67e-24 is proton mass). */
+ (fh is hydrogen fraction by mass; mh is proton mass). */
- float fh = 0.76, ConvertToNumberDensity = DensityUnits/1.67e-24;
+ float fh = 0.76, ConvertToNumberDensity = DensityUnits/mh;
float nElement, nH, LogFraction, Metallicity, Fraction,
deld, delt, logd, logt, dd, dt;
int ilogd, ilogt;
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/Grid_ComputePhotonTimestep.C
--- a/src/enzo/Grid_ComputePhotonTimestep.C
+++ b/src/enzo/Grid_ComputePhotonTimestep.C
@@ -211,7 +211,7 @@
if (G_TotalNumberOfStars > 0 && minStarLifetime < 1e6) {
dtStar = minStarLifetime/NumberOfStepsInLifetime;
} else {
- dtStar = 3.1557e13*mindtNOstars/TimeUnits;
+ dtStar = Myr_s*mindtNOstars/TimeUnits;
}
}
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/Grid_ComputePhotonTimestepHII.C
--- a/src/enzo/Grid_ComputePhotonTimestepHII.C
+++ b/src/enzo/Grid_ComputePhotonTimestepHII.C
@@ -26,6 +26,7 @@
#include "ExternalBoundary.h"
#include "Grid.h"
#include "RadiativeTransferParameters.h"
+#include "phys_constants.h"
#define MAX_CHANGE 0.1
@@ -47,7 +48,6 @@
int i, j, k, dim, index, indexn, size;
float dt, this_dt, sigma_dx, *temperature;
- const double sigmaHI = 6.345e-18, mh = 1.673e-24;
dt = huge_number;
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/Grid_ComputePhotonTimestepTau.C
--- a/src/enzo/Grid_ComputePhotonTimestepTau.C
+++ b/src/enzo/Grid_ComputePhotonTimestepTau.C
@@ -26,6 +26,7 @@
#include "ExternalBoundary.h"
#include "Grid.h"
#include "RadiativeTransferParameters.h"
+#include "phys_constants.h"
#define MAX_CHANGE 0.5
@@ -45,7 +46,6 @@
int i, j, k, dim, index, size;
float dt, this_dt, sigma_dx, *temperature;
- const double sigmaHI = 6.345e-18, mh = 1.673e-24;
dt = huge_number;
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/Grid_ComputeXrayEmissivity.C
--- a/src/enzo/Grid_ComputeXrayEmissivity.C
+++ b/src/enzo/Grid_ComputeXrayEmissivity.C
@@ -30,6 +30,7 @@
#include "fortran.def"
#include "Grid.h"
#include "CosmologyParameters.h"
+#include "phys_constants.h"
/* function prototypes */
@@ -165,7 +166,7 @@
(fh is hydrogen fraction by mass). */
float temp, ne, nH, frac, fh = 0.76;
- float ConvertToNumberDensity = DensityUnits/1.67e-24;
+ float ConvertToNumberDensity = DensityUnits/mh;
float deltemp = (temp2-temp1)/float(NumberOfTemperatureBins-1);
for (i = 0; i < size; i++) {
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/Grid_ConductionBubbleInitialize.C
--- a/src/enzo/Grid_ConductionBubbleInitialize.C
+++ b/src/enzo/Grid_ConductionBubbleInitialize.C
@@ -22,6 +22,7 @@
#include "GridList.h"
#include "ExternalBoundary.h"
#include "Grid.h"
+#include "phys_constants.h"
// Function prototypes
int GetUnits (float *DensityUnits, float *LengthUnits,
@@ -43,7 +44,6 @@
static int ncells;
#define KEV_KELVIN 1.1604e+7
-#define KPC_CGS 3.0857e+21
#define DEFAULT_MU 0.6 // we assume total ionization
// Grid Initializer: all input values are in Enzo internal units _except_
@@ -127,7 +127,7 @@
nofr = new double[ncells];
Tofr = new double[ncells];
- dKdr_cgs = dKdr * KEV_KELVIN / KPC_CGS;
+ dKdr_cgs = dKdr * KEV_KELVIN / kpc_cm;
r_mid_cgs = r_mid * LengthUnits;
r_max_cgs = 2.0*r_mid_cgs;
@@ -147,7 +147,7 @@
// convert 1D arrays into Enzo internal units.
for(i=0; i<ncells; i++){
rad[i] /= LengthUnits; // convert to enzo distance
- nofr[i] *= DEFAULT_MU * 1.67e-24 / DensityUnits; // convert to enzo-unit density (from number density)
+ nofr[i] *= DEFAULT_MU * mh / DensityUnits; // convert to enzo-unit density (from number density)
Tofr[i] /= (TemperatureUnits*(Gamma-1.0)*DEFAULT_MU); // convert from temp to internal energy
}
@@ -331,9 +331,9 @@
- // g*mu*mp/kb;
+ // g*mu*mp/kboltz;
- bunch_of_constants = g*DEFAULT_MU*(1.67e-24)/(1.38e-16);
+ bunch_of_constants = g*DEFAULT_MU*(mh)/(kboltz);
/*
printf("n_mid, bunch_of_constants: %e %e %e\n",n_mid,bunch_of_constants,g);
This diff is so big that we needed to truncate the remainder.
Repository URL: https://bitbucket.org/enzo/enzo-dev/
--
This is a commit notification from bitbucket.org. You are receiving
this because you have the service enabled, addressing the recipient of
this email.
https://bitbucket.org/enzo/enzo-dev/commits/7d8fe75451e0/
Changeset: 7d8fe75451e0
Branch: week-of-code
User: jwise77
Date: 2019-02-21 09:43:17+00:00
Summary: Merged in aemerick/enzo-dev-constants (pull request #417)
Moving all local constants to those in header file
Approved-by: Nathan Goldbaum <ngol...@illinois.edu>
Approved-by: Brian O'Shea <os...@msu.edu>
Approved-by: John Regan <johnanth...@gmail.com>
Affected #: 194 files
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/AMRH5writer.C
--- a/src/enzo/AMRH5writer.C
+++ b/src/enzo/AMRH5writer.C
@@ -24,6 +24,7 @@
#include "typedefs.h"
#include "AMRH5writer.h"
#include "StarParticleData.h"
+#include "phys_constants.h"
AMRHDF5Writer::AMRHDF5Writer() :
h5DataType(H5T_NATIVE_FLOAT),
@@ -151,14 +152,14 @@
error = true;
}
- // The 3.14159 header is to test the endian since this is a binary
+ // The pi header is to test the endian since this is a binary
// file. The root cell spacing is necessary to calculate the dx for
// all other grids.
- float myPi = 3.14159;
char nFields8bit = (char) NumberOfAllFields;
char stag8bit = (char) stag;
int field;
- fwrite(&myPi, sizeof(float), 1, index_file);
+ double test_pi = pi;
+ fwrite(&test_pi, sizeof(float), 1, index_file);
fwrite(&root_dx, sizeof(float), 1, index_file);
fwrite(&nFields8bit, sizeof(char), 1, index_file);
fwrite(&stag8bit, sizeof(char), 1, index_file);
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/AdiabaticExpansionInitialize.C
--- a/src/enzo/AdiabaticExpansionInitialize.C
+++ b/src/enzo/AdiabaticExpansionInitialize.C
@@ -29,6 +29,7 @@
#include "Hierarchy.h"
#include "TopGridData.h"
#include "CosmologyParameters.h"
+#include "phys_constants.h"
/* Set the mean molecular mass as in Grid_ComputeTemperatureField.C */
@@ -135,7 +136,7 @@
ENZO_FAIL("Error in GetUnits.");
}
PressureUnits = DensityUnits * (LengthUnits/TimeUnits)*(LengthUnits/TimeUnits);
- MagneticUnits = sqrt(PressureUnits*4.0*M_PI);
+ MagneticUnits = sqrt(PressureUnits*4.0*pi);
CRUnits = DensityUnits * (LengthUnits/TimeUnits)*(LengthUnits/TimeUnits);
for (int dim = 0; dim < MAX_DIMENSION; dim++)
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/AdjustMustRefineParticlesRefineToLevel.C
--- a/src/enzo/AdjustMustRefineParticlesRefineToLevel.C
+++ b/src/enzo/AdjustMustRefineParticlesRefineToLevel.C
@@ -41,6 +41,7 @@
#include "TopGridData.h"
#include "LevelHierarchy.h"
#include "Star.h"
+#include "phys_constants.h"
int GetUnits(float *DensityUnits, float *LengthUnits,
float *TemperatureUnits, float *TimeUnits,
@@ -61,7 +62,7 @@
ENZO_FAIL("Your setup of assymmetrical topgrid is never tested.\n");
}
- const double pc = 3.086e18, ln2 = 0.69314718;
+ const double ln2 = 0.69314718;
int MustRefineParticlesRefineToLevel_prev = MustRefineParticlesRefineToLevel;
/* Compute Units. */
@@ -80,7 +81,7 @@
MustRefineParticlesRefineToLevel =
(int) floor( log( LengthUnits / MetaData->TopGridDims[0]
- / MustRefineParticlesRefineToLevelAutoAdjust / pc) / ln2);
+ / MustRefineParticlesRefineToLevelAutoAdjust / pc_cm) / ln2);
if (MyProcessorNumber == ROOT_PROCESSOR) {
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/AnalysisBaseClass.C
--- a/src/enzo/AnalysisBaseClass.C
+++ b/src/enzo/AnalysisBaseClass.C
@@ -15,6 +15,7 @@
#include "TopGridData.h"
#include "CosmologyParameters.h"
#include "AnalysisBaseClass.h"
+#include "phys_constants.h"
std::map<HierarchyEntry *, int> OriginalGridID;
#ifdef USE_HDF5_GROUPS
@@ -31,16 +32,8 @@
HierarchyEntry *topgrid,
int level,
FLOAT *Left,
- FLOAT *Right):
-#ifdef USE_BAD_VALUES
- G(6.67e-8),
- MPC_CM(3.0856e+24),
- MSOLAR_G(1.989e+33)
-#else
- G(6.6742e-8),
- MPC_CM(3.0856776e+24),
- MSOLAR_G(1.989e+33)
-#endif
+ FLOAT *Right)
+
{
int i;
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/AnalysisBaseClass.h
--- a/src/enzo/AnalysisBaseClass.h
+++ b/src/enzo/AnalysisBaseClass.h
@@ -33,10 +33,6 @@
protected:
- const float G;
- const float MPC_CM;
- const float MSOLAR_G;
-
TopGridData *MetaData;
HierarchyEntry *TopGrid;
int MaximumAnalysisLevel;
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/CalculateSubtractionParameters.C
--- a/src/enzo/CalculateSubtractionParameters.C
+++ b/src/enzo/CalculateSubtractionParameters.C
@@ -29,6 +29,7 @@
#include "TopGridData.h"
#include "LevelHierarchy.h"
#include "CommunicationUtilities.h"
+#include "phys_constants.h"
int CommunicationBroadcastValue(int *Value, int BroadcastProcessor);
@@ -43,8 +44,7 @@
float &Radius, double &Subtraction)
{
- const double pc = 3.086e18, Msun = 1.989e33, Grav = 6.673e-8, yr = 3.1557e7, Myr = 3.1557e13,
- k_b = 1.38e-16, m_h = 1.673e-24, c = 3.0e10, sigma_T = 6.65e-25, h=0.70;
+ const double h=0.70;
float mdot, AccretedMass, SafetyFactor;
float MassEnclosed = 0, Metallicity = 0, ColdGasMass = 0, OneOverRSquaredSum, AvgVelocity[MAX_DIMENSION];
@@ -129,7 +129,7 @@
density = star_CurrentGrid->ReturnBaryonField(DensNum)[index];
if (MultiSpecies == 0) {
- number_density = density * DensityUnits / (Mu * m_h);
+ number_density = density * DensityUnits / (Mu * mh);
mu = Mu;
} else {
number_density =
@@ -157,7 +157,7 @@
because only integer can be broadcasted at the moment;
anyway, this shouldn't make a huge difference overall. */
- c_s = (int)(sqrt(Gamma * k_b * temperature[index] / (mu * m_h)));
+ c_s = (int)(sqrt(Gamma * kboltz * temperature[index] / (mu * mh)));
fprintf(stdout, "index = %d, temp = %g, mu = %g, density = %g, number_density = %g, c_s = %d\n",
index, temperature[index], mu, density, number_density, c_s);
@@ -182,7 +182,7 @@
/* When calculating Bondi radius, let's not bother to get accurate temperature and mu
just use the MBHAccretionFixedTemperature and default mu */
- c_s = (int)(sqrt(Gamma * k_b * MBHAccretionFixedTemperature / (Mu * m_h)));
+ c_s = (int)(sqrt(Gamma * kboltz * MBHAccretionFixedTemperature / (Mu * mh)));
#endif
@@ -192,14 +192,14 @@
old_mass = (float)(star_mass);
mdot = isnan(star_last_accretion_rate) ? 0.0 : star_last_accretion_rate;
- AccretedMass = mdot * dtForThisStar * TimeUnits; //in Msun
+ AccretedMass = mdot * dtForThisStar * TimeUnits; //in SolarMass
if (MBHAccretionRadius > 0)
- Radius = MBHAccretionRadius * pc / LengthUnits;
+ Radius = MBHAccretionRadius * pc_cm / LengthUnits;
else {
SafetyFactor = -MBHAccretionRadius;
- Radius = SafetyFactor * 2.0 * Grav * old_mass * Msun / (c_s * c_s) / LengthUnits;
+ Radius = SafetyFactor * 2.0 * GravConst * old_mass * SolarMass / (c_s * c_s) / LengthUnits;
}
Radius = min(max(Radius, 4*StarLevelCellWidth), 100*StarLevelCellWidth);
@@ -284,8 +284,8 @@
here Subtraction is the density that should be subtracted.
if one wants to use this, change Grid_SAMFS.C accordingly. */
- Subtraction = -AccretedMass * Msun /
- (4*M_PI/3.0 * pow(Radius*LengthUnits, 3)) / DensityUnits;
+ Subtraction = -AccretedMass * SolarMass /
+ (4*pi/3.0 * pow(Radius*LengthUnits, 3)) / DensityUnits;
#endif
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/ClusterSMBHSumGasMass.C
--- a/src/enzo/ClusterSMBHSumGasMass.C
+++ b/src/enzo/ClusterSMBHSumGasMass.C
@@ -84,7 +84,7 @@
float ColdGasMassMsun=ClusterSMBHColdGasMass*MassUnits/SolarMass;
if (ClusterSMBHCalculateGasMass >1 ){ //2-calculate&remove,3-calculate&remove&re-orient,4-Bondi
if (ClusterSMBHCalculateGasMass == 3){
- ClusterSMBHJetDim = floor(Time*TimeUnits/(1.0e6*3.1557e7*ClusterSMBHJetPrecessionPeriod)); //ClusterSMBHJetPrecessionPeriod is now the Dim changing period.
+ ClusterSMBHJetDim = floor(Time*TimeUnits/(1.0e6*yr_s*ClusterSMBHJetPrecessionPeriod)); //ClusterSMBHJetPrecessionPeriod is now the Dim changing period.
}
if (ColdGasMassMsun > 0.000001) {
ClusterSMBHFeedbackSwitch = TRUE;
@@ -94,7 +94,7 @@
}
else {
ClusterSMBHJetMdot = (ColdGasMassMsun/(ClusterSMBHAccretionTime*1e6))/2.0; // AccretionTime from Myr to yr; reset Mdot, still in Msun/yr. Devide it by 2 because Mdot is for only one jet.
- ClusterSMBHJetEdot = (ClusterSMBHAccretionEpsilon*ClusterSMBHJetMdot * SolarMass/3.1557e7) * POW(clight,2)/1.0e44; //for one jet
+ ClusterSMBHJetEdot = (ClusterSMBHAccretionEpsilon*ClusterSMBHJetMdot * SolarMass/yr_s) * POW(clight,2)/1.0e44; //for one jet
}
}
else
@@ -108,7 +108,7 @@
ClusterSMBHFeedbackSwitch = (ColdGasMassMsun <= ClusterSMBHEnoughColdGas && LastClusterSMBHFeedbackSwitch == FALSE) ? FALSE : TRUE;
if (LastClusterSMBHFeedbackSwitch == FALSE && ClusterSMBHFeedbackSwitch == TRUE) {
- ClusterSMBHStartTime = Time + ClusterSMBHTramp*1.0e6*3.1557e7/TimeUnits;
+ ClusterSMBHStartTime = Time + ClusterSMBHTramp*1.0e6*yr_s/TimeUnits;
if (ClusterSMBHJetPrecessionPeriod < 0.00001 & ClusterSMBHJetPrecessionPeriod > -0.00001) //if precession off (set to 0), change the angle of the jets
ClusterSMBHJetAnglePhi += 0.5;
if (ClusterSMBHJetPrecessionPeriod < -0.00001) //if precession negative, change the jet dimension
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/CommunicationMergeStarParticle.C
--- a/src/enzo/CommunicationMergeStarParticle.C
+++ b/src/enzo/CommunicationMergeStarParticle.C
@@ -33,6 +33,7 @@
#include "TopGridData.h"
#include "Hierarchy.h"
#include "LevelHierarchy.h"
+#include "phys_constants.h"
/* function prototypes */
@@ -148,8 +149,7 @@
GetUnits(&DensityUnits, &LengthUnits, &TemperatureUnits, &TimeUnits,
&VelocityUnits, &MassUnits, 1.0);
- double Msun = 1.989e33;
- if(UsePhysicalUnit) MassUnits = DensityUnits*pow(LengthUnits,3)/Msun;
+ if(UsePhysicalUnit) MassUnits = DensityUnits*pow(LengthUnits,3)/SolarMass;
/*
if (debug) {
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/ConductionBubbleInitialize.C
--- a/src/enzo/ConductionBubbleInitialize.C
+++ b/src/enzo/ConductionBubbleInitialize.C
@@ -29,6 +29,7 @@
#include "Grid.h"
#include "Hierarchy.h"
#include "TopGridData.h"
+#include "phys_constants.h"
int GetUnits(float *DensityUnits, float *LengthUnits,
float *TemperatureUnits, float *TimeUnits,
@@ -108,7 +109,7 @@
total energy (that's not done automatically, since we don't know for sure
what the user is going to do). */
if (HydroMethod == MHD_RK){
- float MagneticUnits = sqrt(DensityUnits*4.0*M_PI)*VelocityUnits;
+ float MagneticUnits = sqrt(DensityUnits*4.0*pi)*VelocityUnits;
ConductionBubbleInitialUniformBField[0] /= MagneticUnits;
ConductionBubbleInitialUniformBField[1] /= MagneticUnits;
ConductionBubbleInitialUniformBField[2] /= MagneticUnits;
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/ConductionCloudInitialize.C
--- a/src/enzo/ConductionCloudInitialize.C
+++ b/src/enzo/ConductionCloudInitialize.C
@@ -25,6 +25,7 @@
#include "Grid.h"
#include "Hierarchy.h"
#include "TopGridData.h"
+#include "phys_constants.h"
int GetUnits(float *DensityUnits, float *LengthUnits,
float *TemperatureUnits, float *TimeUnits,
@@ -123,11 +124,11 @@
return FAIL;
}
- float Boltzmann = 1.38e-16, mu = 0.6, mh=1.67e-24;
+ float mu = 0.6;
if(ConductionCloudTemperature > 1.0){
- ConductionCloudTotalEnergy = (Boltzmann*ConductionCloudTemperature)/((Gamma - 1.0)*mu*mh);
+ ConductionCloudTotalEnergy = (kboltz*ConductionCloudTemperature)/((Gamma - 1.0)*mu*mh);
ConductionCloudTotalEnergy /= (VelocityUnits*VelocityUnits);
printf("ConductionCloudTotalEnergy is %e and ConductionCloudTemperature is %e\n\n",ConductionCloudTotalEnergy, ConductionCloudTemperature);
fflush(stdout);
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/ConductionTestInitialize.C
--- a/src/enzo/ConductionTestInitialize.C
+++ b/src/enzo/ConductionTestInitialize.C
@@ -27,6 +27,7 @@
#include "Grid.h"
#include "Hierarchy.h"
#include "TopGridData.h"
+#include "phys_constants.h"
int GetUnits(float *DensityUnits, float *LengthUnits,
float *TemperatureUnits, float *TimeUnits,
@@ -105,12 +106,12 @@
return FAIL;
}
- float Boltzmann = 1.38e-16, mu = 0.6, mh=1.67e-24;
+ float mu = 0.6;
// User can set
if(ConductionTestTemperature > 1.0){
- ConductionTestTotalEnergy = (Boltzmann*ConductionTestTemperature)/((Gamma - 1.0)*mu*mh);
+ ConductionTestTotalEnergy = (kboltz*ConductionTestTemperature)/((Gamma - 1.0)*mu*mh);
ConductionTestTotalEnergy /= (VelocityUnits*VelocityUnits);
ConductionTestGasEnergy = ConductionTestTotalEnergy;
printf("ConductionTestTotalEnergy is %e and ConductionTestTemperature is %e\n\n",ConductionTestTotalEnergy, ConductionTestTemperature);
@@ -118,7 +119,7 @@
}
if (HydroMethod == MHD_RK){
- float MagneticUnits = sqrt(DensityUnits*4.0*M_PI)*VelocityUnits;
+ float MagneticUnits = sqrt(DensityUnits*4.0*pi)*VelocityUnits;
ConductionTestInitialUniformBField[0] /= MagneticUnits;
ConductionTestInitialUniformBField[1] /= MagneticUnits;
ConductionTestInitialUniformBField[2] /= MagneticUnits;
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/CosmoIonizationInitialize.C
--- a/src/enzo/CosmoIonizationInitialize.C
+++ b/src/enzo/CosmoIonizationInitialize.C
@@ -29,6 +29,7 @@
#include "Hierarchy.h"
#include "LevelHierarchy.h"
#include "TopGridData.h"
+#include "phys_constants.h"
// default constants
#define DEFAULT_MU 0.6 // mean molecular mass
@@ -134,8 +135,7 @@
// convert input temperature to internal energy
RadHydroTemperature = max(RadHydroTemperature,MIN_TEMP); // enforce minimum
- float mp = 1.67262171e-24; // proton mass [g]
- float kb = 1.3806504e-16; // boltzmann constant [erg/K]
+
float nH, HI, HII, nHe, HeI, HeII, HeIII, ne, num_dens, mu;
if (RadHydroChemistry == 1) {
HI = 1.0 - RadHydroInitialFractionHII;
@@ -157,7 +157,7 @@
mu = 1.0/num_dens;
}
// compute the internal energy
- float RadHydroIEnergy = kb*RadHydroTemperature/mu/mp/(Gamma-1.0);
+ float RadHydroIEnergy = kboltz*RadHydroTemperature/mu/mh/(Gamma-1.0);
// set up the grid(s) on this level
HierarchyEntry *Temp = &TopGrid;
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/CosmologyGetUnits.C
--- a/src/enzo/CosmologyGetUnits.C
+++ b/src/enzo/CosmologyGetUnits.C
@@ -39,6 +39,7 @@
#include "ErrorExceptions.h"
#include "macros_and_parameters.h"
#include "CosmologyParameters.h"
+#include "phys_constants.h"
/* function prototypes */
@@ -66,7 +67,7 @@
*DensityUnits = 1.8788e-29*OmegaMatterNow*POW(HubbleConstantNow,2)*
POW(1 + CurrentRedshift,3);
- *LengthUnits = 3.085678e24*ComovingBoxSize/HubbleConstantNow/
+ *LengthUnits = Mpc_cm*ComovingBoxSize/HubbleConstantNow/
(1 + CurrentRedshift);
*TemperatureUnits = 1.81723e6*POW(ComovingBoxSize,2)*OmegaMatterNow*
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/CosmologySimulationInitialize.C
--- a/src/enzo/CosmologySimulationInitialize.C
+++ b/src/enzo/CosmologySimulationInitialize.C
@@ -39,6 +39,7 @@
#include "CosmologyParameters.h"
#include "fortran.def"
#include "CommunicationUtilities.h"
+#include "phys_constants.h"
// Function prototypes
@@ -419,7 +420,7 @@
ENZO_FAIL("Error in GetUnits.");
}
PressureUnits = DensityUnits * VelocityUnits*VelocityUnits;
- MagneticUnits = sqrt(PressureUnits*4.0*M_PI);
+ MagneticUnits = sqrt(PressureUnits*4.0*pi);
for (int dim = 0; dim < MAX_DIMENSION; dim++) {
if (CosmologySimulationInitialUniformBField[dim] != 0.0 && HydroMethod != 4 && HydroMethod != 6)
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/ExternalBoundary_SetWengenCollidingFlowBoundary.C
--- a/src/enzo/ExternalBoundary_SetWengenCollidingFlowBoundary.C
+++ b/src/enzo/ExternalBoundary_SetWengenCollidingFlowBoundary.C
@@ -27,6 +27,7 @@
#include "Grid.h"
#include "MHD2DTestGlobalData.h"
#include "hydro_rk/EOS.h"
+#include "phys_constants.h"
/* function prototypes */
@@ -117,7 +118,7 @@
EOSSoundSpeed*EOSSoundSpeed*rho ;
EOS(pres, rho, eintl, h, cs, dpdrho, dpde, 0, 1); // compute eintl
ramp = 1./(1.+exp(-2/RampWidth*(y-0.5)));
- f = cos(2.*M_PI*x*10.)*exp(-fabs(y-0.5)*10.)*cos(2.*M_PI*x*3);
+ f = cos(2.*pi*x*10.)*exp(-fabs(y-0.5)*10.)*cos(2.*pi*x*3);
vx = f * (vxl + ramp*(vxu-vxl)) ;
vy = vyl + ramp*(vyu - vyl);
bx = (Bxl+ ramp*(Bxu-Bxl)) ;
@@ -152,7 +153,7 @@
EOSSoundSpeed*EOSSoundSpeed*rho ;
EOS(pres, rho, eintl, h, cs, dpdrho, dpde, 0, 1); // compute eintl
ramp = 1./(1.+exp(-2/RampWidth*(y-0.5)));
- f = cos(2.*M_PI*x*10.)*exp(-fabs(y-0.5)*10.)*cos(2.*M_PI*x*3);
+ f = cos(2.*pi*x*10.)*exp(-fabs(y-0.5)*10.)*cos(2.*pi*x*3);
vx = f * (vxl + ramp*(vxu-vxl)) ;
vy = vyl + ramp*(vyu - vyl);
bx = (Bxl+ ramp*(Bxu-Bxl)) ;
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/FOF.C
--- a/src/enzo/FOF.C
+++ b/src/enzo/FOF.C
@@ -26,6 +26,7 @@
#include "macros_and_parameters.h"
#include "typedefs.h"
#include "global_data.h"
+#include "phys_constants.h"
#include "Fluxes.h"
#include "GridList.h"
@@ -148,9 +149,9 @@
void set_units(FOFData &AllVars) {
- AllVars.UnitLength_in_cm = 3.085678e21;
- AllVars.UnitMass_in_g = 1.989e43;
- AllVars.UnitVelocity_in_cm_per_s = 1.0e5;
+ AllVars.UnitLength_in_cm = kpc_cm;
+ AllVars.UnitMass_in_g = SolarMass * 1.0E10;
+ AllVars.UnitVelocity_in_cm_per_s = km_cm;
AllVars.Theta = 0.8; /* opening angle for potential computation */
AllVars.DesDensityNgb = 32;
@@ -1395,7 +1396,7 @@
FILE *in, *out;
int i, Ngroups, *Npart_g;
float *CM, *Mass;
- double SolarMass = UnitMass_in_g / 1.989e33;
+ double Msun = UnitMass_in_g / SolarMass;
float RootBoxSize[3];
if (MyProcessorNumber == 0) {
@@ -1427,7 +1428,7 @@
"Halo", "Npart", "Mass [Msun]", "CM(x)", "CM(y)", "CM(z)");
for (i = 0; i < Ngroups; i++)
fprintf(out, "%8d %12d %12.6g %12.6f %12.6f %12.6f\n",
- i, Npart_g[i], Mass[i]*SolarMass,
+ i, Npart_g[i], Mass[i]*Msun,
CM[3*i+0] / RootBoxSize[0] + leftEdge[0],
CM[3*i+1] / RootBoxSize[1] + leftEdge[1],
CM[3*i+2] / RootBoxSize[2] + leftEdge[2]);
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/FOF_Initialize.C
--- a/src/enzo/FOF_Initialize.C
+++ b/src/enzo/FOF_Initialize.C
@@ -33,6 +33,7 @@
#include "Hierarchy.h"
#include "LevelHierarchy.h"
#include "CosmologyParameters.h"
+#include "phys_constants.h"
#include "FOF_allvars.h"
#include "FOF_nrutil.h"
@@ -82,13 +83,13 @@
D.Time = a / (1 + InitialRedshift);
}
else {
- D.BoxSize = LengthUnits / 3.086e21;
+ D.BoxSize = LengthUnits / kpc_cm;
D.Time = 1.0;
}
// Critical density in units of Msun / kpc^3
D.RhoCritical0 = 1.4775867e31 *
- ((3 * pow(100 * HubbleConstantNow / 3.086e19, 2)) / (8 * M_PI * GRAVITY));
+ ((3.0 * POW(100 * HubbleConstantNow / 3.086e19, 2)) / (8.0 * pi * GRAVITY));
//D.RhoCritical /= pow(D.Time, 3);
// Sometimes MassUnits is infinite (in cgs) when using single
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/FOF_ngbtree.C
--- a/src/enzo/FOF_ngbtree.C
+++ b/src/enzo/FOF_ngbtree.C
@@ -12,7 +12,7 @@
#include "macros_and_parameters.h"
#include "typedefs.h"
#include "global_data.h"
-
+#include "phys_constants.h"
#include "FOF_allvars.h"
#include "FOF_ngbtree.h"
@@ -54,7 +54,7 @@
break;
}
- sr = th->len * pow((3.0/(4*M_PI)*SECFACTOR) * desngb /
+ sr = th->len * pow((3.0/(4*pi)*SECFACTOR) * desngb /
((float)(th->count)), 1.0/3);
} // ENDELSE
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/FOF_properties.C
--- a/src/enzo/FOF_properties.C
+++ b/src/enzo/FOF_properties.C
@@ -6,7 +6,7 @@
#include "macros_and_parameters.h"
#include "typedefs.h"
#include "global_data.h"
-
+#include "phys_constants.h"
#include "FOF_allvars.h"
#include "FOF_proto.h"
@@ -134,7 +134,7 @@
// kpc^3, as is rho_crit.
e2 = D.OmegaLambda + D.Omega * pow(D.Time, -3.0);
// 1/a^3 factor converts comoving radius to proper.
- factor = 1e10 / (4*M_PI/3.0) * pow(D.Time, -3.0);
+ factor = 1e10 / (4*pi/3.0) * pow(D.Time, -3.0);
rho178 = 178.0 * D.RhoCritical0 * e2;
len4 = len/4;
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/FSProb_ComputeOpacityLW.C
--- a/src/enzo/FSProb_ComputeOpacityLW.C
+++ b/src/enzo/FSProb_ComputeOpacityLW.C
@@ -21,12 +21,12 @@
************************************************************************/
#ifdef TRANSFER
#include "FSProb.h"
+#include "phys_constants.h"
int FSProb::ComputeOpacityLW(float *H2Density)
{
const float H2ISigma = 3.71e-18; // cm^-2
- const float mh = 1.673e-24;
double factor = H2ISigma * (DenUnits/mh);
int i, j, k, dim, size;
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/FreeExpansionInitialize.C
--- a/src/enzo/FreeExpansionInitialize.C
+++ b/src/enzo/FreeExpansionInitialize.C
@@ -33,6 +33,7 @@
#include "Grid.h"
#include "Hierarchy.h"
#include "TopGridData.h"
+#include "phys_constants.h"
int GetUnits(float *DensityUnits, float *LengthUnits,
float *TemperatureUnits, float *TimeUnits,
@@ -128,7 +129,7 @@
&VelocityUnits, 0.0);
PressureUnits = DensityUnits * (LengthUnits/TimeUnits)*(LengthUnits/TimeUnits);
- MagneticUnits = sqrt(PressureUnits*4.0*M_PI);
+ MagneticUnits = sqrt(PressureUnits*4.0*pi);
for (dim = 0; dim < MAX_DIMENSION; dim++)
FreeExpansionBField[dim] /= MagneticUnits;
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/FreezeRateData.C
--- a/src/enzo/FreezeRateData.C
+++ b/src/enzo/FreezeRateData.C
@@ -30,6 +30,7 @@
#include "ExternalBoundary.h"
#include "Grid.h"
#include "Hierarchy.h"
+#include "phys_constants.h"
/* function prototypes */
int GetUnits(float *DensityUnits, float *LengthUnits,
@@ -120,9 +121,7 @@
// put the temperature together
float mu = rho/num_density;
- float mp = 1.67262171e-24; // proton mass [g]
- float kb = 1.3806504e-16; // boltzmann constant [erg/K]
- float Temp = max((Gamma-1.0)*mu*mp*eint/kb, 1.0);
+ float Temp = max((Gamma-1.0)*mu*mh*eint/kboltz, 1.0);
// find temperature bin
float lamT = 3.15614e5/Temp;
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/Grid_AddExternalAcceleration.C
--- a/src/enzo/Grid_AddExternalAcceleration.C
+++ b/src/enzo/Grid_AddExternalAcceleration.C
@@ -25,7 +25,8 @@
#include "ExternalBoundary.h"
#include "Grid.h"
#include "SphericalInfall.h"
-
+#include "phys_constants.h"
+
/* function prototypes */
int CosmologyComputeExpansionFactor(FLOAT time, FLOAT *a, FLOAT *dadt);
@@ -39,7 +40,7 @@
if (ProblemType == 24) {
- float Pi = 3.14159, accel;
+ float accel;
FLOAT a = 1, dadt, rcubed, xpos, ypos = 0, zpos = 0;
int i, j, k, n = 0;
@@ -82,7 +83,7 @@
(i.e. 1/a^2 * a = 1/a). */
accel = GravitationalConstant*SphericalInfallFixedMass/
- (4.0*Pi*rcubed*a);
+ (4.0*pi*rcubed*a);
/* Apply force. */
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/Grid_AddFeedbackSphere.C
--- a/src/enzo/Grid_AddFeedbackSphere.C
+++ b/src/enzo/Grid_AddFeedbackSphere.C
@@ -28,6 +28,7 @@
#include "Grid.h"
#include "Hierarchy.h"
#include "CosmologyParameters.h"
+#include "phys_constants.h"
#ifdef CONFIG_BFLOAT_4
#define TOLERANCE 1e-06
@@ -51,8 +52,6 @@
{
const float WhalenMaxVelocity = 35; // km/s
- const double Msun = 1.989e33;
- const double c = 3.0e10;
int dim, i, j, k, index;
int sx, sy, sz;
@@ -117,7 +116,7 @@
if (MetalNum > 0 && SNColourNum > 0 && cstar->type == PopIII)
MetalNum = SNColourNum;
- float BubbleVolume = (4.0 * M_PI / 3.0) * radius * radius * radius;
+ float BubbleVolume = (4.0 * pi / 3.0) * radius * radius * radius;
/***********************************************************************
SUPERNOVAE
@@ -144,7 +143,7 @@
// Correct if the volume with 27 cells is larger than the energy bubble volume
#ifdef UNUSED
float BoxVolume = 27 * CellWidth[0][0] * CellWidth[0][0] * CellWidth[0][0];
- float BubbleVolume = (4.0 * M_PI / 3.0) * radius * radius * radius;
+ float BubbleVolume = (4.0 * pi / 3.0) * radius * radius * radius;
//printf("BoxVolume = %lg, BubbleVolume = %lg\n", BoxVolume, BubbleVolume);
if (BoxVolume > BubbleVolume) {
//printf("Reducing ejecta density by %g\n", BubbleVolume / BoxVolume);
@@ -174,7 +173,7 @@
// printf("grid::AFS: before: cstar->Mass = %lf\n", cstar->Mass);
if (cstar->FeedbackFlag != SUPERNOVA) {
float old_mass = (float)(cstar->Mass);
- cstar->Mass -= EjectaDensity * DensityUnits * BubbleVolume * pow(LengthUnits,3.0) / Msun;
+ cstar->Mass -= EjectaDensity * DensityUnits * BubbleVolume * pow(LengthUnits,3.0) / SolarMass;
float frac = old_mass / cstar->Mass;
cstar->vel[0] *= frac;
cstar->vel[1] *= frac;
@@ -466,7 +465,7 @@
metallicity_inside = 0.0, colour_inside = 0.0, rho_inside, rho_metal_inside, rho_colour_inside;
float m_cell_edge = 0.0, metal_cell_edge = 0.0, colour_cell_edge = 0.0,
metallicity_edge = 0.0, colour_edge = 0.0, rho_jet, rho_metal_jet, rho_colour_jet;
- float L_x, L_y, L_z, L_s, nx_L = 0.0, ny_L = 0.0, nz_L = 0.0, costheta = cos(3.1415926/3.9);
+ float L_x, L_y, L_z, L_s, nx_L = 0.0, ny_L = 0.0, nz_L = 0.0, costheta = cos(pi/3.9);
float EjectaMass, EjectaMetalMass, MBHJetsVelocity;
if (cstar->FeedbackFlag == MBH_JETS) {
@@ -492,7 +491,7 @@
/* Calculate mass that has accumulated thus far; since EjectaDensity is calculated
for isotropic MBH_THERMAL, we use it only to compute EjectaMass, not apply it directly. */
- cstar->NotEjectedMass += EjectaDensity * DensityUnits * BubbleVolume * pow(LengthUnits,3.0) / Msun;
+ cstar->NotEjectedMass += EjectaDensity * DensityUnits * BubbleVolume * pow(LengthUnits,3.0) / SolarMass;
// fprintf(stdout, "d1, d2, d3, i, j, k = %d %d %d / %d %d %d\n",
// GridDimension[0], GridDimension[1], GridDimension[2], i,j,k);
@@ -526,7 +525,7 @@
/* Find ejecta mass */
- EjectaMass = cstar->NotEjectedMass * Msun / DensityUnits / pow(LengthUnits,3.0);
+ EjectaMass = cstar->NotEjectedMass * SolarMass / DensityUnits / pow(LengthUnits,3.0);
EjectaMetalMass = EjectaMass * MBHFeedbackMetalYield;
OutputWhenJetsHaveNotEjected = FALSE;
@@ -561,8 +560,8 @@
//add random noise to theta and phi
srand(time(NULL));
- theta += MaximumNoiseAngle * M_PI / 180.0 * ((2.0*(float)rand()/((float)(RAND_MAX)+(float)(1))) - 1.0);
- phi += MaximumNoiseAngle * M_PI / 180.0 * ((2.0*(float)rand()/((float)(RAND_MAX)+(float)(1))) - 1.0);
+ theta += MaximumNoiseAngle * pi / 180.0 * ((2.0*(float)rand()/((float)(RAND_MAX)+(float)(1))) - 1.0);
+ phi += MaximumNoiseAngle * pi / 180.0 * ((2.0*(float)rand()/((float)(RAND_MAX)+(float)(1))) - 1.0);
//get back to Cartesian coordinate; some tricks needed to preserve the signs of nx_L and ny_L
nx_L = sign(nx_L) * fabs(cos(theta))*sin(phi);
@@ -620,7 +619,7 @@
} // END jj-direction
} // END kk-direction
-// printf("EjectaM in Msun = %g, EjectaM = %g, EjectaMetalM = %g, m_cell_edge = %g, n_cell_edge = %d\n",
+// printf("EjectaM in SolarMass = %g, EjectaM = %g, EjectaMetalM = %g, m_cell_edge = %g, n_cell_edge = %d\n",
// cstar->NotEjectedMass, EjectaMass, EjectaMetalMass, m_cell_edge, n_cell_edge);
/* Calculate the jet density */
@@ -654,13 +653,13 @@
Note that EjectaThermalEnergy is never used; MBHFeedbackEnergyCoupling
should now be calculated considering gravitational redshift (Kim et al. 2010) */
- MBHJetsVelocity = c * sqrt( 2 * MBHFeedbackEnergyCoupling * MBHFeedbackRadiativeEfficiency
+ MBHJetsVelocity = clight * sqrt( 2 * MBHFeedbackEnergyCoupling * MBHFeedbackRadiativeEfficiency
/ MBHFeedbackMassEjectionFraction ) / VelocityUnits;
- if (MBHJetsVelocity * VelocityUnits > 0.99*c) {
+ if (MBHJetsVelocity * VelocityUnits > 0.99*clight) {
ENZO_VFAIL("grid::AddFS: MBHJetsVelocity is ultra-relativistic! (%g/ %g/ %g/ %g c)\n",
MBHFeedbackEnergyCoupling, MBHFeedbackRadiativeEfficiency,
- MBHFeedbackMassEjectionFraction, MBHJetsVelocity * VelocityUnits / c);
+ MBHFeedbackMassEjectionFraction, MBHJetsVelocity * VelocityUnits / clight);
}
/* Finally, add the jet feedback at the edges (outer part of the supercell) */
@@ -757,7 +756,7 @@
// old_mass, cstar->Mass, old_vel1, cstar->vel[1]);
fprintf(stdout, "grid::AddFS: jets injected (EjectaM = %g Ms, JetsVelocity = %g, grid v = %g, rho_jet = %g) along n_L = (%g, %g, %g)\n",
- EjectaMass * DensityUnits * pow(LengthUnits,3.0) / Msun,
+ EjectaMass * DensityUnits * pow(LengthUnits,3.0) / SolarMass,
MBHJetsVelocity, BaryonField[Vel3Num][n_cell_edge-1], rho_jet, nx_L, ny_L, nz_L);
} // END MBH_JETS
@@ -854,7 +853,7 @@
float kph, kheat;
sigma_HI *= (double) TimeUnits / ((double)LengthUnits * (double)LengthUnits)
- / (4.0 * M_PI);
+ / (4.0 * pi);
kph = (float) (Q_HI * sigma_HI);
kheat = kph * deltaE;
@@ -988,7 +987,7 @@
/* Now it's done, unmark. */
//cstar->FeedbackFlag = NO_FEEDBACK;
-
+
return SUCCESS;
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/Grid_AddH2DissociationFromSources.C
--- a/src/enzo/Grid_AddH2DissociationFromSources.C
+++ b/src/enzo/Grid_AddH2DissociationFromSources.C
@@ -29,12 +29,12 @@
#include "Hierarchy.h"
#include "CosmologyParameters.h"
#include "Star.h"
+#include "phys_constants.h"
#ifdef FUTUREAP
#include "ActiveParticle.h"
#include "ActiveParticle_RadiationParticle.h"
#include "ActiveParticle_SmartStar.h"
#endif
-#include "phys_constants.h"
#define THRESHOLD_DENSITY_DB36 1e14
#define THRESHOLD_DENSITY_DB37 5e14
@@ -59,6 +59,7 @@
int DensNum, GENum, TENum, Vel1Num, Vel2Num, Vel3Num, B1Num, B2Num, B3Num;
int DeNum, HINum, HIINum, HeINum, HeIINum, HeIIINum, HMNum, H2INum, H2IINum,
DINum, DIINum, HDINum;
+
double l_char = 0.0, N_H2 = 0.0;
double H2ISigma = 3.71e-18;
@@ -130,7 +131,7 @@
// where the gas is fully molecular. If the solver still breaks,
// then this parameter can be safely increased to ~100 AU or turning
// on H2 self-shielding with RadiationShield = 2.
- float dilutionRadius = 4.848e-6 * pc / (double) LengthUnits; // 1 AU
+ float dilutionRadius = 4.848e-6 * pc_cm / (double) LengthUnits; // 1 AU
float dilRadius2 = dilutionRadius * dilutionRadius;
float LightTravelDist = TimeUnits * clight / LengthUnits;
@@ -196,9 +197,9 @@
int TemperatureField = 0;
/* Pre-compute some quantities to speed things up */
TemperatureField = this->GetTemperatureFieldNumberForH2Shield();
- kdiss_r2 = (float) (LWLuminosity * H2ISigma / (4.0 * M_PI));
- kph_hm = (float) (IRLuminosity * HMSigma / (4.0 * M_PI));
- kdiss_H2II = (float) (H2IILuminosity * H2IISigma / (4.0 * M_PI));
+ kdiss_r2 = (float) (LWLuminosity * H2ISigma / (4.0 * pi));
+ kph_hm = (float) (IRLuminosity * HMSigma / (4.0 * pi));
+ kdiss_H2II = (float) (H2IILuminosity * H2IISigma / (4.0 * pi));
for (k = 0; k < ActiveDims[2]; k++) {
for (j = 0; j < ActiveDims[1]; j++) {
radius2_yz = ddr2[1][j] + ddr2[2][k];
@@ -291,7 +292,8 @@
double radius2, radius2_yz;
- kdiss_r2 = (float) (LWLuminosity * H2ISigma / (4.0 * M_PI));
+ kdiss_r2 = (float) (LWLuminosity * H2ISigma / (4.0 * pi));
+
for (k = 0; k < ActiveDims[2]; k++) {
for (j = 0; j < ActiveDims[1]; j++) {
radius2_yz = ddr2[1][j] + ddr2[2][k];
@@ -365,9 +367,9 @@
int TemperatureField = 0;
/* Pre-compute some quantities to speed things up */
TemperatureField = this->GetTemperatureFieldNumberForH2Shield();
- kdiss_r2 = (float) (LWLuminosity * H2ISigma / (4.0 * M_PI));
- kph_hm = (float) (IRLuminosity * HMSigma / (4.0 * M_PI));
- kdiss_H2II = (float) (H2IILuminosity * H2IISigma / (4.0 * M_PI));
+ kdiss_r2 = (float) (LWLuminosity * H2ISigma / (4.0 * pi));
+ kph_hm = (float) (IRLuminosity * HMSigma / (4.0 * pi));
+ kdiss_H2II = (float) (H2IILuminosity * H2IISigma / (4.0 * pi));
for (k = 0; k < ActiveDims[2]; k++) {
for (j = 0; j < ActiveDims[1]; j++) {
radius2_yz = ddr2[1][j] + ddr2[2][k];
@@ -455,7 +457,7 @@
{
float jeans_length = 0.0;
- jeans_length = 15*kboltz*T/(4.0*M_PI*GravConst*mh*dens*density_units);
+ jeans_length = 15*kboltz*T/(4.0*pi*GravConst*mh*dens*density_units);
return sqrt(jeans_length);
}
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/Grid_AddH2DissociationFromTree.C
--- a/src/enzo/Grid_AddH2DissociationFromTree.C
+++ b/src/enzo/Grid_AddH2DissociationFromTree.C
@@ -23,6 +23,7 @@
#include "Hierarchy.h"
#include "CosmologyParameters.h"
#include "Star.h"
+#include "phys_constants.h"
#define MIN_OPENING_ANGLE 0.2 // 0.2 = arctan(11.3 deg)
@@ -38,8 +39,6 @@
int grid::AddH2DissociationFromTree(void)
{
- const double pc = 3.086e18, clight = 3e10;
-
int i, j, k, index, dim, ci;
FLOAT pos[MAX_DIMENSION];
FLOAT radius2;
@@ -77,7 +76,7 @@
H2ISigma *= (double)TimeUnits / ((double)LengthUnits * (double)LengthUnits);
// Dilution factor (prevent breaking in the rate solver near the star)
- float dilutionRadius = 10.0 * pc / (double) LengthUnits;
+ float dilutionRadius = 10.0 * pc_cm / (double) LengthUnits;
float dilRadius2 = dilutionRadius * dilutionRadius;
// Convert from #/s to RT units
@@ -87,7 +86,7 @@
/* Find sources in the tree that contribute to the cells */
SuperSourceEntry *Leaf;
- float factor = LConv_inv * H2ISigma / (4.0 * M_PI);
+ float factor = LConv_inv * H2ISigma / (4.0 * pi);
float angle;
Leaf = SourceClusteringTree;
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/Grid_AddMagneticSupernovaeToList.C
--- a/src/enzo/Grid_AddMagneticSupernovaeToList.C
+++ b/src/enzo/Grid_AddMagneticSupernovaeToList.C
@@ -56,22 +56,22 @@
// below based on the resolution of the grid at the highest refinement level
if (UseMagneticSupernovaFeedback > 1){
sn_duration = 5.0 * dtFixed;
- MagneticSupernovaDuration = sn_duration * TimeUnits / SecondsPerYear;
+ MagneticSupernovaDuration = sn_duration * TimeUnits / yr_s;
sn_radius = 1.5 * this->CellWidth[0][0];
- MagneticSupernovaRadius = sn_radius * LengthUnits / pc;
+ MagneticSupernovaRadius = sn_radius * LengthUnits / pc_cm;
}
else {
// Converting time from years to seconds, then internal units
- sn_duration = MagneticSupernovaDuration * SecondsPerYear / TimeUnits;
+ sn_duration = MagneticSupernovaDuration * yr_s / TimeUnits;
if(sn_duration < 5.0 * dtFixed) {
printf("WARNING: Magnetic supernova feedback duration is less than the recommended minimum of 5 x dtFixed\n");
- printf("Current dtFixed = %e years \n", dtFixed * TimeUnits / SecondsPerYear);
+ printf("Current dtFixed = %e years \n", dtFixed * TimeUnits / yr_s);
}
// Converting radius from parsecs to cm, then internal units
- sn_radius = MagneticSupernovaRadius * pc / LengthUnits;
+ sn_radius = MagneticSupernovaRadius * pc_cm / LengthUnits;
if(sn_radius < 1.5 * this->CellWidth[0][0]){
printf("WARNING: Magnetic supernova feedback radius is less than the recommended minimum of 1.5 x CellWidth\n");
- printf("Current CellWidth = %e pc \n", this->CellWidth[0][0] * LengthUnits / pc);
+ printf("Current CellWidth = %e pc \n", this->CellWidth[0][0] * LengthUnits / pc_cm);
}
}
@@ -93,7 +93,7 @@
random_u = (float)(mt_random()%32768)/32768.0; // random variable from 0 to 1
random_v = (float)(mt_random()%32768)/32768.0;
- random_phi = 2*M_PI*random_u; // 0 to 2pi
+ random_phi = 2*pi*random_u; // 0 to 2pi
random_theta = acos(2*random_v-1); // 0 to pi
// Setting up randomly oriented magnetic feedback of supernova
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/Grid_AddRadiationImpulse.C
--- a/src/enzo/Grid_AddRadiationImpulse.C
+++ b/src/enzo/Grid_AddRadiationImpulse.C
@@ -23,6 +23,7 @@
#include "Hierarchy.h"
#include "CosmologyParameters.h"
#include "Star.h"
+#include "phys_constants.h"
int FindField(int f, int farray[], int n);
int GetUnits(float *DensityUnits, float *LengthUnits,
@@ -33,8 +34,6 @@
FLOAT BirthTime, FLOAT* pos)
{
- const double pc = 3.086e19, clight = 3e10;
-
int i, j, k, index, dim, FieldNum;
float kdiss_r2;
FLOAT delx, dely, delz, DomainWidth[MAX_DIMENSION];
@@ -62,7 +61,7 @@
DomainWidth[dim] = DomainRightEdge[dim] - DomainLeftEdge[dim];
// Dilution factor (prevent breaking in the rate solver near the star)
- float dilutionRadius = 10.0 * pc / (double) LengthUnits;
+ float dilutionRadius = 10.0 * pc_cm / (double) LengthUnits;
float dilRadius2 = dilutionRadius * dilutionRadius;
float LightTravelDist = TimeUnits * clight / LengthUnits;
@@ -73,7 +72,7 @@
/* Loop over cells */
index = 0;
- kdiss_r2 = (float) (Luminosity * sigma / (4.0 * M_PI));
+ kdiss_r2 = (float) (Luminosity * sigma / (4.0 * pi));
for (k = 0; k < GridDimension[2]; k++) {
delz = fabs(CellLeftEdge[2][k] + 0.5*CellWidth[2][k] - pos[2]);
delz = min(delz, DomainWidth[2]-delz);
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/Grid_CRShockTubesInitializeGrid.C
--- a/src/enzo/Grid_CRShockTubesInitializeGrid.C
+++ b/src/enzo/Grid_CRShockTubesInitializeGrid.C
@@ -10,6 +10,7 @@
#include "ExternalBoundary.h"
#include "Grid.h"
#include "CosmologyParameters.h"
+#include "phys_constants.h"
int GetUnits(float *DensityUnits, float *LengthUnits,
float *TemperatureUnits, float *TimeUnits,
@@ -216,7 +217,7 @@
for (i = 0; i < GridDimension[0]; i++) {
x = CellLeftEdge[0][i] + 0.5*CellWidth[0][i];
x = x - x0; // translate gaussian to center
- BaryonField[iCRD][i] = 1.0/sqrt(4.0*3.14159*CRkappa*t0)
+ BaryonField[iCRD][i] = 1.0/sqrt(4.0*pi*CRkappa*t0)
* PEXP( -x*x/(4.0*CRkappa*t0));
} // end i for
} // end if
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/Grid_ChangeParticleTypeBeforeSN.C
--- a/src/enzo/Grid_ChangeParticleTypeBeforeSN.C
+++ b/src/enzo/Grid_ChangeParticleTypeBeforeSN.C
@@ -23,6 +23,7 @@
#include "GridList.h"
#include "ExternalBoundary.h"
#include "Grid.h"
+#include "phys_constants.h"
int GetUnits(float *DensityUnits, float *LengthUnits,
float *TemperatureUnits, float *TimeUnits,
@@ -40,7 +41,6 @@
if (Stars == NULL)
return SUCCESS;
- const float pc = 3.086e18, mh = 1.673e-24, Msun = 1.989e33;
const float PISNLowerMass = 140.0, PISNUpperMass = 260.0;
const float StartRefineAtTime = 0.99; // Percentage of stellar lifetime
const float EndRefineAtTime = 1.0;
@@ -98,7 +98,7 @@
if (factor < 1) {
Diameter = 2 * BufferZone * PopIIISupernovaRadius *
- (pc/LengthUnits);
+ (pc_cm/LengthUnits);
}
@@ -166,9 +166,6 @@
/* ASSUMES POP III PAIR-INSTABILITY SUPERNOVA FOR ENERGY */
- const float pc = 3.086e18, mh = 1.673e-24, Msun = 1.989e33;
- const float kb = 1.38e-16;
-
FLOAT StartTime, SoundCrossingTime;
float HeliumCoreMass, TransitionTime, ShockVelocity, STradius;
float BlastWaveRadius, SoundSpeed, SNTemperature;
@@ -183,14 +180,14 @@
// With a constant shock velocity, calculate time it takes for
// blastwave to travel from r=0 to r=PopIIISupernovaRadius
- ShockVelocity = 1.165 * sqrt(2.0 * SNenergy / (Mass * Msun));
- StartTime = PopIIISupernovaRadius * pc / ShockVelocity;
+ ShockVelocity = 1.165 * sqrt(2.0 * SNenergy / (Mass * SolarMass));
+ StartTime = PopIIISupernovaRadius * pc_cm / ShockVelocity;
// Because we inject thermal energy, the blastwave is delayed by a
// sound crossing time.
- SNTemperature = min(double(Mass*Msun) / double(mh) / kb, 1e8);
- SoundSpeed = sqrt(kb * SNTemperature / (0.6*mh));
- SoundCrossingTime = PopIIISupernovaRadius * pc / SoundSpeed;
+ SNTemperature = min(double(Mass*SolarMass) / double(mh) / kboltz, 1e8);
+ SoundSpeed = sqrt(kboltz * SNTemperature / (0.6*mh));
+ SoundCrossingTime = PopIIISupernovaRadius * pc_cm / SoundSpeed;
// units in cm
STradius = 3.62e19 * powf(Mass/100., 1./3.) * powf(n0, -1./3.);
@@ -212,7 +209,7 @@
// it slows by a factor of 2 near the transition to the ST phase)
if (Time < TransitionTime) {
- BlastWaveRadius = PopIIISupernovaRadius * pc + ShockVelocity * Time;
+ BlastWaveRadius = PopIIISupernovaRadius * pc_cm + ShockVelocity * Time;
printf("Free expansion\n");
} // ENDIF Free expansion phase
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/Grid_ClusterInitializeGrid.C
--- a/src/enzo/Grid_ClusterInitializeGrid.C
+++ b/src/enzo/Grid_ClusterInitializeGrid.C
@@ -21,6 +21,7 @@
#include "ExternalBoundary.h"
#include "Grid.h"
#include "CosmologyParameters.h"
+#include "phys_constants.h"
int GetUnits(float *DensityUnits, float *LengthUnits,
float *TemperatureUnits, float *TimeUnits,
@@ -95,8 +96,8 @@
/* Set various units. */
- const double Mpc = 3.0856e24, SolarMass = 1.989e33, GravConst = 6.67e-8,
- pi = 3.14159, mh = 1.67e-24, kboltz = 1.381e-16, keV=1.1604e7;
+ const double keV=1.1604e7;
+
float DensityUnits, LengthUnits, TemperatureUnits = 1, TimeUnits,
VelocityUnits, CriticalDensity = 1, BoxLength = 1, mu = 0.6;
double MassUnits = 1;
@@ -110,7 +111,7 @@
BoxLength = ComovingBoxSize*ExpansionFactor/HubbleConstantNow; // in Mpc
} else {
CriticalDensity = 2.78e11*POW(0.74,2); // in Msolar/Mpc^3 for h=0.74
- BoxLength = LengthUnits / 3.086e24;
+ BoxLength = LengthUnits / Mpc_cm;
HubbleConstantNow = 1.0;
OmegaMatterNow = 1.0;
}
@@ -161,32 +162,32 @@
x1 = NFWRadius[i]/SphereCoreRadius[sphere];
NFWDensity[i] = SphereDensity[sphere]/(x1*(1.0+x1)*(1.0+x1)); // DM Density
if (SphereType[sphere]>=6 && SphereType[sphere] <= 8) { //aka 6, 7, 8: Perseus Cluster
- rkpc=NFWRadius[i]*LengthUnits/(1.0e-3*Mpc);
+ rkpc=NFWRadius[i]*LengthUnits/(1.0e-3*Mpc_cm);
/* Initial Temperature */
if (rkpc > 300.0){
NFWTemp[i]=7.0594*1.3*POW((1.0+1.5*NFWRadius[i]/SphereRadius[sphere]),-1.6)*keV; // in K
} else{
- NFWTemp[i]=8.12e7*(1.0+POW(NFWRadius[i]*LengthUnits/(1.0e-3*Mpc*71),3))/(2.3 + pow(NFWRadius[i]*LengthUnits/(1.0e-3*Mpc*71),3));
+ NFWTemp[i]=8.12e7*(1.0+POW(NFWRadius[i]*LengthUnits/(1.0e-3*Mpc_cm*71),3))/(2.3 + pow(NFWRadius[i]*LengthUnits/(1.0e-3*Mpc_cm*71),3));
if (SphereType[sphere]==8)
NFWTemp[i]=8.12e7;
}
/* Set Gravity. NFW Dark Matter, BCG+BH */
Allg[i]=GravConst*PointSourceGravityConstant*SolarMass*
((log(1.0+x1)-x1/(1.0+x1)) /(log(1.0+1.0)-1.0/(1.0+1.0)))/POW(NFWRadius[i]*LengthUnits, 2.0) +
- ClusterSMBHBCG*POW((POW(POW(NFWRadius[i]*LengthUnits/(1.0e-3*Mpc),0.5975)/3.206e-7,0.9) +
- POW(POW(NFWRadius[i]*LengthUnits/(1.0e-3*Mpc), 1.849)/1.861e-6, 0.9)), -1.0/0.9) +
+ ClusterSMBHBCG*POW((POW(POW(NFWRadius[i]*LengthUnits/(1.0e-3*Mpc_cm),0.5975)/3.206e-7,0.9) +
+ POW(POW(NFWRadius[i]*LengthUnits/(1.0e-3*Mpc_cm), 1.849)/1.861e-6, 0.9)), -1.0/0.9) +
GravConst*SolarMass* ClusterSMBHMass / POW(NFWRadius[i]*LengthUnits, 2) ;
}//end Perseus
else{ //Elliptical galaxies
if (SphereType[sphere]==1) { //NGC 4472
- NFWTemp[i]=1.16059e7*(1.17-(1.17-0.6)*exp(-NFWRadius[i]*LengthUnits/(2.0*5*1.0e-3*Mpc)));
+ NFWTemp[i]=1.16059e7*(1.17-(1.17-0.6)*exp(-NFWRadius[i]*LengthUnits/(2.0*5*1.0e-3*Mpc_cm)));
}
if (SphereType[sphere]==2) { //NGC 6482
- NFWTemp[i]=1.16059e7*(0.4+0.4*exp(-NFWRadius[i]*LengthUnits/(2.0*8.0*1.0e-3*Mpc)));
+ NFWTemp[i]=1.16059e7*(0.4+0.4*exp(-NFWRadius[i]*LengthUnits/(2.0*8.0*1.0e-3*Mpc_cm)));
}
Allg[i]=GravConst*PointSourceGravityConstant*SolarMass*
((log(1.0+x1)-x1/(1.0+x1)) /(log(1.0+1.0)-1.0/(1.0+1.0)))/POW(NFWRadius[i]*LengthUnits, 2.0) +
- GravConst*(ClusterSMBHBCG*SolarMass*1.0e11)/POW(NFWRadius[i]*LengthUnits+EllipticalGalaxyRe*1.0e-3*Mpc/1.8153, 2) + //ClusterSMBHBCG is M_* here
+ GravConst*(ClusterSMBHBCG*SolarMass*1.0e11)/POW(NFWRadius[i]*LengthUnits+EllipticalGalaxyRe*1.0e-3*Mpc_cm/1.8153, 2) + //ClusterSMBHBCG is M_* here
GravConst*SolarMass* ClusterSMBHMass / POW(NFWRadius[i]*LengthUnits, 2) ;
}
@@ -197,7 +198,7 @@
}
dpdr_old=dpdr;
if (i>0) {
- NFWPressure[i]= (NFWPressure[i-1] + (NFWRadius[i]-NFWRadius[i-1])*16.0*Mpc*0.5*dpdr_old)/(1.0+(NFWRadius[i]-NFWRadius[i-1])*16.0*Mpc*0.5*Allg[i]/(kboltz*NFWTemp[i]/(mu * mh)));
+ NFWPressure[i]= (NFWPressure[i-1] + (NFWRadius[i]-NFWRadius[i-1])*16.0*Mpc_cm*0.5*dpdr_old)/(1.0+(NFWRadius[i]-NFWRadius[i-1])*16.0*Mpc_cm*0.5*Allg[i]/(kboltz*NFWTemp[i]/(mu * mh)));
GasDensity[i]=NFWPressure[i]/(kboltz * NFWTemp[i]/(mu * mh));
dpdr = -Allg[i]* GasDensity[i];
}
@@ -275,8 +276,8 @@
temp1 = NFWTemp[m] + (NFWTemp[m-1] - NFWTemp[m])*
(r - NFWRadius[m])/(NFWRadius[m-1] - NFWRadius[m]); // in K
if (SphereType[sphere] == 6) {
- gas_dens1 =mh*(0.0192/(1.0+POW(r*LengthUnits/(18.0e-3*Mpc),3.0))+0.046/pow((1.0+pow(r*LengthUnits/(57.0e-3*Mpc), 2.0)), 1.8)+
- 0.0048/POW((1.0+pow(r*LengthUnits/(200.0e-3*Mpc), 2.0)), 1.1))/DensityUnits/0.88;
+ gas_dens1 =mh*(0.0192/(1.0+POW(r*LengthUnits/(18.0e-3*Mpc_cm),3.0))+0.046/pow((1.0+pow(r*LengthUnits/(57.0e-3*Mpc_cm), 2.0)), 1.8)+
+ 0.0048/POW((1.0+pow(r*LengthUnits/(200.0e-3*Mpc_cm), 2.0)), 1.1))/DensityUnits/0.88;
}
else{
gas_dens1 = (GasDensity[m] + (GasDensity[m-1] - GasDensity[m])*
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/Grid_ClusterSMBHEachGridGasMass.C
--- a/src/enzo/Grid_ClusterSMBHEachGridGasMass.C
+++ b/src/enzo/Grid_ClusterSMBHEachGridGasMass.C
@@ -65,7 +65,7 @@
int dim = 0;
float DiskRadius; //ClusterSMBHDiskRadius = 0.5; ClusterSMBHDiskRadiu now a parameter.
- DiskRadius = ClusterSMBHDiskRadius*kpc/LengthUnits; //from kpc to codeunits
+ DiskRadius = ClusterSMBHDiskRadius*kpc_cm/LengthUnits; //from kpc to codeunits
for (dim = 0; dim < GridRank; dim++) {
DiskCenter[dim] = PointSourceGravityPosition[dim];
DiskLeftCorner[dim] = PointSourceGravityPosition[dim]- DiskRadius;
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/Grid_ClusterSMBHFeedback.C
--- a/src/enzo/Grid_ClusterSMBHFeedback.C
+++ b/src/enzo/Grid_ClusterSMBHFeedback.C
@@ -85,7 +85,7 @@
float JetVelocity, FastJetVelocity; // Jet Velocity in km/s (should make parameter)-- gets value from parameter ClusterSMBHJetVelocity
JetScaleRadius = ClusterSMBHJetRadius/2.0; //JetScaleRadius is half the radius of the jet launch region in cellwidths
float DiskRadius; //ClusterSMBHDiskRadius = 0.5; //ClusterSMBHDiskRadiu now a parameter
- DiskRadius = ClusterSMBHDiskRadius*kpc/LengthUnits; //from kpc to codeunits
+ DiskRadius = ClusterSMBHDiskRadius*kpc_cm/LengthUnits; //from kpc to codeunits
for (dim = 0; dim < GridRank; dim++) {
JetCenter[dim] = PointSourceGravityPosition[dim];
@@ -100,8 +100,8 @@
DiskRightCorner[dim] = PointSourceGravityPosition[dim] + DiskRadius;
}
- JetLeftCorner[jet_dim] -= ClusterSMBHJetLaunchOffset*kpc/LengthUnits;
- JetRightCorner[jet_dim] += ClusterSMBHJetLaunchOffset*kpc/LengthUnits; //from kpc to codeunits
+ JetLeftCorner[jet_dim] -= ClusterSMBHJetLaunchOffset*kpc_cm/LengthUnits;
+ JetRightCorner[jet_dim] += ClusterSMBHJetLaunchOffset*kpc_cm/LengthUnits; //from kpc to codeunits
/* Compute indices of jet launch region. */
@@ -143,7 +143,7 @@
} // end: loop over dim
/* Compute mass and momentum to be put into cells in code units if Jet is on this grid. */
- JetMdot = (ClusterSMBHJetMdot*SolarMass/3.1557e7)/(MassUnits/TimeUnits); // from M_sun/yr to code units
+ JetMdot = (ClusterSMBHJetMdot*SolarMass/yr_s)/(MassUnits/TimeUnits); // from M_sun/yr to code units
if (JetOnGrid == true){
float JetNormalization = 0.0, density_normalization, radius, Tramp;
if (jet_dim == 2){
@@ -179,9 +179,9 @@
/* Convert to code units. */
density_normalization = (JetMdot/JetNormalization)*dtFixed/POW(CellWidth[0][0], 3);
- Tramp = ClusterSMBHTramp*1.0e6*3.1557e7/TimeUnits; // from Myr to code units
+ Tramp = ClusterSMBHTramp*1.0e6*yr_s/TimeUnits; // from Myr to code units
- JetVelocity = sqrt((ClusterSMBHJetEdot*1.0e44*ClusterSMBHKineticFraction*2)/(ClusterSMBHJetMdot*SolarMass/3.1557e7))/VelocityUnits;
+ JetVelocity = sqrt((ClusterSMBHJetEdot*1.0e44*ClusterSMBHKineticFraction*2)/(ClusterSMBHJetMdot*SolarMass/yr_s))/VelocityUnits;
JetVelocity *= min((Time-ClusterSMBHStartTime)/Tramp, 1.0); //linear ramp
/* Clip edge of jet launching disk so we don't set cell off the edge of the grid. */
@@ -212,7 +212,7 @@
JetVelocity_x = JetVelocity; // mutiplied by sincos later
JetVelocity_y = JetVelocity; // mutiplied by sincos later
if (ClusterSMBHJetPrecessionPeriod > 0.00001 || ClusterSMBHJetPrecessionPeriod < -0.00001)
- ClusterSMBHJetAnglePhi = Time*2.0/(ClusterSMBHJetPrecessionPeriod*1.0e6*3.1557e7/TimeUnits); // ClusterSMBHJetPrecessionPeriod from Myr to codeunit; *2.0 instead of 2*pi because pi is used later
+ ClusterSMBHJetAnglePhi = Time*2.0/(ClusterSMBHJetPrecessionPeriod*1.0e6*yr_s/TimeUnits); // ClusterSMBHJetPrecessionPeriod from Myr to codeunit; *2.0 instead of 2*pi because pi is used later
}
else {
JetVelocity_z = JetVelocity * ClusterSMBHJetOpenAngleRadius / sqrt(POW(ClusterSMBHJetOpenAngleRadius, 2) + POW(radius, 2));
@@ -273,7 +273,7 @@
JetVelocity_z = JetVelocity; // mutiplied by sincos later
JetVelocity_y = JetVelocity; // mutiplied by sincos later
if (ClusterSMBHJetPrecessionPeriod > 0.00001)
- ClusterSMBHJetAnglePhi = Time*2.0/(ClusterSMBHJetPrecessionPeriod*1.0e6*3.1557e7/TimeUnits); // ClusterSMBHJetPrecessionPeriod from Myr to codeunit; *2.0 instead of 2*pi because pi is used later
+ ClusterSMBHJetAnglePhi = Time*2.0/(ClusterSMBHJetPrecessionPeriod*1.0e6*yr_s/TimeUnits); // ClusterSMBHJetPrecessionPeriod from Myr to codeunit; *2.0 instead of 2*pi because pi is used later
}
else {
JetVelocity_x = JetVelocity * ClusterSMBHJetOpenAngleRadius / sqrt(POW(ClusterSMBHJetOpenAngleRadius, 2) + POW(radius, 2));
@@ -335,7 +335,7 @@
JetVelocity_z = JetVelocity; // mutiplied by sincos later
JetVelocity_x = JetVelocity; // mutiplied by sincos later
if (ClusterSMBHJetPrecessionPeriod > 0.00001)
- ClusterSMBHJetAnglePhi = Time*2.0/(ClusterSMBHJetPrecessionPeriod*1.0e6*3.1557e7/TimeUnits); // ClusterSMBHJetPrecessionPeriod from Myr to codeunit; *2.0 instead of 2*pi because pi is used later
+ ClusterSMBHJetAnglePhi = Time*2.0/(ClusterSMBHJetPrecessionPeriod*1.0e6*yr_s/TimeUnits); // ClusterSMBHJetPrecessionPeriod from Myr to codeunit; *2.0 instead of 2*pi because pi is used later
}
else {
JetVelocity_y = JetVelocity * ClusterSMBHJetOpenAngleRadius / sqrt(POW(ClusterSMBHJetOpenAngleRadius, 2) + POW(radius, 2));
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/Grid_CollapseTestInitializeGrid.C
--- a/src/enzo/Grid_CollapseTestInitializeGrid.C
+++ b/src/enzo/Grid_CollapseTestInitializeGrid.C
@@ -24,6 +24,7 @@
#include "ExternalBoundary.h"
#include "Grid.h"
#include "CosmologyParameters.h"
+#include "phys_constants.h"
#define NTHETA 1000
#define NR 1000
@@ -157,8 +158,6 @@
/* Set various units. */
- const double Mpc = 3.0856e24, SolarMass = 1.989e33, GravConst = 6.67e-8,
- pi = 3.14159, mh = 1.67e-24, kboltz = 1.381e-16, LightSpeed = 2.9979e10;
float DensityUnits, LengthUnits, TemperatureUnits, TimeUnits,
VelocityUnits, CriticalDensity = 1, BoxLength = 1, mu = 0.6;
@@ -168,11 +167,11 @@
if (ComovingCoordinates) {
CosmologyComputeExpansionFactor(Time, &a, &dadt);
ExpansionFactor = a/(1.0+InitialRedshift);
- CriticalDensity = 2.78e11*pow(HubbleConstantNow, 2); // in Msolar/Mpc^3
+ CriticalDensity = 2.78e11*POW(HubbleConstantNow, 2); // in Msolar/Mpc^3
BoxLength = ComovingBoxSize*ExpansionFactor/HubbleConstantNow; // in Mpc
} else {
- CriticalDensity = 2.78e11*pow(0.74,2); // in Msolar/Mpc^3 for h=0.74
- BoxLength = LengthUnits / 3.086e24;
+ CriticalDensity = 2.78e11*POW(0.74,2); // in Msolar/Mpc^3 for h=0.74
+ BoxLength = LengthUnits / Mpc_cm;
HubbleConstantNow = 1.0;
OmegaMatterNow = 1.0;
}
@@ -193,20 +192,20 @@
NFWPressure[0] = 1.0 * kboltz * InitialTemperature / (mu * mh);
FILE *fptr = fopen("NFWProfile.out", "w");
for (i = 0; i < NFW_POINTS; i++) {
- NFWRadius[i] = SphereRadius[sphere]*pow(10, -3*(float(i)/NFW_POINTS));
+ NFWRadius[i] = SphereRadius[sphere]*POW(10, -3*(float(i)/NFW_POINTS));
x1 = NFWRadius[i]/SphereCoreRadius[sphere];
NFWDensity[i] = SphereDensity[sphere]/(x1*(1.0+x1)*(1.0+x1));
NFWMass[i] = 4.0*pi*SphereDensity[sphere]*
- (CriticalDensity/pow(ExpansionFactor, 3)) *
- pow(SphereCoreRadius[sphere]*BoxLength, 3) *
+ (CriticalDensity/POW(ExpansionFactor, 3)) *
+ POW(SphereCoreRadius[sphere]*BoxLength, 3) *
(log(1.0+x1) - x1/(x1+1.0)); // in Msolar
dpdr_old = dpdr;
dpdr = GravConst * NFWMass[i] * SolarMass *
NFWDensity[i] /
- pow(NFWRadius[i]*BoxLength*Mpc, 2);
+ POW(NFWRadius[i]*BoxLength*Mpc_cm, 2);
if (i > 0)
NFWPressure[i] = NFWPressure[i-1] -
- 0.5*(dpdr+dpdr_old)*(NFWRadius[i]-NFWRadius[i-1])*BoxLength*Mpc;
+ 0.5*(dpdr+dpdr_old)*(NFWRadius[i]-NFWRadius[i-1])*BoxLength*Mpc_cm;
NFWTemp[i] = NFWPressure[i]*mu*mh/(kboltz*NFWDensity[i]); // in K
NFWSigma[i] = sqrt(kboltz * NFWTemp[i] / (mu * mh)); // in cm/s
float mean_overdensity = 3.0*SphereDensity[sphere] / (x1*x1*x1) *
@@ -318,7 +317,7 @@
switch (SphereType[sphere]) {
case 1:
- SphereMass = (4*pi/3)*pow((SphereRadius[sphere]*LengthUnits), 3) *
+ SphereMass = (4*pi/3.0)*POW((SphereRadius[sphere]*LengthUnits), 3) *
(SphereDensity[sphere]*DensityUnits);
printf("mass = %"GSYM", lunit = %"GSYM", dunit = %"GSYM", rho = %"GSYM", r = %"GSYM"\n",
SphereMass, LengthUnits, DensityUnits, SphereDensity[sphere],
@@ -332,32 +331,32 @@
case 4:
// Integral of a Gaussian
term1 = 0.25 * sqrt(pi) * 0.3536 * // 2^(-3/2) = 0.3536
- pow((SphereCoreRadius[sphere]), 3) *
+ POW((SphereCoreRadius[sphere]), 3) *
ERF(SphereRadius[sphere]/SphereCoreRadius[sphere]);
term2 = 0.25 * SphereRadius[sphere] *
- pow(SphereCoreRadius[sphere], 2) *
- PEXP(-0.5 * pow((SphereRadius[sphere]/SphereCoreRadius[sphere]), 2));
+ POW(SphereCoreRadius[sphere], 2) *
+ PEXP(-0.5 * POW((SphereRadius[sphere]/SphereCoreRadius[sphere]), 2));
SphereMass = (4*pi*SphereDensity[sphere]*DensityUnits) *
- pow(LengthUnits, 3) * (term1 - term2);
+ POW(LengthUnits, 3) * (term1 - term2);
break;
case 5:
SphereCoreDens = (SphereDensity[sphere]*DensityUnits) *
- pow(SphereCoreRadius[sphere] / SphereRadius[sphere], -2);
- SphereCoreMass = (4*pi/3) *
- pow((SphereCoreRadius[sphere]*LengthUnits),3) * (SphereCoreDens);
+ POW(SphereCoreRadius[sphere] / SphereRadius[sphere], -2);
+ SphereCoreMass = (4*pi/3.0) *
+ POW((SphereCoreRadius[sphere]*LengthUnits),3) * (SphereCoreDens);
SphereMass = (4*pi*SphereDensity[sphere]*DensityUnits *
- pow(SphereRadius[sphere]*LengthUnits, 2) *
+ POW(SphereRadius[sphere]*LengthUnits, 2) *
((SphereRadius[sphere] - SphereCoreRadius[sphere]) *
LengthUnits)) + SphereCoreMass;
break;
case 6:
- SphereMass = pow(VelocitySound[sphere],3) /
- (sqrt(4*pi*pow((GravConst)/(4*pi),3) *
+ SphereMass = POW(VelocitySound[sphere],3) /
+ (sqrt(4*pi*POW((GravConst)/(4*pi),3) *
(SphereDensity[sphere]))) *
q(SphereCutOff[sphere]);
- SphereMass = SphereMass*DensityUnits*pow(LengthUnits,3);
+ SphereMass = SphereMass*DensityUnits*POW(LengthUnits,3);
break;
case 7:
@@ -376,7 +375,7 @@
break;
case 8:
- SphereMass = (4*pi/3)*pow((SphereRadius[sphere]*LengthUnits), 3) *
+ SphereMass = (4*pi/3.0)*POW((SphereRadius[sphere]*LengthUnits), 3) *
(SphereDensity[sphere]*DensityUnits);
ComputeRadialVelocity(SphereDensity[sphere], SphereMass,
SphereRadius[sphere], VelocityUnits,
@@ -402,15 +401,15 @@
BHMass = SphereMass / SolarMass; // in solar masses
SchwarzschildRadius = 2.0 * GravConst * SphereMass /
- (LightSpeed*LightSpeed);
- CavityRadius = 117.0 * SchwarzschildRadius * pow((BHMass/1e6), 0.08);
- InnerDensity = 4.31e-10 * pow((BHMass/1e8), -0.8) *
- pow((CavityRadius/SchwarzschildRadius/1e3), -0.6) / DensityUnits;
- InnerTemperature = 1.7e6 * pow((BHMass/1e6), -0.28);
- InnerScaleHeight = 0.46 * CavityRadius * pow(BHMass/1e6, -0.12) /
+ (clight*clight);
+ CavityRadius = 117.0 * SchwarzschildRadius * POW((BHMass/1e6), 0.08);
+ InnerDensity = 4.31e-10 * POW((BHMass/1e8), -0.8) *
+ POW((CavityRadius/SchwarzschildRadius/1e3), -0.6) / DensityUnits;
+ InnerTemperature = 1.7e6 * POW((BHMass/1e6), -0.28);
+ InnerScaleHeight = 0.46 * CavityRadius * POW(BHMass/1e6, -0.12) /
LengthUnits;
ThickenTransitionRadius = 1.9e3 * SchwarzschildRadius *
- pow((BHMass/1e6), 2.0/21) / LengthUnits;
+ POW((BHMass/1e6), 2.0/21) / LengthUnits;
CavityRadius /= LengthUnits;
printf("cgs: %"GSYM" %"GSYM" %"GSYM" %"GSYM" %"GSYM" %"GSYM"\n", SchwarzschildRadius,
CavityRadius*LengthUnits, InnerDensity*DensityUnits,
@@ -508,9 +507,9 @@
else if (xpos > 0 && ypos < 0)
theta = 2*pi + atan(ypos/xpos);
} else if (xpos == 0 && ypos > 0)
- theta = 3.14159 / 2.0;
+ theta = pi / 2.0;
else if (xpos == 0 && ypos < 0)
- theta = (3*3.14159) / 2.0;
+ theta = (3*pi) / 2.0;
else
theta = 0.0;
@@ -529,7 +528,7 @@
RotVelocity[2] = 2*pi*(xpos*sin(a)) /
SphereRotationalPeriod[sphere];
if (SphereType[sphere] == 7 || SphereType[sphere] == 9)
- GasRotVelocityCorrection = pow(rcyl/SphereRadius[sphere],-1.5);
+ GasRotVelocityCorrection = POW(rcyl/SphereRadius[sphere],-1.5);
for (dim = 0; dim < MAX_DIMENSION; dim++) {
Velocity[dim] = GasRotVelocityCorrection * RotVelocity[dim];
DMVelocity[dim] = DMRotVelocityCorrection * RotVelocity[dim];
@@ -566,10 +565,10 @@
if (SphereType[sphere] == 1)
dens1 = SphereDensity[sphere];
- /* 2) r^-2 power law */
+ /* 2) r^-2 POWer law */
if (SphereType[sphere] == 2)
- dens1 = SphereDensity[sphere]*pow(r/SphereRadius[sphere], -2);
+ dens1 = SphereDensity[sphere]*POW(r/SphereRadius[sphere], -2);
/* 3) NFW profile (use look-up table for temperature and
velocity dispersion)*/
@@ -591,17 +590,17 @@
if (SphereType[sphere] == 4) {
dens1 = SphereDensity[sphere]*
- PEXP(-0.5*pow(r/SphereCoreRadius[sphere], 2));
+ PEXP(-0.5*POW(r/SphereCoreRadius[sphere], 2));
}
- /* 5) r^-2 power law with core radius */
+ /* 5) r^-2 POWer law with core radius */
if (SphereType[sphere] == 5) {
if (r < SphereCoreRadius[sphere]) {
- dens1 = SphereDensity[sphere]*pow(SphereCoreRadius[sphere]/
+ dens1 = SphereDensity[sphere]*POW(SphereCoreRadius[sphere]/
SphereRadius[sphere], -2);
} else {
- dens1 = SphereDensity[sphere]*pow(r/SphereRadius[sphere], -2);
+ dens1 = SphereDensity[sphere]*POW(r/SphereRadius[sphere], -2);
}
}
@@ -638,14 +637,14 @@
/* 9) Circumbinary BH accretion disk (Lippai et al. 2008) */
if (SphereType[sphere] == 9) {
- MidplaneDensity = InnerDensity * pow(rcyl/CavityRadius, -0.6);
+ MidplaneDensity = InnerDensity * POW(rcyl/CavityRadius, -0.6);
MidplaneTemperature =
- InnerTemperature * pow(rcyl/CavityRadius, -0.9);
+ InnerTemperature * POW(rcyl/CavityRadius, -0.9);
if (rcyl < ThickenTransitionRadius)
ScaleHeight = InnerScaleHeight;
else
ScaleHeight = InnerScaleHeight *
- pow(rcyl/ThickenTransitionRadius, 1.05);
+ POW(rcyl/ThickenTransitionRadius, 1.05);
ScaleHeight = max(ScaleHeight, 2*CellWidth[0][0]);
// printf("r=%"FSYM", z=%"FSYM", h=%"GSYM", rho=%"GSYM", T=%"GSYM"\n",
@@ -708,7 +707,7 @@
if (dim == 0)
dens1 = SphereDensity[sphere]*PEXP(-drad/ScaleHeightR)/
- pow(cosh(zheight/max(ScaleHeightz, CellWidth[0][0])), 2);
+ POW(cosh(zheight/max(ScaleHeightz, CellWidth[0][0])), 2);
// if (dens1 < density)
// break;
@@ -723,11 +722,11 @@
float accel = PointSourceGravityConstant;
if (PointSourceGravity == 1)
accel = PointSourceGravityConstant/
- (pow(drad,3) + pow(PointSourceGravityCoreRadius, 3));
+ (POW(drad,3) + POW(PointSourceGravityCoreRadius, 3));
if (PointSourceGravity == 2) {
x1 = drad/PointSourceGravityCoreRadius;
accel = PointSourceGravityConstant*(log(1+x1)-x1/(1+x1))/
- pow(drad, 3);
+ POW(drad, 3);
}
float vel = sqrt(accel);
@@ -748,10 +747,10 @@
} // end: disk
- /* 11) stellar wind, with r^-2 power law density*/
+ /* 11) stellar wind, with r^-2 POWer law density*/
if (SphereType[sphere] == 11) {
radial_velocity = StellarWindSpeed/VelocityUnits;
- dens1 = StellarWindDensity*pow(r/StellarWindRadius, -2);
+ dens1 = StellarWindDensity*POW(r/StellarWindRadius, -2);
Velocity[0] += radial_velocity * xpos / r;
Velocity[1] += radial_velocity * ypos / r;
Velocity[2] += radial_velocity * zpos / r;
@@ -902,13 +901,13 @@
if (MultiSpecies > 1) {
BaryonField[HMNum][n] = CollapseTestInitialFractionHM *
- BaryonField[HIINum][n]* pow(temperature,float(0.88));
+ BaryonField[HIINum][n]* POW(temperature,float(0.88));
BaryonField[H2IINum][n] = CollapseTestInitialFractionH2II *
- 2.0*BaryonField[HIINum][n]* pow(temperature,float(1.8));
+ 2.0*BaryonField[HIINum][n]* POW(temperature,float(1.8));
if (ComovingCoordinates)
BaryonField[H2INum][n] = H2I_Fraction *
- BaryonField[0][n]*CoolData.HydrogenFractionByMass*pow(301.0,5.1)*
- pow(OmegaMatterNow, float(1.5))/
+ BaryonField[0][n]*CoolData.HydrogenFractionByMass*POW(301.0,5.1)*
+ POW(OmegaMatterNow, float(1.5))/
(OmegaMatterNow*BaryonMeanDensity)/
// CosmologySimulationOmegaBaryonNow/
HubbleConstantNow*2.0;
@@ -968,7 +967,7 @@
if (HydroMethod != Zeus_Hydro)
for (dim = 0; dim < GridRank; dim++)
- BaryonField[1][n] += 0.5*pow(BaryonField[ivel+dim][n], 2);
+ BaryonField[1][n] += 0.5*POW(BaryonField[ivel+dim][n], 2);
/* Set particles if being used (generate a number of particle
proportional to density). */
@@ -977,11 +976,11 @@
if (i >= GridStartIndex[0] && i <= GridEndIndex[0] &&
j >= GridStartIndex[1] && j <= GridEndIndex[1] &&
k >= GridStartIndex[2] && k <= GridEndIndex[2] ) {
- ParticleCount += density/pow(float(RefineBy), GridRank*level);
+ ParticleCount += density/POW(float(RefineBy), GridRank*level);
while (ParticleCount > 1) {
if (SetupLoopCount > 0) {
ParticleMass[npart] = ParticleMeanDensity*
- pow(float(RefineBy), GridRank*level);
+ POW(float(RefineBy), GridRank*level);
ParticleNumber[npart] = CollapseTestParticleCount++;
ParticleType[npart] = PARTICLE_TYPE_DARK_MATTER;
@@ -1065,8 +1064,8 @@
double BE(double r)
{
double factor;
- factor = 4.8089e-04*pow(r,5) - 1.0173e-02*pow(r,4) + 7.7899e-02*pow(r,3) -
- 2.3299e-01*pow(r,2) + 1.4721e-02*r + 1.0008e+00;
+ factor = 4.8089e-04*POW(r,5) - 1.0173e-02*POW(r,4) + 7.7899e-02*POW(r,3) -
+ 2.3299e-01*POW(r,2) + 1.4721e-02*r + 1.0008e+00;
return factor;
}
@@ -1075,8 +1074,8 @@
double q(double r)
{
double factor;
- factor = 0.0015970*pow(r,5) - 0.0229113*pow(r,4) + 0.0386709*pow(r,3) +
- 0.7350457*pow(r,2) - 0.5490283*r + 0.0872061;
+ factor = 0.0015970*POW(r,5) - 0.0229113*POW(r,4) + 0.0386709*POW(r,3) +
+ 0.7350457*POW(r,2) - 0.5490283*r + 0.0872061;
return factor;
}
@@ -1092,13 +1091,10 @@
double Maxwellian(double c_tilda, double vel_unit, double mu, double gamma)
{
// Set constants
- double mh = 1.67e-24;
- double kboltz = 1.38e-16;
- double pi = 3.14159;
// Compute temperature in cgs units
double c = c_tilda*vel_unit;
- double T = (pow(c,2)*mu*mh)/(kboltz*gamma);
+ double T = (POW(c,2)*mu*mh)/(kboltz*gamma);
// Compute random Maxwellian velocity
double mean = 0;
@@ -1114,8 +1110,8 @@
double ERF(double x)
{
- return (2.0 / sqrt(M_PI)) *
- (x - pow(x,3)/3.0 + pow(x,5)/10.0 - pow(x,7)/42.0 + pow(x,9)/216.0);
+ return (2.0 / sqrt(pi)) *
+ (x - POW(x,3)/3.0 + POW(x,5)/10.0 - POW(x,7)/42.0 + POW(x,9)/216.0);
}
int ComputeRadialVelocity(float density, double mass, float r_init,
@@ -1124,8 +1120,7 @@
double exterior_rho[], int Npts)
{
- const float theta0 = 0.5*M_PI, theta1 = 1.9*M_PI;
- const float Grav = 6.673e-8, Mpc = 3.086e24, yr = 3.1557e7, Msun = 1.989e33;
+ const float theta0 = 0.5*pi, theta1 = 1.9*pi;
const float delta_i = 0.5;
float dtheta, delta_t[NTHETA], Theta[NTHETA], Lambda[NTHETA], Chi;
float beta[NTHETA], little_d[NTHETA];
@@ -1138,44 +1133,44 @@
dtheta = (theta1 - theta0) / (NTHETA-1.0);
for (i = 0; i < NTHETA; i++) {
Theta[i] = theta0 + i*dtheta;
- Lambda[i] = pow(sin(Theta[i]/2), 2) *
- pow((Theta[i] - sin(Theta[i])) / M_PI, -8.0/9);
- delta_t[i] = 4.5 * pow((Theta[i] - sin(Theta[i])), 2) /
- pow((1.0 - cos(Theta[i])), 3) - 1.0;
+ Lambda[i] = POW(sin(Theta[i]/2), 2) *
+ POW((Theta[i] - sin(Theta[i])) / pi, -8.0/9);
+ delta_t[i] = 4.5 * POW((Theta[i] - sin(Theta[i])), 2) /
+ POW((1.0 - cos(Theta[i])), 3) - 1.0;
beta[i] = sin(0.5*Theta[i]) * sin(0.5*Theta[i]);
little_d[i] = 0.75 * (Theta[i] - sin(Theta[i]));
if (delta_t[i] < density)
iinit = i;
- if (Theta[i] < M_PI) // Turnaround
+ if (Theta[i] < pi) // Turnaround
iturn = i;
}
- t_init = 5.38e8 * yr * pow((1+InitialRedshift) / 10.0, -1.5);
- z_vir = pow((Theta[iinit] - sin(Theta[iinit])) / (2*M_PI), 2.0/3) *
+ t_init = 5.38e8 * yr_s * POW((1+InitialRedshift) / 10.0, -1.5);
+ z_vir = POW((Theta[iinit] - sin(Theta[iinit])) / (2*pi), 2.0/3.0) *
(1 + InitialRedshift) - 1.0;
- t_vir = 5.38e8 * yr * pow((1+z_vir) / 10.0, -1.5);
+ t_vir = 5.38e8 * yr_s * POW((1+z_vir) / 10.0, -1.5);
t_ta = 0.5 * t_vir;
- t_i = t_ta * pow(delta_i, 1.5) / (3*M_PI/4);
- rho_ci = 1.0 / (6 * M_PI * Grav * t_i*t_i);
- r_i = pow(mass / (4*M_PI/3 * rho_ci), 1.0/3);
+ t_i = t_ta * POW(delta_i, 1.5) / (3*pi/4);
+ rho_ci = 1.0 / (6 * pi * GravConst * t_i*t_i);
+ r_i = POW(mass / (4*pi/3.0 * rho_ci), 1.0/3.0);
- dm_rotvel_corr = 2 * pow(2*r_init/r_i, 2);
+ dm_rotvel_corr = 2 * POW(2*r_init/r_i, 2);
// dm_rotvel_corr = 1;
- r_ta_now = pow((3*M_PI/4), -8.0/9) * pow(delta_i, 1.0/3) * r_i
- * pow(t_init / t_i, 8.0/9);
+ r_ta_now = POW((3*pi/4), -8.0/9) * POW(delta_i, 1.0/3.0) * r_i
+ * POW(t_init / t_i, 8.0/9);
for (i = 0; i < NTHETA; i++) {
local_Vr[NTHETA-i-1] = (Lambda[i] * r_ta_now / t_init) *
sin(Theta[i]) * (Theta[i] - sin(Theta[i])) /
- pow(1.0 - cos(Theta[i]), 2);
+ POW(1.0 - cos(Theta[i]), 2);
Vl = local_Vr[NTHETA-i-1] / (r_ta_now / t_init);
// Subtract Hubble flow and convert to code units
- local_Vr[NTHETA-i-1] = (Vl - (2.0/3)*Lambda[i]) * r_ta_now / t_init;
+ local_Vr[NTHETA-i-1] = (Vl - (2.0/3.0)*Lambda[i]) * r_ta_now / t_init;
local_Vr[NTHETA-i-1] /= VelocityUnits;
Chi = 1.0 - 1.5 * Vl / Lambda[i];
- local_rho[NTHETA-i-1] = pow(little_d[i],2) / pow(beta[i],3) /
+ local_rho[NTHETA-i-1] = POW(little_d[i],2) / POW(beta[i],3) /
(1.0 + 3*Chi);
local_radius[NTHETA-i-1] = Lambda[i] * r_ta_now / LengthUnits;
if (local_radius[NTHETA-i-1] > r_init)
@@ -1211,7 +1206,7 @@
exterior_rho[i] = density;
}
// printf("%"ISYM": r = %"GSYM" pc, v_r = %"GSYM" km/s, rho = %"GSYM"\n",
-// i, radius_vr[i]*LengthUnits/3.086e18, Vr[i]*VelocityUnits/1e5,
+// i, radius_vr[i]*LengthUnits/pc, Vr[i]*VelocityUnits/1e5,
// exterior_rho[i]);
}
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/Grid_ComputeAccelerationFieldExternal.C
--- a/src/enzo/Grid_ComputeAccelerationFieldExternal.C
+++ b/src/enzo/Grid_ComputeAccelerationFieldExternal.C
@@ -185,8 +185,8 @@
((log(1.0+x )-x /(1.0+x )) /
(log(1.0+1.0)-1.0/(1.0+1.0))) /
POW(radius*LengthUnits, 2.0) / AccelUnits +
- ClusterSMBHBCG*POW((POW(POW(radius*LengthUnits/(1.0e-3*Mpc),0.5975)/3.206e-7,0.9) +
- POW(POW(radius*LengthUnits/(1.0e-3*Mpc), 1.849)/1.861e-6, 0.9)), -1.0/0.9) / AccelUnits +
+ ClusterSMBHBCG*POW((POW(POW(radius*LengthUnits/(1.0e-3*Mpc_cm),0.5975)/3.206e-7,0.9) +
+ POW(POW(radius*LengthUnits/(1.0e-3*Mpc_cm), 1.849)/1.861e-6, 0.9)), -1.0/0.9) / AccelUnits +
GravConst*SolarMass*ClusterSMBHMass / POW(radius*LengthUnits, 2) / AccelUnits; // + BCG + BH mass
/*Bondi*/
if(ClusterSMBHCalculateGasMass == 4){
@@ -194,8 +194,8 @@
((log(1.0+x )-x /(1.0+x )) /
(log(1.0+1.0)-1.0/(1.0+1.0))) /
POW(radius*LengthUnits, 2.0) / AccelUnits +
- ClusterSMBHBCG*POW((POW(POW(radius*LengthUnits/(1.0e-3*Mpc),0.5975)/3.206e-7,0.9) +
- POW(POW(radius*LengthUnits/(1.0e-3*Mpc), 1.849)/1.861e-6, 0.9)), -1.0/0.9) / AccelUnits +
+ ClusterSMBHBCG*POW((POW(POW(radius*LengthUnits/(1.0e-3*Mpc_cm),0.5975)/3.206e-7,0.9) +
+ POW(POW(radius*LengthUnits/(1.0e-3*Mpc_cm), 1.849)/1.861e-6, 0.9)), -1.0/0.9) / AccelUnits +
GravConst*SolarMass*ClusterSMBHMass/POW(radius*LengthUnits - 2.0*GravConst*SolarMass/POW(clight,2), 2)/ AccelUnits;
}
/*Elliptical Galaxy Fixed Gravity*/
@@ -204,7 +204,7 @@
((log(1.0+x )-x /(1.0+x )) /
(log(1.0+1.0)-1.0/(1.0+1.0))) /
POW(radius*LengthUnits, 2.0) / AccelUnits +
- GravConst*(ClusterSMBHBCG*SolarMass*1.0e11)/POW(radius*LengthUnits+EllipticalGalaxyRe*1.0e-3*Mpc/1.8153, 2)/AccelUnits +
+ GravConst*(ClusterSMBHBCG*SolarMass*1.0e11)/POW(radius*LengthUnits+EllipticalGalaxyRe*1.0e-3*Mpc_cm/1.8153, 2)/AccelUnits +
GravConst*SolarMass*ClusterSMBHMass/POW(radius*LengthUnits - 2.0*GravConst*SolarMass/POW(clight,2), 2)/ AccelUnits;
}
}
@@ -304,8 +304,8 @@
((log(1.0+x )-x /(1.0+x )) /
(log(1.0+1.0)-1.0/(1.0+1.0))) /
POW(radius*LengthUnits, 2.0) / AccelUnits +
- ClusterSMBHBCG*POW((POW(POW(radius*LengthUnits/(1.0e-3*Mpc),0.5975)/3.206e-7,0.9) +
- POW(POW(radius*LengthUnits/(1.0e-3*Mpc), 1.849)/1.861e-6, 0.9)), -1.0/0.9) / AccelUnits +
+ ClusterSMBHBCG*POW((POW(POW(radius*LengthUnits/(1.0e-3*Mpc_cm),0.5975)/3.206e-7,0.9) +
+ POW(POW(radius*LengthUnits/(1.0e-3*Mpc_cm), 1.849)/1.861e-6, 0.9)), -1.0/0.9) / AccelUnits +
GravConst*SolarMass*ClusterSMBHMass / POW(radius*LengthUnits, 2) / AccelUnits; // + BCG + BH mass
/*Bondi*/
if(ClusterSMBHCalculateGasMass == 4){
@@ -313,8 +313,8 @@
((log(1.0+x )-x /(1.0+x )) /
(log(1.0+1.0)-1.0/(1.0+1.0))) /
POW(radius*LengthUnits, 2.0) / AccelUnits +
- ClusterSMBHBCG*POW((POW(POW(radius*LengthUnits/(1.0e-3*Mpc),0.5975)/3.206e-7,0.9) +
- POW(POW(radius*LengthUnits/(1.0e-3*Mpc), 1.849)/1.861e-6, 0.9)), -1.0/0.9) / AccelUnits +
+ ClusterSMBHBCG*POW((POW(POW(radius*LengthUnits/(1.0e-3*Mpc_cm),0.5975)/3.206e-7,0.9) +
+ POW(POW(radius*LengthUnits/(1.0e-3*Mpc_cm), 1.849)/1.861e-6, 0.9)), -1.0/0.9) / AccelUnits +
GravConst*SolarMass*ClusterSMBHMass/POW(radius*LengthUnits - 2.0*GravConst*SolarMass/POW(clight,2), 2)/ AccelUnits;
}
/*Elliptical Galaxy Fixed Gravity*/
@@ -323,7 +323,7 @@
((log(1.0+x )-x /(1.0+x )) /
(log(1.0+1.0)-1.0/(1.0+1.0))) /
POW(radius*LengthUnits, 2.0) / AccelUnits +
- GravConst*(ClusterSMBHBCG*SolarMass*1.0e11)/POW(radius*LengthUnits+EllipticalGalaxyRe*1.0e-3*Mpc/1.8153, 2)/AccelUnits +
+ GravConst*(ClusterSMBHBCG*SolarMass*1.0e11)/POW(radius*LengthUnits+EllipticalGalaxyRe*1.0e-3*Mpc_cm/1.8153, 2)/AccelUnits +
GravConst*SolarMass*ClusterSMBHMass/POW(radius*LengthUnits - 2.0*GravConst*SolarMass/POW(clight,2), 2)/ AccelUnits;
}
}
@@ -427,21 +427,21 @@
rcyl = sqrt(xpos1*xpos1 + ypos1*ypos1 + zpos1*zpos1);
radius = radius*LengthUnits;
rcyl = rcyl*LengthUnits;
- accelsph = (GravConst)*MBulge*SolarMass/POW(radius+rBulge*Mpc,2)
- + pi*GravConst*densDMConst*POW(rDMConst*Mpc,3)/POW(radius,2)
- *(-2.0*atan(radius/rDMConst/Mpc)
- +2.0*log(1.0+radius/rDMConst/Mpc)
- +log(1.0+POW(radius/rDMConst/Mpc,2))
+ accelsph = (GravConst)*MBulge*SolarMass/POW(radius+rBulge*Mpc_cm,2)
+ + pi*GravConst*densDMConst*POW(rDMConst*Mpc_cm,3)/POW(radius,2)
+ *(-2.0*atan(radius/rDMConst/Mpc_cm)
+ +2.0*log(1.0+radius/rDMConst/Mpc_cm)
+ +log(1.0+POW(radius/rDMConst/Mpc_cm,2))
);
accelcylR = GravConst*MSDisk*SolarMass*rcyl/sqrt(POW(POW(rcyl,2)
- + POW(SDiskScaleHeightR*Mpc+sqrt(POW(zheight*LengthUnits,2)
- + POW(SDiskScaleHeightz*Mpc,2)),2),3));
+ + POW(SDiskScaleHeightR*Mpc_cm+sqrt(POW(zheight*LengthUnits,2)
+ + POW(SDiskScaleHeightz*Mpc_cm,2)),2),3));
accelcylz = GravConst*MSDisk*SolarMass/sqrt(POW(zheight*LengthUnits,2)
- + POW(SDiskScaleHeightz*Mpc,2))*zheight*LengthUnits/sqrt(POW(POW(rcyl,2)
- + POW(SDiskScaleHeightR*Mpc+sqrt(POW(zheight*LengthUnits,2)
- + POW(SDiskScaleHeightz*Mpc,2)),2),3))
- *( SDiskScaleHeightR*Mpc+sqrt(POW(zheight*LengthUnits,2)
- + POW(SDiskScaleHeightz*Mpc,2))
+ + POW(SDiskScaleHeightz*Mpc_cm,2))*zheight*LengthUnits/sqrt(POW(POW(rcyl,2)
+ + POW(SDiskScaleHeightR*Mpc_cm+sqrt(POW(zheight*LengthUnits,2)
+ + POW(SDiskScaleHeightz*Mpc_cm,2)),2),3))
+ *( SDiskScaleHeightR*Mpc_cm+sqrt(POW(zheight*LengthUnits,2)
+ + POW(SDiskScaleHeightz*Mpc_cm,2))
)/AccelUnits;
accelsph = (radius ==0.0?0.0:fabs(accelsph )/(radius/LengthUnits)/AccelUnits);
@@ -500,21 +500,21 @@
radius = radius*LengthUnits;
rcyl = rcyl*LengthUnits;
- accelsph = (GravConst)*MBulge*SolarMass/POW(radius+rBulge*Mpc,2)
- + pi*GravConst*densDMConst*POW(rDMConst*Mpc,3)/POW(radius,2)
- * (-2.0*atan(radius/rDMConst/Mpc)
- +2.0*log(1.0+radius/rDMConst/Mpc)
- +log(1.0+POW(radius/rDMConst/Mpc,2))
+ accelsph = (GravConst)*MBulge*SolarMass/POW(radius+rBulge*Mpc_cm,2)
+ + pi*GravConst*densDMConst*POW(rDMConst*Mpc_cm,3)/POW(radius,2)
+ * (-2.0*atan(radius/rDMConst/Mpc_cm)
+ +2.0*log(1.0+radius/rDMConst/Mpc_cm)
+ +log(1.0+POW(radius/rDMConst/Mpc_cm,2))
);
accelcylR = GravConst*MSDisk*SolarMass*rcyl/sqrt(POW(POW(rcyl,2)
- + POW(SDiskScaleHeightR*Mpc+sqrt(POW(zheight*LengthUnits,2)
- + POW(SDiskScaleHeightz*Mpc,2)),2),3));
+ + POW(SDiskScaleHeightR*Mpc_cm+sqrt(POW(zheight*LengthUnits,2)
+ + POW(SDiskScaleHeightz*Mpc_cm,2)),2),3));
accelcylz = GravConst*MSDisk*SolarMass/sqrt(POW(zheight*LengthUnits,2)
- + POW(SDiskScaleHeightz*Mpc,2))*zheight*LengthUnits/sqrt(POW(POW(rcyl,2)
- + POW(SDiskScaleHeightR*Mpc+sqrt(POW(zheight*LengthUnits,2)
- + POW(SDiskScaleHeightz*Mpc,2)),2),3))
- *( SDiskScaleHeightR*Mpc+sqrt(POW(zheight*LengthUnits,2)
- + POW(SDiskScaleHeightz*Mpc,2))
+ + POW(SDiskScaleHeightz*Mpc_cm,2))*zheight*LengthUnits/sqrt(POW(POW(rcyl,2)
+ + POW(SDiskScaleHeightR*Mpc_cm+sqrt(POW(zheight*LengthUnits,2)
+ + POW(SDiskScaleHeightz*Mpc_cm,2)),2),3))
+ *( SDiskScaleHeightR*Mpc_cm+sqrt(POW(zheight*LengthUnits,2)
+ + POW(SDiskScaleHeightz*Mpc_cm,2))
)/AccelUnits;
accelsph = (radius ==0.0?0.0:fabs(accelsph )/(radius/LengthUnits)/AccelUnits);
@@ -561,16 +561,15 @@
FLOAT xc = 0.5, yc = 0.5, zc = 0.5;
double rs = rvir / c;
- double Mvir = 4.0*M_PI*rhoc*POW(rs,3)*(log(1.0+c)-c/(1.0+c));
+ double Mvir = 4.0*pi*rhoc*POW(rs,3)*(log(1.0+c)-c/(1.0+c));
float DensityUnits = 1.0, LengthUnits = 1.0, TemperatureUnits = 1,
TimeUnits = 1.0, VelocityUnits = 1.0;
GetUnits(&DensityUnits, &LengthUnits, &TemperatureUnits,
&TimeUnits, &VelocityUnits, &MassUnits, Time);
double AccelerationUnits = LengthUnits / POW(TimeUnits,2);
- double CGSGravConst = 6.672e-8;
- printf("rhoc=%g, rvir=%g, Mvir=%g\n", rhoc, rvir, Mvir/1.989e33);
+ printf("rhoc=%g, rvir=%g, Mvir=%g\n", rhoc, rvir, Mvir/SolarMass);
FLOAT x, y, z, xpos, ypos, zpos, r;
int n = 0;
@@ -603,12 +602,12 @@
if (r < rvir/LengthUnits) {
x1 = r*LengthUnits/rs;
- M = 4.0*M_PI*rhoc*POW(rs,3)*(log(1.0+x1)-x1/(1.0+x1));
+ M = 4.0*pi*rhoc*POW(rs,3)*(log(1.0+x1)-x1/(1.0+x1));
}
else {
M = Mvir;
}
- g = CGSGravConst*M/POW(r*LengthUnits,2);
+ g = GravConst*M/POW(r*LengthUnits,2);
g /= AccelerationUnits;
if (dim == 0) {
AccelerationField[0][n] += -g*xpos/r;
@@ -636,12 +635,12 @@
if (r < rvir/LengthUnits) {
x1 = r*LengthUnits/rs;
- M = 4.0*M_PI*rhoc*POW(rs,3)*(log(1.0+x1)-x1/(1.0+x1));
+ M = 4.0*pi*rhoc*POW(rs,3)*(log(1.0+x1)-x1/(1.0+x1));
}
else {
M = Mvir;
}
- g = CGSGravConst*M/POW(r*LengthUnits,2);
+ g = GravConst*M/POW(r*LengthUnits,2);
g /= AccelerationUnits;
ParticleAcceleration[0][i] += -g*xpos/r;
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/Grid_ComputeCoolingTime.C
--- a/src/enzo/Grid_ComputeCoolingTime.C
+++ b/src/enzo/Grid_ComputeCoolingTime.C
@@ -25,6 +25,7 @@
#include "Grid.h"
#include "fortran.def"
#include "CosmologyParameters.h"
+#include "phys_constants.h"
/* This parameter controls whether the cooling function recomputes
the metal cooling rates. It is reset by RadiationFieldUpdate. */
@@ -300,8 +301,7 @@
#ifdef TRANSFER
/* unit conversion from Enzo RT units to CGS */
- const float ev2erg = 1.60217653E-12;
- float rtunits = ev2erg / TimeUnits;
+ float rtunits = erg_eV / TimeUnits;
if ( RadiativeTransfer ){
my_fields.RT_HI_ionization_rate = BaryonField[kphHINum];
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/Grid_ComputeElementalDensity.C
--- a/src/enzo/Grid_ComputeElementalDensity.C
+++ b/src/enzo/Grid_ComputeElementalDensity.C
@@ -37,6 +37,7 @@
#include "ExternalBoundary.h"
#include "Grid.h"
#include "CosmologyParameters.h"
+#include "phys_constants.h"
/* function prototypes */
@@ -316,9 +317,9 @@
}
/* Loop over grid and compute emissivity
- (fh is hydrogen fraction by mass; 1.67e-24 is proton mass). */
+ (fh is hydrogen fraction by mass; mh is proton mass). */
- float fh = 0.76, ConvertToNumberDensity = DensityUnits/1.67e-24;
+ float fh = 0.76, ConvertToNumberDensity = DensityUnits/mh;
float nElement, nH, LogFraction, Metallicity, Fraction,
deld, delt, logd, logt, dd, dt;
int ilogd, ilogt;
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/Grid_ComputePhotonTimestep.C
--- a/src/enzo/Grid_ComputePhotonTimestep.C
+++ b/src/enzo/Grid_ComputePhotonTimestep.C
@@ -211,7 +211,7 @@
if (G_TotalNumberOfStars > 0 && minStarLifetime < 1e6) {
dtStar = minStarLifetime/NumberOfStepsInLifetime;
} else {
- dtStar = 3.1557e13*mindtNOstars/TimeUnits;
+ dtStar = Myr_s*mindtNOstars/TimeUnits;
}
}
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/Grid_ComputePhotonTimestepHII.C
--- a/src/enzo/Grid_ComputePhotonTimestepHII.C
+++ b/src/enzo/Grid_ComputePhotonTimestepHII.C
@@ -26,6 +26,7 @@
#include "ExternalBoundary.h"
#include "Grid.h"
#include "RadiativeTransferParameters.h"
+#include "phys_constants.h"
#define MAX_CHANGE 0.1
@@ -47,7 +48,6 @@
int i, j, k, dim, index, indexn, size;
float dt, this_dt, sigma_dx, *temperature;
- const double sigmaHI = 6.345e-18, mh = 1.673e-24;
dt = huge_number;
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/Grid_ComputePhotonTimestepTau.C
--- a/src/enzo/Grid_ComputePhotonTimestepTau.C
+++ b/src/enzo/Grid_ComputePhotonTimestepTau.C
@@ -26,6 +26,7 @@
#include "ExternalBoundary.h"
#include "Grid.h"
#include "RadiativeTransferParameters.h"
+#include "phys_constants.h"
#define MAX_CHANGE 0.5
@@ -45,7 +46,6 @@
int i, j, k, dim, index, size;
float dt, this_dt, sigma_dx, *temperature;
- const double sigmaHI = 6.345e-18, mh = 1.673e-24;
dt = huge_number;
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/Grid_ComputeXrayEmissivity.C
--- a/src/enzo/Grid_ComputeXrayEmissivity.C
+++ b/src/enzo/Grid_ComputeXrayEmissivity.C
@@ -30,6 +30,7 @@
#include "fortran.def"
#include "Grid.h"
#include "CosmologyParameters.h"
+#include "phys_constants.h"
/* function prototypes */
@@ -165,7 +166,7 @@
(fh is hydrogen fraction by mass). */
float temp, ne, nH, frac, fh = 0.76;
- float ConvertToNumberDensity = DensityUnits/1.67e-24;
+ float ConvertToNumberDensity = DensityUnits/mh;
float deltemp = (temp2-temp1)/float(NumberOfTemperatureBins-1);
for (i = 0; i < size; i++) {
diff -r 47eb089a647e -r 7d8fe75451e0 src/enzo/Grid_ConductionBubbleInitialize.C
--- a/src/enzo/Grid_ConductionBubbleInitialize.C
+++ b/src/enzo/Grid_ConductionBubbleInitialize.C
@@ -22,6 +22,7 @@
#include "GridList.h"
#include "ExternalBoundary.h"
#include "Grid.h"
+#include "phys_constants.h"
// Function prototypes
int GetUnits (float *DensityUnits, float *LengthUnits,
@@ -43,7 +44,6 @@
static int ncells;
#define KEV_KELVIN 1.1604e+7
-#define KPC_CGS 3.0857e+21
#define DEFAULT_MU 0.6 // we assume total ionization
// Grid Initializer: all input values are in Enzo internal units _except_
@@ -127,7 +127,7 @@
nofr = new double[ncells];
Tofr = new double[ncells];
- dKdr_cgs = dKdr * KEV_KELVIN / KPC_CGS;
+ dKdr_cgs = dKdr * KEV_KELVIN / kpc_cm;
r_mid_cgs = r_mid * LengthUnits;
r_max_cgs = 2.0*r_mid_cgs;
@@ -147,7 +147,7 @@
// convert 1D arrays into Enzo internal units.
for(i=0; i<ncells; i++){
rad[i] /= LengthUnits; // convert to enzo distance
- nofr[i] *= DEFAULT_MU * 1.67e-24 / DensityUnits; // convert to enzo-unit density (from number density)
+ nofr[i] *= DEFAULT_MU * mh / DensityUnits; // convert to enzo-unit density (from number density)
Tofr[i] /= (TemperatureUnits*(Gamma-1.0)*DEFAULT_MU); // convert from temp to internal energy
}
@@ -331,9 +331,9 @@
- // g*mu*mp/kb;
+ // g*mu*mp/kboltz;
- bunch_of_constants = g*DEFAULT_MU*(1.67e-24)/(1.38e-16);
+ bunch_of_constants = g*DEFAULT_MU*(mh)/(kboltz);
/*
printf("n_mid, bunch_of_constants: %e %e %e\n",n_mid,bunch_of_constants,g);
This diff is so big that we needed to truncate the remainder.
Repository URL: https://bitbucket.org/enzo/enzo-dev/
--
This is a commit notification from bitbucket.org. You are receiving
this because you have the service enabled, addressing the recipient of
this email.
Reply all
Reply to author
Forward
0 new messages