source: ZHANGProjects/ICOSIM/CPP/trunk/source/lattice.cc @ 15

Last change on this file since 15 was 15, checked in by zhangj, 10 years ago

(1) Modify the bug to connect to the Fortran files. (2) Fix the bug of the definition type of the random functions used in Fortran file crystal_dan_FINAL....f.
File size: 49.8 KB
Line 
1#include <iostream>
2#include <vector>
3#include <string>
4#include <cmath>
5#include "lattice.h"
6using namespace std;
7
8int Lattice::count = 0;
9int Lattice::turn = -2;
10int Lattice::choicePart = 1;
11int Lattice::eltOutNber = 1561;
12
13
14
15//the following two constructors are just usefull for the initial tests
16
17Lattice::Lattice(Element element1, Element element2, StandardCollimator stdcolli, FlukaCollimator flukacolli, MagneticCollimator magnetcolli, const vector <int>& ips,  int size, int npcle)
18    : size_res(size), npart(npcle), ips(ips)
19{
20    addElement(&element1);
21    addElement(&element2);
22    addElement(&stdcolli);
23    addElement(&flukacolli);
24    addElement(&magnetcolli);
25    addCollimator(new StandardCollimator(stdcolli));
26    addCollimator(new FlukaCollimator(flukacolli));
27    addCollimator(new MagneticCollimator(magnetcolli));
28
29}
30
31
32Lattice::Lattice(Element start, const vector <int>& ips)
33    : size_res(1), npart(1), ips(ips)
34{
35    addElement(&start);
36}
37
38
39
40Lattice::~Lattice()
41{
42
43    for (int i(0); i < resColli.size(); ++i) {
44        delete resColli[i];
45    }
46    ipcoll.push_back(0);
47    ipcoll.clear();
48}
49
50
51
52void Lattice::setsize_res(int s)
53{
54    size_res = s;
55}
56
57
58void Lattice::addCollimator(Collimator* colli)
59{
60
61    resColli.push_back(colli);
62    cocount.push_back(0);
63}
64
65void Lattice::addElement(Element* elt)
66{
67
68    reseau.push_back(elt);
69
70}
71
72
73double Lattice::time(Particle& p, const double& l, const double& betgam, const int& elt)
74{
75
76    double c(2.99792458e8); //speed of light [m/s]
77    double beta(sqrt(betgam * betgam / (betgam * betgam + 1))); //relativistic beta
78    double gamma(betgam / beta); //relativistic gamma
79    double dist;
80    double dist2;
81
82
83    dist = sqrt(l * l + (p.coordonnees[1][0] - p.coordonnees[0][0]) * (p.coordonnees[1][0] - p.coordonnees[0][0]) + (p.coordonnees[1][2] - p.coordonnees[0][2]) * (p.coordonnees[1][2] - p.coordonnees[0][2]));
84
85    if (elt != reseau.size() - 1) {
86
87        dist2 = sqrt(l * l + (reseau[elt + 1]->XC - reseau[elt]->XC) * (reseau[elt + 1]->XC - reseau[elt]->XC) + (reseau[elt + 1]->YC - reseau[elt]->YC) * (reseau[elt + 1]->YC - reseau[elt]->YC));
88    } else {
89
90        dist2 = sqrt(l * l + (reseau[0]->XC - reseau[elt]->XC) * (reseau[0]->XC - reseau[elt]->XC) + (reseau[0]->YC - reseau[elt]->YC) * (reseau[0]->YC - reseau[elt]->YC));
91    }
92
93    p.dt = (dist - dist2) / (beta * c);
94
95    return (dist / (beta * c));
96}
97
98
99
100void Lattice::outCoord(const Particle& p, const int& indic, const string& fileOut)
101{
102
103    ofstream output;
104
105    int col(15);
106
107    if ((turn == -2) && (indic == 1)) {
108        output.open(fileOut.c_str());
109    } else {
110        output.open(fileOut.c_str(), ios::out | ios::app);
111    }
112
113    if (output.fail()) {
114
115        cerr << "Warning: problem openning the file " << fileOut << "!" << endl;
116    } else {
117
118        output << setw(col) << p.coordonnees[1][0] << setw(col) << p.coordonnees[1][1] << setw(col) << p.coordonnees[1][2] << setw(col) << p.coordonnees[1][3] << setw(col) << p.coordonnees[1][4] << setw(col) << p.coordonnees[1][5] <<  endl;
119
120    }
121
122    output.close();
123
124}
125
126
127void Lattice::outPunct(const int& elt, const Particle& p, const Particle& p2, const double& var, string outputpath)
128{
129
130    ofstream outstream;
131    string file;
132    file = outputpath + "/coordinates_punctual.dat";
133
134    int col(40);
135
136    if (elt == eltOutNber) {
137
138        if ((turn == -2) || (turn == -1)) {
139            outstream.open(file.c_str());
140        } else {
141            outstream.open(file.c_str(), ios::out | ios::app);
142        }
143
144        if (outstream.fail()) {
145            cerr << "Warning: problem openning the file " << file << "!" << endl;
146        } else {
147
148            outstream << setw(col) << p.coordonnees[1][0] << setw(col) << p.coordonnees[1][1] << setw(col) << p.coordonnees[1][2] << setw(col) << p.coordonnees[1][3] << setw(col) << p.coordonnees[1][4] << setw(col) << p.coordonnees[1][5] <<  endl;
149        }
150
151        outstream.close();
152    }
153
154}
155
156void Lattice::outElt(const int& elt, const Particle& p, string outputpath, int& indication)
157{
158
159    ofstream tusors;
160    string file;
161    file = outputpath + "/coordinates_elt.dat";
162
163    int col(35);
164
165    if (elt == eltOutNber) {
166
167        if (((turn == -2) || (turn == -1)) && (indication == 1)) {
168            tusors.open(file.c_str());
169            indication = 0;
170        } else {
171            tusors.open(file.c_str(), ios::out | ios::app);
172        }
173
174        if (tusors.fail()) {
175            cerr << "Warning: problem openning the file " << file << "!" << endl;
176        } else {
177
178            tusors << p.coordonnees[1][0] << setw(col) << p.coordonnees[1][1] << setw(col) << p.coordonnees[1][2] << setw(col) << p.coordonnees[1][3] << setw(col) << p.coordonnees[1][4] << setw(col) << p.coordonnees[1][5] <<  endl;
179
180        }
181
182        tusors.close();
183    }
184}
185
186
187
188void Lattice::outrf(const double& x1, const double& x2, string outputpath)
189{
190
191    ofstream print;
192    string fich;
193    fich = outputpath + "/valuestestrf.dat";
194
195    int col(30);
196
197    if ((turn == -1) || (turn == 0)) {
198        print.open(fich.c_str());
199    } else {
200        print.open(fich.c_str(), ios::out | ios::app);
201    }
202
203    if (print.fail()) {
204        cerr << "Warning: problem openning the file " << fich << "!" << endl;
205    } else {
206
207        print << setw(col) << setprecision(15) << x1 << setw(col) << setprecision(15) << x2 << endl;
208    }
209
210    print.close();
211}
212
213void Lattice::read(vector <Particle>& bunch)
214{
215
216    ifstream lecture;
217    string nomfich;
218    nomfich = "../sample/output_data_120GeV_500pt.txt";
219    //nomfich = "test.txt";
220
221    lecture.open(nomfich.c_str());
222
223    if (lecture.fail()) {
224
225        cerr << "Warning: problem with the file " << nomfich << " !" << endl;
226    } else {
227
228        bunch.clear();
229        double var;
230        int id;
231        Particle p;
232        string phrase;
233
234        lecture >> ws;
235
236        for (int k(0); k < 500; ++k) {
237
238            //getline(lecture, phrase);
239
240            lecture >> var;
241            p.coordonnees[0][0] = var;
242            p.coordonnees[1][0] = var;
243            lecture >> var;
244            p.coordonnees[0][1] = var;
245            p.coordonnees[1][1] = var;
246            lecture >> var;
247            p.coordonnees[0][2] = var;
248            p.coordonnees[1][2] = var;
249            lecture >> var;
250            p.coordonnees[0][3] = var;
251            p.coordonnees[1][3] = var;
252            lecture >> var;
253            p.coordonnees[0][4] = var;
254            p.coordonnees[1][4] = var;
255            lecture >> var;
256            p.coordonnees[0][5] = var;
257            p.coordonnees[1][5] = var;
258            //lecture >> id;
259            p.Ap0 = 1;
260            p.Zp0 = 1;
261
262            bunch.push_back(p);
263            bunch[bunch.size() - 1].inabs = 1;
264        }
265
266        lecture.close();
267    }
268}
269
270
271/*
272 *
273 */
274void Lattice::trackensemblelinearnew(vector <Particle>& bunch, vector <Particle>& bunchhit, const int& nrev, const double& blowup2, const int& blowupperiod)
275{
276
277  /*nip: number of primary collimators in the accelerator; niph: number of primary collimators hit by the particles; nrevhitp: number of turns when the particle hit the primary collimators */
278    int nip=0, niph=0, nrevhitp=0; 
279    vector <int> iph;
280    double blowup=0.0; /* sqrt(blowup2), beam blow up strength*/
281
282    //we only take the primary collimators into account here, as well as the first and the last elements of the lattice
283    nip = ip.size();//the number of primary collimators
284    niph = nip + 1;
285
286    //ip --> ip-1
287
288    iph.push_back(0);
289
290    for (int j(0); j < ip.size(); ++j) {
291        iph.push_back(ip[j] - 1);
292    }
293    iph.push_back(reseau.size() - 1);
294    blowup = sqrt(blowup2);
295
296    for (int q(0); q < bunch.size(); ++q) {
297        bunch[q].in = 1;
298    }
299
300  //  cout <<" "<<endl;
301   //   cout <<"********************************"<<endl;
302    cout << "Initialising LINEAR R-Matrix (solutions of Hill's equations using the Floquet theory, based on the Twiss parameters)." << endl;
303
304    for (int i(0); i < niph; ++i) { //making matrices for the Twiss transform
305
306        long double cx, sx, cy, sy;
307
308        cx = cos(2 * M_PI * (reseau[iph[i + 1]]->MUX - reseau[iph[i]]->MUX));
309        sx = sin(2 * M_PI * (reseau[iph[i + 1]]->MUX - reseau[iph[i]]->MUX));
310        R11X.push_back(sqrt(reseau[iph[i + 1]]->BETX / reseau[iph[i]]->BETX) * (cx + reseau[iph[i]]->ALFX * sx));
311        R12X.push_back(sqrt(reseau[iph[i + 1]]->BETX * reseau[iph[i]]->BETX)*sx);
312        R21X.push_back(((reseau[iph[i]]->ALFX - reseau[iph[i + 1]]->ALFX)*cx - (1 + reseau[iph[i]]->ALFX * reseau[iph[i + 1]]->ALFX)*sx) / sqrt(reseau[iph[i + 1]]->BETX * reseau[iph[i]]->BETX));
313        R22X.push_back(sqrt(reseau[iph[i]]->BETX / reseau[iph[i + 1]]->BETX) * (cx - reseau[iph[i + 1]]->ALFX * sx));
314
315
316        cy = cos(2 * M_PI * (reseau[iph[i + 1]]->MUY - reseau[iph[i]]->MUY));
317        sy = sin(2 * M_PI * (reseau[iph[i + 1]]->MUY - reseau[iph[i]]->MUY));
318        R11Y.push_back(sqrt(reseau[iph[i + 1]]->BETY / reseau[iph[i]]->BETY) * (cy + reseau[iph[i]]->ALFY * sy));
319        R12Y.push_back(sqrt(reseau[iph[i + 1]]->BETY * reseau[iph[i]]->BETY)*sy);
320        R21Y.push_back(((reseau[iph[i]]->ALFY - reseau[iph[i + 1]]->ALFY)*cy - (1 + reseau[iph[i]]->ALFY * reseau[iph[i + 1]]->ALFY)*sy) / sqrt(reseau[iph[i + 1]]->BETY * reseau[iph[i]]->BETY));
321        R22Y.push_back(sqrt(reseau[iph[i]]->BETY / reseau[iph[i + 1]]->BETY) * (cy - reseau[iph[i + 1]]->ALFY * sy));
322    }
323
324    int count(0);
325    int total(bunch.size());
326    for (int k(0); k < nrev; ++k) { //loop over the number of revolution
327        vector <Particle> bunchtemp;
328
329        ++turn;
330
331        for (int w(0); w < bunch.size(); ++w) {
332            bunchtemp.push_back(bunch[w]);
333        }
334        for (int i(0); i < niph; ++i) { //loop through the primary collimators
335            for (int p(0); p < bunch.size(); ++p) { //loop through the particles
336                if (bunch[p].in == 1) { //to assure the particle is still remainding in the accelerator
337
338                    long double pdepth=0.0, pdepth2=0.0;
339
340                    bunch[p].coordonnees[1][0] = R11X[i] * bunch[p].coordonnees[0][0] + R12X[i] * bunch[p].coordonnees[0][1] + (reseau[iph[i + 1]]->DX - R11X[i] * reseau[iph[i]]->DX - R12X[i] * reseau[iph[i]]->DPX) * bunch[p].coordonnees[0][4];
341                    bunch[p].coordonnees[1][1] = R21X[i] * bunch[p].coordonnees[0][0] + R22X[i] * bunch[p].coordonnees[0][1] + (reseau[iph[i + 1]]->DPX - R21X[i] * reseau[iph[i]]->DX - R22X[i] * reseau[iph[i]]->DPX) * bunch[p].coordonnees[0][4];
342                    bunch[p].coordonnees[1][2] = R11Y[i] * bunch[p].coordonnees[0][2] + R12Y[i] * bunch[p].coordonnees[0][3] + (reseau[iph[i + 1]]->DY - R11Y[i] * reseau[iph[i]]->DY - R12Y[i] * reseau[iph[i]]->DPY) * bunch[p].coordonnees[0][4];
343                    bunch[p].coordonnees[1][3] = R21Y[i] * bunch[p].coordonnees[0][2] + R22Y[i] * bunch[p].coordonnees[0][3] + (reseau[iph[i + 1]]->DPY - R21Y[i] * reseau[iph[i]]->DY - R22Y[i] * reseau[iph[i]]->DPY) * bunch[p].coordonnees[0][4];
344
345                    if (i < niph - 1) {
346                        long double lcoll=0.0, sa=0.0, ca=0.0;
347
348                        sa = sin(resColli[ipcoll[i]]->tcang);//sinus of the collimator's angle
349                        lcoll = resColli[ipcoll[i]]->L;//length of the collimator
350                       
351
352                        if (sa == 0) {
353                            pdepth = abs(bunch[p].coordonnees[1][0]) - resColli[ipcoll[i]]->hgap;//impact parameter at the beginning of the collimaator
354                            pdepth2 = abs(bunch[p].coordonnees[1][0] + lcoll * bunch[p].coordonnees[1][1]) - resColli[ipcoll[i]]->hgap2; //impact parameter at the end of the collimator
355
356                        } else {
357
358                            long double xl, xsl;
359                            ca = cos(resColli[ipcoll[i]]->tcang);
360                            xl = bunch[p].coordonnees[1][0] * ca + bunch[p].coordonnees[1][2] * sa;
361                            xsl = bunch[p].coordonnees[1][1] * ca + bunch[p].coordonnees[1][3] * sa;
362                            pdepth = abs(xl) - resColli[ipcoll[i]]->hgap;
363                            pdepth2 = abs(xl + lcoll * xsl) - resColli[ipcoll[i]]->hgap2;
364                        }
365
366                        if ((pdepth <= 0) && (pdepth2 <= 0)) { //we only continue with the particles that have not disappeared
367                            bunch[p].coordonnees[0][0] = bunch[p].coordonnees[1][0];
368                            bunch[p].coordonnees[0][1] = bunch[p].coordonnees[1][1];
369                            bunch[p].coordonnees[0][2] = bunch[p].coordonnees[1][2];
370                            bunch[p].coordonnees[0][3] = bunch[p].coordonnees[1][3];
371                        } else {
372                            bunch[p].in = false;
373                            count = count + 1;
374                            bunchhit.push_back(bunchtemp[p]);
375                            bunch[p].nrevhitp = k + 1;
376                            apdepth.push_back(pdepth);
377                            apdepth2.push_back(pdepth2);
378                            cocount[i] = cocount[i] + 1;
379                        }
380                    }//end (if i<niph)
381                    else {//we only continue with the particles that have not disappeared
382                        bunch[p].coordonnees[0][0] = bunch[p].coordonnees[1][0];
383                        bunch[p].coordonnees[0][1] = bunch[p].coordonnees[1][1];
384                        bunch[p].coordonnees[0][2] = bunch[p].coordonnees[1][2];
385                        bunch[p].coordonnees[0][3] = bunch[p].coordonnees[1][3];
386                    }
387                }//end if in
388
389                //the following part can be uncommented to have an output of the x- and  y-coordinates of a particle in the file coordinates.dat
390
391                /*if(bunch[p].getidentification() == choicePart){
392                  outCoord(bunch[p], i+1, "coordinates.dat");
393                  }*/
394
395            }//end loop over the particles
396        }//end loop over i
397
398       // cout << "Revolution: " << k + 1 << ", number of particles left: " << total - count << endl;
399
400        if (total - count == 0) { //we return if all the particles are gone
401            return;
402        }
403
404        vector <Particle> tempinabs;
405
406        //we continue just with particles that are not lost
407        for (int w(0); w < bunch.size(); ++w) {
408            if (bunch[w].in != 0) {
409                tempinabs.push_back(bunch[w]);
410            }
411        }
412
413        bunch.clear();
414
415        for (int w(0); w < tempinabs.size(); ++w) {
416            bunch.push_back(tempinabs[w]);
417        }
418
419        tempinabs.clear();
420
421        for (int p(0); p < bunch.size(); ++p) {
422            if ((k + 1) % blowupperiod == 0) { //blowup/diffusion
423                if (bunch[p].in == 1) {
424                //    cout << "Blow-up!!" << k+1 << endl;
425                    int im=0;
426                    im = reseau.size() - 1;
427                    bunch[p].coordonnees[0][0] = (bunch[p].coordonnees[0][0] - reseau[im]->DX * bunch[p].coordonnees[0][4]) * blowup + reseau[im]->DX * bunch[p].coordonnees[0][4];
428                    bunch[p].coordonnees[0][1] = (bunch[p].coordonnees[0][1] - reseau[im]->DPX * bunch[p].coordonnees[0][4]) * blowup + reseau[im]->DPX * bunch[p].coordonnees[0][4];
429                    bunch[p].coordonnees[0][2] = (bunch[p].coordonnees[0][2] - reseau[im]->DY * bunch[p].coordonnees[0][4]) * blowup + reseau[im]->DY * bunch[p].coordonnees[0][4];
430                    bunch[p].coordonnees[0][3] = (bunch[p].coordonnees[0][3] - reseau[im]->DPY * bunch[p].coordonnees[0][4]) * blowup + reseau[im]->DPY * bunch[p].coordonnees[0][4];
431                }
432            }
433        }
434
435    }//end loop over k
436
437}
438
439
440void Lattice::trackensemblechrom(vector <Particle>& bunch, const int& irev, const int& i0, const int& im, const double& Apr, const double& Zpr, const double& wecolli, const double& betgam, const int& nonlinflag, const int& scaleorbit, double& attr1, const int& idpart, const int& idelt, const int& outcoord, const string& plotflag, vector <vector <double> >& xco, vector <vector <double> >& yco, string outputpath, int RFflag, int& indication)
441{
442
443    if (idpart >= 0) {
444        if (idpart > bunch.size()) {
445            cerr << "Warning, the particle that you want to spy using IDPART does not exist!" << "  idpart = " << idpart << ",  bunch.size = " << bunch.size() <<endl;
446        } else {
447            choicePart = idpart;
448        }
449    }
450
451    if (idelt >= 0) {
452        if (idelt >= reseau.size()) {
453            cerr << "Warning, the element after which you want to spy using IDELT does not exist!" << endl;
454        } else {
455            eltOutNber = idelt;
456        }
457    }
458
459    ++turn;
460
461    if (nhitcolli.size() < resColli.size()) {
462        for (int k(0); k < resColli.size(); ++k) {
463            nhitcolli.push_back(0);
464        }
465    }
466
467    if (plotflag == "Yes") {
468        xco.clear();
469        yco.clear();
470        vector <double> temp1, temp2;
471        for (int k(0); k < bunch.size(); ++k) {
472            temp1.push_back(bunch[k].coordonnees[0][0]);
473            temp2.push_back(bunch[k].coordonnees[0][2]);
474        }
475
476        xco.push_back(temp1);
477        yco.push_back(temp2);
478
479        temp1.clear();
480        temp2.clear();
481    }
482
483    //we set the revolution number (trackensemblechrom is called two times during the first turn)
484    int rev;
485
486    if (irev == 0) {
487        rev = 1;
488    } else {
489        rev = irev;
490    }
491
492    if (irev == 0) {
493        cout << "Initialising NONLINEAR R-matrix." << endl;
494
495        double K=0.0;
496        double cx=0.0, sx=0.0, cy=0.0, sy=0.0;
497
498        R11X.clear();
499        R12X.clear();
500        R21X.clear();
501        R22X.clear();
502        R11Y.clear();
503        R12Y.clear();
504        R21Y.clear();
505        R22Y.clear();
506        turn = -2;
507
508        Lelem.push_back(0);
509        R11X.push_back(0);
510        R12X.push_back(0);
511        R21X.push_back(0);
512        R22X.push_back(0);
513        R11Y.push_back(0);
514        R12Y.push_back(0);
515        R21Y.push_back(0);
516        R22Y.push_back(0);
517        R11XC.push_back(0);
518        R12XC.push_back(0);
519        R21XC.push_back(0);
520        R22XC.push_back(0);
521        R11YC.push_back(0);
522        R12YC.push_back(0);
523        R21YC.push_back(0);
524        R22YC.push_back(0);
525
526        for (int i(1); i < size_res; ++i) {
527
528            Lelem.push_back(reseau[i]->S - reseau[i - 1]->S);
529           
530            if (reseau[i]->K1L > 0) {
531                if (Lelem[i] == 0) {
532                    Lelem[i] = 0.001;
533                }
534
535                K = reseau[i]->K1L / Lelem[i];
536                cx = cos(sqrt(K) * Lelem[i]);
537                sx = sin(sqrt(K) * Lelem[i]);
538                R11XC.push_back(0.5 * sqrt(K)*Lelem[i]*sx);
539                R12XC.push_back(-0.5 * Lelem[i]*cx + 0.5 / sqrt(K)*sx);
540                R21XC.push_back(0.5 * K * Lelem[i]*cx + 0.5 * sqrt(K)*sx);
541                R22XC.push_back(0.5 * sqrt(K)*Lelem[i]*sx);
542                cy = cosh(sqrt(K) * Lelem[i]);
543                sy = sinh(sqrt(K) * Lelem[i]);
544                R11YC.push_back(-0.5 * sqrt(K)*Lelem[i]*sy);
545                R12YC.push_back(-0.5 * Lelem[i]*cy + 0.5 / sqrt(K)*sy);
546                R21YC.push_back(-0.5 * K * Lelem[i]*cy - 0.5 * sqrt(K)*sy);
547                R22YC.push_back(-0.5 * sqrt(K)*Lelem[i]*sy);
548            } else if (reseau[i]->K1L < 0) {
549                if (Lelem[i] == 0) {
550                    Lelem[i] = 0.001;
551                }
552
553                K = -reseau[i]->K1L / Lelem[i];
554                cx = cosh(sqrt(K) * Lelem[i]);
555                sx = sinh(sqrt(K) * Lelem[i]);
556                R11XC.push_back(-0.5 * sqrt(K)*Lelem[i]*sx);
557                R12XC.push_back(-0.5 * Lelem[i]*cx + 0.5 / sqrt(K)*sx);
558                R21XC.push_back(-0.5 * K * Lelem[i]*cx - 0.5 * sqrt(K)*sx);
559                R22XC.push_back(-0.5 * sqrt(K)*Lelem[i]*sx);
560                cy = cos(sqrt(K) * Lelem[i]);
561                sy = sin(sqrt(K) * Lelem[i]);
562                R11YC.push_back(0.5 * sqrt(K)*Lelem[i]*sy);
563                R12YC.push_back(-0.5 * Lelem[i]*cy + 0.5 / sqrt(K)*sy);
564                R21YC.push_back(0.5 * K * Lelem[i]*cy + 0.5 * sqrt(K)*sy);
565                R22YC.push_back(0.5 * sqrt(K)*Lelem[i]*sy);
566            } else {
567                R11XC.push_back(0);
568                R12XC.push_back(0);
569                R21XC.push_back(0);
570                R22XC.push_back(0);
571                R11YC.push_back(0);
572                R12YC.push_back(0);
573                R21YC.push_back(0);
574                R22YC.push_back(0);
575            }
576
577            cx = cos(2 * M_PI * (reseau[i]->MUX - reseau[i - 1]->MUX));
578            sx = sin(2 * M_PI * (reseau[i]->MUX - reseau[i - 1]->MUX));
579            R11X.push_back(sqrt(reseau[i]->BETX / reseau[i - 1]->BETX) * (cx + reseau[i - 1]->ALFX * sx));
580            R12X.push_back(sqrt(reseau[i]->BETX * reseau[i - 1]->BETX)*sx);
581            R21X.push_back(((reseau[i - 1]->ALFX - reseau[i]->ALFX)*cx - (1 + reseau[i - 1]->ALFX * reseau[i]->ALFX)*sx) / sqrt(reseau[i]->BETX * reseau[i - 1]->BETX));
582            R22X.push_back(sqrt(reseau[i - 1]->BETX / reseau[i]->BETX) * (cx - reseau[i]->ALFX * sx));
583
584
585            cy = cos(2 * M_PI * (reseau[i]->MUY - reseau[i - 1]->MUY));
586            sy = sin(2 * M_PI * (reseau[i]->MUY - reseau[i - 1]->MUY));
587            R11Y.push_back(sqrt(reseau[i]->BETY / reseau[i - 1]->BETY) * (cy + reseau[i - 1]->ALFY * sy));
588            R12Y.push_back(sqrt(reseau[i]->BETY * reseau[i - 1]->BETY)*sy);
589            R21Y.push_back(((reseau[i - 1]->ALFY - reseau[i]->ALFY)*cy - (1 + reseau[i - 1]->ALFY * reseau[i]->ALFY)*sy) / sqrt(reseau[i]->BETY * reseau[i - 1]->BETY));
590            R22Y.push_back(sqrt(reseau[i - 1]->BETY / reseau[i]->BETY) * (cy - reseau[i]->ALFY * sy));
591
592        }
593
594    }
595
596
597    // loop througt the elements in the machine
598
599    for (int i(i0 + 1); i <= im; ++i) { //i is the element number along the ring (similar numbering as s)
600
601        /*if(i == 2822){
602          read(bunch);
603          }*/
604
605        for (int p(0); p < bunch.size(); ++p) { //loop througt the particles
606
607            double dpopeff;
608
609            if (bunch[p].inabs == 1) {
610
611                //cout <<"Element " << i << endl;
612
613                int icop(-1);
614
615                for (int k(0); k < ips.size(); ++k) {
616                    if (i == ips[k]) { //element is a collimator
617                        icop = k;//The element is the 'icop'th collimator
618                    }
619                }
620
621                if (icop != -1) { //the element we consider is a collimator
622
623//cout << "icop = "<<icop << endl;
624
625                    if (resColli[icop]->method == "standard") {
626
627                        resColli[icop]->collipass(bunch[p], dpopeff, scaleorbit, R11X[i], R12X[i], R21X[i], R22X[i], R11Y[i], R12Y[i], R21Y[i], R22Y[i], reseau[i - 1]->DX, reseau[i - 1]->DPX, reseau[i - 1]->DY, reseau[i - 1]->DPY, reseau[i]->S - reseau[i - 1]->S, Apr, Zpr, betgam);
628
629                        if (bunch[p].Ap0 != 0) {
630                            bunch[p].coordonnees[1][5] = bunch[p].coordonnees[0][5] + time(bunch[p], reseau[i]->L, betgam, i);
631                        }
632
633                    }
634#if defined(FLUKA)
635                    else if (resColli[icop]->method == "fluka") {
636
637                        vector <Particle> temp;
638
639                        if (p == 0) { //we send all the bunch at the same time to Fluka only one time per element
640
641                            vector <Particle> bunchend;//bunch after the passage through Fluka
642
643                            for (int y(0); y < bunch.size(); ++y) {
644                                temp.push_back(bunch[y]);
645                            }
646
647                            resColli[icop]->collipassfluka(bunch, bunchend, conn, turn, momentum);
648
649
650                            nhitcolli[icop] = bunch.size() - bunchend.size(); //nhitcolli(icop) gives the number particles getting lost in collimator number icop,
651
652                            int creation(0);
653
654                            bunch.clear();
655
656                            for (int g(0); g < bunchend.size(); ++g) {
657                                if (bunchend[g].getidentification() == 0) {
658                                    bunch.push_back(bunchend[g]);
659                                    break;
660                                }
661                            }
662
663                            for (int g(0); g < bunchend.size(); ++g) {
664                                if (bunchend[g].getidentification() == bunch[bunch.size() - 1].getidentification() + 1) {
665                                    bunch.push_back(bunchend[g]);
666                                } else if (bunchend[g].getidentification() == bunch[bunch.size() - 1].getidentification() + 10000) {
667                                    bunch.push_back(bunchend[g]);
668                                    ++creation;
669                                } else if ((bunchend[g].getidentification() == bunch[bunch.size() - 1].getidentification()) && (bunchend[g].getidentification() != 0)) {
670                                    bunch.push_back(bunchend[g]);
671                                    ++creation;
672                                } else if (bunchend[g].getidentification() == bunch[bunch.size() - 1].getidentification() - 10000 + 1) {
673                                    bunch.push_back(bunchend[g]);
674                                }
675                            }
676
677                            nhitcolli[icop] = nhitcolli[icop] + creation;
678
679                            cout << creation << " particles created." << endl;
680
681                            int uu(1);
682                            int rr0, rr1;
683                            long double tps;
684                            bunch[0].coordonnees[0][5] = temp[0].coordonnees[0][5];
685                            rr0 = bunch[0].getidentification();
686                            bunch[0].setidentification(temp[0].getidentification());
687                            for (int g(1); g < bunch.size(); ++g) {
688                                rr1 = bunch[g].getidentification();
689                                if (bunch[g].getidentification() == rr0 + 1) {
690                                    bunch[g].setidentification(temp[uu].getidentification());
691                                    bunch[g].coordonnees[0][5] = temp[uu].coordonnees[0][5];
692                                    rr0 = rr1;
693                                } else {
694                                    tps = temp[uu - 1].coordonnees[0][5];
695                                    while (bunch[g].getidentification() == rr1) {
696                                        bunch[g].coordonnees[0][5] = tps;
697                                        ++g;
698                                    }
699                                    --g;
700                                    --uu;
701                                }
702                                ++uu;
703                            }
704
705                            for (int g(0); g < bunch.size(); ++g) {
706                                bunch[g].coordonnees[1][0] = bunch[g].coordonnees[0][0];
707                                bunch[g].coordonnees[1][1] = bunch[g].coordonnees[0][1];
708                                bunch[g].coordonnees[1][2] = bunch[g].coordonnees[0][2];
709                                bunch[g].coordonnees[1][3] = bunch[g].coordonnees[0][3];
710                                bunch[g].coordonnees[1][4] = bunch[g].coordonnees[0][4];
711                            }
712
713                        }
714
715                        bunch[p].coordonnees[1][5] = bunch[p].coordonnees[0][5] + time(bunch[p], reseau[i]->L, betgam, i);
716
717                        for (int k(0); k < bunch.size(); ++k) {
718                            bunch[k].inabs = 1;
719                        }
720
721                        temp.clear();
722
723                    }
724#endif
725                    else if (resColli[icop]->method == "magnetic") {
726
727                        if (p == 0) {
728                            resColli[icop]->hgap = resColli[icop]->hgap + resColli[icop]->deltaGap;
729                            resColli[icop]->hgap2 = resColli[icop]->hgap2 + resColli[icop]->deltaGap;
730                        }
731
732
733                        resColli[icop]->collipass(bunch[p], dpopeff, scaleorbit, R11X[i], R12X[i], R21X[i], R22X[i], R11Y[i], R12Y[i], R21Y[i], R22Y[i], reseau[i - 1]->DX, reseau[i - 1]->DPX, reseau[i - 1]->DY, reseau[i - 1]->DPY, reseau[i]->S - reseau[i - 1]->S, Apr, Zpr, betgam);
734
735
736                        if (bunch[p].Ap0 != 0) {
737                            bunch[p].coordonnees[1][5] = bunch[p].coordonnees[0][5] + time(bunch[p], reseau[i]->L, betgam, i);
738                        }
739
740                    } else if (resColli[icop]->method == "crystal") {
741
742                        if (p == 0) { //we send all the bunch at the same time through collipassCrystal
743
744                            int lost(bunch.size());
745
746                            resColli[icop]->collipassCrystal(bunch, betgam, turn, outputpath);
747
748                            nhitcolli[icop] = lost - bunch.size(); //nhitcolli(icop) gives the number particles getting lost in collimator number icop,
749
750                        }
751
752                        bunch[p].coordonnees[1][5] = bunch[p].coordonnees[0][5] + time(bunch[p], reseau[i]->L, betgam, i);
753
754                    } else {
755
756                        cerr << "Error: unknown type of collimator, the method is not good defined!!" << endl;
757                    }
758
759                    //we test if the particle is lost in the preceeding collimator
760                    if (bunch[p].Ap0 == 0) {
761                        nhitcolli[icop] = nhitcolli[icop] + 1;
762                        bunch[p].inabs = 0;
763                    }
764
765                } else {//element not collimator
766
767                    if ((reseau[i]->KEYWORD == "RFCAVITY") && (RFflag == 1)) { //attention: voir jusqu ou aller avec le else (est-ce qu on controle l aperture hit?? pour le moment oui...)
768
769                        //Note thate the phase is taken here to be equal to pi.
770
771                        double period(rfharmonic / freqrf);
772                        double omega;
773                        double phase(0);
774                        double phi;
775                        double beta(sqrt(betgam * betgam / (betgam * betgam + 1))); //relativistic beta
776                        long double c(2.99792458e8);//speed of light [m/s]
777                        long double e(1.60218e-19);//elementary charge [C]
778
779
780                        omega = 2 * M_PI * (freqrf / rfharmonic);
781                        phi = phase + omega * bunch[p].coordonnees[0][5];
782
783                        bunch[p].coordonnees[1][0] = bunch[p].coordonnees[0][0];
784                        bunch[p].coordonnees[1][1] = bunch[p].coordonnees[0][1];
785                        bunch[p].coordonnees[1][2] = bunch[p].coordonnees[0][2];
786                        bunch[p].coordonnees[1][3] = bunch[p].coordonnees[0][3];
787
788                        double attr;
789                        attr = bunch[p].Zp0 * sin(phi) * (rfvoltage) / (beta * momentum);
790                        attr1 =  attr;
791                        bunch[p].coordonnees[1][4] = bunch[p].dpoporiginal + attr;//attention vraiment pas sur des parametre (surtout p dans la formule)
792
793                        bunch[p].coordonnees[1][5] = bunch[p].coordonnees[0][5];
794
795                        //uncomment the following line to have output related to the rf-cavity (cf lattice.h)
796                        //outrf(bunch[p].coordonnees[1][5], phi);
797
798
799                    }
800
801
802                    dpopeff = (bunch[p].Ap0 * Zpr) / (bunch[p].Zp0 * Apr) * (1 + bunch[p].coordonnees[0][4]) - 1;
803
804                    if (Lelem[i] == 0) {
805
806                        bunch[p].coordonnees[1][0] = bunch[p].coordonnees[0][0];
807                        bunch[p].coordonnees[1][1] = bunch[p].coordonnees[0][1];
808                        bunch[p].coordonnees[1][2] = bunch[p].coordonnees[0][2];
809                        bunch[p].coordonnees[1][3] = bunch[p].coordonnees[0][3];
810                        bunch[p].coordonnees[1][4] = bunch[p].coordonnees[0][4];
811                        bunch[p].coordonnees[1][5] = bunch[p].coordonnees[0][5];
812
813                    } else if ((reseau[i]->K1L != 0) && (nonlinflag == 1)) {
814
815                        double R11Xh, R12Xh, R21Xh, R22Xh, R11Yh, R12Yh, R21Yh, R22Yh;
816
817
818                        R11Xh = R11X[i] + R11XC[i] * dpopeff;
819                        R12Xh = R12X[i] + R12XC[i] * dpopeff;
820                        R21Xh = R21X[i] + R21XC[i] * dpopeff;
821                        R22Xh = R22X[i] + R22XC[i] * dpopeff;
822                        R11Yh = R11Y[i] + R11YC[i] * dpopeff;
823                        R12Yh = R12Y[i] + R12YC[i] * dpopeff;
824                        R21Yh = R21Y[i] + R21YC[i] * dpopeff;
825                        R22Yh = R22Y[i] + R22YC[i] * dpopeff;
826
827                        bunch[p].coordonnees[1][0] = R11Xh * bunch[p].coordonnees[0][0] + R12Xh * bunch[p].coordonnees[0][1] + (reseau[i]->DX - R11Xh * reseau[i - 1]->DX - R12Xh * reseau[i - 1]->DPX) * dpopeff;
828                        bunch[p].coordonnees[1][1] = R21Xh * bunch[p].coordonnees[0][0] + R22Xh * bunch[p].coordonnees[0][1] + (reseau[i]->DPX - R21Xh * reseau[i - 1]->DX - R22Xh * reseau[i - 1]->DPX) * dpopeff;
829                        bunch[p].coordonnees[1][2] = R11Yh * bunch[p].coordonnees[0][2] + R12Yh * bunch[p].coordonnees[0][3] + (reseau[i]->DY - R11Yh * reseau[i - 1]->DY - R12Yh * reseau[i - 1]->DPY) * dpopeff;
830                        bunch[p].coordonnees[1][3] = R21Yh * bunch[p].coordonnees[0][2] + R22Yh * bunch[p].coordonnees[0][3] + (reseau[i]->DPY - R21Yh * reseau[i - 1]->DY - R22Yh * reseau[i - 1]->DPY) * dpopeff;
831                        bunch[p].coordonnees[1][4] = bunch[p].coordonnees[0][4];
832
833                        bunch[p].coordonnees[1][5] = bunch[p].coordonnees[0][5] + time(bunch[p], reseau[i]->L, betgam, i);
834
835                    } else if ((reseau[i]->K2L != 0) && (nonlinflag == 1)) {
836                        double Lelemha(Lelem[i] / 2);
837                        double dxha(reseau[i - 1]->DX + reseau[i - 1]->DPX * Lelemha);
838                        double dyha(reseau[i - 1]->DY + reseau[i - 1]->DPY * Lelemha);
839
840                        bunch[p].coordonnees[1][0] = bunch[p].coordonnees[0][0] + Lelemha * bunch[p].coordonnees[0][1] + (dxha - reseau[i - 1]->DX - Lelemha * reseau[i - 1]->DPX) * dpopeff;
841                        bunch[p].coordonnees[1][2] = bunch[p].coordonnees[0][2] + Lelemha * bunch[p].coordonnees[0][3] + (dyha - reseau[i - 1]->DY - Lelemha * reseau[i - 1]->DPY) * dpopeff;
842
843                        bunch[p].coordonnees[1][1] = bunch[p].coordonnees[0][1] - 0.5 * reseau[i]->K2L * (bunch[p].coordonnees[1][0] * bunch[p].coordonnees[1][0] - bunch[p].coordonnees[1][2] * bunch[p].coordonnees[1][2]);
844                        bunch[p].coordonnees[1][3] = bunch[p].coordonnees[0][3] + reseau[i]->K2L * (bunch[p].coordonnees[1][0] * bunch[p].coordonnees[1][2]);
845                        bunch[p].coordonnees[1][0] = bunch[p].coordonnees[1][0] + Lelemha * bunch[p].coordonnees[1][1] + (reseau[i]->DX - dxha - Lelemha * reseau[i - 1]->DPX) * dpopeff;
846                        bunch[p].coordonnees[1][2] = bunch[p].coordonnees[1][2] + Lelemha * bunch[p].coordonnees[1][3] + (reseau[i]->DY - dyha - Lelemha * reseau[i - 1]->DPY) * dpopeff;
847                        bunch[p].coordonnees[1][4] = bunch[p].coordonnees[0][4];
848
849                        bunch[p].coordonnees[1][5] = bunch[p].coordonnees[0][5] + time(bunch[p], reseau[i]->L, betgam, i);
850                    } else {
851
852                        bunch[p].coordonnees[1][0] = R11X[i] * bunch[p].coordonnees[0][0] + R12X[i] * bunch[p].coordonnees[0][1] + (reseau[i]->DX - R11X[i] * reseau[i - 1]->DX - R12X[i] * reseau[i - 1]->DPX) * dpopeff;
853                        bunch[p].coordonnees[1][1] = R21X[i] * bunch[p].coordonnees[0][0] + R22X[i] * bunch[p].coordonnees[0][1] + (reseau[i]->DPX - R21X[i] * reseau[i - 1]->DX - R22X[i] * reseau[i - 1]->DPX) * dpopeff;
854                        bunch[p].coordonnees[1][2] = R11Y[i] * bunch[p].coordonnees[0][2] + R12Y[i] * bunch[p].coordonnees[0][3] + (reseau[i]->DY - R11Y[i] * reseau[i - 1]->DY - R12Y[i] * reseau[i - 1]->DPY) * dpopeff;
855                        bunch[p].coordonnees[1][3] = R21Y[i] * bunch[p].coordonnees[0][2] + R22Y[i] * bunch[p].coordonnees[0][3] + (reseau[i]->DPY - R21Y[i] * reseau[i - 1]->DY - R22Y[i] * reseau[i - 1]->DPY) * dpopeff;
856
857                        if ((reseau[i]->KEYWORD != "RFCAVITY") && (RFflag == 1)) {
858                            bunch[p].coordonnees[1][4] = bunch[p].coordonnees[0][4];
859                        }
860
861                        bunch[p].coordonnees[1][5] = bunch[p].coordonnees[0][5] + time(bunch[p], reseau[i]->L, betgam, i);
862
863                    }
864
865                }
866
867                if ((nonlinflag == 1) && ((icop == -1) || (resColli[icop]->method == "standard") || (resColli[icop]->method == "magnetic")) && (bunch[p].inabs != 0)) {
868
869                    int ixcor(-1), iycor (-1);
870
871                    for (int k(0); k < ixcormag.size(); ++k) {
872                        if (ixcormag[k] == i) {
873                            ixcor = k;
874                        }
875                    }
876
877                    if (ixcor != -1) {
878                        bunch[p].coordonnees[1][0] = bunch[p].coordonnees[1][0] + scaleorbit * (reseau[i]->S - reseau[i - 1]->S) / 2 * xcormag[ixcor] * (1 + dpopeff);
879                        bunch[p].coordonnees[1][1] = bunch[p].coordonnees[1][1] + scaleorbit * xcormag[ixcor] * (1 + dpopeff);
880                    }
881
882                    for (int k(0); k < iycormag.size(); ++k) {
883                        if (iycormag[k] == i) {
884                            iycor = k;
885                        }
886                    }
887
888                    if (iycor != -1) {
889                        bunch[p].coordonnees[1][2] = bunch[p].coordonnees[1][2] + scaleorbit * (reseau[i]->S - reseau[i - 1]->S) / 2 * ycormag[iycor] * (1 + dpopeff);
890                        bunch[p].coordonnees[1][3] = bunch[p].coordonnees[1][3] + scaleorbit * ycormag[iycor] * (1 + dpopeff);
891                    }
892
893                }
894
895
896
897                //checking for aperture hits
898                if (bunch[p].inabs == 1) {
899
900                    if ((icop == -1) || (resColli[icop]->method == "standard") || (resColli[icop]->method == "magnetic") || (resColli[icop]->method == "crystal")) {
901
902                        bool inside(false);
903
904                        if (reseau[i]->APERTYPE == "RECTANGLE") {
905                            if ((abs(bunch[p].coordonnees[1][0]) < reseau[i]->aperx) && (abs(bunch[p].coordonnees[1][2]) < reseau[i]->apery) && (abs(bunch[p].coordonnees[0][0]) < reseau[i]->aperx) && (abs(bunch[p].coordonnees[0][2]) < reseau[i]->apery)) {
906                                inside = true;
907                            }
908                        } else if ((reseau[i]->APERTYPE == "ELLIPSE") || (reseau[i]->APERTYPE == "CIRCLE") || (reseau[i]->APERTYPE == "RECTELLIPSE")) {
909                            if ((((bunch[p].coordonnees[1][0] / reseau[i]->aperx) * (bunch[p].coordonnees[1][0] / reseau[i]->aperx) + (bunch[p].coordonnees[1][2] / reseau[i]->apery) * (bunch[p].coordonnees[1][2] / reseau[i]->apery)) < 1) && (((bunch[p].coordonnees[0][0] / reseau[i]->aperx) * (bunch[p].coordonnees[0][0] / reseau[i]->aperx) + (bunch[p].coordonnees[0][2] / reseau[i]->apery) * (bunch[p].coordonnees[0][2] / reseau[i]->apery)) < 1)) {
910                                inside = true;
911
912                            }
913                        } else {
914                            cerr << "Error for the " << i << "th element: unknown aperture type!" << endl;
915                        }
916
917                        double L;
918
919                        L = reseau[i]->S - reseau[i - 1]->S; //distance to next element in the accelerator
920
921                        double Ap0h, Zp0h, xh0, xh, xhs, yh0, yh, yhs;
922                        int nh(0);//number of lost particles
923
924                        if (inside == false) {
925                            cout <<"The particle is lost at element " << reseau[i]->NAME << endl;
926                            xh0 = bunch[p].coordonnees[0][0];
927                            bunch[p].inabs = 0;
928
929                            if (L == 0) {
930                                nh = nh + 1;
931                                hits.push_back(reseau[i]->S);//saving the value of s where the particles get lost
932                                Aphit.push_back(bunch[p].Ap0);//saving the mass of the lost particles
933                                Zphit.push_back(bunch[p].Zp0);//saving the charge of the lost particles
934                            } else {
935
936                                xh = bunch[p].coordonnees[1][0];
937                                xhs = (xh - xh0) / L;
938
939                                yh0 = bunch[p].coordonnees[0][2];
940                                yh = bunch[p].coordonnees[1][2];
941                                yhs = (yh - yh0) / L;
942
943                                Aphit.push_back(bunch[p].Ap0);
944                                Zphit.push_back(bunch[p].Zp0);
945
946                                double slostx, slosty, splus, sminus;
947
948                                if (reseau[i]->APERTYPE == "RECTANGLE") {
949
950                                    //find hit position in x
951
952                                    if ((abs(xh0) < reseau[i]->aperx) && (abs(xh) < reseau[i]->aperx)) { //both points inside aperture - particle lost in y instead
953                                        slostx = reseau[i]->S;
954
955                                    } else if (abs(xh0) > reseau[i]->aperx) { //first point outside - particle lost already at the entrance of the element
956                                        slostx = reseau[i - 1]->S;
957
958                                    } else {//first point inside, second point outside - do a linear interpolation
959
960                                        splus = (reseau[i]->aperx - xh0) / xhs;
961                                        sminus = (-reseau[i]->aperx - xh0) / xhs;
962
963                                        if (splus > sminus) { //choose the correct point where the straight line hits the aperture, at +-a. The highest s-value is correct
964                                            slostx = splus + reseau[i - 1]->S;
965
966                                        } else {
967                                            slostx = sminus + reseau[i - 1]->S;
968
969                                        }
970                                    }
971
972                                    //find hit position in y
973
974                                    if ((abs(yh0) < reseau[i]->apery) && (abs(yh) < reseau[i]->apery)) { //both points inside aperture - particle lost in x instead
975                                        slosty = reseau[i]->S;
976                                    } else if (abs(yh0) > reseau[i]->apery) { //first point outside - particle lost already at the entrance of the element
977                                        slosty = reseau[i - 1]->S;
978                                    } else {//first point inside, second point outside - do a linear interpolation
979
980                                        splus = (reseau[i]->apery - yh0) / yhs;
981                                        sminus = (-reseau[i]->apery - yh0) / yhs;
982
983                                        if (splus > sminus) { //choose the right point where the straight line hits the aperture, at +-b. The highest s-value is correct
984                                            slosty = splus + reseau[i - 1]->S;
985                                        } else {
986                                            slosty = sminus + reseau[i - 1]->S;
987                                        }
988                                    }
989
990                                    //choose the point where the particle hits first, x or y
991
992                                    if (slostx < slosty) {
993                                        hits.push_back(slostx);
994                                    } else {
995                                        hits.push_back(slosty);
996                                    }
997                                } else if (reseau[i]->APERTYPE == "ELLIPSE") {
998
999                                    double sqrarg;
1000                                    sqrarg = reseau[i]->apery * reseau[i]->apery * xhs * xhs - xhs * xhs * yh0 * yh0 + 2 * xh0 * xhs * yh0 * yhs + reseau[i]->aperx * reseau[i]->aperx * yhs * yhs - xh0 * xh0 * yhs * yhs;
1001                                    if (sqrarg < 0) {
1002                                        cout << "Warning, sqrarg < 0" << endl;
1003                                        hits.push_back(reseau[i - 1]->S);
1004                                    } else {
1005                                        double stry;
1006
1007                                        stry = (-reseau[i]->apery * reseau[i]->apery * xh0 * xhs - reseau[i]->aperx * reseau[i]->aperx * yh0 * yhs - reseau[i]->aperx * reseau[i]->apery * sqrt(sqrarg)) / (reseau[i]->apery * reseau[i]->apery * xhs * xhs + reseau[i]->aperx * reseau[i]->aperx * yhs * yhs);
1008
1009                                        if ((stry > 0) && (stry < L)) {
1010                                            hits.push_back(reseau[i - 1]->S + stry);
1011                                        } else {
1012
1013                                            double stry2;
1014
1015                                            stry2 = (-reseau[i]->apery * reseau[i]->apery * xh0 * xhs - reseau[i]->aperx * reseau[i]->aperx * yh0 * yhs + reseau[i]->aperx * reseau[i]->aperx * sqrt(sqrarg)) / (reseau[i]->apery * reseau[i]->apery * xhs * xhs + reseau[i]->aperx * reseau[i]->aperx * yhs * yhs);
1016
1017                                            if ((stry2 > 0) && (stry2 < L)) {
1018                                                hits.push_back(reseau[i - 1]->S + stry2);
1019                                            } else {
1020                                                hits.push_back(reseau[i - 1]->S);
1021                                            }
1022                                        }
1023                                    }
1024                                }//end if ELLIPSE
1025                                else if (reseau[i]->APERTYPE == "RECTELLIPSE") {
1026
1027                                    double sqrarg;
1028                                    sqrarg = reseau[i]->apery * reseau[i]->apery * xhs * xhs - xhs * xhs * yh0 * yh0 + 2 * xh0 * xhs * yh0 * yhs + reseau[i]->aperx * reseau[i]->aperx * yhs * yhs - xh0 * xh0 * yhs * yhs;
1029                                    if (sqrarg < 0) {
1030                                        cout << "Warning, sqrarg < 0" << endl;
1031                                        hits.push_back(reseau[i - 1]->S);
1032                                    } else {
1033                                        double stry;
1034                                        stry = (-reseau[i]->apery * reseau[i]->apery * xh0 * xhs - reseau[i]->aperx * reseau[i]->aperx * yh0 * yhs - reseau[i]->aperx * reseau[i]->apery * sqrt(sqrarg)) / (reseau[i]->apery * reseau[i]->apery * xhs * xhs + reseau[i]->aperx * reseau[i]->aperx * yhs * yhs);
1035
1036                                        if ((stry > 0) && (stry < L)) {
1037                                            hits.push_back(reseau[i - 1]->S + stry);
1038                                        } else {
1039                                            double stry2;
1040
1041                                            stry2 = (-reseau[i]->apery * reseau[i]->apery * xh0 * xhs - reseau[i]->aperx * reseau[i]->aperx * yh0 * yhs + reseau[i]->aperx * reseau[i]->apery * sqrt(sqrarg)) / (reseau[i]->apery * reseau[i]->apery * xhs * xhs + reseau[i]->aperx * reseau[i]->aperx * yhs * yhs);
1042
1043                                            if ((stry2 > 0) && (stry2 < L)) {
1044                                                hits.push_back(reseau[i - 1]->S + stry2);
1045                                            } else {
1046                                                hits.push_back(reseau[i - 1]->S);
1047                                            }
1048                                        }
1049                                    }
1050                                }//end RECTELLIPSE
1051                            }//end else after if L == 0
1052
1053                        }//end si particle lost
1054                        else {
1055                            //cout << "No losses, all the particles stay in the experience." << endl;
1056                        }
1057                    }//end if icop==vide ||...
1058
1059                }
1060            }
1061
1062            //Different ways to print the coordinates of the particle through the parameters IDPART and OUTCOORD (see manual)
1063
1064            if ((outcoord != 0) && (bunch[p].inabs != 0)) {
1065
1066                if (outcoord == 1) {
1067
1068                    cout << "Particle " <<  p + 1 << " after the passage through element: " << i << " " << reseau[i]->NAME << endl;
1069                    bunch[p].afficheCoordonnees();
1070                    cout << "Turn number " << turn << endl;
1071                    cout << "Massnumber: " << bunch[p].Ap0 << " Chargestate: " << bunch[p].Zp0 << endl;
1072                    cout << "Particle ID: " << bunch[p].getidentification() << endl;
1073                } else if (outcoord == 3) {
1074
1075                    if (bunch[p].getidentification() == choicePart) {
1076                        outCoord(bunch[p], i, outputpath + "/coordinates.dat");
1077                    }
1078                } else if (outcoord == 4) {
1079
1080                    if (bunch[p].getidentification() == choicePart) {
1081                        outPunct(i, bunch[p], bunch[0], attr1, outputpath);
1082                    }
1083                }
1084
1085                else if (outcoord == 2) {
1086                    outElt(i, bunch[p], outputpath, indication);
1087                }
1088            }
1089
1090
1091        }//end loop over particles
1092
1093
1094        //we prepare the particles for the next element
1095        for (int b(0); b < bunch.size(); ++b) {
1096            if (bunch[b].inabs != 0) {
1097                for (int r(0); r < 6; ++r) {
1098                    bunch[b].coordonnees[0][r] = bunch[b].coordonnees[1][r];
1099                }
1100            }
1101        }
1102
1103        if (plotflag == "Yes") {
1104
1105            vector <double> temp1, temp2;
1106
1107            for (int h(0); h < bunch.size(); ++h) {
1108                if (bunch[h].inabs != 0) {
1109                    temp1.push_back(bunch[h].coordonnees[0][0]);
1110                    temp2.push_back(bunch[h].coordonnees[0][2]);
1111                } else {
1112                    temp1.push_back(100);
1113                    temp2.push_back(100);
1114                }
1115            }
1116
1117            xco.push_back(temp1);
1118            yco.push_back(temp2);
1119
1120            temp1.clear();
1121            temp2.clear();
1122        }
1123
1124
1125        //if there are no more particles, we return
1126        if (bunch.size() == 0) {
1127            return;
1128        }
1129
1130    }//end loop elements
1131
1132    vector <Particle> tempinabs;
1133
1134    //we contine just with particles that are not lost
1135    for (int w(0); w < bunch.size(); ++w) {
1136        if (bunch[w].inabs != 0) {
1137            tempinabs.push_back(bunch[w]);
1138        }
1139    }
1140
1141    bunch.clear();
1142
1143    for (int w(0); w < tempinabs.size(); ++w) {
1144        bunch.push_back(tempinabs[w]);
1145    }
1146
1147    tempinabs.clear();
1148
1149}
Note: See TracBrowser for help on using the repository browser.