cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
atmdat_3body.cpp
Go to the documentation of this file.
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 
15 STATIC void blkdata1(void);
16 STATIC double da(double z);
17 
18 static double a2[63],
19  b2[63],
20  x2[63];
21 
22 static double a0[83],
23  x0[83];
24 static realnum b0[83],
25  b1[83];
26 
27 static double a1[83],
28  x1[83];
29 
30 static 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 */
39 STATIC double xmap(double x[],
40  double y[],
41  double x0);
42 
43 /* inverse routine, also below */
44 STATIC double xinvrs(double y,
45  double a,
46  double b,
47  double u,
48  double v,
49  long int *ifail);
50 
51 /* =================================================================== */
52 void 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 
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 /* =================================================================== */
111 STATIC 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 /****************************************************************************** */
470 STATIC 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.26602840e01},
1059  {1, 2.02123032e01},
1060  {1, 1.76353855e01},
1061  {1, 1.47198439e01},
1062  {1, 1.33078461e01},
1063  {1, 1.17181997e01},
1064  {1, 1.04125805e01},
1065  {1, 9.45785904e00},
1066  {1, 8.42799950e00},
1067  {1, 7.62769842e00},
1068  {1, 6.85484743e00},
1069  {1, 6.25903368e00},
1070  {1, 5.75135279e00},
1071  {1, 5.28468180e00},
1072  {1, 4.87669659e00},
1073  {1, 4.57353973e00},
1074  {1, 4.30108690e00},
1075  {1, 4.05412245e00},
1076  {1, 3.83283114e00},
1077  {1, 3.57902861e00},
1078  {1, 3.43705726e00},
1079  {1, 3.26563096e00},
1080  {0L, 0 }
1081  };
1082  for(_i=_r=0L; _i < 79; _i++)
1083  ba1[_i] = (realnum)RC_INI(_rs8); }
1084  { static struct{ long rc; double ini; } _rs9[] = {
1085  {1, 3.07498097e00},
1086  {1, 2.96334076e00},
1087  {1, 2.92890000e00},
1088  {1, 2.89550042e00},
1089  {0L, 0 }
1090  };
1091  for(_i=_r=0L; _i < 4; _i++)
1092  bb1[_i] = (realnum)RC_INI(_rs9); }
1093  { static struct{ long rc; double ini; } _rs10[] = {
1094  {1, -5.46},
1095  {1, -5.51},
1096  {1, -5.49},
1097  {1, -5.30},
1098  {1, -5.29},
1099  {1, -5.28},
1100  {1, -5.37},
1101  {1, -5.33},
1102  {1, -5.38},
1103  {1, -5.55},
1104  {1, -5.55},
1105  {1, -5.55},
1106  {1, -5.55},
1107  {1, -5.55},
1108  {1, -5.55},
1109  {1, -5.55},
1110  {1, -5.55},
1111  {1, -5.55},
1112  {1, -5.55},
1113  {1, -5.38},
1114  {1, -5.19},
1115  {1, -5.14},
1116  {1, -5.27},
1117  {1, -4.93},
1118  {1, -3.64},
1119  {1, -3.68},
1120  {1, -3.74},
1121  {1, -3.78},
1122  {1, -3.82},
1123  {1, -3.88},
1124  {1, -3.40},
1125  {1, -3.41},
1126  {1, -3.32},
1127  {1, -3.32},
1128  {1, -3.32},
1129  {1, -3.32},
1130  {1, -3.31},
1131  {1, -3.31},
1132  {1, -3.29},
1133  {1, -3.29},
1134  {1, -3.27},
1135  {1, -3.16},
1136  {1, -3.14},
1137  {1, -3.11},
1138  {1, -3.10},
1139  {1, -3.07},
1140  {1, -3.03},
1141  {1, -2.99},
1142  {1, -2.96},
1143  {1, -2.91},
1144  {1, -2.87},
1145  {1, -2.81},
1146  {1, -2.78},
1147  {1, -2.72},
1148  {1, -2.66},
1149  {1, -2.61},
1150  {1, -2.56},
1151  {1, -2.51},
1152  {1, -2.47},
1153  {1, -2.42},
1154  {1, -2.35},
1155  {1, -2.31},
1156  {1, -2.26},
1157  {1, -2.21},
1158  {1, -2.17},
1159  {1, -2.12},
1160  {1, -2.08},
1161  {1, -2.03},
1162  {1, -1.99},
1163  {1, -1.95},
1164  {1, -1.91},
1165  {1, -1.87},
1166  {1, -1.84},
1167  {1, -1.81},
1168  {1, -1.78},
1169  {1, -1.75},
1170  {1, -1.71},
1171  {1, -1.69},
1172  {1, -1.66},
1173  {1, -1.62},
1174  {1, -1.60},
1175  {1, -1.60},
1176  {1, -1.60},
1177  {0L, 0 }
1178  };
1179  for(_i=_r=0L; _i < 83; _i++)
1180  x1[_i] = RC_INI(_rs10); }
1181  { static struct{ long rc; double ini; } _rs11[] = {
1182  {1, 20.30049515},
1183  {1, 20.28500366},
1184  {1, 20.25300407},
1185  {1, 20.16626740},
1186  {1, 20.15743256},
1187  {1, 20.11256981},
1188  {1, 20.04818344},
1189  {1, 19.99261856},
1190  {1, 19.94472885},
1191  {1, 19.86478043},
1192  {1, 19.83321571},
1193  {1, 19.80185127},
1194  {1, 19.74884224},
1195  {1, 19.70136070},
1196  {1, 19.65981102},
1197  {1, 19.60598755},
1198  {1, 19.56017494},
1199  {1, 19.52042389},
1200  {1, 19.47429657},
1201  {1, 19.44413757},
1202  {1, 19.40796280},
1203  {1, 19.34819984},
1204  {1, 19.32203293},
1205  {1, 19.27634430},
1206  {1, 19.25627136},
1207  {1, 19.22009087},
1208  {1, 19.18853378},
1209  {1, 19.14809799},
1210  {1, 19.12456703},
1211  {1, 19.08409119},
1212  {1, 19.05431557},
1213  {1, 19.02083015},
1214  {1, 19.00176430},
1215  {1, 18.96817970},
1216  {1, 18.93762589},
1217  {1, 18.91706085},
1218  {1, 18.89299583},
1219  {1, 18.87085915},
1220  {1, 18.85210609},
1221  {1, 18.83035851},
1222  {1, 18.80403900},
1223  {1, 18.78901100},
1224  {1, 18.77099228},
1225  {1, 18.75540161},
1226  {1, 18.74287033},
1227  {1, 18.72928810},
1228  {1, 18.71601868},
1229  {1, 18.70474434},
1230  {1, 18.69515800},
1231  {1, 18.68782425},
1232  {1, 18.68120766},
1233  {1, 18.67630005},
1234  {1, 18.67357445},
1235  {1, 18.67129898},
1236  {1, 18.67042351},
1237  {1, 18.67090988},
1238  {1, 18.67313004},
1239  {1, 18.67636490},
1240  {1, 18.68120003},
1241  {1, 18.68803024},
1242  {1, 18.69487381},
1243  {1, 18.70458412},
1244  {1, 18.71205139},
1245  {0L, 0 }
1246  };
1247  for(_i=_r=0L; _i < 63; _i++)
1248  a2[_i] = RC_INI(_rs11); }
1249  { static struct{ long rc; double ini; } _rs12[] = {
1250  {1, 1.01078403e00},
1251  {1, 1.97956896e00},
1252  {1, 3.14605665e00},
1253  {1, 6.46874905e00},
1254  {1, 3.16406364e01},
1255  {1, 3.74927673e01},
1256  {1, 4.75353088e01},
1257  {1, 5.27809143e01},
1258  {1, 5.86515846e01},
1259  {1, 6.70408707e01},
1260  {1, 1.14904137e02},
1261  {1, 1.03133133e02},
1262  {1, 1.26508728e02},
1263  {1, 1.03827606e02},
1264  {1, 8.79508896e01},
1265  {1, 7.18328934e01},
1266  {1, 6.19807892e01},
1267  {1, 5.51255455e01},
1268  {1, 4.87156143e01},
1269  {1, 4.58579826e01},
1270  {1, 4.19952011e01},
1271  {1, 4.08252220e01},
1272  {1, 3.78243637e01},
1273  {1, 3.34573860e01},
1274  {1, 3.19036102e01},
1275  {1, 2.92026005e01},
1276  {1, 2.74482193e01},
1277  {1, 2.54643936e01},
1278  {1, 2.46636391e01},
1279  {1, 2.33054180e01},
1280  {1, 2.23069897e01},
1281  {1, 2.12891216e01},
1282  {1, 2.06667900e01},
1283  {1, 1.96430798e01},
1284  {1, 1.87381802e01},
1285  {1, 1.76523514e01},
1286  {1, 1.69235287e01},
1287  {1, 1.62647285e01},
1288  {1, 1.56806908e01},
1289  {1, 1.50346069e01},
1290  {1, 1.42240467e01},
1291  {1, 1.37954988e01},
1292  {1, 1.31949224e01},
1293  {1, 1.27211905e01},
1294  {1, 1.22885675e01},
1295  {1, 1.17868662e01},
1296  {1, 1.12577572e01},
1297  {1, 1.08565578e01},
1298  {1, 1.04121590e01},
1299  {1, 1.00410652e01},
1300  {1, 9.64534473e00},
1301  {1, 9.29232121e00},
1302  {1, 8.92519569e00},
1303  {1, 8.60898972e00},
1304  {1, 8.31234550e00},
1305  {1, 8.04089737e00},
1306  {1, 7.74343491e00},
1307  {1, 7.48133039e00},
1308  {1, 7.21957016e00},
1309  {1, 6.94726801e00},
1310  {1, 6.71931219e00},
1311  {1, 6.45107985e00},
1312  {1, 6.28593779e00},
1313  {0L, 0 }
1314  };
1315  for(_i=_r=0L; _i < 63; _i++)
1316  b2[_i] = RC_INI(_rs12); }
1317  { static struct{ long rc; double ini; } _rs13[] = {
1318  {1, -0.43},
1319  {1, -0.75},
1320  {1, -0.93},
1321  {1, -1.20},
1322  {1, -1.78},
1323  {1, -1.85},
1324  {1, -1.95},
1325  {1, -2.00},
1326  {1, -2.05},
1327  {1, -2.12},
1328  {1, -2.34},
1329  {1, -2.31},
1330  {1, -2.42},
1331  {1, -2.36},
1332  {1, -2.31},
1333  {1, -2.25},
1334  {1, -2.21},
1335  {1, -2.18},
1336  {1, -2.15},
1337  {1, -2.14},
1338  {1, -2.12},
1339  {1, -2.14},
1340  {1, -2.12},
1341  {1, -2.09},
1342  {1, -2.08},
1343  {1, -2.06},
1344  {1, -2.05},
1345  {1, -2.04},
1346  {1, -2.04},
1347  {1, -2.04},
1348  {1, -2.04},
1349  {1, -2.04},
1350  {1, -2.04},
1351  {1, -2.04},
1352  {1, -2.04},
1353  {1, -2.04},
1354  {1, -2.04},
1355  {1, -2.04},
1356  {1, -2.04},
1357  {1, -2.04},
1358  {1, -2.04},
1359  {1, -2.04},
1360  {1, -2.04},
1361  {1, -2.04},
1362  {1, -2.04},
1363  {1, -2.04},
1364  {1, -2.04},
1365  {1, -2.04},
1366  {1, -2.04},
1367  {1, -2.04},
1368  {1, -2.04},
1369  {1, -2.04},
1370  {1, -2.04},
1371  {1, -2.04},
1372  {1, -2.04},
1373  {1, -2.04},
1374  {1, -2.04},
1375  {1, -2.04},
1376  {1, -2.04},
1377  {1, -2.04},
1378  {1, -2.04},
1379  {1, -2.04},
1380  {1, -2.04},
1381  {0L, 0 }
1382  };
1383  for(_i=_r=0L; _i < 63; _i++)
1384  x2[_i] = RC_INI(_rs13); }
1385  /*lint +e736 loss of precision in assignment in translated data */
1386 
1387  DEBUG_ENTRY( "blkdata0()" );
1388 }
1389 
1390 /* =================================================================== */
1391 /*xmap mapping function for Cota's 3-body recombination */
1392 STATIC double xmap(double x[],
1393  double y[],
1394  double xmapx0)
1395 {
1396  double a,
1397  b,
1398  c,
1399  xmapx1,
1400  x12m,
1401  x13m,
1402  xmapx2,
1403  x3,
1404  xmap_v,
1405  yinit,
1406  y13m;
1407 
1408  DEBUG_ENTRY( "xmap()" );
1409 
1410  /* PARABOLIC INTERPOLATION.
1411  * */
1412 
1413  yinit = y[0];
1414  xmapx1 = x[0];
1415  xmapx2 = x[1];
1416  x3 = x[2];
1417  x13m = xmapx1 - x3;
1418  x12m = xmapx1 - xmapx2;
1419  y13m = yinit - y[2];
1420  x3 = (xmapx1 + x3)*x13m;
1421  xmapx2 = (xmapx1 + xmapx2)*x12m;
1422  b = ((yinit - y[1])*x3 - y13m*xmapx2)/(x12m*x3 - x13m*xmapx2);
1423  a = (y13m - x13m*b)/x3;
1424  c = yinit - a*xmapx1*xmapx1 - b*xmapx1;
1425 
1426  xmap_v = a*xmapx0*xmapx0 + b*xmapx0 + c;
1427 
1428  return( xmap_v );
1429 }
1430 
1431 /* =================================================================== */
1432 /*xinvrs do inverse function for Cota's three-body recombination */
1433 STATIC double xinvrs(double y,
1434  double a,
1435  double b,
1436  double u,
1437  double v,
1438  long int *ifail)
1439 {
1440  long int i;
1441  double bxu,
1442  dfx,
1443  fx,
1444  fxdfx,
1445  x,
1446  xinvrs_v,
1447  xlog,
1448  xx;
1449  static long itmax = 10;
1450 
1451  DEBUG_ENTRY( "xinvrs()" );
1452 
1453  /* inverts equation of the form :
1454  *
1455  * Y = A + B * X ** U - V * LOG ( X )
1456  * */
1457  *ifail = 0;
1458  xlog = (a - y)/v;
1459  x = pow(10.,xlog);
1460  xx = 0.;
1461 
1462  for( i=0; i < itmax; i++ )
1463  {
1464  bxu = b*pow(x,u);
1465  fx = y - a - bxu + v*xlog;
1466  dfx = v*.4342945 - bxu*u;
1467 
1468  if( dfx != 0. )
1469  {
1470  fxdfx = fabs(fx/dfx);
1471  fxdfx = MIN2(0.2,fxdfx);
1472  xx = x*(1. - sign(fxdfx,fx/dfx));
1473  }
1474  else
1475  {
1476  /* >>chng 96 feb 02 this added in case dfx ever 0
1477  * suggested by Peter van Hoof */
1478  xx = x*(1. - sign(0.2,fx));
1479  }
1480 
1481  if( (fabs(xx-x)/x) < 1.e-4 )
1482  {
1483  xinvrs_v = xx;
1484  return( xinvrs_v );
1485  }
1486  else
1487  {
1488  x = xx;
1489  if( x < 1e-30 )
1490  {
1491  xinvrs_v = 100.;
1492  *ifail = 1;
1493  return( xinvrs_v );
1494  }
1495  xlog = log10(x);
1496  }
1497  }
1498  xinvrs_v = xx;
1499  *ifail = 1;
1500  return( xinvrs_v );
1501 }

Generated for cloudy by doxygen 1.8.3.1