PLUTO
ct_emf.c
Go to the documentation of this file.
1 /* ///////////////////////////////////////////////////////////////////// */
2 /*!
3  \file
4  \brief Store or retrieve the Electromotive Force (EMF).
5 
6  This file provides a database functionality for storing or
7  retrieving EMF components and related information at different
8  points and times in the code.
9 
10  The CT_StoreEMF() function is called immediately after a 1D Riemann
11  solver during the hydro sweeps in order to save Fluxes and
12  characteristic signal velocities into the emf structure for later
13  reuse.
14  The fluxes coming from different sweeps are the different components
15  of the advective part (-v X B) part of the electric field.
16 
17  The function CT_GetEMF() is used to obtain the edge-centered electric
18  field by properly averaging the EMF components previously stored
19  at the zone faces during the 1D sweeps.
20 
21  \author A. Mignone (mignone@ph.unito.it)
22  \date Aug 27, 2014
23 */
24 /* ///////////////////////////////////////////////////////////////////// */
25 #include "pluto.h"
26 
27 static EMF emf;
28 
29 #define eps_UCT_CONTACT 1.e-6
30 #define EX(k,j,i) (vz[k][j][i]*By[k][j][i] - vy[k][j][i]*Bz[k][j][i])
31 #define EY(k,j,i) (vx[k][j][i]*Bz[k][j][i] - vz[k][j][i]*Bx[k][j][i])
32 #define EZ(k,j,i) (vy[k][j][i]*Bx[k][j][i] - vx[k][j][i]*By[k][j][i])
33 
34 /* ********************************************************************* */
35 void CT_StoreEMF (const State_1D *state, int beg, int end, Grid *grid)
36 /*!
37  * Store EMF components and related information available
38  * during 1D sweeps.
39  *
40  * \param [in] state pointer to State_1D structure
41  * \param [in] beg initial index of computation
42  * \param [in] end final index of computation
43  * \param [in] grid pointer to Grid structure;
44  *
45  * \return This function has no return value.
46  *********************************************************************** */
47 {
48  int i, j, k, s;
49 
50 /* ----------------------------------------------------
51  Allocate memory for EMF structure and
52  check for incompatible combinations of algorithms
53  ---------------------------------------------------- */
54 
55  if (emf.ez == NULL){
56 
57  emf.ibeg = emf.iend = 0;
58  emf.jbeg = emf.jend = 0;
59  emf.kbeg = emf.kend = 0;
60 
61  /* -- memory allocation -- */
62 
63  D_EXPAND( ; ,
64  emf.ez = ARRAY_3D(NX3_TOT, NX2_TOT, NX1_TOT, double); ,
65  emf.ex = ARRAY_3D(NX3_TOT, NX2_TOT, NX1_TOT, double);
66  emf.ey = ARRAY_3D(NX3_TOT, NX2_TOT, NX1_TOT, double);
67  )
68 
69  #if CT_EMF_AVERAGE == UCT_CONTACT
70  D_EXPAND(
71  emf.svx = ARRAY_3D(NX3_TOT, NX2_TOT, NX1_TOT, signed char); ,
72  emf.svy = ARRAY_3D(NX3_TOT, NX2_TOT, NX1_TOT, signed char); ,
73  emf.svz = ARRAY_3D(NX3_TOT, NX2_TOT, NX1_TOT, signed char);
74  )
75  #endif
76 
77  D_EXPAND( ; ,
78  emf.ezi = ARRAY_3D(NX3_TOT, NX2_TOT, NX1_TOT, double);
79  emf.ezj = ARRAY_3D(NX3_TOT, NX2_TOT, NX1_TOT, double); ,
80 
81  emf.exj = ARRAY_3D(NX3_TOT, NX2_TOT, NX1_TOT, double);
82  emf.exk = ARRAY_3D(NX3_TOT, NX2_TOT, NX1_TOT, double);
83 
84  emf.eyi = ARRAY_3D(NX3_TOT, NX2_TOT, NX1_TOT, double);
85  emf.eyk = ARRAY_3D(NX3_TOT, NX2_TOT, NX1_TOT, double);
86  )
87 
88  #if CT_EMF_AVERAGE == UCT_HLL
89 
90  D_EXPAND(
91  emf.SxL = ARRAY_3D(NX3_TOT, NX2_TOT, NX1_TOT, double);
92  emf.SxR = ARRAY_3D(NX3_TOT, NX2_TOT, NX1_TOT, double); ,
93 
94  emf.SyL = ARRAY_3D(NX3_TOT, NX2_TOT, NX1_TOT, double);
95  emf.SyR = ARRAY_3D(NX3_TOT, NX2_TOT, NX1_TOT, double); ,
96 
97  emf.SzL = ARRAY_3D(NX3_TOT, NX2_TOT, NX1_TOT, double);
98  emf.SzR = ARRAY_3D(NX3_TOT, NX2_TOT, NX1_TOT, double);
99  )
100  #endif
101  }
102 
103 /* ------------------------------------------------------
104  Store emf components or other necessary 1-D data
105  ------------------------------------------------------ */
106 
107  if (g_dir == IDIR){
108 
109  emf.ibeg = beg; emf.iend = end;
110  for (i = beg; i <= end; i++) {
111 
112  D_EXPAND(emf.ezi[g_k][g_j][i] = -state->flux[i][BX2]; ,
113  ,
114  emf.eyi[g_k][g_j][i] = state->flux[i][BX3]; )
115 
116  #if CT_EMF_AVERAGE == UCT_CONTACT
117  if (state->flux[i][RHO] > eps_UCT_CONTACT) s = 1;
118  else if (state->flux[i][RHO] < -eps_UCT_CONTACT) s = -1;
119  else s = 0;
120 
121  emf.svx[g_k][g_j][i] = s;
122  #endif
123 
124  #if CT_EMF_AVERAGE == UCT_HLL
125  emf.SxL[g_k][g_j][i] = MAX(0.0, -state->SL[i]);
126  emf.SxR[g_k][g_j][i] = MAX(0.0, state->SR[i]);
127  #endif
128 
129  #if CT_EMF_AVERAGE == RIEMANN_2D
130  emf.ezi[g_k][g_j][i] = -2.0*(state->pnt_flx[i][BX2] +
131  state->dff_flx[i][BX2]);
132  #endif
133 
134  }
135 
136  }else if (g_dir == JDIR){
137 
138  emf.jbeg = beg; emf.jend = end;
139  for (j = beg; j <= end; j++) {
140 
141  D_EXPAND( ; ,
142  emf.ezj[g_k][j][g_i] = state->flux[j][BX1]; ,
143  emf.exj[g_k][j][g_i] = -state->flux[j][BX3]; )
144 
145  #if CT_EMF_AVERAGE == UCT_CONTACT
146  if (state->flux[j][RHO] > eps_UCT_CONTACT) s = 1;
147  else if (state->flux[j][RHO] < -eps_UCT_CONTACT) s = -1;
148  else s = 0;
149  emf.svy[g_k][j][g_i] = s;
150  #endif
151 
152  #if CT_EMF_AVERAGE == UCT_HLL
153  emf.SyL[g_k][j][g_i] = MAX(0.0, -state->SL[j]);
154  emf.SyR[g_k][j][g_i] = MAX(0.0, state->SR[j]);
155  #endif
156 
157  #if CT_EMF_AVERAGE == RIEMANN_2D
158 // emf.ezj[g_k][j][g_i] += state->dff_flx[j][BX1];
159  emf.ezj[g_k][j][g_i] = 2.0*(state->pnt_flx[j][BX1] +
160  state->dff_flx[j][BX1]);
161  #endif
162  }
163 
164  }else if (g_dir == KDIR){
165 
166  emf.kbeg = beg; emf.kend = end;
167  for (k = beg; k <= end; k++) {
168 
169  emf.eyk[k][g_j][g_i] = -state->flux[k][BX1];
170  emf.exk[k][g_j][g_i] = state->flux[k][BX2];
171 
172  #if CT_EMF_AVERAGE == UCT_CONTACT
173  if (state->flux[k][RHO] > eps_UCT_CONTACT) s = 1;
174  else if (state->flux[k][RHO] < -eps_UCT_CONTACT) s = -1;
175  else s = 0;
176  emf.svz[k][g_j][g_i] = s;
177  #endif
178 
179  #if CT_EMF_AVERAGE == UCT_HLL
180  emf.SzL[k][g_j][g_i] = MAX(0.0, -state->SL[k]);
181  emf.SzR[k][g_j][g_i] = MAX(0.0, state->SR[k]);
182  #endif
183 
184  }
185  }
186 
187 /* ------------------------------------------------------
188  Store velocity slopes if necessary
189  ------------------------------------------------------ */
190 
191  #if CT_EMF_AVERAGE == UCT_HLL
192  #ifdef CTU
193  if (g_intStage == 2) return;
194  #endif
195 
196  /* -- "end+1" needed to save dvx_dx -- */
197 
198  CT_StoreVelSlopes (&emf, state, beg, end + 1);
199 
200  #endif
201 
202 }
203 
204 /* ********************************************************************* */
205 EMF *CT_GetEMF (const Data *d, Grid *grid)
206 /*!
207  * Retrieve EMF by suitable average of 1D face-centered fluxes.
208  *
209  * \param [in] d
210  * \param [in] grid
211  *
212  * \return a pointer to an edge-centered EMF.
213  *********************************************************************** */
214 {
215  int i, j, k;
216  double ***vx, ***vy, ***vz;
217  double ***Bx, ***By, ***Bz;
218 
219 /* -----------------------------------------------------
220  Return only the resistive part of emf when using
221  super time stepping.
222  ----------------------------------------------------- */
223 
224  #if RESISTIVITY == SUPER_TIME_STEPPING
225  if (g_operatorStep == PARABOLIC_STEP){
226  TOT_LOOP(k,j,i) {
227  #if DIMENSIONS == 3
228  emf.ex[k][j][i] = emf.ey[k][j][i] = 0.0;
229  #endif
230  emf.ez[k][j][i] = 0.0;
231  }
232  CT_AddResistiveEMF(d, grid);
233  return (&emf);
234  }
235  #endif
236 
237 /* -------------------------------------
238  set boundary conditions on
239  face-centered electric fields
240  ------------------------------------- */
241 /*
242  #ifdef CTU
243  if (g_intStage == 2)
244  #endif
245  EMF_BOUNDARY (&emf, grid);
246 */
247 
248 /* ------------------------------------------------------
249  Compute slopes of staggered magnetic fields
250  ------------------------------------------------------ */
251 
252  #if CT_EMF_AVERAGE == UCT_HLL
253  #ifdef CTU
254  if (g_intStage == 1)
255  #endif
256  CT_GetStagSlopes(d->Vs, &emf, grid);
257  #endif
258 
259 /* -----------------------------------------------------
260  Select average
261  ----------------------------------------------------- */
262 
263  #if CT_EMF_AVERAGE == ARITHMETIC || CT_EMF_AVERAGE == RIEMANN_2D
264 
265  CT_EMF_ArithmeticAverage (&emf, 0.25);
266 
267  #elif CT_EMF_AVERAGE == UCT_CONTACT
268 
269  CT_EMF_ArithmeticAverage (&emf, 1.0);
270  CT_EMF_IntegrateToCorner (d, &emf, grid);
271  for (k = emf.kbeg; k <= emf.kend; k++){
272  for (j = emf.jbeg; j <= emf.jend; j++){
273  for (i = emf.ibeg; i <= emf.iend; i++){
274  #if DIMENSIONS == 3
275  emf.ex[k][j][i] *= 0.25;
276  emf.ey[k][j][i] *= 0.25;
277  #endif
278  emf.ez[k][j][i] *= 0.25;
279  }}}
280 
281  #elif CT_EMF_AVERAGE == UCT_HLL
282  #ifdef CTU
283  if (g_intStage == 1) CT_EMF_CMUSCL_Average (d, &emf, grid);
284  else
285  #endif
286  CT_EMF_HLL_Solver (d, &emf, grid);
287 
288  #elif CT_EMF_AVERAGE == UCT0
289 
290 /* -- Subtract cell-centered contribution -- */
291 
292  EXPAND(vx = d->Vc[VX1]; Bx = d->Vc[BX1]; ,
293  vy = d->Vc[VX2]; By = d->Vc[BX2]; ,
294  vz = d->Vc[VX3]; Bz = d->Vc[BX3];)
295 
296  for (k = emf.kbeg; k <= emf.kend + KOFFSET; k++){
297  for (j = emf.jbeg; j <= emf.jend + JOFFSET; j++){
298  for (i = emf.ibeg; i <= emf.iend + IOFFSET; i++){
299  #if DIMENSIONS == 3
300  emf.exj[k][j][i] *= 2.0;
301  emf.exk[k][j][i] *= 2.0;
302  emf.eyi[k][j][i] *= 2.0;
303  emf.eyk[k][j][i] *= 2.0;
304 
305  emf.exj[k][j][i] -= 0.5*(EX(k,j,i) + EX(k,j+1,i));
306  emf.exk[k][j][i] -= 0.5*(EX(k,j,i) + EX(k+1,j,i));
307 
308  emf.eyi[k][j][i] -= 0.5*(EY(k,j,i) + EY(k,j,i+1));
309  emf.eyk[k][j][i] -= 0.5*(EY(k,j,i) + EY(k+1,j,i));
310  #endif
311  emf.ezi[k][j][i] *= 2.0;
312  emf.ezj[k][j][i] *= 2.0;
313  emf.ezi[k][j][i] -= 0.5*(EZ(k,j,i) + EZ(k,j,i+1));
314  emf.ezj[k][j][i] -= 0.5*(EZ(k,j,i) + EZ(k,j+1,i));
315  }}}
316 
317  CT_EMF_ArithmeticAverage (&emf, 0.25);
318 
319  #else
320  print1 ("! CT_GetEMF: unknown EMF average.\n");
321  QUIT_PLUTO(1);
322  #endif
323 
324 
325 /* ------------------------------------------------------
326  Add contributions from resistive terms with explicit
327  time stepping.
328  ------------------------------------------------------ */
329 
330  #if RESISTIVITY == EXPLICIT
331  CT_AddResistiveEMF(d, grid);
332  #endif
333 
334 /* -------------------------------------------------------------
335  Fine Tuning: EMF_USERDEF_BOUNDARY can be used to directly
336  set the edge-centered electric field
337  ------------------------------------------------------------- */
338 /*
339  #ifdef CTU
340  if (step == 2)
341  #endif
342  {
343  int lside[3] = {X1_BEG, X2_BEG, X3_BEG};
344  int rside[3] = {X1_END, X2_END, X3_END};
345  int dir;
346 
347  for (dir = 0; dir < DIMENSIONS; dir++){
348  if (grid[dir].lbound == USERDEF)
349  EMF_USERDEF_BOUNDARY (&emf, lside[dir], EDGE_EMF, grid);
350  if (grid[dir].rbound == USERDEF)
351  EMF_USERDEF_BOUNDARY (&emf, rside[dir], EDGE_EMF, grid);
352  }
353  }
354 */
355 
356  return (&emf);
357 }
358 
359 #if RESISTIVITY != NO
360 /* ********************************************************************* */
361 void CT_AddResistiveEMF (const Data *d, Grid *grid)
362 /*!
363  * Add resistive terms to EMF.
364  *
365  * \param [in] d
366  * \param [in] grid
367  *
368  * \return a pointer to an edge-centered EMF.
369  *********************************************************************** */
370 {
371  int i, j, k;
372  double ***vx, ***vy, ***vz;
373  double ***Bx, ***By, ***Bz;
374  Data_Arr eta;
375 /*
376  #ifdef CTU
377  if (g_intStage == 2) return (&emf);
378  #endif
379 */
380  eta = GetStaggeredEta();
381 
382  for (k = emf.kbeg; k <= emf.kend; k++){
383  for (j = emf.jbeg; j <= emf.jend; j++){
384  for (i = emf.ibeg; i <= emf.iend; i++){
385 
386 /* -------------------------------------------------------------
387  Simple implementation in Cartesian coordinates by computing
388  currents directly here (bypass GetCurrent).
389  -------------------------------------------------------------
390 
391 double ***Bx = d->Vs[BX1s], Jx;
392 double ***By = d->Vs[BX2s], Jy;
393 double ***Bz = d->Vs[BX3s], Jz;
394 double dx = grid[IDIR].dx[i];
395 double dy = grid[JDIR].dx[j];
396 double dz = grid[KDIR].dx[k];
397 
398 #if DIMENSIONS == 3
399 Jx = (Bz[k][j+1][i] - Bz[k][j][i])/dy - (By[k+1][j][i] - By[k][j][i])/dz;
400 Jy = (Bx[k+1][j][i] - Bx[k][j][i])/dz - (Bz[k][j][i+1] - Bz[k][j][i])/dx;
401 emf.ex[k][j][i] += Jx*g_inputParam[ETAX];
402 emf.ey[k][j][i] += Jy*g_inputParam[ETAY];
403 #endif
404 Jz = (By[k][j][i+1] - By[k][j][i])/dx - (Bx[k][j+1][i] - Bx[k][j][i])/dy;
405 emf.ez[k][j][i] += Jz*g_inputParam[ETAZ];
406 */
407  #if DIMENSIONS == 3
408  emf.ex[k][j][i] += eta[IDIR][k][j][i]*d->J[IDIR][k][j][i];
409  emf.ey[k][j][i] += eta[JDIR][k][j][i]*d->J[JDIR][k][j][i];
410  #endif
411  emf.ez[k][j][i] += eta[KDIR][k][j][i]*d->J[KDIR][k][j][i];
412 
413  }}}
414 }
415 #endif /* RESISTIVITY != NO */
416 
417 #undef EX
418 #undef EY
419 #undef EZ
420 #undef eps_UCT_CONTACT
MAX
#define MAX(a, b)
Definition: macros.h:101
emf
static EMF emf
Definition: ct_emf.c:27
CT_StoreEMF
void CT_StoreEMF(const State_1D *state, int beg, int end, Grid *grid)
Definition: ct_emf.c:35
DATA::Vs
double **** Vs
The main four-index data array used for face-centered staggered magnetic fields.
Definition: structs.h:43
CT_GetStagSlopes
void CT_GetStagSlopes(const Data_Arr b, EMF *emf, Grid *grid)
Definition: ct_stag_slopes.c:68
Data_Arr
double **** Data_Arr
Definition: pluto.h:492
DATA::J
double **** J
Electric current defined as curl(B).
Definition: structs.h:54
VX2
#define VX2
Definition: mod_defs.h:29
STATE_1D::flux
double ** flux
upwind flux computed with the Riemann solver
Definition: structs.h:149
eps_UCT_CONTACT
#define eps_UCT_CONTACT
Definition: ct_emf.c:29
ElectroMotiveForce::SxL
double *** SxL
Definition: ct.h:93
ElectroMotiveForce::SzR
double *** SzR
Definition: ct.h:98
print1
void print1(const char *fmt,...)
Definition: amrPluto.cpp:511
ElectroMotiveForce::ex
double *** ex
Definition: ct.h:103
RHO
#define RHO
Definition: mod_defs.h:19
eta
static double *** eta[3]
Definition: res_functions.c:94
ElectroMotiveForce::svy
signed char *** svy
Definition: ct.h:83
NX2_TOT
long int NX2_TOT
Total number of zones in the X2 direction (boundaries included) for the local processor.
Definition: globals.h:57
ElectroMotiveForce::SxR
double *** SxR
Definition: ct.h:94
g_intStage
int g_intStage
Gives the current integration stage of the time stepping method (predictor = 0, 1st corrector = 1...
Definition: globals.h:98
DATA::Vc
double **** Vc
The main four-index data array used for cell-centered primitive variables.
Definition: structs.h:31
CT_GetEMF
EMF * CT_GetEMF(const Data *d, Grid *grid)
Definition: ct_emf.c:205
STATE_1D::SR
double * SR
Rightmost velocity in the Riemann fan at i+1/2.
Definition: structs.h:167
CT_EMF_ArithmeticAverage
void CT_EMF_ArithmeticAverage(const EMF *Z1, const double w)
Compute arithmetic average of EMF at cell edges.
Definition: ct_emf_average.c:13
UCT_CONTACT
#define UCT_CONTACT
Definition: ct.h:21
ARRAY_3D
#define ARRAY_3D(nx, ny, nz, type)
Definition: prototypes.h:172
EZ
#define EZ(k, j, i)
Definition: ct_emf.c:32
VX1
#define VX1
Definition: mod_defs.h:28
KDIR
#define KDIR
Definition: pluto.h:195
CT_EMF_HLL_Solver
void CT_EMF_HLL_Solver(const Data *d, const EMF *emf, Grid *grid)
Solve 2D Riemann problem for induction equation.
Definition: ct_emf_average.c:203
g_i
int g_i
x1 grid index when sweeping along the x2 or x3 direction.
Definition: globals.h:82
CT_EMF_AVERAGE
#define CT_EMF_AVERAGE
Definition: definitions_01.h:48
PARABOLIC_STEP
#define PARABOLIC_STEP
Definition: pluto.h:74
ElectroMotiveForce::ezi
double *** ezi
Ez available at x-faces (i+1/2);.
Definition: ct.h:79
ElectroMotiveForce::ez
double *** ez
Definition: ct.h:105
TOT_LOOP
#define TOT_LOOP(k, j, i)
Definition: macros.h:44
IDIR
#define IDIR
Definition: pluto.h:193
g_dir
int g_dir
Specifies the current sweep or direction of integration.
Definition: globals.h:86
ElectroMotiveForce::exk
double *** exk
Ex available at z-faces (k+1/2);.
Definition: ct.h:76
g_j
int g_j
x2 grid index when sweeping along the x1 or x3 direction.
Definition: globals.h:83
GRID
Definition: structs.h:78
Bx
static double Bx
Definition: hlld.c:62
STATE_1D::SL
double * SL
Leftmost velocity in the Riemann fan at i+1/2.
Definition: structs.h:166
NX3_TOT
long int NX3_TOT
Total number of zones in the X3 direction (boundaries included) for the local processor.
Definition: globals.h:59
j
int j
Definition: analysis.c:2
k
int k
Definition: analysis.c:2
CT_EMF_CMUSCL_Average
void CT_EMF_CMUSCL_Average(const Data *d, const EMF *emf, Grid *grid)
–
Definition: ct_emf_average.c:391
UCT_HLL
#define UCT_HLL
Definition: ct.h:22
ElectroMotiveForce::SyR
double *** SyR
Definition: ct.h:96
g_k
int g_k
x3 grid index when sweeping along the x1 or x2 direction.
Definition: globals.h:84
ElectroMotiveForce::svz
signed char *** svz
Definition: ct.h:83
s
#define s
BX3
#define BX3
Definition: mod_defs.h:27
pluto.h
PLUTO main header file.
ElectroMotiveForce::eyi
double *** eyi
Ey available at x-faces (i+1/2);.
Definition: ct.h:77
ElectroMotiveForce::svx
signed char *** svx
Definition: ct.h:83
CT_StoreVelSlopes
void CT_StoreVelSlopes(EMF *emf, const State_1D *state, int beg, int end)
Definition: ct_stag_slopes.c:5
CT_EMF_IntegrateToCorner
void CT_EMF_IntegrateToCorner(const Data *d, const EMF *emf, Grid *grid)
Definition: ct_emf_average.c:58
D_EXPAND
D_EXPAND(tot/[n]=(double) grid[IDIR].np_int_glob;, tot/[n]=(double) grid[JDIR].np_int_glob;, tot/[n]=(double) grid[KDIR].np_int_glob;)
Definition: analysis.c:27
ElectroMotiveForce::eyk
double *** eyk
Ey available at z-faces (k+1/2);.
Definition: ct.h:78
DATA
Definition: structs.h:30
EY
#define EY(k, j, i)
Definition: ct_emf.c:31
STATE_1D
Definition: structs.h:133
g_operatorStep
int g_operatorStep
Gives the current operator step.
Definition: globals.h:101
i
int i
Definition: analysis.c:2
ElectroMotiveForce::SzL
double *** SzL
Definition: ct.h:97
ElectroMotiveForce::ezj
double *** ezj
Ez available at y-faces (j+1/2);.
Definition: ct.h:80
EX
#define EX(k, j, i)
Definition: ct_emf.c:30
BX1
#define BX1
Definition: mod_defs.h:25
VX3
#define VX3
Definition: mod_defs.h:30
BX2
#define BX2
Definition: mod_defs.h:26
ElectroMotiveForce::ey
double *** ey
Definition: ct.h:104
JDIR
#define JDIR
Definition: pluto.h:194
ElectroMotiveForce::exj
double *** exj
Ex available at y-faces (j+1/2);.
Definition: ct.h:75
QUIT_PLUTO
#define QUIT_PLUTO(e_code)
Definition: macros.h:125
NX1_TOT
long int NX1_TOT
Total number of zones in the X1 direction (boundaries included) for the local processor.
Definition: globals.h:55
ElectroMotiveForce::SyL
double *** SyL
Definition: ct.h:95