← Back to branch summary

~max-mcginley/third-order-lie-map/generator

« back to all changes in this revision

Viewing changes to BASEparticle.cc

  • Committer: Max McGinley
  • Date: 2015-09-07 10:20:10 UTC
  • Revision ID: max.mcginley@stfc.ac.uk-20150907102010-0bot9845lf546xwc
Added a AnalyticLattice class in solenoid_models which uses an equation which has an equivalent MAUS::BTField. Also added function for tracking progress of particles round unit cell

Show diffs side-by-side

added added

removed removed

Lines of Context:
BASEparticle.cc
15
15
 
16
16
 
17
17
namespace MAUS {
18
 
    BASEparticle::BASEparticle() : TransferMatrixGenerator("BASEparticle") {
 
18
    BASEparticle::BASEparticle() : TransferMatrixGenerator("BASEparticle"), usingLattice(false) {
19
19
        getCurrentStream() << "Constructing BASEparticle" << std::endl;
20
20
        matrixGen = false;
21
21
        particlesRun = false;
22
22
        mass = 105.68;
23
 
        refPart = static_cast<MAUSPrimaryGeneratorAction::PGParticle*>(malloc(sizeof(MAUSPrimaryGeneratorAction::PGParticle)));
24
 
        polynomialMatrix = static_cast<Matrix<double>*>(malloc(sizeof(Matrix<double>)));
25
 
        dataOut = NULL;
26
 
        refDataOut = NULL;
27
 
        particles = NULL;
28
 
        thirdOrderMap = NULL;
 
23
        refPart = NULL; polynomialMatrix = NULL; dataOut = NULL; refDataOut = NULL; particles = NULL; thirdOrderMap = NULL;
29
24
        bField = Globals::GetMCFieldConstructor()->GetMagneticField();
30
25
        integrator = new BfieldIntegratorMap(bField);
31
26
        getCurrentStream() << "Memory allocated" << std::endl;
32
27
        setDataOut(Json::Value(Json::nullValue));
33
28
        setRefDataOut(Json::Value(Json::nullValue));
34
29
        setParticles(Json::Value(Json::nullValue));
35
 
        getCurrentStream() << "Values set" << std::endl;
36
 
        
 
30
        _numberOfParticles = 0;
 
31
    }
 
32
    
 
33
    BASEparticle::BASEparticle(AnalyticLattice* lattice) : TransferMatrixGenerator("BASEparticle"), usingLattice(true), _lattice(lattice) {
 
34
        getCurrentStream() << "Constructing BASEparticle using lattice" << std::endl;
 
35
        matrixGen = false;
 
36
        particlesRun = false;
 
37
        mass = 105.68;
 
38
        refPart = NULL; polynomialMatrix = NULL; dataOut = NULL; refDataOut = NULL; particles = NULL; thirdOrderMap = NULL;
 
39
        bField = lattice->getFieldGroup();
 
40
        integrator = new BfieldIntegratorMap(bField);
 
41
        getCurrentStream() << "Memory allocated" << std::endl;
 
42
        setDataOut(Json::Value(Json::nullValue));
 
43
        setRefDataOut(Json::Value(Json::nullValue));
 
44
        setParticles(Json::Value(Json::nullValue));
 
45
        _numberOfParticles = 0;
37
46
    }
38
47
    
39
48
    BASEparticle::~BASEparticle() {
40
 
        free(refPart);
41
 
                free(polynomialMatrix);
42
 
        delete polynomialMatrix;
 
49
        if (refPart != NULL) delete refPart;
 
50
        if (polynomialMatrix != NULL) delete polynomialMatrix;
43
51
        delete integrator;
44
52
    }
45
53
    
 
54
    MAUSPrimaryGeneratorAction::PGParticle BASEparticle::getRefPart() const {if (refPart == NULL) throw std::runtime_error("BASEparticle::getRefPart is NULL"); return *refPart;}
 
55
    Json::Value BASEparticle::getParticles() const {if (particles == NULL) throw std::runtime_error("BASEparticle::getParticles is NULL"); return *particles;}
 
56
    BTField* BASEparticle::getBField() const {if (bField == NULL) throw std::runtime_error("BASEparticle::getBField is NULL"); return bField;}
 
57
    BfieldIntegrator* BASEparticle::getIntegrator() const {if (integrator == NULL) throw std::runtime_error("BASEparticle::getIntegrator is NULL"); return integrator;}
 
58
    Json::Value BASEparticle::getDataOut() const {if (dataOut == NULL) throw std::runtime_error("BASEparticle::getDataOut is NULL"); return *dataOut;}
 
59
    Json::Value BASEparticle::getRefDataOut() const {if (refDataOut == NULL) throw std::runtime_error("BASEparticle::getRefDataOut is NULL"); return *refDataOut;}
 
60
    PolynomialMap BASEparticle::getThirdOrderMap() const {if (thirdOrderMap == NULL) throw std::runtime_error("BASEparticle::getThirdOrderMap is NULL"); return PolynomialMap(*thirdOrderMap);}
 
61
    AnalyticLattice* BASEparticle::getLattice() const {if (usingLattice == false) throw std::runtime_error("BASEparticle::getLattice() is not in usingLattice mode. Use constructor to pass AnalyticLattice* object"); if (_lattice == NULL) throw std::runtime_error("BASEparticle::getLattice is NULL"); return _lattice;}
 
62
    int BASEparticle::getNumberOfParticles() const {return _numberOfParticles;}
 
63
    
46
64
    std::vector<PolynomialMap::PolynomialCoefficient> BASEparticle::get2ndOrderBlanks() {
47
65
        std::vector<PolynomialMap::PolynomialCoefficient> v;
48
66
        for (int i = 0; i < 4; i++) {
57
75
        return v;
58
76
    }
59
77
    
60
 
    void BASEparticle::runParticles() {
61
 
        getCurrentStream() << "Running reference particle" << std::endl;
62
 
        Json::Value refOut = MAUSGeant4Manager::GetInstance()->RunParticle(getRefPart());
63
 
        setRefDataOut(refOut);
64
 
        getCurrentStream() << "Running all particles" << std::endl;
65
 
        Json::Value out = MAUSGeant4Manager::GetInstance()->RunManyParticles(getParticles());
66
 
        setDataOut(out);
67
 
        particlesRun = true;
68
 
        getCurrentStream() << "All particles have been run" << std::endl;
 
78
    void BASEparticle::runParticles(double zOut) {
 
79
        if (usingLattice) {
 
80
                getCurrentStream() << "Running reference particle" << std::endl;
 
81
                trackedReference = new double[8];
 
82
                for (int i = 0; i < 8; i++)
 
83
                        trackedReference[i] = inReference[i];
 
84
                BTTracker::integrate(zOut,trackedReference,getLattice()->getFieldGroup(),BTTracker::z, stepSize, 1.0);
 
85
                getCurrentStream() << "Running all particles" << std::endl;
 
86
                trackedParticles = new double[8*numberOfParticles];
 
87
                for (int i = 0; i < 8*numberOfParticles; i++)
 
88
                        trackedParticles[i] = inParticles[i];
 
89
                for (int i = 0; i < numberOfParticles; i++)
 
90
                        BTTracker::integrate(zOut,trackedParticles + (i*8),getLattice()->getFieldGroup(),BTTracker::z, stepSize, 1.0);
 
91
                particlesRun = true;
 
92
                /*
 
93
                for (int i = 0; i < numberOfParticles; i++) {
 
94
                        getCurrentStream() << "Particle in:\n\t";
 
95
                        for (int j = 0; j < 8; j++)
 
96
                                getCurrentStream() << inParticles[(i*8) + j] << ", ";
 
97
                        getCurrentStream() << std::endl << "Particle out:\n\t";
 
98
                        for (int j = 0; j < 8; j++)
 
99
                                getCurrentStream() << trackedParticles[(i*8) + j] << ", ";
 
100
                        getCurrentStream() << std::endl << std::endl;
 
101
                }
 
102
                */
 
103
                getCurrentStream() << "All particles have been run" << std::endl;
 
104
                }
 
105
                else {
 
106
                        getCurrentStream() << "Running reference particle" << std::endl;
 
107
                Json::Value refOut = MAUSGeant4Manager::GetInstance()->RunParticle(getRefPart());
 
108
                setRefDataOut(refOut);
 
109
                getCurrentStream() << "Running all particles" << std::endl;
 
110
                Json::Value out = MAUSGeant4Manager::GetInstance()->RunManyParticles(getParticles());
 
111
                setDataOut(out);
 
112
                particlesRun = true;
 
113
                getCurrentStream() << "All particles have been run" << std::endl;
 
114
                }
69
115
    }
70
116
    
71
117
    PolynomialMap BASEparticle::buildMap(double zIn, double zOut, bool thirdOrder) {
74
120
                return PolynomialMap(std::vector<PolynomialMap::PolynomialCoefficient>());
75
121
        }
76
122
        
77
 
        getCurrentStream() << "Generating particle in/out data";
78
 
        static int count = 0;
79
 
        getCurrentStream() << "..."  << std::endl;
80
 
        Json::Value refOut = getRefDataOut();
81
 
        double refMom = refOut["primary"]["momentum"]["z"].asDouble();
82
 
        Json::Value out = getDataOut();
83
 
        for(int j = 0; j < out.size(); j++) {
84
 
                loadBar(j,out.size(),50,getCurrentStream());
85
 
                Json::Value vhit = out[j]["virtual_hits"];
86
 
                std::pair<std::vector<double>, std::vector<double> > coords = getVirtualCoordsFromHit(vhit,refMom,zIn,zOut,thirdOrder);
87
 
                partCoordsIn.push_back(coords.first);
88
 
                partCoordsOut.push_back(coords.second);
89
 
        }
90
 
        resetLoadBar();
91
 
        getCurrentStream() << "\n...Generated map" << std::endl;
92
 
        
 
123
        std::vector< std::vector<double> > partCoordsIn;
 
124
        std::vector< std::vector<double> > partCoordsOut;
 
125
        
 
126
        if (usingLattice) {
 
127
                double refMom = trackedReference[7];
 
128
                for(int j = 0; j < getNumberOfParticles(); j++) {
 
129
                        loadBar(j,getNumberOfParticles(),50,getCurrentStream());
 
130
                        std::pair<std::vector<double>, std::vector<double> > coords = getVirtualCoordsFromTracked(inParticles + (j*8), trackedParticles + (j*8),refMom,thirdOrder);
 
131
                        partCoordsIn.push_back(coords.first);
 
132
                        partCoordsOut.push_back(coords.second);
 
133
                }
 
134
                resetLoadBar();
 
135
                getCurrentStream() << "\n...Generated map" << std::endl;
 
136
                }
 
137
                else {
 
138
                        getCurrentStream() << "Generating particle in/out data";
 
139
                static int count = 0;
 
140
                getCurrentStream() << "..."  << std::endl;
 
141
                Json::Value refOut = getRefDataOut();
 
142
                double refMom = refOut["primary"]["momentum"]["z"].asDouble();
 
143
                Json::Value out = getDataOut();
 
144
                for(int j = 0; j < out.size(); j++) {
 
145
                        loadBar(j,out.size(),50,getCurrentStream());
 
146
                        Json::Value vhit = out[j]["virtual_hits"];
 
147
                        std::pair<std::vector<double>, std::vector<double> > coords = getVirtualCoordsFromHit(vhit,refMom,zIn,zOut,thirdOrder);
 
148
                        partCoordsIn.push_back(coords.first);
 
149
                        partCoordsOut.push_back(coords.second);
 
150
                }
 
151
                resetLoadBar();
 
152
                getCurrentStream() << "\n...Generated map" << std::endl;
 
153
                }
93
154
        if(thirdOrder) return *PolynomialMap::ConstrainedPolynomialLeastSquaresFit(partCoordsIn, partCoordsOut, 3,get2ndOrderBlanks());
94
155
        return *PolynomialMap::PolynomialLeastSquaresFit(partCoordsIn, partCoordsOut, 1);
95
156
    }
158
219
        return retval;
159
220
    }
160
221
    
161
 
    void BASEparticle::generateParticles() {
 
222
    void BASEparticle::generateParticles(double zIn) {
162
223
        getCurrentStream() << "Generating " << numberOfParticles << " particles..." << std::endl;
163
 
        generateReferenceParticle();
164
 
        Json::Value retval(Json::arrayValue);
165
 
        for(int i = 0; i < numberOfParticles; i++) {
166
 
                Json::Value val = getRefPart().WriteJson();
167
 
                std::vector<double> offsets = generateRandomOffsets();
168
 
                double x = val["position"]["x"].asDouble() + offsets[0]*sigma;
169
 
                double y = val["position"]["y"].asDouble() + offsets[1]*sigma;
170
 
                double z = val["position"]["z"].asDouble();
171
 
                double px = val["momentum"]["x"].asDouble() + offsets[2]*val["momentum"]["z"].asDouble();
172
 
                double py = val["momentum"]["y"].asDouble() + offsets[3]*val["momentum"]["z"].asDouble();
173
 
                double pz = std::sqrt(val["momentum"]["z"].asDouble()*val["momentum"]["z"].asDouble() - px*px - py*py);
174
 
                val["position"]["x"] = x;
175
 
                val["position"]["y"] = y;
176
 
                //val["position"]["z"] = z;
177
 
                val["momentum"]["x"] = px;
178
 
                val["momentum"]["y"] = py;
179
 
                //val["momentum"]["z"] = pz;
180
 
                Json::Value primary(Json::objectValue);
181
 
                primary["primary"] = val;
182
 
                retval[i] = primary;
183
 
        }
184
 
        setParticles(retval);
185
 
        getCurrentStream() << "Reference particle: " << std::endl << JsonWrapper::JsonToString(getRefPart().WriteJson()) << std::endl;
 
224
        generateReferenceParticle(zIn);
 
225
        if (!usingLattice) {
 
226
                Json::Value retval(Json::arrayValue);
 
227
                for(int i = 0; i < numberOfParticles; i++) {
 
228
                        Json::Value val = getRefPart().WriteJson();
 
229
                        std::vector<double> offsets = generateRandomOffsets();
 
230
                        double x = val["position"]["x"].asDouble() + offsets[0]*sigma;
 
231
                        double y = val["position"]["y"].asDouble() + offsets[1]*sigma;
 
232
                        double px = val["momentum"]["x"].asDouble() + offsets[2]*val["momentum"]["z"].asDouble();
 
233
                        double py = val["momentum"]["y"].asDouble() + offsets[3]*val["momentum"]["z"].asDouble();
 
234
                        double pz = std::sqrt(val["energy"].asDouble()*val["energy"].asDouble() - px*px - py*py);
 
235
                        val["position"]["x"] = x;
 
236
                        val["position"]["y"] = y;
 
237
                        val["position"]["z"] = zIn;
 
238
                        val["momentum"]["x"] = px;
 
239
                        val["momentum"]["y"] = py;
 
240
                        val["momentum"]["z"] = pz;
 
241
                        Json::Value primary(Json::objectValue);
 
242
                        primary["primary"] = val;
 
243
                        retval[i] = primary;
 
244
                }
 
245
                setParticles(retval);
 
246
                _numberOfParticles = numberOfParticles;
 
247
                getCurrentStream() << "Reference particle: " << std::endl << JsonWrapper::JsonToString(getRefPart().WriteJson()) << std::endl;
 
248
        }
 
249
        else {
 
250
                double* parts = new double[8*numberOfParticles];
 
251
                for(int i = 0; i < numberOfParticles; i++) {
 
252
                        std::vector<double> offsets = generateRandomOffsets();
 
253
                        double x = inReference[1] + offsets[0]*sigma;
 
254
                        double y = inReference[2] + offsets[1]*sigma;
 
255
                        double px = inReference[5] + offsets[2]*inReference[7];
 
256
                        double py = inReference[6] + offsets[3]*inReference[7];
 
257
                        double pz = std::sqrt(refEnergy*refEnergy - px*px - py*py);
 
258
                        parts[i*8] = 0.0;
 
259
                        parts[i*8 + 1] = x;
 
260
                        parts[i*8 + 2] = y;
 
261
                        parts[i*8 + 3] = zIn;
 
262
                        parts[i*8 + 4] = refEnergy;
 
263
                        parts[i*8 + 5] = px;
 
264
                        parts[i*8 + 6] = py;
 
265
                        parts[i*8 + 7] = pz;
 
266
                }
 
267
                _numberOfParticles = numberOfParticles;
 
268
                inParticles = parts;
 
269
        }
 
270
        
186
271
    }
187
272
    
188
 
    Matrix<double> BASEparticle::getPolynomialMatrix() {
189
 
        if(!matrixGen) {
190
 
                getCurrentStream() << "Cannot access polynomial matrix - not initialised" << std::endl;
 
273
    Matrix<double> BASEparticle::getPolynomialMatrix() const {
 
274
        if(!matrixGen || polynomialMatrix == NULL) {
 
275
                getCurrentStream() << "Cannot access polynomial matrix - not initialised using setPolynomialMatrix() method" << std::endl;
191
276
                Matrix<double> blank;
192
277
                return blank;
193
278
        }
194
279
        return *polynomialMatrix;
195
280
    }
196
281
    
197
 
    void BASEparticle::generateReferenceParticle() {
198
 
        setRefPart((MAUSGeant4Manager::GetInstance())->GetReferenceParticle());
 
282
    void BASEparticle::generateReferenceParticle(double zIn) {
 
283
        if (usingLattice) {
 
284
                double* ret = new double[8];
 
285
                ret[0] = 0.0; ret[1] = 0.0; ret[2] = 0.0; ret[5] = 0.0; ret[6] = 0.0;
 
286
                ret[3] = zIn;
 
287
                ret[4] = refEnergy;
 
288
                ret[7] = refMomentum;
 
289
                inReference = ret;
 
290
        }
 
291
        else setRefPart((MAUSGeant4Manager::GetInstance())->GetReferenceParticle());
199
292
    }
200
293
    
201
294
    double* BASEparticle::generateRandomCoordinates() {
204
297
        return ret;
205
298
    }
206
299
    
 
300
    void BASEparticle::advanceCell(double* coordsIn, double* coordsOut, double zIn, double zOut) {
 
301
        for (int i = 0; i < 4; i++)
 
302
                coordsOut[i] = 0.0;
 
303
        if (!usingLattice) return;
 
304
        double* properCoords = new double[8];
 
305
        properCoords[0] = 0.;
 
306
        properCoords[1] = coordsIn[0]*sigma;
 
307
        properCoords[2] = coordsIn[1]*sigma;
 
308
        std::vector<double> A1 = integrator->getVectorPotentialAnalytic(properCoords[0],properCoords[1],zIn);
 
309
        properCoords[3] = zIn;
 
310
        properCoords[4] = refEnergy;
 
311
        properCoords[5] = coordsIn[2]*refMomentum - A1[0]*lightSpeed/1.0e-6;
 
312
        properCoords[6] = coordsIn[3]*refMomentum - A1[1]*lightSpeed/1.0e-6;
 
313
        properCoords[7] = std::sqrt(refEnergy*refEnergy - properCoords[5]*properCoords[5] - properCoords[6]*properCoords[6]);
 
314
        BTTracker::integrate(zOut,properCoords,getLattice()->getFieldGroup(),BTTracker::z, stepSize, 1.0);
 
315
        std::vector<double> A2 = integrator->getVectorPotentialAnalytic(properCoords[0],properCoords[1],zOut);
 
316
        coordsOut[0] = properCoords[1]/sigma;
 
317
        coordsOut[1] = properCoords[2]/sigma;
 
318
        coordsOut[2] = (properCoords[5] + A2[0]*lightSpeed/1.0e-6)/sigma;
 
319
        coordsOut[3] = (properCoords[6] + A2[0]*lightSpeed/1.0e-6)/sigma;
 
320
        delete[] properCoords;
 
321
    }
 
322
    
207
323
    /*void BASEparticle::TESTThirdOrderTrack() {
208
324
        TestResult* test = new TestResult("Third Order Tracking Comparison");
209
325
        for (int i = 0; i < 5; i++) {
322
438
        return retval;
323
439
}
324
440
 
 
441
std::pair<std::vector<double>, std::vector<double> > BASEparticle::getVirtualCoordsFromTracked(double* particlesIn, double* particlesOut, double refMom, bool thirdOrder) {
 
442
        //initial kinetic variables
 
443
        double x1 = particlesIn[1];
 
444
        double y1 = particlesIn[2];
 
445
        double px1 = particlesIn[5];
 
446
        double py1 = particlesIn[6];
 
447
        double z1 = particlesIn[3];
 
448
        
 
449
        //final kinetic variables
 
450
        double x2 = particlesOut[1];
 
451
        double y2 = particlesOut[2];
 
452
        double px2 = particlesOut[5];
 
453
        double py2 = particlesOut[6];
 
454
        double z2 = particlesOut[3];
 
455
        
 
456
        //conversion to canonical momenta
 
457
        std::vector<double> A1 = integrator->getVectorPotentialAnalytic(x1,y1,z1);
 
458
        std::vector<double> A2 = integrator->getVectorPotentialAnalytic(x2,y2,z2);
 
459
        
 
460
                //from T m * {fundamental charge unit} -> MeV/c is lightSpeed (m s^-1) / 1,000,000
 
461
        px1 += A1[0] * lightSpeed / 1.0e6;
 
462
        py1 += A1[1] * lightSpeed / 1.0e6;
 
463
        px2 += A2[0] * lightSpeed / 1.0e6;
 
464
        py2 += A2[1] * lightSpeed / 1.0e6;
 
465
        
 
466
        //rescaling to dimensionless MARYLIE variables
 
467
        x1 /= sigma;
 
468
        y1 /= sigma;
 
469
        x2 /= sigma;
 
470
        y2 /= sigma;
 
471
        px1 /= refMom;
 
472
        py1 /= refMom;
 
473
        px2 /= refMom;
 
474
        py2 /= refMom;
 
475
        
 
476
        /*getCurrentStream() << "Particle IN: z = " << step1["position"]["z"] << std::endl;
 
477
        getCurrentStream() << "x = " << x1 << " y = " << y1 << " px = " << px1 << " py = " << py1 << std::endl;
 
478
        getCurrentStream() << "Particle OUT: z = " << step2["position"]["z"] << std::endl;
 
479
        getCurrentStream() << "x = " << x2 << " y = " << y2 << " px = " << px2 << " py = " << py2 << "\n" << std::endl;*/
 
480
        
 
481
        std::vector<double> phaseCoordIn(4);
 
482
        std::vector<double> phaseCoordOut(4);
 
483
       
 
484
        phaseCoordIn[0] = x1;
 
485
        phaseCoordIn[1] = y1;
 
486
        phaseCoordIn[2] = px1;
 
487
        phaseCoordIn[3] = py1;
 
488
 
 
489
        phaseCoordOut[0] = x2;
 
490
        phaseCoordOut[1] = y2;
 
491
        phaseCoordOut[2] = px2;
 
492
        phaseCoordOut[3] = py2;
 
493
        
 
494
        std::pair<std::vector<double>, std::vector<double> > retval(phaseCoordIn, phaseCoordOut);
 
495
        
 
496
        return retval;
 
497
}
 
498
 
325
499
/*
326
500
void BASEparticle::TESTConfig(Json::Value config) {
327
501
        TestResult* test = new TestResult("Json Config");
334
508
}*/
335
509
 
336
510
 
337
 
double BASEparticle::float_tolerance() {
338
 
        return 0.05;
339
 
}
340
 
 
341
 
const char* BASEparticle::getIdentifier() {
342
 
        return "First Order Polynomial Matrix";
 
511
void saveEllipseDataToFile(TransferMatrixGenerator* map, double zIn, double zOut, std::vector<double> initialScales, const char* filename) {
 
512
        std::ostringstream sstream;
 
513
        sstream << lattice_locale << "/" << filename << "/ellipse_data_" << filename << ".txt";
 
514
        std::ofstream stream(sstream.str().c_str(), std::ofstream::out | std::ofstream::trunc);
 
515
        if (!stream.is_open()) {
 
516
                throw InvalidFileException("LieMap::saveStabilityDataToFile",filename);
 
517
        }
 
518
        static const int timesRoundEllipse = 70;
 
519
        
 
520
        
 
521
        for (int i = 0; i < initialScales.size(); i++) {
 
522
                loadBar(i,initialScales.size(),75,getCurrentStream());
 
523
                double* coords = new double[4];
 
524
                coords[0] = initialScales[i];
 
525
                coords[1] = 0.0;
 
526
                coords[2] = initialScales[i];
 
527
                coords[3] = 0.0;
 
528
                double* coordsOut = new double[4];
 
529
                for (int j = 0; j <= timesRoundEllipse; j++) {
 
530
                        map->advanceCell(coords, coordsOut,zIn,zOut);
 
531
                        
 
532
                        if (coordsOut[0] > 1. || coordsOut[2] > 1.) break;
 
533
                        stream << coordsOut[0]*sigma << "\t" << coordsOut[2]*refMomentum << "\n";
 
534
                        
 
535
                        delete[] coords;
 
536
                        coords = coordsOut;
 
537
                        coordsOut = new double[4];
 
538
                }
 
539
                delete[] coordsOut;
 
540
                delete[] coords;
 
541
                stream << "\n";
 
542
        }
 
543
        resetLoadBar();
 
544
        stream.close();
 
545
        std::cerr << "Stability data saved in file " << filename << std::endl;
343
546
}
344
547
 
345
548
 
347
550
 
348
551
Json::Value SetupConfig(int verbose_level, const char* filepath, double zIn) {
349
552
  std::ostringstream sstream;
350
 
  sstream << zIn - 1000.; //start tracking 1m before initialZ value
 
553
  sstream << zIn - 500.; //start tracking 1m before initialZ value
351
554
  std::cerr << "Running with verbose level " << verbose_level << std::endl;
352
555
  Json::Value config(Json::objectValue);
353
556
  config["maximum_module_depth"] = 50;