#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <vector>
#include <cmath>


static double EPSILON=1.e-10;

// facteur  c/ 360. pour calculer (c dphi) / (360.freq) 
// avec freq en Mhz et dphi en degres et résultat en cm: 

static double FACTEUR =  83.3333;  // ameliorer la precision

// contrairement a ce qu'indique la notice, dans parmdesz, les xp et yp 
// sont donnes en radians
static double uniteAngle = 1.0e+3; // passage de radians en milliradians
static double uniteLong = 1.0; // on prend l'unite de long. TRANSPORT : cm


static double EMeV = 0.511;
struct particle
{
  float xx, xxp, begamx,yy,yyp,begamy,z, begamz,phi,wz;
  float phi0, ksi1,ksi2,ksi3;
  int ne,np,ngood,npart;

  int readFromFile(FILE* fp)
  {
    return fscanf(fp, " %e %e %e %e %e %e %e %e %e %e %d %d %d %d %e %e %e %e \n", &xx, &xxp, &begamx,&yy,&yyp,&begamy,&z, &begamz,&phi,&wz,&ne,&np,&ngood,&npart,&phi0, &ksi1,&ksi2,&ksi3);
  }
  void imprim()
  {
    printf( " %e %e %e %e %e %e %e %e %e %e %d %d %d, %d %e %e %e %e \n", xx, xxp, begamx,yy,yyp,begamy,z, begamz,phi,wz,ne,np,ngood,npart,phi0, ksi1,ksi2,ksi3);
  }
};


double readFreqFromParmout(FILE* fp)
{
  float freq=0.0;
  char chaine[80];
  while ( fscanf(fp, " %s ", chaine) > 0 ) 
    {
      if ( strcmp(chaine, "frequency") == 0 ) 
	{
	  printf(" chaone lue :  %s  \n", chaine);
	  fscanf(fp, " %s ", chaine); 
	  fscanf(fp, " %e ", &freq); 
	  break;
	}
    }
  return freq;
};


int main(int argc,char* argv[]) 
{

  double freqRef;
  FILE* filefais;
  FILE* fileParmout;
  FILE* fileFaisTrans;
  std::vector<struct particle> faisceau;
  int k;

  if (argc!=2 ) {
    printf(" usage : recupFaisceau <nom du probleme parmela> \n");
    exit(0);
  }

  printf(" ATTENTION : PROGRAMME recupFaisceauParm NON VERIFIE PRECISEMENT \n"); 

 
  //    printf(" frequence de reference %e Mhz \n", freqRef); 
  char* nomParm = argv[1];
  char nomParmout[132];
  strcpy(nomParmout, nomParm);
  strcat(nomParmout, ".out");

  printf(" probleme parmela :  %s  \n", argv[1]); 
  //
  fileParmout = fopen(nomParmout, "r");
  if ( fileParmout == (FILE*)0 ) 
    {
      printf(" erreur a l'ouverture du fichier %s \n", nomParmout); 
      exit(0);
    }
  else printf(" ouverture du fichier %s \n",  nomParmout ); 
  freqRef= readFreqFromParmout(fileParmout);
  fclose(fileParmout);
  //
  printf(" frequence relue :  %e  \n", freqRef);
  char nomfilefais[132]; 
  strcpy(nomfilefais, nomParm);
  strcat(nomfilefais, ".desz");

  filefais = fopen(nomfilefais, "r");

  if ( filefais == (FILE*)0 ) 
    {
      printf(" erreur a l'ouverture du fichier %s \n",nomfilefais ); 
      exit(0);
    }
  else printf(" ouverture du fichier %s \n", nomfilefais); 

  fileFaisTrans = fopen("faisceauTransp.dat", "w");
  if ( fileFaisTrans == (FILE*)0 ) 
    {
      printf(" erreur a l'ouverture du fichier faisceauTransp.dat \n"); 
      exit(0);
    }

  struct particle partic;
  struct particle* partRefPtr=NULL;
  int numElem=-1;
  while( partic.readFromFile(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;
	}
      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();
	  //	      printf("part. reference \n");
	  //	      partRefPtr->imprim();
	}
    }
  //    for (int k=0; k < faisceau.size(); k++) faisceau[k].imprim(); 
  printf("dernier element %d \n", numElem);
  // pour preserver l'avenir, on definit un centroid (cf. TRANSPORT)
  double xg=0.;
  double xpg=0.;
  double yg=0.;
  double ypg=0.;
  double lg=0.;
  double delg=0.;

  // sigmas TRANSPORT
  double sig11=0.0;
  double sig22=0.0;
  double sig33=0.0;
  double sig44=0.0;
  double sig55=0.0;
  double sig66=0.0;

  // correlations
  double r21=0.0;
  double r31=0.0;
  double r32=0.0;
  double r41=0.0;
  double r42=0.0;
  double r43=0.0;
  double r51=0.0;
  double r52=0.0;
  double r53=0.0;
  double r54=0.0;
  double r61=0.0;
  double r62=0.0;
  double r63=0.0;
  double r64=0.0;
  double r65=0.0;


  double x,xp,y,yp,ll,del;
  // unites : TRANSPORT : cm, mrad, %
  //          PARMELA : phase en degrés


  double gref = partRefPtr->wz/EMeV;
  double PrefMeVsc = sqrt( gref*(gref+2) );


  for (int k=0; k < faisceau.size(); k++)
    //  for (int k=0; k < 10; k++)
    {
      x=faisceau[k].xx;
      xp=faisceau[k].xxp;
      y=faisceau[k].yy;
      yp=faisceau[k].yyp;
      // dephasage par rapport a la reference  
      double dephas = faisceau[k].phi - partRefPtr->phi; // degrés
      double g = faisceau[k].wz/EMeV;
      double betaz = faisceau[k].begamz/(g+1.0);
      // printf( " betaz= %e \n",  betaz);
      // si j'ai bien compris, une phase positive indique une particule en avance dans l'espace
      double deltaz = FACTEUR * betaz * dephas / freqRef;
      //      printf(" correction dz %e \n",  deltaz);
      ll  = deltaz;
      //      printf("ancien x  %e correction %e \n", x, xp * deltaz);
      x += xp * deltaz;
      faisceau[k].xx = x;
      //      printf("ancien y  %e correction %e \n", y, yp * deltaz);
      y += yp * deltaz;
      faisceau[k].yy =y;
      double PMeVsc = sqrt( g*(g+2) );
      del = 100.0 * ( PMeVsc -  PrefMeVsc ) /  PrefMeVsc ; // en %
      //      printf("defaut energie  %e \n", del);
      sig11 += ( x - xg ) * ( x - xg );
      sig22 += ( xp - xpg ) * ( xp - xpg );
      sig33 += ( y - yg ) * ( y - yg );
      sig44 += ( yp - ypg ) * ( yp - ypg );
      sig55 += uniteLong*uniteLong*( ll - lg ) * ( ll - lg );
      sig66 += ( del - delg ) * ( del - delg );
      //


      r21 += ( xp - xpg ) * ( x - xg );
      r31 += ( y - yg ) * ( x - xg );
      r32 += ( y - yg ) * ( xp - xpg );
      r41 += ( yp - ypg ) * ( x - xg );
      r42 += ( yp - ypg ) * ( xp - xpg );
      r43 += ( yp - ypg ) * ( y - yg );
      r51 += ( ll - lg ) * ( x - xg );
      r52 += ( ll - lg ) *  ( xp - xpg );
      r53 += ( ll - lg ) * ( y - yg );
      r54 += ( ll - lg ) * ( yp - ypg );
      r61 += ( del - delg ) * ( x - xg );
      r62 += ( del - delg ) * ( xp - xpg );
      r63 += ( del - delg ) * ( y - yg );
      r64 += ( del - delg ) * ( yp - ypg );
      r65 += ( del - delg ) * ( ll - lg );

    }
  double facmoy = 1.0/double( faisceau.size() );
  sig11 = sqrt(sig11*facmoy);
  sig22 = sqrt(sig22*facmoy);
  sig33 = sqrt(sig33*facmoy);
  sig44 = sqrt(sig44*facmoy);
  sig55 = sqrt(sig55*facmoy);
  sig66 = sqrt(sig66*facmoy);




  //

  r21= r21*facmoy/(sig22*sig11);
  r31= r31*facmoy/(sig33*sig11);
  r32= r32*facmoy/(sig33*sig22);
  r41= r41*facmoy/(sig44*sig11);
  r42= r42*facmoy/(sig44*sig22);
  r43= r43*facmoy/(sig44*sig33);
  r51= r51*facmoy/(sig55*sig11);
  r52= r52*facmoy/(sig55*sig22);
  r53= r53*facmoy/(sig55*sig33);
  r54= r54*facmoy/(sig55*sig44);
  r61= r61*facmoy/(sig66*sig11);
  r62= r62*facmoy/(sig66*sig22);
  r63= r63*facmoy/(sig66*sig33);
  r64= r64*facmoy/(sig66*sig44);
  r65= r65*facmoy/(sig66*sig55);

  // passage aux unites de sortie

  sig11 *= uniteLong;
  sig22 *= uniteAngle;
  sig33 *= uniteLong;
  sig44 *= uniteAngle;
  sig55 *= uniteLong;

  int dummy=0;
  fprintf(fileFaisTrans, " titre du  probleme :  %s  \n", nomParm);
  fprintf(fileFaisTrans, "  %d  \n", dummy);
  fprintf(fileFaisTrans, "  UTRANS  \n");


  fprintf(fileFaisTrans, " BEAM, X=%e, XP=%e, Y=%e, YP=%e, & \n", sig11, sig22, sig33, sig44 );
  fprintf(fileFaisTrans, "  DL=%e, DEL=%e, P0=%e  ; \n", sig55, sig66, 1.0e-3*EMeV*PrefMeVsc );


  fprintf(fileFaisTrans, " CORR, C21=%e, C31=%e, C32=%e, & \n", r21, r31, r32 );
  fprintf(fileFaisTrans, " C41=%e, C42=%e, C43=%e, C51=%e, & \n", r41, r42, r43, r51 );
  fprintf(fileFaisTrans, " C52=%e, C53=%e, C54=%e, C61=%e, & \n", r52, r53, r54, r61 );
  fprintf(fileFaisTrans, " C62=%e, C63=%e, C64=%e, C65=%e ; \n", r62, r63, r64, r65 );


  printf(" le faisceau TRANSPORT est cree sur le fichier  :  faisceauTransp.dat  \n"); 



  fclose(filefais);
  fclose(fileFaisTrans);
}