KADATH
zone_cut.cpp
1 /*
2  Copyright 2017 Philippe Grandclement
3 
4  This file is part of Kadath.
5 
6  Kadath 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 3 of the License, or
9  (at your option) any later version.
10 
11  Kadath 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 Kadath. If not, see <http://www.gnu.org/licenses/>.
18 */
19 
20 #include "headcpp.hpp"
21 #include "utilities.hpp"
22 #include "param.hpp"
23 
24 namespace Kadath {
25 double feta_zero (double eta, const Param& par) {
26  double chi = par.get_double(0) ;
27  double rsa = par.get_double(1) ;
28  double denom = (cosh(eta)-cos(chi)) ;
29  return ((sinh(eta)*sinh(eta)+sin(chi)*sin(chi))/denom/denom-rsa*rsa) ;
30 }
31 
32 
33 double eta_lim_chi (double chi, double rext, double a, double eta_c) {
34  double res ;
35  Param par_func ;
36  par_func.add_double(chi, 0) ;
37  par_func.add_double(rext/a,1) ;
38  double precis = PRECISION ;
39  int nitermax = 500 ;
40  int niter ;
41 
42  if (chi<=PRECISION)
43  res = eta_c ;
44  else
45  res = zerosec(feta_zero, par_func, 0, eta_c, precis, nitermax, niter) ;
46  return res ;
47 }
48 
49 double fchi_zero (double chi, const Param& par) {
50  double eta = par.get_double(0) ;
51  double rsa = par.get_double(1) ;
52  double denom = (cosh(eta)-cos(chi)) ;
53  return ((sinh(eta)*sinh(eta)+sin(chi)*sin(chi))/denom/denom-rsa*rsa) ;
54 }
55 
56 
57 double chi_lim_eta (double eta, double rext, double a, double chi_c) {
58  double res ;
59  Param par_func ;
60  par_func.add_double(eta, 0) ;
61  par_func.add_double(rext/a,1) ;
62  double precis = PRECISION ;
63  int nitermax = 500 ;
64  int niter ;
65  if (fabs(eta)<=PRECISION)
66  res = chi_c ;
67  else
68  res = zerosec(fchi_zero, par_func, 0, chi_c, precis, nitermax, niter) ;
69  return res ;
70 }}
71