/*
*                   NUMERICAL INTEGRATION
*/

#include<iostream.h>
#include<math.h>
#include"lib351.h"

void main(){

	//____DATA_DESCRIPTION____
	bool goout = false;
    ofstream dataOutStr;
	int nCharges;
	float sizeA, x0, y0;
	float x, deltaX;
	float f, v1, v2, v3, v4;
	float lambda, c;
	
	outFileMaker("m349p1c.out", &dataOutStr);
	dataOutStr<<'#'<<"NUMERICAL INTEGRATION"<<endl;

	//_____CYCLE_____
	while (!goout){
		//____DATA_INPUT____
		cout<<"Please, supply values for nCharges, sizeA, x0, y0"<<endl;
		cout<<"number of point charges, nCharges:"<<endl;
		cin>>nCharges;
		cout<<"\nhalf-lenght of the line (-A, A) along x-axis, sizeA:"<<endl;
		cin>>sizeA;
		cout<<"\ncoordinates of point"<<endl;
		cout<<"\nx0:"<<endl;
		cin>>x0;
	  	cout<<"\ny0:"<<endl;
		cin>>y0;
	    dataOutStr<<endl<<endl<<"nCharges = "<<nCharges;
	    dataOutStr<<"; sizeA = "<<sizeA<<"; x0 = "<<x0<<"; y0 = "<<y0<<endl;
	    if(x0 == 0){x0 = float(.00001);}
	    if(y0 == 0){y0 = float(.00001);}
		nCharges = 2 * floor(nCharges/2) + 1;
		//____CALCULATION____
		v1 = log((sqrt(pow(x0-sizeA, 2) + pow(y0, 2)) - x0 + sizeA)/
	  	         (sqrt(pow(x0+sizeA, 2) + pow(y0, 2)) - x0 - sizeA));
		deltaX = 2 * sizeA/(nCharges - 1);
		v2 = 0; v3 = 0; v4 = 0;
		x = -sizeA;
		for (int i = 1; i <= nCharges; i++){
			lambda = 1;
			f = lambda/sqrt(pow(x-x0, 2) + y0*y0);
			c = 1;
			v2 += c * f;
	         	c = 1;
			if(i == 1 || i == nCharges){c = 1./2;}
			v3 += c * f;
			c = 2;
			if(i == 2 * floor(i/2)){c = 4;}
	   		if(i == 1 || i == nCharges){c = 1;}
	   		v4 += c * f;
			x += deltaX;
		}
		v2 *= deltaX;
		v3 *= deltaX;
		v4 *= deltaX/3;
		//_____DATA_OUTPUT_____
		dataOutStr<<"\nv1 ="<<v1;
		dataOutStr<<"\nv2 ="<<v2;
		dataOutStr<<"\nv3 ="<<v3;
	  	dataOutStr<<"\nv4 ="<<v4<<endl;
		goout = query();
	}
};