LORENE
op_d2sdtet2.C
1/*
2 * Copyright (c) 1999-2000 Jean-Alain Marck
3 * Copyright (c) 1999-2001 Eric Gourgoulhon
4 *
5 * This file is part of LORENE.
6 *
7 * LORENE is free software; you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation; either version 2 of the License, or
10 * (at your option) any later version.
11 *
12 * LORENE is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with LORENE; if not, write to the Free Software
19 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 *
21 */
22
23
24
25
26/*
27 * Ensemble des routines de base pour la derivation seconde par rapport a theta
28 * (Utilisation interne)
29 *
30 * void _d2sdtet2_XXXX(Tbl * t, int & b)
31 * t pointeur sur le Tbl d'entree-sortie
32 * b base spectrale
33 *
34 */
35
36/*
37 * $Id: op_d2sdtet2.C,v 1.6 2016/12/05 16:18:07 j_novak Exp $
38 * $Log: op_d2sdtet2.C,v $
39 * Revision 1.6 2016/12/05 16:18:07 j_novak
40 * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
41 *
42 * Revision 1.5 2014/10/13 08:53:24 j_novak
43 * Lorene classes and functions now belong to the namespace Lorene.
44 *
45 * Revision 1.4 2009/10/09 14:00:54 j_novak
46 * New bases T_cos and T_SIN.
47 *
48 * Revision 1.3 2006/03/10 12:45:38 j_novak
49 * Use of C++-style cast.
50 *
51 * Revision 1.2 2004/11/23 15:16:01 m_forot
52 *
53 * Added the bases for the cases without any equatorial symmetry
54 * (T_COSSIN_C, T_COSSIN_S, T_LEG, R_CHEBPI_P, R_CHEBPI_I).
55 *
56 * Revision 1.1.1.1 2001/11/20 15:19:29 e_gourgoulhon
57 * LORENE
58 *
59 * Revision 2.4 2000/10/04 11:50:20 eric
60 * Ajout d' abort() dans le cas non prevu.
61 *
62 * Revision 2.3 2000/02/24 16:40:42 eric
63 * Initialisation a zero du tableau xo avant le calcul.
64 *
65 * Revision 2.2 1999/11/15 16:37:44 eric
66 * Tbl::dim est desormais un Dim_tbl et non plus un Dim_tbl*.
67 *
68 * Revision 2.1 1999/03/01 15:06:00 eric
69 * *** empty log message ***
70 *
71 * Revision 2.0 1999/02/23 11:23:09 hyc
72 * *** empty log message ***
73 *
74 *
75 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/op_d2sdtet2.C,v 1.6 2016/12/05 16:18:07 j_novak Exp $
76 *
77 */
78
79// Headers Lorene
80#include "tbl.h"
81
82// Routine pour les cas non prevus
83//--------------------------------
84namespace Lorene {
85void _d2sdtet2_pas_prevu(Tbl* , int & b) {
86 cout << "Unknown theta basis in Mtbl_cf::d2sdt2() !" << endl ;
87 cout << " basis: " << hex << b << endl ;
88 abort() ;
89}
90
91// cas T_COS
92//----------
93void _d2sdtet2_t_cos(Tbl* tb, int &)
94{
95
96 // Peut-etre rien a faire ?
97 if (tb->get_etat() == ETATZERO) {
98 return ;
99 }
100
101 // Protection
102 assert(tb->get_etat() == ETATQCQ) ;
103
104 // Pour le confort
105 int nr = (tb->dim).dim[0] ; // Nombre
106 int nt = (tb->dim).dim[1] ; // de points
107 int np = (tb->dim).dim[2] ; // physiques REELS
108 np = np - 2 ; // Nombre de points physiques
109
110 // Variables statiques
111 static double* cx = 0 ;
112 static int nt_pre =0 ;
113
114 // Test sur nt pour initialisation eventuelle
115 if (nt > nt_pre) {
116 nt_pre = nt ;
117 cx = reinterpret_cast<double*>(realloc(cx, nt * sizeof(double))) ;
118 for (int i=0 ; i<nt ; i++) {
119 cx[i] = - double(i * i) ;
120 }
121 }
122
123 // pt. sur le tableau de double resultat
124 double* xo = new double[(tb->dim).taille] ;
125
126 // Initialisation a zero :
127 for (int i=0; i<(tb->dim).taille; i++) {
128 xo[i] = 0 ;
129 }
130
131 // On y va...
132 double* xi = tb->t ;
133 double* xci = xi ; // Pointeurs
134 double* xco = xo ; // courants
135
136 // k = 0
137 for (int j=0 ; j<nt ; j++) {
138 for (int i=0 ; i<nr ; i++ ) {
139 *xco = cx[j] * (*xci) ;
140 xci++ ;
141 xco++ ;
142 } // Fin de la boucle sur r
143 } // Fin de la boucle sur theta
144
145 // k = 1
146 xci += nr*nt ;
147 xco += nr*nt ;
148
149 // k >= 2
150 int borne_phi = np + 1 ;
151 if (np == 1) borne_phi = 1 ;
152
153 for (int k=2 ; k<borne_phi ; k++) {
154 for (int j=0 ; j<nt ; j++) {
155 for (int i=0 ; i<nr ; i++ ) {
156 *xco = cx[j] * (*xci) ;
157 xci++ ;
158 xco++ ;
159 } // Fin de la boucle sur r
160 } // Fin de la boucle sur theta
161 } // Fin de la boucle sur phi
162
163 // On remet les choses la ou il faut
164 delete [] tb->t ;
165 tb->t = xo ;
166
167 // base de developpement
168 // inchangee
169}
170
171// cas T_SIN
172//----------
173void _d2sdtet2_t_sin(Tbl* tb, int &)
174{
175
176 // Peut-etre rien a faire ?
177 if (tb->get_etat() == ETATZERO) {
178 return ;
179 }
180
181 // Protection
182 assert(tb->get_etat() == ETATQCQ) ;
183
184 // Pour le confort
185 int nr = (tb->dim).dim[0] ; // Nombre
186 int nt = (tb->dim).dim[1] ; // de points
187 int np = (tb->dim).dim[2] ; // physiques REELS
188 np = np - 2 ; // Nombre de points physiques
189
190 // Variables statiques
191 static double* cx = 0 ;
192 static int nt_pre =0 ;
193
194 // Test sur nt pour initialisation eventuelle
195 if (nt > nt_pre) {
196 nt_pre = nt ;
197 cx = reinterpret_cast<double*>(realloc(cx, nt * sizeof(double))) ;
198 for (int i=0 ; i<nt ; i++) {
199 cx[i] = - double(i * i) ;
200 }
201 }
202
203 // pt. sur le tableau de double resultat
204 double* xo = new double[(tb->dim).taille] ;
205
206 // Initialisation a zero :
207 for (int i=0; i<(tb->dim).taille; i++) {
208 xo[i] = 0 ;
209 }
210
211 // On y va...
212 double* xi = tb->t ;
213 double* xci = xi ; // Pointeurs
214 double* xco = xo ; // courants
215
216 int borne_phi = np + 1 ;
217 if (np == 1) borne_phi = 1 ;
218
219 for (int k=0 ; k< borne_phi ; k++) {
220 for (int j=0 ; j<nt ; j++) {
221 for (int i=0 ; i<nr ; i++ ) {
222 *xco = cx[j] * (*xci) ;
223 xci++ ;
224 xco++ ;
225 } // Fin de la boucle sur r
226 } // Fin de la boucle sur theta
227 } // Fin de la boucle sur phi
228
229 // On remet les choses la ou il faut
230 delete [] tb->t ;
231 tb->t = xo ;
232
233 // base de developpement
234 // inchangee
235}
236
237// cas T_COS_P
238//------------
239void _d2sdtet2_t_cos_p(Tbl* tb, int &)
240{
241
242 // Peut-etre rien a faire ?
243 if (tb->get_etat() == ETATZERO) {
244 return ;
245 }
246
247 // Protection
248 assert(tb->get_etat() == ETATQCQ) ;
249
250 // Pour le confort
251 int nr = (tb->dim).dim[0] ; // Nombre
252 int nt = (tb->dim).dim[1] ; // de points
253 int np = (tb->dim).dim[2] ; // physiques REELS
254 np = np - 2 ; // Nombre de points physiques
255
256 // Variables statiques
257 static double* cx = 0 ;
258 static int nt_pre =0 ;
259
260 // Test sur nt pour initialisation eventuelle
261 if (nt > nt_pre) {
262 nt_pre = nt ;
263 cx = reinterpret_cast<double*>(realloc(cx, nt * sizeof(double))) ;
264 for (int i=0 ; i<nt ; i++) {
265 cx[i] = - (2*i) * (2*i) ;
266 }
267 }
268
269 // pt. sur le tableau de double resultat
270 double* xo = new double[(tb->dim).taille] ;
271
272 // Initialisation a zero :
273 for (int i=0; i<(tb->dim).taille; i++) {
274 xo[i] = 0 ;
275 }
276
277 // On y va...
278 double* xi = tb->t ;
279 double* xci = xi ; // Pointeurs
280 double* xco = xo ; // courants
281
282 // k = 0
283 for (int j=0 ; j<nt ; j++) {
284 for (int i=0 ; i<nr ; i++ ) {
285 *xco = cx[j] * (*xci) ;
286 xci++ ;
287 xco++ ;
288 } // Fin de la boucle sur r
289 } // Fin de la boucle sur theta
290
291 // k = 1
292 xci += nr*nt ;
293 xco += nr*nt ;
294
295 // k >= 2
296 int borne_phi = np + 1 ;
297 if (np == 1) borne_phi = 1 ;
298
299 for (int k=2 ; k<borne_phi ; k++) {
300 for (int j=0 ; j<nt ; j++) {
301 for (int i=0 ; i<nr ; i++ ) {
302 *xco = cx[j] * (*xci) ;
303 xci++ ;
304 xco++ ;
305 } // Fin de la boucle sur r
306 } // Fin de la boucle sur theta
307 } // Fin de la boucle sur phi
308
309 // On remet les choses la ou il faut
310 delete [] tb->t ;
311 tb->t = xo ;
312
313 // base de developpement
314 // inchangee
315}
316
317// cas T_SIN_P
318//------------
319void _d2sdtet2_t_sin_p(Tbl* tb, int &)
320{
321
322 // Peut-etre rien a faire ?
323 if (tb->get_etat() == ETATZERO) {
324 return ;
325 }
326
327 // Protection
328 assert(tb->get_etat() == ETATQCQ) ;
329
330 // Pour le confort
331 int nr = (tb->dim).dim[0] ; // Nombre
332 int nt = (tb->dim).dim[1] ; // de points
333 int np = (tb->dim).dim[2] ; // physiques REELS
334 np = np - 2 ; // Nombre de points physiques
335
336 // Variables statiques
337 static double* cx = 0 ;
338 static int nt_pre =0 ;
339
340 // Test sur nt pour initialisation eventuelle
341 if (nt > nt_pre) {
342 nt_pre = nt ;
343 cx = reinterpret_cast<double*>(realloc(cx, nt * sizeof(double))) ;
344 for (int i=0 ; i<nt ; i++) {
345 cx[i] = - (2*i) * (2*i) ;
346 }
347 }
348
349 // pt. sur le tableau de double resultat
350 double* xo = new double[(tb->dim).taille] ;
351
352 // Initialisation a zero :
353 for (int i=0; i<(tb->dim).taille; i++) {
354 xo[i] = 0 ;
355 }
356
357 // On y va...
358 double* xi = tb->t ;
359 double* xci = xi ; // Pointeurs
360 double* xco = xo ; // courants
361
362 int borne_phi = np + 1 ;
363 if (np == 1) borne_phi = 1 ;
364
365 for (int k=0 ; k< borne_phi ; k++) {
366 for (int j=0 ; j<nt ; j++) {
367 for (int i=0 ; i<nr ; i++ ) {
368 *xco = cx[j] * (*xci) ;
369 xci++ ;
370 xco++ ;
371 } // Fin de la boucle sur r
372 } // Fin de la boucle sur theta
373 } // Fin de la boucle sur phi
374
375 // On remet les choses la ou il faut
376 delete [] tb->t ;
377 tb->t = xo ;
378
379 // base de developpement
380 // inchangee
381}
382
383// cas T_SIN_I
384//------------
385void _d2sdtet2_t_sin_i(Tbl* tb, int &)
386{
387
388 // Peut-etre rien a faire ?
389 if (tb->get_etat() == ETATZERO) {
390 return ;
391 }
392
393 // Protection
394 assert(tb->get_etat() == ETATQCQ) ;
395
396 // Pour le confort
397 int nr = (tb->dim).dim[0] ; // Nombre
398 int nt = (tb->dim).dim[1] ; // de points
399 int np = (tb->dim).dim[2] ; // physiques REELS
400 np = np - 2 ; // Nombre de points physiques
401
402 // Variables statiques
403 static double* cx = 0 ;
404 static int nt_pre =0 ;
405
406 // Test sur nt pour initialisation eventuelle
407 if (nt > nt_pre) {
408 nt_pre = nt ;
409 cx = reinterpret_cast<double*>(realloc(cx, nt * sizeof(double))) ;
410 for (int i=0 ; i<nt ; i++) {
411 cx[i] = - (2*i+1) * (2*i+1) ;
412 }
413 }
414
415 // pt. sur le tableau de double resultat
416 double* xo = new double[(tb->dim).taille] ;
417
418 // Initialisation a zero :
419 for (int i=0; i<(tb->dim).taille; i++) {
420 xo[i] = 0 ;
421 }
422
423 // On y va...
424 double* xi = tb->t ;
425 double* xci = xi ; // Pointeurs
426 double* xco = xo ; // courants
427
428 int borne_phi = np + 1 ;
429 if (np == 1) borne_phi = 1 ;
430
431 for (int k=0 ; k< borne_phi ; k++) {
432 for (int j=0 ; j<nt ; j++) {
433 for (int i=0 ; i<nr ; i++ ) {
434 *xco = cx[j] * (*xci) ;
435 xci++ ;
436 xco++ ;
437 } // Fin de la boucle sur r
438 } // Fin de la boucle sur theta
439 } // Fin de la boucle sur phi
440
441 // On remet les choses la ou il faut
442 delete [] tb->t ;
443 tb->t = xo ;
444
445 // base de developpement
446 // inchangee
447}
448
449// cas T_COS_I
450//------------
451void _d2sdtet2_t_cos_i(Tbl* tb, int &)
452{
453
454 // Peut-etre rien a faire ?
455 if (tb->get_etat() == ETATZERO) {
456 return ;
457 }
458
459 // Protection
460 assert(tb->get_etat() == ETATQCQ) ;
461
462 // Pour le confort
463 int nr = (tb->dim).dim[0] ; // Nombre
464 int nt = (tb->dim).dim[1] ; // de points
465 int np = (tb->dim).dim[2] ; // physiques REELS
466 np = np - 2 ; // Nombre de points physiques
467
468 // Variables statiques
469 static double* cx = 0 ;
470 static int nt_pre =0 ;
471
472 // Test sur nt pour initialisation eventuelle
473 if (nt > nt_pre) {
474 nt_pre = nt ;
475 cx = reinterpret_cast<double*>(realloc(cx, nt * sizeof(double))) ;
476 for (int i=0 ; i<nt ; i++) {
477 cx[i] = - (2*i+1) * (2*i+1) ;
478 }
479 }
480
481 // pt. sur le tableau de double resultat
482 double* xo = new double[(tb->dim).taille] ;
483
484 // Initialisation a zero :
485 for (int i=0; i<(tb->dim).taille; i++) {
486 xo[i] = 0 ;
487 }
488
489 // On y va...
490 double* xi = tb->t ;
491 double* xci = xi ; // Pointeurs
492 double* xco = xo ; // courants
493
494 int borne_phi = np + 1 ;
495 if (np == 1) borne_phi = 1 ;
496
497 for (int k=0 ; k< borne_phi ; k++) {
498 for (int j=0 ; j<nt ; j++) {
499 for (int i=0 ; i<nr ; i++ ) {
500 *xco = cx[j] * (*xci) ;
501 xci++ ;
502 xco++ ;
503 } // Fin de la boucle sur r
504 } // Fin de la boucle sur theta
505 } // Fin de la boucle sur phi
506
507 // On remet les choses la ou il faut
508 delete [] tb->t ;
509 tb->t = xo ;
510
511 // base de developpement
512 // inchangee
513}
514
515// cas T_COSSIN_CP
516//----------------
517void _d2sdtet2_t_cossin_cp(Tbl* tb, int &)
518{
519
520 // Peut-etre rien a faire ?
521 if (tb->get_etat() == ETATZERO) {
522 return ;
523 }
524
525 // Protection
526 assert(tb->get_etat() == ETATQCQ) ;
527
528 // Pour le confort
529 int nr = (tb->dim).dim[0] ; // Nombre
530 int nt = (tb->dim).dim[1] ; // de points
531 int np = (tb->dim).dim[2] ; // physiques REELS
532 np = np - 2 ; // Nombre de points physiques
533
534 // Variables statiques
535 static double* cxp = 0 ;
536 static double* cxi = 0 ;
537 static int nt_pre =0 ;
538
539 // Test sur nt pour initialisation eventuelle
540 if (nt > nt_pre) {
541 nt_pre = nt ;
542 cxp = reinterpret_cast<double*>(realloc(cxp, nt * sizeof(double))) ;
543 cxi = reinterpret_cast<double*>(realloc(cxi, nt * sizeof(double))) ;
544 for (int i=0 ; i<nt ; i++) {
545 cxp[i] = - (2*i) * (2*i) ;
546 cxi[i] = - (2*i+1) * (2*i+1) ;
547 }
548 }
549
550 // pt. sur le tableau de double resultat
551 double* xo = new double[(tb->dim).taille] ;
552
553 // Initialisation a zero :
554 for (int i=0; i<(tb->dim).taille; i++) {
555 xo[i] = 0 ;
556 }
557
558 // On y va...
559 double* xi = tb->t ;
560 double* xci = xi ; // Pointeurs
561 double* xco = xo ; // courants
562
563 // Partie cos(pair)
564 int k ;
565 for (k=0 ; k<np+1 ; k += 4) {
566 for (int m=0 ; m<2 ; m++) {
567 for (int j=0 ; j<nt ; j++) {
568 for (int i=0 ; i<nr ; i++ ) {
569 *xco = cxp[j] * (*xci) ;
570 xci++ ;
571 xco++ ;
572 } // Fin de la boucle sur r
573 } // Fin de la boucle sur theta
574 } // Fin de la boucle intermediaire
575 xci += 2*nr*nt ;
576 xco += 2*nr*nt ;
577 } // Fin de la boucle sur phi
578
579 // Partie sin(impair)
580 xci = xi + 2*nr*nt ;
581 xco = xo + 2*nr*nt ;
582 for (k=2 ; k<np+1 ; k += 4) {
583 for (int m=0 ; m<2 ; m++) {
584 for (int j=0 ; j<nt ; j++) {
585 for (int i=0 ; i<nr ; i++ ) {
586 *xco = cxi[j] * (*xci) ;
587 xci++ ;
588 xco++ ;
589 } // Fin de la boucle sur r
590 } // Fin de la boucle sur theta
591 } // Fin de la boucle intermediaire
592 xci += 2*nr*nt ;
593 xco += 2*nr*nt ;
594 } // Fin de la boucle sur phi
595
596 // On remet les choses la ou il faut
597 delete [] tb->t ;
598 tb->t = xo ;
599
600 // base de developpement
601 // inchangee
602}
603
604// cas T_COSSIN_SP
605//----------------
606void _d2sdtet2_t_cossin_sp(Tbl* tb, int &)
607{
608
609 // Peut-etre rien a faire ?
610 if (tb->get_etat() == ETATZERO) {
611 return ;
612 }
613
614 // Protection
615 assert(tb->get_etat() == ETATQCQ) ;
616
617 // Pour le confort
618 int nr = (tb->dim).dim[0] ; // Nombre
619 int nt = (tb->dim).dim[1] ; // de points
620 int np = (tb->dim).dim[2] ; // physiques REELS
621 np = np - 2 ; // Nombre de points physiques
622
623 // Variables statiques
624 static double* cxp = 0 ;
625 static double* cxi = 0 ;
626 static int nt_pre =0 ;
627
628 // Test sur nt pour initialisation eventuelle
629 if (nt > nt_pre) {
630 nt_pre = nt ;
631 cxp = reinterpret_cast<double*>(realloc(cxp, nt * sizeof(double))) ;
632 cxi = reinterpret_cast<double*>(realloc(cxi, nt * sizeof(double))) ;
633 for (int i=0 ; i<nt ; i++) {
634 cxp[i] = - (2*i) * (2*i) ;
635 cxi[i] = - (2*i+1) * (2*i+1) ;
636 }
637 }
638
639 // pt. sur le tableau de double resultat
640 double* xo = new double[(tb->dim).taille] ;
641
642 // Initialisation a zero :
643 for (int i=0; i<(tb->dim).taille; i++) {
644 xo[i] = 0 ;
645 }
646
647 // On y va...
648 double* xi = tb->t ;
649 double* xci = xi ; // Pointeurs
650 double* xco = xo ; // courants
651
652 // Partie sin(pair)
653 int k ;
654 for (k=0 ; k<np+1 ; k += 4) {
655 for (int m=0 ; m<2 ; m++) {
656 for (int j=0 ; j<nt ; j++) {
657 for (int i=0 ; i<nr ; i++ ) {
658 *xco = cxp[j] * (*xci) ;
659 xci++ ;
660 xco++ ;
661 } // Fin de la boucle sur r
662 } // Fin de la boucle sur theta
663 } // Fin de la boucle intermediaire
664 xci += 2*nr*nt ;
665 xco += 2*nr*nt ;
666 } // Fin de la boucle sur phi
667
668 // Partie cos(impair)
669 xci = xi + 2*nr*nt ;
670 xco = xo + 2*nr*nt ;
671 for (k=2 ; k<np+1 ; k += 4) {
672 for (int m=0 ; m<2 ; m++) {
673 for (int j=0 ; j<nt ; j++) {
674 for (int i=0 ; i<nr ; i++ ) {
675 *xco = cxi[j] * (*xci) ;
676 xci++ ;
677 xco++ ;
678 } // Fin de la boucle sur r
679 } // Fin de la boucle sur theta
680 } // Fin de la boucle intermediaire
681 xci += 2*nr*nt ;
682 xco += 2*nr*nt ;
683 } // Fin de la boucle sur phi
684
685 // On remet les choses la ou il faut
686 delete [] tb->t ;
687 tb->t = xo ;
688
689 // base de developpement
690 // inchangee
691}
692
693// cas T_COSSIN_SI
694//----------------
695void _d2sdtet2_t_cossin_si(Tbl* tb, int &)
696{
697
698 // Peut-etre rien a faire ?
699 if (tb->get_etat() == ETATZERO) {
700 return ;
701 }
702
703 // Protection
704 assert(tb->get_etat() == ETATQCQ) ;
705
706 // Pour le confort
707 int nr = (tb->dim).dim[0] ; // Nombre
708 int nt = (tb->dim).dim[1] ; // de points
709 int np = (tb->dim).dim[2] ; // physiques REELS
710 np = np - 2 ; // Nombre de points physiques
711
712 // Variables statiques
713 static double* cxp = 0 ;
714 static double* cxi = 0 ;
715 static int nt_pre =0 ;
716
717 // Test sur nt pour initialisation eventuelle
718 if (nt > nt_pre) {
719 nt_pre = nt ;
720 cxp = reinterpret_cast<double*>(realloc(cxp, nt * sizeof(double))) ;
721 cxi = reinterpret_cast<double*>(realloc(cxi, nt * sizeof(double))) ;
722 for (int i=0 ; i<nt ; i++) {
723 cxp[i] = - (2*i) * (2*i) ;
724 cxi[i] = - (2*i+1) * (2*i+1) ;
725 }
726 }
727
728 // pt. sur le tableau de double resultat
729 double* xo = new double[(tb->dim).taille] ;
730
731 // Initialisation a zero :
732 for (int i=0; i<(tb->dim).taille; i++) {
733 xo[i] = 0 ;
734 }
735
736 // On y va...
737 double* xi = tb->t ;
738 double* xci = xi ; // Pointeurs
739 double* xco = xo ; // courants
740
741 // Partie sin(impair)
742 int k ;
743 for (k=0 ; k<np+1 ; k += 4) {
744 for (int m=0 ; m<2 ; m++) {
745 for (int j=0 ; j<nt ; j++) {
746 for (int i=0 ; i<nr ; i++ ) {
747 *xco = cxi[j] * (*xci) ;
748 xci++ ;
749 xco++ ;
750 } // Fin de la boucle sur r
751 } // Fin de la boucle sur theta
752 } // Fin de la boucle intermediaire
753 xci += 2*nr*nt ;
754 xco += 2*nr*nt ;
755 } // Fin de la boucle sur phi
756
757 // Partie cos(pair)
758 xci = xi + 2*nr*nt ;
759 xco = xo + 2*nr*nt ;
760 for (k=2 ; k<np+1 ; k += 4) {
761 for (int m=0 ; m<2 ; m++) {
762 for (int j=0 ; j<nt ; j++) {
763 for (int i=0 ; i<nr ; i++ ) {
764 *xco = cxp[j] * (*xci) ;
765 xci++ ;
766 xco++ ;
767 } // Fin de la boucle sur r
768 } // Fin de la boucle sur theta
769 } // Fin de la boucle intermediaire
770 xci += 2*nr*nt ;
771 xco += 2*nr*nt ;
772 } // Fin de la boucle sur phi
773
774 // On remet les choses la ou il faut
775 delete [] tb->t ;
776 tb->t = xo ;
777
778 // base de developpement
779 // inchangee
780}
781
782// cas T_COSSIN_C
783//----------------
784void _d2sdtet2_t_cossin_c(Tbl* tb, int &)
785{
786
787 // Peut-etre rien a faire ?
788 if (tb->get_etat() == ETATZERO) {
789 return ;
790 }
791
792 // Protection
793 assert(tb->get_etat() == ETATQCQ) ;
794
795 // Pour le confort
796 int nr = (tb->dim).dim[0] ; // Nombre
797 int nt = (tb->dim).dim[1] ; // de points
798 int np = (tb->dim).dim[2] ; // physiques REELS
799 np = np - 2 ; // Nombre de points physiques
800
801 // Variables statiques
802 static double* cxp = 0 ;
803 static double* cxi = 0 ;
804 static int nt_pre =0 ;
805
806 // Test sur nt pour initialisation eventuelle
807 if (nt > nt_pre) {
808 nt_pre = nt ;
809 cxp = reinterpret_cast<double*>(realloc(cxp, nt * sizeof(double))) ;
810 cxi = reinterpret_cast<double*>(realloc(cxi, nt * sizeof(double))) ;
811 for (int i=0 ; i<nt ; i++) {
812 cxp[i] = - i*i ;
813 cxi[i] = - i*i ;
814 }
815 }
816
817 // pt. sur le tableau de double resultat
818 double* xo = new double[(tb->dim).taille] ;
819
820 // Initialisation a zero :
821 for (int i=0; i<(tb->dim).taille; i++) {
822 xo[i] = 0 ;
823 }
824
825 // On y va...
826 double* xi = tb->t ;
827 double* xci = xi ; // Pointeurs
828 double* xco = xo ; // courants
829
830 // Partie cos(pair)
831 int k ;
832 for (k=0 ; k<np+1 ; k += 4) {
833 for (int m=0 ; m<2 ; m++) {
834 for (int j=0 ; j<nt ; j++) {
835 for (int i=0 ; i<nr ; i++ ) {
836 *xco = cxp[j] * (*xci) ;
837 xci++ ;
838 xco++ ;
839 } // Fin de la boucle sur r
840 } // Fin de la boucle sur theta
841 } // Fin de la boucle intermediaire
842 xci += 2*nr*nt ;
843 xco += 2*nr*nt ;
844 } // Fin de la boucle sur phi
845
846 // Partie sin(impair)
847 xci = xi + 2*nr*nt ;
848 xco = xo + 2*nr*nt ;
849 for (k=2 ; k<np+1 ; k += 4) {
850 for (int m=0 ; m<2 ; m++) {
851 for (int j=0 ; j<nt ; j++) {
852 for (int i=0 ; i<nr ; i++ ) {
853 *xco = cxi[j] * (*xci) ;
854 xci++ ;
855 xco++ ;
856 } // Fin de la boucle sur r
857 } // Fin de la boucle sur theta
858 } // Fin de la boucle intermediaire
859 xci += 2*nr*nt ;
860 xco += 2*nr*nt ;
861 } // Fin de la boucle sur phi
862
863 // On remet les choses la ou il faut
864 delete [] tb->t ;
865 tb->t = xo ;
866
867 // base de developpement
868 // inchangee
869}
870
871// cas T_COSSIN_S
872//----------------
873void _d2sdtet2_t_cossin_s(Tbl* tb, int &)
874{
875
876 // Peut-etre rien a faire ?
877 if (tb->get_etat() == ETATZERO) {
878 return ;
879 }
880
881 // Protection
882 assert(tb->get_etat() == ETATQCQ) ;
883
884 // Pour le confort
885 int nr = (tb->dim).dim[0] ; // Nombre
886 int nt = (tb->dim).dim[1] ; // de points
887 int np = (tb->dim).dim[2] ; // physiques REELS
888 np = np - 2 ; // Nombre de points physiques
889
890 // Variables statiques
891 static double* cxp = 0 ;
892 static double* cxi = 0 ;
893 static int nt_pre =0 ;
894
895 // Test sur nt pour initialisation eventuelle
896 if (nt > nt_pre) {
897 nt_pre = nt ;
898 cxp = reinterpret_cast<double*>(realloc(cxp, nt * sizeof(double))) ;
899 cxi = reinterpret_cast<double*>(realloc(cxi, nt * sizeof(double))) ;
900 for (int i=0 ; i<nt ; i++) {
901 cxp[i] = - i*i ;
902 cxi[i] = - i*i ;
903 }
904 }
905
906 // pt. sur le tableau de double resultat
907 double* xo = new double[(tb->dim).taille] ;
908
909 // Initialisation a zero :
910 for (int i=0; i<(tb->dim).taille; i++) {
911 xo[i] = 0 ;
912 }
913
914 // On y va...
915 double* xi = tb->t ;
916 double* xci = xi ; // Pointeurs
917 double* xco = xo ; // courants
918
919 // Partie sin(pair)
920 int k ;
921 for (k=0 ; k<np+1 ; k += 4) {
922 for (int m=0 ; m<2 ; m++) {
923 for (int j=0 ; j<nt ; j++) {
924 for (int i=0 ; i<nr ; i++ ) {
925 *xco = cxp[j] * (*xci) ;
926 xci++ ;
927 xco++ ;
928 } // Fin de la boucle sur r
929 } // Fin de la boucle sur theta
930 } // Fin de la boucle intermediaire
931 xci += 2*nr*nt ;
932 xco += 2*nr*nt ;
933 } // Fin de la boucle sur phi
934
935 // Partie cos(impair)
936 xci = xi + 2*nr*nt ;
937 xco = xo + 2*nr*nt ;
938 for (k=2 ; k<np+1 ; k += 4) {
939 for (int m=0 ; m<2 ; m++) {
940 for (int j=0 ; j<nt ; j++) {
941 for (int i=0 ; i<nr ; i++ ) {
942 *xco = cxi[j] * (*xci) ;
943 xci++ ;
944 xco++ ;
945 } // Fin de la boucle sur r
946 } // Fin de la boucle sur theta
947 } // Fin de la boucle intermediaire
948 xci += 2*nr*nt ;
949 xco += 2*nr*nt ;
950 } // Fin de la boucle sur phi
951
952 // On remet les choses la ou il faut
953 delete [] tb->t ;
954 tb->t = xo ;
955
956 // base de developpement
957 // inchangee
958}
959}
Lorene::Tbl
Basic array class.
Definition tbl.h:161
Lorene
Lorene prototypes.
Definition app_hor.h:67