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 /*HCSAR_interp interpolate on collision strengths */ 00004 /*C6cs123 line collision rates for lower levels of hydrogenic carbon, n=1,2,3 */ 00005 /*Ca20cs123 line collision rates for lower levels of hydrogenic calcium, n=1,2,3 */ 00006 /*Hydcs123 Hydrogenic de-excitation collision rates n=1,2,3 */ 00007 /*H1cs123 hydrogen collision data levels involving 1s,2s,2p,3. */ 00008 /*Ne10cs123 line collision rates for lower levels of hydrogenic neon, n=1,2,3 */ 00009 /*He2cs123 line collision strengths for lower levels of helium ion, n=1,2,3, by K Korista */ 00010 /*Fe26cs123 line collision rates for lower levels of hydrogenic iron, n=1,2,3 */ 00011 #include "cddefines.h" 00012 #include "atmdat.h" 00013 #include "dense.h" 00014 #include "helike_cs.h" 00015 #include "hydro_vs_rates.h" 00016 #include "iso.h" 00017 #include "opacity.h" 00018 #include "phycon.h" 00019 #include "physconst.h" 00020 #include "taulines.h" 00021 00022 STATIC double Fe26cs123(long int i, long int j); 00023 STATIC double He2cs123(long int i, long int j); 00024 STATIC double Hydcs123(long int ilow, long int ihigh, long int iz, long int chType); 00025 STATIC double C6cs123(long int i, long int j); 00026 STATIC double Ca20cs123(long int i, long int j); 00027 STATIC double Ne10cs123(long int i, long int j); 00028 00029 STATIC realnum HCSAR_interp( int ipLo , int ipHi ); 00030 STATIC double CS_ThermAve_PR78(long ipISO, long nelem, long nHi, long nLo, double deltaE, double temp ); 00031 STATIC double Therm_ave_coll_str_int_PR78( double EOverKT ); 00032 STATIC double CS_PercivalRichards78( double Ebar ); 00033 00034 static long global_ipISO, global_nelem, global_nHi, global_nLo; 00035 static double kTRyd, global_deltaE; 00036 00037 static const realnum HCSTE[NHCSTE] = {5802.f,11604.f,34812.f,58020.f,116040.f,174060.f,232080.f,290100.f}; 00038 00039 /*HCSAR_interp interpolate on collision strengths */ 00040 STATIC realnum HCSAR_interp( int ipLo , int ipHi ) 00041 { 00042 00043 static int ip=1; 00044 realnum cs; 00045 00046 DEBUG_ENTRY( "HCSAR_interp()" ); 00047 00048 if( ipLo==1 && ipHi==2 ) 00049 { 00050 fprintf(ioQQQ,"HCSAR_interp was called for the 2s-2p transition, which it cannot do\n"); 00051 cdEXIT(EXIT_FAILURE); 00052 } 00053 if( phycon.te <= HCSTE[0] ) 00054 { 00055 cs = t_ADfA::Inst().h_coll_str( ipLo , ipHi , 0 ); 00056 } 00057 else if( phycon.te >= HCSTE[NHCSTE-1] ) 00058 { 00059 cs = t_ADfA::Inst().h_coll_str( ipLo , ipHi , NHCSTE-1 ); 00060 } 00061 else 00062 { 00063 /* the ip index is most likely correct since it points to the last temperature */ 00064 if( (HCSTE[ip-1] >= phycon.te ) || ( phycon.te > HCSTE[ip]) ) 00065 { 00066 /* we must find the temperature in the array */ 00067 for( ip=1; ip<NHCSTE; ++ip ) 00068 { 00069 if( (HCSTE[ip-1] < phycon.te ) && ( phycon.te <= HCSTE[ip]) ) 00070 break; 00071 } 00072 } 00073 /* we now have the index */ 00074 cs = t_ADfA::Inst().h_coll_str( ipLo , ipHi , ip-1 ) + 00075 (t_ADfA::Inst().h_coll_str( ipLo , ipHi , ip ) - t_ADfA::Inst().h_coll_str( ipLo , ipHi , ip-1 ) ) / (HCSTE[ip]-HCSTE[ip-1] ) * 00076 ((realnum)phycon.te - HCSTE[ip-1] ); 00077 if( cs <= 0.) 00078 { 00079 fprintf(ioQQQ," insane cs returned by HCSAR_interp, values are\n"); 00080 fprintf(ioQQQ,"%.3f %.3f \n", t_ADfA::Inst().h_coll_str( ipLo , ipHi , ip-1 ),t_ADfA::Inst().h_coll_str( ipLo , ipHi , ip ) ); 00081 } 00082 } 00083 return(cs); 00084 } 00085 00086 /*Hydcs123 Hydrogenic de-excitation collision strengths between levels n=1,2,3, 00087 * for any charge. routine only called by HydroCSInterp to fill in hydroline arrays 00088 * with collision strengths */ 00089 STATIC double Hydcs123( 00090 /* lower principal quantum number, 1, 2, or 3, in this routine 00091 * 1 is 1s, 2 is 2s, 3 is 2p, and 4 is 3s, 5 is 3p, and 6 is 3d */ 00092 long int ipLow, 00093 /* upper principal quantum nubmer, 2, 3, or 4 */ 00094 long int ipHi, 00095 /* charge, 0 for hydrogen, 1 for helium, etc */ 00096 long int nelem, 00097 /* = 'e' for electron collisions, ='p' for proton */ 00098 long int chType) 00099 { 00100 long int i, 00101 j, 00102 k; 00103 double C, 00104 D, 00105 EE, 00106 expq , 00107 Hydcs123_v, 00108 Ratehigh, 00109 Ratelow, 00110 TeUse, 00111 gLo, 00112 gHi, 00113 q, 00114 rate, 00115 slope, 00116 temp, 00117 temphigh, 00118 templow, 00119 tev, 00120 x, 00121 QuanNLo, 00122 QuanNUp, 00123 Charge, 00124 ChargeSquared, 00125 zhigh, 00126 zlow; 00127 static const double ap[5] = {-2113.113,729.0084,1055.397,854.632,938.9912}; 00128 static const double bp[5] = {-6783.515,-377.7190,724.1936,493.1107,735.7466}; 00129 static const double cp[5] = {-3049.719,226.2320,637.8630,388.5465,554.6369}; 00130 static const double dp[5] = {3514.5153,88.60169,-470.4055,-329.4914,-450.8459}; 00131 static const double ep[5] = {0.005251557,0.009059154,0.008725781,0.009952418,0.01098687}; 00132 static const double ae[5] = {-767.5859,-643.1189,-461.6836,-429.0543,-406.5285}; 00133 static const double be[5] = {-1731.9178,-1442.548,-1055.364,-980.3079,-930.9266}; 00134 static const double ce[5] = {-939.1834,-789.9569,-569.1451,-530.1974,-502.0939}; 00135 static const double de[5] = {927.4773,773.2008,564.3272,524.2944,497.7763}; 00136 static const double ee[5] = {-0.002528027,-0.003793665,-0.002122103,-0.002234207,-0.002317720}; 00137 static const double A[2] = {4.4394,0.0}; 00138 static const double B[2] = {0.8949,0.8879}; 00139 static const double C0[2] = {-0.6012,-0.2474}; 00140 static const double C1[2] = {-3.9710,-3.7562}; 00141 static const double C2[2] = {-4.2176,2.0491}; 00142 static const double D0[2] = {2.930,0.0539}; 00143 static const double D1[2] = {1.7990,3.4009}; 00144 static const double D2[2] = {4.9347,-1.7770}; 00145 00146 DEBUG_ENTRY( "Hydcs123()" ); 00147 /* Hydrogenic de-excitation collision rates n=1,2,3 00148 * >>refer h1 cs Callaway, J. 1983, Phys Let A, 96, 83 00149 * >>refer h1 cs Zygelman, B., & Dalgarno, A. 1987, Phys Rev A, 35, 4085 00150 * for 2p-2s only. 00151 * The fit from Callaway is in nuclear charge for 1s - 2s,2p only. 00152 * For transtions involving level 3, interpolation in Z involving 00153 * the functions He2cs123,C6cs123,Ne10cs123,Ca20cs123, Fe26cs123. 00154 * 00155 * The fits from ZD are for 2p-2s for Z=2,6,12,16,18 only other charges are 00156 * interpolated, both electron and proton rates are included, 00157 * the variable chType is either 'e' or 'p'. 00158 * 00159 * ipLow is the lower level and runs from 1 to 3 (1s, 2s, 2p) 00160 * ipHi is the upper level and runs from 2 to 6 (2s, 2p, 3s, 3p, 3d) */ 00161 00162 /* for Callaway fit: */ 00163 /* for Zygelman and Dalgarno: */ 00164 00165 /* first entry is 2p, then 2s */ 00166 00167 /* fit in nuclear charge Z for 2p-2s collisions in hydrogenic species 00168 * equation is a+bx+cx^2ln(x)+dexp(x)+eln(x)/x^2, where x=te/Z^2 in a.u. 00169 * first are the proton rates: */ 00170 /* these are electron rates: */ 00171 00172 /* following is for charged species */ 00173 /* charge is on scale with He+=1, Li++=2, etc */ 00174 ASSERT( nelem > ipHYDROGEN ); 00175 ASSERT( nelem < LIMELM ); 00176 00177 /* these are the pointers to upper and lower levels. 1=1s, 2=2s, 3=2p, 4=3 */ 00178 ASSERT( ipLow > 0); 00179 ASSERT( ipLow <= 3); 00180 ASSERT( ipHi > 1 ); 00181 ASSERT( ipHi <=6 ); 00182 00183 /* set quantum numbers and stat. weights of the transitions: */ 00184 if( ipHi == 6 ) 00185 { 00186 /* upper is n=3 then set level, stat. weight */ 00187 QuanNUp = 3.; 00188 gHi = 10.; 00189 /* following will be set here even though it is not used in this case, 00190 * to prevent good compilers from falsing on i not set, 00191 * there is assert when used to make sure it is ok */ 00192 i = -1; 00193 } 00194 else if( ipHi == 5 ) 00195 { 00196 /* upper is n=3 then set level, stat. weight */ 00197 QuanNUp = 3.; 00198 gHi = 6.; 00199 /* following will be set here even though it is not used in this case, 00200 * to prevent good compilers from falsing on i not set, 00201 * there is assert when used to make sure it is ok */ 00202 i = -1; 00203 } 00204 else if( ipHi == 4 ) 00205 { 00206 /* upper is n=3 then set level, stat. weight */ 00207 QuanNUp = 3.; 00208 gHi = 2.; 00209 /* following will be set here even though it is not used in this case, 00210 * to prevent good compilers from falsing on i not set, 00211 * there is assert when used to make sure it is ok */ 00212 i = -1; 00213 } 00214 else if( ipHi == 3 ) 00215 { 00216 /* upper is nl=2p then set level, stat. weight */ 00217 QuanNUp = 2.; 00218 gHi = 6.; 00219 /* used to point within vectors defined above */ 00220 i = 0; 00221 } 00222 else if( ipHi == 2 ) 00223 { 00224 /* upper is nl=2s then set level, stat. weight */ 00225 QuanNUp = 2.; 00226 gHi = 2.; 00227 /* used to point within vectors defined above */ 00228 i = 1; 00229 } 00230 else 00231 { 00232 fprintf( ioQQQ, " Insane levels in Hydcs123\n" ); 00233 cdEXIT(EXIT_FAILURE); 00234 } 00235 00236 /* which lower level? */ 00237 if( ipLow == 1 ) 00238 { 00239 /* lower is n=1 then set level, stat. weight */ 00240 QuanNLo = 1.; 00241 gLo = 2.; 00242 } 00243 else if( ipLow == 2 ) 00244 { 00245 /* lower is nl=2s then set level, stat. weight */ 00246 QuanNLo = 2.; 00247 gLo = 2.; 00248 } 00249 else if( ipLow == 3 ) 00250 { 00251 QuanNLo = 2.; 00252 gLo = 6.; 00253 } 00254 else 00255 { 00256 fprintf( ioQQQ, " Insane levels in Hydcs123\n" ); 00257 cdEXIT(EXIT_FAILURE); 00258 } 00259 00260 /* following is the physical charge */ 00261 Charge = (double)(nelem + 1); 00262 /* square of charge */ 00263 ChargeSquared = Charge*Charge; 00264 00265 if( ipLow == 2 && ipHi == 3 ) 00266 { 00267 /*************************************** this is 2p-2s: 00268 * series of if statements determines which entries Charge is between. */ 00269 if( nelem < 1 ) 00270 { 00271 /* this can't happen since routine returned above when ip=1 and 00272 * special atomic hydrogen routine called */ 00273 fprintf( ioQQQ, " insane charge given to Hydcs123\n" ); 00274 cdEXIT(EXIT_FAILURE); 00275 } 00276 else if( nelem == 1 ) 00277 { 00278 zlow = 2.; 00279 j = 1; 00280 zhigh = 2.; 00281 k = 1; 00282 } 00283 /* Li through C */ 00284 else if( nelem <= 5 ) 00285 { 00286 zlow = 2.; 00287 j = 1; 00288 zhigh = 6.; 00289 k = 2; 00290 } 00291 else if( nelem <= 11 ) 00292 { 00293 zlow = 6.; 00294 j = 2; 00295 zhigh = 12.; 00296 k = 3; 00297 } 00298 else if( nelem <= 15 ) 00299 { 00300 zlow = 12.; 00301 j = 3; 00302 zhigh = 16.; 00303 k = 4; 00304 } 00305 else if( nelem <= 17 ) 00306 { 00307 zlow = 16.; 00308 j = 4; 00309 zhigh = 18.; 00310 k = 5; 00311 } 00312 /* following changed to else from else if, 00313 * to prevent false comment in good compilers */ 00314 /*else if( nelem > 18 )*/ 00315 else 00316 { 00317 zlow = 18.; 00318 j = 5; 00319 zhigh = 18.; 00320 k = 5; 00321 } 00322 00323 /* convert Te to a.u./Z^2 00324 * determine rate at the low Charge */ 00325 x = EVRYD/TE1RYD*phycon.te/(27.211396*pow2(zlow)); 00326 TeUse = MIN2(x,0.80); 00327 x = MAX2(0.025,TeUse); 00328 00329 /* what type of collision are we dealing with? */ 00330 if( chType == 'e' ) 00331 { 00332 /* electron collisions */ 00333 Ratelow = ae[j-1] + be[j-1]*x + ce[j-1]*pow2(x)*log(x) + de[j-1]* 00334 exp(x) + ee[j-1]*log(x)/pow2(x); 00335 } 00336 else if( chType == 'p' ) 00337 { 00338 Ratelow = ap[j-1] + bp[j-1]*x + cp[j-1]*pow2(x)*log(x) + dp[j-1]* 00339 exp(x) + ep[j-1]*log(x)/pow2(x); 00340 } 00341 else 00342 { 00343 /* this can't happen */ 00344 fprintf( ioQQQ, " insane collision species given to Hydcs123\n" ); 00345 cdEXIT(EXIT_FAILURE); 00346 } 00347 00348 /* determine rate at the high Charge */ 00349 x = EVRYD/TE1RYD*phycon.te/(27.211396*pow2(zhigh)); 00350 TeUse = MIN2(x,0.80); 00351 x = MAX2(0.025,TeUse); 00352 if( chType == 'e' ) 00353 { 00354 Ratehigh = ae[k-1] + be[k-1]*x + ce[k-1]*pow2(x)*log(x) + 00355 de[k-1]*exp(x) + ee[k-1]*log(x)/pow2(x); 00356 } 00357 else 00358 { 00359 Ratehigh = ap[k-1] + bp[k-1]*x + cp[k-1]*pow2(x)*log(x) + 00360 dp[k-1]*exp(x) + ep[k-1]*log(x)/pow2(x); 00361 } 00362 /* linearly interpolate in charge */ 00363 if( fp_equal( zlow, zhigh ) ) 00364 { 00365 rate = Ratelow; 00366 } 00367 else 00368 { 00369 slope = (Ratehigh - Ratelow)/(zhigh - zlow); 00370 rate = slope*(Charge - zlow) + Ratelow; 00371 } 00372 rate = rate/ChargeSquared/Charge*1.0e-7; 00373 /* must convert to cs and need to know the valid temp range */ 00374 templow = 0.025*27.211396*TE1RYD/EVRYD*ChargeSquared; 00375 temphigh = 0.80*27.211396*TE1RYD/EVRYD*ChargeSquared; 00376 TeUse = MIN2((double)phycon.te,temphigh); 00377 temp = MAX2(TeUse,templow); 00378 Hydcs123_v = rate*gHi*sqrt(temp)/COLL_CONST; 00379 00380 if( chType == 'p' ) 00381 { 00382 /* COLL_CONST is incorrect for protons, correct here */ 00383 Hydcs123_v *= pow( PROTON_MASS/ELECTRON_MASS, 1.5 ); 00384 } 00385 } 00386 else if( ipHi == 4 || ipHi == 5 || ipHi == 6 ) 00387 { 00388 /* n = 3 00389 * for the rates involving n = 3, must do something different. */ 00390 if( nelem < 1 ) 00391 { 00392 fprintf( ioQQQ, " insane charge given to Hydcs123\n" ); 00393 cdEXIT(EXIT_FAILURE); 00394 } 00395 else if( nelem == 1 ) 00396 { 00397 zlow = 2.; 00398 Ratelow = He2cs123(ipLow,ipHi); 00399 zhigh = 2.; 00400 Ratehigh = Ratelow; 00401 } 00402 else if( nelem > 1 && nelem <= 5 ) 00403 { 00404 zlow = 2.; 00405 Ratelow = He2cs123(ipLow,ipHi); 00406 zhigh = 6.; 00407 Ratehigh = C6cs123(ipLow,ipHi); 00408 } 00409 else if( nelem > 5 && nelem <= 9 ) 00410 { 00411 zlow = 6.; 00412 Ratelow = C6cs123(ipLow,ipHi); 00413 zhigh = 10.; 00414 Ratehigh = Ne10cs123(ipLow,ipHi); 00415 } 00416 else if( nelem > 9 && nelem <= 19 ) 00417 { 00418 zlow = 10.; 00419 Ratelow = Ne10cs123(ipLow,ipHi); 00420 zhigh = 20.; 00421 Ratehigh = Ca20cs123(ipLow,ipHi); 00422 } 00423 else if( nelem > 19 && nelem <= 25 ) 00424 { 00425 zlow = 20.; 00426 Ratelow = Ca20cs123(ipLow,ipHi); 00427 zhigh = 26.; 00428 Ratehigh = Fe26cs123(ipLow,ipHi); 00429 } 00430 /*>>>chng 98 dec 17, to else to stop comment from good compilers*/ 00431 /*else if( nelem > 26 )*/ 00432 else 00433 { 00434 Charge = 26.; 00435 zlow = 26.; 00436 Ratelow = Fe26cs123(ipLow,ipHi); 00437 zhigh = 26.; 00438 Ratehigh = Ratelow; 00439 } 00440 00441 /* linearly interpolate */ 00442 if( fp_equal( zlow, zhigh ) ) 00443 { 00444 rate = Ratelow; 00445 } 00446 else 00447 { 00448 slope = (Ratehigh - Ratelow)/(zhigh - zlow); 00449 rate = slope*(Charge - zlow) + Ratelow; 00450 } 00451 00453 Hydcs123_v = rate; 00454 00455 /* the routines C6cs123, Ne10cs123, etc... do not resolve L for n>2 */ 00456 /* dividing by N should roughly recover the original l-resolved data */ 00457 Hydcs123_v /= 3.; 00458 } 00459 else 00460 { 00461 /* this branch 1-2s, 1-2p */ 00462 if( nelem == 1 ) 00463 { 00464 /* this brance for helium, then return */ 00465 Hydcs123_v = He2cs123(ipLow,ipHi); 00466 return( Hydcs123_v ); 00467 } 00468 00469 /* electron temperature in eV */ 00470 tev = phycon.te / EVDEGK; 00471 /* energy in eV for hydrogenic species and these quantum numbers */ 00472 EE = ChargeSquared*EVRYD*(1./QuanNLo/QuanNLo - 1./QuanNUp/QuanNUp); 00473 /* EE/kT for this transion */ 00474 q = EE/tev; 00475 TeUse = MIN2(q,10.); 00476 /* q is now EE/kT but between 1 and 10 */ 00477 q = MAX2(1.,TeUse); 00478 expq = exp(q); 00479 00480 /* i must be 0 or 1 */ 00481 ASSERT( i==0 || i==1 ); 00482 C = C0[i] + C1[i]/Charge + C2[i]/ChargeSquared; 00483 D = D0[i] + D1[i]/Charge + D2[i]/ChargeSquared; 00484 00485 /* following code changed so that ee1 always returns e1, 00486 * orifinal version only returned e1 for x < 1 */ 00487 /* use disabled e1: */ 00488 /*if( q < 1. )*/ 00489 /*{*/ 00490 /*rate = (B[i-1] + D*q)*exp(-q) + (A[i-1] + C*q - D*q*q)**/ 00491 /*ee1(q);*/ 00492 /*}*/ 00493 /*else*/ 00494 /*{*/ 00495 /*rate = (B[i-1] + D*q) + (A[i-1] + C*q - D*q*q)*ee1(q);*/ 00496 /*}*/ 00497 /*rate *= 8.010e-8/2./ChargeSquared/tev*sqrt(tev);*/ 00498 /* convert to de-excitation */ 00499 /*if( q < 1. )*/ 00500 /*{*/ 00501 /*rate = rate*exp(q)*gLo/gHi;*/ 00502 /*}*/ 00503 /*else*/ 00504 /*{*/ 00505 /*rate = rate*gLo/gHi;*/ 00506 /*}*/ 00507 00508 /*>>>chng 98 dec 17, ee1 always returns e1 */ 00509 rate = (B[i] + D*q)/expq + (A[i] + C*q - D*q*q)* 00510 ee1(q); 00511 rate *= 8.010e-8/2./ChargeSquared/tev*sqrt(tev); 00512 /* convert to de-excitation */ 00513 rate *= expq*gLo/gHi; 00514 00515 /* convert to cs */ 00516 Hydcs123_v = rate*gHi*phycon.sqrte/COLL_CONST; 00517 } 00518 return( Hydcs123_v ); 00519 } 00520 00521 /*C6cs123 line collision rates for lower levels of hydrogenic carbon, n=1,2,3 */ 00522 STATIC double C6cs123(long int i, 00523 long int j) 00524 { 00525 double C6cs123_v, 00526 TeUse, 00527 t, 00528 x; 00529 static const double a[3] = {-92.23774,-1631.3878,-6326.4947}; 00530 static const double b[3] = {-11.93818,-218.3341,-849.8927}; 00531 static const double c[3] = {0.07762914,1.50127,5.847452}; 00532 static const double d[3] = {78.401154,1404.8475,5457.9291}; 00533 static const double e[3] = {332.9531,5887.4263,22815.211}; 00534 00535 DEBUG_ENTRY( "C6cs123()" ); 00536 00537 /* These are fits to Table 5 of 00538 * >>refer c6 cs Aggarwal, K.M., & Kingston, A.E. 1991, J Phys B, 24, 4583 00539 * C VI collision rates for 1s-3l, 2s-3l, and 2p-3l, 00540 * principal quantum numbers n and l. 00541 * 00542 * i is the lower level and runs from 1 to 3 (1s, 2s, 2p) 00543 * j is the upper level and runs from 2 to 6 (2s, 2p, 3s, 3p, 3d) 00544 * 1s-2s,2p is not done here. 00545 * check temperature: fits only good between 3.8 < log Te < 6.2 00546 */ 00547 /* arrays for fits of 3 transitions see the code below for key: */ 00548 00549 TeUse = MAX2(phycon.te,6310.); 00550 t = MIN2(TeUse,1.6e6); 00551 x = log10(t); 00552 00553 if( i == 1 && j == 2 ) 00554 { 00555 /* 1s - 2s (first entry) */ 00556 fprintf( ioQQQ, " Carbon VI 2s-1s not done in C6cs123\n" ); 00557 cdEXIT(EXIT_FAILURE); 00558 } 00559 00560 else if( i == 1 && j == 3 ) 00561 { 00562 /* 1s - 2p (second entry) */ 00563 fprintf( ioQQQ, " Carbon VI 2p-1s not done in C6cs123\n" ); 00564 cdEXIT(EXIT_FAILURE); 00565 } 00566 00567 else if( i == 1 && ( j == 4 || j == 5 || j == 6 ) ) 00568 { 00569 /* 1s - 3 (first entry) */ 00570 C6cs123_v = a[0] + b[0]*x + c[0]*pow2(x)*sqrt(x) + d[0]*log(x) + 00571 e[0]*log(x)/pow2(x); 00572 } 00573 else if( i == 2 && ( j == 4 || j == 5 || j == 6 ) ) 00574 { 00575 /* 2s - 3 (second entry) */ 00576 C6cs123_v = a[1] + b[1]*x + c[1]*pow2(x)*sqrt(x) + d[1]*log(x) + 00577 e[1]*log(x)/pow2(x); 00578 } 00579 else if( i == 3 && ( j == 4 || j == 5 || j == 6 ) ) 00580 { 00581 /* 2p - 3s (third entry) */ 00582 C6cs123_v = a[2] + b[2]*x + c[2]*pow2(x)*sqrt(x) + d[2]*log(x) + 00583 e[2]*log(x)/pow2(x); 00584 } 00585 else 00586 { 00587 fprintf( ioQQQ, " insane levels for C VI n=1,2,3 !!!\n" ); 00588 cdEXIT(EXIT_FAILURE); 00589 } 00590 return( C6cs123_v ); 00591 } 00592 00593 /*Ca20cs123 line collision rates for lower levels of hydrogenic calcium, n=1,2,3 */ 00594 STATIC double Ca20cs123(long int i, 00595 long int j) 00596 { 00597 double Ca20cs123_v, 00598 TeUse, 00599 t, 00600 x; 00601 static const double a[3] = {-12.5007,-187.2303,-880.18896}; 00602 static const double b[3] = {-1.438749,-22.17467,-103.1259}; 00603 static const double c[3] = {0.008219688,0.1318711,0.6043752}; 00604 static const double d[3] = {10.116516,153.2650,717.4036}; 00605 static const double e[3] = {45.905343,685.7049,3227.2836}; 00606 00607 DEBUG_ENTRY( "Ca20cs123()" ); 00608 00609 /* 00610 * These are fits to Table 5 of 00611 * >>refer ca20 cs Aggarwal, K.M., & Kingston, A.E. 1992, J Phys B, 25, 751 00612 * Ca XX collision rates for 1s-3l, 2s-3l, and 2p-3l, 00613 * principal quantum numbers n and l. 00614 * 00615 * i is the lower level and runs from 1 to 3 (1s, 2s, 2p) 00616 * j is the upper level and runs from 2 to 6 (2s, 2p, 3s, 3p, 3d) 00617 * 1s-2s,2p is not done here. 00618 * check temperature: fits only good between 5.0 < log Te < 7.2 00619 */ 00620 00621 /* arrays for fits of 3 transitions see the code below for key: */ 00622 00623 TeUse = MAX2(phycon.te,1.0e5); 00624 t = MIN2(TeUse,1.585e7); 00625 x = log10(t); 00626 00627 if( i == 1 && j == 2 ) 00628 { 00629 /* 1s - 2s (first entry) */ 00630 fprintf( ioQQQ, " Ca XX 2s-1s not done in Ca20cs123\n" ); 00631 cdEXIT(EXIT_FAILURE); 00632 } 00633 00634 else if( i == 1 && j == 3 ) 00635 { 00636 /* 1s - 2p (second entry) */ 00637 fprintf( ioQQQ, " Ca XX 2p-1s not done in Ca20cs123\n" ); 00638 cdEXIT(EXIT_FAILURE); 00639 } 00640 00641 else if( i == 1 && ( j == 4 || j == 5 || j == 6 )) 00642 { 00643 /* 1s - 3 (first entry) */ 00644 Ca20cs123_v = a[0] + b[0]*x + c[0]*pow2(x)*sqrt(x) + d[0]*log(x) + 00645 e[0]*log(x)/pow2(x); 00646 } 00647 else if( i == 2 && ( j == 4 || j == 5 || j == 6 )) 00648 { 00649 /* 2s - 3 (second entry) */ 00650 Ca20cs123_v = a[1] + b[1]*x + c[1]*pow2(x)*sqrt(x) + d[1]*log(x) + 00651 e[1]*log(x)/pow2(x); 00652 } 00653 else if( i == 3 && ( j == 4 || j == 5 || j == 6 )) 00654 { 00655 /* 2p - 3s (third entry) */ 00656 Ca20cs123_v = a[2] + b[2]*x + c[2]*pow2(x)*sqrt(x) + d[2]*log(x) + 00657 e[2]*log(x)/pow2(x); 00658 } 00659 else 00660 { 00661 fprintf( ioQQQ, " insane levels for Ca XX n=1,2,3 !!!\n" ); 00662 cdEXIT(EXIT_FAILURE); 00663 } 00664 return( Ca20cs123_v ); 00665 } 00666 00667 /*Ne10cs123 line collision rates for lower levels of hydrogenic neon, n=1,2,3 */ 00668 STATIC double Ne10cs123(long int i, 00669 long int j) 00670 { 00671 double Ne10cs123_v, 00672 TeUse, 00673 t, 00674 x; 00675 static const double a[3] = {3.346644,151.2435,71.7095}; 00676 static const double b[3] = {0.5176036,20.05133,13.1543}; 00677 static const double c[3] = {-0.00408072,-0.1311591,-0.1099238}; 00678 static const double d[3] = {-3.064742,-129.8303,-71.0617}; 00679 static const double e[3] = {-11.87587,-541.8599,-241.2520}; 00680 00681 DEBUG_ENTRY( "Ne10cs123()" ); 00682 00683 /* These are fits to Table 5 of 00684 * >>refer ne10 cs Aggarwal, K.M., & Kingston, A.E. 1991, PhyS, 44, 517 00685 * Ne X collision rates for 1s-3, 2s-3l, and 2p-3l, 00686 * principal quantum numbers n and l. 00687 * 00688 * i is the lower level and runs from 1 to 3 (1s, 2s, 2p) 00689 * j is the upper level and runs from 2 to 6 (2s, 2p, 3s, 3p, 3d) 00690 * 1s-2s,2p is not done here. 00691 * check temperature: fits only good between 3.8 < log Te < 6.2 00692 * */ 00693 /* arrays for fits of 3 transitions see the code below for key: */ 00694 00695 TeUse = MAX2(phycon.te,6310.); 00696 t = MIN2(TeUse,1.6e6); 00697 x = log10(t); 00698 00699 if( i == 1 && j == 2 ) 00700 { 00701 /* 1s - 2s (first entry) */ 00702 fprintf( ioQQQ, " Neon X 2s-1s not done in Ne10cs123\n" ); 00703 cdEXIT(EXIT_FAILURE); 00704 } 00705 00706 else if( i == 1 && j == 3 ) 00707 { 00708 /* 1s - 2p (second entry) */ 00709 fprintf( ioQQQ, " Neon X 2p-1s not done in Ne10cs123\n" ); 00710 cdEXIT(EXIT_FAILURE); 00711 } 00712 00713 else if( i == 1 && ( j == 4 || j == 5 || j == 6 ) ) 00714 { 00715 /* 1s - 3 (first entry) */ 00716 Ne10cs123_v = a[0] + b[0]*x + c[0]*pow2(x)*sqrt(x) + d[0]*log(x) + 00717 e[0]*log(x)/pow2(x); 00718 } 00719 else if( i == 2 && ( j == 4 || j == 5 || j == 6 ) ) 00720 { 00721 /* 2s - 3 (second entry) */ 00722 Ne10cs123_v = a[1] + b[1]*x + c[1]*pow2(x)*sqrt(x) + d[1]*log(x) + 00723 e[1]*log(x)/pow2(x); 00724 } 00725 else if( i == 3 && ( j == 4 || j == 5 || j == 6 ) ) 00726 { 00727 /* 2p - 3s (third entry) */ 00728 Ne10cs123_v = a[2] + b[2]*x + c[2]*pow2(x)*sqrt(x) + d[2]*log(x) + 00729 e[2]*log(x)/pow2(x); 00730 } 00731 else 00732 { 00733 fprintf( ioQQQ, " insane levels for Ne X n=1,2,3 !!!\n" ); 00734 cdEXIT(EXIT_FAILURE); 00735 } 00736 return( Ne10cs123_v ); 00737 } 00738 00739 /*He2cs123 line collision strengths for lower levels of helium ion, n=1,2,3, by K Korista */ 00740 STATIC double He2cs123(long int i, 00741 long int j) 00742 { 00743 double He2cs123_v, 00744 t; 00745 static const double a[11]={0.12176209,0.32916723,0.46546497,0.044501688, 00746 0.040523277,0.5234889,1.4903214,1.4215094,1.0295881,4.769306,9.7226127}; 00747 static const double b[11]={0.039936166,2.9711166e-05,-0.020835863,3.0508137e-04, 00748 -2.004485e-15,4.41475e-06,1.0622666e-05,2.0538877e-06,0.80638448,2.0967075e-06, 00749 7.6089851e-05}; 00750 static const double c[11]={143284.77,0.73158545,-2.159172,0.43254802,2.1338557, 00751 8.9899702e-06,-2.9001451e-12,1.762076e-05,52741.735,-2153.1219,-3.3996921e-11}; 00752 00753 DEBUG_ENTRY( "He2cs123()" ); 00754 00755 /* These are fits to Table 2 00756 * >>refer he2 cs Aggarwal, K.M., Callaway, J., Kingston, A.E., Unnikrishnan, K. 00757 * >>refercon 1992, ApJS, 80, 473 00758 * He II collision rates for 1s-2s, 1s-2p, 1s-3s, 1s-3p, 1s-3d, 2s-3s, 2s-3p, 2s-3d, 00759 * 2p-3s, 2p-3p, and 2p-3d. 00760 * principal quantum numbers n and l. 00761 * 00762 * i is the lower level and runs from 1 to 3 (1s, 2s, 2p) 00763 * j is the upper level and runs from 2 to 6 (2s, 2p, 3s, 3p, 3d) 00764 * check temperature: fits only good between 5,000K and 500,000K 00765 * */ 00766 /* array for fits of 11 transitions see the code below for key: */ 00767 00768 t = phycon.te; 00769 if( t < 5000. ) 00770 { 00771 t = 5000.; 00772 } 00773 else if( t > 5.0e05 ) 00774 { 00775 t = 5.0e05; 00776 } 00777 00778 /**************fits begin here************** 00779 * */ 00780 if( i == 1 && j == 2 ) 00781 { 00782 /* 1s - 2s (first entry) */ 00783 He2cs123_v = a[0] + b[0]*exp(-t/c[0]); 00784 } 00785 else if( i == 1 && j == 3 ) 00786 { 00787 /* 1s - 2p (second entry) */ 00788 He2cs123_v = a[1] + b[1]*pow(t,c[1]); 00789 } 00790 else if( i == 1 && j == 4 ) 00791 { 00792 /* 1s - 3s (third entry) */ 00793 He2cs123_v = a[2] + b[2]*log(t) + c[2]/log(t); 00794 } 00795 else if( i == 1 && j == 5 ) 00796 { 00797 /* 1s - 3p (fourth entry) */ 00798 He2cs123_v = a[3] + b[3]*pow(t,c[3]); 00799 } 00800 else if( i == 1 && j == 6 ) 00801 { 00802 /* 1s - 3d (fifth entry) */ 00803 He2cs123_v = a[4] + b[4]*pow(t,c[4]); 00804 } 00805 else if( i == 2 && j == 4 ) 00806 { 00807 /* 2s - 3s (sixth entry) */ 00808 He2cs123_v = (a[5] + c[5]*t)/(1 + b[5]*t); 00809 } 00810 else if( i == 2 && j == 5 ) 00811 { 00812 /* 2s - 3p (seventh entry) */ 00813 He2cs123_v = a[6] + b[6]*t + c[6]*t*t; 00814 } 00815 else if( i == 2 && j == 6 ) 00816 { 00817 /* 2s - 3d (eighth entry) */ 00818 He2cs123_v = (a[7] + c[7]*t)/(1 + b[7]*t); 00819 } 00820 else if( i == 3 && j == 4 ) 00821 { 00822 /* 2p - 3s (ninth entry) */ 00823 He2cs123_v = a[8] + b[8]*exp(-t/c[8]); 00824 } 00825 else if( i == 3 && j == 5 ) 00826 { 00827 /* 2p - 3p (tenth entry) */ 00828 He2cs123_v = a[9] + b[9]*t + c[9]/t; 00829 } 00830 else if( i == 3 && j == 6 ) 00831 { 00832 /* 2p - 3d (eleventh entry) */ 00833 He2cs123_v = a[10] + b[10]*t + c[10]*t*t; 00834 } 00835 else 00836 { 00837 /**************fits end here************** */ 00838 fprintf( ioQQQ, " insane levels for He II n=1,2,3 !!!\n" ); 00839 cdEXIT(EXIT_FAILURE); 00840 } 00841 return( He2cs123_v ); 00842 } 00843 00844 /*Fe26cs123 line collision rates for lower levels of hydrogenic iron, n=1,2,3 */ 00845 STATIC double Fe26cs123(long int i, 00846 long int j) 00847 { 00848 double Fe26cs123_v, 00849 TeUse, 00850 t, 00851 x; 00852 static const double a[3] = {-4.238398,-238.2599,-1211.5237}; 00853 static const double b[3] = {-0.4448177,-27.06869,-136.7659}; 00854 static const double c[3] = {0.0022861,0.153273,0.7677703}; 00855 static const double d[3] = {3.303775,191.7165,972.3731}; 00856 static const double e[3] = {15.82689,878.1333,4468.696}; 00857 00858 DEBUG_ENTRY( "Fe26cs123()" ); 00859 00860 /* These are fits to Table 5 of 00861 * >>refer fe26 cs Aggarwal, K.M., & Kingston, A.E. 1993, ApJS, 85, 187 00862 * Fe XXVI collision rates for 1s-3, 2s-3, and 2p-3, 00863 * principal quantum numbers n and l. 00864 * 00865 * i is the lower level and runs from 1 to 3 (1s, 2s, 2p) 00866 * j is the upper level and runs from 2 to 6 (2s, 2p, 3s, 3p, 3d) 00867 * 1s-2s,2p is not done here. 00868 * check temperature: fits only good between 5.2 < log Te < 7.2 00869 * */ 00870 /* arrays for fits of 3 transitions see the code below for key: */ 00871 00872 TeUse = MAX2(phycon.te,1.585e5); 00873 t = MIN2(TeUse,1.585e7); 00874 x = log10(t); 00875 00876 if( i == 1 && j == 2 ) 00877 { 00878 /* 1s - 2s (first entry) */ 00879 fprintf( ioQQQ, " Fe XXVI 2s-1s not done in Fe26cs123\n" ); 00880 cdEXIT(EXIT_FAILURE); 00881 } 00882 00883 else if( i == 1 && j == 3 ) 00884 { 00885 /* 1s - 2p (second entry) */ 00886 fprintf( ioQQQ, " Fe XXVI 2p-1s not done in Fe26cs123\n" ); 00887 cdEXIT(EXIT_FAILURE); 00888 } 00889 00890 else if( i == 1 && ( j == 4 || j == 5 || j == 6 ) ) 00891 { 00892 /* 1s - 3 (first entry) */ 00893 Fe26cs123_v = a[0] + b[0]*x + c[0]*pow2(x)*sqrt(x) + d[0]*log(x) + 00894 e[0]*log(x)/pow2(x); 00895 } 00896 else if( i == 2 && ( j == 4 || j == 5 || j == 6 ) ) 00897 { 00898 /* 2s - 3 (second entry) */ 00899 Fe26cs123_v = a[1] + b[1]*x + c[1]*pow2(x)*sqrt(x) + d[1]*log(x) + 00900 e[1]*log(x)/pow2(x); 00901 } 00902 else if( i == 3 && ( j == 4 || j == 5 || j == 6 ) ) 00903 { 00904 /* 2p - 3s (third entry) */ 00905 Fe26cs123_v = a[2] + b[2]*x + c[2]*pow2(x)*sqrt(x) + d[2]*log(x) + 00906 e[2]*log(x)/pow2(x); 00907 } 00908 else 00909 { 00910 fprintf( ioQQQ, " insane levels for Ca XX n=1,2,3 !!!\n" ); 00911 cdEXIT(EXIT_FAILURE); 00912 } 00913 return( Fe26cs123_v ); 00914 } 00915 00916 00917 STATIC double CS_ThermAve_PR78(long ipISO, long nelem, long nHi, long nLo, double deltaE, double temp ) 00918 { 00919 00920 double coll_str; 00921 00922 DEBUG_ENTRY( "CS_ThermAve_PR78()" ); 00923 00924 global_ipISO = ipISO; 00925 global_nelem = nelem; 00926 global_nHi = nHi; 00927 global_nLo = nLo; 00928 global_deltaE = deltaE; 00929 00930 kTRyd = temp / TE1RYD; 00931 00932 if( !iso.lgCS_therm_ave[ipISO] ) 00933 { 00934 /* Must do some thermal averaging for densities greater 00935 * than about 10000 and less than about 1e10, 00936 * because kT gives significantly different results. 00937 * Still, do sparser integration than is done below */ 00938 if( (dense.eden > 10000.) && (dense.eden < 1E10 ) ) 00939 { 00940 coll_str = qg32( 0.0, 6.0, Therm_ave_coll_str_int_PR78); 00941 } 00942 else 00943 { 00944 /* Do NOT average over Maxwellian */ 00945 coll_str = CS_PercivalRichards78( kTRyd ); 00946 } 00947 } 00948 else 00949 { 00950 /* DO average over Maxwellian */ 00951 coll_str = qg32( 0.0, 1.0, Therm_ave_coll_str_int_PR78); 00952 coll_str += qg32( 1.0, 10.0, Therm_ave_coll_str_int_PR78); 00953 } 00954 00955 return coll_str; 00956 } 00957 00958 /* The integrand for calculating the thermal average of collision strengths */ 00959 STATIC double Therm_ave_coll_str_int_PR78( double EOverKT ) 00960 { 00961 double integrand; 00962 00963 DEBUG_ENTRY( "Therm_ave_coll_str_int_PR78()" ); 00964 00965 integrand = exp( -1.*EOverKT ) * CS_PercivalRichards78( EOverKT * kTRyd ); 00966 00967 return integrand; 00968 } 00969 00970 STATIC double CS_PercivalRichards78( double Ebar ) 00971 { 00972 double cross_section, coll_str; 00973 double stat_weight; 00974 double A, D, L, F, G, H; 00975 double np, n, s, Z_3, xPlus, xMinus, y; 00976 long ipISO, nelem; 00977 00978 DEBUG_ENTRY( "CS_PercivalRichards78()" ); 00979 00980 if( Ebar < global_deltaE ) 00981 { 00982 DEBUG_ENTRY( "CS_PercivalRichards78()" ); 00983 return 0.; 00984 } 00985 00986 ipISO = global_ipISO; 00987 nelem = global_nelem; 00988 np = (double)global_nHi; 00989 n = (double)global_nLo; 00990 s = np - n; 00991 ASSERT( s > 0. ); 00992 00993 A = (8./3./s) * pow(np/s/n, 3.) * (0.184 - 0.04 * pow( s, -0.66)) * pow( 1. - 0.2*s/n/np, 1.+2.*s); 00994 00995 Z_3 = (double)(nelem + 1. - ipISO); 00996 00997 D = exp( - Z_3 * Z_3 / n / np / Ebar / Ebar ); 00998 00999 L = log( (1. + 0.53 * Ebar * Ebar * n * np / Z_3 / Z_3) / (1. + 0.4*Ebar) ); 01000 01001 F = pow( 1. - 0.3 * s * D / n /np, 1. + 2.*s ); 01002 01003 G = 0.5* POW3( Ebar * n * n / Z_3 / np ); 01004 01005 xPlus = 2. * Z_3 / ( n * n * Ebar * ( sqrt( 2. - n*n/np/np ) + 1. ) ); 01006 xMinus = 2. * Z_3 / ( n * n * Ebar * ( sqrt( 2. - n*n/np/np ) - 1. ) ); 01007 01008 y = 1. / (1. - D * log ( 18. * s )/ 4. / s); 01009 01010 H = POW2( xMinus) * log( 1. + 2.*xMinus/3. ) / ( 2.*y + 1.5*xMinus ); 01011 H -= POW2( xPlus ) * log( 1. + 2.* xPlus/3. ) / ( 2.*y + 1.5*xPlus ); 01012 01013 /* this is the LHS of equation 1 of PR78 */ 01014 cross_section = (A*D*L + F*G*H); 01015 /* this is the result after solving equation 1 for the cross section */ 01016 cross_section *= PI * POW2( n * n * BOHR_RADIUS_CM / Z_3 ) / Ebar; 01017 01018 if( ipISO == ipH_LIKE ) 01019 stat_weight = 2. * n * n; 01020 else if( ipISO == ipHE_LIKE ) 01021 stat_weight = 4. * n * n; 01022 else 01023 TotalInsanity(); 01024 01025 /* convert to collision strength */ 01026 coll_str = cross_section * stat_weight * Ebar / ( PI * POW2( BOHR_RADIUS_CM ) ); 01027 return coll_str; 01028 } 01029 01030 #if 0 01031 STATIC void TestPercivalRichards( void ) 01032 { 01033 double CStemp; 01034 01035 /* this reproduces Table 1 of PR78 */ 01036 for( long i=0; i<5; i++ ) 01037 { 01038 double Ebar[5] = {0.1, 0.4, 0.8, 1.0, 10.}; 01039 01040 CStemp = CS_PercivalRichards78( 0, 2, 12, 10, Ebar[i] ); 01041 } 01042 01043 /* this reproduces Table 2 of PR78 */ 01044 for( long i=0; i<5; i++ ) 01045 { 01046 double Ebar[5] = {0.1, 0.4, 0.8, 1.0, 10.}; 01047 01048 CStemp = CS_ThermAve_PR78( ipISO, 0, N_(ipHi), N_(ipLo), phycon.te ); 01049 } 01050 01051 return; 01052 } 01053 #endif 01054 01055 realnum HydroCSInterp(long int nelem, 01056 long int ipHi, 01057 long int ipLo, 01058 long int ipCollider ) 01059 { 01060 double CStemp; 01061 long ipISO = ipH_LIKE; 01062 01063 DEBUG_ENTRY( "HydroCSInterp()" ); 01064 01065 /* This set of collision strengths should only be used 01066 * if the Storey and Hummer flag is set */ 01067 if( opac.lgCaseB_HummerStorey ) 01068 { 01069 if( N_(ipLo) == N_(ipHi) ) 01070 { 01071 if( N_(ipHi) <= iso.n_HighestResolved_max[ipH_LIKE][nelem] && 01072 abs( L_(ipLo) - L_(ipHi) ) != 1 ) 01073 { 01074 /* if delta L is not +/- 1, set collision strength to zero. */ 01075 CStemp = 0.; 01076 } 01077 else if( ipCollider != ipELECTRON ) 01078 { 01079 CStemp = CS_l_mixing_PS64( ipH_LIKE, nelem, ipLo, ipHi, ipCollider ); 01080 } 01081 else 01082 CStemp = 0.; 01083 } 01084 01085 else 01086 { 01087 if( N_(ipHi) <= iso.n_HighestResolved_max[ipH_LIKE][nelem] && 01088 abs( L_(ipLo) - L_(ipHi) ) != 1 ) 01089 { 01090 /* if delta L is not +/- 1, set collision strength to zero. */ 01091 CStemp = 0.; 01092 } 01093 else if( ipCollider == ipELECTRON ) 01094 { 01095 CStemp = CS_ThermAve_PR78( ipH_LIKE, nelem, N_(ipHi), N_(ipLo), 01096 Transitions[ipH_LIKE][nelem][ipHi][ipLo].EnergyErg / EN1RYD , phycon.te ); 01097 01098 if( N_(ipHi) <= iso.n_HighestResolved_max[ipH_LIKE][nelem] ) 01099 { 01100 CStemp /= 2.; 01101 } 01102 } 01103 else 01104 CStemp = 0.; 01105 } 01106 } 01107 else 01108 { 01109 /* HCSAR_interp interpolates on a table to return R-matrix collision strengths 01110 * for hydrogen only */ 01111 if( nelem==ipHYDROGEN && ipCollider==ipELECTRON && N_(ipHi) <= 5 && ( N_(ipHi) != N_(ipLo) ) ) 01112 { 01113 CStemp = HCSAR_interp(ipLo,ipHi); 01114 } 01115 else if( nelem==ipHYDROGEN && ipCollider==ipELECTRON && ( N_(ipHi) != N_(ipLo) ) ) 01116 { 01117 CStemp = hydro_vs_deexcit( ipH_LIKE, nelem, ipHi, ipLo, Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->Aul ); 01118 } 01119 else if( nelem>ipHYDROGEN && ipCollider==ipELECTRON && N_(ipHi) <= 3 && N_(ipLo) < 3 ) 01120 { 01121 CStemp = Hydcs123(ipLo + 1,ipHi + 1,nelem,'e'); 01122 } 01123 else if( nelem>ipHYDROGEN && ipCollider==ipPROTON && ipHi==ipH2p && ipLo==ipH2s ) 01124 { 01125 CStemp = Hydcs123(ipLo + 1,ipHi + 1,nelem,'p'); 01126 } 01127 else if( N_(ipLo) == N_(ipHi) ) 01128 { 01129 if( iso.lgCS_Vrinceanu[ipH_LIKE] ) 01130 { 01131 CStemp = CS_l_mixing_VF01(ipH_LIKE, nelem, 01132 StatesElem[ipH_LIKE][nelem][ipLo].n, 01133 StatesElem[ipH_LIKE][nelem][ipLo].l, 01134 StatesElem[ipH_LIKE][nelem][ipHi].l, 01135 StatesElem[ipH_LIKE][nelem][ipLo].S, 01136 phycon.te, 01137 ipCollider ); 01138 } 01139 else 01140 CStemp = CS_l_mixing_PS64( ipH_LIKE, nelem, ipLo, ipHi, ipCollider ); 01141 } 01142 else 01143 { 01144 ASSERT( N_(ipHi) != N_(ipLo) ); 01145 /* highly excited levels */ 01146 CStemp = CS_VS80( ipH_LIKE, nelem, ipHi, ipLo, 01147 Transitions[ipH_LIKE][nelem][ipHi][ipLo].Emis->Aul, 01148 phycon.te, 01149 ipCollider ); 01150 } 01151 } 01152 01153 return (realnum)CStemp; 01154 }