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