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 char 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 
110 namespace Lorene {
111  //-------------------
112  //- Cas Non-Prevu ---
113  //-------------------
114 
115 void 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 
126 void 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
131 int i, j, k ;
132 double* pi = ti ; // pointeur courant sur l'entree
133 double* 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 
185 void 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
190 int i, j, k ;
191 double* pi = ti ; // pointeur courant sur l'entree
192 double* 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 
244 void 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
249 int i, j, k ;
250 double* pi = ti ; // pointeur courant sur l'entree
251 double* 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 
301 void som_r_chebu
302  (double* ti, const int nr, const int nt, const int np, const double x, double* trtp) {
303 
304 // Variables diverses
305 int i, j, k ;
306 double* pi = ti ; // pointeur courant sur l'entree
307 double* 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 
356 void 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
360 int i, j, k ;
361 double* pi = ti ; // pointeur courant sur l'entree
362 double* 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 
414 void 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
418 int i, j, k ;
419 double* pi = ti ; // pointeur courant sur l'entree
420 double* 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 
472 void 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
476 int i, j, k ;
477 double* pi = ti ; // pointeur courant sur l'entree
478 double* 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 
540 void 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
544 int i, j, k ;
545 double* pi = ti ; // pointeur courant sur l'entree
546 double* 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 
607 void 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
612 int i, j, k ;
613 double* pi = ti ; // pointeur courant sur l'entree
614 double* 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 
666 void 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
671 int i, j, k ;
672 double* pi = ti ; // pointeur courant sur l'entree
673 double* 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 
725 void 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
730 int i, j, k ;
731 double* pi = ti ; // pointeur courant sur l'entree
732 double* 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 
783 void 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
788 int i, j, k ;
789 double* pi = ti ; // pointeur courant sur l'entree
790 double* 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