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