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