changa  3.5
 All Classes Files Functions Variables Typedefs Enumerations Friends Macros Groups Pages
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
19 extern "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 
56 typedef 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 
68 typedef 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 
75 typedef struct {
76  double Total;
77 } PERBARYON;
78 
79 typedef struct clDerivsDataStruct clDerivsData;
80 
81 /* Heating Cooling Context */
82 
83 typedef 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 
108 struct 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 
118 COOL *CoolInit( );
119 void CoolFinalize( COOL *cl );
120 clDerivsData *CoolDerivsInit(COOL *cl);
121 void CoolDerivsFinalize(clDerivsData *cld ) ;
122 
123 void clInitConstants( COOL *cl, double dGMPerCcunit, double dComovingGmPerCcUnit,
124  double dErgPerGmUnit, double dSecUnit, double dKpcUnit, COOLPARAM CoolParam);
125 void clInitUV(COOL *cl, int nTableColumns, int nTableRows, double *dTableData );
126 void clInitRatesTable( COOL *cl, double TMin, double TMax, int nTable );
127 void CoolInitRatesTable( COOL *cl, COOLPARAM CoolParam);
128 
129 double clThermalEnergy(double Y_Total, double T );
130 double clTemperature(double Y_Total, double E );
131 
132 double clEdotInstant( COOL *cl, double E, double T, double rho, double r,
133  double *dDeltaTau,
134  double *dEdotHeat, double *dEdotCool);
135 void clIntegrateEnergy(COOL *cl, clDerivsData *clData, double *E,
136  double PdV, double rho, double Y_Total, double radius, double tStep );
137 
138 void clDerivs(double x, const double *y, double *dHeat, double *dCool, void *Data) ;
139 
140 int clJacobn( double x, const double y[], double dfdx[], double *dfdy, void *Data) ;
141 
142 void CoolAddParams( COOLPARAM *CoolParam, PRM );
143 void CoolLogParams( COOLPARAM *CoolParam, FILE *fp );
144 void CoolOutputArray( COOLPARAM *CoolParam, int, int *, char * );
145 
146 #define COOL_ARRAY0_EXT "Y_Tot"
147 double COOL_ARRAY0(COOL *cl_, COOLPARTICLE *cp,double aa);
148 #define COOL_ARRAY0( cl_, cp, aa ) ((cp)->Y_Total)
149 double 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"
153 double COOL_ARRAY1(COOL *cl_, COOLPARTICLE *cp,double aa);
154 #define COOL_ARRAY1( cl_, cp, aa ) ((cp)->dDeltaTau)
155 double 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"
159 double COOL_ARRAY2(COOL *cl_, COOLPARTICLE *cp,double aa);
160 #define COOL_ARRAY2( cl_, cp, aa ) (0)
161 double 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"
165 double COOL_ARRAY3(COOL *cl, COOLPARTICLE *cp, double aa);
166 #define COOL_ARRAY3(cl_, cp, aa ) (0)
167 double COOL_SET_ARRAY3(COOL *cl_, COOLPARTICLE *cp,double aa, double bb_val);
168 #define COOL_SET_ARRAY3( cl_, cp, aa, bb_val ) (0)
169 
170 double 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 
173 double 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 
176 double 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 
179 void CoolPARTICLEtoPERBARYON(COOL *cl_, PERBARYON *Y, COOLPARTICLE *cp);
180 
181 #define CoolPARTICLEtoPERBARYON(cl_, Y, cp) { \
182  (Y)->Total = (cp)->Y_Total; }
183 
184 double 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 */
188 void CoolSetTime( COOL *Cool, double dTime, double z );
189 
190 double CoolCodeTimeToSeconds( COOL *Cool, double dCodeTime );
191 
192 #define CoolCodeTimeToSeconds( Cool, dCodeTime ) ((Cool)->dSecUnit*(dCodeTime))
193 
194 double CoolSecondsToCodeTime( COOL *Cool, double dTime );
195 
196 #define CoolSecondsToCodeTime( Cool, dTime ) ((dTime)/(Cool)->dSecUnit)
197 
198 double CoolCodeEnergyToErgPerGm( COOL *Cool, double dCodeEnergy );
199 
200 #define CoolCodeEnergyToErgPerGm( Cool, dCodeEnergy ) ((Cool)->dErgPerGmUnit*(dCodeEnergy))
201 
202 double CoolErgPerGmToCodeEnergy( COOL *Cool, double dEnergy );
203 
204 #define CoolErgPerGmToCodeEnergy( Cool, dEnergy ) ((Cool)->diErgPerGmUnit*(dEnergy))
205 
206 double CoolCodeWorkToErgPerGmPerSec( COOL *Cool, double dCodeWork );
207 
208 #define CoolCodeWorkToErgPerGmPerSec( Cool, dCodeWork ) ((Cool)->dErgPerGmPerSecUnit*(dCodeWork))
209 
210 double CoolErgPerGmPerSecToCodeWork( COOL *Cool, double dWork );
211 
212 #define CoolErgPerGmPerSecToCodeWork( Cool, dWork ) ((dWork)/(Cool)->dErgPerGmPerSecUnit)
213 
214 double CodeDensityToComovingGmPerCc( COOL *Cool, double dCodeDensity );
215 
216 #define CodeDensityToComovingGmPerCc( Cool, dCodeDensity ) ((Cool)->dComovingGmPerCcUnit*(dCodeDensity))
217 
218 void CoolIntegrateEnergy(COOL *cl, clDerivsData *clData, COOLPARTICLE *cp,
219  double *E, double PdV, double rho, double ZMetal,
220  double tStep );
221 
222 void CoolIntegrateEnergyCode(COOL *cl, clDerivsData *clData, COOLPARTICLE *cp, double *E,
223  double PdV, double rho, double ZMetal, double *r, double tStep );
224 
225 void CoolDefaultParticleData( COOLPARTICLE *cp );
226 
227 void CoolInitEnergyAndParticleData( COOL *cl, COOLPARTICLE *cp, double *E, double dDensity, double dTemp, double fMetals );
228 
229 /* Deprecated */
230 double CoolHeatingRate( COOL *cl, COOLPARTICLE *cp, double E, double dDensity );
231 
232 double CoolEdotInstantCode(COOL *cl, COOLPARTICLE *cp, double ECode,
233  double rhoCode, double ZMetal, double *posCode );
234 
235 void 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 /*
241 double CoolCodePressureOnDensity( COOL *cl, COOLPARTICLE *cp, double uPred, double fDensity, double gammam1 );
242 
243 #define CoolCodePressureOnDensity( cl, cp, uPred, fDensity, gammam1 ) ((gammam1)*(uPred))
244 */
245 
246 void CoolTableReadInfo( COOLPARAM *CoolParam, int cntTable, int *nTableColumns, char *suffix );
247 
248 void CoolTableRead( COOL *Cool, int nData, void *vData);
249 
250 #ifdef __cplusplus
251 }
252 #endif
253 
254 #endif
255 
256 
clDerivsDataStruct::cl
COOL * cl
pointer to cooling context
Definition: cooling_boley.h:110
CoolingParametersStruct
Input parameters for cooling.
Definition: cooling_boley.h:56
prmContext
Object containing the parameter information.
Definition: param.h:38
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
Generated by   doxygen 1.8.5