| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 // *******************************************
26 // *******************************************
27 // |
28 // | G4LowEPComptonModel-- Geant4
29 // | low energy Compton scat
30 // | J. M. C. Brown, Monash Univer
31 // | ## Unpolarised photons
32 // |
33 // |
34 // *******************************************
35 // |
36 // | The following is a Geant4 class to simula
37 // | bound electron Compton scattering. Genera
38 // | based on G4LowEnergyCompton.cc and G4Live
39 // | Algorithms for photon energy, and ejected
40 // | direction taken from:
41 // |
42 // | J. M. C. Brown, M. R. Dimmock, J. E. Gill
43 // | "A low energy bound atomic electron Compt
44 // | for Geant4", NIMB, Vol. 338, 77-88, 2014
45 // |
46 // | The author acknowledges the work of the G
47 // | in developing the following algorithms th
48 // | or adapeted for the present software:
49 // |
50 // | # sampling of photon scattering angle,
51 // | # target element selection in composite
52 // | # target shell selection in element,
53 // | # and sampling of bound electron momentu
54 // |
55 // *******************************************
56 // |
57 // | History:
58 // | --------
59 // |
60 // | Nov. 2011 JMCB - First version
61 // | Feb. 2012 JMCB - Migration to Geant
62 // | Sep. 2012 JMCB - Final fixes for Ge
63 // | Feb. 2013 JMCB - Geant4 9.6 FPE fix
64 // | Jan. 2015 JMCB - Migration to MT (B
65 // | implementation)
66 // | Feb. 2016 JMCB - Geant4 10.2 FPE fi
67 // |
68 // *******************************************
69
70 #include <limits>
71 #include "G4LowEPComptonModel.hh"
72 #include "G4PhysicalConstants.hh"
73 #include "G4SystemOfUnits.hh"
74 #include "G4Electron.hh"
75 #include "G4ParticleChangeForGamma.hh"
76 #include "G4LossTableManager.hh"
77 #include "G4VAtomDeexcitation.hh"
78 #include "G4AtomicShell.hh"
79 #include "G4AutoLock.hh"
80 #include "G4Gamma.hh"
81 #include "G4ShellData.hh"
82 #include "G4DopplerProfile.hh"
83 #include "G4Log.hh"
84 #include "G4Exp.hh"
85
86 //********************************************
87
88 using namespace std;
89 namespace { G4Mutex LowEPComptonModelMutex = G
90
91 G4PhysicsFreeVector* G4LowEPComptonModel::data
92 G4ShellData* G4LowEPComptonModel::shellD
93 G4DopplerProfile* G4LowEPComptonModel::profil
94
95 static const G4double ln10 = G4Log(10.);
96
97 G4LowEPComptonModel::G4LowEPComptonModel(const
98
99 : G4VEmModel(nam),isInitialised(false)
100 {
101 verboseLevel=1 ;
102 // Verbosity scale:
103 // 0 = nothing
104 // 1 = warning for energy non-conservation
105 // 2 = details of energy budget
106 // 3 = calculation of cross sections, file o
107 // 4 = entering in methods
108
109 if( verboseLevel>1 ) {
110 G4cout << "Low energy photon Compton model
111 }
112
113 //Mark this model as "applicable" for atomic
114 SetDeexcitationFlag(true);
115
116 fParticleChange = nullptr;
117 fAtomDeexcitation = nullptr;
118 }
119
120 //********************************************
121
122 G4LowEPComptonModel::~G4LowEPComptonModel()
123 {
124 if(IsMaster()) {
125 delete shellData;
126 shellData = nullptr;
127 delete profileData;
128 profileData = nullptr;
129 }
130 }
131
132 //********************************************
133
134 void G4LowEPComptonModel::Initialise(const G4P
135 const G4D
136 {
137 if (verboseLevel > 1) {
138 G4cout << "Calling G4LowEPComptonModel::In
139 }
140
141 // Initialise element selector
142 if(IsMaster()) {
143
144 // Access to elements
145 const char* path = G4FindDataDir("G4LEDATA
146
147 G4ProductionCutsTable* theCoupleTable =
148 G4ProductionCutsTable::GetProductionCuts
149 G4int numOfCouples = (G4int)theCoupleTable
150
151 for(G4int i=0; i<numOfCouples; ++i) {
152 const G4Material* material =
153 theCoupleTable->GetMaterialCutsCouple(
154 const G4ElementVector* theElementVector
155 std::size_t nelm = material->GetNumberOf
156
157 for (std::size_t j=0; j<nelm; ++j) {
158 G4int Z = G4lrint((*theElementVector)[
159 if(Z < 1) { Z = 1; }
160 else if(Z > maxZ){ Z = maxZ; }
161
162 if( (!data[Z]) ) { ReadData(Z, path);
163 }
164 }
165
166 // For Doppler broadening
167 if(!shellData) {
168 shellData = new G4ShellData();
169 shellData->SetOccupancyData();
170 G4String file = "/doppler/shell-doppler"
171 shellData->LoadData(file);
172 }
173 if(!profileData) { profileData = new G4Dop
174
175 InitialiseElementSelectors(particle, cuts)
176 }
177
178 if (verboseLevel > 2) {
179 G4cout << "Loaded cross section files" <<
180 }
181
182 if( verboseLevel>1 ) {
183 G4cout << "G4LowEPComptonModel is initiali
184 << "Energy range: "
185 << LowEnergyLimit() / eV << " eV -
186 << HighEnergyLimit() / GeV << " GeV
187 << G4endl;
188 }
189
190 if(isInitialised) { return; }
191
192 fParticleChange = GetParticleChangeForGamma(
193 fAtomDeexcitation = G4LossTableManager::Ins
194 isInitialised = true;
195 }
196
197 //********************************************
198
199 void G4LowEPComptonModel::InitialiseLocal(cons
200
201 {
202 SetElementSelectors(masterModel->GetElementS
203 }
204
205 //********************************************
206
207 void G4LowEPComptonModel::ReadData(std::size_t
208 {
209 if (verboseLevel > 1)
210 {
211 G4cout << "G4LowEPComptonModel::ReadData()
212 << G4endl;
213 }
214 if(data[Z]) { return; }
215 const char* datadir = path;
216 if(!datadir)
217 {
218 datadir = G4FindDataDir("G4LEDATA");
219 if(!datadir)
220 {
221 G4Exception("G4LowEPComptonModel::ReadDa
222 "em0006",FatalException,
223 "Environment variable G4LEDA
224 return;
225 }
226 }
227
228 data[Z] = new G4PhysicsFreeVector();
229
230 std::ostringstream ost;
231 ost << datadir << "/livermore/comp/ce-cs-" <
232 std::ifstream fin(ost.str().c_str());
233
234 if( !fin.is_open())
235 {
236 G4ExceptionDescription ed;
237 ed << "G4LowEPComptonModel data file <"
238 << "> is not opened!" << G4endl;
239 G4Exception("G4LowEPComptonModel::ReadData
240 "em0003",FatalException,
241 ed,"G4LEDATA version should be
242 return;
243 } else {
244 if(verboseLevel > 3) {
245 G4cout << "File " << ost.str()
246 << " is opened by G4LowEPCompto
247 }
248 data[Z]->Retrieve(fin, true);
249 data[Z]->ScaleVector(MeV, MeV*barn);
250 }
251 fin.close();
252 }
253
254 //********************************************
255
256
257 G4double
258 G4LowEPComptonModel::ComputeCrossSectionPerAto
259
260
261
262 {
263 if (verboseLevel > 3) {
264 G4cout << "G4LowEPComptonModel::ComputeCro
265 << G4endl;
266 }
267 G4double cs = 0.0;
268
269 if (GammaEnergy < LowEnergyLimit()) { return
270
271 G4int intZ = G4lrint(Z);
272 if(intZ < 1 || intZ > maxZ) { return cs; }
273
274 G4PhysicsFreeVector* pv = data[intZ];
275
276 // if element was not initialised
277 // do initialisation safely for MT mode
278 if(!pv)
279 {
280 InitialiseForElement(0, intZ);
281 pv = data[intZ];
282 if(!pv) { return cs; }
283 }
284
285 G4int n = G4int(pv->GetVectorLength() - 1);
286 G4double e1 = pv->Energy(0);
287 G4double e2 = pv->Energy(n);
288
289 if(GammaEnergy <= e1) { cs = GammaEnerg
290 else if(GammaEnergy <= e2) { cs = pv->Value(
291 else if(GammaEnergy > e2) { cs = pv->Value(
292
293 return cs;
294 }
295
296 //********************************************
297
298 void G4LowEPComptonModel::SampleSecondaries(st
299
300
301
302 {
303
304 // The scattered gamma energy is sampled acc
305 // then accepted or rejected depending on th
306 // by factor from Klein - Nishina formula.
307 // Expression of the angular distribution as
308 // angular and energy distribution and Scatt
309 // D. E. Cullen "A simple model of photon tr
310 // Phys. Res. B 101 (1995). Method of sampli
311 // data are interpolated while in the articl
312 // Reference to the article is from J. Stepa
313 // and Electron Interaction Data for GEANT i
314 // TeV (draft).
315 // The random number techniques of Butcher &
316 // (Nucl Phys 20(1960),15).
317
318 G4double photonEnergy0 = aDynamicGamma->GetK
319
320 if (verboseLevel > 3) {
321 G4cout << "G4LowEPComptonModel::SampleSeco
322 << photonEnergy0/MeV << " in " << c
323 << G4endl;
324 }
325 // do nothing below the threshold
326 // should never get here because the XS is z
327 if (photonEnergy0 < LowEnergyLimit())
328 return ;
329
330 G4double e0m = photonEnergy0 / electron_mass
331 G4ParticleMomentum photonDirection0 = aDynam
332
333 // Select randomly one element in the curren
334 const G4ParticleDefinition* particle = aDyn
335 const G4Element* elm = SelectRandomAtom(coup
336 G4int Z = (G4int)elm->GetZ();
337
338 G4double LowEPCepsilon0 = 1. / (1. + 2. * e0
339 G4double LowEPCepsilon0Sq = LowEPCepsilon0 *
340 G4double alpha1 = -std::log(LowEPCepsilon0);
341 G4double alpha2 = 0.5 * (1. - LowEPCepsilon0
342
343 G4double wlPhoton = h_Planck*c_light/photonE
344
345 // Sample the energy of the scattered photon
346 G4double LowEPCepsilon;
347 G4double LowEPCepsilonSq;
348 G4double oneCosT;
349 G4double sinT2;
350 G4double gReject;
351
352 if (verboseLevel > 3) {
353 G4cout << "Started loop to sample gamma en
354 }
355
356 do
357 {
358 if ( alpha1/(alpha1+alpha2) > G4UniformR
359 {
360 LowEPCepsilon = G4Exp(-alpha1 * G4Un
361 LowEPCepsilonSq = LowEPCepsilon * Lo
362 }
363 else
364 {
365 LowEPCepsilonSq = LowEPCepsilon0Sq +
366 LowEPCepsilon = std::sqrt(LowEPCepsi
367 }
368
369 oneCosT = (1. - LowEPCepsilon) / ( LowEP
370 sinT2 = oneCosT * (2. - oneCosT);
371 G4double x = std::sqrt(oneCosT/2.) / (wl
372 G4double scatteringFunction = ComputeSca
373 gReject = (1. - LowEPCepsilon * sinT2 /
374
375 } while(gReject < G4UniformRand()*Z);
376
377 G4double cosTheta = 1. - oneCosT;
378 G4double sinTheta = std::sqrt(sinT2);
379 G4double phi = twopi * G4UniformRand();
380 G4double dirx = sinTheta * std::cos(phi);
381 G4double diry = sinTheta * std::sin(phi);
382 G4double dirz = cosTheta ;
383
384 // Scatter photon energy and Compton electro
385 // J. M. C. Brown, M. R. Dimmock, J. E. Gill
386 // "A low energy bound atomic electron Compt
387 // NIMB, Vol. 338, 77-88, 2014.
388
389 // Set constants and initialize scattering p
390 const G4double vel_c = c_light / (m/s);
391 const G4double momentum_au_to_nat = halfpi*
392 const G4double e_mass_kg = electron_mass_c2
393
394 const G4int maxDopplerIterations = 1000;
395 G4double bindingE = 0.;
396 G4double pEIncident = photonEnergy0 ;
397 G4double pERecoil = -1.;
398 G4double eERecoil = -1.;
399 G4double e_alpha =0.;
400 G4double e_beta = 0.;
401
402 G4double CE_emission_flag = 0.;
403 G4double ePAU = -1;
404 G4int shellIdx = 0;
405 G4double u_temp = 0;
406 G4double cosPhiE =0;
407 G4double sinThetaE =0;
408 G4double cosThetaE =0;
409 G4int iteration = 0;
410
411 if (verboseLevel > 3) {
412 G4cout << "Started loop to sample photon e
413 }
414
415 do{
416 // *************************************
417 // | Determine scatter photon energy
418 // *************************************
419 do
420 {
421 iteration++;
422 // *****************************************
423 // | Sample bound electron information
424 // *****************************************
425
426 // Select shell based on shell occupancy
427 shellIdx = shellData->SelectRandomShell(Z);
428 bindingE = shellData->BindingEnergy(Z,shellI
429
430 // Randomly sample bound electron momentum (
431 ePAU = profileData->RandomSelectMomentum(Z,s
432
433 // Convert to SI units
434 G4double ePSI = ePAU * momentum_au_to_nat;
435
436 //Calculate bound electron velocity and norm
437 u_temp = sqrt( ((ePSI*ePSI)*(vel_c*vel_c)) /
438
439 // Sample incident electron direction, amorp
440 e_alpha = pi*G4UniformRand();
441 e_beta = twopi*G4UniformRand();
442
443 // Total energy of system
444
445 G4double eEIncident = electron_mass_c2 / sqr
446 G4double systemE = eEIncident + pEIncident;
447 G4double gamma_temp = 1.0 / sqrt( 1 - (u_tem
448 G4double numerator = gamma_temp*electron_mas
449 G4double subdenom1 = u_temp*cosTheta*std::c
450 G4double subdenom2 = u_temp*sinTheta*std::si
451 G4double denominator = (1.0 - cosTheta) + (
452 pERecoil = (numerator/denominator);
453 eERecoil = systemE - pERecoil;
454 CE_emission_flag = pEIncident - pERecoil;
455 } while ( (iteration <= maxDopplerIterat
456
457 // End of recalculation of photon energy w
458
459 // ***************************************
460 // | Determine ejected Compton electro
461 // ***************************************
462
463 // Calculate velocity of ejected Compton e
464
465 G4double a_temp = eERecoil / electron_mass
466 G4double u_p_temp = sqrt(1 - (1 / (a_temp*
467
468 // Coefficients and terms from simulatenou
469 G4double sinAlpha = std::sin(e_alpha);
470 G4double cosAlpha = std::cos(e_alpha);
471 G4double sinBeta = std::sin(e_beta);
472 G4double cosBeta = std::cos(e_beta);
473
474 G4double gamma = 1.0 / sqrt(1 - (u_temp*u_
475 G4double gamma_p = 1.0 / sqrt(1 - (u_p_tem
476
477 G4double var_A = pERecoil*u_p_temp*sinThet
478 G4double var_B = u_p_temp* (pERecoil*cosTh
479 G4double var_C = (pERecoil-pEIncident) - (
480
481 G4double var_D1 = gamma*electron_mass_c2*p
482 G4double var_D2 = (1 - (u_temp*cosTheta*co
483 G4double var_D3 = ((electron_mass_c2*elect
484 G4double var_D = var_D1*var_D2 + var_D3;
485
486 G4double var_E1 = ((gamma*gamma_p)*(electr
487 G4double var_E2 = gamma_p*electron_mass_c2
488 G4double var_E = var_E1 - var_E2;
489
490 G4double var_F1 = ((gamma*gamma_p)*(electr
491 G4double var_F2 = (gamma_p*electron_mass_c
492 G4double var_F = var_F1 - var_F2;
493
494 G4double var_G = (gamma*gamma_p)*(electron
495
496 // Two equations form a quadratic form of
497 // Coefficents and solution to quadratic
498 G4double var_W1 = (var_F*var_B - var_E*var
499 G4double var_W2 = (var_G*var_G)*(var_A*var
500 G4double var_W = var_W1 + var_W2;
501
502 G4double var_Y = 2.0*(((var_A*var_D-var_F*
503
504 G4double var_Z1 = (var_A*var_D - var_F*var
505 G4double var_Z2 = (var_G*var_G)*(var_C*var
506 G4double var_Z = var_Z1 + var_Z2;
507 G4double diff1 = var_Y*var_Y;
508 G4double diff2 = 4*var_W*var_Z;
509 G4double diff = diff1 - diff2;
510
511 // Check if diff is less than zero, if so
512 //Determine number of digits (in decimal b
513 G4double g4d_order = G4double(numeric_limi
514 G4double g4d_limit = std::pow(10.,-g4d_ord
515 //Confirm that diff less than zero is due
516 //than 10^(-g4d_order), then set diff to z
517
518 if ((diff < 0.0) && (abs(diff / diff1) < g
519 {
520 diff = 0.0;
521 }
522
523
524 // Plus and minus of quadratic
525 G4double X_p = (-var_Y + sqrt (diff))/(2*v
526 G4double X_m = (-var_Y - sqrt (diff))/(2*v
527
528 // Floating point precision protection
529 // Check if X_p and X_m are greater than o
530 // Issue due to propagation of FPE and onl
531 if(X_p >1){X_p=1;} if(X_p<-1){X_p=-1;}
532 if(X_m >1){X_m=1;} if(X_m<-1){X_m=-1;}
533 // End of FP protection
534
535 G4double ThetaE = 0.;
536 // Randomly sample one of the two possible
537 G4double sol_select = G4UniformRand();
538
539 if (sol_select < 0.5)
540 {
541 ThetaE = std::acos(X_p);
542 }
543 if (sol_select > 0.5)
544 {
545 ThetaE = std::acos(X_m);
546 }
547 cosThetaE = std::cos(ThetaE);
548 sinThetaE = std::sin(ThetaE);
549 G4double Theta = std::acos(cosTheta);
550
551 //Calculate electron Phi
552 G4double iSinThetaE = std::sqrt(1+std::tan
553 G4double iSinTheta = std::sqrt(1+std::tan(
554 G4double ivar_A = iSinTheta/ (pERecoil*u_p
555 // Trigs
556 cosPhiE = (var_C - var_B*cosThetaE)*(ivar_
557
558 // End of calculation of ejection Compton
559 //Fix for floating point errors
560
561 } while ( (iteration <= maxDopplerIterations
562
563 // Revert to original if maximum number of i
564 if (iteration >= maxDopplerIterations)
565 {
566 pERecoil = photonEnergy0 ;
567 bindingE = 0.;
568 dirx=0.0;
569 diry=0.0;
570 dirz=1.0;
571 }
572
573 // Set "scattered" photon direction and ener
574 G4ThreeVector photonDirection1(dirx,diry,dir
575 photonDirection1.rotateUz(photonDirection0);
576 fParticleChange->ProposeMomentumDirection(ph
577
578 if (pERecoil > 0.)
579 {
580 fParticleChange->SetProposedKineticEnergy
581
582 // Set ejected Compton electron direction
583 G4double PhiE = std::acos(cosPhiE);
584 G4double eDirX = sinThetaE * std::cos(phi
585 G4double eDirY = sinThetaE * std::sin(phi
586 G4double eDirZ = cosThetaE;
587
588 G4double eKineticEnergy = pEIncident - pE
589
590 G4ThreeVector eDirection(eDirX,eDirY,eDir
591 eDirection.rotateUz(photonDirection0);
592 G4DynamicParticle* dp = new G4DynamicPart
593
594 fvect->push_back(dp);
595
596 }
597 else
598 {
599 fParticleChange->SetProposedKineticEnerg
600 fParticleChange->ProposeTrackStatus(fSto
601 }
602
603 // sample deexcitation
604 //
605 if (verboseLevel > 3) {
606 G4cout << "Started atomic de-excitation "
607 }
608
609 if(fAtomDeexcitation && iteration < maxDoppl
610 G4int index = couple->GetIndex();
611 if(fAtomDeexcitation->CheckDeexcitationAct
612 std::size_t nbefore = fvect->size();
613 G4AtomicShellEnumerator as = G4AtomicShe
614 const G4AtomicShell* shell = fAtomDeexci
615 fAtomDeexcitation->GenerateParticles(fve
616 std::size_t nafter = fvect->size();
617 if(nafter > nbefore) {
618 for (std::size_t i=nbefore; i<nafter;
619 //Check if there is enough residual
620 if (bindingE >= ((*fvect)[i])->GetKi
621 {
622 //Ok, this is a valid secondary:
623 bindingE -= ((*fvect)[i])->GetKin
624 }
625 else
626 {
627 //Invalid secondary: not enough e
628 //Keep its energy in the local de
629 delete (*fvect)[i];
630 (*fvect)[i]=nullptr;
631 }
632 }
633 }
634 }
635 }
636
637 //This should never happen
638 if(bindingE < 0.0)
639 G4Exception("G4LowEPComptonModel::SampleS
640 "em2051",FatalException,"Nega
641
642 fParticleChange->ProposeLocalEnergyDeposit(b
643 }
644
645 //********************************************
646
647 G4double
648 G4LowEPComptonModel::ComputeScatteringFunction
649 {
650 G4double value = Z;
651 if (x <= ScatFuncFitParam[Z][2]) {
652
653 G4double lgq = G4Log(x)/ln10;
654
655 if (lgq < ScatFuncFitParam[Z][1]) {
656 value = ScatFuncFitParam[Z][3] + lgq*Sca
657 } else {
658 value = ScatFuncFitParam[Z][5] + lgq*Sca
659 lgq*lgq*ScatFuncFitParam[Z][7] + lgq*l
660 }
661 value = G4Exp(value*ln10);
662 }
663 return value;
664 }
665
666
667 //********************************************
668
669 void
670 G4LowEPComptonModel::InitialiseForElement(cons
671
672 {
673 G4AutoLock l(&LowEPComptonModelMutex);
674 if(!data[Z]) { ReadData(Z); }
675 l.unlock();
676 }
677
678 //********************************************
679
680 //Fitting data to compute scattering function
681
682 const G4double G4LowEPComptonModel::ScatFuncFi
683 { 0, 0., 0., 0., 0.,
684 { 1, 6.673, 1.49968E+08, -14.352, 1.999, -143
685 { 2, 6.500, 2.50035E+08, -14.215, 1.970, -53.
686 { 3, 6.551, 3.99945E+08, -13.555, 1.993, -62.
687 { 4, 6.500, 5.00035E+08, -13.746, 1.998, -127
688 { 5, 6.500, 5.99791E+08, -13.800, 1.998, -131
689 { 6, 6.708, 6.99842E+08, -13.885, 1.999, -128
690 { 7, 6.685, 7.99834E+08, -13.885, 2.000, -131
691 { 8, 6.669, 7.99834E+08, -13.962, 2.001, -128
692 { 9, 6.711, 7.99834E+08, -13.999, 2.000, -122
693 { 10, 6.702, 7.99834E+08, -14.044, 1.999, -110
694 { 11, 6.425, 1.00000E+09, -13.423, 1.993, -41.
695 { 12, 6.542, 1.00000E+09, -13.389, 1.997, -53.
696 { 13, 6.570, 1.49968E+09, -13.401, 1.997, -66.
697 { 14, 6.364, 1.49968E+09, -13.452, 1.999, -78.
698 { 15, 6.500, 1.49968E+09, -13.488, 1.998, -85.
699 { 16, 6.500, 1.49968E+09, -13.532, 1.998, -93.
700 { 17, 6.500, 1.49968E+09, -13.584, 2.000, -98.
701 { 18, 6.500, 1.49968E+09, -13.618, 1.999, -100
702 { 19, 6.500, 1.99986E+09, -13.185, 1.992, -53.
703 { 20, 6.490, 1.99986E+09, -13.123, 1.993, -52.
704 { 21, 6.498, 1.99986E+09, -13.157, 1.994, -55.
705 { 22, 6.495, 1.99986E+09, -13.183, 1.994, -57.
706 { 23, 6.487, 1.99986E+09, -13.216, 1.995, -58.
707 { 24, 6.500, 1.99986E+09, -13.330, 1.997, -62.
708 { 25, 6.488, 1.99986E+09, -13.277, 1.997, -58.
709 { 26, 6.500, 5.00035E+09, -13.292, 1.997, -61.
710 { 27, 6.500, 5.00035E+09, -13.321, 1.998, -61.
711 { 28, 6.500, 5.00035E+09, -13.340, 1.998, -62.
712 { 29, 6.500, 5.00035E+09, -13.439, 1.998, -67.
713 { 30, 6.500, 5.00035E+09, -13.383, 1.999, -62.
714 { 31, 6.500, 5.00035E+09, -13.349, 1.997, -61.
715 { 32, 6.500, 5.00035E+09, -13.373, 1.999, -63.
716 { 33, 6.500, 5.00035E+09, -13.395, 1.999, -65.
717 { 34, 6.500, 5.00035E+09, -13.417, 1.999, -69.
718 { 35, 6.500, 5.00035E+09, -13.442, 2.000, -72.
719 { 36, 6.500, 5.00035E+09, -13.451, 1.998, -74.
720 { 37, 6.500, 5.00035E+09, -13.082, 1.991, -46.
721 { 38, 6.465, 5.00035E+09, -13.022, 1.993, -41.
722 { 39, 6.492, 5.00035E+09, -13.043, 1.994, -44.
723 { 40, 6.499, 5.00035E+09, -13.064, 1.994, -47.
724 { 41, 6.384, 5.00035E+09, -13.156, 1.996, -53.
725 { 42, 6.500, 5.00035E+09, -13.176, 1.996, -54.
726 { 43, 6.500, 5.00035E+09, -13.133, 1.997, -51.
727 { 44, 6.500, 5.00035E+09, -13.220, 1.996, -58.
728 { 45, 6.500, 5.00035E+09, -13.246, 1.998, -59.
729 { 46, 6.500, 5.00035E+09, -13.407, 1.999, -72.
730 { 47, 6.500, 5.00035E+09, -13.277, 1.998, -60.
731 { 48, 6.500, 5.00035E+09, -13.222, 1.998, -56.
732 { 49, 6.500, 5.00035E+09, -13.199, 1.997, -56.
733 { 50, 6.500, 5.00035E+09, -13.215, 1.998, -58.
734 { 51, 6.500, 5.00035E+09, -13.230, 1.998, -60.
735 { 52, 6.500, 7.99834E+09, -13.246, 1.998, -63.
736 { 53, 6.500, 5.00035E+09, -13.262, 1.998, -66.
737 { 54, 6.500, 7.99834E+09, -13.279, 1.998, -67.
738 { 55, 6.500, 5.00035E+09, -12.951, 1.990, -45.
739 { 56, 6.425, 5.00035E+09, -12.882, 1.992, -39.
740 { 57, 6.466, 2.82488E+09, -12.903, 1.992, -38.
741 { 58, 6.451, 5.00035E+09, -12.915, 1.993, -41.
742 { 59, 6.434, 5.00035E+09, -12.914, 1.993, -40.
743 { 60, 6.444, 5.00035E+09, -12.922, 1.992, -39.
744 { 61, 6.414, 7.99834E+09, -12.930, 1.993, -42.
745 { 62, 6.420, 7.99834E+09, -12.938, 1.992, -42.
746 { 63, 6.416, 7.99834E+09, -12.946, 1.993, -42.
747 { 64, 6.443, 7.99834E+09, -12.963, 1.993, -43.
748 { 65, 6.449, 7.99834E+09, -12.973, 1.993, -43.
749 { 66, 6.419, 7.99834E+09, -12.966, 1.993, -42.
750 { 67, 6.406, 7.99834E+09, -12.976, 1.993, -42.
751 { 68, 6.424, 7.99834E+09, -12.986, 1.993, -41.
752 { 69, 6.417, 7.99834E+09, -12.989, 1.993, -42.
753 { 70, 6.405, 7.99834E+09, -13.000, 1.994, -42.
754 { 71, 6.449, 7.99834E+09, -13.015, 1.994, -42.
755 { 72, 6.465, 7.99834E+09, -13.030, 1.994, -43.
756 { 73, 6.447, 7.99834E+09, -13.048, 1.996, -44.
757 { 74, 6.452, 7.99834E+09, -13.073, 1.997, -45.
758 { 75, 6.432, 7.99834E+09, -13.082, 1.997, -46.
759 { 76, 6.439, 7.99834E+09, -13.100, 1.997, -47.
760 { 77, 6.432, 7.99834E+09, -13.110, 1.997, -48.
761 { 78, 6.500, 7.99834E+09, -13.185, 1.997, -53.
762 { 79, 6.500, 7.99834E+09, -13.200, 1.997, -53.
763 { 80, 6.500, 7.99834E+09, -13.156, 1.998, -49.
764 { 81, 6.500, 7.99834E+09, -13.128, 1.997, -49.
765 { 82, 6.500, 7.99834E+09, -13.134, 1.997, -51.
766 { 83, 6.500, 7.99834E+09, -13.148, 1.998, -52.
767 { 84, 6.500, 7.99834E+09, -13.161, 1.998, -54.
768 { 85, 6.500, 7.99834E+09, -13.175, 1.998, -56.
769 { 86, 6.500, 7.99834E+09, -13.189, 1.998, -57.
770 { 87, 6.500, 7.99834E+09, -12.885, 1.990, -39.
771 { 88, 6.417, 7.99834E+09, -12.816, 1.991, -34.
772 { 89, 6.442, 7.99834E+09, -12.831, 1.992, -36.
773 { 90, 6.463, 7.99834E+09, -12.850, 1.993, -37.
774 { 91, 6.447, 7.99834E+09, -12.852, 1.993, -37.
775 { 92, 6.439, 7.99834E+09, -12.858, 1.993, -37.
776 { 93, 6.437, 1.00000E+10, -12.866, 1.993, -39.
777 { 94, 6.432, 7.99834E+09, -12.862, 1.993, -37.
778 { 95, 6.435, 7.99834E+09, -12.869, 1.993, -37.
779 { 96, 6.449, 1.00000E+10, -12.886, 1.993, -39.
780 { 97, 6.446, 1.00000E+10, -12.892, 1.993, -40.
781 { 98, 6.421, 1.00000E+10, -12.887, 1.993, -39.
782 { 99, 6.414, 1.00000E+10, -12.894, 1.993, -39.
783 {100, 6.412, 1.00000E+10, -12.900, 1.993, -39.
784 };
785
786
787
788