LORENE
op_primr.C
1/*
2 * Computation of primitive in a single domain
3 *
4 *
5 */
6
7/*
8 * Copyright (c) 2004 Eric Gourgoulhon.
9 *
10 * This file is part of LORENE.
11 *
12 * LORENE is free software; you can redistribute it and/or modify
13 * it under the terms of the GNU General Public License version 2
14 * as published by the Free Software Foundation.
15 *
16 * LORENE is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 * GNU General Public License for more details.
20 *
21 * You should have received a copy of the GNU General Public License
22 * along with LORENE; if not, write to the Free Software
23 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
24 *
25 */
26
27
28
29/*
30 * $Id: op_primr.C,v 1.13 2017/02/24 14:55:55 j_novak Exp $
31 * $Log: op_primr.C,v $
32 * Revision 1.13 2017/02/24 14:55:55 j_novak
33 * Corrected error in the primitive formula
34 *
35 * Revision 1.12 2017/02/22 17:11:33 j_novak
36 * Addition of new Legendre basis.
37 *
38 * Revision 1.11 2016/12/05 16:18:08 j_novak
39 * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
40 *
41 * Revision 1.10 2014/10/13 08:53:26 j_novak
42 * Lorene classes and functions now belong to the namespace Lorene.
43 *
44 * Revision 1.9 2014/10/06 15:16:06 j_novak
45 * Modified #include directives to use c++ syntax.
46 *
47 * Revision 1.8 2013/04/25 15:46:06 j_novak
48 * Added special treatment in the case np = 1, for type_p = NONSYM.
49 *
50 * Revision 1.7 2007/12/21 13:59:02 j_novak
51 * Suppression of call to pow(-1, something).
52 *
53 * Revision 1.6 2007/12/11 15:28:18 jl_cornou
54 * Jacobi(0,2) polynomials partially implemented
55 *
56 * Revision 1.5 2006/05/17 15:01:16 j_novak
57 * Treatment of the case nr = 1 and R_CHEB
58 *
59 * Revision 1.4 2004/11/23 15:16:01 m_forot
60 *
61 * Added the bases for the cases without any equatorial symmetry
62 * (T_COSSIN_C, T_COSSIN_S, T_LEG, R_CHEBPI_P, R_CHEBPI_I).
63 *
64 * Revision 1.3 2004/10/12 09:58:24 j_novak
65 * Better memory management.
66 *
67 * Revision 1.2 2004/06/14 15:24:57 e_gourgoulhon
68 * First operationnal version (tested).
69 *
70 * Revision 1.1 2004/06/13 21:33:13 e_gourgoulhon
71 * First version.
72 *
73 *
74 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/op_primr.C,v 1.13 2017/02/24 14:55:55 j_novak Exp $
75 *
76 */
77
78
79// C headers
80#include <cstdlib>
81#include <cmath>
82
83// Lorene headers
84#include "tbl.h"
85
86// Unexpected case
87//----------------
88namespace Lorene {
89void _primr_pas_prevu(const Tbl&, int bin, const Tbl&, Tbl&, int&, Tbl& ) {
90
91 cout << "Unexpected basis in primr : basis = " << hex << bin << endl ;
92 abort() ;
93
94}
95
96// case R_CHEB
97//------------
98void _primr_r_cheb(const Tbl& tin, int bin, const Tbl& valm1, Tbl& tout,
99 int& bout, Tbl& valp1) {
100
101 assert(tin.dim == tout.dim) ;
102
103 // Output spectral basis
104 bout = bin ;
105
106 // Number of coefficients
107 int nr = tin.get_dim(0) ;
108 int nt = tin.get_dim(1) ;
109 int np = tin.get_dim(2) - 2 ;
110 int borne_phi = np + 1 ;
111 if (np == 1) borne_phi = 1 ;
112
113 // Case of a zero input or pure angular grid
114 // -----------------------------------------
115 if ((tin.get_etat() == ETATZERO)||(nr == 1)) {
116 if (valm1.get_etat() == ETATZERO) {
117 tout.set_etat_zero() ;
118 valp1.set_etat_zero() ;
119 return ;
120 }
121 else {
122 assert(valm1.get_etat() == ETATQCQ) ;
123 tout.set_etat_qcq() ;
124 valp1.set_etat_qcq() ;
125 double* xco = tout.t ;
126 for (int k=0 ; k< borne_phi ; k++) {
127 if (k==1) { // jump over the coefficient of sin(0*phi)
128 xco += nr*nt ;
129 }
130 else {
131 for (int j=0 ; j<nt ; j++) {
132 xco[0] = valm1(k,j) ; // constant value = boundary value
133 for (int i=1; i<nr; i++) xco[i] = 0 ;
134 valp1.set(k,j) = xco[0] ;
135 xco += nr ;
136 }
137 }
138 }
139 return ;
140 }
141 }
142
143 // Case of a non-zero input
144 // ------------------------
145
146 assert(tin.get_etat() == ETATQCQ ) ;
147 tout.annule_hard() ;
148 valp1.annule_hard() ;
149
150 const double* xci = tin.t ;
151 double* xco = tout.t ;
152
153 for (int k=0 ; k< borne_phi ; k++) {
154 if (k==1) { // jump over the coefficient of sin(0*phi)
155 xci += nr*nt ;
156 xco += nr*nt ;
157 }
158 else {
159 for (int j=0 ; j<nt ; j++) {
160
161 xco[1] = xci[0] - 0.5 * xci[2] ; // special case i = 1
162
163 for (int i=2; i<nr-2; i++) {
164 xco[i] = (xci[i-1] - xci[i+1]) / double(2*i) ;
165 }
166
167 xco[nr-2] = xci[nr-3] / double(2*nr - 4) ;
168 xco[nr-1] = xci[nr-2] / double(2*nr - 2) ;
169
170 // Determination of the T_0 coefficient by matching with
171 // provided value at xi = - 1 :
172 double som = - xco[1] ;
173 for (int i=2; i<nr; i+=2) som += xco[i] ;
174 for (int i=3; i<nr; i+=2) som -= xco[i] ;
175 xco[0] = valm1(k,j) - som ;
176
177 // Value of primitive at xi = + 1 :
178 som = xco[0] ;
179 for (int i=1; i<nr; i++) som += xco[i] ;
180 valp1.set(k,j) = som ;
181
182 xci += nr ;
183 xco += nr ;
184 } // end of theta loop
185 }
186 } // end of phi loop
187
188}
189
190
191
192// case R_CHEBP
193//-------------
194void _primr_r_chebp(const Tbl& tin, int bin, const Tbl&, Tbl& tout,
195 int& bout, Tbl& valp1) {
196
197 assert(tin.dim == tout.dim) ;
198
199 // Output spectral basis
200 int base_t = bin & MSQ_T ;
201 int base_p = bin & MSQ_P ;
202 bout = base_p | base_t | R_CHEBI ;
203
204 // Number of coefficients
205 int nr = tin.get_dim(0) ;
206 int nt = tin.get_dim(1) ;
207 int np = tin.get_dim(2) - 2 ;
208 int borne_phi = np + 1 ;
209 if (np == 1) borne_phi = 1 ;
210
211 // Case of a zero input
212 // --------------------
213 if (tin.get_etat() == ETATZERO) {
214 tout.set_etat_zero() ;
215 valp1.set_etat_zero() ;
216 return ;
217 }
218
219 // Case of a non-zero input
220 // ------------------------
221
222 assert(tin.get_etat() == ETATQCQ ) ;
223 tout.annule_hard() ;
224 valp1.annule_hard() ;
225
226 const double* xci = tin.t ;
227 double* xco = tout.t ;
228
229 for (int k=0 ; k< borne_phi ; k++) {
230 if (k==1) { // jump over the coefficient of sin(0*phi)
231 xci += nr*nt ;
232 xco += nr*nt ;
233 }
234 else {
235 for (int j=0 ; j<nt ; j++) {
236
237 xco[0] = xci[0] - 0.5*xci[1] ; // special case i = 0
238
239 for (int i=1; i<nr-2; i++) {
240 xco[i] = (xci[i] - xci[i+1]) / double(4*i+2) ;
241 }
242
243 xco[nr-2] = xci[nr-2] / double(4*nr - 6) ;
244 xco[nr-1] = 0 ;
245
246 // Value of primitive at xi = + 1 :
247 double som = xco[0] ;
248 for (int i=1; i<nr; i++) som += xco[i] ;
249 valp1.set(k,j) = som ;
250
251 xci += nr ;
252 xco += nr ;
253 } // end of theta loop
254 }
255 } // end of phi loop
256
257}
258
259
260// case R_CHEBI
261//-------------
262void _primr_r_chebi(const Tbl& tin, int bin, const Tbl& val0, Tbl& tout,
263 int& bout, Tbl& valp1) {
264
265 assert(tin.dim == tout.dim) ;
266
267 // Output spectral basis
268 int base_t = bin & MSQ_T ;
269 int base_p = bin & MSQ_P ;
270 bout = base_p | base_t | R_CHEBP ;
271
272 // Number of coefficients
273 int nr = tin.get_dim(0) ;
274 int nt = tin.get_dim(1) ;
275 int np = tin.get_dim(2) - 2 ;
276 int borne_phi = np + 1 ;
277 if (np == 1) borne_phi = 1 ;
278
279
280 // Case of a zero input
281 // --------------------
282 if (tin.get_etat() == ETATZERO) {
283 if (val0.get_etat() == ETATZERO) {
284 tout.set_etat_zero() ;
285 valp1.set_etat_zero() ;
286 return ;
287 }
288 else {
289 assert(val0.get_etat() == ETATQCQ) ;
290 tout.annule_hard() ;
291 valp1.annule_hard() ;
292 double* xco = tout.t ;
293 for (int k=0 ; k< borne_phi ; k++) {
294 if (k==1) { // jump over the coefficient of sin(0*phi)
295 xco += nr*nt ;
296 }
297 else {
298 for (int j=0 ; j<nt ; j++) {
299 xco[0] = val0(k,j) ; // constant value = boundary value
300 for (int i=1; i<nr; i++) xco[i] = 0 ;
301 valp1.set(k,j) = xco[0] ;
302 xco += nr ;
303 }
304 }
305 }
306 return ;
307 }
308 }
309
310 // Case of a non-zero input
311 // ------------------------
312
313 assert(tin.get_etat() == ETATQCQ ) ;
314 tout.annule_hard() ;
315 valp1.annule_hard() ;
316
317 const double* xci = tin.t ;
318 double* xco = tout.t ;
319
320 for (int k=0 ; k< borne_phi ; k++) {
321 if (k==1) { // jump over the coefficient of sin(0*phi)
322 xci += nr*nt ;
323 xco += nr*nt ;
324 }
325 else {
326 for (int j=0 ; j<nt ; j++) {
327
328 for (int i=1; i<nr-1; i++) {
329 xco[i] = (xci[i-1] - xci[i]) / double(4*i) ;
330 }
331
332 xco[nr-1] = xci[nr-2] / double(4*nr - 4) ;
333
334 // Determination of the T_0 coefficient by maching with
335 // provided value at xi = 0 :
336 double som = - xco[1] ;
337 for (int i=2; i<nr; i+=2) som += xco[i] ;
338 for (int i=3; i<nr; i+=2) som -= xco[i] ;
339 xco[0] = val0(k,j) - som ;
340
341 // Value of primitive at xi = + 1 :
342 som = xco[0] ;
343 for (int i=1; i<nr; i++) som += xco[i] ;
344 valp1.set(k,j) = som ;
345
346 xci += nr ;
347 xco += nr ;
348 } // end of theta loop
349 }
350 } // end of phi loop
351
352}
353
354
355
356// case R_CHEBPIM_P
357//-----------------
358void _primr_r_chebpim_p(const Tbl& tin, int bin, const Tbl& val0, Tbl& tout,
359 int& bout, Tbl& valp1) {
360
361 assert(tin.dim == tout.dim) ;
362
363 // Output spectral basis
364 int base_t = bin & MSQ_T ;
365 int base_p = bin & MSQ_P ;
366 bout = base_p | base_t | R_CHEBPIM_I ;
367
368 // Number of coefficients
369 int nr = tin.get_dim(0) ;
370 int nt = tin.get_dim(1) ;
371 int np = tin.get_dim(2) - 2 ;
372 int borne_phi = np + 1 ;
373 if (np == 1) borne_phi = 1 ;
374
375
376 // Case of a zero input
377 // --------------------
378 if (tin.get_etat() == ETATZERO) {
379 if (val0.get_etat() == ETATZERO) {
380 tout.set_etat_zero() ;
381 valp1.set_etat_zero() ;
382 return ;
383 }
384 else {
385 assert(val0.get_etat() == ETATQCQ) ;
386 tout.annule_hard() ;
387 valp1.annule_hard() ;
388 double* xco = tout.t ;
389
390 // m even part
391 for (int k=0 ; k<borne_phi ; k += 4) {
392 int auxiliaire = (k==np) ? 1 : 2 ; // to avoid the last coef
393 for (int kmod=0 ; kmod<auxiliaire ; kmod++) {
394 if ((k==0) && (kmod == 1)) { // jump over the coefficient of sin(0*phi)
395 xco += nr*nt ;
396 }
397 else {
398 for (int j=0 ; j<nt ; j++) {
399 assert( val0(k+kmod,j) == double(0) ) ;
400 for (int i=0; i<nr; i++) xco[i] = 0 ;
401 valp1.set(k+kmod,j) = 0. ;
402 xco += nr ;
403 }
404 }
405 }
406 xco += 2*nr*nt ; // next even m
407 }
408
409 // m odd part
410 xco = tout.t + 2*nr*nt ;
411 for (int k=2 ; k<borne_phi ; k += 4) {
412 int auxiliaire = (k==np) ? 1 : 2 ; // to avoid the last coef
413 for (int kmod=0 ; kmod<auxiliaire ; kmod++) {
414 for (int j=0 ; j<nt ; j++) {
415 xco[0] = val0(k+kmod,j) ; // constant value = boundary value
416 for (int i=1; i<nr; i++) xco[i] = 0 ;
417 valp1.set(k+kmod,j) = xco[0] ;
418 xco += nr ;
419 }
420 }
421 xco += 2*nr*nt ; // next odd m
422 }
423 return ;
424 }
425 }
426
427 // Case of a non-zero input
428 // ------------------------
429
430 assert(tin.get_etat() == ETATQCQ ) ;
431 tout.annule_hard() ;
432 valp1.annule_hard() ;
433
434 const double* xci = tin.t ;
435 double* xco = tout.t ;
436
437 // m even part
438 // -----------
439 for (int k=0 ; k<borne_phi ; k += 4) {
440 int auxiliaire = (k==np) ? 1 : 2 ; // to avoid the last coef
441 for (int kmod=0 ; kmod<auxiliaire ; kmod++) {
442 if ((k==0) && (kmod == 1)) { // jump over the coefficient of sin(0*phi)
443 xci += nr*nt ;
444 xco += nr*nt ;
445 }
446 else {
447 for (int j=0 ; j<nt ; j++) {
448 xco[0] = xci[0] - 0.5*xci[1] ; // special case i = 0
449
450 for (int i=1; i<nr-2; i++) {
451 xco[i] = (xci[i] - xci[i+1]) / double(4*i+2) ;
452 }
453
454 xco[nr-2] = xci[nr-2] / double(4*nr - 6) ;
455 xco[nr-1] = 0 ;
456
457 // Value of primitive at xi = + 1 :
458 double som = xco[0] ;
459 for (int i=1; i<nr; i++) som += xco[i] ;
460 valp1.set(k+kmod,j) = som ;
461
462 xci += nr ;
463 xco += nr ;
464 } // end of theta loop
465
466 }
467 }
468 xci += 2*nr*nt ; // next even m
469 xco += 2*nr*nt ; //
470 }
471
472 // m odd part
473 // ----------
474 xci = tin.t + 2*nr*nt ;
475 xco = tout.t + 2*nr*nt ;
476 for (int k=2 ; k<borne_phi ; k += 4) {
477 int auxiliaire = (k==np) ? 1 : 2 ; // to avoid the last coef
478 for (int kmod=0 ; kmod<auxiliaire ; kmod++) {
479 for (int j=0 ; j<nt ; j++) {
480
481 for (int i=1; i<nr-1; i++) {
482 xco[i] = (xci[i-1] - xci[i]) / double(4*i) ;
483 }
484
485 xco[nr-1] = xci[nr-2] / double(4*nr - 4) ;
486
487 // Determination of the T_0 coefficient by maching with
488 // provided value at xi = 0 :
489 double som = - xco[1] ;
490 for (int i=2; i<nr; i+=2) som += xco[i] ;
491 for (int i=3; i<nr; i+=2) som -= xco[i] ;
492 xco[0] = val0(k+kmod,j) - som ;
493
494 // Value of primitive at xi = + 1 :
495 som = xco[0] ;
496 for (int i=1; i<nr; i++) som += xco[i] ;
497 valp1.set(k+kmod,j) = som ;
498
499 xci += nr ;
500 xco += nr ;
501 } // end of theta loop
502 }
503 xci += 2*nr*nt ; // next odd m
504 xco += 2*nr*nt ; //
505 }
506
507
508}
509
510
511
512// case R_CHEBPIM_I
513//-----------------
514void _primr_r_chebpim_i(const Tbl& tin, int bin, const Tbl& val0, Tbl& tout,
515 int& bout, Tbl& valp1) {
516
517 assert(tin.dim == tout.dim) ;
518
519 // Output spectral basis
520 int base_t = bin & MSQ_T ;
521 int base_p = bin & MSQ_P ;
522 bout = base_p | base_t | R_CHEBPIM_P ;
523
524 // Number of coefficients
525 int nr = tin.get_dim(0) ;
526 int nt = tin.get_dim(1) ;
527 int np = tin.get_dim(2) - 2 ;
528 int borne_phi = np + 1 ;
529 if (np == 1) borne_phi = 1 ;
530
531 // Case of a zero input
532 // --------------------
533 if (tin.get_etat() == ETATZERO) {
534 if (val0.get_etat() == ETATZERO) {
535 tout.set_etat_zero() ;
536 valp1.set_etat_zero() ;
537 return ;
538 }
539 else {
540 assert(val0.get_etat() == ETATQCQ) ;
541 tout.annule_hard() ;
542 valp1.annule_hard() ;
543 double* xco = tout.t ;
544
545 // m odd part
546 for (int k=0 ; k<borne_phi ; k += 4) {
547 int auxiliaire = (k==np) ? 1 : 2 ; // to avoid the last coef
548 for (int kmod=0 ; kmod<auxiliaire ; kmod++) {
549 if ((k==0) && (kmod == 1)) { // jump over the coefficient of sin(0*phi)
550 xco += nr*nt ;
551 }
552 else {
553 for (int j=0 ; j<nt ; j++) {
554 xco[0] = val0(k+kmod,j) ; // constant value = boundary value
555 for (int i=1; i<nr; i++) xco[i] = 0 ;
556 valp1.set(k+kmod,j) = xco[0] ;
557 xco += nr ;
558 }
559 }
560 }
561 xco += 2*nr*nt ; // next odd m
562 }
563
564 // m even part
565 xco = tout.t + 2*nr*nt ;
566 for (int k=2 ; k<borne_phi ; k += 4) {
567 int auxiliaire = (k==np) ? 1 : 2 ; // to avoid the last coef
568 for (int kmod=0 ; kmod<auxiliaire ; kmod++) {
569 for (int j=0 ; j<nt ; j++) {
570 assert( val0(k+kmod,j) == double(0) ) ;
571 for (int i=0; i<nr; i++) xco[i] = 0 ;
572 valp1.set(k+kmod,j) = 0. ;
573 xco += nr ;
574 }
575 }
576 xco += 2*nr*nt ; // next even m
577 }
578 return ;
579 }
580 }
581
582 // Case of a non-zero input
583 // ------------------------
584
585 assert(tin.get_etat() == ETATQCQ ) ;
586 tout.annule_hard() ;
587 valp1.annule_hard() ;
588
589 const double* xci = tin.t ;
590 double* xco = tout.t ;
591
592 // m odd part
593 // ----------
594 for (int k=0 ; k<borne_phi ; k += 4) {
595 int auxiliaire = (k==np) ? 1 : 2 ; // to avoid the last coef
596 for (int kmod=0 ; kmod<auxiliaire ; kmod++) {
597 if ((k==0) && (kmod == 1)) { // jump over the coefficient of sin(0*phi)
598 xci += nr*nt ;
599 xco += nr*nt ;
600 }
601 else {
602 for (int j=0 ; j<nt ; j++) {
603
604 for (int i=1; i<nr-1; i++) {
605 xco[i] = (xci[i-1] - xci[i]) / double(4*i) ;
606 }
607
608 xco[nr-1] = xci[nr-2] / double(4*nr - 4) ;
609
610 // Determination of the T_0 coefficient by maching with
611 // provided value at xi = 0 :
612 double som = - xco[1] ;
613 for (int i=2; i<nr; i+=2) som += xco[i] ;
614 for (int i=3; i<nr; i+=2) som -= xco[i] ;
615 xco[0] = val0(k+kmod,j) - som ;
616
617 // Value of primitive at xi = + 1 :
618 som = xco[0] ;
619 for (int i=1; i<nr; i++) som += xco[i] ;
620 valp1.set(k+kmod,j) = som ;
621
622 xci += nr ;
623 xco += nr ;
624 } // end of theta loop
625 }
626 }
627 xci += 2*nr*nt ; // next odd m
628 xco += 2*nr*nt ; //
629 }
630
631 // m even part
632 // -----------
633 xci = tin.t + 2*nr*nt ;
634 xco = tout.t + 2*nr*nt ;
635 for (int k=2 ; k<borne_phi ; k += 4) {
636 int auxiliaire = (k==np) ? 1 : 2 ; // to avoid the last coef
637 for (int kmod=0 ; kmod<auxiliaire ; kmod++) {
638
639 for (int j=0 ; j<nt ; j++) {
640 xco[0] = xci[0] - 0.5*xci[1] ; // special case i = 0
641
642 for (int i=1; i<nr-2; i++) {
643 xco[i] = (xci[i] - xci[i+1]) / double(4*i+2) ;
644 }
645
646 xco[nr-2] = xci[nr-2] / double(4*nr - 6) ;
647 xco[nr-1] = 0 ;
648
649 // Value of primitive at xi = + 1 :
650 double som = xco[0] ;
651 for (int i=1; i<nr; i++) som += xco[i] ;
652 valp1.set(k+kmod,j) = som ;
653
654 xci += nr ;
655 xco += nr ;
656 } // end of theta loop
657
658 }
659 xci += 2*nr*nt ; // next even m
660 xco += 2*nr*nt ; //
661 }
662
663
664}
665
666// case R_CHEBPI_P
667//-------------
668void _primr_r_chebpi_p(const Tbl& tin, int bin, const Tbl& val0, Tbl& tout,
669 int& bout, Tbl& valp1) {
670
671 assert(tin.dim == tout.dim) ;
672
673 // Output spectral basis
674 int base_t = bin & MSQ_T ;
675 int base_p = bin & MSQ_P ;
676 bout = base_p | base_t | R_CHEBPI_I ;
677
678 // Number of coefficients
679 int nr = tin.get_dim(0) ;
680 int nt = tin.get_dim(1) ;
681 int np = tin.get_dim(2) - 2 ;
682 int borne_phi = np + 1 ;
683 if (np == 1) borne_phi = 1 ;
684
685
686 // Case of a zero input
687 // --------------------
688 if (tin.get_etat() == ETATZERO) {
689 if (val0.get_etat() == ETATZERO) {
690 tout.set_etat_zero() ;
691 valp1.set_etat_zero() ;
692 return ;
693 }
694 else {
695 assert(val0.get_etat() == ETATQCQ) ;
696 tout.annule_hard() ;
697 valp1.annule_hard() ;
698 double* xco = tout.t ;
699 for (int k=0 ; k< borne_phi ; k++) {
700 if (k==1) { // jump over the coefficient of sin(0*phi)
701 xco += nr*nt ;
702 }
703 else {
704 for (int j=0 ; j<nt ; j++) {
705 int l = j%2;
706 if(l==0){
707 for (int i=0; i<nr; i++) xco[i] = 0 ;
708 valp1.set(k,j) = 0. ;
709 } else {
710 xco[0] = val0(k,j) ; // constant value = boundary value
711 for (int i=1; i<nr; i++) xco[i] = 0 ;
712 valp1.set(k,j) = xco[0] ;
713 }
714 xco += nr ;
715 }
716 }
717 }
718 return ;
719 }
720 }
721
722 // Case of a non-zero input
723 // ------------------------
724
725 assert(tin.get_etat() == ETATQCQ ) ;
726 tout.annule_hard() ;
727 valp1.annule_hard() ;
728
729 const double* xci = tin.t ;
730 double* xco = tout.t ;
731
732 for (int k=0 ; k< borne_phi ; k++) {
733 if (k==1) { // jump over the coefficient of sin(0*phi)
734 xci += nr*nt ;
735 xco += nr*nt ;
736 }
737 else {
738 for (int j=0 ; j<nt ; j++) {
739
740 int l = j%2;
741 if(l==0){
742 xco[0] = xci[0] - 0.5*xci[1] ; // special case i = 0
743
744 for (int i=1; i<nr-2; i++) {
745 xco[i] = (xci[i] - xci[i+1]) / double(4*i+2) ;
746 }
747
748 xco[nr-2] = xci[nr-2] / double(4*nr - 6) ;
749 xco[nr-1] = 0 ;
750
751 // Value of primitive at xi = + 1 :
752 double som = xco[0] ;
753 for (int i=1; i<nr; i++) som += xco[i] ;
754 valp1.set(k,j) = som ;
755 } else {
756 for (int i=1; i<nr-1; i++) {
757 xco[i] = (xci[i-1] - xci[i]) / double(4*i) ;
758 }
759
760 xco[nr-1] = xci[nr-2] / double(4*nr - 4) ;
761
762 // Determination of the T_0 coefficient by maching with
763 // provided value at xi = 0 :
764 double som = - xco[1] ;
765 for (int i=2; i<nr; i+=2) som += xco[i] ;
766 for (int i=3; i<nr; i+=2) som -= xco[i] ;
767 xco[0] = val0(k,j) - som ;
768
769 // Value of primitive at xi = + 1 :
770 som = xco[0] ;
771 for (int i=1; i<nr; i++) som += xco[i] ;
772 valp1.set(k,j) = som ;
773 }
774 xci += nr ;
775 xco += nr ;
776 } // end of theta loop
777 }
778 } // end of phi loop
779
780}
781// case R_CHEBPI_I
782//-------------
783void _primr_r_chebpi_i(const Tbl& tin, int bin, const Tbl& val0, Tbl& tout,
784 int& bout, Tbl& valp1) {
785
786 assert(tin.dim == tout.dim) ;
787
788 // Output spectral basis
789 int base_t = bin & MSQ_T ;
790 int base_p = bin & MSQ_P ;
791 bout = base_p | base_t | R_CHEBPI_P ;
792
793 // Number of coefficients
794 int nr = tin.get_dim(0) ;
795 int nt = tin.get_dim(1) ;
796 int np = tin.get_dim(2) - 2 ;
797 int borne_phi = np + 1 ;
798 if (np == 1) borne_phi = 1 ;
799
800
801 // Case of a zero input
802 // --------------------
803 if (tin.get_etat() == ETATZERO) {
804 if (val0.get_etat() == ETATZERO) {
805 tout.set_etat_zero() ;
806 valp1.set_etat_zero() ;
807 return ;
808 }
809 else {
810 assert(val0.get_etat() == ETATQCQ) ;
811 tout.annule_hard() ;
812 valp1.annule_hard() ;
813 double* xco = tout.t ;
814 for (int k=0 ; k< borne_phi ; k++) {
815 if (k==1) { // jump over the coefficient of sin(0*phi)
816 xco += nr*nt ;
817 }
818 else {
819 for (int j=0 ; j<nt ; j++) {
820 int l = j%2;
821 if(l==1){
822 for (int i=0; i<nr; i++) xco[i] = 0 ;
823 valp1.set(k,j) = 0. ;
824 } else {
825 xco[0] = val0(k,j) ; // constant value = boundary value
826 for (int i=1; i<nr; i++) xco[i] = 0 ;
827 valp1.set(k,j) = xco[0] ;
828 }
829 xco += nr ;
830 }
831 }
832 }
833 return ;
834 }
835 }
836
837 // Case of a non-zero input
838 // ------------------------
839
840 assert(tin.get_etat() == ETATQCQ ) ;
841 tout.annule_hard() ;
842 valp1.annule_hard() ;
843
844 const double* xci = tin.t ;
845 double* xco = tout.t ;
846
847 for (int k=0 ; k< borne_phi ; k++) {
848 if (k==1) { // jump over the coefficient of sin(0*phi)
849 xci += nr*nt ;
850 xco += nr*nt ;
851 }
852 else {
853 for (int j=0 ; j<nt ; j++) {
854
855 int l = j%2;
856 if(l==1){
857 xco[0] = xci[0] - 0.5*xci[1] ; // special case i = 0
858
859 for (int i=1; i<nr-2; i++) {
860 xco[i] = (xci[i] - xci[i+1]) / double(4*i+2) ;
861 }
862
863 xco[nr-2] = xci[nr-2] / double(4*nr - 6) ;
864 xco[nr-1] = 0 ;
865
866 // Value of primitive at xi = + 1 :
867 double som = xco[0] ;
868 for (int i=1; i<nr; i++) som += xco[i] ;
869 valp1.set(k,j) = som ;
870 } else {
871 for (int i=1; i<nr-1; i++) {
872 xco[i] = (xci[i-1] - xci[i]) / double(4*i) ;
873 }
874
875 xco[nr-1] = xci[nr-2] / double(4*nr - 4) ;
876
877 // Determination of the T_0 coefficient by maching with
878 // provided value at xi = 0 :
879 double som = - xco[1] ;
880 for (int i=2; i<nr; i+=2) som += xco[i] ;
881 for (int i=3; i<nr; i+=2) som -= xco[i] ;
882 xco[0] = val0(k,j) - som ;
883
884 // Value of primitive at xi = + 1 :
885 som = xco[0] ;
886 for (int i=1; i<nr; i++) som += xco[i] ;
887 valp1.set(k,j) = som ;
888 }
889 xci += nr ;
890 xco += nr ;
891 } // end of theta loop
892 }
893 } // end of phi loop
894
895}
896
897
898// case R_JACO02
899//------------
900void _primr_r_jaco02(const Tbl& tin, int bin, const Tbl& valm1, Tbl& tout,
901 int& bout, Tbl& valp1) {
902
903 assert(tin.dim == tout.dim) ;
904
905 // Output spectral basis
906 bout = bin ;
907
908 // Number of coefficients
909 int nr = tin.get_dim(0) ;
910 int nt = tin.get_dim(1) ;
911 int np = tin.get_dim(2) - 2 ;
912 int borne_phi = np + 1 ;
913 if (np == 1) borne_phi = 1 ;
914
915 // Case of a zero input or pure angular grid
916 // -----------------------------------------
917 if ((tin.get_etat() == ETATZERO)||(nr == 1)) {
918 if (valm1.get_etat() == ETATZERO) {
919 tout.set_etat_zero() ;
920 valp1.set_etat_zero() ;
921 return ;
922 }
923 else {
924 assert(valm1.get_etat() == ETATQCQ) ;
925 tout.set_etat_qcq() ;
926 valp1.set_etat_qcq() ;
927 double* xco = tout.t ;
928 for (int k=0 ; k< borne_phi ; k++) {
929 if (k==1) { // jump over the coefficient of sin(0*phi)
930 xco += nr*nt ;
931 }
932 else {
933 for (int j=0 ; j<nt ; j++) {
934 xco[0] = valm1(k,j) ; // constant value = boundary value
935 for (int i=1; i<nr; i++) xco[i] = 0 ;
936 valp1.set(k,j) = xco[0] ;
937 xco += nr ;
938 }
939 }
940 }
941 return ;
942 }
943 }
944
945 // Case of a non-zero input
946 // ------------------------
947
948 assert(tin.get_etat() == ETATQCQ ) ;
949 tout.annule_hard() ;
950 valp1.annule_hard() ;
951
952 const double* xci = tin.t ;
953 double* xco = tout.t ;
954
955 for (int k=0 ; k< borne_phi ; k++) {
956 if (k==1) { // jump over the coefficient of sin(0*phi)
957 xci += nr*nt ;
958 xco += nr*nt ;
959 }
960 else {
961 for (int j=0 ; j<nt ; j++) {
962
963 for (int i=1; i<nr-1; i++) {
964 xco[i] = (i+2)/double((i+1)*(2*i+1))*xci[i-1] - xci[i]/double((i+1)*(i+2)) - (i+1)/double((i+2)*(2*i+5))*xci[i+1] ;
965 }
966 xco[nr-1] = (nr+1)/double((nr)*(2*nr-1))*xci[nr-2] - xci[nr-1]/double((nr)*(nr+1));
967
968 // Determination of the J_0 coefficient by matching with
969 // provided value at xi = - 1 :
970
971 double som = -3*xco[1] ;
972 for (int i=2; i<nr; i++) {
973 int signe = (i%2 == 0 ? 1 : -1) ;
974 som += xco[i]*signe*(i+1)*(i+2)/double(2) ;
975 }
976 xco[0] = valm1(k,j) - som ;
977
978 // Value of primitive at xi = + 1 :
979 som = xco[0] ;
980 for (int i=1; i<nr; i++) som += xco[i] ;
981 valp1.set(k,j) = som ;
982
983 xci += nr ;
984 xco += nr ;
985 } // end of theta loop
986 }
987 } // end of phi loop
988
989}
990
991
992// case R_LEG
993//-----------
994void _primr_r_leg(const Tbl& tin, int bin, const Tbl& valm1, Tbl& tout,
995 int& bout, Tbl& valp1) {
996
997 assert(tin.dim == tout.dim) ;
998
999 // Output spectral basis
1000 bout = bin ;
1001
1002 // Number of coefficients
1003 int nr = tin.get_dim(0) ;
1004 int nt = tin.get_dim(1) ;
1005 int np = tin.get_dim(2) - 2 ;
1006 int borne_phi = np + 1 ;
1007 if (np == 1) borne_phi = 1 ;
1008
1009 // Case of a zero input or pure angular grid
1010 // -----------------------------------------
1011 if ((tin.get_etat() == ETATZERO)||(nr == 1)) {
1012 if (valm1.get_etat() == ETATZERO) {
1013 tout.set_etat_zero() ;
1014 valp1.set_etat_zero() ;
1015 return ;
1016 }
1017 else {
1018 assert(valm1.get_etat() == ETATQCQ) ;
1019 tout.set_etat_qcq() ;
1020 valp1.set_etat_qcq() ;
1021 double* xco = tout.t ;
1022 for (int k=0 ; k< borne_phi ; k++) {
1023 if (k==1) { // jump over the coefficient of sin(0*phi)
1024 xco += nr*nt ;
1025 }
1026 else {
1027 for (int j=0 ; j<nt ; j++) {
1028 xco[0] = valm1(k,j) ; // constant value = boundary value
1029 for (int i=1; i<nr; i++) xco[i] = 0 ;
1030 valp1.set(k,j) = xco[0] ;
1031 xco += nr ;
1032 }
1033 }
1034 }
1035 return ;
1036 }
1037 }
1038
1039 // Case of a non-zero input
1040 // ------------------------
1041
1042 assert(tin.get_etat() == ETATQCQ ) ;
1043 tout.annule_hard() ;
1044 valp1.annule_hard() ;
1045
1046 const double* xci = tin.t ;
1047 double* xco = tout.t ;
1048
1049 for (int k=0 ; k< borne_phi ; k++) {
1050 if (k==1) { // jump over the coefficient of sin(0*phi)
1051 xci += nr*nt ;
1052 xco += nr*nt ;
1053 }
1054 else {
1055 for (int j=0 ; j<nt ; j++) {
1056
1057 for (int i=1; i<nr-2; i++) {
1058 xco[i] = xci[i-1] / double(2*i-1) - xci[i+1] / double(2*i+3) ;
1059 }
1060
1061 xco[nr-2] = xci[nr-3] / double(2*nr - 5) ;
1062 xco[nr-1] = xci[nr-2] / double(2*nr - 3) ;
1063
1064 // Determination of the T_0 coefficient by matching with
1065 // provided value at xi = - 1 :
1066 double som = - xco[1] ;
1067 for (int i=2; i<nr; i+=2) som += xco[i] ;
1068 for (int i=3; i<nr; i+=2) som -= xco[i] ;
1069 xco[0] = valm1(k,j) - som ;
1070
1071 // Value of primitive at xi = + 1 :
1072 som = xco[0] ;
1073 for (int i=1; i<nr; i++) som += xco[i] ;
1074 valp1.set(k,j) = som ;
1075
1076 xci += nr ;
1077 xco += nr ;
1078 } // end of theta loop
1079 }
1080 } // end of phi loop
1081
1082}
1083
1084
1085
1086// case R_LEGP
1087//-------------
1088void _primr_r_legp(const Tbl& tin, int bin, const Tbl&, Tbl& tout,
1089 int& bout, Tbl& valp1) {
1090
1091 assert(tin.dim == tout.dim) ;
1092
1093 // Output spectral basis
1094 int base_t = bin & MSQ_T ;
1095 int base_p = bin & MSQ_P ;
1096 bout = base_p | base_t | R_LEGI ;
1097
1098 // Number of coefficients
1099 int nr = tin.get_dim(0) ;
1100 int nt = tin.get_dim(1) ;
1101 int np = tin.get_dim(2) - 2 ;
1102 int borne_phi = np + 1 ;
1103 if (np == 1) borne_phi = 1 ;
1104
1105 // Case of a zero input
1106 // --------------------
1107 if ( (tin.get_etat() == ETATZERO) || (nr == 1) ){
1108 tout.set_etat_zero() ;
1109 valp1.set_etat_zero() ;
1110 return ;
1111 }
1112
1113 // Case of a non-zero input
1114 // ------------------------
1115
1116 assert(tin.get_etat() == ETATQCQ ) ;
1117 tout.annule_hard() ;
1118 valp1.annule_hard() ;
1119
1120 const double* xci = tin.t ;
1121 double* xco = tout.t ;
1122
1123 for (int k=0 ; k< borne_phi ; k++) {
1124 if (k==1) { // jump over the coefficient of sin(0*phi)
1125 xci += nr*nt ;
1126 xco += nr*nt ;
1127 }
1128 else {
1129 for (int j=0 ; j<nt ; j++) {
1130
1131 for (int i=0; i<nr-2; i++) {
1132 xco[i] = xci[i]/ double(4*i+1) - xci[i+1]/double(4*i+5) ;
1133 }
1134
1135 xco[nr-2] = xci[nr-2] / double(4*nr - 7) ;
1136 xco[nr-1] = 0 ;
1137
1138 // Value of primitive at xi = + 1 :
1139 double som = xco[0] ;
1140 for (int i=1; i<nr; i++) som += xco[i] ;
1141 valp1.set(k,j) = som ;
1142
1143 xci += nr ;
1144 xco += nr ;
1145 } // end of theta loop
1146 }
1147 } // end of phi loop
1148
1149}
1150
1151
1152// case R_LEGI
1153//------------
1154void _primr_r_legi(const Tbl& tin, int bin, const Tbl& val0, Tbl& tout,
1155 int& bout, Tbl& valp1) {
1156
1157 assert(tin.dim == tout.dim) ;
1158
1159 // Output spectral basis
1160 int base_t = bin & MSQ_T ;
1161 int base_p = bin & MSQ_P ;
1162 bout = base_p | base_t | R_LEGP ;
1163
1164 // Number of coefficients
1165 int nr = tin.get_dim(0) ;
1166 int nt = tin.get_dim(1) ;
1167 int np = tin.get_dim(2) - 2 ;
1168 int borne_phi = np + 1 ;
1169 if (np == 1) borne_phi = 1 ;
1170
1171
1172 // Case of a zero input
1173 // --------------------
1174 if ( (tin.get_etat() == ETATZERO) || (nr == 1) ){
1175 if (val0.get_etat() == ETATZERO) {
1176 tout.set_etat_zero() ;
1177 valp1.set_etat_zero() ;
1178 return ;
1179 }
1180 else {
1181 assert(val0.get_etat() == ETATQCQ) ;
1182 tout.annule_hard() ;
1183 valp1.annule_hard() ;
1184 double* xco = tout.t ;
1185 for (int k=0 ; k< borne_phi ; k++) {
1186 if (k==1) { // jump over the coefficient of sin(0*phi)
1187 xco += nr*nt ;
1188 }
1189 else {
1190 for (int j=0 ; j<nt ; j++) {
1191 xco[0] = val0(k,j) ; // constant value = boundary value
1192 for (int i=1; i<nr; i++) xco[i] = 0 ;
1193 valp1.set(k,j) = xco[0] ;
1194 xco += nr ;
1195 }
1196 }
1197 }
1198 return ;
1199 }
1200 }
1201
1202 // Case of a non-zero input
1203 // ------------------------
1204
1205 assert(tin.get_etat() == ETATQCQ ) ;
1206 tout.annule_hard() ;
1207 valp1.annule_hard() ;
1208
1209 const double* xci = tin.t ;
1210 double* xco = tout.t ;
1211
1212 for (int k=0 ; k< borne_phi ; k++) {
1213 if (k==1) { // jump over the coefficient of sin(0*phi)
1214 xci += nr*nt ;
1215 xco += nr*nt ;
1216 }
1217 else {
1218 for (int j=0 ; j<nt ; j++) {
1219
1220 for (int i=1; i<nr-1; i++) {
1221 xco[i] = xci[i-1]/ double(4*i-1) - xci[i]/double(4*i+3) ;
1222 }
1223
1224 xco[nr-1] = xci[nr-2] / double(4*nr - 5) ;
1225
1226 // Determination of the T_0 coefficient by matching with
1227 // provided value at xi = 0 :
1228 double val = -0.5 ;
1229 double som = val*xco[1] ;
1230 for (int i=2; i<nr; i++) {
1231 val *= -double(2*i-1) / double(2*i) ;
1232 som += val*xco[i] ;
1233 }
1234 xco[0] = val0(k,j) - som ;
1235
1236 // Value of primitive at xi = + 1 :
1237 som = xco[0] ;
1238 for (int i=1; i<nr; i++) som += xco[i] ;
1239 valp1.set(k,j) = som ;
1240
1241 xci += nr ;
1242 xco += nr ;
1243 } // end of theta loop
1244 }
1245 } // end of phi loop
1246
1247}
1248
1249
1250
1251
1252}
Lorene::Tbl
Basic array class.
Definition tbl.h:161
R_LEGP
#define R_LEGP
base de Legendre paire (rare) seulement
Definition type_parite.h:184
R_LEGI
#define R_LEGI
base de Legendre impaire (rare) seulement
Definition type_parite.h:186
R_CHEBI
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition type_parite.h:170
R_CHEBPIM_I
#define R_CHEBPIM_I
Cheb. pair-impair suivant m, impair pour m=0.
Definition type_parite.h:178
R_CHEBPI_I
#define R_CHEBPI_I
Cheb. pair-impair suivant l impair pour l=0.
Definition type_parite.h:174
R_CHEBPIM_P
#define R_CHEBPIM_P
Cheb. pair-impair suivant m, pair pour m=0.
Definition type_parite.h:176
MSQ_T
#define MSQ_T
Extraction de l'info sur Theta.
Definition type_parite.h:154
R_CHEBP
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition type_parite.h:168
MSQ_P
#define MSQ_P
Extraction de l'info sur Phi.
Definition type_parite.h:156
R_CHEBPI_P
#define R_CHEBPI_P
Cheb. pair-impair suivant l pair pour l=0.
Definition type_parite.h:172
Lorene
Lorene prototypes.
Definition app_hor.h:67