LORENE
spheroid.C
1 /*
2  * Definition of methods for the class Spheroid and its subclass App_hor
3  *
4  */
5 
6 /*
7  * Copyright (c) 2006 Nicolas Vasset, Jerome Novak & Jose-Luis Jaramillo
8  *
9  * This file is part of LORENE.
10  *
11  * LORENE is free software; you can redistribute it and/or modify
12  * it under the terms of the GNU General Public License version 2
13  * as published by the Free Software Foundation.
14  *
15  * LORENE is distributed in the hope that it will be useful,
16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18  * GNU General Public License for more details.
19  *
20  * You should have received a copy of the GNU General Public License
21  * along with LORENE; if not, write to the Free Software
22  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
23  *
24  */
25 
26 
27 
28 /*
29  * $Id: spheroid.C,v 1.22 2016/12/05 16:17:44 j_novak Exp $
30  * $Log: spheroid.C,v $
31  * Revision 1.22 2016/12/05 16:17:44 j_novak
32  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
33  *
34  * Revision 1.21 2014/10/13 08:52:38 j_novak
35  * Lorene classes and functions now belong to the namespace Lorene.
36  *
37  * Revision 1.20 2014/10/06 15:12:56 j_novak
38  * Modified #include directives to use c++ syntax.
39  *
40  * Revision 1.19 2012/05/18 16:27:05 j_novak
41  * ... and another memory leak
42  *
43  * Revision 1.18 2012/05/18 15:19:14 j_novak
44  * Corrected a memory leak
45  *
46  * Revision 1.17 2012/05/11 14:11:24 j_novak
47  * Modifications to deal with the accretion of a scalar field
48  *
49  * Revision 1.16 2009/12/01 08:44:24 j_novak
50  * Added the missing operator=().
51  *
52  * Revision 1.15 2008/12/10 13:55:55 jl_jaramillo
53  * versions developed at Meudon in November 2008
54  *
55  * Revision 1.14 2008/11/17 08:30:01 jl_jaramillo
56  * Nicolas's changes on multipoles. Sign correction in outgoing shear expression
57  *
58  * Revision 1.13 2008/11/12 15:17:47 n_vasset
59  * Bug-hunting, and new definition for computation of the Ricci scalar
60  * (instead of Ricci tensor previously)
61  *
62  * Revision 1.12 2008/07/09 08:47:33 n_vasset
63  * new version for multipole calculation. Function zeta implemented.
64  * Revision 1.11 2008/06/04 12:31:23 n_vasset
65  * New functions multipole_mass and multipole_angu. first version.
66  *
67  * Revision 1.9 2006/09/07 08:39:45 j_novak
68  * Minor changes.
69  *
70  * Revision 1.8 2006/07/03 10:13:48 n_vasset
71  * More efficient method for calculation of ricci tensor. Adding of flag issphere
72  *
73  * Revision 1.4 2006/06/01 11:47:50 j_novak
74  * Memory error hunt.
75  *
76  * Revision 1.3 2006/06/01 08:18:16 n_vasset
77  * Further implementation of Spheroid class definitions
78  *
79  * $Header: /cvsroot/Lorene/C++/Source/App_hor/spheroid.C,v 1.22 2016/12/05 16:17:44 j_novak Exp $
80  *
81  */
82 
83 // C headers
84 #include <cmath>
85 
86 // Lorene headers
87 #include "metric.h"
88 #include "spheroid.h"
89 #include "proto.h"
90 #include "param.h"
91 
92 //---------------//
93 // Constructors //
94 //--------------//
95 
96 namespace Lorene {
97 Spheroid::Spheroid(const Map_af& map, double radius):
98  h_surf(map),
99  jac2d(map, 2, COV, map.get_bvect_spher()),
100  proj(map, 2, COV, map.get_bvect_spher()),
101  qq(map, COV, map.get_bvect_spher()),
102  ss (map, COV, map.get_bvect_spher()),
103  ephi(map, CON, map.get_bvect_spher()),
104  qab(map.flat_met_spher()),
105  ricci(map),
106  hh(map, COV, map.get_bvect_spher()),
107  trk(map),
108  ll(map, COV, map.get_bvect_spher()),
109  jj(map, COV, map.get_bvect_spher()),
110  fff(map),
111  ggg(map),
112  zeta(map),
113  issphere(true)
114 {
115 
116  // Itbl type (2) ;
117  // type.set(1) = CON ;
118  // type.set(2) = COV ;
119  // Tensor proj(map, 2, type, map.get_bvect_spher());
120 
121 
122  assert(radius > 1.e-15) ;
123  assert(map.get_mg()->get_nzone() == 1) ; // one domain
124  assert(map.get_mg()->get_nr(0) == 1) ; // angular grid
125  assert(map.get_mg()->get_type_r(0) == FIN) ; //considered as a shell
126 
127 
128  // Setting of real index types for jacobian and projector (first contravariant, second covariant)
129  jac2d.set_index_type(0) = CON ;
130  proj.set_index_type(0) = CON ;
131 
132  jac2d.set_etat_zero() ;
133 
134 
135  h_surf = radius ;
136  ss.set_etat_zero();
139  hh.set_etat_zero() ;
140  for (int i=1; i<=3; i++)
141  hh.set(i,i) = 2./radius ;
142  trk.set_etat_zero() ;
143  ll.set_etat_zero() ;
144  jj.set_etat_zero() ;
145  fff.set_etat_zero();
146  ggg.set_etat_zero();
147  zeta.set_etat_zero();
148  set_der_0x0() ;
149 
150 
151 }
152 
153 Spheroid::Spheroid(const Scalar& h_in, const Metric& gamij, const Sym_tensor& Kij):
154  h_surf(h_in),
155  jac2d(h_in.get_mp(),2, COV, h_in.get_mp().get_bvect_spher()),
156  proj(h_in.get_mp(),2, COV, h_in.get_mp().get_bvect_spher()),
157  qq(h_in.get_mp(), COV, h_in.get_mp().get_bvect_spher()),
158  ss (h_in.get_mp(), COV, h_in.get_mp().get_bvect_spher()),
159  ephi(h_in.get_mp(), CON, h_in.get_mp().get_bvect_spher()),
160  qab(h_in.get_mp().flat_met_spher()),
161  ricci(h_in.get_mp()),
162  hh(h_in.get_mp(), COV, h_in.get_mp().get_bvect_spher()),
163  trk(h_in.get_mp()),
164  ll(h_in.get_mp(), COV, h_in.get_mp().get_bvect_spher()),
165  jj(h_in.get_mp(), COV, h_in.get_mp().get_bvect_spher()),
166  fff(h_in.get_mp()),
167  ggg(h_in.get_mp()),
168  zeta(h_in.get_mp()),
169  issphere(true) {
170 
171  set_der_0x0() ;
172  const Map& map = h_in.get_mp() ; // The 2-d 1-domain map
173 
174  const Map& map3 = Kij.get_mp(); // The 3-d map
175  const Mg3d& gri2d = *map.get_mg() ;
176 
177  const Map_radial* mp_rad = dynamic_cast<const Map_radial*>(&Kij.get_mp()) ;
178  assert(mp_rad != 0x0) ;
179 
180  const Mg3d& gri3d = *map3.get_mg();
181  assert(&gri2d == Kij.get_mp().get_mg()->get_angu_1dom()) ;
182 
183  int np = gri2d.get_np(0) ;
184  int nt = gri2d.get_nt(0) ;
185 
186  int nr3 = gri3d.get_nr(0) ;
187  int nt3 = gri3d.get_nt(0) ;
188  int np3 = gri3d.get_np(0) ;
189  int nz = gri3d.get_nzone() ;
190  assert( nt == nt3 ) ;
191  assert ( np == np3 );
192  assert ( nr3 != 1 );
193 
194  Param pipo ;
195  double xi = 0. ;
196  int lz = 0 ;
197 
198  if(nz >2){
199  lz =1;
200  }
201 
202  // Setting of real index types forjacobian and projector(first contravariant, other covariant)
203  proj.set_index_type(0) = CON;
204  jac2d.set_index_type(0) = CON;
205 
206  // Copy of the 2-dimensional h_surf to a 3_d h_surf (calculus commodity, in order to match grids)
207  Scalar h_surf3 (map3);
208 
209  h_surf3.allocate_all();
210  h_surf3.std_spectral_base();
211  for (int f= 0; f<nz; f++)
212  for (int k=0; k<np3; k++)
213  for (int j=0; j<nt3; j++) {
214  for (int l=0; l<nr3; l++) {
215  h_surf3.set_grid_point(f,k,j,l) = h_surf.val_grid_point(0,k,j,0);
216  }
217  }
218  if (nz >2){
219  h_surf3.annule_domain(0);
220  h_surf3.annule_domain(nz - 1);
221  }
222  // h_surf3.std_spectral_base();
223 
224  /* Computation of the jacobian projector linked to the mapping from the
225  spheroid to a coordinate sphere. All quantities will then be calculated
226  as from a real coordinate sphere
227  */
228  Itbl type (2);
229  type.set(0) = CON ;
230  type.set(1) = COV ;
231  Tensor jac (Kij.get_mp(),2,type, Kij.get_mp().get_bvect_spher());
232  jac.set_etat_zero();
233  jac.std_spectral_base();
234  jac.set(1,1) = 1. ;
235  jac.set(2,2)= 1. ;
236  jac.set(3,3) = 1. ;
237  jac.set(1,2) = -h_surf3.srdsdt() ;
238  jac.set(1,3) = -h_surf3.srstdsdp() ;
239  jac.std_spectral_base() ;
240 
241  // Copy on the 2-d grid
242  jac2d.allocate_all() ;
244  for (int l=1; l<4; l++)
245  for (int m=1; m<4; m++)
246  for (int k=0; k<np; k++)
247  for (int j=0; j<nt; j++) {
248  mp_rad->val_lx_jk((h_surf.val_grid_point(0, k, j, 0))*1.0000000000001,
249  j, k, pipo, lz, xi) ;
250  jac2d.set(l,m).set_grid_point(0, k, j, 0) =
251  jac(l,m).get_spectral_va().val_point_jk(lz, xi, j, k) ;
252  }
253 
254  // Inverse jacobian (on 3-d grid)
255  Tensor jac_inv = jac ;
256  jac_inv.set(1,2) = - jac_inv(1,2);
257  jac_inv.set(1,3) = - jac_inv(1,3) ;
258 
259  // Scalar field which annulation characterizes the 2-surface
260  Scalar carac = Kij.get_mp().r - h_surf3;
261  carac.std_spectral_base();
262  // Computation of the normal vector (covariant form) on both grids
263  Vector ss3d= carac.derive_cov(gamij) ;
264  Vector ss3dcon= carac.derive_con(gamij) ;
265  Scalar ssnorm = contract (ss3d.up(0, gamij), 0, ss3d, 0);
266  ssnorm.std_spectral_base() ;
267  ss3d = ss3d / sqrt (ssnorm) ;
268  ss3dcon = ss3dcon / sqrt (ssnorm) ;
269  if (nz >2){
270  ss3d.annule_domain(0);
271  ss3dcon.annule_domain(0);
272  ss3d.annule_domain(nz-1);
273  ss3dcon.annule_domain(nz -1);
274  }
275  ss3d.std_spectral_base();
276  ss3dcon.std_spectral_base();
277 
278 
279  // Provisory handling of dzpuis problems
280  if (nz >2){
281  h_surf3.annule_domain(nz-1);
282  }
283  for (int ii=1; ii <=3; ii++){
284  ss3d.set(ii).dec_dzpuis(ss3d(ii).get_dzpuis());
285  ss3dcon.set(ii).dec_dzpuis(ss3dcon(ii).get_dzpuis());
286  }
287  for (int ii=1; ii <=3; ii++)
288  for (int jjy = 1; jjy <=3; jjy ++){
289  jac_inv.set(ii, jjy).dec_dzpuis(jac_inv(ii, jjy).get_dzpuis());
290  jac.set(ii, jjy).dec_dzpuis(jac(ii, jjy).get_dzpuis());
291  }
292  // End of dzpuis handling.
293 
294  Sym_tensor sxss3d = ss3d * ss3d ;
295  Sym_tensor sxss3dcon = ss3dcon * ss3dcon ;
296  Vector ss3 (Kij.get_mp(), COV, Kij.get_mp().get_bvect_spher()) ;
297  Vector ss3con(Kij.get_mp(), CON, Kij.get_mp().get_bvect_spher()) ;
298  ss3 = contract(jac_inv, 0, ss3d, 0);
299  ss.allocate_all() ;
301 
302  for (int l=1; l<4; l++)
303  for (int k=0; k<np; k++)
304  for (int j=0; j<nt; j++) {
305  mp_rad->val_lx_jk(h_surf.val_grid_point(0, k, j, 0)*1.0000000000001, j, k, pipo,
306  lz, xi) ;
307  ss.set(l).set_grid_point(0, k, j, 0) =
308  ss3(l).get_spectral_va().val_point_jk(lz, xi, j, k) ;
309  }
310 
311  // The vector field associated with Kiling conformal symmetry
312 
313  Vector ephi3(gamij.get_mp(), CON, gamij.get_mp().get_bvect_spher());
314  ephi3.set(1).annule_hard();
315  ephi3.set(2).annule_hard();
316  Scalar ephi33(gamij.get_mp()); ephi33 = 1.; ephi33.std_spectral_base();
317  ephi33.mult_r(); ephi33.mult_sint();
318  ephi3.set(3) = ephi33;
319  ephi3.std_spectral_base();
320  ephi3 = contract (jac, 1, ephi3,0);
321  ephi.allocate_all() ;
323 
324  for (int l=1; l<4; l++)
325  for (int k=0; k<np; k++)
326  for (int j=0; j<nt; j++) {
327  mp_rad->val_lx_jk(h_surf.val_grid_point(0, k, j, 0)*1.000000000000001,
328  j, k, pipo, lz, xi) ;
329  ephi.set(l).set_grid_point(0, k, j, 0) =
330  ephi3(l).get_spectral_va().val_point_jk(lz, xi, j, k) ;
331  }
332 
333  // Computation of the 3-d projector on the 2-sphere
334  Tensor proj3d(Kij.get_mp(),2, type, Kij.get_mp().get_bvect_spher());
335  Tensor proj_prov = gamij.con().down(1, gamij) - ss3dcon*ss3d;
336  proj.allocate_all();
338  proj3d.allocate_all();
339  proj3d = contract(jac, 1, contract( jac_inv, 0, proj_prov , 1), 1 );
340  proj3d.std_spectral_base();
341 
342  for (int l=1; l<4; l++)
343  for (int m=1; m<4; m++)
344  for (int k=0; k<np; k++)
345  for (int j=0; j<nt; j++) {
346  mp_rad->val_lx_jk(h_surf.val_grid_point(0, k, j, 0)*1.00000000000001,
347  j, k, pipo, lz, xi) ;
348  proj.set(l,m).set_grid_point(0, k, j, 0) =
349  proj3d(l,m).get_spectral_va().val_point_jk(lz, xi, j, k) ;
350  }
351 
352  /* Computation of the metric linked to the 2-surface (linked to covariant form
353  of the degenerated 2-metric).
354  It will be designed as a block-diagonal 3-metric, with 1 for
355  the first coordinate and the concerned 2-d metric as a
356  second block */
357 
358  Sym_tensor qq3d (Kij.get_mp(), COV, Kij.get_mp().get_bvect_spher());
359  Sym_tensor qab2 (h_in.get_mp(), COV, h_in.get_mp().get_bvect_spher());
360  qq3d.set_etat_zero();
361  qq3d.set(1,1) = 1;
362  qq3d.set(2,2)= gamij.cov()(1,1) * (h_surf3.srdsdt())* (h_surf3.srdsdt())
363  + 2* gamij.cov()(1,2)* (h_surf3.srdsdt()) + gamij.cov()(2,2);
364  qq3d.set(3,3)= gamij.cov()(1,1)* (h_surf3.srstdsdp())*(h_surf3.srstdsdp())
365  +2*gamij.cov()(1,3)* (h_surf3.srstdsdp()) +gamij.cov()(3,3);
366  qq3d.set(2,3)= gamij.cov()(1,1)* (h_surf3.srdsdt())* (h_surf3.srstdsdp())+
367  gamij.cov()(1,2)* (h_surf3.srstdsdp())+gamij.cov()(1,3)* (h_surf3.srdsdt())
368  + gamij.cov()(2,3) ;
369  qq3d.set(3,2)= gamij.cov()(1,1)* (h_surf3.srdsdt())* (h_surf3.srstdsdp())+
370  gamij.cov()(1,2)* (h_surf3.srstdsdp())+gamij.cov()(1,3)* (h_surf3.srdsdt())
371  + gamij.cov()(2,3) ;
372  qq3d.std_spectral_base();
373  qab2.allocate_all() ;
374  qab2.std_spectral_base();
375  for (int l=1; l<4; l++)
376  for (int m=1; m<4; m++)
377  for (int k=0; k<np; k++)
378  for (int j=0; j<nt; j++) {
379  mp_rad->val_lx_jk(h_surf.val_grid_point(0, k, j, 0)*1.000000000000001,
380  j, k, pipo, lz, xi) ;
381  qab2.set(l,m).set_grid_point(0, k, j, 0) =
382  qq3d(l,m).get_spectral_va().val_point_jk(lz, xi, j, k) ;
383  }
384  qab= qab2;
385 
386  // Computation of the degenerated 3d degenerated covariant metric on the 2-surface
387  Sym_tensor qqq = contract(jac_inv, 0, contract( jac_inv, 0, (gamij.cov()
388  - ss3d * ss3d) , 0), 1) ;
389  qq.allocate_all() ;
391  for (int l=1; l<4; l++)
392  for (int m=1; m<4; m++)
393  for (int k=0; k<np; k++)
394  for (int j=0; j<nt; j++) {
395  mp_rad->val_lx_jk(h_surf.val_grid_point(0, k, j, 0)*1.000000000000001,
396  j, k, pipo, lz, xi) ;
397  qq.set(l,m).set_grid_point(0, k, j, 0) =
398  qqq(l,m).get_spectral_va().val_point_jk(lz, xi, j, k) ;
399  }
400 
401  // Computation of the trace of the extrinsic curvature of 3-slice
402  Scalar trk3d = Kij.trace(gamij) ;
403  trk.allocate_all() ;
404  for (int k=0; k<np; k++)
405  for (int j=0; j<nt; j++) {
406  mp_rad->val_lx_jk(h_surf.val_grid_point(0, k, j, 0)*1.00000000000001,
407  j, k, pipo, lz, xi) ;
408  trk.set_grid_point(0, k, j, 0) =
409  trk3d.get_spectral_va().val_point_jk(lz, xi, j, k) ;
410  }
411 
412  // Computation of the normalization factor of the outgoing null vector.
413  Scalar fff3d (map3);
414  fff3d = 1. ;
415  fff.allocate_all() ;
416  fff3d.std_spectral_base();
417  for (int k=0; k<np; k++)
418  for (int j=0; j<nt; j++) {
419  mp_rad->val_lx_jk(h_surf.val_grid_point(0, k, j, 0)*1.000000000000001, j, k, pipo,
420  lz, xi) ;
421  fff.set_grid_point(0, k, j, 0) =
422  fff3d.get_spectral_va().val_point_jk(lz, xi, j, k) ;
423  }
424 
425  // Computation of the normalization factor of the ingoing null vector.
426  Scalar ggg3d (map3);
427  ggg3d = 1. ;
428  ggg.allocate_all() ;
429  ggg3d.std_spectral_base();
430  for (int k=0; k<np; k++)
431  for (int j=0; j<nt; j++) {
432  mp_rad->val_lx_jk(h_surf.val_grid_point(0, k, j, 0)*1.00000000000001, j, k, pipo,
433  lz, xi) ;
434  ggg.set_grid_point(0, k, j, 0) =
435  ggg3d.get_spectral_va().val_point_jk(lz, xi, j, k) ;
436  }
437 
438  //Computation of function zeta, changing spheroid coordinates to get a round measure.
439  Scalar ftilde = sqrt_q();
440  double rayon = sqrt(area()/(4.*M_PI));
441  ftilde = -ftilde/(rayon*rayon);
442  ftilde = ftilde*h_surf*h_surf;
443  ftilde.set_spectral_va().ylm();
444 
445  Base_val base = ftilde.get_spectral_base() ;
446 
447  Mtbl_cf *coefftilde = ftilde.get_spectral_va().c_cf;
448  int nombre = 2*nt; // ### Doubled in SYM base!!
449  double *a_tilde = new double[nombre];
450 
451  lz = 0; // Now we work with 2d map associated with sqrt(q)
452 
453  for (int k=0; k<np+1; k++)
454  for (int j=0; j<nt; j++) {
455  int l_q, m_q, base_r ;
456  base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
457  if (nullite_plm(j, nt, k, np, base) == 1) {
458  donne_lm(1, lz, j, k, base, m_q, l_q,base_r) ;
459  if (m_q ==0) {
460  a_tilde[l_q] = (*coefftilde)(0, k, j, 0);
461  }
462  }
463  }
464 
465  Scalar zeta2(map); zeta2= 3.; zeta2.std_spectral_base(); zeta2.mult_cost();
466  zeta2.set_spectral_va().ylm();
467  Base_val base2 = zeta2.get_spectral_base() ;
468  Mtbl_cf *dzeta = zeta2.set_spectral_va().c_cf;
469 
470  for (int k=0; k<np+1; k++)
471  for (int j=0; j<nt; j++) {
472  int l_q2, m_q2, base_r2 ;
473  base.give_quant_numbers(lz, k, j, m_q2, l_q2, base_r2) ;
474  if (nullite_plm(j, nt, k, np, base2) == 1) {
475  donne_lm(1, lz, j, k, base2, m_q2, l_q2,base_r2) ;
476  if (m_q2 ==0){
477  if(l_q2 ==0){
478  (*dzeta).set(0,k,j,0) = a_tilde[0] + (a_tilde[1]/3.)*sqrt(2.*1. + 1.);
479  }
480  if (l_q2 >0){
481  (*dzeta).set(0,k,j,0) =
482  (a_tilde[l_q2 +1]/(2.*l_q2 + 3.))*sqrt((2.*(l_q2 +1.)+1.)/(2.*l_q2 + 1.))
483  - (a_tilde[l_q2 -1]/(2.*l_q2 - 1.))
484  *sqrt((2.*(l_q2 - 1.)+1.)/(2.*l_q2 + 1.));
485  }
486  }
487  }
488  }
489  zeta2.set_spectral_va().coef();
490  zeta = zeta2;
491 
492  /* Computation of the tangent part of the extrinsic curvature of
493  * the 2 surface embedded in the 3 slice. The reference vector
494  used is the vector field s */
495 
496  Vector ll3d = contract( proj_prov, 0, contract(Kij, 1, ss3dcon, 0), 0) ;
497  Vector ll3 = contract( jac_inv, 0, ll3d, 0) ;
498 
499  ll.allocate_all() ;
501  for (int l=1; l<4; l++)
502  for (int k=0; k<np; k++)
503  for (int j=0; j<nt; j++) {
504  mp_rad->val_lx_jk(h_surf.val_grid_point(0, k, j, 0)*1.000000000001, j, k, pipo,
505  lz, xi) ;
506  ll.set(l).set_grid_point(0, k, j, 0) =
507  ll3(l).get_spectral_va().val_point_jk(lz, xi, j, k) ;
508  }
509 
510  /* computation of the tangent components of the extrinsic curvature
511  *of the 3-slice
512  *(extracted from the curvature of the timeslice)
513  * Note: this is not the actual 2d_ extrinsic curvature, but the
514  *tangent part of the time-slice extrinsic curvature */
515 
516  Tensor jj3d = Kij - ss3d*ll3d - ll3d*ss3d
517  - contract(Kij, 0 , 1, sxss3dcon , 0, 1)* sxss3d ;
518  Tensor jj3 =contract(jac_inv, 0 , contract(jac_inv,0 , jj3d,1),1);
519  jj.allocate_all() ;
521  for (int l=1; l<4; l++)
522  for (int m=1; m<4; m++)
523  for (int k=0; k<np; k++)
524  for (int j=0; j<nt; j++) {
525  mp_rad->val_lx_jk(h_surf.val_grid_point(0, k, j, 0)*1.000000000000001,
526  j, k, pipo, lz, xi) ;
527  jj.set(l,m).set_grid_point(0, k, j, 0) =
528  jj3(l,m).get_spectral_va().val_point_jk(lz, xi, j, k) ;
529  }
530 
531  // Computation of 2d extrinsic curvature in the 3-slice
532 
533  Tensor hh3d = contract(proj_prov, 0,
534  contract(proj_prov, 0,ss3d.derive_cov(gamij),1), 1 ) ;
535  Tensor hh3 =contract(jac_inv, 0 , contract(jac_inv,0 , hh3d,1),1);
536  hh.allocate_all() ;
538  for (int l=1; l<4; l++)
539  for (int m=1; m<4; m++)
540  for (int k=0; k<np; k++)
541  for (int j=0; j<nt; j++) {
542  mp_rad->val_lx_jk(h_surf.val_grid_point(0, k, j, 0)*1.000000000000001,
543  j, k, pipo, lz, xi) ;
544  hh.set(l,m).set_grid_point(0, k, j, 0) =
545  hh3(l,m).get_spectral_va().val_point_jk(lz, xi, j, k) ;
546  }
547 
548  // Computation of 2d ricci scalar
549 
550  Tensor hh3dupdown = hh3d.up(0, gamij);
551  Scalar ricciscal3 = gamij.ricci().trace(gamij);
552  if (nz>2){
553  // Ricci scalar on the 3-surface.
554  ricciscal3.annule_domain(nz-1); ricciscal3.std_spectral_base();
555  }
556  Tensor hh3dupup = hh3dupdown.up(1,gamij);
557  if (nz>2){
558  hh3dupup.annule_domain(nz-1); hh3dupup.std_spectral_base();
559  }
560 
561  Scalar ricci22 = ricciscal3
562  - 2.*contract(contract(gamij.ricci(), 1, ss3dcon, 0),0, ss3dcon, 0);
563  if (nz >2){
564  ricci22.annule_domain(nz-1);
565  ricci22.std_spectral_base();
566  }
567  ricci22 += (hh3d.trace(gamij)*hh3d.trace(gamij))
568  - contract(contract(hh3dupup,0, hh3d,0),0,1);
569  if (nz >2){
570  ricci22.annule_domain(nz-1);
571  }
572  ricci22.std_spectral_base();
575  for (int k=0; k<np; k++)
576  for (int j=0; j<nt; j++) {
577  mp_rad->val_lx_jk(h_surf.val_grid_point(0, k, j, 0)*1.000000000000001, j, k, pipo,
578  lz, xi) ;
579  ricci.set_grid_point(0, k, j, 0) =
580  ricci22.get_spectral_va().val_point_jk(lz, xi, j, k) ;
581  }
582 
583  del_deriv() ; //### to be checked...
584  set_der_0x0() ;
585  delete [] a_tilde ;
586 }
587 
588 
589 
590 
591 
592 //Copy constructor//
593 
594 Spheroid::Spheroid (const Spheroid &sph_in) :h_surf(sph_in.h_surf),
595  jac2d(sph_in.jac2d),
596  proj(sph_in.proj),
597  qq(sph_in.qq),
598  ss (sph_in.ss),
599  ephi (sph_in.ephi),
600  qab(sph_in.qab),
601  ricci(sph_in.ricci),
602  hh(sph_in.hh),
603  trk(sph_in.trk),
604  ll(sph_in.ll),
605  jj(sph_in.jj),
606  fff(sph_in.fff),
607  ggg(sph_in.ggg),
608  zeta(sph_in.zeta),
609  issphere(sph_in.issphere)
610 
611 {
612  set_der_0x0() ;
613 
614 }
615 
616 // Assignment to another Spheroid
617 void Spheroid::operator=(const Spheroid& sph_in)
618 {
619 
620  h_surf = sph_in.h_surf ;
621  jac2d = sph_in.jac2d ;
622  proj = sph_in.proj ;
623  qq = sph_in.qq ;
624  ss = sph_in.ss ;
625  ephi = sph_in.ephi ;
626  qab = sph_in.qab ;
627  ricci = sph_in.ricci ;
628  hh = sph_in.hh ;
629  trk = sph_in.trk ;
630  ll = sph_in.ll ;
631  jj = sph_in.jj ;
632  fff = sph_in.fff ;
633  ggg = sph_in.ggg ;
634  zeta = sph_in.zeta ;
635  issphere = sph_in.issphere ;
636 
637  del_deriv() ; // Deletes all derived quantities
638 
639 }
640 
641 //------------//
642 //Destructor //
643 //-----------//
644 
646 {
647  del_deriv() ;
648 }
649 
650 // -----------------//
651 // Memory management//
652 //------------------//
653 void Spheroid::del_deriv() const {
654  if (p_sqrt_q != 0x0) delete p_sqrt_q ;
655  if (p_area != 0x0) delete p_area ;
656  if (p_angu_mom != 0x0) delete p_angu_mom ;
657  if (p_mass != 0x0) delete p_mass ;
658  if (p_multipole_mass != 0x0) delete p_multipole_mass ;
659  if (p_multipole_angu != 0x0) delete p_multipole_angu ;
660  if (p_epsilon_A_minus_one != 0x0) delete p_epsilon_A_minus_one ;
661  if (p_epsilon_P_minus_one != 0x0) delete p_epsilon_P_minus_one ;
662  if (p_theta_plus != 0x0) delete p_theta_plus ;
663  if (p_theta_minus != 0x0) delete p_theta_minus ;
664  if (p_shear != 0x0) delete p_shear ;
665  if (p_delta != 0x0) delete p_delta ;
666  set_der_0x0() ;
667 }
668 
669 void Spheroid::set_der_0x0() const {
670  p_sqrt_q = 0x0 ;
671  p_area = 0x0 ;
672  p_angu_mom = 0x0 ;
673  p_mass = 0x0 ;
674  p_multipole_mass = 0x0;
675  p_multipole_angu = 0x0;
676  p_epsilon_A_minus_one = 0x0;
677  p_epsilon_P_minus_one = 0x0;
678  p_theta_plus = 0x0 ;
679  p_theta_minus = 0x0 ;
680  p_shear = 0x0 ;
681  p_delta = 0x0;
682 
683 }
684 
685 
686 
687 //---------//
688 //Accessors//
689 //---------//
690 
691 
692 
693 
694 
695 // Computation of the 2-dimensional Jacobian amplitude for the surface
696 const Scalar& Spheroid::sqrt_q() const {
697  if (p_sqrt_q == 0x0) {
698  p_sqrt_q = new Scalar(sqrt((get_qq()(2,2)*get_qq()(3,3))- (get_qq()(2,3)*get_qq()(2,3)))) ;
700  }
701  return *p_sqrt_q ;
702 }
703 
704 
705 
706 
707 
708 // Computation of the 2-dimensional area of the surface
709 
710 
711 double Spheroid::area() const {
712  if (p_area == 0x0) {
713  const Map_af& mp_ang = dynamic_cast<const Map_af&>(h_surf.get_mp()) ;
714  p_area = new double (mp_ang.integrale_surface((sqrt_q()) * h_surf *h_surf, 1)) ;
715  }
716  return *p_area;
717 
718 }
719 
720 
721 
722 // Computation of the angular momentum of the surface (G is set to be identically one)
723 
724 double Spheroid::angu_mom() const {
725  if (p_angu_mom == 0x0) {
726  const Map_af& mp_ang = dynamic_cast<const Map_af&>(h_surf.get_mp()) ;
727  Vector phi(mp_ang, CON, mp_ang.get_bvect_spher());
728  phi = get_ephi();
729  Scalar tmp = contract(ll,0, contract (jac2d, 1,phi,0), 0 );
730  p_angu_mom = new double (mp_ang.integrale_surface((sqrt_q()*h_surf*h_surf*tmp),1)) ;
731  *p_angu_mom = *p_angu_mom /(8. * M_PI) ;
732  }
733 
734  return *p_angu_mom;
735 
736 }
737 
738 
739 double Spheroid::mass() const {
740  if (p_mass == 0x0) {
741  double rayon = sqrt(area()/(4.*M_PI));
742  p_mass = new double ((1/(2.*rayon))*sqrt(rayon*rayon*rayon*rayon + 4.*angu_mom()*angu_mom()));
743 
744  }
745  return *p_mass;
746 
747 }
748 
749 
750 
751 double Spheroid::multipole_mass(const int order) const{
752  const Map_af& mp_ang = dynamic_cast<const Map_af&>(h_surf.get_mp()) ;
753  double rayon = sqrt(area()/(4.*M_PI));
754  // Multiplicative factor before integral.
755  double factor = mass()/(8. * M_PI); // To check later
756  if (order >0)
757  { for (int compte=0; compte <=order -1; compte++)
758  factor = factor*rayon;
759  }
760  // Calculus of legendre polynomial of order n, as function of cos(theta)
761  Scalar Pn(mp_ang); Pn=1; Pn.std_spectral_base(); Pn.set_spectral_va().ylm();
762  if (order >0)
763  { Pn = Pn*zeta;
764  }
765  if (order >1)
766  { Scalar Pnold(mp_ang); Pnold = 1; Pnold.std_spectral_base(); Pnold.set_spectral_va().ylm();
767 
768  for (int nn=1; nn<order; nn++){
769 
770  Scalar Pnnew = (2.*nn +1.)*Pn;
771  Pnnew = Pnnew*zeta;
772  Pnnew = Pnnew - nn*Pnold;
773  Pnnew = Pnnew/(double(nn) + 1.);
774 
775  Pnold = Pn;
776  Pn = Pnnew;
777 
778  }
779  }
780 
781  // Calculus of Ricci Scalar over the surface
782  Scalar ricciscal(mp_ang);
783  ricciscal = get_ricci();
784  ricciscal.set_spectral_va().ylm();
785 
786  Scalar rayyon = h_surf;
787  rayyon.std_spectral_base();
788  rayyon.set_spectral_va().ylm();
789 
790  Scalar sqq = sqrt_q();
791  Scalar integrande = sqq * rayyon *rayyon*ricciscal*Pn; integrande.std_spectral_base();
792 
793 
794  p_multipole_mass = new double (factor*mp_ang.integrale_surface(integrande, 1)) ;
795 
796 
797 
798 
799  return *p_multipole_mass;
800 }
801 
802 
803 
804 double Spheroid::multipole_angu(const int order) const{
805 
806  assert (order >=1) ;
807  const Map_af& mp_ang = dynamic_cast<const Map_af&>(h_surf.get_mp()) ;
808  Vector phi(mp_ang, CON, mp_ang.get_bvect_spher());
809  phi = get_ephi();
810  double rayon = sqrt(area()/(4.*M_PI));
811 
812  double factor = 1./(8. * M_PI);
813  if (order >1)
814  { for (int compte=0; compte <=order -2; compte++)
815  factor = factor*rayon;
816  }
817 
818  // Calculus of legendre polynomial of order n, as function of cos(theta)
819  Scalar Pn(mp_ang); Pn=1; Pn.std_spectral_base(); Pn.set_spectral_va().ylm();
820  Scalar dPn = Pn;
821 
822  Pn = Pn*zeta;
823 
824  if (order >1)
825 
826  { Scalar Pnold(mp_ang); Pnold = 1; Pnold.std_spectral_base(); Pnold.set_spectral_va().ylm();
827 
828  for (int nn=1; nn<order; nn++){
829 
830  Scalar Pnnew = (2.*nn +1.)*Pn;
831  Pnnew = Pnnew*zeta;
832  Pnnew = Pnnew - nn*Pnold;
833  Pnnew = Pnnew/(double(nn) + 1.);
834 
835  Pnold = Pn;
836  Pn = Pnnew; // Pn is now P(n+1)
837 
838  }
839 
840  // Calculus of functional derivative of order N legendre polynomial.
841 
842  dPn = Pn* zeta; dPn = dPn - Pnold; dPn = double(order)*dPn;
843 
844  Scalar quotient(mp_ang); quotient = 1.; quotient.std_spectral_base();
845  quotient = quotient*zeta*zeta; quotient = quotient -1.;
846 
847  dPn = dPn/quotient;
848 
849  }
850 
851  // Computation of the multipole;
852  Scalar tmp = contract(ll,0, contract (jac2d, 1,phi,0), 0 ); tmp.std_spectral_base();
853  Scalar tmp2 = (sqrt_q()) * h_surf *h_surf*tmp*dPn; tmp2.std_spectral_base();
854 
855 
856 
857  p_multipole_angu = new double (factor*mp_ang.integrale_surface(tmp2, 1)) ;
858 
859 
860  return *p_multipole_angu;
861 
862 }
863 
864 
865 // Computation of the refined Penrose parameter for axisymmetric spacetimes, and its difference wrt one.
866 
868  if (p_epsilon_A_minus_one == 0x0) {
869  assert (pow(mass(),4) - pow (angu_mom(),2) > 0.);
870  p_epsilon_A_minus_one = new double(area()/(8.*M_PI*(mass()*mass() + sqrt(pow(mass(),4) - pow (angu_mom(),2)))) - 1.);
871  }
872 
873  return *p_epsilon_A_minus_one;
874 
875 }
876 
877 // Computation of the classical Penrose parameter, and its difference wrt one.
878 // To use in replacement of epsilon_A_minus_one when the computed spacetime is not axisymmetric.
880  if (p_epsilon_P_minus_one == 0x0) {
881  assert (pow(mass(),4) - pow (angu_mom(),2) > 0.);
882  p_epsilon_P_minus_one = new double(area()/(16.*M_PI*mass()*mass()) - 1.);
883  }
884 
885  return *p_epsilon_P_minus_one;
886 
887 }
888 
889 
890 // Outgoing null expansion of 2-surface
891 
892 const Scalar &Spheroid::theta_plus() const {
893 
894  if (p_theta_plus == 0x0) {
895  p_theta_plus = new Scalar(fff*(hh.trace(qab) - jj.trace(qab))) ;
898  }
899 
900  return *p_theta_plus;
901 }
902 
903 
904 
905 
906 
907 
908 
909 
910 // ingoing null expansion of 2-surface
911 
913 
914  if (p_theta_minus == 0x0) {
915  p_theta_minus = new Scalar(ggg*(-hh.trace(qab) - jj.trace(qab))) ;
917  }
918 
919  return *p_theta_minus;
920 
921 }
922 
923 
924 
925 
926 //outer null-oriented shear of 2-surface
927 
928 const Sym_tensor& Spheroid::shear() const {
929 
930  if (p_shear == 0x0) {
931  p_shear = new Sym_tensor( fff*(hh - jj) - (qab.cov()/2) *(hh.trace(qab) - jj.trace(qab))) ;
932 // This is associated with the null vector: "l = n + s";
933 // For a null vector, "l = f (n + s)", multiply by the global factor "f" (e.g. the lapse "N" for "l=N(n+s)".
934 
935 
937  }
938 
939  return *p_shear;
940 
941 }
942 
943 
944 
945 
946 
947 
948 
949 
950 
951 
952 
953 //-------------------------------------------//
954 // Covariant flat derivative, returning a pointer.//
955 //-------------------------------------------//
956 
958 
959  // Notations: suffix 0 in name <=> input tensor
960  // suffix 1 in name <=> output tensor
961 
962  int valence0 = uu.get_valence() ;
963  int valence1 = valence0 + 1 ;
964  int valence1m1 = valence1 - 1 ; // same as valence0, but introduced for
965  // the sake of clarity
966 
967 
968  // Protections
969  // -----------
970  if (valence0 >= 1) {
971 
972  }
973 
974  // Creation of the result (pointer)
975  // --------------------------------
976  Tensor *resu ;
977 
978  // If uu is a Scalar, the result is a vector
979  if (valence0 == 0) {
980  resu = new Vector(uu.get_mp(), COV, uu.get_mp().get_bvect_spher()) ;
981  }
982  else {
983 
984  // Type of indices of the result :
985  Itbl tipe(valence1) ;
986  const Itbl& tipeuu = uu.get_index_type() ;
987  for (int id = 0; id<valence0; id++) {
988  tipe.set(id) = tipeuu(id) ; // First indices = same as uu
989  }
990  tipe.set(valence1m1) = COV ; // last index is the derivation index
991 
992  // if uu is a Tensor_sym, the result is also a Tensor_sym:
993  const Tensor* puu = &uu ;
994  const Tensor_sym* puus = dynamic_cast<const Tensor_sym*>(puu) ;
995  if ( puus != 0x0 ) { // the input tensor is symmetric
996  resu = new Tensor_sym(uu.get_mp(), valence1, tipe, *uu.get_triad(),
997  puus->sym_index1(), puus->sym_index2()) ;
998  }
999  else {
1000  resu = new Tensor(uu.get_mp(), valence1, tipe, *uu.get_triad()) ; // no symmetry
1001  }
1002 
1003  }
1004 
1005  int ncomp1 = resu->get_n_comp() ;
1006 
1007  Itbl ind1(valence1) ; // working Itbl to store the indices of resu
1008  Itbl ind0(valence0) ; // working Itbl to store the indices of uu
1009  Itbl ind(valence0) ; // working Itbl to store the indices of uu
1010 
1011  Scalar tmp(uu.get_mp()) ; // working scalar
1012 
1013  // Determination of the dzpuis parameter of the result --> dz_resu
1014  // --------------------------------------------------
1015 
1016  int dz_resu = 0;
1017 
1018  // (We only work here on a non-compactified shell) //
1019 
1020 
1021 
1022 
1023  // Loop on all the components of the output tensor
1024  // -----------------------------------------------
1025  /* Note: we have here preserved all the non-useful terms in this case(typically christoffel symbols) for the sake of understandng what's going on...
1026  */
1027 
1028 
1029  for (int ic=0; ic<ncomp1; ic++) {
1030 
1031  // indices corresponding to the component no. ic in the output tensor
1032  ind1 = resu->indices(ic) ;
1033 
1034  // Component no. ic:
1035  Scalar& cresu = resu->set(ind1) ;
1036 
1037  // Indices of the input tensor
1038  for (int id = 0; id < valence0; id++) {
1039  ind0.set(id) = ind1(id) ;
1040  }
1041 
1042  // Value of last index (derivation index)
1043  int k = ind1(valence1m1) ;
1044 
1045  switch (k) {
1046 
1047  case 1 : { // Derivation index = r
1048  //---------------------
1049 
1050  cresu = 0; //(uu(ind0)).dsdr() ; // d/dr
1051 
1052  // all the connection coefficients Gamma^i_{jk} are zero for k=1
1053  break ;
1054  }
1055 
1056  case 2 : { // Derivation index = theta
1057  //-------------------------
1058 
1059  cresu = (uu(ind0)).srdsdt() ; // 1/r d/dtheta
1060 
1061  // Loop on all the indices of uu
1062  for (int id=0; id<valence0; id++) {
1063 
1064  switch ( ind0(id) ) {
1065 
1066  case 1 : { // Gamma^r_{l theta} V^l
1067  // or -Gamma^l_{r theta} V_l
1068  /* ind = ind0 ;
1069  ind.set(id) = 2 ; // l = theta
1070 
1071  // Division by r :
1072  tmp = uu(ind) ;
1073  tmp.div_r_dzpuis(dz_resu) ;
1074 
1075  cresu -= tmp ; */
1076  break ;
1077  }
1078 
1079  case 2 : { // Gamma^theta_{l theta} V^l
1080  // or -Gamma^l_{theta theta} V_l
1081  /* ind = ind0 ;
1082  ind.set(id) = 1 ; // l = r
1083  tmp = uu(ind) ;
1084  tmp.div_r_dzpuis(dz_resu) ;
1085 
1086  cresu += tmp ; */
1087  break ;
1088  }
1089 
1090  case 3 : { // Gamma^phi_{l theta} V^l
1091  // or -Gamma^l_{phi theta} V_l
1092  break ;
1093  }
1094 
1095  default : {
1096  cerr << "Connection_fspher::p_derive_cov : index problem ! "
1097  << endl ;
1098  abort() ;
1099  }
1100  }
1101 
1102  }
1103  break ;
1104  }
1105 
1106 
1107  case 3 : { // Derivation index = phi
1108  //-----------------------
1109 
1110  cresu = (uu(ind0)).srstdsdp() ; // 1/(r sin(theta)) d/dphi
1111 
1112  // Loop on all the indices of uu
1113  for (int id=0; id<valence0; id++) {
1114 
1115  switch ( ind0(id) ) {
1116 
1117  case 1 : { // Gamma^r_{l phi} V^l
1118  // or -Gamma^l_{r phi} V_l
1119  /* ind = ind0 ;
1120  ind.set(id) = 3 ; // l = phi
1121  tmp = uu(ind) ;
1122  tmp.div_r_dzpuis(dz_resu) ;
1123 
1124  cresu -= tmp ; */
1125  break ;
1126  }
1127 
1128  case 2 : { // Gamma^theta_{l phi} V^l
1129  // or -Gamma^l_{theta phi} V_l
1130  ind = ind0 ;
1131  ind.set(id) = 3 ; // l = phi
1132  tmp = uu(ind) ;
1133  tmp.div_r_dzpuis(dz_resu) ;
1134 
1135  tmp.div_tant() ; // division by tan(theta)
1136 
1137  cresu -= tmp ;
1138  break ;
1139  }
1140 
1141  case 3 : { // Gamma^phi_{l phi} V^l
1142  // or -Gamma^l_{phi phi} V_l
1143 
1144  ind = ind0 ;
1145  // ind.set(id) = 1 ; // l = r
1146  // tmp = uu(ind) ;
1147  // tmp.div_r_dzpuis(dz_resu) ;
1148 
1149  // cresu += tmp ;
1150 
1151  ind.set(id) = 2 ; // l = theta
1152  tmp = uu(ind) ;
1153  tmp.div_r_dzpuis(dz_resu) ;
1154 
1155  tmp.div_tant() ; // division by tan(theta)
1156 
1157  cresu += tmp ;
1158  break ;
1159  }
1160 
1161  default : {
1162  cerr << "Connection_fspher::p_derive_cov : index problem ! \n"
1163  << endl ;
1164  abort() ;
1165  }
1166  }
1167 
1168  }
1169 
1170  break ;
1171  }
1172 
1173  default : {
1174  cerr << "Connection_fspher::p_derive_cov : index problem ! \n" ;
1175  abort() ;
1176  }
1177 
1178  } // End of switch on the derivation index
1179 
1180 
1181  } // End of loop on all the components of the output tensor
1182 
1183  // C'est fini !
1184  // -----------
1185  return *resu ;
1186 
1187 }
1188 
1189 void Spheroid::sauve(FILE* ) const {
1190 
1191  cout << "c'est pas fait!" << endl ;
1192  return ;
1193 
1194 }
1195 
1196 
1197 
1198 
1199 
1200 
1201 
1202 
1203 // Computation of the delta coefficients
1204 
1205 
1206 const Tensor& Spheroid::delta() const {
1207 
1208  if (p_delta == 0x0) {
1209 
1210  Tensor christoflat(qab.get_mp(), 3, COV, qab.get_mp().get_bvect_spher());
1211  christoflat.set_index_type(0) = CON;
1212  christoflat.set_etat_zero() ;
1213 
1214  // assert(flat_met != 0x0) ;
1215  Tensor dgam = derive_cov2dflat(qab.cov()) ;
1216 
1217  for (int k=1; k<=3; k++) {
1218  for (int i=1; i<=3; i++) {
1219  for (int j=1; j<=i; j++) {
1220  Scalar& cc= christoflat.set(k,i,j);
1221  for (int l=1; l<=3; l++) {
1222  cc += qab.con()(k,l) * (
1223  dgam(l,j,i) + dgam(i,l,j) - dgam(i,j,l) ) ;
1224 
1225  }
1226  cc = 0.5 * cc ;
1227  }
1228  }
1229  }
1230 
1231  p_delta = new Tensor (christoflat) ;
1232 
1233  }
1234  return *p_delta;
1235 }
1236 
1237 
1238 
1239 
1240 
1241 
1242 // Computation of global derivative on 2-sphere
1244 
1245  if(uu.get_valence()>=1){
1246  int nbboucle = uu.get_valence();
1247  Tensor resu = derive_cov2dflat(uu);
1248  for (int y=1; y<=nbboucle; y++){
1249 
1250  int df = uu.get_index_type(y-1);
1251  if (df == COV) {
1252  resu -= contract(delta(),0, uu, y-1);
1253  }
1254  else {resu += contract(delta(),1, uu, y-1);}
1255 
1256  return resu;
1257  }
1258  }
1259  else return derive_cov2dflat(uu);
1260 
1261  return derive_cov2dflat(uu); // to avoid warnings...
1262 }
1263 
1264 
1265 
1266 
1267 
1268 
1269 
1270 
1271 
1272 // // COmputation of the ricci tensor
1273 
1274 // const Sym_tensor& Spheroid::ricci() const {
1275 
1276 // if (p_ricci == 0x0) { // a new computation is necessary
1277 
1278 // assert( issphere == true ) ;
1279 // Sym_tensor riccia(h_surf.get_mp(), CON, h_surf.get_mp().get_bvect_spher()) ;
1280 // riccia.set_etat_zero();
1281 
1282 // const Tensor& d_delta = derive_cov2dflat(delta()) ;
1283 
1284 // for (int i=1; i<=3; i++) {
1285 
1286 // int jmax = 3 ;
1287 
1288 // for (int j=1; j<=jmax; j++) {
1289 
1290 // Scalar tmp1(h_surf.get_mp()) ;
1291 // tmp1.set_etat_zero() ;
1292 // for (int k=1; k<=3; k++) {
1293 // tmp1 += d_delta(k,i,j,k) ;
1294 // }
1295 
1296 // Scalar tmp2(h_surf.get_mp()) ;
1297 // tmp2.set_etat_zero() ;
1298 // for (int k=1; k<=3; k++) {
1299 // tmp2 += d_delta(k,i,k,j) ;
1300 // }
1301 
1302 // Scalar tmp3(h_surf.get_mp()) ;
1303 // tmp3.set_etat_zero() ;
1304 // for (int k=1; k<=3; k++) {
1305 // for (int m=1; m<=3; m++) {
1306 // tmp3 += delta()(k,k,m) * delta()(m,i,j) ;
1307 // }
1308 // }
1309 // tmp3.dec_dzpuis() ; // dzpuis 4 -> 3
1310 
1311 // Scalar tmp4(h_surf.get_mp()) ;
1312 // tmp4.set_etat_zero() ;
1313 // for (int k=1; k<=3; k++) {
1314 // for (int m=1; m<=3; m++) {
1315 // tmp4 += delta()(k,j,m) * delta()(m,i,k) ;
1316 // }
1317 // }
1318 // tmp4.dec_dzpuis() ; // dzpuis 4 -> 3
1319 
1320 
1321 // riccia.set(i,j) = tmp1 - tmp2 + tmp3 - tmp4 ;
1322 
1323 
1324 // }
1325 // }
1326 // /* Note: Here we must take into account the fact that a round metric on a spheroid doesn't give zero as "flat" ricci part. Then a diagonal scalar term must be added.
1327 // WARNING: this only works with "round" horizons!! */
1328 
1329 // double rayon = sqrt(area()/(4.*M_PI));
1330 // Scalar rayon2 = h_surf;
1331 // rayon2 = rayon;
1332 // rayon2.std_spectral_base();
1333 
1334 // for (int hi=1; hi<=3; hi++){
1335 
1336 // riccia.set(hi,hi) += 2/(rayon2 * rayon2) ; // Plutot 1/hsurf^2, non?
1337 // }
1338 // p_ricci = new Sym_tensor(riccia);
1339 // }
1340 
1341 // return *p_ricci ;
1342 
1343 // }
1344 
1345 
1346 
1347 
1348 // COmputation of the ricci tensor
1349 
1350 // const Sym_tensor& Spheroid::ricci() const {
1351 
1352 // if (p_ricci == 0x0) { // a new computation is necessary
1353 // Sym_tensor riccia(h_surf.get_mp(), CON, h_surf.get_mp().get_bvect_spher()) ;
1354 // Sym_tensor ricci3 = gamij.ricci();
1355 
1356 // Sym_tensor ricci3-2d(h_surf.get_mp(), COV, h_surf.get_mp().get_bvect_spher());
1357 // ricci3-2d.allocate_all() ;
1358 // ricci3-2d.std_spectral_base();
1359 // for (int l=1; l<4; l++)
1360 // for (int m=1; m<4; m++)
1361 // for (int k=0; k<np; k++)
1362 // for (int j=0; j<nt; j++)
1363 // {
1364 // mp_rad->val_lx_jk(h_surf.val_grid_point(0, k, j, 0)*1.000000000001, j, k, pipo,
1365 // lz, xi) ;
1366 // ricci3-2d.set(l,m).set_grid_point(0, k, j, 0) =
1367 // ricci3(l,m).get_spectral_va().val_point_jk(lz, xi, j, k) ;
1368 
1369 // }
1370 
1371 
1372 
1373 // }
1374 // return *p_ricci ;
1375 
1376 // }
1377 
1378 
1379 
1380 
1381 
1382 }
Sym_tensor * p_shear
The shear tensor.
Definition: spheroid.h:167
Scalar h_surf
The location of the 2-surface as r = h_surf .
Definition: spheroid.h:91
const Base_val & get_spectral_base() const
Returns the spectral bases of the Valeur va.
Definition: scalar.h:1328
void annule_domain(int l)
Sets the Tensor to zero in a given domain.
Definition: tensor.C:675
const Scalar & theta_plus() const
Computes the outgoing null expansion .
Definition: spheroid.C:892
Metric for tensor calculation.
Definition: metric.h:90
double * p_epsilon_P_minus_one
Refined Penrose parameter, difference wrt one.
Definition: spheroid.h:164
double multipole_mass(const int order) const
Computes the mass multipole of a given order for the spheroid, assumed to be spherical.
Definition: spheroid.C:751
double * p_area
The area of the 2-surface.
Definition: spheroid.h:158
double angu_mom() const
Computes the angular momentum with respect to a divergence-free vector field tangent to the 2-surface...
Definition: spheroid.C:724
virtual const Sym_tensor & con() const
Read-only access to the contravariant representation.
Definition: metric.C:293
Tensor jac2d
The jacobian of the adaptation of coordinates (contravariant/covariant representation) ...
Definition: spheroid.h:96
const Scalar & sqrt_q() const
Computes the normal vector field to the 2-surface.
Definition: spheroid.C:696
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
Definition: valeur.h:312
int & set(int i)
Read/write of a particular element (index i ) (1D case)
Definition: itbl.h:247
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.
Definition: spheroid.C:653
double * p_angu_mom
The angular momentum.
Definition: spheroid.h:159
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition: grilles.h:479
void coef() const
Computes the coeffcients of *this.
Definition: valeur_coef.C:151
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:223
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: scalar.C:330
double multipole_angu(const int order) const
Computes the angular multipole of a given order for the spheroid, assumed to be spherical.
Definition: spheroid.C:804
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 )
Definition: tensor.h:1162
const Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition: map.h:801
Lorene prototypes.
Definition: app_hor.h:67
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.
Definition: valeur_ylm.C:141
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition: map.h:783
Tensor field of valence 0 (or component of a tensorial field).
Definition: scalar.h:393
Base class for coordinate mappings.
Definition: map.h:688
Vector ephi
The conformal Killing vector field on the 2-surface (set to by default to the axial vector associated...
Definition: spheroid.h:113
int get_n_comp() const
Returns the number of stored components.
Definition: tensor.h:885
int sym_index2() const
Number of the second symmetric index (id_sym1 < id_sym2 < valence )
Definition: tensor.h:1167
Basic integer array class.
Definition: itbl.h:122
const Sym_tensor & shear() const
Computes the shear of the 2-surface .
Definition: spheroid.C:928
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field. ...
Definition: scalar.C:790
const Sym_tensor & ricci() const
Returns the Ricci tensor (given by the Connection p_connect )
Definition: metric.C:341
virtual void allocate_all()
Sets the logical state to ETATQCQ (ordinary state) and performs the memory allocation of all the elem...
Definition: scalar.C:371
const Vector & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of *this with respect to some metric , by raising the index of...
Definition: scalar_deriv.C:402
Tensor field of valence 1.
Definition: vector.h:188
Scalar * p_theta_minus
Null ingoing expansion.
Definition: spheroid.h:166
const Map & get_mp() const
Returns the mapping.
Definition: metric.h:202
double val_grid_point(int l, int k, int j, int i) const
Returns the value of the field at a specified grid point.
Definition: scalar.h:643
virtual void std_spectral_base()
Sets the standard spectal bases of decomposition for each component.
Definition: vector.C:322
void annule_hard()
Sets the Scalar to zero in a hard way.
Definition: scalar.C:386
Spheroidal 2-surfaces embedded in a time-slice of the 3+1 formalism.
Definition: spheroid.h:84
Spheroid(const Map_af &map, double radius)
The delta tensorial fields linked to Christoffel symbols.
Definition: spheroid.C:97
const Vector & get_ephi() const
Returns the conformal Killing symmetry vector on the 2-surface.
Definition: spheroid.h:253
Tensor derive_cov2dflat(const Tensor &uu) const
Computes the round covariant derivative on the spheroid.
Definition: spheroid.C:957
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition: tensor.h:879
Tensor derive_cov2d(const Tensor &uu) const
Computes the total covariant derivative on the spheroid.
Definition: spheroid.C:1243
int get_index_type(int i) const
Gives the type (covariant or contravariant) of the index number i .
Definition: tensor.h:899
virtual void allocate_all()
Performs the memory allocation of all the elements, down to the double arrays of the Tbl s...
Definition: tensor.C:517
bool issphere
Flag to know whether the horizon is geometrically round or distorted.
Definition: spheroid.h:151
Scalar trk
Trace K of the extrinsic curvature of the 3-slice.
Definition: spheroid.h:124
Parameter storage.
Definition: param.h:125
Base class for pure radial mappings.
Definition: map.h:1557
Vector ll
Normal-tangent component of the extrinsic curvature of the 3-slice.
Definition: spheroid.h:129
double * p_multipole_angu
Angular momentum multipole for the spheroid.
Definition: spheroid.h:162
void set_der_0x0() const
Sets to 0x0 all the pointers on derived quantities.
Definition: spheroid.C:669
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:465
void operator=(const Spheroid &)
Assignment to another Spheroid.
Definition: spheroid.C:617
virtual void sauve(FILE *) const
Save in a file.
Definition: spheroid.C:1189
double epsilon_P_minus_one() const
Computation of the classical Penrose parameter, and its difference wrt one.
Definition: spheroid.C:879
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:903
void mult_cost()
Multiplication by .
const Scalar & get_ricci() const
Returns the 2-ricci scalar .
Definition: spheroid.h:229
double area() const
Computes the area of the 2-surface.
Definition: spheroid.C:711
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
Tensor handling.
Definition: tensor.h:294
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.
Definition: spheroid.h:134
int get_valence() const
Returns the valence.
Definition: tensor.h:882
int & set_index_type(int i)
Sets the type of the index number i .
Definition: tensor.h:922
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.
Definition: grilles.h:469
virtual const Sym_tensor & cov() const
Read-only access to the covariant representation.
Definition: metric.C:283
Scalar ggg
Normalization function for the ingoing null vector k.
Definition: spheroid.h:143
Multi-domain grid.
Definition: grilles.h:279
Bases of the spectral expansions.
Definition: base_val.h:325
double & set_grid_point(int l, int k, int j, int i)
Setting the value of the field at a given grid point.
Definition: scalar.h:690
Affine radial mapping.
Definition: map.h:2048
const Scalar & srdsdt() const
Returns of *this .
Definition: scalar_deriv.C:145
const Scalar & theta_minus() const
Computes the ingoing null expansion .
Definition: spheroid.C:912
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition: tensor.C:663
Sym_tensor hh
The ricci scalar on the 2-surface.
Definition: spheroid.h:122
const Sym_tensor & get_qq() const
returns the 3-d degenerate 2-metric
Definition: spheroid.h:235
Scalar fff
Normalization function for the outgoing null vector l.
Definition: spheroid.h:138
Coefficients storage for the multi-domain spectral method.
Definition: mtbl_cf.h:196
double * p_mass
Mass defined from angular momentum.
Definition: spheroid.h:160
virtual void set_etat_zero()
Sets the logical state of all components to ETATZERO (zero state).
Definition: tensor.C:506
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition: grilles.h:474
Valeur & set_spectral_va()
Returns va (read/write version)
Definition: scalar.h:610
Scalar & set(int)
Read/write access to a component.
Definition: vector.C:302
Symmetric tensors (with respect to two of their arguments).
Definition: tensor.h:1050
Vector ss
The adapted normal vector field to spheroid in the 3-slice.
Definition: spheroid.h:108
double * p_multipole_mass
Mass multipole for the spheroid.
Definition: spheroid.h:161
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
Definition: grilles.h:491
Tensor proj
The 3-d projector on the 2-surface (contravariant-covariant form).
Definition: spheroid.h:100
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this
Definition: scalar_deriv.C:390
const Map & get_mp() const
Returns the mapping.
Definition: tensor.h:874
const Mg3d * get_angu_1dom() const
Returns the pointer on the associated mono-domain angular grid.
Definition: mg3d.C:625
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.
Definition: tensor.C:935
Scalar * p_theta_plus
Classical Penrose parameter, difference wrt one.
Definition: spheroid.h:165
Scalar ricci
Induced metric on the 2-surface .
Definition: spheroid.h:117
Class intended to describe valence-2 symmetric tensors.
Definition: sym_tensor.h:226
double epsilon_A_minus_one() const
Computation of the refined Penrose parameter for axisymmetric spacetimes, and its difference wrt one...
Definition: spheroid.C:867
Scalar * p_sqrt_q
Surface element .
Definition: spheroid.h:157
const Valeur & get_spectral_va() const
Returns va (read only version)
Definition: scalar.h:607
Sym_tensor qq
The 3-d covariant degenerated 2-metric on the surface.
Definition: spheroid.h:104
virtual Itbl indices(int pos) const
Returns the indices of a component given by its position in the array cmp .
Definition: tensor.C:548
const Scalar & srstdsdp() const
Returns of *this .
Definition: scalar_deriv.C:177
double mass() const
Computes the mass as defined from the calculus of angular momentum, done with respect to a divergence...
Definition: spheroid.C:739
Tensor trace(int ind1, int ind2) const
Trace on two different type indices.
Coord r
r coordinate centered on the grid
Definition: map.h:736
virtual ~Spheroid()
Destructor.
Definition: spheroid.C:645
const Tensor & delta() const
Computes the delta coefficients for covariant derivative.
Definition: spheroid.C:1206
Tensor down(int ind, const Metric &gam) const
Computes a new tensor by lowering an index of *this.