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 /*atmdat_3body derive three-body recombination coefficients */ 00004 /*da interpolate on three body recombination by Steve Cota */ 00005 #include "cddefines.h" 00006 #include "ionbal.h" 00007 #include "phycon.h" 00008 #include "dense.h" 00009 #include "trace.h" 00010 #include "punch.h" 00011 #include "atmdat.h" 00012 00013 #define MAXZ 28 00014 00015 STATIC void blkdata1(void); 00016 STATIC double da(double z); 00017 00018 static double a2[63], 00019 b2[63], 00020 x2[63]; 00021 00022 static double a0[83], 00023 x0[83]; 00024 static realnum b0[83], 00025 b1[83]; 00026 00027 static double a1[83], 00028 x1[83]; 00029 00030 static double tz[83], 00031 zlog7[28], 00032 zlog2[28]; 00033 00034 #define RC_INI(rs) (rs[_r].rc>1 ? DEC_RC_(rs) : (rs[_r].rc==1 ? INC_NDX_(rs) : rs[_r].ini )) 00035 #define DEC_RC_(rs) (rs[_r].rc--,rs[_r].ini) 00036 #define INC_NDX_(rs) (_r++,rs[_r-1].ini) 00037 00038 /* this "mapping function" occurs below */ 00039 STATIC double xmap(double x[], 00040 double y[], 00041 double x0); 00042 00043 /* inverse routine, also below */ 00044 STATIC double xinvrs(double y, 00045 double a, 00046 double b, 00047 double u, 00048 double v, 00049 long int *ifail); 00050 00051 /* =================================================================== */ 00052 void atmdat_3body(void) 00053 { 00054 long int i, 00055 iup; 00056 00057 DEBUG_ENTRY( "atmdat_3body()" ); 00058 00059 00060 if( ionbal.lgNoCota ) 00061 { 00062 for( i=0; i < LIMELM; i++ ) 00063 { 00064 ionbal.CotaRate[i] = 0.; 00065 } 00066 atmdat.nsbig = 0; 00067 return; 00068 } 00069 00070 if( atmdat.nsbig == 0 ) 00071 { 00072 /* steve cota only defined things up to 28 */ 00073 iup = MIN2(28,LIMELM); 00074 } 00075 else 00076 { 00077 iup = MIN3( LIMELM , atmdat.nsbig , 28 ); 00078 } 00079 00080 for( i=0; i < iup; i++ ) 00081 { 00082 ionbal.CotaRate[i] = (realnum)da((double)(i+1)); 00083 } 00084 00085 atmdat.nsbig = 0; 00086 00087 if( trace.lgTrace && trace.lgTrace3Bod ) 00088 { 00089 fprintf( ioQQQ, " 3BOD rate:" ); 00090 for( i=1; i <= 14; i++ ) 00091 { 00092 fprintf( ioQQQ, "%8.1e", ionbal.CotaRate[i-1] ); 00093 } 00094 fprintf( ioQQQ, "\n" ); 00095 } 00096 00097 if( punch.lgioRecom ) 00098 { 00099 /* option to punch coefficients */ 00100 fprintf( punch.ioRecom, " 3-body rec coef vs charge \n" ); 00101 for( i=0; i < iup; i++ ) 00102 { 00103 fprintf( punch.ioRecom, "%3ld%10.2e\n", i+1, ionbal.CotaRate[i] ); 00104 } 00105 fprintf( punch.ioRecom, "\n"); 00106 } 00107 return; 00108 } 00109 00110 /* =================================================================== */ 00111 STATIC double da(double z) 00112 { 00113 /*lint -e736 loss of precision in assignment in translated data */ 00114 long int jfail, 00115 nt, 00116 nt0, 00117 nt1, 00118 nz; 00119 00120 static bool lgCalled=false; 00121 double a, 00122 alogn, 00123 alognc, 00124 alogt, 00125 b, 00126 c, 00127 d, 00128 da_v, 00129 expp, 00130 u, 00131 v, 00132 x, 00133 xnc, 00134 y, 00135 zlog; 00136 double yya[3], 00137 xx[3], 00138 yyb[3], 00139 yyx[3], 00140 yyy[3]; 00141 00142 /* alogte is base 10 log of temperature and elec density*/ 00143 double alogte , alogne; 00144 00145 DEBUG_ENTRY( "da()" ); 00146 00147 /* WRITTEN BY S. A. COTA, 2/87 00148 * */ 00149 00150 /* MAXZ IS THE MAXIMUM EFFECTIVE NUCLEAR CHARGE ( = IONIC CHARGE + 1 ) 00151 * WHICH THE DIMENSION STATEMENTS ACCOMODATE. 00152 * 00153 * IT IS USED ONLY FOR THE ARRAY ZLOG7 ( = 7 * LOG ( Z ) ) 00154 * AND THE ARRAY ZLOG2 ( = 2 * LOG ( Z ) ) . THESE ARRAYS 00155 * CONTAIN EASILY CALCULATED VALUES, WHICH HAVE BEEN STORED 00156 * TO SAVE TIME IN EXECUTION. 00157 * 00158 * IF MAXZ IS EXCEEDED, THIS PROGRAM SIMPLY CALCULATES THE 00159 * LOGS INSTEAD OF LOOKING THEM UP. 00160 * 00161 * */ 00162 00163 if( !lgCalled ) 00164 { 00165 lgCalled = true; 00166 blkdata1(); 00167 } 00168 00169 /*begin sanity check */ 00170 ASSERT( zlog7[1] > 0. ); 00171 00172 nz = (long)(z + .1); 00173 alogne = log10(dense.eden); 00174 alogte = log10(phycon.te); 00175 if( nz > MAXZ ) 00176 { 00177 zlog = log10(z); 00178 alogt = alogte - 2.*zlog; 00179 alogn = alogne - 7.*zlog; 00180 } 00181 else 00182 { 00183 alogt = alogte - zlog2[nz-1]; 00184 alogn = alogne - zlog7[nz-1]; 00185 } 00186 00187 /* CHECK IF PARAMETERS ARE WITHIN BOUNDS. IF NOT, INCREMENT 00188 * APPROPRIATE ERROR COUNTER AND SET TO BOUNDARY IF 00189 * NECESSARY: 00190 * 00191 * DEFINITION OF ERROR COUNTERS: 00192 * 00193 * ILT : LOW T 00194 * ILTLN : LOW T , LOW N 00195 * ILTHN : LOW T , HIGH N 00196 * IHTHN : HIGH T , HIGH N 00197 * */ 00198 if( alogt < 0. ) 00199 { 00200 ionbal.ilt += 1; 00201 alogt = 0.; 00202 } 00203 00204 if( alogt <= 2.1760913 ) 00205 { 00206 if( alogn < (3.5*alogt - 8.) ) 00207 { 00208 ionbal.iltln += 1; 00209 } 00210 else if( alogn > (3.5*alogt - 2.) ) 00211 { 00212 ionbal.ilthn += 1; 00213 alogn = 3.5*alogt - 2.; 00214 } 00215 00216 } 00217 else if( alogt <= 2.4771213 ) 00218 { 00219 if( alogn > 9. ) 00220 { 00221 ionbal.ilthn += 1; 00222 alogn = 9.; 00223 } 00224 00225 } 00226 else if( alogt <= 5.1139434 ) 00227 { 00228 if( alogn > 13. ) 00229 { 00230 ionbal.ihthn += 1; 00231 alogn = 13.; 00232 } 00233 00234 } 00235 else 00236 { 00237 da_v = 0.; 00238 return( da_v ); 00239 } 00240 00241 /* LOCATE POSITION IN ARRAYS */ 00242 if( alogt <= 2. ) 00243 { 00244 nt = (long)(9.9657843*alogt + 1.); 00245 } 00246 else 00247 { 00248 nt = (long)(19.931568*alogt - 19.); 00249 } 00250 nt = MIN2(83,nt); 00251 nt = MAX2(1,nt); 00252 00253 /* CENTER UP SINCE ARRAY VALUES ARE ROUNDED */ 00254 if( fabs(alogt-tz[nt-1]) >= fabs(alogt-tz[MIN2(83,nt+1)-1]) ) 00255 { 00256 nt = MIN2(83,nt+1); 00257 } 00258 else if( fabs(alogt-tz[nt-1]) > fabs(alogt-tz[MAX2(1,nt-1)-1]) ) 00259 { 00260 nt = MAX2(1,nt-1); 00261 } 00262 00263 /* CHECK IF INTERPOLATION IS NEEDED AND PROCEED IF NOT.*/ 00264 if( fabs(alogt-tz[nt-1]) < 0.00001 ) 00265 { 00266 if( z != 1.0 ) 00267 { 00268 c = a1[nt-1]; 00269 d = b1[nt-1]; 00270 u = x1[nt-1]; 00271 v = 8.90; 00272 } 00273 else 00274 { 00275 nt = MAX2(21,nt); 00276 nt = MIN2(83,nt); 00277 c = a2[nt-(21)]; 00278 d = b2[nt-(21)]; 00279 u = x2[nt-(21)]; 00280 v = 9.40; 00281 } 00282 00283 xnc = xinvrs(alogn,c,d,u,v,&jfail); 00284 if( xnc <= 0. || jfail != 0 ) 00285 { 00286 ionbal.ifail = 1; 00287 jfail = 1; 00288 da_v = 0.; 00289 return( da_v ); 00290 } 00291 alognc = log10(xnc); 00292 00293 a = a0[nt-1]; 00294 b = b0[nt-1]; 00295 x = -2.45; 00296 y = x0[nt-1]; 00297 nt0 = nt - 1; 00298 00299 /* IF INTERPOLATION WAS REQUIRED, 00300 * CHECK THAT NT IS NOT ON THE EDGE OF A DISCONTINUITY, 00301 * EITHER AT END OF ARRAYS OR WITHIN THEM, 00302 * WHERE VALUES CHANGE ABRUPTLY. 00303 * */ 00304 } 00305 else 00306 { 00307 if( (nt <= 21) && (z == 1.0) ) 00308 { 00309 nt = 22; 00310 } 00311 else if( nt <= 1 ) 00312 { 00313 nt = 2; 00314 } 00315 else if( nt >= 83 ) 00316 { 00317 nt = 82; 00318 } 00319 else if( nt == 24 ) 00320 { 00321 if( alogt <= 2.1760913 ) 00322 { 00323 nt = 23; 00324 } 00325 else 00326 { 00327 nt = 25; 00328 } 00329 } 00330 else if( nt == 30 ) 00331 { 00332 if( alogt <= 2.4771213 ) 00333 { 00334 nt = 29; 00335 } 00336 else 00337 { 00338 nt = 31; 00339 } 00340 } 00341 00342 nt0 = nt - 1; 00343 nt1 = nt + 1; 00344 xx[0] = tz[nt0-1]; 00345 xx[1] = tz[nt-1]; 00346 xx[2] = tz[nt1-1]; 00347 00348 if( z != 1.0 ) 00349 { 00350 if( nt0 == 24 ) 00351 { 00352 yya[0] = 17.2880135; 00353 yyb[0] = 6.93410742e03; 00354 yyx[0] = -3.75; 00355 } 00356 else if( nt0 == 30 ) 00357 { 00358 yya[0] = 17.4317989; 00359 yyb[0] = 1.36653821e03; 00360 yyx[0] = -3.40; 00361 } 00362 else 00363 { 00364 yya[0] = a1[nt0-1]; 00365 yyb[0] = b1[nt0-1]; 00366 yyx[0] = x1[nt0-1]; 00367 } 00368 00369 yya[1] = a1[nt-1]; 00370 yya[2] = a1[nt1-1]; 00371 c = xmap(xx,yya,alogt); 00372 yyb[1] = b1[nt-1]; 00373 yyb[2] = b1[nt1-1]; 00374 d = xmap(xx,yyb,alogt); 00375 yyx[1] = x1[nt-1]; 00376 yyx[2] = x1[nt1-1]; 00377 u = xmap(xx,yyx,alogt); 00378 v = 8.90; 00379 00380 } 00381 else 00382 { 00383 if( nt0 == 24 ) 00384 { 00385 yya[0] = 20.1895161; 00386 yyb[0] = 2.25774918e01; 00387 yyx[0] = -1.66; 00388 } 00389 else if( nt0 == 30 ) 00390 { 00391 yya[0] = 19.8647804; 00392 yyb[0] = 6.70408707e02; 00393 yyx[0] = -2.12; 00394 } 00395 else 00396 { 00397 yya[0] = a2[nt0-(21)]; 00398 yyb[0] = b2[nt0-(21)]; 00399 yyx[0] = x2[nt0-(21)]; 00400 } 00401 00402 yya[1] = a2[nt-(21)]; 00403 yya[2] = a2[nt1-(21)]; 00404 c = xmap(xx,yya,alogt); 00405 yyb[1] = b2[nt-(21)]; 00406 yyb[2] = b2[nt1-(21)]; 00407 d = xmap(xx,yyb,alogt); 00408 yyx[1] = x2[nt-(21)]; 00409 yyx[2] = x2[nt1-(21)]; 00410 u = xmap(xx,yyx,alogt); 00411 v = 9.40; 00412 } 00413 00414 xnc = xinvrs(alogn,c,d,u,v,&jfail); 00415 if( xnc <= 0. || jfail != 0 ) 00416 { 00417 ionbal.ifail = 1; 00418 jfail = 1; 00419 da_v = 0.; 00420 return( da_v ); 00421 } 00422 alognc = log10(xnc); 00423 00424 if( nt0 == 24 ) 00425 { 00426 yya[0] = -8.04963875; 00427 yyb[0] = 1.07205127e03; 00428 yyy[0] = 2.05; 00429 } 00430 else if( nt0 == 30 ) 00431 { 00432 yya[0] = -8.54721069; 00433 yyb[0] = 4.70450195e02; 00434 yyy[0] = 2.05; 00435 } 00436 else 00437 { 00438 yya[0] = a0[nt0-1]; 00439 yyb[0] = b0[nt0-1]; 00440 yyy[0] = x0[nt0-1]; 00441 } 00442 00443 yya[1] = a0[nt-1]; 00444 yya[2] = a0[nt1-1]; 00445 a = xmap(xx,yya,alogt); 00446 yyb[1] = b0[nt-1]; 00447 yyb[2] = b0[nt1-1]; 00448 b = xmap(xx,yyb,alogt); 00449 x = -2.45; 00450 yyy[1] = x0[nt-1]; 00451 yyy[2] = x0[nt1-1]; 00452 y = xmap(xx,yyy,alogt); 00453 } 00454 00455 expp = a - y*alognc + b*pow(xnc,x); 00456 if( expp < 37 ) 00457 { 00458 da_v = z*pow(10.,expp); 00459 } 00460 else 00461 { 00462 da_v = 0.; 00463 } 00464 ionbal.ifail += jfail; 00465 00466 return( da_v ); 00467 } 00468 00469 /****************************************************************************** */ 00470 STATIC void blkdata1(void) 00471 { 00472 /*block data with Steve Cota's 3-body recombination coefficients */ 00473 00474 /* data for function da. 00475 * 00476 * S. A. COTA, 2/1987 00477 * */ 00478 00479 long int _i, 00480 _r; 00481 realnum *const ba0 = (realnum*)b0; 00482 realnum *const ba1 = (realnum*)b1; 00483 realnum *const bb0 = (realnum*)((char*)(b0 + 79)); 00484 realnum *const bb1 = (realnum*)((char*)(b1 + 79)); 00485 00486 /* to fix all the conversion errors, change realnum ini to double ini, 00487 * but chech that results still ok */ 00488 { static struct{ long rc; double ini; } _rs0[] = { 00489 {1, 0.00000e00}, 00490 {1, 2.10721e00}, 00491 {1, 3.33985e00}, 00492 {1, 4.21442e00}, 00493 {1, 4.89279e00}, 00494 {1, 5.44706e00}, 00495 {1, 5.91569e00}, 00496 {1, 6.32163e00}, 00497 {1, 6.67970e00}, 00498 {1, 7.00000e00}, 00499 {1, 7.28975e00}, 00500 {1, 7.55427e00}, 00501 {1, 7.79760e00}, 00502 {1, 8.02290e00}, 00503 {1, 8.23264e00}, 00504 {1, 8.42884e00}, 00505 {1, 8.61314e00}, 00506 {1, 8.78691e00}, 00507 {1, 8.95128e00}, 00508 {1, 9.10721e00}, 00509 {1, 9.25554e00}, 00510 {1, 9.39696e00}, 00511 {1, 9.53209e00}, 00512 {1, 9.66148e00}, 00513 {1, 9.78558e00}, 00514 {1, 9.90481e00}, 00515 {1, 10.01954e00}, 00516 {1, 10.13010e00}, 00517 {0L, 0} 00518 }; 00519 for(_i=_r=0L; _i < 28; _i++) 00520 zlog7[_i] = RC_INI(_rs0); } 00521 { static struct{ long rc; double ini; } _rs1[] = { 00522 {1, 0.00000e00}, 00523 {1, 6.02060e-01}, 00524 {1, 9.54243e-01}, 00525 {1, 1.20412e00}, 00526 {1, 1.39794e00}, 00527 {1, 1.55630e00}, 00528 {1, 1.69020e00}, 00529 {1, 1.80618e00}, 00530 {1, 1.90849e00}, 00531 {1, 2.00000e00}, 00532 {1, 2.08279e00}, 00533 {1, 2.15836e00}, 00534 {1, 2.22789e00}, 00535 {1, 2.29226e00}, 00536 {1, 2.35218e00}, 00537 {1, 2.40824e00}, 00538 {1, 2.46090e00}, 00539 {1, 2.51055e00}, 00540 {1, 2.55751e00}, 00541 {1, 2.60206e00}, 00542 {1, 2.64444e00}, 00543 {1, 2.68485e00}, 00544 {1, 2.72346e00}, 00545 {1, 2.76042e00}, 00546 {1, 2.79588e00}, 00547 {1, 2.82995e00}, 00548 {1, 2.86272e00}, 00549 {1, 2.89431e00}, 00550 {0L, 0} 00551 }; 00552 for(_i=_r=0L; _i < 28; _i++) 00553 zlog2[_i] = RC_INI(_rs1); } 00554 { static struct{ long rc; double ini; } _rs2[] = { 00555 {1, 0.}, 00556 {1, 0.09691}, 00557 {1, 0.17609}, 00558 {1, 0.30103}, 00559 {1, 0.39794}, 00560 {1, 0.47712}, 00561 {1, 0.60206}, 00562 {1, 0.69897}, 00563 {1, 0.77815}, 00564 {1, 0.90309}, 00565 {1, 1.00000}, 00566 {1, 1.07918}, 00567 {1, 1.20412}, 00568 {1, 1.30103}, 00569 {1, 1.39794}, 00570 {1, 1.47712}, 00571 {1, 1.60206}, 00572 {1, 1.69897}, 00573 {1, 1.77815}, 00574 {1, 1.90309}, 00575 {1, 2.00000}, 00576 {1, 2.06070}, 00577 {1, 2.09691}, 00578 {1, 2.17609}, 00579 {1, 2.20412}, 00580 {1, 2.24304}, 00581 {1, 2.30103}, 00582 {1, 2.35218}, 00583 {1, 2.39794}, 00584 {1, 2.47712}, 00585 {1, 2.51188}, 00586 {1, 2.54407}, 00587 {1, 2.60206}, 00588 {1, 2.65321}, 00589 {1, 2.69897}, 00590 {1, 2.75967}, 00591 {1, 2.81291}, 00592 {1, 2.86034}, 00593 {1, 2.91645}, 00594 {1, 2.95424}, 00595 {1, 3.00000}, 00596 {1, 3.07918}, 00597 {1, 3.11394}, 00598 {1, 3.17609}, 00599 {1, 3.20412}, 00600 {1, 3.25527}, 00601 {1, 3.30103}, 00602 {1, 3.36173}, 00603 {1, 3.39794}, 00604 {1, 3.46240}, 00605 {1, 3.51188}, 00606 {1, 3.56820}, 00607 {1, 3.60206}, 00608 {1, 3.66276}, 00609 {1, 3.72016}, 00610 {1, 3.76343}, 00611 {1, 3.81291}, 00612 {1, 3.86034}, 00613 {1, 3.90309}, 00614 {1, 3.95424}, 00615 {1, 4.02119}, 00616 {1, 4.06070}, 00617 {1, 4.11394}, 00618 {1, 4.16137}, 00619 {1, 4.20412}, 00620 {1, 4.25527}, 00621 {1, 4.31175}, 00622 {1, 4.36173}, 00623 {1, 4.41497}, 00624 {1, 4.46240}, 00625 {1, 4.51521}, 00626 {1, 4.56526}, 00627 {1, 4.61542}, 00628 {1, 4.66605}, 00629 {1, 4.71600}, 00630 {1, 4.76343}, 00631 {1, 4.81624}, 00632 {1, 4.86629}, 00633 {1, 4.91645}, 00634 {1, 4.96614}, 00635 {1, 5.02119}, 00636 {1, 5.06726}, 00637 {1, 5.11394}, 00638 {0L, 0 } 00639 }; 00640 for(_i=_r=0L; _i < 83; _i++) 00641 tz[_i] = RC_INI(_rs2); } 00642 { static struct{ long rc; double ini; } _rs3[] = { 00643 {1, -4.31396484}, 00644 {1, -4.56640625}, 00645 {1, -4.74560547}, 00646 {1, -4.98535156}, 00647 {1, -5.15373850}, 00648 {1, -5.28123093}, 00649 {1, -5.48215008}, 00650 {1, -5.63811255}, 00651 {1, -5.76573515}, 00652 {1, -5.96755028}, 00653 {1, -6.12449837}, 00654 {1, -6.25304174}, 00655 {1, -6.45615673}, 00656 {1, -6.61384058}, 00657 {1, -6.77161551}, 00658 {1, -6.90069818}, 00659 {1, -7.10470295}, 00660 {1, -7.26322412}, 00661 {1, -7.39289951}, 00662 {1, -7.59792519}, 00663 {1, -7.75725508}, 00664 {1, -7.85722494}, 00665 {1, -7.91697407}, 00666 {1, -8.04758644}, 00667 {1, -8.09447479}, 00668 {1, -8.15859795}, 00669 {1, -8.25424385}, 00670 {1, -8.33880615}, 00671 {1, -8.41452408}, 00672 {1, -8.54581165}, 00673 {1, -8.60400581}, 00674 {1, -8.65751839}, 00675 {1, -8.75414848}, 00676 {1, -8.83946800}, 00677 {1, -8.91589737}, 00678 {1, -9.01741695}, 00679 {1, -9.10663033}, 00680 {1, -9.18621922}, 00681 {1, -9.28059292}, 00682 {1, -9.34430218}, 00683 {1, -9.42154408}, 00684 {1, -9.55562973}, 00685 {1, -9.61459446}, 00686 {1, -9.72023010}, 00687 {1, -9.76802444}, 00688 {1, -9.85540199}, 00689 {1, -9.93374062}, 00690 {1, -10.03800774}, 00691 {1, -10.10044670}, 00692 {1, -10.21178055}, 00693 {1, -10.29757786}, 00694 {1, -10.39561272}, 00695 {1, -10.45469666}, 00696 {1, -10.56102180}, 00697 {1, -10.66205502}, 00698 {1, -10.73780537}, 00699 {1, -10.82557774}, 00700 {1, -10.91007328}, 00701 {1, -10.98659325}, 00702 {1, -11.07857418}, 00703 {1, -11.19975281}, 00704 {1, -11.27170753}, 00705 {1, -11.36930943}, 00706 {1, -11.45675945}, 00707 {1, -11.53620148}, 00708 {1, -11.63198853}, 00709 {1, -11.73875237}, 00710 {1, -11.83400822}, 00711 {1, -11.93677044}, 00712 {1, -12.02933311}, 00713 {1, -12.13374519}, 00714 {1, -12.23410702}, 00715 {1, -12.33664989}, 00716 {1, -12.44163322}, 00717 {1, -12.54730415}, 00718 {1, -12.64975739}, 00719 {1, -12.76682186}, 00720 {1, -12.88185978}, 00721 {1, -13.00052643}, 00722 {1, -13.12289810}, 00723 {1, -13.26689529}, 00724 {1, -13.39390945}, 00725 {1, -30.00000000}, 00726 {0L, 0 } 00727 }; 00728 for(_i=_r=0L; _i < 83; _i++) 00729 a0[_i] = RC_INI(_rs3); } 00730 { static struct{ long rc; double ini; } _rs4[] = { 00731 {1, 4.53776000e05}, 00732 {1, 3.48304000e05}, 00733 {1, 2.80224000e05}, 00734 {1, 1.98128000e05}, 00735 {1, 1.51219797e05}, 00736 {1, 1.21113266e05}, 00737 {1, 8.52812109e04}, 00738 {1, 6.49598125e04}, 00739 {1, 5.20075781e04}, 00740 {1, 3.66190977e04}, 00741 {1, 2.79060723e04}, 00742 {1, 2.23634102e04}, 00743 {1, 1.57683135e04}, 00744 {1, 1.20284307e04}, 00745 {1, 9.17755273e03}, 00746 {1, 7.36044873e03}, 00747 {1, 5.19871680e03}, 00748 {1, 3.97240796e03}, 00749 {1, 3.18934326e03}, 00750 {1, 2.25737622e03}, 00751 {1, 1.72767114e03}, 00752 {1, 1.46202722e03}, 00753 {1, 1.32456628e03}, 00754 {1, 1.06499792e03}, 00755 {1, 9.92735291e02}, 00756 {1, 8.91604858e02}, 00757 {1, 7.59411560e02}, 00758 {1, 6.59120056e02}, 00759 {1, 5.80688965e02}, 00760 {1, 4.66602264e02}, 00761 {1, 4.27612854e02}, 00762 {1, 3.91531494e02}, 00763 {1, 3.34516968e02}, 00764 {1, 2.91021820e02}, 00765 {1, 2.56853912e02}, 00766 {1, 2.17598007e02}, 00767 {1, 1.88145462e02}, 00768 {1, 1.65329865e02}, 00769 {1, 1.41960342e02}, 00770 {1, 1.28181686e02}, 00771 {1, 1.13336761e02}, 00772 {1, 9.17785034e01}, 00773 {1, 8.36242981e01}, 00774 {1, 7.08843536e01}, 00775 {1, 6.58346100e01}, 00776 {1, 5.75790634e01}, 00777 {1, 5.11293755e01}, 00778 {1, 4.37563019e01}, 00779 {1, 3.99226875e01}, 00780 {1, 3.39562836e01}, 00781 {1, 3.00413170e01}, 00782 {1, 2.61871891e01}, 00783 {1, 2.41310368e01}, 00784 {1, 2.08853607e01}, 00785 {1, 1.82632275e01}, 00786 {1, 1.60007000e01}, 00787 {1, 1.42617064e01}, 00788 {1, 1.27951088e01}, 00789 {1, 1.16221066e01}, 00790 {1, 1.03779335e01}, 00791 {1, 8.97864914e00}, 00792 {1, 8.25593281e00}, 00793 {1, 7.39339924e00}, 00794 {1, 6.70784378e00}, 00795 {1, 6.16084862e00}, 00796 {1, 5.57818031e00}, 00797 {1, 5.01341105e00}, 00798 {1, 4.55679178e00}, 00799 {1, 4.13692093e00}, 00800 {1, 3.80004382e00}, 00801 {1, 3.46328306e00}, 00802 {1, 3.17340493e00}, 00803 {1, 2.93525696e00}, 00804 {1, 2.69083858e00}, 00805 {1, 2.46588683e00}, 00806 {1, 2.26083040e00}, 00807 {1, 2.04337358e00}, 00808 {1, 1.89027369e00}, 00809 {1, 1.69208312e00}, 00810 {0L, 0 } 00811 }; 00812 for(_i=_r=0L; _i < 79; _i++) 00813 ba0[_i] = (realnum)RC_INI(_rs4); } 00814 { static struct{ long rc; double ini; } _rs5[] = { 00815 {1, 1.48992336e00}, 00816 {1, 1.32466662e00}, 00817 {1, 1.10697949e00}, 00818 {1, 9.29813862e-01}, 00819 {0L, 0 } 00820 }; 00821 for(_i=_r=0L; _i < 4; _i++) 00822 bb0[_i] = (realnum)RC_INI(_rs5); } 00823 { static struct{ long rc; double ini; } _rs6[] = { 00824 {1, 2.12597656}, 00825 {1, 2.08984375}, 00826 {1, 2.06958008}, 00827 {1, 2.05444336}, 00828 {1, 2.05}, 00829 {1, 2.05}, 00830 {1, 2.05}, 00831 {1, 2.05}, 00832 {1, 2.05}, 00833 {1, 2.05}, 00834 {1, 2.05}, 00835 {1, 2.05}, 00836 {1, 2.05}, 00837 {1, 2.05}, 00838 {1, 2.05}, 00839 {1, 2.05}, 00840 {1, 2.05}, 00841 {1, 2.05}, 00842 {1, 2.05}, 00843 {1, 2.05}, 00844 {1, 2.05}, 00845 {1, 2.05}, 00846 {1, 2.05}, 00847 {1, 2.05}, 00848 {1, 2.05}, 00849 {1, 2.05}, 00850 {1, 2.05}, 00851 {1, 2.05}, 00852 {1, 2.05}, 00853 {1, 2.05}, 00854 {1, 2.05}, 00855 {1, 2.05}, 00856 {1, 2.05}, 00857 {1, 2.05}, 00858 {1, 2.05}, 00859 {1, 2.05}, 00860 {1, 2.05}, 00861 {1, 2.05}, 00862 {1, 2.05}, 00863 {1, 2.05}, 00864 {1, 2.05}, 00865 {1, 2.05}, 00866 {1, 2.05}, 00867 {1, 2.05}, 00868 {1, 2.05}, 00869 {1, 2.05}, 00870 {1, 2.05}, 00871 {1, 2.05}, 00872 {1, 2.05}, 00873 {1, 2.05}, 00874 {1, 2.05}, 00875 {1, 2.05}, 00876 {1, 2.05}, 00877 {1, 2.05}, 00878 {1, 2.05}, 00879 {1, 2.05}, 00880 {1, 2.05}, 00881 {1, 2.05}, 00882 {1, 2.05}, 00883 {1, 2.05}, 00884 {1, 2.05}, 00885 {1, 2.05}, 00886 {1, 2.05}, 00887 {1, 2.05}, 00888 {1, 2.05}, 00889 {1, 2.05}, 00890 {1, 2.05}, 00891 {1, 2.05}, 00892 {1, 2.05}, 00893 {1, 2.05}, 00894 {1, 2.05}, 00895 {1, 2.05}, 00896 {1, 2.05}, 00897 {1, 2.05}, 00898 {1, 2.05}, 00899 {1, 2.05}, 00900 {1, 2.05}, 00901 {1, 2.05}, 00902 {1, 2.05}, 00903 {1, 2.05}, 00904 {1, 2.05}, 00905 {1, 2.05}, 00906 {1, 2.05}, 00907 {0L, 0 } 00908 }; 00909 for(_i=_r=0L; _i < 83; _i++) 00910 x0[_i] = RC_INI(_rs6); } 00911 00912 { static struct{ long rc; double ini; } _rs7[] = { 00913 {1, 16.23337936}, 00914 {1, 16.27946854}, 00915 {1, 16.31696320}, 00916 {1, 16.37597656}, 00917 {1, 16.42210960}, 00918 {1, 16.45996284}, 00919 {1, 16.51994896}, 00920 {1, 16.56644440}, 00921 {1, 16.60460854}, 00922 {1, 16.66510773}, 00923 {1, 16.71198654}, 00924 {1, 16.75038719}, 00925 {1, 16.81106949}, 00926 {1, 16.85778809}, 00927 {1, 16.90416527}, 00928 {1, 16.94209099}, 00929 {1, 17.00195694}, 00930 {1, 17.04838943}, 00931 {1, 17.08633804}, 00932 {1, 17.14627838}, 00933 {1, 17.19270515}, 00934 {1, 17.22186279}, 00935 {1, 17.23933601}, 00936 {1, 17.27728271}, 00937 {1, 17.30161858}, 00938 {1, 17.32085800}, 00939 {1, 17.34928894}, 00940 {1, 17.37349129}, 00941 {1, 17.39528084}, 00942 {1, 17.43282318}, 00943 {1, 17.44827652}, 00944 {1, 17.46357536}, 00945 {1, 17.49082375}, 00946 {1, 17.51517677}, 00947 {1, 17.53697205}, 00948 {1, 17.56587219}, 00949 {1, 17.59125519}, 00950 {1, 17.61410332}, 00951 {1, 17.64081383}, 00952 {1, 17.65900803}, 00953 {1, 17.68086433}, 00954 {1, 17.71843529}, 00955 {1, 17.73512840}, 00956 {1, 17.76512146}, 00957 {1, 17.77873421}, 00958 {1, 17.80340767}, 00959 {1, 17.82530022}, 00960 {1, 17.85470963}, 00961 {1, 17.87210464}, 00962 {1, 17.90334511}, 00963 {1, 17.92751503}, 00964 {1, 17.95458603}, 00965 {1, 17.97117233}, 00966 {1, 18.00062943}, 00967 {1, 18.02842712}, 00968 {1, 18.04934502}, 00969 {1, 18.07340050}, 00970 {1, 18.09639168}, 00971 {1, 18.11732864}, 00972 {1, 18.14218903}, 00973 {1, 18.17465591}, 00974 {1, 18.19370079}, 00975 {1, 18.21962166}, 00976 {1, 18.24237251}, 00977 {1, 18.26305962}, 00978 {1, 18.28767967}, 00979 {1, 18.31531525}, 00980 {1, 18.33900452}, 00981 {1, 18.36478043}, 00982 {1, 18.38741112}, 00983 {1, 18.41271973}, 00984 {1, 18.43644333}, 00985 {1, 18.46075630}, 00986 {1, 18.48509216}, 00987 {1, 18.50897980}, 00988 {1, 18.53143501}, 00989 {1, 18.55570030}, 00990 {1, 18.58008003}, 00991 {1, 18.60348320}, 00992 {1, 18.62536430}, 00993 {1, 18.65199852}, 00994 {1, 18.67623520}, 00995 {1, 18.70072174}, 00996 {0L, 0 } 00997 }; 00998 for(_i=_r=0L; _i < 83; _i++) 00999 a1[_i] = RC_INI(_rs7); } 01000 { static struct{ long rc; double ini; } _rs8[] = { 01001 {1, 1.09462866e10}, 01002 {1, 9.32986675e09}, 01003 {1, 6.15947008e09}, 01004 {1, 1.54486170e09}, 01005 {1, 1.00812454e09}, 01006 {1, 7.00559552e08}, 01007 {1, 6.25999232e08}, 01008 {1, 3.50779968e08}, 01009 {1, 3.11956288e08}, 01010 {1, 3.74866016e08}, 01011 {1, 2.47019424e08}, 01012 {1, 1.73169776e08}, 01013 {1, 1.01753168e08}, 01014 {1, 6.81861920e07}, 01015 {1, 4.61764000e07}, 01016 {1, 3.31671360e07}, 01017 {1, 2.03160540e07}, 01018 {1, 1.40249480e07}, 01019 {1, 1.02577860e07}, 01020 {1, 3.53822650e06}, 01021 {1, 1.32563388e06}, 01022 {1, 9.14284688e05}, 01023 {1, 1.25230388e06}, 01024 {1, 3.17865156e05}, 01025 {1, 4.76750244e03}, 01026 {1, 4.81107031e03}, 01027 {1, 4.88406152e03}, 01028 {1, 4.80611279e03}, 01029 {1, 4.78843652e03}, 01030 {1, 4.65988477e03}, 01031 {1, 1.26723059e03}, 01032 {1, 1.20825342e03}, 01033 {1, 8.66052612e02}, 01034 {1, 7.76661316e02}, 01035 {1, 7.05279358e02}, 01036 {1, 6.21722656e02}, 01037 {1, 5.46207581e02}, 01038 {1, 4.96247742e02}, 01039 {1, 4.26340118e02}, 01040 {1, 3.96090424e02}, 01041 {1, 3.48429657e02}, 01042 {1, 2.37949142e02}, 01043 {1, 2.14678406e02}, 01044 {1, 1.81019180e02}, 01045 {1, 1.68923676e02}, 01046 {1, 1.45979385e02}, 01047 {1, 1.25311127e02}, 01048 {1, 1.05205528e02}, 01049 {1, 9.39378357e01}, 01050 {1, 7.75339966e01}, 01051 {1, 6.68987427e01}, 01052 {1, 5.53580055e01}, 01053 {1, 5.00100212e01}, 01054 {1, 4.14198608e01}, 01055 {1, 3.46289063e01}, 01056 {1, 3.00775223e01}, 01057 {1, 2.60294399e01}, 01058 {1, 2.26602840e01}, 01059 {1, 2.02123032e01}, 01060 {1, 1.76353855e01}, 01061 {1, 1.47198439e01}, 01062 {1, 1.33078461e01}, 01063 {1, 1.17181997e01}, 01064 {1, 1.04125805e01}, 01065 {1, 9.45785904e00}, 01066 {1, 8.42799950e00}, 01067 {1, 7.62769842e00}, 01068 {1, 6.85484743e00}, 01069 {1, 6.25903368e00}, 01070 {1, 5.75135279e00}, 01071 {1, 5.28468180e00}, 01072 {1, 4.87669659e00}, 01073 {1, 4.57353973e00}, 01074 {1, 4.30108690e00}, 01075 {1, 4.05412245e00}, 01076 {1, 3.83283114e00}, 01077 {1, 3.57902861e00}, 01078 {1, 3.43705726e00}, 01079 {1, 3.26563096e00}, 01080 {0L, 0 } 01081 }; 01082 for(_i=_r=0L; _i < 79; _i++) 01083 ba1[_i] = (realnum)RC_INI(_rs8); } 01084 { static struct{ long rc; double ini; } _rs9[] = { 01085 {1, 3.07498097e00}, 01086 {1, 2.96334076e00}, 01087 {1, 2.92890000e00}, 01088 {1, 2.89550042e00}, 01089 {0L, 0 } 01090 }; 01091 for(_i=_r=0L; _i < 4; _i++) 01092 bb1[_i] = (realnum)RC_INI(_rs9); } 01093 { static struct{ long rc; double ini; } _rs10[] = { 01094 {1, -5.46}, 01095 {1, -5.51}, 01096 {1, -5.49}, 01097 {1, -5.30}, 01098 {1, -5.29}, 01099 {1, -5.28}, 01100 {1, -5.37}, 01101 {1, -5.33}, 01102 {1, -5.38}, 01103 {1, -5.55}, 01104 {1, -5.55}, 01105 {1, -5.55}, 01106 {1, -5.55}, 01107 {1, -5.55}, 01108 {1, -5.55}, 01109 {1, -5.55}, 01110 {1, -5.55}, 01111 {1, -5.55}, 01112 {1, -5.55}, 01113 {1, -5.38}, 01114 {1, -5.19}, 01115 {1, -5.14}, 01116 {1, -5.27}, 01117 {1, -4.93}, 01118 {1, -3.64}, 01119 {1, -3.68}, 01120 {1, -3.74}, 01121 {1, -3.78}, 01122 {1, -3.82}, 01123 {1, -3.88}, 01124 {1, -3.40}, 01125 {1, -3.41}, 01126 {1, -3.32}, 01127 {1, -3.32}, 01128 {1, -3.32}, 01129 {1, -3.32}, 01130 {1, -3.31}, 01131 {1, -3.31}, 01132 {1, -3.29}, 01133 {1, -3.29}, 01134 {1, -3.27}, 01135 {1, -3.16}, 01136 {1, -3.14}, 01137 {1, -3.11}, 01138 {1, -3.10}, 01139 {1, -3.07}, 01140 {1, -3.03}, 01141 {1, -2.99}, 01142 {1, -2.96}, 01143 {1, -2.91}, 01144 {1, -2.87}, 01145 {1, -2.81}, 01146 {1, -2.78}, 01147 {1, -2.72}, 01148 {1, -2.66}, 01149 {1, -2.61}, 01150 {1, -2.56}, 01151 {1, -2.51}, 01152 {1, -2.47}, 01153 {1, -2.42}, 01154 {1, -2.35}, 01155 {1, -2.31}, 01156 {1, -2.26}, 01157 {1, -2.21}, 01158 {1, -2.17}, 01159 {1, -2.12}, 01160 {1, -2.08}, 01161 {1, -2.03}, 01162 {1, -1.99}, 01163 {1, -1.95}, 01164 {1, -1.91}, 01165 {1, -1.87}, 01166 {1, -1.84}, 01167 {1, -1.81}, 01168 {1, -1.78}, 01169 {1, -1.75}, 01170 {1, -1.71}, 01171 {1, -1.69}, 01172 {1, -1.66}, 01173 {1, -1.62}, 01174 {1, -1.60}, 01175 {1, -1.60}, 01176 {1, -1.60}, 01177 {0L, 0 } 01178 }; 01179 for(_i=_r=0L; _i < 83; _i++) 01180 x1[_i] = RC_INI(_rs10); } 01181 { static struct{ long rc; double ini; } _rs11[] = { 01182 {1, 20.30049515}, 01183 {1, 20.28500366}, 01184 {1, 20.25300407}, 01185 {1, 20.16626740}, 01186 {1, 20.15743256}, 01187 {1, 20.11256981}, 01188 {1, 20.04818344}, 01189 {1, 19.99261856}, 01190 {1, 19.94472885}, 01191 {1, 19.86478043}, 01192 {1, 19.83321571}, 01193 {1, 19.80185127}, 01194 {1, 19.74884224}, 01195 {1, 19.70136070}, 01196 {1, 19.65981102}, 01197 {1, 19.60598755}, 01198 {1, 19.56017494}, 01199 {1, 19.52042389}, 01200 {1, 19.47429657}, 01201 {1, 19.44413757}, 01202 {1, 19.40796280}, 01203 {1, 19.34819984}, 01204 {1, 19.32203293}, 01205 {1, 19.27634430}, 01206 {1, 19.25627136}, 01207 {1, 19.22009087}, 01208 {1, 19.18853378}, 01209 {1, 19.14809799}, 01210 {1, 19.12456703}, 01211 {1, 19.08409119}, 01212 {1, 19.05431557}, 01213 {1, 19.02083015}, 01214 {1, 19.00176430}, 01215 {1, 18.96817970}, 01216 {1, 18.93762589}, 01217 {1, 18.91706085}, 01218 {1, 18.89299583}, 01219 {1, 18.87085915}, 01220 {1, 18.85210609}, 01221 {1, 18.83035851}, 01222 {1, 18.80403900}, 01223 {1, 18.78901100}, 01224 {1, 18.77099228}, 01225 {1, 18.75540161}, 01226 {1, 18.74287033}, 01227 {1, 18.72928810}, 01228 {1, 18.71601868}, 01229 {1, 18.70474434}, 01230 {1, 18.69515800}, 01231 {1, 18.68782425}, 01232 {1, 18.68120766}, 01233 {1, 18.67630005}, 01234 {1, 18.67357445}, 01235 {1, 18.67129898}, 01236 {1, 18.67042351}, 01237 {1, 18.67090988}, 01238 {1, 18.67313004}, 01239 {1, 18.67636490}, 01240 {1, 18.68120003}, 01241 {1, 18.68803024}, 01242 {1, 18.69487381}, 01243 {1, 18.70458412}, 01244 {1, 18.71205139}, 01245 {0L, 0 } 01246 }; 01247 for(_i=_r=0L; _i < 63; _i++) 01248 a2[_i] = RC_INI(_rs11); } 01249 { static struct{ long rc; double ini; } _rs12[] = { 01250 {1, 1.01078403e00}, 01251 {1, 1.97956896e00}, 01252 {1, 3.14605665e00}, 01253 {1, 6.46874905e00}, 01254 {1, 3.16406364e01}, 01255 {1, 3.74927673e01}, 01256 {1, 4.75353088e01}, 01257 {1, 5.27809143e01}, 01258 {1, 5.86515846e01}, 01259 {1, 6.70408707e01}, 01260 {1, 1.14904137e02}, 01261 {1, 1.03133133e02}, 01262 {1, 1.26508728e02}, 01263 {1, 1.03827606e02}, 01264 {1, 8.79508896e01}, 01265 {1, 7.18328934e01}, 01266 {1, 6.19807892e01}, 01267 {1, 5.51255455e01}, 01268 {1, 4.87156143e01}, 01269 {1, 4.58579826e01}, 01270 {1, 4.19952011e01}, 01271 {1, 4.08252220e01}, 01272 {1, 3.78243637e01}, 01273 {1, 3.34573860e01}, 01274 {1, 3.19036102e01}, 01275 {1, 2.92026005e01}, 01276 {1, 2.74482193e01}, 01277 {1, 2.54643936e01}, 01278 {1, 2.46636391e01}, 01279 {1, 2.33054180e01}, 01280 {1, 2.23069897e01}, 01281 {1, 2.12891216e01}, 01282 {1, 2.06667900e01}, 01283 {1, 1.96430798e01}, 01284 {1, 1.87381802e01}, 01285 {1, 1.76523514e01}, 01286 {1, 1.69235287e01}, 01287 {1, 1.62647285e01}, 01288 {1, 1.56806908e01}, 01289 {1, 1.50346069e01}, 01290 {1, 1.42240467e01}, 01291 {1, 1.37954988e01}, 01292 {1, 1.31949224e01}, 01293 {1, 1.27211905e01}, 01294 {1, 1.22885675e01}, 01295 {1, 1.17868662e01}, 01296 {1, 1.12577572e01}, 01297 {1, 1.08565578e01}, 01298 {1, 1.04121590e01}, 01299 {1, 1.00410652e01}, 01300 {1, 9.64534473e00}, 01301 {1, 9.29232121e00}, 01302 {1, 8.92519569e00}, 01303 {1, 8.60898972e00}, 01304 {1, 8.31234550e00}, 01305 {1, 8.04089737e00}, 01306 {1, 7.74343491e00}, 01307 {1, 7.48133039e00}, 01308 {1, 7.21957016e00}, 01309 {1, 6.94726801e00}, 01310 {1, 6.71931219e00}, 01311 {1, 6.45107985e00}, 01312 {1, 6.28593779e00}, 01313 {0L, 0 } 01314 }; 01315 for(_i=_r=0L; _i < 63; _i++) 01316 b2[_i] = RC_INI(_rs12); } 01317 { static struct{ long rc; double ini; } _rs13[] = { 01318 {1, -0.43}, 01319 {1, -0.75}, 01320 {1, -0.93}, 01321 {1, -1.20}, 01322 {1, -1.78}, 01323 {1, -1.85}, 01324 {1, -1.95}, 01325 {1, -2.00}, 01326 {1, -2.05}, 01327 {1, -2.12}, 01328 {1, -2.34}, 01329 {1, -2.31}, 01330 {1, -2.42}, 01331 {1, -2.36}, 01332 {1, -2.31}, 01333 {1, -2.25}, 01334 {1, -2.21}, 01335 {1, -2.18}, 01336 {1, -2.15}, 01337 {1, -2.14}, 01338 {1, -2.12}, 01339 {1, -2.14}, 01340 {1, -2.12}, 01341 {1, -2.09}, 01342 {1, -2.08}, 01343 {1, -2.06}, 01344 {1, -2.05}, 01345 {1, -2.04}, 01346 {1, -2.04}, 01347 {1, -2.04}, 01348 {1, -2.04}, 01349 {1, -2.04}, 01350 {1, -2.04}, 01351 {1, -2.04}, 01352 {1, -2.04}, 01353 {1, -2.04}, 01354 {1, -2.04}, 01355 {1, -2.04}, 01356 {1, -2.04}, 01357 {1, -2.04}, 01358 {1, -2.04}, 01359 {1, -2.04}, 01360 {1, -2.04}, 01361 {1, -2.04}, 01362 {1, -2.04}, 01363 {1, -2.04}, 01364 {1, -2.04}, 01365 {1, -2.04}, 01366 {1, -2.04}, 01367 {1, -2.04}, 01368 {1, -2.04}, 01369 {1, -2.04}, 01370 {1, -2.04}, 01371 {1, -2.04}, 01372 {1, -2.04}, 01373 {1, -2.04}, 01374 {1, -2.04}, 01375 {1, -2.04}, 01376 {1, -2.04}, 01377 {1, -2.04}, 01378 {1, -2.04}, 01379 {1, -2.04}, 01380 {1, -2.04}, 01381 {0L, 0 } 01382 }; 01383 for(_i=_r=0L; _i < 63; _i++) 01384 x2[_i] = RC_INI(_rs13); } 01385 /*lint +e736 loss of precision in assignment in translated data */ 01386 01387 DEBUG_ENTRY( "blkdata0()" ); 01388 } 01389 01390 /* =================================================================== */ 01391 /*xmap mapping function for Cota's 3-body recombination */ 01392 STATIC double xmap(double x[], 01393 double y[], 01394 double xmapx0) 01395 { 01396 double a, 01397 b, 01398 c, 01399 xmapx1, 01400 x12m, 01401 x13m, 01402 xmapx2, 01403 x3, 01404 xmap_v, 01405 yinit, 01406 y13m; 01407 01408 DEBUG_ENTRY( "xmap()" ); 01409 01410 /* PARABOLIC INTERPOLATION. 01411 * */ 01412 01413 yinit = y[0]; 01414 xmapx1 = x[0]; 01415 xmapx2 = x[1]; 01416 x3 = x[2]; 01417 x13m = xmapx1 - x3; 01418 x12m = xmapx1 - xmapx2; 01419 y13m = yinit - y[2]; 01420 x3 = (xmapx1 + x3)*x13m; 01421 xmapx2 = (xmapx1 + xmapx2)*x12m; 01422 b = ((yinit - y[1])*x3 - y13m*xmapx2)/(x12m*x3 - x13m*xmapx2); 01423 a = (y13m - x13m*b)/x3; 01424 c = yinit - a*xmapx1*xmapx1 - b*xmapx1; 01425 01426 xmap_v = a*xmapx0*xmapx0 + b*xmapx0 + c; 01427 01428 return( xmap_v ); 01429 } 01430 01431 /* =================================================================== */ 01432 /*xinvrs do inverse function for Cota's three-body recombination */ 01433 STATIC double xinvrs(double y, 01434 double a, 01435 double b, 01436 double u, 01437 double v, 01438 long int *ifail) 01439 { 01440 long int i; 01441 double bxu, 01442 dfx, 01443 fx, 01444 fxdfx, 01445 x, 01446 xinvrs_v, 01447 xlog, 01448 xx; 01449 static long itmax = 10; 01450 01451 DEBUG_ENTRY( "xinvrs()" ); 01452 01453 /* inverts equation of the form : 01454 * 01455 * Y = A + B * X ** U - V * LOG ( X ) 01456 * */ 01457 *ifail = 0; 01458 xlog = (a - y)/v; 01459 x = pow(10.,xlog); 01460 xx = 0.; 01461 01462 for( i=0; i < itmax; i++ ) 01463 { 01464 bxu = b*pow(x,u); 01465 fx = y - a - bxu + v*xlog; 01466 dfx = v*.4342945 - bxu*u; 01467 01468 if( dfx != 0. ) 01469 { 01470 fxdfx = fabs(fx/dfx); 01471 fxdfx = MIN2(0.2,fxdfx); 01472 xx = x*(1. - sign(fxdfx,fx/dfx)); 01473 } 01474 else 01475 { 01476 /* >>chng 96 feb 02 this added in case dfx ever 0 01477 * suggested by Peter van Hoof */ 01478 xx = x*(1. - sign(0.2,fx)); 01479 } 01480 01481 if( (fabs(xx-x)/x) < 1.e-4 ) 01482 { 01483 xinvrs_v = xx; 01484 return( xinvrs_v ); 01485 } 01486 else 01487 { 01488 x = xx; 01489 if( x < 1e-30 ) 01490 { 01491 xinvrs_v = 100.; 01492 *ifail = 1; 01493 return( xinvrs_v ); 01494 } 01495 xlog = log10(x); 01496 } 01497 } 01498 xinvrs_v = xx; 01499 *ifail = 1; 01500 return( xinvrs_v ); 01501 }