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

Last change on this file since 342 was 342, checked in by touze, 11 years ago

nvx element snapshot
File size: 9.9 KB
Line 
1#include "softwareParmela.h"
2#include "abstractElement.h"
3#include "parmelaParticle.h"
4#include "mathematicalConstants.h"
5#include "PhysicalConstants.h"
6#include "dataManager.h"
7#include "mixedTools.h"
8
9softwareParmela::softwareParmela() : abstractSoftware()
10{;}
11
12softwareParmela::softwareParmela(string inputFileName, globalParameters* globals, dataManager* dt) : abstractSoftware(inputFileName, globals, dt)
13{;}
14
15bool softwareParmela::createInputFile( particleBeam* beamBefore, unsigned int numeroDeb, unsigned int numeroFin, string workingDir)
16{
17  unsigned int k;
18
19  if ( numeroDeb < 1 || numeroFin > dataManager_->getBeamLineSize() ) {
20    dataManager_->consoleMessage(" softwareParmela::createInputFile : numero of element out of limits " );
21    cerr << " numero of element out of limits " << endl;
22    return false;
23  }
24
25  ofstream outfile;
26  string name = workingDir + inputFileName_;
27  outfile.open(name.c_str(), ios::out);
28  if (!outfile) {
29    dataManager_->consoleMessage(" softwareParmela::createInputFile : error opening output stream " );
30    cerr << " softwareParmela::createInputFile : error opening output stream " << name << endl;
31    return false;
32  }
33
34  abstractElement* elPtr;
35  double initalKineticEnergy = 0.0;
36  elPtr = dataManager_->getElementPointerFromNumero(numeroDeb);
37  bool there_is_rfgun = ( elPtr->getNomdElement().getElementType() == RFgun ); 
38  if ( !there_is_rfgun ) {
39    if ( !beamToParmela(workingDir, beamBefore ) ) return false;
40    initalKineticEnergy = beamBefore->referenceKineticEnergyMeV();
41        // les elements de parmela sont indexes de 1 à max, s'in n'y a pas de rfgun
42    offsetNumElem_ = numeroDeb -1;
43  }
44  else {
45    elPtr->setPhaseStep( globParamPtr_->getIntegrationStep() );
46    initalKineticEnergy = elPtr->getInitialKineticEnergy();
47        // les elements de parmela sont indexes de 0 à max, s'il y a un rfgun
48    offsetNumElem_ = numeroDeb;
49  }
50
51 
52  outfile << "TITLE" << endl;
53  outfile << " titre provisoire " << endl;
54  outfile << "RUN /n0=1 /ip=999 /freq=" << globParamPtr_->getFrequency() << "  /z0=0.0 /W0=" << initalKineticEnergy << "  /itype=1" << endl;
55  outfile << "OUTPUT 0" << endl;
56  unsigned int premier = numeroDeb ;
57  if ( there_is_rfgun ) {
58      outfile << dataManager_->getElementPointerFromNumero(numeroDeb)->parmelaOutputFlow() << endl;
59      premier++;
60  } else {
61    outfile << "INPUT 0 /NP=" << beamBefore->getNbParticles() << endl;
62  }
63
64 
65  for ( k = premier; k <= numeroFin; k++)
66    {
67      elPtr = dataManager_->getElementPointerFromNumero(k);
68      if(elPtr->getNomdElement().getElementType() == snapshot) continue;
69      outfile << elPtr->parmelaOutputFlow() << endl;
70    }
71 
72  outfile << "ZOUT" << endl;
73  outfile << "START /wt=0.0 /dwt=" << globParamPtr_->getIntegrationStep() << "  /nsteps=" << globParamPtr_->getNbSteps() << "  nsc=" << globParamPtr_->getScPeriod() << "  /nout=10" << endl;
74  outfile << "END" << endl;
75  outfile.close();
76  dataManager_->consoleMessage("fichier input termine pour PARMELA");
77  return true;
78}
79
80
81bool softwareParmela::execute(unsigned int numeroDeb,unsigned int numeroFin,string workingDir)
82{
83  bool ExecuteStatus = true;
84 
85  ostringstream sortie;
86  sortie << " EXECUTION DE PARMELA DE l'ELEMENT " << numeroDeb << " A L'ELEMENT " << numeroFin << endl;
87
88  string parmelaJob = workingDir + "parmela";
89  parmelaJob += string("   ");
90  parmelaJob += workingDir;
91 
92  string resultOfRun;
93  bool success = launchJob(parmelaJob,resultOfRun);
94  sortie << resultOfRun << endl;
95  if ( !success ) {
96    sortie << " launching of parmela failed " << endl;
97    ExecuteStatus = false;
98  } else {
99    sortie << " successful launching of parmela " << endl;
100    cout << " execution parmela MARCHE " << endl;
101    string::size_type nn = (resultOfRun).find("normal");
102    if ( nn == string::npos ) 
103      {
104        sortie << " abnormal exit of parmela " << endl;
105        ExecuteStatus = false;
106      }
107  }
108
109  dataManager_->consoleMessage(sortie.str());
110  return ExecuteStatus;
111}
112
113bool softwareParmela::buildBeamAfterElements(unsigned int numeroDeb,unsigned int numeroFin, vector<particleBeam>& beams, string workingDir) 
114{
115  bool result = true;
116  cout << "debug:: debut " << numeroDeb << ", fin " << numeroFin << endl;
117
118  // index du premier element de parmela
119  int id= numeroDeb-offsetNumElem_;
120
121  for(unsigned k = numeroDeb; k <= numeroFin; k++) 
122    {
123      abstractElement* elem = dataManager_->getElementPointerFromNumero(k);
124
125      // si l'element est un snapshot, recuperer la sortie precedente
126      if(elem->getNomdElement().getElementType() == snapshot) {
127        int avantDernier = beams.size() - 1;
128        cout<<"["<<k<<"] : ecrit sur fichier le contenu de beam["<<avantDernier<<"]"<<endl;
129        continue;
130      }
131
132      // sinon c'est un element de parmela d'index id
133      beams.push_back(particleBeam());
134      vector<double> centroid;
135      bareParticle refPart;
136      vector<bareParticle> particles;
137      double freq= globParamPtr_->getFrequency();
138
139      if(!beamFromParmela(workingDir,id,freq,centroid,refPart,particles)) 
140        {
141          if(elem->is_accepted_by_software(nomDeLogiciel::parmela) == warning) {
142            int avantDernier = beams.size() -2;
143            beams.back() = beams.at(avantDernier);
144          } else {
145            dataManager_->consoleMessage("softwareParmela::buildBeamAfterElements : reading parmdesz failed");
146            result = false;
147            break;
148          }
149        } else {
150        beams.back().setWithParticles(centroid,refPart,particles); 
151      }
152
153      // l'element de parmela suivant
154      id++;
155    }
156
157  return result;
158}
159
160bool softwareParmela::beamFromParmela(string workingDir,unsigned numeroParmel, double referencefrequency, vector<double>& centroid, bareParticle& refPart,vector<bareParticle>& particles ) 
161{
162  string nomfilefais = workingDir + "parmdesz";
163  cout << " nom fichier desz : " << nomfilefais << endl;
164  FILE *filefais = fopen(nomfilefais.c_str(), "r");
165
166  if ( filefais == (FILE*)0 ) {
167    dataManager_->consoleMessage(" beamFromParmela() erreur a l'ouverture du fichier 'parmdesz'");
168    cerr << " beamFromParmela() erreur a l'ouverture du fichier" << nomfilefais  << endl; 
169    return false;
170  } else 
171    cout << " beamFromParmela() : ouverture du fichier " << nomfilefais << endl; 
172 
173  parmelaParticle partic;
174  std::vector<parmelaParticle> faisceau;
175  int testNombrePartRef =0;
176  double phaseRef = 0.0;
177
178  while( partic.readFromParmelaFile(filefais) > 0 ) 
179    {
180
181      if ( partic.ne == (int)numeroParmel ) {
182        faisceau.push_back(partic);
183        if ( partic.np == 1 ) {
184          // en principe on est sur la particule de reference
185          if ( fabs(partic.xx) > EPSILON || fabs(partic.yy) > EPSILON || fabs(partic.xxp) > EPSILON  || fabs(partic.yyp) > EPSILON) {
186            printf(" ATTENTION part. reference douteuse  \n");
187            partic.imprim();
188          }
189          phaseRef = partic.phi;
190          TRIDVECTOR  posRef(partic.xx,partic.yy,0.0);
191          TRIDVECTOR betagammaRef(partic.xxp*partic.begamz, partic.yyp*partic.begamz, partic.begamz);
192          refPart = bareParticle(posRef, betagammaRef);
193          testNombrePartRef++;
194          if ( testNombrePartRef != 1 ) {
195            dataManager_->consoleMessage(" beamFromParmela : nombre de part. de ref different de 1 :");
196            cerr << " nombre de part. de ref different de 1 : " << testNombrePartRef << " !! " << endl;
197            return false;
198          }
199        }
200      }
201    } //while
202
203  if ( faisceau.size() == 0) 
204    {
205      string stringNum = mixedTools::intToString( (int)numeroParmel );
206      dataManager_->consoleMessage("beamFromParmela echec lecture  element numero relatif  parmela : " + stringNum);
207      cerr << " beamFromParmela echec lecture  element numero relatif  parmela " << numeroParmel << endl;
208      return false;
209    }
210 
211  // facteur  c/ 360. pour calculer (c dphi) / (360.freq)
212  // avec freq en Mhz et dphi en degres et résultat en cm:
213  double FACTEUR =  83.3333;  // ameliorer la precision
214
215  // pour l'instant on choisit un centroid nul;
216  centroid.clear();
217  centroid = vector<double>(6,0.0);
218
219  particles.clear();
220  particles.resize(faisceau.size(), bareParticle());
221  double x,xp,y,yp;
222  double betagammaz;
223  double betaz;
224  double deltaz;
225  double dephas;
226  double g;
227  TRIDVECTOR  pos;
228  TRIDVECTOR betagamma;
229  // contrairement a ce qu'indique la notice PARMELA, dans parmdesz, les xp et yp
230  // sont donnes en radians
231  for (unsigned k = 0; k < faisceau.size(); k++) {
232    x= faisceau.at(k).xx;
233    xp=faisceau.at(k).xxp;
234    y= faisceau.at(k).yy;
235    yp=faisceau.at(k).yyp;
236    // dephasage par rapport a la reference 
237    dephas = faisceau.at(k).phi - phaseRef; // degrés
238    g = faisceau.at(k).wz/ERESTMeV;
239    betagammaz = faisceau.at(k).begamz;
240    betaz = betagammaz/(g+1.0);
241    deltaz = FACTEUR * betaz * dephas / referencefrequency;
242    x += xp * deltaz;
243    y += yp * deltaz;
244    pos.setComponents(x,y,deltaz);
245    betagamma.setComponents(xp*betagammaz, yp*betagammaz, betagammaz);
246    particles.at(k) = bareParticle(pos,betagamma);
247  }
248 
249  return true;
250}
251
252
253// sauvegarde d'un 'particleBeam' sur un fichier parmela, en guise d'INPUT
254// pour l'instant de nom standard 'parin.input0'
255bool softwareParmela::beamToParmela(string workingDir, particleBeam* beam ) {
256  if ( !beam->particleRepresentationOk() ) {
257    dataManager_->consoleMessage("softwareParmela::beamToParmela : beam not in particles form : not yet programmed");
258    cout << " softwareParmela::beamToParmela : beam not in particles form : not yet programmed " << endl;
259    return false;
260  }
261  ofstream outfile;
262  string name = workingDir + "parin.input0";
263  outfile.open(name.c_str(), ios::out);
264  if (!outfile) {
265    dataManager_->consoleMessage(" softwareParmela::beamToParmela : error opening output stream ");
266    cerr << " softwareParmela::beamToParmela : error opening output stream " << name << endl;
267    return false;
268  }
269
270  const vector<bareParticle>& partic = beam->getParticleVector();
271  unsigned k;
272  double weight = 1.0;
273  double xx,yy,zz;
274  double begamx, begamy, begamz;
275  for ( k=0; k < partic.size(); k++) {
276    partic.at(k).getPosition().getComponents(xx,yy,zz);
277    partic.at(k).getBetaGamma().getComponents(begamx, begamy, begamz);
278    outfile << xx << " " << begamx << " " <<  yy << " " << begamy << " " << zz << " " << begamz  << " " << weight << endl;
279  }
280  outfile.close();
281  return true;
282}
283
284
Note: See TracBrowser for help on using the repository browser.