26char 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++) {
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());
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) ;
315 ephi3.
set(3) = ephi33;
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;
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) ;
360 + 2* gamij.
cov()(1,2)* (h_surf3.
srdsdt()) + gamij.
cov()(2,2);
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) ;
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);
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.));
498 for (
int l=1; l<4; l++)
499 for (
int k=0; k<np; k++)
500 for (
int j=0; j<nt; j++) {
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++) {
523 j, k, pipo, lz, xi) ;
525 jj3(l,m).get_spectral_va().val_point_jk(lz, xi, j, k) ;
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++) {
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++) {
606 issphere(sph_in.issphere)
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) ) {
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) {
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++){
Bases of the spectral expansions.
void give_quant_numbers(int, int, int, int &, int &, int &) const
Computes the various quantum numbers and 1d radial base.
Basic integer array class.
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 .
Base class for pure radial mappings.
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...
Base class for coordinate mappings.
Coord r
r coordinate centered on the grid
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Metric for tensor calculation.
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
const Map & get_mp() const
Returns the mapping.
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
const Sym_tensor & ricci() const
Returns the Ricci tensor (given by the Connection p_connect )
const Mg3d * get_angu_1dom() const
Returns the pointer on the associated mono-domain angular grid.
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.
Coefficients storage for the multi-domain spectral method.
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
void div_r_dzpuis(int ced_mult_r)
Division by r everywhere but with the output flag dzpuis set to ced_mult_r .
void mult_cost()
Multiplication by .
const Scalar & srstdsdp() const
Returns of *this .
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field.
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 set_etat_zero()
Sets the logical state to ETATZERO (zero).
void div_tant()
Division by .
Valeur & set_spectral_va()
Returns va (read/write version)
const Valeur & get_spectral_va() const
Returns va (read only version)
void annule_hard()
Sets the Scalar to zero in a hard way.
void mult_sint()
Multiplication by .
virtual void allocate_all()
Sets the logical state to ETATQCQ (ordinary state) and performs the memory allocation of all the elem...
const Base_val & get_spectral_base() const
Returns the spectral bases of the Valeur va
double & set_grid_point(int l, int k, int j, int i)
Setting the value of the field at a given grid point.
const Vector & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of *this with respect to some metric , by raising the index of...
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...
void mult_r()
Multiplication by r everywhere; dzpuis is not changed.
Spheroidal 2-surfaces embedded in a time-slice of the 3+1 formalism.
Tensor derive_cov2dflat(const Tensor &uu) const
Computes the round covariant derivative on the spheroid.
virtual ~Spheroid()
Destructor.
const Scalar & sqrt_q() const
Computes the normal vector field to the 2-surface.
Tensor proj
The 3-d projector on the 2-surface (contravariant-covariant form).
Sym_tensor qq
The 3-d covariant degenerated 2-metric on the surface.
const Scalar & theta_minus() const
Computes the ingoing null expansion .
virtual void sauve(FILE *) const
Save in a file.
Vector ll
Normal-tangent component of the extrinsic curvature of the 3-slice.
double * p_multipole_angu
Angular momentum multipole for the spheroid.
double * p_mass
Mass defined from angular momentum.
Scalar trk
Trace K of the extrinsic curvature of the 3-slice.
Scalar * p_theta_plus
Classical Penrose parameter, difference wrt one.
Scalar ggg
Normalization function for the ingoing null vector k.
void operator=(const Spheroid &)
Assignment to another Spheroid.
double epsilon_A_minus_one() const
Computation of the refined Penrose parameter for axisymmetric spacetimes, and its difference wrt one.
void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
Spheroid(const Map_af &map, double radius)
The delta tensorial fields linked to Christoffel symbols.
const Tensor & delta() const
Computes the delta coefficients for covariant derivative.
double multipole_angu(const int order) const
Computes the angular multipole of a given order for the spheroid, assumed to be spherical.
bool issphere
Flag to know whether the horizon is geometrically round or distorted.
Tensor jac2d
The jacobian of the adaptation of coordinates (contravariant/covariant representation)
double mass() const
Computes the mass as defined from the calculus of angular momentum, done with respect to a divergence...
Sym_tensor jj
Tangent components of the extrinsic curvature of the 3-slice.
double * p_area
The area of the 2-surface.
Scalar fff
Normalization function for the outgoing null vector l.
double area() const
Computes the area of the 2-surface.
double multipole_mass(const int order) const
Computes the mass multipole of a given order for the spheroid, assumed to be spherical.
Scalar h_surf
The location of the 2-surface as r = h_surf .
Sym_tensor hh
The ricci scalar on the 2-surface.
double * p_epsilon_P_minus_one
Refined Penrose parameter, difference wrt one.
Scalar * p_sqrt_q
Surface element .
Vector ss
The adapted normal vector field to spheroid in the 3-slice.
virtual void del_deriv() const
Deletes all the derived quantities.
const Scalar & get_ricci() const
Returns the 2-ricci scalar .
const Sym_tensor & shear() const
Computes the shear of the 2-surface .
Scalar * p_theta_minus
Null ingoing expansion.
double angu_mom() const
Computes the angular momentum with respect to a divergence-free vector field tangent to the 2-surface...
const Sym_tensor & get_qq() const
returns the 3-d degenerate 2-metric
Vector ephi
The conformal Killing vector field on the 2-surface (set to by default to the axial vector associated...
const Vector & get_ephi() const
Returns the conformal Killing symmetry vector on the 2-surface.
double * p_multipole_mass
Mass multipole for the spheroid.
double * p_angu_mom
The angular momentum.
Tensor derive_cov2d(const Tensor &uu) const
Computes the total covariant derivative on the spheroid.
Scalar ricci
Induced metric on the 2-surface .
Sym_tensor * p_shear
The shear tensor.
double epsilon_P_minus_one() const
Computation of the classical Penrose parameter, and its difference wrt one.
const Scalar & theta_plus() const
Computes the outgoing null expansion .
Class intended to describe valence-2 symmetric tensors.
Symmetric tensors (with respect to two of their arguments).
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...
void ylm()
Computes the coefficients of *this.
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
void coef() const
Computes the coeffcients of *this.
Tensor field of valence 1.
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Scalar & set(int)
Read/write access to a component.
Cmp sqrt(const Cmp &)
Square root.
Cmp pow(const Cmp &, int)
Power .
const Map & get_mp() const
Returns the mapping.
int sym_index1() const
Number of the first symmetric index (0<= id_sym1 < valence )
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.
int sym_index2() const
Number of the second symmetric index (id_sym1 < id_sym2 < valence )
int get_valence() const
Returns the valence.
virtual Itbl indices(int pos) const
Returns the indices of a component given by its position in the array cmp .
int & set_index_type(int i)
Sets the type of the index number i .
int get_n_comp() const
Returns the number of stored components.
Tensor down(int ind, const Metric &gam) const
Computes a new tensor by lowering an index of *this.
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
const Tensor & derive_cov(const Metric &gam) const
Returns the covariant derivative of this with respect to some metric .
Tensor up(int ind, const Metric &gam) const
Computes a new tensor by raising an index of *this.
Tensor trace(int ind1, int ind2) const
Trace on two different type indices.
virtual void set_etat_zero()
Sets the logical state of all components to ETATZERO (zero state).
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .