LORENE
et_bin_bhns_extr_max.C
1/*
2 * Method of class Et_bin_bhns_extr to search the position of the maximum
3 * enthalpy
4 *
5 * (see file et_bin_bhns_extr.h for documentation).
6 *
7 */
8
9/*
10 * Copyright (c) 2004 Keisuke Taniguchi
11 *
12 * This file is part of LORENE.
13 *
14 * LORENE is free software; you can redistribute it and/or modify
15 * it under the terms of the GNU General Public License version 2
16 * as published by the Free Software Foundation.
17 *
18 * LORENE is distributed in the hope that it will be useful,
19 * but WITHOUT ANY WARRANTY; without even the implied warranty of
20 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 * GNU General Public License for more details.
22 *
23 * You should have received a copy of the GNU General Public License
24 * along with LORENE; if not, write to the Free Software
25 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
26 *
27 */
28
29char et_bin_bhns_extr_max_C[] = "$Header: /cvsroot/Lorene/C++/Source/Etoile/et_bin_bhns_extr_max.C,v 1.3 2014/10/13 08:52:55 j_novak Exp $" ;
30
31/*
32 * $Id: et_bin_bhns_extr_max.C,v 1.3 2014/10/13 08:52:55 j_novak Exp $
33 * $Log: et_bin_bhns_extr_max.C,v $
34 * Revision 1.3 2014/10/13 08:52:55 j_novak
35 * Lorene classes and functions now belong to the namespace Lorene.
36 *
37 * Revision 1.2 2014/10/06 15:13:07 j_novak
38 * Modified #include directives to use c++ syntax.
39 *
40 * Revision 1.1 2004/11/30 20:50:24 k_taniguchi
41 * *** empty log message ***
42 *
43 *
44 * $Header: /cvsroot/Lorene/C++/Source/Etoile/et_bin_bhns_extr_max.C,v 1.3 2014/10/13 08:52:55 j_novak Exp $
45 *
46 */
47
48// C headers
49#include <cmath>
50
51// Lorene headers
52#include "et_bin_bhns_extr.h"
53//#include "utilitaires.h"
54
55namespace Lorene {
56void Et_bin_bhns_extr::ent_max_search(double& xx, double& yy) const {
57
58 //------------------------------------------------------------------//
59 // Calculate the derivative of the enthalpy field //
60 //------------------------------------------------------------------//
61
62 const Tenseur& dent = ent.gradient() ;
63
64 double xxp = 0. ; // Position of the x-coordinate, initialized to zero
65 double yyp = 0. ; // Position of the y-coordinate, initialized to zero
66 double xtmp, ytmp ;
67 int mm, nn ; // Number of steps to the x and y directions
68 double rr = 0. ; // r coordinate, initialized to zero
69 double pp = M_PI/2. ; // phi coordinate, initialized to M_PI/2.
70 double dval_x ; // Direction of dent(0) (1 or -1)
71 double dval_y ; // Direction of dent(1) (1 or -1)
72 double ss ;
73
74 while ( fabs(dent(0).val_point(rr, M_PI/2., pp)) > 1.e-15 ||
75 fabs(dent(1).val_point(rr, M_PI/2., pp)) > 5.e-15) {
76
77 double dx = 1. ; // Step interval to the x-direction, initialized to 1
78 double dy = 1. ; // Step interval to the y-direction, initialized to 1
79 double diff_dent_x ;
80 double diff_dent_y ;
81
82 while ( dy > 1.e-15 ) {
83
84 diff_dent_y = 1. ;
85 nn = 0 ;
86 dy = 0.1 * dy ;
87
88 rr = sqrt( xxp*xxp + yyp*yyp ) ;
89
90 if ( xxp == 0. ) {
91 pp = M_PI/2. ; // There is a possibility of (pp = 1.5*M_PI)
92 }
93 else {
94 pp = acos( xxp / rr ) ;
95 }
96
97 dval_y = dent(1).val_point(rr, M_PI/2., pp) ;
98
99 if ( dval_y > 0. ) {
100 ss = 1. ;
101 }
102 else {
103 ss = -1. ;
104 }
105
106 while( diff_dent_y > 1.e-15 ) {
107
108 nn++ ;
109 ytmp = yyp + ss * nn * dy ;
110
111 rr = sqrt( xxp*xxp + ytmp*ytmp ) ;
112
113 if ( xxp == 0. ) {
114 if ( ss > 0. ) {
115 pp = M_PI/2. ;
116 }
117 else {
118 pp = 1.5*M_PI ;
119 }
120 }
121 else {
122 pp = acos( xxp / rr ) ;
123 }
124
125 diff_dent_y = ss * dent(1).val_point(rr, M_PI/2., pp) ;
126
127 }
128 yyp += ss * (nn - 1) * dy ;
129
130 }
131
132 while ( dx > 1.e-15 ) {
133
134 diff_dent_x = 1. ;
135 mm = 0 ;
136 dx = 0.1 * dx ;
137
138 rr = sqrt( xxp*xxp + yyp*yyp ) ;
139
140 if ( xxp == 0. ) {
141 pp = M_PI/2. ; // There is a possibility of (pp = 1.5*M_PI)
142 }
143 else {
144 pp = acos( xxp / rr ) ;
145 }
146
147 dval_x = dent(0).val_point(rr, M_PI/2., pp) ;
148
149 if ( dval_x > 0. ) {
150 ss = 1. ;
151 }
152 else {
153 ss = -1. ;
154 }
155
156 while( diff_dent_x > 1.e-15 ) {
157
158 mm++ ;
159 xtmp = xxp + ss * mm * dx ;
160
161 rr = sqrt( xtmp*xtmp + yyp*yyp ) ;
162
163 if ( xtmp == 0. ) {
164 if ( ss > 0. ) {
165 pp = M_PI/2. ;
166 }
167 else {
168 pp = 1.5*M_PI ;
169 }
170 }
171 else {
172 pp = acos( xtmp / rr ) ;
173 }
174
175 diff_dent_x = ss * dent(0).val_point(rr, M_PI/2., pp) ;
176
177 }
178 xxp += ss * (mm - 1) * dx ;
179
180 }
181 }
182
183 xx = xxp ;
184 yy = yyp ;
185
186}
187}
void ent_max_search(double &xx, double &yy) const
Searches the position of the maximum enthalpy.
Tenseur ent
Log-enthalpy (relativistic case) or specific enthalpy (Newtonian case)
Definition etoile.h:457
Tensor handling *** DEPRECATED : use class Tensor instead ***.
Definition tenseur.h:301
const Tenseur & gradient() const
Returns the gradient of *this (Cartesian coordinates)
Definition tenseur.C:1542
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:220
Cmp acos(const Cmp &)
Arccosine.
Definition cmp_math.C:169
Lorene prototypes.
Definition app_hor.h:64