LORENE
strot_dirac_equilibrium.C
1 /*
2  * Function Star_rot_Dirac::equilibrium
3  *
4  * (see file star_rot_dirac.h for documentation).
5  *
6  */
7 
8 /*
9  * Copyright (c) 2005 Lap-Ming Lin & Jerome Novak
10  *
11  * This file is part of LORENE.
12  *
13  * LORENE is free software; you can redistribute it and/or modify
14  * it under the terms of the GNU General Public License version 2
15  * as published by the Free Software Foundation.
16  *
17  * LORENE is distributed in the hope that it will be useful,
18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20  * GNU General Public License for more details.
21  *
22  * You should have received a copy of the GNU General Public License
23  * along with LORENE; if not, write to the Free Software
24  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25  *
26  */
27 
28 
29 
30 /*
31  * $Id: strot_dirac_equilibrium.C,v 1.15 2023/12/20 10:16:18 j_novak Exp $
32  * $Log: strot_dirac_equilibrium.C,v $
33  * Revision 1.15 2023/12/20 10:16:18 j_novak
34  * New options in Star_rot_Dirac/rotstar_dirac: a flag to choose between GR, CFC and CFC+ A^{ij}_{TT}=0 (as described n Cordero_Carrion et al. 2009).
35  *
36  * Revision 1.14 2016/12/05 16:18:15 j_novak
37  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
38  *
39  * Revision 1.13 2014/10/13 08:53:40 j_novak
40  * Lorene classes and functions now belong to the namespace Lorene.
41  *
42  * Revision 1.12 2014/10/06 15:13:18 j_novak
43  * Modified #include directives to use c++ syntax.
44  *
45  * Revision 1.11 2009/10/26 10:54:33 j_novak
46  * Added the case of a NONSYM base in theta.
47  *
48  * Revision 1.10 2008/02/25 10:40:52 j_novak
49  * Added the flag mer_hij to control the step from which the equation for h^{ij}
50  * is being solved.
51  *
52  * Revision 1.9 2006/03/14 15:18:21 lm_lin
53  *
54  * Add convergence to a given baryon mass.
55  *
56  * Revision 1.8 2005/11/28 14:45:16 j_novak
57  * Improved solution of the Poisson tensor equation in the case of a transverse
58  * tensor.
59  *
60  * Revision 1.7 2005/09/16 14:04:49 j_novak
61  * The equation for hij is now solved only for mer > mer_fix_omega. It uses the
62  * Poisson solver of the class Sym_tensor_trans.
63  *
64  * Revision 1.6 2005/04/20 14:26:29 j_novak
65  * Removed some outputs.
66  *
67  * Revision 1.5 2005/04/05 09:24:05 j_novak
68  * minor modifs
69  *
70  * Revision 1.4 2005/03/10 09:39:19 j_novak
71  * The order of resolution has been changed in the iteration step.
72  *
73  * Revision 1.3 2005/02/17 17:30:09 f_limousin
74  * Change the name of some quantities to be consistent with other classes
75  * (for instance nnn is changed to nn, shift to beta, beta to lnq...)
76  *
77  * Revision 1.2 2005/02/09 13:36:42 lm_lin
78  *
79  * Calculate GRV2 during iterations.
80  *
81  * Revision 1.1 2005/01/31 08:51:48 j_novak
82  * New files for rotating stars in Dirac gauge (still under developement).
83  *
84  *
85  * $Header: /cvsroot/Lorene/C++/Source/Star/strot_dirac_equilibrium.C,v 1.15 2023/12/20 10:16:18 j_novak Exp $
86  *
87  */
88 
89 
90 // C headers
91 #include <cmath>
92 #include <cassert>
93 
94 // Lorene headers
95 #include "star_rot_dirac.h"
96 
97 #include "utilitaires.h"
98 #include "unites.h"
99 
100 namespace Lorene {
101 void Star_rot_Dirac::equilibrium(double ent_c, double omega0,
102  double fact_omega, int , const Tbl& ent_limit,
103  const Itbl& icontrol, const Tbl& control,
104  double mbar_wanted, double aexp_mass, Tbl& diff){
105 
106 
107  // Fundamental constants and units
108  // --------------------------------
109  using namespace Unites ;
110 
111  // For the display
112  // ---------------
113  char display_bold[]="x[1m" ; display_bold[0] = 27 ;
114  char display_normal[] = "x[0m" ; display_normal[0] = 27 ;
115 
116 
117  // Grid parameters
118  // ----------------
119 
120  const Mg3d* mg = mp.get_mg() ;
121  int nz = mg->get_nzone() ; // total number of domains
122  int nzm1 = nz - 1 ;
123 
124  // Index of the point at phi=0, theta=pi/2 at the surface of the star:
125  int type_t = mg->get_type_t() ;
126  assert( ( type_t == SYM) || (type_t == NONSYM) ) ;
127  int l_b = nzet - 1 ;
128  int i_b = mg->get_nr(l_b) - 1 ;
129  int j_b = (type_t == SYM ? mg->get_nt(l_b) - 1 : mg->get_nt(l_b)/2 ) ;
130  int k_b = 0 ;
131 
132  // Value of the enthalpy defining the surface of the star
133  double ent_b = ent_limit(nzet-1) ;
134 
135  // Parameters to control the iteration
136  // -----------------------------------
137 
138  int mer_max = icontrol(0) ;
139  int mer_rot = icontrol(1) ;
140  int mer_change_omega = icontrol(2) ;
141  int mer_fix_omega = icontrol(3) ;
142  int mer_mass = icontrol(4) ;
143  int delta_mer_kep = icontrol(5) ;
144  int mer_hij = icontrol(6) ;
145 
146  // Protections:
147  if (mer_change_omega < mer_rot) {
148  cout << "Star_rot_Dirac::equilibrium: mer_change_omega < mer_rot !"
149  << endl ;
150  cout << " mer_change_omega = " << mer_change_omega << endl ;
151  cout << " mer_rot = " << mer_rot << endl ;
152  abort() ;
153  }
154  if (mer_fix_omega < mer_change_omega) {
155  cout << "Star_rot_Dirac::equilibrium: mer_fix_omega < mer_change_omega !"
156  << endl ;
157  cout << " mer_fix_omega = " << mer_fix_omega << endl ;
158  cout << " mer_change_omega = " << mer_change_omega << endl ;
159  abort() ;
160  }
161 
162  // In order to converge to a given baryon mass, shall the central
163  // enthalpy be varied or Omega ?
164  bool change_ent = true ;
165  if (mer_mass < 0) {
166  change_ent = false ;
167  mer_mass = abs(mer_mass) ;
168  }
169 
170 
171  double precis = control(0) ;
172  double omega_ini = control(1) ;
173  double relax = control(2) ;
174  double relax_prev = double(1) - relax ;
175 
176  // Error indicators
177  // ----------------
178 
179  diff.annule_hard() ;
180  double& diff_ent = diff.set(0) ;
181 
182  double alpha_r = 1 ;
183 
184  // Initializations
185  // ---------------
186 
187  // Initial angular velocities
188  omega = 0 ;
189 
190  double accrois_omega = (omega0 - omega_ini) /
191  double(mer_fix_omega - mer_change_omega) ;
192 
193 
194  update_metric() ; //update of the metric quantities
195 
196  equation_of_state() ; // update of the density, pressure,...etc
197 
198  hydro_euler() ; //update of the hydro quantities relative to the
199  // Eulerian observer
200 
201  // Quantities at the previous step :
202  Scalar ent_prev = ent ;
203  Scalar logn_prev = logn ;
204  Scalar qqq_prev = qqq ;
205  // Vector beta_prev = beta ;
206  // Sym_tensor_trans hh_prev = hh ;
207 
208  // Output files
209  // -------------
210 
211  ofstream fichconv("convergence.d") ; // Output file for diff_ent
212  fichconv << "# diff_ent GRV2 max_triax vit_triax" << endl ;
213 
214  ofstream fichfreq("frequency.d") ; // Output file for omega
215  fichfreq << "# f [Hz]" << endl ;
216 
217  ofstream fichevol("evolution.d") ; // Output file for various quantities
218  fichevol <<
219  "# |dH/dr_eq/dH/dr_pole| r_pole/r_eq ent_c"
220  << endl ;
221 
222  diff_ent = 1 ;
223  double err_grv2 = 1 ;
224 
225 
226 
227 //=========================================================================
228 // Start of iteration
229 //=========================================================================
230 
231  for(int mer=0 ; (diff_ent > precis) && (mer<mer_max) ; mer++) {
232 
233  cout << "-----------------------------------------------" << endl ;
234  cout << "step: " << mer << endl ;
235  cout << "diff_ent = " << display_bold << diff_ent << display_normal
236  << endl ;
237  cout << "err_grv2 = " << err_grv2 << endl ;
238  fichconv << mer ;
239  fichfreq << mer ;
240  fichevol << mer ;
241 
242 
243  // switch on rotation
244  if (mer >= mer_rot) {
245 
246  if (mer < mer_change_omega) {
247  omega = omega_ini ;
248  }
249  else {
250  if (mer <= mer_fix_omega) {
251  omega = omega_ini + accrois_omega *
252  (mer - mer_change_omega) ;
253  }
254  }
255 
256 
257  }
258 
259 
260  //---------------------------------------------------//
261  // Resolution of the Poisson equation for logn //
262  // Note: ln_f is due to the fluid part //
263  // ln_q is due to the quadratic metric part //
264  //---------------------------------------------------//
265 
266  Scalar ln_f_new(mp) ;
267  Scalar ln_q_new(mp) ;
268 
269  solve_logn_f( ln_f_new ) ;
270  solve_logn_q( ln_q_new ) ;
271 
272  ln_f_new.std_spectral_base() ;
273  ln_q_new.std_spectral_base() ;
274 
275 
276  //--------------------------------------------------//
277  // Resolution of the Poisson equation for shift //
278  //--------------------------------------------------//
279 
280  Vector beta_new(mp, CON, mp.get_bvect_spher()) ;
281 
282  solve_shift( beta_new ) ;
283 
284  //------------------------------------
285  // Determination of the fluid velocity
286  //------------------------------------
287 
288  if (mer > mer_fix_omega + delta_mer_kep) {
289 
290  omega *= fact_omega ; // Increase of the angular velocity if
291  } // fact_omega != 1
292 
293  bool omega_trop_grand = false ;
294  bool kepler = true ;
295 
296  while ( kepler ) {
297 
298  // Possible decrease of Omega to ensure a velocity < c
299 
300  bool superlum = true ;
301 
302  while ( superlum ){
303 
304  // New fluid velocity :
305  //
306 
307  u_euler.set(1).set_etat_zero() ;
308  u_euler.set(2).set_etat_zero() ;
309 
310  u_euler.set(3) = omega ;
311  u_euler.set(3).annule(nzet,nzm1) ; // nzet is defined in class Star
313  u_euler.set(3).mult_rsint() ;
314  u_euler.set(3) += beta(3) ;
315  u_euler.set(3).annule(nzet,nzm1) ;
316 
317  u_euler = u_euler / nn ;
318 
319 
320  // v2 (square of norm of u_euler)
321  // -------------------------------
322 
323  v2 = contract(contract(gamma.cov(), 0, u_euler, 0), 0, u_euler, 0) ;
324 
325  // Is the new velocity larger than c in the equatorial plane ?
326 
327  superlum = false ;
328 
329  for (int l=0; l<nzet; l++) {
330  for (int i=0; i<mg->get_nr(l); i++) {
331 
332  double u2 = v2.val_grid_point(l, 0, j_b, i) ;
333  if (u2 >= 1.) { // superluminal velocity
334  superlum = true ;
335  cout << "U > c for l, i : " << l << " " << i
336  << " U = " << sqrt( u2 ) << endl ;
337  }
338  }
339  }
340  if ( superlum ) {
341  cout << "**** VELOCITY OF LIGHT REACHED ****" << endl ;
342  omega /= fact_omega ; // Decrease of Omega
343  cout << "New rotation frequency : "
344  << omega/(2.*M_PI) * f_unit << " Hz" << endl ;
345  omega_trop_grand = true ;
346  }
347  } // end of while ( superlum )
348 
349 
350  // New computation of U (this time is not superluminal)
351  // as well as of gam_euler, ener_euler,...etc
352 
353  hydro_euler() ;
354 
355 
356 
357  //--------------------------------//
358  // First integral of motion //
359  //--------------------------------//
360 
361  Scalar mlngamma(mp) ; // -log( gam_euler )
362 
363  mlngamma = - log( gam_euler ) ;
364 
365  // Equatorial values of various potentials :
366  double ln_f_b = ln_f_new.val_grid_point(l_b, k_b, j_b, i_b) ;
367  double ln_q_b = ln_q_new.val_grid_point(l_b, k_b, j_b, i_b) ;
368  double mlngamma_b = mlngamma.val_grid_point(l_b, k_b, j_b, i_b) ;
369 
370  // Central values of various potentials :
371  double ln_f_c = ln_f_new.val_grid_point(0,0,0,0) ;
372  double ln_q_c = ln_q_new.val_grid_point(0,0,0,0) ;
373  double mlngamma_c = 0 ;
374 
375  // Scale factor to ensure that the (log of) enthalpy is equal to
376  // ent_b at the equator
377  double alpha_r2 = ( ent_c - ent_b + mlngamma_c - mlngamma_b
378  + ln_q_c - ln_q_b) / ( ln_f_b - ln_f_c ) ;
379 
380  alpha_r = sqrt(alpha_r2) ;
381 
382  cout << "alpha_r = " << alpha_r << endl ;
383 
384  // Rescaling of the grid (no adaptation!)
385  //---------------------------------------
386  mp.homothetie(alpha_r) ;
387 
388  // Readjustment of logn :
389  // -----------------------
390 
391  logn = alpha_r2 * ln_f_new + ln_q_new ;
392 
393  double logn_c = logn.val_grid_point(0,0,0,0) ;
394 
395  // First integral of motion -> (log of) enthalpy in all space
396  // ----------------------------------------------------------
397 
398  ent = (ent_c + logn_c + mlngamma_c) - logn - mlngamma ;
399 
400 
401  // --------------------------------------------------------------
402  // Test: is the enthalpy negative somewhere in the equatorial plane
403  // inside the star?
404  // --------------------------------------------------------
405 
406  kepler = false ;
407  for (int l=0; l<nzet; l++) {
408  int imax = mg->get_nr(l) - 1 ;
409  if (l == l_b) imax-- ; // The surface point is skipped
410  for (int i=0; i<imax; i++) {
411  if ( ent.val_grid_point(l, 0, j_b, i) < 0. ) {
412  kepler = true ;
413  cout << "ent < 0 for l, i : " << l << " " << i
414  << " ent = " << ent.val_grid_point(l, 0, j_b, i) << endl ;
415  }
416  }
417  }
418 
419  if ( kepler ) {
420  cout << "**** KEPLERIAN VELOCITY REACHED ****" << endl ;
421  omega /= fact_omega ; // Omega is decreased
422  cout << "New rotation frequency : "
423  << omega/(2.*M_PI) * f_unit << " Hz" << endl ;
424  omega_trop_grand = true ;
425  }
426 
427  } // End of while ( kepler )
428 
429  if ( omega_trop_grand ) { // fact_omega is decreased for the
430  // next step
431  fact_omega = sqrt( fact_omega ) ;
432  cout << "**** New fact_omega : " << fact_omega << endl ;
433  }
434 
435 
436 //---------------------------------
437  // Equation of state
438  //---------------------------------
439 
440  equation_of_state() ; // computes new values for nbar (n), ener (e),
441  // and press (p) from the new ent (H)
442 
443  hydro_euler() ;
444 
445  //---------------------------------------------//
446  // Resolution of the Poisson equation for qqq //
447  //---------------------------------------------//
448 
449  Scalar q_new(mp) ;
450 
451  solve_qqq( q_new ) ;
452 
453  q_new.std_spectral_base() ;
454 
455  //----------------------------------------------//
456  // Resolution of the Poisson equation for hh //
457  //----------------------------------------------//
458 
459  Sym_tensor_trans hij_new(mp, mp.get_bvect_spher(), flat) ;
460 
461  if ((mer > mer_hij) && (relat_type == 1))
462  solve_hij( hij_new ) ;
463  else
464  hij_new.set_etat_zero() ;
465 
466  hh = hij_new ;
467  qqq = q_new ;
468  beta = beta_new ;
469 
470  //---------------------------------------
471  // Calculate error of the GRV2 identity
472  //---------------------------------------
473 
474  err_grv2 = grv2() ;
475 
476 
477  //--------------------------------------
478  // Relaxations on some quantities....?
479  //
480  // ** On logn and qqq?
481  //--------------------------------------
482 
483  if (mer >= 10) {
484  logn = relax * logn + relax_prev * logn_prev ;
485 
486  qqq = relax * qqq + relax_prev * qqq_prev ;
487 
488  }
489 
490  // Update of the metric quantities :
491 
492  update_metric() ;
493 
494  //---------------------------
495  // Informations display
496  // More to come later......
497  //---------------------------
498 
499  // partial_display(cout) ; // What is partial_display(cout) ?
500  fichfreq << " " << omega / (2*M_PI) * f_unit ;
501  fichevol << " " << ent_c ;
502 
503 
504  //-----------------------------------------
505  // Convergence towards a given baryon mass
506  //-----------------------------------------
507 
508  if (mer > mer_mass) {
509 
510  double xx ;
511  if (mbar_wanted > 0.) {
512  xx = mass_b() / mbar_wanted - 1. ;
513  cout << "Discrep. baryon mass <-> wanted bar. mass : " << xx
514  << endl ;
515  }
516  else{
517  xx = mass_g() / fabs(mbar_wanted) - 1. ;
518  cout << "Discrep. grav. mass <-> wanted grav. mass : " << xx
519  << endl ;
520  }
521  double xprog = ( mer > 2*mer_mass) ? 1. :
522  double(mer-mer_mass)/double(mer_mass) ;
523  xx *= xprog ;
524  double ax = .5 * ( 2. + xx ) / (1. + xx ) ;
525  double fact = pow(ax, aexp_mass) ;
526  cout << " xprog, xx, ax, fact : " << xprog << " " <<
527  xx << " " << ax << " " << fact << endl ;
528 
529  if ( change_ent ) {
530  ent_c *= fact ;
531  }
532  else {
533  if (mer%4 == 0) omega *= fact ;
534  }
535  }
536 
537 
538  //-----------------------------------------------------------
539  // Relative change in enthalpy with respect to previous step
540  // ** Check: Is diffrel(ent, ent_prev) ok?
541  //-----------------------------------------------------------
542 
543  Tbl diff_ent_tbl = diffrel( ent, ent_prev ) ;
544  diff_ent = diff_ent_tbl(0) ;
545  for (int l=1; l<nzet; l++) {
546  diff_ent += diff_ent_tbl(l) ;
547  }
548  diff_ent /= nzet ;
549 
550  fichconv << " " << log10( fabs(diff_ent) + 1.e-16 ) ;
551  fichconv << " " << log10( fabs(err_grv2) + 1.e-16 ) ;
552 
553  //------------------------------
554  // Recycling for the next step
555  //------------------------------
556 
557  ent_prev = ent ;
558  logn_prev = logn ;
559  qqq_prev = qqq ;
560 
561  fichconv << endl ;
562  fichfreq << endl ;
563  fichevol << endl ;
564  fichconv.flush() ;
565  fichfreq.flush() ;
566  fichevol.flush() ;
567 
568 
569  } // End of main loop
570 
571  //=================================================
572  // End of iteration
573  //=================================================
574 
575  fichconv.close() ;
576  fichfreq.close() ;
577  fichevol.close() ;
578 
579 
580 }
581 }
void solve_logn_f(Scalar &ln_f_new) const
Solution of the two scalar Poisson equations for rotating stars in Dirac gauge.
Cmp log(const Cmp &)
Neperian logarithm.
Definition: cmp_math.C:299
double omega
Rotation angular velocity ([f_unit] )
Map & mp
Mapping associated with the star.
Definition: star.h:180
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:223
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: scalar.C:330
virtual void annule(int l_min, int l_max)
Sets the Scalar to zero in several domains.
Definition: scalar.C:397
Metric gamma
3-metric
Definition: star.h:235
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition: map.h:801
Lorene prototypes.
Definition: app_hor.h:67
Standard units of space, time and mass.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:783
double & set(int i)
Read/write of a particular element (index i) (1D case)
Definition: tbl.h:301
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
Sym_tensor_trans hh
is defined by .
int relat_type
Relativistic flag.
int get_type_t() const
Returns the type of sampling in the direction: SYM : : symmetry with respect to the equatorial pl...
Definition: grilles.h:502
Basic integer array class.
Definition: itbl.h:122
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field. ...
Definition: scalar.C:790
Scalar ent
Log-enthalpy.
Definition: star.h:190
Vector beta
Shift vector.
Definition: star.h:228
Tensor field of valence 1.
Definition: vector.h:188
void update_metric()
Computes metric quantities from known potentials.
Tbl diffrel(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (norme version).
Definition: cmp_math.C:507
virtual void hydro_euler()
Computes the hydrodynamical quantities relative to the Eulerian observer from those in the fluid fram...
double val_grid_point(int l, int k, int j, int i) const
Returns the value of the field at a specified grid point.
Definition: scalar.h:643
int nzet
Number of domains of *mp occupied by the star.
Definition: star.h:183
Scalar gam_euler
Lorentz factor between the fluid and Eulerian observers.
Definition: star.h:204
virtual double mass_g() const
Gravitational mass.
virtual void homothetie(double lambda)=0
Sets a new radial scale.
virtual double grv2() const
Error on the virial identity GRV2.
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
Vector u_euler
Fluid 3-velocity with respect to the Eulerian observer.
Definition: star.h:207
void mult_rsint()
Multiplication by everywhere; dzpuis is not changed.
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Scalar logn
Logarithm of the lapse N .
Definition: star.h:222
virtual double mass_b() const
Baryonic mass.
Transverse symmetric tensors of rank 2.
Definition: sym_tensor.h:611
void equilibrium(double ent_c, double omega0, double fact_omega, int nzadapt, const Tbl &ent_limit, const Itbl &icontrol, const Tbl &control, double mbar_wanted, double aexp_mass, Tbl &diff)
Computes an equilibrium configuration.
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:469
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Definition: metric.C:283
Multi-domain grid.
Definition: grilles.h:279
void solve_hij(Sym_tensor_trans &hij_new) const
Solution of the tensor Poisson equation for rotating stars in Dirac gauge.
Scalar nn
Lapse function N .
Definition: star.h:225
Cmp log10(const Cmp &)
Basis 10 logarithm.
Definition: cmp_math.C:325
Cmp abs(const Cmp &)
Absolute value.
Definition: cmp_math.C:413
void solve_shift(Vector &shift_new) const
Solution of the shift equation for rotating stars in Dirac gauge.
void solve_logn_q(Scalar &ln_q_new) const
Solution of the two scalar Poisson equations for rotating stars in Dirac gauge.
void solve_qqq(Scalar &q_new) const
Solution of the two scalar Poisson equations for rotating stars in Dirac gauge.
Basic array class.
Definition: tbl.h:164
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:474
Scalar & set(int)
Read/write access to a component.
Definition: vector.C:302
void equation_of_state()
Computes the proper baryon and energy density, as well as pressure from the enthalpy.
Definition: star.C:465
const Metric_flat & flat
flat metric (spherical components)
void annule_hard()
Sets the Tbl to zero in a hard way.
Definition: tbl.C:375