36 return ((0.075*
pow2(x) - 1./6.)*
pow2(x) + 1.)*x;
37 else if( abs(x) <= 1./sqrt(DBL_EPSILON) )
40 return -log(sqrt(1. +
pow2(x)) - x);
42 return log(sqrt(1. +
pow2(x)) + x);
54 #define FREE_CHECK(PTR) { ASSERT( PTR != NULL ); free( PTR ); PTR = NULL; }
55 #define FREE_SAFE(PTR) { if( PTR != NULL ) free( PTR ); PTR = NULL; }
103 static const double AC0 = 3.e-9;
104 static const double AC1G = 4.e-8;
105 static const double AC2G = 7.e-8;
192 inline void Yfunc(
long,
long,
double,
double,
double,
double,
double,
double*,
double*,
203 inline double y2pa(
double,
double,
long,
double*);
205 inline double y2s(
double,
double,
long,
double*);
212 double*,
double*,
double*,
bool);
229 double*,
double*,
double*,
double*);
241 long int nelem, ion, ion_to;
333 for( nelem=0; nelem <
LIMELM; nelem++ )
335 for( ion=0; ion <= nelem+1; ion++ )
337 for( ion_to=0; ion_to <= nelem+1; ion_to++ )
366 for( nd=0; nd <
gv.
nBin; nd++ )
399 for( nd=0; nd <
gv.
nBin; nd++ )
433 fprintf(
ioQQQ,
" The code has run out of grain bins; increase NDUST and recompile.\n" );
445 for( ns=0; ns <
NSHL; ns++ )
449 for( nz=0; nz <
NCHS; nz++ )
498 for( nd=0; nd <
gv.
nBin; nd++ )
508 for( nz=0; nz <
NCHS; nz++ )
514 for( ns=0; ns <
NSHL; ns++ )
516 if(
gv.
bin[nd]->
sd[ns] != NULL )
530 for( nelem=0; nelem <
LIMELM; nelem++ )
559 for( nd=0; nd <
NDUST; nd++ )
562 for( nelem=0; nelem <
LIMELM; nelem++ )
598 fprintf(
ioQQQ,
" GrainsInit called.\n" );
626 for( nelem=0; nelem <
LIMELM; nelem++ )
656 fprintf(
ioQQQ,
" GrainsInit exits.\n" );
677 for( nd=0; nd <
gv.
nBin; nd++ )
679 double help,
atoms,p_rad,ThresInf,ThresInfVal,Emin,d[5];
680 long low1,low2,low3,lowm;
688 fprintf(
ioQQQ,
" Grain work function for %s has insane value: %.4e\n",
711 #ifndef IGNORE_QUANTUM_HEATING
731 for( nz=0; nz <
NCHS; nz++ )
765 help = ceil(-(1.2*
POW2(help)+3.9*help+0.2)/1.44);
772 help = ceil(-(0.7*
POW2(help)+2.5*help+0.8)/1.44);
783 fprintf(
ioQQQ,
" GrainsInit detected unknown Zmin type: %d\n" , zcase );
789 ASSERT( help > (
double)(LONG_MIN+1) );
803 while( low3-low2 > 1 )
805 lowm = (low2+low3)/2;
867 double Ethres = ( ns == 0 ) ? ThresInfVal :
gv.
bin[nd]->
sd[ns]->
ionPot;
877 sptr->p.reserve( len );
880 sptr->y01.reserve( len );
887 ASSERT( sptr->AvNr == NULL );
888 sptr->AvNr = (
double*)
MALLOC((
size_t)sptr->nData*
sizeof(double));
889 ASSERT( sptr->Ener == NULL );
890 sptr->Ener = (
double*)
MALLOC((
size_t)sptr->nData*
sizeof(double));
891 sptr->y01A.resize( sptr->nData );
893 for(
long n=0; n < sptr->nData; n++ )
898 sptr->y01A[n].reserve( len );
919 p_rad = 1./(1.+exp(20.-atoms));
942 for( nelem=0; nelem <
LIMELM; nelem++ )
947 for( nd=0; nd <
gv.
nBin; nd++ )
949 for( nelem=0; nelem <
LIMELM; nelem++ )
975 fprintf(
ioQQQ,
" There are %ld grain types turned on.\n",
gv.
nBin );
977 fprintf(
ioQQQ,
" grain depletion factors, dstfactor*GrainMetal=" );
978 for( nd=0; nd <
gv.
nBin; nd++ )
982 fprintf(
ioQQQ,
"\n" );
984 fprintf(
ioQQQ,
" nChrg =" );
985 for( nd=0; nd <
gv.
nBin; nd++ )
989 fprintf(
ioQQQ,
"\n" );
991 fprintf(
ioQQQ,
" lowest charge (e) =" );
992 for( nd=0; nd <
gv.
nBin; nd++ )
996 fprintf(
ioQQQ,
"\n" );
998 fprintf(
ioQQQ,
" lgDustFunc flag for user requested custom depth dependence:" );
999 for( nd=0; nd <
gv.
nBin; nd++ )
1003 fprintf(
ioQQQ,
"\n" );
1005 fprintf(
ioQQQ,
" Quantum heating flag:" );
1006 for( nd=0; nd <
gv.
nBin; nd++ )
1010 fprintf(
ioQQQ,
"\n" );
1013 fprintf(
ioQQQ,
" NU(Ryd), Abs cross sec per proton\n" );
1015 fprintf(
ioQQQ,
" Ryd " );
1016 for( nd=0; nd <
gv.
nBin; nd++ )
1020 fprintf(
ioQQQ,
"\n" );
1025 for( nd=0; nd <
gv.
nBin; nd++ )
1029 fprintf(
ioQQQ,
"\n" );
1032 fprintf(
ioQQQ,
" NU(Ryd), Sct cross sec per proton\n" );
1034 fprintf(
ioQQQ,
" Ryd " );
1035 for( nd=0; nd <
gv.
nBin; nd++ )
1039 fprintf(
ioQQQ,
"\n" );
1044 for( nd=0; nd <
gv.
nBin; nd++ )
1048 fprintf(
ioQQQ,
"\n" );
1051 fprintf(
ioQQQ,
" NU(Ryd), Q abs\n" );
1053 fprintf(
ioQQQ,
" Ryd " );
1054 for( nd=0; nd <
gv.
nBin; nd++ )
1058 fprintf(
ioQQQ,
"\n" );
1063 for( nd=0; nd <
gv.
nBin; nd++ )
1067 fprintf(
ioQQQ,
"\n" );
1070 fprintf(
ioQQQ,
" NU(Ryd), Q sct\n" );
1072 fprintf(
ioQQQ,
" Ryd " );
1073 for( nd=0; nd <
gv.
nBin; nd++ )
1077 fprintf(
ioQQQ,
"\n" );
1082 for( nd=0; nd <
gv.
nBin; nd++ )
1086 fprintf(
ioQQQ,
"\n" );
1089 fprintf(
ioQQQ,
" NU(Ryd), asymmetry factor\n" );
1091 fprintf(
ioQQQ,
" Ryd " );
1092 for( nd=0; nd <
gv.
nBin; nd++ )
1096 fprintf(
ioQQQ,
"\n" );
1101 for( nd=0; nd <
gv.
nBin; nd++ )
1105 fprintf(
ioQQQ,
"\n" );
1108 fprintf(
ioQQQ,
" GrainsInit exits.\n" );
1122 static char chFile[] =
"auger_spectrum.dat";
1126 sscanf( chString,
"%ld", &version );
1129 fprintf(
ioQQQ,
" File %s has wrong version number\n", chFile );
1130 fprintf(
ioQQQ,
" please update your installation...\n" );
1142 res = sscanf( chString,
"%ld", &nelem );
1153 res = sscanf( chString,
"%lu", &ptr->
nSubShell );
1166 res = sscanf( chString,
"%lu", &ss );
1167 ASSERT( res == 1 && ns == ss );
1170 res = sscanf( chString,
"%le", &ptr->
IonThres[ns] );
1175 res = sscanf( chString,
"%lu", &ptr->
nData[ns] );
1181 for(
unsigned long n=0; n < ptr->
nData[ns]; n++ )
1184 res = sscanf(chString,
"%le %le",&ptr->
Energy[ns][n],&ptr->
AvNumber[ns][n]);
1204 long i, ipLo, nelem;
1222 for( i=ipLo; i < ipEnd; i++ )
1251 temp[i] +=
gv.
bin[nd]->
sd[ns]->
p[i];
1262 temp[i] +=
gv.
bin[nd]->
sd[ns]->
p[i];
1267 for( i=ipBegin; i < ipEnd && !
gv.
lgWD01; i++ )
1278 for( i=ipLo; i < ipEnd; i++ )
1281 gv.
bin[nd]->
sd[ns]->
p[i] /= temp[i];
1283 gv.
bin[nd]->
sd[ns]->
p[i] = ( ns == 0 ) ? 1.f : 0.f;
1294 ipLo =
max( sptr->
ipLo, ipBegin );
1299 for( i=ipLo; i < ipEnd; i++ )
1301 double elec_en,yzero,yone;
1304 yzero =
y0psa( nd, ns, i, elec_en );
1308 yone =
y1psa( nd, i, elec_en );
1325 for( n=0; n < sptr->
nData; n++ )
1327 sptr->
y01A[n].realloc( ipEnd );
1329 for( i=ipLo; i < ipEnd; i++ )
1331 double yzero = sptr->
AvNr[n] *
y0psa( nd, ns, i, sptr->
Ener[n] );
1335 double yone =
y1psa( nd, i, sptr->
Ener[n] );
1357 fprintf(
ioQQQ,
" Could not read from %s\n", chFile );
1359 fprintf(
ioQQQ,
" EOF reached\n");
1363 while( chLine[0] ==
'#' );
1366 str = strstr(chLine,
"#");
1385 fprintf(
ioQQQ,
" InitEmissivities starts\n" );
1386 fprintf(
ioQQQ,
" ND Tdust Emis BB Check 4pi*a^2*<Q>\n" );
1395 for( nd=0; nd <
gv.
nBin; nd++ )
1407 for( nd=0; nd <
gv.
nBin; nd++ )
1418 mul = sqrt((
double)(NDEMS-NTOP));
1420 mul = (double)NTOP + sqrt((
double)
NDEMS);
1424 for( nd=0; nd <
gv.
nBin; nd++ )
1435 for( i=0; i < 2000; i++ )
1438 z = pow(10.,-40. + (
double)i/50.);
1443 printf(
" input %.6e temp %.6e output %.6e rel. diff. %.6e\n",z,exp(y),x,(x-z)/z);
1482 x = 0.999*log(DBL_MAX);
1493 ExpM1 = arg*(1. + arg/2.);
1498 ExpM1 = exp(
MIN2(x,arg)) - 1.;
1512 if( Planck1/integral1 < DBL_EPSILON && Planck2/integral2 < DBL_EPSILON )
1515 integral1 += Planck1;
1516 integral2 += Planck2;
1522 fprintf(
ioQQQ,
" %4ld %11.4e %11.4e %11.4e %11.4e\n", nd, tdust,
1523 integral2, integral1/4./5.67051e-5/
powi(tdust,4), integral2*4./integral1 );
1526 ASSERT( integral2 > 0. );
1538 for( nz=0; nz <
NCHS; nz++ )
1558 for( nz=0; nz <
NCHS; nz++ )
1567 for( nz=0; nz <
NCHS; nz++ )
1580 double GrnStdDpth_v;
1614 ASSERT( GrnStdDpth_v > 0. && GrnStdDpth_v <= 1.000001 );
1616 return GrnStdDpth_v;
1628 static double tesave=-1.f;
1629 static long int nzonesave=-1;
1654 fprintf(
ioQQQ,
" grain5 calling GrainChargeTemp\n");
1669 long nelem, ion, ion_to, nd;
1674 for( nd=0; nd <
gv.
nBin; nd++ )
1740 for( nelem=0; nelem <
LIMELM; nelem++ )
1742 for( ion=0; ion <= nelem+1; ion++ )
1744 for( ion_to=0; ion_to <= nelem+1; ion_to++ )
1782 bool lgAnyNegCharge =
false;
1800 static long int oldZone = -1;
1801 static double oldTe = -DBL_MAX,
1827 for( nelem=0; nelem <
LIMELM; nelem++ )
1829 for( ion=0; ion <= nelem+1; ion++ )
1831 for( ion_to=0; ion_to <= nelem+1; ion_to++ )
1843 for( nd=0; nd <
gv.
nBin; nd++ )
1846 double ChTdBracketLo = 0., ChTdBracketHi = -DBL_MAX;
1865 fprintf(
ioQQQ,
" >>GrainChargeTemp starting grain %s\n",
1871 for( i=0; i < relax*CT_LOOP_MAX && delta >
TOLER; ++i )
1875 double TdBracketLo = 0., TdBracketHi = -DBL_MAX;
1876 double ThresEst = 0.;
1906 fprintf(
ioQQQ,
" >>loop %ld BracketLo %.6e BracketHi %.6e",
1907 j, TdBracketLo, TdBracketHi );
1919 fprintf(
ioQQQ,
" converged\n" );
1925 TdBracketHi = oldTemp2;
1927 TdBracketLo = oldTemp2;
1937 if( TdBracketHi > TdBracketLo )
1941 if( ( j >= 2 && TdBracketLo > 0. ) ||
1947 fprintf(
ioQQQ,
" bisection\n" );
1952 fprintf(
ioQQQ,
" iteration\n" );
1958 fprintf(
ioQQQ,
" iteration\n" );
1968 ChTdBracketHi = oldtemp;
1970 ChTdBracketLo = oldtemp;
1980 fprintf(
ioQQQ,
" PROBLEM temperature of grain species %s (Tg=%.3eK) not converged\n",
2005 delta = ( delta < 0.9*log(DBL_MAX) ) ?
2006 ThermRatio*fabs(
POW2(ratio)*exp(delta)-1.) : DBL_MAX;
2012 strncpy( which,
"iteration",
sizeof(which) );
2022 if( ChTdBracketHi > ChTdBracketLo )
2026 if( ( i >= 2 && ChTdBracketLo > 0. ) ||
2032 strncpy( which,
"bisection",
sizeof(which) );
2039 fprintf(
ioQQQ,
" >>GrainChargeTemp finds delta=%.4e, ", delta );
2040 fprintf(
ioQQQ,
" old/new temp=%.5e %.5e, ", oldtemp,
gv.
bin[nd]->
tedust );
2042 fprintf(
ioQQQ,
"doing another %s\n", which );
2044 fprintf(
ioQQQ,
"converged\n" );
2049 fprintf(
ioQQQ,
" PROBLEM charge/temperature not converged for %s zone %.2f\n",
2057 lgAnyNegCharge =
true;
2107 enum {DEBUG_LOC=
false};
2111 fprintf(
ioQQQ,
" DEBUG grn chr nz\t%.2f\teden\t%.3e\tnd\t%li",
2115 fprintf(
ioQQQ,
"\tne\t%.2e\tAveDustZ\t%.2e\t%.2e\t%.2e\t%.2e",
2120 fprintf(
ioQQQ,
"\n");
2126 fprintf(
ioQQQ,
" %s Pot %.5e Thermal %.5e GasCoolColl %.5e" ,
2128 fprintf(
ioQQQ,
" GasPEHeat %.5e GasThermHeat %.5e ChemHeat %.5e\n\n" ,
2200 printf(
"wd test: proton fraction %.5e Total DustZ %.6f heating/cooling rate %.5e %.5e\n",
2209 enum {DEBUG_LOC=
false};
2213 fprintf(
ioQQQ,
" %.2f Grain surface charge transfer rates\n",
fnzone );
2215 for( nelem=0; nelem <
LIMELM; nelem++ )
2222 for( ion_to=0; ion_to <= nelem+1; ion_to++ )
2226 printf(
" %ld->%ld %.2e", ion, ion_to,
2236 fprintf(
ioQQQ,
" %.2f Grain contribution to electron density %.2e\n",
2239 fprintf(
ioQQQ,
" Grain electons: " );
2240 for( nd=0; nd <
gv.
nBin; nd++ )
2248 fprintf(
ioQQQ,
" %.2e", sum );
2250 fprintf(
ioQQQ,
"\n" );
2252 fprintf(
ioQQQ,
" Grain potentials:" );
2253 for( nd=0; nd <
gv.
nBin; nd++ )
2257 fprintf(
ioQQQ,
"\n" );
2259 fprintf(
ioQQQ,
" Grain temperatures:" );
2260 for( nd=0; nd <
gv.
nBin; nd++ )
2264 fprintf(
ioQQQ,
"\n" );
2299 netloss0 = -DBL_MAX,
2300 netloss1 = -DBL_MAX,
2315 for( nz=0; nz <
NCHS; nz++ )
2380 power = (int)(log(step)/log((
double)stride0));
2381 power =
MAX2(power,0);
2382 xxx =
powi((
double)stride0,power);
2392 UpdatePot( nd, Zlo, stride, rate_up, rate_dn );
2397 netloss0 = rate_up[0] - rate_dn[0];
2400 if( netloss0*netloss1 <= 0. )
2413 UpdatePot( nd, Zlo, stride, rate_up, rate_dn );
2416 if( netloss0*netloss1 > 0. ) {
2417 fprintf(
ioQQQ,
" insanity: could not bracket grain charge for %s\n",
gv.
bin[nd]->
chDstLab );
2427 netloss0 = rate_up[0] - rate_dn[0];
2430 netloss1 = rate_up[nz+1] - rate_dn[nz+1];
2432 if( netloss0*netloss1 <= 0. )
2438 netloss0 = netloss1;
2440 UpdatePot( nd, Zlo, stride, rate_up, rate_dn );
2443 ASSERT( netloss0*netloss1 <= 0. );
2451 for( i=0; i < loopMax; i++ )
2453 GetFracPop( nd, Zlo, rate_up, rate_dn, &newZlo );
2462 UpdatePot( nd, Zlo, 1, rate_up, rate_dn );
2466 fprintf(
ioQQQ,
" insanity: could not converge grain charge for %s\n",
gv.
bin[nd]->
chDstLab );
2487 *ThermRatio = ( crate > 0. ) ? csum3/crate : 0.;
2489 ASSERT( *ThermRatio >= 0. );
2495 fprintf(
ioQQQ,
"\n" );
2497 crate = csum1a = csum1b = csum2 = csum3 = 0.;
2508 fprintf(
ioQQQ,
" ElecEm rate1a=%.4e, rate1b=%.4e, ", csum1a, csum1b );
2509 fprintf(
ioQQQ,
"rate2=%.4e, rate3=%.4e, sum=%.4e\n", csum2, csum3, crate );
2512 fprintf(
ioQQQ,
" rate1a/sum=%.4e, rate1b/sum=%.4e, ", csum1a/crate, csum1b/crate );
2513 fprintf(
ioQQQ,
"rate2/sum=%.4e, rate3/sum=%.4e\n", csum2/crate, csum3/crate );
2516 crate = csum1 = csum2 = 0.;
2524 fprintf(
ioQQQ,
" ElecRc rate1=%.4e, rate2=%.4e, sum=%.4e\n", csum1, csum2, crate );
2527 fprintf(
ioQQQ,
" rate1/sum=%.4e, rate2/sum=%.4e\n", csum1/crate, csum2/crate );
2530 fprintf(
ioQQQ,
" Charging rates:" );
2533 fprintf(
ioQQQ,
" Zg %ld up %.4e dn %.4e",
2536 fprintf(
ioQQQ,
"\n" );
2543 fprintf(
ioQQQ,
"\n" );
2545 fprintf(
ioQQQ,
" >Grain potential:%12.12s %.6fV",
2549 fprintf(
ioQQQ,
" Thres[%ld]: %.4f V ThresVal[%ld]: %.4f V",
2553 fprintf(
ioQQQ,
"\n" );
2584 rate = *sum1 + *sum2;
2604 #ifndef IGNORE_GRAIN_ION_COLLISIONS
2605 for( ion=0; ion <=
LIMELM; ion++ )
2607 double CollisionRateAll = 0.;
2609 for( nelem=
MAX2(ion-1,0); nelem <
LIMELM; nelem++ )
2620 if( CollisionRateAll > 0. )
2626 *sum2 += CollisionRateAll*eta;
2631 rate = *sum1 + *sum2;
2637 ASSERT( *sum1 >= 0. && *sum2 >= 0. );
2672 rate = *sum1a + *sum1b + *sum2 + *sum3;
2718 # ifndef IGNORE_GRAIN_ION_COLLISIONS
2719 for( ion=0; ion <=
LIMELM; ion++ )
2721 double CollisionRateAll = 0.;
2723 for( nelem=
MAX2(ion-1,0); nelem <
LIMELM; nelem++ )
2734 if( CollisionRateAll > 0. )
2740 *sum2 += CollisionRateAll*eta;
2752 rate = *sum1a + *sum1b + *sum2 + *sum3;
2759 ASSERT( *sum1a >= 0. && *sum1b >= 0. && *sum2 >= 0. && *sum3 >= 0. );
2808 *eta = (1. - nu/tau)*(1. + sqrt(2./(tau - 2.*nu)));
2809 *xi = (1. - nu/(2.*tau))*(1. + 1./sqrt(tau - nu));
2813 *eta = 1. + sqrt(
PI/(2.*tau));
2814 *xi = 1. + 0.75*sqrt(
PI/(2.*tau));
2818 double theta_nu =
ThetaNu(nu);
2820 double xxx = 1. + 1./sqrt(4.*tau+3.*nu);
2821 *eta =
POW2(xxx)*exp(-theta_nu/tau);
2823 *xi = (1. + nu/(2.*tau))*(1. + 1./sqrt(3./(2.*tau)+3.*nu))*exp(-theta_nu/tau);
2828 xxx = 0.25*pow(nu/tau,0.75)/(pow(nu/tau,0.75) + pow((25.+3.*nu)/5.,0.75)) +
2829 (1. + 0.75*sqrt(
PI/(2.*tau)))/(1. + sqrt(
PI/(2.*tau)));
2830 *xi = (
MIN2(xxx,1.) + theta_nu/(2.*tau))*(*eta);
2834 ASSERT( *eta >= 0. && *xi >= 0. );
2853 const double REL_TOLER = 10.*DBL_EPSILON;
2856 double xi_nu = 1. + 1./sqrt(3.*nu);
2857 double xi_nu2 =
POW2(xi_nu);
2863 double fnu = 2.*xi_nu2 - 1. - nu*xi_nu*
POW2(xi_nu2 - 1.);
2864 double dfdxi = 4.*xi_nu - nu*((5.*xi_nu2 - 6.)*xi_nu2 + 1.);
2866 xi_nu2 =
POW2(xi_nu);
2867 xxx = fabs(old-xi_nu);
2868 }
while( xxx > REL_TOLER*xi_nu );
2870 theta_nu = nu/xi_nu - 1./(2.*xi_nu2*(xi_nu2-1.));
2909 fprintf(
ioQQQ,
" %ld/%ld", Zlo, stride );
2924 Zg = Zlo + nz*stride;
2928 for( zz=0; zz <
NCHS-1; zz++ )
2942 for( zz=ind-1; zz >= nz; zz-- )
2960 ASSERT( rate_up[nz] >= 0. && rate_dn[nz] >= 0. );
2973 HighEnergy =
MAX2(HighEnergy,
3013 for( i=0; i < 2; i++ )
3018 sum = pop[i][0] = 1.;
3022 if( rate_dn[nz] > 10.*rate_up[nz-1]/sqrt(DBL_MAX) )
3024 pop[i][j] = pop[i][j-1]*rate_up[nz-1]/rate_dn[nz];
3029 for( k=0; k < j; k++ )
3037 if( pop[i][j] > sqrt(DBL_MAX) )
3039 for( k=0; k <= j; k++ )
3041 pop[i][k] /= DBL_MAX/10.;
3051 netloss[i] += pop[i][j]*(rate_up[nz] - rate_dn[nz]);
3057 if( netloss[0]*netloss[1] > 0. )
3058 *newZlo = ( netloss[1] > 0. ) ? Zlo + 1 : Zlo - 1;
3069 ( *newZlo == Zlo && pop[1][
gv.
bin[nd]->
nChrg-2] < DBL_EPSILON ) ) )
3080 printf(
" fnzone %.2f nd %ld Zlo %ld newZlo %ld netloss %.4e %.4e nChrg %ld lgRedo %d\n",
3081 fnzone, nd, Zlo, *newZlo, netloss[0], netloss[1],
gv.
bin[nd]->
nChrg, lgRedo );
3088 fprintf(
ioQQQ,
" could not converge charge state populations for %s\n",
gv.
bin[nd]->
chDstLab );
3093 if( *newZlo == Zlo )
3095 double frac0, frac1;
3097 double test1, test2, test3,
x1,
x2;
3101 frac0 = netloss[1]/(netloss[1]-netloss[0]);
3102 frac1 = -netloss[0]/(netloss[1]-netloss[0]);
3112 test1 = test2 = test3 = 0.;
3120 x1 = fabs(test1-1.);
3121 x2 = fabs(test3/test2);
3160 const double INV_DELTA_E =
EVRYD/3.;
3161 const double CS_PDT = 1.2e-17;
3199 long len = nfill + 10 - ipLo;
3230 for( i=
max(ipLo,ipStart); i < nfill; i++ )
3234 double Yp1,Ys1,Ehp,Ehs,yp,ya,ys,eyp,eya,eys;
3235 double yzero,yone,
Ehi,Ehi_band,Wcorr,Eel;
3239 eyp = eya = eys = 0.;
3242 Ehi = Ehi_band = anu[i] - ThresInfVal - Emin;
3243 Wcorr = ThresInfVal + Emin - GrainPot;
3244 Eel = anu[i] - DustWorkFcn;
3245 yzero =
y0b( nd, nz, i );
3246 yone = sptr->
y01[i];
3247 Yfunc(nd,nz,yzero*yone,sptr->
p[i],Elo,Ehi,Eel,&Yp1,&Ys1,&Ehp,&Ehs);
3255 for( ns=1; ns < nsmax; ns++ )
3259 if( i < sptr->ipLo )
3262 Ehi = Ehi_band + Wcorr - sptr->
ionPot;
3263 Eel = anu[i] - sptr->
ionPot;
3264 Yfunc(nd,nz,sptr->
y01[i],sptr->
p[i],Elo,Ehi,Eel,&Yp1,&Ys1,&Ehp,&Ehs);
3271 long nmax = sptr->
nData;
3272 for( n=0; n < nmax; n++ )
3274 double max = sptr->
AvNr[n]*sptr->
p[i];
3275 Ehi = sptr->
Ener[n] - GrainPot;
3276 Eel = sptr->
Ener[n];
3277 Yfunc(nd,nz,sptr->
y01A[n][i],max,Elo,Ehi,Eel,&Yp1,&Ys1,&Ehp,&Ehs);
3306 len = nfill + 10 - ipLo;
3323 c1 = -CS_PDT*(double)Zg;
3327 for( i=
max(ipLo,ipStart); i < nfill; i++ )
3329 double x = (anu[i] - ThresInf)*INV_DELTA_E;
3330 double cs = c1*x/
POW2(1.+(1./3.)*
POW2(x));
3354 for( i=
max(ipLo,ipStart); i < nfill; i++ )
3356 double cs1,cs2,cs_tot,cool1,cool2,ehat1,ehat2;
3359 PE_init(nd,nz,i,&cs1,&cs2,&cs_tot,&cool1,&cool2,&ehat1,&ehat2);
3361 gv.
bin[nd]->
chrg[nz]->
fac1[i] = cs_tot*anu[i]-cs1*cool1-cs2*cool2;
3428 # if defined( WD_TEST2 ) || defined( IGNORE_THERMIONIC )
3455 y2pr =
y2pa( Elo, Ehi, Zg, Ehp );
3462 *Yp = y2pr*
min(y01,maxval);
3464 y2sec =
y2s( Elo, Ehi, Zg, Ehs );
3467 else if( pcase ==
PE_SIL )
3471 fprintf(
ioQQQ,
" Yfunc: unknown type for PE effect: %d\n" , pcase );
3477 *Ys = y2sec*f3*
min(y01,maxval);
3500 yzero =
y0b01( nd, nz, i );
3505 if( Eph <= 20./
EVRYD )
3506 yzero =
y0b01( nd, nz, i );
3507 else if( Eph < 50./
EVRYD )
3509 double y0a =
y0b01( nd, nz, i );
3512 double frac = log(Eph*(
EVRYD/20.))*1.0913566679372915;
3514 yzero = y0a * exp(log(y0b/y0a)*frac);
3542 yzero = xv/((1./9.e-3) + (3.7e-2/9.e-3)*xv);
3546 yzero = xv/(2.+10.*xv);
3549 fprintf(
ioQQQ,
" y0b01: unknown type for PE effect: %d\n" , pcase );
3564 double yzero, leola;
3577 yzero =
gv.
bin[nd]->
sd[ns]->
p[i]*leola*(1. - leola*log(1.+1./leola));
3580 double x = 1./leola;
3581 yzero =
gv.
bin[nd]->
sd[ns]->
p[i]*(((-1./5.*x+1./4.)*x-1./3.)*x+1./2.);
3594 double alpha, beta, af, bf, yone;
3600 bf =
pow2(beta) - 2.*beta + 2. - 2.*exp(-beta);
3602 bf = ((1./60.*beta - 1./12.)*beta + 1./3.)*
pow3(beta);
3606 af =
pow2(alpha) - 2.*alpha + 2. - 2.*exp(-alpha);
3608 af = ((1./60.*alpha - 1./12.)*alpha + 1./3.)*
pow3(alpha);
3610 yone =
pow2(beta/alpha)*af/bf;
3632 *Ehp = 0.5*Ehi*(1.-2.*x)/(1.-3.*x);
3634 ytwo = ( abs(x) > 1e-4 ) ? (1.-3.*x)/
pow3(1.-x) : 1. - (3. + 8.*x)*x*x;
3635 ASSERT( *Ehp > 0. && *Ehp <= Ehi && ytwo > 0. && ytwo <= 1. );
3647 *Ehp = 0.5*(Elo+
Ehi);
3649 ASSERT( *Ehp >= Elo && *Ehp <= Ehi );
3682 double x2 = x*x, x3 = x2*x, x4 = x3*x, x5 = x4*x;
3683 double yh2 = yh*yh, yh3 = yh2*yh, yh4 = yh3*yh, yh5 = yh4*yh;
3684 double help1 = 2.*x-yh;
3685 double help2 = (6.*x3-15.*yh*x2+12.*yh2*x-3.*yh3)/4.;
3686 double help3 = (22.*x5-95.*yh*x4+164.*yh2*x3-141.*yh3*x2+60.*yh4*x-10.*yh5)/16.;
3687 N0 = yh*(help1 - help2 + help3)/x2;
3689 help1 = (3.*x-yh)/3.;
3690 help2 = (15.*x3-25.*yh*x2+15.*yh2*x-3.*yh3)/20.;
3691 help3 = (1155.*x5-3325.*yh*x4+4305.*yh2*x3-2961.*yh3*x2+1050.*yh4*x-150.*yh5)/1680.;
3692 E0 =
ETILDE*yh2*(help1 - help2 + help3)/x2;
3696 double sR0 = (1. + yl*yl);
3697 double sqR0 = sqrt(sR0);
3698 double sqRh = sqrt(1. + x*x);
3699 double alpha = sqRh/(sqRh - 1.);
3703 double z = yh*(yh - 2.*yl)/sR0;
3704 N0 = ((((7./256.*z-5./128.)*z+1./16.)*z-1./8.)*z+1./2.)*z/(sqRh-1.);
3706 double yl2 = yl*yl, yl3 = yl2*yl, yl4 = yl3*yl;
3707 double help1 = yl/2.;
3708 double help2 = (2.*yl2-1.)/3.;
3709 double help3 = (6.*yl3-9.*yl)/8.;
3710 double help4 = (8.*yl4-24.*yl2+3.)/10.;
3712 E0 = -alpha*Ehi*(((help4*h + help3)*h + help2)*h + help1)*h/sqR0;
3716 N0 = alpha*(1./sqR0 - 1./sqRh);
3717 E0 = alpha*
ETILDE*(
ASINH(x*sqR0 + yl*sqRh) - yh/sqRh);
3720 ASSERT( N0 > 0. && N0 <= 1. );
3724 ASSERT( *Ehs > 0. && *Ehs <= Ehi );
3744 double sqRh = sqrt(1. + x2);
3745 double alpha = sqRh/(sqRh - 1.);
3751 *Ehs = Ehi - (Ehi-
Elo)*((-37./840.*x2 + 1./10.)*x2 + 1./3.);
3754 ASSERT( *Ehs >= Elo && *Ehs <= Ehi );
3778 for( nelem=
LIMELM-1; nelem >= 0; nelem-- )
3782 for( ion=nelem+1; ion >= 0; ion-- )
3787 high =
MAX2(high,ion);
3798 bool lgAllIonStages)
3816 for( i=1; i <=
LIMELM; i++ )
3821 phi_s_up[i] = -DBL_MAX;
3827 for( nelem=0; nelem <
LIMELM; nelem++ )
3831 for( ion=0; ion <= nelem+1; ion++ )
3855 double *ThresInfVal,
3857 double *ThresSurfVal,
3860 bool lgUseTunnelCorr)
3899 fprintf(
ioQQQ,
" GetPotValues detected unknown type for ionization pot: %d\n" , pcase );
3905 IP_v =
MAX2(IP,IP_v);
3910 double help = fabs(dZg+1);
3913 if( lgUseTunnelCorr )
3924 *ThresInf = IP - *Emin;
3925 *ThresInfVal = IP_v - *Emin;
3926 *ThresSurf = *ThresInf;
3927 *ThresSurfVal = *ThresInfVal;
3933 *ThresInfVal = IP_v;
3934 *ThresSurf = *ThresInf - dstpot;
3935 *ThresSurfVal = *ThresInfVal - dstpot;
3950 const double phi_s_up[],
3951 const double phi_s_dn[],
3969 phi_s = phi_s_up[0];
3977 *ChemEn -= (
realnum)(phi_s - phi_s_up[0]);
3980 phi_s = phi_s_up[save-ion];
3992 phi_s = phi_s_dn[0];
4000 *ChemEn += (
realnum)(phi_s - phi_s_dn[0]);
4005 phi_s = phi_s_dn[ion-save];
4033 # ifndef IGNORE_GRAIN_ION_COLLISIONS
4039 double fac1 = gptr->
FracPop*fac0;
4044 for( ion=0; ion <=
LIMELM; ion++ )
4047 double eta, fac2, xi;
4058 for( nelem=
MAX2(0,ion-1); nelem <
LIMELM; nelem++ )
4082 for( nelem=0; nelem <
LIMELM; nelem++ )
4088 for( nd=0; nd <
gv.
nBin; nd++ )
4105 for( nelem=0; nelem <
LIMELM; nelem++ )
4118 bool lgChangeCS_PDT;
4141 for( nd=0; nd <
gv.
nBin; nd++ )
4165 if( gptr->
DustZ <= -1 && i >= ipLo )
4191 double EhatThermionic,
4241 double bolflux1 = 0.;
4245 bool lgReEvaluate1 = gptr->
hcon1 < 0.;
4246 bool lgReEvaluate2 = gptr->
hots1 < 0.;
4247 bool lgReEvaluate = lgReEvaluate1 || lgReEvaluate2;
4253 for( i=0; i < loopmax; i++ )
4283 gptr->
hcon1 = hcon1;
4287 hcon1 = gptr->
hcon1;
4305 if( gptr->
DustZ <= -1 )
4309 gptr->
hots1 = hots1;
4315 hots1 = gptr->
hots1;
4338 ASSERT( hcon1 >= 0. && hots1 >= 0. && hla1 >= 0. && bolflux1 >= 0. && pe1 >= 0. );
4350 fprintf(
ioQQQ,
" Zg %ld bolflux: %.4e\n", gptr->
DustZ,
4468 fprintf(
ioQQQ,
" >GrainTemperature finds %s Tdst %.5e hcon %.4e ",
4470 fprintf(
ioQQQ,
"hots %.4e dcheat %.4e GrainCoolTherm %.4e\n",
4509 *ehat1 = gptr->
ehat[i];
4530 if( gptr->
DustZ <= -1 )
4535 ASSERT( *ehat1 > 0. && *cool1 > 0. );
4545 if( gptr->
DustZ <= -1 && i >= ipLo2 )
4554 ASSERT( *ehat2 > 0. && *cool2 > 0. );
4577 double Accommodation,
4578 CollisionRateElectr,
4605 const double H2_FORMATION_GRAIN_HEATING[
H2_TOP] = { 0.20, 0.4, 1.72 };
4636 double ChemEn1 = 0.;
4642 for( ion=0; ion <=
LIMELM; ion++ )
4651 CollisionRateIon = 0.;
4653 CoolPotentialGas = 0.;
4654 HeatRecombination = 0.;
4661 for( nelem=
MAX2(0,ion-1); nelem <
LIMELM; nelem++ )
4665 double CollisionRateOne;
4670 #if defined( IGNORE_GRAIN_ION_COLLISIONS )
4672 #elif defined( WD_TEST2 )
4673 Stick = ( ion == gptr->
RecomZ0[nelem][ion] ) ?
4676 Stick = ( ion == gptr->
RecomZ0[nelem][ion] ) ?
4683 CollisionRateIon += CollisionRateOne;
4692 if( ion >= gptr->
RecomZ0[nelem][ion] )
4694 CoolPotential += CollisionRateOne * (double)ion *
4696 CoolPotentialGas += CollisionRateOne *
4697 (
double)gptr->
RecomZ0[nelem][ion] *
4702 CoolPotential += CollisionRateOne * (double)ion *
4704 CoolPotentialGas += CollisionRateOne *
4705 (
double)gptr->
RecomZ0[nelem][ion] *
4712 HeatRecombination += CollisionRateOne *
4714 HeatChem += CollisionRateOne * gptr->
ChemEn[nelem][ion];
4728 CoolPotential *= eta*
EN1RYD;
4729 CoolPotentialGas *= eta*
EN1RYD;
4731 HeatRecombination *= eta*
EN1RYD;
4737 Heat1 += HeatCollisions - CoolPotential + HeatRecombination - CoolEmitted;
4742 Cool1 += HeatCollisions - CoolEmitted - CoolPotentialGas;
4744 ChemEn1 += HeatChem;
4748 HeatIons += gptr->
FracPop*Heat1;
4752 fprintf(
ioQQQ,
" Zg %ld ions heat/cool: %.4e %.4e\n", gptr->
DustZ,
4767 CollisionRateElectr = Stick*
dense.
eden*ve;
4778 CoolPotential = CollisionRateElectr * (double)ion*gptr->
PotSurfInc*eta*
EN1RYD;
4786 HeatCollisions = 0.;
4788 HeatRecombination = 0.;
4796 CoolBounce = CollisionRateElectr *
4798 CoolBounce =
MAX2(CoolBounce,0.);
4803 HeatElectrons = HeatCollisions-CoolPotential+HeatRecombination+HeatBounce-CoolBounce;
4804 Heat1 += HeatElectrons;
4806 CoolElectrons = HeatCollisions+HeatBounce-CoolBounce;
4807 Cool1 += CoolElectrons;
4811 fprintf(
ioQQQ,
" Zg %ld electrons heat/cool: %.4e %.4e\n", gptr->
DustZ,
4834 HeatTot += gptr->
FracPop*Heat1;
4837 CoolTot += gptr->
FracPop*Cool1;
4852 #ifndef IGNORE_GRAIN_ION_COLLISIONS
4862 H2_FORMATION_GRAIN_HEATING[ipH2]*
EN1EV;
4866 CollisionRateMol = 0.;
4873 #ifndef IGNORE_GRAIN_ION_COLLISIONS
4877 CollisionRateMol = 0.;
4884 HeatMolecules = HeatCollisions - CoolEmitted +
gv.
bin[nd]->
ChemEnH2;
4885 HeatTot += HeatMolecules;
4888 CoolMolecules = HeatCollisions - CoolEmitted;
4889 CoolTot += CoolMolecules;
4915 fprintf(
ioQQQ,
" molecules heat/cool: %.4e %.4e heatcor: %.4e\n",
4980 for( nd=0; nd <
gv.
nBin; nd++ )
5008 phi2lm =
POW2(psi)*alam;
5011 for( loop = 0; loop < 50 && fabs(corr-1.) > 0.001; loop++ )
5017 g0 = 1.5045*si*sqrt(1.+0.4418*si*si);
5018 g2 = si/(1.329 +
POW3(si));
5026 g0 = 1.5045*si*sqrt(1.+0.4418*si*si);
5027 g2 = si/(1.329 +
POW3(si));
5028 fdrag += fac*
dense.
eden*(g0 + phi2lm*g2);
5032 g0 = 1.5045*si*sqrt(1.+0.4418*si*si);
5037 g0 = 1.5045*si*sqrt(1.+0.4418*si*si);
5038 g2 = si/(1.329 +
POW3(si));
5048 corr = sqrt(volmom/fdrag);
5060 fprintf(
ioQQQ,
" %2ld new drift velocity:%10.2e momentum absorbed:%10.2e\n",
5079 double GrnVryDpth_v;
5107 ASSERT( GrnVryDpth_v > 0. && GrnVryDpth_v <= 1.000001 );
5108 return GrnVryDpth_v;