#include "particleBeam.h"
#include "mathematicalConstants.h"
#include "PhysicalConstants.h"
//#include <string>
#include <stdio.h>

#include "environmentVariables.h"


bool particleBeam::setFromParmela(unsigned numeroElement, double referencefrequency)
{
  unsigned int k;
  FILE* filefais;
  string nomfilefais = WORKINGAREA + "parmdesz";
  cout << " nom fichier desz : " << nomfilefais << endl;
  filefais = fopen(nomfilefais.c_str(), "r");

  if ( filefais == (FILE*)0 ) 
    {
      cerr << " particleBeam::setFromParmela() erreur a l'ouverture du fichier" << nomfilefais  << endl;; 
      return false;
    }
  else cout << " particleBeam::setFromParmela() : ouverture du fichier " << nomfilefais << endl; 

  struct particle partic;
  struct particle* partRefPtr=NULL;

  std::vector<struct particle> faisceau;

  //  int numElem=-1;
  partRefPtr=NULL;
  cout << " particleBeam::setFromParmela : numeroElement = " << numeroElement << endl;
  numeroElement--;
  while( partic.readFromParmelaFile(filefais) > 0 )

    {
      // // on ne va conserver que les resultats a la sortie du dernier element
      // if ( partic.ne != numElem )
      // 	{
      // 	  faisceau.clear();
      // 	  partRefPtr=NULL;
      // 	  numElem = partic.ne;
      // 	}
      //      cout << " partic.ne = " << partic.ne << endl;
      if ( partic.ne == numeroElement )
	{
	  //	  cout << " trouve particule " << endl;
	  faisceau.push_back(partic);
	  if ( partic.np == 1 )
	    {
	      if ( fabs(partic.xx) > EPSILON || fabs(partic.yy) > EPSILON || fabs(partic.xxp) > EPSILON  || fabs(partic.yyp) > EPSILON) 
		{
		  printf(" ATTENTION part. reference douteuse  \n");
		  partic.imprim();
		}
	      partRefPtr=&faisceau.back();
	  //	  cout << " indice de la part. ref. " << faisceau.size()-1 << endl;
	  //	      printf("part. reference \n");
	  //	      partRefPtr->imprim();
	    }
	}
    }

  if ( faisceau.size() == 0) return false;

  //  printf("dernier element %d \n", numElem);
  // facteur  c/ 360. pour calculer (c dphi) / (360.freq) 
  // avec freq en Mhz et dphi en degres et résultat en cm: 
  double FACTEUR =  83.3333;  // ameliorer la precision
  double x,xp,y,yp;
  double betagammaz;
  // contrairement a ce qu'indique la notice PARMELA, dans parmdesz, les xp et yp 
  // sont donnes en radians
  // particule de reference
  
  x=partRefPtr->xx;
  xp=partRefPtr->xxp;
  y=partRefPtr->yy;
  yp=partRefPtr->yyp;
  betagammaz = partRefPtr->begamz;
  TRIDVECTOR  posRef(x,y,0.0);
  TRIDVECTOR betagammaRef(xp*betagammaz, yp*betagammaz, betagammaz);

  referenceParticle_ = bareParticle(posRef, betagammaRef);


  // pour l'instant on choisit un centroid nul;
  centroid_ = vector<double>(6,0.0);

  goodPartic_.clear();
  for ( k=0; k < faisceau.size(); k++)
    //  for (int k=0; k < 10; k++)
    {
      x=faisceau.at(k).xx;
      xp=faisceau.at(k).xxp;
      y=faisceau.at(k).yy;
      yp=faisceau.at(k).yyp;
      // dephasage par rapport a la reference  
      double dephas = faisceau.at(k).phi - partRefPtr->phi; // degrés
      double g = faisceau.at(k).wz/ERESTMeV;
      betagammaz = faisceau.at(k).begamz;
      double betaz = betagammaz/(g+1.0);
      double deltaz = FACTEUR * betaz * dephas / referencefrequency;
      x += xp * deltaz;
      y += yp * deltaz;
      TRIDVECTOR  pos(x,y,deltaz);
      TRIDVECTOR betagamma(xp*betagammaz, yp*betagammaz, betagammaz);
      bareParticle bp(pos,betagamma);
      goodPartic_.push_back(bp);
    }
  particleRepresentationOk_ = true;
  //  buildMomentRepresentation();
  return true;
}


void particleBeam::getVariance(double& varx, double& vary, double& varz) const
{

  double x,y,z;
  double xav = 0.;
  double yav = 0.;
  double zav = 0.;
  double xavsq = 0.;
  double yavsq = 0.;
  double zavsq = 0.;


  TRIDVECTOR pos;

  unsigned int k;

  for ( k = 0 ; k < goodPartic_.size(); k++)
    {
      pos = goodPartic_.at(k).getPosition();
      pos.getComponents(x,y,z);
      //      partic_[k].getXYZ(x,y,z);
      xav += x;
      xavsq += x*x;
      yav += y;
      yavsq += y*y;
      zav += z;
      zavsq += z*z;
    }

  double aginv = double (goodPartic_.size());
  aginv = 1.0/aginv;

  varx =  aginv * ( xavsq - xav*xav*aginv );
  vary =  aginv * ( yavsq - yav*yav*aginv );
  varz =  aginv * ( zavsq - zav*zav*aginv );

}


void particleBeam::printAllXYZ() const
{
  cout << " dump du faisceau : " << endl;

  cout <<  goodPartic_.size() << " particules " << endl;
  unsigned int k;
  for ( k = 0 ; k < goodPartic_.size(); k++)
    {
      double xx,yy,zz;
      goodPartic_.at(k).getPosition().getComponents(xx,yy,zz);
      cout << " part. numero " << k << "  x= " << xx << " y= " << yy  << " z= " << zz << endl;
    }
}



void particleBeam::Zrange(double& zmin, double& zmax) const
{
  double z;
  zmin = GRAND;
  zmax = -zmin;

  unsigned int k;
  for ( k = 0 ; k < goodPartic_.size(); k++)
    {
      z = goodPartic_.at(k).getZ();
      if ( z < zmin ) zmin = z;
      else if ( z > zmax) zmax = z;	    
    }
}



string particleBeam::FileOutputFlow() const
{
  ostringstream sortie;
  unsigned int k;
  for ( k = 0 ; k < goodPartic_.size(); k++)
    {
      sortie << goodPartic_.at(k).FileOutputFlow() << endl;
    }
  sortie << endl;
  return sortie.str();
}

bool particleBeam::FileInput( ifstream& ifs)
{
  bool test = true;
  string dum1, dum2;
  double dummy;
  if ( !( ifs >> dum1 >> dum2 >> dummy) ) return false;
 
  bareParticle pp;
  while ( pp.FileInput(ifs) )
    {
      addParticle( pp);
    }
  return test;
}

void particleBeam::buildMomentRepresentation()
{

  unsigned k,j,m;
  double auxj, auxm;

  if ( !particleRepresentationOk_)
    {
      cerr << " particleBeam::buildMomentRepresentation() vecteur de particules invalide" << endl;
    }

  cout << " buildMomentRepresentation " << endl;
  //  printAllXYZ();

  double gref = referenceParticle_.getGamma() - 1.0;
  double P_reference_MeV_sur_c = sqrt( gref*(gref+2) );

  // initialisation des moments
  for ( j = 0; j < 6; j++)
    {
      for (m=0; m <= j; m++) 
	{
	  ( rij_transportMoments_.at(j) ).at(m) = 0.0; // element r_jm
	}
    }

  // accumulation
  for (k=0; k < goodPartic_.size(); k++)
    {
      bareParticle bp = goodPartic_.at(k);
      double gamma = bp.getGamma();
      TRIDVECTOR pos = bp.getPosition();
      TRIDVECTOR begam= bp.getBetaGamma();
      double begamz = begam.getComponent(2);
      double g = gamma -1.0;
      double PMeVsc = sqrt( g*(g+2) );
      double del = 100.0 * ( PMeVsc -  P_reference_MeV_sur_c ) / P_reference_MeV_sur_c ; // en %

      vector<double> part(6);
      part[0] = pos.getComponent(0);
      part[1] = begam.getComponent(0)/begamz;
      part[2] = pos.getComponent(1);
      part[3] = begam.getComponent(1)/begamz;

      part[4] = pos.getComponent(2);

      part[5] = del;

      //      cout << " buildMomentRepresentation part. " << k << " x= " << part[0] << " xp= " << part[1] << " y= " << part[2] << " yp= " << part[3] << endl;

      for ( j = 0; j < 6; j++)
      	{
	  auxj = part.at(j) - centroid_.at(j);
      	  for (m=0; m <= j; m++) 
      	    {
	      auxm = part.at(m) - centroid_.at(m);
	      ( rij_transportMoments_.at(j) ).at(m) += auxj*auxm;
	      //	      cout << " j= " << j << " m= " << m << " rjm= " << ( rij_transportMoments_.at(j) ).at(m) << endl;
      	    }
      	}
    }


  // moyenne
  double facmoy = 1.0/double( goodPartic_.size() );
  for ( j = 0; j < 6; j++)
    {
      ( rij_transportMoments_.at(j) ).at(j) = sqrt(( rij_transportMoments_.at(j) ).at(j) * facmoy );
    }

  for ( j = 0; j < 6; j++)
    {
      auxj =  ( rij_transportMoments_.at(j) ).at(j);
      for (m=0; m < j; m++) 
	{
	  auxm = ( rij_transportMoments_.at(m) ).at(m);
	  ( rij_transportMoments_.at(j) ).at(m) *= facmoy/(auxj * auxm);
	}
    }
    
  // les longueurs sont en cm
  // les angles en radians, on passe en mrad;

  double uniteAngle = 1.0e+3;
  ( rij_transportMoments_.at(1) ).at(1)  *= uniteAngle;

  ( rij_transportMoments_.at(3) ).at(3)  *= uniteAngle;
   
  P0Transport_ = 1.0e-3*ERESTMeV*P_reference_MeV_sur_c;

  // cout << " impression des moments " << endl;
  // for ( j = 0; j < 6; j++)
  //   {
  //     for (m=0; m <= j; m++) 
  // 	{
  // 	      cout  << ( rij_transportMoments_.at(j) ).at(m) << " ";
  // 	}
  //     cout << endl;
  //   }

  momentRepresentationOk_ = true;
}

bool  particleBeam::setFromTransport(ifstream& inp, unsigned nblignes)
  {
    unsigned k;
       unsigned j,m;
    unsigned  nl0;
    //    cout << " particleBeam : lecture resultats transport nblignes= " << nblignes << endl;
      string buf;
       nl0 = nblignes - 7;
       //       cout << " setFromTransport nblignes= " << nblignes << endl;
       for (k=0; k < nl0; k++) 
       	 {
	   getline(inp, buf);
	   // cout  << " ligne " << k+1 << " buf= " << buf << endl;
	 }

       readTransportMoments(inp);

  cout << " impression des moments " << endl;
  for ( j = 0; j < 6; j++)
    {
      for (m=0; m <= j; m++) 
  	{
  	      cout  << ( rij_transportMoments_.at(j) ).at(m) << " ";
  	}
      cout << endl;
    }

     momentRepresentationOk_ = true;
     return true;
  }



void particleBeam::readTransportMoments(ifstream& inp)
{
        unsigned j,m;
      string  bidString;
       double bidon;

  // initialisation des moments
  for ( j = 0; j < 6; j++)
    {
      for (m=0; m <= j; m++) 
	{
	  ( rij_transportMoments_.at(j) ).at(m) = 0.0; // element r_jm
	}
    }

        inp >> bidon >>  bidString >>  bidon >>  ( rij_transportMoments_.at(0) ).at(0) >> bidString;
        inp >> bidon >> ( rij_transportMoments_.at(1) ).at(1) >>  bidString >> ( rij_transportMoments_.at(1) ).at(0);
        inp >> bidon >> ( rij_transportMoments_.at(2) ).at(2) >>  bidString >> ( rij_transportMoments_.at(2) ).at(0)  >> ( rij_transportMoments_.at(2) ).at(1);
        inp >> bidon >> ( rij_transportMoments_.at(3) ).at(3) >>  bidString >> ( rij_transportMoments_.at(3) ).at(0)  >> ( rij_transportMoments_.at(3) ).at(1) >> ( rij_transportMoments_.at(3) ).at(2);

        inp >> bidon >> ( rij_transportMoments_.at(4) ).at(4) >>  bidString >> ( rij_transportMoments_.at(4) ).at(0)  >> ( rij_transportMoments_.at(4) ).at(1) >> ( rij_transportMoments_.at(4) ).at(2) >> ( rij_transportMoments_.at(4) ).at(3);

        inp >> bidon >> ( rij_transportMoments_.at(5) ).at(5) >>  bidString >> ( rij_transportMoments_.at(5) ).at(0)  >> ( rij_transportMoments_.at(5) ).at(1) >> ( rij_transportMoments_.at(5) ).at(2) >> ( rij_transportMoments_.at(5) ).at(3) >> ( rij_transportMoments_.at(5) ).at(4);

}



void particleBeam::donneesDessinEllipseXxp(vector<double>& xcor, vector<double>& ycor)
{
  if ( !momentRepresentationOk_ ) return;

  xcor.clear();
  ycor.clear();

  double xm = ( rij_transportMoments_.at(0) ).at(0);
  double ym = ( rij_transportMoments_.at(1) ).at(1);
  double r  = ( rij_transportMoments_.at(1) ).at(0);

  //  cout << " r= " << r << endl;

  double rac = sqrt(1 - r*r);
  double alpha = -r / rac;
  double beta = xm / ( ym * rac);
  //  double gamma = ym / ( xm * rac );
  double epsil = xm * ym * rac;


  int nbintv = 50;
  double pas = 2.0 * xm / nbintv;
  double fac1 = -1.0 / ( beta * beta);
  double fac2 = epsil/beta;
  double fac3 = -alpha/beta;
  int k;
  double x,y;
  double aux;
  for ( k=0; k < nbintv; k++)
    {
      x = -xm + k*pas;
      aux = fac1 * x * x + fac2;
      //     cout << " aux2= " << aux << endl;
      if ( aux <= 0.0 )
	{
	  aux = 0.0;
	}
      else aux = sqrt(aux);
      
      //        y = fac3*x;
            y = fac3*x + aux;
      xcor.push_back(x);
      ycor.push_back(y);
    }

  for ( k=0; k <= nbintv; k++)
    {
      x = xm - k*pas;
      aux =  fac1 * x * x + fac2;
      if ( aux <= 0.0 ) 
	{
	  aux = 0.0;
	}
      else aux = sqrt(aux);
      //   y = fac3*x;
            y = fac3*x - aux;
      xcor.push_back(x);
      ycor.push_back(y);
    }
}