changa 3.5
Loading...
Searching...
No Matches
cooling_boley.h
1#ifndef COOLING_BOLEY_HINCLUDED
2#define COOLING_BOLEY_HINCLUDED
3
4/*
5 * Planetary disk cooling code as described in Boley 2009, ApJ 695,L53
6 * and Boley et al 2010, Icarus 207, 509. The cooling is calculated
7 * from \grad \cdot F \sim (T^4 - T^4_{irr})/(\Delta \tau + 1/\Delta
8 * \tau). This approximates an irradiated disk.
9 *
10 * This implementation is taken from the GASOLINE implementation by
11 * Hayfield.
12 */
13
14/* Global consts */
15
16#include "param.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))
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 * Work around for Dec ev6 flawed
31 * treatment of sub-normal numbers
32 */
33#define CL_MAX_NEG_EXP_ARG -500.
34
35#define CL_NMAXBYTETABLE 56000
36
37/* The opacity tables are expected to have a total of CL_TABRC2
38 * entries arranged in CL_TABRC rows of length CL_TABRC. Each row is
39 * for a given pressure, with temperature varying along the column.
40 * The pressures are expected to range from log(P(cgs)) = CL_TABPMIN
41 * to CL_TABPMAX. The temperatures are expected to range from log(T
42 * (K)) = CL_TABTMIN to CL_TABTMAX. */
43
44#define CL_TABRC 100
45#define CL_TABRCM2 ((CL_TABRC)-2)
46#define CL_TABRC2 10000
47
48#define CL_TABPMIN -12.0
49#define CL_TABPMAX 9.0
50#define CL_TABDP (((CL_TABPMAX) - (CL_TABPMIN))/(CL_TABRC-1.0))
51
52#define CL_TABTMIN 0.5
53#define CL_TABTMAX 7.0
54#define CL_TABDT (((CL_TABTMAX) - (CL_TABTMIN))/(CL_TABRC-1.0))
55
56typedef struct CoolingParametersStruct {
57 double dParam1;
58 double dRhoCutoff;
59 double dParam3;
60 double dParam4;
61 double Y_Total;
62 double dCoolingTmin;
63 double dCoolingTmax;
64 char achRossName[256];
65 char achPlckName[256];
66 } COOLPARAM;
67
68typedef struct CoolingParticleStruct {
69 double Y_Total;
70 double mrho; /* (m/rho)^(1/3), i.e V^{1/3} in Boley 2009 */
71 double dDeltaTau; /* Optical depth across particle */
72 } COOLPARTICLE;
73
75typedef struct {
76 double Total;
77} PERBARYON;
78
79typedef struct clDerivsDataStruct clDerivsData;
80
81/* Heating Cooling Context */
82
83typedef struct CoolingPKDStruct {
84 double z; /* Redshift */
85 double dTime;
86
87 double dGmPerCcUnit;
88 double dComovingGmPerCcUnit;
89 double dErgPerGmUnit;
90 double dSecUnit;
91 double dErgPerGmPerSecUnit;
92 double diErgPerGmUnit;
93 double dKpcUnit;
94
95 double dParam1;
96 double dRhoCutoff;
97 double dParam3;
98 double dParam4;
99 double Y_Total;
100 double Tmin;
101 double Tmax;
102 clDerivsData *DerivsData;
103 double ** rossTab, ** plckTab, t4i;
104
105 int nTableRead; /* Internal Tables read from Files */
106} COOL;
107
108struct clDerivsDataStruct {
109 void *IntegratorContext;
110 COOL *cl;
111 double rho,PdV,E,T,Y_Total,rFactor;
112 double dDeltaTau; /* optical depth across particle */
113 double dlnE;
114 int its; /* Debug */
115};
116
117
118COOL *CoolInit( );
119void CoolFinalize( COOL *cl );
120clDerivsData *CoolDerivsInit(COOL *cl);
121void CoolDerivsFinalize(clDerivsData *cld ) ;
122
123void clInitConstants( COOL *cl, double dGMPerCcunit, double dComovingGmPerCcUnit,
124 double dErgPerGmUnit, double dSecUnit, double dKpcUnit, COOLPARAM CoolParam);
125void clInitUV(COOL *cl, int nTableColumns, int nTableRows, double *dTableData );
126void clInitRatesTable( COOL *cl, double TMin, double TMax, int nTable );
127void CoolInitRatesTable( COOL *cl, COOLPARAM CoolParam);
128
129double clThermalEnergy(double Y_Total, double T );
130double clTemperature(double Y_Total, double E );
131
132double clEdotInstant( COOL *cl, double E, double T, double rho, double r,
133 double *dDeltaTau,
134 double *dEdotHeat, double *dEdotCool);
135void clIntegrateEnergy(COOL *cl, clDerivsData *clData, double *E,
136 double PdV, double rho, double Y_Total, double radius, double tStep );
137
138void clDerivs(double x, const double *y, double *dHeat, double *dCool, void *Data) ;
139
140int clJacobn( double x, const double y[], double dfdx[], double *dfdy, void *Data) ;
141
142void CoolAddParams( COOLPARAM *CoolParam, PRM );
143void CoolLogParams( COOLPARAM *CoolParam, FILE *fp );
144void CoolOutputArray( COOLPARAM *CoolParam, int, int *, char * );
145
146#define COOL_ARRAY0_EXT "Y_Tot"
147double COOL_ARRAY0(COOL *cl_, COOLPARTICLE *cp,double aa);
148#define COOL_ARRAY0( cl_, cp, aa ) ((cp)->Y_Total)
149double COOL_SET_ARRAY0(COOL *cl_, COOLPARTICLE *cp,double aa, double bb_val);
150#define COOL_SET_ARRAY0( cl_, cp, aa, bb_val ) ((cp)->Y_Total = (bb_val))
151
152#define COOL_ARRAY1_EXT "dtau"
153double COOL_ARRAY1(COOL *cl_, COOLPARTICLE *cp,double aa);
154#define COOL_ARRAY1( cl_, cp, aa ) ((cp)->dDeltaTau)
155double COOL_SET_ARRAY1(COOL *cl_, COOLPARTICLE *cp,double aa, double bb_val);
156#define COOL_SET_ARRAY1( cl_, cp, aa, bb_val ) ((cp)->dDeltaTau = (bb_val))
157
158#define COOL_ARRAY2_EXT "NA"
159double COOL_ARRAY2(COOL *cl_, COOLPARTICLE *cp,double aa);
160#define COOL_ARRAY2( cl_, cp, aa ) (0)
161double COOL_SET_ARRAY2(COOL *cl_, COOLPARTICLE *cp,double aa, double bb_val);
162#define COOL_SET_ARRAY2( cl_, cp, aa, bb_val ) (0)
163
164#define COOL_ARRAY3_EXT "NA3"
165double COOL_ARRAY3(COOL *cl, COOLPARTICLE *cp, double aa);
166#define COOL_ARRAY3(cl_, cp, aa ) (0)
167double COOL_SET_ARRAY3(COOL *cl_, COOLPARTICLE *cp,double aa, double bb_val);
168#define COOL_SET_ARRAY3( cl_, cp, aa, bb_val ) (0)
169
170double COOL_EDOT( COOL *cl_, COOLPARTICLE *cp_, double ECode_, double rhoCode_, double ZMetal_, double *posCode_ );
171#define COOL_EDOT( cl_, cp_, ECode_, rhoCode_, ZMetal_, posCode_) (CoolCodeWorkToErgPerGmPerSec( cl_, CoolEdotInstantCode( cl_, cp_, ECode_, rhoCode_, ZMetal_, posCode_ )))
172
173double COOL_COOLING( COOL *cl_, COOLPARTICLE *cp_, double ECode_, double rhoCode_, double ZMetal_, double *posCode_ );
174#define COOL_COOLING( cl_, cp_, ECode_, rhoCode_, ZMetal_, posCode_) (CoolCodeWorkToErgPerGmPerSec( cl_, CoolEdotInstantCode( cl_, cp_, ECode_, rhoCode_, ZMetal_, posCode_ )))
175
176double COOL_HEATING( COOL *cl_, COOLPARTICLE *cp_, double ECode_, double rhoCode_, double ZMetal_, double *posCode_ );
177#define COOL_HEATING( cl_, cp_, ECode_, rhoCode_, ZMetal_, posCode_) (0)
178
179void CoolPARTICLEtoPERBARYON(COOL *cl_, PERBARYON *Y, COOLPARTICLE *cp);
180
181#define CoolPARTICLEtoPERBARYON(cl_, Y, cp) { \
182 (Y)->Total = (cp)->Y_Total; }
183
184double CoolCodeEnergyToTemperature( COOL *Cool, COOLPARTICLE *cp, double E,
185 double fMetals );
186
187/* Note: nod to cosmology (z parameter) unavoidable unless we want to access cosmo.[ch] from here */
188void CoolSetTime( COOL *Cool, double dTime, double z );
189
190double CoolCodeTimeToSeconds( COOL *Cool, double dCodeTime );
191
192#define CoolCodeTimeToSeconds( Cool, dCodeTime ) ((Cool)->dSecUnit*(dCodeTime))
193
194double CoolSecondsToCodeTime( COOL *Cool, double dTime );
195
196#define CoolSecondsToCodeTime( Cool, dTime ) ((dTime)/(Cool)->dSecUnit)
197
198double CoolCodeEnergyToErgPerGm( COOL *Cool, double dCodeEnergy );
199
200#define CoolCodeEnergyToErgPerGm( Cool, dCodeEnergy ) ((Cool)->dErgPerGmUnit*(dCodeEnergy))
201
202double CoolErgPerGmToCodeEnergy( COOL *Cool, double dEnergy );
203
204#define CoolErgPerGmToCodeEnergy( Cool, dEnergy ) ((Cool)->diErgPerGmUnit*(dEnergy))
205
206double CoolCodeWorkToErgPerGmPerSec( COOL *Cool, double dCodeWork );
207
208#define CoolCodeWorkToErgPerGmPerSec( Cool, dCodeWork ) ((Cool)->dErgPerGmPerSecUnit*(dCodeWork))
209
210double CoolErgPerGmPerSecToCodeWork( COOL *Cool, double dWork );
211
212#define CoolErgPerGmPerSecToCodeWork( Cool, dWork ) ((dWork)/(Cool)->dErgPerGmPerSecUnit)
213
214double CodeDensityToComovingGmPerCc( COOL *Cool, double dCodeDensity );
215
216#define CodeDensityToComovingGmPerCc( Cool, dCodeDensity ) ((Cool)->dComovingGmPerCcUnit*(dCodeDensity))
217
218void CoolIntegrateEnergy(COOL *cl, clDerivsData *clData, COOLPARTICLE *cp,
219 double *E, double PdV, double rho, double ZMetal,
220 double tStep );
221
222void CoolIntegrateEnergyCode(COOL *cl, clDerivsData *clData, COOLPARTICLE *cp, double *E,
223 double PdV, double rho, double ZMetal, double *r, double tStep );
224
225void CoolDefaultParticleData( COOLPARTICLE *cp );
226
227void CoolInitEnergyAndParticleData( COOL *cl, COOLPARTICLE *cp, double *E, double dDensity, double dTemp, double fMetals );
228
229/* Deprecated */
230double CoolHeatingRate( COOL *cl, COOLPARTICLE *cp, double E, double dDensity );
231
232double CoolEdotInstantCode(COOL *cl, COOLPARTICLE *cp, double ECode,
233 double rhoCode, double ZMetal, double *posCode );
234
235void CoolCodePressureOnDensitySoundSpeed( COOL *cl, COOLPARTICLE *cp, double uPred, double fDensity, double gamma, double gammam1, double *PoverRho, double *c );
236#define CoolCodePressureOnDensitySoundSpeed( cl__, cp__, uPred__, fDensity__, gamma__, gammam1__, PoverRho__, c__ ) { \
237 *(PoverRho__) = ((gammam1__)*(uPred__)); \
238 *(c__) = sqrt((gamma__)*(*(PoverRho__))); }
239
240/*
241double CoolCodePressureOnDensity( COOL *cl, COOLPARTICLE *cp, double uPred, double fDensity, double gammam1 );
242
243#define CoolCodePressureOnDensity( cl, cp, uPred, fDensity, gammam1 ) ((gammam1)*(uPred))
244*/
245
246void CoolTableReadInfo( COOLPARAM *CoolParam, int cntTable, int *nTableColumns, char *suffix );
247
248void CoolTableRead( COOL *Cool, int nData, void *vData);
249
250#ifdef __cplusplus
251}
252#endif
253
254#endif
255
256
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
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