/*
* plot_algorithms.h -- plot function templates, for build use only
*
* This file is part of ePiX, a C++ library for creating high-quality
* figures in LaTeX
*
* Andrew D. Hwang <rot 13 nujnat at zngupf dot ubylpebff dot rqh>
*
* Version 1.1.9
* Last Change: August 01, 2007
*/
/*
* Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2007
* Andrew D. Hwang <rot 13 nujnat at zngupf dot ubylpebff dot rqh>
* Department of Mathematics and Computer Science
* College of the Holy Cross
* Worcester, MA, 01610-2395, USA
*/
/*
* ePiX is free software; you can redistribute it and/or modify it
* under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* ePiX is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
* License for more details.
*
* You should have received a copy of the GNU General Public License
* along with ePiX; if not, write to the Free Software Foundation, Inc.,
* 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*/
/*
* This file provides plotting templates, with T a class
*
* void plot_map_dom(const T&, const domain&)
* void plot_map_domlist(const T&, const domain_list&)
* void plot_function(const T&, double, double, int)
* void euler_plot(const T&, const P&, double, double, int)
* P euler_flow(const T&, P start, double, int)
*/
#ifndef EPIX_PLOT_ALGO
#define EPIX_PLOT_ALGO
#include <cmath>
#include <list>
#include "constants.h"
#include "triples.h"
#include "functions.h"
#include "edge_data.h"
#include "paint_style.h"
#include "domain.h"
#include "path_data.h"
namespace ePiX {
typedef std::list<domain>::const_iterator dolci;
const double INF(EPIX_INFTY);
// Are args endpoints of a plottable edge?
static bool is_valid(const P& tail, const P& head)
{
return tail.is_valid() && head.is_valid();
}
template<class T> void plot_map_dom(const T& map, const domain& R)
{
// only attempt to fill if plotting a curve
const bool is_curve(R.dim() == 1);
const bool filling(is_curve && the_paint_style().fill_flag());
// max summation indices
unsigned int i_max((R.dx1() > 0) ? R.coarse_n1() : 0);
unsigned int j_max((R.dx2() > 0) ? R.coarse_n2() : 0);
unsigned int k_max((R.dx3() > 0) ? R.coarse_n3() : 0);
if (R.dx1() > 0) // there's something to draw
{
// number of intervals in subdivision
unsigned int count1(R.fine_n1());
const double st1(R.dx1());
const double st2(R.step2());
const double st3(R.step3());
for (unsigned int j=0; j <= j_max; ++j)
for (unsigned int k=0; k <= k_max; ++k)
{
std::list<edge3d> data1;
const P init(R.corner1() + (j*st2)*E_2 + (k*st3)*E_3);
for (unsigned int i=0; i < count1; ++i)
{
P curr(map(init + P(i*st1,0,0)));
P next(map(init + P((i+1)*st1,0,0)));
data1.push_back(edge3d(curr, next, is_valid(curr, next)));
}
// close path if initial and terminal points coincide
bool closed(is_curve &&
map(init) == map(init + P(count1*st1,0,0)));
path_data path1(data1, closed && filling, filling);
path1.draw();
}
}
if (R.dx2() > 0)
{
unsigned int count2(R.fine_n2());
const double st1(R.step1());
const double st2(R.dx2());
const double st3(R.step3());
for (unsigned int i=0; i <= i_max; ++i)
for (unsigned int k=0; k <= k_max; ++k)
{
std::list<edge3d> data2;
const P init(R.corner1() + (i*st1)*E_1 + (k*st3)*E_3);
for (unsigned int j=0; j < count2; ++j)
{
P curr(map(init + P(0,j*st2,0)));
P next(map(init + P(0,(j+1)*st2,0)));
data2.push_back(edge3d(curr, next, is_valid(curr, next)));
}
// close path if initial and terminal points coincide
bool closed(is_curve &&
map(init) == map(init + P(0,count2*st2,0)));
path_data path2(data2, closed && filling, filling);
path2.draw();
}
}
if (R.dx3() > 0)
{
unsigned int count3(R.fine_n3());
const double st1(R.step1());
const double st2(R.step2());
const double st3(R.dx3());
for (unsigned int i=0; i <= i_max; ++i)
for (unsigned int j=0; j <= j_max; ++j)
{
std::list<edge3d> data3;
const P init(R.corner1() + (i*st1)*E_1 + (j*st2)*E_2);
for (unsigned int k=0; k < count3; ++k)
{
P curr(map(init + P(0,0,k*st3)));
P next(map(init + P(0,0,(k+1)*st3)));
data3.push_back(edge3d(curr, next, is_valid(curr, next)));
}
// close path if initial and terminal points coincide
bool closed(is_curve &&
map(init) == map(init + P(0,0,count3*st3)));
path_data path3(data3, closed && filling, filling);
path3.draw();
}
}
}; // end of plot_map_dom
// plot over a list of domains
template<class T> void plot_map_domlist(const T& map, const domain_list& R)
{
for (dolci p=R.m_list.begin(); p != R.m_list.end(); ++p)
plot_map_dom(map, *p);
}
// paths
template<class T>void plot_function(const T& f, double t1, double t2,
unsigned int n)
{
plot_map_dom(f, domain(t1, t2, n));
}
// Solutions of ODE systems
template<class VF> void euler_plot(const VF& field, P curr,
double t_min, double t_max,
unsigned int num_pts)
{
std::list<edge3d> data;
const double dt(t_max/(num_pts*EPIX_ITERATIONS));
const double dseek(t_min/(num_pts*EPIX_ITERATIONS));
if (fabs(t_min/num_pts) > EPIX_EPSILON) // seek beginning of path
for (unsigned int i=0; i <= num_pts*EPIX_ITERATIONS; ++i)
curr += dseek*field(curr);
// curr = "start, flowed by t_min"
P next(curr);
for (unsigned int i=0; i <= num_pts*EPIX_ITERATIONS; ++i)
{
next += dt*field(next); // Euler's method
if (i%EPIX_ITERATIONS == 0)
{
data.push_back(edge3d(curr, next, is_valid(curr, next)));
curr = next;
}
}
path_data temp(data, false, false);
temp.draw();
} // end of euler_plot
// flow x0 under field for specified time; pass x0 by value
template<class VF> P euler_flow(const VF& field, P x0,
double t_max, unsigned int num_pts=0)
{
if (num_pts == 0) // use "sensible" default; hardwired constant 4
num_pts = 4*(1 + (unsigned int)ceil(fabs(t_max)));
const double dt(t_max/(num_pts*EPIX_ITERATIONS));
for (unsigned int i=0; i <= num_pts*EPIX_ITERATIONS; ++i)
x0 += dt*field(x0);
return x0;
}
} // end of namespace
#endif /* EPIX_PLOT_ALGO */