1
// (C) Copyright John Maddock 2008.
2
// Use, modification and distribution are subject to the
3
// Boost Software License, Version 1.0. (See accompanying file
4
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
7
# pragma warning (disable : 4127) // conditional expression is constant
8
# pragma warning (disable : 4180) // qualifier applied to function type has no meaning; ignored
9
# pragma warning (disable : 4503) // decorated name length exceeded, name was truncated
10
# pragma warning (disable : 4512) // assignment operator could not be generated
11
# pragma warning (disable : 4224) // nonstandard extension used : formal parameter 'function_ptr' was previously defined as a type
14
#include <boost/math/special_functions.hpp>
15
#include <boost/math/tools/roots.hpp>
16
#include <boost/function.hpp>
17
#include <boost/bind.hpp>
22
#include <boost/svg_plot/svg_2d_plot.hpp>
23
#include <boost/svg_plot/show_2d_settings.hpp>
25
class function_arity1_plotter
28
function_arity1_plotter() : m_min_x(0), m_max_x(0), m_min_y(0), m_max_y(0), m_has_legend(false) {}
30
void add(boost::function<double(double)> f, double a, double b, const std::string& name)
35
// Now set our x-axis limits:
37
if(m_max_x == m_min_x)
49
m_points.push_back(std::pair<std::string, std::map<double,double> >(name, std::map<double,double>()));
50
std::map<double,double>& points = m_points.rbegin()->second;
51
double interval = (b - a) / 200;
52
for(double x = a; x <= b; x += interval)
55
// Evaluate the function, set the Y axis limits
56
// if needed and then store the pair of points:
59
if((m_min_y == m_max_y) && (m_min_y == 0))
60
m_min_y = m_max_y = y;
69
void add(const std::map<double, double>& m, const std::string& name)
73
m_points.push_back(std::pair<std::string, std::map<double,double> >(name, m));
74
std::map<double, double>::const_iterator i = m.begin();
78
if((m_min_x == m_min_y) && (m_min_y == 0))
80
m_min_x = m_max_x = i->first;
82
if(i->first < m_min_x)
86
if(i->first > m_max_x)
91
if((m_min_y == m_max_y) && (m_min_y == 0))
93
m_min_y = m_max_y = i->second;
95
if(i->second < m_min_y)
99
if(i->second > m_max_y)
108
void plot(const std::string& title, const std::string& file,
109
const std::string& x_lable = std::string(),
110
const std::string& y_lable = std::string())
112
using namespace boost::svg;
114
static const svg_color colors[5] =
124
plot.image_x_size(600);
125
plot.image_y_size(400);
126
plot.copyright_holder("John Maddock").copyright_date("2008").boost_license_on(true);
127
plot.coord_precision(4); // Could be 3 for smaller plots?
128
plot.title(title).title_font_size(20).title_on(true);
129
plot.legend_on(m_has_legend);
131
double x_delta = (m_max_x - m_min_x) / 50;
132
double y_delta = (m_max_y - m_min_y) / 50;
133
plot.x_range(m_min_x, m_max_x + x_delta)
134
.y_range(m_min_y, m_max_y + y_delta);
135
plot.x_label_on(true).x_label(x_lable);
136
plot.y_label_on(true).y_label(y_lable);
137
plot.y_major_grid_on(false).x_major_grid_on(false);
138
plot.x_num_minor_ticks(3);
139
plot.y_num_minor_ticks(3);
141
// Work out axis tick intervals:
143
double l = std::floor(std::log10((m_max_x - m_min_x) / 10) + 0.5);
144
double interval = std::pow(10.0, (int)l);
145
if(((m_max_x - m_min_x) / interval) > 10)
147
plot.x_major_interval(interval);
148
l = std::floor(std::log10((m_max_y - m_min_y) / 10) + 0.5);
149
interval = std::pow(10.0, (int)l);
150
if(((m_max_y - m_min_y) / interval) > 10)
152
plot.y_major_interval(interval);
153
plot.plot_window_on(true);
154
plot.plot_border_color(lightslategray)
155
.background_border_color(lightslategray)
156
.legend_border_color(lightslategray)
157
.legend_background_color(white);
161
for(std::list<std::pair<std::string, std::map<double,double> > >::const_iterator i = m_points.begin();
162
i != m_points.end(); ++i)
164
plot.plot(i->second, i->first)
166
.line_color(colors[color_index])
171
color_index = color_index % (sizeof(colors)/sizeof(colors[0]));
179
m_min_x = m_min_y = m_max_x = m_max_y = 0;
180
m_has_legend = false;
184
std::list<std::pair<std::string, std::map<double, double> > > m_points;
185
double m_min_x, m_max_x, m_min_y, m_max_y;
190
struct location_finder
192
location_finder(F _f, double t, double x0) : f(_f), target(t), x_off(x0){}
194
double operator()(double x)
198
return f(x + x_off) - target;
200
catch(const std::overflow_error&)
202
return boost::math::tools::max_value<double>();
204
catch(const std::domain_error&)
206
if(x + x_off == x_off)
207
return f(x_off + boost::math::tools::epsilon<double>() * x_off);
219
double find_end_point(F f, double x0, double target, bool rising, double x_off = 0)
221
boost::math::tools::eps_tolerance<double> tol(50);
222
boost::uintmax_t max_iter = 1000;
223
return boost::math::tools::bracket_and_solve_root(
224
location_finder<F>(f, target, x_off),
232
double sqrt1pm1(double x)
234
return boost::math::sqrt1pm1(x);
237
double lbeta(double a, double b)
239
return std::log(boost::math::beta(a, b));
244
function_arity1_plotter plot;
246
double (*f2)(double, double);
247
double (*f2u)(unsigned, double);
248
double (*f2i)(int, double);
249
double (*f3)(double, double, double);
250
double (*f4)(double, double, double, double);
252
f = boost::math::zeta;
253
plot.add(f, 1 + find_end_point(f, 0.1, 40.0, false, 1.0), 10, "");
254
plot.add(f, -20, 1 + find_end_point(f, -0.1, -40.0, false, 1.0), "");
255
plot.plot("Zeta Function Over [-20,10]", "zeta1.svg", "z", "zeta(z)");
258
plot.add(f, -14, 0, "");
259
plot.plot("Zeta Function Over [-14,0]", "zeta2.svg", "z", "zeta(z)");
261
f = boost::math::tgamma;
262
double max_val = f(6);
264
plot.add(f, find_end_point(f, 0.1, max_val, false), 6, "");
265
plot.add(f, -1 + find_end_point(f, 0.1, -max_val, true, -1), find_end_point(f, -0.1, -max_val, false), "");
266
plot.add(f, -2 + find_end_point(f, 0.1, max_val, false, -2), -1 + find_end_point(f, -0.1, max_val, true, -1), "");
267
plot.add(f, -3 + find_end_point(f, 0.1, -max_val, true, -3), -2 + find_end_point(f, -0.1, -max_val, false, -2), "");
268
plot.add(f, -4 + find_end_point(f, 0.1, max_val, false, -4), -3 + find_end_point(f, -0.1, max_val, true, -3), "");
269
plot.plot("tgamma", "tgamma.svg", "z", "tgamma(z)");
271
f = boost::math::lgamma;
274
plot.add(f, find_end_point(f, 0.1, max_val, false), 10, "");
275
plot.add(f, -1 + find_end_point(f, 0.1, max_val, false, -1), find_end_point(f, -0.1, max_val, true), "");
276
plot.add(f, -2 + find_end_point(f, 0.1, max_val, false, -2), -1 + find_end_point(f, -0.1, max_val, true, -1), "");
277
plot.add(f, -3 + find_end_point(f, 0.1, max_val, false, -3), -2 + find_end_point(f, -0.1, max_val, true, -2), "");
278
plot.add(f, -4 + find_end_point(f, 0.1, max_val, false, -4), -3 + find_end_point(f, -0.1, max_val, true, -3), "");
279
plot.add(f, -5 + find_end_point(f, 0.1, max_val, false, -5), -4 + find_end_point(f, -0.1, max_val, true, -4), "");
280
plot.plot("lgamma", "lgamma.svg", "z", "lgamma(z)");
282
f = boost::math::digamma;
285
plot.add(f, find_end_point(f, 0.1, -max_val, true), 10, "");
286
plot.add(f, -1 + find_end_point(f, 0.1, -max_val, true, -1), find_end_point(f, -0.1, max_val, true), "");
287
plot.add(f, -2 + find_end_point(f, 0.1, -max_val, true, -2), -1 + find_end_point(f, -0.1, max_val, true, -1), "");
288
plot.add(f, -3 + find_end_point(f, 0.1, -max_val, true, -3), -2 + find_end_point(f, -0.1, max_val, true, -2), "");
289
plot.add(f, -4 + find_end_point(f, 0.1, -max_val, true, -4), -3 + find_end_point(f, -0.1, max_val, true, -3), "");
290
plot.plot("digamma", "digamma.svg", "z", "digamma(z)");
292
f = boost::math::erf;
294
plot.add(f, -3, 3, "");
295
plot.plot("erf", "erf.svg", "z", "erf(z)");
296
f = boost::math::erfc;
298
plot.add(f, -3, 3, "");
299
plot.plot("erfc", "erfc.svg", "z", "erfc(z)");
301
f = boost::math::erf_inv;
303
plot.add(f, -1 + find_end_point(f, 0.1, -3, true, -1), 1 + find_end_point(f, -0.1, 3, true, 1), "");
304
plot.plot("erf_inv", "erf_inv.svg", "z", "erf_inv(z)");
305
f = boost::math::erfc_inv;
307
plot.add(f, find_end_point(f, 0.1, 3, false), 2 + find_end_point(f, -0.1, -3, false, 2), "");
308
plot.plot("erfc_inv", "erfc_inv.svg", "z", "erfc_inv(z)");
310
f = boost::math::log1p;
312
plot.add(f, -1 + find_end_point(f, 0.1, -10, true, -1), 10, "");
313
plot.plot("log1p", "log1p.svg", "z", "log1p(z)");
315
f = boost::math::expm1;
317
plot.add(f, -4, 2, "");
318
plot.plot("expm1", "expm1.svg", "z", "expm1(z)");
320
f = boost::math::cbrt;
322
plot.add(f, -10, 10, "");
323
plot.plot("cbrt", "cbrt.svg", "z", "cbrt(z)");
327
plot.add(f, -1 + find_end_point(f, 0.1, -10, true, -1), 5, "");
328
plot.plot("sqrt1pm1", "sqrt1pm1.svg", "z", "sqrt1pm1(z)");
330
f2 = boost::math::powm1;
332
plot.add(boost::bind(f2, 0.0001, _1), find_end_point(boost::bind(f2, 0.0001, _1), -1, 10, false), 5, "a=0.0001");
333
plot.add(boost::bind(f2, 0.001, _1), find_end_point(boost::bind(f2, 0.001, _1), -1, 10, false), 5, "a=0.001");
334
plot.add(boost::bind(f2, 0.01, _1), find_end_point(boost::bind(f2, 0.01, _1), -1, 10, false), 5, "a=0.01");
335
plot.add(boost::bind(f2, 0.1, _1), find_end_point(boost::bind(f2, 0.1, _1), -1, 10, false), 5, "a=0.1");
336
plot.add(boost::bind(f2, 0.75, _1), -5, 5, "a=0.75");
337
plot.add(boost::bind(f2, 1.25, _1), -5, 5, "a=1.25");
338
plot.plot("powm1", "powm1.svg", "z", "powm1(a, z)");
340
f = boost::math::sinc_pi;
342
plot.add(f, -10, 10, "");
343
plot.plot("sinc_pi", "sinc_pi.svg", "z", "sinc_pi(z)");
345
f = boost::math::sinhc_pi;
347
plot.add(f, -5, 5, "");
348
plot.plot("sinhc_pi", "sinhc_pi.svg", "z", "sinhc_pi(z)");
350
f = boost::math::acosh;
352
plot.add(f, 1, 10, "");
353
plot.plot("acosh", "acosh.svg", "z", "acosh(z)");
355
f = boost::math::asinh;
357
plot.add(f, -10, 10, "");
358
plot.plot("asinh", "asinh.svg", "z", "asinh(z)");
360
f = boost::math::atanh;
362
plot.add(f, -1 + find_end_point(f, 0.1, -5, true, -1), 1 + find_end_point(f, -0.1, 5, true, 1), "");
363
plot.plot("atanh", "atanh.svg", "z", "atanh(z)");
365
f2 = boost::math::tgamma_delta_ratio;
367
plot.add(boost::bind(f2, _1, -0.5), 1, 40, "delta = -0.5");
368
plot.add(boost::bind(f2, _1, -0.2), 1, 40, "delta = -0.2");
369
plot.add(boost::bind(f2, _1, -0.1), 1, 40, "delta = -0.1");
370
plot.add(boost::bind(f2, _1, 0.1), 1, 40, "delta = 0.1");
371
plot.add(boost::bind(f2, _1, 0.2), 1, 40, "delta = 0.2");
372
plot.add(boost::bind(f2, _1, 0.5), 1, 40, "delta = 0.5");
373
plot.add(boost::bind(f2, _1, 1.0), 1, 40, "delta = 1.0");
374
plot.plot("tgamma_delta_ratio", "tgamma_delta_ratio.svg", "z", "tgamma_delta_ratio(delta, z)");
376
f2 = boost::math::gamma_p;
378
plot.add(boost::bind(f2, 0.5, _1), 0, 20, "a = 0.5");
379
plot.add(boost::bind(f2, 1.0, _1), 0, 20, "a = 1.0");
380
plot.add(boost::bind(f2, 5.0, _1), 0, 20, "a = 5.0");
381
plot.add(boost::bind(f2, 10.0, _1), 0, 20, "a = 10.0");
382
plot.plot("gamma_p", "gamma_p.svg", "z", "gamma_p(a, z)");
384
f2 = boost::math::gamma_q;
386
plot.add(boost::bind(f2, 0.5, _1), 0, 20, "a = 0.5");
387
plot.add(boost::bind(f2, 1.0, _1), 0, 20, "a = 1.0");
388
plot.add(boost::bind(f2, 5.0, _1), 0, 20, "a = 5.0");
389
plot.add(boost::bind(f2, 10.0, _1), 0, 20, "a = 10.0");
390
plot.plot("gamma_q", "gamma_q.svg", "z", "gamma_q(a, z)");
394
plot.add(boost::bind(f2, 0.5, _1), 0.00001, 5, "a = 0.5");
395
plot.add(boost::bind(f2, 1.0, _1), 0.00001, 5, "a = 1.0");
396
plot.add(boost::bind(f2, 5.0, _1), 0.00001, 5, "a = 5.0");
397
plot.add(boost::bind(f2, 10.0, _1), 0.00001, 5, "a = 10.0");
398
plot.plot("beta", "beta.svg", "z", "log(beta(a, z))");
400
f = boost::math::expint;
403
plot.add(f, find_end_point(f, 0.1, -max_val, true), 4, "");
404
plot.add(f, -3, find_end_point(f, -0.1, -max_val, false), "");
405
plot.plot("Exponential Integral Ei", "expint_i.svg", "z", "expint(z)");
407
f2u = boost::math::expint;
410
plot.add(boost::bind(f2u, 1, _1), find_end_point(boost::bind(f2u, 1, _1), 0.1, max_val, false), 2, "n = 1 ");
411
plot.add(boost::bind(f2u, 2, _1), find_end_point(boost::bind(f2u, 2, _1), 0.1, max_val, false), 2, "n = 2 ");
412
plot.add(boost::bind(f2u, 3, _1), 0, 2, "n = 3 ");
413
plot.add(boost::bind(f2u, 4, _1), 0, 2, "n = 4 ");
414
plot.plot("Exponential Integral En", "expint2.svg", "z", "expint(n, z)");
416
f3 = boost::math::ibeta;
418
plot.add(boost::bind(f3, 9, 1, _1), 0, 1, "a = 9, b = 1");
419
plot.add(boost::bind(f3, 7, 2, _1), 0, 1, "a = 7, b = 2");
420
plot.add(boost::bind(f3, 5, 5, _1), 0, 1, "a = 5, b = 5");
421
plot.add(boost::bind(f3, 2, 7, _1), 0, 1, "a = 2, b = 7");
422
plot.add(boost::bind(f3, 1, 9, _1), 0, 1, "a = 1, b = 9");
423
plot.plot("ibeta", "ibeta.svg", "z", "ibeta(a, b, z)");
425
f2i = boost::math::legendre_p;
427
plot.add(boost::bind(f2i, 1, _1), -1, 1, "l = 1");
428
plot.add(boost::bind(f2i, 2, _1), -1, 1, "l = 2");
429
plot.add(boost::bind(f2i, 3, _1), -1, 1, "l = 3");
430
plot.add(boost::bind(f2i, 4, _1), -1, 1, "l = 4");
431
plot.add(boost::bind(f2i, 5, _1), -1, 1, "l = 5");
432
plot.plot("Legendre Polynomials", "legendre_p.svg", "x", "legendre_p(l, x)");
434
f2u = boost::math::legendre_q;
436
plot.add(boost::bind(f2u, 1, _1), -0.95, 0.95, "l = 1");
437
plot.add(boost::bind(f2u, 2, _1), -0.95, 0.95, "l = 2");
438
plot.add(boost::bind(f2u, 3, _1), -0.95, 0.95, "l = 3");
439
plot.add(boost::bind(f2u, 4, _1), -0.95, 0.95, "l = 4");
440
plot.add(boost::bind(f2u, 5, _1), -0.95, 0.95, "l = 5");
441
plot.plot("Legendre Polynomials of the Second Kind", "legendre_q.svg", "x", "legendre_q(l, x)");
443
f2u = boost::math::laguerre;
445
plot.add(boost::bind(f2u, 0, _1), -5, 10, "n = 0");
446
plot.add(boost::bind(f2u, 1, _1), -5, 10, "n = 1");
447
plot.add(boost::bind(f2u, 2, _1),
448
find_end_point(boost::bind(f2u, 2, _1), -2, 20, false),
449
find_end_point(boost::bind(f2u, 2, _1), 4, 20, true),
451
plot.add(boost::bind(f2u, 3, _1),
452
find_end_point(boost::bind(f2u, 3, _1), -2, 20, false),
453
8 + find_end_point(boost::bind(f2u, 3, _1), 1, 20, false, 8),
455
plot.add(boost::bind(f2u, 4, _1),
456
find_end_point(boost::bind(f2u, 4, _1), -2, 20, false),
457
8 + find_end_point(boost::bind(f2u, 4, _1), 1, 20, true, 8),
459
plot.add(boost::bind(f2u, 5, _1),
460
find_end_point(boost::bind(f2u, 5, _1), -2, 20, false),
461
8 + find_end_point(boost::bind(f2u, 5, _1), 1, 20, true, 8),
463
plot.plot("Laguerre Polynomials", "laguerre.svg", "x", "laguerre(n, x)");
465
f2u = boost::math::hermite;
467
plot.add(boost::bind(f2u, 0, _1), -1.8, 1.8, "n = 0");
468
plot.add(boost::bind(f2u, 1, _1), -1.8, 1.8, "n = 1");
469
plot.add(boost::bind(f2u, 2, _1), -1.8, 1.8, "n = 2");
470
plot.add(boost::bind(f2u, 3, _1), -1.8, 1.8, "n = 3");
471
plot.add(boost::bind(f2u, 4, _1), -1.8, 1.8, "n = 4");
472
plot.plot("Hermite Polynomials", "hermite.svg", "x", "hermite(n, x)");
474
f2 = boost::math::cyl_bessel_j;
476
plot.add(boost::bind(f2, 0, _1), -20, 20, "v = 0");
477
plot.add(boost::bind(f2, 1, _1), -20, 20, "v = 1");
478
plot.add(boost::bind(f2, 2, _1), -20, 20, "v = 2");
479
plot.add(boost::bind(f2, 3, _1), -20, 20, "v = 3");
480
plot.add(boost::bind(f2, 4, _1), -20, 20, "v = 4");
481
plot.plot("Bessel J", "cyl_bessel_j.svg", "x", "cyl_bessel_j(v, x)");
483
f2 = boost::math::cyl_neumann;
485
plot.add(boost::bind(f2, 0, _1), find_end_point(boost::bind(f2, 0, _1), 0.1, -5, true), 20, "v = 0");
486
plot.add(boost::bind(f2, 1, _1), find_end_point(boost::bind(f2, 1, _1), 0.1, -5, true), 20, "v = 1");
487
plot.add(boost::bind(f2, 2, _1), find_end_point(boost::bind(f2, 2, _1), 0.1, -5, true), 20, "v = 2");
488
plot.add(boost::bind(f2, 3, _1), find_end_point(boost::bind(f2, 3, _1), 0.1, -5, true), 20, "v = 3");
489
plot.add(boost::bind(f2, 4, _1), find_end_point(boost::bind(f2, 4, _1), 0.1, -5, true), 20, "v = 4");
490
plot.plot("Bessel Y", "cyl_neumann.svg", "x", "cyl_neumann(v, x)");
492
f2 = boost::math::cyl_bessel_i;
494
plot.add(boost::bind(f2, 0, _1), find_end_point(boost::bind(f2, 0, _1), -0.1, 20, false), find_end_point(boost::bind(f2, 0, _1), 0.1, 20, true), "v = 0");
495
plot.add(boost::bind(f2, 2, _1), find_end_point(boost::bind(f2, 2, _1), -0.1, 20, false), find_end_point(boost::bind(f2, 2, _1), 0.1, 20, true), "v = 2");
496
plot.add(boost::bind(f2, 5, _1), find_end_point(boost::bind(f2, 5, _1), -0.1, -20, true), find_end_point(boost::bind(f2, 5, _1), 0.1, 20, true), "v = 5");
497
plot.add(boost::bind(f2, 7, _1), find_end_point(boost::bind(f2, 7, _1), -0.1, -20, true), find_end_point(boost::bind(f2, 7, _1), 0.1, 20, true), "v = 7");
498
plot.add(boost::bind(f2, 10, _1), find_end_point(boost::bind(f2, 10, _1), -0.1, 20, false), find_end_point(boost::bind(f2, 10, _1), 0.1, 20, true), "v = 10");
499
plot.plot("Bessel I", "cyl_bessel_i.svg", "x", "cyl_bessel_i(v, x)");
501
f2 = boost::math::cyl_bessel_k;
503
plot.add(boost::bind(f2, 0, _1), find_end_point(boost::bind(f2, 0, _1), 0.1, 10, false), 10, "v = 0");
504
plot.add(boost::bind(f2, 2, _1), find_end_point(boost::bind(f2, 2, _1), 0.1, 10, false), 10, "v = 2");
505
plot.add(boost::bind(f2, 5, _1), find_end_point(boost::bind(f2, 5, _1), 0.1, 10, false), 10, "v = 5");
506
plot.add(boost::bind(f2, 7, _1), find_end_point(boost::bind(f2, 7, _1), 0.1, 10, false), 10, "v = 7");
507
plot.add(boost::bind(f2, 10, _1), find_end_point(boost::bind(f2, 10, _1), 0.1, 10, false), 10, "v = 10");
508
plot.plot("Bessel K", "cyl_bessel_k.svg", "x", "cyl_bessel_k(v, x)");
510
f2u = boost::math::sph_bessel;
512
plot.add(boost::bind(f2u, 0, _1), 0, 20, "v = 0");
513
plot.add(boost::bind(f2u, 2, _1), 0, 20, "v = 2");
514
plot.add(boost::bind(f2u, 5, _1), 0, 20, "v = 5");
515
plot.add(boost::bind(f2u, 7, _1), 0, 20, "v = 7");
516
plot.add(boost::bind(f2u, 10, _1), 0, 20, "v = 10");
517
plot.plot("Bessel j", "sph_bessel.svg", "x", "sph_bessel(v, x)");
519
f2u = boost::math::sph_neumann;
521
plot.add(boost::bind(f2u, 0, _1), find_end_point(boost::bind(f2u, 0, _1), 0.1, -5, true), 20, "v = 0");
522
plot.add(boost::bind(f2u, 2, _1), find_end_point(boost::bind(f2u, 2, _1), 0.1, -5, true), 20, "v = 2");
523
plot.add(boost::bind(f2u, 5, _1), find_end_point(boost::bind(f2u, 5, _1), 0.1, -5, true), 20, "v = 5");
524
plot.add(boost::bind(f2u, 7, _1), find_end_point(boost::bind(f2u, 7, _1), 0.1, -5, true), 20, "v = 7");
525
plot.add(boost::bind(f2u, 10, _1), find_end_point(boost::bind(f2u, 10, _1), 0.1, -5, true), 20, "v = 10");
526
plot.plot("Bessel y", "sph_neumann.svg", "x", "sph_neumann(v, x)");
528
f4 = boost::math::ellint_rj;
530
plot.add(boost::bind(f4, _1, _1, _1, _1), find_end_point(boost::bind(f4, _1, _1, _1, _1), 0.1, 10, false), 4, "RJ");
531
f3 = boost::math::ellint_rf;
532
plot.add(boost::bind(f3, _1, _1, _1), find_end_point(boost::bind(f3, _1, _1, _1), 0.1, 10, false), 4, "RF");
533
plot.plot("Elliptic Integrals", "ellint_carlson.svg", "x", "");
535
f2 = boost::math::ellint_1;
537
plot.add(boost::bind(f2, _1, 0.5), -0.9, 0.9, "φ=0.5");
538
plot.add(boost::bind(f2, _1, 0.75), -0.9, 0.9, "φ=0.75");
539
plot.add(boost::bind(f2, _1, 1.25), -0.9, 0.9, "φ=1.25");
540
plot.add(boost::bind(f2, _1, boost::math::constants::pi<double>() / 2), -0.9, 0.9, "φ=π/2");
541
plot.plot("Elliptic Of the First Kind", "ellint_1.svg", "k", "ellint_1(k, phi)");
543
f2 = boost::math::ellint_2;
545
plot.add(boost::bind(f2, _1, 0.5), -1, 1, "φ=0.5");
546
plot.add(boost::bind(f2, _1, 0.75), -1, 1, "φ=0.75");
547
plot.add(boost::bind(f2, _1, 1.25), -1, 1, "φ=1.25");
548
plot.add(boost::bind(f2, _1, boost::math::constants::pi<double>() / 2), -1, 1, "φ=π/2");
549
plot.plot("Elliptic Of the Second Kind", "ellint_2.svg", "k", "ellint_2(k, phi)");
551
f3 = boost::math::ellint_3;
553
plot.add(boost::bind(f3, _1, 0, 1.25), -1, 1, "n=0 φ=1.25");
554
plot.add(boost::bind(f3, _1, 0.5, 1.25), -1, 1, "n=0.5 φ=1.25");
555
plot.add(boost::bind(f3, _1, 0.25, boost::math::constants::pi<double>() / 2),
557
boost::bind(f3, _1, 0.25, boost::math::constants::pi<double>() / 2),
558
0.5, 4, false, -1) - 1,
560
boost::bind(f3, _1, 0.25, boost::math::constants::pi<double>() / 2),
561
-0.5, 4, true, 1) + 1, "n=0.25 φ=π/2");
562
plot.add(boost::bind(f3, _1, 0.75, boost::math::constants::pi<double>() / 2),
564
boost::bind(f3, _1, 0.75, boost::math::constants::pi<double>() / 2),
565
0.5, 4, false, -1) - 1,
567
boost::bind(f3, _1, 0.75, boost::math::constants::pi<double>() / 2),
568
-0.5, 4, true, 1) + 1, "n=0.75 φ=π/2");
569
plot.plot("Elliptic Of the Third Kind", "ellint_3.svg", "k", "ellint_3(k, n, phi)");