changa 3.5
Loading...
Searching...
No Matches
cooling_metal_H2.h
1
2#ifndef COOLING_METAL_HINCLUDED
3#define COOLING_METAL_HINCLUDED
4
5/*
6 * Cooling code for cosmology simulations.
7 * Originally written by James Wadsley and Sijing Shen, McMaster
8 * University for GASOLINE.
9 */
10
11/* Global consts */
12/*#if defined(COOLDEBUG)
13#include "mdl.h"
14#endif*/
15#include "param.h"
16#include "rpc/xdr.h"
17
18#ifdef __cplusplus
19extern "C" {
20#endif
21
22#include "stiff.h"
23
24/* Constants */
25#define CL_B_gm (6.022e23*(938.7830/931.494))/*Avegadro's Number * Mass_Hydrogen/Energy_AMU */
26#define CL_k_Boltzmann 1.38066e-16
27#define CL_eV_erg 1.60219e-12
28#define CL_eV_per_K (CL_k_Boltzmann/CL_eV_erg)
29/*
30#define CL_RT_FLOAT float
31#define CL_RT_MIN 1e-38
32#define CL_RT_MIN FLT_MIN
33*/
34
35#define CL_RT_FLOAT double
36#define CL_RT_MIN 1e-100
37
38/*
39#define CL_RT_MIN DBL_MIN
40*/
41/*
42 * Work around for Dec ev6 flawed
43 * treatment of sub-normal numbers
44 */
45#define CL_MAX_NEG_EXP_ARG -500.
46
47#define CL_NMAXBYTETABLE 56000
48#define MU_METAL 17.6003
49#define ZSOLAR 0.0130215
50
51typedef struct CoolingParametersStruct {
52 int bIonNonEqm;
53 int nCoolingTable;
54 int bUV;
55 int bMetal;
56 char *CoolInFile;
57 int bUVTableUsesTime;
58 int bDoIonOutput;
59 int bLowTCool;
60 int bSelfShield;
61 int bShieldHI; /* Set to true if dust shields HI;*/
62 double dMassFracHelium;
63 double dCoolingTmin;
64 double dCoolingTmax;
65 double dClump;
66 double dLymanWernerFrac; /* Fraction of Lyman Werner radiation that escapes birth cloud. 0.5 is a good value.*/
67#ifdef SHIELDSF
68 double dMaxTShield; // Temperature ceiling for dust shielding.
69#endif
70} COOLPARAM;
71
72typedef struct CoolingParticleStruct {
73 double f_HI,f_HeI,f_HeII;
74 double f_H2; /* Abundance of ions */
75 double dLymanWerner; /* Flux of Lyman Werner radiation at the gas particle */
76} COOLPARTICLE;
77
78typedef struct {
79 double e,Total;
80 double HI,HII,HeI,HeII,HeIII;
81 double H2;
82} PERBARYON;
83
84typedef struct {
85 double zTime;
86
87 double Rate_Phot_HI;
88 double Rate_Phot_HeI;
89 double Rate_Phot_HeII;
90 double Rate_Phot_H2_cosmo; /* Dissociating radiation from the cosmic background for H2*/
91
92 double Heat_Phot_HI;
93 double Heat_Phot_HeI;
94 double Heat_Phot_HeII;
95 double Heat_Phot_H2;
96} UVSPECTRUM;
97
98typedef struct {
99 double Rate_Phot_HI;
100 double Rate_Phot_HeI;
101 double Rate_Phot_HeII;
102 double Rate_Phot_H2_cosmo;
103
104 double Heat_Phot_HI;
105 double Heat_Phot_HeI;
106 double Heat_Phot_HeII;
107 double Heat_Phot_H2;
108
109 double Cool_Coll_HI;
110 double Cool_Coll_HeI;
111 double Cool_Coll_HeII;
112 double Cool_Diel_HeII;
113 double Cool_Coll_H2;
114
115 double Cool_Comp;
116 double Tcmb;
117 double Cool_LowTFactor;
118
119} RATES_NO_T;
120
121typedef struct {
122 CL_RT_FLOAT Rate_Coll_HI;
123 CL_RT_FLOAT Rate_Coll_HeI;
124 CL_RT_FLOAT Rate_Coll_HeII;
125 CL_RT_FLOAT Rate_Coll_e_H2;
126 CL_RT_FLOAT Rate_Coll_HI_H2;
127 CL_RT_FLOAT Rate_Coll_H2_H2;
128 CL_RT_FLOAT Rate_Coll_Hm_e; /*gas phase formation of H2 */
129 CL_RT_FLOAT Rate_Coll_HI_e; /*--------------------*/
130 CL_RT_FLOAT Rate_Coll_HII_H2; /*--------------------*/
131 CL_RT_FLOAT Rate_Coll_Hm_HII; /*-------------------- */
132 CL_RT_FLOAT Rate_HI_e; /*-------------------- */
133 CL_RT_FLOAT Rate_HI_Hm; /*gas phase formation of H2 */
134 CL_RT_FLOAT Rate_Radr_HII;
135 CL_RT_FLOAT Rate_Radr_HeII;
136 CL_RT_FLOAT Rate_Radr_HeIII;
137 CL_RT_FLOAT Rate_Diel_HeII;
138 CL_RT_FLOAT Rate_Chtr_HeII;
139
140 CL_RT_FLOAT Cool_Brem_1;
141 CL_RT_FLOAT Cool_Brem_2;
142 CL_RT_FLOAT Cool_Radr_HII;
143 CL_RT_FLOAT Cool_Radr_HeII;
144 CL_RT_FLOAT Cool_Radr_HeIII;
145 CL_RT_FLOAT Cool_Line_HI;
146 CL_RT_FLOAT Cool_Line_HeI;
147 CL_RT_FLOAT Cool_Line_HeII;
148 CL_RT_FLOAT Cool_Line_H2_H;
149 CL_RT_FLOAT Cool_Line_H2_H2;
150 CL_RT_FLOAT Cool_Line_H2_He;
151 CL_RT_FLOAT Cool_Line_H2_e;
152 CL_RT_FLOAT Cool_Line_H2_HII;
153 CL_RT_FLOAT Cool_LowT;
154} RATES_T;
155
156typedef struct clDerivsDataStruct clDerivsData;
157
158/* Heating Cooling Context */
159
160typedef struct CoolingPKDStruct {
161 double z; /* Redshift */
162 double dTime;
163 /* Rates independent of Temperature */
164 RATES_NO_T R;
165 /* Table for Temperature dependent rates */
166 int nTable;
167 double TMin;
168 double TMax;
169 double TlnMin;
170 double TlnMax;
171 double rDeltaTln;
172 RATES_T *RT;
173
174 int bMetal;
175 int nzMetalTable;
176 int nnHMetalTable;
177 int nTMetalTable;
178 double MetalTMin;
179 double MetalTMax;
180 double MetalTlogMin;
181 double MetalTlogMax;
182 double rDeltaTlog;
183 double MetalnHMin;
184 double MetalnHMax;
185 double MetalnHlogMin;
186 double MetalnHlogMax;
187 double rDeltanHlog;
188 double MetalzMin;
189 double MetalzMax;
190 double rDeltaz;
191 float ***MetalCoolln;
192 float ***MetalHeatln;
193 double *Rate_DustForm_H2;
194
195 int nTableRead; /* number of Tables read from files */
196
197 int bUV;
198 int nUV;
199 UVSPECTRUM *UV;
200 int bUVTableUsesTime;
201 int bUVTableLinear;
202 int bLowTCool;
203 int bSelfShield;
204
205 int bShieldHI;
206 double dClump; /* Subgrid clumping factor for determining rate of H2 formation on dust. 10 is a good value*/
207 double dLymanWernerFrac; /* Set to true to determine age of star particle from mass compared to formation mass when calculating LW radiation. Useful in running ICs which already have stars*/
208 double dGmPerCcUnit;
209 double dComovingGmPerCcUnit;
210 double dExpand; /*cosmological expansion factor*/
211 double dErgPerGmUnit;
212 double dSecUnit;
213 double dErgPerGmPerSecUnit;
214 double diErgPerGmUnit;
215 double dKpcUnit;
216#ifdef SHIELDSF
217 double dMaxTShield; /* Temperature ceiling for dust shielding */
218#endif
219 double dMsolUnit;
220 double dMassFracHelium;
221
222/* Diagnostic */
223 int its;
224#if defined(COOLDEBUG)
225 /* MDL mdl; *//* For diag/debug outputs */
226 /*struct particle *p;*/ /* particle pointer needed for SN feedback */
227 int iOrder;
228#endif
229} COOL;
230
231typedef struct {
232 double T, Tln;
233 double Coll_HI;
234 double Coll_HeI;
235 double Coll_HeII;
236 double Coll_e_H2;
237 double Coll_HI_H2;
238 double Coll_H2_H2;
239 double Coll_Hm_e; /*gas phase formation of H2 */
240 double Coll_Hm_HII; /*------------------- */
241 double Coll_HI_e; /*------------------- */
242 double Coll_HII_H2; /*--------------------- */
243 double HI_e; /*---------------------- */
244 double HI_Hm; /*gas phase formation of H2 */
245 double Radr_HII;
246 double Radr_HeII;
247 double Diel_HeII;
248 double Chtr_HeII;
249 double Totr_HeII;
250 double Radr_HeIII;
251 double Cool_Metal;
252 double Heat_Metal;
253
254 double Phot_HI;
255 double Phot_HeI;
256 double Phot_HeII;
257 double Phot_H2; /*Photon dissociation of H2*/
258 double DustForm_H2; /* Formation of H2 on dust */
259 double CorreLength; /* The correlation length of subgrid turbulence, used when calculating shielding*/
260 double LymanWernerCode;
261} RATE;
262
263typedef struct {
264 double compton;
265 double bremHII;
266 double bremHeII;
267 double bremHeIII;
268 double radrecHII;
269 double radrecHeII;
270 double radrecHeIII;
271 double collionHI;
272 double collionHeI;
273 double collionHeII;
274 double collion_e_H2;
275 double collion_H_H2;
276 double collion_H2_H2;
277 double collion_HII_H2;
278 double dielrecHeII;
279 double lineHI;
280 double lineHeI;
281 double lineHeII;
282 double lineH2;
283 double lowT;
284 double NetMetalCool;
285} COOL_ERGPERSPERGM;
286
287
288struct clDerivsDataStruct {
289 STIFF *IntegratorContext;
290 COOL *cl;
291 double rho,ExternalHeating,E,ZMetal,dLymanWerner, columnL;
292/* double Y_H, Y_He; */ /* will be needed -- also for temperature , Y_MetalIon, Y_eMetal */
293 RATE Rate;
294 PERBARYON Y;
295 double Y_H, Y_He, Y_eMax;
296 double Y_Total0, Y_Total1;
297 double dlnE;
298 int its; /* Debug */
299 int bCool;
300};
301
302COOL *CoolInit( );
303void CoolFinalize( COOL *cl );
304clDerivsData *CoolDerivsInit(COOL *cl);
305void CoolDerivsFinalize(clDerivsData *cld ) ;
306
307void clInitConstants( COOL *cl, double dGMPerCcunit, double dComovingGmPerCcUnit,
308 double dErgPerGmUnit, double dSecUnit, double dKpcUnit, COOLPARAM CoolParam);
309void clInitUV(COOL *cl, int nTableColumns, int nTableRows, double *dTableData );
310void clInitRatesTable( COOL *cl, double TMin, double TMax, int nTable );
311void clReadMetalTable(COOL *cl, COOLPARAM clParam);
312void clRateMetalTable(COOL *cl, RATE *Rate, double T, double rho, double Y_H, double ZMetal);
313void clHHeTotal(COOL *cl, double ZMetal);
314void CoolInitRatesTable( COOL *cl, COOLPARAM CoolParam);
315
316void clRatesTableError( COOL *cl );
317void clRatesRedshift( COOL *cl, double z, double dTime );
318double clHeatTotal ( COOL *cl, PERBARYON *Y, RATE *Rate, double rho, double ZMetal );
319void clRates( COOL *cl, RATE *Rate, double T, double rho, double ZMetal, double columnL, double Rate_Phot_H2_stellar);
320double clCoolTotal( COOL *cl, PERBARYON *Y, RATE *Rate, double rho, double ZMetal );
321COOL_ERGPERSPERGM clTestCool ( COOL *cl, PERBARYON *Y, RATE *Rate, double rho );
322void clPrintCool( COOL *cl, PERBARYON *Y, RATE *Rate, double rho );
323void clPrintCoolFile( COOL *cl, PERBARYON *Y, RATE *Rate, double rho, double ZMetal, FILE *fp );
324
325void clAbunds( COOL *cl, PERBARYON *Y, RATE *Rate, double rho, double ZMetal);
326double clThermalEnergy( double Y_Total, double T );
327double clTemperature( double Y_Total, double E );
328double clSelfShield (double yH2, double h);
329double clDustShield (double yHI, double yH2, double z, double h);
330double clRateCollHI( double T );
331double clRateCollHeI( double T );
332double clRateCollHeII( double T );
333double clRateColl_e_H2( double T );
334double clRateColl_HI_H2( double T );
335double clRateColl_H2_H2( double T );
336double clRateColl_HII_H2(double T);
337double clRateColl_Hm_e(double T);
338double clRateColl_HI_e(double T);
339double clRateColl_Hm_HII(double T);
340double clRateHI_e(double T);
341double clRateHI_Hm(double T);
342double clRateRadrHII( double T );
343double clRateRadrHeII( double T );
344double clRateDielHeII( double T );
345double clRateChtrHeII(double T);
346double clRateRadrHeIII( double T );
347double clCoolBrem1( double T );
348double clCoolBrem2( double T );
349double clCoolRadrHII( double T );
350double clCoolRadrHeII( double T );
351double clCoolRadrHeIII( double T );
352double clCoolLineHI( double T );
353double clCoolLineHeI( double T );
354double clCoolLineHeII( double T );
355double clCoolLineH2_table( double T );
356double clCoolLineH2_HI( double T );
357double clCoolLineH2_H2( double T );
358double clCoolLineH2_He( double T );
359double clCoolLineH2_e( double T );
360double clCoolLineH2_HII( double T );
361double clCoolLowT( double T );
362double clRateDustFormH2(double z, double clump);
363double clEdotInstant ( COOL *cl, PERBARYON *Y, RATE *Rate, double rho,
364 double ZMetal, double *dEdotHeat, double *dEdotCool );
365 void clIntegrateEnergy(COOL *cl, clDerivsData *clData, PERBARYON *Y, double *E,
366 double ExternalHeating, double rho, double ZMetal, double dt, double columnL, double dLymanWerner );
367 void clIntegrateEnergyDEBUG(COOL *cl, clDerivsData *clData, PERBARYON *Y, double *E,
368 double ExternalHeating, double rho, double ZMetal, double dt );
369
370
371void clDerivs(double x, const double *y, double *yheat,
372 double *ycool, void *Data) ;
373
374void CoolAddParams( COOLPARAM *CoolParam, PRM );
375void CoolLogParams( COOLPARAM *CoolParam, FILE *fp );
376void CoolOutputArray( COOLPARAM *CoolParam, int, int *, char * );
377
378#define COOL_ARRAY0_EXT "HI"
379double COOL_ARRAY0(COOL *cl, COOLPARTICLE *cp, double ZMetal);
380double COOL_SET_ARRAY0(COOL *cl, COOLPARTICLE *cp,double ZMetal, double val);
381
382#define COOL_ARRAY1_EXT "HeI"
383double COOL_ARRAY1(COOL *cl, COOLPARTICLE *cp, double ZMetal);
384double COOL_SET_ARRAY1(COOL *cl, COOLPARTICLE *cp,double ZMetal, double val);
385
386#define COOL_ARRAY2_EXT "HeII"
387double COOL_ARRAY2(COOL *cl, COOLPARTICLE *cp, double ZMetal);
388double COOL_SET_ARRAY2(COOL *cl, COOLPARTICLE *cp,double ZMetal, double val);
389
390#define COOL_ARRAY3_EXT "H2"
391double COOL_ARRAY3(COOL *cl, COOLPARTICLE *cp, double ZMetal);
392double COOL_SET_ARRAY3(COOL *cl, COOLPARTICLE *cp,double ZMetal, double val);
393
394double COOL_EDOT( COOL *cl_, COOLPARTICLE *cp_, double ECode_, double rhoCode_, double ZMetal_, double *posCode_, double columnL_ );
395#define COOL_EDOT( cl_, cp_, ECode_, rhoCode_, ZMetal_, posCode_, columnL_) (CoolCodeWorkToErgPerGmPerSec( cl_, CoolEdotInstantCode( cl_, cp_, ECode_, rhoCode_, ZMetal_, posCode_ , columnL_)))
396
397double COOL_COOLING( COOL *cl_, COOLPARTICLE *cp_, double ECode_, double rhoCode_, double ZMetal_, double *posCode_ , double columnL_);
398#define COOL_COOLING( cl_, cp_, ECode_, rhoCode_, ZMetal_, posCode_, columnL_) (CoolCodeWorkToErgPerGmPerSec( cl_, CoolCoolingCode( cl_, cp_, ECode_, rhoCode_, ZMetal_, posCode_ , columnL_)))
399
400double COOL_HEATING( COOL *cl_, COOLPARTICLE *cp_, double ECode_, double rhoCode_, double ZMetal_, double *posCode_, double columnL_ );
401#define COOL_HEATING( cl_, cp_, ECode_, rhoCode_, ZMetal_, posCode_, columnL_) (CoolCodeWorkToErgPerGmPerSec( cl_, CoolHeatingCode( cl_, cp_, ECode_, rhoCode_, ZMetal_, posCode_ , columnL_)))
402
403void clSetAbundanceTotals(COOL *cl, double ZMetal, double *Y_H, double *Y_He, double *Y_eMAX);
404void CoolPARTICLEtoPERBARYON(COOL *cl_, PERBARYON *Y, COOLPARTICLE *cp, double ZMetal);
405void CoolPERBARYONtoPARTICLE(COOL *cl_, PERBARYON *Y, COOLPARTICLE *cp, double ZMetal);
406
407double CoolLymanWerner(double dAge);
408
409double CoolEnergyToTemperature( COOL *Cool, COOLPARTICLE *cp, double E, double ZMetal);
410double CoolCodeEnergyToTemperature( COOL *Cool, COOLPARTICLE *cp, double E, double ZMetal);
411
412/* Note: nod to cosmology (z parameter) unavoidable unless we want to access cosmo.[ch] from here */
413void CoolSetTime( COOL *Cool, double dTime, double z );
414
415double CoolCodeTimeToSeconds( COOL *Cool, double dCodeTime );
416
417#define CoolCodeTimeToSeconds( Cool, dCodeTime ) ((Cool)->dSecUnit*(dCodeTime))
418
419double CoolSecondsToCodeTime( COOL *Cool, double dTime );
420
421#define CoolSecondsToCodeTime( Cool, dTime ) ((dTime)/(Cool)->dSecUnit)
422
423double CoolCodeEnergyToErgPerGm( COOL *Cool, double dCodeEnergy );
424
425#define CoolCodeEnergyToErgPerGm( Cool, dCodeEnergy ) ((Cool)->dErgPerGmUnit*(dCodeEnergy))
426
427double CoolErgPerGmToCodeEnergy( COOL *Cool, double dEnergy );
428
429#define CoolErgPerGmToCodeEnergy( Cool, dEnergy ) ((Cool)->diErgPerGmUnit*(dEnergy))
430
431double CoolCodeWorkToErgPerGmPerSec( COOL *Cool, double dCodeWork );
432
433#define CoolCodeWorkToErgPerGmPerSec( Cool, dCodeWork ) ((Cool)->dErgPerGmPerSecUnit*(dCodeWork))
434
435double CoolErgPerGmPerSecToCodeWork( COOL *Cool, double dWork );
436
437#define CoolErgPerGmPerSecToCodeWork( Cool, dWork ) ((dWork)/(Cool)->dErgPerGmPerSecUnit)
438
439double CodeDensityToComovingGmPerCc( COOL *Cool, double dCodeDensity );
440
441#define CodeDensityToComovingGmPerCc( Cool, dCodeDensity ) ((Cool)->dComovingGmPerCcUnit*(dCodeDensity))
442
443void CoolIntegrateEnergy(COOL *cl, clDerivsData *cData, COOLPARTICLE *cp, double *E,
444 double ExternalHeating, double rho, double ZMetal, double tStep, double columnL );
445
446void CoolIntegrateEnergyCode(COOL *cl, clDerivsData *cData, COOLPARTICLE *cp, double *E,
447 double ExternalHeating, double rho, double ZMetal, double *r, double tStep, double columnL );
448
449void CoolDefaultParticleData( COOLPARTICLE *cp );
450
451void CoolInitEnergyAndParticleData( COOL *cl, COOLPARTICLE *cp, double *E, double dDensity, double dTemp, double ZMetal);
452
453/* Deprecated */
454double CoolHeatingRate( COOL *cl, COOLPARTICLE *cp, double E, double dDensity, double ZMetal, double columnL);
455
456double CoolEdotInstantCode(COOL *cl, COOLPARTICLE *cp, double ECode,
457 double rhoCode, double ZMetal, double *posCode, double columnL );
458double CoolCoolingCode(COOL *cl, COOLPARTICLE *cp, double ECode,
459 double rhoCode, double ZMetal, double *posCode, double columnL );
460double CoolHeatingCode(COOL *cl, COOLPARTICLE *cp, double ECode,
461 double rhoCode, double ZMetal, double *posCode, double columnL );
462
463void CoolCodePressureOnDensitySoundSpeed( COOL *cl, COOLPARTICLE *cp, double uPred, double fDensity, double gamma, double gammam1, double *PoverRho, double *c );
464
465/* Note: gamma should be 5/3 for this to be consistent! */
466#define CoolCodePressureOnDensitySoundSpeed( cl__, cp__, uPred__, fDensity__, gamma__, gammam1__, PoverRho__, c__ ) { \
467 *(PoverRho__) = ((5./3.-1)*(uPred__)); \
468 *(c__) = sqrt((5./3.)*(*(PoverRho__))); }
469
470/*
471double CoolCodePressureOnDensity( COOL *cl, COOLPARTICLE *cp, double uPred, double fDensity, double gammam1 );
472
473#define CoolCodePressureOnDensity( cl, cp, uPred, fDensity, gammam1 ) ((gammam1)*(uPred))
474*/
475
476void CoolTableReadInfo( COOLPARAM *CoolParam, int cntTable, int *nTableColumns, char *suffix );
477
478void CoolTableRead( COOL *Cool, int nData, void *vData);
479
480#ifdef __cplusplus
481}
482#endif
483
484#endif
485
486
COOL_ERGPERSPERGM
return structure for clTestCool()
Definition cooling_cosmo.h:187
CoolingPKDStruct
Heating/Cooling context: parameters and tables.
Definition cooling_boley.h:83
CoolingParametersStruct
Input parameters for cooling.
Definition cooling_boley.h:56
CoolingParticleStruct
per-particle cooling data
Definition cooling_boley.h:68
PERBARYON
abundance of various species in particles/baryon
Definition cooling_boley.h:75
RATES_NO_T
structure to hold Temperature independent cooling and heating rates
Definition cooling_cosmo.h:86
RATES_T
structure to hold Temperature dependent cooling rates
Definition cooling_cosmo.h:106
RATE
Rate information for a given particle.
Definition cooling_cosmo.h:170
UVSPECTRUM
photoionization and heating rates from a uniform UV background
Definition cooling_cosmo.h:72
clDerivsDataStruct
context for calculating cooling derivatives
Definition cooling_boley.h:108
clDerivsDataStruct::cl
COOL * cl
pointer to cooling context
Definition cooling_boley.h:110
Generated by doxygen 1.13.2