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