root/branches/newmole/source/atmdat_3body.cpp

Revision 1739, 27.8 kB (checked in by rjrw, 12 months ago)

Merged from trunk r1700:1738

  • Property svn:eol-style set to native
Line 
1/* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
2 * others.  For conditions of distribution and use see copyright notice in license.txt */
3/*atmdat_3body derive three-body recombination coefficients */
4/*da interpolate on three body recombination by Steve Cota */
5#include "cddefines.h"
6#include "ionbal.h"
7#include "phycon.h"
8#include "dense.h"
9#include "trace.h"
10#include "punch.h"
11#include "atmdat.h"
12
13#define MAXZ    28
14
15STATIC void blkdata1(void);
16STATIC double da(double z);
17
18static double a2[63],
19          b2[63],
20          x2[63];
21
22static double a0[83],
23          x0[83];
24static realnum b0[83],
25          b1[83];
26
27static double a1[83],
28          x1[83];
29
30static double tz[83],
31          zlog7[28],
32          zlog2[28];
33
34#define RC_INI(rs)      (rs[_r].rc>1 ? DEC_RC_(rs) : (rs[_r].rc==1 ? INC_NDX_(rs) : rs[_r].ini ))
35#define DEC_RC_(rs)     (rs[_r].rc--,rs[_r].ini)
36#define INC_NDX_(rs)    (_r++,rs[_r-1].ini)
37
38/* this "mapping function" occurs below */
39STATIC double xmap(double x[],
40          double y[],
41          double x0);
42
43/* inverse routine, also below */
44STATIC double xinvrs(double y,
45          double a,
46          double b,
47          double u,
48          double v,
49          long int *ifail);
50
51/* =================================================================== */
52void atmdat_3body(void)
53{
54        long int i,
55          iup;
56
57        DEBUG_ENTRY( "atmdat_3body()" );
58
59
60        if( ionbal.lgNoCota )
61        {
62                for( i=0; i < LIMELM; i++ )
63                {
64                        ionbal.CotaRate[i] = 0.;
65                }
66                atmdat.nsbig = 0;
67                return;
68        }
69
70        if( atmdat.nsbig == 0 )
71        {
72                /* steve cota only defined things up to 28 */
73                iup = MIN2(28,LIMELM);
74        }
75        else
76        {
77                iup = MIN3( LIMELM , atmdat.nsbig , 28 );
78        }
79
80        for( i=0; i < iup; i++ )
81        {
82                ionbal.CotaRate[i] = (realnum)da((double)(i+1));
83        }
84
85        atmdat.nsbig = 0;
86
87        if( trace.lgTrace && trace.lgTrace3Bod )
88        {
89                fprintf( ioQQQ, "     3BOD rate:" );
90                for( i=1; i <= 14; i++ )
91                {
92                        fprintf( ioQQQ, "%8.1e", ionbal.CotaRate[i-1] );
93                }
94                fprintf( ioQQQ, "\n" );
95        }
96
97        if( punch.lgioRecom )
98        {
99                /* option to punch coefficients */
100                fprintf( punch.ioRecom, " 3-body rec coef vs charge \n" );
101                for( i=0; i < iup; i++ )
102                {
103                        fprintf( punch.ioRecom, "%3ld%10.2e\n", i+1, ionbal.CotaRate[i] );
104                }
105                fprintf( punch.ioRecom, "\n");
106        }
107        return;
108}
109
110/* =================================================================== */
111STATIC double da(double z)
112{
113        /*lint -e736 loss of precision in assignment in translated data */
114        long int jfail,
115          nt,
116          nt0,
117          nt1,
118          nz;
119
120        static bool lgCalled=false;
121        double a,
122          alogn,
123          alognc,
124          alogt,
125          b,
126          c,
127          d,
128          da_v,
129          expp,
130          u,
131          v,
132          x,
133          xnc,
134          y,
135          zlog;
136        double yya[3],
137          xx[3],
138          yyb[3],
139          yyx[3],
140          yyy[3];
141
142        /* alogte is base 10 log of temperature and elec density*/
143        double alogte , alogne;
144
145        DEBUG_ENTRY( "da()" );
146
147        /* WRITTEN BY S. A. COTA, 2/87
148         * */
149
150        /* MAXZ IS THE MAXIMUM EFFECTIVE NUCLEAR CHARGE ( = IONIC CHARGE + 1 )
151         * WHICH THE DIMENSION STATEMENTS ACCOMODATE. 
152         *
153         * IT IS USED ONLY FOR THE ARRAY ZLOG7 ( = 7 * LOG ( Z ) )
154         * AND THE ARRAY ZLOG2 ( = 2 * LOG ( Z ) ) .   THESE ARRAYS
155         * CONTAIN EASILY CALCULATED VALUES, WHICH HAVE BEEN STORED
156         * TO SAVE TIME IN EXECUTION.
157         *
158         * IF MAXZ IS EXCEEDED, THIS PROGRAM SIMPLY CALCULATES THE
159         * LOGS INSTEAD OF LOOKING THEM UP.
160         *
161         * */
162
163        if( !lgCalled )
164        {
165                lgCalled = true;
166                blkdata1();
167        }
168
169        /*begin sanity check */
170        ASSERT( zlog7[1] > 0. );
171
172        nz = (long)(z + .1);
173        alogne = log10(dense.eden);
174        alogte = log10(phycon.te);
175        if( nz > MAXZ )
176        {
177                zlog = log10(z);
178                alogt = alogte - 2.*zlog;
179                alogn = alogne - 7.*zlog;
180        }
181        else
182        {
183                alogt = alogte - zlog2[nz-1];
184                alogn = alogne - zlog7[nz-1];
185        }
186
187        /* CHECK IF PARAMETERS ARE WITHIN BOUNDS.  IF NOT, INCREMENT
188         * APPROPRIATE ERROR COUNTER AND SET TO BOUNDARY IF
189         * NECESSARY:
190         *
191         * DEFINITION OF ERROR COUNTERS:
192         *
193         * ILT    : LOW T         
194         * ILTLN  : LOW T , LOW  N
195         * ILTHN  : LOW T , HIGH N
196         * IHTHN  : HIGH T , HIGH N
197         * */
198        if( alogt < 0. )
199        {
200                ionbal.ilt += 1;
201                alogt = 0.;
202        }
203
204        if( alogt <= 2.1760913 )
205        {
206                if( alogn < (3.5*alogt - 8.) )
207                {
208                        ionbal.iltln += 1;
209                }
210                else if( alogn > (3.5*alogt - 2.) )
211                {
212                        ionbal.ilthn += 1;
213                        alogn = 3.5*alogt - 2.;
214                }
215
216        }
217        else if( alogt <= 2.4771213 )
218        {
219                if( alogn > 9. )
220                {
221                        ionbal.ilthn += 1;
222                        alogn = 9.;
223                }
224
225        }
226        else if( alogt <= 5.1139434 )
227        {
228                if( alogn > 13. )
229                {
230                        ionbal.ihthn += 1;
231                        alogn = 13.;
232                }
233
234        }
235        else
236        {
237                da_v = 0.;
238                return( da_v );
239        }
240
241        /* LOCATE POSITION IN ARRAYS */
242        if( alogt <= 2. )
243        {
244                nt = (long)(9.9657843*alogt + 1.);
245        }
246        else
247        {
248                nt = (long)(19.931568*alogt - 19.);
249        }
250        nt = MIN2(83,nt);
251        nt = MAX2(1,nt);
252
253        /* CENTER UP SINCE ARRAY VALUES ARE ROUNDED */
254        if( fabs(alogt-tz[nt-1]) >= fabs(alogt-tz[MIN2(83,nt+1)-1]) )
255        {
256                nt = MIN2(83,nt+1);
257        }
258        else if( fabs(alogt-tz[nt-1]) > fabs(alogt-tz[MAX2(1,nt-1)-1]) )
259        {
260                nt = MAX2(1,nt-1);
261        }
262
263        /* CHECK IF INTERPOLATION IS NEEDED AND PROCEED IF NOT.*/
264        if( fabs(alogt-tz[nt-1]) < 0.00001 )
265        {
266                if( z != 1.0 )
267                {
268                        c = a1[nt-1];
269                        d = b1[nt-1];
270                        u = x1[nt-1];
271                        v = 8.90;
272                }
273                else
274                {
275                        nt = MAX2(21,nt);
276                        nt = MIN2(83,nt);
277                        c = a2[nt-(21)];
278                        d = b2[nt-(21)];
279                        u = x2[nt-(21)];
280                        v = 9.40;
281                }
282
283                xnc = xinvrs(alogn,c,d,u,v,&jfail);
284                if( xnc <= 0. || jfail != 0 )
285                {
286                        ionbal.ifail = 1;
287                        jfail = 1;
288                        da_v = 0.;
289                        return( da_v );
290                }
291                alognc = log10(xnc);
292
293                a = a0[nt-1];
294                b = b0[nt-1];
295                x = -2.45;
296                y = x0[nt-1];
297                nt0 = nt - 1;
298
299                /* IF INTERPOLATION WAS REQUIRED,
300                 * CHECK THAT NT IS NOT ON THE EDGE OF A DISCONTINUITY,
301                 * EITHER AT END OF ARRAYS OR WITHIN THEM,
302                 * WHERE VALUES CHANGE ABRUPTLY.
303                 * */
304        }
305        else
306        {
307                if( (nt <= 21) && (z == 1.0) )
308                {
309                        nt = 22;
310                }
311                else if( nt <= 1 )
312                {
313                        nt = 2;
314                }
315                else if( nt >= 83 )
316                {
317                        nt = 82;
318                }
319                else if( nt == 24 )
320                {
321                        if( alogt <= 2.1760913 )
322                        {
323                                nt = 23;
324                        }
325                        else
326                        {
327                                nt = 25;
328                        }
329                }
330                else if( nt == 30 )
331                {
332                        if( alogt <= 2.4771213 )
333                        {
334                                nt = 29;
335                        }
336                        else
337                        {
338                                nt = 31;
339                        }
340                }
341
342                nt0 = nt - 1;
343                nt1 = nt + 1;
344                xx[0] = tz[nt0-1];
345                xx[1] = tz[nt-1];
346                xx[2] = tz[nt1-1];
347
348                if( z != 1.0 )
349                {
350                        if( nt0 == 24 )
351                        {
352                                yya[0] = 17.2880135;
353                                yyb[0] = 6.93410742e03;
354                                yyx[0] = -3.75;
355                        }
356                        else if( nt0 == 30 )
357                        {
358                                yya[0] = 17.4317989;
359                                yyb[0] = 1.36653821e03;
360                                yyx[0] = -3.40;
361                        }
362                        else
363                        {
364                                yya[0] = a1[nt0-1];
365                                yyb[0] = b1[nt0-1];
366                                yyx[0] = x1[nt0-1];
367                        }
368
369                        yya[1] = a1[nt-1];
370                        yya[2] = a1[nt1-1];
371                        c = xmap(xx,yya,alogt);
372                        yyb[1] = b1[nt-1];
373                        yyb[2] = b1[nt1-1];
374                        d = xmap(xx,yyb,alogt);
375                        yyx[1] = x1[nt-1];
376                        yyx[2] = x1[nt1-1];
377                        u = xmap(xx,yyx,alogt);
378                        v = 8.90;
379
380                }
381                else
382                {
383                        if( nt0 == 24 )
384                        {
385                                yya[0] = 20.1895161;
386                                yyb[0] = 2.25774918e01;
387                                yyx[0] = -1.66;
388                        }
389                        else if( nt0 == 30 )
390                        {
391                                yya[0] = 19.8647804;
392                                yyb[0] = 6.70408707e02;
393                                yyx[0] = -2.12;
394                        }
395                        else
396                        {
397                                yya[0] = a2[nt0-(21)];
398                                yyb[0] = b2[nt0-(21)];
399                                yyx[0] = x2[nt0-(21)];
400                        }
401
402                        yya[1] = a2[nt-(21)];
403                        yya[2] = a2[nt1-(21)];
404                        c = xmap(xx,yya,alogt);
405                        yyb[1] = b2[nt-(21)];
406                        yyb[2] = b2[nt1-(21)];
407                        d = xmap(xx,yyb,alogt);
408                        yyx[1] = x2[nt-(21)];
409                        yyx[2] = x2[nt1-(21)];
410                        u = xmap(xx,yyx,alogt);
411                        v = 9.40;
412                }
413
414                xnc = xinvrs(alogn,c,d,u,v,&jfail);
415                if( xnc <= 0. || jfail != 0 )
416                {
417                        ionbal.ifail = 1;
418                        jfail = 1;
419                        da_v = 0.;
420                        return( da_v );
421                }
422                alognc = log10(xnc);
423
424                if( nt0 == 24 )
425                {
426                        yya[0] = -8.04963875;
427                        yyb[0] = 1.07205127e03;
428                        yyy[0] = 2.05;
429                }
430                else if( nt0 == 30 )
431                {
432                        yya[0] = -8.54721069;
433                        yyb[0] = 4.70450195e02;
434                        yyy[0] = 2.05;
435                }
436                else
437                {
438                        yya[0] = a0[nt0-1];
439                        yyb[0] = b0[nt0-1];
440                        yyy[0] = x0[nt0-1];
441                }
442
443                yya[1] = a0[nt-1];
444                yya[2] = a0[nt1-1];
445                a = xmap(xx,yya,alogt);
446                yyb[1] = b0[nt-1];
447                yyb[2] = b0[nt1-1];
448                b = xmap(xx,yyb,alogt);
449                x = -2.45;
450                yyy[1] = x0[nt-1];
451                yyy[2] = x0[nt1-1];
452                y = xmap(xx,yyy,alogt);
453        }
454
455        expp = a - y*alognc + b*pow(xnc,x);
456        if( expp < 37 )
457        {
458                da_v = z*pow(10.,expp);
459        }
460        else
461        {
462                da_v = 0.;
463        }
464        ionbal.ifail += jfail;
465
466        return( da_v );
467}
468
469/****************************************************************************** */
470STATIC void blkdata1(void)
471{
472        /*block data with Steve Cota's 3-body recombination coefficients */
473
474        /* data for function  da.
475         *
476         * S. A. COTA, 2/1987
477         * */
478
479        long int _i,
480          _r;
481        realnum *const ba0 = (realnum*)b0;
482        realnum *const ba1 = (realnum*)b1;
483        realnum *const bb0 = (realnum*)((char*)(b0 + 79));
484        realnum *const bb1 = (realnum*)((char*)(b1 + 79));
485
486                /* to fix all the conversion errors, change realnum ini to double ini,
487                 * but chech that results still ok */
488                { static struct{ long rc; double ini; } _rs0[] = {
489                        {1,     0.00000e00},
490                        {1,     2.10721e00},
491                        {1,     3.33985e00},
492                        {1,     4.21442e00},
493                        {1,     4.89279e00},
494                        {1,     5.44706e00},
495                        {1,     5.91569e00},
496                        {1,     6.32163e00},
497                        {1,     6.67970e00},
498                        {1,     7.00000e00},
499                        {1,     7.28975e00},
500                        {1,     7.55427e00},
501                        {1,     7.79760e00},
502                        {1,     8.02290e00},
503                        {1,     8.23264e00},
504                        {1,     8.42884e00},
505                        {1,     8.61314e00},
506                        {1,     8.78691e00},
507                        {1,     8.95128e00},
508                        {1,     9.10721e00},
509                        {1,     9.25554e00},
510                        {1,     9.39696e00},
511                        {1,     9.53209e00},
512                        {1,     9.66148e00},
513                        {1,     9.78558e00},
514                        {1,     9.90481e00},
515                        {1,     10.01954e00},
516                        {1,     10.13010e00},
517                        {0L, 0}
518                };
519                for(_i=_r=0L; _i < 28; _i++)
520                        zlog7[_i] = RC_INI(_rs0); }
521                { static struct{ long rc; double ini; } _rs1[] = {
522                        {1,     0.00000e00},
523                        {1,     6.02060e-01},
524                        {1,     9.54243e-01},
525                        {1,     1.20412e00},
526                        {1,     1.39794e00},
527                        {1,     1.55630e00},
528                        {1,     1.69020e00},
529                        {1,     1.80618e00},
530                        {1,     1.90849e00},
531                        {1,     2.00000e00},
532                        {1,     2.08279e00},
533                        {1,     2.15836e00},
534                        {1,     2.22789e00},
535                        {1,     2.29226e00},
536                        {1,     2.35218e00},
537                        {1,     2.40824e00},
538                        {1,     2.46090e00},
539                        {1,     2.51055e00},
540                        {1,     2.55751e00},
541                        {1,     2.60206e00},
542                        {1,     2.64444e00},
543                        {1,     2.68485e00},
544                        {1,     2.72346e00},
545                        {1,     2.76042e00},
546                        {1,     2.79588e00},
547                        {1,     2.82995e00},
548                        {1,     2.86272e00},
549                        {1,     2.89431e00},
550                        {0L, 0}
551                };
552                for(_i=_r=0L; _i < 28; _i++)
553                        zlog2[_i] = RC_INI(_rs1); }
554                { static struct{ long rc; double ini; } _rs2[] = {
555                        {1,     0.},
556                        {1,     0.09691},
557                        {1,     0.17609},
558                        {1,     0.30103},
559                        {1,     0.39794},
560                        {1,     0.47712},
561                        {1,     0.60206},
562                        {1,     0.69897},
563                        {1,     0.77815},
564                        {1,     0.90309},
565                        {1,     1.00000},
566                        {1,     1.07918},
567                        {1,     1.20412},
568                        {1,     1.30103},
569                        {1,     1.39794},
570                        {1,     1.47712},
571                        {1,     1.60206},
572                        {1,     1.69897},
573                        {1,     1.77815},
574                        {1,     1.90309},
575                        {1,     2.00000},
576                        {1,     2.06070},
577                        {1,     2.09691},
578                        {1,     2.17609},
579                        {1,     2.20412},
580                        {1,     2.24304},
581                        {1,     2.30103},
582                        {1,     2.35218},
583                        {1,     2.39794},
584                        {1,     2.47712},
585                        {1,     2.51188},
586                        {1,     2.54407},
587                        {1,     2.60206},
588                        {1,     2.65321},
589                        {1,     2.69897},
590                        {1,     2.75967},
591                        {1,     2.81291},
592                        {1,     2.86034},
593                        {1,     2.91645},
594                        {1,     2.95424},
595                        {1,     3.00000},
596                        {1,     3.07918},
597                        {1,     3.11394},
598                        {1,     3.17609},
599                        {1,     3.20412},
600                        {1,     3.25527},
601                        {1,     3.30103},
602                        {1,     3.36173},
603                        {1,     3.39794},
604                        {1,     3.46240},
605                        {1,     3.51188},
606                        {1,     3.56820},
607                        {1,     3.60206},
608                        {1,     3.66276},
609                        {1,     3.72016},
610                        {1,     3.76343},
611                        {1,     3.81291},
612                        {1,     3.86034},
613                        {1,     3.90309},
614                        {1,     3.95424},
615                        {1,     4.02119},
616                        {1,     4.06070},
617                        {1,     4.11394},
618                        {1,     4.16137},
619                        {1,     4.20412},
620                        {1,     4.25527},
621                        {1,     4.31175},
622                        {1,     4.36173},
623                        {1,     4.41497},
624                        {1,     4.46240},
625                        {1,     4.51521},
626                        {1,     4.56526},
627                        {1,     4.61542},
628                        {1,     4.66605},
629                        {1,     4.71600},
630                        {1,     4.76343},
631                        {1,     4.81624},
632                        {1,     4.86629},
633                        {1,     4.91645},
634                        {1,     4.96614},
635                        {1,     5.02119},
636                        {1,     5.06726},
637                        {1,     5.11394},
638                        {0L, 0 }
639                        };
640                for(_i=_r=0L; _i < 83; _i++)
641                        tz[_i] = RC_INI(_rs2); }
642                { static struct{ long rc; double ini; } _rs3[] = {
643                        {1,     -4.31396484},
644                        {1,     -4.56640625},
645                        {1,     -4.74560547},
646                        {1,     -4.98535156},
647                        {1,     -5.15373850},
648                        {1,     -5.28123093},
649                        {1,     -5.48215008},
650                        {1,     -5.63811255},
651                        {1,     -5.76573515},
652                        {1,     -5.96755028},
653                        {1,     -6.12449837},
654                        {1,     -6.25304174},
655                        {1,     -6.45615673},
656                        {1,     -6.61384058},
657                        {1,     -6.77161551},
658                        {1,     -6.90069818},
659                        {1,     -7.10470295},
660                        {1,     -7.26322412},
661                        {1,     -7.39289951},
662                        {1,     -7.59792519},
663                        {1,     -7.75725508},
664                        {1,     -7.85722494},
665                        {1,     -7.91697407},
666                        {1,     -8.04758644},
667                        {1,     -8.09447479},
668                        {1,     -8.15859795},
669                        {1,     -8.25424385},
670                        {1,     -8.33880615},
671                        {1,     -8.41452408},
672                        {1,     -8.54581165},
673                        {1,     -8.60400581},
674                        {1,     -8.65751839},
675                        {1,     -8.75414848},
676                        {1,     -8.83946800},
677                        {1,     -8.91589737},
678                        {1,     -9.01741695},
679                        {1,     -9.10663033},
680                        {1,     -9.18621922},
681                        {1,     -9.28059292},
682                        {1,     -9.34430218},
683                        {1,     -9.42154408},
684                        {1,     -9.55562973},
685                        {1,     -9.61459446},
686                        {1,     -9.72023010},
687                        {1,     -9.76802444},
688                        {1,     -9.85540199},
689                        {1,     -9.93374062},
690                        {1,     -10.03800774},
691                        {1,     -10.10044670},
692                        {1,     -10.21178055},
693                        {1,     -10.29757786},
694                        {1,     -10.39561272},
695                        {1,     -10.45469666},
696                        {1,     -10.56102180},
697                        {1,     -10.66205502},
698                        {1,     -10.73780537},
699                        {1,     -10.82557774},
700                        {1,     -10.91007328},
701                        {1,     -10.98659325},
702                        {1,     -11.07857418},
703                        {1,     -11.19975281},
704                        {1,     -11.27170753},
705                        {1,     -11.36930943},
706                        {1,     -11.45675945},
707                        {1,     -11.53620148},
708                        {1,     -11.63198853},
709                        {1,     -11.73875237},
710                        {1,     -11.83400822},
711                        {1,     -11.93677044},
712                        {1,     -12.02933311},
713                        {1,     -12.13374519},
714                        {1,     -12.23410702},
715                        {1,     -12.33664989},
716                        {1,     -12.44163322},
717                        {1,     -12.54730415},
718                        {1,     -12.64975739},
719                        {1,     -12.76682186},
720                        {1,     -12.88185978},
721                        {1,     -13.00052643},
722                        {1,     -13.12289810},
723                        {1,     -13.26689529},
724                        {1,     -13.39390945},
725                        {1,     -30.00000000},
726                        {0L, 0 }
727                        };
728                for(_i=_r=0L; _i < 83; _i++)
729                        a0[_i] = RC_INI(_rs3); }
730                { static struct{ long rc; double ini; } _rs4[] = {
731                        {1,     4.53776000e05},
732                        {1,     3.48304000e05},
733                        {1,     2.80224000e05},
734                        {1,     1.98128000e05},
735                        {1,     1.51219797e05},
736                        {1,     1.21113266e05},
737                        {1,     8.52812109e04},
738                        {1,     6.49598125e04},
739                        {1,     5.20075781e04},
740                        {1,     3.66190977e04},
741                        {1,     2.79060723e04},
742                        {1,     2.23634102e04},
743                        {1,     1.57683135e04},
744                        {1,     1.20284307e04},
745                        {1,     9.17755273e03},
746                        {1,     7.36044873e03},
747                        {1,     5.19871680e03},
748                        {1,     3.97240796e03},
749                        {1,     3.18934326e03},
750                        {1,     2.25737622e03},
751                        {1,     1.72767114e03},
752                        {1,     1.46202722e03},
753                        {1,     1.32456628e03},
754                        {1,     1.06499792e03},
755                        {1,     9.92735291e02},
756                        {1,     8.91604858e02},
757                        {1,     7.59411560e02},
758                        {1,     6.59120056e02},
759                        {1,     5.80688965e02},
760                        {1,     4.66602264e02},
761                        {1,     4.27612854e02},
762                        {1,     3.91531494e02},
763                        {1,     3.34516968e02},
764                        {1,     2.91021820e02},
765                        {1,     2.56853912e02},
766                        {1,     2.17598007e02},
767                        {1,     1.88145462e02},
768                        {1,     1.65329865e02},
769                        {1,     1.41960342e02},
770                        {1,     1.28181686e02},
771                        {1,     1.13336761e02},
772                        {1,     9.17785034e01},
773                        {1,     8.36242981e01},
774                        {1,     7.08843536e01},
775                        {1,     6.58346100e01},
776                        {1,     5.75790634e01},
777                        {1,     5.11293755e01},
778                        {1,     4.37563019e01},
779                        {1,     3.99226875e01},
780                        {1,     3.39562836e01},
781                        {1,     3.00413170e01},
782                        {1,     2.61871891e01},
783                        {1,     2.41310368e01},
784                        {1,     2.08853607e01},
785                        {1,     1.82632275e01},
786                        {1,     1.60007000e01},
787                        {1,     1.42617064e01},
788                        {1,     1.27951088e01},
789                        {1,     1.16221066e01},
790                        {1,     1.03779335e01},
791                        {1,     8.97864914e00},
792                        {1,     8.25593281e00},
793                        {1,     7.39339924e00},
794                        {1,     6.70784378e00},
795                        {1,     6.16084862e00},
796                        {1,     5.57818031e00},
797                        {1,     5.01341105e00},
798                        {1,     4.55679178e00},
799                        {1,     4.13692093e00},
800                        {1,     3.80004382e00},
801                        {1,     3.46328306e00},
802                        {1,     3.17340493e00},
803                        {1,     2.93525696e00},
804                        {1,     2.69083858e00},
805                        {1,     2.46588683e00},
806                        {1,     2.26083040e00},
807                        {1,     2.04337358e00},
808                        {1,     1.89027369e00},
809                        {1,     1.69208312e00},
810                        {0L, 0 }
811                };
812                for(_i=_r=0L; _i < 79; _i++)
813                        ba0[_i] = (realnum)RC_INI(_rs4); }
814                { static struct{ long rc; double ini; } _rs5[] = {
815                        {1,     1.48992336e00},
816                        {1,     1.32466662e00},
817                        {1,     1.10697949e00},
818                        {1,     9.29813862e-01},
819                        {0L, 0 }
820                };
821                for(_i=_r=0L; _i < 4; _i++)
822                        bb0[_i] = (realnum)RC_INI(_rs5); }
823                { static struct{ long rc; double ini; } _rs6[] = {
824                        {1,     2.12597656},
825                        {1,     2.08984375},
826                        {1,     2.06958008},
827                        {1,     2.05444336},
828                        {1,     2.05},
829                        {1,     2.05},
830                        {1,     2.05},
831                        {1,     2.05},
832                        {1,     2.05},
833                        {1,     2.05},
834                        {1,     2.05},
835                        {1,     2.05},
836                        {1,     2.05},
837                        {1,     2.05},
838                        {1,     2.05},
839                        {1,     2.05},
840                        {1,     2.05},
841                        {1,     2.05},
842                        {1,     2.05},
843                        {1,     2.05},
844                        {1,     2.05},
845                        {1,     2.05},
846                        {1,     2.05},
847                        {1,     2.05},
848                        {1,     2.05},
849                        {1,     2.05},
850                        {1,     2.05},
851                        {1,     2.05},
852                        {1,     2.05},
853                        {1,     2.05},
854                        {1,     2.05},
855                        {1,     2.05},
856                        {1,     2.05},
857                        {1,     2.05},
858                        {1,     2.05},
859                        {1,     2.05},
860                        {1,     2.05},
861                        {1,     2.05},
862                        {1,     2.05},
863                        {1,     2.05},
864                        {1,     2.05},
865                        {1,     2.05},
866                        {1,     2.05},
867                        {1,     2.05},
868                        {1,     2.05},
869                        {1,     2.05},
870                        {1,     2.05},
871                        {1,     2.05},
872                        {1,     2.05},
873                        {1,     2.05},
874                        {1,     2.05},
875                        {1,     2.05},
876                        {1,     2.05},
877                        {1,     2.05},
878                        {1,     2.05},
879                        {1,     2.05},
880                        {1,     2.05},
881                        {1,     2.05},
882                        {1,     2.05},
883                        {1,     2.05},
884                        {1,     2.05},
885                        {1,     2.05},
886                        {1,     2.05},
887                        {1,     2.05},
888                        {1,     2.05},
889                        {1,     2.05},
890                        {1,     2.05},
891                        {1,     2.05},
892                        {1,     2.05},
893                        {1,     2.05},
894                        {1,     2.05},
895                        {1,     2.05},
896                        {1,     2.05},
897                        {1,     2.05},
898                        {1,     2.05},
899                        {1,     2.05},
900                        {1,     2.05},
901                        {1,     2.05},
902                        {1,     2.05},
903                        {1,     2.05},
904                        {1,     2.05},
905                        {1,     2.05},
906                        {1,     2.05},
907                        {0L, 0 }
908                };
909                for(_i=_r=0L; _i < 83; _i++)
910                        x0[_i] = RC_INI(_rs6); }
911
912                { static struct{ long rc; double ini; } _rs7[] = {
913                        {1,     16.23337936},
914                        {1,     16.27946854},
915                        {1,     16.31696320},
916                        {1,     16.37597656},
917                        {1,     16.42210960},
918                        {1,     16.45996284},
919                        {1,     16.51994896},
920                        {1,     16.56644440},
921                        {1,     16.60460854},
922                        {1,     16.66510773},
923                        {1,     16.71198654},
924                        {1,     16.75038719},
925                        {1,     16.81106949},
926                        {1,     16.85778809},
927                        {1,     16.90416527},
928                        {1,     16.94209099},
929                        {1,     17.00195694},
930                        {1,     17.04838943},
931                        {1,     17.08633804},
932                        {1,     17.14627838},
933                        {1,     17.19270515},
934                        {1,     17.22186279},
935                        {1,     17.23933601},
936                        {1,     17.27728271},
937                        {1,     17.30161858},
938                        {1,     17.32085800},
939                        {1,     17.34928894},
940                        {1,     17.37349129},
941                        {1,     17.39528084},
942                        {1,     17.43282318},
943                        {1,     17.44827652},
944                        {1,     17.46357536},
945                        {1,     17.49082375},
946                        {1,     17.51517677},
947                        {1,     17.53697205},
948                        {1,     17.56587219},
949                        {1,     17.59125519},
950                        {1,     17.61410332},
951                        {1,     17.64081383},
952                        {1,     17.65900803},
953                        {1,     17.68086433},
954                        {1,     17.71843529},
955                        {1,     17.73512840},
956                        {1,     17.76512146},
957                        {1,     17.77873421},
958                        {1,     17.80340767},
959                        {1,     17.82530022},
960                        {1,     17.85470963},
961                        {1,     17.87210464},
962                        {1,     17.90334511},
963                        {1,     17.92751503},
964                        {1,     17.95458603},
965                        {1,     17.97117233},
966                        {1,     18.00062943},
967                        {1,     18.02842712},
968                        {1,     18.04934502},
969                        {1,     18.07340050},
970                        {1,     18.09639168},
971                        {1,     18.11732864},
972                        {1,     18.14218903},
973                        {1,     18.17465591},
974                        {1,     18.19370079},
975                        {1,     18.21962166},
976                        {1,     18.24237251},
977                        {1,     18.26305962},
978                        {1,     18.28767967},
979                        {1,     18.31531525},
980                        {1,     18.33900452},
981                        {1,     18.36478043},
982                        {1,     18.38741112},
983                        {1,     18.41271973},
984                        {1,     18.43644333},
985                        {1,     18.46075630},
986                        {1,     18.48509216},
987                        {1,     18.50897980},
988                        {1,     18.53143501},
989                        {1,     18.55570030},
990                        {1,     18.58008003},
991                        {1,     18.60348320},
992                        {1,     18.62536430},
993                        {1,     18.65199852},
994                        {1,     18.67623520},
995                        {1,     18.70072174},
996                        {0L, 0 }
997                };
998                for(_i=_r=0L; _i < 83; _i++)
999                        a1[_i] = RC_INI(_rs7); }
1000                { static struct{ long rc; double ini; } _rs8[] = {
1001                        {1,     1.09462866e10},
1002                        {1,     9.32986675e09},
1003                        {1,     6.15947008e09},
1004                        {1,     1.54486170e09},
1005                        {1,     1.00812454e09},
1006                        {1,     7.00559552e08},
1007                        {1,     6.25999232e08},
1008                        {1,     3.50779968e08},
1009                        {1,     3.11956288e08},
1010                        {1,     3.74866016e08},
1011                        {1,     2.47019424e08},
1012                        {1,     1.73169776e08},
1013                        {1,     1.01753168e08},
1014                        {1,     6.81861920e07},
1015                        {1,     4.61764000e07},
1016                        {1,     3.31671360e07},
1017                        {1,     2.03160540e07},
1018                        {1,     1.40249480e07},
1019                        {1,     1.02577860e07},
1020                        {1,     3.53822650e06},
1021                        {1,     1.32563388e06},
1022                        {1,     9.14284688e05},
1023                        {1,     1.25230388e06},
1024                        {1,     3.17865156e05},
1025                        {1,     4.76750244e03},
1026                        {1,     4.81107031e03},
1027                        {1,     4.88406152e03},
1028                        {1,     4.80611279e03},
1029                        {1,     4.78843652e03},
1030                        {1,     4.65988477e03},
1031                        {1,     1.26723059e03},
1032                        {1,     1.20825342e03},
1033                        {1,     8.66052612e02},
1034                        {1,     7.76661316e02},
1035                        {1,     7.05279358e02},
1036                        {1,     6.21722656e02},
1037                        {1,     5.46207581e02},
1038                        {1,     4.96247742e02},
1039                        {1,     4.26340118e02},
1040                        {1,     3.96090424e02},
1041                        {1,     3.48429657e02},
1042                        {1,     2.37949142e02},
1043                        {1,     2.14678406e02},
1044                        {1,     1.81019180e02},
1045                        {1,     1.68923676e02},
1046                        {1,     1.45979385e02},
1047                        {1,     1.25311127e02},
1048                        {1,     1.05205528e02},
1049                        {1,     9.39378357e01},
1050                        {1,     7.75339966e01},
1051                        {1,     6.68987427e01},
1052                        {1,     5.53580055e01},
1053                        {1,     5.00100212e01},
1054                        {1,     4.14198608e01},
1055                        {1,     3.46289063e01},
1056                        {1,     3.00775223e01},
1057                        {1,     2.60294399e01},
1058                        {1,     2.26602840