LORENE
sx_1d.C
1/*
2 * Copyright (c) 2006 Jerome Novak
3 *
4 * This file is part of LORENE.
5 *
6 * LORENE is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * LORENE is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with LORENE; if not, write to the Free Software
18 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 *
20 */
21
22
23char sx_1d_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/sx_1d.C,v 1.4 2015/03/05 08:49:32 j_novak Exp $" ;
24
25/*
26 * $Id: sx_1d.C,v 1.4 2015/03/05 08:49:32 j_novak Exp $
27 * $Log: sx_1d.C,v $
28 * Revision 1.4 2015/03/05 08:49:32 j_novak
29 * Implemented operators with Legendre bases.
30 *
31 * Revision 1.3 2014/10/13 08:53:27 j_novak
32 * Lorene classes and functions now belong to the namespace Lorene.
33 *
34 * Revision 1.2 2014/10/06 15:16:07 j_novak
35 * Modified #include directives to use c++ syntax.
36 *
37 * Revision 1.1 2006/04/10 15:19:20 j_novak
38 * New definition of 1D operators dsdx and sx in the nucleus (bases R_CHEBP and
39 * R_CHEBI).
40 *
41 *
42 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/sx_1d.C,v 1.4 2015/03/05 08:49:32 j_novak Exp $
43 *
44 */
45
46
47 // Includes
48#include <cstdlib>
49
50#include "headcpp.h"
51#include "type_parite.h"
52#include "proto.h"
53
54
55/*
56 * 1/x operator (division by x)
57 *
58 * Only for the bases R_CHEBP and R_CHEBI
59 *
60 * Input :
61 *
62 * int nr : number of coefficients
63 * tb : the array of coefficients of the input function
64 *
65 * Output :
66 * tb : the array of coefficients of the result
67 */
68
69
70 //----------------------------
71 // Bases not implemented -----
72 //----------------------------
73
74namespace Lorene {
75void _sx_1d_pas_prevu(int , double* , double *) {
76 cout << "sx_1d : base not implemented..." << endl ;
77 abort() ;
78 exit(-1) ;
79}
80
81
82 //-------------------
83 // case R_CHEBP ---
84 //-------------------
85
86void _sx_1d_r_chebp (int nr, double* tb, double* res) {
87
88 double som ;
89 int sign = 1 ;
90 res[nr-1] = 0 ;
91 som = 2 * sign * tb[nr-1] ;
92 res[nr-2] = som ;
93 for (int i=nr-3 ; i>=0 ; i--) {
94 sign = - sign ;
95 som += 2 * sign * tb[i+1] ;
96 res[i] = (i%2 == 0 ? -1 : 1) * som ;
97 }
98
99}
100
101
102 //-------------------
103 // case R_CHEBI ---
104 //-------------------
105
106void _sx_1d_r_chebi (int nr, double* tb, double* res) {
107
108 double som ;
109 int sign = 1 ;
110 res[nr-1] = 0 ;
111 som = 2 * sign * tb[nr-2] ;
112 res[nr-2] = som ;
113 for (int i=nr-3 ; i>=0 ; i--) {
114 sign = - sign ;
115 som += 2 * sign * tb[i] ;
116 res[i] = (i%2 == 0 ? -1 : 1) * som ;
117 }
118 res[0] *= 0.5 ;
119}
120
121 //-------------------
122 // case R_CHEBU ---
123 //-------------------
124
125void _sx_1d_r_chebu (int nr, double* tb, double* res) {
126
127 for (int i=0; i<nr; i++)
128 res[i] = tb[i] ;
129
130 sxm1_1d_cheb(nr, res) ;
131
132}
133
134 //------------------
135 // case R_LEGP ---
136 //------------------
137
138void _sx_1d_r_legp (int nr, double* tb, double* res) {
139
140 double term = 0 ;
141 res[nr-1] = 0 ;
142 for (int i=nr-2 ; i>=0 ; i--) {
143 term += tb[i+1] ;
144 res[i] = double(4*i+3)/double(2*i+2)*term ;
145 term *= -double(2*i+1)/double(2*i+2) ;
146 }
147
148}
149
150
151 //------------------
152 // case R_LEGI ---
153 //------------------
154
155void _sx_1d_r_legi (int nr, double* tb, double* res) {
156
157 double term = 0 ;
158 res[nr-1] = 0 ;
159 for (int i=nr-2; i>=0; i--) {
160 term += tb[i] ;
161 res[i] = double(4*i+1)/double(2*i+1)*term ;
162 term *= -double(2*i)/double(2*i+1) ;
163 }
164}
165
166
167
168 // ----------------------
169 // The calling function
170 //-----------------------
171
172void sx_1d(int nr, double **tb, int base_r)
173{
174
175 static void (*sx_1d[MAX_BASE])(int, double *, double *) ;
176 static int nap = 0 ;
177
178 // Premier appel
179 if (nap==0) {
180 nap = 1 ;
181 for (int i=0 ; i<MAX_BASE ; i++) {
182 sx_1d[i] = _sx_1d_pas_prevu ;
183 }
184 // Les routines existantes
185 sx_1d[R_CHEBP >> TRA_R] = _sx_1d_r_chebp ;
186 sx_1d[R_CHEBI >> TRA_R] = _sx_1d_r_chebi ;
187 sx_1d[R_LEGP >> TRA_R] = _sx_1d_r_legp ;
188 sx_1d[R_LEGI >> TRA_R] = _sx_1d_r_legi ;
189 sx_1d[R_CHEBU >> TRA_R] = _sx_1d_r_chebu ;
190 }
191
192
193 double *result = new double[nr] ;
194 sx_1d[base_r](nr, *tb, result) ;
195 delete [] (*tb) ;
196 (*tb) = result ;
197}
198
199}
#define R_LEGP
base de Legendre paire (rare) seulement
#define MAX_BASE
Nombre max. de bases differentes.
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
#define R_LEGI
base de Legendre impaire (rare) seulement
#define R_CHEBI
base de Cheb. impaire (rare) seulement
#define TRA_R
Translation en R, used for a bitwise shift (in hex)
#define R_CHEBP
base de Cheb. paire (rare) seulement
Lorene prototypes.
Definition app_hor.h:64