#include <cmath>
#include <iostream>
#include <fstream>
#include <sys/stat.h>
#include <sys/types.h>
#include <sstream>

#include "vec3.h"
#include "curve.h"

using namespace std;

vec3 integrand(vec3 s, vec3 ds)
{
    cout << "point " << s << "ds = " << ds << endl;

    return ds;
}

vec3 point;

/* dl x r / (|r| ^ 2) */
vec3 dB(vec3 s, vec3 ds)
{
    vec3 r = s - point;

    scalar r2 = r.magnitudeSquared();

    vec3 rnorm = r / std::sqrt(r2);

    return ds.cross(rnorm) / r2;
}

/* dl * r / (|r| ^ 2) */
vec3 dE(vec3 s, vec3 ds)
{
    vec3 r = point - s;

    scalar r2 = r.magnitudeSquared();

    vec3 rnorm = r / std::sqrt(r2);

    return rnorm * ds.magnitude() / r2;
}

/* A current or charge distribution */
struct Entity {
    enum { CHARGE, CURRENT } type;
    
};

Arc loop(vec3(0, 0, 0), vec3(0, 1, 0), vec3(1, 0, 0), M_PI * 2);
//Spiral loop(vec3(0, 0, 0), vec3(0, 1, 0), vec3(1, 0, 0), M_PI * 2 * 10, 1);
//LineSegment loop(vec3(-.1, .1, 0), vec3(.1, .1, 0));
//Toroid loop;

ostream *dump_ofs = NULL;

vec3 dump(vec3 s, vec3 ds)
{
    *dump_ofs << s << " " << ds << endl;
    //cout << "Magn: " << ds.magnitude() << endl;

    return 0;
}

void dump_path(ostream &out, const Curve *c)
{
    dump_ofs = &out;
    c->integrate(dump, 1e-2);
}

/* dump the field (gnuplot format) at z = 0 */
/* requires x0 < x1, y0 < y1 */

const scalar U0 = 4e-7 * M_PI;
const scalar I = 1;
const scalar Q = 1;
const scalar C = 299792458;
const scalar E0 = 1 / ( U0 * C * C );
const scalar K_E = 1 / (4 * M_PI * E0);

void dump_field(scalar x0, scalar y0, scalar z0, scalar x1, scalar y1, scalar z1, scalar delta)
{
    for(scalar z = z0; z <= z1; z += delta)
        for(scalar y = y0; y <= y1; y += delta)
            for(scalar x = x0; x <= x1; x += delta)
            {
                point = vec3(x, y, z);

                vec3 B = loop.integrate(dE, 1e-1) * Q * K_E;

                if(B.magnitude() > 1e-8)
                {
                    B=B.normalize() / 10;
                    cout << point[0] << " " << point[1] << " " << point[2] << " ";
                    cout << B[0] << " " << B[1] << " " << B[2] << endl;
                }
            }
}

void dump_fieldline(ostream &out, vec3 x, scalar len)
{
    point = x;
    scalar delta = .1;
    while(len > 0)
    {
        out << point << endl;
        
        vec3 B = loop.integrate(dB, 1e-1) * U0 * I;

        point += delta * B;
        len -= delta;
    }
}

void dump_values(vec3 start, vec3 del, int times)
{
    point = start;
    while(times--)
    {
        
        
        point += del;
    }    
}

int main(int argc, char *argv[])
{
    //if(argc != 2)
        //    return 1;
    //loop = Toroid(vec3(0, 0, 0), vec3(1, 0, 0), vec3(0, 0, 1), 1 * M_PI, .1, atof(argv[1]));
    
    
    //LineSegment wire(vec3(0, -100, 0), vec3(0, 100, 0));

    //std::cout << "length = " << loop.integrate(integrand, 1e-2) << endl;

    //vec3 dx = .01;

    //point = 0;

    //scalar I = 1;

    //for(int i = 0; i < 1000; i++, point += dx)
    //std::cout << point[0] << " " << U0 / ( 4 * M_PI ) * loop.integrate(dB, 1e-2)[0] << endl;
    
   dump_field(-3, -3, 0,
              3, 3, 0,
              .1);
    //dump_field(0,0,0,0,0,0);

    stringstream ss;
    ss << "curve.fld";
    ofstream ofs(ss.str());

    dump_path(ofs, (Curve*)&loop);
    
    ofs.close();
    
    
#if 0
    mkdir("field", 0755);

    for(scalar y = -1.5; y <= 1.5; y += .1)
    {
        stringstream ss;
        ss << "field/" << y << ".fld";
        ofstream ofs(ss.str());

        dump_fieldline(ofs, vec3(0, y, 0), 10);

        ofs.close();
    }
#endif
}