LORENE
regularise_shift.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 regularise_shift_C[] = "$Header: /cvsroot/Lorene/C++/Source/Bhole_binaire/regularise_shift.C,v 1.7 2014/10/13 08:52:41 j_novak Exp $" ;
24
25/*
26 * $Id: regularise_shift.C,v 1.7 2014/10/13 08:52:41 j_novak Exp $
27 * $Log: regularise_shift.C,v $
28 * Revision 1.7 2014/10/13 08:52:41 j_novak
29 * Lorene classes and functions now belong to the namespace Lorene.
30 *
31 * Revision 1.6 2014/10/06 15:12:59 j_novak
32 * Modified #include directives to use c++ syntax.
33 *
34 * Revision 1.5 2006/04/27 09:12:32 p_grandclement
35 * First try at irrotational black holes
36 *
37 * Revision 1.4 2005/08/29 15:10:14 p_grandclement
38 * Addition of things needed :
39 * 1) For BBH with different masses
40 * 2) Provisory files for the mixted binaries (Bh and NS) : THIS IS NOT
41 * WORKING YET !!!
42 *
43 * Revision 1.3 2003/10/03 15:58:44 j_novak
44 * Cleaning of some headers
45 *
46 * Revision 1.2 2003/02/13 16:40:25 p_grandclement
47 * Addition of various things for the Bin_ns_bh project, non of them being
48 * completely tested
49 *
50 * Revision 1.1.1.1 2001/11/20 15:19:28 e_gourgoulhon
51 * LORENE
52 *
53 * Revision 2.6 2001/05/07 09:11:48 phil
54 * *** empty log message ***
55 *
56 * Revision 2.5 2001/01/29 14:30:48 phil
57 * ajout type rotation
58 *
59 * Revision 2.4 2000/11/02 10:18:05 phil
60 * modification du degre de l;a fonction\
61 * de correction. On passe de 2 a 3
62 *
63 * Revision 2.3 2000/10/31 10:04:28 phil
64 * correction importation
65 *
66 * Revision 2.2 2000/10/26 12:34:17 phil
67 * *** empty log message ***
68 *
69 * Revision 2.1 2000/10/26 12:31:55 phil
70 * correction orientation pour calcul du shift total
71 *
72 * Revision 2.0 2000/10/19 10:08:03 phil
73 * *** empty log message ***
74 *
75 *
76 * $Header: /cvsroot/Lorene/C++/Source/Bhole_binaire/regularise_shift.C,v 1.7 2014/10/13 08:52:41 j_novak Exp $
77 *
78 */
79
80
81//Standard
82#include <cstdlib>
83#include <cmath>
84
85//Lorene
86#include "nbr_spx.h"
87#include "tenseur.h"
88
89namespace Lorene {
90double regle (Tenseur& shift_auto, const Tenseur& shift_comp, double omega, double omega_local) {
91
92 Tenseur shift_old (shift_auto) ;
93
94 double orientation = shift_auto.get_mp()->get_rot_phi() ;
95 assert ((orientation==0) || (orientation == M_PI)) ;
96 double orientation_autre = shift_comp.get_mp()->get_rot_phi() ;
97 assert ((orientation_autre==0) || (orientation_autre == M_PI)) ;
98
99 int alignes = (orientation == orientation_autre) ? 1 : -1 ;
100
101 // Cas triades identiques
102 if (*shift_comp.get_triad() == *shift_auto.get_triad())
103 alignes = 1 ;
104
105 int nz = shift_auto.get_mp()->get_mg()->get_nzone() ;
106 int np = shift_auto.get_mp()->get_mg()->get_np(1) ;
107 int nt = shift_auto.get_mp()->get_mg()->get_nt(1) ;
108 int nr = shift_auto.get_mp()->get_mg()->get_nr(1) ;
109
110 // On minimise la valeur de la derivee de B sur R :
111 Tenseur shift_tot (*shift_auto.get_mp(), 1, CON, *shift_auto.get_triad()) ;
112 shift_tot.set_etat_qcq() ;
113 shift_tot.set(0).import_asymy (alignes*shift_comp(0)) ;
114 shift_tot.set(1).import_symy (alignes*shift_comp(1)) ;
115 shift_tot.set(2).import_asymy (shift_comp(2)) ;
116
117 shift_tot = shift_tot + shift_auto ;
118
119 double indic = (orientation == 0) ? 1 : -1 ;
120
121 Mtbl xa_mtbl (shift_tot.get_mp()->get_mg()) ;
122 xa_mtbl = shift_tot.get_mp()->xa ;
123 Mtbl ya_mtbl (shift_tot.get_mp()->get_mg()) ;
124 ya_mtbl = shift_tot.get_mp()->ya ;
125
126 Tenseur tbi (shift_tot) ;
127 if (omega != 0) {
128 for (int i=0 ; i<3 ; i++) {
129 tbi.set(i).va.coef_i() ;
130 tbi.set(i).va.set_etat_c_qcq() ;
131 }
132
133 tbi.set(0) = *shift_tot(0).va.c - indic *omega * ya_mtbl(0,0,0,0) - indic*omega_local* shift_tot.get_mp()->y ;
134 tbi.set(1) = *shift_tot(1).va.c + indic *omega * xa_mtbl(0,0,0,0) + indic*omega_local* shift_tot.get_mp()->x ;
135 tbi.set_std_base() ;
136 tbi.set(0).annule(nz-1) ;
137 tbi.set(1).annule(nz-1) ;
138 }
139
140 Tenseur derive_r (*shift_auto.get_mp(), 1, CON, *shift_auto.get_triad()) ;
141 derive_r.set_etat_qcq() ;
142 for (int i=0 ; i<3 ; i++) {
143 if (tbi(i).get_etat() != ETATZERO)
144 derive_r.set(i) = tbi(i).dsdr() ;
145 else
146 derive_r.set(i).set_etat_zero() ;
147 }
148
149 // On enleve un fonction pour rendre Kij regulier !
150
151 Valeur val_hor (shift_auto.get_mp()->get_mg()) ;
152 Valeur fonction_radiale (shift_auto.get_mp()->get_mg()) ;
153 Cmp enleve (shift_auto.get_mp()) ;
154
155 double erreur = 0 ;
156 for (int comp=0 ; comp<3 ; comp++)
157 {
158 val_hor.annule_hard() ; // Pour initialiser les tableaux
159 for (int k=0 ; k<np ; k++)
160 for (int j=0 ; j<nt ; j++)
161 for (int i=0 ; i<nr ; i++)
162 val_hor.set(1, k, j, i) = derive_r(comp) (1, k, j, 0) ;
163
164 double r_0 = shift_auto.get_mp()->val_r (1, -1, 0, 0) ;
165 double r_1 = shift_auto.get_mp()->val_r (1, 1, 0, 0) ;
166
167 fonction_radiale = pow(r_1-shift_auto.get_mp()->r, 3.)*
168 (shift_auto.get_mp()->r-r_0)/pow(r_1-r_0, 3.) ;
169 fonction_radiale.annule(0) ;
170 fonction_radiale.annule(2, nz-1) ;
171
172 enleve = fonction_radiale * val_hor ;
173 enleve.va.base = shift_auto(comp).va.base ;
174
175 if (norme(enleve)(1) != 0)
176 shift_auto.set(comp) = shift_auto(comp) - enleve ;
177 if (norme(shift_auto(comp))(1) > 1e-5) {
178 Tbl diff (diffrelmax (shift_auto(comp), shift_old(comp))) ;
179 if (erreur < diff(1))
180 erreur = diff(1) ;
181 }
182 }
183 return erreur ;
184}
185}
void annule(int l)
Sets the Cmp to zero in a given domain.
Definition cmp.C:348
Tbl norme(const Cmp &)
Sums of the absolute values of all the values of the Cmp in each domain.
Definition cmp_math.C:481
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
Tbl diffrelmax(const Cmp &a, const Cmp &b)
Relative difference between two Cmp (max version).
Definition cmp_math.C:539
Lorene prototypes.
Definition app_hor.h:64