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