28char star_QI_global_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Compobj/star_QI_global.C,v 1.3 2014/10/13 08:52:49 j_novak Exp $" ;
133 - dlogn(1)*dlogn(1) - dlogn(2)*dlogn(2) - dlogn(3)*dlogn(3) ;
169 - dlogn(1)*dlogn(1) - dlogn(2)*dlogn(2) - dlogn(3)*dlogn(3)
170 + 0.5 * ( dalpha(1)*dbet(1) + dalpha(2)*dbet(2) + dalpha(3)*dbet(3) ) ;
172 Scalar aa = alpha - 0.5 * bet ;
178 cout <<
"Star_rot::grv3: the mapping does not belong"
179 <<
" to the class Map_radial !" << endl ;
186 vdaadt = vdaadt.
ssint() ;
199 vtemp = (mpr->
xsr) * vtemp ;
207 source =
bbb * source + 0.5 * temp ;
223 source = qpig *
a_car *
bbb * s_euler ;
232 *ost <<
"Star_rot::grv3 : gravitational term : " << int_grav
234 *ost <<
"Star_rot::grv3 : matter term : " << int_mat
238 p_grv3 =
new double( (int_grav + int_mat) / int_mat ) ;
274 - dlogn(1)*dlogn_bet(1) - dlogn(2)*dlogn_bet(2) - dlogn(3)*dlogn_bet(3) ;
305 source = 0.5 * source - 1.5 * temp ;
326 cout <<
"Star_rot::lambda_grv2: the mapping of sou_m does not"
327 << endl <<
" belong to the class Map_radial !" << endl ;
331 assert( &sou_q.
get_mp() == mprad ) ;
341 double theta0 = M_PI / 2 ;
346 for (
int l=0 ; l<nz ; l++) {
347 double rmax = mprad->
val_r(l,
double(1), theta0, phi0) ;
350 double rmin = mprad->
val_r(l,
double(0), theta0, phi0) ;
357 double rmin = mprad->
val_r(l,
double(-1), theta0, phi0) ;
358 mpaff.
set_alpha(
double(.5) * (rmax - rmin), l ) ;
359 mpaff.
set_beta(
double(.5) * (rmax + rmin), l) ;
364 double rmin = mprad->
val_r(l,
double(-1), theta0, phi0) ;
365 double umax = double(1) / rmin ;
366 double umin = double(1) / rmax ;
367 mpaff.
set_alpha(
double(.5) * (umin - umax), l) ;
368 mpaff.
set_beta(
double(.5) * (umin + umax), l) ;
373 cout <<
"Star_rot::lambda_grv2: unknown type_r ! " << endl ;
390 for (
int l=0; l<nz; l++) {
394 *(jac.
t[l]) = *(jac.
t[l]) / (a1*a1) ;
401 assert( jac.
t[l]->
get_etat() == ETATQCQ ) ;
402 double* tjac = jac.
t[l]->
t ;
404 for (
int k=0; k<mg->
get_np(l); k++) {
405 for (
int j=0; j<mg->
get_nt(l); j++) {
406 for (
int i=0; i<mg->
get_nr(l); i++) {
408 (a1 * (a1 * xi[i] + b1) ) ;
419 *(jac.
t[l]) = - *(jac.
t[l]) / (a1*a1) ;
424 cout <<
"Star_rot::lambda_grv2: unknown type_r ! " << endl ;
438 if ( sou_m.
get_etat() == ETATZERO ) {
448 if ( sou_q.
get_etat() == ETATZERO ) {
471 int* ndl =
new int[nz+4] ;
473 for (
int l=0; l<nz; l++) {
474 ndl[1+l] = mg->
get_nr(l) ;
483 for (
int l=0; l<nz ; l++) {
484 nrmax = ( ndl[1+l] > nrmax ) ? ndl[1+l] : nrmax ;
486 int ndr = nrmax + 5 ;
493 double* erre =
new double [nz*ndr] ;
495 for (
int l=0; l<nz; l++) {
498 for (
int i=0; i<ndl[1+l]; i++) {
500 erre[ ndr*l + i ] = a1 * xi + b1 ;
508 int ndrtp = ndr*ndt*ndp ;
509 int taille = ndrtp*nz ;
511 double* tsou_m =
new double[ taille ] ;
512 double* tsou_q =
new double[ taille ] ;
515 for (
int i=0; i<taille; i++) {
523 for (
int l=0; l<nz; l++) {
524 for (
int k=0; k<np1; k++) {
525 for (
int j=0; j<nt; j++) {
526 for (
int i=0; i<mg->
get_nr(l); i++) {
527 double xx = s_m(l, k, j, i) ;
528 tsou_m[ndrtp*l + ndrt*k + ndr*j + i] = xx ;
530 tsou_m[ndrtp*l + ndrt*k + ndr*(nt2-1-j) + i] = xx ;
539 for (
int l=0; l<nz; l++) {
540 for (
int k=0; k<np1; k++) {
541 for (
int j=0; j<nt; j++) {
542 for (
int i=0; i<mg->
get_nr(l); i++) {
543 double xx = s_q(l, k, j, i) ;
544 tsou_q[ndrtp*l + ndrt*k + ndr*j + i] = xx ;
546 tsou_q[ndrtp*l + ndrt*k + ndr*(nt2-1-j) + i] = xx ;
556 double int_m, int_q ;
557 F77_integrale2d(ndl, &ndr, &ndt, &ndp, erre, tsou_m, &int_m) ;
558 F77_integrale2d(ndl, &ndr, &ndt, &ndp, erre, tsou_q, &int_q) ;
572 if ( int_q !=
double(0) ) {
573 lambda = - int_m / int_q ;
Scalar b_car
Square of the metric factor B.
Scalar bbb
Metric factor B.
Scalar a_car
Square of the metric factor A.
double * p_angu_mom
Angular momentum.
Vector mom_euler
Total 3-momentum density in the Eulerian frame.
Sym_tensor stress_euler
Stress tensor with respect to the Eulerian observer.
Scalar ener_euler
Total energy density E in the Eulerian frame.
Scalar nn
Lapse function N .
Map & mp
Mapping describing the coordinate system (r,theta,phi)
double * x
Array of values of at the nr collocation points.
const double * get_beta() const
Returns the pointer on the array beta.
const double * get_alpha() const
Returns the pointer on the array alpha.
void set_beta(double beta0, int l)
Modifies the value of in domain no. l.
void set_alpha(double alpha0, int l)
Modifies the value of in domain no. l.
Base class for pure radial mappings.
Coord xsr
in the nucleus; \ 1/R in the non-compactified shells; \ in the compactified outer domain.
Coord dxdr
in the nucleus and in the non-compactified shells; \ in the compactified outer domain.
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
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.
const Metric_flat & flat_met_spher() const
Returns the flat metric associated with the spherical coordinates and with components expressed in th...
const Grille3d * get_grille3d(int l) const
Returns a pointer on the 3D mono-grid for domain no. l.
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
int get_nzone() const
Returns the number of domains.
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
Tbl ** t
Array (size nzone ) of pointers on the Tbl 's.
Tensor field of valence 0 (or component of a tensorial field).
const Scalar & srdsdt() const
Returns of *this .
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this
int get_dzpuis() const
Returns dzpuis.
double integrale() const
Computes the integral over all space of *this .
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field.
virtual void inc_dzpuis(int inc=1)
Increases by inc units the value of dzpuis and changes accordingly the values of the Scalar in the co...
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...
const Scalar & dsdr() const
Returns of *this .
Valeur & set_spectral_va()
Returns va (read/write version)
const Valeur & get_spectral_va() const
Returns va (read only version)
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
void set_dzpuis(int)
Modifies the dzpuis flag.
void mult_r()
Multiplication by r everywhere; dzpuis is not changed.
virtual double grv2() const
Error on the virial identity GRV2.
double * p_grv2
Error on the virial identity GRV2.
Scalar logn
Logarithm of the lapse N .
double * p_mass_g
Gravitational mass (ADM mass as a volume integral)
virtual double grv3(ostream *ost=0x0) const
Error on the virial identity GRV3.
virtual double mom_quad() const
Quadrupole moment.
double * p_grv3
Error on the virial identity GRV3.
static double lambda_grv2(const Scalar &sou_m, const Scalar &sou_q)
Computes the coefficient which ensures that the GRV2 virial identity is satisfied.
virtual double angu_mom() const
Angular momentum.
virtual double mass_g() const
Gravitational mass.
Scalar tnphi
Component of the shift vector.
double * p_mom_quad
Quadrupole moment
Scalar dzeta
Metric potential .
int get_etat() const
Gives the logical state.
double * t
The array of double.
Values and coefficients of a (real-value) function.
const Valeur & mult_ct() const
Returns applied to *this.
const Valeur & sx() const
Returns (r -sampling = RARE ) \ Id (r sampling = FIN ) \ (r -sampling = UNSURR )
int get_etat() const
Returns the logical state.
Mtbl * c
Values of the function at the points of the multi-grid
const Valeur & ssint() const
Returns of *this.
void coef_i() const
Computes the physical value of *this.
Tensor field of valence 1.
Cmp log(const Cmp &)
Neperian logarithm.
const Map & get_mp() const
Returns the mapping.
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Tensor trace(int ind1, int ind2) const
Trace on two different type indices.
Standard units of space, time and mass.