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

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

Initial import
File size: 49.2 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
272void Lattice::trackensemblelinearnew(vector <Particle>& bunch, vector <Particle>& bunchhit, const int& nrev, const double& blowup2, const int& blowupperiod)
273{
274
275    int nip, niph, nrevhitp;
276    vector <int> iph;
277    double blowup;
278
279    //we only take the primary collimators into account here, as well as the first and the last elements of the lattice
280    nip = ip.size();//the number of primary collimators
281    niph = nip + 1;
282
283    //ip --> ip-1
284
285    iph.push_back(0);
286
287    for (int j(0); j < ip.size(); ++j) {
288        iph.push_back(ip[j] - 1);
289    }
290    iph.push_back(reseau.size() - 1);
291    blowup = sqrt(blowup2);
292
293    for (int q(0); q < bunch.size(); ++q) {
294        bunch[q].in = 1;
295    }
296
297    cout << "Initialising R-Matrix" << endl;
298
299    for (int i(0); i < niph; ++i) { //making matrices for the Twiss transform
300
301        long double cx, sx, cy, sy;
302
303        cx = cos(2 * M_PI * (reseau[iph[i + 1]]->MUX - reseau[iph[i]]->MUX));
304        sx = sin(2 * M_PI * (reseau[iph[i + 1]]->MUX - reseau[iph[i]]->MUX));
305        R11X.push_back(sqrt(reseau[iph[i + 1]]->BETX / reseau[iph[i]]->BETX) * (cx + reseau[iph[i]]->ALFX * sx));
306        R12X.push_back(sqrt(reseau[iph[i + 1]]->BETX * reseau[iph[i]]->BETX)*sx);
307        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));
308        R22X.push_back(sqrt(reseau[iph[i]]->BETX / reseau[iph[i + 1]]->BETX) * (cx - reseau[iph[i + 1]]->ALFX * sx));
309
310
311        cy = cos(2 * M_PI * (reseau[iph[i + 1]]->MUY - reseau[iph[i]]->MUY));
312        sy = sin(2 * M_PI * (reseau[iph[i + 1]]->MUY - reseau[iph[i]]->MUY));
313        R11Y.push_back(sqrt(reseau[iph[i + 1]]->BETY / reseau[iph[i]]->BETY) * (cy + reseau[iph[i]]->ALFY * sy));
314        R12Y.push_back(sqrt(reseau[iph[i + 1]]->BETY * reseau[iph[i]]->BETY)*sy);
315        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));
316        R22Y.push_back(sqrt(reseau[iph[i]]->BETY / reseau[iph[i + 1]]->BETY) * (cy - reseau[iph[i + 1]]->ALFY * sy));
317    }
318
319    int count(0);
320    int total(bunch.size());
321    for (int k(0); k < nrev; ++k) { //loop over the number of revolution
322        vector <Particle> bunchtemp;
323
324        ++turn;
325
326        for (int w(0); w < bunch.size(); ++w) {
327            bunchtemp.push_back(bunch[w]);
328        }
329        for (int i(0); i < niph; ++i) { //loop through the primary collimators
330            for (int p(0); p < bunch.size(); ++p) { //loop through the particles
331                if (bunch[p].in == 1) { //to assure the particle is still remainding in the accelerator
332
333                    long double pdepth, pdepth2;
334
335                    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];
336                    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];
337                    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];
338                    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];
339
340                    if (i < niph - 1) {
341                        long double lcoll, sa, ca;
342
343                        sa = sin(resColli[ipcoll[i]]->tcang);//sinus of the collimator's angle
344                        lcoll = resColli[ipcoll[i]]->L;//length of the collimator
345
346                        if (sa == 0) {
347                            pdepth = abs(bunch[p].coordonnees[1][0]) - resColli[ipcoll[i]]->hgap;//impact parameter at the beginning of the collimaator
348                            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
349
350                        } else {
351
352                            long double xl, xsl;
353                            ca = cos(resColli[ipcoll[i]]->tcang);
354                            xl = bunch[p].coordonnees[1][0] * ca + bunch[p].coordonnees[1][2] * sa;
355                            xsl = bunch[p].coordonnees[1][1] * ca + bunch[p].coordonnees[1][3] * sa;
356                            pdepth = abs(xl) - resColli[ipcoll[i]]->hgap;
357                            pdepth2 = abs(xl + lcoll * xsl) - resColli[ipcoll[i]]->hgap2;
358                        }
359
360                        if ((pdepth <= 0) && (pdepth2 <= 0)) { //we only continue with the particles that have not disappeared
361                            bunch[p].coordonnees[0][0] = bunch[p].coordonnees[1][0];
362                            bunch[p].coordonnees[0][1] = bunch[p].coordonnees[1][1];
363                            bunch[p].coordonnees[0][2] = bunch[p].coordonnees[1][2];
364                            bunch[p].coordonnees[0][3] = bunch[p].coordonnees[1][3];
365                        } else {
366                            bunch[p].in = false;
367                            count = count + 1;
368                            bunchhit.push_back(bunchtemp[p]);
369                            bunch[p].nrevhitp = k + 1;
370                            apdepth.push_back(pdepth);
371                            apdepth2.push_back(pdepth2);
372                            cocount[i] = cocount[i] + 1;
373                        }
374                    }//end (if i<niph)
375                    else {//we only continue with the particles that have not disappeared
376                        bunch[p].coordonnees[0][0] = bunch[p].coordonnees[1][0];
377                        bunch[p].coordonnees[0][1] = bunch[p].coordonnees[1][1];
378                        bunch[p].coordonnees[0][2] = bunch[p].coordonnees[1][2];
379                        bunch[p].coordonnees[0][3] = bunch[p].coordonnees[1][3];
380                    }
381                }//end if in
382
383                //the following part can be uncommented to have an output of the x- and  y-coordinates of a particle in the file coordinates.dat
384
385                /*if(bunch[p].getidentification() == choicePart){
386                  outCoord(bunch[p], i+1, "coordinates.dat");
387                  }*/
388
389            }//end loop over the particles
390        }//end loop over i
391
392        cout << "Revolution: " << k + 1 << ", number of particles left: " << total - count << endl;
393
394        if (total - count == 0) { //we return if all the particles are gone
395            return;
396        }
397
398        vector <Particle> tempinabs;
399
400        //we continue just with particles that are not lost
401        for (int w(0); w < bunch.size(); ++w) {
402            if (bunch[w].in != 0) {
403                tempinabs.push_back(bunch[w]);
404            }
405        }
406
407        bunch.clear();
408
409        for (int w(0); w < tempinabs.size(); ++w) {
410            bunch.push_back(tempinabs[w]);
411        }
412
413        tempinabs.clear();
414
415        for (int p(0); p < bunch.size(); ++p) {
416            if ((k + 1) % blowupperiod == 0) { //blowup/diffusion
417                if (bunch[p].in == 1) {
418                    //cout << "Blow-up!!" << k+1 << endl;
419                    int im;
420                    im = reseau.size() - 1;
421                    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];
422                    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];
423                    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];
424                    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];
425                }
426            }
427        }
428
429    }//end loop over k
430
431}
432
433
434void 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)
435{
436
437    if (idpart >= 0) {
438        if (idpart >= bunch.size()) {
439            cerr << "Warning, the particle that you want to spy using IDPART does not exist!" << endl;
440        } else {
441            choicePart = idpart;
442        }
443    }
444
445    if (idelt >= 0) {
446        if (idelt >= reseau.size()) {
447            cerr << "Warning, the element after which you want to spy using IDELT does not exist!" << endl;
448        } else {
449            eltOutNber = idelt;
450        }
451    }
452
453    ++turn;
454
455    if (nhitcolli.size() < resColli.size()) {
456        for (int k(0); k < resColli.size(); ++k) {
457            nhitcolli.push_back(0);
458        }
459    }
460
461    if (plotflag == "Yes") {
462        xco.clear();
463        yco.clear();
464        vector <double> temp1, temp2;
465        for (int k(0); k < bunch.size(); ++k) {
466            temp1.push_back(bunch[k].coordonnees[0][0]);
467            temp2.push_back(bunch[k].coordonnees[0][2]);
468        }
469
470        xco.push_back(temp1);
471        yco.push_back(temp2);
472
473        temp1.clear();
474        temp2.clear();
475    }
476
477    //we set the revolution number (trackensemblechrom is called two times during the first turn)
478    int rev;
479
480    if (irev == 0) {
481        rev = 1;
482    } else {
483        rev = irev;
484    }
485
486    if (irev == 0) {
487        cout << "Initialising R-matrix." << endl;
488
489        double K;
490        double cx, sx, cy, sy;
491
492        R11X.clear();
493        R12X.clear();
494        R21X.clear();
495        R22X.clear();
496        R11Y.clear();
497        R12Y.clear();
498        R21Y.clear();
499        R22Y.clear();
500        turn = -2;
501
502        Lelem.push_back(0);
503        R11X.push_back(0);
504        R12X.push_back(0);
505        R21X.push_back(0);
506        R22X.push_back(0);
507        R11Y.push_back(0);
508        R12Y.push_back(0);
509        R21Y.push_back(0);
510        R22Y.push_back(0);
511        R11XC.push_back(0);
512        R12XC.push_back(0);
513        R21XC.push_back(0);
514        R22XC.push_back(0);
515        R11YC.push_back(0);
516        R12YC.push_back(0);
517        R21YC.push_back(0);
518        R22YC.push_back(0);
519
520        for (int i(1); i < size_res; ++i) {
521
522            Lelem.push_back(reseau[i]->S - reseau[i - 1]->S);
523            if (reseau[i]->K1L > 0) {
524                if (Lelem[i] == 0) {
525                    Lelem[i] = 0.001;
526                }
527
528                K = reseau[i]->K1L / Lelem[i];
529                cx = cos(sqrt(K) * Lelem[i]);
530                sx = sin(sqrt(K) * Lelem[i]);
531                R11XC.push_back(0.5 * sqrt(K)*Lelem[i]*sx);
532                R12XC.push_back(-0.5 * Lelem[i]*cx + 0.5 / sqrt(K)*sx);
533                R21XC.push_back(0.5 * K * Lelem[i]*cx + 0.5 * sqrt(K)*sx);
534                R22XC.push_back(0.5 * sqrt(K)*Lelem[i]*sx);
535                cy = cosh(sqrt(K) * Lelem[i]);
536                sy = sinh(sqrt(K) * Lelem[i]);
537                R11YC.push_back(-0.5 * sqrt(K)*Lelem[i]*sy);
538                R12YC.push_back(-0.5 * Lelem[i]*cy + 0.5 / sqrt(K)*sy);
539                R21YC.push_back(-0.5 * K * Lelem[i]*cy - 0.5 * sqrt(K)*sy);
540                R22YC.push_back(-0.5 * sqrt(K)*Lelem[i]*sy);
541            } else if (reseau[i]->K1L < 0) {
542                if (Lelem[i] == 0) {
543                    Lelem[i] = 0.001;
544                }
545
546                K = -reseau[i]->K1L / Lelem[i];
547                cx = cosh(sqrt(K) * Lelem[i]);
548                sx = sinh(sqrt(K) * Lelem[i]);
549                R11XC.push_back(-0.5 * sqrt(K)*Lelem[i]*sx);
550                R12XC.push_back(-0.5 * Lelem[i]*cx + 0.5 / sqrt(K)*sx);
551                R21XC.push_back(-0.5 * K * Lelem[i]*cx - 0.5 * sqrt(K)*sx);
552                R22XC.push_back(-0.5 * sqrt(K)*Lelem[i]*sx);
553                cy = cos(sqrt(K) * Lelem[i]);
554                sy = sin(sqrt(K) * Lelem[i]);
555                R11YC.push_back(0.5 * sqrt(K)*Lelem[i]*sy);
556                R12YC.push_back(-0.5 * Lelem[i]*cy + 0.5 / sqrt(K)*sy);
557                R21YC.push_back(0.5 * K * Lelem[i]*cy + 0.5 * sqrt(K)*sy);
558                R22YC.push_back(0.5 * sqrt(K)*Lelem[i]*sy);
559            } else {
560                R11XC.push_back(0);
561                R12XC.push_back(0);
562                R21XC.push_back(0);
563                R22XC.push_back(0);
564                R11YC.push_back(0);
565                R12YC.push_back(0);
566                R21YC.push_back(0);
567                R22YC.push_back(0);
568            }
569
570            cx = cos(2 * M_PI * (reseau[i]->MUX - reseau[i - 1]->MUX));
571            sx = sin(2 * M_PI * (reseau[i]->MUX - reseau[i - 1]->MUX));
572            R11X.push_back(sqrt(reseau[i]->BETX / reseau[i - 1]->BETX) * (cx + reseau[i - 1]->ALFX * sx));
573            R12X.push_back(sqrt(reseau[i]->BETX * reseau[i - 1]->BETX)*sx);
574            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));
575            R22X.push_back(sqrt(reseau[i - 1]->BETX / reseau[i]->BETX) * (cx - reseau[i]->ALFX * sx));
576
577
578            cy = cos(2 * M_PI * (reseau[i]->MUY - reseau[i - 1]->MUY));
579            sy = sin(2 * M_PI * (reseau[i]->MUY - reseau[i - 1]->MUY));
580            R11Y.push_back(sqrt(reseau[i]->BETY / reseau[i - 1]->BETY) * (cy + reseau[i - 1]->ALFY * sy));
581            R12Y.push_back(sqrt(reseau[i]->BETY * reseau[i - 1]->BETY)*sy);
582            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));
583            R22Y.push_back(sqrt(reseau[i - 1]->BETY / reseau[i]->BETY) * (cy - reseau[i]->ALFY * sy));
584
585        }
586
587    }
588
589
590    // loop througt the elements in the machine
591
592    for (int i(i0 + 1); i <= im; ++i) { //i is the element number along the ring (similar numbering as s)
593
594        /*if(i == 2822){
595          read(bunch);
596          }*/
597
598        for (int p(0); p < bunch.size(); ++p) { //loop througt the particles
599
600            double dpopeff;
601
602            if (bunch[p].inabs == 1) {
603
604                //cout <<"Element " << i << endl;
605
606                int icop(-1);
607
608                for (int k(0); k < ips.size(); ++k) {
609                    if (i == ips[k]) { //element is a collimator
610                        icop = k;//The element is the 'icop'th collimator
611                    }
612                }
613
614                if (icop != -1) { //the element we consider is a collimator
615
616                    if (resColli[icop]->method == "standard") {
617
618                        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);
619
620                        if (bunch[p].Ap0 != 0) {
621                            bunch[p].coordonnees[1][5] = bunch[p].coordonnees[0][5] + time(bunch[p], reseau[i]->L, betgam, i);
622                        }
623
624                    }
625#if defined(FLUKA)
626                    else if (resColli[icop]->method == "fluka") {
627
628                        vector <Particle> temp;
629
630                        if (p == 0) { //we send all the bunch at the same time to Fluka only one time per element
631
632                            vector <Particle> bunchend;//bunch after the passage through Fluka
633
634                            for (int y(0); y < bunch.size(); ++y) {
635                                temp.push_back(bunch[y]);
636                            }
637
638                            resColli[icop]->collipassfluka(bunch, bunchend, conn, turn, momentum);
639
640
641                            nhitcolli[icop] = bunch.size() - bunchend.size(); //nhitcolli(icop) gives the number particles getting lost in collimator number icop,
642
643                            int creation(0);
644
645                            bunch.clear();
646
647                            for (int g(0); g < bunchend.size(); ++g) {
648                                if (bunchend[g].getidentification() == 0) {
649                                    bunch.push_back(bunchend[g]);
650                                    break;
651                                }
652                            }
653
654                            for (int g(0); g < bunchend.size(); ++g) {
655                                if (bunchend[g].getidentification() == bunch[bunch.size() - 1].getidentification() + 1) {
656                                    bunch.push_back(bunchend[g]);
657                                } else if (bunchend[g].getidentification() == bunch[bunch.size() - 1].getidentification() + 10000) {
658                                    bunch.push_back(bunchend[g]);
659                                    ++creation;
660                                } else if ((bunchend[g].getidentification() == bunch[bunch.size() - 1].getidentification()) && (bunchend[g].getidentification() != 0)) {
661                                    bunch.push_back(bunchend[g]);
662                                    ++creation;
663                                } else if (bunchend[g].getidentification() == bunch[bunch.size() - 1].getidentification() - 10000 + 1) {
664                                    bunch.push_back(bunchend[g]);
665                                }
666                            }
667
668                            nhitcolli[icop] = nhitcolli[icop] + creation;
669
670                            cout << creation << " particles created." << endl;
671
672                            int uu(1);
673                            int rr0, rr1;
674                            long double tps;
675                            bunch[0].coordonnees[0][5] = temp[0].coordonnees[0][5];
676                            rr0 = bunch[0].getidentification();
677                            bunch[0].setidentification(temp[0].getidentification());
678                            for (int g(1); g < bunch.size(); ++g) {
679                                rr1 = bunch[g].getidentification();
680                                if (bunch[g].getidentification() == rr0 + 1) {
681                                    bunch[g].setidentification(temp[uu].getidentification());
682                                    bunch[g].coordonnees[0][5] = temp[uu].coordonnees[0][5];
683                                    rr0 = rr1;
684                                } else {
685                                    tps = temp[uu - 1].coordonnees[0][5];
686                                    while (bunch[g].getidentification() == rr1) {
687                                        bunch[g].coordonnees[0][5] = tps;
688                                        ++g;
689                                    }
690                                    --g;
691                                    --uu;
692                                }
693                                ++uu;
694                            }
695
696                            for (int g(0); g < bunch.size(); ++g) {
697                                bunch[g].coordonnees[1][0] = bunch[g].coordonnees[0][0];
698                                bunch[g].coordonnees[1][1] = bunch[g].coordonnees[0][1];
699                                bunch[g].coordonnees[1][2] = bunch[g].coordonnees[0][2];
700                                bunch[g].coordonnees[1][3] = bunch[g].coordonnees[0][3];
701                                bunch[g].coordonnees[1][4] = bunch[g].coordonnees[0][4];
702                            }
703
704                        }
705
706                        bunch[p].coordonnees[1][5] = bunch[p].coordonnees[0][5] + time(bunch[p], reseau[i]->L, betgam, i);
707
708                        for (int k(0); k < bunch.size(); ++k) {
709                            bunch[k].inabs = 1;
710                        }
711
712                        temp.clear();
713
714                    }
715#endif
716                    else if (resColli[icop]->method == "magnetic") {
717
718                        if (p == 0) {
719                            resColli[icop]->hgap = resColli[icop]->hgap + resColli[icop]->deltaGap;
720                            resColli[icop]->hgap2 = resColli[icop]->hgap2 + resColli[icop]->deltaGap;
721                        }
722
723
724                        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);
725
726
727                        if (bunch[p].Ap0 != 0) {
728                            bunch[p].coordonnees[1][5] = bunch[p].coordonnees[0][5] + time(bunch[p], reseau[i]->L, betgam, i);
729                        }
730
731                    } else if (resColli[icop]->method == "crystal") {
732
733                        if (p == 0) { //we send all the bunch at the same time through collipassCrystal
734
735                            int lost(bunch.size());
736
737                            resColli[icop]->collipassCrystal(bunch, betgam, turn, outputpath);
738
739                            nhitcolli[icop] = lost - bunch.size(); //nhitcolli(icop) gives the number particles getting lost in collimator number icop,
740
741                        }
742
743                        bunch[p].coordonnees[1][5] = bunch[p].coordonnees[0][5] + time(bunch[p], reseau[i]->L, betgam, i);
744
745                    } else {
746
747                        cerr << "Error: unknown type of collimator, the method is not good defined!!" << endl;
748                    }
749
750                    //we test if the particle is lost in the preceeding collimator
751                    if (bunch[p].Ap0 == 0) {
752                        nhitcolli[icop] = nhitcolli[icop] + 1;
753                        bunch[p].inabs = 0;
754                    }
755
756                } else {//element not collimator
757
758                    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...)
759
760                        //Note thate the phase is taken here to be equal to pi.
761
762                        double period(rfharmonic / freqrf);
763                        double omega;
764                        double phase(0);
765                        double phi;
766                        double beta(sqrt(betgam * betgam / (betgam * betgam + 1))); //relativistic beta
767                        long double c(2.99792458e8);//speed of light [m/s]
768                        long double e(1.60218e-19);//elementary charge [C]
769
770
771                        omega = 2 * M_PI * (freqrf / rfharmonic);
772                        phi = phase + omega * bunch[p].coordonnees[0][5];
773
774                        bunch[p].coordonnees[1][0] = bunch[p].coordonnees[0][0];
775                        bunch[p].coordonnees[1][1] = bunch[p].coordonnees[0][1];
776                        bunch[p].coordonnees[1][2] = bunch[p].coordonnees[0][2];
777                        bunch[p].coordonnees[1][3] = bunch[p].coordonnees[0][3];
778
779                        double attr;
780                        attr = bunch[p].Zp0 * sin(phi) * (rfvoltage) / (beta * momentum);
781                        attr1 =  attr;
782                        bunch[p].coordonnees[1][4] = bunch[p].dpoporiginal + attr;//attention vraiment pas sur des parametre (surtout p dans la formule)
783
784                        bunch[p].coordonnees[1][5] = bunch[p].coordonnees[0][5];
785
786                        //uncomment the following line to have output related to the rf-cavity (cf lattice.h)
787                        //outrf(bunch[p].coordonnees[1][5], phi);
788
789
790                    }
791
792
793                    dpopeff = (bunch[p].Ap0 * Zpr) / (bunch[p].Zp0 * Apr) * (1 + bunch[p].coordonnees[0][4]) - 1;
794
795                    if (Lelem[i] == 0) {
796
797                        bunch[p].coordonnees[1][0] = bunch[p].coordonnees[0][0];
798                        bunch[p].coordonnees[1][1] = bunch[p].coordonnees[0][1];
799                        bunch[p].coordonnees[1][2] = bunch[p].coordonnees[0][2];
800                        bunch[p].coordonnees[1][3] = bunch[p].coordonnees[0][3];
801                        bunch[p].coordonnees[1][4] = bunch[p].coordonnees[0][4];
802                        bunch[p].coordonnees[1][5] = bunch[p].coordonnees[0][5];
803
804                    } else if ((reseau[i]->K1L != 0) && (nonlinflag == 1)) {
805
806                        double R11Xh, R12Xh, R21Xh, R22Xh, R11Yh, R12Yh, R21Yh, R22Yh;
807
808
809                        R11Xh = R11X[i] + R11XC[i] * dpopeff;
810                        R12Xh = R12X[i] + R12XC[i] * dpopeff;
811                        R21Xh = R21X[i] + R21XC[i] * dpopeff;
812                        R22Xh = R22X[i] + R22XC[i] * dpopeff;
813                        R11Yh = R11Y[i] + R11YC[i] * dpopeff;
814                        R12Yh = R12Y[i] + R12YC[i] * dpopeff;
815                        R21Yh = R21Y[i] + R21YC[i] * dpopeff;
816                        R22Yh = R22Y[i] + R22YC[i] * dpopeff;
817
818                        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;
819                        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;
820                        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;
821                        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;
822                        bunch[p].coordonnees[1][4] = bunch[p].coordonnees[0][4];
823
824                        bunch[p].coordonnees[1][5] = bunch[p].coordonnees[0][5] + time(bunch[p], reseau[i]->L, betgam, i);
825
826                    } else if ((reseau[i]->K2L != 0) && (nonlinflag == 1)) {
827                        double Lelemha(Lelem[i] / 2);
828                        double dxha(reseau[i - 1]->DX + reseau[i - 1]->DPX * Lelemha);
829                        double dyha(reseau[i - 1]->DY + reseau[i - 1]->DPY * Lelemha);
830
831                        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;
832                        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;
833
834                        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]);
835                        bunch[p].coordonnees[1][3] = bunch[p].coordonnees[0][3] + reseau[i]->K2L * (bunch[p].coordonnees[1][0] * bunch[p].coordonnees[1][2]);
836                        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;
837                        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;
838                        bunch[p].coordonnees[1][4] = bunch[p].coordonnees[0][4];
839
840                        bunch[p].coordonnees[1][5] = bunch[p].coordonnees[0][5] + time(bunch[p], reseau[i]->L, betgam, i);
841                    } else {
842
843                        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;
844                        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;
845                        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;
846                        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;
847
848                        if ((reseau[i]->KEYWORD != "RFCAVITY") && (RFflag == 1)) {
849                            bunch[p].coordonnees[1][4] = bunch[p].coordonnees[0][4];
850                        }
851
852                        bunch[p].coordonnees[1][5] = bunch[p].coordonnees[0][5] + time(bunch[p], reseau[i]->L, betgam, i);
853
854                    }
855
856                }
857
858                if ((nonlinflag == 1) && ((icop == -1) || (resColli[icop]->method == "standard") || (resColli[icop]->method == "magnetic")) && (bunch[p].inabs != 0)) {
859
860                    int ixcor(-1), iycor (-1);
861
862                    for (int k(0); k < ixcormag.size(); ++k) {
863                        if (ixcormag[k] == i) {
864                            ixcor = k;
865                        }
866                    }
867
868                    if (ixcor != -1) {
869                        bunch[p].coordonnees[1][0] = bunch[p].coordonnees[1][0] + scaleorbit * (reseau[i]->S - reseau[i - 1]->S) / 2 * xcormag[ixcor] * (1 + dpopeff);
870                        bunch[p].coordonnees[1][1] = bunch[p].coordonnees[1][1] + scaleorbit * xcormag[ixcor] * (1 + dpopeff);
871                    }
872
873                    for (int k(0); k < iycormag.size(); ++k) {
874                        if (iycormag[k] == i) {
875                            iycor = k;
876                        }
877                    }
878
879                    if (iycor != -1) {
880                        bunch[p].coordonnees[1][2] = bunch[p].coordonnees[1][2] + scaleorbit * (reseau[i]->S - reseau[i - 1]->S) / 2 * ycormag[iycor] * (1 + dpopeff);
881                        bunch[p].coordonnees[1][3] = bunch[p].coordonnees[1][3] + scaleorbit * ycormag[iycor] * (1 + dpopeff);
882                    }
883
884                }
885
886
887
888                //checking for aperture hits
889                if (bunch[p].inabs == 1) {
890
891                    if ((icop == -1) || (resColli[icop]->method == "standard") || (resColli[icop]->method == "magnetic") || (resColli[icop]->method == "crystal")) {
892
893                        bool inside(false);
894
895                        if (reseau[i]->APERTYPE == "RECTANGLE") {
896                            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)) {
897                                inside = true;
898                            }
899                        } else if ((reseau[i]->APERTYPE == "ELLIPSE") || (reseau[i]->APERTYPE == "CIRCLE") || (reseau[i]->APERTYPE == "RECTELLIPSE")) {
900                            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)) {
901                                inside = true;
902
903                            }
904                        } else {
905                            cerr << "Error for the " << i << "th element: unknown aperture type!" << endl;
906                        }
907
908                        double L;
909
910                        L = reseau[i]->S - reseau[i - 1]->S; //distance to next element in the accelerator
911
912                        double Ap0h, Zp0h, xh0, xh, xhs, yh0, yh, yhs;
913                        int nh(0);//number of lost particles
914
915                        if (inside == false) {
916                            //cout <<"The particle is lost at element " << reseau[i]->NAME << endl;
917                            xh0 = bunch[p].coordonnees[0][0];
918                            bunch[p].inabs = 0;
919
920                            if (L == 0) {
921                                nh = nh + 1;
922                                hits.push_back(reseau[i]->S);//saving the value of s where the particles get lost
923                                Aphit.push_back(bunch[p].Ap0);//saving the mass of the lost particles
924                                Zphit.push_back(bunch[p].Zp0);//saving the charge of the lost particles
925                            } else {
926
927                                xh = bunch[p].coordonnees[1][0];
928                                xhs = (xh - xh0) / L;
929
930                                yh0 = bunch[p].coordonnees[0][2];
931                                yh = bunch[p].coordonnees[1][2];
932                                yhs = (yh - yh0) / L;
933
934                                Aphit.push_back(bunch[p].Ap0);
935                                Zphit.push_back(bunch[p].Zp0);
936
937                                double slostx, slosty, splus, sminus;
938
939                                if (reseau[i]->APERTYPE == "RECTANGLE") {
940
941                                    //find hit position in x
942
943                                    if ((abs(xh0) < reseau[i]->aperx) && (abs(xh) < reseau[i]->aperx)) { //both points inside aperture - particle lost in y instead
944                                        slostx = reseau[i]->S;
945
946                                    } else if (abs(xh0) > reseau[i]->aperx) { //first point outside - particle lost already at the entrance of the element
947                                        slostx = reseau[i - 1]->S;
948
949                                    } else {//first point inside, second point outside - do a linear interpolation
950
951                                        splus = (reseau[i]->aperx - xh0) / xhs;
952                                        sminus = (-reseau[i]->aperx - xh0) / xhs;
953
954                                        if (splus > sminus) { //choose the correct point where the straight line hits the aperture, at +-a. The highest s-value is correct
955                                            slostx = splus + reseau[i - 1]->S;
956
957                                        } else {
958                                            slostx = sminus + reseau[i - 1]->S;
959
960                                        }
961                                    }
962
963                                    //find hit position in y
964
965                                    if ((abs(yh0) < reseau[i]->apery) && (abs(yh) < reseau[i]->apery)) { //both points inside aperture - particle lost in x instead
966                                        slosty = reseau[i]->S;
967                                    } else if (abs(yh0) > reseau[i]->apery) { //first point outside - particle lost already at the entrance of the element
968                                        slosty = reseau[i - 1]->S;
969                                    } else {//first point inside, second point outside - do a linear interpolation
970
971                                        splus = (reseau[i]->apery - yh0) / yhs;
972                                        sminus = (-reseau[i]->apery - yh0) / yhs;
973
974                                        if (splus > sminus) { //choose the right point where the straight line hits the aperture, at +-b. The highest s-value is correct
975                                            slosty = splus + reseau[i - 1]->S;
976                                        } else {
977                                            slosty = sminus + reseau[i - 1]->S;
978                                        }
979                                    }
980
981                                    //choose the point where the particle hits first, x or y
982
983                                    if (slostx < slosty) {
984                                        hits.push_back(slostx);
985                                    } else {
986                                        hits.push_back(slosty);
987                                    }
988                                } else if (reseau[i]->APERTYPE == "ELLIPSE") {
989
990                                    double sqrarg;
991                                    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;
992                                    if (sqrarg < 0) {
993                                        cout << "Warning, sqrarg < 0" << endl;
994                                        hits.push_back(reseau[i - 1]->S);
995                                    } else {
996                                        double stry;
997
998                                        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);
999
1000                                        if ((stry > 0) && (stry < L)) {
1001                                            hits.push_back(reseau[i - 1]->S + stry);
1002                                        } else {
1003
1004                                            double stry2;
1005
1006                                            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);
1007
1008                                            if ((stry2 > 0) && (stry2 < L)) {
1009                                                hits.push_back(reseau[i - 1]->S + stry2);
1010                                            } else {
1011                                                hits.push_back(reseau[i - 1]->S);
1012                                            }
1013                                        }
1014                                    }
1015                                }//end if ELLIPSE
1016                                else if (reseau[i]->APERTYPE == "RECTELLIPSE") {
1017
1018                                    double sqrarg;
1019                                    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;
1020                                    if (sqrarg < 0) {
1021                                        cout << "Warning, sqrarg < 0" << endl;
1022                                        hits.push_back(reseau[i - 1]->S);
1023                                    } else {
1024                                        double stry;
1025                                        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);
1026
1027                                        if ((stry > 0) && (stry < L)) {
1028                                            hits.push_back(reseau[i - 1]->S + stry);
1029                                        } else {
1030                                            double stry2;
1031
1032                                            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);
1033
1034                                            if ((stry2 > 0) && (stry2 < L)) {
1035                                                hits.push_back(reseau[i - 1]->S + stry2);
1036                                            } else {
1037                                                hits.push_back(reseau[i - 1]->S);
1038                                            }
1039                                        }
1040                                    }
1041                                }//end RECTELLIPSE
1042                            }//end else after if L == 0
1043
1044                        }//end si particle lost
1045                        else {
1046                            //cout << "No losses, all the particles stay in the experience." << endl;
1047                        }
1048                    }//end if icop==vide ||...
1049
1050                }
1051            }
1052
1053            //Different ways to print the coordinates of the particle through the parameters IDPART and OUTCOORD (see manual)
1054
1055            if ((outcoord != 0) && (bunch[p].inabs != 0)) {
1056
1057                if (outcoord == 1) {
1058
1059                    cout << "Particle " <<  p + 1 << " after the passage through element: " << i << " " << reseau[i]->NAME << endl;
1060                    bunch[p].afficheCoordonnees();
1061                    cout << "Turn number " << turn << endl;
1062                    cout << "Massnumber: " << bunch[p].Ap0 << " Chargestate: " << bunch[p].Zp0 << endl;
1063                    cout << "Particle ID: " << bunch[p].getidentification() << endl;
1064                } else if (outcoord == 3) {
1065
1066                    if (bunch[p].getidentification() == choicePart) {
1067                        outCoord(bunch[p], i, outputpath + "/coordinates.dat");
1068                    }
1069                } else if (outcoord == 4) {
1070
1071                    if (bunch[p].getidentification() == choicePart) {
1072                        outPunct(i, bunch[p], bunch[0], attr1, outputpath);
1073                    }
1074                }
1075
1076                else if (outcoord == 2) {
1077                    outElt(i, bunch[p], outputpath, indication);
1078                }
1079            }
1080
1081
1082        }//end loop over particles
1083
1084
1085        //we prepare the particles for the next element
1086        for (int b(0); b < bunch.size(); ++b) {
1087            if (bunch[b].inabs != 0) {
1088                for (int r(0); r < 6; ++r) {
1089                    bunch[b].coordonnees[0][r] = bunch[b].coordonnees[1][r];
1090                }
1091            }
1092        }
1093
1094        if (plotflag == "Yes") {
1095
1096            vector <double> temp1, temp2;
1097
1098            for (int h(0); h < bunch.size(); ++h) {
1099                if (bunch[h].inabs != 0) {
1100                    temp1.push_back(bunch[h].coordonnees[0][0]);
1101                    temp2.push_back(bunch[h].coordonnees[0][2]);
1102                } else {
1103                    temp1.push_back(100);
1104                    temp2.push_back(100);
1105                }
1106            }
1107
1108            xco.push_back(temp1);
1109            yco.push_back(temp2);
1110
1111            temp1.clear();
1112            temp2.clear();
1113        }
1114
1115
1116        //if there are no more particles, we return
1117        if (bunch.size() == 0) {
1118            return;
1119        }
1120
1121    }//end loop elements
1122
1123    vector <Particle> tempinabs;
1124
1125    //we contine just with particles that are not lost
1126    for (int w(0); w < bunch.size(); ++w) {
1127        if (bunch[w].inabs != 0) {
1128            tempinabs.push_back(bunch[w]);
1129        }
1130    }
1131
1132    bunch.clear();
1133
1134    for (int w(0); w < tempinabs.size(); ++w) {
1135        bunch.push_back(tempinabs[w]);
1136    }
1137
1138    tempinabs.clear();
1139
1140}
Note: See TracBrowser for help on using the repository browser.