/****************************************************************************
**
** This file is part of the LibreCAD project, a 2D CAD program
**
** Copyright (C) 2011-2015 Dongxu Li (dongxuli2011@gmail.com)
** Copyright (C) 2010 R. van Twisk (librecad@rvt.dds.nl)
** Copyright (C) 2001-2003 RibbonSoft. All rights reserved.
**
**
** This file may be distributed and/or modified under the terms of the
** GNU General Public License version 2 as published by the Free Software
** Foundation and appearing in the file gpl-2.0.txt included in the
** packaging of this file.
**
** This program 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 program; if not, write to the Free Software
** Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
**
** This copyright notice MUST APPEAR in all copies of the script!
**
**********************************************************************/

#include "rs_ellipse.h"

#include "rs_circle.h"
#include "rs_line.h"
#include "rs_graphic.h"
#include "rs_graphicview.h"
#include "rs_painter.h"
#include "rs_information.h"
#include "rs_linetypepattern.h"
#include "rs_math.h"
#include  "lc_quadratic.h"
#include "rs_painterqt.h"
#include "rs_debug.h"

#ifdef EMU_C99
#include "emu_c99.h" /* C99 math */
#endif
// Workaround for Qt bug: https://bugreports.qt-project.org/browse/QTBUG-22829
// TODO: the Q_MOC_RUN detection shouldn't be necessary after this Qt bug is resolved
#ifndef Q_MOC_RUN
#include <boost/version.hpp>
#include <boost/math/tools/roots.hpp>
#include <boost/math/special_functions/ellint_2.hpp>
#if BOOST_VERSION > 104500
#include <boost/math/tools/tuple.hpp>
#endif
#endif

namespace{
//functor to solve for distance, used by snapDistance
class EllipseDistanceFunctor
{
public:
	EllipseDistanceFunctor(RS_Ellipse const* ellipse, double const& target) :
		distance{target}
	  , e{ellipse}
	  , ra{e->getMajorRadius()}
	  , k2{1.- e->getRatio()*e->getRatio()}
	  , k2ra{k2 * ra}
	{
	}
	void setDistance(const double& target){
		distance=target;
	}
#if BOOST_VERSION > 104500
	boost::math::tuple<double, double, double> operator()(double const& z) const {
#else
	boost::fusion::tuple<double, double, double> operator()(double const& z) const {
#endif
	double const cz=cos(z);
	double const sz=sin(z);
		//delta amplitude
	double const d=sqrt(1-k2*sz*sz);
		// return f(x), f'(x) and f''(x)
#if BOOST_VERSION > 104500
	return boost::math::make_tuple(
#else
	return boost::fusion::make_tuple(
#endif
					e->getEllipseLength(z)-distance,
					ra*d,
					k2ra*sz*cz/d
					);
	}

private:
	double distance;
	RS_Ellipse const* const e;
	const double ra;
	const double k2;
	const double k2ra;
};

/**
 * @brief getNearestDistHelper find end point after trimmed by amount
 * @param e ellipse which is not reversed, assume ratio (a/b) >= 1
 * @param trimAmount the length of the trimmed is increased by this amount
 * @param coord current mouse position
 * @param dist if this pointer is not nullptr, save the distance from the new
 * end point to mouse position coord
 * @return the new end point of the trimmed. Only one end of the entity is
 *  trimmed
 */
RS_Vector getNearestDistHelper(RS_Ellipse const& e,
							   double trimAmount,
							   RS_Vector const& coord,
							   double* dist = nullptr)
{
	double const x1 = e.getAngle1();

	double const guess= x1 + M_PI;
	int const digits=std::numeric_limits<double>::digits;

	double const wholeLength = e.getEllipseLength(0, 0); // start/end angle 0 is used for whole ellipses

	double trimmed = e.getLength() + trimAmount;

	// choose the end to trim by the mouse position coord
	bool const trimEnd = coord.squaredTo(e.getStartpoint()) <= coord.squaredTo(e.getEndpoint());

	if (trimEnd)
		trimmed = trimAmount > 0 ? wholeLength - trimAmount : - trimAmount;

	//solve equation of the distance by second order newton_raphson
	EllipseDistanceFunctor X(&e, trimmed);
	using namespace boost::math::tools;
	double const sol =
			halley_iterate<EllipseDistanceFunctor,double>(X,
														  guess,
														  x1,
														  x1 + 2 * M_PI - RS_TOLERANCE_ANGLE,
														  digits);

	RS_Vector const vp = e.getEllipsePoint(sol);
	if (dist) *dist = vp.distanceTo(coord);
	return vp;
}
}

std::ostream& operator << (std::ostream& os, const RS_EllipseData& ed) {
	os << "(" << ed.center <<
		  " " << ed.majorP <<
		  " " << ed.ratio <<
		  " " << ed.angle1 <<
		  "," << ed.angle2 <<
		  ")";
	return os;
}

/**
 * Constructor.
 */
RS_Ellipse::RS_Ellipse(RS_EntityContainer* parent,
                       const RS_EllipseData& d)
	:RS_AtomicEntity(parent)
	,data(d) {
    //calculateEndpoints();
    calculateBorders();
}

RS_Entity* RS_Ellipse::clone() const {
	RS_Ellipse* e = new RS_Ellipse(*this);
	e->initId();
	return e;
}

/**
 * Calculates the boundary box of this ellipse.
  * @author Dongxu Li
 */
void RS_Ellipse::calculateBorders() {

    RS_Vector startpoint = getStartpoint();
    RS_Vector endpoint = getEndpoint();

    double minX = std::min(startpoint.x, endpoint.x);
    double minY = std::min(startpoint.y, endpoint.y);
    double maxX = std::max(startpoint.x, endpoint.x);
    double maxY = std::max(startpoint.y, endpoint.y);

    RS_Vector vp;

    double amin,amax,a;
//      x range
//    vp.set(radius1*cos(angle),radius2*sin(angle));
    vp.set(getMajorP().x,getRatio()*getMajorP().y);
    a=vp.angle();

    amin=RS_Math::correctAngle(getAngle1()+a); // to the range of 0 to 2*M_PI
    amax=RS_Math::correctAngle(getAngle2()+a); // to the range of 0 to 2*M_PI
    if( RS_Math::isAngleBetween(M_PI,amin,amax,isReversed()) ) {
        minX= data.center.x-vp.magnitude();
    }

    if( RS_Math::isAngleBetween(2.*M_PI,amin,amax,isReversed()) ) {
        maxX= data.center.x+vp.magnitude();
    }
    //y range
    vp.set(getMajorP().y, -getRatio()*getMajorP().x);
    a=vp.angle();
    amin=RS_Math::correctAngle(getAngle1()+a); // to the range of 0 to 2*M_PI
    amax=RS_Math::correctAngle(getAngle2()+a); // to the range of 0 to 2*M_PI
    if( RS_Math::isAngleBetween(M_PI,amin,amax,isReversed()) ) {
        minY= data.center.y-vp.magnitude();
    }
    if( RS_Math::isAngleBetween(2.*M_PI,amin,amax,isReversed()) ) {
        maxY= data.center.y+vp.magnitude();
    }

    minV.set(minX, minY);
	maxV.set(maxX, maxY);
}


/**
  * return the foci of ellipse
  *
  * @author Dongxu Li
  */

RS_VectorSolutions RS_Ellipse::getFoci() const {
	RS_Ellipse e=*this;
	if(getRatio()>1.)
		e.switchMajorMinor();
	RS_Vector vp(e.getMajorP()*sqrt(1.-e.getRatio()*e.getRatio()));
	return RS_VectorSolutions({getCenter()+vp, getCenter()-vp});
}

RS_VectorSolutions RS_Ellipse::getRefPoints() const
{
    RS_VectorSolutions ret;
	if(isEllipticArc()){
        //no start/end point for whole ellipse
        ret.push_back(getStartpoint());
        ret.push_back(getEndpoint());
    }
    ret.push_back(data.center);
    ret.push_back(getFoci());
    ret.push_back(getMajorPoint());
    ret.push_back(getMinorPoint());
    return ret;
}



RS_Vector RS_Ellipse::getNearestEndpoint(const RS_Vector& coord, double* dist)const {
    double dist1, dist2;
    RS_Vector startpoint = getStartpoint();
    RS_Vector endpoint = getEndpoint();

    dist1 = (startpoint-coord).squared();
    dist2 = (endpoint-coord).squared();

    if (dist2<dist1) {
		if (dist) {
            *dist = sqrt(dist2);
        }
        return endpoint;
    } else {
		if (dist) {
            *dist = sqrt(dist1);
        }
        return startpoint;
    }
}

/**
  *find the tangential points from a given point, i.e., the tangent lines should pass
  * the given point and tangential points
  *
  * \author Dongxu Li
  */
RS_VectorSolutions RS_Ellipse::getTangentPoint(const RS_Vector& point) const {
    RS_Vector point2(point);
    point2.move(-getCenter());
    RS_Vector aV(-getAngle());
    point2.rotate(aV);
    RS_VectorSolutions sol;
    double a=getMajorRadius();
    if(a<RS_TOLERANCE || getRatio()<RS_TOLERANCE) return sol;
	RS_Circle c(nullptr, RS_CircleData(RS_Vector(0.,0.),a));
    point2.y /=getRatio();
    sol=c.getTangentPoint(point2);
    sol.scale(RS_Vector(1.,getRatio()));
    aV.y *= -1.;
    sol.rotate(aV);
    sol.move(getCenter());
    return sol;
}

RS_Vector RS_Ellipse::getTangentDirection(const RS_Vector& point) const {
    RS_Vector vp(point-getCenter());
    RS_Vector aV(-getAngle());
    vp.rotate(aV);
	vp.y /= getRatio();
    double a=getMajorRadius();
    if(a<RS_TOLERANCE || getRatio()<RS_TOLERANCE) return RS_Vector(false);
	RS_Circle c(nullptr, RS_CircleData(RS_Vector(0.,0.),a));
    RS_Vector direction=c.getTangentDirection(vp);
    direction.y *= getRatio();
    aV.y *= -1.;
    direction.rotate(aV);
    return direction;
}

/**
  * find total length of the ellipse (arc)
  *
  * \author: Dongxu Li
  */
double RS_Ellipse::getLength() const
{
		RS_Ellipse e(nullptr, data);
        //switch major/minor axis, because we need the ratio smaller than one
    if(e.getRatio()>1.)  e.switchMajorMinor();
    if(e.isReversed()) {
        e.setReversed(false);
        std::swap(e.data.angle1,e.data.angle2);
    }
    return e.getEllipseLength(e.data.angle1,e.data.angle2);
}

/**
//Ellipse must have ratio<1, and not reversed
*@ x1, ellipse angle
*@ x2, ellipse angle
//@return the arc length between ellipse angle x1, x2
* \author Dongxu Li
**/
double RS_Ellipse::getEllipseLength(double x1, double x2) const
{
    double a(getMajorRadius()),k(getRatio());
    k= 1-k*k;//elliptic modulus, or eccentricity
//    std::cout<<"1, angle1="<<x1/M_PI<<" angle2="<<x2/M_PI<<std::endl;
//    if(isReversed())  std::swap(x1,x2);
    x1=RS_Math::correctAngle(x1);
    x2=RS_Math::correctAngle(x2);
//    std::cout<<"2, angle1="<<x1/M_PI<<" angle2="<<x2/M_PI<<std::endl;
    if(x2 < x1+RS_TOLERANCE_ANGLE) x2 += 2.*M_PI;
    double ret;
//    std::cout<<"3, angle1="<<x1/M_PI<<" angle2="<<x2/M_PI<<std::endl;
    if( x2 >= M_PI) {
        // the complete elliptic integral
        ret=  (static_cast< int>((x2+RS_TOLERANCE_ANGLE)/M_PI) -
               (static_cast<int>((x1+RS_TOLERANCE_ANGLE)/M_PI)
                ))*2;
//        std::cout<<"Adding "<<ret<<" of E("<<k<<")\n";
        ret*=boost::math::ellint_2<double>(k);
    } else {
        ret=0.;
    }
    x1=fmod(x1,M_PI);
    x2=fmod(x2,M_PI);
    if( fabs(x2-x1)>RS_TOLERANCE_ANGLE)  {
        ret += RS_Math::ellipticIntegral_2(k,x2)-RS_Math::ellipticIntegral_2(k,x1);
    }
    return a*ret;
}

/**
  * arc length from start point (angle1)
  */
double RS_Ellipse::getEllipseLength(double x2) const
{
    return getEllipseLength(getAngle1(),x2);
}


/**
  * get the point on the ellipse arc and with arc distance from the start point
  * the distance is expected to be within 0 and getLength()
  * using Newton-Raphson from boost
  *
  *@author: Dongxu Li
  */

RS_Vector RS_Ellipse::getNearestDist(double distance,
                                     const RS_Vector& coord,
									 double* dist) const{
//    RS_DEBUG->print("RS_Ellipse::getNearestDist() begin\n");
	if( ! isEllipticArc() ) {
        // both angles being 0, whole ellipse
        // no end points for whole ellipse, therefore, no snap by distance from end points.
		return {};
    }
	RS_Ellipse e(nullptr,data);
	if(e.getRatio()>1.) e.switchMajorMinor();
    if(e.isReversed()) {
        std::swap(e.data.angle1,e.data.angle2);
        e.setReversed(false);
    }

	if(e.getMajorRadius() < RS_TOLERANCE)
		return {}; //ellipse too small

    if(getRatio()<RS_TOLERANCE) {
        //treat the ellipse as a line
		RS_Line line{e.minV,e.maxV};
		return line.getNearestDist(distance, coord, dist);
    }
    double x1=e.getAngle1();
    double x2=e.getAngle2();
	if(x2 < x1+RS_TOLERANCE_ANGLE) x2 += 2 * M_PI;
	double const l0=e.getEllipseLength(x1,x2); // the getEllipseLength() function only defined for proper e
//    distance=fabs(distance);
	if(distance > l0+RS_TOLERANCE) return {}; // can not trim more than the current length
	if(distance > l0-RS_TOLERANCE) return(getNearestEndpoint(coord,dist)); // trim to zero length
	return getNearestDistHelper(e, distance, coord, dist);
}


/**
  * switch the major/minor axis naming
  *
  * \author: Dongxu Li
  */
bool RS_Ellipse::switchMajorMinor(void)
//switch naming of major/minor, return true if success
{
    if (fabs(data.ratio) < RS_TOLERANCE) return false;
    RS_Vector vp_start=getStartpoint();
    RS_Vector vp_end=getEndpoint();
    RS_Vector vp=getMajorP();
    setMajorP(RS_Vector(- data.ratio*vp.y, data.ratio*vp.x)); //direction pi/2 relative to old MajorP;
    setRatio(1./data.ratio);
	if( isEllipticArc() )  {
        //only reset start/end points for ellipse arcs, i.e., angle1 angle2 are not both zero
        setAngle1(getEllipseAngle(vp_start));
        setAngle2(getEllipseAngle(vp_end));
    }
    return true;
}

/**
 * @return Start point of the entity.
 */
RS_Vector  RS_Ellipse::getStartpoint() const {
	if(isEllipticArc()) return getEllipsePoint(data.angle1);
	return RS_Vector(false);
}

/**
 * @return End point of the entity.
 */
RS_Vector  RS_Ellipse::getEndpoint() const {
	if(isEllipticArc()) return getEllipsePoint(data.angle2);
	return RS_Vector(false);
}

/**
 * @return Ellipse point by ellipse angle
 */
RS_Vector  RS_Ellipse::getEllipsePoint(const double& a) const {
    RS_Vector p(a);
    double ra=getMajorRadius();
    p.scale(RS_Vector(ra,ra*getRatio()));
    p.rotate(getAngle());
    p.move(getCenter());
    return p;
}

/** \brief implemented using an analytical algorithm
* find nearest point on ellipse to a given point
*
* @author Dongxu Li <dongxuli2011@gmail.com>
*/

RS_Vector RS_Ellipse::getNearestPointOnEntity(const RS_Vector& coord,
        bool onEntity, double* dist, RS_Entity** entity)const
{

    RS_DEBUG->print("RS_Ellipse::getNearestPointOnEntity");
    RS_Vector ret(false);

    if( ! coord.valid ) {
		if ( dist ) *dist=RS_MAXDOUBLE;
        return ret;

    }

	if (entity) {
        *entity = const_cast<RS_Ellipse*>(this);
    }
    ret=coord;
    ret.move(-getCenter());
    ret.rotate(-getAngle());
    double x=ret.x,y=ret.y;
    double a=getMajorRadius();
    double b=a*getRatio();
    //std::cout<<"(a= "<<a<<" b= "<<b<<" x= "<<x<<" y= "<<y<<" )\n";
    //std::cout<<"finding minimum for ("<<x<<"-"<<a<<"*cos(t))^2+("<<y<<"-"<<b<<"*sin(t))^2\n";
    double twoa2b2=2*(a*a-b*b);
    double twoax=2*a*x;
    double twoby=2*b*y;
    double a0=twoa2b2*twoa2b2;
    std::vector<double> ce(4,0.);
    std::vector<double> roots(0,0.);

    //need to handle a=b
    if(a0 > RS_TOLERANCE2 ) { // a != b , ellipse
        ce[0]=-2.*twoax/twoa2b2;
        ce[1]= (twoax*twoax+twoby*twoby)/a0-1.;
        ce[2]= - ce[0];
        ce[3]= -twoax*twoax/a0;
        //std::cout<<"1::find cosine, variable c, solve(c^4 +("<<ce[0]<<")*c^3+("<<ce[1]<<")*c^2+("<<ce[2]<<")*c+("<<ce[3]<<")=0,c)\n";
        roots=RS_Math::quarticSolver(ce);
    } else {//a=b, quadratic equation for circle
        a0=twoby/twoax;
        roots.push_back(sqrt(1./(1.+a0*a0)));
        roots.push_back(-roots[0]);
    }
    if(roots.size()==0) {
        //this should not happen
        std::cout<<"(a= "<<a<<" b= "<<b<<" x= "<<x<<" y= "<<y<<" )\n";
        std::cout<<"finding minimum for ("<<x<<"-"<<a<<"*cos(t))^2+("<<y<<"-"<<b<<"*sin(t))^2\n";
        std::cout<<"2::find cosine, variable c, solve(c^4 +("<<ce[0]<<")*c^3+("<<ce[1]<<")*c^2+("<<ce[2]<<")*c+("<<ce[3]<<")=0,c)\n";
        std::cout<<ce[0]<<' '<<ce[1]<<' '<<ce[2]<<' '<<ce[3]<<std::endl;
        std::cerr<<"RS_Math::RS_Ellipse::getNearestPointOnEntity() finds no root from quartic, this should not happen\n";
        return RS_Vector(coord); // better not to return invalid: return RS_Vector(false);
    }

//    RS_Vector vp2(false);
	double d,dDistance(RS_MAXDOUBLE*RS_MAXDOUBLE);
    //double ea;
    for(size_t i=0; i<roots.size(); i++) {
        //I don't understand the reason yet, but I can do without checking whether sine/cosine are valid
        //if ( fabs(roots[i])>1.) continue;
		double const s=twoby*roots[i]/(twoax-twoa2b2*roots[i]); //sine
        //if (fabs(s) > 1. ) continue;
		double const d2=twoa2b2+(twoax-2.*roots[i]*twoa2b2)*roots[i]+twoby*s;
        if (d2<0) continue; // fartherest
        RS_Vector vp3;
        vp3.set(a*roots[i],b*s);
        d=(vp3-ret).squared();
//        std::cout<<i<<" Checking: cos= "<<roots[i]<<" sin= "<<s<<" angle= "<<atan2(roots[i],s)<<" ds2= "<<d<<" d="<<d2<<std::endl;
        if( ret.valid && d>dDistance) continue;
        ret=vp3;
        dDistance=d;
//			ea=atan2(roots[i],s);
    }
    if( ! ret.valid ) {
        //this should not happen
//        std::cout<<ce[0]<<' '<<ce[1]<<' '<<ce[2]<<' '<<ce[3]<<std::endl;
//        std::cout<<"(x,y)=( "<<x<<" , "<<y<<" ) a= "<<a<<" b= "<<b<<" sine= "<<s<<" d2= "<<d2<<" dist= "<<d<<std::endl;
//        std::cout<<"RS_Ellipse::getNearestPointOnEntity() finds no minimum, this should not happen\n";
        RS_DEBUG->print(RS_Debug::D_ERROR,"RS_Ellipse::getNearestPointOnEntity() finds no minimum, this should not happen\n");
    }
	if (dist) {
        *dist = sqrt(dDistance);
    }
    ret.rotate(getAngle());
    ret.move(getCenter());
//    ret=vp2;
    if (onEntity) {
        if (!RS_Math::isAngleBetween(getEllipseAngle(ret), getAngle1(), getAngle2(), isReversed())) { // not on entity, use the nearest endpoint
               //std::cout<<"not on ellipse, ( "<<getAngle1()<<" "<<getEllipseAngle(ret)<<" "<<getAngle2()<<" ) reversed= "<<isReversed()<<"\n";
            ret=getNearestEndpoint(coord,dist);
        }
    }

//    if(! ret.valid) {
//        std::cout<<"RS_Ellipse::getNearestOnEntity() returns invalid by mistake. This should not happen!"<<std::endl;
//    }
    return ret;
}




/**
 * @param tolerance Tolerance.
 *
 * @retval true if the given point is on this entity.
 * @retval false otherwise
 */
bool RS_Ellipse::isPointOnEntity(const RS_Vector& coord,
                                 double tolerance) const {
    double t=fabs(tolerance);
    double a=getMajorRadius();
    double b=a*getRatio();
    RS_Vector vp((coord - getCenter()).rotate(-getAngle()));
    if ( a<RS_TOLERANCE ) {
        //radius treated as zero
        if(fabs(vp.x)<RS_TOLERANCE && fabs(vp.y) < b) return true;
        return false;
    }
    if ( b<RS_TOLERANCE ) {
        //radius treated as zero
        if (fabs(vp.y)<RS_TOLERANCE && fabs(vp.x) < a) return true;
        return false;
    }
    vp.scale(RS_Vector(1./a,1./b));

    if (fabs(vp.squared()-1.) > t) return false;
	return RS_Math::isAngleBetween(vp.angle(),getAngle1(),getAngle2(),isReversed());
}



RS_Vector RS_Ellipse::getNearestCenter(const RS_Vector& coord,
									   double* dist) const{
    RS_Vector   vCenter = data.center;
    double      distCenter = coord.distanceTo(data.center);

    RS_VectorSolutions  vsFoci = getFoci();
    if( 2 == vsFoci.getNumber()) {
		RS_Vector const& vFocus1 = vsFoci.get(0);
		RS_Vector const& vFocus2 = vsFoci.get(1);

        double distFocus1 = coord.distanceTo(vFocus1);
        double distFocus2 = coord.distanceTo(vFocus2);

        /* if (distFocus1 < distCenter) is true
         * then (distFocus1 < distFocus2) must be true too
         * and vice versa
         * no need to check this */
        if( distFocus1 < distCenter) {
            vCenter = vFocus1;
            distCenter = distFocus1;
        }
        else if( distFocus2 < distCenter) {
            vCenter = vFocus2;
            distCenter = distFocus2;
        }
    }

	if (dist) {
        *dist = distCenter;
    }
    return vCenter;
}

/**
//create Ellipse with axes in x-/y- directions from 4 points
*
*
*@author Dongxu Li
*/
bool	RS_Ellipse::createFrom4P(const RS_VectorSolutions& sol)
{
    if (sol.getNumber() != 4 ) return (false); //only do 4 points
	std::vector<std::vector<double> > mt;
	std::vector<double> dn;
    int mSize(4);
    mt.resize(mSize);
    for(int i=0;i<mSize;i++) {//form the linear equation, c0 x^2 + c1 x + c2 y^2 + c3 y = 1
        mt[i].resize(mSize+1);
        mt[i][0]=sol.get(i).x * sol.get(i).x;
        mt[i][1]=sol.get(i).x ;
        mt[i][2]=sol.get(i).y * sol.get(i).y;
        mt[i][3]=sol.get(i).y ;
        mt[i][4]=1.;
    }
    if ( ! RS_Math::linearSolver(mt,dn)) return false;
    double d(1.+0.25*(dn[1]*dn[1]/dn[0]+dn[3]*dn[3]/dn[2]));
    if(fabs(dn[0])<RS_TOLERANCE15
            ||fabs(dn[2])<RS_TOLERANCE15
            ||d/dn[0]<RS_TOLERANCE15
            ||d/dn[2]<RS_TOLERANCE15
            ) {
        //ellipse not defined
        return false;
    }
    data.center.set(-0.5*dn[1]/dn[0],-0.5*dn[3]/dn[2]); // center
    d=sqrt(d/dn[0]);
    data.majorP.set(d,0.);
    data.ratio=sqrt(dn[0]/dn[2]);
    data.angle1=0.;
    data.angle2=0.;
//    DEBUG_HEADER
//    std::cout<<"center="<<data.center;
//    std::cout<<"majorP="<<data.majorP;
//    std::cout<<"ratio="<<data.ratio;
//    std::cout<<"successful"<<std::endl;
    return true;

}

/**
//create Ellipse with center and 3 points
*
*
*@author Dongxu Li
*/
bool	RS_Ellipse::createFromCenter3Points(const RS_VectorSolutions& sol) {
    if(sol.getNumber()<3) return false; //need one center and 3 points on ellipse
	std::vector<std::vector<double> > mt;
    int mSize(sol.getNumber() -1);
    if( (sol.get(mSize) - sol.get(mSize-1)).squared() < RS_TOLERANCE15 ) {
        //remove the last point
        mSize--;
    }

    mt.resize(mSize);
	std::vector<double> dn(mSize);
    switch(mSize){
    case 2:
        for(int i=0;i<mSize;i++){//form the linear equation
            mt[i].resize(mSize+1);
            RS_Vector vp(sol.get(i+1)-sol.get(0)); //the first vector is center
            mt[i][0]=vp.x*vp.x;
            mt[i][1]=vp.y*vp.y;
            mt[i][2]=1.;
        }
        if ( ! RS_Math::linearSolver(mt,dn) ) return false;
        if( dn[0]<RS_TOLERANCE15 || dn[1]<RS_TOLERANCE15) return false;
        setMajorP(RS_Vector(1./sqrt(dn[0]),0.));
        setRatio(sqrt(dn[0]/dn[1]));
        setAngle1(0.);
        setAngle2(0.);
        setCenter(sol.get(0));
        return true;

    case 3:
        for(int i=0;i<mSize;i++){//form the linear equation
            mt[i].resize(mSize+1);
            RS_Vector vp(sol.get(i+1)-sol.get(0)); //the first vector is center
            mt[i][0]=vp.x*vp.x;
            mt[i][1]=vp.x*vp.y;
            mt[i][2]=vp.y*vp.y;
            mt[i][3]=1.;
        }
        if ( ! RS_Math::linearSolver(mt,dn) ) return false;
        setCenter(sol.get(0));
		return createFromQuadratic(dn);
    default:
        return false;
    }
    return false; // only for compiler warning
}

/** \brief create from quadratic form:
  * dn[0] x^2 + dn[1] xy + dn[2] y^2 =1
  * keep the ellipse center before calling this function
  *
  *@author: Dongxu Li
  */
bool RS_Ellipse::createFromQuadratic(const std::vector<double>& dn){
	RS_DEBUG->print("RS_Ellipse::createFromQuadratic() begin\n");
	if(dn.size()!=3) return false;
//	if(fabs(dn[0]) <RS_TOLERANCE2 || fabs(dn[2])<RS_TOLERANCE2) return false; //invalid quadratic form

	//eigenvalues and eigenvectors of quadratic form
    // (dn[0] 0.5*dn[1])
	// (0.5*dn[1] dn[2])
	double a=dn[0];
	const double c=dn[1];
	double b=dn[2];

	//Eigen system
	const double d = a - b;
	const double s=hypot(d,c);
	// { a>b, d>0
	// eigenvalue: ( a+b - s)/2, eigenvector: ( -c, d + s)
	// eigenvalue: ( a+b + s)/2, eigenvector: ( d + s, c)
	// }
	// { a<b, d<0
	// eigenvalue: ( a+b - s)/2, eigenvector: ( s-d,-c)
	// eigenvalue: ( a+b + s)/2, eigenvector: ( c, s-d)
	// }

	// eigenvalues are required to be positive for ellipses
	if(s >= a+b ) return false;
	if(a>=b) {
		setMajorP(RS_Vector(atan2(d+s, -c))/sqrt(0.5*(a+b-s)));
	}else{
		setMajorP(RS_Vector(atan2(-c, s-d))/sqrt(0.5*(a+b-s)));
	}
	setRatio(sqrt((a+b-s)/(a+b+s)));

	// start/end angle at 0. means a whole ellipse, instead of an elliptic arc
    setAngle1(0.);
	setAngle2(0.);

	RS_DEBUG->print("RS_Ellipse::createFromQuadratic(): successful\n");
	return true;
}

bool RS_Ellipse::createFromQuadratic(const LC_Quadratic& q){
	if (!q.isQuadratic()) return false;
	auto  const& mQ=q.getQuad();
	double const& a=mQ(0,0);
	double const& c=2.*mQ(0,1);
	double const& b=mQ(1,1);
	auto  const& mL=q.getLinear();
	double const& d=mL(0);
	double const& e=mL(1);
	double determinant=c*c-4.*a*b;
	if(determinant>= -DBL_EPSILON) return false;
	// find center of quadratic
	// 2 A x + C y = D
	// C x   + 2 B y = E
	// x = (2BD - EC)/( 4AB - C^2)
	// y = (2AE - DC)/(4AB - C^2)
	const RS_Vector eCenter=RS_Vector(2.*b*d - e*c, 2.*a*e - d*c)/determinant;
	//generate centered quadratic
	LC_Quadratic qCentered=q;
	qCentered.move(-eCenter);
	if(qCentered.constTerm()>= -DBL_EPSILON) return false;
	const auto& mq2=qCentered.getQuad();
	const double factor=-1./qCentered.constTerm();
	//quadratic terms
	if(!createFromQuadratic({mq2(0,0)*factor, 2.*mq2(0,1)*factor, mq2(1,1)*factor})) return false;

	//move back to center
	move(eCenter);
	return true;
}

/**
//create Ellipse inscribed in a quadrilateral
*
*algorithm: http://chrisjones.id.au/Ellipses/ellipse.html
*finding the tangential points and ellipse center
*
*@author Dongxu Li
*/
bool	RS_Ellipse::createInscribeQuadrilateral(const std::vector<RS_Line*>& lines)
{
	if(lines.size() != 4) return false; //only do 4 lines
	std::vector<std::unique_ptr<RS_Line> > quad(4);
	{ //form quadrilateral from intersections
		RS_EntityContainer c0(nullptr, false);
		for(RS_Line*const p: lines){//copy the line pointers
			c0.addEntity(p);
		}
		RS_VectorSolutions const& s0=RS_Information::createQuadrilateral(c0);
		if(s0.size()!=4) return false;
		for(size_t i=0; i<4; ++i){
			quad[i].reset(new RS_Line{s0[i], s0[(i+1)%4]});
		}
	}

	//center of original square projected, intersection of diagonal
	RS_Vector centerProjection;
	{
		std::vector<RS_Line> diagonal;
		diagonal.emplace_back(quad[0]->getStartpoint(), quad[1]->getEndpoint());
		diagonal.emplace_back(quad[1]->getStartpoint(), quad[2]->getEndpoint());
		RS_VectorSolutions const& sol=RS_Information::getIntersectionLineLine( & diagonal[0],& diagonal[1]);
		if(sol.getNumber()==0) {//this should not happen
			//        RS_DEBUG->print(RS_Debug::D_WARNING, "RS_Ellipse::createInscribeQuadrilateral(): can not locate projection Center");
			RS_DEBUG->print("RS_Ellipse::createInscribeQuadrilateral(): can not locate projection Center");
			return false;
		}
		centerProjection=sol.get(0);
	}
	//        std::cout<<"RS_Ellipse::createInscribe(): centerProjection="<<centerProjection<<std::endl;

	std::vector<RS_Vector> tangent;//holds the tangential points on edges, in the order of edges: 1 3 2 0
	int parallel=0;
	int parallel_index=0;
	for(int i=0;i<=1;++i) {
		RS_VectorSolutions const& sol1=RS_Information::getIntersectionLineLine(quad[i].get(), quad[(i+2)%4].get());
		RS_Vector direction;
		if(sol1.getNumber()==0) {
			direction=quad[i]->getEndpoint()-quad[i]->getStartpoint();
			++parallel;
			parallel_index=i;
		}else{
			direction=sol1.get(0)-centerProjection;
		}
		//                std::cout<<"Direction: "<<direction<<std::endl;
		RS_Line l(centerProjection, centerProjection+direction);
		for(int k=1;k<=3;k+=2){
			RS_VectorSolutions sol2=RS_Information::getIntersectionLineLine(&l, quad[(i+k)%4].get());
			if(sol2.size()) tangent.push_back(sol2.get(0));
		}
	}

	if(tangent.size()<3) return false;

	//find ellipse center by projection
	RS_Vector ellipseCenter;
	{
		RS_Line cl0(quad[1]->getEndpoint(),(tangent[0]+tangent[2])*0.5);
		RS_Line cl1(quad[2]->getEndpoint(),(tangent[1]+tangent[2])*0.5);
		RS_VectorSolutions const& sol=RS_Information::getIntersection(&cl0, &cl1,false);
		if(sol.getNumber()==0){
			//this should not happen
			//        RS_DEBUG->print(RS_Debug::D_WARNING, "RS_Ellipse::createInscribeQuadrilateral(): can not locate Ellipse Center");
			RS_DEBUG->print("RS_Ellipse::createInscribeQuadrilateral(): can not locate Ellipse Center");
			return false;
		}
		ellipseCenter=sol.get(0);
	}
	//	qDebug()<<"parallel="<<parallel;
	if(parallel==1){
		RS_DEBUG->print("RS_Ellipse::createInscribeQuadrilateral(): trapezoid detected\n");
		//trapezoid
		RS_Line* l0=quad[parallel_index].get();
		RS_Line* l1=quad[(parallel_index+2)%4].get();
		RS_Vector centerPoint=(l0->getMiddlePoint()+l1->getMiddlePoint())*0.5;
		//not symmetric, no inscribed ellipse
		if( fabs(centerPoint.distanceTo(l0->getStartpoint()) - centerPoint.distanceTo(l0->getEndpoint()))>RS_TOLERANCE)
			return false;
		//symmetric
		RS_DEBUG->print("RS_Ellipse::createInscribeQuadrilateral(): symmetric trapezoid detected\n");
		double d=l0->getDistanceToPoint(centerPoint);
		double l=((l0->getLength()+l1->getLength()))*0.25;
		double k= 4.*d/fabs(l0->getLength()-l1->getLength());
		double theta=d/(l*k);
		if(theta>=1. || d<RS_TOLERANCE) {
			RS_DEBUG->print("RS_Ellipse::createInscribeQuadrilateral(): this should not happen\n");
			return false;
		}
		theta=asin(theta);

		//major axis
		double a=d/(k*tan(theta));
		setCenter(RS_Vector(0., 0.));
		setMajorP(RS_Vector(a, 0.));
		setRatio(d/a);
		rotate(l0->getAngle1());
		setCenter(centerPoint);
		return true;

	}
	//    double ratio;
	//        std::cout<<"dn="<<dn[0]<<' '<<dn[1]<<' '<<dn[2]<<std::endl;
	std::vector<double> dn(3);
	RS_Vector angleVector(false);

	for(size_t i=0;i<tangent.size();i++) {
		tangent[i] -= ellipseCenter;//relative to ellipse center
	}
	std::vector<std::vector<double> > mt;
	mt.clear();
	const double symTolerance=20.*RS_TOLERANCE;
	for(const RS_Vector& vp: tangent){
		//form the linear equation
		// need to remove duplicated {x^2, xy, y^2} terms due to symmetry (x => -x, y=> -y)
		// i.e. rotation of 180 degrees around ellipse center
		//		std::cout<<"point  : "<<vp<<std::endl;
		std::vector<double> mtRow;
		mtRow.push_back(vp.x*vp.x);
		mtRow.push_back(vp.x*vp.y);
		mtRow.push_back(vp.y*vp.y);
		const double l= hypot(hypot(mtRow[0], mtRow[1]), mtRow[2]);
		bool addRow(true);
		for(const auto& v: mt){
			RS_Vector const dv{v[0] - mtRow[0], v[1] - mtRow[1], v[2] - mtRow[2]};
			if( dv.magnitude() < symTolerance*l){
				//symmetric
				addRow=false;
				break;
			}
		}
		if(addRow) {
			mtRow.push_back(1.);
			mt.push_back(mtRow);
		}
	}
	//    std::cout<<"mt.size()="<<mt.size()<<std::endl;
	switch(mt.size()){
	case 2:{// the quadrilateral is a parallelogram
		RS_DEBUG->print("RS_Ellipse::createInscribeQuadrilateral(): parallelogram detected\n");

		//fixme, need to handle degenerate case better
		//        double angle(center.angleTo(tangent[0]));
		RS_Vector majorP(tangent[0]);
		double dx(majorP.magnitude());
		if(dx<RS_TOLERANCE2) return false; //refuse to return zero size ellipse
		angleVector.set(majorP.x/dx,-majorP.y/dx);
		for(size_t i=0;i<tangent.size();i++)tangent[i].rotate(angleVector);

		RS_Vector minorP(tangent[2]);
		double dy2(minorP.squared());
		if(fabs(minorP.y)<RS_TOLERANCE || dy2<RS_TOLERANCE2) return false; //refuse to return zero size ellipse
		// y'= y
		// x'= x-y/tan
		// reverse scale
		// y=y'
		// x=x'+y' tan
		//
		double ia2=1./(dx*dx);
		double ib2=1./(minorP.y*minorP.y);
		//ellipse scaled:drawi
		// ia2*x'^2+ib2*y'^2=1
		// ia2*(x-y*minor.x/minor.y)^2+ib2*y^2=1
		// ia2*x^2 -2*ia2*minor.x/minor.y xy + ia2*minor.x^2*ib2 y^2 + ib2*y^2 =1
		dn[0]=ia2;
		dn[1]=-2.*ia2*minorP.x/minorP.y;
		dn[2]=ib2*ia2*minorP.x*minorP.x+ib2;
	}
		break;
	case 4:
		mt.pop_back(); //only 3 points needed to form the qudratic form
		if ( ! RS_Math::linearSolver(mt,dn) ) return false;
		break;
	default:
		RS_DEBUG->print(RS_Debug::D_WARNING,"No inscribed ellipse for non isosceles trapezoid");
		return false; //invalid quadrilateral
	}

	if(! createFromQuadratic(dn)) return false;
	setCenter(ellipseCenter);

	if(angleVector.valid) {//need to rotate back, for the parallelogram case
		angleVector.y *= -1.;
		rotate(ellipseCenter,angleVector);
	}
	return true;

}

/**
 * a naive implementation of middle point
 * to accurately locate the middle point from arc length is possible by using elliptic integral to find the total arc length, then, using elliptic function to find the half length point
 */
RS_Vector RS_Ellipse::getMiddlePoint()const{
    return getNearestMiddle(getCenter());
}
/**
  * get Nearest equidistant point
  *
  *@author: Dongxu Li
  */
RS_Vector RS_Ellipse::getNearestMiddle(const RS_Vector& coord,
                                       double* dist,
                                       int middlePoints
                                       ) const{
    RS_DEBUG->print("RS_Ellpse::getNearestMiddle(): begin\n");
	if ( ! isEllipticArc() ) {
        //no middle point for whole ellipse, angle1=angle2=0
		if (dist) {
            *dist = RS_MAXDOUBLE;
        }
        return RS_Vector(false);
    }
    double ra(getMajorRadius());
    double rb(getRatio()*ra);
    if ( ra < RS_TOLERANCE || rb < RS_TOLERANCE ) {
        //zero radius, return the center
        RS_Vector vp(getCenter());
		if (dist) {
            *dist = vp.distanceTo(coord);
        }
        return vp;
    }
    double amin=getCenter().angleTo(getStartpoint());
    double amax=getCenter().angleTo(getEndpoint());
    if(isReversed()) {
        std::swap(amin,amax);
    }
    double da=fmod(amax-amin+2.*M_PI, 2.*M_PI);
    if ( da < RS_TOLERANCE ) {
        da = 2.*M_PI; //whole ellipse
    }
    RS_Vector vp(getNearestPointOnEntity(coord,true,dist));
    double a=getCenter().angleTo(vp);
    int counts(middlePoints + 1);
    int i( static_cast<int>(fmod(a-amin+2.*M_PI,2.*M_PI)/da*counts+0.5));
    if(!i) i++; // remove end points
    if(i==counts) i--;
    a=amin + da*(double(i)/double(counts))-getAngle();
    vp.set(a);
    RS_Vector vp2(vp);
    vp2.scale( RS_Vector(1./ra,1./rb));
    vp.scale(1./vp2.magnitude());
    vp.rotate(getAngle());
    vp.move(getCenter());

	if (dist) {
        *dist = vp.distanceTo(coord);
    }
    //RS_DEBUG->print("RS_Ellipse::getNearestMiddle: angle1=%g, angle2=%g, middle=%g\n",amin,amax,a);
    RS_DEBUG->print("RS_Ellpse::getNearestMiddle(): end\n");
    return vp;
}

/**
  * get the tangential point of a tangential line orthogonal to a given line
  *@ normal, the given line
  *@ onEntity, should the tangential be required to on entity of the elliptic arc
  *@ coord, current cursor position
  *
  *@author: Dongxu Li
  */

RS_Vector RS_Ellipse::getNearestOrthTan(const RS_Vector& coord,
                                        const RS_Line& normal,
										bool onEntity ) const
{
    if ( !coord.valid ) {
        return RS_Vector(false);
    }
    RS_Vector direction=normal.getEndpoint() - normal.getStartpoint();
	if (direction.squared()< RS_TOLERANCE15) {
        //undefined direction
        return RS_Vector(false);
    }
    //scale to ellipse angle
    RS_Vector aV(-getAngle());
    direction.rotate(aV);
    double angle=direction.scale(RS_Vector(1.,getRatio())).angle();
    double ra(getMajorRadius());
    direction.set(ra*cos(angle),getRatio()*ra*sin(angle));//relative to center
	std::vector<RS_Vector> sol;
    for(int i=0;i<2;i++){
        if(!onEntity ||
                RS_Math::isAngleBetween(angle,getAngle1(),getAngle2(),isReversed())) {
            if(i){
				sol.push_back(- direction);
            }else{
				sol.push_back(direction);
            }
        }
        angle=RS_Math::correctAngle(angle+M_PI);
    }
    if(sol.size()<1) return RS_Vector(false);
    aV.y*=-1.;
	for(auto& v: sol){
		v.rotate(aV);
	}
    RS_Vector vp;
	switch(sol.size()) {
    case 0:
        return RS_Vector(false);
    case 2:
        if( RS_Vector::dotP(sol[1],coord-getCenter())>0.) {
            vp=sol[1];
            break;
        }
        // fall-through
    default:
        vp=sol[0];
        break;
    }
    return getCenter() + vp;
}


void RS_Ellipse::move(const RS_Vector& offset) {
    data.center.move(offset);
    //calculateEndpoints();
    //    minV.move(offset);
    //    maxV.move(offset);
    moveBorders(offset);
}



void RS_Ellipse::rotate(const RS_Vector& center, const double& angle) {
    RS_Vector angleVector(angle);
    data.center.rotate(center, angleVector);
    data.majorP.rotate(angleVector);
    //calculateEndpoints();
    calculateBorders();
}
void RS_Ellipse::rotate(const RS_Vector& center, const RS_Vector& angleVector) {
    data.center.rotate(center, angleVector);
    data.majorP.rotate(angleVector);
    //calculateEndpoints();
    calculateBorders();
}

void RS_Ellipse::rotate( const double& angle) {//rotate around origin
    RS_Vector aV(angle);
    data.center.rotate(aV);
    data.majorP.rotate(aV);
    calculateBorders();
}
void RS_Ellipse::rotate( const RS_Vector& angleVector) {//rotate around origin
    data.center.rotate(angleVector);
    data.majorP.rotate(angleVector);
    //calculateEndpoints();
    calculateBorders();
}

/**
 * make sure angleLength() is not more than 2*M_PI
 */
void RS_Ellipse::correctAngles() {
        double *pa1= & data.angle1;
        double *pa2= & data.angle2;
        if (isReversed()) std::swap(pa1,pa2);
        *pa2 = *pa1 + fmod(*pa2 - *pa1, 2.*M_PI);
        if ( fabs(data.angle1 - data.angle2) < RS_TOLERANCE_ANGLE ) *pa2 += 2.*M_PI;
}

void RS_Ellipse::moveStartpoint(const RS_Vector& pos) {
    data.angle1 = getEllipseAngle(pos);
    //data.angle1 = data.center.angleTo(pos);
    //calculateEndpoints();
    correctAngles(); // make sure angleLength is no more than 2*M_PI
    calculateBorders();
}



void RS_Ellipse::moveEndpoint(const RS_Vector& pos) {
    data.angle2 = getEllipseAngle(pos);
    //data.angle2 = data.center.angleTo(pos);
    //calculateEndpoints();
    correctAngles(); // make sure angleLength is no more than 2*M_PI
    calculateBorders();
}


RS2::Ending RS_Ellipse::getTrimPoint(const RS_Vector& trimCoord,
                                     const RS_Vector& /*trimPoint*/) {

    //double angEl = getEllipseAngle(trimPoint);
    double angM = getEllipseAngle(trimCoord);
    if (RS_Math::getAngleDifference(angM, data.angle1,isReversed()) > RS_Math::getAngleDifference(data.angle2,angM,isReversed())) {
        return RS2::EndingStart;
    } else {
        return RS2::EndingEnd;
    }
}

RS_Vector RS_Ellipse::prepareTrim(const RS_Vector& trimCoord,
                                  const RS_VectorSolutions& trimSol) {
//special trimming for ellipse arc
        RS_DEBUG->print("RS_Ellipse::prepareTrim()");
    if( ! trimSol.hasValid() ) return (RS_Vector(false));
    if( trimSol.getNumber() == 1 ) return (trimSol.get(0));
    double am=getEllipseAngle(trimCoord);
	std::vector<double> ias;
    double ia(0.),ia2(0.);
    RS_Vector is,is2;
	for(size_t ii=0; ii<trimSol.getNumber(); ++ii) { //find closest according ellipse angle
		ias.push_back(getEllipseAngle(trimSol.get(ii)));
        if( !ii ||  fabs( remainder( ias[ii] - am, 2*M_PI)) < fabs( remainder( ia -am, 2*M_PI)) ) {
            ia = ias[ii];
            is = trimSol.get(ii);
        }
    }
    std::sort(ias.begin(),ias.end());
	for(size_t ii=0; ii<trimSol.getNumber(); ++ii) { //find segment to include trimCoord
        if ( ! RS_Math::isSameDirection(ia,ias[ii],RS_TOLERANCE)) continue;
        if( RS_Math::isAngleBetween(am,ias[(ii+trimSol.getNumber()-1)% trimSol.getNumber()],ia,false))  {
            ia2=ias[(ii+trimSol.getNumber()-1)% trimSol.getNumber()];
        } else {
            ia2=ias[(ii+1)% trimSol.getNumber()];
        }
        break;
    }
	for(const RS_Vector& vp: trimSol) { //find segment to include trimCoord
		if ( ! RS_Math::isSameDirection(ia2,getEllipseAngle(vp),RS_TOLERANCE)) continue;
		is2=vp;
        break;
    }
    if(RS_Math::isSameDirection(getAngle1(),getAngle2(),RS_TOLERANCE_ANGLE)
            ||  RS_Math::isSameDirection(ia2,ia,RS_TOLERANCE) ) {
        //whole ellipse
        if( !RS_Math::isAngleBetween(am,ia,ia2,isReversed())) {
            std::swap(ia,ia2);
            std::swap(is,is2);
        }
        setAngle1(ia);
        setAngle2(ia2);
        double da1=fabs(remainder(getAngle1()-am,2*M_PI));
        double da2=fabs(remainder(getAngle2()-am,2*M_PI));
        if(da2<da1) {
            std::swap(is,is2);
        }

    } else {
        double dia=fabs(remainder(ia-am,2*M_PI));
        double dia2=fabs(remainder(ia2-am,2*M_PI));
        double ai_min=std::min(dia,dia2);
        double da1=fabs(remainder(getAngle1()-am,2*M_PI));
        double da2=fabs(remainder(getAngle2()-am,2*M_PI));
        double da_min=std::min(da1,da2);
        if( da_min < ai_min ) {
            //trimming one end of arc
            bool irev= RS_Math::isAngleBetween(am,ia2,ia, isReversed()) ;
            if ( RS_Math::isAngleBetween(ia,getAngle1(),getAngle2(), isReversed()) &&
                    RS_Math::isAngleBetween(ia2,getAngle1(),getAngle2(), isReversed()) ) { //
                if(irev) {
                    setAngle2(ia);
                    setAngle1(ia2);
                } else {
                    setAngle1(ia);
                    setAngle2(ia2);
                }
                da1=fabs(remainder(getAngle1()-am,2*M_PI));
                da2=fabs(remainder(getAngle2()-am,2*M_PI));
            }
            if( ((da1 < da2) && (RS_Math::isAngleBetween(ia2,ia,getAngle1(),isReversed()))) ||
                    ((da1 > da2) && (RS_Math::isAngleBetween(ia2,getAngle2(),ia,isReversed())))
              ) {
                std::swap(is,is2);
                //std::cout<<"reset: angle1="<<getAngle1()<<" angle2="<<getAngle2()<<" am="<< am<<" is="<<getEllipseAngle(is)<<" ia2="<<ia2<<std::endl;
            }
        } else {
            //choose intersection as new end
            if( dia > dia2) {
                std::swap(is,is2);
                std::swap(ia,ia2);
            }
            if(RS_Math::isAngleBetween(ia,getAngle1(),getAngle2(),isReversed())) {
                if(RS_Math::isAngleBetween(am,getAngle1(),ia,isReversed())) {
                    setAngle2(ia);
                } else {
                    setAngle1(ia);
                }
            }
        }
    }
    return is;
}

double RS_Ellipse::getEllipseAngle(const RS_Vector& pos) const {
    RS_Vector m = pos-data.center;
    m.rotate(-data.majorP.angle());
    m.x *= data.ratio;
    return m.angle();
}

const RS_EllipseData& RS_Ellipse::getData() const
{
	return data;
}



/* Dongxu Li's Version, 19 Aug 2011
 * scale an ellipse
 * Find the eigen vectors and eigen values by optimization
 * original ellipse equation,
 * x= a cos t
 * y= b sin t
 * rotated by angle,
 *
 * x = a cos t cos (angle) - b sin t sin(angle)
 * y = a cos t sin (angle) + b sin t cos(angle)
 * scaled by ( kx, ky),
 * x *= kx
 * y *= ky
 * find the maximum and minimum of x^2 + y^2,
 */
void RS_Ellipse::scale(const RS_Vector& center, const RS_Vector& factor) {
    RS_Vector vpStart;
    RS_Vector vpEnd;
	if(isEllipticArc()){
        //only handle start/end points for ellipse arc
        vpStart=getStartpoint().scale(center,factor);
        vpEnd=getEndpoint().scale(center,factor);
    }
    data.center.scale(center, factor);
    RS_Vector vp1(getMajorP());
    double a(vp1.magnitude());
    if(a<RS_TOLERANCE) return; //ellipse too small
    vp1 *= 1./a;
    double ct=vp1.x;
    double ct2 = ct*ct; // cos^2 angle
    double st=vp1.y;
    double st2=1.0 - ct2; // sin^2 angle
    double kx2= factor.x * factor.x;
    double ky2= factor.y * factor.y;
//    double a=getMajorRadius();
    double b=getRatio()*a;
    double cA=0.5*a*a*(kx2*ct2+ky2*st2);
    double cB=0.5*b*b*(kx2*st2+ky2*ct2);
    double cC=a*b*ct*st*(ky2-kx2);
    if (factor.x < 0)
        setReversed(!isReversed());
    if (factor.y < 0)
        setReversed(!isReversed());
    RS_Vector vp(cA-cB,cC);
    vp1.set(a,b);
    vp1.scale(RS_Vector(0.5*vp.angle()));
    vp1.rotate(RS_Vector(ct,st));
    vp1.scale(factor);
    setMajorP(vp1);
    a=cA+cB;
    b=vp.magnitude();
    setRatio( sqrt((a - b)/(a + b) ));
	if( isEllipticArc() ) {
        //only reset start/end points for ellipse arcs, i.e., angle1 angle2 are not both zero
        setAngle1(getEllipseAngle(vpStart));
        setAngle2(getEllipseAngle(vpEnd));
    }
    correctAngles();//avoid extra 2.*M_PI in angles
    //calculateEndpoints();
    scaleBorders(center,factor);
// calculateBorders();

}


/**
 * is the Ellipse an Arc
 * @return false, if both angle1/angle2 are zero
 *
 *@author: Dongxu Li
 */
bool RS_Ellipse::isEllipticArc() const{
#ifndef EMU_C99
    using std::isnormal;
#endif
    return /*std::*/isnormal(getAngle1()) || /*std::*/isnormal(getAngle2());
}
/**
 * mirror by the axis of the line axisPoint1 and axisPoint2
 *
 *@author: Dongxu Li
 */
void RS_Ellipse::mirror(const RS_Vector& axisPoint1, const RS_Vector& axisPoint2) {
    RS_Vector center=getCenter();
    RS_Vector majorp = center + getMajorP();
    RS_Vector startpoint,endpoint;
	if( isEllipticArc() )  {
        startpoint = getStartpoint();
        endpoint = getEndpoint();
    }

    center.mirror(axisPoint1, axisPoint2);
    majorp.mirror(axisPoint1, axisPoint2);

    setCenter(center);
    setReversed(!isReversed());
    setMajorP(majorp - center);
	if( isEllipticArc() )  {
        //only reset start/end points for ellipse arcs, i.e., angle1 angle2 are not both zero
        startpoint.mirror(axisPoint1, axisPoint2);
        endpoint.mirror(axisPoint1, axisPoint2);
        setAngle1( getEllipseAngle(startpoint));
        setAngle2( getEllipseAngle(endpoint));
    }
    correctAngles();//avoid extra 2.*M_PI in angles
    calculateBorders();
}

/**
  * get direction1 and direction2
  * get the tangent pointing outside at end points
  *
  *@author: Dongxu Li
  */
//getDirection1 for start point
double RS_Ellipse::getDirection1() const {
    RS_Vector vp;
    if (isReversed()){
        vp.set(sin(getAngle1()), -getRatio()*cos(getAngle1()));
    } else {
        vp.set(-sin(getAngle1()), getRatio()*cos(getAngle1()));
    }
    return vp.angle()+getAngle();
}

//getDirection2 for end point
double RS_Ellipse::getDirection2() const {
    RS_Vector vp;
    if (isReversed()){
        vp.set(-sin(getAngle2()), getRatio()*cos(getAngle2()));
    } else {
        vp.set(sin(getAngle2()), -getRatio()*cos(getAngle2()));
    }
    return vp.angle()+getAngle();
}

void RS_Ellipse::moveRef(const RS_Vector& ref, const RS_Vector& offset) {
	if(isEllipticArc()){
        RS_Vector startpoint = getStartpoint();
        RS_Vector endpoint = getEndpoint();

        //    if (ref.distanceTo(startpoint)<1.0e-4) {
		//instead of
		if ((ref-startpoint).squared()<RS_TOLERANCE_ANGLE) {
            moveStartpoint(startpoint+offset);
            correctAngles();//avoid extra 2.*M_PI in angles
            return;
        }
		if ((ref-endpoint).squared()<RS_TOLERANCE_ANGLE) {
            moveEndpoint(endpoint+offset);
            correctAngles();//avoid extra 2.*M_PI in angles
            return;
        }
    }
	if ((ref-getCenter()).squared()<RS_TOLERANCE_ANGLE) {
        //move center
        setCenter(getCenter()+offset);
        return;
    }

    if(data.ratio>1.) switchMajorMinor();
	auto foci=getFoci();
	for(size_t i=0; i< 2 ; i++){
		if ((ref-foci.at(i)).squared()<RS_TOLERANCE_ANGLE) {
			auto focusNew=foci.at(i) + offset;
            //move focus
			auto center = getCenter() + offset*0.5;
            RS_Vector majorP;
            if(getMajorP().dotP( foci.at(i) - getCenter()) >= 0.){
                majorP = focusNew - center;
            }else{
                majorP = center - focusNew;
            }
            double d=getMajorP().magnitude();
            double c=0.5*focusNew.distanceTo(foci.at(1-i));
            double k=majorP.magnitude();
            if(k<RS_TOLERANCE2 || d < RS_TOLERANCE ||
                    c >= d - RS_TOLERANCE) return;
            //            DEBUG_HEADER
            //            std::cout<<__func__<<" : moving focus";
            majorP *= d/k;
            setCenter(center);
            setMajorP(majorP);
            setRatio(sqrt(d*d-c*c)/d);
            correctAngles();//avoid extra 2.*M_PI in angles
            if(data.ratio>1.) switchMajorMinor();
            return;
        }
    }

    //move major/minor points
	if ((ref-getMajorPoint()).squared()<RS_TOLERANCE_ANGLE) {
        RS_Vector majorP=getMajorP()+offset;
        double r=majorP.magnitude();
        if(r<RS_TOLERANCE) return;
        double ratio = getRatio()*getMajorRadius()/r;
        setMajorP(majorP);
        setRatio(ratio);
        if(data.ratio>1.) switchMajorMinor();
        return;
    }
	if ((ref-getMinorPoint()).squared()<RS_TOLERANCE_ANGLE) {
        RS_Vector minorP=getMinorPoint() + offset;
        double r2=getMajorP().squared();
        if(r2<RS_TOLERANCE2) return;
         RS_Vector projected= getCenter() +
                getMajorP()*getMajorP().dotP(minorP-getCenter())/r2;
        double r=(minorP - projected).magnitude();
        if(r<RS_TOLERANCE) return;
        double ratio = getRatio()*r/getMinorRadius();
        setRatio(ratio);
        if(data.ratio>1.) switchMajorMinor();
        return;
    }
}

/** whether the entity's bounding box intersects with visible portion of graphic view
//fix me, need to handle overlay container separately
*/
bool RS_Ellipse::isVisibleInWindow(RS_GraphicView* view) const
{
    RS_Vector vpMin(view->toGraph(0,view->getHeight()));
    RS_Vector vpMax(view->toGraph(view->getWidth(),0));
    //viewport
    QRectF visualRect(vpMin.x,vpMin.y,vpMax.x-vpMin.x, vpMax.y-vpMin.y);
    QPolygonF visualBox(visualRect);
	std::vector<RS_Vector> vps;
    for(unsigned short i=0;i<4;i++){
        const QPointF& vp(visualBox.at(i));
		vps.push_back(RS_Vector(vp.x(),vp.y()));
    }
    //check for intersection points with viewport
    for(unsigned short i=0;i<4;i++){
		RS_Line line{vps.at(i),vps.at((i+1)%4)};
		RS_Ellipse e0(nullptr, getData());
        if( RS_Information::getIntersection(&e0, &line, true).size()>0) return true;
    }
    //is startpoint within viewport
    QRectF ellipseRect(minV.x, minV.y, maxV.x - minV.x, maxV.y - minV.y);
    return ellipseRect.intersects(visualRect);
}

/** return the equation of the entity
for quadratic,

return a vector contains:
m0 x^2 + m1 xy + m2 y^2 + m3 x + m4 y + m5 =0

for linear:
m0 x + m1 y + m2 =0
**/
LC_Quadratic RS_Ellipse::getQuadratic() const
{
    std::vector<double> ce(6,0.);
    ce[0]=data.majorP.squared();
    ce[2]= data.ratio*data.ratio*ce[0];
    if(ce[0]<RS_TOLERANCE2 || ce[2]<RS_TOLERANCE2){
        return LC_Quadratic();
    }
    ce[0]=1./ce[0];
    ce[2]=1./ce[2];
    ce[5]=-1.;
    LC_Quadratic ret(ce);
    ret.rotate(getAngle());
    ret.move(data.center);
    return ret;
}

/**
 * @brief areaLineIntegral, line integral for contour area calculation by Green's Theorem
 * Contour Area =\oint x dy
 * @return line integral \oint x dy along the entity
 * \oint x dy = Cx y + \frac{1}{4}((a^{2}+b^{2})sin(2a)cos^{2}(t)-ab(2sin^{2}(a)sin(2t)-2t-sin(2t)))
 *@author Dongxu Li
 */
double RS_Ellipse::areaLineIntegral() const
{
    const double a=getMajorRadius();
    const double b=getMinorRadius();
	if(!isEllipticArc())
        return M_PI*a*b;
    const double ab=a*b;
    const double r2=a*a+b*b;
    const double& cx=data.center.x;
    const double aE=getAngle();
    const double& a0=data.angle1;
    const double& a1=data.angle2;
    const double fStart=cx*getStartpoint().y+0.25*r2*sin(2.*aE)*cos(a0)*cos(a0)-0.25*ab*(2.*sin(aE)*sin(aE)*sin(2.*a0)-sin(2.*a0));
    const double fEnd=cx*getEndpoint().y+0.25*r2*sin(2.*aE)*cos(a1)*cos(a1)-0.25*ab*(2.*sin(aE)*sin(aE)*sin(2.*a1)-sin(2.*a1));
    return (isReversed()?fStart-fEnd:fEnd-fStart) + 0.5*ab*getAngleLength();
}

bool RS_Ellipse::isReversed() const {
	return data.reversed;
}

void RS_Ellipse::setReversed(bool r) {
	data.reversed = r;
}

double RS_Ellipse::getAngle() const {
	return data.majorP.angle();
}

double RS_Ellipse::getAngle1() const {
	return data.angle1;
}

void RS_Ellipse::setAngle1(double a1) {
	data.angle1 = a1;
}

double RS_Ellipse::getAngle2() const {
	return data.angle2;
}

void RS_Ellipse::setAngle2(double a2) {
	data.angle2 = a2;
}

RS_Vector RS_Ellipse::getCenter() const {
	return data.center;
}

void RS_Ellipse::setCenter(const RS_Vector& c) {
	data.center = c;
}


const RS_Vector& RS_Ellipse::getMajorP() const {
	return data.majorP;
}

void RS_Ellipse::setMajorP(const RS_Vector& p) {
	data.majorP = p;
}

double RS_Ellipse::getRatio() const {
	return data.ratio;
}

void RS_Ellipse::setRatio(double r) {
	data.ratio = r;
}

double RS_Ellipse::getAngleLength() const {
    double ret;
    if (isReversed()) {
        ret= RS_Math::correctAngle(data.angle1-data.angle2);
    } else {
        ret= RS_Math::correctAngle(data.angle2-data.angle1);
    }
    if(ret<RS_TOLERANCE_ANGLE) ret=2.*M_PI;
    return ret;
}


double RS_Ellipse::getMajorRadius() const {
	return data.majorP.magnitude();
}

RS_Vector RS_Ellipse::getMajorPoint() const{
	return data.center + data.majorP;
}

RS_Vector RS_Ellipse::getMinorPoint() const{
	return data.center +
			RS_Vector(-data.majorP.y, data.majorP.x)*data.ratio;
}

double RS_Ellipse::getMinorRadius() const {
	return data.majorP.magnitude()*data.ratio;
}

void RS_Ellipse::draw(RS_Painter* painter, RS_GraphicView* view, double& patternOffset) {
	if(isEllipticArc()==false){
        RS_Ellipse arc(*this);
        arc.setAngle2(2.*M_PI);
        arc.setReversed(false);
        arc.draw(painter,view,patternOffset);
        return;
    }
    //only draw the visible portion of line
    RS_Vector vpMin(view->toGraph(0,view->getHeight()));
    RS_Vector vpMax(view->toGraph(view->getWidth(),0));
    QPolygonF visualBox(QRectF(vpMin.x,vpMin.y,vpMax.x-vpMin.x, vpMax.y-vpMin.y));

    RS_Vector vpStart(isReversed()?getEndpoint():getStartpoint());
    RS_Vector vpEnd(isReversed()?getStartpoint():getEndpoint());

	std::vector<RS_Vector> vertex(0);
    for(unsigned short i=0;i<4;i++){
        const QPointF& vp(visualBox.at(i));
		vertex.emplace_back(vp.x(),vp.y());
    }
    /** angles at cross points */
	std::vector<double> crossPoints(0);

    double baseAngle=isReversed()?getAngle2():getAngle1();
    for(unsigned short i=0;i<4;i++){
		RS_Line line{vertex.at(i),vertex.at((i+1)%4)};
		auto vpIts=RS_Information::getIntersection(
                    static_cast<RS_Entity*>(this), &line, true);
//    std::cout<<"vpIts.size()="<<vpIts.size()<<std::endl;
        if( vpIts.size()==0) continue;
		for(const RS_Vector& vp: vpIts){
			auto ap1=getTangentDirection(vp).angle();
			auto ap2=line.getTangentDirection(vp).angle();
            //ignore tangent points, because the arc doesn't cross over
            if( fabs( remainder(ap2 - ap1, M_PI) ) < RS_TOLERANCE_ANGLE) continue;
            crossPoints.push_back(
                        RS_Math::getAngleDifference(baseAngle, getEllipseAngle(vp))
                        );
        }
    }
    if(vpStart.isInWindowOrdered(vpMin, vpMax)) crossPoints.push_back(0.);
    if(vpEnd.isInWindowOrdered(vpMin, vpMax)) crossPoints.push_back(
                RS_Math::getAngleDifference(baseAngle,isReversed()?getAngle1():getAngle2())
                );


    //sorting
    std::sort(crossPoints.begin(),crossPoints.end());
    //draw visible
//    DEBUG_HEADER
//    std::cout<<"crossPoints.size()="<<crossPoints.size()<<std::endl;
    RS_Ellipse arc(*this);
    arc.setSelected(isSelected());
    arc.setPen(getPen());
    arc.setReversed(false);
    if( crossPoints.size() >= 2) {
		for(size_t i=1;i<crossPoints.size();i+=2){
			arc.setAngle1(baseAngle+crossPoints[i-1]);
			arc.setAngle2(baseAngle+crossPoints[i]);
            arc.drawVisible(painter,view,patternOffset);
        }
        return;
    }
    //a workaround for buggy equation solver, can be removed, when line-ellipse equation solver is reliable
    if(isVisibleInWindow(view)){
        arc.drawVisible(painter,view,patternOffset);
    }
}

/** directly draw the arc, assuming the whole arc is within visible window */
void RS_Ellipse::drawVisible(RS_Painter* painter, RS_GraphicView* view, double& /*patternOffset*/) {
//    std::cout<<"RS_Ellipse::drawVisible(): begin\n";
//    std::cout<<*this<<std::endl;
	if (!(painter && view)) return;

    //visible in graphic view
	if(!isVisibleInWindow(view)) return;
    double ra(getMajorRadius()*view->getFactor().x);
    double rb(getRatio()*ra);
	if(std::min(ra, rb) < RS_TOLERANCE) {//ellipse too small
        painter->drawLine(view->toGui(minV),view->toGui(maxV));
        return;
    }

    bool drawAsSelected = isSelected() && !(view->isPrinting() || view->isPrintPreview());

    double mAngle=getAngle();
    RS_Vector cp(view->toGui(getCenter()));
	if (!drawAsSelected && (
             getPen().getLineType()==RS2::SolidLine ||
             view->getDrawingMode()==RS2::ModePreview)) {
        painter->drawEllipse(cp,
                             ra, rb,
                             mAngle,
                             getAngle1(), getAngle2(),
                             isReversed());
        return;
    }

    // Pattern:
	const RS_LineTypePattern* pat = nullptr;
	if (drawAsSelected) {
		pat = &RS_LineTypePattern::patternSelected;
	}
	else {
		pat = view->getPattern(getPen().getLineType());
	}

	if (!pat) {
        RS_DEBUG->print(RS_Debug::D_WARNING, "Invalid pattern for Ellipse");
        return;
    }

    // Pen to draw pattern is always solid:
    RS_Pen pen = painter->getPen();
    pen.setLineType(RS2::SolidLine);
    double a1(RS_Math::correctAngle(getAngle1()));
    double a2(RS_Math::correctAngle(getAngle2()));
    if (isReversed()) std::swap(a1,a2);
	if(a2 <a1+RS_TOLERANCE_ANGLE) a2 += 2.*M_PI;
    painter->setPen(pen);
	if(pat->num <= 0){
		RS_DEBUG->print(RS_Debug::D_WARNING,"Invalid pattern when drawing ellipse");
		painter->drawEllipse(cp, ra, rb, mAngle, a1, a2, false);
		return;
	}

	std::vector<double> ds(pat->num, 0.);

	double dpmm=static_cast<RS_PainterQt*>(painter)->getDpmm();
	for (size_t i = 0; i < pat->num; i++) {
		ds[i]= dpmm * pat->pattern[i]; //pattern length
		if(fabs(ds[i]) < 1.)
			ds[i] = copysign(1., ds[i]);
	}

    double curA(a1);
    bool notDone(true);

	//draw patterned ellipse
	for (size_t i=0; notDone; i = (i + 1) % pat->num) {
		double nextA = curA + std::abs(ds[i])/
				RS_Vector(ra*sin(curA), rb*cos(curA)).magnitude();
		if (nextA > a2){
			nextA = a2;
			notDone = false;
        }
		if (ds[i] > 0.)
			painter->drawEllipse(cp, ra, rb, mAngle, curA, nextA, false);

        curA=nextA;
    }

}



/**
 * Dumps the point's data to stdout.
 */
std::ostream& operator << (std::ostream& os, const RS_Ellipse& a) {
    os << " Ellipse: " << a.data << "\n";
    return os;
}