source: PSPA/Interface_Web/trunk/pspaWT/sources/controler/src/softwareParmela.cc @ 298

#include "softwareParmela.h"
#include "abstractElement.h"
#include "parmelaParticle.h"
#include "mathematicalConstants.h"
#include "PhysicalConstants.h"


softwareParmela::softwareParmela() : abstractSoftware()
{
  ;
}

softwareParmela::softwareParmela(globalParameters* globals, elementsCollection* beamLine) : abstractSoftware(globals, beamLine)
{
  ;
}

bool softwareParmela::createInputFile(particleBeam* beamBefore, unsigned int numeroDeb, unsigned int numeroFin, string workingDir)
{
  unsigned int k;

  if ( numeroDeb < 1 || numeroFin > beamLinePtr_->size() ) {
    cerr << " index of element out of limits " << endl;
    return false;
  }


  unsigned indexDeb = numeroDeb - 1;
  unsigned indexFin = numeroFin - 1;

  ofstream outfile;
  string name = workingDir + "parmin";
  outfile.open(name.c_str(), ios::out);
  if (!outfile) {
    cerr << " error opening output stream " << name << endl;
    return false;
  }

  abstractElement* elPtr;
  double initalKineticEnergy = 0.0;
  elPtr = beamLinePtr_->getElementPointerFromIndex(indexDeb);
  if ( elPtr->getNomdElement().getElementType() != RFgun ) {
    cerr << " dataManager::createInputFileParmela : the first element should be rfgun" << endl;
    return false;
  }
  else {
    elPtr->setPhaseStep( globParamPtr_->getIntegrationStep() );
    initalKineticEnergy = elPtr->getInitialKineticEnergy();
  }


  outfile << "TITLE" << endl;
  outfile << " titre provisoire " << endl;
  outfile << "RUN /n0=1 /ip=999 /freq=" << globParamPtr_->getFrequency() << "  /z0=0.0 /W0=" << initalKineticEnergy << "  /itype=1" << endl;
  outfile << "OUTPUT 0" << endl;
    
  for ( k = indexDeb; k <= indexFin; k++)
    {
      elPtr = beamLinePtr_->getElementPointerFromIndex(k);
      outfile << elPtr->parmelaOutputFlow() << endl;
    }

  outfile << "ZOUT" << endl;
  outfile << "START /wt=0.0 /dwt=" << globParamPtr_->getIntegrationStep() << "  /nsteps=" << globParamPtr_->getNbSteps() << "  nsc=" << globParamPtr_->getScPeriod() << "  /nout=10" << endl;
  outfile << "END" << endl;
  outfile.close();
  return true;
}


bool  softwareParmela::execute(unsigned int numeroDeb,unsigned int numeroFin,string workingDir,string& resul)
{
  ostringstream sortie;
  bool ExecuteStatus = true;
  resul.clear();
  // if ( !createInputFile(NULL,numeroDeb,numeroFin, workingDir) )
  //   {
  //     sortie << " error creating parmela input file "  << endl;
  //     resul = sortie.str();
  //     return false;
  //   }

  sortie << " EXECUTION DE PARMELA DE l'ELEMENT " << numeroDeb << " A L'ELEMENT " << numeroFin << endl;

  char buf[132];
  string parmelaJob = workingDir + "parmela";
  parmelaJob += string("   ");
  parmelaJob += workingDir;
  //  cout << " job parmela= " << parmelaJob << endl;

  string resultOfRun;
  bool success = launchJob(parmelaJob, resultOfRun);
  sortie << resultOfRun << endl;
  if ( !success) {
    sortie << " launching of parmela failed " << endl;
    ExecuteStatus = false;
  }
  else {
    cout << " execution parmela MARCHE " << endl;
    sortie << resultOfRun;
    string::size_type nn = (resultOfRun).find("normal");
    if ( nn == string::npos ) 
      {
        sortie << " abnormal exit of parmela " << endl;
        ExecuteStatus = false;
      }
    // else
    //   {
    //  ExecuteStatus = buildBeamAfterElements(numeroDeb,numeroFin, beamAfterElement, workingDir);
    //  if ( !ExecuteStatus ) {
    //            sortie << " reading parmdesz  failed " << endl;
    //  }

    //   }
  }

  resul =  sortie.str();  
  return ExecuteStatus;
}


bool  softwareParmela::buildBeamAfterElements(unsigned int numeroDeb,unsigned int numeroFin, vector<particleBeam>& beams, string workingDir) {
  bool result = true;
        unsigned k;
        for ( k= numeroDeb; k <= numeroFin; k++)
          {
            beams.push_back(particleBeam());
            vector<double> centroid;
            bareParticle refPart;
            vector<bareParticle> particles;
            if (!beamFromParmela(workingDir,k, globParamPtr_->getFrequency(), centroid, refPart,particles ))
              {
                abstractElement* elem = beamLinePtr_->getElementPointerFromNumero(k);
                if ( elem->is_accepted_by_software(nomDeLogiciel::parmela) == warning) {
                  int avantDernier = beams.size() -2;
                  beams.back() = beams.at(avantDernier);
                } else {
                  // sortie << " reading parmdesz  failed " << endl;
                  result = false;
                  break;
                }
              }
            else {
              beams.back().setWithParticles(centroid, refPart,particles); 
            }
          }
        return result;
}




bool softwareParmela::beamFromParmela(string workingDir,unsigned numeroElement, double referencefrequency, vector<double>& centroid, bareParticle& refPart,vector<bareParticle>& particles ) {
  unsigned  k;
  FILE* filefais;
  string nomfilefais = workingDir + "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; 

  parmelaParticle partic;
  std::vector<parmelaParticle> faisceau;

  cout << " particleBeam::setFromParmela : numeroElement = " << numeroElement << endl;
  unsigned indexElement = numeroElement-1;
 



  int testNombrePartRef =0;
  double phaseRef;

  while( partic.readFromParmelaFile(filefais) > 0 ) {
    if ( partic.ne == (int)indexElement )
      {
        faisceau.push_back(partic);

        if ( partic.np == 1 ) {
          // en principe on est sur la particule de reference
          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();
          }
          phaseRef = partic.phi;
          TRIDVECTOR  posRef(partic.xx,partic.yy,0.0);
          TRIDVECTOR betagammaRef(partic.xxp*partic.begamz, partic.yyp*partic.begamz, partic.begamz);
          refPart = bareParticle(posRef, betagammaRef);
          testNombrePartRef++;
          if ( testNombrePartRef != 1 ) {
            cerr << " TROP DE PART. DE REF : " << testNombrePartRef << " !! " << endl;
            return false;
          }
        }
      }
  }

  if ( faisceau.size() == 0) 
    {
      cerr << " particleBeam::setFromParmela echec lecture  element " << numeroElement << endl;
      return false;
    }
  
  // 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



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

  particles.clear();
  particles.resize(faisceau.size(), bareParticle());
  double x,xp,y,yp;
  double betagammaz;
  double betaz;
  double deltaz;
  double dephas;
  double g;
  TRIDVECTOR  pos;
  TRIDVECTOR betagamma;
  // contrairement a ce qu'indique la notice PARMELA, dans parmdesz, les xp et yp 
  // sont donnes en radians
  for ( k=0; k < faisceau.size(); 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  
    dephas = faisceau.at(k).phi - phaseRef; // degrés
    g = faisceau.at(k).wz/ERESTMeV;
    betagammaz = faisceau.at(k).begamz;
    betaz = betagammaz/(g+1.0);
    deltaz = FACTEUR * betaz * dephas / referencefrequency;
    x += xp * deltaz;
    y += yp * deltaz;
    pos.setComponents(x,y,deltaz);
    betagamma.setComponents(xp*betagammaz, yp*betagammaz, betagammaz);
    particles.at(k) = bareParticle(pos,betagamma);
  }
  return true;
}