42 #if defined(PRINT_DR) || defined(PRINT_RR)
43 static const char FILE_NAME_OUT[] =
"array.out";
61 int nAtomicNumberCScale,
63 int n_core_e_before_recomb )
66 double RateCoefficient, sum;
70 ASSERT( nAtomicNumberCScale>=0 && nAtomicNumberCScale<
LIMELM );
72 if( nAtomicNumberCScale==
ipIRON && n_core_e_before_recomb>=12 &&
73 n_core_e_before_recomb<=18 )
84 {5.636e-4, 7.390e-3, 3.635e-2, 1.693e-1, 3.315e-2, 2.288e-1, 7.316e-2, 0.},
85 {1.090e-3, 7.801e-3, 1.132e-2, 4.740e-2, 1.990e-1, 3.379e-2, 1.140e-1, 1.250e-1},
86 {3.266e-3, 7.637e-3, 1.005e-2, 2.527e-2, 6.389e-2, 1.564e-1, 0., 0.},
87 {1.074e-3, 6.080e-3, 1.887e-2, 2.540e-2, 7.580e-2, 2.773e-1, 0., 0.},
88 {9.073e-4, 3.777e-3, 1.027e-2, 3.321e-2, 8.529e-2, 2.778e-1, 0., 0.},
89 {5.335e-4, 1.827e-3, 4.851e-3, 2.710e-2, 8.226e-2, 3.147e-1, 0., 0.},
90 {7.421e-4, 2.526e-3, 4.605e-3, 1.489e-2, 5.891e-2, 2.318e-1, 0., 0.}
96 {3.628e3, 2.432e4, 1.226e5, 4.351e5, 1.411e6, 6.589e6, 1.030e7, 0},
97 {1.246e3, 1.063e4, 4.719e4, 1.952e5, 5.637e5, 2.248e6, 7.202e6, 3.999e9},
98 {1.242e3, 1.001e4, 4.466e4, 1.497e5, 3.919e5, 6.853e5, 0. , 0.},
99 {1.387e3, 1.048e4, 3.955e4, 1.461e5, 4.010e5, 7.208e5, 0. , 0.},
100 {1.525e3, 1.071e4, 4.033e4, 1.564e5, 4.196e5, 7.580e5, 0. , 0.},
101 {2.032e3, 1.018e4, 4.638e4, 1.698e5, 4.499e5, 7.880e5, 0. , 0.},
102 {3.468e3, 1.353e4, 3.690e4, 1.957e5, 4.630e5, 8.202e5, 0. , 0.}
105 long int nion = n_core_e_before_recomb - 12;
106 ASSERT( nion>=0 && nion <=6 );
113 sum += (cFe_q[nion][i] *
sexp( EFe_q[nion][i]/
phycon.
te));
119 return RateCoefficient;
123 else if( nAtomicNumberCScale < n_core_e_before_recomb )
125 RateCoefficient = -2;
128 else if( nAtomicNumberCScale >=
LIMELM )
130 RateCoefficient = -2;
135 RateCoefficient = -1;
143 for(i=0; i<
nDRFitPar[nAtomicNumberCScale][n_core_e_before_recomb]; ++i )
145 sum += (
DRFitParPart1[nAtomicNumberCScale][n_core_e_before_recomb][i] *
155 RateCoefficient = -99;
158 return RateCoefficient;
167 int nAtomicNumberCScale,
169 int n_core_e_before_recomb )
171 double RateCoefficient;
176 ASSERT( nAtomicNumberCScale>=0 && nAtomicNumberCScale<
LIMELM );
178 if( nAtomicNumberCScale==
ipIRON &&
179 n_core_e_before_recomb>=12 && n_core_e_before_recomb<=18 )
188 double parFeq[7][6] ={
189 {1.179e-9 , 0.7096, 4.508e2, 3.393e7, 0.0154, 3.977e6},
190 {1.050e-9 , 0.6939, 4.568e2, 3.987e7, 0.0066, 5.451e5},
191 {9.832e-10, 0.7146, 3.597e2, 3.808e7, 0.0045, 3.952e5},
192 {8.303e-10, 0.7156, 3.531e2, 3.554e7, 0.0132, 2.951e5},
193 {1.052e-9 , 0.7370, 1.639e2, 2.924e7, 0.0224, 4.291e5},
194 {1.338e-9 , 0.7495, 7.242e1, 2.453e7, 0.0404, 4.199e5},
195 {1.263e-9 , 0.7532, 5.209e1, 2.169e7, 0.0421, 2.917e5}
200 long int nion = n_core_e_before_recomb - 12;
201 ASSERT( nion>=0 && nion <=6 );
204 B = parFeq[nion][1] + parFeq[nion][4]*exp(temp);
205 D = sqrt(
phycon.
te/parFeq[nion][2]);
206 F = sqrt(
phycon.
te/parFeq[nion][3]);
207 RateCoefficient = parFeq[nion][0]/(D*pow((1.+D),(1.-B))*pow((1.+F),(1.+B)));
209 return RateCoefficient;
213 else if( nAtomicNumberCScale < n_core_e_before_recomb )
215 RateCoefficient = -2;
218 else if( nAtomicNumberCScale >=
LIMELM )
220 RateCoefficient = -2;
225 RateCoefficient = -1;
236 B =
RRFitPar[nAtomicNumberCScale][n_core_e_before_recomb][1] +
237 RRFitPar[nAtomicNumberCScale][n_core_e_before_recomb][4]*exp(temp);
238 D = sqrt(
phycon.
te/
RRFitPar[nAtomicNumberCScale][n_core_e_before_recomb][2]);
239 F = sqrt(
phycon.
te/
RRFitPar[nAtomicNumberCScale][n_core_e_before_recomb][3]);
240 RateCoefficient =
RRFitPar[nAtomicNumberCScale][n_core_e_before_recomb][0]/(D*pow((1.+D),(1.-B))*pow((1.+F),(1.+B)));
245 RateCoefficient = -99;
247 return RateCoefficient;
260 int NuclearCharge=-1, NumberElectrons=-1;
264 int M_state, W_state,
267 const int NBLOCK = 2;
268 int data_begin_line[NBLOCK];
271 const char* chFilename;
275 #define BIGGEST_INDEX_TO_USE 103
283 long INDX=0,INDP=0,
N=0,
S=0,L=0,J=0,maxINDX=0,loopindex=0,max_N_of_data=-1;
286 static int nCalled = 0;
288 const char* cdDATAFILE[] =
292 "nrb00#h_he1ic12.dat",
293 "nrb00#h_li2ic12.dat",
294 "nrb00#h_be3ic12.dat",
295 "nrb00#h_b4ic12.dat",
296 "nrb00#h_c5ic12.dat",
297 "nrb00#h_n6ic12.dat",
298 "nrb00#h_o7ic12.dat",
299 "nrb00#h_f8ic12.dat",
300 "nrb00#h_ne9ic12.dat",
301 "nrb00#h_na10ic12.dat",
302 "nrb00#h_mg11ic12.dat",
303 "nrb00#h_al12ic12.dat",
304 "nrb00#h_si13ic12.dat",
305 "nrb00#h_p14ic12.dat",
306 "nrb00#h_s15ic12.dat",
307 "nrb00#h_cl16ic12.dat",
308 "nrb00#h_ar17ic12.dat",
309 "nrb00#h_k18ic12.dat",
310 "nrb00#h_ca19ic12.dat",
311 "nrb00#h_sc20ic12.dat",
312 "nrb00#h_ti21ic12.dat",
313 "nrb00#h_v22ic12.dat",
314 "nrb00#h_cr23ic12.dat",
315 "nrb00#h_mn24ic12.dat",
316 "nrb00#h_fe25ic12.dat",
317 "nrb00#h_co26ic12.dat",
318 "nrb00#h_ni27ic12.dat",
319 "nrb00#h_cu28ic12.dat",
320 "nrb00#h_zn29ic12.dat"
333 # if defined(PRINT_DR) || defined(PRINT_RR)
343 ioDATA=
open_data( cdDATAFILE[nelem],
"r" );
349 TheirIndexToOurIndex[i] = -1;
356 if(
nMatch(
"INDX INDP ",
string) )
361 fprintf(
ioQQQ,
" Badnell data file appears to be corrupted.\n");
368 if( strcmp(
string,
"\n")==0 )
376 if( INDX >= BIGGEST_INDEX_TO_USE )
383 ASSERT( INDX < BIGGEST_INDEX_TO_USE );
390 if( (i1=
nMatch(
"1S1 ",
string)) > 0 )
401 if( (i1=
nMatch(
" (",
string)) > 0 )
421 ASSERT( J <= ( L + (
int)((
S+1)/2) ) &&
422 J >= ( L - (
int)((
S+1)/2) ) && J >= 0 );
437 if(
N==2 && L==1 &&
S==3 )
440 TheirIndexToOurIndex[INDX] = 3;
442 TheirIndexToOurIndex[INDX] = 4;
446 ASSERT( TheirIndexToOurIndex[INDX] == 5 );
449 max_N_of_data =
MAX2( max_N_of_data,
N );
468 ASSERT( maxINDX < BIGGEST_INDEX_TO_USE );
477 if(
nMatch(
"INDX TE= ",
string) )
484 fprintf(
ioQQQ,
" Badnell data file appears to be corrupted.\n");
492 if(
nMatch(
"PRTF",
string) || INDX >= maxINDX || INDX<0 )
500 for(loopindex=0;loopindex<10;loopindex++)
504 TheirIndexToOurIndex[INDX] > 0 )
513 fprintf(
ioQQQ,
" Badnell data file appears to be corrupted.\n");
519 for(loopindex=10;loopindex<19;loopindex++)
523 TheirIndexToOurIndex[INDX] > 0 )
536 ASSERT( maxINDX < BIGGEST_INDEX_TO_USE );
537 ASSERT( max_N_of_data > 0 );
546 for( i=TheirIndexToOurIndex[maxINDX]+1;
556 for(loopindex=0;loopindex<19;loopindex++)
572 for(loopindex=0;loopindex<19;loopindex++)
584 for( i=0; i<NBLOCK; ++i )
587 data_begin_line[i] = INT_MIN;
590 chFilename =
"badnell_dr.dat";
606 data_begin_line[number] = count;
607 ASSERT( number < NBLOCK );
626 for( nelem=0; nelem<
LIMELM; nelem++ )
639 for(
long ion=0; ion<nelem+1; ++ion )
667 fseek(ioDATA, 0, SEEK_SET);
671 fprintf(
ioQQQ,
" DISASTER PROBLEM Badnell_rec_init could not read first line of badnell_dr.dat.\n");
677 if( (chs = strchr(chLine,
')'))==NULL )
680 fprintf(
ioQQQ,
" DISASTER PROBLEM Badnell_rec_init data file incorrect format.\n");
685 sscanf(chs,
"%4i%2i%2i",&yr, &mo, &dy);
687 int dr_yr = 2007, dr_mo = 10, dr_dy = 27;
688 if((yr != dr_yr) || (mo != dr_mo) || (dy != dr_dy))
691 "DISASTER PROBLEM Badnell_rec_init The version of %s I found (%i %i %i) is not the current version (%i %i %i).\n",
692 chFilename, yr, mo, dy, dr_yr, dr_mo, dr_dy);
693 fprintf(
ioQQQ,
" The first line of the file is the following\n %s\n", chLine );
700 length_of_line = (int)strlen(chLine);
703 if( count > data_begin_line[0] && count < data_begin_line[1] && length_of_line >3 )
710 sscanf(chLine,
"%i%i%i%i%lf%lf%lf%lf%lf%lf%lf%lf%lf",
711 &NuclearCharge, &NumberElectrons, &M_state, &W_state, &par_C[0], &par_C[1], &par_C[2],
712 &par_C[3], &par_C[4], &par_C[5], &par_C[6], &par_C[7], &par_C[8]);
715 long int NuclearChargeM1 = NuclearCharge-1;
717 if(M_state == 1 && NuclearChargeM1 < LIMELM )
720 ASSERT( NumberElectrons < LIMELM );
721 ASSERT( NuclearChargeM1 < LIMELM );
725 nDRFitPar[NuclearChargeM1][NumberElectrons] = 9;
726 for( i=8; i>=0; i-- )
729 --
nDRFitPar[NuclearChargeM1][NumberElectrons];
736 DRFitParPart1[NuclearChargeM1][NumberElectrons][i] = par_C[i];
742 fseek(ioDATA, 0, SEEK_SET);
747 length_of_line = (int)strlen(chLine);
748 if( count > data_begin_line[1] && length_of_line > 3 )
756 sscanf(chLine,
"%i%i%i%i%lf%lf%lf%lf%lf%lf%lf%lf%lf",
757 &NuclearCharge, &NumberElectrons, &M_state, &W_state, &par_E[0], &par_E[1], &par_E[2],
758 &par_E[3], &par_E[4], &par_E[5], &par_E[6], &par_E[7], &par_E[8]);
760 long int NuclearChargeM1 = NuclearCharge-1;
762 if(M_state == 1 && NuclearChargeM1<LIMELM)
764 ASSERT( NumberElectrons < LIMELM );
765 ASSERT( NuclearChargeM1 < LIMELM );
769 nDRFitPar[NuclearChargeM1][NumberElectrons] = 9;
770 for( i=8; i>=0; i-- )
773 --
nDRFitPar[NuclearChargeM1][NumberElectrons];
779 for( i=0; i<
nDRFitPar[NuclearChargeM1][NumberElectrons]; i++ )
780 DRFitParPart2[NuclearChargeM1][NumberElectrons][i] = par_E[i];
789 for( nelem=0; nelem<
LIMELM; nelem++ )
791 for(
int ion=0; ion<nelem+1;++ion )
795 fprintf(ofp,
"%i %i %e %e %e %e %e %e %e %e %e\n",
804 for( nelem=0; nelem<
LIMELM; nelem++ )
806 for(
int ion=0; ion<nelem+1; ion++ )
810 fprintf(ofp,
"%i %i %e %e %e %e %e %e %e %e %e\n",
824 bool lgDRBadnellBothDefined =
true;
825 for( nelem=0; nelem<
LIMELM; nelem++ )
827 for(
int ion=0; ion<nelem+1; ion++ )
833 fprintf(
ioQQQ,
"DR %i, RR %i: %c %c\n", nelem, ion,
836 fprintf(
ioQQQ,
"PROBLEM ion_recomb_Badnell first and second half of Badnell DR not consistent.\n");
837 lgDRBadnellBothDefined =
false;
842 if( !lgDRBadnellBothDefined )
846 "DISASTER PROBLEM The DR data files are corrupted - part 1 and 2 do not agree.\n");
847 fprintf(
ioQQQ,
" Start again with a fresh copy of the data directory\n" );
852 chFilename =
"badnell_rr.dat";
859 fprintf(
ioQQQ,
" DISASTER PROBLEM Badnell_rec_init could not read first line of badnell_rr.dat.\n");
863 if( (chs = strchr(chLine,
')'))==NULL )
866 fprintf(
ioQQQ,
" DISASTER PROBLEM Badnell_rec_init data file incorrect format.\n");
870 sscanf(chs,
"%4i%2i%2i", &yr, &mo, &dy);
871 int rr_yr = 2007, rr_mo = 10, rr_dy = 27;
872 if((yr != rr_yr)||(mo != rr_mo)||(dy != rr_dy))
874 fprintf(
ioQQQ,
"DISASTER PROBLEM The version of %s I found (%i %i %i) is not the current version (%i %i %i).\n",
875 chFilename, yr, mo, dy, rr_yr, rr_mo, rr_dy);
876 fprintf(
ioQQQ,
" The line was as follows:\n %s\n", chLine );
890 sscanf(chLine,
"%i%i%i%i%lf%lf%lf%lf%lf%lf",
891 &NuclearCharge, &NumberElectrons, &M_state, &W_state, &temp_par[0], &temp_par[1],
892 &temp_par[2], &temp_par[3], &temp_par[4], &temp_par[5]);
893 long NuclearChargeM1 = NuclearCharge-1;
895 if(M_state == 1 && NuclearChargeM1<LIMELM)
897 ASSERT( NuclearChargeM1 < LIMELM );
898 ASSERT( NumberElectrons <= LIMELM );
903 RRFitPar[NuclearChargeM1][NumberElectrons][i] = temp_par[i];
911 for( nelem=0; nelem<
LIMELM; nelem++ )
913 for( ion=0; ion<nelem+1; ion++ )
917 fprintf(ofp,
"%i %i %e %e %e %e %e %e\n",
918 nelem, ion,
RRFitPar[nelem][ion][0],
926 fprintf(ofp,
"total lines are %i ", count);
934 enum {DEBUG_LOC=
false};
940 fprintf(
ioQQQ,
"\nDEBUG rr rec\t%i",nelem);
941 for( ion=0; ion<=nelem; ++ion )
946 fprintf(
ioQQQ,
"DEBUG dr rec\t%i",nelem);
947 for( ion=0; ion<=nelem; ++ion )
962 static double TeUsed = -1;
970 if( fabs(
phycon.
te/TeUsed - 1. ) < 0.001 )
976 for( ion=0; ion<
LIMELM; ++ion )
987 for( ion=0; ion<nelem+1; ++ion )
989 long int n_bnd_elec_before_recom ,
990 n_bnd_elec_after_recom;
992 n_bnd_elec_before_recom = nelem-ion;
993 n_bnd_elec_after_recom = nelem-ion+1;
994 # define DR2SMALL 1e-15
1002 n_bnd_elec_before_recom )) >= 0. )
1025 n_bnd_elec_before_recom )) >= 0. )
1043 n_bnd_elec_after_recom ,
1049 double Fe_Gu_c[9][6] = {
1050 { 2.50507e-11, 5.60226e-11, 1.85001e-10, 3.57495e-9, 1.66321e-7, 0. },
1051 { 9.19610e-11, 2.92460e-10, 1.02120e-9, 1.14852e-8, 3.25418e-7, 0. },
1052 { 9.02625e-11, 6.22962e-10, 5.77545e-9, 1.78847e-8, 3.40610e-7, 0. },
1053 { 9.04286e-12, 9.68148e-10, 4.83636e-9, 2.48159e-8, 3.96815e-7, 0. },
1054 { 6.77873e-10, 1.47252e-9, 5.31202e-9, 2.54793e-8, 3.47407e-7, 0. },
1055 { 1.29742e-9, 4.10172e-9, 1.23605e-8, 2.33615e-8, 2.97261e-7, 0. },
1056 { 8.78027e-10, 2.31680e-9, 3.49333e-9, 1.16927e-8, 8.18537e-8, 1.54740e-7 },
1057 { 2.23178e-10, 1.87313e-9, 2.86171e-9, 1.38575e-8, 1.17803e-7, 1.06251e-7 },
1058 { 2.17263e-10, 7.35929e-10, 2.81276e-9, 1.32411e-8, 1.15761e-7, 4.80389e-8 }
1062 { 8.30501e-2, 8.52897e-1, 3.40225e0, 2.23053e1, 6.80367e1, 0. },
1063 { 1.44392e-1, 9.23999e-1, 5.45498e0, 2.04301e1, 7.06112e1, 0. },
1064 { 5.79132e-2, 1.27852e0, 3.22439e0, 1.79602e1, 6.96277e1, 0. },
1065 { 1.02421e-1, 1.79393e0, 4.83226e0, 1.91117e1, 6.80858e1, 0. },
1066 { 1.24630e-1, 6.86045e-1, 3.09611e0, 1.44023e1, 6.42820e1, 0. },
1067 { 1.34459e-1, 6.63028e-1, 2.61753e0, 1.30392e1, 6.10222e1, 0. },
1068 { 7.79748e-2, 5.35522e-1, 1.88407e0, 8.38459e0, 3.38613e1, 7.89706e1 },
1069 { 8.83019e-2, 6.12756e-1, 2.36035e0, 9.61736e0, 3.64467e1, 8.72406e1 },
1070 { 1.51322e-1, 5.63155e-1, 2.57013e0, 9.08166e0, 3.69528e1, 1.08067e2 }
1073 double s_c[5][2] = {
1074 {0.1565e-3, 0.1617e-2},
1075 {0.3026e-3, 0.2076e-1},
1076 {0.3177e-1, 0.6309e-3},
1077 {0.2464e-1, 0.5553e-3},
1078 {0.1924e-1, 0.5111e-3}
1082 {0.1157e6, 0.1868e6},
1083 {0.9662e5, 0.1998e6},
1084 {0.1928e6, 0.6126e5},
1085 {0.1806e6, 0.3142e5},
1086 {0.1519e6, 0.4978e5}
1101 double fe14_c[6]={7.07E-04,7.18E-03,2.67E-02,3.15E-02,1.62E-01,5.37E-04};
1105 double fe14_E[6]={4.12E-01,2.06E+00,1.03E+01,2.20E+01,4.22E+01,3.41E+03};
1108 for( i=0; i<6; i++ )
1126 double fe13_c[10]={ 3.55e-4, 2.40e-3, 7.83e-3, 1.10e-2, 3.30e-2, 1.45e-1, 8.50e-2, 2.59e-2, 8.93e-3, 9.80e-3 },
1127 fe13_E[10]={ 2.19e-2, 1.79e-1, 7.53e-1, 2.21e0, 9.57e0, 3.09e1, 6.37e1, 2.19e2, 1.50e3, 7.86e3 };
1132 for( i=0; i<10; i++ )
1155 for( ion=0; ion<9; ion++ )
1161 for( i=0; i<6; i++ )
1178 for( ion=0; ion<
LIMELM; ++ion )
1190 for( ion=0; ion<5; ion++ )
1200 for( i=0; i<2; i++ )
1240 fprintf(
ioQQQ,
"DEBUG Badnell recombination RR, then DR, T=%.3e\n",
phycon.
te );
1243 fprintf(
ioQQQ,
"nelem=%li %s, RR then DR\n",
1245 for( ion=0; ion<nelem+1; ++ion )
1253 fprintf(
ioQQQ,
" %.2e", -1. );
1256 fprintf(
ioQQQ,
"\n" );
1257 for( ion=0; ion<nelem+1; ++ion )
1265 fprintf(
ioQQQ,
" %.2e", -1. );
1268 fprintf(
ioQQQ,
"\n\n" );
1271 fprintf(
ioQQQ,
"mean DR recombination ion mean stddev stddev/mean \n" );
1272 for( ion=0; ion<
LIMELM; ++ion )
1276 for( nelem=ion; nelem<
LIMELM; ++nelem )
1282 stddev = sqrt( stddev /
MAX2( 1 ,
nsumrec[ion]-1 ) );
1283 fprintf(
ioQQQ,
" %2li %.2e %.2e %.2e\n",