/*
    Hyperbolic Tessellations
    Copyright (C) 2007 Dmitry Brant
    http://dmitrybrant.com

    Based largely on code from
    http://www.plunk.org/~hatch/HyperbolicApplet/
*/
//---------------------------------------------------------------------------
#ifndef IsometryH
#define IsometryH

#include "mycomplex.h"


class Isometry{
public:
    complex T, P;
    int R;

    Isometry(const complex& _T, const complex& _P, int _R=1):T(_T),P(_P),R(_R) { }
    //------------------------------------------------------------

    Isometry():
    T(0.,0.), P(0.,0.), R(1)
    { }
    //------------------------------------------------------------

    Isometry(const Isometry& from):T(from.T),P(from.P),R(from.R) { }
    //------------------------------------------------------------

    inline
    Isometry& operator=(const Isometry& from){
        T = from.T;
        P = from.P;
        R = from.R;
        return *this;
    }

    inline
    void apply(complex& z, complex& result){
        if (R < 0)
            result = (T*(~z) + P) / ((~P)*T*(~z) + 1.0);
        else
            result = (T*z + P) / ((~P)*T*z + 1.0);
    }
    //------------------------------------------------------------

    inline
    Isometry operator*(const Isometry& _F0){
        Isometry F0(_F0);
        if (this->R < 0){
            F0.T = ~F0.T;
            F0.P = ~F0.P;
        }
        complex denom = this->T*F0.P*(~this->P) + 1.0;
        Isometry result;
        result.R = this->R * F0.R;
        result.P = (this->T*F0.P + this->P) / denom;
        result.T = (F0.T*this->T + F0.T*(~F0.P)*this->P) / denom;
        // renormalize, as recommended in the paper.
        // use the fact that sqrt(1+eps) is approximately 1+eps/2 for small eps.
        double invLenT = 1.0 / ((1+result.T[0]*result.T[0]+result.T[1]*result.T[1])*0.5);
        result.T[0] *= invLenT;
        result.T[1] *= invLenT;
        return result;
    }

    //---------------------------------------------------------


    inline
    static Isometry pow(Isometry& F, int e){
        if (e > 0){
            Isometry temp = pow(F, e/2);
            if ((e % 2) == 0)
                return (temp*temp);
            else
                return ((F*temp)*temp);
        }
		else if (e < 0){
			Isometry temp = F.inverse();
			return pow(temp, -e);
        }
        else{
            return Isometry(one,zero,1);
        }
    }
    //------------------------------------------------------------


    inline
    Isometry inverse(){
        Isometry F0(one,-P), F1(~T,zero), F2(one,zero,R);
        return F2 * F1 * F0;
    }
    //------------------------------------------------------------


    inline
    static Isometry pureRotation(double angle){
        return Isometry(complex(cos(angle),sin(angle)), zero, 1);

    }
    //------------------------------------------------------------

    // the simple isometry that takes 0,0 to x1,y1, and takes -x1,-y1 to 0,0
    inline
    static Isometry pureTranslation(double x1, double y1){
        return Isometry(one, complex(x1,y1), 1);
    }
    //------------------------------------------------------------

    inline
    static Isometry pureTranslation(complex& p0, complex& p1){
        complex A = p1-p0;
        complex B = p1*p0;
        double denom = 1.0 - (B[0]*B[0] + B[1]*B[1]);
        return Isometry(one, complex((A[0]*(1+B[0])+A[1]*B[1])/denom,
                                                   (A[1]*(1-B[0])+A[0]*B[1])/denom),
                                                    1);
    }
    //------------------------------------------------------------

    inline
    static Isometry reflectionAcrossLine(complex& l0, complex& l1){
        Isometry takel0ToOrigin = Isometry::pureTranslation(-l0[0],-l0[1]);
        Isometry takeOriginTol0 = Isometry::pureTranslation(l0[0],l0[1]);
        complex l1_;
        takel0ToOrigin.apply(l1, l1_);
        double angle = atan2(l1_[1], l1_[0]);
        Isometry rotatel1_ToXAxis = Isometry::pureRotation(-angle);
        Isometry rotateXAxisTol1_ = Isometry::pureRotation(angle);
        Isometry reflectAboutXAxis(one, zero, -1);
        Isometry result = takel0ToOrigin;
        result = rotatel1_ToXAxis * result;
        result = reflectAboutXAxis * result;
        result = rotateXAxisTol1_ * result;
        result = takeOriginTol0 * result;
        return result;
    }
    //------------------------------------------------------------

};

//-----------------------------------------------------------




inline

bool operator==(const Isometry& first, const Isometry& other){

    // XXX clean up-- this is error-prone!
    double eps=1e-6;
    return (first.R == other.R
        && first.T[0]-other.T[0] <= eps
        && other.T[0]-first.T[0] <= eps
        && first.T[1]-other.T[1] <= eps
        && other.T[1]-first.T[1] <= eps
        && first.P[0]-other.P[0] <= eps
        && other.P[0]-first.P[0] <= eps
        && first.P[1]-other.P[1] <= eps
        && other.P[1]-first.P[1] <= eps);
}
//----------------------------------------------------------





//---------------------------------------------------------------------------
#endif