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
26char 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
93namespace Lorene {
94Spheroid::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
130
131
132 h_surf = radius ;
136 hh.set_etat_zero() ;
137 for (int i=1; i<=3; i++)
138 hh.set(i,i) = 2./radius ;
140 ll.set_etat_zero() ;
141 jj.set_etat_zero() ;
144 zeta.set_etat_zero();
145 set_der_0x0() ;
146
147
148}
149
150Spheroid::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
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;
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
591Spheroid::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
614void 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//------------------//
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 ;
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
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;
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
693const 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
708double 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
721double 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
736double 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
748double 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
801double 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
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
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
1186void 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
1203const 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 Base_vect_spher & get_bvect_spher() const
Returns the orthonormal vectorial basis associated with the coordinates of the mapping.
Definition map.h:783
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
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
const Map & get_mp() const
Returns the mapping.
Definition metric.h:202
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
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 Scalar & srdsdt() const
Returns of *this .
const Vector & derive_cov(const Metric &gam) const
Returns the gradient (1-form = covariant vector) of *this
void div_r_dzpuis(int ced_mult_r)
Division by r everywhere but with the output flag dzpuis set to ced_mult_r .
void mult_cost()
Multiplication by .
const Scalar & srstdsdp() const
Returns of *this .
virtual void std_spectral_base()
Sets the spectral bases of the Valeur va to the standard ones for a scalar field.
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
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition scalar.C:324
void div_tant()
Division by .
Valeur & set_spectral_va()
Returns va (read/write version)
Definition scalar.h:604
const Valeur & get_spectral_va() const
Returns va (read only version)
Definition scalar.h:601
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 Base_val & get_spectral_base() const
Returns the spectral bases of the Valeur va
Definition scalar.h:1294
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
const Vector & derive_con(const Metric &gam) const
Returns the "contravariant" derivative of *this with respect to some metric , by raising the index of...
virtual void dec_dzpuis(int dec=1)
Decreases by dec units the value of dzpuis and changes accordingly the values of the Scalar in the co...
void mult_r()
Multiplication by r everywhere; dzpuis is not changed.
Spheroidal 2-surfaces embedded in a time-slice of the 3+1 formalism.
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 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 Scalar & get_ricci() const
Returns the 2-ricci scalar .
Definition spheroid.h:229
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
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
const Sym_tensor & get_qq() const
returns the 3-d degenerate 2-metric
Definition spheroid.h:235
Vector ephi
The conformal Killing vector field on the 2-surface (set to by default to the axial vector associated...
Definition spheroid.h:113
const Vector & get_ephi() const
Returns the conformal Killing symmetry vector on the 2-surface.
Definition spheroid.h:253
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
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.
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 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 & set_index_type(int i)
Sets the type of the index number i .
Definition tensor.h:909
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.
const Base_vect * get_triad() const
Returns the vectorial basis (triad) on which the components are defined.
Definition tensor.h:866
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
Tenseur contract(const Tenseur &, int id1, int id2)
Self contraction of two indices of a Tenseur .
Lorene prototypes.
Definition app_hor.h:64