LORENE
solh_helmholtz_minus.C
1/*
2 * Copyright (c) 1999-2001 Philippe Grandclement
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 solh_helmholtz_minusC[] = "$Header $" ;
24
25/*
26 * $Id: solh_helmholtz_minus.C,v 1.9 2014/10/13 08:53:31 j_novak Exp $
27 * $Log: solh_helmholtz_minus.C,v $
28 * Revision 1.9 2014/10/13 08:53:31 j_novak
29 * Lorene classes and functions now belong to the namespace Lorene.
30 *
31 * Revision 1.8 2014/10/06 15:16:10 j_novak
32 * Modified #include directives to use c++ syntax.
33 *
34 * Revision 1.7 2008/07/10 11:20:33 p_grandclement
35 * mistake fixed in solh_helmholtz_minus
36 *
37 * Revision 1.6 2008/07/08 11:45:28 p_grandclement
38 * Add helmholtz_minus in the nucleus
39 *
40 * Revision 1.5 2008/02/18 13:53:43 j_novak
41 * Removal of special indentation instructions.
42 *
43 * Revision 1.4 2004/08/24 10:11:12 p_grandclement
44 * Correction of the includes of gsl
45 *
46 * Revision 1.3 2004/08/24 09:14:44 p_grandclement
47 * Addition of some new operators, like Poisson in 2d... It now requieres the
48 * GSL library to work.
49 *
50 * Also, the way a variable change is stored by a Param_elliptic is changed and
51 * no longer uses Change_var but rather 2 Scalars. The codes using that feature
52 * will requiere some modification. (It should concern only the ones about monopoles)
53 *
54 * Revision 1.2 2004/01/15 09:15:37 p_grandclement
55 * Modification and addition of the Helmholtz operators
56 *
57 * Revision 1.1 2003/12/11 14:48:49 p_grandclement
58 * Addition of ALL (and that is a lot !) the files needed for the general elliptic solver ... UNDER DEVELOPEMENT...
59 *
60 *
61 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/solh_helmholtz_minus.C,v 1.9 2014/10/13 08:53:31 j_novak Exp $
62 *
63 */
64
65//fichiers includes
66#include <cstdio>
67#include <cstdlib>
68#include <cmath>
69#include <gsl/gsl_sf_bessel.h>
70
71#include "proto.h"
72#include "matrice.h"
73#include "type_parite.h"
74
75
76 //------------------------------------
77 // Routine pour les cas non prevus --
78 //------------------------------------
79namespace Lorene {
80Tbl _solh_helmholtz_minus_pas_prevu (int, int, double, double, double) {
81
82 cout << "Homogeneous solution not implemented in hemlholtz_minus : "<< endl ;
83 abort() ;
84 exit(-1) ;
85 Tbl res(1) ;
86 return res;
87}
88
89
90
91 //-------------------
92 //-- R_CHEB ------
93 //-------------------
94
95Tbl _solh_helmholtz_minus_r_cheb (int n, int lq, double alpha, double beta,
96 double masse) {
97
98 assert (masse > 0) ;
99
100 Tbl res(2,n) ;
101 res.set_etat_qcq() ;
102 double* coloc = new double[n] ;
103
104 int * deg = new int[3] ;
105 deg[0] = 1 ;
106 deg[1] = 1 ;
107 deg[2] = n ;
108
109 // First SH
110 for (int i=0 ; i<n ; i++){
111 double air = alpha*(-cos(M_PI*i/(n-1))) + beta ;
112 coloc[i] = gsl_sf_bessel_il_scaled (lq, masse*air)/exp(-masse*air) ;
113 }
114
115 cfrcheb(deg, deg, coloc, deg, coloc) ;
116 for (int i=0 ; i<n ;i++)
117 res.set(0,i) = coloc[i] ;
118
119 // Second SH
120 for (int i=0 ; i<n ; i++){
121 double air = alpha*(-cos(M_PI*i/(n-1))) + beta ;
122 coloc[i] = gsl_sf_bessel_kl_scaled (lq, masse*air) / exp(masse*air) ;
123 }
124
125 cfrcheb(deg, deg, coloc, deg, coloc) ;
126 for (int i=0 ; i<n ;i++)
127 res.set(1,i) = coloc[i] ;
128
129 delete [] deg ;
130 delete [] coloc ;
131 return res ;
132}
133 //-------------------
134 //-- R_CHEBP ------
135 //-------------------
136
137Tbl _solh_helmholtz_minus_r_chebp (int n, int lq, double alpha, double,
138 double masse) {
139
140 assert (masse > 0) ;
141
142 Tbl res(n) ;
143 res.set_etat_qcq() ;
144 double* coloc = new double[n] ;
145
146 int * deg = new int[3] ;
147 deg[0] = 1 ;
148 deg[1] = 1 ;
149 deg[2] = n ;
150
151 // First SH
152 for (int i=0 ; i<n ; i++){
153 double air = alpha*(sin(M_PI*i/2./(n-1))) ;
154 coloc[i] = gsl_sf_bessel_il_scaled (lq, masse*air)/exp(-masse*air) ;
155 }
156
157 cfrchebp(deg, deg, coloc, deg, coloc) ;
158 for (int i=0 ; i<n ;i++)
159 res.set(i) = coloc[i] ;
160
161 delete [] deg ;
162 delete [] coloc ;
163 return res ;
164}
165 //-------------------
166 //-- R_CHEBI ------
167 //-------------------
168
169Tbl _solh_helmholtz_minus_r_chebi (int n, int lq, double alpha, double,
170 double masse) {
171
172 assert (masse > 0) ;
173
174 Tbl res(n) ;
175 res.set_etat_qcq() ;
176 double* coloc = new double[n] ;
177
178 int * deg = new int[3] ;
179 deg[0] = 1 ;
180 deg[1] = 1 ;
181 deg[2] = n ;
182
183 // First SH
184 for (int i=0 ; i<n ; i++){
185 double air = alpha*(sin(M_PI*i/2./(n-1))) ;
186 coloc[i] = gsl_sf_bessel_il_scaled (lq, masse*air)/exp(-masse*air) ;
187 }
188
189 cfrchebi(deg, deg, coloc, deg, coloc) ;
190 for (int i=0 ; i<n ;i++)
191 res.set(i) = coloc[i] ;
192
193 delete [] deg ;
194 delete [] coloc ;
195 return res ;
196}
197
198 //-------------------
199 //-- R_CHEBU ------
200 //-------------------
201
202Tbl _solh_helmholtz_minus_r_chebu (int n, int lq,
203 double alpha, double, double masse) {
204
205 assert (masse > 0) ;
206
207 Tbl res(n) ;
208 res.set_etat_qcq() ;
209 double* coloc = new double[n] ;
210
211 int * deg = new int[3] ;
212 deg[0] = 1 ;
213 deg[1] = 1 ;
214 deg[2] = n ;
215
216 for (int i=0 ; i<n-1 ; i++){
217 double air = 1./(alpha*(-1-cos(M_PI*i/(n-1)))) ;
218 coloc[i] = gsl_sf_bessel_kl_scaled (lq, masse*air) / exp(masse*air) ;
219 }
220 coloc[n-1] = 0 ;
221
222 cfrcheb(deg, deg, coloc, deg, coloc) ;
223 for (int i=0 ; i<n ;i++)
224 res.set(i) = coloc[i] ;
225
226 delete [] deg ;
227 delete [] coloc ;
228 return res ;
229}
230
231
232 //-------------------
233 //-- Fonction ----
234 //-------------------
235
236
237Tbl solh_helmholtz_minus (int n, int lq, double alpha, double beta,
238 double masse, int base_r) {
239
240 // Routines de derivation
241 static Tbl (*solh_helmholtz_minus[MAX_BASE])(int, int, double, double, double) ;
242 static int nap = 0 ;
243
244 // Premier appel
245 if (nap==0) {
246 nap = 1 ;
247 for (int i=0 ; i<MAX_BASE ; i++) {
248 solh_helmholtz_minus[i] = _solh_helmholtz_minus_pas_prevu ;
249 }
250 // Les routines existantes
251 solh_helmholtz_minus[R_CHEB >> TRA_R] = _solh_helmholtz_minus_r_cheb ;
252 solh_helmholtz_minus[R_CHEBU >> TRA_R] = _solh_helmholtz_minus_r_chebu ;
253 solh_helmholtz_minus[R_CHEBP >> TRA_R] = _solh_helmholtz_minus_r_chebp ;
254 solh_helmholtz_minus[R_CHEBI >> TRA_R] = _solh_helmholtz_minus_r_chebi ;
255 }
256
257 Tbl res(solh_helmholtz_minus[base_r](n, lq, alpha, beta, masse)) ;
258 return res ;
259}
260}
Cmp exp(const Cmp &)
Exponential.
Definition cmp_math.C:270
Cmp sin(const Cmp &)
Sine.
Definition cmp_math.C:69
Cmp cos(const Cmp &)
Cosine.
Definition cmp_math.C:94
#define MAX_BASE
Nombre max. de bases differentes.
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
#define R_CHEBI
base de Cheb. impaire (rare) seulement
#define TRA_R
Translation en R, used for a bitwise shift (in hex)
#define R_CHEB
base de Chebychev ordinaire (fin)
#define R_CHEBP
base de Cheb. paire (rare) seulement
Lorene prototypes.
Definition app_hor.h:64