#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;
	ifstream fileBxy;
	ofstream fileDF, fileDP; // poner las novedades para estos ficheros.

	// 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]=dP; // 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;
}

