LORENE
val_dern_1d.C
1/*
2 * Copyright (c) 2004 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 val_dern_1d_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/val_dern_1d.C,v 1.2 2014/10/13 08:53:27 j_novak Exp $" ;
24
25/*
26 * $Id: val_dern_1d.C,v 1.2 2014/10/13 08:53:27 j_novak Exp $
27 * $Log: val_dern_1d.C,v $
28 * Revision 1.2 2014/10/13 08:53:27 j_novak
29 * Lorene classes and functions now belong to the namespace Lorene.
30 *
31 * Revision 1.1 2004/02/17 09:21:39 j_novak
32 * New functions for calculating values of the derivatives of a function
33 * using its Chebyshev coefficients.
34 *
35 *
36 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/val_dern_1d.C,v 1.2 2014/10/13 08:53:27 j_novak Exp $
37 *
38 */
39
40#include "type_parite.h"
41#include "tbl.h"
42
43/*
44 * Functions computing value of f^(n) at boundaries of the interval [-1, 1],
45 * using the Chebyshev expansion of f. Note: n=0 works too.
46 *
47 * Input : 1-dimensional Tbl containing the Chebyshev coefficients of f.
48 * int base : base of spectral expansion.
49 *
50 * Output : double : the value of the n-th derivative of f at x=+/- 1.
51 *
52 */
53
54
55namespace Lorene {
56double val1_dern_1d(int n, const Tbl& tb, int base_r)
57{
58
59 //This function should be OK for any radial base
60 assert ( (base_r == R_CHEB) || (base_r == R_CHEBI) || (base_r == R_CHEBP) ||
61 (base_r == R_CHEBU) ) ;
62
63 assert (n>=0) ;
64 assert (tb.get_ndim() == 1) ;
65 int nr = tb.get_dim(0) ;
66
67 double resu = 0. ;
68
69 int n_ini = ( (base_r == R_CHEBP) || (base_r == R_CHEBI) ) ? n / 2 : n ;
70
71 double *tbi = &tb.t[n_ini] ;
72 for (int i=n_ini; i<nr; i++) {
73 double fact = 1. ;
74 int ii = i ;
75 if (base_r == R_CHEBP) ii *= 2 ;
76 if (base_r == R_CHEBI) ii = 2*i + 1 ;
77 for (int j=0; j<n; j++)
78 fact *= double(ii*ii - j*j)/double(2*j + 1) ;
79 resu += fact * (*tbi) ;
80 tbi++ ;
81 }
82
83 return resu ;
84}
85
86double valm1_dern_1d(int n, const Tbl& tb, int base_r)
87{
88
89 //This function should be OK for any radial base
90 assert ( (base_r == R_CHEB) || (base_r == R_CHEBI) || (base_r == R_CHEBP) ||
91 (base_r == R_CHEBU) ) ;
92
93 assert (n>=0) ;
94 assert (tb.get_ndim() == 1) ;
95 int nr = tb.get_dim(0) ;
96
97 double resu = 0. ;
98 double parite, fac ;
99 int n_ini ;
100 switch (base_r) {
101 case R_CHEBP:
102 n_ini = n / 2 ;
103 parite = 1 ;
104 fac = (n%2 == 0 ? 1 : -1) ;
105 break ;
106 case R_CHEBI:
107 n_ini = n / 2 ;
108 fac = (n%2 == 0 ? -1 : 1) ;
109 parite = 1 ;
110 break ;
111 default:
112 n_ini = n ;
113 parite = -1 ;
114 fac = 1 ;
115 break ;
116 }
117 double *tbi = &tb.t[n_ini] ;
118
119 for (int i=n_ini; i<nr; i++) {
120 double fact = fac ;
121 int ii = i ;
122 if (base_r == R_CHEBP) ii *= 2 ;
123 if (base_r == R_CHEBI) ii = 2*i + 1 ;
124 for (int j=0; j<n; j++)
125 fact *= double(ii*ii - j*j)/double(2*j + 1) ;
126 resu += fact * (*tbi) ;
127 fac *= parite ;
128 tbi++ ;
129 }
130
131 return resu ;
132}
133}
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
#define R_CHEBI
base de Cheb. impaire (rare) seulement
#define R_CHEB
base de Chebychev ordinaire (fin)
#define R_CHEBP
base de Cheb. paire (rare) seulement
Lorene prototypes.
Definition app_hor.h:64