LORENE
som_r.C
1/*
2 * Copyright (c) 1999-2000 Jean-Alain Marck
3 * Copyright (c) 1999-2001 Philippe Grandclement
4 * Copyright (c) 1999-2002 Eric Gourgoulhon
5 *
6 * This file is part of LORENE.
7 *
8 * LORENE is free software; you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation; either version 2 of the License, or
11 * (at your option) any later version.
12 *
13 * LORENE is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public License
19 * along with LORENE; if not, write to the Free Software
20 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
21 *
22 */
23
24
25char som_r_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/som_r.C,v 1.11 2014/10/13 08:53:27 j_novak Exp $" ;
26
27/*
28 * Ensemble des routine pour la sommation directe en r
29 *
30 * SYNOPSYS:
31 * double som_r_XX
32 * (double* ti, int nr, int nt, int np, double x, double* trtp)
33 *
34 * x est l'argument du polynome de Chebychev: x in [0, 1] ou x in [-1, 1]
35 *
36 * ATTENTION: np est suppose etre le nombre de points (sans le +2)
37 * on suppose que trtp tient compte de ca.
38 */
39
40
41/*
42 * $Id: som_r.C,v 1.11 2014/10/13 08:53:27 j_novak Exp $
43 * $Log: som_r.C,v $
44 * Revision 1.11 2014/10/13 08:53:27 j_novak
45 * Lorene classes and functions now belong to the namespace Lorene.
46 *
47 * Revision 1.10 2014/10/06 15:16:07 j_novak
48 * Modified #include directives to use c++ syntax.
49 *
50 * Revision 1.9 2013/06/13 14:18:18 j_novak
51 * Inclusion of new bases R_LEG, R_LEGP and R_LEGI.
52 *
53 * Revision 1.8 2008/08/27 08:50:10 jl_cornou
54 * Added Jacobi(0,2) polynomials
55 *
56 * Revision 1.7 2007/12/11 15:28:18 jl_cornou
57 * Jacobi(0,2) polynomials partially implemented
58 *
59 * Revision 1.6 2006/05/17 13:15:19 j_novak
60 * Added a test for the pure angular grid case (nr=1), in the shell.
61 *
62 * Revision 1.5 2005/02/18 13:14:17 j_novak
63 * Changing of malloc/free to new/delete + suppression of some unused variables
64 * (trying to avoid compilation warnings).
65 *
66 * Revision 1.4 2004/11/23 15:16:02 m_forot
67 *
68 * Added the bases for the cases without any equatorial symmetry
69 * (T_COSSIN_C, T_COSSIN_S, T_LEG, R_CHEBPI_P, R_CHEBPI_I).
70 *
71 * Revision 1.3 2002/10/16 14:36:59 j_novak
72 * Reorganization of #include instructions of standard C++, in order to
73 * use experimental version 3 of gcc.
74 *
75 * Revision 1.2 2002/05/05 16:20:40 e_gourgoulhon
76 * Error message (in unknwon basis case) in English
77 * Added the basis T_COSSIN_SP
78 *
79 * Revision 1.1.1.1 2001/11/20 15:19:29 e_gourgoulhon
80 * LORENE
81 *
82 * Revision 2.4 2000/03/06 09:34:21 eric
83 * Suppression des #include inutiles.
84 *
85 * Revision 2.3 1999/04/28 12:11:15 phil
86 * changements de sommations
87 *
88 * Revision 2.2 1999/04/26 14:26:31 phil
89 * remplacement des malloc en new
90 *
91 * Revision 2.1 1999/04/13 14:44:06 phil
92 * ajout de som_r_chebi
93 *
94 * Revision 2.0 1999/04/12 15:42:33 phil
95 * *** empty log message ***
96 *
97 * Revision 1.1 1999/04/12 15:40:25 phil
98 * Initial revision
99 *
100 *
101 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/som_r.C,v 1.11 2014/10/13 08:53:27 j_novak Exp $
102 *
103 */
104// Headers C
105#include <cstdlib>
106
107#include "headcpp.h"
108#include "proto.h"
109
110namespace Lorene {
111 //-------------------
112 //- Cas Non-Prevu ---
113 //-------------------
114
115void som_r_pas_prevu
116 (double*, const int, const int, const int, const double, double*) {
117 cout << "Mtbl_cf::val_point: r basis not implemented yet !"
118 << endl ;
119 abort() ;
120}
121
122 //----------------
123 //- Cas R_CHEB ---
124 //----------------
125
126void som_r_cheb
127 (double* ti, const int nr, const int nt, const int np, const double x,
128 double* trtp) {
129
130// Variables diverses
131int i, j, k ;
132double* pi = ti ; // pointeur courant sur l'entree
133double* po = trtp ; // pointeur courant sur la sortie
134
135 // Valeurs des polynomes de Chebyshev au point x demande
136 double* cheb = new double [nr] ;
137 cheb[0] = 1. ;
138 if (nr > 1) cheb[1] = x ;
139 for (i=2; i<nr; i++) {
140 cheb[i] = 2*x* cheb[i-1] - cheb[i-2] ;
141 }
142
143 // Sommation pour le premier phi, k=0
144 for (j=0 ; j<nt ; j++) {
145 *po = 0 ;
146 // Sommation sur r
147 for (i=0 ; i<nr ; i++) {
148 *po += (*pi) * cheb[i] ;
149 pi++ ; // R suivant
150 }
151 po++ ; // Theta suivant
152 }
153
154 if (np > 1) {
155
156 // On saute le deuxieme phi (sin(0)), k=1
157 pi += nr*nt ;
158 po += nt ;
159
160 // Sommation sur les phi suivants (pour k=2,...,np)
161 for (k=2 ; k<np+1 ; k++) {
162 for (j=0 ; j<nt ; j++) {
163 // Sommation sur r
164 *po = 0 ;
165 for (i=0 ; i<nr ; i++) {
166 *po += (*pi) * cheb[i] ;
167 pi++ ; // R suivant
168 }
169 po++ ; // Theta suivant
170 }
171 }
172
173 } // fin du cas np > 1
174
175 // Menage
176 delete [] cheb ;
177
178}
179
180
181 //-----------------
182 //- Cas R_CHEBP ---
183 //-----------------
184
185void som_r_chebp
186 (double* ti, const int nr, const int nt, const int np, const double x,
187 double* trtp) {
188
189// Variables diverses
190int i, j, k ;
191double* pi = ti ; // pointeur courant sur l'entree
192double* po = trtp ; // pointeur courant sur la sortie
193
194 // Valeurs des polynomes de Chebyshev au point x demande
195 double* cheb = new double [nr] ;
196 cheb[0] = 1. ;
197 double t2im1 = x ;
198 for (i=1; i<nr; i++) {
199 cheb[i] = 2*x* t2im1 - cheb[i-1] ;
200 t2im1 = 2*x* cheb[i] - t2im1 ;
201 }
202
203 // Sommation pour le premier phi, k=0
204 for (j=0 ; j<nt ; j++) {
205 *po = 0 ;
206 // Sommation sur r
207 for (i=0 ; i<nr ; i++) {
208 *po += (*pi) * cheb[i] ;
209 pi++ ; // R suivant
210 }
211 po++ ; // Theta suivant
212 }
213
214 if (np > 1) {
215
216 // On saute le deuxieme phi (sin(0)), k=1
217 pi += nr*nt ;
218 po += nt ;
219
220 // Sommation sur les phi suivants (pour k=2,...,np)
221 for (k=2 ; k<np+1 ; k++) {
222 for (j=0 ; j<nt ; j++) {
223 // Sommation sur r
224 *po = 0 ;
225 for (i=0 ; i<nr ; i++) {
226 *po += (*pi) * cheb[i] ;
227 pi++ ; // R suivant
228 }
229 po++ ; // Theta suivant
230 }
231 }
232
233 } // fin du cas np > 1
234 // Menage
235 delete [] cheb ;
236
237}
238
239
240 //-----------------
241 //- Cas R_CHEBI ---
242 //-----------------
243
244void som_r_chebi
245 (double* ti, const int nr, const int nt, const int np, const double x,
246 double* trtp) {
247
248// Variables diverses
249int i, j, k ;
250double* pi = ti ; // pointeur courant sur l'entree
251double* po = trtp ; // pointeur courant sur la sortie
252
253 // Valeurs des polynomes de Chebyshev au point x demande
254 double* cheb = new double [nr] ;
255 cheb[0] = x ;
256 double t2im1 = 1. ;
257 for (i=1; i<nr; i++) {
258 t2im1 = 2*x* cheb[i-1] - t2im1 ;
259 cheb[i] = 2*x* t2im1 - cheb[i-1] ;
260 }
261
262 // Sommation pour le premier phi, k=0
263 for (j=0 ; j<nt ; j++) {
264 *po = 0 ;
265 // Sommation sur r
266 for (i=0 ; i<nr ; i++) {
267 *po += (*pi) * cheb[i] ;
268 pi++ ; // R suivant
269 }
270 po++ ; // Theta suivant
271 }
272
273 if (np > 1) {
274
275 // On saute le deuxieme phi (sin(0)), k=1
276 pi += nr*nt ;
277 po += nt ;
278
279 // Sommation sur les phi suivants (pour k=2,...,np)
280 for (k=2 ; k<np+1 ; k++) {
281 for (j=0 ; j<nt ; j++) {
282 // Sommation sur r
283 *po = 0 ;
284 for (i=0 ; i<nr ; i++) {
285 *po += (*pi) * cheb[i] ;
286 pi++ ; // R suivant
287 }
288 po++ ; // Theta suivant
289 }
290 }
291
292 } // fin du cas np > 1
293 // Menage
294 delete [] cheb ;
295
296}
297 //-----------------
298 //- Cas R_CHEBU ---
299 //-----------------
300
301void som_r_chebu
302 (double* ti, const int nr, const int nt, const int np, const double x, double* trtp) {
303
304// Variables diverses
305int i, j, k ;
306double* pi = ti ; // pointeur courant sur l'entree
307double* po = trtp ; // pointeur courant sur la sortie
308
309 // Valeurs des polynomes de Chebyshev au point x demande
310 double* cheb = new double [nr] ;
311 cheb[0] = 1. ;
312 cheb[1] = x ;
313 for (i=2; i<nr; i++) {
314 cheb[i] = 2*x* cheb[i-1] - cheb[i-2] ;
315 }
316
317 // Sommation pour le premier phi, k=0
318 for (j=0 ; j<nt ; j++) {
319 *po = 0 ;
320 // Sommation sur r
321 for (i=0 ; i<nr ; i++) {
322 *po += (*pi) * cheb[i] ;
323 pi++ ; // R suivant
324 }
325 po++ ; // Theta suivant
326 }
327
328 if (np > 1) {
329
330 // On saute le deuxieme phi (sin(0)), k=1
331 pi += nr*nt ;
332 po += nt ;
333
334 // Sommation sur les phi suivants (pour k=2,...,np)
335 for (k=2 ; k<np+1 ; k++) {
336 for (j=0 ; j<nt ; j++) {
337 // Sommation sur r
338 *po = 0 ;
339 for (i=0 ; i<nr ; i++) {
340 *po += (*pi) * cheb[i] ;
341 pi++ ; // R suivant
342 }
343 po++ ; // Theta suivant
344 }
345 }
346
347 } // fin du cas np > 1
348 // Menage
349 delete [] cheb ;
350}
351
352 //----------------------
353 // Cas R_CHEBPIM_P ---
354 //---------------------
355
356void som_r_chebpim_p
357 (double* ti, const int nr, const int nt, const int np, const double x, double* trtp) {
358
359// Variables diverses
360int i, j, k ;
361double* pi = ti ; // pointeur courant sur l'entree
362double* po = trtp ; // pointeur courant sur la sortie
363
364 // Valeurs des polynomes de Chebyshev au point x demande
365 double* cheb = new double [2*nr] ;
366 cheb[0] = 1. ;
367 cheb[1] = x ;
368 for (i=2 ; i<2*nr ; i++) {
369 cheb[i] = 2*x* cheb[i-1] - cheb[i-2] ;
370 }
371
372 // Sommation pour le premier phi, k=0
373 int m = 0;
374 for (j=0 ; j<nt ; j++) {
375 *po = 0 ;
376 // Sommation sur r
377 for (i=0 ; i<nr ; i++) {
378 *po += (*pi) * cheb[2*i] ;
379 pi++ ; // R suivant
380 }
381 po++ ; // Theta suivant
382 }
383
384 if (np > 1) {
385
386 // On saute le deuxieme phi (sin(0)), k=1
387 pi += nr*nt ;
388 po += nt ;
389
390 // Sommation sur les phi suivants (pour k=2,...,np)
391 for (k=2 ; k<np+1 ; k++) {
392 m = (k/2) % 2 ; // parite: 0 <-> 1
393 for (j=0 ; j<nt ; j++) {
394 // Sommation sur r
395 *po = 0 ;
396 for (i=0 ; i<nr ; i++) {
397 *po += (*pi) * cheb[2*i + m] ;
398 pi++ ; // R suivant
399 }
400 po++ ; // Theta suivant
401 }
402 }
403
404 } // fin du cas np > 1
405 // Menage
406 delete [] cheb ;
407
408}
409
410 //----------------------
411 //- Cas R_CHEBPIM_I ----
412 //----------------------
413
414void som_r_chebpim_i
415 (double* ti, const int nr, const int nt, const int np, const double x, double* trtp) {
416
417// Variables diverses
418int i, j, k ;
419double* pi = ti ; // pointeur courant sur l'entree
420double* po = trtp ; // pointeur courant sur la sortie
421
422 // Valeurs des polynomes de Chebyshev au point x demande
423 double* cheb = new double [2*nr] ;
424 cheb[0] = 1. ;
425 cheb[1] = x ;
426 for (i=2 ; i<2*nr ; i++) {
427 cheb[i] = 2*x* cheb[i-1] - cheb[i-2] ;
428 }
429
430 // Sommation pour le premier phi, k=0
431 int m = 0;
432 for (j=0 ; j<nt ; j++) {
433 *po = 0 ;
434 // Sommation sur r
435 for (i=0 ; i<nr ; i++) {
436 *po += (*pi) * cheb[2*i+1] ;
437 pi++ ; // R suivant
438 }
439 po++ ; // Theta suivant
440 }
441
442 if (np > 1) {
443
444 // On saute le deuxieme phi (sin(0)), k=1
445 pi += nr*nt ;
446 po += nt ;
447
448 // Sommation sur les phi suivants (pour k=2,...,np)
449 for (k=2 ; k<np+1 ; k++) {
450 m = (k/2) % 2 ; // parite: 0 <-> 1
451 for (j=0 ; j<nt ; j++) {
452 // Sommation sur r
453 *po = 0 ;
454 for (i=0 ; i<nr ; i++) {
455 *po += (*pi) * cheb[2*i + 1 - m] ;
456 pi++ ; // R suivant
457 }
458 po++ ; // Theta suivant
459 }
460 }
461
462 } // fin du cas np > 1
463 // Menage
464 delete [] cheb ;
465
466}
467
468 //----------------------
469 // Cas R_CHEBPI_P ---
470 //---------------------
471
472void som_r_chebpi_p
473 (double* ti, const int nr, const int nt, const int np, const double x, double* trtp) {
474
475// Variables diverses
476int i, j, k ;
477double* pi = ti ; // pointeur courant sur l'entree
478double* po = trtp ; // pointeur courant sur la sortie
479
480 // Valeurs des polynomes de Chebyshev au point x demande
481 double* cheb = new double [2*nr] ;
482 cheb[0] = 1. ;
483 cheb[1] = x ;
484 for (i=2 ; i<2*nr ; i++) {
485 cheb[i] = 2*x* cheb[i-1] - cheb[i-2] ;
486 }
487
488 // Sommation pour le premier phi, k=0
489 for (j=0 ; j<nt ; j++) {
490 int l = j%2;
491 if(l==0){
492 *po = 0 ;
493 // Sommation sur r
494 for (i=0 ; i<nr ; i++) {
495 *po += (*pi) * cheb[2*i] ;
496 pi++ ; // R suivant
497 }
498 po++ ; // Theta suivant
499 } else {
500 *po = 0 ;
501 // Sommation sur r
502 for (i=0 ; i<nr ; i++) {
503 *po += (*pi) * cheb[2*i+1] ;
504 pi++ ; // R suivant
505 }
506 po++ ; // Theta suivant
507 }
508 }
509
510 if (np > 1) {
511
512 // On saute le deuxieme phi (sin(0)), k=1
513 pi += nr*nt ;
514 po += nt ;
515
516 // Sommation sur les phi suivants (pour k=2,...,np)
517 for (k=2 ; k<np+1 ; k++) {
518 for (j=0 ; j<nt ; j++) {
519 int l = j% 2 ; // parite: 0 <-> 1
520 // Sommation sur r
521 *po = 0 ;
522 for (i=0 ; i<nr ; i++) {
523 *po += (*pi) * cheb[2*i + l] ;
524 pi++ ; // R suivant
525 }
526 po++ ; // Theta suivant
527 }
528 }
529
530 } // fin du cas np > 1
531 // Menage
532 delete [] cheb ;
533
534}
535
536 //----------------------
537 //- Cas R_CHEBPI_I ----
538 //----------------------
539
540void som_r_chebpi_i
541 (double* ti, const int nr, const int nt, const int np, const double x, double* trtp) {
542
543// Variables diverses
544int i, j, k ;
545double* pi = ti ; // pointeur courant sur l'entree
546double* po = trtp ; // pointeur courant sur la sortie
547
548 // Valeurs des polynomes de Chebyshev au point x demande
549 double* cheb = new double [2*nr] ;
550 cheb[0] = 1. ;
551 cheb[1] = x ;
552 for (i=2 ; i<2*nr ; i++) {
553 cheb[i] = 2*x* cheb[i-1] - cheb[i-2] ;
554 }
555
556 // Sommation pour le premier phi, k=0
557 for (j=0 ; j<nt ; j++) {
558 int l = j%2 ;
559 if(l==1){
560 *po = 0 ;
561 // Sommation sur r
562 for (i=0 ; i<nr ; i++) {
563 *po += (*pi) * cheb[2*i] ;
564 pi++ ; // R suivant
565 }
566 po++ ; // Theta suivant
567 } else {
568 *po = 0 ;
569 // Sommation sur r
570 for (i=0 ; i<nr ; i++) {
571 *po += (*pi) * cheb[2*i+1] ;
572 pi++ ; // R suivant
573 }
574 po++ ; // Theta suivant
575 }
576 }
577
578 if (np > 1) {
579
580 // On saute le deuxieme phi (sin(0)), k=1
581 pi += nr*nt ;
582 po += nt ;
583
584 // Sommation sur les phi suivants (pour k=2,...,np)
585 for (k=2 ; k<np+1 ; k++) {
586 for (j=0 ; j<nt ; j++) {
587 int l = j % 2 ; // parite: 0 <-> 1
588 // Sommation sur r
589 *po = 0 ;
590 for (i=0 ; i<nr ; i++) {
591 *po += (*pi) * cheb[2*i + 1 - l] ;
592 pi++ ; // R suivant
593 }
594 po++ ; // Theta suivant
595 }
596 }
597
598 } // fin du cas np > 1
599 // Menage
600 delete [] cheb ;
601
602}
603 //----------------
604 //- Cas R_LEG ---
605 //----------------
606
607void som_r_leg
608 (double* ti, const int nr, const int nt, const int np, const double x,
609 double* trtp) {
610
611// Variables diverses
612int i, j, k ;
613double* pi = ti ; // pointeur courant sur l'entree
614double* po = trtp ; // pointeur courant sur la sortie
615
616 // Valeurs des polynomes de Legendre au point x demande
617 double* leg = new double [nr] ;
618 leg[0] = 1. ;
619 if (nr > 1) leg[1] = x ;
620 for (i=2; i<nr; i++) {
621 leg[i] = (double(2*i-1)*x* leg[i-1] - double(i-1)*leg[i-2]) / double(i) ;
622 }
623
624 // Sommation pour le premier phi, k=0
625 for (j=0 ; j<nt ; j++) {
626 *po = 0 ;
627 // Sommation sur r
628 for (i=0 ; i<nr ; i++) {
629 *po += (*pi) * leg[i] ;
630 pi++ ; // R suivant
631 }
632 po++ ; // Theta suivant
633 }
634
635 if (np > 1) {
636
637 // On saute le deuxieme phi (sin(0)), k=1
638 pi += nr*nt ;
639 po += nt ;
640
641 // Sommation sur les phi suivants (pour k=2,...,np)
642 for (k=2 ; k<np+1 ; k++) {
643 for (j=0 ; j<nt ; j++) {
644 // Sommation sur r
645 *po = 0 ;
646 for (i=0 ; i<nr ; i++) {
647 *po += (*pi) * leg[i] ;
648 pi++ ; // R suivant
649 }
650 po++ ; // Theta suivant
651 }
652 }
653
654 } // fin du cas np > 1
655
656 // Menage
657 delete [] leg ;
658
659}
660
661
662 //-----------------
663 //- Cas R_LEGP ---
664 //-----------------
665
666void som_r_legp
667 (double* ti, const int nr, const int nt, const int np, const double x,
668 double* trtp) {
669
670// Variables diverses
671int i, j, k ;
672double* pi = ti ; // pointeur courant sur l'entree
673double* po = trtp ; // pointeur courant sur la sortie
674
675 // Valeurs des polynomes de Legendre au point x demande
676 double* leg = new double [nr] ;
677 leg[0] = 1. ;
678 double p2im1 = x ;
679 for (i=1; i<nr; i++) {
680 leg[i] = (double(4*i-1)*x* p2im1 - double(2*i-1)*leg[i-1]) / double(2*i) ;
681 p2im1 = (double(4*i+1)*x* leg[i] - double(2*i)*p2im1) / double(2*i+1) ;
682 }
683
684 // Sommation pour le premier phi, k=0
685 for (j=0 ; j<nt ; j++) {
686 *po = 0 ;
687 // Sommation sur r
688 for (i=0 ; i<nr ; i++) {
689 *po += (*pi) * leg[i] ;
690 pi++ ; // R suivant
691 }
692 po++ ; // Theta suivant
693 }
694
695 if (np > 1) {
696
697 // On saute le deuxieme phi (sin(0)), k=1
698 pi += nr*nt ;
699 po += nt ;
700
701 // Sommation sur les phi suivants (pour k=2,...,np)
702 for (k=2 ; k<np+1 ; k++) {
703 for (j=0 ; j<nt ; j++) {
704 // Sommation sur r
705 *po = 0 ;
706 for (i=0 ; i<nr ; i++) {
707 *po += (*pi) * leg[i] ;
708 pi++ ; // R suivant
709 }
710 po++ ; // Theta suivant
711 }
712 }
713
714 } // fin du cas np > 1
715 // Menage
716 delete [] leg ;
717
718}
719
720
721 //-----------------
722 //- Cas R_LEGI ---
723 //-----------------
724
725void som_r_legi
726 (double* ti, const int nr, const int nt, const int np, const double x,
727 double* trtp) {
728
729// Variables diverses
730int i, j, k ;
731double* pi = ti ; // pointeur courant sur l'entree
732double* po = trtp ; // pointeur courant sur la sortie
733
734 // Valeurs des polynomes de Legendre au point x demande
735 double* leg = new double [nr] ;
736 leg[0] = x ;
737 double p2im1 = 1. ;
738 for (i=1; i<nr; i++) {
739 p2im1 = (double(4*i-1)*x* leg[i-1] - double(2*i-1)*p2im1) / double(2*i) ;
740 leg[i] = (double(4*i+1)*x* p2im1 - double(2*i)*leg[i-1]) / double(2*i+1) ;
741 }
742
743 // Sommation pour le premier phi, k=0
744 for (j=0 ; j<nt ; j++) {
745 *po = 0 ;
746 // Sommation sur r
747 for (i=0 ; i<nr ; i++) {
748 *po += (*pi) * leg[i] ;
749 pi++ ; // R suivant
750 }
751 po++ ; // Theta suivant
752 }
753
754 if (np > 1) {
755
756 // On saute le deuxieme phi (sin(0)), k=1
757 pi += nr*nt ;
758 po += nt ;
759
760 // Sommation sur les phi suivants (pour k=2,...,np)
761 for (k=2 ; k<np+1 ; k++) {
762 for (j=0 ; j<nt ; j++) {
763 // Sommation sur r
764 *po = 0 ;
765 for (i=0 ; i<nr ; i++) {
766 *po += (*pi) * leg[i] ;
767 pi++ ; // R suivant
768 }
769 po++ ; // Theta suivant
770 }
771 }
772
773 } // fin du cas np > 1
774 // Menage
775 delete [] leg ;
776
777}
778
779 //----------------
780 //- Cas R_JACO02 -
781 //----------------
782
783void som_r_jaco02
784 (double* ti, const int nr, const int nt, const int np, const double x,
785 double* trtp) {
786
787// Variables diverses
788int i, j, k ;
789double* pi = ti ; // pointeur courant sur l'entree
790double* po = trtp ; // pointeur courant sur la sortie
791
792 // Valeurs des polynomes de Jacobi(0,2) au point x demande
793 double* jaco = jacobi(nr-1,x) ;
794
795 // Sommation pour le premier phi, k=0
796 for (j=0 ; j<nt ; j++) {
797 *po = 0 ;
798 // Sommation sur r
799 for (i=0 ; i<nr ; i++) {
800 *po += (*pi) * jaco[i] ;
801 pi++ ; // R suivant
802 }
803 po++ ; // Theta suivant
804 }
805
806 if (np > 1) {
807
808 // On saute le deuxieme phi (sin(0)), k=1
809 pi += nr*nt ;
810 po += nt ;
811
812 // Sommation sur les phi suivants (pour k=2,...,np)
813 for (k=2 ; k<np+1 ; k++) {
814 for (j=0 ; j<nt ; j++) {
815 // Sommation sur r
816 *po = 0 ;
817 for (i=0 ; i<nr ; i++) {
818 *po += (*pi) * jaco[i] ;
819 pi++ ; // R suivant
820 }
821 po++ ; // Theta suivant
822 }
823 }
824
825 } // fin du cas np > 1
826
827 // Menage
828 delete [] jaco ;
829
830}
831}
Lorene prototypes.
Definition app_hor.h:64