LORENE
sym_tensor_aux.C
1/*
2 * Methods of class Sym_tensor linked with auxiliary members (eta, mu, W, X...)
3 *
4 * (see file sym_tensor.h for documentation)
5 *
6 */
7
8/*
9 * Copyright (c) 2005 Jerome Novak
10 *
11 * This file is part of LORENE.
12 *
13 * LORENE is free software; you can redistribute it and/or modify
14 * it under the terms of the GNU General Public License as published by
15 * the Free Software Foundation; either version 2 of the License, or
16 * (at your option) any later version.
17 *
18 * LORENE is distributed in the hope that it will be useful,
19 * but WITHOUT ANY WARRANTY; without even the implied warranty of
20 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 * GNU General Public License for more details.
22 *
23 * You should have received a copy of the GNU General Public License
24 * along with LORENE; if not, write to the Free Software
25 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
26 *
27 */
28
29
30char sym_tensor__aux_C[] = "$Header: /cvsroot/Lorene/C++/Source/Tensor/sym_tensor_aux.C,v 1.15 2014/10/13 08:53:43 j_novak Exp $" ;
31
32/*
33 * $Id: sym_tensor_aux.C,v 1.15 2014/10/13 08:53:43 j_novak Exp $
34 * $Log: sym_tensor_aux.C,v $
35 * Revision 1.15 2014/10/13 08:53:43 j_novak
36 * Lorene classes and functions now belong to the namespace Lorene.
37 *
38 * Revision 1.14 2014/10/06 15:13:19 j_novak
39 * Modified #include directives to use c++ syntax.
40 *
41 * Revision 1.13 2007/11/27 15:49:51 n_vasset
42 * new function compute_tilde_C in class sym_tensor
43 *
44 * Revision 1.12 2006/10/24 14:52:38 j_novak
45 * The Mtbl corresponding to the physical space is destroyed after the
46 * calculation of tilde_B(tt), to get the updated result.
47 *
48 * Revision 1.11 2006/10/24 13:03:19 j_novak
49 * New methods for the solution of the tensor wave equation. Perhaps, first
50 * operational version...
51 *
52 * Revision 1.10 2006/08/31 12:12:43 j_novak
53 * Correction of a small mistake in a phi loop.
54 *
55 * Revision 1.9 2006/06/28 07:48:26 j_novak
56 * Better treatment of some null cases.
57 *
58 * Revision 1.8 2006/06/12 11:40:07 j_novak
59 * Better memory and spectral base handling for Sym_tensor::compute_tilde_B.
60 *
61 * Revision 1.7 2006/06/12 07:42:29 j_novak
62 * Fields A and tilde{B} are defined only for l>1.
63 *
64 * Revision 1.6 2006/06/12 07:27:20 j_novak
65 * New members concerning A and tilde{B}, dealing with the transverse part of the
66 * Sym_tensor.
67 *
68 * Revision 1.5 2005/09/07 16:47:43 j_novak
69 * Removed method Sym_tensor_trans::T_from_det_one
70 * Modified Sym_tensor::set_auxiliary, so that it takes eta/r and mu/r as
71 * arguments.
72 * Modified Sym_tensor_trans::set_hrr_mu.
73 * Added new protected method Sym_tensor_trans::solve_hrr
74 *
75 * Revision 1.4 2005/04/05 15:38:08 j_novak
76 * Changed the formulas for W and X. There was an error before...
77 *
78 * Revision 1.3 2005/04/05 09:22:15 j_novak
79 * Use of the right formula with poisson_angu(2) for the determination of W and
80 * X.
81 *
82 * Revision 1.2 2005/04/04 15:25:24 j_novak
83 * Added new members www, xxx, ttt and the associated methods.
84 *
85 * Revision 1.1 2005/04/01 14:39:57 j_novak
86 * Methods dealing with auxilliary derived members of the symmetric
87 * tensor (eta, mu, W, X, etc ...).
88 *
89 *
90 *
91 * $Header: /cvsroot/Lorene/C++/Source/Tensor/sym_tensor_aux.C,v 1.15 2014/10/13 08:53:43 j_novak Exp $
92 *
93 */
94
95// Headers C
96#include <cstdlib>
97#include <cassert>
98#include <cmath>
99
100// Headers Lorene
101#include "metric.h"
102#include "proto.h"
103
104
105 //--------------//
106 // eta //
107 //--------------//
108
109
110namespace Lorene {
112
113 if (p_eta == 0x0) { // a new computation is necessary
114
115 // All this has a meaning only for spherical components:
116 assert(dynamic_cast<const Base_vect_spher*>(triad) != 0x0) ;
117
118 // eta is computed from its definition (148) of Bonazzola et al. (2004)
119 int dzp = operator()(1,1).get_dzpuis() ;
120 int dzp_resu = ((dzp == 0) ? 0 : dzp-1) ;
121
123 source_eta.div_tant() ;
124 source_eta += operator()(1,2).dsdt() + operator()(1,3).stdsdp() ;
125 source_eta.mult_r_dzpuis(dzp_resu) ;
126
127 // Resolution of the angular Poisson equation for eta
128 // --------------------------------------------------
129 if (dynamic_cast<const Map_af*>(mp) != 0x0) {
130 p_eta = new Scalar( source_eta.poisson_angu() ) ;
131 }
132 else {
133 Scalar resu (*mp) ;
134 resu = 0. ;
135 mp->poisson_angu(source_eta, *par, resu) ;
136 p_eta = new Scalar( resu ) ;
137 }
138
139 }
140
141 return *p_eta ;
142
143}
144
145
146 //--------------//
147 // mu //
148 //--------------//
149
150
152
153 if (p_mu == 0x0) { // a new computation is necessary
154
155 // All this has a meaning only for spherical components:
156 assert(dynamic_cast<const Base_vect_spher*>(triad) != 0x0) ;
157
158 Scalar source_mu = operator()(1,3) ; // h^{r ph}
159 source_mu.div_tant() ; // h^{r ph} / tan(th)
160
161 // dh^{r ph}/dth + h^{r ph}/tan(th) - 1/sin(th) dh^{r th}/dphi
162 source_mu += operator()(1,3).dsdt() - operator()(1,2).stdsdp() ;
163
164 // Multiplication by r
165 int dzp = operator()(1,2).get_dzpuis() ;
166 int dzp_resu = ((dzp == 0) ? 0 : dzp-1) ;
167 source_mu.mult_r_dzpuis(dzp_resu) ;
168
169 // Resolution of the angular Poisson equation for mu
170 // --------------------------------------------------
171 if (dynamic_cast<const Map_af*>(mp) != 0x0) {
172 p_mu = new Scalar( source_mu.poisson_angu() ) ;
173 }
174 else {
175 Scalar resu (*mp) ;
176 resu = 0. ;
177 mp->poisson_angu(source_mu, *par, resu) ;
178 p_mu = new Scalar( resu ) ;
179 }
180 }
181 return *p_mu ;
182
183}
184
185 //-------------//
186 // T //
187 //-------------//
188
189
190const Scalar& Sym_tensor::ttt() const {
191
192 if (p_ttt == 0x0) { // a new computation is necessary
193
194 // All this has a meaning only for spherical components:
195 assert(dynamic_cast<const Base_vect_spher*>(triad) != 0x0) ;
196
197 p_ttt = new Scalar( operator()(2,2) + operator()(3,3) ) ;
198
199 }
200 return *p_ttt ;
201
202}
203
204 //------------//
205 // W //
206 //------------//
207
208
209const Scalar& Sym_tensor::www() const {
210
211 if (p_www == 0x0) { // a new computation is necessary
212
213 // All this has a meaning only for spherical components:
214 assert(dynamic_cast<const Base_vect_spher*>(triad) != 0x0) ;
215
216 Scalar ppp = 0.5*(operator()(2,2) - operator()(3,3)) ;
217 Scalar tmp = ppp.dsdt() ;
218 tmp.div_tant() ;
219 Scalar source_w = ppp.dsdt().dsdt() +3*tmp - ppp.stdsdp().stdsdp() -2*ppp ;
220 tmp = operator()(2,3) ;
221 tmp.div_tant() ;
222 tmp += operator()(2,3).dsdt() ;
223 source_w += 2*tmp.stdsdp() ;
224
225 // Resolution of the angular Poisson equation for W
226 p_www = new Scalar( source_w.poisson_angu(2).poisson_angu() ) ;
227
228 }
229
230 return *p_www ;
231
232}
233
234
235 //------------//
236 // X //
237 //------------//
238
239
240const Scalar& Sym_tensor::xxx() const {
241
242 if (p_xxx == 0x0) { // a new computation is necessary
243
244 // All this has a meaning only for spherical components:
245 assert(dynamic_cast<const Base_vect_spher*>(triad) != 0x0) ;
246
247 const Scalar& htp = operator()(2,3) ;
248 Scalar tmp = htp.dsdt() ;
249 tmp.div_tant() ;
250 Scalar source_x = htp.dsdt().dsdt() + 3*tmp - htp.stdsdp().stdsdp() -2*htp ;
251 Scalar ppp = 0.5*(operator()(2,2) - operator()(3,3)) ;
252 tmp = ppp ;
253 tmp.div_tant() ;
254 tmp += ppp.dsdt() ;
255 source_x -= 2*tmp.stdsdp() ;
256
257 // Resolution of the angular Poisson equation for W
258 p_xxx = new Scalar( source_x.poisson_angu(2).poisson_angu() ) ;
259
260 }
261
262 return *p_xxx ;
263
264}
265
267 const Scalar& mu_over_r, const Scalar& w_in,
268 const Scalar& x_in, const Scalar& t_in ) {
269
270 // All this has a meaning only for spherical components:
271 assert(dynamic_cast<const Base_vect_spher*>(triad) != 0x0) ;
272 int dzp = trr.get_dzpuis() ;
273 int dzeta = (dzp == 0 ? 0 : dzp - 1) ;
274
275 assert(eta_over_r.check_dzpuis(dzp)) ;
276 assert(mu_over_r.check_dzpuis(dzp)) ;
277 assert(w_in.check_dzpuis(dzp)) ;
278 assert(x_in.check_dzpuis(dzp)) ;
279 assert(t_in.check_dzpuis(dzp)) ;
280
281 assert(&w_in != p_www) ;
282 assert(&x_in != p_xxx) ;
283 assert(&t_in != p_ttt) ;
284
285 set(1,1) = trr ;
286 set(1,2) = eta_over_r.dsdt() - mu_over_r.stdsdp() ;
287 // set(1,2).div_r_dzpuis(dzp) ;
288 set(1,3) = eta_over_r.stdsdp() + mu_over_r.dsdt() ;
289 // set(1,3).div_r_dzpuis(dzp) ;
290 Scalar tmp = w_in.lapang() ;
291 tmp.set_spectral_va().ylm_i() ;
292 Scalar ppp = 2*w_in.dsdt().dsdt() - tmp - 2*x_in.stdsdp().dsdt() ;
293 tmp = x_in.lapang() ;
294 tmp.set_spectral_va().ylm_i() ;
295 set(2,3) = 2*x_in.dsdt().dsdt() - tmp + 2*w_in.stdsdp().dsdt() ;
296 set(2,2) = 0.5*t_in + ppp ;
297 set(3,3) = 0.5*t_in - ppp ;
298
299 // Deleting old derived quantities ...
300 del_deriv() ;
301
302 // .. and affecting new ones.
303 p_eta = new Scalar(eta_over_r) ; p_eta->mult_r_dzpuis(dzeta) ;
304 p_mu = new Scalar(mu_over_r) ; p_mu->mult_r_dzpuis(dzeta) ;
305 p_www = new Scalar(w_in) ;
306 p_xxx = new Scalar(x_in) ;
307 p_ttt = new Scalar(t_in) ;
308
309}
310
311 //------------//
312 // A //
313 //------------//
314
315
317
318 if (p_aaa == 0x0) { // a new computation is necessary
319
320 // All this has a meaning only for spherical components:
321 assert(dynamic_cast<const Base_vect_spher*>(triad) != 0x0) ;
322
323 int dzp = xxx().get_dzpuis() ;
324 int dzp_resu = ((dzp == 0) ? 2 : dzp+1) ;
325
326 Scalar source_mu = operator()(1,3) ; // h^{r ph}
327 source_mu.div_tant() ; // h^{r ph} / tan(th)
328
329 // dh^{r ph}/dth + h^{r ph}/tan(th) - 1/sin(th) dh^{r th}/dphi
330 source_mu += operator()(1,3).dsdt() - operator()(1,2).stdsdp() ;
331
332 // Resolution of the angular Poisson equation for mu / r
333 // -----------------------------------------------------
334 Scalar tilde_mu(*mp) ;
335 tilde_mu = 0 ;
336
337 if (dynamic_cast<const Map_af*>(mp) != 0x0) {
338 tilde_mu = source_mu.poisson_angu() ;
339 }
340 else {
341 assert(par != 0x0) ;
342 mp->poisson_angu(source_mu, *par, tilde_mu) ;
343 }
344 tilde_mu.div_r_dzpuis(dzp_resu) ;
345
346 p_aaa = new Scalar( xxx().dsdr() - tilde_mu) ;
347 p_aaa->annule_l(0, 1, output_ylm) ;
348 }
349
351 else p_aaa->set_spectral_va().ylm_i() ;
352
353 return *p_aaa ;
354
355}
356
357 //--------------//
358 // tilde(B) //
359 //--------------//
360
361
363
364 if (p_tilde_b == 0x0) { // a new computation is necessary
365
366 // All this has a meaning only for spherical components:
367 assert(dynamic_cast<const Base_vect_spher*>(triad) != 0x0) ;
368
369 int dzp = operator()(1,1).get_dzpuis() ;
370 int dzp_resu = ((dzp == 0) ? 2 : dzp+1) ;
371
372 p_tilde_b = new Scalar(*mp) ;
374
376 source_eta.div_tant() ;
377 source_eta += operator()(1,2).dsdt() + operator()(1,3).stdsdp() ;
378
379 // Resolution of the angular Poisson equation for eta / r
380 // ------------------------------------------------------
382 tilde_eta = 0 ;
383
384 if (dynamic_cast<const Map_af*>(mp) != 0x0) {
385 tilde_eta = source_eta.poisson_angu() ;
386 }
387 else {
388 assert (par != 0x0) ;
389 mp->poisson_angu(source_eta, *par, tilde_eta) ;
390 }
391
392 Scalar dwdr = www().dsdr() ;
393 dwdr.set_spectral_va().ylm() ;
394 Scalar wsr = www() ;
395 wsr.div_r_dzpuis(dzp_resu) ;
396 wsr.set_spectral_va().ylm() ;
398 etasr2.div_r_dzpuis(dzp_resu) ;
399 etasr2.set_spectral_va().ylm() ;
400 Scalar dtdr = ttt().dsdr() ;
401 dtdr.set_spectral_va().ylm() ;
402 Scalar tsr = ttt() ;
403 tsr.div_r_dzpuis(dzp_resu) ;
404 tsr.set_spectral_va().ylm() ;
405 Scalar hrrsr = operator()(1,1) ;
406 hrrsr.div_r_dzpuis(dzp_resu) ;
407 hrrsr.set_spectral_va().ylm() ;
408
409 int nz = mp->get_mg()->get_nzone() ;
410
411 Base_val base(nz) ;
412
413 if (wsr.get_etat() != ETATZERO) {
414 base = wsr.get_spectral_base() ;
415 }
416 else {
417 if (etasr2.get_etat() != ETATZERO) {
418 base = etasr2.get_spectral_base() ;
419 }
420 else {
421 if (tsr.get_etat() != ETATZERO) {
422 base = tsr.get_spectral_base() ;
423 }
424 else {
425 if (hrrsr.get_etat() != ETATZERO) {
426 base = hrrsr.get_spectral_base() ;
427 }
428 else { //Everything is zero...
430 return *p_tilde_b ;
431 }
432 }
433 }
434 }
435
439
440 if (dwdr.get_etat() == ETATZERO) dwdr.annule_hard() ;
441 if (wsr.get_etat() == ETATZERO) wsr.annule_hard() ;
442 if (etasr2.get_etat() == ETATZERO) etasr2.annule_hard() ;
443 if (dtdr.get_etat() == ETATZERO) dtdr.annule_hard() ;
444 if (tsr.get_etat() == ETATZERO) tsr.annule_hard() ;
445 if (hrrsr.get_etat() == ETATZERO) hrrsr.annule_hard() ;
446
447 int m_q, l_q, base_r ;
448 for (int lz=0; lz<nz; lz++) {
449 int np = mp->get_mg()->get_np(lz) ;
450 int nt = mp->get_mg()->get_nt(lz) ;
451 int nr = mp->get_mg()->get_nr(lz) ;
452 for (int k=0; k<np+1; k++)
453 for (int j=0; j<nt; j++) {
454 base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
455 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q > 1))
456 {
457 for (int i=0; i<nr; i++)
459 = (l_q+2)*(*dwdr.get_spectral_va().c_cf)(lz, k, j, i)
460 + l_q*(l_q+2)*(*wsr.get_spectral_va().c_cf)(lz, k, j, i)
461 - 2*(*etasr2.get_spectral_va().c_cf)(lz, k, j, i)
462 + 0.5*double(l_q+2)/double(l_q+1)*(*tsr.get_spectral_va().c_cf)(lz, k, j, i)
463 + 0.5/double(l_q+1)*(*dtdr.get_spectral_va().c_cf)(lz, k, j, i)
464 - 1./double(l_q+1)*(*hrrsr.get_spectral_va().c_cf)(lz, k, j, i) ;
465 }
466 }
467 }
469 } //End of p_tilde_b == 0x0
470
473
474 return *p_tilde_b ;
475
476}
477
479
480 Scalar resu = compute_tilde_B(true, par) ;
481 if (resu.get_etat() == ETATZERO) return resu ;
482 else {
483 assert(resu.get_etat() == ETATQCQ) ;
484 assert(resu.get_spectral_va().c_cf != 0x0) ;
485 int dzp = operator()(1,1).get_dzpuis() ;
486 int dzp_resu = ((dzp == 0) ? 2 : dzp+1) ;
487
488 Scalar hsr = operator()(1,1) + ttt() ;
489 if (hsr.get_etat() == ETATZERO) return resu ;
490 Scalar dhdr = hsr.dsdr() ;
491 dhdr.set_spectral_va().ylm() ;
492 hsr.div_r_dzpuis(dzp_resu) ;
493 hsr.set_spectral_va().ylm() ;
494
495 int nz = mp->get_mg()->get_nzone() ;
496
497 const Base_val& base = resu.get_spectral_base() ;
498
499 if (dhdr.get_etat() == ETATZERO) dhdr.annule_hard() ;
500
501 int m_q, l_q, base_r ;
502 for (int lz=0; lz<nz; lz++) {
503 int np = mp->get_mg()->get_np(lz) ;
504 int nt = mp->get_mg()->get_nt(lz) ;
505 int nr = mp->get_mg()->get_nr(lz) ;
506 for (int k=0; k<np+1; k++)
507 for (int j=0; j<nt; j++) {
508 base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
509 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q > 1))
510 {
511 for (int i=0; i<nr; i++)
512 resu.set_spectral_va().c_cf->set(lz, k, j, i)
513 -= 0.5*(*hsr.get_spectral_va().c_cf)(lz, k, j, i)
514 + 0.5/double(l_q+1)*(*dhdr.get_spectral_va().c_cf)(lz, k, j, i);
515 }
516 }
517 }
518 resu.set_dzpuis(dzp_resu) ;
519 if (resu.set_spectral_va().c != 0x0)
520 delete resu.set_spectral_va().c ;
521 resu.set_spectral_va().c = 0x0 ;
522 } //End of resu != 0
523
524 if (output_ylm) resu.set_spectral_va().ylm() ;
525 else resu.set_spectral_va().ylm_i() ;
526
527 return resu ;
528
529}
530
532 tras) const {
533
534 // All this has a meaning only for spherical components:
535 assert(dynamic_cast<const Base_vect_spher*>(triad) != 0x0) ;
536
537 Scalar resu = tbtt ;
538 if (resu.get_etat() == ETATZERO) {
539 if (tras.get_etat() == ETATZERO) return resu ;
540 else {
541 resu.annule_hard() ;
542 Base_val base = tras.get_spectral_base() ;
543 base.mult_x() ;
544 resu.set_spectral_base(base) ;
545 }
546 }
547 resu.set_spectral_va().ylm() ;
548 int dzp = operator()(1,1).get_dzpuis() ;
549 int dzp_resu = ((dzp == 0) ? 2 : dzp+1) ;
550
551 Scalar hsr = tras ;
552 if (hsr.get_etat() == ETATZERO) return resu ;
553 else {
554 Scalar dhdr = hsr.dsdr() ;
555 if (dhdr.get_etat() == ETATZERO) dhdr.annule_hard() ;
556 dhdr.set_spectral_va().ylm() ;
557 hsr.div_r_dzpuis(dzp_resu) ;
558 hsr.set_spectral_va().ylm() ;
559
560 int nz = mp->get_mg()->get_nzone() ;
561 const Base_val& base = resu.get_spectral_base() ;
562
563 int m_q, l_q, base_r ;
564 for (int lz=0; lz<nz; lz++) {
565 if ((*resu.get_spectral_va().c_cf)(lz).get_etat() == ETATZERO)
566 resu.get_spectral_va().c_cf->set(lz).annule_hard() ;
567 int np = mp->get_mg()->get_np(lz) ;
568 int nt = mp->get_mg()->get_nt(lz) ;
569 int nr = mp->get_mg()->get_nr(lz) ;
570 for (int k=0; k<np+1; k++)
571 for (int j=0; j<nt; j++) {
572 base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
573 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q > 1))
574 {
575 for (int i=0; i<nr; i++)
576 resu.set_spectral_va().c_cf->set(lz, k, j, i)
577 += 0.5*(*hsr.get_spectral_va().c_cf)(lz, k, j, i)
578 + 0.5/double(l_q+1)*(*dhdr.get_spectral_va().c_cf)(lz, k, j, i);
579 }
580 }
581 }
582 resu.set_dzpuis(dzp_resu) ;
583 if (resu.set_spectral_va().c != 0x0)
584 delete resu.set_spectral_va().c ;
585 resu.set_spectral_va().c = 0x0 ;
586 } //End of trace != 0
587 return resu ;
588}
589
590
591 //--------------//
592 // tilde(C) //
593 //--------------//
594
595
597
598 if (p_tilde_c == 0x0) { // a new computation is necessary
599
600 // All this has a meaning only for spherical components:
601 assert(dynamic_cast<const Base_vect_spher*>(triad) != 0x0) ;
602
603 int dzp = operator()(1,1).get_dzpuis() ;
604 int dzp_resu = ((dzp == 0) ? 2 : dzp+1) ;
605
606 p_tilde_c = new Scalar(*mp) ;
608
610 source_eta.div_tant() ;
611 source_eta += operator()(1,2).dsdt() + operator()(1,3).stdsdp() ;
612
613 // Resolution of the angular Poisson equation for eta / r
614 // ------------------------------------------------------
616 tilde_eta = 0 ;
617
618 if (dynamic_cast<const Map_af*>(mp) != 0x0) {
619 tilde_eta = source_eta.poisson_angu() ;
620 }
621 else {
622 assert (par != 0x0) ;
623 mp->poisson_angu(source_eta, *par, tilde_eta) ;
624 }
625
626 Scalar dwdr = www().dsdr() ;
627 dwdr.set_spectral_va().ylm() ;
628 Scalar wsr = www() ;
629 wsr.div_r_dzpuis(dzp_resu) ;
630 wsr.set_spectral_va().ylm() ;
632 etasr2.div_r_dzpuis(dzp_resu) ;
633 etasr2.set_spectral_va().ylm() ;
634 Scalar dtdr = ttt().dsdr() ;
635 dtdr.set_spectral_va().ylm() ;
636 Scalar tsr = ttt() ;
637 tsr.div_r_dzpuis(dzp_resu) ;
638 tsr.set_spectral_va().ylm() ;
639 Scalar hrrsr = operator()(1,1) ;
640 hrrsr.div_r_dzpuis(dzp_resu) ;
641 hrrsr.set_spectral_va().ylm() ;
642
643 int nz = mp->get_mg()->get_nzone() ;
644
645 Base_val base(nz) ;
646
647 if (wsr.get_etat() != ETATZERO) {
648 base = wsr.get_spectral_base() ;
649 }
650 else {
651 if (etasr2.get_etat() != ETATZERO) {
652 base = etasr2.get_spectral_base() ;
653 }
654 else {
655 if (tsr.get_etat() != ETATZERO) {
656 base = tsr.get_spectral_base() ;
657 }
658 else {
659 if (hrrsr.get_etat() != ETATZERO) {
660 base = hrrsr.get_spectral_base() ;
661 }
662 else { //Everything is zero...
664 return *p_tilde_c ;
665 }
666 }
667 }
668 }
669
673
674 if (dwdr.get_etat() == ETATZERO) dwdr.annule_hard() ;
675 if (wsr.get_etat() == ETATZERO) wsr.annule_hard() ;
676 if (etasr2.get_etat() == ETATZERO) etasr2.annule_hard() ;
677 if (dtdr.get_etat() == ETATZERO) dtdr.annule_hard() ;
678 if (tsr.get_etat() == ETATZERO) tsr.annule_hard() ;
679 if (hrrsr.get_etat() == ETATZERO) hrrsr.annule_hard() ;
680
681 int m_q, l_q, base_r ;
682 for (int lz=0; lz<nz; lz++) {
683 int np = mp->get_mg()->get_np(lz) ;
684 int nt = mp->get_mg()->get_nt(lz) ;
685 int nr = mp->get_mg()->get_nr(lz) ;
686 for (int k=0; k<np+1; k++)
687 for (int j=0; j<nt; j++) {
688 base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
689 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q > 1))
690 {
691 for (int i=0; i<nr; i++)
693 = - (l_q - 1)*(*dwdr.get_spectral_va().c_cf)(lz, k, j, i)
694 + (l_q + 1)*(l_q - 1)*(*wsr.get_spectral_va().c_cf)(lz, k, j, i)
695 - 2*(*etasr2.get_spectral_va().c_cf)(lz, k, j, i)
696 + 0.5*double(l_q-1)/double(l_q)*(*tsr.get_spectral_va().c_cf)(lz, k, j, i)
697 - 0.5/double(l_q)*(*dtdr.get_spectral_va().c_cf)(lz, k, j, i)
698 + 1./double(l_q)*(*hrrsr.get_spectral_va().c_cf)(lz, k, j, i) ;
699 }
700 }
701 }
703 } //End of p_tilde_c == 0x0
704
707
708 return *p_tilde_c ;
709
710}
711}
Bases of the spectral expansions.
Definition base_val.h:322
void mult_x()
The basis is transformed as with a multiplication by .
void give_quant_numbers(int, int, int, int &, int &, int &) const
Computes the various quantum numbers and 1d radial base.
Spherical orthonormal vectorial bases (triads).
Definition base_vect.h:308
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
Affine radial mapping.
Definition map.h:2027
Tbl & set(int l)
Read/write of the Tbl containing the coefficients in a given domain.
Definition mtbl_cf.h:294
void annule_hard()
Sets the Mtbl_cf to zero in a hard way.
Definition mtbl_cf.C:312
Parameter storage.
Definition param.h:125
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
void annule_l(int l_min, int l_max, bool ylm_output=false)
Sets all the multipolar components between l_min and l_max to zero.
const Scalar & dsdt() const
Returns of *this .
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition scalar.C:353
const Scalar & stdsdp() const
Returns of *this .
virtual void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition scalar.C:324
Valeur & set_spectral_va()
Returns va (read/write version)
Definition scalar.h:604
void set_dzpuis(int)
Modifies the dzpuis flag.
Definition scalar.C:808
void mult_r_dzpuis(int ced_mult_r)
Multiplication by r everywhere but with the output flag dzpuis set to ced_mult_r .
void set_spectral_base(const Base_val &)
Sets the spectral bases of the Valeur va
Definition scalar.C:797
const Scalar & xxx() const
Gives the field X (see member p_xxx ).
const Scalar & compute_tilde_C(bool output_ylm=true, Param *par=0x0) const
Gives the field (see member p_tilde_c ).
Scalar * p_ttt
Field T defined as .
Definition sym_tensor.h:315
Scalar * p_aaa
Field A defined from X and insensitive to the longitudinal part of the Sym_tensor (only for ).
Definition sym_tensor.h:322
const Scalar & compute_A(bool output_ylm=true, Param *par=0x0) const
Gives the field A (see member p_aaa ).
Scalar * p_tilde_b
Field defined from and h insensitive to the longitudinal part of the Sym_tensor.
Definition sym_tensor.h:334
Scalar * p_mu
Field such that the components of the tensor are written (has only meaning with spherical component...
Definition sym_tensor.h:274
Scalar * p_tilde_c
Field defined from and h insensitive to the longitudinal part of the Sym_tensor.
Definition sym_tensor.h:346
const Scalar & ttt() const
Gives the field T (see member p_ttt ).
virtual const Scalar & eta(Param *par=0x0) const
Gives the field (see member p_eta ).
const Scalar & www() const
Gives the field W (see member p_www ).
Scalar get_tilde_B_from_TT_trace(const Scalar &tilde_B_tt_in, const Scalar &trace) const
Computes (see Sym_tensor::p_tilde_b ) from its transverse-traceless part and the trace.
const Scalar & mu(Param *par=0x0) const
Gives the field (see member p_mu ).
virtual void del_deriv() const
Deletes the derived quantities.
Definition sym_tensor.C:286
Scalar compute_tilde_B_tt(bool output_ylm=true, Param *par=0x0) const
Gives the field (see member p_tilde_b ) associated with the TT-part of the Sym_tensor .
Scalar * p_eta
Field such that the components of the tensor are written (has only meaning with spherical component...
Definition sym_tensor.h:260
Scalar * p_xxx
Field X such that the components and of the tensor are written (has only meaning with spherical com...
Definition sym_tensor.h:312
void set_auxiliary(const Scalar &trr, const Scalar &eta_over_r, const Scalar &mu_over_r, const Scalar &www, const Scalar &xxx, const Scalar &ttt)
Assigns the component and the derived members p_eta , p_mu , p_www, p_xxx and p_ttt ,...
Scalar * p_www
Field W such that the components and of the tensor are written (has only meaning with spherical com...
Definition sym_tensor.h:293
const Scalar & compute_tilde_B(bool output_ylm=true, Param *par=0x0) const
Gives the field (see member p_tilde_b ).
void set_etat_cf_qcq()
Sets the logical state to ETATQCQ (ordinary state) for values in the configuration space (Mtbl_cf c_c...
Definition valeur.C:712
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 ylm_i()
Inverse of ylm()
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition tensor.h:295
const Scalar & operator()(const Itbl &ind) const
Returns the value of a component (read-only version).
Definition tensor.C:798
Scalar & set(const Itbl &ind)
Returns the value of a component (read/write version).
Definition tensor.C:654
const Base_vect * triad
Vectorial basis (triad) with respect to which the tensor components are defined.
Definition tensor.h:303
Lorene prototypes.
Definition app_hor.h:64