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);
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));
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)
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 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
00188
00189
00190
00191
00192
00193
00194
00195
00196
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
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
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
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
00300
00301
00302
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
00473
00474
00475
00476
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
00487
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
01386
01387 DEBUG_ENTRY( "blkdata0()" );
01388 }
01389
01390
01391
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
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
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
01454
01455
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
01477
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 }