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 
23 char 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 
61 namespace Lorene {
62 double 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 ;
89 return result ;
90 }
91 
92 double 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 ;
119 return result ;
120 }
121 
122 double 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  }
138 return result ;
139 }
140 
141 double 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  }
159 return result ;
160 }
161 
162 void 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) =
189 serie_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) =
215 serie_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 }
Lorene
Lorene prototypes.
Definition: app_hor.h:64
Lorene::sqrt
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:220