forked from N-BodyShop/gasoline
-
Notifications
You must be signed in to change notification settings - Fork 0
/
cooling_metal.h
483 lines (386 loc) · 14.5 KB
/
cooling_metal.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
#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/types.h>
#include <rpc/xdr.h>
#include <sys/param.h> /* for MAXPATHLEN */
#ifndef MAXPATHLEN
#define MAXPATHLEN 256
#endif
/* 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;
double dMassFracHelium;
double dCoolingTmin;
double dCoolingTmax;
double dPhotoelectricHeating;
double dPhotoelectricnMin;
double dPhotoelectricScaleLength;
double dPhotoelectricInnerRadius;
double dzTimeClampUV;
double dMetalCoolFactor;
} COOLPARAM;
typedef struct CoolingParticleStruct {
FLOAT f_HI,f_HeI,f_HeII;
} COOLPARTICLE;
typedef struct {
double e,Total;
double HI,HII,HeI,HeII,HeIII;
} PERBARYON;
typedef struct {
double zTime;
double Rate_Phot_HI;
double Rate_Phot_HeI;
double Rate_Phot_HeII;
double Heat_Phot_HI;
double Heat_Phot_HeI;
double Heat_Phot_HeII;
} UVSPECTRUM;
typedef struct {
double Rate_Phot_HI;
double Rate_Phot_HeI;
double Rate_Phot_HeII;
double Heat_Phot_HI;
double Heat_Phot_HeI;
double Heat_Phot_HeII;
double Heat_Photoelectric;
double nMin_Photoelectric;
double Cool_Coll_HI;
double Cool_Coll_HeI;
double Cool_Coll_HeII;
double Cool_Diel_HeII;
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_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_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;
int nTableRead; /* number of Tables read from files */
int bUV;
int nUV;
UVSPECTRUM *UV;
int bUVTableUsesTime;
int bUVTableLinear;
int bLowTCool;
int bSelfShield;
double dPhotoelectricFactor;
double dPhotoelectricScaleLength;
double dPhotoelectricInnerRadius;
double dzTimeClampUV;
double dMetalCoolFactor;
double dGmPerCcUnit;
double dComovingGmPerCcUnit;
double dErgPerGmUnit;
double dSecUnit;
double dErgPerGmPerSecUnit;
double diErgPerGmUnit;
double dKpcUnit;
double dMassFracHelium;
void *DerivsData;
/* Diagnostic */
int its;
#if defined(COOLDEBUG)
MDL mdl; /* For diag/debug outputs */
struct particle *p; /* particle pointer NEVER TO BE USED EXCEPT FOR DEBUG */
#endif
} COOL;
typedef struct {
double T, Tln;
double Coll_HI;
double Coll_HeI;
double Coll_HeII;
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 Heat_Photoelectric;
double nMin_Photoelectric;
double Phot_HI;
double Phot_HeI;
double Phot_HeII;
} RATE;
typedef struct {
double compton;
double bremHII;
double bremHeII;
double bremHeIII;
double radrecHII;
double radrecHeII;
double radrecHeIII;
double collionHI;
double collionHeI;
double collionHeII;
double dielrecHeII;
double lineHI;
double lineHeI;
double lineHeII;
double lowT;
double NetMetalCool;
} COOL_ERGPERSPERGM;
typedef struct clDerivsDataStruct {
void *IntegratorContext;
COOL *cl;
double rho,ExternalHeating,E,ZMetal;
/* 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;
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, 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 );
void clRates( COOL *cl, RATE *Rate, double T, double rho);
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, 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 clRateCollHI( double T );
double clRateCollHeI( double T );
double clRateCollHeII( 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 clCoolLowT( double T );
double clEdotInstant ( COOL *cl, PERBARYON *Y, RATE *Rate, double rho,
double ZMetal, double *dEdotHeat, double *EdotCool );
void clSetPhotoelectricFactor(COOL *cl, double rkpc);
void clIntegrateEnergy(COOL *cl, PERBARYON *Y, double *E,
double ExternalHeating, double rho, double ZMetal, double rkpc, double dt );
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);
void COOL_IN_ARRAY0(COOL *cl, COOLPARTICLE *cp, double ZMetal, double Data);
#define COOL_ARRAY1_EXT "HeI"
FLOAT COOL_ARRAY1(COOL *cl, COOLPARTICLE *cp, double ZMetal);
void COOL_IN_ARRAY1(COOL *cl, COOLPARTICLE *cp, double ZMetal, double Data);
#define COOL_ARRAY2_EXT "HeII"
FLOAT COOL_ARRAY2(COOL *cl, COOLPARTICLE *cp, double ZMetal);
void COOL_IN_ARRAY2(COOL *cl, COOLPARTICLE *cp, double ZMetal, double Data);
#define COOL_ARRAY3_EXT "dummy"
#define COOL_ARRAY3(x,y,z) 0
#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_ );
#define COOL_EDOT( cl_, cp_, ECode_, rhoCode_, ZMetal_, posCode_) (CoolCodeWorkToErgPerGmPerSec( cl_, CoolEdotInstantCode( cl_, cp_, ECode_, rhoCode_, ZMetal_, posCode_ )))
FLOAT COOL_COOLING( COOL *cl_, COOLPARTICLE *cp_, double ECode_, double rhoCode_, double ZMetal_, double *posCode_ );
#define COOL_COOLING( cl_, cp_, ECode_, rhoCode_, ZMetal_, posCode_) (CoolCodeWorkToErgPerGmPerSec( cl_, CoolCoolingCode( cl_, cp_, ECode_, rhoCode_, ZMetal_, posCode_ )))
FLOAT COOL_HEATING( COOL *cl_, COOLPARTICLE *cp_, double ECode_, double rhoCode_, double ZMetal_, double *posCode_ );
#define COOL_HEATING( cl_, cp_, ECode_, rhoCode_, ZMetal_, posCode_) (CoolCodeWorkToErgPerGmPerSec( cl_, CoolHeatingCode( cl_, cp_, ECode_, rhoCode_, ZMetal_, posCode_ )))
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 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 *rp, double tStep );
void CoolIntegrateEnergyCode(COOL *cl, COOLPARTICLE *cp, double *E,
double ExternalHeating, double rho, double ZMetal, double *r, double tStep );
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 rkpc);
double CoolEdotInstantCode(COOL *cl, COOLPARTICLE *cp, double ECode,
double rhoCode, double ZMetal, double *posCode );
double CoolCoolingCode(COOL *cl, COOLPARTICLE *cp, double ECode,
double rhoCode, double ZMetal, double *posCode );
double CoolHeatingCode(COOL *cl, COOLPARTICLE *cp, double ECode,
double rhoCode, double ZMetal, double *posCode );
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 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