LORENE
lindquist.C
1/*
2 * Copyright (c) 2000-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 lindquist_C[] = "$Header: /cvsroot/Lorene/C++/Source/Bhole_binaire/lindquist.C,v 1.4 2014/10/13 08:52:40 j_novak Exp $" ;
24
25/*
26 * $Id: lindquist.C,v 1.4 2014/10/13 08:52:40 j_novak Exp $
27 * $Log: lindquist.C,v $
28 * Revision 1.4 2014/10/13 08:52:40 j_novak
29 * Lorene classes and functions now belong to the namespace Lorene.
30 *
31 * Revision 1.3 2014/10/06 15:12:58 j_novak
32 * Modified #include directives to use c++ syntax.
33 *
34 * Revision 1.2 2002/10/16 14:36:33 j_novak
35 * Reorganization of #include instructions of standard C++, in order to
36 * use experimental version 3 of gcc.
37 *
38 * Revision 1.1.1.1 2001/11/20 15:19:28 e_gourgoulhon
39 * LORENE
40 *
41 * Revision 2.0 2000/12/13 15:42:57 phil
42 * *** empty log message ***
43 *
44 *
45 * $Header: /cvsroot/Lorene/C++/Source/Bhole_binaire/lindquist.C,v 1.4 2014/10/13 08:52:40 j_novak Exp $
46 *
47 */
48
49
50//standard
51#include <cstdlib>
52#include <cmath>
53
54// LORENE
55#include "type_parite.h"
56#include "nbr_spx.h"
57#include "proto.h"
58#include "coord.h"
59#include "tenseur.h"
60
61namespace Lorene {
62double serie_lindquist_plus (double rayon, double distance, double xa, double ya,
63 double za, double precision, double itemax) {
64
65
66 double result = 0.5 ;
67 double c_n = rayon ;
68 double d_n = distance/2. ;
69
70 int indic = 1 ;
71 int conte=0 ;
72 while ((indic == 1) && (conte <= itemax)) {
73
74 double norme_plus = sqrt ((xa+d_n)*(xa+d_n)+ya*ya+za*za) ;
75
76 double terme = c_n * 1./norme_plus ;
77 if (fabs(terme/result) < precision)
78 indic = -1 ;
79
80 result = result + terme ;
81
82 c_n *= rayon/(distance/2. + d_n) ;
83 d_n = distance/2. - rayon*rayon/(distance/2.+d_n) ;
84 conte ++ ;
85 }
86
87 if (conte > itemax)
88 result = 1 ;
89return result ;
90}
91
92double serie_lindquist_moins (double rayon, double distance, double xa, double ya,
93 double za, double precision, double itemax) {
94
95
96 double result = 0.5 ;
97 double c_n = rayon ;
98 double d_n = distance/2. ;
99
100 int indic = 1 ;
101 int conte=0 ;
102 while ((indic == 1) && (conte <= itemax)) {
103
104 double norme_plus = sqrt ((xa-d_n)*(xa-d_n)+ya*ya+za*za) ;
105
106 double terme = c_n * 1./norme_plus ;
107 if (fabs(terme/result) < precision)
108 indic = -1 ;
109
110 result = result + terme ;
111
112 c_n *= rayon/(distance/2. + d_n) ;
113 d_n = distance/2. - rayon*rayon/(distance/2.+d_n) ;
114 conte ++ ;
115 }
116
117 if (conte > itemax)
118 result = 1 ;
119return result ;
120}
121
122double adm_serie (double rayon, double distance, double precision) {
123
124 double result = 0 ;
125 double c_n = rayon ;
126 double d_n = distance/2. ;
127
128 int indic =1 ;
129 while (indic == 1) {
130
131 result += 4*c_n ;
132 if (fabs(c_n/result) < precision)
133 indic = -1 ;
134
135 c_n *= rayon/(distance/2. + d_n) ;
136 d_n = distance/2. - rayon*rayon/(distance/2.+d_n) ;
137 }
138return result ;
139}
140
141double bare_serie (double rayon, double distance, double precision) {
142
143 double result = 0 ;
144 double c_n = rayon ;
145 double d_n = distance/2. ;
146
147 int indic =1 ;
148 int n = 1 ;
149 while (indic == 1) {
150
151 result += 4*c_n*n ;
152 if (fabs(c_n/result) < precision)
153 indic = -1 ;
154
155 c_n *= rayon/(distance/2. + d_n) ;
156 d_n = distance/2. - rayon*rayon/(distance/2.+d_n) ;
157 n++ ;
158 }
159return result ;
160}
161
162void set_lindquist (Cmp& psi_un, Cmp& psi_deux, double rayon, double precision) {
163
164 // Pour les allocations !
165 psi_un = 1 ;
166 psi_deux = 1 ;
167
168 double distance = psi_un.get_mp()->get_ori_x()-psi_deux.get_mp()->get_ori_x() ;
169
170 // On regarde psi 1 :
171 Mtbl xa_mtbl_un (psi_un.get_mp()->xa) ;
172 Mtbl ya_mtbl_un (psi_un.get_mp()->ya) ;
173 Mtbl za_mtbl_un (psi_un.get_mp()->za) ;
174
175 int nz_un = psi_un.get_mp()->get_mg()->get_nzone() ;
176 for (int l=1 ; l<nz_un ; l++) {
177 int np = psi_un.get_mp()->get_mg()->get_np (l) ;
178 int nt = psi_un.get_mp()->get_mg()->get_nt (l) ;
179 int nr = psi_un.get_mp()->get_mg()->get_nr (l) ;
180 double xa, ya, za ;
181 for (int k=0 ; k<np ; k++)
182 for (int j=0 ; j<nt ; j++)
183 for (int i=0 ; i<nr ; i++) {
184 xa = xa_mtbl_un (l, k, j, i) ;
185 ya = ya_mtbl_un (l, k, j, i) ;
186 za = za_mtbl_un (l, k, j, i) ;
187
188 psi_un.set(l, k, j, i) =
189serie_lindquist_moins (rayon, distance, xa, ya, za, precision, 30) ;
190 }
191 }
192
193 psi_un.set_val_inf (0.5) ;
194 psi_un.std_base_scal() ;
195
196 // On regarde psi 2 :
197 Mtbl xa_mtbl_deux (psi_deux.get_mp()->xa) ;
198 Mtbl ya_mtbl_deux (psi_deux.get_mp()->ya) ;
199 Mtbl za_mtbl_deux (psi_deux.get_mp()->za) ;
200
201 int nz_deux = psi_deux.get_mp()->get_mg()->get_nzone() ;
202 for (int l=1 ; l<nz_deux ; l++) {
203 int np = psi_deux.get_mp()->get_mg()->get_np (l) ;
204 int nt = psi_deux.get_mp()->get_mg()->get_nt (l) ;
205 int nr = psi_deux.get_mp()->get_mg()->get_nr (l) ;
206 double xa, ya, za ;
207 for (int k=0 ; k<np ; k++)
208 for (int j=0 ; j<nt ; j++)
209 for (int i=0 ; i<nr ; i++) {
210 xa = xa_mtbl_deux (l, k, j, i) ;
211 ya = ya_mtbl_deux (l, k, j, i) ;
212 za = za_mtbl_deux (l, k, j, i) ;
213
214 psi_deux.set(l, k, j, i) =
215serie_lindquist_plus (rayon, distance, xa, ya, za, precision, 30) ;
216 }
217 }
218 psi_deux.set_val_inf (0.5) ;
219 psi_deux.std_base_scal() ;
220}
221}
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:220
Lorene prototypes.
Definition app_hor.h:64