1 /*
2 [auto_generated]
3 libs/numeric/odeint/examples/black_hole.cpp
4
5 [begin_description]
6 This example shows how the __float128 from gcc libquadmath can be used with odeint.
7 [end_description]
8
9 Copyright 2012 Karsten Ahnert
10 Copyright 2012 Lee Hodgkinson
11 Copyright 2012 Mario Mulansky
12
13 Distributed under the Boost Software License, Version 1.0.
14 (See accompanying file LICENSE_1_0.txt or
15 copy at http://www.boost.org/LICENSE_1_0.txt)
16 */
17
18 #include <cstdlib>
19 #include <cmath>
20 #include <iostream>
21 #include <iterator>
22 #include <utility>
23 #include <algorithm>
24 #include <cassert>
25 #include <vector>
26 #include <complex>
27
28 extern "C" {
29 #include <quadmath.h>
30 }
31
32 const __float128 zero =strtoflt128 ("0.0", NULL);
33
34 namespace std {
35
abs(__float128 x)36 inline __float128 abs( __float128 x )
37 {
38 return fabsq( x );
39 }
40
sqrt(__float128 x)41 inline __float128 sqrt( __float128 x )
42 {
43 return sqrtq( x );
44 }
45
pow(__float128 x,__float128 y)46 inline __float128 pow( __float128 x , __float128 y )
47 {
48 return powq( x , y );
49 }
50
abs(std::complex<__float128> x)51 inline __float128 abs( std::complex< __float128 > x )
52 {
53 return sqrtq( x.real() * x.real() + x.imag() * x.imag() );
54 }
55
pow(std::complex<__float128> x,__float128 y)56 inline std::complex< __float128 > pow( std::complex< __float128> x , __float128 y )
57 {
58 __float128 r = pow( abs(x) , y );
59 __float128 phi = atanq( x.imag() / x.real() );
60 return std::complex< __float128 >( r * cosq( y * phi ) , r * sinq( y * phi ) );
61 }
62 }
63
operator <<(std::ostream & os,const __float128 & f)64 inline std::ostream& operator<< (std::ostream& os, const __float128& f) {
65
66 char* y = new char[1000];
67 quadmath_snprintf(y, 1000, "%.30Qg", f) ;
68 os.precision(30);
69 os<<y;
70 delete[] y;
71 return os;
72 }
73
74
75 #include <boost/array.hpp>
76 #include <boost/range/algorithm.hpp>
77 #include <boost/range/adaptor/filtered.hpp>
78 #include <boost/range/numeric.hpp>
79 #include <boost/numeric/odeint.hpp>
80
81
82
83 using namespace boost::numeric::odeint;
84 using namespace std;
85
86 typedef __float128 my_float;
87 typedef std::vector< std::complex < my_float > > state_type;
88
89 struct radMod
90 {
91 my_float m_om;
92 my_float m_l;
93
radModradMod94 radMod( my_float om , my_float l )
95 : m_om( om ) , m_l( l ) { }
96
operator ()radMod97 void operator()( const state_type &x , state_type &dxdt , my_float r ) const
98 {
99
100 dxdt[0] = x[1];
101 dxdt[1] = -(2*(r-1)/(r*(r-2)))*x[1]-((m_om*m_om*r*r/((r-2)*(r-2)))-(m_l*(m_l+1)/(r*(r-2))))*x[0];
102 }
103 };
104
105
106
107
108
109
110
main(int argc,char ** argv)111 int main( int argc , char **argv )
112 {
113
114
115 state_type x(2);
116
117 my_float re0 = strtoflt128( "-0.00008944230755601224204687038354994353820468" , NULL );
118 my_float im0 = strtoflt128( "0.00004472229441850588228136889483397204368247" , NULL );
119 my_float re1 = strtoflt128( "-4.464175354293244250869336196695966076150E-6 " , NULL );
120 my_float im1 = strtoflt128( "-8.950483248390306670770345406051469584488E-6" , NULL );
121
122 x[0] = complex< my_float >( re0 ,im0 );
123 x[1] = complex< my_float >( re1 ,im1 );
124
125 const my_float dt =strtoflt128 ("-0.001", NULL);
126 const my_float start =strtoflt128 ("10000.0", NULL);
127 const my_float end =strtoflt128 ("9990.0", NULL);
128 const my_float omega =strtoflt128 ("2.0", NULL);
129 const my_float ell =strtoflt128 ("1.0", NULL);
130
131
132
133 my_float abs_err = strtoflt128( "1.0E-15" , NULL ) , rel_err = strtoflt128( "1.0E-10" , NULL );
134 my_float a_x = strtoflt128( "1.0" , NULL ) , a_dxdt = strtoflt128( "1.0" , NULL );
135
136 typedef runge_kutta_dopri5< state_type, my_float > dopri5_type;
137 typedef controlled_runge_kutta< dopri5_type > controlled_dopri5_type;
138 typedef dense_output_runge_kutta< controlled_dopri5_type > dense_output_dopri5_type;
139
140 dense_output_dopri5_type dopri5( controlled_dopri5_type( default_error_checker< my_float >( abs_err , rel_err , a_x , a_dxdt ) ) );
141
142 std::for_each( make_adaptive_time_iterator_begin(dopri5 , radMod(omega , ell) , x , start , end , dt) ,
143 make_adaptive_time_iterator_end(dopri5 , radMod(omega , ell) , x ) ,
144 []( const std::pair< state_type&, my_float > &x ) {
145 std::cout << x.second << ", " << x.first[0].real() << "\n"; }
146 );
147
148
149
150 return 0;
151 }
152