36 double total, mtotal; \
37 for( i=0;i<N_H_MOLEC;i++) { \
39 for( j=0;j<N_H_MOLEC;j++) { \
40 total += c[i][j]*nprot[j]; \
42 if( fabs(total) > 1e-6*fabs(c[i][i]*nprot[i])) { \
43 fprintf(ioQQQ,"PROBLEM Subtotal1 %c %.2e\n",a,fabs(total)/fabs(c[i][i]*nprot[i])); \
44 fprintf(ioQQQ,"Species %li Total %g Diag %g\n",i,total,c[i][i]*nprot[i]); \
47 total = mtotal = 0.;for( j=0;j<N_H_MOLEC;j++) { total += bvec[j]*nprot[j]; mtotal += fabs(bvec[j]*nprot[j]); }\
48 if( fabs(total) > 1e-30 && fabs(total) > 1e-10*rtot) { \
49 fprintf(ioQQQ,"PROBLEM Subtotal2 %c %.2e\n",a,fabs(total)/mtotal); \
50 fprintf(ioQQQ,"RHS Total %g cf %g\n",total,mtotal); \
51 } else if( a == '.' && fabs(total) > 1e-7*mtotal) { \
52 fprintf(ioQQQ,"PROBLEM zone %li Hmole RHS conservation error %.2e of %.2e\n",nzone,total,mtotal); \
53 fprintf(ioQQQ,"(may be due to high rate equilibrium reactions)\n"); \
113 # define BIGERROR 1e-4
123 enum {DEBUG_LOC=
false};
125 if( DEBUG_LOC && (
nzone>140) )
127 fprintf(
ioQQQ,
"DEBUG hmole in \t%.2f\t%.5e\t%.5e",
133 fprintf(
ioQQQ,
"\n" );
176 enum {DEBUG_LOC=
false};
180 fprintf(
ioQQQ,
"DEBUG hmole out\t%i\t%.2f\t%.5e\t%.5e",
186 "\terror\t%.4e\tH0\t%.4e\tH+\t%.4e\tsink\t%.4e\t%.4e\tsour\t%.4e\t%.4e\tion\t%.4e\trec\t%.4e",
199 fprintf(
ioQQQ,
"\n" );
212 fprintf(
ioQQQ,
" PROBLEM hmole, zone %li: %d iters, %d bad; final error: %g lgSearch %i\n",
267 else if( te < 1200. )
271 else if( te < 3800. )
276 else if( te <= 1.4e4 )
316 static double teused=-1;
421 enum {DEBUG_LOC=
false};
426 fprintf(
ioQQQ,
"ph2lteee\t%.2e\t%.1e\t%.1e\n",
496 enum {DEBUG_LOC=
false};
514 enum {DEBUG_LOC=
false};
516 if( DEBUG_LOC &&
nzone>400)
625 double sqrtx = sqrt(1.+x);
629 double fshield = 0.965/
POW2(1.+x/b5) + 0.035/sqrtx *
682 fshield_H2g, fshield_H2s,
684 static double a_H2g, a_H2s,
717 a_H2g = 0.06 * log_G0_face + 1.32;
718 a_H2s = 0.375 * log_G0_face + 2.125;
720 e1_H2g = -0.05 * log_G0_face + 2.25;
721 e1_H2s = -0.125 * log_G0_face + 2.625;
723 e2_H2g = -0.005 * log_G0_face + 0.625;
725 b_H2g = -4.0e-3 * log_G0_face + 3.2e-2;
732 k_f_H2s =
MAX2(0.1,2.375 * log_G0_face - 1.875 );
735 k_H2g_to_H2s =
MAX2(1.,-1.75 * log_G0_face + 11.25);
750 fshield_H2g = 0.965/pow(1.+x_H2g/b5,e1_H2g) + b_H2g/pow(1.+x_H2g/b5,e2_H2g);
751 fshield_H2s = 0.965/pow(1.+x_H2s/b5,e1_H2s);
916 enum {DEBUG_LOC=
false};
920 fprintf(
ioQQQ,
" Solomon H2 dest rates: TH85 %.2e BD96 %.2e Big %.2e excit rates: TH85 %.2e Big %.2e\n",