37 STATIC double TempInterp(
double* TempArray ,
double* ValueArray,
long NumElements );
43 DEBUG_ENTRY(
"iso_radrecomb_from_cross_section()" );
62 long ipFirstCollapsed, LastN=0L, ThisN=1L, ipLevel;
63 double topoff, TotMinusExplicitResolved,
64 TotRRThisN=0., TotRRLastN=0., Total_DR_Added=0.;
65 double RecExplictLevels, TotalRadRecomb, RecCollapsed;
67 {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
68 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
69 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},
70 {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
71 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
72 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.}};
80 RecExplictLevels = 0.;
89 for( ipLevel=0; ipLevel<ipFirstCollapsed; ++ipLevel )
169 if(
N_(ipLevel) == ThisN )
178 TotRRLastN = TotRRThisN;
184 enum {DEBUG_LOC=
false};
186 static long RUNONCE =
false;
188 if( !RUNONCE && DEBUG_LOC )
190 static long FIRSTTIME =
true;
194 fprintf(
ioQQQ,
"Sum of Radiative Recombination at current iso, nelem, temp = %li %li %.2f\n",
199 fprintf(
ioQQQ,
"%li\t%.2e\n",LastN,TotRRLastN );
210 TotalRadRecomb = RecCollapsed + RecExplictLevels;
252 TotMinusExplicitResolved = TotalRadRecomb - RecExplictLevels;
254 topoff = TotMinusExplicitResolved - RecCollapsed;
260 fabs( topoff/TotalRadRecomb ) > 0.01 )
262 fprintf(
ioQQQ,
" PROBLEM negative RR topoff for %li, rel error was %.2e, temp was %.2f\n",
263 nelem, topoff/TotalRadRecomb,
phycon.
te );
273 if( Total_DR_Added > TotalRadRecomb/100. )
277 fprintf(
ioQQQ,
" PROBLEM negative DR topoff for %li, rel error was %.2e, temp was %.2f\n",
344 fprintf(
ioQQQ,
" iso_radiative_recomb trace ipISO=%3ld Z=%3ld\n", ipISO, nelem );
347 fprintf(
ioQQQ,
" iso_radiative_recomb recomb effic" );
352 fprintf(
ioQQQ,
"\n" );
355 fprintf(
ioQQQ,
" iso_radiative_recomb recomb net effic" );
360 fprintf(
ioQQQ,
"\n" );
363 fprintf(
ioQQQ,
" iso_radiative_recomb in optic dep" );
368 fprintf(
ioQQQ,
"\n" );
371 fprintf(
ioQQQ,
" iso_radiative_recomb out op depth" );
376 fprintf(
ioQQQ,
"\n" );
379 fprintf(
ioQQQ,
" iso_radiative_recomb rad rec coef " );
384 fprintf(
ioQQQ,
"\n" );
389 fprintf(
ioQQQ,
" iso_radiative_recomb total rec coef" );
391 fprintf(
ioQQQ,
" case A=" );
394 fprintf(
ioQQQ,
" case B=");
396 fprintf(
ioQQQ,
"\n" );
409 " PROBLEM iso_radiative_recomb non-positive recombination coefficient for ipISO=%3ld Z=%3ld lev n=%3ld rec=%11.2e te=%11.2e\n",
464 enum {DEBUG_LOC=
false};
470 fprintf(
ioQQQ,
"Effective recombination, ipISO=%li, nelem=%li, Te = %e\n", ipISO, nelem,
phycon.
te );
471 fprintf(
ioQQQ,
"N\tL\tS\tRadEffec\tLifetime\n" );
473 for(
long i=0; i<maxPrt; i++ )
475 fprintf(
ioQQQ,
"%li\t%li\t%li\t%e\t%e\n",
N_(i),
L_(i),
S_(i),
479 fprintf(
ioQQQ,
"\n" );
486 dprintf(
ioQQQ,
"ipHi\tipLo\tWL\tEmiss\tSigmaEmiss\tRadEffec\tSigRadEff\tBrRat\tSigBrRat\n" );
508 for(
long ipLo = 0; ipLo < ipHi; ipLo++ )
510 if( ((
L_(ipLo) ==
L_(ipHi) + 1 ) || (
L_(ipHi) ==
L_(ipLo) + 1 )) )
515 double sigma_emiss = 0., SigmaBranchRatio = 0.;
517 if( ( emissivity > 2.E-29 ) && ( wavelength < 1.E6 ) && (
N_(ipHi)<=5) )
519 SigmaBranchRatio =
iso.
BranchRatio[ipISO][nelem][ipHi][ipLo] * sqrt(
520 pow( (
double)
iso.
Error[ipISO][nelem][ipHi][ipLo][
IPRAD], 2. ) +
523 sigma_emiss =
EN1RYD * EnergyInRydbergs * sqrt(
525 pow( SigmaBranchRatio *
iso.
RadEffec[ipISO][nelem][ipHi], 2.) );
530 fprintf(
ioQQQ,
"\t%e\t%e\t%e\t%e\t%e\t%e\n",
565 rate = pow( 10. , rate );
574 double RecombRelError ,
595 RecombRelError = ( RecombInterp - RecombCalc ) /
MAX2( RecombInterp , RecombCalc );
597 return RecombRelError;
609 for(
long ipISO=0; ipISO<
NISO; ipISO++ )
616 for(
long nelem=ipISO; nelem <
LIMELM; ++nelem )
618 long int MaxLevels, maxN;
638 RRCoef[ipISO][nelem] = (
double**)
MALLOC(
sizeof(
double*)*(unsigned)(MaxLevels) );
640 for(
long ipLo=0; ipLo < MaxLevels;++ipLo )
642 RRCoef[ipISO][nelem][ipLo] = (
double*)
MALLOC(
sizeof(
double)*(unsigned)N_ISO_TE_RECOMB );
656 TeRRCoef[N_ISO_TE_RECOMB-1] += 0.01f;
665 for(
long i = 0; i<
NISO; i++ )
676 double RadRecombReturn;
677 long int i, i1, i2, i3, i4, i5;
678 long int ipLo, nelem;
681 char chLine[chLine_LENGTH];
683 const char* chFilename[
NISO] =
684 {
"h_iso_recomb.dat",
"he_iso_recomb.dat" };
706 fprintf(
ioQQQ,
" iso_recomb_setup opening %s:", chFilename[ipISO] );
711 ioDATA =
open_data( chFilename[ipISO],
"r" );
715 fprintf(
ioQQQ,
" Defaulting to on-the-fly computation, ");
716 fprintf(
ioQQQ,
" but the code runs much faster if you compile %s!\n", chFilename[ipISO]);
722 for( nelem = ipISO; nelem <
LIMELM; nelem++ )
743 RRCoef[ipISO][nelem][ipLo][i] = log10(RadRecombReturn);
767 fprintf(
ioQQQ,
" iso_recomb_setup could not read first line of %s.\n", chFilename[ipISO]);
778 " iso_recomb_setup: the version of %s is not the current version.\n", chFilename[ipISO] );
780 " iso_recomb_setup: I expected to find the numbers %i %li %li %i and got %li %li %li %li instead.\n" ,
786 fprintf(
ioQQQ,
"Here is the line image:\n==%s==\n", chLine );
788 " iso_recomb_setup: please recompile the data file with the COMPile RECOmb COEFficient H-LIke [or HE-Like] command.\n" );
794 for( nelem = ipISO; nelem <
LIMELM; nelem++ )
796 for( ipLo=0; ipLo <=
NumLevRecomb[ipISO][nelem]; ipLo++ )
802 fprintf(
ioQQQ,
" iso_recomb_setup could not read line %li of %s.\n", i5,
811 if( i1!=nelem || i2!=ipLo )
813 fprintf(
ioQQQ,
" iso_recomb_setup detected insanity in %s.\n", chFilename[ipISO]);
815 " iso_recomb_setup: please recompile the data file with the COMPile RECOmb COEFficient H-LIke [or HE-Like] command.\n" );
822 double ThisCoef =
FFmtRead(chLine,&i3,chLine_LENGTH,&lgEOL);
832 RRCoef[ipISO][nelem][ipLo][i] = ThisCoef;
837 fprintf(
ioQQQ,
" iso_recomb_setup detected insanity in %s.\n", chFilename[ipISO]);
839 " iso_recomb_setup: please recompile the data file with the COMPile RECOmb COEFficient H-LIke [or HE-Like] command.\n" );
864 fprintf(
ioQQQ,
" iso_recomb_setup could not read last line of %s.\n", chFilename[ipISO]);
876 " iso_recomb_setup: the version of %s is not the current version.\n", chFilename[ipISO] );
878 " iso_recomb_setup: I expected to find the numbers %i %li %li %i and got %li %li %li %li instead.\n" ,
884 fprintf(
ioQQQ,
"Here is the line image:\n==%s==\n", chLine );
886 " iso_recomb_setup: please recompile the data file with the COMPile RECOmb COEFficient H-LIke [or HE-Like] command.\n" );
907 fprintf(ioRECOMB,
"%i\t%li\t%li\t%i\t%s isoelectronic sequence recomb data, created by COMPile RECOmb COEFficient H-LIke [or HE-Like] command, with %li %s levels, %li ion levels, and %i temperatures.\n",
918 for( nelem = ipISO; nelem <
LIMELM; nelem++ )
931 fprintf(ioRECOMB,
"%li\t%li", nelem, ipLo );
938 RRCoef[ipISO][nelem][ipLo][i] = log10(RadRecombReturn);
939 fprintf(ioRECOMB,
"\t%f",
RRCoef[ipISO][nelem][ipLo][i] );
941 fprintf(ioRECOMB,
"\n" );
947 fprintf(ioRECOMB,
"%li\t%li", nelem, NumLevRecomb[ipISO][nelem] );
958 fprintf(ioRECOMB,
"\t%f", log10(
TotalRecomb[ipISO][nelem][i] ) );
960 fprintf(ioRECOMB,
"\n" );
963 fprintf(ioRECOMB,
"%i\t%li\t%li\t%i\t%s isoelectronic sequence recomb data, created by COMPile RECOmb COEFficient [H-LIke/HE-Like] command, with %li %s levels, %li ion levels, and %i temperatures.\n",
976 fprintf(
ioQQQ,
"iso_recomb_setup: compilation complete, %s created.\n", chFilename[ipISO] );
977 fprintf(
ioQQQ,
"The compilation is completed successfully.\n");
991 1.00000, 1.30103, 1.69897, 2.00000, 2.30103, 2.69897, 3.00000,
992 3.30103, 3.69897, 4.00000, 4.30103, 4.69897, 5.00000, 5.30103,
993 5.69897, 6.00000, 6.30103, 6.69897, 7.00000 };
1004 TeDRCoef[i] = Te_over_Z1_Squared[i] + 2. * log10( (
double) nelem );
1022 pow( 10., 1.5* (TeDRCoef[NUM_DR_TEMPS-1] -
phycon.
alogte ) ) ;
1029 ASSERT( (ipTe >=0) && (ipTe < NUM_DR_TEMPS-1) );
1041 (TeDRCoef[ipTe+1]-TeDRCoef[ipTe]);
1043 rate = pow( 10., rate );
1052 ASSERT( rate >= 0. && rate < 1.0e-12 );
1061 static long int ipTe=-1;
1075 fprintf(
ioQQQ,
" TempInterp called with te out of bounds \n");
1085 while( (
phycon.
alogte < TempArray[ipTe] ) && ipTe > 0)
1093 while( (
phycon.
alogte > TempArray[ipTe+1] ) && ipTe < NumElements-1)
1097 ASSERT( (ipTe >=0) && (ipTe < NumElements-1) );
1101 && (
phycon.
alogte <= TempArray[ipTe+1] ) && ( ipTe < NumElements-1 ) );
1103 if( ValueArray[ipTe+1] == 0. && ValueArray[ipTe] == 0. )
1110 const int ORDER = 3;
1111 i0 =
max(
min(ipTe-ORDER/2,NumElements-ORDER-1),0);