26 char spheroid_C[] =
"$Header: /cvsroot/Lorene/C++/Source/App_hor/spheroid.C,v 1.21 2014/10/13 08:52:38 j_novak Exp $" ;
96 jac2d(map, 2, COV, map.get_bvect_spher()),
97 proj(map, 2, COV, map.get_bvect_spher()),
98 qq(map, COV, map.get_bvect_spher()),
99 ss (map, COV, map.get_bvect_spher()),
100 ephi(map, CON, map.get_bvect_spher()),
101 qab(map.flat_met_spher()),
103 hh(map, COV, map.get_bvect_spher()),
105 ll(map, COV, map.get_bvect_spher()),
106 jj(map, COV, map.get_bvect_spher()),
119 assert(radius > 1.e-15) ;
137 for (
int i=1; i<=3; i++)
138 hh.
set(i,i) = 2./radius ;
152 jac2d(h_in.get_mp(),2, COV, h_in.get_mp().get_bvect_spher()),
153 proj(h_in.get_mp(),2, COV, h_in.get_mp().get_bvect_spher()),
154 qq(h_in.get_mp(), COV, h_in.get_mp().get_bvect_spher()),
155 ss (h_in.get_mp(), COV, h_in.get_mp().get_bvect_spher()),
156 ephi(h_in.get_mp(), CON, h_in.get_mp().get_bvect_spher()),
157 qab(h_in.get_mp().flat_met_spher()),
158 ricci(h_in.get_mp()),
159 hh(h_in.get_mp(), COV, h_in.get_mp().get_bvect_spher()),
161 ll(h_in.get_mp(), COV, h_in.get_mp().get_bvect_spher()),
162 jj(h_in.get_mp(), COV, h_in.get_mp().get_bvect_spher()),
175 assert(mp_rad != 0x0) ;
180 int np = gri2d.
get_np(0) ;
181 int nt = gri2d.
get_nt(0) ;
183 int nr3 = gri3d.
get_nr(0) ;
184 int nt3 = gri3d.
get_nt(0) ;
185 int np3 = gri3d.
get_np(0) ;
187 assert( nt == nt3 ) ;
188 assert ( np == np3 );
208 for (
int f= 0; f<nz; f++)
209 for (
int k=0; k<np3; k++)
210 for (
int j=0; j<nt3; j++) {
211 for (
int l=0; l<nr3; l++) {
230 jac.std_spectral_base();
234 jac.set(1,2) = -h_surf3.
srdsdt() ;
241 for (
int l=1; l<4; l++)
242 for (
int m=1; m<4; m++)
243 for (
int k=0; k<np; k++)
244 for (
int j=0; j<nt; j++) {
246 j, k, pipo, lz, xi) ;
248 jac(l,m).get_spectral_va().val_point_jk(lz, xi, j, k) ;
253 jac_inv.
set(1,2) = - jac_inv(1,2);
254 jac_inv.
set(1,3) = - jac_inv(1,3) ;
264 ss3d = ss3d /
sqrt (ssnorm) ;
265 ss3dcon = ss3dcon /
sqrt (ssnorm) ;
280 for (
int ii=1; ii <=3; ii++){
284 for (
int ii=1; ii <=3; ii++)
285 for (
int jjy = 1; jjy <=3; jjy ++){
286 jac_inv.
set(ii, jjy).
dec_dzpuis(jac_inv(ii, jjy).get_dzpuis());
287 jac.set(ii, jjy).dec_dzpuis(jac(ii, jjy).get_dzpuis());
295 ss3 =
contract(jac_inv, 0, ss3d, 0);
299 for (
int l=1; l<4; l++)
300 for (
int k=0; k<np; k++)
301 for (
int j=0; j<nt; j++) {
305 ss3(l).get_spectral_va().val_point_jk(lz, xi, j, k) ;
312 ephi3.set(2).annule_hard();
314 ephi33.mult_r(); ephi33.mult_sint();
315 ephi3.set(3) = ephi33;
316 ephi3.std_spectral_base();
321 for (
int l=1; l<4; l++)
322 for (
int k=0; k<np; k++)
323 for (
int j=0; j<nt; j++) {
325 j, k, pipo, lz, xi) ;
327 ephi3(l).get_spectral_va().val_point_jk(lz, xi, j, k) ;
332 Tensor proj_prov = gamij.
con().
down(1, gamij) - ss3dcon*ss3d;
335 proj3d.allocate_all();
337 proj3d.std_spectral_base();
339 for (
int l=1; l<4; l++)
340 for (
int m=1; m<4; m++)
341 for (
int k=0; k<np; k++)
342 for (
int j=0; j<nt; j++) {
344 j, k, pipo, lz, xi) ;
346 proj3d(l,m).get_spectral_va().val_point_jk(lz, xi, j, k) ;
359 qq3d.set(2,2)= gamij.
cov()(1,1) * (h_surf3.
srdsdt())* (h_surf3.
srdsdt())
360 + 2* gamij.
cov()(1,2)* (h_surf3.
srdsdt()) + gamij.
cov()(2,2);
370 qab2.allocate_all() ;
371 qab2.std_spectral_base();
372 for (
int l=1; l<4; l++)
373 for (
int m=1; m<4; m++)
374 for (
int k=0; k<np; k++)
375 for (
int j=0; j<nt; j++) {
377 j, k, pipo, lz, xi) ;
378 qab2.set(l,m).set_grid_point(0, k, j, 0) =
379 qq3d(l,m).get_spectral_va().val_point_jk(lz, xi, j, k) ;
385 - ss3d * ss3d) , 0), 1) ;
388 for (
int l=1; l<4; l++)
389 for (
int m=1; m<4; m++)
390 for (
int k=0; k<np; k++)
391 for (
int j=0; j<nt; j++) {
393 j, k, pipo, lz, xi) ;
395 qqq(l,m).get_spectral_va().val_point_jk(lz, xi, j, k) ;
401 for (
int k=0; k<np; k++)
402 for (
int j=0; j<nt; j++) {
404 j, k, pipo, lz, xi) ;
414 for (
int k=0; k<np; k++)
415 for (
int j=0; j<nt; j++) {
427 for (
int k=0; k<np; k++)
428 for (
int j=0; j<nt; j++) {
437 double rayon =
sqrt(
area()/(4.*M_PI));
438 ftilde = -ftilde/(rayon*rayon);
446 double *a_tilde =
new double[nombre];
450 for (
int k=0; k<np+1; k++)
451 for (
int j=0; j<nt; j++) {
452 int l_q, m_q, base_r ;
454 if (nullite_plm(j, nt, k, np, base) == 1) {
455 donne_lm(1, lz, j, k, base, m_q, l_q,base_r) ;
457 a_tilde[l_q] = (*coefftilde)(0, k, j, 0);
462 Scalar zeta2(map); zeta2= 3.; zeta2.std_spectral_base(); zeta2.mult_cost();
463 zeta2.set_spectral_va().ylm();
464 Base_val base2 = zeta2.get_spectral_base() ;
465 Mtbl_cf *dzeta = zeta2.set_spectral_va().c_cf;
467 for (
int k=0; k<np+1; k++)
468 for (
int j=0; j<nt; j++) {
469 int l_q2, m_q2, base_r2 ;
471 if (nullite_plm(j, nt, k, np, base2) == 1) {
472 donne_lm(1, lz, j, k, base2, m_q2, l_q2,base_r2) ;
475 (*dzeta).set(0,k,j,0) = a_tilde[0] + (a_tilde[1]/3.)*
sqrt(2.*1. + 1.);
478 (*dzeta).set(0,k,j,0) =
479 (a_tilde[l_q2 +1]/(2.*l_q2 + 3.))*
sqrt((2.*(l_q2 +1.)+1.)/(2.*l_q2 + 1.))
480 - (a_tilde[l_q2 -1]/(2.*l_q2 - 1.))
481 *
sqrt((2.*(l_q2 - 1.)+1.)/(2.*l_q2 + 1.));
486 zeta2.set_spectral_va().coef();
498 for (
int l=1; l<4; l++)
499 for (
int k=0; k<np; k++)
500 for (
int j=0; j<nt; j++) {
501 mp_rad->
val_lx_jk(h_surf.val_grid_point(0, k, j, 0)*1.000000000001, j, k, pipo,
504 ll3(l).get_spectral_va().val_point_jk(lz, xi, j, k) ;
513 Tensor jj3d = Kij - ss3d*ll3d - ll3d*ss3d
514 -
contract(Kij, 0 , 1, sxss3dcon , 0, 1)* sxss3d ;
518 for (
int l=1; l<4; l++)
519 for (
int m=1; m<4; m++)
520 for (
int k=0; k<np; k++)
521 for (
int j=0; j<nt; j++) {
522 mp_rad->
val_lx_jk(h_surf.val_grid_point(0, k, j, 0)*1.000000000000001,
523 j, k, pipo, lz, xi) ;
525 jj3(l,m).get_spectral_va().val_point_jk(lz, xi, j, k) ;
531 contract(proj_prov, 0,ss3d.derive_cov(gamij),1), 1 ) ;
535 for (
int l=1; l<4; l++)
536 for (
int m=1; m<4; m++)
537 for (
int k=0; k<np; k++)
538 for (
int j=0; j<nt; j++) {
539 mp_rad->
val_lx_jk(h_surf.val_grid_point(0, k, j, 0)*1.000000000000001,
540 j, k, pipo, lz, xi) ;
542 hh3(l,m).get_spectral_va().val_point_jk(lz, xi, j, k) ;
547 Tensor hh3dupdown = hh3d.
up(0, gamij);
553 Tensor hh3dupup = hh3dupdown.
up(1,gamij);
558 Scalar ricci22 = ricciscal3
564 ricci22 += (hh3d.
trace(gamij)*hh3d.
trace(gamij))
572 for (
int k=0; k<np; k++)
573 for (
int j=0; j<nt; j++) {
574 mp_rad->
val_lx_jk(h_surf.val_grid_point(0, k, j, 0)*1.000000000000001, j, k, pipo,
657 if (p_epsilon_A_minus_one != 0x0)
delete p_epsilon_A_minus_one ;
662 if (p_delta != 0x0)
delete p_delta ;
673 p_epsilon_A_minus_one = 0x0;
738 double rayon =
sqrt(
area()/(4.*M_PI));
750 double rayon =
sqrt(
area()/(4.*M_PI));
752 double factor =
mass()/(8. * M_PI);
754 {
for (
int compte=0; compte <=order -1; compte++)
755 factor = factor*rayon;
765 for (
int nn=1; nn<order; nn++){
767 Scalar Pnnew = (2.*nn +1.)*Pn;
769 Pnnew = Pnnew - nn*Pnold;
770 Pnnew = Pnnew/(double(nn) + 1.);
807 double rayon =
sqrt(
area()/(4.*M_PI));
809 double factor = 1./(8. * M_PI);
811 {
for (
int compte=0; compte <=order -2; compte++)
812 factor = factor*rayon;
825 for (
int nn=1; nn<order; nn++){
827 Scalar Pnnew = (2.*nn +1.)*Pn;
829 Pnnew = Pnnew - nn*Pnold;
830 Pnnew = Pnnew/(double(nn) + 1.);
839 dPn = Pn* zeta; dPn = dPn - Pnold; dPn = double(order)*dPn;
842 quotient = quotient*zeta*zeta; quotient = quotient -1.;
865 if (p_epsilon_A_minus_one == 0x0) {
870 return *p_epsilon_A_minus_one;
960 int valence1 = valence0 + 1 ;
961 int valence1m1 = valence1 - 1 ;
982 Itbl tipe(valence1) ;
984 for (
int id = 0;
id<valence0;
id++) {
985 tipe.
set(
id) = tipeuu(
id) ;
987 tipe.
set(valence1m1) = COV ;
1004 Itbl ind1(valence1) ;
1005 Itbl ind0(valence0) ;
1006 Itbl ind(valence0) ;
1026 for (
int ic=0; ic<ncomp1; ic++) {
1035 for (
int id = 0;
id < valence0;
id++) {
1036 ind0.
set(
id) = ind1(
id) ;
1040 int k = ind1(valence1m1) ;
1056 cresu = (uu(ind0)).srdsdt() ;
1059 for (
int id=0;
id<valence0;
id++) {
1061 switch ( ind0(
id) ) {
1093 cerr <<
"Connection_fspher::p_derive_cov : index problem ! " 1107 cresu = (uu(ind0)).srstdsdp() ;
1110 for (
int id=0;
id<valence0;
id++) {
1112 switch ( ind0(
id) ) {
1130 tmp.div_r_dzpuis(dz_resu) ;
1150 tmp.div_r_dzpuis(dz_resu) ;
1159 cerr <<
"Connection_fspher::p_derive_cov : index problem ! \n" 1171 cerr <<
"Connection_fspher::p_derive_cov : index problem ! \n" ;
1188 cout <<
"c'est pas fait!" << endl ;
1205 if (p_delta == 0x0) {
1209 christoflat.set_etat_zero() ;
1214 for (
int k=1; k<=3; k++) {
1215 for (
int i=1; i<=3; i++) {
1216 for (
int j=1; j<=i; j++) {
1218 for (
int l=1; l<=3; l++) {
1219 cc += qab.
con()(k,l) * (
1220 dgam(l,j,i) + dgam(i,l,j) - dgam(i,j,l) ) ;
1228 p_delta =
new Tensor (christoflat) ;
1245 for (
int y=1; y<=nbboucle; y++){
Sym_tensor * p_shear
The shear tensor.
Scalar h_surf
The location of the 2-surface as r = h_surf .
const Base_val & get_spectral_base() const
Returns the spectral bases of the Valeur va.
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
const Scalar & theta_plus() const
Computes the outgoing null expansion .
Metric for tensor calculation.
double * p_epsilon_P_minus_one
Refined Penrose parameter, difference wrt one.
double multipole_mass(const int order) const
Computes the mass multipole of a given order for the spheroid, assumed to be spherical.
double * p_area
The area of the 2-surface.
double angu_mom() const
Computes the angular momentum with respect to a divergence-free vector field tangent to the 2-surface...
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Tensor jac2d
The jacobian of the adaptation of coordinates (contravariant/covariant representation) ...
const Scalar & sqrt_q() const
Computes the normal vector field to the 2-surface.
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
int & set(int i)
Read/write of a particular element (index i ) (1D case)
double integrale_surface(const Cmp &ci, double rayon) const
Performs the surface integration of ci on the sphere of radius rayon .
virtual void del_deriv() const
Deletes all the derived quantities.
double * p_angu_mom
The angular momentum.
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Cmp sqrt(const Cmp &)
Square root.
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
double multipole_angu(const int order) const
Computes the angular multipole of a given order for the spheroid, assumed to be spherical.
void give_quant_numbers(int, int, int, int &, int &, int &) const
Computes the various quantum numbers and 1d radial base.
int sym_index1() const
Number of the first symmetric index (0<= id_sym1 < valence )
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Tensor up(int ind, const Metric &gam) const
Computes a new tensor by raising an index of *this.
void ylm()
Computes the coefficients of *this.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Tensor field of valence 0 (or component of a tensorial field).
Base class for coordinate mappings.
Vector ephi
The conformal Killing vector field on the 2-surface (set to by default to the axial vector associated...
int get_n_comp() const
Returns the number of stored components.
int sym_index2() const
Number of the second symmetric index (id_sym1 < id_sym2 < valence )
Basic integer array class.
const Sym_tensor & shear() const
Computes the shear of the 2-surface .
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field. ...
const Sym_tensor & ricci() const
Returns the Ricci tensor (given by the Connection p_connect )
virtual void allocate_all()
Sets the logical state to ETATQCQ (ordinary state) and performs the memory allocation of all the elem...
const Vector & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of *this with respect to some metric , by raising the index of...
Tensor field of valence 1.
Scalar * p_theta_minus
Null ingoing expansion.
const Map & get_mp() const
Returns the mapping.
double val_grid_point(int l, int k, int j, int i) const
Returns the value of the field at a specified grid point.
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
void annule_hard()
Sets the Scalar to zero in a hard way.
Spheroidal 2-surfaces embedded in a time-slice of the 3+1 formalism.
Spheroid(const Map_af &map, double radius)
The delta tensorial fields linked to Christoffel symbols.
const Vector & get_ephi() const
Returns the conformal Killing symmetry vector on the 2-surface.
Tensor derive_cov2dflat(const Tensor &uu) const
Computes the round covariant derivative on the spheroid.
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Tensor derive_cov2d(const Tensor &uu) const
Computes the total covariant derivative on the spheroid.
int get_index_type(int i) const
Gives the type (covariant or contravariant) of the index number i .
virtual void allocate_all()
Performs the memory allocation of all the elements, down to the double arrays of the Tbl s...
bool issphere
Flag to know whether the horizon is geometrically round or distorted.
Scalar trk
Trace K of the extrinsic curvature of the 3-slice.
Base class for pure radial mappings.
Vector ll
Normal-tangent component of the extrinsic curvature of the 3-slice.
double * p_multipole_angu
Angular momentum multipole for the spheroid.
void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
int get_nzone() const
Returns the number of domains.
void operator=(const Spheroid &)
Assignment to another Spheroid.
virtual void sauve(FILE *) const
Save in a file.
double epsilon_P_minus_one() const
Computation of the classical Penrose parameter, and its difference wrt one.
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...
const Scalar & get_ricci() const
Returns the 2-ricci scalar .
double area() const
Computes the area of the 2-surface.
Cmp pow(const Cmp &, int)
Power .
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Sym_tensor jj
Tangent components of the extrinsic curvature of the 3-slice.
int get_valence() const
Returns the valence.
int & set_index_type(int i)
Sets the type of the index number i .
virtual void val_lx_jk(double rr, int j, int k, const Param &par, int &l, double &xi) const =0
Computes the domain index l and the value of corresponding to a point of arbitrary r but collocation...
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Scalar ggg
Normalization function for the ingoing null vector k.
Bases of the spectral expansions.
double & set_grid_point(int l, int k, int j, int i)
Setting the value of the field at a given grid point.
const Scalar & srdsdt() const
Returns of *this .
const Scalar & theta_minus() const
Computes the ingoing null expansion .
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Sym_tensor hh
The ricci scalar on the 2-surface.
const Sym_tensor & get_qq() const
returns the 3-d degenerate 2-metric
Scalar fff
Normalization function for the outgoing null vector l.
Coefficients storage for the multi-domain spectral method.
double * p_mass
Mass defined from angular momentum.
virtual void set_etat_zero()
Sets the logical state of all components to ETATZERO (zero state).
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Valeur & set_spectral_va()
Returns va (read/write version)
Scalar & set(int)
Read/write access to a component.
Symmetric tensors (with respect to two of their arguments).
Vector ss
The adapted normal vector field to spheroid in the 3-slice.
double * p_multipole_mass
Mass multipole for the spheroid.
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
Tensor proj
The 3-d projector on the 2-surface (contravariant-covariant form).
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this.
const Map & get_mp() const
Returns the mapping.
const Mg3d * get_angu_1dom() const
Returns the pointer on the associated mono-domain angular grid.
virtual void dec_dzpuis(int dec=1)
Decreases by dec units the value of dzpuis and changes accordingly the values of the Scalar in the co...
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Scalar * p_theta_plus
Classical Penrose parameter, difference wrt one.
Scalar ricci
Induced metric on the 2-surface .
Class intended to describe valence-2 symmetric tensors.
double epsilon_A_minus_one() const
Computation of the refined Penrose parameter for axisymmetric spacetimes, and its difference wrt one...
Scalar * p_sqrt_q
Surface element .
const Valeur & get_spectral_va() const
Returns va (read only version)
Sym_tensor qq
The 3-d covariant degenerated 2-metric on the surface.
virtual Itbl indices(int pos) const
Returns the indices of a component given by its position in the array cmp .
const Scalar & srstdsdp() const
Returns of *this .
double mass() const
Computes the mass as defined from the calculus of angular momentum, done with respect to a divergence...
Tensor trace(int ind1, int ind2) const
Trace on two different type indices.
Coord r
r coordinate centered on the grid
virtual ~Spheroid()
Destructor.
const Tensor & delta() const
Computes the delta coefficients for covariant derivative.
Tensor down(int ind, const Metric &gam) const
Computes a new tensor by lowering an index of *this.