m_twotorus.h

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////////
///
/// @file m_twotorus.h
///
/// Two circles at a torus - a simple test map with two invariant circles.
///
/// This map defines a simple test mapping with a formula given by Jacek.
///
/// @author Pawel Pilarczyk
///
/////////////////////////////////////////////////////////////////////////////
 
// Copyright (C) 1997-2014 by Pawel Pilarczyk.
//
// This file is part of my research software package. This 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 3 of the License, or (at your option) any later version.
//
// This software 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 this software; see the file "license.txt". If not,
// please, see <https://www.gnu.org/licenses/>.
 
// Started on January 20, 2012. Last revision: January 20, 2012.
 
 
#ifndef _CMGRAPHS_TWOTORUS_H_
#define _CMGRAPHS_TWOTORUS_H_
 
 
// some standard header files
#include <cmath>
 
// include local header files
#include "maptype.h"
#include "typeintv.h"
 
// for debugging
#include "chomp/system/textfile.h"
 
 
// --------------------------------------------------
// ------------ a two circles model map -------------
// --------------------------------------------------
 
/// This class defines a test mapping by a formula from Jacek Szybowski.
/// There is only one parameter to the map. If the left value is <= 1
/// then the first map is selected, otherwise the second map is selected.
class MapTwoTorus: public MapType
{
public:
        /// The default constructor.
        MapTwoTorus ();
 
        /// Computes the image of a box whose left and right coordinates
        /// are given. Fills in the images with the left and right
        /// coordinates of the image box.
        void compute (const double *xleft, const double *xright,
                double *yleft, double *yright, int dim,
                const spcCoord *coord, int subdiv) const;
 
private:
        /// Computes the core map. Replaces 'x' and 'y' with their images.
        void computeMap0 (IntervalType &x, IntervalType &y) const;
 
}; /* class MapTwoTorus */
 
// --------------------------------------------------
 
inline MapTwoTorus::MapTwoTorus ()
{
        testIntervals ();
        return;
} /* MapTwoTorus::MapTwoTorus */
 
inline void MapTwoTorus::computeMap0 (IntervalType &x, IntervalType &y) const
{
        using chomp::homology::sbug;
 
        y = 16 * power (y - x / 2 - 1.0 / 4, 3) + 1.0 / 4 + x / 2;
 
        return;
} /* MapTwoTorus::computeMap0 */
 
inline void MapTwoTorus::compute (const double *xleft, const double *xright,
        double *yleft, double *yright, int dim, const spcCoord *, int) const
{
        using chomp::homology::sbug;
 
        // this formula is programmed for phase space of dim 2 only
        if (dim != 2)
                throw "Mapping for two circles on a torus: dim = 2 only!";
 
        // debug: show the values
//      sbug << "input = [" << xleft [0] << "," << xright [0] <<
//              "] x [" << xleft [1] << "," << xright [1] << "]\n";
 
        // determine the map choice: 0, 1
        // corresponding to Example 3, f_1 or f_2
        int whichMap = static_cast<int> (std::floor (getLeftParam (0)));
 
        // the first map is very simple, so apply the formula immediately
        if (whichMap == 0)
        {
                IntervalType y (xleft [1], xright [1]);
                y = 4 * power (y - 1.0 / 2, 3) + 1.0 / 2;
                yleft [0] = xleft [0];
                yright [0] = xright [0];
                yleft [1] = y. leftBound ();
                yright [1] = y. rightBound ();
                return;
        }
 
        // get the intervals for the x and y variables
        IntervalType x (xleft [0], xright [0]);
        IntervalType y (xleft [1], xright [1]);
 
        // check if the box touches the inside and the outside the region D
        bool inside (((x / 2). leftBound () <= y. rightBound ()) &&
                (y. leftBound () <= (x / 2 + 1.0 / 2). rightBound ()));
        bool outside (((x / 2). rightBound () > y. leftBound ()) ||
                (y. rightBound () > (x / 2 + 1.0 / 2). leftBound ()));
 
        // compute the image in the region D if necessary
        IntervalType xInD (x);
        IntervalType yInD (y);
        if (inside)
        {
                computeMap0 (xInD, yInD);
        }
 
        // compute the image outside the region D if necessary
        IntervalType xOutD (x);
        IntervalType yOutD (y);
        if (outside)
        {
                // prepare y + 1/2
                IntervalType yHalf (yOutD + 1.0 / 2);
 
                // if it sticks out of the [0,1] interval, shift it backwards
                if (yHalf. leftBound () >= 1.0)
                        yHalf -= 1.0;
 
                // compute f_0 (x, y + 1/2); it may be necessary to split y
                if (yHalf. rightBound () <= 1.0)
                {
                        IntervalType xCopy (xOutD);
                        computeMap0 (xCopy, yHalf);
                }
                else
                {
                        IntervalType xCopy1 (xOutD);
                        IntervalType y1 (yHalf. leftBound (), 1.0);
                        computeMap0 (xCopy1, y1);
                        IntervalType xCopy2 (xOutD);
                        IntervalType y2 (0.0, (yHalf - 1.0). rightBound ());
                        computeMap0 (xCopy2, y2);
                        yHalf = intervalHull (y1, y2);
                }
 
                // add 1/2 to the result
                yOutD = yHalf + IntervalType (1.0 / 2, 1.0 / 2);
 
                // if it sticks out of the [0,1] interval, shift it backwards
                if (yOutD. leftBound () >= 1.0)
                        yOutD -= 1.0;
        }
 
        // compute the resulting interval
        if (!outside)
        {
                x = xInD;
                y = yInD;
        }
        else if (!inside)
        {
                x = xOutD;
                y = yOutD;
        }
        else
        {
                x = intervalHull (xInD, xOutD);
                y = intervalHull (yInD, yOutD);
        }
 
        // copy the result values to the output variables
        yleft [0] = x. leftBound ();
        yright [0] = x. rightBound ();
        yleft [1] = y. leftBound ();
        yright [1] = y. rightBound ();
 
        // swich the rounding direction to the neutral one
        resetRounding ();
        return;
} /* MapTwoTorus::compute */
 
 
#endif // _CMGRAPHS_TWOTORUS_H_