changa  3.5
 All Classes Files Functions Variables Typedefs Enumerations Friends Macros Groups Pages
cooling_cosmo.h
1 #ifndef COOLING_COSMO_HINCLUDED
2 #define COOLING_COSMO_HINCLUDED
3 
4 /*
5  * Cooling code for cosmology simulations.
6  * Originally written by James Wadsley, McMaster University for
7  * GASOLINE.
8  */
9 
10 /* Global consts */
11 
12 #include "param.h"
13 
14 #ifdef __cplusplus
15 extern "C" {
16 #endif
17 
18 #include "stiff.h"
19 
20 /* Constants */
21 #define CL_B_gm (6.022e23*(938.7830/931.494))
22 #define CL_k_Boltzmann 1.38066e-16
23 #define CL_eV_erg 1.60219e-12
24 #define CL_eV_per_K (CL_k_Boltzmann/CL_eV_erg)
25 /*
26 #define CL_RT_FLOAT float
27 #define CL_RT_MIN 1e-38
28 #define CL_RT_MIN FLT_MIN
29 */
30 
31 #define CL_RT_FLOAT double
32 #define CL_RT_MIN 1e-100
33 
34 /*
35 #define CL_RT_MIN DBL_MIN
36 */
37 /*
38  * Work around for Dec ev6 flawed
39  * treatment of sub-normal numbers
40  */
41 #define CL_MAX_NEG_EXP_ARG -500.
42 
43 #define CL_NMAXBYTETABLE 56000
44 
46 typedef struct CoolingParametersStruct {
47  int bIonNonEqm;
48  int nCoolingTable;
49  int bUV;
50  int bUVTableUsesTime;
51  int bDoIonOutput;
52  int bLowTCool;
53  int bSelfShield;
54  double dMassFracHelium;
55  double dCoolingTmin;
56  double dCoolingTmax;
57  } COOLPARAM;
58 
60 typedef struct CoolingParticleStruct {
61  double Y_HI,Y_HeI,Y_HeII; /* Abundance of ions */
62  double dLymanWerner; /* Flux of Lyman Werner radiation at the gas particle. Temporary CC */
63  } COOLPARTICLE;
64 
66 typedef struct {
67  double e,Total;
68  double HI,HII,HeI,HeII,HeIII;
69 } PERBARYON;
70 
72 typedef struct {
73  double zTime;
74 
75  double Rate_Phot_HI;
76  double Rate_Phot_HeI;
77  double Rate_Phot_HeII;
78 
79  double Heat_Phot_HI;
80  double Heat_Phot_HeI;
81  double Heat_Phot_HeII;
82 } UVSPECTRUM;
83 
86 typedef struct {
87  double Rate_Phot_HI;
88  double Rate_Phot_HeI;
89  double Rate_Phot_HeII;
90 
91  double Heat_Phot_HI;
92  double Heat_Phot_HeI;
93  double Heat_Phot_HeII;
94 
95  double Cool_Coll_HI;
96  double Cool_Coll_HeI;
97  double Cool_Coll_HeII;
98  double Cool_Diel_HeII;
99  double Cool_Comp;
100  double Tcmb;
101  double Cool_LowTFactor;
102 
103 } RATES_NO_T;
104 
106 typedef struct {
107  CL_RT_FLOAT Rate_Coll_HI;
108  CL_RT_FLOAT Rate_Coll_HeI;
109  CL_RT_FLOAT Rate_Coll_HeII;
110  CL_RT_FLOAT Rate_Radr_HII;
111  CL_RT_FLOAT Rate_Radr_HeII;
112  CL_RT_FLOAT Rate_Radr_HeIII;
113  CL_RT_FLOAT Rate_Diel_HeII;
114 
115  CL_RT_FLOAT Cool_Brem_1;
116  CL_RT_FLOAT Cool_Brem_2;
117  CL_RT_FLOAT Cool_Radr_HII;
118  CL_RT_FLOAT Cool_Radr_HeII;
119  CL_RT_FLOAT Cool_Radr_HeIII;
120  CL_RT_FLOAT Cool_Line_HI;
121  CL_RT_FLOAT Cool_Line_HeI;
122  CL_RT_FLOAT Cool_Line_HeII;
123  CL_RT_FLOAT Cool_LowT;
124 
125 } RATES_T;
126 
127 typedef struct clDerivsDataStruct clDerivsData;
128 
131 typedef struct CoolingPKDStruct {
132  double z; /* Redshift */
133  double dTime;
134  /* Rates independent of Temperature */
135  RATES_NO_T R;
136  /* Table for Temperature dependent rates */
137  int nTable;
138  double TMin;
139  double TMax;
140  double TlnMin;
141  double TlnMax;
142  double rDeltaTln;
143  RATES_T *RT;
144 
145  int nTableRead; /* number of Tables read from files */
146 
147  int bUV;
148  int nUV;
149  UVSPECTRUM *UV;
150  int bUVTableUsesTime;
151  int bUVTableLinear;
152  int bLowTCool;
153  int bSelfShield;
154 
155  double dGmPerCcUnit;
156  double dComovingGmPerCcUnit;
157  double dErgPerGmUnit;
158  double dSecUnit;
159  double dErgPerGmPerSecUnit;
160  double diErgPerGmUnit;
161  double dKpcUnit;
162  double Y_H;
163  double Y_He;
164  double Y_eMAX;
165 /* Diagnostic */
166  int its;
167 } COOL;
168 
170 typedef struct {
171  double T, Tln;
172  double Coll_HI;
173  double Coll_HeI;
174  double Coll_HeII;
175  double Radr_HII;
176  double Radr_HeII;
177  double Diel_HeII;
178  double Totr_HeII;
179  double Radr_HeIII;
180 
181  double Phot_HI;
182  double Phot_HeI;
183  double Phot_HeII;
184 } RATE;
185 
187 typedef struct {
188  double compton;
189  double bremHII;
190  double bremHeII;
191  double bremHeIII;
192  double radrecHII;
193  double radrecHeII;
194  double radrecHeIII;
195  double collionHI;
196  double collionHeI;
197  double collionHeII;
198  double dielrecHeII;
199  double lineHI;
200  double lineHeI;
201  double lineHeII;
202  double lowT;
203 } COOL_ERGPERSPERGM;
204 
205 
207 struct clDerivsDataStruct {
208  STIFF *IntegratorContext;
209  COOL *cl;
210  double rho,ExternalHeating,E,ZMetal;
211  RATE Rate;
212  PERBARYON Y;
213  double Y_Total0, Y_Total1;
214  int its; /* Debug */
215 };
216 
217 
218 COOL *CoolInit( );
219 void CoolFinalize( COOL *cl );
220 clDerivsData *CoolDerivsInit(COOL *cl);
221 void CoolDerivsFinalize(clDerivsData *cld ) ;
222 
223 void clInitConstants( COOL *cl, double dGMPerCcunit, double dComovingGmPerCcUnit,
224  double dErgPerGmUnit, double dSecUnit, double dKpcUnit, COOLPARAM CoolParam);
225 void clInitUV(COOL *cl, int nTableColumns, int nTableRows, double *dTableData );
226 void clInitRatesTable( COOL *cl, double TMin, double TMax, int nTable );
227 void CoolInitRatesTable( COOL *cl, COOLPARAM CoolParam);
228 
229 void clRatesTableError( COOL *cl );
230 void clRatesRedshift( COOL *cl, double z, double dTime );
231 double clHeatTotal ( COOL *cl, PERBARYON *Y, RATE *Rate );
232 void clRates( COOL *cl, RATE *Rate, double T, double rho );
233 double clCoolTotal( COOL *cl, PERBARYON *Y, RATE *Rate, double rho, double ZMetal );
234 COOL_ERGPERSPERGM clTestCool ( COOL *cl, PERBARYON *Y, RATE *Rate, double rho );
235 void clPrintCool( COOL *cl, PERBARYON *Y, RATE *Rate, double rho );
236 void clPrintCoolFile( COOL *cl, PERBARYON *Y, RATE *Rate, double rho, FILE *fp );
237 
238 void clAbunds( COOL *cl, PERBARYON *Y, RATE *Rate, double rho);
239 double clThermalEnergy( double Y_Total, double T );
240 double clTemperature( double Y_Total, double E );
241 double clRateCollHI( double T );
242 double clRateCollHeI( double T );
243 double clRateCollHeII( double T );
244 double clRateRadrHII( double T );
245 double clRateRadrHeII( double T );
246 double clRateDielHeII( double T );
247 double clRateRadrHeIII( double T );
248 double clCoolBrem1( double T );
249 double clCoolBrem2( double T );
250 double clCoolRadrHII( double T );
251 double clCoolRadrHeII( double T );
252 double clCoolRadrHeIII( double T );
253 double clCoolLineHI( double T );
254 double clCoolLineHeI( double T );
255 double clCoolLineHeII( double T );
256 double clCoolLowT( double T );
257 
258 double clEdotInstant ( COOL *cl, PERBARYON *Y, RATE *Rate, double rho,
259  double ZMetal, double *EdotHeat, double *EdotCool );
260 void clIntegrateEnergy(COOL *cl, clDerivsData *clData, PERBARYON *Y, double *E,
261  double ExternalHeating, double rho, double ZMetal, double dt );
262 void clIntegrateEnergyDEBUG(COOL *cl, PERBARYON *Y, double *E,
263  double ExternalHeating, double rho, double ZMetal, double dt );
264 
265 
266 void clDerivs(double x, const double *y, double *dHeat, double *dCool,
267  void *Data) ;
268 
269 int clJacobn( double x, const double y[], double dfdx[], double *dfdy, void *Data) ;
270 
271 void CoolAddParams( COOLPARAM *CoolParam, PRM );
272 void CoolLogParams( COOLPARAM *CoolParam, FILE *fp );
273 void CoolOutputArray( COOLPARAM *CoolParam, int, int *, char * );
274 
275 #define COOL_ARRAY0_EXT "HI"
276 double COOL_ARRAY0(COOL *cl_, COOLPARTICLE *cp,double aa);
277 #define COOL_ARRAY0( cl_, cp, aa ) ((cp)->Y_HI)
278 double COOL_SET_ARRAY0(COOL *cl_, COOLPARTICLE *cp,double aa, double bb_val);
279 #define COOL_SET_ARRAY0( cl_, cp, aa, bb_val ) ((cp)->Y_HI = (bb_val))
280 
281 #define COOL_ARRAY1_EXT "HeI"
282 double COOL_ARRAY1(COOL *cl_, COOLPARTICLE *cp,double aa);
283 #define COOL_ARRAY1( cl_, cp, aa ) ((cp)->Y_HeI)
284 double COOL_SET_ARRAY1(COOL *cl_, COOLPARTICLE *cp,double aa, double bb_val);
285 #define COOL_SET_ARRAY1( cl_, cp, aa, bb_val ) ((cp)->Y_HeI = (bb_val))
286 
287 #define COOL_ARRAY2_EXT "HeII"
288 double COOL_ARRAY2(COOL *cl_, COOLPARTICLE *cp,double aa);
289 #define COOL_ARRAY2( cl_, cp, aa ) ((cp)->Y_HeII)
290 double COOL_SET_ARRAY2(COOL *cl_, COOLPARTICLE *cp,double aa, double bb_val);
291 #define COOL_SET_ARRAY2( cl_, cp, aa, bb_val ) ((cp)->Y_HeII = (bb_val))
292 
293 #define COOL_ARRAY3_EXT "H2"
294 double COOL_ARRAY3(COOL *cl, COOLPARTICLE *cp, double ZMetal);
295 #define COOL_ARRAY3(cl_, cp, aa ) (0)
296 double COOL_SET_ARRAY3(COOL *cl_, COOLPARTICLE *cp,double aa, double bb_val);
297 #define COOL_SET_ARRAY3( cl_, cp, aa, bb_val ) (0)
298 
299 double COOL_EDOT( COOL *cl_, COOLPARTICLE *cp_, double ECode_, double rhoCode_, double ZMetal_, double *posCode_ );
300 #define COOL_EDOT( cl_, cp_, ECode_, rhoCode_, ZMetal_, posCode_) (CoolCodeWorkToErgPerGmPerSec( cl_, CoolEdotInstantCode( cl_, cp_, ECode_, rhoCode_, ZMetal_, posCode_ )))
301 
302 double COOL_COOLING( COOL *cl_, COOLPARTICLE *cp_, double ECode_, double rhoCode_, double ZMetal_, double *posCode_ );
303 #define COOL_COOLING( cl_, cp_, ECode_, rhoCode_, ZMetal_, posCode_) (CoolCodeWorkToErgPerGmPerSec( cl_, CoolCoolingCode( cl_, cp_, ECode_, rhoCode_, ZMetal_, posCode_ )))
304 
305 double COOL_HEATING( COOL *cl_, COOLPARTICLE *cp_, double ECode_, double rhoCode_, double ZMetal_, double *posCode_ );
306 #define COOL_HEATING( cl_, cp_, ECode_, rhoCode_, ZMetal_, posCode_) (CoolCodeWorkToErgPerGmPerSec( cl_, CoolHeatingCode( cl_, cp_, ECode_, rhoCode_, ZMetal_, posCode_ )))
307 
308 void CoolPARTICLEtoPERBARYON(COOL *cl_, PERBARYON *Y, COOLPARTICLE *cp, double HTotal, double HeTotal);
309 
310 #define CoolPARTICLEtoPERBARYON(cl_, Y, cp) { \
311  (Y)->HI = (cp)->Y_HI; \
312  (Y)->HII = (cl_)->Y_H - (Y)->HI; \
313  (Y)->HeI = (cp)->Y_HeI; \
314  (Y)->HeII = (cp)->Y_HeII; \
315  (Y)->HeIII = (cl_)->Y_He - (Y)->HeI - (Y)->HeII; \
316  (Y)->e = (Y)->HII + (Y)->HeII + 2*(Y)->HeIII; \
317  (Y)->Total = (Y)->e + (cl_)->Y_H + (cl_)->Y_He; }
318 
319 void CoolPERBARYONtoPARTICLE(COOL *cl_, PERBARYON *Y, COOLPARTICLE *cp);
320 
321 #define CoolPERBARYONtoPARTICLE(cl_, Y, cp) { \
322  (cp)->Y_HI = (Y)->HI; \
323  (cp)->Y_HeI = (Y)->HeI; \
324  (cp)->Y_HeII = (Y)->HeII; }
325 
326 
327 double CoolEnergyToTemperature( COOL *Cool, COOLPARTICLE *cp, double E, double fMetal );
328 double CoolCodeEnergyToTemperature( COOL *Cool, COOLPARTICLE *cp, double E, double fMetal );
329 
330 /* Note: nod to cosmology (z parameter) unavoidable unless we want to access cosmo.[ch] from here */
331 void CoolSetTime( COOL *Cool, double dTime, double z );
332 
333 double CoolCodeTimeToSeconds( COOL *Cool, double dCodeTime );
334 
335 #define CoolCodeTimeToSeconds( Cool, dCodeTime ) ((Cool)->dSecUnit*(dCodeTime))
336 
337 double CoolSecondsToCodeTime( COOL *Cool, double dTime );
338 
339 #define CoolSecondsToCodeTime( Cool, dTime ) ((dTime)/(Cool)->dSecUnit)
340 
341 double CoolCodeEnergyToErgPerGm( COOL *Cool, double dCodeEnergy );
342 
343 #define CoolCodeEnergyToErgPerGm( Cool, dCodeEnergy ) ((Cool)->dErgPerGmUnit*(dCodeEnergy))
344 
345 double CoolErgPerGmToCodeEnergy( COOL *Cool, double dEnergy );
346 
347 #define CoolErgPerGmToCodeEnergy( Cool, dEnergy ) ((Cool)->diErgPerGmUnit*(dEnergy))
348 
349 double CoolCodeWorkToErgPerGmPerSec( COOL *Cool, double dCodeWork );
350 
351 #define CoolCodeWorkToErgPerGmPerSec( Cool, dCodeWork ) ((Cool)->dErgPerGmPerSecUnit*(dCodeWork))
352 
353 double CoolErgPerGmPerSecToCodeWork( COOL *Cool, double dWork );
354 
355 #define CoolErgPerGmPerSecToCodeWork( Cool, dWork ) ((dWork)/(Cool)->dErgPerGmPerSecUnit)
356 
357 double CodeDensityToComovingGmPerCc( COOL *Cool, double dCodeDensity );
358 
359 #define CodeDensityToComovingGmPerCc( Cool, dCodeDensity ) ((Cool)->dComovingGmPerCcUnit*(dCodeDensity))
360 
361 void CoolIntegrateEnergy(COOL *cl, clDerivsData *cData, COOLPARTICLE *cp, double *E,
362  double ExternalHeating, double rho, double ZMetal, double tStep );
363 
364 void CoolIntegrateEnergyCode(COOL *cl, clDerivsData *cData, COOLPARTICLE *cp, double *E,
365  double ExternalHeating, double rho, double ZMetal, double *r, double tStep );
366 
367 void CoolDefaultParticleData( COOLPARTICLE *cp );
368 
369 void CoolInitEnergyAndParticleData( COOL *cl, COOLPARTICLE *cp, double *E, double dDensity, double dTemp, double fMetal );
370 
371 /* Deprecated */
372 double CoolHeatingRate( COOL *cl, COOLPARTICLE *cp, double E, double dDensity, double ZMetal );
373 
374 double CoolEdotInstantCode(COOL *cl, COOLPARTICLE *cp, double ECode,
375  double rhoCode, double ZMetal, double *posCode );
376 double CoolCoolingCode(COOL *cl, COOLPARTICLE *cp, double ECode,
377  double rhoCode, double ZMetal, double *posCode );
378 double CoolHeatingCode(COOL *cl, COOLPARTICLE *cp, double ECode,
379  double rhoCode, double ZMetal, double *posCode );
380 
381 void CoolCodePressureOnDensitySoundSpeed( COOL *cl, COOLPARTICLE *cp, double uPred, double fDensity, double gamma, double gammam1, double *PoverRho, double *c );
382 
383 /* Note: gamma should be 5/3 for this to be consistent! */
384 #define CoolCodePressureOnDensitySoundSpeed( cl__, cp__, uPred__, fDensity__, gamma__, gammam1__, PoverRho__, c__ ) { \
385  *(PoverRho__) = ((5./3.-1)*(uPred__)); \
386  *(c__) = sqrt((5./3.)*(*(PoverRho__))); }
387 
388 /*
389 double CoolCodePressureOnDensity( COOL *cl, COOLPARTICLE *cp, double uPred, double fDensity, double gammam1 );
390 
391 #define CoolCodePressureOnDensity( cl, cp, uPred, fDensity, gammam1 ) ((gammam1)*(uPred))
392 */
393 
394 void CoolTableReadInfo( COOLPARAM *CoolParam, int cntTable, int *nTableColumns, char *suffix );
395 
396 void CoolTableRead( COOL *Cool, int nData, void *vData);
397 
398 #ifdef __cplusplus
399 }
400 #endif
401 
402 #endif
403 
404 
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
clDerivsDataStruct::IntegratorContext
STIFF * IntegratorContext
Context for diff. eq. integrator.
Definition: cooling_cosmo.h:208
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