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