#include <fstream>
#include <iostream>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <map>
using namespace std;

#define PI 3.14159265358979323846264338327950288419716939937510
#define VACIO 0.00000125663706144

int main(int argc, char** args)
{
	char ignorarC;
	double ignorarD,cenX,cenY,cooX,cooY,radio,theta,Bx,By,t11,t22,t12,dfx,dfy,dP,Rx,Ry,
	*ffunction,*ftheta,H,SumEven=0,SumOdd=0,SUM=0;
	map<double,double> thetaForce;
	map<double,double>::iterator it;
	
	/*Archivo de entrada fileBxy esta generado con GMSH. Se obtiene a partir de
	una distribucion de campo magnetico. Se genera una circunferencia con
	el plugin cutparametric donde se mostrara los valores del campo en cada punto,
	teniendo cuidado de hacer el barrido desde -3,1415 hasta 3,1415.
	Imprimimos la vista obtenida en un archivo RAW, que sera el que utilice este
	algoritmo.*/
	ifstream fileBxy;
	ofstream fileDF; // archivo de salida visualizable en GMSH.

	// Leyendo -x -y -f de la linea de comandos
	cout << "Usar: max.bin -x coordenada -y coordenada -f archivo\n\n";
	for(int i=0; i<argc; i++ ){
		if(strcmp(args[i],"-x")==0){
			cenX=atof(args[i+1]);
		}
		if(strcmp(args[i],"-y")==0){
			cenY=atof(args[i+1]);
		}
		if(strcmp(args[i],"-f")==0){
			fileBxy.open(args[i+1], ios::in);
			if( !fileBxy.is_open() ){ cout<<"maxwel: Error al leer \"-f\"\n";exit(0);}
		}
	}

	fileDF.open("fuerzas.txt", ios::out);
		if( !fileDF.is_open() ){ cout<<"maxwel: Error al crear archivo.\n";exit(0);}


// comenzamos a escribir en los archivos
	fileDF << "View \"f\" {\n";

	// leyendo el archivo "gmsh raw" y calculando dfx, dfy y dpar para cada punto
	while (!fileBxy.eof( )){
		fileBxy >> cooX >> cooY >> ignorarD >> Bx >> By >> ignorarD;

		Rx=cooX-cenX;
		Ry=cooY-cenY;

		t11=0.5*(Bx*Bx-By*By);
		t22=0.5*(By*By-Bx*Bx);
		t12=Bx*By;

		// Ojo, que el barrido de la circunferencia de puntos en el gmsh hay que hacerlos desde
		// -3.1415 hasta 3.1415, para que el incremento de angulo sea cte en todo el barrido
		// calculando theta y radio corregidos respecto al centro del circulo de valores
		theta=atan2(Ry,Rx);
		radio=sqrt(Rx*Rx+Ry*Ry);

		// calculando el valor del interior de la integral para cada paso de angulo
		dfx = (t11*radio*cos(theta) + t12*radio*sin(theta))/VACIO;
		dfy = (t12*radio*cos(theta) + t22*radio*sin(theta))/VACIO;
		dP = dfy*Rx - dfx*Ry;
		thetaForce[theta]=dfx; // con lo que se calculara la fuerza en el eje X. podría haber sido tambien sobre dfy para calcular la fuerza en la direccion Y

		fileDF << "VP(" << cooX << ',' << cooY << ",0){";
		fileDF << dfx << ',' << dfy << ",0};" << endl;

	}

	fileDF << "};";

	fileDF.close();

	fileBxy.close();

	// pasando la variable map a puntero a doble, para mayor control sobre los indices

	ffunction=new double[thetaForce.size()];
	ftheta=new double[thetaForce.size()];

	int i=0;
	for ( it=thetaForce.begin() ; it != thetaForce.end(); it++,i++ ){
		ftheta[i]=(*it).first;
		ffunction[i]=(*it).second;
	}
	thetaForce.clear();
	i=i-1; // el ultimo valor es repetido
	H=ftheta[i-1]-ftheta[i-2]; //variacion de los dos ultimos valores de theta. Solo es para  saber el incremento entre medidas


 // METODO DE SIMPSON //
///////////////////////

//	cout << "Pares:\n"; (2,4,6,8...n-2)
    for ( int K = 2; K <= i-2; K=K+2 ) {
         SumEven = SumEven + ffunction[K];
    }

//	cout << "Impares:\n"; (1,3,5,7...n-1)
    for ( int K = 1; K <= i-1; K=K+2 ) {
         SumOdd = SumOdd + ffunction[K];
    }

	SUM = H * ( ffunction[0] + ffunction[i] + 2*SumEven + 4*SumOdd ) / 3;

	cout << "Mediante SIMPSON:\nR=" << radio << endl;
	cout << "i: " << i << " ; H: " << H << endl;
	cout << "ffunction[0]: " << ffunction[0] << "; ffunction[i]: " << ffunction[i] << endl;; 
	cout << "2*SumEven: " << 2*SumEven << "; 4*SumOdd: " << 4*SumOdd << endl;
	cout << "Valor de la integral definida: " << SUM << endl;


 // METODO DE TRAPEZOIDAL //
///////////////////////////

	SUM=0;
	for ( int K = 1; K <= i-1; K++ ) {
		SUM = SUM + ffunction[K];
	}
	SUM = H * ( ffunction[0] + ffunction[i] + 2*SUM )/2;

	cout << "Mediante TRAEZOIDES:\nR=" << radio << endl;
	cout << "i: " << i << " ; H: " << H << endl;
	cout << "ffunction[0]: " << ffunction[0] << "; ffunction[i]: " << ffunction[i] << endl; 
	cout << "Valor de la integral definida: " << SUM << endl;

	return 0;
}

