LORENE
op_lapang.C
1/*
2 * Copyright (c) 1999-2001 Eric Gourgoulhon
3 *
4 * This file is part of LORENE.
5 *
6 * LORENE is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * LORENE is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with LORENE; if not, write to the Free Software
18 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 *
20 */
21
22
23char op_lapang_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/op_lapang.C,v 1.9 2014/10/13 08:53:25 j_novak Exp $" ;
24
25/*
26 * Ensemble des routines de base pour le calcul du laplacien angulaire,
27 * c'est-a-dire de l'operateur
28 *
29 * d^2/dtheta^2 + cos(theta)/sin(theta) d/dtheta + 1/sin(theta) d^2/dphi^2
30 *
31 * (Utilisation interne)
32 *
33 * void _lapang_XXXX(Mtbl_cf * mt, int l)
34 * mt pointeur sur le Mtbl_cf d'entree-sortie
35 * l indice de la zone ou l'on doit effectuer le calcul
36 *
37 */
38
39/*
40 * $Id: op_lapang.C,v 1.9 2014/10/13 08:53:25 j_novak Exp $
41 * $Log: op_lapang.C,v $
42 * Revision 1.9 2014/10/13 08:53:25 j_novak
43 * Lorene classes and functions now belong to the namespace Lorene.
44 *
45 * Revision 1.8 2014/10/06 15:16:06 j_novak
46 * Modified #include directives to use c++ syntax.
47 *
48 * Revision 1.7 2009/10/23 12:55:04 j_novak
49 * New base T_LEG_MI
50 *
51 * Revision 1.6 2009/10/13 19:45:01 j_novak
52 * New base T_LEG_MP.
53 *
54 * Revision 1.5 2005/05/18 07:47:36 j_novak
55 * Corrected an error for the T_LEG_II base (ll was set to 2j+1 instead of 2j for
56 * sin(phi)).
57 *
58 * Revision 1.4 2004/12/17 13:20:55 m_forot
59 * Add T_LEG base
60 *
61 * Revision 1.3 2003/09/16 12:11:59 j_novak
62 * Added the base T_LEG_II.
63 *
64 * Revision 1.2 2002/10/16 14:36:58 j_novak
65 * Reorganization of #include instructions of standard C++, in order to
66 * use experimental version 3 of gcc.
67 *
68 * Revision 1.1.1.1 2001/11/20 15:19:29 e_gourgoulhon
69 * LORENE
70 *
71 * Revision 2.2 2000/11/14 15:09:08 eric
72 * Traitement du cas np=1 dans T_LEG_PI
73 *
74 * Revision 2.1 2000/10/04 14:54:59 eric
75 * Ajout des bases T_LEG_IP et T_LEG_PI.
76 *
77 * Revision 2.0 1999/04/26 16:42:04 phil
78 * *** empty log message ***
79 *
80 *
81 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/op_lapang.C,v 1.9 2014/10/13 08:53:25 j_novak Exp $
82 *
83 */
84
85// Headers C
86#include <cstdlib>
87
88// Headers Lorene
89#include "mtbl_cf.h"
90#include "grilles.h"
91#include "type_parite.h"
92
93
94 //--------------------------------------
95 // Routine pour les cas non prevus ----
96 //------------------------------------
97
98namespace Lorene {
99void _lapang_pas_prevu(Mtbl_cf* mt, int l) {
100 cout << "Unknwon theta basis in the operator Mtbl_cf::lapang() !" << endl ;
101 cout << " basis : " << hex << (mt->base).b[l] << endl ;
102 abort () ;
103}
104 //---------------
105 // cas T_LEG --
106 //---------------
107
108void _lapang_t_leg(Mtbl_cf* mt, int l)
109{
110
111 Tbl* tb = mt->t[l] ; // pt. sur tbl de travail
112
113 // Peut-etre rien a faire ?
114 if (tb->get_etat() == ETATZERO) {
115 return ;
116 }
117
118 int k, j, i ;
119 // Pour le confort
120 int nr = mt->get_mg()->get_nr(l) ; // Nombre
121 int nt = mt->get_mg()->get_nt(l) ; // de points
122 int np = mt->get_mg()->get_np(l) ; // physiques
123
124 int np1 = ( np == 1 ) ? 1 : np+1 ;
125
126 double* tuu = tb->t ;
127
128 // k = 0 :
129
130 for (j=0 ; j<nt ; j++) {
131 int ll = j ;
132 double xl = - ll*(ll+1) ;
133 for (i=0 ; i<nr ; i++) {
134 tuu[i] *= xl ;
135 } // Fin de boucle sur r
136 tuu += nr ;
137 } // Fin de boucle sur theta
138
139 // On saute k = 1 :
140 tuu += nt*nr ;
141
142 // k=2,...
143 for (k=2 ; k<np1 ; k++) {
144 int m = k/2 ;
145 tuu += (m/2)*nr ;
146 for (j=m/2 ; j<nt ; j++) {
147 int ll = j;
148 double xl = - ll*(ll+1) ;
149 for (i=0 ; i<nr ; i++) {
150 tuu[i] *= xl ;
151 } // Fin de boucle sur r
152 tuu += nr ;
153 } // Fin de boucle sur theta
154 } // Fin de boucle sur phi
155
156 // base de developpement inchangee
157}
158
159 //---------------
160 // cas T_LEG_P --
161 //---------------
162
163void _lapang_t_leg_p(Mtbl_cf* mt, int l)
164{
165
166 Tbl* tb = mt->t[l] ; // pt. sur tbl de travail
167
168 // Peut-etre rien a faire ?
169 if (tb->get_etat() == ETATZERO) {
170 return ;
171 }
172
173 int k, j, i ;
174 // Pour le confort
175 int nr = mt->get_mg()->get_nr(l) ; // Nombre
176 int nt = mt->get_mg()->get_nt(l) ; // de points
177 int np = mt->get_mg()->get_np(l) ; // physiques
178
179 int np1 = ( np == 1 ) ? 1 : np+1 ;
180
181 double* tuu = tb->t ;
182
183 // k = 0 :
184
185 for (j=0 ; j<nt ; j++) {
186 int ll = 2*j ;
187 double xl = - ll*(ll+1) ;
188 for (i=0 ; i<nr ; i++) {
189 tuu[i] *= xl ;
190 } // Fin de boucle sur r
191 tuu += nr ;
192 } // Fin de boucle sur theta
193
194 // On saute k = 1 :
195 tuu += nt*nr ;
196
197 // k=2,...
198 for (k=2 ; k<np1 ; k++) {
199 int m = k/2 ;
200 tuu += (m/2)*nr ;
201 for (j=m/2 ; j<nt ; j++) {
202 int ll = 2*j + (m%2) ;
203 double xl = - ll*(ll+1) ;
204 for (i=0 ; i<nr ; i++) {
205 tuu[i] *= xl ;
206 } // Fin de boucle sur r
207 tuu += nr ;
208 } // Fin de boucle sur theta
209 } // Fin de boucle sur phi
210
211 // base de developpement inchangee
212}
213
214 //------------------
215 // cas T_LEG_PP --
216 //----------------
217
218void _lapang_t_leg_pp(Mtbl_cf* mt, int l)
219{
220
221 Tbl* tb = mt->t[l] ; // pt. sur tbl de travail
222
223 // Peut-etre rien a faire ?
224 if (tb->get_etat() == ETATZERO) {
225 return ;
226 }
227
228 int k, j, i ;
229 // Pour le confort
230 int nr = mt->get_mg()->get_nr(l) ; // Nombre
231 int nt = mt->get_mg()->get_nt(l) ; // de points
232 int np = mt->get_mg()->get_np(l) ; // physiques
233
234 int np1 = ( np == 1 ) ? 1 : np+1 ;
235
236 double* tuu = tb->t ;
237
238 // k = 0 :
239
240 for (j=0 ; j<nt ; j++) {
241 int ll = 2*j ;
242 double xl = - ll*(ll+1) ;
243 for (i=0 ; i<nr ; i++) {
244 tuu[i] *= xl ;
245 } // Fin de boucle sur r
246 tuu += nr ;
247 } // Fin de boucle sur theta
248
249 // On saute k = 1 :
250 tuu += nt*nr ;
251
252 // k=2,...
253 for (k=2 ; k<np1 ; k++) {
254 int m = 2*(k/2);
255 tuu += (m/2)*nr ;
256 for (j=m/2 ; j<nt ; j++) {
257 int ll = 2*j ;
258 double xl = - ll*(ll+1) ;
259 for (i=0 ; i<nr ; i++) {
260 tuu[i] *= xl ;
261 } // Fin de boucle sur r
262 tuu += nr ;
263 } // Fin de boucle sur theta
264 } // Fin de boucle sur phi
265
266 // base de developpement inchangee
267}
268
269 //----------------
270 // cas T_LEG_I --
271 //---------------
272
273void _lapang_t_leg_i(Mtbl_cf* mt, int l)
274{
275
276 Tbl* tb = mt->t[l] ; // pt. sur tbl de travail
277
278 // Peut-etre rien a faire ?
279 if (tb->get_etat() == ETATZERO) {
280 return ;
281 }
282
283 int k, j, i ;
284 // Pour le confort
285 int nr = mt->get_mg()->get_nr(l) ; // Nombre
286 int nt = mt->get_mg()->get_nt(l) ; // de points
287 int np = mt->get_mg()->get_np(l) ; // physiques
288
289 int np1 = ( np == 1 ) ? 1 : np+1 ;
290
291 double* tuu = tb->t ;
292
293 // k = 0 :
294
295 for (j=0 ; j<nt-1 ; j++) {
296 int ll = 2*j+1 ;
297 double xl = - ll*(ll+1) ;
298 for (i=0 ; i<nr ; i++) {
299 tuu[i] *= xl ;
300 } // Fin de boucle sur r
301 tuu += nr ;
302 } // Fin de boucle sur theta
303 tuu += nr ; // On saute j=nt-1
304
305 // On saute k = 1 :
306 tuu += nt*nr ;
307
308 // k=2,...
309 for (k=2 ; k<np1 ; k++) {
310 int m = k/2 ;
311 tuu += ((m+1)/2)*nr ;
312 for (j=(m+1)/2 ; j<nt-1 ; j++) {
313 int ll = 2*j + ((m+1)%2) ;
314 double xl = - ll*(ll+1) ;
315 for (i=0 ; i<nr ; i++) {
316 tuu[i] *= xl ;
317 } // Fin de boucle sur r
318 tuu += nr ;
319 } // Fin de boucle sur theta
320 tuu += nr ; // On saute j=nt-1
321 } // Fin de boucle sur phi
322
323 // base de developpement inchangee
324}
325
326 //------------------
327 // cas T_LEG_IP --
328 //----------------
329
330void _lapang_t_leg_ip(Mtbl_cf* mt, int l)
331{
332
333 Tbl* tb = mt->t[l] ; // pt. sur tbl de travail
334
335 // Peut-etre rien a faire ?
336 if (tb->get_etat() == ETATZERO) {
337 return ;
338 }
339
340 int k, j, i ;
341 // Pour le confort
342 int nr = mt->get_mg()->get_nr(l) ; // Nombre
343 int nt = mt->get_mg()->get_nt(l) ; // de points
344 int np = mt->get_mg()->get_np(l) ; // physiques
345
346 int np1 = ( np == 1 ) ? 1 : np+1 ;
347
348 double* tuu = tb->t ;
349
350 // k = 0 :
351
352 for (j=0 ; j<nt-1 ; j++) {
353 int ll = 2*j+1 ;
354 double xl = - ll*(ll+1) ;
355 for (i=0 ; i<nr ; i++) {
356 tuu[i] *= xl ;
357 } // Fin de boucle sur r
358 tuu += nr ;
359 } // Fin de boucle sur theta
360 tuu += nr ; // On saute j=nt-1
361
362 // On saute k = 1 :
363 tuu += nt*nr ;
364
365 // k=2,...
366 for (k=2 ; k<np1 ; k++) {
367 int m = 2*(k/2);
368 tuu += (m/2)*nr ;
369 for (j=m/2 ; j<nt-1 ; j++) {
370 int ll = 2*j+1 ;
371 double xl = - ll*(ll+1) ;
372 for (i=0 ; i<nr ; i++) {
373 tuu[i] *= xl ;
374 } // Fin de boucle sur r
375 tuu += nr ;
376 } // Fin de boucle sur theta
377 tuu += nr ; // On saute j=nt-1
378 } // Fin de boucle sur phi
379
380//## Verif
381 assert (tuu == tb->t + (np+1)*nt*nr) ;
382
383 // base de developpement inchangee
384}
385
386 //------------------
387 // cas T_LEG_PI --
388 //----------------
389
390void _lapang_t_leg_pi(Mtbl_cf* mt, int l)
391{
392
393 Tbl* tb = mt->t[l] ; // pt. sur tbl de travail
394
395 // Peut-etre rien a faire ?
396 if (tb->get_etat() == ETATZERO) {
397 return ;
398 }
399
400 int k, j, i ;
401 // Pour le confort
402 int nr = mt->get_mg()->get_nr(l) ; // Nombre
403 int nt = mt->get_mg()->get_nt(l) ; // de points
404 int np = mt->get_mg()->get_np(l) ; // physiques
405
406 double* tuu = tb->t ;
407
408 // k = 0 : cos(phi)
409 // -----
410
411 for (j=0 ; j<nt-1 ; j++) {
412 int ll = 2*j+1 ;
413 double xl = - ll*(ll+1) ;
414 for (i=0 ; i<nr ; i++) {
415 tuu[i] *= xl ;
416 } // Fin de boucle sur r
417 tuu += nr ;
418 } // Fin de boucle sur theta
419 tuu += nr ; // On saute j=nt-1
420
421 if (np==1) {
422 return ;
423 }
424
425 // k = 1 : on saute
426 // -----
427 tuu += nt*nr ;
428
429 // k = 2 : sin(phi)
430 // ------
431 for (j=0 ; j<nt-1 ; j++) {
432 int ll = 2*j+1 ;
433 double xl = - ll*(ll+1) ;
434 for (i=0 ; i<nr ; i++) {
435 tuu[i] *= xl ;
436 } // Fin de boucle sur r
437 tuu += nr ;
438 } // Fin de boucle sur theta
439 tuu += nr ; // On saute j=nt-1
440
441 // 3 <= k <= np
442 // ------------
443 for (k=3 ; k<np+1 ; k++) {
444 int m = (k%2 == 0) ? k-1 : k ;
445 tuu += (m-1)/2*nr ;
446 for (j=(m-1)/2 ; j<nt-1 ; j++) {
447 int ll = 2*j+1 ;
448 double xl = - ll*(ll+1) ;
449 for (i=0 ; i<nr ; i++) {
450 tuu[i] *= xl ;
451 } // Fin de boucle sur r
452 tuu += nr ;
453 } // Fin de boucle sur theta
454 tuu += nr ; // On saute j=nt-1
455 } // Fin de boucle sur phi
456
457//## Verif
458 assert (tuu == tb->t + (np+1)*nt*nr) ;
459
460 // base de developpement inchangee
461}
462
463 //------------------
464 // cas T_LEG_II --
465 //----------------
466
467void _lapang_t_leg_ii(Mtbl_cf* mt, int l)
468{
469
470 Tbl* tb = mt->t[l] ; // pt. sur tbl de travail
471
472 // Peut-etre rien a faire ?
473 if (tb->get_etat() == ETATZERO) {
474 return ;
475 }
476
477 int k, j, i ;
478 // Pour le confort
479 int nr = mt->get_mg()->get_nr(l) ; // Nombre
480 int nt = mt->get_mg()->get_nt(l) ; // de points
481 int np = mt->get_mg()->get_np(l) ; // physiques
482
483 double* tuu = tb->t ;
484
485 // k = 0 : cos(phi)
486 // -----
487
488 for (j=0 ; j<nt-1 ; j++) {
489 int ll = 2*j ;
490 double xl = - ll*(ll+1) ;
491 for (i=0 ; i<nr ; i++) {
492 tuu[i] *= xl ;
493 } // Fin de boucle sur r
494 tuu += nr ;
495 } // Fin de boucle sur theta
496 tuu += nr ; // On saute j=nt-1
497
498 if (np==1) {
499 return ;
500 }
501
502 // k = 1 : on saute
503 // -----
504 tuu += nt*nr ;
505
506 // k = 2 : sin(phi)
507 // ------
508 for (j=0 ; j<nt-1 ; j++) {
509 int ll = 2*j ;
510 double xl = - ll*(ll+1) ;
511 for (i=0 ; i<nr ; i++) {
512 tuu[i] *= xl ;
513 } // Fin de boucle sur r
514 tuu += nr ;
515 } // Fin de boucle sur theta
516 tuu += nr ; // On saute j=nt-1
517
518 // 3 <= k <= np
519 // ------------
520 for (k=3 ; k<np+1 ; k++) {
521 int m = (k%2 == 0) ? k-1 : k ;
522 tuu += (m+1)/2*nr ;
523 for (j=(m+1)/2 ; j<nt-1 ; j++) {
524 int ll = 2*j ;
525 double xl = - ll*(ll+1) ;
526 for (i=0 ; i<nr ; i++) {
527 tuu[i] *= xl ;
528 } // Fin de boucle sur r
529 tuu += nr ;
530 } // Fin de boucle sur theta
531 tuu += nr ; // On saute j=nt-1
532 } // Fin de boucle sur phi
533
534//## Verif
535 assert (tuu == tb->t + (np+1)*nt*nr) ;
536
537 // base de developpement inchangee
538}
539
540 //----------------
541 // cas T_LEG_MP --
542 //----------------
543
544void _lapang_t_leg_mp(Mtbl_cf* mt, int l)
545{
546
547 Tbl* tb = mt->t[l] ; // pt. sur tbl de travail
548
549 // Peut-etre rien a faire ?
550 if (tb->get_etat() == ETATZERO) {
551 return ;
552 }
553
554 int k, j, i ;
555 // Pour le confort
556 int nr = mt->get_mg()->get_nr(l) ; // Nombre
557 int nt = mt->get_mg()->get_nt(l) ; // de points
558 int np = mt->get_mg()->get_np(l) ; // physiques
559
560 int np1 = ( np == 1 ) ? 1 : np+1 ;
561
562 double* tuu = tb->t ;
563
564 // k = 0 :
565
566 for (j=0 ; j<nt ; j++) {
567 int ll = j ;
568 double xl = - ll*(ll+1) ;
569 for (i=0 ; i<nr ; i++) {
570 tuu[i] *= xl ;
571 } // Fin de boucle sur r
572 tuu += nr ;
573 } // Fin de boucle sur theta
574
575 // On saute k = 1 :
576 tuu += nt*nr ;
577
578 // k=2,...
579 for (k=2 ; k<np1 ; k++) {
580 int m = 2*(k/2);
581 tuu += m*nr ;
582 for (j=m ; j<nt ; j++) {
583 int ll = j ;
584 double xl = - ll*(ll+1) ;
585 for (i=0 ; i<nr ; i++) {
586 tuu[i] *= xl ;
587 } // Fin de boucle sur r
588 tuu += nr ;
589 } // Fin de boucle sur theta
590 } // Fin de boucle sur phi
591
592 // base de developpement inchangee
593}
594 //----------------
595 // cas T_LEG_MI --
596 //----------------
597
598void _lapang_t_leg_mi(Mtbl_cf* mt, int l)
599{
600
601 Tbl* tb = mt->t[l] ; // pt. sur tbl de travail
602
603 // Peut-etre rien a faire ?
604 if (tb->get_etat() == ETATZERO) {
605 return ;
606 }
607
608 int k, j, i ;
609 // Pour le confort
610 int nr = mt->get_mg()->get_nr(l) ; // Nombre
611 int nt = mt->get_mg()->get_nt(l) ; // de points
612 int np = mt->get_mg()->get_np(l) ; // physiques
613
614 int np1 = ( np == 1 ) ? 1 : np+1 ;
615
616 double* tuu = tb->t ;
617
618 // k = 0 :
619
620 for (j=0 ; j<nt ; j++) {
621 int ll = j ;
622 double xl = - ll*(ll+1) ;
623 for (i=0 ; i<nr ; i++) {
624 tuu[i] *= xl ;
625 } // Fin de boucle sur r
626 tuu += nr ;
627 } // Fin de boucle sur theta
628
629 // On saute k = 1 :
630 tuu += nt*nr ;
631
632 // k=2,...
633 for (k=2 ; k<np1 ; k++) {
634 int m = 2*((k-1)/2) + 1;
635 tuu += m*nr ;
636 for (j=m ; j<nt ; j++) {
637 int ll = j ;
638 double xl = - ll*(ll+1) ;
639 for (i=0 ; i<nr ; i++) {
640 tuu[i] *= xl ;
641 } // Fin de boucle sur r
642 tuu += nr ;
643 } // Fin de boucle sur theta
644 } // Fin de boucle sur phi
645
646 // base de developpement inchangee
647}
648}
Lorene prototypes.
Definition app_hor.h:64