StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StFtpcMagboltz1.cc
1 
2 // $Id: StFtpcMagboltz1.cc,v 1.2 2003/09/02 17:58:15 perev Exp $
3 //
4 // $Log: StFtpcMagboltz1.cc,v $
5 // Revision 1.2 2003/09/02 17:58:15 perev
6 // gcc 3.2 updates + WarnOff
7 //
8 // Revision 1.1 2000/12/20 08:44:02 jcs
9 // Replace pam/ftpc/fmg with maker
10 //
11 
12 #include "StFtpcMagboltz1.hh"
13 #include <math.h>
14 #include <stdio.h>
15 
16 StFtpcMagboltz1::StFtpcMagboltz1()
17 {
18  /* Table of constant values */
19 c__1 = 1;
20 c__4 = 4;
21 c__0 = 0;
22 c__3 = 3;
23 c_b12 = 10.;
24 }
25 
26 StFtpcMagboltz1::~StFtpcMagboltz1()
27 {
28 
29 }
30 
31 int StFtpcMagboltz1::magboltz_(float *e_magni__, float *b_magni__,
32 float *b_ang__, float *press, float *p_ar__,
33  float *p_co2__, float *p_ne__, float *p_he__, float *temper,
34 float *vdr, float *psiang, float *efin)
35 {
36  /* System generated locals */
37  int i__1;
38  double d__1;
39 
40  /* Local variables */
41  static double aold, anew;
42  static int last, itot, iout, j, l, m;
43  static double aconv, anorm;
44  static int ist;
45 
46 
47 /* ---------------------------------------------------------------------
48 */
49 /* FILL CALLPARS COMMON BLOCK WITH CALLUP PARAMETERS */
50  callpars_1.e_magnitude__ = *e_magni__;
51  printf("setting e_magnitude to %f\n", callpars_1.e_magnitude__);
52  callpars_1.b_magnitude__ = *b_magni__;
53  callpars_1.b_angle__ = *b_ang__;
54  callpars_1.pressure = *press;
55  callpars_1.perc_ar__ = *p_ar__;
56  callpars_1.perc_co2__ = *p_co2__;
57  callpars_1.perc_ne__ = *p_ne__;
58  callpars_1.perc_he__ = *p_he__;
59  callpars_1.temperature = *temper;
60  inpt_1.efinal = *efin;
61  last = 0;
62  setup_(&last);
63 L777:
64  mixer_();
65  prnter_();
66  bfield_(&inpt_1.nstep1);
67 
68 /* BACKWARD PROLONGATION AND GAUSS-SEIDEL ITERATION */
69 
70  i__1 = inpt_1.lhigh;
71  for (l = 1; l <= i__1; ++l) {
72  aold = (float)1.;
73  itot = 0;
74  iout = 0;
75  f0calc_(&c__0, &c__0, &anew, &l);
76  goto L100;
77 L5:
78  f0calc_(&c__1, &c__0, &anew, &l);
79 L100:
80  ++itot;
81  ++iout;
82  anorm = (d__1 = (float)1. - anew / aold, fabs(d__1));
83  printf("NEW SUM = %f OLD SUM = %f FRACTIONAL DIFFERENCE = %f ITERATION NUMBER = %d\n", anew, aold, anorm, itot);
84  aold = anew;
85  if (iout == inpt_1.nout) {
86  output_(&c__0);
87  }
88  if (iout == inpt_1.nout) {
89  iout = 0;
90  }
91  if (itot > inpt_1.itmax) {
92  goto L20;
93  }
94  if (anorm > inpt_1.conv) {
95  goto L5;
96  }
97  stepph_(&l);
98  if (inpt_1.idlong == 1) {
99  steppr_(&c__0, &l);
100  }
101  output_(&l);
102  if (inpt_1.i2type == 1) {
103  goto L50;
104  }
105  if (inpt_1.nitalp == 1 && inpt_1.alpha != (float)0.) {
106  goto L30;
107  }
108  goto L60;
109 
110 /* INCLUDE COLLISIONS OF 2ND KIND */
111 
112 L50:
113  itot = 0;
114  iout = 0;
115 L10:
116  f0calc_(&c__1, &inpt_1.i2type, &anew, &l);
117  stepph_(&l);
118  ++iout;
119  ++itot;
120  anorm = (d__1 = (float)1. - anew / aold, fabs(d__1));
121  printf("NEW SUM = %f OLD SUM = %f FRACTIONAL DIFFERENCE = %f ITERATION NUMBER = %d\n", anew, aold, anorm, itot);
122  aold = anew;
123  if (iout == inpt_1.nout) {
124  output_(&c__0);
125  }
126  if (iout == inpt_1.nout) {
127  iout = 0;
128  }
129  if (itot > inpt_1.itmax) {
130  goto L20;
131  }
132  if (anorm >= inpt_1.conv) {
133  goto L10;
134  }
135  if (inpt_1.idlong == 1) {
136  steppr_(&inpt_1.i2type, &l);
137  }
138  output_(&c__3);
139  goto L60;
140 
141 /* INCLUDE NET IONISATION AND ATTACHMENT TO ALL ORDERS */
142 
143 L30:
144  itot = 0;
145  ist = 0;
146 L40:
147  f0calc_(&c__1, &inpt_1.i2type, &anew, &l);
148  stepph_(&l);
149  ++itot;
150  anorm = (d__1 = (float)1. - anew / aold, fabs(d__1));
151  if (anorm < inpt_1.conv) {
152  printf("NEW SUM = %f OLD SUM = %f FRACTIONAL DIFFERENCE = %f ITERATION NUMBER = %d\n", anew, aold, anorm, itot);
153  }
154  aold = anew;
155  if (itot > inpt_1.itmax) {
156  goto L20;
157  }
158  if (anorm >= inpt_1.conv) {
159  goto L40;
160  }
161 
162 /* ITERATE ON NET IONISATION */
163 
164  nalpha_();
165  ++ist;
166  if (inpt_1.alpold == (float)0.) {
167  inpt_1.alpold = 1e-10;
168  }
169  if (inpt_1.alpnew < (float)700.) {
170  goto L876;
171  }
172  if (ist == 1) {
173  inpt_1.alpnew *= (float).5;
174  }
175  if (ist == 1) {
176  inpt_1.alpnax *= (float).5;
177  }
178  if (ist == 1) {
179  inpt_1.alpnay *= (float).5;
180  }
181  if (ist == 1) {
182  inpt_1.alpnaz *= (float).5;
183  }
184 L876:
185  aconv = (d__1 = (float)1. - inpt_1.alpnew / inpt_1.alpold,
186 fabs(d__1));
187  if (aconv >= inpt_1.conalp) {
188  goto L40;
189  }
190  if (inpt_1.idlong == 1) {
191  steppr_(&inpt_1.i2type, &l);
192  }
193  output_(&c__4);
194 
195 /* CALCULATE HIGHER MULTIPOLE TERM IN F-DISTRIBUTION TO FIRST ORDER
196 */
197 
198 /* REDUCED OUTPUT */
199 L60:
200  if (l == 2) {
201  goto L2;
202  }
203  if (l == 1) {
204  goto L3000;
205  }
206 /* REDUCED OUTPUT */
207  fncalc_(&l);
208  stepph_(&l);
209  if (inpt_1.idlong == 1) {
210  steppr_(&inpt_1.i2type, &l);
211  }
212  m = l + 4;
213 /* MK CALL OUTPUT(M) */
214  for (j = 1; j <= 2002; ++j) {
215  f2c_1.f2[j - 1] = (float)0.;
216 /* L433: */
217  f2c_1.df2[j - 1] = (float)0.;
218  }
219 
220 L3000:
221  ;
222  }
223  goto L2;
224 L20:
225  printf("NUMBER OF ITERATIONS GREATER THAN ALLOWED = %d\n", itot);
226 /* ********** MK new 99-routine ********** */
227 L2:
228  if (f0c_1.f[inpt_1.nstep - 1] > (float)1e-7) {
229  printf("Jump 1: %f\n", f0c_1.f[inpt_1.nstep - 1]);
230  inpt_1.efinal *= (float)1.1;
231 /* increase by 10% */
232  last = 20;
233  setup_(&last);
234  last = 0;
235  goto L777;
236  } else if (f0c_1.f[inpt_1.nstep - 1] < (float)1e-10) {
237  printf("Jump 1: %f\n", f0c_1.f[inpt_1.nstep - 1]);
238  inpt_1.efinal *= (float).9;
239 /* decrease by 10% */
240  last = 20;
241  setup_(&last);
242  last = 0;
243  goto L777;
244  } else {
245  printf("********************\n");
246  d__1 = mk_1.velz / (float)1e6;
247  printf("Drift velocity == %f\n", d__1);
248  printf("Longitudinal diff = %f\n",mk_1.sig_long__);
249  printf("Transvers_xx diff = %f\n", mk_1.sig_tranx__);
250  printf("Transvers_yy diff = %f\n", mk_1.sig_trany__);
251  printf("Lorentz angle == %f\n",mk_1.angle);
252  printf("Drift Field = %f\n", mk_1.amk_emag__);
253  printf("********************\n");
254  *vdr = mk_1.velz / (float)1e6;
255  *psiang = mk_1.angle;
256  *efin = inpt_1.efinal;
257  }
258  return 0;
259 } /* magboltz_ */
260 
261 
262 int StFtpcMagboltz1::bfield_(int *nstep1)
263 {
264  /* System generated locals */
265  int i__1;
266  double d__1, d__2;
267 
268  /* Local variables */
269  static double bcos, bsin, emag2;
270  static int j;
271  static double en, pi, atheta, ef2, sincos, cos2, sin2;
272 
273  pi =(double) acos(-1.);
274  atheta = mag_1.btheta * pi / 180.;
275  bsin = sin(atheta);
276  bcos = cos(atheta);
277 
278 /* REMOVE ROUNDING ERRORS IN SIN AND COS AT 90 DEGREES. */
279 
280  if (mag_1.btheta == (float)90.) {
281  bsin = (float)1.;
282  }
283  if (mag_1.btheta == (float)90.) {
284  bcos = (float)0.;
285  }
286  sin2 = bsin * bsin;
287  cos2 = bcos * bcos;
288  sincos = bsin * bcos;
289  emag2 = mag_1.emag * mag_1.emag;
290  i__1 = *nstep1 + 1;
291  for (j = 1; j <= i__1; ++j) {
292  en = mix2_1.e[j - 1];
293  if (en == (float)0.) {
294  en = (float)1e-4;
295  }
296 /* Computing 2nd power */
297  d__1 = mag_1.bmag;
298 /* Computing 2nd power */
299  d__2 = mix2_1.qtot[j - 1];
300  mag_1.denom[j - 1] = cnsts_1.echarg * (d__1 * d__1) / (cnsts_1.emass *
301  2e6 * en * (d__2 * d__2));
302  mag_1.sod[j - 1] = ::sqrt(mag_1.denom[j - 1]) * bsin;
303  mag_1.cod2[j - 1] = mag_1.denom[j - 1] * cos2;
304  mag_1.sod2[j - 1] = mag_1.denom[j - 1] * sin2;
305  mag_1.scd[j - 1] = mag_1.denom[j - 1] * sincos;
306  mag_1.cod2[j - 1] += (float)1.;
307  mag_1.sod2[j - 1] += (float)1.;
308  mag_1.denom[j - 1] += (float)1.;
309  mag_1.sod[j - 1] /= mag_1.denom[j - 1];
310  mag_1.sod2[j - 1] /= mag_1.denom[j - 1];
311  mag_1.cod2[j - 1] /= mag_1.denom[j - 1];
312  mag_1.scd[j - 1] /= mag_1.denom[j - 1];
313  ef2 = emag2 * cos2 + emag2 * sin2 / mag_1.denom[j - 1];
314  mag_1.ef[j - 1] = ef2 / mag_1.emag;
315  mag_1.qe[j - 1] = mag_1.ef[j - 1] / mix2_1.qtot[j - 1];
316  mag_1.qef[j - 1] = en / (mix2_1.qtot[j - 1] * (float)3.);
317  mag_1.qfemag[j - 1] = mag_1.qef[j - 1] * mag_1.emag;
318  mag_1.qeef[j - 1] = mag_1.qef[j - 1] * mag_1.ef[j - 1];
319 /* L100: */
320  mag_1.qeeef[j - 1] = ef2 * mag_1.qef[j - 1];
321  }
322  return 0;
323 } /* bfield_ */
324 
325 int StFtpcMagboltz1::f0calc_(int *iback, int *icon, double *anew, int *l)
326 {
327  /* System generated locals */
328  int i__1, i__2;
329 
330 
331  /* Local variables */
332  static int mion;
333  static double asum;
334  static double fmus, ssum, f1sum, f2sum;
335  static int i__, j, m;
336  static double s2sum, z__, f2sum1, f2sum2, fnsum, fnssm1, fnssm2, ai,
337  sfb, f2t1, f2t2, ext, sum;
338  static int min1, min2, min3, min4;
339 
340 
341 
342 
343 /* STARTING VALUES */
344 
345 
346  printf("starting f0calc\n");
347  if (*iback == 1) {
348  goto L10;
349  }
350  if (*l != 1) {
351  goto L10;
352  }
353  f0c_1.df[inpt_1.nstep + 1] = (float)0.;
354  f0c_1.f[inpt_1.nstep1 - 1] = 1e-7;
355 L10:
356  ssum = (float)0.;
357  f2sum = (float)0.;
358  fnssm1 = (float)0.;
359  fnssm2 = (float)0.;
360  z__ = (float)0.;
361  i__1 = inpt_1.nstep;
362  for (m = 1; m <= i__1; ++m) {
363  i__ = inpt_1.nstep - m + 2;
364  ai = (double) i__;
365  i__2 = inpt_1.ngas;
366  for (j = 1; j <= i__2; ++j) {
367  mion = mix3_1.lion[j - 1] + i__;
368  if (mion >= inpt_1.nstep1) {
369  goto L60;
370  }
371  z__ = z__ + f0c_1.f[mion - 1] * mix1_1.qion[j + (mion << 2) - 5] /
372  ai + (f0c_1.f[mion] * mix1_1.qion[j + ((mion + 1) << 2) -
373  5] - f0c_1.f[mion - 1] * mix1_1.qion[j + (mion << 2) - 5])
374  * mix3_1.alion[j - 1] / ai;
375 L60:
376  ;
377  }
378 
379  asum = f0c_1.f[i__ - 1] * (mag_1.qef[i__ - 1] * inpt_1.alpnaz *
380  inpt_1.alpnew + mag_1.qeef[i__ - 1] * inpt_1.alpnew *
381  f0c_1.df[i__ - 1]);
382  sum = z__ - f0c_1.f[i__ - 1] * mix1_1.qsum[i__ - 1] + asum;
383 
384  if (mix3_1.nin1 == 0) {
385  goto L610;
386  }
387  i__2 = mix3_1.nin1;
388  for (j = 1; j <= i__2; ++j) {
389  min1 = mix3_1.lin1[j - 1] + i__;
390  if (min1 >= inpt_1.nstep1) {
391  goto L61;
392  }
393  sum = sum + f0c_1.f[min1 - 1] * mix1_1.qin1[j + min1 * 24 - 25] +
394  (f0c_1.f[min1] * mix1_1.qin1[j + (min1 + 1) * 24 - 25] -
395  f0c_1.f[min1 - 1] * mix1_1.qin1[j + min1 * 24 - 25]) *
396  mix3_1.alin1[j - 1];
397 L61:
398  ;
399  }
400 L610:
401  if (mix3_1.nin2 == 0) {
402  goto L620;
403  }
404  i__2 = mix3_1.nin2;
405  for (j = 1; j <= i__2; ++j) {
406  min2 = mix3_1.lin2[j - 1] + i__;
407  if (min2 >= inpt_1.nstep1) {
408  goto L62;
409  }
410  sum = sum + f0c_1.f[min2 - 1] * mix1_1.qin2[j + min2 * 24 - 25] +
411  (f0c_1.f[min2] * mix1_1.qin2[j + (min2 + 1) * 24 - 25] -
412  f0c_1.f[min2 - 1] * mix1_1.qin2[j + min2 * 24 - 25]) *
413  mix3_1.alin2[j - 1];
414 L62:
415  ;
416  }
417 L620:
418  if (mix3_1.nin3 == 0) {
419  goto L630;
420  }
421  i__2 = mix3_1.nin3;
422  for (j = 1; j <= i__2; ++j) {
423  min3 = mix3_1.lin3[j - 1] + i__;
424  if (min3 >= inpt_1.nstep1) {
425  goto L63;
426  }
427  sum = sum + f0c_1.f[min3 - 1] * mix1_1.qin3[j + min3 * 24 - 25] +
428  (f0c_1.f[min3] * mix1_1.qin3[j + (min3 + 1) * 24 - 25] -
429  f0c_1.f[min3 - 1] * mix1_1.qin3[j + min3 * 24 - 25]) *
430  mix3_1.alin3[j - 1];
431 L63:
432  ;
433  }
434 L630:
435  if (mix3_1.nin4 == 0) {
436  goto L640;
437  }
438  i__2 = mix3_1.nin4;
439  for (j = 1; j <= i__2; ++j) {
440  min4 = mix3_1.lin4[j - 1] + i__;
441  if (min4 >= inpt_1.nstep1) {
442  goto L64;
443  }
444  sum = sum + f0c_1.f[min4 - 1] * mix1_1.qin4[j + min4 * 24 - 25] +
445  (f0c_1.f[min4] * mix1_1.qin4[j + (min4 + 1) * 24 - 25] -
446  f0c_1.f[min4 - 1] * mix1_1.qin4[j + min4 * 24 - 25]) *
447  mix3_1.alin4[j - 1];
448 L64:
449  ;
450  }
451 L640:
452 
453 /* COLLISIONS OF 2ND KIND */
454 
455  if (*icon != 1) {
456  goto L90;
457  }
458  type2_(&s2sum, &i__, &inpt_1.nstep1);
459  sum += s2sum;
460 L90:
461  fnsum = sum - asum;
462  fnsum *= inpt_1.estep;
463  fnssm1 += fnsum;
464  sum *= inpt_1.estep;
465  ssum += sum;
466  f2sum1 = inpt_1.alpnaz * (float).4 * inpt_1.alpnew * f2c_1.f2[i__ - 1]
467  * mag_1.qef[j - 1] * inpt_1.estep;
468  f2sum2 = inpt_1.alpnew * (float).4 * mag_1.qeef[i__ - 1] * (f2c_1.df2[
469  i__ - 1] + f2c_1.f2[i__ - 1] * (float)1.5 / mix2_1.e[i__ - 1])
470  * inpt_1.estep;
471  f2sum = f2sum + f2sum1 + f2sum2;
472  f2t1 = mag_1.qeeef[i__ - 1] * (float).4 * (f2c_1.df2[i__ - 1] +
473  f2c_1.f2[i__ - 1] * (float)1.5 / mix2_1.e[i__ - 1]);
474  f2t2 = mag_1.qfemag[i__ - 1] * (float).4 * inpt_1.alpnaz * f2c_1.f2[
475  i__ - 1];
476  if (inpt_1.isfb == 1) {
477  sfb = mag_1.qeeef[i__ - 1] * f0c_1.f[i__ - 1] + ((mix1_1.qelm[i__
478  - 1] + mag_1.qfemag[i__ - 1] * inpt_1.alpnaz) * f0c_1.f[
479  i__ - 1] - ssum) * inpt_1.akt;
480  }
481  if (inpt_1.isfb == 0) {
482  sfb = mag_1.qeeef[i__ - 1] * f0c_1.f[i__ - 1] + mix1_1.qelm[i__ -
483  1] * f0c_1.f[i__ - 1] * inpt_1.akt;
484  }
485  f0c_1.df[i__ - 1] = (ssum - f0c_1.f[i__ - 1] * (mix1_1.qelm[i__ - 1]
486  + mag_1.qfemag[i__ - 1] * inpt_1.alpnaz) + f2sum - f2t1 -
487  f2t2) / sfb;
488 
489  ext = mix1_1.qelm[i__ - 1] * f0c_1.f[i__ - 1] + mix1_1.qelm[i__ - 1] *
490  f0c_1.f[i__ - 1] * f0c_1.df[i__ - 1] * inpt_1.akt;
491  f1c_1.f1[i__ - 1] = (fnssm2 + ext - fnssm1) * (float)3. / (mag_1.emag
492  * mix2_1.e[i__ - 1]);
493  f1sum = inpt_1.alpnew * (mix2_1.e[i__ - 1] * (float)2. * f1c_1.f1[i__
494  - 1] - mix2_1.e[i__] * f1c_1.f1[i__]) * inpt_1.estep / (float)
495  3.;
496  fnssm2 += f1sum;
497  if (i__ >= inpt_1.nstep1 - 2) {
498  goto L623;
499  }
500  f1c_1.df1[i__ - 1] = (f1c_1.f1[i__ + 1] - f1c_1.f1[i__ - 1]) /
501  inpt_1.estep - f1c_1.df1[i__ + 1];
502  f1c_1.df1[i__] = (f1c_1.f1[i__ + 1] - f1c_1.f1[i__ - 1]) / (
503  inpt_1.estep * (float)2.);
504  goto L624;
505 L623:
506  if (i__ == inpt_1.nstep1) {
507  f1c_1.df1[i__ - 1] = -f1c_1.f1[inpt_1.nstep1 - 1] / (inpt_1.estep
508  * (float)2.);
509  }
510  if (i__ == inpt_1.nstep1 - 1) {
511  f1c_1.df1[i__ - 1] = (f1c_1.f1[inpt_1.nstep1 - 1] - f1c_1.f1[
512  inpt_1.nstep1 - 2]) / inpt_1.estep;
513  }
514  if (i__ == inpt_1.nstep1 - 2) {
515  f1c_1.df1[i__ - 1] = (f1c_1.f1[inpt_1.nstep1 - 2] - f1c_1.f1[
516  inpt_1.nstep1 - 3]) / inpt_1.estep;
517  }
518 
519 L624:
520  if (*l == 1) {
521  goto L629;
522  }
523  f2c_1.f2[i__ - 1] = mag_1.qe[i__ - 1] * (float)-2. * (f1c_1.df1[i__ -
524  1] - f1c_1.f1[i__ - 1] / (mix2_1.e[i__ - 1] * (float)2.)) / (
525  float)3. - inpt_1.alpnaz * (float)2. * f1c_1.f1[i__ - 1] / (
526  mix2_1.qtot[i__ - 1] * (float)3.);
527  f2c_1.df2[i__ - 1] = f0c_1.df[i__ - 1] * (float)-2.5 * f0c_1.f[i__ -
528  1] - inpt_1.alpnaz * (float)2.5 * f0c_1.f[i__ - 1] / mag_1.ef[
529  i__ - 1] - f2c_1.f2[i__ - 1] * (float)1.5 / mix2_1.e[i__ - 1]
530  - inpt_1.alpnaz * f2c_1.f2[i__ - 1] / mag_1.ef[i__ - 1] -
531  f1c_1.f1[i__ - 1] * (float)2.5 / mag_1.qe[i__ - 1];
532  f2c_1.df2[inpt_1.nstep1] = f2c_1.df2[inpt_1.nstep1 - 1];
533  f2c_1.f2[i__ - 2] = f2c_1.f2[i__ - 1] - inpt_1.estep * (f2c_1.df2[i__
534  - 1] * (float)1.5 - f2c_1.df2[i__] * (float).5);
535  f2c_1.df2[i__ - 2] = f2c_1.df2[i__ - 1] * (float)2. - f2c_1.df2[i__];
536  goto L296;
537 L629:
538  if (*iback == 1) {
539  goto L1000;
540  }
541 
542 
543 /* BACKWARD PROLONGATION */
544 
545 
546 L296:
547  f0c_1.f[i__ - 2] = f0c_1.f[i__ - 1] * exp(-(f0c_1.df[i__ - 1] * (
548  float)3. - f0c_1.df[i__]) * (float).5 * inpt_1.estep);
549 L1000:
550  ;
551  }
552  f0c_1.df[0] = (float)0.;
553 
554 
555 /* GAUSS-SEIDEL ITERATION */
556 
557 
558  f0c_1.f[1] = exp((f0c_1.df[0] + f0c_1.df[1]) * (float).5 * inpt_1.estep);
559  f0c_1.f[0] = f0c_1.f[1] * exp(-(f0c_1.df[1] * (float)3. - f0c_1.df[2]) * (
560  float).5 * inpt_1.estep);
561  f0c_1.f[2] = exp(inpt_1.estep * (f0c_1.df[0] + f0c_1.df[1] * (float)4. +
562  f0c_1.df[2]) / (float)3.);
563  i__1 = inpt_1.nstep1;
564  for (i__ = 5; i__ <= i__1; i__ += 2) {
565 /* L300: */
566  f0c_1.f[i__ - 1] = f0c_1.f[i__ - 3] * exp(inpt_1.estep * (f0c_1.df[
567  i__ - 3] + f0c_1.df[i__ - 2] * (float)4. + f0c_1.df[i__ - 1])
568  / (float)3.);
569  }
570  i__1 = inpt_1.nstep;
571  for (i__ = 4; i__ <= i__1; i__ += 2) {
572 /* L310: */
573  f0c_1.f[i__ - 1] = f0c_1.f[i__ - 3] * exp(inpt_1.estep * (f0c_1.df[
574  i__ - 3] + f0c_1.df[i__ - 2] * (float)4. + f0c_1.df[i__ - 1])
575  / (float)3.);
576  }
577 
578 
579 
580 /* NORMALIZATION OF PROBABILITY DISTRIBUTION */
581 
582 
583  i__1 = inpt_1.nstep1;
584  for (j = 1; j <= i__1; ++j) {
585 /* L400: */
586  sint_1.simf[j - 1] = f0c_1.f[j - 1] * mix2_1.eroot[j - 1];
587  }
588  simp_(anew);
589  fmus = (float)1. / *anew;
590  i__1 = inpt_1.nstep1;
591  for (j = 1; j <= i__1; ++j) {
592  f0c_1.f[j - 1] *= fmus;
593 /* L410: */
594  f0c_1.df0[j - 1] = f0c_1.df[j - 1] * f0c_1.f[j - 1];
595  }
596 
597  printf("done with f0calc %f\n", f0c_1.f[inpt_1.nstep - 1]);
598  return 0;
599 
600 } /* f0calc_ */
601 
602 int StFtpcMagboltz1::fncalc_(int *lmax)
603 {
604  /* System generated locals */
605  int i__1;
606 
607  /* Local variables */
608  static int j;
609 
610  if (*lmax == 2) {
611  goto L300;
612  }
613  f2c_1.f2[inpt_1.nstep1] = (float)0.;
614  i__1 = inpt_1.nstep1;
615  for (j = 2; j <= i__1; ++j) {
616 /* L100: */
617  f2c_1.f2[j - 1] = -(mag_1.ef[j - 1] * (float)2. * (mix2_1.e[j - 1] *
618  f1c_1.df1[j - 1] - f1c_1.f1[j - 1] / (float)2.) / (
619  mix2_1.qtot[j - 1] * (float)3. * mix2_1.e[j - 1]) + f1c_1.f1[
620  j - 1] * (float)2. * inpt_1.alpnaz / (mix2_1.qtot[j - 1] * (
621  float)3.));
622  }
623  f2c_1.f2[0] = (float)0.;
624  i__1 = inpt_1.nstep;
625  for (j = 2; j <= i__1; ++j) {
626 /* L150: */
627  f2c_1.df2[j - 1] = (f2c_1.f2[j] - f2c_1.f2[j - 2]) / (inpt_1.estep * (
628  float)2.);
629  }
630  f2c_1.df2[inpt_1.nstep1 - 1] = (f2c_1.f2[inpt_1.nstep1 - 1] - f2c_1.f2[
631  inpt_1.nstep - 1]) / inpt_1.estep;
632  f2c_1.df2[0] = f2c_1.f2[1] / inpt_1.estep;
633  return 0;
634 L300:
635  f3c_1.f3[inpt_1.nstep1] = (float)0.;
636  i__1 = inpt_1.nstep1;
637  for (j = 2; j <= i__1; ++j) {
638 /* L310: */
639  f3c_1.f3[j - 1] = -mag_1.qe[j - 1] * (float).6 * (f2c_1.df2[j - 1] -
640  f2c_1.f2[j - 1] / mix2_1.e[j - 1]) - inpt_1.alpnaz * (float)
641  .6 * f2c_1.f2[j - 1] / mix2_1.qtot[j - 1];
642  }
643  f3c_1.f3[0] = (float)0.;
644  i__1 = inpt_1.nstep;
645  for (j = 2; j <= i__1; ++j) {
646 /* L350: */
647  f3c_1.df3[j - 1] = (f3c_1.f3[j] - f3c_1.f3[j - 2]) / (inpt_1.estep * (
648  float)2.);
649  }
650  f3c_1.df3[inpt_1.nstep1 - 1] = (f3c_1.f3[inpt_1.nstep1 - 1] - f3c_1.f3[
651  inpt_1.nstep - 1]) / inpt_1.estep;
652  f3c_1.df3[0] = f3c_1.f3[1] / inpt_1.estep;
653  return 0;
654 } /* fncalc_ */
655 
656 int StFtpcMagboltz1::g0calc_(int *icon, double *gfinal, double *eg0sum, int *lmax)
657 {
658  /* System generated locals */
659  int i__1, i__2;
660 
661  /* Local variables */
662  static double anew;
663  static int mion;
664  static double asum, bsum, ssum, g1sum;
665  static int i__, j;
666  static double s2sum;
667  static int m;
668  static double z__, gssum, gnsum, ai;
669  static double gnssum, vel, ext, sum;
670  static int min1, min2, min3, min4;
671  static double sum1, sum2, sum3, sum4;
672 
673  inpt_1.alpnew = (float)0.;
674  inpt_1.alpnaz = (float)0.;
675  i__1 = inpt_1.nstep1;
676  for (j = 1; j <= i__1; ++j) {
677 /* L10: */
678  sint_1.simf[j - 1] = mix2_1.e[j - 1] * f1c_1.f1[j - 1];
679  }
680  simp_(&sum);
681  vel = sum * mag_1.eovm / (float)3.;
682 
683 /* STARTING VALUES */
684 
685  g0c_1.g[inpt_1.nstep1 - 1] = *gfinal;
686  g2c_1.g2[inpt_1.nstep1 - 1] = (float)0.;
687  g2c_1.dg2[inpt_1.nstep1 - 1] = (float)0.;
688  ssum = (float)0.;
689  sum4 = (float)0.;
690  gnssum = (float)0.;
691  gssum = (float)0.;
692  z__ = (float)0.;
693  i__1 = inpt_1.nstep;
694  for (m = 1; m <= i__1; ++m) {
695  i__ = inpt_1.nstep - m + 2;
696  ai = (double) i__;
697  i__2 = inpt_1.ngas;
698  for (j = 1; j <= i__2; ++j) {
699  mion = mix3_1.lion[j - 1] + i__;
700  if (mion >= inpt_1.nstep1) {
701  goto L60;
702  }
703  z__ = z__ + g0c_1.g[mion - 1] * mix1_1.qion[j + (mion << 2) - 5] /
704  ai + (g0c_1.g[mion] * mix1_1.qion[j + ((mion + 1) << 2) -
705  5] - g0c_1.g[mion - 1] * mix1_1.qion[j + (mion << 2) - 5])
706  * mix3_1.alion[j - 1] / ai;
707 L60:
708  ;
709  }
710 
711  asum = g0c_1.g[i__ - 1] * (mag_1.qef[i__ - 1] * inpt_1.alpnaz *
712  inpt_1.alpnew + mag_1.qeef[i__ - 1] * inpt_1.alpnew *
713  g0c_1.dg[i__ - 1]);
714  sum = z__ - g0c_1.g[i__ - 1] * mix1_1.qsum[i__ - 1] + asum;
715  if (*lmax == 1) {
716  goto L601;
717  }
718  bsum = inpt_1.alpnew * (inpt_1.alpnaz * mag_1.qef[i__ - 1] * g2c_1.g2[
719  i__ - 1] + mag_1.qeef[i__ - 1] * (g2c_1.dg2[i__ - 1] +
720  g2c_1.g2[i__ - 1] * (float)1.5 / mix2_1.e[i__ - 1]));
721  sum += bsum * (float).4;
722 
723 L601:
724  if (mix3_1.nin1 == 0) {
725  goto L610;
726  }
727  i__2 = mix3_1.nin1;
728  for (j = 1; j <= i__2; ++j) {
729  min1 = mix3_1.lin1[j - 1] + i__;
730  if (min1 >= inpt_1.nstep1) {
731  goto L61;
732  }
733  sum = sum + g0c_1.g[min1 - 1] * mix1_1.qin1[j + min1 * 24 - 25] +
734  (g0c_1.g[min1] * mix1_1.qin1[j + (min1 + 1) * 24 - 25] -
735  g0c_1.g[min1 - 1] * mix1_1.qin1[j + min1 * 24 - 25]) *
736  mix3_1.alin1[j - 1];
737 L61:
738  ;
739  }
740 L610:
741  if (mix3_1.nin2 == 0) {
742  goto L620;
743  }
744  i__2 = mix3_1.nin2;
745  for (j = 1; j <= i__2; ++j) {
746  min2 = mix3_1.lin2[j - 1] + i__;
747  if (min2 >= inpt_1.nstep1) {
748  goto L62;
749  }
750  sum = sum + g0c_1.g[min2 - 1] * mix1_1.qin2[j + min2 * 24 - 25] +
751  (g0c_1.g[min2] * mix1_1.qin2[j + (min2 + 1) * 24 - 25] -
752  g0c_1.g[min2 - 1] * mix1_1.qin2[j + min2 * 24 - 25]) *
753  mix3_1.alin2[j - 1];
754 L62:
755  ;
756  }
757 L620:
758  if (mix3_1.nin3 == 0) {
759  goto L630;
760  }
761  i__2 = mix3_1.nin3;
762  for (j = 1; j <= i__2; ++j) {
763  min3 = mix3_1.lin3[j - 1] + i__;
764  if (min3 >= inpt_1.nstep1) {
765  goto L63;
766  }
767  sum = sum + g0c_1.g[min3 - 1] * mix1_1.qin3[j + min3 * 24 - 25] +
768  (g0c_1.g[min3] * mix1_1.qin3[j + (min3 + 1) * 24 - 25] -
769  g0c_1.g[min3 - 1] * mix1_1.qin3[j + min3 * 24 - 25]) *
770  mix3_1.alin3[j - 1];
771 L63:
772  ;
773  }
774 L630:
775  if (mix3_1.nin4 == 0) {
776  goto L640;
777  }
778  i__2 = mix3_1.nin4;
779  for (j = 1; j <= i__2; ++j) {
780  min4 = mix3_1.lin4[j - 1] + i__;
781  if (min4 >= inpt_1.nstep1) {
782  goto L64;
783  }
784  sum = sum + g0c_1.g[min4 - 1] * mix1_1.qin4[j + min4 * 24 - 25] +
785  (g0c_1.g[min4] * mix1_1.qin4[j + (min4 + 1) * 24 - 25] -
786  g0c_1.g[min4 - 1] * mix1_1.qin4[j + min4 * 24 - 25]) *
787  mix3_1.alin4[j - 1];
788 L64:
789  ;
790  }
791 L640:
792  sum1 = -mag_1.qef[i__ - 1] * inpt_1.alpnew * (f0c_1.f[i__ - 1] +
793  f2c_1.f2[i__ - 1] * (float).4) + inpt_1.alpnew * mix2_1.eroot[
794  i__ - 1] * vel * f1c_1.f1[i__ - 1] / (mag_1.eovm * (float)3. *
795  mix2_1.qtot[i__ - 1]) + mix2_1.e[i__ - 1] * f1c_1.f1[i__ - 1]
796  / (float)3. - mix2_1.eroot[i__ - 1] * vel * f0c_1.f[i__ - 1]
797  / mag_1.eovm;
798  gnsum = sum - asum;
799  if (*lmax == 1) {
800  goto L641;
801  }
802  gnsum -= bsum * (float).4;
803 L641:
804  gnsum = gnsum + mix2_1.e[i__ - 1] * f1c_1.f1[i__ - 1] / (float)3. -
805  mix2_1.eroot[i__ - 1] * vel * f0c_1.f[i__ - 1] / mag_1.eovm;
806  sum += sum1;
807 
808 /* COLLISIONS OF 2ND KIND */
809 
810  if (*icon != 1) {
811  goto L90;
812  }
813  type2g_(&s2sum, &i__, &inpt_1.nstep1);
814  gnsum += s2sum;
815  sum += s2sum;
816 L90:
817  sum *= inpt_1.estep;
818  gnsum *= inpt_1.estep;
819  gnssum += gnsum;
820  ssum += sum;
821  sum2 = mag_1.qfemag[i__ - 1] * (f0c_1.f[i__ - 1] + f2c_1.f2[i__ - 1] *
822  (float).4 - vel * f1c_1.f1[i__ - 1] / (mag_1.eovm *
823  mix2_1.eroot[i__ - 1]));
824  sum3 = -g0c_1.g[i__ - 1] * (mix1_1.qelm[i__ - 1] + mag_1.qfemag[i__ -
825  1] * inpt_1.alpnaz);
826  if (*lmax == 1) {
827  goto L901;
828  }
829  sum4 = (g2c_1.g2[i__ - 1] * inpt_1.alpnaz * mag_1.qfemag[i__ - 1] +
830  mag_1.qeeef[i__ - 1] * (g2c_1.dg2[i__ - 1] + g2c_1.g2[i__ - 1]
831  * (float)1.5 / mix2_1.e[i__ - 1])) * (float).4;
832 L901:
833  g0c_1.dg0[i__ - 1] = (ssum + sum2 + sum3 - sum4) / (mag_1.qeeef[i__ -
834  1] + mix1_1.qelm[i__ - 1] * inpt_1.akt);
835 
836  ext = g0c_1.g[i__ - 1] * mix1_1.qelm[i__ - 1] + mix1_1.qelm[i__ - 1] *
837  g0c_1.dg0[i__ - 1] * inpt_1.akt;
838  g1c_1.g1[i__ - 1] = (ext + gssum - gnssum) * (float)3. / (mag_1.emag *
839  mix2_1.e[i__ - 1]);
840  g1sum = inpt_1.alpnew * (mix2_1.e[i__ - 1] * (float)2. * g1c_1.g1[i__
841  - 1] - mix2_1.e[i__] * g1c_1.g1[i__]) * inpt_1.estep / (float)
842  3.;
843  gssum += g1sum;
844  if (i__ >= inpt_1.nstep1 - 2) {
845  goto L623;
846  }
847  g1c_1.dg1[i__ - 1] = (g1c_1.g1[i__ + 1] - g1c_1.g1[i__ - 1]) /
848  inpt_1.estep - g1c_1.dg1[i__ + 1];
849  g1c_1.dg1[i__] = (g1c_1.g1[i__ + 1] - g1c_1.g1[i__ - 1]) / (
850  inpt_1.estep * (float)2.);
851  goto L624;
852 L623:
853  g1c_1.dg1[i__ - 1] = (g1c_1.g1[i__] - g1c_1.g1[i__ - 1]) /
854  inpt_1.estep;
855  if (i__ == inpt_1.nstep1) {
856  g1c_1.dg1[i__ - 1] = -g1c_1.g1[i__ - 1] / (inpt_1.estep * (float)
857  2.);
858  }
859 
860 L624:
861  if (*lmax == 1) {
862  goto L903;
863  }
864 
865 /* PN, 07.10.99: make continuation readable: */
866  g2c_1.g2[i__ - 1] = f1c_1.f1[i__ - 1] * (float)2. / (float)3. - vel *
867  f2c_1.f2[i__ - 1] / (mix2_1.eroot[i__ - 1] * mag_1.eovm) -
868  mag_1.ef[i__ - 1] * (float)2. * (g1c_1.dg1[i__ - 1] -
869  g1c_1.g1[i__ - 1] / (mix2_1.e[i__ - 1] * (float)2.)) / (float)
870  3. - inpt_1.alpnaz * (float)2. * g1c_1.g1[i__ - 1] / (float)
871  3. + f3c_1.f3[i__ - 1] * (float)3. / (float)7.;
872  g2c_1.g2[i__ - 1] /= mix2_1.qtot[i__ - 1];
873  g2c_1.dg2[i__ - 1] = g0c_1.dg0[i__ - 1] * (float)-2.5 - g2c_1.g2[i__
874  - 1] * (float)1.5 / mix2_1.e[i__ - 1] - inpt_1.alpnaz * (
875  float)2.5 * (g0c_1.g[i__ - 1] + g2c_1.g2[i__ - 1] * (float).4)
876  / mag_1.ef[i__ - 1] - g1c_1.g1[i__ - 1] * (float)2.5 /
877  mag_1.qe[i__ - 1] + (f0c_1.f[i__ - 1] + f2c_1.f2[i__ - 1] * (
878  float).4) * (float)2.5 / mag_1.ef[i__ - 1] - vel * (float)2.5
879  * f1c_1.f1[i__ - 1] / (mix2_1.eroot[i__ - 1] * mag_1.eovm *
880  mag_1.ef[i__ - 1]);
881  g2c_1.dg2[inpt_1.nstep1] = g2c_1.dg2[inpt_1.nstep1 - 1];
882  g2c_1.g2[i__ - 2] = g2c_1.g2[i__ - 1] - inpt_1.estep * (g2c_1.dg2[i__
883  - 1] * (float)1.5 - g2c_1.dg2[i__] * (float).5);
884  g2c_1.dg2[i__ - 2] = g2c_1.dg2[i__ - 1] * (float)2. - g2c_1.dg2[i__];
885 L903:
886  g0c_1.g[i__ - 2] = g0c_1.g[i__ - 1] - inpt_1.estep * (g0c_1.dg0[i__ -
887  1] * (float)1.5 - g0c_1.dg0[i__] * (float).5);
888  g0c_1.dg[i__ - 1] = g0c_1.dg0[i__ - 1] / g0c_1.g[i__ - 1];
889 /* L1000: */
890  }
891  i__1 = inpt_1.nstep1;
892  for (j = 1; j <= i__1; ++j) {
893 /* L400: */
894  sint_1.simf[j - 1] = g0c_1.g[j - 1] * mix2_1.eroot[j - 1];
895  }
896  simp_(&anew);
897  *eg0sum = anew;
898  return 0;
899 
900 } /* g0calc_ */
901 
902 /* Argon gas */
903 int StFtpcMagboltz1::gas1_(double *q, double *qin, double *q2ro, int *nin,
904 int *n2ro, double *e, double *ein, double *e2ro, char *name__,
905  double *virial, int *monte, int name_len)
906 {
907  /* Initialized data */
908 
909  static double xen[34] = { 1.,1.2,1.5,1.7,2.,2.5,3.,4.,4.9,5.,6.,6.67,
910  7.,8.,8.71,9.,10.,11.,12.,13.,13.6,14.,15.,16.,16.5,18.,20.,30.,
911  30.6,50.,54.4,100.,400.,1e3 };
912  static double yxsec[34] = { 1.39,1.66,2.05,2.33,2.7,3.43,4.15,5.65,
913  7.26,7.46,9.32,10.6,11.3,13.1,14.1,14.4,15.4,15.8,15.8,15.4,15.1,
914  14.8,14.1,13.2,13.,11.4,10.2,6.13,6.01,4.17,3.97,2.71,1.3,1. };
915  static double xeni[54] = { 15.7,16.,16.5,17.,17.5,18.,18.5,19.,19.5,
916  20.,20.5,21.,21.5,22.,22.5,23.,23.5,24.,24.5,25.,25.5,26.,28.,30.,
917  32.,34.,36.,38.,40.,50.,60.,70.,80.,90.,100.,110.,120.,130.,140.,
918  150.,160.,180.,200.,250.,300.,350.,400.,450.,500.,600.,700.,800.,
919  900.,1e3 };
920  static double yxeni[54] = { -.2,.306,.825,1.126,1.326,1.468,1.577,
921  1.663,1.737,1.797,1.853,1.896,1.933,1.97,1.997,2.024,2.048,2.071,
922  2.094,2.115,2.132,2.148,2.204,2.256,2.293,2.325,2.351,2.368,2.379,
923  2.404,2.424,2.443,2.454,2.456,2.455,2.452,2.448,2.441,2.436,2.429,
924  2.419,2.401,2.379,2.337,2.296,2.258,2.225,2.19,2.164,2.115,2.065,
925  2.027,1.994,1.961 };
926  static double xin[15] = { 11.55,13.,13.2,13.4,14.,16.,20.,30.,40.,50.,
927  80.,100.,200.,500.,1e3 };
928  static double yxsin[15] = { 0.,.057,.075,.072,.096,.17,.18,.21,.24,
929  .28,.22,.17,.11,.06,.04 };
930  static double yxpin[15] = { 0.,0.,.01,.03,.06,.17,.35,.45,.44,.41,.3,
931  .27,.18,.09,.05 };
932  static double yxdin[15] = { 0.,0.,0.,0.,0.,.05,.11,.22,.26,.28,.35,
933  .35,.26,.13,.08 };
934 
935  /* System generated locals */
936  int i__1, i__2;
937  double d__1;
938 
939  /* Local variables */
940  static double apol;
941  static int lmax;
942  static double sumi, a, b;
943  static int i__, j, ndata;
944  static double e1, aa, dd, ff, ak, en;
945  static int nidata, nxdata;
946  static double an0, an1, an2, api, sum;
947 
948  /* Parameter adjustments */
949  --e2ro;
950  --ein;
951  --e;
952  q2ro -= 9;
953  qin -= 25;
954  q -= 7;
955 
956  /* Function Body */
957  sprintf(name__, " ARGON 1988 ");
958 /* ---------------------------------------------------------------- */
959 /* MULTI-TERM CROSS-SECTION. */
960 /* FOR PURE ARGON: */
961 /* ACCURACY OF DERIVED VELOCITY AND DIFFUSION COEFFICIENTS 0.5% BELOW */
962 /* 3000VOLTS . BELOW 20000VOLTS ACCURACY 1.0%. IONISATION COEFFICIENT */
963 /* AND DRIFT VELOCITY ACCURACY BETTER THAN 10% BELOW 500,000VOLTS. */
964 /* ----------------------------------------------------------------- */
965 
966 /* PARAMETERS OF PHASE SHIFT ANALYSIS. */
967 
968  apol = (float)11.08;
969  lmax = 100;
970  aa = (float)-1.488;
971  dd = (float)65.4;
972  ff = (float)-84.3;
973  e1 = (float).883;
974  api = (float)3.141592654;
975 
976  *nin = 3;
977  *n2ro = 0;
978  e[1] = (float)0.;
979  e[2] = cnsts_1.emass * (float)2. / (cnsts_1.amu * (float)39.948);
980  e[3] = (float)15.7;
981  e[4] = (float)0.;
982  e[5] = (float)0.;
983  e[6] = (float)0.;
984  ein[1] = (float)11.55;
985  ein[2] = (float)13.;
986  ein[3] = (float)14.;
987  en = -inpt_1.estep;
988  if (*monte == 1) {
989  en = -inpt_1.estep / (float)2.;
990  }
991  i__1 = inpt_1.nstep1 + 1;
992  for (i__ = 1; i__ <= i__1; ++i__) {
993  en += inpt_1.estep;
994  if (en > (float)1.) {
995  goto L100;
996  }
997  if (en == (float)0.) {
998  q[i__ * 6 + 2] = (float)7.79e-16;
999  }
1000  if (en == (float)0.) {
1001  goto L200;
1002  }
1003  ak = ::sqrt(en / inpt_1.ary);
1004  an0 = -aa * ak * (apol * (float)4. / (float)3. * ak * ak * ::log(ak) + (
1005  float)1.) - api * apol / (float)3. * ak * ak + dd * ak * ak *
1006  ak + ff * ak * ak * ak * ak;
1007  an1 = api / (float)15. * apol * ak * ak * ((float)1. - ::sqrt(en / e1));
1008  an2 = api * apol * ak * ak / (float)105.;
1009  an0 = atan(an0);
1010  an1 = atan(an1);
1011  an2 = atan(an2);
1012 /* Computing 2nd power */
1013  d__1 = sin(an0 - an1);
1014  sum = d__1 * d__1;
1015 /* Computing 2nd power */
1016  d__1 = sin(an1 - an2);
1017  sum += d__1 * d__1 * (float)2.;
1018  i__2 = lmax - 1;
1019  for (j = 2; j <= i__2; ++j) {
1020  sumi = (float)6. / ((j * (float)2. + (float)5.) * (j * (float)2.
1021  + (float)3.) * (j * (float)2. + (float)1.) * (j * (float)
1022  2. - (float)1.));
1023 /* Computing 2nd power */
1024  d__1 = sin(atan(api * apol * ak * ak * sumi));
1025  sum += (j + (float)1.) * (d__1 * d__1);
1026 /* L10: */
1027  }
1028  q[i__ * 6 + 2] = sum * (float)4. * cnsts_1.pir2 / (ak * ak);
1029  goto L200;
1030 L100:
1031  ndata = 34;
1032  i__2 = ndata;
1033  for (j = 2; j <= i__2; ++j) {
1034  if (en <= xen[j - 1]) {
1035  goto L160;
1036  }
1037 /* L150: */
1038  }
1039  j = ndata;
1040 L160:
1041  a = (yxsec[j - 1] - yxsec[j - 2]) / (xen[j - 1] - xen[j - 2]);
1042  b = (xen[j - 2] * yxsec[j - 1] - xen[j - 1] * yxsec[j - 2]) / (xen[j
1043  - 2] - xen[j - 1]);
1044  q[i__ * 6 + 2] = (a * en + b) * (float)1e-16;
1045 L200:
1046  q[i__ * 6 + 3] = (float)0.;
1047  if (en <= e[3]) {
1048  goto L230;
1049  }
1050  nidata = 54;
1051  i__2 = nidata;
1052  for (j = 2; j <= i__2; ++j) {
1053  if (en <= xeni[j - 1]) {
1054  goto L220;
1055  }
1056 /* L210: */
1057  }
1058  j = nidata;
1059 L220:
1060  a = (yxeni[j - 1] - yxeni[j - 2]) / (xeni[j - 1] - xeni[j - 2]);
1061  b = (xeni[j - 2] * yxeni[j - 1] - xeni[j - 1] * yxeni[j - 2]) / (xeni[
1062  j - 2] - xeni[j - 1]);
1063  d__1 = a * en + b;
1064  q[i__ * 6 + 3] = ::pow(c_b12, d__1) * (float)1e-18;
1065 L230:
1066  q[i__ * 6 + 4] = (float)0.;
1067  q[i__ * 6 + 5] = (float)0.;
1068  q[i__ * 6 + 6] = (float)0.;
1069 
1070  qin[i__ * 24 + 1] = (float)0.;
1071  qin[i__ * 24 + 2] = (float)0.;
1072  qin[i__ * 24 + 3] = (float)0.;
1073  if (en <= ein[1]) {
1074  goto L400;
1075  }
1076  nxdata = 15;
1077  i__2 = nxdata;
1078  for (j = 2; j <= i__2; ++j) {
1079  if (en <= xin[j - 1]) {
1080  goto L320;
1081  }
1082 /* L310: */
1083  }
1084  j = nxdata;
1085 L320:
1086  a = (yxsin[j - 1] - yxsin[j - 2]) / (xin[j - 1] - xin[j - 2]);
1087  b = (xin[j - 2] * yxsin[j - 1] - xin[j - 1] * yxsin[j - 2]) / (xin[j
1088  - 2] - xin[j - 1]);
1089  qin[i__ * 24 + 1] = (a * en + b) * (float)1e-16;
1090  if (en <= ein[2]) {
1091  goto L400;
1092  }
1093  a = (yxpin[j - 1] - yxpin[j - 2]) / (xin[j - 1] - xin[j - 2]);
1094  b = (xin[j - 2] * yxpin[j - 1] - xin[j - 1] * yxpin[j - 2]) / (xin[j
1095  - 2] - xin[j - 1]);
1096  qin[i__ * 24 + 2] = (a * en + b) * (float)1e-16;
1097  if (en <= ein[3]) {
1098  goto L400;
1099  }
1100  a = (yxdin[j - 1] - yxdin[j - 2]) / (xin[j - 1] - xin[j - 2]);
1101  b = (xin[j - 2] * yxdin[j - 1] - xin[j - 1] * yxdin[j - 2]) / (xin[j
1102  - 2] - xin[j - 1]);
1103  qin[i__ * 24 + 3] = (a * en + b) * (float)1e-16;
1104 L400:
1105  q[i__ * 6 + 1] = q[i__ * 6 + 2] + q[i__ * 6 + 3] + qin[i__ * 24 + 1]
1106  + qin[i__ * 24 + 2] + qin[i__ * 24 + 3];
1107 /* L900: */
1108  }
1109 /* SAVE COMPUTE TIME */
1110  if (inpt_1.efinal <= ein[3]) {
1111  *nin = 2;
1112  }
1113  if (inpt_1.efinal <= ein[2]) {
1114  *nin = 1;
1115  }
1116  if (inpt_1.efinal <= ein[1]) {
1117  *nin = 0;
1118  }
1119 
1120  return 0;
1121 } /* gas1_ */
1122 
1123 
1124 /* CO2 */
1125 int StFtpcMagboltz1::gas2_(double *q, double *qin, double *q2ro, int *nin,
1126 int *n2ro, double *e, double *ein, double *e2ro, char *name__,
1127  double *virial, int *monte, int name_len)
1128 {
1129  /* Initialized data */
1130 
1131  static double xmom[54] = { 0.,.001,.002,.003,.005,.007,.0085,.01,.015,
1132  .02,.03,.04,.05,.07,.1,.12,.15,.17,.2,.25,.3,.35,.4,.5,.7,1.,1.2,
1133  1.3,1.5,1.7,1.9,2.1,2.2,2.5,2.8,3.,3.3,3.6,4.,4.5,5.,6.,7.,8.,10.,
1134  12.,15.,17.,20.,25.,30.,50.,75.,100. };
1135  static double yvib4[24] = { 0.,0.,.95,1.7,1.85,2.,2.15,2.2,2.15,2.,
1136  1.85,1.55,1.23,1.,.76,.64,.49,.44,.41,.48,.26,.135,.1,0. };
1137  static double xvib5[12] = { 0.,.339,1.5,1.95,2.5,3.,3.56,4.1,4.5,5.06,
1138  6.,100. };
1139  static double yvib5[12] = { 0.,0.,0.,.07,.2,.41,.66,.34,.155,0.,0.,0.
1140  };
1141  static double xvib6[12] = { 0.,.422,1.5,1.95,2.5,3.,3.56,4.1,4.5,5.06,
1142  6.,100. };
1143  static double yvib6[12] = { 0.,0.,0.,0.,0.,.105,.225,.1,0.,0.,0.,0. };
1144  static double xvib7[12] = { 0.,.505,1.5,1.95,2.5,3.,3.56,4.1,4.5,5.06,
1145  6.,100. };
1146  static double yvib7[12] = { 0.,0.,0.,0.,0.,.156,.33,.156,0.,0.,0.,0. }
1147  ;
1148  static double xexc1[7] = { 0.,2.5,3.,3.6,4.1,4.5,100. };
1149  static double yexc1[7] = { 0.,0.,.18,.25,.18,0.,0. };
1150  static double xatt[12] = { 0.,3.85,4.3,4.5,5.1,6.6,7.2,8.2,8.4,8.9,
1151  9.7,100. };
1152  static double ymom[54] = { 600.,540.,380.,307.,237.,200.,182.,170.,
1153  138.,120.,97.,85.,76.,63.,50.,44.,39.,34.,29.,24.,18.,15.,13.,10.,
1154  7.1,5.2,4.8,4.7,4.65,4.65,4.85,5.05,5.2,6.4,7.6,9.,11.5,14.,15.2,
1155  14.8,12.7,10.,10.,10.8,12.1,13.1,14.4,15.,15.8,16.,15.8,12.6,9.5,
1156  8. };
1157  static double yatt[12] = { 0.,0.,.0014,.0014,0.,0.,7e-4,.0045,.0042,
1158  .001,0.,0. };
1159  static double xexc2[6] = { 0.,7.,8.,8.5,11.,100. };
1160  static double yexc2[6] = { 0.,0.,.6,.6,0.,0. };
1161  static double xexc3[10] = { 0.,10.5,12.,12.7,13.5,15.,17.,20.,40.,
1162  100. };
1163  static double yexc3[10] = { 0.,0.,.69,.73,.78,.88,1.04,1.24,3.6,6.3 };
1164  static double xion[12] = { 0.,13.3,14.5,15.,16.,18.,20.,30.,40.,50.,
1165  70.,100. };
1166  static double yion[12] = { 0.,0.,.06,.104,.188,.359,.532,1.63,2.28,
1167  2.79,3.43,3.79 };
1168  static double xvib1[38] = { 0.,.083,.0844,.0862,.0932,.1035,.1208,
1169  .1382,.1726,.207,.275,.345,.5,.7,.9,1.1,1.4,1.6,1.8,2.3,2.6,3.,
1170  3.2,3.4,3.6,3.8,4.,4.2,4.6,5.1,5.5,6.,7.,8.,10.,20.,50.,100. };
1171  static double yvib1[38] = { 0.,0.,.85,1.16,1.85,2.3,2.6,2.68,2.54,2.2,
1172  1.72,1.43,1.08,.8,.66,.57,.45,.42,.44,.7,.93,1.34,1.58,1.75,1.8,
1173  1.79,1.7,1.52,1.05,.57,.51,.5,.48,.45,.2,0.,0.,0. };
1174  static double xvib2[28] = { 0.,.167,.172,.18,.2,.25,.5,1.,1.5,1.9,2.,
1175  2.25,2.5,3.,3.2,3.4,3.55,3.7,3.9,4.1,4.5,4.9,5.2,6.,8.,10.,20.,
1176  100. };
1177  static double yvib2[28] = { 0.,0.,.3,.33,.35,.325,.117,.05,.04,.06,
1178  .08,.2,.4,1.28,1.57,1.77,1.78,1.75,1.6,1.28,.88,.39,.33,.27,.25,
1179  .21,0.,0. };
1180  static double xvib3[12] = { 0.,.252,1.5,1.95,2.5,3.,3.56,4.1,4.5,5.06,
1181  6.,100. };
1182  static double yvib3[12] = { 0.,0.,0.,0.,0.,.32,.54,.34,.16,.044,0.,0.
1183  };
1184  static double xvib4[24] = { 0.,.291,.3,.31,.32,.33,.35,.38,.4,.45,.5,
1185  .6,.8,1.,1.5,2.,3.,4.5,6.,8.,10.,25.,30.,100. };
1186 
1187  /* System generated locals */
1188  int i__1, i__2;
1189 
1190  /* Local variables */
1191  static int nion, nmom, natt, nexc1, nvib1, nvib2, nvib3, nvib4, nvib5,
1192  nvib6, nvib7, nexc2, nexc3, i__, j;
1193  static double a, b, en;
1194 
1195  /* Parameter adjustments */
1196  --e2ro;
1197  --ein;
1198  --e;
1199  q2ro -= 9;
1200  qin -= 25;
1201  q -= 7;
1202 
1203  /* Function Body */
1204 
1205  sprintf(name__, "C02 TEST # 2 ");
1206  *nin = 10;
1207  *n2ro = 0;
1208  nmom = 54;
1209  nvib1 = 38;
1210  nvib2 = 28;
1211  nvib3 = 12;
1212  nvib4 = 24;
1213  nvib5 = 12;
1214  nvib6 = 12;
1215  nvib7 = 12;
1216  nexc1 = 7;
1217  natt = 12;
1218  nexc2 = 6;
1219  nexc3 = 10;
1220  nion = 12;
1221  e[1] = (float)0.;
1222  e[2] = cnsts_1.emass * (float)2. / (cnsts_1.amu * (float)44.0098);
1223  e[3] = (float)13.3;
1224  e[4] = (float)3.85;
1225  e[5] = (float)0.;
1226  e[6] = (float)0.;
1227  ein[1] = (float).083;
1228  ein[2] = (float).167;
1229  ein[3] = (float).252;
1230  ein[4] = (float).291;
1231  ein[5] = (float).339;
1232  ein[6] = (float).422;
1233  ein[7] = (float).505;
1234  ein[8] = (float)2.5;
1235  ein[9] = (float)7.;
1236  ein[10] = (float)10.5;
1237  en = -inpt_1.estep;
1238  if (*monte == 1) {
1239  en = -inpt_1.estep / (float)2.;
1240  }
1241  i__1 = inpt_1.nstep1 + 1;
1242  for (i__ = 1; i__ <= i__1; ++i__) {
1243  en += inpt_1.estep;
1244 
1245  i__2 = nmom;
1246  for (j = 2; j <= i__2; ++j) {
1247  if (en <= xmom[j - 1]) {
1248  goto L150;
1249  }
1250 /* L100: */
1251  }
1252  j = nmom;
1253 L150:
1254  a = (ymom[j - 1] - ymom[j - 2]) / (xmom[j - 1] - xmom[j - 2]);
1255  b = (xmom[j - 2] * ymom[j - 1] - xmom[j - 1] * ymom[j - 2]) / (xmom[j
1256  - 2] - xmom[j - 1]);
1257  q[i__ * 6 + 2] = (a * en + b) * (float)1e-16;
1258 
1259  qin[i__ * 24 + 1] = (float)0.;
1260  if (en <= ein[1]) {
1261  goto L260;
1262  }
1263  i__2 = nvib1;
1264  for (j = 2; j <= i__2; ++j) {
1265  if (en <= xvib1[j - 1]) {
1266  goto L250;
1267  }
1268 /* L200: */
1269  }
1270  j = nvib1;
1271 L250:
1272  a = (yvib1[j - 1] - yvib1[j - 2]) / (xvib1[j - 1] - xvib1[j - 2]);
1273  b = (xvib1[j - 2] * yvib1[j - 1] - xvib1[j - 1] * yvib1[j - 2]) / (
1274  xvib1[j - 2] - xvib1[j - 1]);
1275  qin[i__ * 24 + 1] = (a * en + b) * (float)1e-16;
1276 
1277 L260:
1278  qin[i__ * 24 + 2] = (float)0.;
1279  if (en <= ein[2]) {
1280  goto L360;
1281  }
1282  i__2 = nvib2;
1283  for (j = 2; j <= i__2; ++j) {
1284  if (en <= xvib2[j - 1]) {
1285  goto L350;
1286  }
1287 /* L300: */
1288  }
1289  j = nvib2;
1290 L350:
1291  a = (yvib2[j - 1] - yvib2[j - 2]) / (xvib2[j - 1] - xvib2[j - 2]);
1292  b = (xvib2[j - 2] * yvib2[j - 1] - xvib2[j - 1] * yvib2[j - 2]) / (
1293  xvib2[j - 2] - xvib2[j - 1]);
1294  qin[i__ * 24 + 2] = (a * en + b) * (float)1e-16;
1295 
1296 L360:
1297  qin[i__ * 24 + 3] = (float)0.;
1298  if (en <= ein[3]) {
1299  goto L460;
1300  }
1301  i__2 = nvib3;
1302  for (j = 2; j <= i__2; ++j) {
1303  if (en <= xvib3[j - 1]) {
1304  goto L450;
1305  }
1306 /* L400: */
1307  }
1308  j = nvib3;
1309 L450:
1310  a = (yvib3[j - 1] - yvib3[j - 2]) / (xvib3[j - 1] - xvib3[j - 2]);
1311  b = (xvib3[j - 2] * yvib3[j - 1] - xvib3[j - 1] * yvib3[j - 2]) / (
1312  xvib3[j - 2] - xvib3[j - 1]);
1313  qin[i__ * 24 + 3] = (a * en + b) * (float)1e-16;
1314 
1315 L460:
1316  qin[i__ * 24 + 4] = (float)0.;
1317  if (en <= ein[4]) {
1318  goto L560;
1319  }
1320  i__2 = nvib4;
1321  for (j = 2; j <= i__2; ++j) {
1322  if (en <= xvib4[j - 1]) {
1323  goto L550;
1324  }
1325 /* L500: */
1326  }
1327  j = nvib4;
1328 L550:
1329  a = (yvib4[j - 1] - yvib4[j - 2]) / (xvib4[j - 1] - xvib4[j - 2]);
1330  b = (xvib4[j - 2] * yvib4[j - 1] - xvib4[j - 1] * yvib4[j - 2]) / (
1331  xvib4[j - 2] - xvib4[j - 1]);
1332  qin[i__ * 24 + 4] = (a * en + b) * (float)1e-16;
1333 
1334 L560:
1335  qin[i__ * 24 + 5] = (float)0.;
1336  if (en <= ein[5]) {
1337  goto L660;
1338  }
1339  i__2 = nvib5;
1340  for (j = 2; j <= i__2; ++j) {
1341  if (en <= xvib5[j - 1]) {
1342  goto L650;
1343  }
1344 /* L600: */
1345  }
1346  j = nvib5;
1347 L650:
1348  a = (yvib5[j - 1] - yvib5[j - 2]) / (xvib5[j - 1] - xvib5[j - 2]);
1349  b = (xvib5[j - 2] * yvib5[j - 1] - xvib5[j - 1] * yvib5[j - 2]) / (
1350  xvib5[j - 2] - xvib5[j - 1]);
1351  qin[i__ * 24 + 5] = (a * en + b) * (float)1e-16;
1352 
1353 L660:
1354  qin[i__ * 24 + 6] = (float)0.;
1355  if (en <= ein[6]) {
1356  goto L760;
1357  }
1358  i__2 = nvib6;
1359  for (j = 2; j <= i__2; ++j) {
1360  if (en <= xvib6[j - 1]) {
1361  goto L750;
1362  }
1363 /* L700: */
1364  }
1365  j = nvib6;
1366 L750:
1367  a = (yvib6[j - 1] - yvib6[j - 2]) / (xvib6[j - 1] - xvib6[j - 2]);
1368  b = (xvib6[j - 2] * yvib6[j - 1] - xvib6[j - 1] * yvib6[j - 2]) / (
1369  xvib6[j - 2] - xvib6[j - 1]);
1370  qin[i__ * 24 + 6] = (a * en + b) * (float)1e-16;
1371 
1372 L760:
1373  qin[i__ * 24 + 7] = (float)0.;
1374  if (en <= ein[7]) {
1375  goto L860;
1376  }
1377  i__2 = nvib7;
1378  for (j = 2; j <= i__2; ++j) {
1379  if (en <= xvib7[j - 1]) {
1380  goto L850;
1381  }
1382 /* L800: */
1383  }
1384  j = nvib7;
1385 L850:
1386  a = (yvib7[j - 1] - yvib7[j - 2]) / (xvib7[j - 1] - xvib7[j - 2]);
1387  b = (xvib7[j - 2] * yvib7[j - 1] - xvib7[j - 1] * yvib7[j - 2]) / (
1388  xvib7[j - 2] - xvib7[j - 1]);
1389  qin[i__ * 24 + 7] = (a * en + b) * (float)1e-16;
1390 
1391 L860:
1392  qin[i__ * 24 + 8] = (float)0.;
1393  if (en <= ein[8]) {
1394  goto L960;
1395  }
1396  i__2 = nexc1;
1397  for (j = 2; j <= i__2; ++j) {
1398  if (en <= xexc1[j - 1]) {
1399  goto L950;
1400  }
1401 /* L900: */
1402  }
1403  j = nexc1;
1404 L950:
1405  a = (yexc1[j - 1] - yexc1[j - 2]) / (xexc1[j - 1] - xexc1[j - 2]);
1406  b = (xexc1[j - 2] * yexc1[j - 1] - xexc1[j - 1] * yexc1[j - 2]) / (
1407  xexc1[j - 2] - xexc1[j - 1]);
1408  qin[i__ * 24 + 8] = (a * en + b) * (float)1e-16;
1409 
1410 L960:
1411  q[i__ * 6 + 4] = (float)0.;
1412  if (en <= e[4]) {
1413  goto L1060;
1414  }
1415  i__2 = natt;
1416  for (j = 2; j <= i__2; ++j) {
1417  if (en <= xatt[j - 1]) {
1418  goto L1050;
1419  }
1420 /* L1000: */
1421  }
1422  j = natt;
1423 L1050:
1424  a = (yatt[j - 1] - yatt[j - 2]) / (xatt[j - 1] - xatt[j - 2]);
1425  b = (xatt[j - 2] * yatt[j - 1] - xatt[j - 1] * yatt[j - 2]) / (xatt[j
1426  - 2] - xatt[j - 1]);
1427  q[i__ * 6 + 4] = (a * en + b) * (float)1e-16;
1428 
1429 L1060:
1430  qin[i__ * 24 + 9] = (float)0.;
1431  if (en <= ein[9]) {
1432  goto L1160;
1433  }
1434  i__2 = nexc2;
1435  for (j = 2; j <= i__2; ++j) {
1436  if (en <= xexc2[j - 1]) {
1437  goto L1150;
1438  }
1439 /* L1100: */
1440  }
1441  j = nexc2;
1442 L1150:
1443  a = (yexc2[j - 1] - yexc2[j - 2]) / (xexc2[j - 1] - xexc2[j - 2]);
1444  b = (xexc2[j - 2] * yexc2[j - 1] - xexc2[j - 1] * yexc2[j - 2]) / (
1445  xexc2[j - 2] - xexc2[j - 1]);
1446  qin[i__ * 24 + 9] = (a * en + b) * (float)1e-16;
1447 
1448 L1160:
1449  qin[i__ * 24 + 10] = (float)0.;
1450  if (en <= ein[10]) {
1451  goto L1260;
1452  }
1453  i__2 = nexc3;
1454  for (j = 2; j <= i__2; ++j) {
1455  if (en <= xexc3[j - 1]) {
1456  goto L1250;
1457  }
1458 /* L1200: */
1459  }
1460  j = nexc3;
1461 L1250:
1462  a = (yexc3[j - 1] - yexc3[j - 2]) / (xexc3[j - 1] - xexc3[j - 2]);
1463  b = (xexc3[j - 2] * yexc3[j - 1] - xexc3[j - 1] * yexc3[j - 2]) / (
1464  xexc3[j - 2] - xexc3[j - 1]);
1465  qin[i__ * 24 + 10] = (a * en + b) * (float)1e-16;
1466 
1467 L1260:
1468  q[i__ * 6 + 3] = (float)0.;
1469  if (en <= e[3]) {
1470  goto L1360;
1471  }
1472  i__2 = nion;
1473  for (j = 2; j <= i__2; ++j) {
1474  if (en <= xion[j - 1]) {
1475  goto L1350;
1476  }
1477 /* L1300: */
1478  }
1479  j = nion;
1480 L1350:
1481  a = (yion[j - 1] - yion[j - 2]) / (xion[j - 1] - xion[j - 2]);
1482  b = (xion[j - 2] * yion[j - 1] - xion[j - 1] * yion[j - 2]) / (xion[j
1483  - 2] - xion[j - 1]);
1484  q[i__ * 6 + 3] = (a * en + b) * (float)1e-16;
1485 
1486 L1360:
1487  q[i__ * 6 + 1] = q[i__ * 6 + 2] + q[i__ * 6 + 3] + q[i__ * 6 + 4];
1488  q[i__ * 6 + 2] = q[i__ * 6 + 2] - qin[i__ * 24 + 1] - qin[i__ * 24 +
1489  2] - qin[i__ * 24 + 3] - qin[i__ * 24 + 4] - qin[i__ * 24 + 5]
1490  - qin[i__ * 24 + 6] - qin[i__ * 24 + 7] - qin[i__ * 24 + 8]
1491  - qin[i__ * 24 + 9] - qin[i__ * 24 + 10];
1492  q[i__ * 6 + 5] = (float)0.;
1493  q[i__ * 6 + 6] = (float)0.;
1494 /* L9000: */
1495  }
1496 
1497 /* SAVE ON COMPUTING TIME */
1498 
1499  if (inpt_1.efinal < ein[10]) {
1500  *nin = 9;
1501  }
1502  if (inpt_1.efinal < ein[9]) {
1503  *nin = 8;
1504  }
1505  if (inpt_1.efinal < ein[8]) {
1506  *nin = 7;
1507  }
1508  if (inpt_1.efinal < ein[7]) {
1509  *nin = 6;
1510  }
1511  if (inpt_1.efinal < ein[6]) {
1512  *nin = 5;
1513  }
1514  if (inpt_1.efinal < ein[5]) {
1515  *nin = 4;
1516  }
1517  if (inpt_1.efinal < ein[4]) {
1518  *nin = 3;
1519  }
1520  if (inpt_1.efinal < ein[3]) {
1521  *nin = 2;
1522  }
1523  if (inpt_1.efinal < ein[2]) {
1524  *nin = 1;
1525  }
1526  if (inpt_1.efinal < ein[1]) {
1527  *nin = 0;
1528  }
1529  return 0;
1530 } /* gas2_ */
1531 
1532 /* Ne */
1533 int StFtpcMagboltz1::gas3_(double *q, double *qin, double *q2ro, int *nin,
1534 int *n2ro, double *e, double *ein, double *e2ro, char *name__,
1535  double *virial, int *monte, int name_len)
1536 {
1537  /* Initialized data */
1538 
1539  static double xen[37] = { 1.,1.2,1.5,1.8,2.,2.5,3.,4.,5.,6.,7.,8.,
1540  8.71,9.,10.,11.,13.6,15.,16.5,19.6,20.,30.,30.6,50.,54.4,65.,77.,
1541  100.,130.,150.,170.,200.,300.,400.,600.,800.,1e3 };
1542  static double yxsec[37] = { 1.62,1.69,1.75,1.79,1.82,1.86,1.91,1.98,
1543  2.07,2.14,2.21,2.29,2.35,2.37,2.44,2.51,2.66,2.71,2.76,2.83,2.84,
1544  2.83,2.82,2.45,2.36,2.18,1.97,1.56,1.22,1.05,.91,.76,.52,.42,.33,
1545  .27,.25 };
1546  static double xion[35] = { 21.56,22.,22.5,23.,23.5,24.,24.5,25.,25.5,
1547  26.,27.,28.,29.,30.,32.,34.,36.,40.,45.,50.,55.,60.,70.,80.,90.,
1548  100.,120.,140.,175.,200.,300.,400.,600.,800.,1e3 };
1549  static double yion[35] = { 0.,.0033,.0089,.0146,.02,.026,.032,.038,
1550  .044,.05,.063,.076,.089,.102,.128,.154,.179,.228,.282,.338,.391,
1551  .435,.514,.58,.63,.67,.72,.76,.78,.78,.72,.63,.53,.44,.39 };
1552  static double xexc[35] = { 16.615,16.78,16.97,17.3,18.4,18.7,18.8,
1553  19.8,20.,21.,22.,24.,26.,28.,30.,35.,40.,50.,60.,70.,80.,90.,100.,
1554  120.,150.,200.,250.,300.,400.,500.,600.,700.,800.,900.,1e3 };
1555  static double yexc[35] = { 0.,.0034,.0185,.012,.0181,.0349,.028,.05,
1556  .0523,.0732,.0923,.123,.143,.162,.176,.195,.205,.207,.203,.195,
1557  .187,.179,.171,.157,.138,.117,.102,.091,.075,.064,.057,.051,.047,
1558  .043,.04 };
1559 
1560  /* System generated locals */
1561  int i__1, i__2;
1562  double d__1;
1563 
1564  /* Local variables */
1565  static double apol;
1566  static int nexc, lmax, nion;
1567  static double sumi, a, b;
1568  static int i__, j, ndata;
1569  static double a1, b1, a2, aa, dd, ff, ak, en, an0, an1, an2, api, sum;
1570 
1571  /* Parameter adjustments */
1572  --e2ro;
1573  --ein;
1574  --e;
1575  q2ro -= 9;
1576  qin -= 25;
1577  q -= 7;
1578 
1579  /* Function Body */
1580 
1581  sprintf(name__, " NEON 92 ");
1582  *nin = 1;
1583  *n2ro = 0;
1584  ndata = 37;
1585  nion = 35;
1586  nexc = 35;
1587  e[1] = (float)0.;
1588  e[2] = cnsts_1.emass * (float)2. / (cnsts_1.amu * (float)20.179);
1589  e[3] = (float)21.56;
1590  e[4] = (float)0.;
1591  e[5] = (float)0.;
1592  e[6] = (float)0.;
1593  ein[1] = (float)16.615;
1594  apol = (float)2.672;
1595  lmax = 100;
1596  aa = (float).2135;
1597  dd = (float)3.86;
1598  ff = (float)-2.656;
1599  a1 = (float)1.846;
1600  b1 = (float)3.29;
1601  a2 = (float)-.037;
1602  api = (float)3.141592654;
1603  en = -inpt_1.estep;
1604  if (*monte == 1) {
1605  en = -inpt_1.estep / (float)2.;
1606  }
1607  i__1 = inpt_1.nstep1 + 1;
1608  for (i__ = 1; i__ <= i__1; ++i__) {
1609  en += inpt_1.estep;
1610  if (en > (float)1.) {
1611  goto L100;
1612  }
1613  if (en == (float)0.) {
1614  q[i__ * 6 + 2] = (float)1.61e-17;
1615  }
1616  if (en == (float)0.) {
1617  goto L200;
1618  }
1619  ak = ::sqrt(en / inpt_1.ary);
1620  an0 = -aa * ak * (apol * (float)4. / (float)3. * ak * ak * ::log(ak) + (
1621  float)1.) - api * apol / (float)3. * ak * ak + dd * ak * ak *
1622  ak + ff * ak * ak * ak * ak;
1623  an1 = (ak * (float).56 * ak - a1 * ak * ak * ak) / (b1 * ak * ak + (
1624  float)1.);
1625  an2 = ak * (float).08 * ak - a2 * ak * ak * ak * ak * ak;
1626 /* Computing 2nd power */
1627  d__1 = sin(an0 - an1);
1628  sum = d__1 * d__1;
1629 /* Computing 2nd power */
1630  d__1 = sin(an1 - an2);
1631  sum += d__1 * d__1 * (float)2.;
1632  i__2 = lmax - 1;
1633  for (j = 2; j <= i__2; ++j) {
1634  sumi = (float)6. / ((j * (float)2. + (float)5.) * (j * (float)2.
1635  + (float)3.) * (j * (float)2. + (float)1.) * (j * (float)
1636  2. - (float)1.));
1637 /* Computing 2nd power */
1638  d__1 = sin(api * apol * ak * ak * sumi);
1639  sum += (j + (float)1.) * (d__1 * d__1);
1640 /* L10: */
1641  }
1642  q[i__ * 6 + 2] = sum * (float)4. * cnsts_1.pir2 / (ak * ak);
1643  goto L200;
1644 L100:
1645  i__2 = ndata;
1646  for (j = 2; j <= i__2; ++j) {
1647  if (en <= xen[j - 1]) {
1648  goto L160;
1649  }
1650 /* L150: */
1651  }
1652  j = ndata;
1653 L160:
1654  a = (yxsec[j - 1] - yxsec[j - 2]) / (xen[j - 1] - xen[j - 2]);
1655  b = (xen[j - 2] * yxsec[j - 1] - xen[j - 1] * yxsec[j - 2]) / (xen[j
1656  - 2] - xen[j - 1]);
1657  q[i__ * 6 + 2] = (a * en + b) * (float)1e-16;
1658 L200:
1659  q[i__ * 6 + 3] = (float)0.;
1660  if (en <= e[3]) {
1661  goto L230;
1662  }
1663  i__2 = nion;
1664  for (j = 2; j <= i__2; ++j) {
1665  if (en <= xion[j - 1]) {
1666  goto L220;
1667  }
1668 /* L210: */
1669  }
1670  j = nion;
1671 L220:
1672  a = (yion[j - 1] - yion[j - 2]) / (xion[j - 1] - xion[j - 2]);
1673  b = (xion[j - 2] * yion[j - 1] - xion[j - 1] * yion[j - 2]) / (xion[j
1674  - 2] - xion[j - 1]);
1675  q[i__ * 6 + 3] = (a * en + b) * (float)1e-16;
1676 L230:
1677  q[i__ * 6 + 4] = (float)0.;
1678  q[i__ * 6 + 5] = (float)0.;
1679  q[i__ * 6 + 6] = (float)0.;
1680 
1681  qin[i__ * 24 + 1] = (float)0.;
1682  if (en <= ein[1]) {
1683  goto L370;
1684  }
1685  i__2 = nexc;
1686  for (j = 2; j <= i__2; ++j) {
1687  if (en <= xexc[j - 1]) {
1688  goto L360;
1689  }
1690 /* L350: */
1691  }
1692  j = nexc;
1693 L360:
1694  a = (yexc[j - 1] - yexc[j - 2]) / (xexc[j - 1] - xexc[j - 2]);
1695  b = (xexc[j - 2] * yexc[j - 1] - xexc[j - 1] * yexc[j - 2]) / (xexc[j
1696  - 2] - xexc[j - 1]);
1697  qin[i__ * 24 + 1] = (a * en + b) * (float)1e-16;
1698 L370:
1699  q[i__ * 6 + 1] = q[i__ * 6 + 2] + q[i__ * 6 + 3] + qin[i__ * 24 + 1];
1700 /* L900: */
1701  }
1702  if (inpt_1.efinal < ein[1]) {
1703  *nin = 0;
1704  }
1705  return 0;
1706 } /* gas3_ */
1707 /* He */
1708 int StFtpcMagboltz1::gas4_(double *q, double *qin, double *q2ro, int *nin,
1709 int *n2ro, double *e, double *ein, double *e2ro, char *name__,
1710  double *virial, int *monte, int name_len)
1711 {
1712  /* Initialized data */
1713 
1714  static double xen[63] = { 0.,.008,.009,.01,.013,.017,.02,.025,.03,.04,
1715  .05,.06,.07,.08,.09,.1,.12,.15,.18,.2,.25,.3,.4,.5,.6,.7,.8,.9,1.,
1716  1.2,1.5,1.8,2.,2.5,3.,4.,5.,6.,7.,8.,9.,10.,11.,12.,13.6,16.5,18.,
1717  20.,25.,30.,40.,50.,60.,70.,75.,80.,90.,100.,150.,200.,400.,600.,
1718  1e3 };
1719  static double yxsec[63] = { 4.9,5.18,5.19,5.21,5.26,5.31,5.35,5.41,
1720  5.46,5.54,5.62,5.68,5.74,5.79,5.83,5.86,5.94,6.04,6.12,6.16,6.27,
1721  6.35,6.49,6.59,6.66,6.73,6.77,6.82,6.85,6.91,6.96,6.98,6.99,6.96,
1722  6.89,6.62,6.31,6.,5.68,5.35,5.03,4.72,4.44,4.15,3.83,3.25,2.99,
1723  2.58,1.95,1.51,.98,.7,.5,.4,.34,.31,.25,.21,.104,.063,.02,.01,
1724  .0035 };
1725  static double xion[28] = { 24.59,25.,25.5,26.,26.5,27.,28.,29.,30.,
1726  32.,34.,36.,38.,40.,45.,50.,55.,60.,70.,80.,100.,120.,150.,175.,
1727  200.,300.,500.,1e3 };
1728  static double yion[28] = { 0.,.0052,.0113,.0175,.0236,.03,.043,.055,
1729  .067,.092,.114,.135,.155,.172,.21,.243,.271,.29,.321,.344,.366,
1730  .373,.37,.359,.347,.297,.224,.141 };
1731  static double xexc[21] = { 19.82,20.,25.,30.,40.,50.,60.,70.,75.,80.,
1732  90.,100.,150.,200.,300.,400.,500.,600.,700.,800.,1e3 };
1733  static double yexc[21] = { 0.,.03,.1,.163,.187,.191,.205,.194,.192,
1734  .189,.187,.185,.165,.146,.121,.099,.087,.076,.068,.061,.051 };
1735 
1736  /* System generated locals */
1737  int i__1, i__2;
1738 
1739  /* Local variables */
1740  static int nexc, nion;
1741  static double a, b;
1742  static int i__, j, ndata;
1743  static double en;
1744 
1745  /* Parameter adjustments */
1746  --e2ro;
1747  --ein;
1748  --e;
1749  q2ro -= 9;
1750  qin -= 25;
1751  q -= 7;
1752 
1753  /* Function Body */
1754 
1755  sprintf(name__, " HELIUM4 ");
1756 /* -------------------------------------------------------------------- */
1757 /* HELIUM 4 BEST KNOWN GAS USED AS STANDARD ACCURACY BETTER THAN 0.2% */
1758 /* AT ALL FIELDS. */
1759 /* -------------------------------------------------------------------- */
1760  *nin = 1;
1761  *n2ro = 0;
1762  ndata = 63;
1763  nion = 28;
1764  nexc = 21;
1765  e[1] = (float)0.;
1766  e[2] = cnsts_1.emass * (float)2. / (cnsts_1.amu * (float)4.0026);
1767  e[3] = (float)24.59;
1768  e[4] = (float)0.;
1769  e[5] = (float)0.;
1770  e[6] = (float)0.;
1771  ein[1] = (float)19.82;
1772  en = -inpt_1.estep;
1773  if (*monte == 1) {
1774  en = -inpt_1.estep / (float)2.;
1775  }
1776  i__1 = inpt_1.nstep1 + 1;
1777  for (i__ = 1; i__ <= i__1; ++i__) {
1778  en += inpt_1.estep;
1779  i__2 = ndata;
1780  for (j = 2; j <= i__2; ++j) {
1781  if (en <= xen[j - 1]) {
1782  goto L20;
1783  }
1784 /* L10: */
1785  }
1786  j = ndata;
1787 L20:
1788  a = (yxsec[j - 1] - yxsec[j - 2]) / (xen[j - 1] - xen[j - 2]);
1789  b = (xen[j - 2] * yxsec[j - 1] - xen[j - 1] * yxsec[j - 2]) / (xen[j
1790  - 2] - xen[j - 1]);
1791  q[i__ * 6 + 2] = (a * en + b) * (float)1e-16;
1792 
1793  q[i__ * 6 + 3] = (float)0.;
1794  if (en <= e[3]) {
1795  goto L200;
1796  }
1797  i__2 = nion;
1798  for (j = 2; j <= i__2; ++j) {
1799  if (en <= xion[j - 1]) {
1800  goto L120;
1801  }
1802 /* L110: */
1803  }
1804  j = nion;
1805 L120:
1806  a = (yion[j - 1] - yion[j - 2]) / (xion[j - 1] - xion[j - 2]);
1807  b = (xion[j - 2] * yion[j - 1] - xion[j - 1] * yion[j - 2]) / (xion[j
1808  - 2] - xion[j - 1]);
1809  q[i__ * 6 + 3] = (a * en + b) * (float)1e-16;
1810 
1811 L200:
1812  q[i__ * 6 + 4] = (float)0.;
1813  q[i__ * 6 + 5] = (float)0.;
1814  q[i__ * 6 + 6] = (float)0.;
1815 
1816  qin[i__ * 24 + 1] = (float)0.;
1817  if (en <= ein[1]) {
1818  goto L600;
1819  }
1820  i__2 = nexc;
1821  for (j = 2; j <= i__2; ++j) {
1822  if (en <= xexc[j - 1]) {
1823  goto L520;
1824  }
1825 /* L510: */
1826  }
1827  j = nexc;
1828 L520:
1829  a = (yexc[j - 1] - yexc[j - 2]) / (xexc[j - 1] - xexc[j - 2]);
1830  b = (xexc[j - 2] * yexc[j - 1] - xexc[j - 1] * yexc[j - 2]) / (xexc[j
1831  - 2] - xexc[j - 1]);
1832  qin[i__ * 24 + 1] = (a * en + b) * (float)1e-16;
1833  qin[i__ * 24 + 1] = qin[i__ * 24 + 1];
1834 L600:
1835 
1836  q[i__ * 6 + 1] = q[i__ * 6 + 2] + qin[i__ * 24 + 1] + q[i__ * 6 + 3];
1837 /* L900: */
1838  }
1839 /* SAVE COMPUTE TIME */
1840  if (inpt_1.efinal <= ein[1]) {
1841  *nin = 0;
1842  }
1843 
1844  return 0;
1845 } /* gas4_ */
1846 
1847 int StFtpcMagboltz1::h1calc_(int *l, double *dhfnal, double *dxx, double *dhfrst)
1848 {
1849  /* System generated locals */
1850  int i__1;
1851 
1852  /* Local variables */
1853  static double ssum0;
1854  static int i__, j, m;
1855  static double termb, termt, sum, sum0, sum1, sum2, sum3, sum4, sum5;
1856 
1857  if (*l == 2) {
1858  goto L300;
1859  }
1860  i__1 = inpt_1.nstep1;
1861  for (i__ = 1; i__ <= i__1; ++i__) {
1862 /* L1000: */
1863  h1c_1.h1[i__ - 1] = (f0c_1.f[i__ - 1] - f2c_1.f2[i__ - 1] * (float).2)
1864  / mix2_1.qtot[i__ - 1];
1865  }
1866  i__1 = inpt_1.nstep1;
1867  for (j = 2; j <= i__1; ++j) {
1868 /* L200: */
1869  h1c_1.dh1[j - 1] = (h1c_1.h1[j] - h1c_1.h1[j - 2]) / (inpt_1.estep * (
1870  float)2.);
1871  }
1872  h1c_1.dh1[0] = h1c_1.h1[0] / inpt_1.estep;
1873  return 0;
1874 L300:
1875  ssum0 = (float)0.;
1876  for (j = 1; j <= 2002; ++j) {
1877  h1c_1.dh1[j - 1] = (float)0.;
1878 /* L344: */
1879  h1c_1.h1[j - 1] = (float)0.;
1880  }
1881  h1c_1.h1[inpt_1.nstep1 - 1] = (float)0.;
1882  h1c_1.dh1[inpt_1.nstep1 - 1] = -(*dhfnal);
1883  h1c_1.h1[inpt_1.nstep - 1] = -inpt_1.estep * h1c_1.dh1[inpt_1.nstep1 - 1];
1884  h1c_1.dh1[inpt_1.nstep - 1] = h1c_1.dh1[inpt_1.nstep1 - 1];
1885  i__1 = inpt_1.nstep - 1;
1886  for (m = 1; m <= i__1; ++m) {
1887  i__ = inpt_1.nstep - m + 1;
1888  sum0 = h1c_1.h1[i__ - 1];
1889  sum1 = (f0c_1.f[i__ - 1] - f2c_1.f2[i__ - 1] * (float).2) /
1890  mix2_1.qtot[i__ - 1];
1891  sum2 = mag_1.qe[i__ - 1] * (float).3 * f1c_1.f1[i__ - 1] / (mix2_1.e[
1892  i__ - 1] * mix2_1.qtot[i__ - 1]);
1893  sum3 = mag_1.qe[i__ - 1] * (float).3 * mag_1.qe[i__ - 1] * h1c_1.dh1[
1894  i__ - 1] / mix2_1.e[i__ - 1];
1895  sum4 = mag_1.qe[i__ - 1] * (float).15 * mag_1.qe[i__ - 1] * h1c_1.h1[
1896  i__ - 1] / (mix2_1.e[i__ - 1] * mix2_1.e[i__ - 1]);
1897  sum5 = mag_1.qe[i__ - 1] * (float)1.8 * f3c_1.f3[i__ - 1] / (mix2_1.e[
1898  i__ - 1] * mix2_1.qtot[i__ - 1] * (float)14.);
1899  sum = sum1 - sum0 - sum2 + sum3 - sum4 + sum5;
1900  sum *= inpt_1.estep;
1901  ssum0 += sum;
1902  termt = mag_1.qe[i__ - 1] * (float).2 * (f1c_1.f1[i__ - 1] * (float)
1903  1. / mix2_1.qtot[i__ - 1] + mag_1.qe[i__ - 1] * h1c_1.h1[i__
1904  - 1] / (mix2_1.e[i__ - 1] * (float)2.) - f3c_1.f3[i__ - 1] * (
1905  float)3. / (mix2_1.qtot[i__ - 1] * (float)7.));
1906  termb = mag_1.qe[i__ - 1] * (float).2 * mag_1.qe[i__ - 1];
1907  h1c_1.dh1[i__ - 1] = (ssum0 + termt) / termb;
1908  h1c_1.dh1[i__ - 1] -= *dhfnal;
1909  h1c_1.h1[i__ - 2] = h1c_1.h1[i__ - 1] - inpt_1.estep * (h1c_1.dh1[i__
1910  - 1] * (float)1.5 - h1c_1.dh1[i__] * (float).5);
1911  h1c_1.dh1[i__ - 2] = h1c_1.dh1[i__ - 1] * (float)2. - h1c_1.dh1[i__];
1912 /* L500: */
1913  }
1914  *dhfrst = h1c_1.dh1[0];
1915  i__1 = inpt_1.nstep1;
1916  for (j = 1; j <= i__1; ++j) {
1917 /* L125: */
1918  sint_1.simf[j - 1] = mix2_1.e[j - 1] * h1c_1.h1[j - 1];
1919  }
1920  simp_(&sum);
1921  *dxx = sum * mag_1.eovm / (float)3.;
1922  return 0;
1923 } /* h1calc_ */
1924 
1925 int StFtpcMagboltz1::mixer_()
1926 {
1927  /* System generated locals */
1928  int i__1, i__2;
1929 
1930  /* Local variables */
1931  static double a2ro1, aion, a2ro2, e2ro1[3], e2ro2[3], e2ro3[3], e2ro4[
1932  3], eion[4], qatt[8008] /* was [4][2002] */, a2ro3, a2ro4;
1933  static int i__, j, k;
1934  static double e1[6], e2[6], e3[6], e4[6], qdrot[2002], q1[12012]
1935  /* was [6][2002] */, q2[12012] /* was [6][2002] */, q3[12012]
1936  /* was [6][2002] */, q4[12012] /* was [6][2002] */, qqrot[
1937  2002], eg, en, ei1[24], ei2[24], ei3[24], ei4[24], virial1,
1938  virial2, virial3, virial4, ain1, ain2, ain3;
1939  static double ain4;
1940 
1941 
1942 
1943 
1944 /* ---------------------------------------------------------------------
1945 */
1946 
1947 /* SUBROUTINE MIXER FILLS ARRAYS USED IN SUBROUTINE F0CALC */
1948 /* CAN HAVE A MIXTURE OF UP TO 4 GASES */
1949 
1950 
1951 /* ---------------------------------------------------------------------
1952 */
1953 
1954  mix3_1.nin1 = 0;
1955  mix3_1.nin2 = 0;
1956  mix3_1.nin3 = 0;
1957  mix3_1.nin4 = 0;
1958  mix4_1.n2ro1 = 0;
1959  mix4_1.n2ro2 = 0;
1960  mix4_1.n2ro3 = 0;
1961  mix4_1.n2ro4 = 0;
1962  sprintf(names_1.name1, "000000000000000");
1963  sprintf(names_1.name2, "000000000000000");
1964  sprintf(names_1.name3, "000000000000000");
1965  sprintf(names_1.name4, "000000000000000");
1966  for (j = 1; j <= 6; ++j) {
1967  for (i__ = 1; i__ <= 2002; ++i__) {
1968  q1[j + i__ * 6 - 7] = (float)0.;
1969  q2[j + i__ * 6 - 7] = (float)0.;
1970  q3[j + i__ * 6 - 7] = (float)0.;
1971  q4[j + i__ * 6 - 7] = (float)0.;
1972 /* L10: */
1973  }
1974  e1[j - 1] = (float)0.;
1975  e2[j - 1] = (float)0.;
1976  e3[j - 1] = (float)0.;
1977 /* L20: */
1978  e4[j - 1] = (float)0.;
1979  }
1980 
1981 /* CALL GAS CROSS-SECTIONS */
1982 
1983  gas1_(q1, mix1_1.qin1, mix4_1.q2ro1, &mix3_1.nin1, &mix4_1.n2ro1, e1, ei1,
1984  e2ro1, names_1.name1, &virial1, &c__0, 15L);
1985  if (inpt_1.ngas == 1) {
1986  goto L100;
1987  }
1988  gas2_(q2, mix1_1.qin2, mix4_1.q2ro2, &mix3_1.nin2, &mix4_1.n2ro2, e2, ei2,
1989  e2ro2, names_1.name2, &virial2, &c__0, 15L);
1990  if (inpt_1.ngas == 2) {
1991  goto L100;
1992  }
1993  gas3_(q3, mix1_1.qin3, mix4_1.q2ro3, &mix3_1.nin3, &mix4_1.n2ro3, e3, ei3,
1994  e2ro3, names_1.name3, &virial3, &c__0, 15L);
1995  if (inpt_1.ngas == 3) {
1996  goto L100;
1997  }
1998  gas4_(q4, mix1_1.qin4, mix4_1.q2ro4, &mix3_1.nin4, &mix4_1.n2ro4, e4, ei4,
1999  e2ro4, names_1.name4, &virial4, &c__0, 15L);
2000 L100:
2001 
2002 /* ------------------------------------------------------------------- */
2003 /* CORRECTION FOR NUMBER DENSITY DUE TO SECOND VIRIAL COEFFICIENT TO */
2004 /* BE ENTERED HERE (NOT YET IMPLEMENTED) . IMPORTANT FOR HIGH PRESSURE. */
2005 /* -------------------------------------------------------------------- */
2006 
2007  i__1 = inpt_1.nstep1 + 1;
2008  for (i__ = 1; i__ <= i__1; ++i__) {
2009  mix2_1.qtot[i__ - 1] = ratio_1.an1 * q1[i__ * 6 - 6] + ratio_1.an2 *
2010  q2[i__ * 6 - 6] + ratio_1.an3 * q3[i__ * 6 - 6] + ratio_1.an4
2011  * q4[i__ * 6 - 6];
2012  mix2_1.qel[i__ - 1] = ratio_1.an1 * q1[i__ * 6 - 5] + ratio_1.an2 *
2013  q2[i__ * 6 - 5] + ratio_1.an3 * q3[i__ * 6 - 5] + ratio_1.an4
2014  * q4[i__ * 6 - 5];
2015  mix1_1.qelm[i__ - 1] = ratio_1.an1 * q1[i__ * 6 - 5] * e1[1] +
2016  ratio_1.an2 * q2[i__ * 6 - 5] * e2[1] + ratio_1.an3 * q3[i__ *
2017  6 - 5] * e3[1] + ratio_1.an4 * q4[i__ * 6 - 5] * e4[1];
2018 
2019  qqrot[i__ - 1] = ratio_1.an1 * q1[i__ * 6 - 2] * e1[4] + ratio_1.an2 *
2020  q2[i__ * 6 - 2] * e2[4] + ratio_1.an3 * q3[i__ * 6 - 2] * e3[
2021  4] + ratio_1.an4 * q4[i__ * 6 - 2] * e4[4];
2022  qqrot[i__ - 1] *= (float)4.;
2023 
2024  en = mix2_1.e[i__ - 1] * (float)2.;
2025  if (mix2_1.e[i__ - 1] == (float)0.) {
2026  en = (float)1e-8;
2027  }
2028  qdrot[i__ - 1] = (float)0.;
2029  if (e1[5] == (float)0.) {
2030  goto L170;
2031  }
2032  qdrot[i__ - 1] = ratio_1.an1 * q1[i__ * 6 - 1] * (float)2. * e1[5] *
2033  inpt_1.ary * ::log(en / ::sqrt(e1[5] * inpt_1.akt));
2034 L170:
2035  if (e2[5] == (float)0.) {
2036  goto L171;
2037  }
2038  qdrot[i__ - 1] += ratio_1.an2 * q2[i__ * 6 - 1] * (float)2. * e2[5] *
2039  inpt_1.ary * ::log(en / ::sqrt(e2[5] * inpt_1.akt));
2040 L171:
2041  if (e3[5] == (float)0.) {
2042  goto L172;
2043  }
2044  qdrot[i__ - 1] += ratio_1.an3 * q3[i__ * 6 - 1] * (float)2. * e3[5] *
2045  inpt_1.ary * ::log(en / ::sqrt(e3[5] * inpt_1.akt));
2046 L172:
2047  if (e4[5] == (float)0.) {
2048  goto L173;
2049  }
2050  qdrot[i__ - 1] += ratio_1.an4 * q4[i__ * 6 - 1] * (float)2. * e4[5] *
2051  inpt_1.ary * ::log(en / ::sqrt(e4[5] * inpt_1.akt));
2052 
2053 L173:
2054  mix1_1.qelm[i__ - 1] = mix1_1.qelm[i__ - 1] * mix2_1.e[i__ - 1] *
2055  mix2_1.e[i__ - 1] + qqrot[i__ - 1] * mix2_1.e[i__ - 1] +
2056  qdrot[i__ - 1];
2057 
2058  mix1_1.qion[(i__ << 2) - 4] = q1[i__ * 6 - 4] * ratio_1.an1 *
2059  mix2_1.e[i__ - 1];
2060  mix1_1.qion[(i__ << 2) - 3] = q2[i__ * 6 - 4] * ratio_1.an2 *
2061  mix2_1.e[i__ - 1];
2062  mix1_1.qion[(i__ << 2) - 2] = q3[i__ * 6 - 4] * ratio_1.an3 *
2063  mix2_1.e[i__ - 1];
2064  mix1_1.qion[(i__ << 2) - 1] = q4[i__ * 6 - 4] * ratio_1.an4 *
2065  mix2_1.e[i__ - 1];
2066  qatt[(i__ << 2) - 4] = q1[i__ * 6 - 3] * ratio_1.an1 * mix2_1.e[i__ -
2067  1];
2068  qatt[(i__ << 2) - 3] = q2[i__ * 6 - 3] * ratio_1.an2 * mix2_1.e[i__ -
2069  1];
2070  qatt[(i__ << 2) - 2] = q3[i__ * 6 - 3] * ratio_1.an3 * mix2_1.e[i__ -
2071  1];
2072  qatt[(i__ << 2) - 1] = q4[i__ * 6 - 3] * ratio_1.an4 * mix2_1.e[i__ -
2073  1];
2074 
2075  if (mix3_1.nin1 == 0) {
2076  goto L1810;
2077  }
2078  i__2 = mix3_1.nin1;
2079  for (j = 1; j <= i__2; ++j) {
2080 /* L181: */
2081  mix1_1.qin1[j + i__ * 24 - 25] = mix1_1.qin1[j + i__ * 24 - 25] *
2082  ratio_1.an1 * mix2_1.e[i__ - 1];
2083  }
2084 L1810:
2085  if (mix3_1.nin2 == 0) {
2086  goto L1820;
2087  }
2088  i__2 = mix3_1.nin2;
2089  for (j = 1; j <= i__2; ++j) {
2090 /* L182: */
2091  mix1_1.qin2[j + i__ * 24 - 25] = mix1_1.qin2[j + i__ * 24 - 25] *
2092  ratio_1.an2 * mix2_1.e[i__ - 1];
2093  }
2094 L1820:
2095  if (mix3_1.nin3 == 0) {
2096  goto L1830;
2097  }
2098  i__2 = mix3_1.nin3;
2099  for (j = 1; j <= i__2; ++j) {
2100 /* L183: */
2101  mix1_1.qin3[j + i__ * 24 - 25] = mix1_1.qin3[j + i__ * 24 - 25] *
2102  ratio_1.an3 * mix2_1.e[i__ - 1];
2103  }
2104 L1830:
2105  if (mix3_1.nin4 == 0) {
2106  goto L1840;
2107  }
2108  i__2 = mix3_1.nin4;
2109  for (j = 1; j <= i__2; ++j) {
2110 /* L184: */
2111  mix1_1.qin4[j + i__ * 24 - 25] = mix1_1.qin4[j + i__ * 24 - 25] *
2112  ratio_1.an4 * mix2_1.e[i__ - 1];
2113  }
2114 
2115 L1840:
2116  for (k = 1; k <= 2; ++k) {
2117  if (mix4_1.n2ro1 == 0) {
2118  goto L1910;
2119  }
2120  i__2 = mix4_1.n2ro1;
2121  for (j = 1; j <= i__2; ++j) {
2122 /* L191: */
2123  mix4_1.q2ro1[k + ((j + i__ * 3) << 1) - 9] = mix4_1.q2ro1[k + (
2124  (j + i__ * 3) << 1) - 9] * ratio_1.an1 * mix2_1.e[i__ -
2125  1];
2126  }
2127 L1910:
2128  if (mix4_1.n2ro2 == 0) {
2129  goto L1920;
2130  }
2131  i__2 = mix4_1.n2ro2;
2132  for (j = 1; j <= i__2; ++j) {
2133 /* L192: */
2134  mix4_1.q2ro2[k + ((j + i__ * 3) << 1) - 9] = mix4_1.q2ro2[k + (
2135  (j + i__ * 3) << 1) - 9] * ratio_1.an2 * mix2_1.e[i__ -
2136  1];
2137  }
2138 L1920:
2139  if (mix4_1.n2ro3 == 0) {
2140  goto L1930;
2141  }
2142  i__2 = mix4_1.n2ro3;
2143  for (j = 1; j <= i__2; ++j) {
2144 /* L193: */
2145  mix4_1.q2ro3[k + ((j + i__ * 3) << 1) - 9] = mix4_1.q2ro3[k + (
2146  (j + i__ * 3) << 1) - 9] * ratio_1.an3 * mix2_1.e[i__ -
2147  1];
2148  }
2149 L1930:
2150  if (mix4_1.n2ro4 == 0) {
2151  goto L195;
2152  }
2153  i__2 = mix4_1.n2ro4;
2154  for (j = 1; j <= i__2; ++j) {
2155 /* L194: */
2156  mix4_1.q2ro4[k + ((j + i__ * 3) << 1) - 9] = mix4_1.q2ro4[k + (
2157  (j + i__ * 3) << 1) - 9] * ratio_1.an4 * mix2_1.e[i__ -
2158  1];
2159  }
2160 L195:
2161  ;
2162  }
2163 
2164 
2165  mix2_1.qrel[i__ - 1] = (float)0.;
2166  mix1_1.qsatt[i__ - 1] = (float)0.;
2167  mix1_1.qsum[i__ - 1] = (float)0.;
2168  i__2 = inpt_1.ngas;
2169  for (j = 1; j <= i__2; ++j) {
2170  mix1_1.qsum[i__ - 1] = mix1_1.qsum[i__ - 1] + mix1_1.qion[j + (
2171  i__ << 2) - 5] + qatt[j + (i__ << 2) - 5];
2172  mix1_1.qsatt[i__ - 1] += qatt[j + (i__ << 2) - 5];
2173 /* L90: */
2174  mix2_1.qrel[i__ - 1] = mix2_1.qrel[i__ - 1] + mix1_1.qion[j + (
2175  i__ << 2) - 5] - qatt[j + (i__ << 2) - 5];
2176  }
2177 
2178  if (mix3_1.nin1 == 0) {
2179  goto L910;
2180  }
2181  i__2 = mix3_1.nin1;
2182  for (j = 1; j <= i__2; ++j) {
2183 /* L91: */
2184  mix1_1.qsum[i__ - 1] += mix1_1.qin1[j + i__ * 24 - 25];
2185  }
2186 L910:
2187  if (mix3_1.nin2 == 0) {
2188  goto L920;
2189  }
2190  i__2 = mix3_1.nin2;
2191  for (j = 1; j <= i__2; ++j) {
2192 /* L92: */
2193  mix1_1.qsum[i__ - 1] += mix1_1.qin2[j + i__ * 24 - 25];
2194  }
2195 L920:
2196  if (mix3_1.nin3 == 0) {
2197  goto L930;
2198  }
2199  i__2 = mix3_1.nin3;
2200  for (j = 1; j <= i__2; ++j) {
2201 /* L93: */
2202  mix1_1.qsum[i__ - 1] += mix1_1.qin3[j + i__ * 24 - 25];
2203  }
2204 L930:
2205  if (mix3_1.nin4 == 0) {
2206  goto L940;
2207  }
2208  i__2 = mix3_1.nin4;
2209  for (j = 1; j <= i__2; ++j) {
2210 /* L94: */
2211  mix1_1.qsum[i__ - 1] += mix1_1.qin4[j + i__ * 24 - 25];
2212  }
2213 L940:
2214  eg = mix2_1.e[i__ - 1];
2215  if (mix2_1.e[i__ - 1] == (float)0.) {
2216  eg = (float)1e-7;
2217  }
2218  mix2_1.qinel[i__ - 1] = mix1_1.qsum[i__ - 1] / eg;
2219 
2220 /* L200: */
2221  }
2222 
2223  eion[0] = e1[2];
2224  eion[1] = e2[2];
2225  eion[2] = e3[2];
2226  eion[3] = e4[2];
2227  i__1 = inpt_1.ngas;
2228  for (j = 1; j <= i__1; ++j) {
2229  aion = eion[j - 1] / inpt_1.estep;
2230  mix3_1.lion[j - 1] = (int) aion;
2231 /* L300: */
2232  mix3_1.alion[j - 1] = aion - mix3_1.lion[j - 1];
2233  }
2234 
2235  if (mix3_1.nin1 == 0) {
2236  goto L4000;
2237  }
2238  i__1 = mix3_1.nin1;
2239  for (j = 1; j <= i__1; ++j) {
2240  ain1 = ei1[j - 1] / inpt_1.estep;
2241  mix3_1.lin1[j - 1] = (int) ain1;
2242 /* L400: */
2243  mix3_1.alin1[j - 1] = ain1 - mix3_1.lin1[j - 1];
2244  }
2245 L4000:
2246  if (mix3_1.nin2 == 0) {
2247  goto L4100;
2248  }
2249  i__1 = mix3_1.nin2;
2250  for (j = 1; j <= i__1; ++j) {
2251  ain2 = ei2[j - 1] / inpt_1.estep;
2252  mix3_1.lin2[j - 1] = (int) ain2;
2253 /* L410: */
2254  mix3_1.alin2[j - 1] = ain2 - mix3_1.lin2[j - 1];
2255  }
2256 L4100:
2257  if (mix3_1.nin3 == 0) {
2258  goto L4200;
2259  }
2260  i__1 = mix3_1.nin3;
2261  for (j = 1; j <= i__1; ++j) {
2262  ain3 = ei3[j - 1] / inpt_1.estep;
2263  mix3_1.lin3[j - 1] = (int) ain3;
2264 /* L420: */
2265  mix3_1.alin3[j - 1] = ain3 - mix3_1.lin3[j - 1];
2266  }
2267 L4200:
2268  if (mix3_1.nin4 == 0) {
2269  goto L4300;
2270  }
2271  i__1 = mix3_1.nin4;
2272  for (j = 1; j <= i__1; ++j) {
2273  ain4 = ei4[j - 1] / inpt_1.estep;
2274  mix3_1.lin4[j - 1] = (int) ain4;
2275 /* L430: */
2276  mix3_1.alin4[j - 1] = ain4 - mix3_1.lin4[j - 1];
2277  }
2278 L4300:
2279 
2280  if (mix4_1.n2ro1 == 0) {
2281  goto L5000;
2282  }
2283  i__1 = mix4_1.n2ro1;
2284  for (j = 1; j <= i__1; ++j) {
2285  a2ro1 = e2ro1[j - 1] / inpt_1.estep;
2286  mix4_1.l2ro1[j - 1] = (int) a2ro1;
2287 /* L500: */
2288  mix4_1.al2ro1[j - 1] = a2ro1 - mix4_1.l2ro1[j - 1];
2289  }
2290 L5000:
2291  if (mix4_1.n2ro2 == 0) {
2292  goto L5100;
2293  }
2294  i__1 = mix4_1.n2ro2;
2295  for (j = 1; j <= i__1; ++j) {
2296  a2ro2 = e2ro2[j - 1] / inpt_1.estep;
2297  mix4_1.l2ro2[j - 1] = (int) a2ro2;
2298 /* L510: */
2299  mix4_1.al2ro2[j - 1] = a2ro2 - mix4_1.l2ro2[j - 1];
2300  }
2301 L5100:
2302  if (mix4_1.n2ro3 == 0) {
2303  goto L5200;
2304  }
2305  i__1 = mix4_1.n2ro3;
2306  for (j = 1; j <= i__1; ++j) {
2307  a2ro3 = e2ro3[j - 1] / inpt_1.estep;
2308  mix4_1.l2ro3[j - 1] = (int) a2ro3;
2309 /* L520: */
2310  mix4_1.al2ro3[j - 1] = a2ro3 - mix4_1.l2ro3[j - 1];
2311  }
2312 L5200:
2313  if (mix4_1.n2ro4 == 0) {
2314  goto L5300;
2315  }
2316  i__1 = mix4_1.n2ro4;
2317  for (j = 1; j <= i__1; ++j) {
2318  a2ro4 = e2ro4[j - 1] / inpt_1.estep;
2319  mix4_1.l2ro4[j - 1] = (int) a2ro4;
2320 /* L530: */
2321  mix4_1.al2ro4[j - 1] = a2ro4 - mix4_1.l2ro4[j - 1];
2322  }
2323 L5300:
2324  return 0;
2325 } /* mixer_ */
2326 
2327 int StFtpcMagboltz1::nalpha_()
2328 {
2329  /* System generated locals */
2330  int i__1;
2331  double d__1;
2332 
2333  /* Local variables */
2334  static double velx, vely, velz, vtot, velx1, vely1, velz1;
2335  static int j;
2336  static double vtot1, ratei, cjk, dss, dxx, sum, dyy, dzz;
2337 
2338 
2339  inpt_1.alpold = inpt_1.alpnew;
2340  inpt_1.alpoax = inpt_1.alpnax;
2341  inpt_1.alpoay = inpt_1.alpnay;
2342  inpt_1.alpoaz = inpt_1.alpnaz;
2343  i__1 = inpt_1.nstep1;
2344  for (j = 1; j <= i__1; ++j) {
2345 /* L50: */
2346  sint_1.simf[j - 1] = mix2_1.e[j - 1] * f0c_1.f[j - 1] / (mag_1.denom[
2347  j - 1] * mix2_1.qtot[j - 1]);
2348  }
2349  simp_(&sum);
2350  dxx = sum * mag_1.eovm / (float)3.;
2351  i__1 = inpt_1.nstep1;
2352  for (j = 1; j <= i__1; ++j) {
2353 /* L60: */
2354  sint_1.simf[j - 1] = mix2_1.e[j - 1] * f0c_1.f[j - 1] * mag_1.sod2[j
2355  - 1] / mix2_1.qtot[j - 1];
2356  }
2357  simp_(&sum);
2358  dyy = sum * mag_1.eovm / (float)3.;
2359  i__1 = inpt_1.nstep1;
2360  for (j = 1; j <= i__1; ++j) {
2361 /* L70: */
2362  sint_1.simf[j - 1] = mix2_1.e[j - 1] * f0c_1.f[j - 1] * mag_1.cod2[j
2363  - 1] / mix2_1.qtot[j - 1];
2364  }
2365  simp_(&sum);
2366  dzz = sum * mag_1.eovm / (float)3.;
2367  dss = (dzz + dyy + dxx) / (float)3.;
2368  i__1 = inpt_1.nstep1;
2369  for (j = 1; j <= i__1; ++j) {
2370 /* L100: */
2371  sint_1.simf[j - 1] = mag_1.cod2[j - 1] * mag_1.qfemag[j - 1] *
2372  f0c_1.df[j - 1] * f0c_1.f[j - 1];
2373  }
2374  simp_(&sum);
2375  velz1 = -sum * mag_1.eovm;
2376  i__1 = inpt_1.nstep1;
2377  for (j = 1; j <= i__1; ++j) {
2378 /* L200: */
2379  sint_1.simf[j - 1] = mag_1.scd[j - 1] * mag_1.qfemag[j - 1] *
2380  f0c_1.df[j - 1] * f0c_1.f[j - 1];
2381  }
2382  simp_(&sum);
2383  vely1 = -sum * mag_1.eovm;
2384  i__1 = inpt_1.nstep1;
2385  for (j = 1; j <= i__1; ++j) {
2386 /* L300: */
2387  sint_1.simf[j - 1] = mag_1.sod[j - 1] * mag_1.qfemag[j - 1] *
2388  f0c_1.df[j - 1] * f0c_1.f[j - 1];
2389  }
2390  simp_(&sum);
2391  velx1 = -sum * mag_1.eovm;
2392  vtot1 = ::sqrt(velx1 * velx1 + vely1 * vely1 + velz1 * velz1);
2393  i__1 = inpt_1.nstep1;
2394  for (j = 1; j <= i__1; ++j) {
2395 /* L10: */
2396  sint_1.simf[j - 1] = mix2_1.qrel[j - 1] * f0c_1.f[j - 1];
2397  }
2398  simp_(&sum);
2399  ratei = sum * mag_1.eovm;
2400 /* Computing 2nd power */
2401  d__1 = vtot1 / (dss * (float)2.);
2402  cjk = d__1 * d__1 - ratei / dss;
2403  if (cjk < (float)0.) {
2404  goto L15;
2405  }
2406  inpt_1.alpnew = vtot1 / (dss * (float)2.) - ::sqrt(cjk);
2407  goto L16;
2408 L15:
2409  inpt_1.alpnew = ratei / vtot1;
2410 L16:
2411  i__1 = inpt_1.nstep1;
2412  for (j = 1; j <= i__1; ++j) {
2413 /* L20: */
2414  sint_1.simf[j - 1] = mag_1.cod2[j - 1] * (mag_1.qfemag[j - 1] *
2415  f0c_1.df[j - 1] - mag_1.qef[j - 1] * inpt_1.alpnew) * f0c_1.f[
2416  j - 1];
2417  }
2418  simp_(&sum);
2419  velz = -sum * mag_1.eovm;
2420  i__1 = inpt_1.nstep1;
2421  for (j = 1; j <= i__1; ++j) {
2422 /* L30: */
2423  sint_1.simf[j - 1] = mag_1.scd[j - 1] * (mag_1.qfemag[j - 1] *
2424  f0c_1.df[j - 1] - mag_1.qef[j - 1] * inpt_1.alpnew) * f0c_1.f[
2425  j - 1];
2426  }
2427  simp_(&sum);
2428  vely = -sum * mag_1.eovm;
2429  i__1 = inpt_1.nstep1;
2430  for (j = 1; j <= i__1; ++j) {
2431 /* L40: */
2432  sint_1.simf[j - 1] = mag_1.sod[j - 1] * (mag_1.qfemag[j - 1] *
2433  f0c_1.df[j - 1] - mag_1.qef[j - 1] * inpt_1.alpnew) * f0c_1.f[
2434  j - 1];
2435  }
2436  simp_(&sum);
2437  velx = -sum * mag_1.eovm;
2438  vtot = ::sqrt(velx * velx + vely * vely + velz * velz);
2439  if (cjk < (float)0.) {
2440  printf("WARNING NALPHA DID NOT CONVERGE (SUBROUTINE NALPHA)\n");
2441  }
2442  inpt_1.alpnax = velx1 * inpt_1.alpnew / vtot1;
2443  inpt_1.alpnay = vely1 * inpt_1.alpnew / vtot1;
2444  inpt_1.alpnaz = velz1 * inpt_1.alpnew / vtot1;
2445  printf("NEW ALPHA =%f / CM. OLD ALPHA =%f /CM. VELZ =%f VTOT =%f CM./SEC.\n",inpt_1.alpnew,inpt_1.alpold,velz,vtot);
2446  return 0;
2447 } /* nalpha_ */
2448 
2449 
2450 int StFtpcMagboltz1::output_(int *n)
2451 {
2452 
2453  /* System generated locals */
2454  int i__1;
2455  double d__1;
2456 
2457  /* Local variables */
2458  static double faci;
2459  static double velx, vely, vtot, velx1, vely1, velz1;
2460  static int j;
2461  static double vtot1, vtot2, fudge, ratei, dovmb, rteel, vtot12, dl,
2462  ri, wm, select, erootf[2002], colrte, rteinl, ratatt, alpatt,
2463  dlovmb, fac, cjk, dss, dxx, sum, dyy, dyz, dzz, dxz;
2464 
2465 
2466 
2467  i__1 = inpt_1.nstep1;
2468  for (j = 1; j <= i__1; ++j) {
2469 /* L3: */
2470  sint_1.simf[j - 1] = mag_1.cod2[j - 1] * (mag_1.qfemag[j - 1] *
2471  f0c_1.df[j - 1] * f0c_1.f[j - 1] + mag_1.emag * (float).4 * (
2472  mix2_1.e[j - 1] * f2c_1.df2[j - 1] + f2c_1.f2[j - 1] * (float)
2473  1.5) / (mix2_1.qtot[j - 1] * (float)3.));
2474  }
2475  simp_(&sum);
2476  velz1 = -sum * mag_1.eovm;
2477  if (*n == 0 && inpt_1.i2type == 1) {
2478  printf("COLLISIONS OF 2ND KIND INCLUDED.\n");
2479  }
2480  if (*n == 0) {
2481  printf("INTERMEDIATE VALUE FOR DRIFT VELOCITY = %f\n",velz1);
2482  }
2483  if (*n == 0) {
2484  return 0;
2485  }
2486  i__1 = inpt_1.nstep1;
2487  for (j = 1; j <= i__1; ++j) {
2488 /* L5: */
2489  sint_1.simf[j - 1] = mix2_1.e[j - 1] * f0c_1.f[j - 1] / (mag_1.denom[
2490  j - 1] * mix2_1.qtot[j - 1]);
2491  }
2492  simp_(&sum);
2493  dxx = sum * mag_1.eovm / (float)3.;
2494  i__1 = inpt_1.nstep1;
2495  for (j = 1; j <= i__1; ++j) {
2496 /* L6: */
2497  sint_1.simf[j - 1] = mix2_1.e[j - 1] * f0c_1.f[j - 1] * mag_1.sod2[j
2498  - 1] / mix2_1.qtot[j - 1];
2499  }
2500  simp_(&sum);
2501  dyy = sum * mag_1.eovm / (float)3.;
2502  i__1 = inpt_1.nstep1;
2503  for (j = 1; j <= i__1; ++j) {
2504 /* L7: */
2505  sint_1.simf[j - 1] = mix2_1.e[j - 1] * f0c_1.f[j - 1] * mag_1.cod2[j
2506  - 1] / mix2_1.qtot[j - 1];
2507  }
2508  simp_(&sum);
2509  dzz = sum * mag_1.eovm / (float)3.;
2510  dss = (dzz + dyy + dxx) / (float)3.;
2511  i__1 = inpt_1.nstep1;
2512  for (j = 1; j <= i__1; ++j) {
2513 /* L10: */
2514  sint_1.simf[j - 1] = mix2_1.qrel[j - 1] * f0c_1.f[j - 1];
2515  }
2516  simp_(&sum);
2517  ratei = sum * mag_1.eovm;
2518  i__1 = inpt_1.nstep1;
2519  for (j = 1; j <= i__1; ++j) {
2520 /* L11: */
2521  sint_1.simf[j - 1] = mag_1.cod2[j - 1] * (mag_1.qfemag[j - 1] *
2522  f0c_1.df[j - 1] * f0c_1.f[j - 1] + mag_1.emag * (float).4 * (
2523  mix2_1.e[j - 1] * f2c_1.df2[j - 1] + f2c_1.f2[j - 1] * (float)
2524  1.5) / (mix2_1.qtot[j - 1] * (float)3.));
2525  }
2526  simp_(&sum);
2527  velz1 = -sum * mag_1.eovm;
2528  i__1 = inpt_1.nstep1;
2529  for (j = 1; j <= i__1; ++j) {
2530 /* L12: */
2531  sint_1.simf[j - 1] = mag_1.scd[j - 1] * (mag_1.qfemag[j - 1] *
2532  f0c_1.df[j - 1] * f0c_1.f[j - 1] + mag_1.emag * (float).4 * (
2533  mix2_1.e[j - 1] * f2c_1.df2[j - 1] + f2c_1.f2[j - 1] * (float)
2534  1.5) / (mix2_1.qtot[j - 1] * (float)3.));
2535  }
2536  simp_(&sum);
2537  vely1 = -sum * mag_1.eovm;
2538  i__1 = inpt_1.nstep1;
2539  for (j = 1; j <= i__1; ++j) {
2540 /* L13: */
2541  sint_1.simf[j - 1] = mag_1.sod[j - 1] * (mag_1.qfemag[j - 1] *
2542  f0c_1.df[j - 1] * f0c_1.f[j - 1] + mag_1.emag * (float).4 * (
2543  mix2_1.e[j - 1] * f2c_1.df2[j - 1] + f2c_1.f2[j - 1] * (float)
2544  1.5) / (mix2_1.qtot[j - 1] * (float)3.));
2545  }
2546  simp_(&sum);
2547  velx1 = -sum * mag_1.eovm;
2548  vtot12 = velx1 * velx1 + vely1 * vely1 + velz1 * velz1;
2549  vtot1 = ::sqrt(vtot12);
2550 /* Computing 2nd power */
2551  d__1 = vtot1 / (dss * (float)2.);
2552  cjk = d__1 * d__1 - ratei / dss;
2553  if (cjk < (float)0.) {
2554  printf("WARNING ALPHA DID NOT CONVERGE USED ALPHA = RATE/VEL0CITY\n");
2555  }
2556  if (cjk < (float)0. || ratei == (float)0.) {
2557  goto L15;
2558  }
2559  inpt_1.alpha = vtot1 / (dss * (float)2.) - ::sqrt(cjk);
2560  goto L16;
2561 L15:
2562  inpt_1.alpha = ratei / vtot1;
2563 L16:
2564  i__1 = inpt_1.nstep1;
2565  for (j = 1; j <= i__1; ++j) {
2566 /* L20: */
2567  sint_1.simf[j - 1] = mag_1.cod2[j - 1] * (mag_1.qfemag[j - 1] *
2568  f0c_1.df[j - 1] * f0c_1.f[j - 1] + mag_1.emag * (float).4 * (
2569  mix2_1.e[j - 1] * f2c_1.df2[j - 1] + f2c_1.f2[j - 1] * (float)
2570  1.5) / (mix2_1.qtot[j - 1] * (float)3.) - mag_1.qef[j - 1] *
2571  inpt_1.alpha * (f0c_1.f[j - 1] + f2c_1.f2[j - 1] * (float).4))
2572  ;
2573  }
2574  simp_(&sum);
2575  mk_1.velz = -sum * mag_1.eovm;
2576  i__1 = inpt_1.nstep1;
2577  for (j = 1; j <= i__1; ++j) {
2578 /* L30: */
2579  sint_1.simf[j - 1] = mag_1.scd[j - 1] * (mag_1.qfemag[j - 1] *
2580  f0c_1.df[j - 1] * f0c_1.f[j - 1] + mag_1.emag * (float).4 * (
2581  mix2_1.e[j - 1] * f2c_1.df2[j - 1] + f2c_1.f2[j - 1] * (float)
2582  1.5) / (mix2_1.qtot[j - 1] * (float)3.) - mag_1.qef[j - 1] *
2583  inpt_1.alpha * (f0c_1.f[j - 1] + f2c_1.f2[j - 1] * (float).4))
2584  ;
2585  }
2586  simp_(&sum);
2587  vely = -sum * mag_1.eovm;
2588  i__1 = inpt_1.nstep1;
2589  for (j = 1; j <= i__1; ++j) {
2590 /* L40: */
2591  sint_1.simf[j - 1] = mag_1.sod[j - 1] * (mag_1.qfemag[j - 1] *
2592  f0c_1.df[j - 1] * f0c_1.f[j - 1] + mag_1.emag * (float).4 * (
2593  mix2_1.e[j - 1] * f2c_1.df2[j - 1] + f2c_1.f2[j - 1] * (float)
2594  1.5) / (mix2_1.qtot[j - 1] * (float)3.) - mag_1.qef[j - 1] *
2595  inpt_1.alpha * (f0c_1.f[j - 1] + f2c_1.f2[j - 1] * (float).4))
2596  ;
2597  }
2598  simp_(&sum);
2599  velx = -sum * mag_1.eovm;
2600  vtot2 = velx * velx + vely * vely + mk_1.velz * mk_1.velz;
2601  vtot = ::sqrt(vtot2);
2602  ri = inpt_1.alpha * vtot1;
2603  i__1 = inpt_1.nstep1;
2604  for (j = 1; j <= i__1; ++j) {
2605 /* L50: */
2606  sint_1.simf[j - 1] = mix2_1.e[j - 1] * h1c_1.h1[j - 1] / mag_1.denom[
2607  j - 1];
2608  }
2609  simp_(&sum);
2610  dxx = sum * mag_1.eovm / (float)3.;
2611  i__1 = inpt_1.nstep1;
2612  for (j = 1; j <= i__1; ++j) {
2613 /* L60: */
2614  sint_1.simf[j - 1] = mix2_1.e[j - 1] * h1c_1.h1[j - 1] * mag_1.sod2[j
2615  - 1];
2616  }
2617  simp_(&sum);
2618  dyy = sum * mag_1.eovm / (float)3.;
2619  i__1 = inpt_1.nstep1;
2620  for (j = 1; j <= i__1; ++j) {
2621 /* L70: */
2622  sint_1.simf[j - 1] = mix2_1.e[j - 1] * h1c_1.h1[j - 1] * mag_1.cod2[j
2623  - 1];
2624  }
2625  simp_(&sum);
2626  dzz = sum * mag_1.eovm / (float)3.;
2627  i__1 = inpt_1.nstep1;
2628  for (j = 1; j <= i__1; ++j) {
2629 /* L80: */
2630  sint_1.simf[j - 1] = mix2_1.e[j - 1] * h1c_1.h1[j - 1] * (float)2. *
2631  mag_1.scd[j - 1];
2632  }
2633  simp_(&sum);
2634  dyz = sum * mag_1.eovm / (float)3.;
2635  i__1 = inpt_1.nstep1;
2636  for (j = 1; j <= i__1; ++j) {
2637 /* L90: */
2638  sint_1.simf[j - 1] = mix2_1.e[j - 1] * h1c_1.h1[j - 1] * (float)2. *
2639  mag_1.sod[j - 1];
2640  }
2641  simp_(&sum);
2642  dxz = sum * mag_1.eovm / (float)3.;
2643  i__1 = inpt_1.nstep1;
2644  for (j = 1; j <= i__1; ++j) {
2645 /* L100: */
2646  sint_1.simf[j - 1] = mix2_1.e[j - 1] * mix2_1.eroot[j - 1] * f0c_1.f[
2647  j - 1];
2648  }
2649  simp_(&sum);
2650  mk_1.emean = sum;
2651  i__1 = inpt_1.nstep1;
2652  for (j = 1; j <= i__1; ++j) {
2653 /* L110: */
2654  sint_1.simf[j - 1] = mix2_1.eroot[j - 1] * f0c_1.f[j - 1] /
2655  mag_1.denom[j - 1];
2656  }
2657  simp_(&sum);
2658  fac = sum;
2659  i__1 = inpt_1.nstep1;
2660  for (j = 1; j <= i__1; ++j) {
2661 /* L120: */
2662  sint_1.simf[j - 1] = mix2_1.qtot[j - 1] * mix2_1.e[j - 1] * f0c_1.f[j
2663  - 1];
2664  }
2665  simp_(&sum);
2666  colrte = sum * mag_1.eovm;
2667 /* Computing 2nd power */
2668  d__1 = mag_1.wb / colrte;
2669  faci = (float)1. / (d__1 * d__1 + (float)1.);
2670  i__1 = inpt_1.nstep1;
2671  for (j = 1; j <= i__1; ++j) {
2672 /* L130: */
2673  erootf[j - 1] = mix2_1.eroot[j - 1] * f0c_1.f[j - 1];
2674  }
2675  i__1 = inpt_1.nstep1;
2676  for (j = 1; j <= i__1; ++j) {
2677 /* L140: */
2678  sint_1.simf[j - 1] = mix2_1.qinel[j - 1] * mix2_1.e[j - 1] * f0c_1.f[
2679  j - 1];
2680  }
2681  simp_(&sum);
2682  rteinl = sum * mag_1.eovm;
2683  i__1 = inpt_1.nstep1;
2684  for (j = 1; j <= i__1; ++j) {
2685 /* L150: */
2686  sint_1.simf[j - 1] = mix2_1.qel[j - 1] * mix2_1.e[j - 1] * f0c_1.f[j
2687  - 1];
2688  }
2689  simp_(&sum);
2690  rteel = sum * mag_1.eovm;
2691  i__1 = inpt_1.nstep1;
2692  for (j = 1; j <= i__1; ++j) {
2693 /* L160: */
2694  sint_1.simf[j - 1] = g1c_1.g1[j - 1] * mix2_1.e[j - 1];
2695  }
2696  simp_(&sum);
2697  dl = sum * mag_1.eovm / (float)3.;
2698 
2699  i__1 = inpt_1.nstep1;
2700  for (j = 1; j <= i__1; ++j) {
2701 /* L161: */
2702  sint_1.simf[j - 1] = f0c_1.f[j - 1] * (mix1_1.qin1[j * 24 - 22] +
2703  mix1_1.qin1[j * 24 - 23]);
2704  }
2705  simp_(&sum);
2706  select = sum * mag_1.eovm;
2707  printf("SELECTED INELASTIC FREQUENCY= %f\n", select);
2708  i__1 = inpt_1.nstep1;
2709  for (j = 1; j <= i__1; ++j) {
2710 /* L170: */
2711  sint_1.simf[j - 1] = f0c_1.f[j - 1] * mix1_1.qsatt[j - 1];
2712  }
2713  simp_(&sum);
2714  ratatt = sum * mag_1.eovm;
2715  alpatt = ratatt / vtot1;
2716  wm = (float)0.;
2717  if (mag_1.bmag > (float)0.) {
2718  wm = velx * mag_1.emag * (float)1e5 / (mk_1.velz * mag_1.bmag);
2719  }
2720  dovmb = (dzz + dyy + dxx) * mag_1.emag / (vtot1 * (float)3.);
2721  dlovmb = dl * mag_1.emag / vtot1;
2722  mk_1.angle = atan(velx / mk_1.velz) * (float)180. / acos(-1.);
2723  printf("-----------------------------------------------------------------------------------------------------------------------");
2724  if (*n == 1) {
2725  printf("FINAL VALUES WITHOUT COLLISIONS OF 2ND KIND. LORENTZ SOLUTION (L=1).\n");
2726  }
2727  if (*n == 2) {
2728  printf("FINAL VALUES WITH L=2 TERM FULLY INCLUDED IN CALCULATION.\n");
2729  }
2730  if (*n == 3) {
2731  printf("FINAL VALUES WITH COLLISIONS OF 2ND KIND. \n");
2732  }
2733  if (*n == 4) {
2734  printf("FINAL VALUES CONVERGED ON ALPHA TO ALL ORDERS. \n");
2735  }
2736  if (*n == 5) {
2737  printf("FINAL VALUES WITH HIGHER MULTIPOLE. (L=2 TERM) HIGHER TERM ONLY TO FIRST ORDER .\n");
2738  }
2739  if (*n == 6) {
2740  printf("FINAL VALUES WITH HIGHER MULTIPOLE. (L=3 TERM) HIGHER TERM ONLY TO FIRST ORDER .\n");
2741  }
2742  printf("N.B. LORENTZ ANGLE = VELX/VELZ =%f DEGREES.\n", mk_1.angle);
2743  printf("VELZ =%f VELY =%f VELX =%f VTOT =%f CM./SEC. ( WITH IONISATION )\n",mk_1.velz,vely,velx,vtot);
2744  printf("VELZ1 =%f VELY1 =%f VELX1 =%f VTOT1 =%f CM./SEC. ( WITHOUT IONISATION )\n",velz1,vely1,velx1,vtot1);
2745  printf(" DZZ =%f DYY =%f DXX =%f DYZ =%f DXZ =%f CM.**2/SEC.",dzz,dyy,dxx,dyz,dxz);
2746  printf("MEAN ELECTRON ENERGY =%f EV. DIFFUSIONMOBILITY =%f EV. = CHARACTERISTIC ENERGY.\n",mk_1.emean,dovmb);
2747  printf("NET ( IONISATION - ATTACHMENT ) RATE =%f / CM. =%f /SEC. ATTACHMENT RATE =%f / CM.\n",inpt_1.alpha,ratei,alpatt);
2748  printf("CYCLOTRON FREQUENCY =%f RADS./SEC. TOTAL COLLISION FREQUENCY =%f /SEC. 1/(1+(CYCLOTRON FREQ./COLL.FREQ.)**2) = K FACTOR =%f TRUE K FACTOR =%f\n",mag_1.wb,colrte,faci,fac);
2749  printf("MAGNETIC DRIFT VELOCITY =%f CM/SEC. (ONLY VALID FOR E AND B AT 90 DEGREES TO EACH OTHER.)\n",wm);
2750  printf("INELASTIC COLLISION FREQUENCY =%f /SEC. ELASTIC COLLISION FREQUENCY =%f /SEC.\n",rteinl,rteel);
2751  printf(" F0 =%f AT FINAL ENERGY OF %f EV.\n",f0c_1.f[inpt_1.nstep1 - 1],inpt_1.efinal);
2752  printf("LONGITUDINAL DIFFUSION =%f CM.**2/SEC. = %f EV. \n",dl,dlovmb);
2753 /* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
2754 /* convert long. diff. to um/::sqrt(cm) by W.Gong */
2755 /* also applying the correction */
2756  fudge = dl / dzz;
2757  mk_1.sig_long__ = ::sqrt(dl * (float)2. * (float)1e6 / mk_1.velz) * (float)
2758  10.;
2759  mk_1.sig_tranx__ = ::sqrt(dxx * (float)2. * (float)1e6 / mk_1.velz) * (
2760  float)10.;
2761  mk_1.sig_trany__ = ::sqrt(dyy * (float)2. * (float)1e6 / mk_1.velz) * (
2762  float)10.;
2763 /* L9921: */
2764  mk_1.btheta_mk__ = mag_1.btheta;
2765  mk_1.bmag_mk__ = mag_1.bmag;
2766  printf("PULSED TOWNSEND OR TOF IONISATION RATE=%f /SEC.\n",ri);
2767 /* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
2768 /* write out the relevant quantities into fort.10 by W.Gong */
2769 /* not needed for STAR version, output to table */
2770 /* s_wsfe(&io___58); */
2771 /* do_fio(&c__1, (char *)&mag_1.emag, (int)sizeof(double)); */
2772 /* d__1 = mk_1.velz / (float)1e6; */
2773 /* do_fio(&c__1, (char *)&d__1, (int)sizeof(double)); */
2774 /* do_fio(&c__1, (char *)&mk_1.sig_long__, (int)sizeof(double)); */
2775 /* do_fio(&c__1, (char *)&mk_1.sig_tranx__, (int)sizeof(double)); */
2776 /* do_fio(&c__1, (char *)&mk_1.sig_trany__, (int)sizeof(double)); */
2777 /* do_fio(&c__1, (char *)&mk_1.angle, (int)sizeof(double)); */
2778 /* do_fio(&c__1, (char *)&mk_1.emean, (int)sizeof(double)); */
2779 /* do_fio(&c__1, (char *)&mag_1.bmag, (int)sizeof(double)); */
2780 /* do_fio(&c__1, (char *)&mag_1.btheta, (int)sizeof(double)); */
2781 /* do_fio(&c__1, (char *)&f0c_1.f[inpt_1.nstep1 - 1], (int)sizeof( */
2782 /* double)); */
2783 /* do_fio(&c__1, (char *)&inpt_1.efinal, (int)sizeof(double)); */
2784 /* e_wsfe(); */
2785 /* cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc */
2786 /* ********** MK ********** */
2787  mk_1.amk_emag__ = mag_1.emag;
2788  printf("MK-EMAG %f\n", mk_1.amk_emag__);
2789 /* ************************ */
2790  if (*n < 5) {
2791  return 0;
2792  }
2793  return 0;
2794 } /* output_ */
2795 
2796 
2797 int StFtpcMagboltz1::prnter_()
2798 {
2799  printf("BOLTZMAN SOLUTION FOR MIXTURE OF %d GASES.\n------------------------------------------------------\n",inpt_1.ngas);
2800  printf(" GASES REQUESTED = %s %s %s %s\n",gasn_1.ngas1,gasn_1.ngas2,gasn_1.ngas3,gasn_1.ngas4);
2801  printf(" GASES USED = %s %s %s %s\n",names_1.name1,names_1.name2,names_1.name3,names_1.name4);
2802  printf(" PERCENTAGE USED = %f %f %f %f\n",ratio_1.frac1,ratio_1.frac2,ratio_1.frac3,ratio_1.frac4);
2803  printf("GAS TEMPERATURE =%f DEGREES CENTIGRADE. GAS PRESSURE = %f TORR.\n",inpt_1.tempc,inpt_1.torr);
2804  printf(" INTEGRATION FROM 0.0 TO %f EV. IN %d STEPS. CONVERGENCE ERROR LESS THAN %f\n",inpt_1.efinal,inpt_1.nstep,inpt_1.conv);
2805  if (inpt_1.i2type == 1) {
2806  printf("COLLISIONS OF SECOND TYPE INCLUDED.\n");
2807  }
2808  if (inpt_1.idlong == 1) {
2809  goto L71;
2810  }
2811  printf("PRINTOUT EVERY %d ITERATIONS. MAXIMUM NUMBER OF ITERATIONS ALLOWED =%d\nLONGITUDINAL DIFF. NOT CALCULATED.\n",inpt_1.nout,inpt_1.itmax);
2812  goto L74;
2813 L71:
2814  printf("PRINTOUT EVERY %d ITERATIONS. MAXIMUM NUMBER OF ITERATIONS ALLOWED =%d\nLONGITUDINAL DIFFUSION CALCULATED.\n",inpt_1.nout,inpt_1.itmax);
2815 L74:
2816  printf(" ELECTRIC FIELD =%f VOLTS/CM. MAGNETIC FIELD =%f KGAUSS. \nANGLE BETWEEN MAGNETIC AND ELECTRIC FIELD = %f DEGREES.\n",mag_1.emag,mag_1.bmag,mag_1.btheta);
2817  if (inpt_1.isfb == 0) {
2818  printf(" STANDARD PARAMETERISATION OF KT TERM.\n");
2819  }
2820  if (inpt_1.isfb == 1) {
2821  printf(" INELASTIC LEVELS IN KT TERM.\n");
2822  }
2823  printf(" N.B. MULTIPLY DIFFUSION COEFFICIENTS BELOW BY :- DZZ=DZZ*FUDGE, DYZ=DYZ*SQRT(FUDGE), DXZ=DXZ*SQRT(FUDGE). USE (FUDGE) FACTOR = RATIO OF LONGITUDINAL/TRANSVERSE DIFFUSION COEFFICIENTS IN ZERO MAGNETIC FIELD.\n");
2824  return 0;
2825 } /* prnter_ */
2826 
2827 int StFtpcMagboltz1::setup_(int *last)
2828 {
2829  /* Local variables */
2830  static double corr;
2831  static int i__, j;
2832  static double boltz, aj, alosch, awb;
2833 
2834 
2835 /* PHYSICAL CONSTANTS 1987 UPDATE OF TAYLOR AND COHEN */
2836 /* NEW VALUE OF BOLTZ 1988 */
2837 
2838  inpt_1.ary = (float)13.6056981;
2839  cnsts_1.pir2 = 8.79735669e-17;
2840  cnsts_1.echarg = 1.60217733e-19;
2841  cnsts_1.emass = 9.1093897e-31;
2842  cnsts_1.amu = 1.6605402e-27;
2843  boltz = 8.617343e-5;
2844  awb = 17588196200.;
2845  alosch = 2.686763e19;
2846  mag_1.eovm = ::sqrt(cnsts_1.echarg * (float)2. / cnsts_1.emass) * (float)
2847  100.;
2848 
2849 
2850 /* READ IN OUTPUT CONTROL AND INTEGRATION DATA */
2851 
2852 
2853 /* INTEGRATE TO 0.1% ACCURACY */
2854  inpt_1.conv = (float)5e-4;
2855  inpt_1.nstep = 2000;
2856 /* CALCULATE ALPHA TO 1% ACCURACY */
2857  inpt_1.conalp = (float).01;
2858  inpt_1.nout = 10;
2859  inpt_1.itmax = 1200;
2860  inpt_1.isfb = 0;
2861  inpt_1.idbug = 0;
2862 /* ********** hh ********** */
2863  if (*last < 10) {
2864  inpt_1.ngas = 4;
2865  inpt_1.i2type = 0;
2866  inpt_1.nitalp = 0;
2867  inpt_1.idlong = 1;
2868  inpt_1.lhigh = 1;
2869  }
2870 /* ********** hh ********** */
2871  printf("************************\n");
2872  printf("Here comes the final E: %f\n", inpt_1.efinal);
2873  printf("************************\n");
2874  if (inpt_1.ngas == 0) {
2875  goto L99;
2876  }
2877  inpt_1.nstep1 = inpt_1.nstep + 1;
2878  inpt_1.estep = inpt_1.efinal / inpt_1.nstep;
2879  for (i__ = 1; i__ <= 2002; ++i__) {
2880  j = i__ - 1;
2881  aj = (double) j;
2882  mix2_1.e[i__ - 1] = aj * inpt_1.estep;
2883 /* L10: */
2884  mix2_1.eroot[i__ - 1] = ::sqrt(mix2_1.e[i__ - 1]);
2885  }
2886  mix2_1.eroot[0] = 1e-10;
2887 
2888 /* GAS PARAMETERS */
2889 
2890 /* ********** hh ********** */
2891  sprintf(gasn_1.ngas1, "Ar");
2892  sprintf(gasn_1.ngas1, "CO2");
2893  sprintf(gasn_1.ngas1, "Ne");
2894  sprintf(gasn_1.ngas1, "He");
2895  ratio_1.frac1 = callpars_1.perc_ar__;
2896  ratio_1.frac2 = callpars_1.perc_co2__;
2897  ratio_1.frac3 = callpars_1.perc_ne__;
2898  ratio_1.frac4 = callpars_1.perc_he__;
2899  inpt_1.tempc = callpars_1.temperature;
2900  inpt_1.torr = callpars_1.pressure;
2901 /* ********** hh ********** */
2902  corr = inpt_1.torr * (float)2.7316 / ((inpt_1.tempc + (float)273.16) * (
2903  float)760.);
2904  ratio_1.an1 = ratio_1.frac1 * corr * alosch;
2905  ratio_1.an2 = ratio_1.frac2 * corr * alosch;
2906  ratio_1.an3 = ratio_1.frac3 * corr * alosch;
2907  ratio_1.an4 = ratio_1.frac4 * corr * alosch;
2908  inpt_1.akt = (inpt_1.tempc + (float)273.16) * boltz;
2909  ratio_1.an = corr * (float)100. * alosch;
2910 
2911 /* FIELD VALUES */
2912 
2913  mag_1.emag = callpars_1.e_magnitude__;
2914  mag_1.bmag = callpars_1.b_magnitude__;
2915  mag_1.btheta = callpars_1.b_angle__;
2916 /* HH end *************************** */
2917  *last = 0;
2918 /* ********** MK ********** */
2919  mag_1.wb = awb * mag_1.bmag;
2920  for (j = 1; j <= 2002; ++j) {
2921  f1c_1.f1[j - 1] = (float)0.;
2922  f1c_1.df1[j - 1] = (float)0.;
2923  f2c_1.f2[j - 1] = (float)0.;
2924  f2c_1.df2[j - 1] = (float)0.;
2925  f3c_1.f3[j - 1] = (float)0.;
2926  f3c_1.df3[j - 1] = (float)0.;
2927  h1c_1.h1[j - 1] = (float)0.;
2928  h1c_1.dh1[j - 1] = (float)0.;
2929  g2c_1.g2[j - 1] = (float)0.;
2930  g2c_1.dg2[j - 1] = (float)0.;
2931  g1c_1.g1[j - 1] = (float)0.;
2932  g1c_1.dg1[j - 1] = (float)0.;
2933  g0c_1.g[j - 1] = (float)0.;
2934  g0c_1.dg0[j - 1] = (float)0.;
2935  g0c_1.dg[j - 1] = (float)0.;
2936  f0c_1.f[j - 1] = (float)0.;
2937  f0c_1.df0[j - 1] = (float)0.;
2938 /* L6: */
2939  f0c_1.df[j - 1] = (float)0.;
2940  }
2941  inpt_1.alpnew = (float)0.;
2942  inpt_1.alpold = (float)0.;
2943  inpt_1.alpnax = (float)0.;
2944  inpt_1.alpnay = (float)0.;
2945  inpt_1.alpnaz = (float)0.;
2946  inpt_1.alpoax = (float)0.;
2947  inpt_1.alpoay = (float)0.;
2948  inpt_1.alpoaz = (float)0.;
2949  inpt_1.alpha = (float)0.;
2950  return 0;
2951 L99:
2952  *last = 1;
2953  return 0;
2954 } /* setup_ */
2955 
2956 int StFtpcMagboltz1::simp_(double *sum)
2957 {
2958  /* System generated locals */
2959  int i__1;
2960 
2961  /* Local variables */
2962  static double aodd, even;
2963  static int i__, n0;
2964 
2965  aodd = (float)0.;
2966  even = (float)0.;
2967  i__1 = inpt_1.nstep;
2968  for (i__ = 2; i__ <= i__1; i__ += 2) {
2969 /* L11: */
2970  aodd += sint_1.simf[i__ - 1];
2971  }
2972  n0 = inpt_1.nstep - 1;
2973  i__1 = n0;
2974  for (i__ = 3; i__ <= i__1; i__ += 2) {
2975 /* L12: */
2976  even += sint_1.simf[i__ - 1];
2977  }
2978  *sum = inpt_1.estep * (sint_1.simf[0] + sint_1.simf[inpt_1.nstep] + aodd *
2979  (float)4. + even * (float)2.) / (float)3.;
2980  return 0;
2981 } /* simp_ */
2982 int StFtpcMagboltz1::stepph_(int *l)
2983 {
2984  /* System generated locals */
2985  double d__1;
2986  int c__2=2;
2987  /* Local variables */
2988  static double frac;
2989  static int k;
2990  static double accur, dxold;
2991  static int kstep;
2992  static double dhfrs1, dhfrs2, h01, dhfnal, dhlast, dhstep, dhfrst;
2993  static int kprint;
2994  static double dum, dxx;
2995 
2996  if (*l == 2) {
2997  goto L21;
2998  }
2999  h1calc_(&c__1, &dum, &dum, &dum);
3000  return 0;
3001 /* ACCURACY 3.0% */
3002 L21:
3003  accur = (float).02;
3004  dxold = (float)0.;
3005  dhlast = (float)0.;
3006  kprint = 0;
3007  kstep = 0;
3008 /* INITIALLY USE WIDE STEPS */
3009  dhstep = (float)100.;
3010  dhfnal = (float)1.0000000000000011e-10;
3011 L146:
3012  h1calc_(&c__2, &dhfnal, &dxx, &dhfrst);
3013 L10:
3014  h01 = dhfnal;
3015  dhfrs1 = dhfrst;
3016  ++kstep;
3017  if (kstep > 9) {
3018  goto L160;
3019  }
3020  dhfnal *= dhstep;
3021  h1calc_(&c__2, &dhfnal, &dxx, &dhfrst);
3022  dhfrs2 = dhfrst;
3023  if (dhfrs1 / dhfrs2 <= (float)0.) {
3024  goto L20;
3025  }
3026  if (kstep > 1) {
3027  goto L10;
3028  }
3029  if (fabs(dhfrs1) > fabs(dhfrs2)) {
3030  goto L10;
3031  }
3032  dhstep = (float)1. / dhstep;
3033  goto L10;
3034 /* SMALLER STEPS CLOSER TO MINIMUM */
3035 L20:
3036  if (dhstep > (float)1.) {
3037  goto L55;
3038  }
3039 L25:
3040  dhstep = h01 * (float).5;
3041 L30:
3042  dhfnal = h01;
3043  dhfnal -= dhstep;
3044  ++kstep;
3045  h1calc_(&c__2, &dhfnal, &dxx, &dhfrst);
3046  if (dhlast == dhfrst) {
3047  goto L60;
3048  }
3049  dhlast = dhfrst;
3050  dhfrs2 = dhfrst;
3051  frac = (d__1 = (dxx - dxold) / dxx, fabs(d__1));
3052  if (frac < accur) {
3053  goto L70;
3054  }
3055  dxold = dxx;
3056 /* L46: */
3057  if (dhfrs1 / dhfrs2 <= (float)0.) {
3058  goto L50;
3059  }
3060  h01 = dhfnal;
3061  dhfrs1 = dhfrst;
3062 L50:
3063  dhstep *= (float).5;
3064  goto L30;
3065 L55:
3066  h01 = dhfnal;
3067  dhfrs1 = dhfrs2;
3068  goto L25;
3069 L60:
3070  printf(" WARNING TRANSVERSE DIFFUSION DID NOT CONVERGE IN %d ITERATIONS . LAST VALUE DXX=%f\n REMEDY:INCREASE COMPUTER PRECISION. IF FAIL AT FLOAT*8 THEN REDUCE UPPER INTEGRATION ENERGY LIMIT.\n",kstep,dxx);
3071  for (k = 1; k <= 2002; ++k) {
3072 /* L91: */
3073  h1c_1.h1[k - 1] = (float)0.;
3074  }
3075  return 0;
3076 L70:
3077  printf("TRANSVERSE DIFFUSION CONVERGED AFTER %d ITERATIONS.\n",kstep);
3078  return 0;
3079 L160:
3080  if (dhfnal < (float)0.) {
3081  goto L876;
3082  }
3083  dhfnal = (float)-1.0000000000000011e-10;
3084  dhstep = (float)100.;
3085  kstep = 0;
3086  goto L146;
3087 L876:
3088  printf(" WARNING TRANSVERSE DIFFUSION DID NOT APPROACH CONVGENCE.\n");
3089  return 0;
3090 } /* stepph_ */
3091 
3092 int StFtpcMagboltz1::steppr_(int *itype, int *lmax)
3093 {
3094 
3095  /* System generated locals */
3096  int i__1;
3097  double d__1;
3098 
3099  /* Local variables */
3100  static double frac;
3101  static int j, k;
3102  static double accur, dlold, gstep;
3103  static int kstep;
3104  static double eg0sum, egsum1, egsum2, g01, dl, gfinal, sum;
3105 
3106 /* ACCURACY 0.5% */
3107  accur = (float).003;
3108  dlold = (float)0.;
3109  kstep = 0;
3110 /* INITIALLY USE WIDE STEPS */
3111  gstep = (float).01;
3112  gfinal = (float).010000000000000002;
3113  g0calc_(&c__0, &gfinal, &eg0sum, lmax);
3114 L10:
3115  g01 = gfinal;
3116  egsum1 = eg0sum;
3117  ++kstep;
3118  if (kstep > 10) {
3119  goto L876;
3120  }
3121  gfinal *= gstep;
3122  g0calc_(&c__0, &gfinal, &eg0sum, lmax);
3123  egsum2 = eg0sum;
3124  if (egsum1 / egsum2 <= (float)0.) {
3125  goto L25;
3126  }
3127  goto L10;
3128 /* SMALLER STEPS CLOSER TO MINIMUM */
3129 L25:
3130  gstep = g01;
3131 L30:
3132  gstep *= (float).5;
3133  gfinal = g01;
3134  gfinal -= gstep;
3135  ++kstep;
3136  g0calc_(itype, &gfinal, &eg0sum, lmax);
3137  if (kstep > 40) {
3138  goto L60;
3139  }
3140  i__1 = inpt_1.nstep1;
3141  for (j = 1; j <= i__1; ++j) {
3142 /* L161: */
3143  sint_1.simf[j - 1] = g1c_1.g1[j - 1] * mix2_1.e[j - 1];
3144  }
3145  simp_(&sum);
3146  dl = sum * mag_1.eovm / (float)3.;
3147  egsum2 = eg0sum;
3148  frac = (d__1 = (dl - dlold) / dl, fabs(d__1));
3149  if (frac < accur) {
3150  goto L70;
3151  }
3152  dlold = dl;
3153  if (egsum1 / egsum2 <= (float)0.) {
3154  goto L30;
3155  }
3156  g01 = gfinal;
3157  egsum1 = eg0sum;
3158  goto L30;
3159 L60:
3160  printf(" WARNING LONGITUDINAL DIFFUSION DID NOT CONVERGE IN %d ITERATIONS . LAST VALUE DL =%f\n",kstep,dl);
3161  for (k = 1; k <= 2002; ++k) {
3162 /* L823: */
3163  g1c_1.g1[k - 1] = (float)0.;
3164  }
3165  return 0;
3166 L70:
3167  printf("LONGITUDINAL DIFFUSION CONVERGED AFTER %d ITERATIONS.\n",kstep);
3168  return 0;
3169 L876:
3170  printf(" WARNING LONGITUDINAL DIFFUSION DID NOT APPROACH CONVERGENCE.\n");
3171  return 0;
3172 } /* steppr_ */
3173 
3174 
3175 int StFtpcMagboltz1::type2_(double *s2sum, int *i__, int *nstep1)
3176 {
3177  /* System generated locals */
3178  int i__1;
3179 
3180  /* Local variables */
3181  static double s1sum;
3182  static int j, ldn, lup;
3183 
3184 
3185  s1sum = (float)0.;
3186  *s2sum = (float)0.;
3187  if (mix4_1.n2ro1 == 0) {
3188  goto L1000;
3189  }
3190  i__1 = mix4_1.n2ro1;
3191  for (j = 1; j <= i__1; ++j) {
3192  s1sum += f0c_1.f[*i__ - 1] * (mix4_1.q2ro1[((j + *i__ * 3) << 1) - 8] +
3193  mix4_1.q2ro1[((j + *i__ * 3) << 1) - 7]);
3194  lup = *i__ + mix4_1.l2ro1[j - 1];
3195  if (lup >= *nstep1) {
3196  goto L50;
3197  }
3198  *s2sum = *s2sum + f0c_1.f[lup - 1] * mix4_1.q2ro1[((j + lup * 3) << 1)
3199  - 8] + (f0c_1.f[lup] * mix4_1.q2ro1[((j + (lup + 1) * 3) << 1)
3200  - 8] - f0c_1.f[lup - 1] * mix4_1.q2ro1[((j + lup * 3) << 1) - 8]
3201  ) * mix4_1.al2ro1[j - 1];
3202 L50:
3203  ;
3204  }
3205  i__1 = mix4_1.n2ro1;
3206  for (j = 2; j <= i__1; ++j) {
3207  ldn = *i__ - mix4_1.l2ro1[j - 2];
3208  if (ldn < 1) {
3209  goto L100;
3210  }
3211  *s2sum = *s2sum + f0c_1.f[ldn - 1] * mix4_1.q2ro1[((j + ldn * 3) << 1)
3212  - 7] + (f0c_1.f[ldn] * mix4_1.q2ro1[((j + (ldn + 1) * 3) << 1)
3213  - 7] - f0c_1.f[ldn - 1] * mix4_1.q2ro1[((j + ldn * 3) << 1) - 7]
3214  ) * mix4_1.al2ro1[j - 1];
3215 L100:
3216  ;
3217  }
3218 L1000:
3219  if (mix4_1.n2ro2 == 0) {
3220  goto L2001;
3221  }
3222  i__1 = mix4_1.n2ro2;
3223  for (j = 1; j <= i__1; ++j) {
3224  s1sum += f0c_1.f[*i__ - 1] * (mix4_1.q2ro2[((j + *i__ * 3) << 1) - 8] +
3225  mix4_1.q2ro2[((j + *i__ * 3) << 1) - 7]);
3226  lup = *i__ + mix4_1.l2ro2[j - 1];
3227  if (lup >= *nstep1) {
3228  goto L150;
3229  }
3230  *s2sum = *s2sum + f0c_1.f[lup - 1] * mix4_1.q2ro2[((j + lup * 3) << 1)
3231  - 8] + (f0c_1.f[lup] * mix4_1.q2ro2[((j + (lup + 1) * 3) << 1)
3232  - 8] - f0c_1.f[lup - 1] * mix4_1.q2ro2[((j + lup * 3) << 1) - 8]
3233  ) * mix4_1.al2ro2[j - 1];
3234 L150:
3235  ;
3236  }
3237  i__1 = mix4_1.n2ro2;
3238  for (j = 2; j <= i__1; ++j) {
3239  ldn = *i__ - mix4_1.l2ro2[j - 2];
3240  if (ldn < 1) {
3241  goto L200;
3242  }
3243  *s2sum = *s2sum + f0c_1.f[ldn - 1] * mix4_1.q2ro2[((j + ldn * 3) << 1)
3244  - 7] + (f0c_1.f[ldn] * mix4_1.q2ro2[((j + (ldn + 1) * 3) << 1)
3245  - 7] - f0c_1.f[ldn - 1] * mix4_1.q2ro2[((j + ldn * 3) << 1) - 7]
3246  ) * mix4_1.al2ro2[j - 1];
3247 L200:
3248  ;
3249  }
3250 L2001:
3251  if (mix4_1.n2ro3 == 0) {
3252  goto L3000;
3253  }
3254  i__1 = mix4_1.n2ro3;
3255  for (j = 1; j <= i__1; ++j) {
3256  s1sum += f0c_1.f[*i__ - 1] * (mix4_1.q2ro3[((j + *i__ * 3) << 1) - 8] +
3257  mix4_1.q2ro3[((j + *i__ * 3) << 1) - 7]);
3258  lup = *i__ + mix4_1.l2ro3[j - 1];
3259  if (lup >= *nstep1) {
3260  goto L250;
3261  }
3262  *s2sum = *s2sum + f0c_1.f[lup - 1] * mix4_1.q2ro3[((j + lup * 3) << 1)
3263  - 8] + (f0c_1.f[lup] * mix4_1.q2ro3[((j + (lup + 1) * 3) << 1)
3264  - 8] - f0c_1.f[lup - 1] * mix4_1.q2ro3[((j + lup * 3) << 1) - 8]
3265  ) * mix4_1.al2ro3[j - 1];
3266 L250:
3267  ;
3268  }
3269  i__1 = mix4_1.n2ro3;
3270  for (j = 2; j <= i__1; ++j) {
3271  ldn = *i__ - mix4_1.l2ro3[j - 2];
3272  if (ldn < 1) {
3273  goto L300;
3274  }
3275  *s2sum = *s2sum + f0c_1.f[ldn - 1] * mix4_1.q2ro3[((j + ldn * 3) << 1)
3276  - 7] + (f0c_1.f[ldn] * mix4_1.q2ro3[((j + (ldn + 1) * 3) << 1)
3277  - 7] - f0c_1.f[ldn - 1] * mix4_1.q2ro3[((j + ldn * 3) << 1) - 7]
3278  ) * mix4_1.al2ro3[j - 1];
3279 L300:
3280  ;
3281  }
3282 L3000:
3283  if (mix4_1.n2ro4 == 0) {
3284  goto L4000;
3285  }
3286  i__1 = mix4_1.n2ro4;
3287  for (j = 1; j <= i__1; ++j) {
3288  s1sum += f0c_1.f[*i__ - 1] * (mix4_1.q2ro4[((j + *i__ * 3) << 1) - 8] +
3289  mix4_1.q2ro4[((j + *i__ * 3) << 1) - 7]);
3290  lup = *i__ + mix4_1.l2ro4[j - 1];
3291  if (lup >= *nstep1) {
3292  goto L350;
3293  }
3294  *s2sum = *s2sum + f0c_1.f[lup - 1] * mix4_1.q2ro4[((j + lup * 3) << 1)
3295  - 8] + (f0c_1.f[lup] * mix4_1.q2ro4[((j + (lup + 1) * 3) << 1)
3296  - 8] - f0c_1.f[lup - 1] * mix4_1.q2ro4[((j + lup * 3) << 1) - 8]
3297  ) * mix4_1.al2ro4[j - 1];
3298 L350:
3299  ;
3300  }
3301  i__1 = mix4_1.n2ro4;
3302  for (j = 2; j <= i__1; ++j) {
3303  ldn = *i__ - mix4_1.l2ro4[j - 2];
3304  if (ldn < 1) {
3305  goto L400;
3306  }
3307  *s2sum = *s2sum + f0c_1.f[ldn - 1] * mix4_1.q2ro4[((j + ldn * 3) << 1)
3308  - 7] + (f0c_1.f[ldn] * mix4_1.q2ro4[((j + (ldn + 1) * 3) << 1)
3309  - 7] - f0c_1.f[ldn - 1] * mix4_1.q2ro4[((j + ldn * 3) << 1) - 7]
3310  ) * mix4_1.al2ro4[j - 1];
3311 L400:
3312  ;
3313  }
3314 L4000:
3315  *s2sum -= s1sum;
3316  return 0;
3317 } /* type2_ */
3318 
3319 int StFtpcMagboltz1::type2g_(double *s2sum, int *i__, int *nstep1)
3320 {
3321  /* System generated locals */
3322  int i__1;
3323 
3324  /* Local variables */
3325  static double s1sum;
3326  static int j, ldn, lup;
3327 
3328 
3329  s1sum = (float)0.;
3330  *s2sum = (float)0.;
3331  if (mix4_1.n2ro1 == 0) {
3332  goto L1000;
3333  }
3334  i__1 = mix4_1.n2ro1;
3335  for (j = 1; j <= i__1; ++j) {
3336  s1sum += g0c_1.g[*i__ - 1] * (mix4_1.q2ro1[((j + *i__ * 3) << 1) - 8] +
3337  mix4_1.q2ro1[((j + *i__ * 3) << 1) - 7]);
3338  lup = *i__ + mix4_1.l2ro1[j - 1];
3339  if (lup >= *nstep1) {
3340  goto L50;
3341  }
3342  *s2sum = *s2sum + g0c_1.g[lup - 1] * mix4_1.q2ro1[((j + lup * 3) << 1)
3343  - 8] + (g0c_1.g[lup] * mix4_1.q2ro1[((j + (lup + 1) * 3) << 1)
3344  - 8] - g0c_1.g[lup - 1] * mix4_1.q2ro1[((j + lup * 3) << 1) - 8]
3345  ) * mix4_1.al2ro1[j - 1];
3346 L50:
3347  ;
3348  }
3349  i__1 = mix4_1.n2ro1;
3350  for (j = 2; j <= i__1; ++j) {
3351  ldn = *i__ - mix4_1.l2ro1[j - 2];
3352  if (ldn < 1) {
3353  goto L100;
3354  }
3355  *s2sum = *s2sum + g0c_1.g[ldn - 1] * mix4_1.q2ro1[((j + ldn * 3) << 1)
3356  - 7] + (g0c_1.g[ldn] * mix4_1.q2ro1[((j + (ldn + 1) * 3) << 1)
3357  - 7] - g0c_1.g[ldn - 1] * mix4_1.q2ro1[((j + ldn * 3) << 1) - 7]
3358  ) * mix4_1.al2ro1[j - 1];
3359 L100:
3360  ;
3361  }
3362 L1000:
3363  if (mix4_1.n2ro2 == 0) {
3364  goto L2001;
3365  }
3366  i__1 = mix4_1.n2ro2;
3367  for (j = 1; j <= i__1; ++j) {
3368  s1sum += g0c_1.g[*i__ - 1] * (mix4_1.q2ro2[((j + *i__ * 3) << 1) - 8] +
3369  mix4_1.q2ro2[((j + *i__ * 3) << 1) - 7]);
3370  lup = *i__ + mix4_1.l2ro2[j - 1];
3371  if (lup >= *nstep1) {
3372  goto L150;
3373  }
3374  *s2sum = *s2sum + g0c_1.g[lup - 1] * mix4_1.q2ro2[((j + lup * 3) << 1)
3375  - 8] + (g0c_1.g[lup] * mix4_1.q2ro2[((j + (lup + 1) * 3) << 1)
3376  - 8] - g0c_1.g[lup - 1] * mix4_1.q2ro2[((j + lup * 3) << 1) - 8]
3377  ) * mix4_1.al2ro2[j - 1];
3378 L150:
3379  ;
3380  }
3381  i__1 = mix4_1.n2ro2;
3382  for (j = 2; j <= i__1; ++j) {
3383  ldn = *i__ - mix4_1.l2ro2[j - 2];
3384  if (ldn < 1) {
3385  goto L200;
3386  }
3387  *s2sum = *s2sum + g0c_1.g[ldn - 1] * mix4_1.q2ro2[((j + ldn * 3) << 1)
3388  - 7] + (g0c_1.g[ldn] * mix4_1.q2ro2[((j + (ldn + 1) * 3 )<< 1)
3389  - 7] - g0c_1.g[ldn - 1] * mix4_1.q2ro2[((j + ldn * 3) << 1) - 7]
3390  ) * mix4_1.al2ro2[j - 1];
3391 L200:
3392  ;
3393  }
3394 L2001:
3395  if (mix4_1.n2ro3 == 0) {
3396  goto L3000;
3397  }
3398  i__1 = mix4_1.n2ro3;
3399  for (j = 1; j <= i__1; ++j) {
3400  s1sum += g0c_1.g[*i__ - 1] * (mix4_1.q2ro3[((j + *i__ * 3) << 1) - 8] +
3401  mix4_1.q2ro3[((j + *i__ * 3) << 1) - 7]);
3402  lup = *i__ + mix4_1.l2ro3[j - 1];
3403  if (lup >= *nstep1) {
3404  goto L250;
3405  }
3406  *s2sum = *s2sum + g0c_1.g[lup - 1] * mix4_1.q2ro3[((j + lup * 3) << 1)
3407  - 8] + (g0c_1.g[lup] * mix4_1.q2ro3[((j + (lup + 1) * 3) << 1)
3408  - 8] - g0c_1.g[lup - 1] * mix4_1.q2ro3[((j + lup * 3) << 1) - 8]
3409  ) * mix4_1.al2ro3[j - 1];
3410 L250:
3411  ;
3412  }
3413  i__1 = mix4_1.n2ro3;
3414  for (j = 2; j <= i__1; ++j) {
3415  ldn = *i__ - mix4_1.l2ro3[j - 2];
3416  if (ldn < 1) {
3417  goto L300;
3418  }
3419  *s2sum = *s2sum + g0c_1.g[ldn - 1] * mix4_1.q2ro3[((j + ldn * 3) << 1)
3420  - 7] + (g0c_1.g[ldn] * mix4_1.q2ro3[((j + (ldn + 1) * 3) << 1)
3421  - 7] - g0c_1.g[ldn - 1] * mix4_1.q2ro3[((j + ldn * 3) << 1) - 7]
3422  ) * mix4_1.al2ro3[j - 1];
3423 L300:
3424  ;
3425  }
3426 L3000:
3427  if (mix4_1.n2ro4 == 0) {
3428  goto L4000;
3429  }
3430  i__1 = mix4_1.n2ro4;
3431  for (j = 1; j <= i__1; ++j) {
3432  s1sum += g0c_1.g[*i__ - 1] * (mix4_1.q2ro4[((j + *i__ * 3) << 1) - 8] +
3433  mix4_1.q2ro4[((j + *i__ * 3) << 1) - 7]);
3434  lup = *i__ + mix4_1.l2ro4[j - 1];
3435  if (lup >= *nstep1) {
3436  goto L350;
3437  }
3438  *s2sum = *s2sum + g0c_1.g[lup - 1] * mix4_1.q2ro4[((j + lup * 3) << 1)
3439  - 8] + (g0c_1.g[lup] * mix4_1.q2ro4[((j + (lup + 1) * 3) << 1)
3440  - 8] - g0c_1.g[lup - 1] * mix4_1.q2ro4[((j + lup * 3) << 1) - 8]
3441  ) * mix4_1.al2ro4[j - 1];
3442 L350:
3443  ;
3444  }
3445  i__1 = mix4_1.n2ro4;
3446  for (j = 2; j <= i__1; ++j) {
3447  ldn = *i__ - mix4_1.l2ro4[j - 2];
3448  if (ldn < 1) {
3449  goto L400;
3450  }
3451  *s2sum = *s2sum + g0c_1.g[ldn - 1] * mix4_1.q2ro4[((j + ldn * 3) << 1)
3452  - 7] + (g0c_1.g[ldn] * mix4_1.q2ro4[((j + (ldn + 1) * 3) << 1)
3453  - 7] - g0c_1.g[ldn - 1] * mix4_1.q2ro4[((j + ldn * 3) << 1) - 7]
3454  ) * mix4_1.al2ro4[j - 1];
3455 L400:
3456  ;
3457  }
3458 L4000:
3459  *s2sum -= s1sum;
3460  return 0;
3461 } /* type2g_ */
3462 
3463 
3464 
3465 
3466 
3467 
3468 
3469 
3470 
3471 
3472 
3473 
3474 
3475 
3476 
Generated by   doxygen 1.8.5