Changeset 2127
- Timestamp:
- 06/25/08 11:17:40 (2 months ago)
- Location:
- trunk/source
- Files:
-
- 25 modified
-
abundances.cpp (modified) (1 diff)
-
assert_results.cpp (modified) (1 diff)
-
cont_createmesh.cpp (modified) (1 diff)
-
cont_createpointers.cpp (modified) (3 diffs)
-
conv_init_solution.cpp (modified) (2 diffs)
-
grid_do.cpp (modified) (1 diff)
-
helike_cs.cpp (modified) (12 diffs)
-
helike_cs.h (modified) (1 diff)
-
helike_einsta.cpp (modified) (20 diffs)
-
helike_recom.cpp (modified) (10 diffs)
-
helike_recom.h (modified) (1 diff)
-
hydrocollid.cpp (modified) (2 diffs)
-
hydroeinsta.cpp (modified) (5 diffs)
-
hydrogenic.h (modified) (1 diff)
-
hydro_recom.cpp (modified) (4 diffs)
-
hydro_vs_rates.cpp (modified) (9 diffs)
-
hydro_vs_rates.h (modified) (2 diffs)
-
init_coreload.cpp (modified) (1 diff)
-
iso.h (modified) (3 diffs)
-
iso_collide.cpp (modified) (7 diffs)
-
iso_error.cpp (modified) (3 diffs)
-
iso_radiative_recomb.cpp (modified) (6 diffs)
-
iso_solve.cpp (modified) (2 diffs)
-
nemala2.cpp (modified) (2 diffs)
-
parse_atom_iso.cpp (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/source/abundances.cpp
r1771 r2127 269 269 270 270 /* print warning if temperature set below default but C > O */ 271 if( dense. gas_phase[ipCARBON]/SDIV( dense.gas_phase[ipOXYGEN]) >= 1. )271 if( dense.lgElmtOn[ipOXYGEN] && dense.gas_phase[ipCARBON]/SDIV( dense.gas_phase[ipOXYGEN]) >= 1. ) 272 272 { 273 273 fprintf( ioQQQ, "\n >>> \n" -
trunk/source/assert_results.cpp
r2012 r2127 1943 1943 1944 1944 int iCase = 1; 1945 if( strcmp( "Ca A", chAssertLineLabel[i] ) == 0 ) 1946 iCase = 0; 1947 else if( strcmp( "Ca B", chAssertLineLabel[i] ) == 0 ) 1948 iCase = 1; 1949 else if( strcmp( "+Col", chAssertLineLabel[i] ) == 0 ) 1950 iCase = 1; 1951 else 1952 TotalInsanity(); 1953 1945 1954 RelError[i] = 0.; 1946 1955 long nHighestPrinted = StatesElem[nISOCaseB][nelemCaseB][iso.numPrintLevels[nISOCaseB][nelemCaseB]-1].n; -
trunk/source/cont_createmesh.cpp
r2120 r2127 462 462 * 463 463 * rfield.fine_opac_nelem is the most massive (hence sharpest line) 464 * we will worry about. By default this is ir ion but can be changed464 * we will worry about. By default this is iron but can be changed 465 465 * with SET FINE CONTINUUM command 466 466 * -
trunk/source/cont_createpointers.cpp
r1926 r2127 632 632 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO ) 633 633 { 634 /* do remaining part of the he iso sequence*/634 /* do remaining part of the iso sequences */ 635 635 for( nelem=2; nelem < LIMELM; nelem++ ) 636 636 { … … 639 639 /* generate label for this ion */ 640 640 sprintf( chLab, "%2s%2ld",elementnames.chElementSym[nelem], nelem+1-ipISO ); 641 /* Lya itself is the only transition below n=3 - we explicitly do not 642 * want to generate pointers for 2s-1s or 2p-2s */ 643 /* array index for continuum edges for levels in He-like ions */ 641 /* array index for continuum edges */ 644 642 iso.ipIsoLevNIonCon[ipISO][nelem][0] = 645 643 ipContEnergy(iso.xIsoLevNIonRyd[ipISO][nelem][0],chLab); … … 647 645 for( ipHi=1; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ ) 648 646 { 649 /* array index for continuum edges for levels in He-like ions*/647 /* array index for continuum edges */ 650 648 iso.ipIsoLevNIonCon[ipISO][nelem][ipHi] = ipContEnergy(iso.xIsoLevNIonRyd[ipISO][nelem][ipHi],chLab); 651 649 652 /* define all he-likeline pointers */650 /* define all line pointers */ 653 651 for( ipLo=0; ipLo < ipHi; ipLo++ ) 654 652 { -
trunk/source/conv_init_solution.cpp
r2119 r2127 20 20 #include "pressure.h" 21 21 #include "conv.h" 22 #include "hmi.h" 22 23 23 24 /* derivative of net cooling wrt temperature to check on sign oscillations */ … … 383 384 TeFirst = MAX2( 4000. , thermal.ConstTemp ); 384 385 /* true allow ionization stage trimming, false blocks it */ 386 387 /* override this if molecule deliberately disabled */ 388 if( hmi.lgNoH2Mole ) 389 TeFirst = thermal.ConstTemp; 385 390 } 386 391 else if( thermal.lgTeHigh ) -
trunk/source/grid_do.cpp
r1894 r2127 310 310 311 311 fprintf( ioQQQ, "\n %s\n", chLine ); 312 fprintf( ioQQQ, " Initial increment is %6.3f, the limits are%10.2e to %10.2e\n",312 fprintf( ioQQQ, " Initial increment is %6.3f, the limits are%10.2e to %10.2e\n", 313 313 optimize.vincr[i], optimize.varang[i][0], optimize.varang[i][1] ); 314 314 } -
trunk/source/helike_cs.cpp
r1960 r2127 225 225 } 226 226 227 /* set 15% errors on all collision rates. */228 iso_put_error(ipHE_LIKE, nelem, ipHi, ipLo, IPCOLLIS, 0.15f );229 230 227 ASSERT( cs >= 0.f ); 231 228 … … 451 448 /* >>refer He cs Pengelly, R.M., & Seaton, M.J., 1964, MNRAS, 127, 165 */ 452 449 /* statistical weights included */ 453 cs = (realnum)CS_l_mixing_PS64( ipHE_LIKE, nelem, ipLo, ipHi, Collider); 450 cs = (realnum)CS_l_mixing_PS64( 451 nelem, 452 StatesElem[ipHE_LIKE][nelem][ipLo].lifetime, 453 nelem+1.-ipHE_LIKE, 454 StatesElem[ipHE_LIKE][nelem][ipLo].n, 455 StatesElem[ipHE_LIKE][nelem][ipLo].l, 456 StatesElem[ipHE_LIKE][nelem][ipHi].g, 457 Collider); 454 458 } 455 459 else … … 589 593 ASSERT( cs >= 0.f ); 590 594 591 /*iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPCOLLIS,-1);*/592 593 595 return(cs); 594 596 … … 797 799 /* >>refer He cs Pengelly, R.M., & Seaton, M.J., 1964, MNRAS, 127, 165 */ 798 800 /* statistical weights included */ 799 cs = (realnum)CS_l_mixing_PS64( ipHE_LIKE, nelem, ipLo, ipHi, Collider); 801 cs = (realnum)CS_l_mixing_PS64( 802 nelem, 803 StatesElem[ipHE_LIKE][nelem][ipLo].lifetime, 804 nelem+1.-ipHE_LIKE, 805 StatesElem[ipHE_LIKE][nelem][ipLo].n, 806 StatesElem[ipHE_LIKE][nelem][ipLo].l, 807 StatesElem[ipHE_LIKE][nelem][ipHi].g, 808 Collider); 800 809 } 801 810 else … … 857 866 858 867 ASSERT( cs >= 0.f ); 859 860 /*iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPCOLLIS,-1);*/861 868 862 869 return(cs); … … 903 910 coll_str += qg32( 1.0, 10.0, S62_Therm_ave_coll_str); 904 911 ASSERT( coll_str > 0. ); 905 906 /*iso_put_error(ipISO,nelem,ipHi,ipLo,ipCollis,-1);*/907 912 908 913 return coll_str; … … 1024 1029 /*CS_l_mixing_PS64 - find rate for l-mixing collisions by protons, for neutrals */ 1025 1030 double CS_l_mixing_PS64( 1026 long ipISO, 1027 long nelem /* the chemical element, 1 for He */, 1028 long ipLo /* lower level, 0 for ground */, 1029 long ipHi /* upper level, 0 for ground */, 1031 long nelem, /* the chemical element, 1 for He */ 1032 double tau, /* the radiative lifetime of the initial level. */ 1033 double target_charge, 1034 long n, 1035 long l, 1036 double gHi, 1030 1037 long Collider) 1031 1038 { … … 1037 1044 bestfactor, factorpart, 1038 1045 reduced_mass, reduced_mass_2_emass, 1039 rate , tau;1046 rate; 1040 1047 const double ChargIncoming = ColliderCharge[Collider]; 1041 1048 1042 1049 DEBUG_ENTRY( "CS_l_mixing_PS64()" ); 1043 1050 1044 ASSERT( ipHi > ipLo );1045 1051 /* In this routine, two different cutoff radii are calculated, and as per PS64, 1046 1052 * the least of these is selected. We take the least positive result. */ … … 1051 1057 /* this mass always appears relative to the electron mass, so define it that way */ 1052 1058 reduced_mass_2_emass = reduced_mass / ELECTRON_MASS; 1053 1054 /* This is the lifetime of ipLo. */1055 tau = StatesElem[ipISO][nelem][ipLo].lifetime;1056 1057 if( ipISO == ipH_LIKE && Transitions[ipISO][nelem][ipLo][ipH1s].ipCont > 0 )1058 {1059 tau = (1./tau) + Transitions[ipISO][nelem][ipLo][ipH1s].Emis->Aul;1060 tau = 1./tau;1061 }1062 1059 1063 1060 /* equation 46 of PS64 */ … … 1074 1071 1075 1072 /* this is equation 44 of PS64 */ 1076 Dnl = POW2( ChargIncoming / (double)(nelem + 1. - ipISO)) * 6. * POW2( (double)StatesElem[ipISO][nelem][ipLo].n ) * 1077 ( POW2((double)StatesElem[ipISO][nelem][ipLo].n) - 1078 POW2((double)StatesElem[ipISO][nelem][ipLo].l) - StatesElem[ipISO][nelem][ipLo].l - 1); 1073 Dnl = POW2( ChargIncoming / target_charge) * 6. * POW2( (double)n) * 1074 ( POW2((double)n) - POW2((double)l) - l - 1); 1079 1075 1080 1076 ASSERT( Dnl > 0. ); … … 1110 1106 /* this is the TOTAL rate from nl to nl+/-1. So assume we can 1111 1107 * divide by two to get the rate nl to either nl+1 or nl-1. */ 1112 if( StatesElem[ipISO][nelem][ipLo].l > 0 )1108 if( l > 0 ) 1113 1109 rate /=2; 1114 1110 … … 1117 1113 * for electron colliders and is off by reduced_mass_2_emass^-1.5 */ 1118 1114 cs = rate / ( COLL_CONST * pow(reduced_mass_2_emass, -1.5) ) * 1119 phycon.sqrte * StatesElem[ipISO][nelem][ipHi].g;1115 phycon.sqrte * gHi; 1120 1116 1121 1117 ASSERT( cs > 0. ); 1122 1123 /*iso_put_error(ipISO,nelem,ipHi,ipLo,ipCollis,-1);*/1124 1118 1125 1119 return cs; -
trunk/source/helike_cs.h
r1732 r2127 76 76 77 77 /**CS_l_mixing_PS64 Collision treatment based on Pengelly and Seaton 1964 78 \param ipISO 79 \param nelem the chemical element, 1 for He 80 \param ipLo lower level, 0 for ground 81 \param ipHi upper level, 0 for ground 78 \param nelem, the chemical element, 1 for He 79 \param tau, 80 \param target_charge, 81 \param n, 82 \param l, 83 \param gHi, 82 84 \param Collider 83 85 */ 84 86 double CS_l_mixing_PS64( 85 long int ipISO,86 87 long int nelem, 87 long int ipLo , 88 long int ipHi , 89 long int Collider ); 88 double tau, 89 double target_charge, 90 long int n, 91 long int l, 92 double gHi, 93 long int Collider); 90 94 91 95 /**CS_l_mixing_VF01 Collision treatment based on Vrinceanu and Flannery 2001 -
trunk/source/helike_einsta.cpp
r1997 r2127 253 253 fixit(); // adding the 2-photon decay 2^3S - 1^1S may be important in early universe 254 254 A = ForbiddenHe[ipHi - 1]; 255 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f); 255 iso_put_error(ipHE_LIKE,nelem,ipHe2s3S ,ipHe1s1S,IPRAD, 0.01f, 0.20f); 256 iso_put_error(ipHE_LIKE,nelem,ipHe2s1S ,ipHe1s1S,IPRAD, 0.01f, 0.05f); 257 iso_put_error(ipHE_LIKE,nelem,ipHe2p3P0,ipHe1s1S,IPRAD, 0.00f, 0.00f); 258 iso_put_error(ipHE_LIKE,nelem,ipHe2p3P1,ipHe1s1S,IPRAD, 0.01f, 0.05f); 259 iso_put_error(ipHE_LIKE,nelem,ipHe2p3P2,ipHe1s1S,IPRAD, 0.01f, 0.01f); 256 260 } 257 261 else … … 293 297 TotalInsanity(); 294 298 } 295 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f); 296 } 297 /* For now just just put 1% error for forbidden lines. */ 298 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,.01f); 299 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f); 300 } 301 299 302 return A; 300 303 } … … 313 316 A = 8.0E-3 * exp(9.283/sqrt((double)N_(ipLo))) * pow((double)nelem,9.091) / 314 317 pow((double)N_(ipHi),2.877); 315 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f );318 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f); 316 319 } 317 320 … … 328 331 A *= (2.*(ipLo-3)+1.0)/3.0; 329 332 } 330 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f );333 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f); 331 334 } 332 335 … … 344 347 /* 20% error is based on difference between Savukov, Labzowsky, and Johnson (2005) 345 348 * and Lach and Pachucki (2001) for the helium transition. */ 346 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.2f );349 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.2f,0.2f); 347 350 } 348 351 else … … 351 354 * the above values to 10% or better. */ 352 355 A= 7.22E-9*pow((double)nelem, 9.33); 353 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.3f );356 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.3f,0.3f); 354 357 } 355 358 } … … 361 364 { 362 365 A = 1.549; 363 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f );366 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f); 364 367 } 365 368 else … … 369 372 * >>refercon astro-ph 0205163 */ 370 373 A= 0.1834*pow((double)nelem, 6.5735); 371 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f );374 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f); 372 375 } 373 376 } … … 383 386 * See discussion "Energy order within 2 3P" near the top of helike.c */ 384 387 A = 5.78E-12; 385 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f );388 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f); 386 389 387 390 } … … 396 399 * See discussion "Energy order within 2 3P" near the top of helike.c */ 397 400 A = 3.61E-15; 398 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f );401 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f); 399 402 } 400 403 … … 403 406 /* Current transition is not supported. */ 404 407 A = iso.SmallA; 405 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f); 406 } 407 408 /* For now just put 1% error for forbidden lines. */ 409 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,.01f); 408 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f); 409 } 410 410 411 411 ASSERT( A > 0.); … … 488 488 Aul = TransProbs[nelem][ipHi][ipLo]; 489 489 490 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.000 5f);490 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.0001f,0.002f); 491 491 } 492 492 … … 502 502 Aul = (1.59208e10) / pow(Eff_nupper,3.0); 503 503 ASSERT( Aul > 0.); 504 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.00 5f);504 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.0002f,0.002f); 505 505 } 506 506 … … 511 511 Aul = H_Einstein_A(nHi ,lHi , nLo , lLo , nelem); 512 512 ASSERT( Aul > 0.); 513 514 if( lHi + lLo >= 7 ) 515 { 516 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.001f); 517 } 518 else 519 { 520 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.01f); 521 } 513 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.006f,0.04f); 522 514 } 523 515 /* These fits come from extrapolations of Drake's oscillator strengths … … 616 608 617 609 ASSERT( Aul > 0. ); 618 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.0 1f);610 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.0002f,0.002f); 619 611 } 620 612 else … … 633 625 634 626 ASSERT( Aul > 0. ); 635 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.0 3f);627 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.01f,0.07f); 636 628 } 637 629 } … … 652 644 653 645 Aul = TransProbs[nelem][ipHi][ipLo]; 654 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f );646 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f); 655 647 } 656 648 … … 784 776 785 777 } 786 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f );778 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.1f,0.1f); 787 779 } 788 789 /* for now just give ions some a 5% error across the board */790 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.05f);791 780 } 792 781 } … … 1005 994 /* Neither upper nor lower is resolved...Aul's are easy. */ 1006 995 Aul = HydroEinstA( N_(ipLo), N_(ipHi) )*z4; 1007 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.001f );996 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.001f,0.001f); 1008 997 1009 998 ASSERT( Aul > 0.); … … 1034 1023 iso.CachedAs[ipISO][nelem][ N_(ipHi)-iso.n_HighestResolved_max[ipISO][nelem]-1 ][ ipLo ][1] = 0.f; 1035 1024 1036 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.01f );1025 iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPRAD,0.01f,0.01f); 1037 1026 ASSERT( Aul > 0.); 1038 1027 } -
trunk/source/helike_recom.cpp
r1960 r2127 4 4 /*cross_section - calculates the photoionization cross_section for a given level and photon energy*/ 5 5 /*radrecomb - calculates radiative recombination coefficients. */ 6 /*RecomInt - Integral in milne relation. Called by qg32. */7 6 /*He_cross_section returns cross section (cm^-2), 8 7 * given EgammaRyd, the photon energy in Ryd, … … 31 30 STATIC double X1Int( double u ); 32 31 STATIC double X2Int( double u ); 33 34 STATIC double cross_section(double EgammaRyd); 35 36 static double RecomInt(double EE); 32 STATIC double cross_section(double EgammaRyd, long nelem, long ipLev); 37 33 38 34 #if 0 … … 46 42 #endif 47 43 48 static double kTRyd,EthRyd,Xn_S59; 49 static long int ipLev,globalZ; 44 static double EthRyd,Xn_S59; 50 45 51 46 /* Here are several structures containing fits to cross-sections calculated in the following: */ … … 1512 1507 1513 1508 /* These corrections ensure helium cross-sections agree with the highly accurate 1514 * values given by Hummer and Storey (1998). The index is [ipLev - 1] */1509 * values given by Hummer and Storey (1998). The index is [ipLevel - 1] */ 1515 1510 /* >>refer He PCS Hummer, D.G., & Storey, P.J. 1998, MNRAS 297, 1073 */ 1516 1511 double PreferredThres[NUM_PREF_TCS] = { … … 1527 1522 13.54,26.58,55.71,-1.00,-1.00,-1.00,-1.00,-1.00,-1.00,58.07}; 1528 1523 1529 globalZ = nelem; 1530 1531 ipLev = ipLevel; 1532 1533 EthRyd = iso.xIsoLevNIonRyd[ipHE_LIKE][nelem][ipLev]; 1524 EthRyd = iso.xIsoLevNIonRyd[ipHE_LIKE][nelem][ipLevel]; 1534 1525 1535 1526 /* make sure this routine not called within collapsed high levels */ … … 1537 1528 1538 1529 /* this is a resolved level */ 1539 cs = (1.e-18)*cross_section( EgammaRyd );1530 cs = (1.e-18)*cross_section( EgammaRyd, nelem, ipLevel ); 1540 1531 1541 1532 /* Scale to preferred threshold values */ … … 1546 1537 if( ipCurrentLevel!=ipLevel ) 1547 1538 { 1548 ThreshCS = cross_section( EthRyd );1539 ThreshCS = cross_section( EthRyd, nelem, ipLevel ); 1549 1540 &