#include <iostream.h>;
#include <fstream.h>;
#include <mapping.h>;
#include <stdlib.h>;
#include <math.h>;


//------------------------------------------------------------------
//definition of Sequence class functions
//------------------------------------------------------------------


ChargeSet::RenumCharges()
{
	int num = ChargesNum;
	const Link * pc = TList<LCharge>::GetHead();
	while ( pc != 0 )
	{
		((LCharge*)pc)->SetId(num);
		num--;
		pc = ( pc -> Next() );
	}


}



ChargeSet::AddCharge(Charge qq)
{
	TList<LCharge>::Add();
	ChargesNum++;
	LCharge * ptr = TList<LCharge>::Get();
	(ptr->x) = qq.x; (ptr->y) = qq.y; (ptr->z) = qq.z; (ptr->q) = qq.q;
	(ptr->Id) = ChargesNum;
}

ChargeSet::RemoveCharge(int n)
{
	TList<LCharge>::Remove(n);
	ChargesNum--;
	RenumCharges();
}

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

double ChargeSet::Potential(double x, double y, double z)
{
//GAMMA electricity constant
 double GAMMA = 8987551787.3681764;

   double PT = 0;
   const Link * cp = TList<LCharge>::GetHead();

	while ( cp != 0 )
	{
	 double cx = (((LCharge*)cp) ->x);
	 double cy = (((LCharge*)cp) ->y);
	 double cz = (((LCharge*)cp) ->z);
	 double cq = (((LCharge*)cp) ->q);

	 double delta = sqrt( pow( (x-cx), 2)
			     +pow( (y-cy), 2)
			     +pow( (z-cz), 2 ));

	 PT = PT + (GAMMA*cq)/(delta);

		cp = ( cp -> Next() );
	}
return PT;
}

////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////
void ElectrodeArray::SetInternalModel()
{
//Fukuda mapping

double pi = 3.1415;
double x1;
double fi;
double dx=4.5;  //vertical distance beetwen electrodes
//double dfi=(double)(2*pi/Meshy());


double xFactor = 0;
double y1,z1;
double xStart=-18;



// mapping designations   A   B      C       D       E
double angles[13] =      {0, pi/8, 2*pi/8, 3*pi/8, 4*pi/8,
//                           F       G       H       I
			   5*pi/8, 6*pi/8, 7*pi/8,   pi,
//                           J       K       L       M
			   6*pi/5, 7*pi/5, 8*pi/5, 9*pi/5};

////////////
Electrode * ArrPtr = Ptr();
//////////

x1=xStart+3; //start -x- for electrodes setting

int HorizCounter=0;
int VertCounter=0;

///start of main loop
int PointCounter = 0;

for (int i=0; i<7; i++) //x direction for breast, vertical for voltages
{
VertCounter++;

  if (x1 - xStart<2) {xFactor = (2-x1+xStart)/2;}
  if (x1 - xStart>=2  && x1 - xStart<11) {xFactor = (11-(x1 - xStart))/9;}
  if (x1 - xStart>=11 && x1 - xStart<18) {xFactor = (18-(x1 - xStart))/7;}
  if (x1 - xStart>=18 && x1 - xStart<25) {xFactor = (25-(x1 - xStart))/7;}
  if (x1 - xStart>=25 && x1 - xStart<37) {xFactor = (37-(x1 - xStart))/12;}


  HorizCounter=0;
  for (int j=0; j<13; j++) //y direction for breast, horizontal for voltages
      {
       HorizCounter++;
       fi = angles[j];
      //5 - 6 section
      if (x1-xStart<2) {
y1 = 22*cos(fi)-6*cos(fi)*cos(fi)*cos(fi)
     -xFactor*(-6*cos(fi)*cos(fi)*cos(fi) + 14.75*cos(fi));
z1 = 7*sin(fi) -1 - xFactor*(sin(fi) );
		 }

      //4 - 5 section
      if (x1 - xStart>=2 && x1 - xStart<11) {
y1 = 22*cos(fi)-6*cos(fi)*cos(fi)*cos(fi);
z1 = 11.75*sin(fi) - xFactor*(3.75*sin(fi)+1);
		 }

      //3 - 4 section
      if (x1-xStart>=11 && x1-xStart<18) {
y1 = 21*cos(fi)-6*cos(fi)*cos(fi)*cos(fi) + xFactor*cos(fi);
z1 = 12.5*sin(fi) - xFactor*0.75*sin(fi);
		 }


      //2 - 3 section
      if (x1 - xStart>=18 && x1 - xStart<25) {
y1 = 19.5*cos(fi) - 6*cos(fi)*cos(fi)*cos(fi) - cos(fi)*sin(fi)
   + xFactor*(cos(fi)*sin(fi) + 1.5*cos(fi));
z1 = 11.5*sin(fi) - 1 + xFactor*(sin(fi) +1);
		 }



      //1 - 2 section
      if (x1 - xStart>=25 && x1 - xStart<37) {
      y1 = 12.5*cos(fi)
       + xFactor*(7*cos(fi) - 6*cos(fi)*cos(fi)*cos(fi) - cos(fi)*sin(fi));
      z1 = 9*sin(fi) + xFactor*(2.5*sin(fi)-1);
		 }


       ArrPtr[PointCounter].x = x1;
       ArrPtr[PointCounter].y = y1;
       ArrPtr[PointCounter].z = z1;
       ArrPtr[PointCounter].hIdSet(1+j);
       ArrPtr[PointCounter].vIdSet(7-i);
       ArrPtr[PointCounter].SetStatus(1);

       PointCounter++;
      }
x1=x1+dx;
}
//seting inactive (status = 0) additional 4 electrodes
ArrPtr[0].SetStatus(0);
ArrPtr[8].SetStatus(0);
ArrPtr[13].SetStatus(0);
ArrPtr[21].SetStatus(0);
}
////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////


//------------------------------------------------------------------
//definition of Mapping class functions
//------------------------------------------------------------------

Mapping::RenumSequences()
{
	int num = SequenceNum;
	const Link * pc = TList<LSequence>::GetHead();
	while ( pc != 0 )
	{
		((LSequence*)pc)->SetId(num);
		num--;
		pc = ( pc -> Next() );
	}


}


Mapping::AddSequence()   //adding next empty Sequence
{
	TList<LSequence>::Add();
	SequenceNum++;
	LSequence * ptr = TList<LSequence>::Get();
	ptr->SetId(SequenceNum);
	SetCurSequence(SequenceNum);
}


Mapping::RemoveSequence()
{
	int next = CurSequence - 1;
	TList<LSequence>::Remove(CurSequence);
	SequenceNum--;
	RenumSequences();
	SetCurSequence(next); // ????????????????????????
}


void Mapping::SetCurSequence(int id)
{
if (  (1 <= id)&(id<=SequenceNum) ){CurSequence = id;}
}


Sequence * Mapping::GetCurSequence()
{
   const Link * pc = GetHead();
   while ( pc != 0 )
     {
      if (   ((LSequence*)pc)->Id() == CurSequence )
	 {   return (LSequence*)pc;   }
      pc = ( pc -> Next() );
     }
  return 0;
}


Mapping::CalculateCurSequence()
{

Electrode * eArrPtr = ElectrodeArray::Ptr();
Sequence  * cseq = GetCurSequence();
Voltage   * vArrPtr = cseq -> Ptr();

for(int i = 0; i < ElectrodeArray::GetSize(); i++)
   {
    double mx = 0.01*(eArrPtr[i].x);
    double my = 0.01*(eArrPtr[i].y);  //scaling in [m]
    double mz = 0.01*(eArrPtr[i].z);

    double volt =  ( cseq ->Potential(mx, my, mz) - GND );

    vArrPtr[i].USet(1000000*volt); // scaling in [muV]

    vArrPtr[i].hIdSet(eArrPtr[i].hId());
    vArrPtr[i].vIdSet(eArrPtr[i].vId());
    vArrPtr[i].SetStatus(1);

   //seting inactive (status = 0) additional 4 voltages
   vArrPtr[0].SetStatus(0);
   vArrPtr[8].SetStatus(0);
   vArrPtr[13].SetStatus(0);
   vArrPtr[21].SetStatus(0);


   }
}
//////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////

void Manager::SaveElectrodes(const Mapping& mapp,
		 char *filename,char *format)
{
Electrode * ArrPtr = mapp.ElectrodeArray::Ptr();
ofstream out(filename);

for(int i = 0; i< mapp.ElectrodeArray::GetSize(); i++)
   {
    if( ArrPtr[i].GetStatus() != 0 )
      {
	 out<<"POINT"<<endl;
	 out<<(ArrPtr[i]).x<<","
	     <<(ArrPtr[i]).y<<","
	     <<(ArrPtr[i]).z<<endl;
      }
   }
}




void Manager::SaveCurSequence(const Mapping& mapp,
		 char *filename,char *format)
{
Sequence  * cseq = mapp.GetCurSequence();
Voltage   * ArrPtr = cseq -> Ptr();

ofstream out(filename);

//out<<"Sequence"<<endl;

for(int i = 0; i< mapp.ElectrodeArray::GetSize(); i++)
   {
    if( ArrPtr[i].GetStatus() != 0 )
      {
	 //out<<"POINT"<<endl;
	 out<<(ArrPtr[i]).hId()<<","
	    <<(ArrPtr[i]).vId()<<","
	    <<(ArrPtr[i]).U()<<endl;
      }
   }
}

///////////////////////////////////////////////////////





Mapping mapp;
Sequence seq;
main()
{

Charge q1;
q1.x=0;
q1.y=0;
q1.z=0;
q1.q=1*10e-15;

Charge q2;
q2.x=-0.05;
q2.y=0;
q2.z=0;
q2.q=-1*10e-15;

Charge q3;
q3.x=-0.025;
q3.y=0;
q3.z=0;
q3.q=1;

mapp.AddSequence();
mapp.AddSequence();

mapp.SetCurSequence(1);
mapp.AddCharge(q1);
mapp.AddCharge(q2);
//mapp.RemoveCharge(1);

mapp.SetCurSequence(1);
mapp.CalculateCurSequence();

mapp.Out();

//mapp.RemoveSequence();

//mapp.Out();
//mapp.ElectrodeArray::Out();

//ChargeSet cs;
//cout<<endl<<"adding"<<endl;
//cs.AddCharge(q1);
//cs.AddCharge(q2);
//cs.AddCharge(q3);
//cs.AddCharge(q3);
//cout<<endl<<"size"<<cs.Size();
//cs.Out();
//cs.RemoveCharge(4);
//cs.RemoveCharge(2);
//cs.RemoveCharge(1);
//cs.Out();
//cs.RemoveCharge(1);
//cs.Out();

//seq.AddCharge(q1);
//seq.AddCharge(q1);
//seq.AddCharge(q1);
//seq.Out();

Manager m;
m.SaveElectrodes(mapp,"pizda.scr","ACAD");
//m.SaveCurSequence(mapp,"s.scr","SURFER");
m.SaveCurSequence(mapp,"s.dat","SURFER");
cout<<endl<<"ok";


}