#include <iostream.h>;
#include <fstream.h>;
#include <math.h>;
//constants

//GAMMA electricity constant
const double GAMMA = 8987551787.3681764;
//e elementary charge
const double e = 1.602*pow(10,-19);


//n = number of electric charges
const int n = 2;

const float Q[n] = { pow(10,-11), pow(10,-11) };

//----------------------
/*const float gridx = 0.01;
const float gridy = 0.01;
*/
//const double scale = 0.5*gridx;

const float X[n] = { 0, 0 };

const float Y[n] = { 0, 0.055 };

const float Z[n] = { 0 ,0 };

float z = 0;





/////////////////////////////////////////////////////////////////////
//------------------------------------------------------------------
//definition for potential of n electric charges PT( , , )

double PT(float XP, float YP, float ZP)
{
  double PT = 0;
  for (int i=0; i<n; i++)
	 {
double
	  deltai =  pow(  (XP-X[i]), 2) + pow( (YP-Y[i]), 2) + pow( (ZP-Z[i]), 2);
	  PT = PT + GAMMA*Q[i]/(sqrt(deltai));
	 }
return PT;
}
//------------------------------------------------------------------
//
//definition for ai function of n electric charges a(i, , , )
double a(int i,float XP, float YP, float ZP)
{
double
deltai =  sqrt(pow(  (XP-X[i]), 2) + pow( (YP-Y[i]), 2) + pow( (ZP-Z[i]), 2));
double aa = GAMMA*Q[i]/( pow(deltai,3) );
return aa;
}


//------------------------------------------------------------------
//definition for magnitude "E" of n electric charges Mag( , , )

double Mag[3];    //Magnitude vector
double MagVrs[3]; //Magnitude versor
GetMag(float XP, float YP, float ZP)
{
  double magx = 0;
  double magy = 0;
  double magz = 0;
  for (int i=0; i<n; i++)
	 {
		double aaa = a(i,XP,YP,ZP);
		magx = magx + aaa*(XP-X[i]);
		magy = magy + aaa*(YP-Y[i]);
		magz = magz + aaa*(ZP-Z[i]);
	 }
  Mag[0] = magx;
  Mag[1] = magy;
  Mag[2] = magz;

double ModMag = sqrt( Mag[0]*Mag[0] + Mag[1]*Mag[1] + Mag[2]*Mag[2] );


  MagVrs[0] = magx/ModMag;
  MagVrs[1] = magy/ModMag;
  MagVrs[2] = magz/ModMag;


}






//------------------------------------------------------------------

void OutFileLine(char* filename, int cnum , int steps, int accur)
{

ofstream out(filename);

/*  out << endl << "Current lines for charges";
  for (int i=0; i<n; i++)
		{

		 out << endl << "Q["<<i<<"] ="<<Q[i]
			  <<"  ("<<X[i]<<","<<Y[i] <<","<< Z[i]<<")";
		}
		out <<endl
		<<endl<<"gridx = "<<gridx
		<<endl<<"gridy = "<<gridy

		<<endl<<"start point"
		<<endl<<"near Q[ "<<cnum<<"]"
		<<endl<<"z value = "<<z
		<<endl;
*/

float xstart, ystart, zstart;
int spnum = 50;

  for (int j = 0; j < spnum; j++)
	{
	double reduce =0.01;
	xstart = X[cnum] + reduce*cos(2*3.14*j/(spnum));
	ystart = Y[cnum] + reduce*sin(2*3.14*j/(spnum));
	zstart = Z[cnum];
//	+ 0.01*gridx;

//		out<<"LINE";
		out<<"3DPOLY";
	for (int jj=0; jj < steps; jj++)
		  {

			GetMag(xstart, ystart, zstart);

if ( (10*floor(jj/10) - jj) == 0 ) {out <<endl<<xstart<<","<<ystart<<","<<zstart;}
//			out <<endl<<xstart<<","<<ystart<<","<<zstart;

		xstart = xstart + (Q[cnum]/fabs(Q[cnum]))*MagVrs[0]/accur;
		ystart = ystart + (Q[cnum]/fabs(Q[cnum]))*MagVrs[1]/accur;
		zstart = zstart + (Q[cnum]/fabs(Q[cnum]))*MagVrs[2]/accur;


		 }
out <<endl<<endl;
	}

}

//---------------------------------------------------------------



//void OutSingleLine(char* filename, int cnum , int steps, int accur,
//double xstart, double ystart, double zstart)

void OutSingleLine(char* filename, int steps, int accur,
double xstart, double ystart, double zstart)

{

ofstream out(filename);


//float xstart, ystart, zstart;


		out<<"3DPOLY";
	for (int jj=0; jj < steps; jj++)
		  {

			GetMag(xstart, ystart, zstart);

if ( (10*floor(jj/10) - jj) == 0 ) {out <<endl<<xstart<<","<<ystart<<","<<zstart;}
//			out <<endl<<xstart<<","<<ystart<<","<<zstart;

//		xstart = xstart + (Q[cnum]/fabs(Q[cnum]))*MagVrs[0]/accur;
//		ystart = ystart + (Q[cnum]/fabs(Q[cnum]))*MagVrs[1]/accur;
//		zstart = zstart + (Q[cnum]/fabs(Q[cnum]))*MagVrs[2]/accur;

		xstart = xstart + MagVrs[0]/accur;
		ystart = ystart + MagVrs[1]/accur;
		zstart = zstart + MagVrs[2]/accur;

		 }
out <<endl<<endl;


}

//---------------------------------------------------------------



void OutSingleIzoPotLine(char* filename, int steps, int accur,
double xstart, double ystart, double zstart)

{

ofstream out(filename);


//float xstart, ystart, zstart;


		out<<"3DPOLY";
	for (int jj=0; jj < steps; jj++)
		  {

			GetMag(xstart, ystart, zstart);

if ( (10*floor(jj/10) - jj) == 0 ) {out <<endl<<xstart<<","<<ystart<<","<<zstart;}
//			out <<endl<<xstart<<","<<ystart<<","<<zstart;

//		xstart = xstart + (Q[cnum]/fabs(Q[cnum]))*MagVrs[0]/accur;
//		ystart = ystart + (Q[cnum]/fabs(Q[cnum]))*MagVrs[1]/accur;
//		zstart = zstart + (Q[cnum]/fabs(Q[cnum]))*MagVrs[2]/accur;

		xstart = xstart - MagVrs[1]/accur;
		ystart = ystart + MagVrs[0]/accur;
		zstart = zstart + MagVrs[2]/accur;

		 }
out <<endl<<endl;


}

//---------------------------------------------------------------

void GetPositions(char* filename)
{
ofstream out(filename);

  out << "Positions of charges"<< endl ;
  for (int i=0; i<n; i++)
		{
  for (int ii=i+1; ii<n; ii++)
		{
		 out <<endl<<"LINE"<<endl
			  <<X[i]<<","<<Y[i] <<","<< Z[i]<<endl
			  <<X[ii]<<","<<Y[ii] <<","<< Z[ii]<<endl;
		}
		}
}

//---------------------------------------------------------------

//////////////////////////////////////////////////////////////////
main()
{



//OutFileLine(char* filename, int cnum , int steps, int accur)
//OutFileLine("l0.scr",0, 500,10000);

//void OutSingleLine(char* filename, int cnum , int steps, int accur,
//double xstart, double ystart, double zstart)

OutSingleLine("ls.scr", 500,10000,
0.01, 0.01, 0);


OutSingleIzoPotLine("li.scr", 5000,100000,
0, 0.06, 0);

};