LORENE
et_bfrot_global.C
1 /*
2  * Copyright (c) 2001 Jerome Novak
3  *
4  * This file is part of LORENE.
5  *
6  * LORENE is free software; you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation; either version 2 of the License, or
9  * (at your option) any later version.
10  *
11  * LORENE is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with LORENE; if not, write to the Free Software
18  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19  *
20  */
21 
22 
23 char et_bfrot_global_C[] = "$Header: /cvsroot/Lorene/C++/Source/Etoile/et_bfrot_global.C,v 1.17 2015/06/10 14:39:17 a_sourie Exp $" ;
24 
25 /*
26  * $Id: et_bfrot_global.C,v 1.17 2015/06/10 14:39:17 a_sourie Exp $
27  * $Log: et_bfrot_global.C,v $
28  * Revision 1.17 2015/06/10 14:39:17 a_sourie
29  * New class Eos_bf_tabul for tabulated 2-fluid EoSs and associated functions for the computation of rotating stars with such EoSs.
30  *
31  * Revision 1.16 2014/10/13 08:52:54 j_novak
32  * Lorene classes and functions now belong to the namespace Lorene.
33  *
34  * Revision 1.15 2014/10/06 15:13:07 j_novak
35  * Modified #include directives to use c++ syntax.
36  *
37  * Revision 1.14 2004/09/03 13:55:07 j_novak
38  * Use of enerps_euler instead of ener_euler in the calculation of mom_angu()
39  *
40  * Revision 1.13 2004/03/25 10:29:03 j_novak
41  * All LORENE's units are now defined in the namespace Unites (in file unites.h).
42  *
43  * Revision 1.12 2003/11/20 14:01:26 r_prix
44  * changed member names to better conform to Lorene coding standards:
45  * J_euler -> j_euler, EpS_euler -> enerps_euler, Delta_car -> delta_car
46  *
47  * Revision 1.11 2003/11/18 18:38:11 r_prix
48  * use of new member EpS_euler: matter sources in equilibrium() and global quantities
49  * no longer distinguish Newtonian/relativistic, as all terms should have the right limit...
50  *
51  * Revision 1.10 2003/11/13 12:13:15 r_prix
52  * *) removed all uses of Etoile-type specific quantities: u_euler, press
53  * and exclusively use ener_euler, s_euler, sphph_euler, J_euler instead
54  * *) this also fixes a bug in previous expression for mass_g() in 2-fluid case
55  *
56  * NOTE: some current Newtonian expressions are still completely broken..
57  *
58  * Revision 1.9 2003/09/17 08:27:50 j_novak
59  * New methods: mass_b1() and mass_b2().
60  *
61  * Revision 1.8 2003/02/07 10:47:43 j_novak
62  * The possibility of having gamma5 xor gamma6 =0 has been introduced for
63  * tests. The corresponding parameter files have been added.
64  *
65  * Revision 1.7 2002/10/16 14:36:36 j_novak
66  * Reorganization of #include instructions of standard C++, in order to
67  * use experimental version 3 of gcc.
68  *
69  * Revision 1.6 2002/10/14 14:20:08 j_novak
70  * Error corrected for angu_mom()
71  *
72  * Revision 1.5 2002/05/10 09:26:52 j_novak
73  * Added new class Et_rot_mag for magnetized rotating neutron stars (under development)
74  *
75  * Revision 1.4 2002/04/05 09:09:36 j_novak
76  * The inversion of the EOS for 2-fluids polytrope has been modified.
77  * Some errors in the determination of the surface were corrected.
78  *
79  * Revision 1.3 2002/01/08 14:43:53 j_novak
80  * better determination of surfaces for 2-fluid stars
81  *
82  * Revision 1.2 2002/01/03 15:30:28 j_novak
83  * Some comments modified.
84  *
85  * Revision 1.1.1.1 2001/11/20 15:19:28 e_gourgoulhon
86  * LORENE
87  *
88  * Revision 1.4 2001/08/31 15:07:12 novak
89  * Retour a la version 1.2, sans la routine prolonge_c1
90  *
91  * Revision 1.2 2001/08/28 15:32:17 novak
92  * The surface is now defined by the baryonic density
93  *
94  * Revision 1.1 2001/06/22 15:39:46 novak
95  * Initial revision
96  *
97  *
98  * $Header: /cvsroot/Lorene/C++/Source/Etoile/et_bfrot_global.C,v 1.17 2015/06/10 14:39:17 a_sourie Exp $
99  *
100  */
101 
102 
103 // Headers C
104 #include <cstdlib>
105 #include <cmath>
106 
107 // Headers Lorene
108 #include "et_rot_bifluid.h"
109 #include "unites.h"
110 
111 //--------------------------//
112 // Baryon mass //
113 //--------------------------//
114 
115 namespace Lorene {
116 double Et_rot_bifluid::mass_b1() const {
117 
118  if (p_mass_b1 == 0x0) { // a new computation is required
119 
120  assert(nbar.get_etat() == ETATQCQ);
121  Cmp dens1 = a_car() * bbb() * gam_euler() * eos.get_m1()* nbar();
122  dens1.std_base_scal() ;
123 
124  p_mass_b1 = new double( dens1.integrale() ) ;
125 
126  }
127 
128  return *p_mass_b1 ;
129 
130 }
131 
132 double Et_rot_bifluid::mass_b2() const {
133 
134  if (p_mass_b2 == 0x0) { // a new computation is required
135  assert(nbar2.get_etat() == ETATQCQ);
136 
137  Cmp dens2 = a_car() * bbb() * gam_euler2() * eos.get_m2() * nbar2();
138  dens2.std_base_scal() ;
139 
140  p_mass_b2 = new double( dens2.integrale() ) ;
141 
142  }
143 
144  return *p_mass_b2 ;
145 
146 }
147 
148 double Et_rot_bifluid::mass_b() const {
149 
150  if (p_mass_b == 0x0)
151  p_mass_b = new double(mass_b1() + mass_b2() ) ;
152 
153  return *p_mass_b;
154 
155 }
156 
157 
158 //----------------------------//
159 // Gravitational mass //
160 //----------------------------//
161 
162 double Et_rot_bifluid::mass_g() const {
163 
164  if (p_mass_g == 0x0) { // a new computation is required
165 
166  Tenseur vtmp (j_euler);
167  vtmp.change_triad( mp.get_bvect_spher() ); // change to spherical triad
168 
169  Tenseur tJphi (mp);
170  tJphi = vtmp(2); // get the phi tetrad-component: J^ph r sin(th)
171 
172  // relativistic AND Newtonian: enerps_euler has right limit and tnphi->0
173  Tenseur source = nnn * enerps_euler + 2 * b_car * tnphi * tJphi;
174  source = a_car * bbb * source ;
175 
176  source.set_std_base() ;
177 
178  p_mass_g = new double( source().integrale() ) ;
179 
180  }
181 
182  return *p_mass_g ;
183 
184 }
185 
186 //----------------------------//
187 // Angular momentum //
188 //----------------------------//
189 
190 double Et_rot_bifluid::angu_mom() const {
191 
192  if (p_angu_mom == 0x0) { // a new computation is required
193 
194  Cmp dens(mp) ;
195 
196  // this should work for both relativistic AND Newtonian,
197  // provided j_euler has the right limit...
198  Tenseur tmp = j_euler;
199  tmp.change_triad( mp.get_bvect_spher() ) ;
200  dens = tmp(2) ;
201  dens.mult_rsint() ;
202 
203  dens = a_car() * b_car()* bbb() * dens ;
204 
205  dens.std_base_scal() ;
206 
207  p_angu_mom = new double( dens.integrale() ) ;
208 
209  }
210 
211  return *p_angu_mom ;
212 
213 }
214 
215 
217 
218  if (p_angu_mom_1 == 0x0) { // a new computation is required
219 
220  Cmp dens(mp) ;
221 
222  // this should work for both relativistic AND Newtonian,
223  // provided j_euler has the right limit...
224  Tenseur tmp = j_euler1;
225  tmp.change_triad( mp.get_bvect_spher() ) ;
226  dens = tmp(2) ;
227  dens.mult_rsint() ;
228 
229  dens = a_car() * b_car()* bbb() * dens ;
230 
231  dens.std_base_scal() ;
232 
233  p_angu_mom_1 = new double( dens.integrale() ) ;
234 
235  }
236 
237  return *p_angu_mom_1 ;
238 
239 }
240 
242 
243  if (p_angu_mom_2 == 0x0) { // a new computation is required
244 
245  Cmp dens(mp) ;
246 
247  // this should work for both relativistic AND Newtonian,
248  // provided j_euler has the right limit...
249  Tenseur tmp = j_euler2;
250  tmp.change_triad( mp.get_bvect_spher() ) ;
251  dens = tmp(2) ;
252  dens.mult_rsint() ;
253 
254  dens = a_car() * b_car()* bbb() * dens ;
255 
256  dens.std_base_scal() ;
257 
258  p_angu_mom_2 = new double( dens.integrale() ) ;
259 
260  }
261 
262  return *p_angu_mom_2 ;
263 
264 }
265 
267 
268  if (p_angu_mom_1_part1_1 == 0x0) { // a new computation is required
269 
270  Cmp dens(mp) ;
271 
272  // this should work for both relativistic AND Newtonian,
273  // provided j_euler has the right limit...
274  Tenseur tmp = j_euler11_1;
275  tmp.change_triad( mp.get_bvect_spher() ) ;
276  dens = tmp(2) ;
277  dens.mult_rsint() ;
278 
279  dens = a_car() * b_car()* bbb() * dens ;
280 
281  dens.std_base_scal() ;
282 
283  p_angu_mom_1_part1_1 = new double( dens.integrale() ) ;
284 
285  }
286 
287  return *p_angu_mom_1_part1_1 ;
288 
289 }
290 
292 
293  if (p_angu_mom_2_part1_1 == 0x0) { // a new computation is required
294 
295  Cmp dens(mp) ;
296 
297  // this should work for both relativistic AND Newtonian,
298  // provided j_euler has the right limit...
299  Tenseur tmp = j_euler21_1;
300  tmp.change_triad( mp.get_bvect_spher() ) ;
301  dens = tmp(2) ;
302  dens.mult_rsint() ;
303 
304  dens = a_car() * b_car()* bbb() * dens ;
305 
306  dens.std_base_scal() ;
307 
308  p_angu_mom_2_part1_1 = new double( dens.integrale() ) ;
309 
310  }
311 
312  return *p_angu_mom_2_part1_1 ;
313 
314 }
315 
317 
318  if (p_angu_mom_1_part2_1 == 0x0) { // a new computation is required
319 
320  Cmp dens(mp) ;
321 
322  // this should work for both relativistic AND Newtonian,
323  // provided j_euler has the right limit...
324  Tenseur tmp = j_euler12_1;
325  tmp.change_triad( mp.get_bvect_spher() ) ;
326  dens = tmp(2) ;
327  dens.mult_rsint() ;
328 
329  dens = a_car() * b_car()* bbb() * dens ;
330 
331  dens.std_base_scal() ;
332 
333  p_angu_mom_1_part2_1 = new double( dens.integrale() ) ;
334 
335  }
336 
337  return *p_angu_mom_1_part2_1 ;
338 
339 }
340 
342 
343  if (p_angu_mom_2_part2_1 == 0x0) { // a new computation is required
344 
345  Cmp dens(mp) ;
346 
347  // this should work for both relativistic AND Newtonian,
348  // provided j_euler has the right limit...
349  Tenseur tmp = j_euler22_1;
350  tmp.change_triad( mp.get_bvect_spher() ) ;
351  dens = tmp(2) ;
352  dens.mult_rsint() ;
353 
354  dens = a_car() * b_car()* bbb() * dens ;
355 
356  dens.std_base_scal() ;
357 
358  p_angu_mom_2_part2_1 = new double( dens.integrale() ) ;
359 
360  }
361 
362  return *p_angu_mom_2_part2_1 ;
363 
364 }
365 
367 
368  if (p_angu_mom_1_part1_2 == 0x0) { // a new computation is required
369 
370  Cmp dens(mp) ;
371 
372  // this should work for both relativistic AND Newtonian,
373  // provided j_euler has the right limit...
374  Tenseur tmp = j_euler11_2;
375  tmp.change_triad( mp.get_bvect_spher() ) ;
376  dens = tmp(2) ;
377  dens.mult_rsint() ;
378 
379  dens = a_car() * b_car()* bbb() * dens ;
380 
381  dens.std_base_scal() ;
382 
383  p_angu_mom_1_part1_2 = new double( dens.integrale() ) ;
384 
385  }
386 
387  return *p_angu_mom_1_part1_2 ;
388 
389 }
390 
392 
393  if (p_angu_mom_2_part1_2 == 0x0) { // a new computation is required
394 
395  Cmp dens(mp) ;
396 
397  // this should work for both relativistic AND Newtonian,
398  // provided j_euler has the right limit...
399  Tenseur tmp = j_euler21_2;
400  tmp.change_triad( mp.get_bvect_spher() ) ;
401  dens = tmp(2) ;
402  dens.mult_rsint() ;
403 
404  dens = a_car() * b_car()* bbb() * dens ;
405 
406  dens.std_base_scal() ;
407 
408  p_angu_mom_2_part1_2 = new double( dens.integrale() ) ;
409 
410  }
411 
412  return *p_angu_mom_2_part1_2 ;
413 
414 }
415 
417 
418  if (p_angu_mom_1_part2_2 == 0x0) { // a new computation is required
419 
420  Cmp dens(mp) ;
421 
422  // this should work for both relativistic AND Newtonian,
423  // provided j_euler has the right limit...
424  Tenseur tmp = j_euler12_2;
425  tmp.change_triad( mp.get_bvect_spher() ) ;
426  dens = tmp(2) ;
427  dens.mult_rsint() ;
428 
429  dens = a_car() * b_car()* bbb() * dens ;
430 
431  dens.std_base_scal() ;
432 
433  p_angu_mom_1_part2_2 = new double( dens.integrale() ) ;
434 
435  }
436 
437  return *p_angu_mom_1_part2_2 ;
438 
439 }
440 
442 
443  if (p_angu_mom_2_part2_2 == 0x0) { // a new computation is required
444 
445  Cmp dens(mp) ;
446 
447  // this should work for both relativistic AND Newtonian,
448  // provided j_euler has the right limit...
449  Tenseur tmp = j_euler22_2;
450  tmp.change_triad( mp.get_bvect_spher() ) ;
451  dens = tmp(2) ;
452  dens.mult_rsint() ;
453 
454  dens = a_car() * b_car()* bbb() * dens ;
455 
456  dens.std_base_scal() ;
457 
458  p_angu_mom_2_part2_2 = new double( dens.integrale() ) ;
459 
460  }
461 
462  return *p_angu_mom_2_part2_2 ;
463 
464 }
465 
466 
467 //----------------------------//
468 // GRV2 //
469 //----------------------------//
470 
471 double Et_rot_bifluid::grv2() const {
472 
473  using namespace Unites ;
474 
475  if (p_grv2 == 0x0) { // a new computation is required
476 
477  Tenseur sou_m(mp);
478  // should work for both relativistic AND Newtonian, provided sphph_euler has right limit..
479  sou_m = 2 * qpig * a_car * sphph_euler ;
480 
482 
483  p_grv2 = new double( double(1) - lambda_grv2(sou_m(), sou_q()) ) ;
484 
485  }
486 
487  return *p_grv2 ;
488 
489 }
490 
491 
492 //----------------------------//
493 // GRV3 //
494 //----------------------------//
495 
496 double Et_rot_bifluid::grv3(ostream* ost) const {
497 
498  using namespace Unites ;
499 
500  if (p_grv3 == 0x0) { // a new computation is required
501 
502  Tenseur source(mp) ;
503 
504  // Gravitational term [cf. Eq. (43) of Gourgoulhon & Bonazzola
505  // ------------------ Class. Quantum Grav. 11, 443 (1994)]
506 
507  if (relativistic) {
508  Tenseur alpha = dzeta - logn ;
509  Tenseur beta = log( bbb ) ;
510  beta.set_std_base() ;
511 
512  source = 0.75 * ak_car - flat_scalar_prod(logn.gradient_spher(), logn.gradient_spher() )
513  + 0.5 * flat_scalar_prod(alpha.gradient_spher(), beta.gradient_spher() ) ;
514 
515  Cmp aa = alpha() - 0.5 * beta() ;
516  Cmp daadt = aa.srdsdt() ; // 1/r d/dth
517 
518  // What follows is valid only for a mapping of class Map_radial :
519  const Map_radial* mpr = dynamic_cast<const Map_radial*>(&mp) ;
520  if (mpr == 0x0) {
521  cout << "Et_rot_bifluid::grv3: the mapping does not belong"
522  << " to the class Map_radial !" << endl ;
523  abort() ;
524  }
525 
526  // Computation of 1/tan(theta) * 1/r daa/dtheta
527  if (daadt.get_etat() == ETATQCQ) {
528  Valeur& vdaadt = daadt.va ;
529  vdaadt = vdaadt.ssint() ; // division by sin(theta)
530  vdaadt = vdaadt.mult_ct() ; // multiplication by cos(theta)
531  }
532 
533  Cmp temp = aa.dsdr() + daadt ;
534  temp = ( bbb() - a_car()/bbb() ) * temp ;
535  temp.std_base_scal() ;
536 
537  // Division by r
538  Valeur& vtemp = temp.va ;
539  vtemp = vtemp.sx() ; // division by xi in the nucleus
540  // Id in the shells
541  // division by xi-1 in the ZEC
542  vtemp = (mpr->xsr) * vtemp ; // multiplication by xi/r in the nucleus
543  // by 1/r in the shells
544  // by r(xi-1) in the ZEC
545 
546  // In the ZEC, a multiplication by r has been performed instead
547  // of the division:
548  temp.set_dzpuis( temp.get_dzpuis() + 2 ) ;
549 
550  source = bbb() * source() + 0.5 * temp ;
551 
552  }
553  else{
554  source = - 0.5 * flat_scalar_prod(logn.gradient_spher(),logn.gradient_spher() ) ;
555  }
556 
557  source.set_std_base() ;
558 
559  double int_grav = source().integrale() ;
560 
561  // Matter term
562  // -----------
563 
564  // should work for relativistic AND Newtonian, provided s_euler has the right limit...
565  source = qpig * a_car * bbb * s_euler ;
566 
567  source.set_std_base() ;
568 
569  double int_mat = source().integrale() ;
570 
571  // Virial error
572  // ------------
573  if (ost != 0x0) {
574  *ost << "Et_rot_bifluid::grv3 : gravitational term : " << int_grav
575  << endl ;
576  *ost << "Et_rot_bifluid::grv3 : matter term : " << int_mat
577  << endl ;
578  }
579 
580  p_grv3 = new double( (int_grav + int_mat) / int_mat ) ;
581 
582  }
583 
584  return *p_grv3 ;
585 
586 }
587 
588 
589 //----------------------------//
590 // R_circ //
591 //----------------------------//
592 
593 double Et_rot_bifluid::r_circ2() const {
594 
595  if (p_r_circ2 == 0x0) { // a new computation is required
596 
597  // Index of the point at phi=0, theta=pi/2 at the surface of the star:
598  const Mg3d* mg = mp.get_mg() ;
599  assert(mg->get_type_t() == SYM) ;
600  int l_b = nzet - 1 ;
601  int i_b = mg->get_nr(l_b) - 1 ;
602  int j_b = mg->get_nt(l_b) - 1 ;
603  int k_b = 0 ;
604 
605  p_r_circ2 = new double( bbb()(l_b, k_b, j_b, i_b) * ray_eq2() ) ;
606 
607  }
608 
609  return *p_r_circ2 ;
610 
611 }
612 
613 
614 //----------------------------//
615 // Flattening //
616 //----------------------------//
617 
618 double Et_rot_bifluid::aplat2() const {
619 
620  if (p_aplat2 == 0x0) { // a new computation is required
621 
622  p_aplat2 = new double( ray_pole2() / ray_eq2() ) ;
623 
624  }
625 
626  return *p_aplat2 ;
627 
628 }
629 
630 
631 
632 //----------------------------//
633 // Surface area //
634 //----------------------------//
635 
636  double Et_rot_bifluid::area2() const {
637 
638  if (p_area2 == 0x0) { // a new computation is required
639  const Mg3d& mg = *(mp.get_mg()) ;
640  int np = mg.get_np(0) ;
641  int nt = mg.get_nt(0) ;
642  assert(np == 1) ; //Only valid for axisymmetric configurations
643 
644  const Map_radial* mp_rad = dynamic_cast<const Map_radial*>(&mp) ;
645  assert(mp_rad != 0x0) ;
646 
647  Valeur va_r(mg.get_angu()) ;
648  va_r.annule_hard() ;
649  Itbl lsurf2 = l_surf2() ;
650  Tbl xisurf2 = xi_surf2() ;
651 
652  for (int k=0; k<np; k++) {
653  for (int j=0; j<nt; j++) {
654  int l_star2 = lsurf2(k, j) ;
655  double xi_star2 = xisurf2(k, j) ;
656 
657  va_r.set(0, k, j, 0) = mp_rad->val_r_jk(l_star2, xi_star2, j, k) ;
658  }
659  }
660  va_r.std_base_scal() ;
661 
662  Valeur integ(mg.get_angu()) ;
663  Valeur dr = va_r.dsdt() ;
664  integ = sqrt(va_r*va_r + dr*dr) ;
665  Cmp aaaa = get_a_car()() ;
666  Valeur a2 = aaaa.va ; a2.std_base_scal() ;
667  Cmp bbbb = get_bbb()() ;
668  Valeur b = bbbb.va ; b.std_base_scal() ;
669  for (int k=0; k<np; k++) {
670  for (int j=0; j<nt; j++) {
671  integ.set(0, k, j, 0) *= sqrt(a2.val_point_jk(lsurf2(k, j), xisurf2(k, j), j, k))
672  * b.val_point_jk(lsurf2(k, j), xisurf2(k, j), j, k) * va_r(0, k, j, 0) ;
673  }
674  }
675  integ.std_base_scal() ;
676  Valeur integ2 = integ.mult_st() ;
677  double surftot = 0. ;
678  for (int j=0; j<nt; j++) {
679  surftot += (*integ2.c_cf)(0, 0, j, 0) / double(2*j+1) ;
680  }
681 
682  p_area2 = new double( 4*M_PI*surftot) ;
683 
684  }
685 
686  return *p_area2 ;
687 
688 }
689 
691 
692  return sqrt(area2()/(4*M_PI)) ;
693 
694  }
695 
696 
697 
698 
699 //----------------------------//
700 // Quadrupole moment //
701 //----------------------------//
702 
703 double Et_rot_bifluid::mom_quad() const {
704 
705  using namespace Unites ;
706 
707  if (p_mom_quad == 0x0) { // a new computation is required
708 
709  double b = mom_quad_Bo() / ( mass_g() * mass_g() ) ;
710  p_mom_quad = new double( mom_quad_old() - 4./3. * ( 1./4. + b ) * pow(mass_g(), 3) * ggrav * ggrav ) ;
711 
712  }
713 
714  return *p_mom_quad ;
715 
716 }
717 
718 
720 
721  using namespace Unites ;
722 
723  if (p_mom_quad_Bo == 0x0) { // a new computation is required
724 
725  Cmp dens(mp) ;
726 
727  dens = press() ;
728  dens = a_car() * bbb() * nnn() * dens ;
729  dens.mult_rsint() ;
730  dens.std_base_scal() ;
731 
732  p_mom_quad_Bo = new double( - 32. * dens.integrale() / qpig ) ;
733 
734  }
735 
736  return *p_mom_quad_Bo ;
737 
738 }
739 
740 
741 
743 
744  using namespace Unites ;
745 
746  if (p_mom_quad_old == 0x0) { // a new computation is required
747 
748  // Source for of the Poisson equation for nu
749  // -----------------------------------------
750 
751  Tenseur source(mp) ;
752 
753  if (relativistic) {
754  Tenseur beta = log(bbb) ;
755  beta.set_std_base() ;
756  source = qpig * a_car * enerps_euler
758  logn.gradient_spher() + beta.gradient_spher()) ;
759  }
760  else {
761  source = qpig * (nbar + nbar2);
762  }
763  source.set_std_base() ;
764 
765  // Multiplication by -r^2 P_2(cos(theta))
766  // [cf Eq.(7) of Salgado et al. Astron. Astrophys. 291, 155 (1994) ]
767  // ------------------------------------------------------------------
768 
769  // Multiplication by r^2 :
770  // ----------------------
771  Cmp& csource = source.set() ;
772  csource.mult_r() ;
773  csource.mult_r() ;
774  if (csource.check_dzpuis(2)) {
775  csource.inc2_dzpuis() ;
776  }
777 
778  // Muliplication by cos^2(theta) :
779  // -----------------------------
780  Cmp temp = csource ;
781 
782  // What follows is valid only for a mapping of class Map_radial :
783  assert( dynamic_cast<const Map_radial*>(&mp) != 0x0 ) ;
784 
785  if (temp.get_etat() == ETATQCQ) {
786  Valeur& vtemp = temp.va ;
787  vtemp = vtemp.mult_ct() ; // multiplication by cos(theta)
788  vtemp = vtemp.mult_ct() ; // multiplication by cos(theta)
789  }
790 
791  // Muliplication by -P_2(cos(theta)) :
792  // ----------------------------------
793  source = 0.5 * source() - 1.5 * temp ;
794 
795  // Final result
796  // ------------
797 
798  p_mom_quad_old = new double(- source().integrale() / qpig ) ;
799 
800  }
801 
802  return *p_mom_quad_old ;
803 
804 }
805 
806 
807 //--------------------------//
808 // Stellar surface //
809 //--------------------------//
810 
811 const Itbl& Et_rot_bifluid::l_surf() const {
812 
813  if (p_l_surf == 0x0) { // a new computation is required
814 
815  assert(p_xi_surf == 0x0) ; // consistency check
816 
817  int np = mp.get_mg()->get_np(0) ;
818  int nt = mp.get_mg()->get_nt(0) ;
819 
820  p_l_surf = new Itbl(np, nt) ;
821  p_xi_surf = new Tbl(np, nt) ;
822 
823  double nb0 = 0 ; // definition of the surface
824  double precis = 1.e-15 ;
825  int nitermax = 100 ;
826  int niter ;
827 
828  // Cmp defining the surface of the star (via the density fields)
829  //
830  Cmp surf(mp) ;
831  surf = -0.2*nbar()(0,0,0,0) ;
832  surf.annule(0, nzet-1) ;
833  surf += nbar() ; ;
834  surf = prolonge_c1(surf, nzet) ;
835 
836  (surf.va).equipot(nb0, nzet, precis, nitermax, niter, *p_l_surf,
837  *p_xi_surf) ;
838 
839  }
840 
841  return *p_l_surf ;
842 
843 }
845 
846  if (p_l_surf2 == 0x0) { // a new computation is required
847 
848  assert(p_xi_surf2 == 0x0) ; // consistency check
849 
850  int np = mp.get_mg()->get_np(0) ;
851  int nt = mp.get_mg()->get_nt(0) ;
852 
853  p_l_surf2 = new Itbl(np, nt) ;
854  p_xi_surf2 = new Tbl(np, nt) ;
855 
856  double nb0 = 0 ; // definition of the surface
857  double precis = 1.e-15 ;
858  int nitermax = 100 ;
859  int niter ;
860 
861  // Cmp defining the surface of the star (via the density fields)
862  //
863  Cmp surf2(mp) ;
864  surf2 = -0.2*nbar2()(0,0,0,0) ;
865  surf2.annule(0, nzet-1) ;
866  surf2 += nbar2() ; ;
867  surf2 = prolonge_c1(surf2, nzet) ;
868 
869  (surf2.va).equipot(nb0, nzet, precis, nitermax, niter, *p_l_surf2,
870  *p_xi_surf2) ;
871 
872  }
873 
874  return *p_l_surf2 ;
875 
876 }
877 
879 
880  if (p_xi_surf2 == 0x0) { // a new computation is required
881 
882  assert(p_l_surf2 == 0x0) ; // consistency check
883 
884  l_surf2() ; // the computation is done by l_surf2()
885 
886  }
887 
888  return *p_xi_surf2 ;
889 
890 }
891 
892 //--------------------------//
893 // Coordinate radii //
894 //--------------------------//
895 
896 double Et_rot_bifluid::ray_eq2() const {
897 
898  if (p_ray_eq2 == 0x0) { // a new computation is required
899 
900  const Mg3d& mg = *(mp.get_mg()) ;
901 
902  int type_t = mg.get_type_t() ;
903  int nt = mg.get_nt(0) ;
904 
905  if ( type_t == SYM ) {
906  assert( ( mg.get_type_p() == SYM) || (mg.get_type_p() == NONSYM) ) ;
907  int k = 0 ;
908  int j = nt-1 ;
909  int l = l_surf2()(k, j) ;
910  double xi = xi_surf2()(k, j) ;
911  double theta = M_PI / 2 ;
912  double phi = 0 ;
913 
914  p_ray_eq2 = new double( mp.val_r(l, xi, theta, phi) ) ;
915 
916  }
917  else {
918  cout << "Et_rot_bifluid::ray_eq2 : the case type_t = " << type_t
919  << " is not contemplated yet !" << endl ;
920  abort() ;
921  }
922 
923  }
924 
925  return *p_ray_eq2 ;
926 
927 }
928 
929 
931 
932  if (p_ray_eq2_pis2 == 0x0) { // a new computation is required
933 
934  const Mg3d& mg = *(mp.get_mg()) ;
935 
936  int type_t = mg.get_type_t() ;
937  int type_p = mg.get_type_p() ;
938  int nt = mg.get_nt(0) ;
939  int np = mg.get_np(0) ;
940 
941  if ( type_t == SYM ) {
942 
943  int j = nt-1 ;
944  double theta = M_PI / 2 ;
945  double phi = M_PI / 2 ;
946 
947  switch (type_p) {
948 
949  case SYM : {
950  int k = np / 2 ;
951  int l = l_surf2()(k, j) ;
952  double xi = xi_surf2()(k, j) ;
953  p_ray_eq2_pis2 = new double( mp.val_r(l, xi, theta, phi) ) ;
954  break ;
955  }
956 
957  case NONSYM : {
958  assert( np % 4 == 0 ) ;
959  int k = np / 4 ;
960  int l = l_surf2()(k, j) ;
961  double xi = xi_surf2()(k, j) ;
962  p_ray_eq2_pis2 = new double( mp.val_r(l, xi, theta, phi) ) ;
963  break ;
964  }
965 
966  default : {
967  cout << "Et_rot_bifluid::ray_eq2_pis2 : the case type_p = "
968  << type_p << " is not contemplated yet !" << endl ;
969  abort() ;
970  }
971  }
972 
973  }
974  else {
975  cout << "Et_rot_bifluid::ray_eq2_pis2 : the case type_t = " << type_t
976  << " is not contemplated yet !" << endl ;
977  abort() ;
978  }
979 
980  }
981 
982  return *p_ray_eq2_pis2 ;
983 
984 }
985 
986 
988 
989  if (p_ray_eq2_pi == 0x0) { // a new computation is required
990 
991  const Mg3d& mg = *(mp.get_mg()) ;
992 
993  int type_t = mg.get_type_t() ;
994  int type_p = mg.get_type_p() ;
995  int nt = mg.get_nt(0) ;
996  int np = mg.get_np(0) ;
997 
998  if ( type_t == SYM ) {
999 
1000  switch (type_p) {
1001 
1002  case SYM : {
1003  p_ray_eq2_pi = new double( ray_eq2() ) ;
1004  break ;
1005  }
1006 
1007  case NONSYM : {
1008  int k = np / 2 ;
1009  int j = nt-1 ;
1010  int l = l_surf2()(k, j) ;
1011  double xi = xi_surf2()(k, j) ;
1012  double theta = M_PI / 2 ;
1013  double phi = M_PI ;
1014 
1015  p_ray_eq2_pi = new double( mp.val_r(l, xi, theta, phi) ) ;
1016  break ;
1017  }
1018 
1019  default : {
1020 
1021  cout << "Et_rot_bifluid::ray_eq2_pi : the case type_t = " << type_t
1022  << " and type_p = " << type_p << endl ;
1023  cout << " is not contemplated yet !" << endl ;
1024  abort() ;
1025  break ;
1026  }
1027  }
1028  }
1029 
1030  }
1031 
1032  return *p_ray_eq2_pi ;
1033 
1034 }
1035 
1037 
1038  if (p_ray_pole2 == 0x0) { // a new computation is required
1039 
1040  assert( ((mp.get_mg())->get_type_t() == SYM)
1041  || ((mp.get_mg())->get_type_t() == NONSYM) ) ;
1042 
1043  int k = 0 ;
1044  int j = 0 ;
1045  int l = l_surf2()(k, j) ;
1046  double xi = xi_surf2()(k, j) ;
1047  double theta = 0 ;
1048  double phi = 0 ;
1049 
1050  p_ray_pole2 = new double( mp.val_r(l, xi, theta, phi) ) ;
1051 
1052  }
1053 
1054  return *p_ray_pole2 ;
1055 
1056 }
1057 
1058 
1059 
1060 }
int get_type_p() const
Returns the type of sampling in the direction: SYM : : symmetry with respect to the transformatio...
Definition: grilles.h:495
const Cmp & dsdr() const
Returns of *this .
Definition: cmp_deriv.C:84
double get_m2() const
Return the individual particule mass .
Definition: eos_bifluid.h:262
double * p_angu_mom_2_part1_2
To compute Ip (2nd version)
double * p_mass_b2
Baryon mass of fluid 2.
double * p_angu_mom_1_part1_1
To compute In (1st version)
double * p_mom_quad_old
Part of the quadrupole moment.
Definition: etoile.h:1642
const Valeur & dsdt() const
Returns of *this.
Definition: valeur_dsdt.C:112
Cmp log(const Cmp &)
Neperian logarithm.
Definition: cmp_math.C:296
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition: valeur.h:302
double * p_angu_mom_1_part1_2
To compute In (2nd version)
double * p_mom_quad_Bo
Part of the quadrupole moment.
Definition: etoile.h:1643
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:446
virtual double mass_g() const
Gravitational mass.
Tenseur tnphi
Component of the shift vector.
Definition: etoile.h:1515
const Tenseur & get_a_car() const
Returns the total conformal factor .
Definition: etoile.h:733
virtual double angu_mom_1_part2_2() const
To compute Xn (2nd version)
Tenseur j_euler2
To compute Jp.
const Tenseur & gradient_spher() const
Returns the gradient of *this (Spherical coordinates) (scalar field only).
Definition: tenseur.C:1548
Tenseur j_euler1
To compute Jn.
virtual double mass_b() const
Total Baryon mass.
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:462
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:220
double mass_b1() const
Baryon mass of fluid 1.
void annule(int l)
Sets the Cmp to zero in a given domain.
Definition: cmp.C:348
void set_std_base()
Set the standard spectal basis of decomposition for each component.
Definition: tenseur.C:1170
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition: map.h:783
double * p_r_circ2
Circumferential radius of fluid no.2.
virtual double angu_mom_1_part2_1() const
To compute Xn (1st version)
Tbl * p_xi_surf2
Description of the surface of fluid 2: 2-D Tbl containing the values of the radial coordinate on the...
Lorene prototypes.
Definition: app_hor.h:64
Standard units of space, time and mass.
void annule_hard()
Sets the Valeur to zero in a hard way.
Definition: valeur.C:723
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:765
int get_etat() const
Returns the logical state.
Definition: cmp.h:899
Tenseur nnn
Total lapse function.
Definition: etoile.h:509
Tenseur s_euler
Trace of the stress tensor in the Eulerian frame.
Definition: etoile.h:468
double * p_angu_mom_1
Angular momentum of fluid 1.
Tenseur flat_scalar_prod(const Tenseur &t1, const Tenseur &t2)
Scalar product of two Tenseur when the metric is : performs the contraction of the last index of t1 w...
double * p_mom_quad
Quadrupole moment.
Definition: etoile.h:1641
const Cmp & srdsdt() const
Returns of *this .
Definition: cmp_deriv.C:105
Tenseur b_car
Square of the metric factor B.
Definition: etoile.h:1507
int get_type_t() const
Returns the type of sampling in the direction: SYM : : symmetry with respect to the equatorial pl...
Definition: grilles.h:485
Basic integer array class.
Definition: itbl.h:122
const Valeur & sx() const
Returns (r -sampling = RARE ) \ Id (r sampling = FIN ) \ (r -sampling = UNSURR ) ...
Definition: valeur_sx.C:110
Values and coefficients of a (real-value) function.
Definition: valeur.h:287
double * p_ray_eq2
Coordinate radius at , .
double * p_angu_mom_1_part2_1
To compute Xn (1st version)
Tenseur press
Fluid pressure.
Definition: etoile.h:461
virtual double mom_quad_Bo() const
Part of the quadrupole moment.
const Itbl & l_surf2() const
Description of the surface of fluid 2: returns a 2-D Itbl containing the values of the domain index l...
void std_base_scal()
Sets the bases for spectral expansions (member base ) to the standard ones for a scalar.
Definition: valeur.C:824
Tenseur j_euler
Total angular momentum (flat-space!) 3-vector , which is related to of the "3+1" decomposition...
void mult_r()
Multiplication by r everywhere.
Definition: cmp_r_manip.C:91
double * p_mass_b1
Baryon mass of fluid 1.
virtual double val_r(int l, double xi, double theta, double pphi) const =0
Returns the value of the radial coordinate r for a given in a given domain.
Tenseur gam_euler2
Lorentz factor between the fluid 2 and Eulerian observers.
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition: tenseur.C:824
double * p_grv3
Error on the virial identity GRV3.
Definition: etoile.h:1634
void change_triad(const Base_vect &new_triad)
Sets a new vectorial basis (triad) of decomposition and modifies the components accordingly.
Definition: tenseur.C:668
Tenseur j_euler11_1
To compute In (1st version)
double ray_eq2_pi() const
Coordinate radius for fluid 2 at , [r_unit].
const Eos_bifluid & eos
Equation of state for two-fluids model.
double * p_angu_mom_2_part1_1
To compute Ip (1st version)
const Valeur & ssint() const
Returns of *this.
Definition: valeur_ssint.C:112
Tenseur nbar
Baryon density in the fluid frame.
Definition: etoile.h:459
const Mg3d * get_angu() const
Returns the pointer on the associated angular grid.
Definition: mg3d.C:473
Tenseur gam_euler
Lorentz factor between the fluid and Eulerian observers.
Definition: etoile.h:471
virtual double grv3(ostream *ost=0x0) const
Error on the virial identity GRV3.
Base class for pure radial mappings.
Definition: map.h:1536
virtual double angu_mom_2() const
Angular momentum of fluid 2.
void mult_rsint()
Multiplication by .
Definition: cmp_r_manip.C:116
virtual double angu_mom_1_part1_1() const
To compute In (1st version)
Map & mp
Mapping associated with the star.
Definition: etoile.h:429
double * p_area2
Surface area of fluid no.2.
const Valeur & mult_st() const
Returns applied to *this.
virtual double val_r_jk(int l, double xi, int j, int k) const =0
Returns the value of the radial coordinate r for a given and a given collocation point in in a give...
double integrale() const
Computes the integral over all space of *this .
Definition: cmp_integ.C:55
Tbl * p_xi_surf
Description of the stellar surface: 2-D Tbl containing the values of the radial coordinate on the su...
Definition: etoile.h:545
const Tbl & xi_surf2() const
Description of the surface of fluid 2: returns a 2-D Tbl containing the values of the radial coordina...
double ray_eq2_pis2() const
Coordinate radius for fluid 2 at , [r_unit].
Tenseur ak_car
Scalar .
Definition: etoile.h:1586
int get_etat() const
Returns the logical state.
Definition: tenseur.h:707
static double lambda_grv2(const Cmp &sou_m, const Cmp &sou_q)
Computes the coefficient which ensures that the GRV2 virial identity is satisfied.
virtual double angu_mom_1_part1_2() const
To compute In (2nd version)
Tenseur enerps_euler
the combination : useful because in the Newtonian limit .
double val_point_jk(int l, double x, int j, int k) const
Computes the value of the field represented by *this at an arbitrary point in , but collocation point...
Definition: valeur.C:900
Coord xsr
in the nucleus; \ 1/R in the non-compactified shells; \ in the compactified outer domain...
Definition: map.h:1549
Tenseur bbb
Metric factor B.
Definition: etoile.h:1504
double * p_ray_eq2_pi
Coordinate radius at , .
virtual double angu_mom_2_part2_2() const
To compute Xp (2nd version)
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
virtual const Itbl & l_surf() const
Description of the surface of fluid 1: returns a 2-D Itbl containing the values of the domain index l...
virtual double mom_quad_old() const
Part of the quadrupole moment.
virtual double grv2() const
Error on the virial identity GRV2.
void inc2_dzpuis()
Increases by 2 the value of dzpuis and changes accordingly the values of the Cmp in the external comp...
Definition: cmp_r_manip.C:192
double * p_angu_mom_2_part2_2
To compute Xp (2nd version)
double * p_mass_g
Gravitational mass.
Definition: etoile.h:548
virtual double mean_radius2() const
Mean radius for fluid 2.
void std_base_scal()
Sets the spectral bases of the Valeur va to the standard ones for a scalar.
Definition: cmp.C:644
double * p_angu_mom_2_part2_1
To compute Xp (1st version)
int nzet
Number of domains of *mp occupied by the star.
Definition: etoile.h:432
double ray_pole2() const
Coordinate radius for fluid 2 at [r_unit].
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition: grilles.h:452
virtual double angu_mom_1() const
Angular momentum of fluid 1.
Tenseur a_car
Total conformal factor .
Definition: etoile.h:515
double * p_angu_mom_2
Angular momentum of fluid 2.
virtual double angu_mom_2_part2_1() const
To compute Xp (1st version)
bool relativistic
Indicator of relativity: true for a relativistic star, false for a Newtonian one. ...
Definition: etoile.h:437
Multi-domain grid.
Definition: grilles.h:273
Tenseur sphph_euler
The component of the stress tensor .
virtual double mom_quad() const
Quadrupole moment.
int get_dzpuis() const
Returns dzpuis.
Definition: cmp.h:903
double * p_ray_eq2_pis2
Coordinate radius at , .
double * p_ray_pole2
Coordinate radius at .
virtual double r_circ2() const
Circumferential radius for fluid 2.
Tenseur nbar2
Baryon density in the fluid frame, for fluid 2.
virtual double angu_mom_2_part1_1() const
To compute Ip (1st version)
bool check_dzpuis(int dzi) const
Returns false if the last domain is compactified and *this is not zero in this domain and dzpuis is n...
Definition: cmp.C:715
const Valeur & mult_ct() const
Returns applied to *this.
Tenseur & logn
Metric potential = logn_auto.
Definition: etoile.h:1521
void set_dzpuis(int)
Set a value to dzpuis.
Definition: cmp.C:654
Basic array class.
Definition: tbl.h:161
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:457
virtual double angu_mom() const
Angular momentum.
double mass_b2() const
Baryon mass of fluid 2.
double * p_grv2
Error on the virial identity GRV2.
Definition: etoile.h:1633
const Tenseur & get_bbb() const
Returns the metric factor B.
Definition: etoile.h:1712
virtual double angu_mom_2_part1_2() const
To compute Ip (2nd version)
double ray_eq2() const
Coordinate radius for fluid 2 at , [r_unit].
Tenseur & dzeta
Metric potential = beta_auto.
Definition: etoile.h:1534
double * p_angu_mom_1_part2_2
To compute Xn (2nd version)
virtual double aplat2() const
Flatening r_pole/r_eq for fluid 2.
double * p_aplat2
Flatening r_pole/r_eq of fluid no.2.
Valeur va
The numerical value of the Cmp.
Definition: cmp.h:464
double * p_mass_b
Baryon mass.
Definition: etoile.h:547
double get_m1() const
Return the individual particule mass .
Definition: eos_bifluid.h:256
virtual double area2() const
Surface area for fluid 2.
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition: tenseur.h:301
Tenseur j_euler12_1
To compute Ip (1st version)
double * p_angu_mom
Angular momentum.
Definition: etoile.h:1631
Itbl * p_l_surf2
Description of the surface of fluid 2: 2-D Itbl containing the values of the domain index l on the su...
Itbl * p_l_surf
Description of the stellar surface: 2-D Itbl containing the values of the domain index l on the surfa...
Definition: etoile.h:539