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
23char 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
115namespace Lorene {
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
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
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
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
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;
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;
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;
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;
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;
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;
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;
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;
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;
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;
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;
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
471double 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
496double 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
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
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
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
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
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
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}
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
void mult_rsint()
Multiplication by .
int get_etat() const
Returns the logical state.
Definition cmp.h:899
Valeur va
The numerical value of the Cmp
Definition cmp.h:464
void std_base_scal()
Sets the spectral bases of the Valeur va to the standard ones for a scalar.
Definition cmp.C:644
void annule(int l)
Sets the Cmp to zero in a given domain.
Definition cmp.C:348
int get_dzpuis() const
Returns dzpuis.
Definition cmp.h:903
void mult_r()
Multiplication by r everywhere.
Definition cmp_r_manip.C:91
void inc2_dzpuis()
Increases by 2 the value of dzpuis and changes accordingly the values of the Cmp in the external comp...
void set_dzpuis(int)
Set a value to dzpuis.
Definition cmp.C:654
double integrale() const
Computes the integral over all space of *this .
Definition cmp_integ.C:55
const Cmp & srdsdt() const
Returns of *this .
Definition cmp_deriv.C:105
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 Cmp & dsdr() const
Returns of *this .
Definition cmp_deriv.C:84
double get_m1() const
Return the individual particule mass
double get_m2() const
Return the individual particule mass
double * p_angu_mom_1_part2_2
To compute Xn (2nd version)
Tenseur gam_euler2
Lorentz factor between the fluid 2 and Eulerian observers
Tenseur enerps_euler
the combination : useful because in the Newtonian limit .
double * p_mass_b1
Baryon mass of fluid 1.
virtual double angu_mom_1_part1_1() const
To compute In (1st version)
double * p_mass_b2
Baryon mass of fluid 2.
virtual double grv3(ostream *ost=0x0) const
Error on the virial identity GRV3.
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...
Tenseur j_euler12_1
To compute Ip (1st version)
virtual double mom_quad_Bo() const
Part of the quadrupole moment.
double * p_angu_mom_1_part1_2
To compute In (2nd version)
Tenseur j_euler2
To compute Jp.
Tenseur sphph_euler
The component of the stress tensor .
double * p_angu_mom_2_part2_1
To compute Xp (1st version)
virtual double area2() const
Surface area for fluid 2.
virtual double angu_mom_2_part1_2() const
To compute Ip (2nd version)
Tbl * p_xi_surf2
Description of the surface of fluid 2: 2-D Tbl containing the values of the radial coordinate on the...
virtual double angu_mom_2_part1_1() const
To compute Ip (1st version)
double * p_angu_mom_1
Angular momentum of fluid 1.
virtual double angu_mom_2() const
Angular momentum of fluid 2.
Tenseur j_euler11_1
To compute In (1st version)
virtual double angu_mom_2_part2_2() const
To compute Xp (2nd version)
double * p_ray_eq2
Coordinate radius at , .
double * p_angu_mom_1_part2_1
To compute Xn (1st version)
double ray_pole2() const
Coordinate radius for fluid 2 at [r_unit].
virtual double mass_b() const
Total Baryon mass.
Tenseur j_euler1
To compute Jn.
Tenseur j_euler
Total angular momentum (flat-space!) 3-vector , which is related to of the "3+1" decomposition,...
const Tbl & xi_surf2() const
Description of the surface of fluid 2: returns a 2-D Tbl containing the values of the radial coordi...
double * p_ray_pole2
Coordinate radius at .
virtual double mass_g() const
Gravitational mass.
virtual double angu_mom_2_part2_1() const
To compute Xp (1st version)
double * p_aplat2
Flatening r_pole/r_eq of fluid no.2.
double * p_angu_mom_2_part2_2
To compute Xp (2nd version)
double * p_area2
Surface area of fluid no.2.
double * p_angu_mom_2
Angular momentum of fluid 2.
const Itbl & l_surf2() const
Description of the surface of fluid 2: returns a 2-D Itbl containing the values of the domain index...
const Eos_bifluid & eos
Equation of state for two-fluids model.
double * p_angu_mom_2_part1_2
To compute Ip (2nd version)
double mass_b1() const
Baryon mass of fluid 1.
virtual double angu_mom_1_part1_2() const
To compute In (2nd version)
virtual double angu_mom_1_part2_1() const
To compute Xn (1st version)
double ray_eq2() const
Coordinate radius for fluid 2 at , [r_unit].
virtual double r_circ2() const
Circumferential radius for fluid 2.
double mass_b2() const
Baryon mass of fluid 2.
virtual double mom_quad() const
Quadrupole moment.
virtual double aplat2() const
Flatening r_pole/r_eq for fluid 2.
virtual double grv2() const
Error on the virial identity GRV2.
double * p_angu_mom_1_part1_1
To compute In (1st version)
double * p_r_circ2
Circumferential radius of fluid no.2.
double * p_ray_eq2_pi
Coordinate radius at , .
virtual double angu_mom_1_part2_2() const
To compute Xn (2nd version)
double * p_angu_mom_2_part1_1
To compute Ip (1st version)
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...
Tenseur nbar2
Baryon density in the fluid frame, for fluid 2.
virtual double angu_mom_1() const
Angular momentum of fluid 1.
virtual double angu_mom() const
Angular momentum.
double * p_ray_eq2_pis2
Coordinate radius at , .
double ray_eq2_pis2() const
Coordinate radius for fluid 2 at , [r_unit].
double ray_eq2_pi() const
Coordinate radius for fluid 2 at , [r_unit].
virtual double mean_radius2() const
Mean radius for fluid 2.
virtual double mom_quad_old() const
Part of the quadrupole moment.
Tenseur & logn
Metric potential = logn_auto.
Definition etoile.h:1521
double * p_mom_quad_old
Part of the quadrupole moment.
Definition etoile.h:1642
const Tenseur & get_bbb() const
Returns the metric factor B.
Definition etoile.h:1712
static double lambda_grv2(const Cmp &sou_m, const Cmp &sou_q)
Computes the coefficient which ensures that the GRV2 virial identity is satisfied.
Tenseur bbb
Metric factor B.
Definition etoile.h:1504
Tenseur ak_car
Scalar .
Definition etoile.h:1586
Tenseur & dzeta
Metric potential = beta_auto.
Definition etoile.h:1534
double * p_grv3
Error on the virial identity GRV3.
Definition etoile.h:1634
double * p_grv2
Error on the virial identity GRV2.
Definition etoile.h:1633
double * p_mom_quad_Bo
Part of the quadrupole moment.
Definition etoile.h:1643
double * p_mom_quad
Quadrupole moment.
Definition etoile.h:1641
double * p_angu_mom
Angular momentum.
Definition etoile.h:1631
Tenseur b_car
Square of the metric factor B.
Definition etoile.h:1507
Tenseur tnphi
Component of the shift vector.
Definition etoile.h:1515
int nzet
Number of domains of *mp occupied by the star.
Definition etoile.h:432
double * p_mass_b
Baryon mass.
Definition etoile.h:547
double * p_mass_g
Gravitational mass.
Definition etoile.h:548
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
Tenseur nnn
Total lapse function.
Definition etoile.h:509
Tenseur nbar
Baryon density in the fluid frame.
Definition etoile.h:459
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
Tenseur gam_euler
Lorentz factor between the fluid and Eulerian observers.
Definition etoile.h:471
Map & mp
Mapping associated with the star.
Definition etoile.h:429
const Tenseur & get_a_car() const
Returns the total conformal factor .
Definition etoile.h:733
Tenseur press
Fluid pressure.
Definition etoile.h:461
bool relativistic
Indicator of relativity: true for a relativistic star, false for a Newtonian one.
Definition etoile.h:437
Tenseur s_euler
Trace of the stress tensor in the Eulerian frame.
Definition etoile.h:468
Tenseur a_car
Total conformal factor .
Definition etoile.h:515
Basic integer array class.
Definition itbl.h:122
Base class for pure radial mappings.
Definition map.h:1536
Coord xsr
in the nucleus; \ 1/R in the non-compactified shells; \ in the compactified outer domain.
Definition map.h:1549
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...
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition map.h:783
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.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
Multi-domain grid.
Definition grilles.h:273
const Mg3d * get_angu() const
Returns the pointer on the associated angular grid.
Definition mg3d.C:473
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
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition grilles.h:462
int get_type_p() const
Returns the type of sampling in the direction: SYM : : symmetry with respect to the transformatio...
Definition grilles.h:495
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition grilles.h:457
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition grilles.h:452
Basic array class.
Definition tbl.h:161
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition tenseur.h:301
Cmp & set()
Read/write for a scalar (see also operator=(const Cmp&) ).
Definition tenseur.C:824
void set_std_base()
Set the standard spectal basis of decomposition for each component.
Definition tenseur.C:1170
const Tenseur & gradient_spher() const
Returns the gradient of *this (Spherical coordinates) (scalar field only).
Definition tenseur.C:1548
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
int get_etat() const
Returns the logical state.
Definition tenseur.h:707
Values and coefficients of a (real-value) function.
Definition valeur.h:287
const Valeur & mult_ct() const
Returns applied to *this.
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
const Valeur & sx() const
Returns (r -sampling = RARE ) \ Id (r sampling = FIN ) \ (r -sampling = UNSURR )
Definition valeur_sx.C:110
const Valeur & dsdt() const
Returns of *this.
Tbl & set(int l)
Read/write of the value in a given domain (configuration space).
Definition valeur.h:363
const Valeur & ssint() const
Returns of *this.
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition valeur.h:302
const Valeur & mult_st() const
Returns applied to *this.
void std_base_scal()
Sets the bases for spectral expansions (member base ) to the standard ones for a scalar.
Definition valeur.C:824
void annule_hard()
Sets the Valeur to zero in a hard way.
Definition valeur.C:723
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:220
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
Cmp log(const Cmp &)
Neperian logarithm.
Definition cmp_math.C:296
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...
Lorene prototypes.
Definition app_hor.h:64
Standard units of space, time and mass.