23char lindquist_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Bhole_binaire/lindquist.C,v 1.4 2014/10/13 08:52:40 j_novak Exp $" ;
55#include "type_parite.h"
62double serie_lindquist_plus (
double rayon,
double distance,
double xa,
double ya,
63 double za,
double precision,
double itemax) {
68 double d_n = distance/2. ;
72 while ((indic == 1) && (conte <= itemax)) {
74 double norme_plus =
sqrt ((xa+d_n)*(xa+d_n)+ya*ya+za*za) ;
76 double terme = c_n * 1./norme_plus ;
77 if (fabs(terme/result) < precision)
80 result = result + terme ;
82 c_n *= rayon/(distance/2. + d_n) ;
83 d_n = distance/2. - rayon*rayon/(distance/2.+d_n) ;
92double serie_lindquist_moins (
double rayon,
double distance,
double xa,
double ya,
93 double za,
double precision,
double itemax) {
98 double d_n = distance/2. ;
102 while ((indic == 1) && (conte <= itemax)) {
104 double norme_plus =
sqrt ((xa-d_n)*(xa-d_n)+ya*ya+za*za) ;
106 double terme = c_n * 1./norme_plus ;
107 if (fabs(terme/result) < precision)
110 result = result + terme ;
112 c_n *= rayon/(distance/2. + d_n) ;
113 d_n = distance/2. - rayon*rayon/(distance/2.+d_n) ;
122double adm_serie (
double rayon,
double distance,
double precision) {
126 double d_n = distance/2. ;
132 if (fabs(c_n/result) < precision)
135 c_n *= rayon/(distance/2. + d_n) ;
136 d_n = distance/2. - rayon*rayon/(distance/2.+d_n) ;
141double bare_serie (
double rayon,
double distance,
double precision) {
145 double d_n = distance/2. ;
152 if (fabs(c_n/result) < precision)
155 c_n *= rayon/(distance/2. + d_n) ;
156 d_n = distance/2. - rayon*rayon/(distance/2.+d_n) ;
162void set_lindquist (Cmp& psi_un, Cmp& psi_deux,
double rayon,
double precision) {
168 double distance = psi_un.get_mp()->get_ori_x()-psi_deux.get_mp()->get_ori_x() ;
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) ;
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) ;
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) ;
188 psi_un.set(l, k, j, i) =
189serie_lindquist_moins (rayon, distance, xa, ya, za, precision, 30) ;
193 psi_un.set_val_inf (0.5) ;
194 psi_un.std_base_scal() ;
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) ;
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) ;
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) ;
214 psi_deux.set(l, k, j, i) =
215serie_lindquist_plus (rayon, distance, xa, ya, za, precision, 30) ;
218 psi_deux.set_val_inf (0.5) ;
219 psi_deux.std_base_scal() ;
Cmp sqrt(const Cmp &)
Square root.