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 
113 namespace Lorene {
114  //-------------------
115  //- Cas Non-Prevu ---
116  //-------------------
117 
118 void 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 
129 void 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
134 int i, j, k ;
135 double* pi = ti ; // pointeur courant sur l'entree
136 double* 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 
188 void 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
193 int i, j, k ;
194 double* pi = ti ; // pointeur courant sur l'entree
195 double* 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 
247 void 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
252 int i, j, k ;
253 double* pi = ti ; // pointeur courant sur l'entree
254 double* 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 
304 void som_r_chebu
305  (double* ti, const int nr, const int nt, const int np, const double x, double* trtp) {
306 
307 // Variables diverses
308 int i, j, k ;
309 double* pi = ti ; // pointeur courant sur l'entree
310 double* 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 
359 void 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
363 int i, j, k ;
364 double* pi = ti ; // pointeur courant sur l'entree
365 double* 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 
417 void 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
421 int i, j, k ;
422 double* pi = ti ; // pointeur courant sur l'entree
423 double* 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 
475 void 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
479 int i, j, k ;
480 double* pi = ti ; // pointeur courant sur l'entree
481 double* 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 
543 void 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
547 int i, j, k ;
548 double* pi = ti ; // pointeur courant sur l'entree
549 double* 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 
610 void 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
615 int i, j, k ;
616 double* pi = ti ; // pointeur courant sur l'entree
617 double* 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 
669 void 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
674 int i, j, k ;
675 double* pi = ti ; // pointeur courant sur l'entree
676 double* 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 
728 void 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
733 int i, j, k ;
734 double* pi = ti ; // pointeur courant sur l'entree
735 double* 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 
786 void 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
791 int i, j, k ;
792 double* pi = ti ; // pointeur courant sur l'entree
793 double* 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