m_twocirc.h

Go to the documentation of this file.

/////////////////////////////////////////////////////////////////////////////
///
/// @file m_twocirc.h
///
/// Two circles - 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 May 21, 2011. Last revision: January 17, 2012.
 
 
#ifndef _CMGRAPHS_TWOCIRC_H_
#define _CMGRAPHS_TWOCIRC_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 MapTwoCirc: public MapType
{
public:
        /// The default constructor.
        MapTwoCirc ();
 
        /// 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, int root) const;
 
}; /* class MapTwoCirc */
 
// --------------------------------------------------
 
inline MapTwoCirc::MapTwoCirc ()
{
        testIntervals ();
        return;
} /* MapTwoCirc::MapTwoCirc */
 
inline void MapTwoCirc::computeMap0 (IntervalType &x, IntervalType &y,
        int root) const
{
        using chomp::homology::sbug;
 
        // apply the formula if at least one variable does not contain 0
        if (((x. leftBound () > 0) || (x. rightBound () < 0)) ||
                ((y. leftBound () > 0) || (y. rightBound () < 0)))
        {
                IntervalType denomTmp = power (x, 2) + power (y, 2);
                IntervalType denominator = (root == 2) ?
                        sqrt (denomTmp) : power (denomTmp,
                        IntervalType (1.0 / root, 1.0 / root));
                if ((denominator. leftBound () <= 0) &&
                        (denominator. rightBound () >= 0))
                {
                        throw "Denominator contains zero.";
                }
                x /= denominator;
                y /= denominator;
                return;
        }
 
        // if x == 0
        if ((x. leftBound () == 0) && (x. rightBound () == 0))
        {
                x = IntervalType (0.0, 0.0);
                if (root == 2)
                {
                        y = IntervalType (-1.0, 1.0);
                        return;
                }
                double yMin = y. leftBound ();
                if (yMin != 0)
                {
                        y = (IntervalType (yMin, yMin) /
                                power (IntervalType (-yMin, -yMin),
                                IntervalType (2.0 / root, 2.0 / root)));
                }
                double yMax = y. rightBound ();
                if (yMax != 0)
                {
                        y = (IntervalType (yMax, yMax) /
                                power (IntervalType (yMax, yMax),
                                IntervalType (2.0 / root, 2.0 / root)));
                }
                y = IntervalType (yMin, yMax);
                return;
        }
 
        // if y == 0
        if ((y. leftBound () == 0) && (y. rightBound () == 0))
        {
                computeMap0 (y, x, root);
                return;
        }
 
        // calculate the endpoints separately if both variables contain 0
        double xMaxAbs ((-(x. leftBound ()) > x. rightBound ()) ?
                -(x. leftBound ()) : x. rightBound ());
        double yMaxAbs ((-(y. leftBound ()) > y. rightBound ()) ?
                -(y. leftBound ()) : y. rightBound ());
        IntervalType xMaxAbs2 (power (IntervalType (xMaxAbs, xMaxAbs), 2));
        IntervalType yMaxAbs2 (power (IntervalType (yMaxAbs, yMaxAbs), 2));
 
        // calculate xLeft
        double xLeft (0);
        if (x. leftBound () == 0)
        {
                xLeft = (root == 2) ? -1 : 0;
        }
        else
        {
                IntervalType numerator (x. leftBound (), x. leftBound ());
                IntervalType fraction2 (2.0 / root, 2.0 / root);
                IntervalType denominator (power (-numerator, fraction2));
                xLeft = (numerator / denominator). leftBound ();
        }
 
        // calculate yLeft
        double yLeft (0);
        if (y. leftBound () == 0)
        {
                yLeft = (root == 2) ? -1 : 0;
        }
        else
        {
                IntervalType numerator (y. leftBound (), y. leftBound ());
                IntervalType fraction2 (2.0 / root, 2.0 / root);
                IntervalType denominator (power (-numerator, fraction2));
                yLeft = (numerator / denominator). leftBound ();
        }
 
        // calculate xRight
        double xRight (0);
        if (x. rightBound () == 0)
        {
                xRight = (root == 2) ? 1 : 0;
        }
        else
        {
                IntervalType numerator (x. rightBound (), x. rightBound ());
                IntervalType fraction2 (2.0 / root, 2.0 / root);
                IntervalType denominator (power (numerator, fraction2));
                xRight = (numerator / denominator). leftBound ();
        }
 
        // calculate yRight
        double yRight (0);
        if (y. rightBound () == 0)
        {
                yRight = (root == 2) ? 1 : 0;
        }
        else
        {
                IntervalType numerator (y. rightBound (), y. rightBound ());
                IntervalType fraction2 (2.0 / root, 2.0 / root);
                IntervalType denominator (power (numerator, fraction2));
                yRight = (numerator / denominator). leftBound ();
        }
 
        x = IntervalType (xLeft, xRight);
        y = IntervalType (yLeft, yRight);
 
        return;
} /* MapTwoCirc::computeMap0 */
 
inline void MapTwoCirc::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 "Test mapping for two circles: 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, 2 or 3,
        // corresponding to Example 1, f_1 or f_2, and Example 2, f_1 or f_2
        int whichMap = static_cast<int> (std::floor (getLeftParam (0)));
 
        // 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 is inside or outside the disk or radius 2
        IntervalType radius2 = power (x, 2) + power (y, 2);
 
        // compute the image of the auxiliary base map
        IntervalType x4 (x);
        IntervalType y4 (y);
        computeMap0 (x4, y4, 4);
        IntervalType x2 (x);
        IntervalType y2 (y);
        if (whichMap > 1)
                computeMap0 (x2, y2, 2);
 
        // determine the coefficients
        double coefInt = 0;
        double coefExt = 0;
        IntervalType xInt (x4);
        IntervalType yInt (y4);
        IntervalType xExt (x4);
        IntervalType yExt (y4);
        switch (whichMap)
        {
        case 0:
                coefInt = 1.0;
                coefExt = std::sqrt (2);
                break;
        case 1:
                coefInt = std::sqrt (2);
                coefExt = std::sqrt (3);
                break;
        case 2:
                coefInt = 1.0;
                coefExt = std::sqrt (2);
                xInt = x2;
                yInt = y2;
                break;
        case 3:
                coefInt = std::sqrt (2);
                coefExt = 3;
                xExt = x2;
                yExt = y2;
                break;
        default:
                throw "Wrong parameter value - unknown map requested.";
        }
 
        // multiply the result by the right coefficients
        if (radius2. rightBound () < 4)
        {
                x = coefInt * xInt;
                y = coefInt * yInt;
        }
        else if (radius2. leftBound () > 4)
        {
                x = coefExt * xExt;
                y = coefExt * yExt;
        }
        else
        {
                x = intervalHull (coefInt * xInt, coefExt * xExt);
                y = intervalHull (coefInt * yInt, coefExt * yExt);
        }
 
        // copy the result values
        yleft [0] = x. leftBound ();
        yright [0] = x. rightBound ();
        yleft [1] = y. leftBound ();
        yright [1] = y. rightBound ();
 
        // debug: show the values
/*      sbug << "output = [" << yleft [0] << "," << yright [0] <<
                "] x [" << yleft [1] << "," << yright [1] << "]\n";
        double exp0 = (yright [0] - yleft [0]) / (xright [0] - xleft [0]);
        double exp1 = (yright [1] - yleft [1]) / (xright [1] - xleft [1]);
        sbug << "expansion: " << exp0 << " x " << exp1 << "\n";
*/
        // swich the rounding direction to the neutral one
        resetRounding ();
        return;
} /* MapTwoCirc::compute */
 
 
#endif // _CMGRAPHS_TWOCIRC_H_