#include "bigfloat.h"
#include "ringforce.h"

struct BigPoint
{
    void Init(const BigFloat& x, const BigFloat& y, const BigFloat& z)
    {
        _x.Copy(x);
        _y.Copy(y);
        _z.Copy(z);
    }

    void Copy(BigPoint p)
    {
        Init(p._x, p._y, p._z);
    }


    void Sub(BigPoint p)
    {
        _x.Sub(p._x);
        _y.Sub(p._y);
        _z.Sub(p._z);
    }

    BigFloat Distance(BigPoint p)
    {
        BigFloat d;
        BigFloat t;
        d.Add(t.Copy(_x).Sub(p._x).Mult(t));
        d.Add(t.Copy(_y).Sub(p._y).Mult(t));
        d.Add(t.Copy(_z).Sub(p._z).Mult(t));
        return d.Sqrt();
    }

    void Show()
    {
        printf("%.18g %.18g %.18g", _x.ToDouble(), _y.ToDouble(), _z.ToDouble());
    }
    
    BigFloat _x;
    BigFloat _y;
    BigFloat _z;
};


class BigRing
{
public:
    // build a ring passing through p, in the (x,y) plane and z=0, mass 1, with nMoons moons in it
    BigRing(s8 nMoons, bool odd, BigPoint p)
    {
        _nMoons = nMoons;
        _odd = odd;
        _ring = new BigPoint[_nMoons];
        BigFloat x2(p._x);
        x2.Mult(p._x);
        BigFloat y2(p._y);
        y2.Mult(p._y);
        BigFloat scale(x2);
        scale.Add(y2);
        scale.Sqrt();
        // make all the points for a ring of 
        for (s8 i=0; i<nMoons; i++)
        {
            BigFloat x(2*i+odd);
            x.Mult(2).Mult(BigFloat::Pi()).Div(2*nMoons);
            BigFloat y(x);
            y.Sin().Mult(scale);
            x.Cos().Mult(scale);
            _ring[i].Init(x, y, p._z);
        }
    }

    ~BigRing()
    {
        delete[] _ring;
    }

    // find the force exerted by a unit ring of unit mass on a point exactly by nbody calculations
    BigPoint NBodyForce(const BigPoint& p) const
    {
        // find the acceleration it exerts on the point of reference
        BigPoint force;
        for (s8 i=0; i<_nMoons; ++i)
        {
            // do not measure the attraction of a point to itself
            if (_ring[i]._x.Compare(p._x) == 0 &&
                _ring[i]._y.Compare(p._y) == 0 &&
                _ring[i]._z.Compare(p._z) == 0)
                continue;
            
            BigPoint q(_ring[i]);
            q.Sub(p);
            BigFloat t;
            BigFloat d;
            d.Add(t.Copy(q._x).Mult(t));
            d.Add(t.Copy(q._y).Mult(t));
            d.Add(t.Copy(q._z).Mult(t));
            t.Copy(d);
            t.Mult(d.Sqrt());
            force._x.Add(q._x.Div(t));
            force._y.Add(q._y.Div(t));
            force._z.Add(q._z.Div(t));
        }
        BigFloat scale(_nMoons);
        scale.Invert();
        force._x.Mult(scale);
        force._y.Mult(scale);
        force._z.Mult(scale);
        return force;
    }

    
private:
    s8 _nMoons;
    BigPoint* _ring;
    bool _odd;
};



void Driver(s8 nMoons)
{
    // the ring
    BigPoint br;
    br.Init(BigFloat(3), BigFloat(4), BigFloat(2));
    
    // the point being acted upon
    BigPoint bp;
    bp.Init(BigFloat(3), BigFloat(4), BigFloat(2));
    
    // make all the points for a ring of nMoons moons
    BigRing oddRing(nMoons, true, br);
    BigRing evenRing(nMoons, false, br);

    // find the acceleration the ring through br exerts on the point bp
    BigPoint accelEven = oddRing.NBodyForce(bp);
    BigPoint accelOdd = evenRing.NBodyForce(bp);

    // report the result
    accelOdd.Show();
    printf("  diff %g\n", accelEven.Distance(accelOdd).ToDouble());

    Point r(br._x.ToDouble(), br._y.ToDouble(), br._z.ToDouble());
    Point p(bp._x.ToDouble(), bp._y.ToDouble(), bp._z.ToDouble());
    Point a = RingAttraction(r, p);

    printf("%.18g %.18g %.18g\n", a._x, a._y, a._z);
    printf("diff %.18g %.18g %.18g\n",
           a._x - accelOdd._x.ToDouble(), a._y - accelOdd._y.ToDouble(), a._z - accelOdd._z.ToDouble());
}

int main(int argc, char** argv)
{
    BigFloatCache::Init();
    
    ASSERT(argc == 2);
    s8 nMoons;
    sscanf(argv[1], "%ld", &nMoons);
    printf("nmoons: %d\n", nMoons);

    try
    {
        Driver(nMoons);
    }
    catch( const std::exception & ex )
    {
        fprintf(stderr, "%s\n", ex.what());
    }
    return 0;
}