cloudy
trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /*HeCollid evaluate collisional rates */ 00004 /*HeCSInterp interpolate on He1 collision strengths */ 00005 /*AtomCSInterp do the atom */ 00006 /*IonCSInterp do the ions */ 00007 /*CS_l_mixing_PS64 - find rate for l-mixing collisions by protons, for neutrals */ 00008 #include "cddefines.h" 00009 #include "atmdat.h" 00010 #include "conv.h" 00011 #include "dense.h" 00012 #include "helike.h" 00013 #include "helike_cs.h" 00014 #include "hydro_vs_rates.h" 00015 #include "iso.h" 00016 #include "opacity.h" 00017 #include "phycon.h" 00018 #include "physconst.h" 00019 #include "rfield.h" 00020 #include "taulines.h" 00021 #include "thirdparty.h" 00022 #include "trace.h" 00023 00024 /* returns thermally-averaged Seaton 62 collision strength. */ 00025 STATIC double S62_Therm_ave_coll_str( double EProjectile_eV ); 00026 00027 /* all of these are used in the calculation of Stark collision strengths 00028 * following the algorithms in Vrinceanu & Flannery (2001). */ 00029 STATIC double Therm_ave_coll_str_int_VF01( double EProjectileRyd); 00030 STATIC double collision_strength_VF01( double velOrEner, bool lgParamIsRedVel ); 00031 STATIC double L_mix_integrand_VF01( double alpha ); 00032 STATIC double StarkCollTransProb_VF01( long int n, long int l, long int lp, double alpha, double deltaPhi); 00033 00034 static long int global_ipISO, global_n, global_l, global_l_prime, global_s, global_z, global_Collider; 00035 static double global_bmax, global_red_vel, global_an, global_collider_charge; 00036 static double global_I_energy_eV, global_deltaE, global_temp, global_osc_str, global_stat_weight; 00037 static double kTRyd; 00038 00039 /* These are masses relative to the proton mass of the electron, proton, he+, and alpha particle. */ 00040 static double ColliderMass[4] = {ELECTRON_MASS/PROTON_MASS, 1.0, 4.0, 4.0}; 00041 static double ColliderCharge[4] = {1.0, 1.0, 1.0, 2.0}; 00042 00043 void HeCollidSetup( void ) 00044 { 00045 /* this must be longer than data path string, set in path.h/cpu.cpp */ 00046 long i, i1, j, nelem, ipHi, ipLo; 00047 bool lgEOL, lgHIT; 00048 FILE *ioDATA; 00049 long ipISO = ipHE_LIKE; 00050 00051 # define chLine_LENGTH 1000 00052 char chLine[chLine_LENGTH]; 00053 00054 DEBUG_ENTRY( "HeCollidSetup()" ); 00055 00056 /* get the collision strength data for the He 1 lines */ 00057 if( trace.lgTrace ) 00058 fprintf( ioQQQ," HeCollidSetup opening he1_cs.dat:"); 00059 00060 ioDATA = open_data( "he1_cs.dat", "r" ); 00061 00062 /* check that magic number is ok */ 00063 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 00064 { 00065 fprintf( ioQQQ, " HeCollidSetup could not read first line of he1_cs.dat.\n"); 00066 cdEXIT(EXIT_FAILURE); 00067 } 00068 i = 1; 00069 i1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL); 00070 /*i2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL); 00071 i3 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);*/ 00072 00073 /* the following is to check the numbers that appear at the start of he1_cs.dat */ 00074 if( i1 !=COLLISMAGIC ) 00075 { 00076 fprintf( ioQQQ, 00077 " HeCollidSetup: the version of he1_cs.dat is not the current version.\n" ); 00078 fprintf( ioQQQ, 00079 " HeCollidSetup: I expected to find the number %i and got %li instead.\n" , 00080 COLLISMAGIC, i1 ); 00081 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine ); 00082 cdEXIT(EXIT_FAILURE); 00083 } 00084 00085 /* get the array of temperatures */ 00086 lgHIT = false; 00087 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL ) 00088 { 00089 /* only look at lines without '#' in first col */ 00090 if( chLine[0] != '#') 00091 { 00092 lgHIT = true; 00093 iso.nCS[ipISO] = 0; 00094 lgEOL = false; 00095 i = 1; 00096 while( !lgEOL && iso.nCS[ipISO] < HE1CSARRAY) 00097 { 00098 iso.CSTemp[ipISO][iso.nCS[ipISO]] = FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL); 00099 ++iso.nCS[ipISO]; 00100 } 00101 break; 00102 } 00103 } 00104 --iso.nCS[ipISO]; 00105 if( !lgHIT ) 00106 { 00107 fprintf( ioQQQ, " HeCollidSetup could not find line in CS temperatures in he1_cs.dat.\n"); 00108 cdEXIT(EXIT_FAILURE); 00109 } 00110 00111 /* create space for array of CS values, if not already done */ 00112 iso.HeCS.reserve( NISO ); 00113 iso.HeCS.reserve( ipISO, LIMELM ); 00114 00115 for( nelem=ipHELIUM; nelem<LIMELM; ++nelem ) 00116 { 00120 if( nelem != ipHELIUM ) 00121 continue; 00122 00123 /* only grab core for elements that are turned on */ 00124 if( nelem == ipHELIUM || dense.lgElmtOn[nelem]) 00125 { 00126 iso.HeCS.reserve( ipISO, nelem, iso.numLevels_max[ipHE_LIKE][nelem]-iso.nCollapsed_max[ipHE_LIKE][nelem] ); 00127 00128 for( ipHi=ipHe2s3S; ipHi < iso.numLevels_max[ipHE_LIKE][nelem] - iso.nCollapsed_max[ipHE_LIKE][nelem];++ipHi ) 00129 { 00130 iso.HeCS.reserve( ipISO, nelem, ipHi, ipHi ); 00131 00132 for( ipLo=ipHe1s1S; ipLo<ipHi; ++ipLo ) 00133 iso.HeCS.reserve( ipISO, nelem, ipHi, ipLo, iso.nCS[ipISO] ); 00134 } 00135 } 00136 } 00137 00138 iso.HeCS.alloc(); 00139 iso.HeCS.zero(); 00140 00141 /* now read in the CS data */ 00142 lgHIT = false; 00143 nelem = ipHELIUM; 00144 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL ) 00145 { 00146 char *chTemp; 00147 /* only look at lines without '#' in first col */ 00148 if( (chLine[0] == ' ') || (chLine[0]=='\n') ) 00149 break; 00150 if( chLine[0] != '#') 00151 { 00152 lgHIT = true; 00153 00154 /* get lower and upper level index */ 00155 j = 1; 00156 ipLo = (long)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL); 00157 ipHi = (long)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL); 00158 ASSERT( ipLo < ipHi ); 00159 if( ipHi >= iso.numLevels_max[ipHE_LIKE][nelem] - iso.nCollapsed_max[ipHE_LIKE][nelem] ) 00160 continue; 00161 else 00162 { 00163 chTemp = chLine; 00164 /* skip over 4 tabs to start of cs data */ 00165 for( i=0; i<3; ++i ) 00166 { 00167 if( (chTemp = strchr( chTemp, '\t' )) == NULL ) 00168 { 00169 fprintf( ioQQQ, " HeCollidSetup no init cs\n" ); 00170 cdEXIT(EXIT_FAILURE); 00171 } 00172 ++chTemp; 00173 } 00174 00175 for( i=0; i<iso.nCS[ipISO]; ++i ) 00176 { 00177 double a; 00178 if( (chTemp = strchr( chTemp, '\t' )) == NULL ) 00179 { 00180 fprintf( ioQQQ, " HeCollidSetup not scan cs, current indices: %li %li\n", ipHi, ipLo ); 00181 cdEXIT(EXIT_FAILURE); 00182 } 00183 ++chTemp; 00184 sscanf( chTemp , "%le" , &a ); 00185 iso.HeCS[ipISO][nelem][ipHi][ipLo][i] = (realnum)a; 00186 } 00187 } 00188 } 00189 } 00190 00191 /* close the data file */ 00192 fclose( ioDATA ); 00193 00194 return; 00195 } 00196 00197 /* Choose either AtomCSInterp or IonCSInterp */ 00198 realnum HeCSInterp(long int nelem, 00199 long int ipHi, 00200 long int ipLo, 00201 long int Collider ) 00202 { 00203 realnum cs, factor1; 00204 00205 /* This variable is for diagnostic purposes: 00206 * a string used in the output to designate where each cs comes from. */ 00207 const char *where = " "; 00208 00209 DEBUG_ENTRY( "HeCSInterp()" ); 00210 00211 if( !iso.lgColl_excite[ipHE_LIKE] ) 00212 { 00213 return (realnum)1E-10; 00214 } 00215 00216 if( nelem == ipHELIUM ) 00217 { 00218 /* do for helium */ 00219 cs = AtomCSInterp( nelem, ipHi , ipLo, &factor1, &where, Collider ); 00220 } 00221 else 00222 { 00223 /* get collision strengths for an ion */ 00224 cs = IonCSInterp( nelem, ipHi , ipLo, &factor1, &where, Collider ); 00225 } 00226 00227 /* set 15% errors on all collision rates. */ 00228 iso_put_error(ipHE_LIKE, nelem, ipHi, ipLo, IPCOLLIS, 0.15f ); 00229 00230 ASSERT( cs >= 0.f ); 00231 00232 /* in many cases the correction factor for split states has already been made, 00233 * if not then factor is still negative */ 00234 /* Remove the second test here when IonCSInterp is up to par with AtomCSInterp */ 00235 ASSERT( factor1 >= 0.f || nelem!=ipHELIUM ); 00236 if( factor1 < 0.f ) 00237 { 00238 ASSERT( iso.lgCS_Vriens[ipHE_LIKE] ); 00239 00240 factor1 = 1.f; 00241 } 00242 00243 /* take factor into account, usually 1, ratio of stat weights if within 2 3P 00244 * and with collisions from collapsed to resolved levels */ 00245 cs *= factor1; 00246 00247 { 00248 /*@-redef@*/ 00249 enum {DEBUG_LOC=false}; 00250 /*@+redef@*/ 00251 00252 if( DEBUG_LOC && ( nelem==ipOXYGEN ) && (cs > 0.f) && (StatesElem[ipHE_LIKE][nelem][ipHi].n == 2) 00253 && ( StatesElem[ipHE_LIKE][nelem][ipLo].n <= 2 ) ) 00254 fprintf(ioQQQ,"%li\t%li\t%li\t-\t%li\t%li\t%li\t%.2e\t%s\n", 00255 StatesElem[ipHE_LIKE][nelem][ipLo].n , StatesElem[ipHE_LIKE][nelem][ipLo].S , 00256 StatesElem[ipHE_LIKE][nelem][ipLo].l , StatesElem[ipHE_LIKE][nelem][ipHi].n , 00257 StatesElem[ipHE_LIKE][nelem][ipHi].S , StatesElem[ipHE_LIKE][nelem][ipHi].l , cs,where); 00258 } 00259 00260 return MAX2(cs, 1.e-10f); 00261 } 00262 00263 realnum AtomCSInterp(long int nelem, 00264 long int ipHi, 00265 long int ipLo, 00266 realnum *factor1, 00267 const char **where, 00268 long int Collider ) 00269 { 00270 long ipArray; 00271 realnum cs; 00272 00273 DEBUG_ENTRY( "AtomCSInterp()" ); 00274 00275 ASSERT( nelem == ipHELIUM ); 00276 00277 /* init values, better be positive when we exit */ 00278 cs = -1.f; 00279 00280 /* this may be used for splitting up the collision strength within 2 3P when 00281 * the lower level is withint 2 3P, and for collisions between resolved and collapsed levels. 00282 * It may be set somewhere in routine, so set to negative value here as flag saying not set */ 00283 *factor1 = -1.f; 00284 00285 /* for most of the helium iso sequence, the order of the J levels within 2 3P 00286 * in increasing energy, is 0 1 2 - the exception is atomic helium itself, 00287 * which is swapped, 2 1 0 */ 00288 00289 /* this branch is for upper and lower levels within 2p3P */ 00290 if( ipLo >= ipHe2p3P0 && ipHi <= ipHe2p3P2 && Collider==ipELECTRON ) 00291 { 00292 *factor1 = 1.f; 00293 /* atomic helium, use Berrington private comm cs */ 00294 00295 /* >>refer he1 cs Berrington, Keith, 2001, private communication - email follows 00296 > Dear Gary, 00297 > I could not find any literature on the He fine-structure 00298 > problem (but I didn't look very hard, so there may be 00299 > something somewhere). However, I did a quick R-matrix run to 00300 > see what the magnitude of the collision strengths are... At 00301 > 1000K, I get the effective collision strength for 2^3P J=0-1, 00302 > 0-2, 1-2 as 0.8,0.7,2.7; for 10000K I get 1.2, 2.1, 6.0 00303 */ 00304 /* indices are the same and correct, only thing special is that energies are in inverted order...was right first time. */ 00305 if( ipLo == ipHe2p3P0 && ipHi == ipHe2p3P1 ) 00306 { 00307 cs = 1.2f; 00308 } 00309 else if( ipLo == ipHe2p3P0 && ipHi == ipHe2p3P2 ) 00310 { 00311 cs = 2.1f; 00312 } 00313 else if( ipLo == ipHe2p3P1 && ipHi == ipHe2p3P2 ) 00314 { 00315 cs = 6.0f; 00316 } 00317 else 00318 { 00319 cs = 1.0f; 00320 TotalInsanity(); 00321 } 00322 00323 *where = "Berr "; 00324 /* statistical weights included */ 00325 } 00326 /* >>chng 02 feb 25, Bray data should come first since it is the best we have. */ 00327 /* this branch is the Bray et al data, for n <= 5, where quantal calcs exist 00328 * must exclude ipLo >= ipHe2p1P because they give no numbers for those */ 00329 else if( StatesElem[ipHE_LIKE][nelem][ipHi].n <= 5 && 00330 ( ipHi < iso.numLevels_max[ipHE_LIKE][nelem] - iso.nCollapsed_max[ipHE_LIKE][nelem] ) && 00331 iso.HeCS[ipHE_LIKE][nelem][ipHi][ipLo][0] >= 0.f && Collider== ipELECTRON ) 00332 { 00333 ASSERT( *factor1 == -1.f ); 00334 ASSERT( ipLo < ipHi ); 00335 ASSERT( ipHe2p3P0 == 3 ); 00336 00337 /* ipLo is within 2^3P */ 00338 if( ipLo >= ipHe2p3P0 && ipLo <= ipHe2p3P2 ) 00339 { 00340 /* *factor1 is ratio of statistical weights of level to term */ 00341 00342 /* ipHe2p3P0, ipHe2p3P1, ipHe2p3P2 have indices 3,4,and 5, but j=0,1,and 2. */ 00343 *factor1 = (2.f*((realnum)ipLo-3.f)+1.f) / 9.f; 00344 00345 /* ipHi must be above ipHe2p3P2 since 2p3Pj->2p3Pk taken care of above */ 00346 ASSERT( ipHi > ipHe2p3P2 ); 00347 } 00348 /* ipHi is within 2^3P */ 00349 else if( ipHi >= ipHe2p3P0 && ipHi <= ipHe2p3P2 ) 00350 { 00351 ASSERT( ipLo < ipHe2p3P0 ); 00352 00353 *factor1 = (2.f*((realnum)ipHi-3.f)+1.f) / 9.f; 00354 } 00355 /* neither are within 2^3P...no splitting necessary */ 00356 else 00357 { 00358 *factor1 = 1.f; 00359 } 00360 00361 /* SOME OF THESE ARE NOT N-CHANGING! */ 00362 /* Must be careful about turning each one on or off. */ 00363 00364 /* this is the case where we have quantal calculations */ 00365 /* >>refer He1 cs Bray, I., Burgess, A., Fursa, D.V., & Tully, J.A., 2000, A&AS, 146, 481-498 */ 00366 /* check whether we are outside temperature array bounds, 00367 * and use extreme value if we are */ 00368 if( phycon.alogte <= iso.CSTemp[ipHE_LIKE][0] ) 00369 { 00370 cs = iso.HeCS[ipHE_LIKE][nelem][ipHi][ipLo][0]; 00371 } 00372 else if( phycon.alogte >= iso.CSTemp[ipHE_LIKE][iso.nCS[ipHE_LIKE]-1] ) 00373 { 00374 cs = iso.HeCS[ipHE_LIKE][nelem][ipHi][ipLo][iso.nCS[ipHE_LIKE]-1]; 00375 } 00376 else 00377 { 00378 realnum flow; 00379 00380 /* find which array element within the cs vs temp array */ 00381 ipArray = (long)((phycon.alogte - iso.CSTemp[ipHE_LIKE][0])/(iso.CSTemp[ipHE_LIKE][1]-iso.CSTemp[ipHE_LIKE][0])); 00382 ASSERT( ipArray < iso.nCS[ipHE_LIKE] ); 00383 ASSERT( ipArray >= 0 ); 00384 /* when taking the average, this is the fraction from the lower temperature value */ 00385 flow = (realnum)( (phycon.alogte - iso.CSTemp[ipHE_LIKE][ipArray])/ 00386 (iso.CSTemp[ipHE_LIKE][ipArray+1]-iso.CSTemp[ipHE_LIKE][ipArray])); 00387 ASSERT( (flow >= 0.f) && (flow <= 1.f) ); 00388 00389 cs = iso.HeCS[ipHE_LIKE][nelem][ipHi][ipLo][ipArray] * (1.f-flow) + 00390 iso.HeCS[ipHE_LIKE][nelem][ipHi][ipLo][ipArray+1] * flow; 00391 } 00392 00393 *where = "Bray "; 00394 00395 /* options to kill collisional excitation and/or l-mixing */ 00396 if( StatesElem[ipHE_LIKE][nelem][ipHi].n == StatesElem[ipHE_LIKE][nelem][ipLo].n ) 00397 /* iso.lgColl_l_mixing turned off with atom he-like l-mixing collisions off command */ 00398 cs *= (realnum)iso.lgColl_l_mixing[ipHE_LIKE]; 00399 else 00400 { 00401 /* iso.lgColl_excite turned off with atom he-like collisional excitation off command */ 00402 cs *= (realnum)iso.lgColl_excite[ipHE_LIKE]; 00403 } 00404 00405 ASSERT( cs >= 0.f ); 00406 /* statistical weights included */ 00407 } 00408 /* this branch, n-same, l-changing collision, and not case of He with quantal data */ 00409 else if( (StatesElem[ipHE_LIKE][nelem][ipHi].n == StatesElem[ipHE_LIKE][nelem][ipLo].n ) && 00410 (StatesElem[ipHE_LIKE][nelem][ipHi].S == StatesElem[ipHE_LIKE][nelem][ipLo].S ) ) 00411 { 00412 ASSERT( *factor1 == -1.f ); 00413 *factor1 = 1.f; 00414 00415 /* ASSERT( StatesElem[ipHE_LIKE][nelem][ipHi].n >= 3 ); */ 00416 ASSERT( StatesElem[ipHE_LIKE][nelem][ipHi].n <= iso.n_HighestResolved_max[ipHE_LIKE][nelem] ); 00417 00418 if( (StatesElem[ipHE_LIKE][nelem][ipLo].l <=2) && 00419 abs(StatesElem[ipHE_LIKE][nelem][ipHi].l - StatesElem[ipHE_LIKE][nelem][ipLo].l)== 1 ) 00420 { 00421 /* Use the method given in 00422 * >>refer He CS Seaton, M. J. 1962, Proc. Phys. Soc. 79, 1105 00423 * statistical weights included */ 00424 cs = (realnum)CS_l_mixing_S62(ipHE_LIKE, nelem, ipLo, ipHi, (double)phycon.te, Collider); 00425 } 00426 else if( iso.lgCS_Vrinceanu[ipHE_LIKE] ) 00427 { 00428 if( StatesElem[ipHE_LIKE][nelem][ipLo].l >=3 && 00429 StatesElem[ipHE_LIKE][nelem][ipHi].l >=3 ) 00430 { 00431 /* Use the method given in 00432 * >>refer He CS Vrinceanu, D. \& Flannery, M. R. 2001, PhysRevA 63, 032701 00433 * statistical weights included */ 00434 cs = (realnum)CS_l_mixing_VF01( ipHE_LIKE, 00435 nelem, 00436 StatesElem[ipHE_LIKE][nelem][ipLo].n, 00437 StatesElem[ipHE_LIKE][nelem][ipLo].l, 00438 StatesElem[ipHE_LIKE][nelem][ipHi].l, 00439 StatesElem[ipHE_LIKE][nelem][ipLo].S, 00440 (double)phycon.te, 00441 Collider ); 00442 } 00443 else 00444 { 00445 cs = 0.f; 00446 } 00447 } 00448 /* this branch, l changing by one */ 00449 else if( abs(StatesElem[ipHE_LIKE][nelem][ipHi].l - StatesElem[ipHE_LIKE][nelem][ipLo].l)== 1) 00450 { 00451 /* >>refer He cs Pengelly, R.M., & Seaton, M.J., 1964, MNRAS, 127, 165 */ 00452 /* statistical weights included */ 00453 cs = (realnum)CS_l_mixing_PS64( ipHE_LIKE, nelem, ipLo, ipHi, Collider); 00454 } 00455 else 00456 { 00457 /* l changes by more than 1, but same-n collision */ 00458 cs = 0.f; 00459 } 00460 00461 /* ipLo is within 2^3P */ 00462 if( ipLo >= ipHe2p3P0 && ipLo <= ipHe2p3P2 ) 00463 { 00464 *factor1 = (2.f*((realnum)ipLo-3.f)+1.f) / 9.f; 00465 } 00466 00467 /* ipHi is within 2^3P */ 00468 if( ipHi >= ipHe2p3P0 && ipHi <= ipHe2p3P2 ) 00469 { 00470 *factor1 = (2.f*((realnum)ipHi-3.f)+1.f) / 9.f; 00471 } 00472 00473 *where = "lmix "; 00474 cs *= (realnum)iso.lgColl_l_mixing[ipHE_LIKE]; 00475 } 00476 00477 /* this is an atomic n-changing collision with no quantal calculations */ 00478 else if( StatesElem[ipHE_LIKE][nelem][ipHi].n != StatesElem[ipHE_LIKE][nelem][ipLo].n ) 00479 { 00480 ASSERT( *factor1 == -1.f ); 00481 /* this is an atomic n-changing collision with no quantal calculations */ 00482 /* gbar g-bar goes here */ 00483 00484 /* >>chng 02 jul 25, add option for fits to quantal cs data */ 00485 if( iso.lgCS_Vriens[ipHE_LIKE] ) 00486 { 00487 /* >>refer He CS Vriens, L., & Smeets, A.H.M. 1980, Phys Rev A 22, 940 00488 * statistical weight IS included in the routine */ 00489 cs = (realnum)CS_VS80( 00490 ipHE_LIKE, 00491 nelem, 00492 ipHi, 00493 ipLo, 00494 Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->Aul, 00495 phycon.te, 00496 Collider ); 00497 00498 *factor1 = 1.f; 00499 *where = "Vriens"; 00500 } 00501 else if( iso.lgCS_None[ipHE_LIKE] ) 00502 { 00503 cs = 0.f; 00504 *factor1 = 1.f; 00505 *where = "no gb"; 00506 } 00507 else if( iso.nCS_new[ipHE_LIKE] ) 00508 { 00509 *factor1 = 1.f; 00510 /* Don't know if stat weights are included in this, but they're probably 00511 * wrong anyway since they are based in part on the former (incorrect) 00512 * implementation of Vriens and Smeets rates */ 00513 00514 /* two different fits, allowed and forbidden */ 00515 if( Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->Aul > 1. ) 00516 { 00517 /* permitted lines - large A */ 00518 double x = 00519 log10(MAX2(34.7,Transitions[ipHE_LIKE][nelem][ipHi][ipLo].EnergyWN)); 00520 00521 if( iso.nCS_new[ipHE_LIKE] == 1 ) 00522 { 00523 /* this is the broken power law fit, passing through both quantal 00524 * calcs at high energy and asymptotically goes to VS at low energies */ 00525 if( x < 4.5 ) 00526 { 00527 /* low energy fit for permitted transitions */ 00528 cs = (realnum)pow( 10. , -1.45*x + 6.75); 00529 } 00530 else 00531 { 00532 /* higher energy fit for permitted transitions */ 00533 cs = (realnum)pow( 10. , -3.33*x+15.15); 00534 } 00535 } 00536 else if( iso.nCS_new[ipHE_LIKE] == 2 ) 00537 { 00538 /* single parallel fit for permitted transitions, runs parallel to VS */ 00539 cs = (realnum)pow( 10. , -2.3*x+10.3); 00540 } 00541 else 00542 TotalInsanity(); 00543 } 00544 else 00545 { 00546 /* fit for forbidden transitions */ 00547 if( Transitions[ipHE_LIKE][nelem][ipHi][ipLo].EnergyWN < 25119.f ) 00548 { 00549 cs = 0.631f; 00550 } 00551 else 00552 { 00553 cs = (realnum)pow(10., 00554 -3.*log10(Transitions[ipHE_LIKE][nelem][ipHi][ipLo].EnergyWN)+12.8); 00555 } 00556 } 00557 00558 *where = "newgb"; 00559 } 00560 else 00561 TotalInsanity(); 00562 00563 /* ipLo is within 2^3P */ 00564 if( ipLo >= ipHe2p3P0 && ipLo <= ipHe2p3P2 ) 00565 { 00566 *factor1 = (2.f*((realnum)ipLo-3.f)+1.f) / 9.f; 00567 } 00568 00569 /* ipHi is within 2^3P */ 00570 if( ipHi >= ipHe2p3P0 && ipHi <= ipHe2p3P2 ) 00571 { 00572 *factor1 = (2.f*((realnum)ipHi-3.f)+1.f) / 9.f; 00573 } 00574 00575 /* options to turn off collisions */ 00576 cs *= (realnum)iso.lgColl_excite[ipHE_LIKE]; 00577 00578 } 00579 else 00580 { 00581 /* If spin changing collisions are prohibited in the l-mixing routine, 00582 * they will fall here, and will have been assigned no collision strength. 00583 * Assign zero for now. */ 00584 ASSERT( StatesElem[ipHE_LIKE][nelem][ipHi].S != StatesElem[ipHE_LIKE][nelem][ipLo].S ); 00585 cs = 0.f; 00586 *factor1 = 1.f; 00587 } 00588 00589 ASSERT( cs >= 0.f ); 00590 00591 /*iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPCOLLIS,-1);*/ 00592 00593 return(cs); 00594 00595 } 00596 00597 /* IonCSInterp interpolate on collision strengths for element nelem */ 00598 realnum IonCSInterp( long nelem , long ipHi , long ipLo, realnum *factor1, const char **where, long Collider ) 00599 { 00600 realnum cs; 00601 00602 DEBUG_ENTRY( "IonCSInterp()" ); 00603 00604 ASSERT( nelem > ipHELIUM ); 00605 ASSERT( nelem < LIMELM ); 00606 00607 /* init values, better be positive when we exit */ 00608 cs = -1.f; 00609 00610 /* this may be used for splitting up the collision strength for collisions between 00611 * resolved and collapsed levels. It may be set somewhere in routine, so set to 00612 * negative value here as flag saying not set */ 00613 *factor1 = -1.f; 00614 00615 00616 /* >>chng 02 mar 04, the approximation here for transitions within 2p3P was not needed, 00617 * because the Zhang data give these transitions. They are of the same order, but are 00618 * specific to the three transitions */ 00619 00620 /* this branch is ground to n=2 or from n=2 to n=2, for ions only */ 00621 /*>>refer Helike CS Zhang, Honglin, & Sampson, Douglas H. 1987, ApJS 63, 487 */ 00622 if( StatesElem[ipHE_LIKE][nelem][ipHi].n==2 00623 && StatesElem[ipHE_LIKE][nelem][ipLo].n<=2 && Collider==ipELECTRON) 00624 { 00625 *where = "Zhang"; 00626 *factor1 = 1.; 00627 00628 /* Collisions from gound */ 00629 if( ipLo == ipHe1s1S ) 00630 { 00631 switch( ipHi ) 00632 { 00633 case 1: /* to 2tripS */ 00634 cs = 0.25f/POW2(nelem+1.f); 00635 break; 00636 case 2: /* to 2singS */ 00637 cs = 0.4f/POW2(nelem+1.f); 00638 break; 00639 case 3: /* to 2tripP0 */ 00640 cs = 0.15f/POW2(nelem+1.f); 00641 break; 00642 case 4: /* to 2tripP1 */ 00643 cs = 0.45f/POW2(nelem+1.f); 00644 break; 00645 case 5: /* to 2tripP2 */ 00646 cs = 0.75f/POW2(nelem+1.f); 00647 break; 00648 case 6: /* to 2singP */ 00649 cs = 1.3f/POW2(nelem+1.f); 00650 break; 00651 default: 00652 TotalInsanity(); 00653 break; 00654 } 00655 cs *= (realnum)iso.lgColl_excite[ipHE_LIKE]; 00656 } 00657 /* collisions from 2tripS to n=2 */ 00658 else if( ipLo == ipHe2s3S ) 00659 { 00660 switch( ipHi ) 00661 { 00662 case 2: /* to 2singS */ 00663 cs = 2.75f/POW2(nelem+1.f); 00664 break; 00665 case 3: /* to 2tripP0 */ 00666 cs = 60.f/POW2(nelem+1.f); 00667 break; 00668 case 4: /* to 2tripP1 */ 00669 cs = 180.f/POW2(nelem+1.f); 00670 break; 00671 case 5: /* to 2tripP2 */ 00672 cs = 300.f/POW2(nelem+1.f); 00673 break; 00674 case 6: /* to 2singP */ 00675 cs = 5.8f/POW2(nelem+1.f); 00676 break; 00677 default: 00678 TotalInsanity(); 00679 break; 00680 } 00681 cs *= (realnum)iso.lgColl_l_mixing[ipHE_LIKE]; 00682 } 00683 /* collisions from 2singS to n=2 */ 00684 else if( ipLo == ipHe2s1S ) 00685 { 00686 switch( ipHi ) 00687 { 00688 case 3: /* to 2tripP0 */ 00689 cs = 0.56f/POW2(nelem+1.f); 00690 break; 00691 case 4: /* to 2tripP1 */ 00692 cs = 1.74f/POW2(nelem+1.f); 00693 break; 00694 case 5: /* to 2tripP2 */ 00695 cs = 2.81f/POW2(nelem+1.f); 00696 break; 00697 case 6: /* to 2singP */ 00698 cs = 190.f/POW2(nelem+1.f); 00699 break; 00700 default: 00701 TotalInsanity(); 00702 break; 00703 } 00704 cs *= (realnum)iso.lgColl_l_mixing[ipHE_LIKE]; 00705 } 00706 /* collisions from 2tripP0 to n=2 */ 00707 else if( ipLo == ipHe2p3P0 ) 00708 { 00709 switch( ipHi ) 00710 { 00711 case 4: /* to 2tripP1 */ 00712 cs = 8.1f/POW2(nelem+1.f); 00713 break; 00714 case 5: /* to 2tripP2 */ 00715 cs = 8.2f/POW2(nelem+1.f); 00716 break; 00717 case 6: /* to 2singP */ 00718 cs = 3.9f/POW2(nelem+1.f); 00719 break; 00720 default: 00721 TotalInsanity(); 00722 break; 00723 } 00724 cs *= (realnum)iso.lgColl_l_mixing[ipHE_LIKE]; 00725 } 00726 /* collisions from 2tripP1 to n=2 */ 00727 else if( ipLo == ipHe2p3P1 ) 00728 { 00729 switch( ipHi ) 00730 { 00731 case 5: /* to 2tripP2 */ 00732 cs = 30.f/POW2(nelem+1.f); 00733 break; 00734 case 6: /* to 2singP */ 00735 cs = 11.7f/POW2(nelem+1.f); 00736 break; 00737 default: 00738 TotalInsanity(); 00739 break; 00740 } 00741 cs *= (realnum)iso.lgColl_l_mixing[ipHE_LIKE]; 00742 } 00743 /* collisions from 2tripP2 to n=2 */ 00744 else if( ipLo == ipHe2p3P2 ) 00745 { 00746 /* to 2singP */ 00747 cs = 19.5f/POW2(nelem+1.f) * (realnum)iso.lgColl_l_mixing[ipHE_LIKE]; 00748 } 00749 else 00750 TotalInsanity(); 00751 00752 /* statistical weights included */ 00753 } 00754 00755 /* this branch, n-same, l-changing collisions */ 00756 else if( (StatesElem[ipHE_LIKE][nelem][ipHi].n == StatesElem[ipHE_LIKE][nelem][ipLo].n ) && 00757 (StatesElem[ipHE_LIKE][nelem][ipHi].S == StatesElem[ipHE_LIKE][nelem][ipLo].S ) ) 00758 { 00759 ASSERT( *factor1 == -1.f ); 00760 *factor1 = 1.f; 00761 00762 ASSERT( StatesElem[ipHE_LIKE][nelem][ipHi].n <= iso.n_HighestResolved_max[ipHE_LIKE][nelem] ); 00763 00764 if( (StatesElem[ipHE_LIKE][nelem][ipLo].l <=2) && 00765 abs(StatesElem[ipHE_LIKE][nelem][ipHi].l - StatesElem[ipHE_LIKE][nelem][ipLo].l)== 1 ) 00766 { 00767 /* Use the method given in 00768 * >>refer He CS Seaton, M. J. 1962, Proc. Phys. Soc. 79, 1105 00769 * statistical weights included */ 00770 cs = (realnum)CS_l_mixing_S62(ipHE_LIKE, nelem, ipLo, ipHi, (double)phycon.te, Collider); 00771 } 00772 else if( iso.lgCS_Vrinceanu[ipHE_LIKE] ) 00773 { 00774 if( StatesElem[ipHE_LIKE][nelem][ipLo].l >=3 && 00775 StatesElem[ipHE_LIKE][nelem][ipHi].l >=3 ) 00776 { 00777 /* Use the method given in 00778 * >>refer He CS Vrinceanu, D. \& Flannery, M. R. 2001, PhysRevA 63, 032701 00779 * statistical weights included */ 00780 cs = (realnum)CS_l_mixing_VF01( ipHE_LIKE, 00781 nelem, 00782 StatesElem[ipHE_LIKE][nelem][ipLo].n, 00783 StatesElem[ipHE_LIKE][nelem][ipLo].l, 00784 StatesElem[ipHE_LIKE][nelem][ipHi].l, 00785 StatesElem[ipHE_LIKE][nelem][ipLo].S, 00786 (double)phycon.te, 00787 Collider ); 00788 } 00789 else 00790 { 00791 cs = 0.f; 00792 } 00793 } 00794 /* this branch, l changing by one */ 00795 else if( abs(StatesElem[ipHE_LIKE][nelem][ipHi].l - StatesElem[ipHE_LIKE][nelem][ipLo].l)== 1) 00796 { 00797 /* >>refer He cs Pengelly, R.M., & Seaton, M.J., 1964, MNRAS, 127, 165 */ 00798 /* statistical weights included */ 00799 cs = (realnum)CS_l_mixing_PS64( ipHE_LIKE, nelem, ipLo, ipHi, Collider); 00800 } 00801 else 00802 { 00803 /* l changes by more than 1, but same-n collision */ 00804 cs = 0.f; 00805 } 00806 00807 /* ipHi is within 2^3P, do not need to split on ipLo. */ 00808 if( ipHi >= ipHe2p3P0 && ipHi <= ipHe2p3P2 ) 00809 { 00810 *factor1 = (2.f*((realnum)ipHi-3.f)+1.f) / 9.f; 00811 } 00812 00813 *where = "lmix "; 00814 cs *= (realnum)iso.lgColl_l_mixing[ipHE_LIKE]; 00815 } 00816 00817 /* this branch, n changing collisions for ions */ 00818 else if( StatesElem[ipHE_LIKE][nelem][ipHi].n != StatesElem[ipHE_LIKE][nelem][ipLo].n ) 00819 { 00820 if( iso.lgCS_Vriens[ipHE_LIKE] ) 00821 { 00822 /* this is the default branch */ 00823 /* >>refer He CS Vriens, L., & Smeets, A.H.M. 1980, Phys Rev A 22, 940 00824 * statistical weight is NOT included in the routine */ 00825 cs = (realnum)CS_VS80( 00826 ipHE_LIKE, 00827 nelem, 00828 ipHi, 00829 ipLo, 00830 Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->Aul, 00831 phycon.te, 00832 Collider ); 00833 00834 *factor1 = 1.f; 00835 *where = "Vriens"; 00836 } 00837 else 00838 { 00839 /* \todo 2 this branch and above now do the same thing. Change something. */ 00840 /* this branch is for testing and reached with command ATOM HE-LIKE COLLISIONS VRIENS OFF */ 00841 fixit(); /* use Percival and Richards here */ 00842 00843 cs = 0.f; 00844 *where = "hydro"; 00845 } 00846 } 00847 else 00848 { 00849 /* what's left are deltaN=0, spin changing collisions. 00850 * These have not been accounted for. */ 00851 /* Make sure what comes here is what we think it is. */ 00852 ASSERT( StatesElem[ipHE_LIKE][nelem][ipHi].n == StatesElem[ipHE_LIKE][nelem][ipLo].n ); 00853 ASSERT( StatesElem[ipHE_LIKE][nelem][ipHi].S != StatesElem[ipHE_LIKE][nelem][ipLo].S ); 00854 cs = 0.f; 00855 *where = "spin "; 00856 } 00857 00858 ASSERT( cs >= 0.f ); 00859 00860 /*iso_put_error(ipHE_LIKE,nelem,ipHi,ipLo,IPCOLLIS,-1);*/ 00861 00862 return(cs); 00863 } 00864 00865 00866 /*CS_l_mixing_S62 - find rate for l-mixing collisions by protons, for neutrals */ 00867 /* The S62 stands for Seaton 1962 */ 00868 double CS_l_mixing_S62( 00869 long ipISO, 00870 long nelem /* the chemical element, 1 for He */, 00871 long ipLo /* lower level, 0 for ground */, 00872 long ipHi /* upper level, 0 for ground */, 00873 double temp, 00874 long Collider) 00875 { 00876 /* >>refer He l-mixing Seaton, M.J., 1962, Proc. Phys. Soc. */ 00877 double coll_str, deltaE; 00878 00879 DEBUG_ENTRY( "CS_l_mixing_S62()" ); 00880 00881 if( Transitions[ipISO][nelem][ipHi][ipLo].ipCont <= 0 ) 00882 { 00883 return 0.; 00884 } 00885 00886 global_temp = temp; 00887 global_stat_weight = StatesElem[ipISO][nelem][ipLo].g; 00888 /* >>chng 05 sep 06, RP - update energies of excited states */ 00889 /* global_deltaE = EVRYD*(iso.xIsoLevNIonRyd[ipISO][nelem][ipLo] - 00890 iso.xIsoLevNIonRyd[ipISO][nelem][ipHi]); */ 00891 global_deltaE = Transitions[ipISO][nelem][ipHi][ipLo].EnergyErg/EN1EV; 00892 deltaE = global_deltaE; 00893 global_I_energy_eV = EVRYD*iso.xIsoLevNIonRyd[ipISO][nelem][0]; 00894 global_Collider = Collider; 00895 00896 ASSERT( TRANS_PROB_CONST*POW2(deltaE/WAVNRYD/EVRYD) > 0. ); 00897 00898 global_osc_str = Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul/ 00899 (TRANS_PROB_CONST*POW2(deltaE/WAVNRYD/EVRYD)); 00900 00901 /* This returns a thermally averaged collision strength */ 00902 coll_str = qg32( 0.0, 1.0, S62_Therm_ave_coll_str); 00903 coll_str += qg32( 1.0, 10.0, S62_Therm_ave_coll_str); 00904 ASSERT( coll_str > 0. ); 00905 00906 /*iso_put_error(ipISO,nelem,ipHi,ipLo,ipCollis,-1);*/ 00907 00908 return coll_str; 00909 } 00910 00911 /* The integrand for calculating the thermal average of collision strengths */ 00912 STATIC double S62_Therm_ave_coll_str( double proj_energy_overKT ) 00913 { 00914 00915 double integrand, cross_section, /*Rnot,*/ osc_strength, coll_str, zOverB2; 00916 double deltaE, /* betanot, */ betaone, zeta_of_betaone, /* cs1, */ cs2, temp, stat_weight; 00917 double I_energy_eV, Dubya, proj_energy; 00918 long i, Collider; 00919 00920 DEBUG_ENTRY( "S62_Therm_ave_coll_str()" ); 00921 00922 /* projectile energy in eV */ 00923 proj_energy = proj_energy_overKT * phycon.te / EVDEGK; 00924 00925 deltaE = global_deltaE; 00926 osc_strength = global_osc_str; 00927 temp = (double)global_temp; 00928 stat_weight = global_stat_weight; 00929 I_energy_eV = global_I_energy_eV; 00930 Collider = global_Collider; 00931 00932 /* Rnot = 1.1229*H_BAR/sqrt(ELECTRON_MASS*deltaE*EN1EV)/Bohr_radius; in units of Bohr_radius */ 00933 00934 proj_energy *= ColliderMass[ipELECTRON]/ColliderMass[Collider]; 00935 /* The deltaE here is to make sure that the collider has no less 00936 * than the energy difference between the initial and final levels. */ 00937 proj_energy += deltaE; 00938 Dubya = 0.5*(2.*proj_energy-deltaE); 00939 ASSERT( Dubya > 0. ); 00940 00941 /* betanot = sqrt(proj_energy/I_energy_eV)*(deltaE/2./Dubya)*Rnot; */ 00942 00943 ASSERT( I_energy_eV > 0. ); 00944 ASSERT( osc_strength > 0. ); 00945 00946 /* from equation 33 */ 00947 zOverB2 = 0.5*POW2(Dubya/deltaE)*deltaE/I_energy_eV/osc_strength; 00948 00949 ASSERT( zOverB2 > 0. ); 00950 00951 if( zOverB2 > 100. ) 00952 { 00953 betaone = sqrt( 1./zOverB2 ); 00954 } 00955 else if( zOverB2 < 0.54 ) 00956 { 00957 /* Low betaone approximation */ 00958 betaone = (1./3.)*(log(PI)-log(zOverB2)+1.28); 00959 if( betaone > 2.38 ) 00960 { 00961 /* average with this over approximation */ 00962 betaone += 0.5*(log(PI)-log(zOverB2)); 00963 betaone *= 0.5; 00964 } 00965 } 00966 else 00967 { 00968 long ip_zOverB2 = 0; 00969 double zetaOVERbeta2[101] = { 00970 1.030E+02,9.840E+01,9.402E+01,8.983E+01,8.583E+01,8.200E+01,7.835E+01,7.485E+01, 00971 7.151E+01,6.832E+01,6.527E+01,6.236E+01,5.957E+01,5.691E+01,5.436E+01,5.193E+01, 00972 4.961E+01,4.738E+01,4.526E+01,4.323E+01,4.129E+01,3.943E+01,3.766E+01,3.596E+01, 00973 3.434E+01,3.279E+01,3.131E+01,2.989E+01,2.854E+01,2.724E+01,2.601E+01,2.482E+01, 00974 2.369E+01,2.261E+01,2.158E+01,2.059E+01,1.964E+01,1.874E+01,1.787E+01,1.705E+01, 00975 1.626E+01,1.550E+01,1.478E+01,1.409E+01,1.343E+01,1.280E+01,1.219E+01,1.162E+01, 00976 1.107E+01,1.054E+01,1.004E+01,9.554E+00,9.094E+00,8.655E+00,8.234E+00,7.833E+00, 00977 7.449E+00,7.082E+00,6.732E+00,6.397E+00,6.078E+00,5.772E+00,5.481E+00,5.202E+00, 00978 4.937E+00,4.683E+00,4.441E+00,4.210E+00,3.989E+00,3.779E+00,3.578E+00,3.387E+00, 00979 3.204E+00,3.031E+00,2.865E+00,2.707E+00,2.557E+00,2.414E+00,2.277E+00,2.148E+00, 00980 2.024E+00,1.907E+00,1.795E+00,1.689E+00,1.589E+00,1.493E+00,1.402E+00,1.316E+00, 00981 1.235E+00,1.157E+00,1.084E+00,1.015E+00,9.491E-01,8.870E-01,8.283E-01,7.729E-01, 00982 7.206E-01,6.712E-01,6.247E-01,5.808E-01,5.396E-01}; 00983 00984 ASSERT( zOverB2 >= zetaOVERbeta2[100] ); 00985 00986 /* find beta in the table */ 00987 for( i=0; i< 100; ++i ) 00988 { 00989 if( ( zOverB2 < zetaOVERbeta2[i] ) && ( zOverB2 >= zetaOVERbeta2[i+1] ) ) 00990 { 00991 /* found the temperature - use it */ 00992 ip_zOverB2 = i; 00993 break; 00994 } 00995 } 00996 00997 ASSERT( (ip_zOverB2 >=0) && (ip_zOverB2 < 100) ); 00998 00999 betaone = (zOverB2 - zetaOVERbeta2[ip_zOverB2]) / 01000 (zetaOVERbeta2[ip_zOverB2+1] - zetaOVERbeta2[ip_zOverB2]) * 01001 (pow(10., ((double)ip_zOverB2+1.)/100. - 1.) - pow(10., ((double)ip_zOverB2)/100. - 1.)) + 01002 pow(10., (double)ip_zOverB2/100. - 1.); 01003 01004 } 01005 01006 zeta_of_betaone = zOverB2 * POW2(betaone); 01007 01008 /* cs1 = betanot * bessel_k0(betanot) * bessel_k1(betanot); */ 01009 cs2 = 0.5*zeta_of_betaone + betaone * bessel_k0(betaone) * bessel_k1(betaone); 01010 01011 /* cross_section = MIN2(cs1, cs2); */ 01012 cross_section = cs2; 01013 01014 /* cross section in units of PI * a_o^2 */ 01015 cross_section *= 8. * (I_energy_eV/deltaE) * osc_strength * (I_energy_eV/proj_energy); 01016 01017 /* convert to collision strength */ 01018 coll_str = cross_section * stat_weight * ( proj_energy/EVRYD ); 01019 01020 integrand = exp( -1.*(proj_energy-deltaE)*EVDEGK/temp ) * coll_str; 01021 return integrand; 01022 } 01023 01024 /*CS_l_mixing_PS64 - find rate for l-mixing collisions by protons, for neutrals */ 01025 double CS_l_mixing_PS64( 01026 long ipISO, 01027 long nelem /* the chemical element, 1 for He */, 01028 long ipLo /* lower level, 0 for ground */, 01029 long ipHi /* upper level, 0 for ground */, 01030 long Collider) 01031 { 01032 /* >>refer H-like l-mixing Pengelly, R.M., & Seaton, M.J., 1964, MNRAS, 127, 165 */ 01033 /* >>refer He-like l-mixing Pengelly, R.M., & Seaton, M.J., 1964, MNRAS, 127, 165 */ 01034 double cs, Dnl, 01035 TwoLogDebye, TwoLogRc1, 01036 factor1, factor2, 01037 bestfactor, factorpart, 01038 reduced_mass, reduced_mass_2_emass, 01039 rate, tau; 01040 const double ChargIncoming = ColliderCharge[Collider]; 01041 01042 DEBUG_ENTRY( "CS_l_mixing_PS64()" ); 01043 01044 ASSERT( ipHi > ipLo ); 01045 /* In this routine, two different cutoff radii are calculated, and as per PS64, 01046 * the least of these is selected. We take the least positive result. */ 01047 01048 /* This reduced mass is in grams. */ 01049 reduced_mass = dense.AtomicWeight[nelem]*ColliderMass[Collider]/ 01050 (dense.AtomicWeight[nelem]+ColliderMass[Collider])*ATOMIC_MASS_UNIT; 01051 /* this mass always appears relative to the electron mass, so define it that way */ 01052 reduced_mass_2_emass = reduced_mass / ELECTRON_MASS; 01053 01054 /* This is the lifetime of ipLo. */ 01055 tau = StatesElem[ipISO][nelem][ipLo].lifetime; 01056 01057 if( ipISO == ipH_LIKE && Transitions[ipISO][nelem][ipLo][ipH1s].ipCont > 0 ) 01058 { 01059 tau = (1./tau) + Transitions[ipISO][nelem][ipLo][ipH1s].Emis->Aul; 01060 tau = 1./tau; 01061 } 01062 01063 /* equation 46 of PS64 */ 01064 /* min on density added to prevent this from becoming large and negative 01065 * at very high n_e - Pengelly & Seaton did not apply this above 1e12 cm^-3 */ 01066 /* This is 2 times the log of the Debye radius. */ 01067 //TwoLogDebye = 1.68 + log10( phycon.te / MIN2(1e11 , dense.eden ) ); 01068 /* Brocklehurst (1971, equation 3.40) has 1.181 instead of 1.68. This appears to be due to 01069 * Pengelly and Seaton neglecting one of the two factors of PI in their Equation 33 */ 01070 TwoLogDebye = 1.181 + log10( phycon.te / MIN2(1e10 , dense.eden ) ); 01071 01072 /* This is w times the log of cutoff = 0.72v(tau), where tau is the lifetime of the level nl. */ 01073 TwoLogRc1 = 10.95 + log10( phycon.te * tau * tau / reduced_mass_2_emass ); 01074 01075 /* this is equation 44 of PS64 */ 01076 Dnl = POW2( ChargIncoming / (double)(nelem + 1. - ipISO)) * 6. * POW2( (double)StatesElem[ipISO][nelem][ipLo].n ) * 01077 ( POW2((double)StatesElem[ipISO][nelem][ipLo].n) - 01078 POW2((double)StatesElem[ipISO][nelem][ipLo].l) - StatesElem[ipISO][nelem][ipLo].l - 1); 01079 01080 ASSERT( Dnl > 0. ); 01081 ASSERT( phycon.te / Dnl / reduced_mass_2_emass > 0. ); 01082 01083 factorpart = (11.54 + log10( phycon.te / Dnl / reduced_mass_2_emass ) ); 01084 01085 if( (factor1 = factorpart + TwoLogDebye) <= 0.) 01086 factor1 = BIGDOUBLE; 01087 if( (factor2 = factorpart + TwoLogRc1) <= 0.) 01088 factor2 = BIGDOUBLE; 01089 01090 /* Now we find the least positive result. */ 01091 bestfactor = MIN2(factor1,factor2); 01092 01093 ASSERT( bestfactor > 0. ); 01094 01095 /* If both factors are bigger than 100, toss out the result, and return SMALLFLOAT. */ 01096 if( bestfactor > 100. ) 01097 { 01098 return SMALLFLOAT; 01099 } 01100 01101 /* This is the rate coefficient. Units: cm^3 s-1 */ 01102 rate = 9.93e-6 * sqrt( reduced_mass_2_emass ) * Dnl / phycon.sqrte * bestfactor; 01103 01104 /***** NB NB NB NB 01105 * Brocklehurst (1971) has a factor of electron density in the rate coefficient (equation 3.38). 01106 * This is NOT a proper rate, as can be seen in his equations 2.2 and 2.4. This differs 01107 * from the formulism given by PS64. */ 01108 //rate *= dense.eden; 01109 01110 /* this is the TOTAL rate from nl to nl+/-1. So assume we can 01111 * divide by two to get the rate nl to either nl+1 or nl-1. */ 01112 if( StatesElem[ipISO][nelem][ipLo].l > 0 ) 01113 rate /=2; 01114 01115 /* convert rate to collision strength */ 01116 /* NB - the term in parentheses corrects for the fact that COLL_CONST is only appropriate 01117 * for electron colliders and is off by reduced_mass_2_emass^-1.5 */ 01118 cs = rate / ( COLL_CONST * pow(reduced_mass_2_emass, -1.5) ) * 01119 phycon.sqrte * StatesElem[ipISO][nelem][ipHi].g; 01120 01121 ASSERT( cs > 0. ); 01122 01123 /*iso_put_error(ipISO,nelem,ipHi,ipLo,ipCollis,-1);*/ 01124 01125 return cs; 01126 } 01127 01128 /*CS_l_mixing - find collision strength for l-mixing collisions for neutrals */ 01129 /* The VF stands for Vrinceanu & Flannery 2001 */ 01130 /* >>refer He l-mixing Vrinceanu, D. & Flannery, M. R. 2001, PhysRevA 63, 032701 */ 01131 /* >>refer He l-mixing Hezel, T. P., Burkhardt, C. E., Ciocca, M., He, L-W., */ 01132 /* >>refercon Leventhal, J. J. 1992, Am. J. Phys. 60, 329 */ 01133 double CS_l_mixing_VF01(long int ipISO, 01134 long int nelem, 01135 long int n, 01136 long int l, 01137 long int lp, 01138 long int s, 01139 double temp, 01140 long int Collider ) 01141 { 01142 01143 double coll_str; 01144 01145 DEBUG_ENTRY( "CS_l_mixing_VF01()" ); 01146 01147 global_ipISO = ipISO; 01148 global_z = nelem; 01149 global_n = n; 01150 global_l = l; 01151 global_l_prime = lp; 01152 global_s = s; 01153 global_temp = temp; 01154 global_Collider = Collider; 01155 global_collider_charge = ColliderCharge[Collider]; 01156 ASSERT( global_collider_charge > 0. ); 01157 01158 /* no need to do this for h-like */ 01159 if( ipISO > ipH_LIKE ) 01160 { 01161 ASSERT( l != 0 ); 01162 ASSERT( lp != 0 ); 01163 } 01164 01165 kTRyd = temp / TE1RYD; 01166 if( !iso.lgCS_therm_ave[ipISO] ) 01167 { 01168 /* Must do some thermal averaging for densities greater 01169 * than about 10000 and less than about 1e10, 01170 * because kT gives significantly different results. 01171 * Still, do sparser integration than is done below */ 01172 if( (dense.eden > 10000.) && (dense.eden < 1E10 ) ) 01173 { 01174 coll_str = qg32( 0.0, 6.0, Therm_ave_coll_str_int_VF01); 01175 } 01176 else 01177 { 01178 /* Do NOT average over Maxwellian */ 01179 coll_str = collision_strength_VF01( kTRyd, false ); 01180 } 01181 } 01182 else 01183 { 01184 /* DO average over Maxwellian */ 01185 coll_str = qg32( 0.0, 1.0, Therm_ave_coll_str_int_VF01); 01186 coll_str += qg32( 1.0, 10.0, Therm_ave_coll_str_int_VF01); 01187 } 01188 01189 return coll_str; 01190 } 01191 01192 /* The integrand for calculating the thermal average of collision strengths */ 01193 STATIC double Therm_ave_coll_str_int_VF01( double EOverKT ) 01194 { 01195 double integrand; 01196 01197 DEBUG_ENTRY( "Therm_ave_coll_str_int_VF01()" ); 01198 01199 integrand = exp( -1.*EOverKT ) * collision_strength_VF01( EOverKT * kTRyd, false ); 01200 01201 return integrand; 01202 } 01203 01204 STATIC double collision_strength_VF01( double velOrEner, bool lgParamIsRedVel ) 01205 { 01206 01207 double cross_section, coll_str, RMSv, aveRadius, reduced_vel, E_Proj_Ryd; 01208 double ConstantFactors, reduced_mass, CSIntegral, stat_weight; 01209 double ColliderCharge = global_collider_charge; 01210 double quantum_defect1, quantum_defect2, omega_qd1, omega_qd2, omega_qd; 01211 double reduced_b_max, reduced_b_min, alphamax, alphamin, step, alpha1 /*, alpha2*/; 01212 01213 long ipISO, nelem, n, l, lp, s, ipLo, ipHi, Collider = global_Collider; 01214 01215 DEBUG_ENTRY( "collision_strength_VF01()" ); 01216 01217 nelem = global_z; 01218 n = global_n; 01219 ASSERT( n > 0 ); 01220 l = global_l; 01221 lp = global_l_prime; 01222 s = global_s; 01223 ipISO = global_ipISO; 01224 01225 /* >>chng 06 may 30, move these up from below. Also ipHi needs lp not l. */ 01226 ipLo = iso.QuantumNumbers2Index[ipISO][nelem][n][l][s]; 01227 ipHi = iso.QuantumNumbers2Index[ipISO][nelem][n][lp][s]; 01228 01229 stat_weight = StatesElem[ipISO][nelem][ipLo].g; 01230 01231 /* no need to do this for h-like */ 01232 if( ipISO > ipH_LIKE ) 01233 { 01234 /* these shut up the lint, already done above */ 01235 ASSERT( l > 0 ); 01236 ASSERT( lp > 0 ); 01237 } 01238 01239 /* This reduced mass is in grams. */ 01240 reduced_mass = dense.AtomicWeight[nelem]*ColliderMass[Collider]/ 01241 (dense.AtomicWeight[nelem]+ColliderMass[Collider])*ATOMIC_MASS_UNIT; 01242 ASSERT( reduced_mass > 0. ); 01243 01244 /* this is root mean squared velocity */ 01245 /* use this as projectile velocity for thermally-averaged cross-section */ 01246 aveRadius = (BOHR_RADIUS_CM/((double)nelem+1.-(double)ipISO))*POW2((double)n); 01247 ASSERT( aveRadius < 1.e-4 ); 01248 /* >>chng 05 jul 14, as per exchange with Ryan Porter & Peter van Hoof, avoid 01249 * roundoff error and give ability to go beyond zinc */ 01250 /*ASSERT( aveRadius >= BOHR_RADIUS_CM );*/ 01251 ASSERT( aveRadius > 3.9/LIMELM * BOHR_RADIUS_CM ); 01252 global_an = aveRadius; 01253 01254 /* vn = n * H_BAR/ m / r = Z * e^2 / n / H_BAR 01255 * where Z is the effective charge. */ 01256 RMSv = ((double)nelem+1.-(double)ipISO)*POW2(ELEM_CHARGE_ESU)/(double)n/H_BAR; 01257 ASSERT( RMSv > 0. ); 01258 01259 ASSERT( ColliderMass[Collider] > 0. ); 01260 01261 if( lgParamIsRedVel ) 01262 { 01263 /* velOrEner is a reduced velocity */ 01264 reduced_vel = velOrEner; 01265 /* The proton mass is needed here because the ColliderMass array is 01266 * expressed in units of the proton mass, but here we need absolute mass. */ 01267 E_Proj_Ryd = 0.5 * POW2( reduced_vel * RMSv ) * ColliderMass[Collider] * 01268 PROTON_MASS / EN1RYD; 01269 } 01270 else 01271 { 01272 /* velOrEner is a projectile energy in Rydbergs */ 01273 E_Proj_Ryd = velOrEner; 01274 reduced_vel = sqrt( 2.*E_Proj_Ryd*EN1RYD/ColliderMass[Collider]/PROTON_MASS )/RMSv; 01275 } 01276 01277 /* put limits on the reduced velocity. These limits should be more than fair. */ 01278 ASSERT( reduced_vel > 1.e-10 ); 01279 ASSERT( reduced_vel < 1.e10 ); 01280 01281 global_red_vel = reduced_vel; 01282 01283 /* Factors outside integral */ 01284 ConstantFactors = 4.5*PI*POW2(ColliderCharge*aveRadius/reduced_vel); 01285 01286 /* Reduced here means in units of aveRadius: */ 01287 reduced_b_min = 1.5 * ColliderCharge / reduced_vel; 01288 ASSERT( reduced_b_min > 1.e-10 ); 01289 ASSERT( reduced_b_min < 1.e10 ); 01290 01291 if( ipISO == ipH_LIKE ) 01292 { 01293 /* Debye radius: appears to be too large, results in 1/v^2 variation. */ 01294 reduced_b_max = sqrt( BOLTZMANN*global_temp/ColliderCharge/dense.eden ) 01295 / (PI2*ELEM_CHARGE_ESU)/aveRadius; 01296 } 01297 else if( ipISO == ipHE_LIKE ) 01298 { 01299 quantum_defect1 = (double)n- (double)nelem/sqrt(iso.xIsoLevNIonRyd[ipISO][nelem][ipLo]); 01300 quantum_defect2 = (double)n- (double)nelem/sqrt(iso.xIsoLevNIonRyd[ipISO][nelem][ipHi]); 01301 01302 /* The magnitude of each quantum defect must be between zero and one. */ 01303 ASSERT( fabs(quantum_defect1) < 1.0 ); 01304 ASSERT( fabs(quantum_defect1) > 0.0 ); 01305 ASSERT( fabs(quantum_defect2) < 1.0 ); 01306 ASSERT( fabs(quantum_defect2) > 0.0 ); 01307 01308 /* The quantum defect precession frequencies */ 01309 omega_qd1 = fabs( 5.* quantum_defect1 * (1.-0.6*POW2((double)l/(double)n)) / POW3( (double)n ) / (double)l ); 01310 /* >>chng 06 may 30, this needs lp not l. */ 01311 omega_qd2 = fabs( 5.* quantum_defect2 * (1.-0.6*POW2((double)lp/(double)n)) / POW3( (double)n ) / (double)lp ); 01312 /* Take the average for the two levels, for reciprocity. */ 01313 omega_qd = 0.5*( omega_qd1 + omega_qd2 ); 01314 01315 ASSERT( omega_qd > 0. ); 01316 01317 reduced_b_max = sqrt( 1.5 * ColliderCharge * n / omega_qd )/aveRadius; 01318 } 01319 else 01320 /* rethink this before using on other iso sequences. */ 01321 TotalInsanity(); 01322 01323 reduced_b_max = MAX2( reduced_b_max, reduced_b_min ); 01324 ASSERT( reduced_b_max > 0. ); 01325 global_bmax = reduced_b_max*aveRadius; 01326 alphamin = 1.5*ColliderCharge/(reduced_vel * reduced_b_max); 01327 alphamax = 1.5*ColliderCharge/(reduced_vel * reduced_b_min); 01328 01329 ASSERT( alphamin > 0. ); 01330 ASSERT( alphamax > 0. ); 01331 01332 alphamin = MAX2( alphamin, 1.e-30 ); 01333 alphamax = MAX2( alphamax, 1.e-20 ); 01334 01335 CSIntegral = 0.; 01336 01337 if( alphamax > alphamin ) 01338 { 01339 01340 step = (alphamax - alphamin)/5.; 01341 alpha1 = alphamin; 01342 CSIntegral += qg32( alpha1, alpha1+step, L_mix_integrand_VF01); 01343 CSIntegral += qg32( alpha1+step, alpha1+4.*step, L_mix_integrand_VF01); 01344 } 01345 01346 /* Calculate cross section */ 01347 cross_section = ConstantFactors * CSIntegral; 01348 01349 /* convert to collision strength */ 01350 coll_str = cross_section * stat_weight / PI / BOHR_RADIUS_CM / BOHR_RADIUS_CM; 01351 coll_str *= E_Proj_Ryd; 01352 01353 coll_str = MAX2( SMALLFLOAT, coll_str); 01354 01355 return coll_str; 01356 } 01357 01358 STATIC double L_mix_integrand_VF01( double alpha ) 01359 { 01360 double integrand, deltaPhi, b; 01361 01362 DEBUG_ENTRY( "L_mix_integrand_VF01()" ); 01363 01364 ASSERT( alpha >= 1.e-30 ); 01365 ASSERT( global_bmax > 0. ); 01366 ASSERT( global_red_vel > 0. ); 01367 01368 /* >>refer He l-mixing Kazansky, A. K. & Ostrovsky, V. N. 1996, JPhysB: At. Mol. Opt. Phys. 29, 3651 */ 01369 b = 1.5*global_collider_charge*global_an/(global_red_vel * alpha); 01370 /* deltaPhi is the effective angle swept out by the projector as viewed by the target. */ 01371 if( b < global_bmax ) 01372 { 01373 deltaPhi = -1.*PI + 2.*asin(b/global_bmax); 01374 } 01375 else 01376 { 01377 deltaPhi = 0.; 01378 } 01379 integrand = 1./alpha/alpha/alpha; 01380 integrand *= StarkCollTransProb_VF01( global_n, global_l, global_l_prime, alpha, deltaPhi ); 01381 01382 return integrand; 01383 } 01384 01385 STATIC double StarkCollTransProb_VF01( long n, long l, long lp, double alpha, double deltaPhi) 01386 { 01387 double probability; 01388 double cosU1, cosU2, sinU1, sinU2, cosChiOver2, sinChiOver2, cosChi, A, B; 01389 01390 DEBUG_ENTRY( "StarkCollTransProb_VF01()" ); 01391 01392 ASSERT( alpha > 0. ); 01393 01394 /* These are defined on page 11 of VF01 */ 01395 cosU1 = 2.*POW2((double)l/(double)n) - 1.; 01396 cosU2 = 2.*POW2((double)lp/(double)n) - 1.; 01397 01398 sinU1 = sqrt( 1. - cosU1*cosU1 ); 01399 sinU2 = sqrt( 1. - cosU2*cosU2 ); 01400 01401 01402 cosChiOver2 = (1. + alpha*alpha*cos( sqrt(1.+alpha*alpha) * deltaPhi ) )/(1.+alpha*alpha); 01403 sinChiOver2 = sqrt( 1. - cosChiOver2*cosChiOver2 ); 01404 cosChi = 2. * POW2( cosChiOver2 ) - 1.; 01405 01406 if( l == 0 ) 01407 { 01408 if( -1.*cosU2 - cosChi < 0. ) 01409 { 01410 probability = 0.; 01411 } 01412 else 01413 { 01414 /* Here is the initial state S case */ 01415 ASSERT( sinChiOver2 > 0. ); 01416 ASSERT( sinChiOver2*sinChiOver2 > POW2((double)lp/(double)n) ); 01417 /* This is equation 35 of VF01. There is a factor of hbar missing in the denominator of the 01418 * paper, but it's okay if you use atomic units (hbar = 1). */ 01419 probability = (double)lp/(POW2((double)n)*sinChiOver2*sqrt( POW2(sinChiOver2) - POW2((double)lp/(double)n) ) ); 01420 } 01421 } 01422 else 01423 { 01424 double OneMinusCosChi = 1. - cosChi; 01425 01426 if( OneMinusCosChi == 0. ) 01427 { 01428 double hold = sin( deltaPhi / 2. ); 01429 /* From approximation at bottom of page 10 of VF01. */ 01430 OneMinusCosChi = 8. * alpha * alpha * POW2( hold ); 01431 } 01432 01433 if( OneMinusCosChi == 0. ) 01434 { 01435 probability = 0.; 01436 } 01437 else 01438 { 01439 /* Here is the general case */ 01440 A = (cosU1*cosU2 - sinU1*sinU2 - cosChi)/OneMinusCosChi; 01441 B = (cosU1*cosU2 + sinU1*sinU2 - cosChi)/OneMinusCosChi; 01442 01443 ASSERT( B > A ); 01444 01445 /* These are the three cases of equation 34. */ 01446 if( B <= 0. ) 01447 { 01448 probability = 0.; 01449 } 01450 else 01451 { 01452 ASSERT( POW2( sinChiOver2 ) > 0. ); 01453 01454 probability = 2.*lp/(PI* /*H_BAR* */ n*n*POW2( sinChiOver2 )); 01455 01456 if( A < 0. ) 01457 { 01458 probability *= ellpk( -A/(B-A) ); 01459 probability /= sqrt( B - A ); 01460 } 01461 else 01462 { 01463 probability *= ellpk( A/B ); 01464 probability /= sqrt( B ); 01465 } 01466 } 01467 } 01468 01469 } 01470 01471 return probability; 01472 01473 } 01474 01475 #if 0 01476 STATIC void updateHeColl(long int nelem) 01477 { 01478 long int ipLo , ipHi, i; 01479 01480 /* collision strengths are assumed to be roughly constant for small changes in temperature 01481 * and are not recalculated as often as other data. This array stores the last temperature 01482 * at which collision strengths were evaluated for each species of the isoelectronic sequence. */ 01483 static double TeUsedForCS[LIMELM]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0., 01484 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0}; 01485 01486 DEBUG_ENTRY("updateHeColl()"); 01487 /* Calculate initial value of collision strengths, OR 01488 * Reevaluate collision strengths if the temperature has changed by more than 15%. */ 01489 01490 /* This was the case in the initial coding -- if the number/order of 01491 * colliders changes, must make sure all the auxiliary information 01492 * is correct. In the long term the code should be robust to 01493 * removing these conditions */ 01494 ASSERT(ipELECTRON == 0 && ipPROTON == 1 && ipHE_PLUS == 2); 01495 01496 if( (TeUsedForCS[nelem] == 0.) || 01497 ( TeUsedForCS[nelem]/phycon.te > 1.15 ) || 01498 ( TeUsedForCS[nelem]/phycon.te < 0.85 ) ) 01499 { 01500 for( ipHi=ipHe2s3S; ipHi <iso.numLevels_max[ipHE_LIKE][nelem]; ipHi++ ) 01501 { 01502 for( ipLo=ipHe1s1S; ipLo < ipHi; ipLo++ ) 01503 { 01504 for (i=0;i<=ipNCOLLIDER;i++) 01505 { 01506 /* Should remove this limit when collider data is complete */ 01507 if (i > ipHE_PLUS) 01508 continue; 01509 01510 /* collsion strengths due to impact */ 01511 Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Coll.csi[i] = 01512 HeCSInterp( nelem , ipHi , ipLo, i ); 01513 01514 /* check for sanity */ 01515 ASSERT( Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Coll.csi[i] >= 0. ); 01516 } 01517 } 01518 } 01519 /* Update temperature at which collision strengths were evaluated. */ 01520 TeUsedForCS[nelem] = phycon.te; 01521 } 01522 } 01523 #endif 01524