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