PLUTO
Macros | Functions
zeta_tables.c File Reference
#include "pluto.h"
Include dependency graph for zeta_tables.c:

Go to the source code of this file.

Macros

#define THETA_V   6140.0
 
#define THETA_R   85.5
 
#define ORTHO_PARA_MODE   1
 

Functions

static void MakeZetaTables (double *, double *, int)
 
void GetFuncDum (double T, double *funcdum_val)
 

Macro Definition Documentation

#define ORTHO_PARA_MODE   1

Definition at line 13 of file zeta_tables.c.

#define THETA_R   85.5

Definition at line 7 of file zeta_tables.c.

#define THETA_V   6140.0

Definition at line 6 of file zeta_tables.c.

Function Documentation

void GetFuncDum ( double  T,
double *  funcdum_val 
)

Interpolate the value of a function of zetaR from the table that is used to estimate the value of EH2.

Parameters
[in]TValue of temperature in kelvin.
[out]*funcdum_valPointer to the value of function of zetaR
Returns
This function has no return value.

Reference: D'Angelo, G. et al, ApJ 778 2013.

Definition at line 16 of file zeta_tables.c.

29 {
30  int klo, khi, kmid;
31  int nsteps = 5000;
32  double mu, Tmid, dT;
33  static double *lnT, *funcdum;
34  double y, dy;
35  int indx;
36 
37 /* -------------------------------------------
38  Make table on first call
39  ------------------------------------------- */
40 
41  if (lnT == NULL){
42  lnT = ARRAY_1D(nsteps, double);
43  funcdum = ARRAY_1D(nsteps, double);
44  MakeZetaTables(lnT, funcdum, nsteps);
45  }
46  y = log(T);
47 
48 /* -------------------------------------------------
49  Since the table has regular spacing in log T,
50  we divide by the increment to find the nearest
51  node in the table.
52  ------------------------------------------------- */
53 
54  if (y > lnT[nsteps-2]) {
55  *funcdum_val = funcdum[nsteps-2];
56  } else if (y < lnT[0]) {
57  *funcdum_val = funcdum[0];
58  } else{
59  dy = lnT[1] - lnT[0];
60  indx = floor((y - lnT[0])/dy);
61 
62  if (indx >= nsteps || indx < 0){
63  print1 ("! GetFuncDum: indx out of range, indx = %d\n",indx);
64  print1 ("! T = %12.6e\n",T);
65  QUIT_PLUTO(1);
66  }
67  *funcdum_val = (funcdum[indx]*(lnT[indx+1] - y)
68  + funcdum[indx+1]*(y - lnT[indx]))/dy;
69  }
70 }
Sph_disk.T
tuple T
Definition: Sph_disk.py:33
print1
void print1(const char *fmt,...)
Definition: amrPluto.cpp:511
MakeZetaTables
static void MakeZetaTables(double *, double *, int)
Definition: zeta_tables.c:73
ARRAY_1D
#define ARRAY_1D(nx, type)
Definition: prototypes.h:170
QUIT_PLUTO
#define QUIT_PLUTO(e_code)
Definition: macros.h:125

Here is the call graph for this function:

Here is the caller graph for this function:

void MakeZetaTables ( double *  lnT,
double *  funcdum,
int  nsteps 
)
static

Compute tables from iterative summation involed to estimate partition function of parahydrogen ( $ \zeta_P$) and orthohydrogen ( $ \zeta_O$) and their derivatives with respect to temperature. Then further estimate $\zeta_R$ and finally the function of $ \zeta_R$ that goes into the estimation of EH2. The partition function and its derivative can be written as

\[ \left\{\begin{array}{l} \zeta = \DS \sum_i a_ie^{-b_i/T} \\ \noalign{\medskip} \zeta' = \DS \frac{1}{T^2}\sum_i a_ib_ie^{-b_i/T} \end{array}\right. \]

where i=0,2,4,... is even for parahydrogen while i=1,3,5,... is odd for orthoyhdrogen. One can see that in the limit of zero temperatures, $ \zeta_P \to 1, \, \zeta'_P \to 0$ while $\zeta_O\to 0,\,\zeta'_O\to0$. Then $ \zeta'_O/\zeta_O$ is ill-defined in the low temperature limit since it becomes 0/0. To avoid this, we rewrite zetaO by extracting the first term from the summation:

\[ \left\{\begin{array}{l} \zeta_O = \DS e^{-b_1/T}\sum_i a_ie^{-\Delta b_i/T} \\ \noalign{\medskip} \zeta'_O = \DS \frac{1}{T^2}e^{-b_1/T} \left[ b_1 \sum_i a_ie^{-\Delta b_i/T} +\sum_i a_i\Delta b_ie^{-\Delta b_i/T}\right] \end{array}\right. \]

where $\Delta b_i = b_i-b_1$. Taking the ratio:

\[ \frac{\zeta'_O}{\zeta_O} = \frac{1}{T^2} \frac{ b_1 \sum_i a_ie^{-\Delta b_i/T} + \sum_i a_i\Delta b_ie^{-\Delta b_i/T}} {\sum_i a_ie^{-\Delta b_i/T}} \qquad\Longrightarrow\qquad \frac{\zeta'_O}{\zeta_O} - \frac{b_1}{T^2}= \frac{1}{T^2} \frac{\sum_i a_i\Delta b_ie^{-\Delta b_i/T}} {\sum_i a_ie^{-\Delta b_i/T}} \to 0 \quad{\rm as}\quad T\to 0 \]

The expression on the right (that appears in the computation of $\zeta'_R$) is now well-behaved and it reproduces the low temperature limit correctly. You can convince yourself by playing with the following MAPLE script

1 restart;
2 a[0] := 1; a[1] := 3; a[2] := 5; a[3] := 7;
3 b[0] := 0; b[1] := 2*theta; b[2] := 6*theta; b[3] := 12*theta;
4 
5 zeta[P] := a[0]*exp(-b[0]*y) + a[2]*exp(-b[2]*y);
6 dzeta[P] := y^2*(a[0]*b[0]*exp(-b[0]*y) + a[2]*b[2]*exp(-b[2]*y));
7 
8 zeta[O] := a[1]*exp(-b[1]*y) + a[3]*exp(-b[3]*y);
9 dzeta[O] := y^2*(a[1]*b[1]*exp(-b[1]*y) + a[3]*b[3]*exp(-b[3]*y));
10 rP := dzeta[P]/zeta[P];
11 rO := dzeta[O]/zeta[O];
12 
13 rOminus := rO - b[1]*y^2;
14 simplify(rOminus);
Parameters
[in]lnTArray of Logrithmic values of Gas temperatures from 0.01 K to 10^12 K.
[in]nstepsNumber of equal spacings in T.
[out]funcdumThe function of zetaR that goes into EH2.
Returns
This function has no return value

Reference:

  • D'Angelo, G. et. al ApJ 778, 2013 (Eq. 23)

Definition at line 73 of file zeta_tables.c.

149 {
150  int i,j;
151  double Temp0 = 0.01*T_CUT_RHOE;
152  double Tmax = 1.0e12;
153  double dy = log(Tmax/Temp0)*(1./nsteps);
154  double dT = Temp0*exp(dy);
155  double T, a, b, b1, scrh, inv_T2;
156  double zetaP, dzetaP, zetaO, dzetaO, zetaR, dzetaR;
157  double dum1, dum2, dum3;
158  double alpha, beta, gamma;
159  double dzO_zO_m, db, sum1, sum2;
160 
161  print1 ("> MakeZetaTables(): generating Zeta tables...\n");
162 
163  if (ORTHO_PARA_MODE == 0){
164  alpha = 1.0; beta = 0.0; gamma = 0.0;
165  }else if(ORTHO_PARA_MODE == 2){
166  alpha = 0.25; beta = 0.75; gamma = 0.0;
167  }else{
168  alpha = 1.0; beta = 0.0; gamma = 1.0;
169  }
170 
171  b1 = 2.0*THETA_R;
172  for(j = 0; j < nsteps; j++){
173  T = Temp0*exp(j*dy);
174  inv_T2 = 1.0/(T*T);
175  zetaO = zetaP = dzetaP = dzetaO = 0.0;
176  dzO_zO_m = sum1 = sum2 = 0.0;
177  for(i = 0; i <= 10000; i++){
178  a = 2*i + 1;
179  b = i*(i + 1)*THETA_R;
180  if (i%2 == 0){
181  scrh = a*exp(-b/T);
182  zetaP += scrh;
183  dzetaP += scrh*b;
184  }else{
185  db = b - b1;
186  scrh = a*exp(-db/T);
187  sum1 += scrh;
188  sum2 += scrh*db;
189  }
190  }
191  dzetaP *= inv_T2;
192 
193  zetaO = exp(-b1/T)*sum1;
194  dzetaO = exp(-b1/T)*(b1*sum1 + sum2)*inv_T2;
195 
196  dzO_zO_m = sum2/sum1*inv_T2; /* = zeta'(O)/zeta(O) - 2*theta/T^2 */
197 
198  lnT[j] = log(T);
199 
200  /* -----------------------------------------
201  Compute table
202  ----------------------------------------- */
203 
204  scrh = zetaO*exp(2.0*THETA_R/T);
205 
206  zetaR = pow(zetaP,alpha)*pow(scrh,beta) + 3.0*gamma*zetaO;
207  dzetaR = (zetaR - 3.0*gamma*zetaO)*(alpha*(dzetaP/zetaP) +
208  beta*dzO_zO_m) + 3.0*gamma*dzetaO;
209  dum1 = THETA_V/T;
210  dum2 = dum1*exp(-dum1)/(1.0 - exp(-dum1));
211  dum3 = (T/zetaR)*dzetaR;
212  funcdum[j] = 1.5 + dum2 + dum3;
213  }
214 }
alpha
static double alpha
Definition: init.c:3
a
static double a
Definition: init.c:135
Sph_disk.T
tuple T
Definition: Sph_disk.py:33
T_CUT_RHOE
#define T_CUT_RHOE
Sets the lowest cut-off temperature used in the PVTE_LAW equation of state.
Definition: eos.h:69
print1
void print1(const char *fmt,...)
Definition: amrPluto.cpp:511
configure.scrh
tuple scrh
Definition: configure.py:200
THETA_R
#define THETA_R
Definition: zeta_tables.c:7
THETA_V
#define THETA_V
Definition: zeta_tables.c:6
j
int j
Definition: analysis.c:2
i
int i
Definition: analysis.c:2
ORTHO_PARA_MODE
#define ORTHO_PARA_MODE
Definition: zeta_tables.c:13

Here is the call graph for this function:

Here is the caller graph for this function: