31 #define FREE_CHECK(PTR) { ASSERT( PTR != NULL ); free( PTR ); PTR = NULL; }
32 #define FREE_SAFE(PTR) { if( PTR != NULL ) free( PTR ); PTR = NULL; }
71 static const double AC0 = 3.e-9;
72 static const double AC1G = 4.e-8;
73 static const double AC2G = 7.e-8;
76 static const double ETILDE = 2.*SQRT2/EVRYD;
89 return ELEM_CHARGE/EVRYD/
gv.
bin[nd]->Capacity;
111 if( e <=
gv.
bin[nd]->le_thres )
114 return 3.e-6*
gv.
bin[nd]->eec*
powpq(e*EVRYD*1.e-3,3,2);
173 double*,
double*,
double*,
bool);
189 STATIC void PE_init(
size_t,
long,
long,
double*,
double*,
double*,
190 double*,
double*,
double*,
double*);
255 for(
unsigned int ns=0; ns <
sd.size(); ns++ )
259 for(
int nz=0; nz <
NCHS; nz++ )
296 for(
int i=0; i < 3; ++i )
340 for(
int nz=0; nz <
NCHS; nz++ )
346 for(
size_t nd=0; nd <
bin.size(); nd++ )
350 for(
int nelem=0; nelem <
LIMELM; nelem++ )
368 for(
int nelem=0; nelem <
LIMELM; nelem++ )
464 for(
int nelem=0; nelem <
LIMELM; nelem++ )
466 for(
int ion=0; ion <= nelem+1; ion++ )
468 for(
int ion_to=0; ion_to <= nelem+1; ion_to++ )
510 for(
size_t nd=0; nd <
gv.
bin.size(); nd++ )
512 gv.
bin[nd]->ZloSave =
gv.
bin[nd]->chrg[0]->DustZ;
513 gv.
bin[nd]->qtmin = (
gv.
bin[nd]->qtmin_zone1 > 0. ) ?
514 gv.
bin[nd]->qtmin_zone1 : DBL_MAX;
515 gv.
bin[nd]->avdust = 0.;
516 gv.
bin[nd]->avdpot = 0.;
517 gv.
bin[nd]->avdft = 0.;
518 gv.
bin[nd]->avDGRatio = 0.;
519 gv.
bin[nd]->TeGrainMax = -1.f;
520 gv.
bin[nd]->lgEverQHeat =
false;
521 gv.
bin[nd]->QHeatFailures = 0L;
522 gv.
bin[nd]->lgQHTooWide =
false;
523 gv.
bin[nd]->lgPAHsInIonizedRegion =
false;
539 for(
size_t nd=0; nd <
gv.
bin.size(); nd++ )
541 gv.
bin[nd]->lgIterStart =
true;
583 for( nelem=0; nelem <
LIMELM; nelem++ )
629 for(
size_t nd=0; nd <
gv.
bin.size(); nd++ )
631 double help,
atoms,p_rad,ThresInf,ThresInfVal,Emin,d[5];
632 long low1,low2,low3,lowm;
640 fprintf(
ioQQQ,
" Grain work function for %s has insane value: %.4e\n",
641 gv.
bin[nd]->chDstLab,
gv.
bin[nd]->DustWorkFcn );
648 gv.
bin[nd]->lgQHeat =
true;
654 gv.
bin[nd]->lgQHeat =
false;
660 gv.
bin[nd]->lgQHeat =
false;
665 gv.
bin[nd]->lgQHTooWide =
false;
666 gv.
bin[nd]->qtmin = DBL_MAX;
669 gv.
bin[nd]->lgIterStart =
true;
677 gv.
bin[nd]->dstAbund = -FLT_MAX;
679 gv.
bin[nd]->GrnDpth = 1.f;
681 gv.
bin[nd]->qtmin_zone1 = -1.;
687 for(
long nz=0; nz <
NCHS; nz++ )
700 help =
gv.
bin[nd]->AvRadius*1.e7;
701 help = ceil(-(1.2*
POW2(help)+3.9*help+0.2)/1.44);
705 help =
gv.
bin[nd]->AvRadius*1.e7;
706 help = ceil(-(0.7*
POW2(help)+2.5*help+0.8)/1.44);
709 fprintf(
ioQQQ,
" GrainsInit detected unknown Zmin type: %d\n" , zcase );
714 ASSERT( help > (
double)(LONG_MIN+1) );
729 while( low3-low2 > 1 )
731 lowm = (low2+low3)/2;
746 gv.
bin[nd]->ZloSave = LONG_MIN;
747 for(
int i=0; i < 3; ++i )
749 gv.
bin[nd]->AveZUsed[i] = -DBL_MAX;
750 gv.
bin[nd]->TeUsed[i] = -DBL_MAX;
759 gv.
bin[nd]->sd[ns]->nelem = -1;
760 gv.
bin[nd]->sd[ns]->ns = -1;
761 gv.
bin[nd]->sd[ns]->ionPot =
gv.
bin[nd]->DustWorkFcn;
766 if(
gv.
bin[nd]->elmAbund[nelem] > 0. )
770 fprintf(
ioQQQ,
" Grain Auger data are missing for element %s\n",
772 fprintf(
ioQQQ,
" Please include the NO GRAIN X-RAY TREATMENT command "
773 "to disable the Auger treatment in grains.\n" );
783 gv.
bin[nd]->sd[ns]->nelem = nelem;
784 gv.
bin[nd]->sd[ns]->ns = j;
792 for( ns=0; ns <
gv.
bin[nd]->sd.size(); ns++ )
795 double Ethres = ( ns == 0 ) ? ThresInfVal :
gv.
bin[nd]->sd[ns]->ionPot;
804 sptr->p.reserve( len );
807 sptr->y01.reserve( len );
814 sptr->AvNr.resize( sptr->nData );
815 sptr->Ener.resize( sptr->nData );
816 sptr->y01A.resize( sptr->nData );
818 for(
long n=0; n < sptr->nData; n++ )
823 sptr->y01A[n].reserve( len );
842 atoms =
gv.
bin[nd]->AvVol*
gv.
bin[nd]->dustp[0]/ATOMIC_MASS_UNIT/
gv.
bin[nd]->atomWeight;
843 p_rad = 1./(1.+exp(20.-atoms));
844 gv.
bin[nd]->StickElecNeg =
gv.
bin[nd]->StickElecPos*p_rad;
854 gv.
bin[nd]->cnv_CM3_pH = 1./
gv.
bin[nd]->cnv_H_pCM3;
856 gv.
bin[nd]->cnv_CM3_pGR =
gv.
bin[nd]->cnv_H_pGR/
gv.
bin[nd]->cnv_H_pCM3;
857 gv.
bin[nd]->cnv_GR_pCM3 = 1./
gv.
bin[nd]->cnv_CM3_pGR;
865 for( nelem=0; nelem <
LIMELM; nelem++ )
870 for(
size_t nd=0; nd <
gv.
bin.size(); nd++ )
872 for( nelem=0; nelem <
LIMELM; nelem++ )
896 fprintf(
ioQQQ,
" There are %ld grain types turned on.\n", (
unsigned long)
gv.
bin.size() );
898 fprintf(
ioQQQ,
" grain depletion factors, dstfactor*GrainMetal=" );
899 for(
size_t nd=0; nd <
gv.
bin.size(); nd++ )
906 for(
size_t nd=0; nd <
gv.
bin.size(); nd++ )
913 for(
size_t nd=0; nd <
gv.
bin.size(); nd++ )
919 fprintf(
ioQQQ,
" nDustFunc flag for user requested custom depth dependence:" );
920 for(
size_t nd=0; nd <
gv.
bin.size(); nd++ )
927 for(
size_t nd=0; nd <
gv.
bin.size(); nd++ )
934 fprintf(
ioQQQ,
" NU(Ryd), Abs cross sec per proton\n" );
937 for(
size_t nd=0; nd <
gv.
bin.size(); nd++ )
946 for(
size_t nd=0; nd <
gv.
bin.size(); nd++ )
953 fprintf(
ioQQQ,
" NU(Ryd), Sct cross sec per proton\n" );
956 for(
size_t nd=0; nd <
gv.
bin.size(); nd++ )
965 for(
size_t nd=0; nd <
gv.
bin.size(); nd++ )
975 for(
size_t nd=0; nd <
gv.
bin.size(); nd++ )
984 for(
size_t nd=0; nd <
gv.
bin.size(); nd++ )
994 for(
size_t nd=0; nd <
gv.
bin.size(); nd++ )
1003 for(
size_t nd=0; nd <
gv.
bin.size(); nd++ )
1013 for(
size_t nd=0; nd <
gv.
bin.size(); nd++ )
1022 for(
size_t nd=0; nd <
gv.
bin.size(); nd++ )
1043 static const char chFile[] =
"auger_spectrum.dat";
1047 sscanf( chString,
"%ld", &version );
1050 fprintf(
ioQQQ,
" File %s has wrong version number\n", chFile );
1051 fprintf(
ioQQQ,
" please update your installation...\n" );
1062 if( sscanf( chString,
"%ld", &nelem ) != 1 )
1076 if( sscanf( chString,
"%u", &ptr->
nSubShell ) != 1 )
1093 if( sscanf( chString,
"%u", &ss ) != 1 )
1101 if( sscanf( chString,
"%le", &ptr->
IonThres[ns] ) != 1 )
1109 if( sscanf( chString,
"%u", &ptr->
nData[ns] ) != 1 )
1119 for(
unsigned int n=0; n < ptr->
nData[ns]; n++ )
1122 if( sscanf(chString,
"%le %le",&ptr->
Energy[ns][n],&ptr->
AvNumber[ns][n]) != 2 )
1127 ptr->
Energy[ns][n] /= EVRYD;
1146 long i, ipLo, nelem;
1153 double norm =
gv.
bin[nd]->cnv_H_pGR/
gv.
bin[nd]->AvVol;
1156 for( ns=0; ns <
gv.
bin[nd]->sd.size(); ns++ )
1158 ipLo =
max(
gv.
bin[nd]->sd[ns]->ipLo, ipBegin );
1160 gv.
bin[nd]->sd[ns]->p.realloc( ipEnd );
1164 for( i=ipLo; i < ipEnd; i++ )
1175 gv.
bin[nd]->sd[ns]->p[i] = 0.;
1179 if(
gv.
bin[nd]->elmAbund[nelem] == 0. )
1188 gv.
bin[nd]->sd[ns]->p[i] +=
1193 temp[i] +=
gv.
bin[nd]->sd[ns]->p[i];
1198 nelem =
gv.
bin[nd]->sd[ns]->nelem;
1200 nsh =
gv.
bin[nd]->sd[ns]->ns;
1202 gv.
bin[nd]->sd[ns]->p[i] =
1204 temp[i] +=
gv.
bin[nd]->sd[ns]->p[i];
1209 for( i=ipBegin; i < ipEnd && !
gv.
lgWD01; i++ )
1213 gv.
bin[nd]->inv_att_len[i] = temp[i];
1216 for( ns=0; ns <
gv.
bin[nd]->sd.size(); ns++ )
1218 ipLo =
max(
gv.
bin[nd]->sd[ns]->ipLo, ipBegin );
1220 for( i=ipLo; i < ipEnd; i++ )
1223 gv.
bin[nd]->sd[ns]->p[i] /= temp[i];
1225 gv.
bin[nd]->sd[ns]->p[i] = ( ns == 0 ) ? 1.f : 0.f;
1231 for( ns=0; ns <
gv.
bin[nd]->sd.size(); ns++ )
1236 ipLo =
max( sptr->
ipLo, ipBegin );
1241 for( i=ipLo; i < ipEnd; i++ )
1243 double elec_en,yzero,yone;
1246 yzero =
y0psa( nd, ns, i, elec_en );
1250 yone =
y1psa( nd, i, elec_en );
1267 for( n=0; n < sptr->
nData; n++ )
1269 sptr->
y01A[n].realloc( ipEnd );
1271 for( i=ipLo; i < ipEnd; i++ )
1273 double yzero = sptr->
AvNr[n] *
y0psa( nd, ns, i, sptr->
Ener[n] );
1277 double yone =
y1psa( nd, i, sptr->
Ener[n] );
1305 while( chLine[0] ==
'#' );
1326 fprintf(
ioQQQ,
" ND Tdust Emis BB Check 4pi*a^2*<Q>\n" );
1335 for(
size_t nd=0; nd <
gv.
bin.size(); nd++ )
1347 for(
size_t nd=0; nd <
gv.
bin.size(); nd++ )
1359 mul = sqrt((
double)(NDEMS-NTOP));
1361 mul = (double)NTOP + sqrt((
double)
NDEMS);
1366 for(
size_t nd=0; nd <
gv.
bin.size(); nd++ )
1377 for( i=0; i < 2000; i++ )
1380 z =
exp10(-40. + (
double)i/50.);
1385 printf(
" input %.6e temp %.6e output %.6e rel. diff. %.6e\n",z,exp(y),x,(x-z)/z);
1411 double TDustRyg = TE1RYD/tdust;
1412 double x = 0.999*log(DBL_MAX);
1419 if( arg[i] < 0.9999e-5 )
1421 if( arg[i] > x+0.0001 )
1427 vexp( arg.ptr0(), expval.
ptr0(), mini, maxi );
1429 double integral1 = 0.;
1430 double integral2 = 0.;
1432 for(
long i=0; i < maxi; i++ )
1436 if( arg[i] < 1.e-5 )
1439 ExpM1 = arg[i]*(1. + arg[i]/2.);
1444 ExpM1 = expval[i] - 1.;
1447 double Planck1 = PI4*2.*HPLANCK/
POW2(SPEEDLIGHT)*
POW2(FR1RYD)*
POW2(FR1RYD)*
1449 double Planck2 = Planck1*
gv.
bin[nd]->dstab1[i];
1458 if( Planck1/integral1 < DBL_EPSILON && Planck2/integral2 < DBL_EPSILON )
1461 integral1 += Planck1;
1462 integral2 += Planck2;
1468 fprintf(
ioQQQ,
" %4ld %11.4e %11.4e %11.4e %11.4e\n", (
unsigned long)nd, tdust,
1469 integral2, integral1/4./5.67051e-5/
powi(tdust,4), integral2*4./integral1 );
1472 ASSERT( integral2 > 0. );
1484 for( nz=0; nz <
NCHS; nz++ )
1486 gv.
bin[nd]->chrg[nz]->RSum1 = -DBL_MAX;
1487 gv.
bin[nd]->chrg[nz]->RSum2 = -DBL_MAX;
1488 gv.
bin[nd]->chrg[nz]->ESum1a = -DBL_MAX;
1489 gv.
bin[nd]->chrg[nz]->ESum1b = -DBL_MAX;
1490 gv.
bin[nd]->chrg[nz]->ESum2 = -DBL_MAX;
1492 gv.
bin[nd]->chrg[nz]->hots1 = -DBL_MAX;
1493 gv.
bin[nd]->chrg[nz]->bolflux1 = -DBL_MAX;
1494 gv.
bin[nd]->chrg[nz]->pe1 = -DBL_MAX;
1497 gv.
bin[nd]->chrg[nz]->ThermRate = -DBL_MAX;
1498 gv.
bin[nd]->chrg[nz]->GrainHeat = -DBL_MAX;
1499 gv.
bin[nd]->chrg[nz]->HeatingRate2 = -DBL_MAX;
1504 for( nz=0; nz <
NCHS; nz++ )
1506 memset(
gv.
bin[nd]->chrg[nz]->eta, 0, (
LIMELM+2)*
sizeof(
double) );
1507 memset(
gv.
bin[nd]->chrg[nz]->xi, 0, (
LIMELM+2)*
sizeof(
double) );
1513 for( nz=0; nz <
NCHS; nz++ )
1515 gv.
bin[nd]->chrg[nz]->hcon1 = -DBL_MAX;
1526 double GrnStdDpth_v;
1590 GrnStdDpth_v =
max(1.e-10,GrnStdDpth_v);
1592 return GrnStdDpth_v;
1604 static double tesave = -1.;
1605 static long int nzonesave = -1;
1643 long nelem, ion, ion_to;
1648 for(
size_t nd=0; nd <
gv.
bin.size(); nd++ )
1653 gv.
bin[nd]->lgPAHsInIonizedRegion =
false;
1655 for( nz=0; nz <
gv.
bin[nd]->nChrg; nz++ )
1657 gv.
bin[nd]->chrg[nz]->DustZ = nz;
1658 gv.
bin[nd]->chrg[nz]->FracPop = ( nz == 0 ) ? 1. : 0.;
1659 gv.
bin[nd]->chrg[nz]->nfill = 0;
1660 gv.
bin[nd]->chrg[nz]->tedust = 100.f;
1663 gv.
bin[nd]->AveDustZ = 0.;
1666 gv.
bin[nd]->tedust = 100.f;
1667 gv.
bin[nd]->TeGrainMax = 100.;
1670 gv.
bin[nd]->BolFlux = 0.;
1671 gv.
bin[nd]->GrainCoolTherm = 0.;
1672 gv.
bin[nd]->GasHeatPhotoEl = 0.;
1673 gv.
bin[nd]->GrainHeat = 0.;
1674 gv.
bin[nd]->GrainHeatColl = 0.;
1675 gv.
bin[nd]->ChemEn = 0.;
1676 gv.
bin[nd]->ChemEnH2 = 0.;
1677 gv.
bin[nd]->thermionic = 0.;
1679 gv.
bin[nd]->lgUseQHeat =
false;
1680 gv.
bin[nd]->lgEverQHeat =
false;
1681 gv.
bin[nd]->QHeatFailures = 0;
1683 gv.
bin[nd]->DustDftVel = 0.;
1686 gv.
bin[nd]->avdft = 0.f;
1688 gv.
bin[nd]->avDGRatio = -1.f;
1712 for( nelem=0; nelem <
LIMELM; nelem++ )
1714 for( ion=0; ion <= nelem+1; ion++ )
1716 for( ion_to=0; ion_to <= nelem+1; ion_to++ )
1770 static long int oldZone = -1;
1771 static double oldTe = -DBL_MAX,
1797 for( nelem=0; nelem <
LIMELM; nelem++ )
1799 for( ion=0; ion <= nelem+1; ion++ )
1801 for( ion_to=0; ion_to <= nelem+1; ion_to++ )
1811 for(
size_t nd=0; nd <
gv.
bin.size(); nd++ )
1814 double ChTdBracketLo = 0., ChTdBracketHi = -DBL_MAX;
1822 gv.
bin[nd]->lgPAHsInIonizedRegion =
true;
1833 fprintf(
ioQQQ,
" >>GrainChargeTemp starting grain %s\n",
1834 gv.
bin[nd]->chDstLab );
1839 for( i=0; i < relax*CT_LOOP_MAX && delta >
TOLER; ++i )
1843 double TdBracketLo = 0., TdBracketHi = -DBL_MAX;
1844 double ThresEst = 0.;
1845 oldtemp =
gv.
bin[nd]->tedust;
1860 gv.
bin[nd]->lgTdustConverged =
false;
1863 double oldTemp2 =
gv.
bin[nd]->tedust;
1864 double oldHeat2 =
gv.
bin[nd]->GrainHeat;
1865 double oldCool =
gv.
bin[nd]->GrainGasCool;
1870 gv.
bin[nd]->GrainGasCool = dccool;
1874 fprintf(
ioQQQ,
" >>loop %ld BracketLo %.6e BracketHi %.6e",
1875 j, TdBracketLo, TdBracketHi );
1885 gv.
bin[nd]->lgTdustConverged =
true;
1892 if(
gv.
bin[nd]->tedust < oldTemp2 )
1893 TdBracketHi = oldTemp2;
1895 TdBracketLo = oldTemp2;
1905 if( TdBracketHi > TdBracketLo )
1909 if( ( j >= 2 && TdBracketLo > 0. ) ||
1910 gv.
bin[nd]->tedust <= TdBracketLo ||
1911 gv.
bin[nd]->tedust >= TdBracketHi )
1913 gv.
bin[nd]->tedust = (
realnum)(0.5*(TdBracketLo + TdBracketHi));
1932 if(
gv.
bin[nd]->lgTdustConverged )
1935 if(
gv.
bin[nd]->tedust < oldtemp )
1936 ChTdBracketHi = oldtemp;
1938 ChTdBracketLo = oldtemp;
1943 double y, x = log(
gv.
bin[nd]->tedust);
1946 gv.
bin[nd]->GrainHeat = exp(y)*
gv.
bin[nd]->cnv_H_pCM3;
1948 fprintf(
ioQQQ,
" PROBLEM temperature of grain species %s (Tg=%.3eK) not converged\n",
1966 ratio =
gv.
bin[nd]->tedust/oldtemp;
1967 for( nz=0; nz <
gv.
bin[nd]->nChrg; nz++ )
1969 ThresEst +=
gv.
bin[nd]->chrg[nz]->FracPop*
gv.
bin[nd]->chrg[nz]->ThresInf;
1971 delta = ThresEst*TE1RYD/
gv.
bin[nd]->tedust*(ratio - 1.);
1973 delta = ( delta < 0.9*log(DBL_MAX) ) ?
1974 ThermRatio*fabs(
POW2(ratio)*exp(delta)-1.) : DBL_MAX;
1980 which =
"iteration";
1990 if( ChTdBracketHi > ChTdBracketLo )
1994 if( ( i >= 2 && ChTdBracketLo > 0. ) ||
1995 gv.
bin[nd]->tedust <= ChTdBracketLo ||
1996 gv.
bin[nd]->tedust >= ChTdBracketHi )
1998 gv.
bin[nd]->tedust = (
realnum)(0.5*(ChTdBracketLo + ChTdBracketHi));
2000 which =
"bisection";
2007 fprintf(
ioQQQ,
" >>GrainChargeTemp finds delta=%.4e, ", delta );
2017 fprintf(
ioQQQ,
" PROBLEM charge/temperature not converged for %s zone %.2f\n",
2060 for( nz=0; nz <
gv.
bin[nd]->nChrg; nz++ )
2062 one +=
gv.
bin[nd]->chrg[nz]->FracPop*(double)
gv.
bin[nd]->chrg[nz]->DustZ*
2063 gv.
bin[nd]->cnv_GR_pCM3;
2069 enum {DEBUG_LOC=
false};
2073 fprintf(
ioQQQ,
" DEBUG grn chr nz\t%.2f\teden\t%.3e\tnd\t%li",
2077 fprintf(
ioQQQ,
"\tne\t%.2e\tAveDustZ\t%.2e\t%.2e\t%.2e\t%.2e",
2079 gv.
bin[nd]->AveDustZ,
2080 gv.
bin[nd]->chrg[0]->FracPop,(
double)
gv.
bin[nd]->chrg[0]->DustZ,
2081 gv.
bin[nd]->cnv_GR_pCM3);
2088 fprintf(
ioQQQ,
" %s Pot %.5e Thermal %.5e GasCoolColl %.5e" ,
2089 gv.
bin[nd]->chDstLab,
gv.
bin[nd]->dstpot,
gv.
bin[nd]->GrainHeat, dccool );
2090 fprintf(
ioQQQ,
" GasPEHeat %.5e GasThermHeat %.5e ChemHeat %.5e\n\n" ,
2091 gv.
bin[nd]->GasHeatPhotoEl,
gv.
bin[nd]->thermionic,
gv.
bin[nd]->ChemEn );
2153 printf(
"wd test: proton fraction %.5e Total DustZ %.6f heating/cooling rate %.5e %.5e\n",
2162 enum {DEBUG_LOC=
false};
2168 for( nelem=0; nelem <
LIMELM; nelem++ )
2175 for( ion_to=0; ion_to <= nelem+1; ion_to++ )
2179 printf(
" %ld->%ld %.2e", ion, ion_to,
2189 fprintf(
ioQQQ,
" %.2f Grain contribution to electron density %.2e\n",
2193 for(
size_t nd=0; nd <
gv.
bin.size(); nd++ )
2196 for( nz=0; nz <
gv.
bin[nd]->nChrg; nz++ )
2198 sum +=
gv.
bin[nd]->chrg[nz]->FracPop*(double)
gv.
bin[nd]->chrg[nz]->DustZ*
2199 gv.
bin[nd]->cnv_GR_pCM3;
2206 for(
size_t nd=0; nd <
gv.
bin.size(); nd++ )
2213 for(
size_t nd=0; nd <
gv.
bin.size(); nd++ )
2251 netloss0 = -DBL_MAX,
2252 netloss1 = -DBL_MAX,
2267 for( nz=0; nz <
NCHS; nz++ )
2269 gv.
bin[nd]->chrg[nz]->FracPop = -DBL_MAX;
2320 backup =
gv.
bin[nd]->nChrg;
2323 stride0 =
gv.
bin[nd]->nChrg-1;
2326 if(
gv.
bin[nd]->lgIterStart )
2336 double anuav =
safe_div( sum2, sum1, 0. );
2347 Zlo = -long(
nint(step/4.));
2349 power = (int)(log(step)/log((
double)stride0));
2350 power =
MAX2(power,0);
2351 double xxx =
powi((
double)stride0,power);
2358 Zlo =
gv.
bin[nd]->ZloSave;
2369 double d2ZdTe2 = 0.;
2371 double y0 =
gv.
bin[nd]->AveZUsed[0], y1 =
gv.
bin[nd]->AveZUsed[1];
2372 double dZdTe = (y1 - y0)/(
x1 - x0);
2373 if(
gv.
bin[nd]->TeUsed[2] > 0. )
2375 double x2 =
gv.
bin[nd]->TeUsed[2];
2376 double y2 =
gv.
bin[nd]->AveZUsed[2];
2377 double deriv2 = (y2 - y1)/(x2 -
x1);
2378 d2ZdTe2 = 2.*(deriv2 - dZdTe)/(x2 - x0);
2382 double step2 = d2ZdTe2*
pow2(dTe)/2.;
2383 if( abs(step2) < abs(step) )
2385 power = (int)(log(0.16*fabs(step))/log((
double)stride0));
2386 power =
MAX2(power,0);
2387 double xxx =
powi((
double)stride0,power);
2389 Zlo = (long)floor(
gv.
bin[nd]->AveDustZ + 0.92*step);
2395 Zlo =
gv.
bin[nd]->chrg[0]->DustZ;
2398 Zlo =
max(Zlo,
gv.
bin[nd]->LowestZg);
2399 UpdatePot( nd, Zlo, stride, rate_up, rate_dn );
2404 netloss0 = rate_up[0] - rate_dn[0];
2405 netloss1 = rate_up[
gv.
bin[nd]->nChrg-1] - rate_dn[
gv.
bin[nd]->nChrg-1];
2407 if( netloss0*netloss1 <= 0. )
2411 Zlo += (
gv.
bin[nd]->nChrg-1)*stride;
2417 Zlo -= (
gv.
bin[nd]->nChrg-1)*stride;
2420 UpdatePot( nd, Zlo, stride, rate_up, rate_dn );
2423 if( netloss0*netloss1 > 0. ) {
2424 fprintf(
ioQQQ,
" insanity: could not bracket grain charge for %s\n",
gv.
bin[nd]->chDstLab );
2434 netloss0 = rate_up[0] - rate_dn[0];
2435 for( nz=0; nz <
gv.
bin[nd]->nChrg-1; nz++ )
2437 netloss1 = rate_up[nz+1] - rate_dn[nz+1];
2439 if( netloss0*netloss1 <= 0. )
2441 Zlo =
gv.
bin[nd]->chrg[nz]->DustZ;
2445 netloss0 = netloss1;
2447 UpdatePot( nd, Zlo, stride, rate_up, rate_dn );
2450 ASSERT( netloss0*netloss1 <= 0. );
2452 gv.
bin[nd]->nChrg = backup;
2455 loopMax = (
gv.
bin[nd]->lgIterStart ) ? 4*
gv.
bin[nd]->nChrg : 2*
gv.
bin[nd]->nChrg;
2458 for( i=0; i < loopMax; i++ )
2460 GetFracPop( nd, Zlo, rate_up, rate_dn, &newZlo );
2469 UpdatePot( nd, Zlo, 1, rate_up, rate_dn );
2473 fprintf(
ioQQQ,
" insanity: could not converge grain charge for %s\n",
gv.
bin[nd]->chDstLab );
2478 gv.
bin[nd]->AveDustZ = 0.;
2480 for( nz=0; nz <
gv.
bin[nd]->nChrg; nz++ )
2484 gv.
bin[nd]->AveDustZ +=
gv.
bin[nd]->chrg[nz]->FracPop*
gv.
bin[nd]->chrg[nz]->DustZ;
2487 csum3 +=
gv.
bin[nd]->chrg[nz]->FracPop*d[3];
2490 *ThermRatio = ( crate > 0. ) ? csum3/crate : 0.;
2492 ASSERT( *ThermRatio >= 0. );
2494 gv.
bin[nd]->lgIterStart =
false;
2498 for(
int i=2; i > 0; i-- )
2500 gv.
bin[nd]->AveZUsed[i] =
gv.
bin[nd]->AveZUsed[i-1];
2501 gv.
bin[nd]->TeUsed[i] =
gv.
bin[nd]->TeUsed[i-1];
2503 gv.
bin[nd]->AveZUsed[0] =
gv.
bin[nd]->AveDustZ;
2513 crate = csum1a = csum1b = csum2 = csum3 = 0.;
2514 for( nz=0; nz <
gv.
bin[nd]->nChrg; nz++ )
2516 crate +=
gv.
bin[nd]->chrg[nz]->FracPop*
2518 csum1a +=
gv.
bin[nd]->chrg[nz]->FracPop*d[0];
2519 csum1b +=
gv.
bin[nd]->chrg[nz]->FracPop*d[1];
2520 csum2 +=
gv.
bin[nd]->chrg[nz]->FracPop*d[2];
2521 csum3 +=
gv.
bin[nd]->chrg[nz]->FracPop*d[3];
2524 fprintf(
ioQQQ,
" ElecEm rate1a=%.4e, rate1b=%.4e, ", csum1a, csum1b );
2525 fprintf(
ioQQQ,
"rate2=%.4e, rate3=%.4e, sum=%.4e\n", csum2, csum3, crate );
2528 fprintf(
ioQQQ,
" rate1a/sum=%.4e, rate1b/sum=%.4e, ", csum1a/crate, csum1b/crate );
2529 fprintf(
ioQQQ,
"rate2/sum=%.4e, rate3/sum=%.4e\n", csum2/crate, csum3/crate );
2532 crate = csum1 = csum2 = 0.;
2533 for( nz=0; nz <
gv.
bin[nd]->nChrg; nz++ )
2536 csum1 +=
gv.
bin[nd]->chrg[nz]->FracPop*d[0];
2537 csum2 +=
gv.
bin[nd]->chrg[nz]->FracPop*d[1];
2540 fprintf(
ioQQQ,
" ElecRc rate1=%.4e, rate2=%.4e, sum=%.4e\n", csum1, csum2, crate );
2543 fprintf(
ioQQQ,
" rate1/sum=%.4e, rate2/sum=%.4e\n", csum1/crate, csum2/crate );
2547 for( nz=0; nz <
gv.
bin[nd]->nChrg; nz++ )
2550 gv.
bin[nd]->chrg[nz]->DustZ, rate_up[nz], rate_dn[nz] );
2554 fprintf(
ioQQQ,
" FracPop: fnzone %.2f nd %ld AveZg %.5e",
2555 fnzone, (
unsigned long)nd,
gv.
bin[nd]->AveDustZ );
2556 for( nz=0; nz <
gv.
bin[nd]->nChrg; nz++ )
2563 gv.
bin[nd]->chDstLab,
gv.
bin[nd]->dstpot*EVRYD );
2564 for( nz=0; nz <
gv.
bin[nd]->nChrg; nz++ )
2566 fprintf(
ioQQQ,
" Thres[%ld]: %.4f V ThresVal[%ld]: %.4f V",
2567 gv.
bin[nd]->chrg[nz]->DustZ,
gv.
bin[nd]->chrg[nz]->ThresInf*EVRYD,
2568 gv.
bin[nd]->chrg[nz]->DustZ,
gv.
bin[nd]->chrg[nz]->ThresInfVal*EVRYD );
2597 if(
gv.
bin[nd]->chrg[nz]->RSum1 >= 0. )
2599 *sum1 =
gv.
bin[nd]->chrg[nz]->RSum1;
2600 *sum2 =
gv.
bin[nd]->chrg[nz]->RSum2;
2601 rate = *sum1 + *sum2;
2609 ve = sqrt(8.*BOLTZMANN/PI/ELECTRON_MASS*
phycon.
te);
2611 Stick = (
gv.
bin[nd]->chrg[nz]->DustZ <= -1 ) ?
gv.
bin[nd]->StickElecNeg :
gv.
bin[nd]->StickElecPos;
2616 *sum1 = (
gv.
bin[nd]->chrg[nz]->DustZ >
gv.
bin[nd]->LowestZg ) ? Stick*
dense.
eden*ve*eta : 0.;
2621 #ifndef IGNORE_GRAIN_ION_COLLISIONS
2622 for( ion=0; ion <=
LIMELM; ion++ )
2624 double CollisionRateAll = 0.;
2626 for( nelem=
MAX2(ion-1,0); nelem <
LIMELM; nelem++ )
2629 gv.
bin[nd]->chrg[nz]->RecomZ0[nelem][ion] > ion )
2634 (double)(
gv.
bin[nd]->chrg[nz]->RecomZ0[nelem][ion]-ion);
2638 if( CollisionRateAll > 0. )
2644 *sum2 += CollisionRateAll*eta;
2649 rate = *sum1 + *sum2;
2652 gv.
bin[nd]->chrg[nz]->RSum1 = *sum1;
2653 gv.
bin[nd]->chrg[nz]->RSum2 = *sum2;
2655 ASSERT( *sum1 >= 0. && *sum2 >= 0. );
2675 if(
gv.
bin[nd]->chrg[nz]->ESum1a >= 0. )
2677 *sum1a =
gv.
bin[nd]->chrg[nz]->ESum1a;
2678 *sum1b =
gv.
bin[nd]->chrg[nz]->ESum1b;
2679 *sum2 =
gv.
bin[nd]->chrg[nz]->ESum2;
2681 *sum3 = 4.*
gv.
bin[nd]->chrg[nz]->ThermRate;
2682 return *sum1a + *sum1b + *sum2 + *sum3;
2706 *sum1a /=
gv.
bin[nd]->IntArea/4.;
2709 if( gptr->
DustZ <= -1 )
2718 *sum1b /=
gv.
bin[nd]->IntArea/4.;
2723 # ifndef IGNORE_GRAIN_ION_COLLISIONS
2724 for(
long ion=0; ion <=
LIMELM; ion++ )
2726 double CollisionRateAll = 0.;
2727 for(
long nelem=
MAX2(ion-1,0); nelem <
LIMELM; nelem++ )
2730 ion > gptr->
RecomZ0[nelem][ion] )
2735 (double)(ion-gptr->
RecomZ0[nelem][ion]);
2740 if( CollisionRateAll > 0. )
2747 *sum2 += CollisionRateAll*eta;
2755 *sum3 = 4.*
gv.
bin[nd]->chrg[nz]->ThermRate;
2762 gptr->
ESum2 = *sum2;
2764 ASSERT( *sum1a >= 0. && *sum1b >= 0. && *sum2 >= 0. && *sum3 >= 0. );
2765 return *sum1a + *sum1b + *sum2 + *sum3;
2785 if(
gv.
bin[nd]->chrg[nz]->eta[ind] > 0. )
2787 *eta =
gv.
bin[nd]->chrg[nz]->eta[ind];
2788 *xi =
gv.
bin[nd]->chrg[nz]->xi[ind];
2807 double nu = (double)
gv.
bin[nd]->chrg[nz]->DustZ/(
double)ion;
2808 double tau =
gv.
bin[nd]->Capacity*BOLTZMANN*
phycon.
te*1.e-7/
POW2((
double)ion*ELEM_CHARGE);
2811 *eta = (1. - nu/tau)*(1. + sqrt(2./(tau - 2.*nu)));
2812 *xi = (1. - nu/(2.*tau))*(1. + 1./sqrt(tau - nu));
2816 *eta = 1. + sqrt(PI/(2.*tau));
2817 *xi = 1. + 0.75*sqrt(PI/(2.*tau));
2821 double theta_nu =
ThetaNu(nu);
2823 double xxx = 1. + 1./sqrt(4.*tau+3.*nu);
2824 *eta =
POW2(xxx)*exp(-theta_nu/tau);
2826 *xi = (1. + nu/(2.*tau))*(1. + 1./sqrt(3./(2.*tau)+3.*nu))*exp(-theta_nu/tau);
2831 xxx = 0.25*
powpq(nu/tau,3,4)/(
powpq(nu/tau,3,4) +
powpq((25.+3.*nu)/5.,3,4)) +
2832 (1. + 0.75*sqrt(PI/(2.*tau)))/(1. + sqrt(PI/(2.*tau)));
2833 *xi = (
MIN2(xxx,1.) + theta_nu/(2.*tau))*(*eta);
2837 ASSERT( *eta >= 0. && *xi >= 0. );
2840 gv.
bin[nd]->chrg[nz]->eta[ind] = *eta;
2841 gv.
bin[nd]->chrg[nz]->xi[ind] = *xi;
2855 const double REL_TOLER = 10.*DBL_EPSILON;
2858 double xi_nu = 1. + 1./sqrt(3.*nu);
2859 double xi_nu2 =
POW2(xi_nu);
2865 double fnu = 2.*xi_nu2 - 1. - nu*xi_nu*
POW2(xi_nu2 - 1.);
2866 double dfdxi = 4.*xi_nu - nu*((5.*xi_nu2 - 6.)*xi_nu2 + 1.);
2868 xi_nu2 =
POW2(xi_nu);
2869 xxx = fabs(old-xi_nu);
2870 }
while( xxx > REL_TOLER*xi_nu );
2872 theta_nu = nu/xi_nu - 1./(2.*xi_nu2*(xi_nu2-1.));
2920 for( nz=0; nz <
gv.
bin[nd]->nChrg; nz++ )
2926 Zg = Zlo + nz*stride;
2930 for( zz=0; zz <
NCHS-1; zz++ )
2932 if(
gv.
bin[nd]->chrg[zz]->DustZ == Zg )
2941 ptr =
gv.
bin[nd]->chrg[ind];
2944 for( zz=ind-1; zz >= nz; zz-- )
2945 gv.
bin[nd]->chrg[zz+1] =
gv.
bin[nd]->chrg[zz];
2947 gv.
bin[nd]->chrg[nz] = ptr;
2949 if(
gv.
bin[nd]->chrg[nz]->DustZ != Zg )
2962 ASSERT( rate_up[nz] >= 0. && rate_dn[nz] >= 0. );
2970 BoltzFac = (-log(
CONSERV_TOL) + 8.)*BOLTZMANN/EN1RYD;
2972 for( nz=0; nz <
gv.
bin[nd]->nChrg; nz++ )
2975 HighEnergy =
MAX2(HighEnergy,
2999 double netloss[2], pop[2][
NCHU-1];
3010 for(
long i=0; i < 2; i++ )
3012 double sum = pop[i][0] = 1.;
3013 long nlim =
gv.
bin[nd]->nChrg-1;
3014 for(
long j=1; j < nlim; j++ )
3017 if( rate_dn[nz] > 10.*rate_up[nz-1]/sqrt(DBL_MAX) )
3019 pop[i][j] = pop[i][j-1]*rate_up[nz-1]/rate_dn[nz];
3024 for(
long k=0; k < j; k++ )
3032 if( pop[i][j] > sqrt(DBL_MAX) )
3034 for(
long k=0; k <= j; k++ )
3036 pop[i][k] /= DBL_MAX/10.;
3042 for(
long j=0; j < nlim; j++ )
3046 netloss[i] += pop[i][j]*(rate_up[nz] - rate_dn[nz]);
3052 if( netloss[0]*netloss[1] > 0. )
3053 *newZlo = ( netloss[1] > 0. ) ? Zlo + 1 : Zlo - 1;
3062 if(
gv.
bin[nd]->nChrg > 2 &&
3063 ( *newZlo <
gv.
bin[nd]->LowestZg ||
3064 ( *newZlo == Zlo && pop[1][
gv.
bin[nd]->nChrg-2] < DBL_EPSILON ) ) )
3066 gv.
bin[nd]->nChrg--;
3075 printf(
" fnzone %.2f nd %ld Zlo %ld newZlo %ld netloss %.4e %.4e nChrg %ld lgRedo %d\n",
3076 fnzone, nd, Zlo, *newZlo, netloss[0], netloss[1],
gv.
bin[nd]->nChrg, lgRedo );
3081 if( *newZlo <
gv.
bin[nd]->LowestZg )
3083 fprintf(
ioQQQ,
" could not converge charge state populations for %s\n",
gv.
bin[nd]->chDstLab );
3088 if( *newZlo == Zlo )
3091 double frac0 = netloss[1]/(netloss[1]-netloss[0]);
3092 double frac1 = -netloss[0]/(netloss[1]-netloss[0]);
3094 gv.
bin[nd]->chrg[0]->FracPop = frac0*pop[0][0];
3095 gv.
bin[nd]->chrg[
gv.
bin[nd]->nChrg-1]->FracPop = frac1*pop[1][
gv.
bin[nd]->nChrg-2];
3096 for(
long nz=1; nz <
gv.
bin[nd]->nChrg-1; nz++ )
3098 gv.
bin[nd]->chrg[nz]->FracPop = frac0*pop[0][nz] + frac1*pop[1][nz-1];
3102 double test1 = 0., test2 = 0., test3 = 0.;
3103 for(
long nz=0; nz <
gv.
bin[nd]->nChrg; nz++ )
3106 test1 +=
gv.
bin[nd]->chrg[nz]->FracPop;
3107 test2 +=
gv.
bin[nd]->chrg[nz]->FracPop*rate_up[nz];
3108 test3 +=
gv.
bin[nd]->chrg[nz]->FracPop*(rate_up[nz]-rate_dn[nz]);
3110 double x1 = fabs(test1-1.);
3111 double x2 = fabs(
safe_div(test3, test2, 0.) );
3112 ASSERT(
MAX2(x1,x2) < 10.*sqrt((
double)
gv.
bin[nd]->nChrg)*DBL_EPSILON );
3138 const double INV_DELTA_E = EVRYD/3.;
3139 const double CS_PDT = 1.2e-17;
3149 gv.
bin[nd]->chrg[nz]->DustZ = Zg;
3152 memset(
gv.
bin[nd]->chrg[nz]->eta, 0, (
LIMELM+2)*
sizeof(
double) );
3153 memset(
gv.
bin[nd]->chrg[nz]->xi, 0, (
LIMELM+2)*
sizeof(
double) );
3156 &
gv.
bin[nd]->chrg[nz]->ThresSurf,&
gv.
bin[nd]->chrg[nz]->ThresSurfVal,
3170 long ipLo =
gv.
bin[nd]->chrg[nz]->ipThresInfVal;
3174 long nfill =
gv.
bin[nd]->chrg[nz]->nfill;
3180 gv.
bin[nd]->chrg[nz]->yhat.reserve( len );
3181 gv.
bin[nd]->chrg[nz]->yhat_primary.reserve( len );
3182 gv.
bin[nd]->chrg[nz]->ehat.reserve( len );
3183 gv.
bin[nd]->chrg[nz]->yhat.alloc( ipLo, nfill );
3184 gv.
bin[nd]->chrg[nz]->yhat_primary.alloc( ipLo, nfill );
3185 gv.
bin[nd]->chrg[nz]->ehat.alloc( ipLo, nfill );
3189 gv.
bin[nd]->chrg[nz]->yhat.realloc( nfill );
3190 gv.
bin[nd]->chrg[nz]->yhat_primary.realloc( nfill );
3191 gv.
bin[nd]->chrg[nz]->ehat.realloc( nfill );
3199 double DustWorkFcn =
gv.
bin[nd]->DustWorkFcn;
3200 double Elo = -
gv.
bin[nd]->chrg[nz]->PotSurf;
3201 double ThresInfVal =
gv.
bin[nd]->chrg[nz]->ThresInfVal;
3202 double Emin =
gv.
bin[nd]->chrg[nz]->Emin;
3203 double Wcorr = ThresInfVal + Emin - GrainPot;
3212 long ipBegin =
max(ipLo,ipStart);
3214 y0b(nd, nz, yzero.
ptr0(), ipBegin, nfill);
3217 avx_ptr<double> Ehi(ipBegin, nfill), Ehi_band(ipBegin, nfill), Eel(ipBegin, nfill);
3219 for(
long i=ipBegin; i < nfill; i++ )
3221 Ehi[i] = anu[i] - ThresInfVal - Emin;
3222 Ehi_band[i] = Ehi[i];
3223 Eel[i] = anu[i] - DustWorkFcn;
3226 avx_ptr<realnum> yp(ipBegin, nfill), ya(ipBegin, nfill), ys(ipBegin, nfill), eyp(ipBegin, nfill),
3227 eya(ipBegin, nfill), eys(ipBegin, nfill), Yp1(ipBegin, nfill), Ys1(ipBegin, nfill),
3228 Ehp(ipBegin, nfill), Ehs(ipBegin, nfill);
3231 Yp1.ptr0(), Ys1.ptr0(), Ehp.ptr0(), Ehs.
ptr0(), ipBegin, nfill);
3233 if( nfill > ipBegin )
3235 memset( ya.data(), 0, size_t((nfill-ipBegin)*
sizeof(ya[ipBegin])) );
3236 memset( eya.data(), 0, size_t((nfill-ipBegin)*
sizeof(eya[ipBegin])) );
3240 for(
long i=ipBegin; i < nfill; i++ )
3243 eyp[i] = Yp1[i]*Ehp[i];
3245 for(
long i=ipBegin; i < nfill; i++ )
3248 eys[i] = Ys1[i]*Ehs[i];
3252 unsigned int nsmax =
gv.
bin[nd]->sd.size();
3253 for(
unsigned int ns=1; ns < nsmax; ns++ )
3255 sptr =
gv.
bin[nd]->sd[ns];
3257 long ipBeginShell =
max(ipBegin, sptr->
ipLo);
3258 double ionPot = sptr->
ionPot;
3259 for(
long i=ipBeginShell; i < nfill; i++ )
3261 Ehi[i] = Ehi_band[i] + Wcorr - ionPot;
3262 Eel[i] = anu[i] - ionPot;
3265 Yfunc(nd, nz, sptr->
y01.
ptr0(), NULL, sptr->
p.
ptr0(), Elo, Ehi.ptr0(), Eel.ptr0(),
3266 Yp1.ptr0(), Ys1.ptr0(), Ehp.ptr0(), Ehs.
ptr0(), ipBeginShell, nfill);
3269 for(
long i=ipBeginShell; i < nfill; i++ )
3272 eyp[i] += Yp1[i]*Ehp[i];
3274 for(
long i=ipBeginShell; i < nfill; i++ )
3277 eys[i] += Ys1[i]*Ehs[i];
3281 long nmax = sptr->
nData;
3282 for(
long n=0; n < nmax; n++ )
3285 double AvNr = sptr->
AvNr[n];
3286 double Ener = sptr->
Ener[n];
3287 double Ehic = Ener - GrainPot;
3289 for(
long i=ipBeginShell; i < nfill; i++ )
3290 max[i] = AvNr*sptr->
p[i];
3292 Yfunc(nd, nz, sptr->
y01A[n].ptr0(), max.
ptr0(), Elo, Ehic, Eelc,
3293 Yp1.ptr0(), Ys1.ptr0(), Ehp.ptr0(), Ehs.
ptr0(), ipBeginShell, nfill);
3296 for(
long i=ipBeginShell; i < nfill; i++ )
3299 eya[i] += Yp1[i]*Ehp[i];
3301 for(
long i=ipBeginShell; i < nfill; i++ )
3304 eys[i] += Ys1[i]*Ehs[i];
3310 for(
long i=ipBegin; i < nfill; i++ )
3312 gv.
bin[nd]->chrg[nz]->yhat[i] = (
realnum)(yp[i] + ya[i] + ys[i]);
3314 gv.
bin[nd]->chrg[nz]->ehat[i] = (
gv.
bin[nd]->chrg[nz]->yhat[i] > 0. ) ?
3315 (
realnum)((eyp[i] + eya[i] + eys[i])/
gv.
bin[nd]->chrg[nz]->yhat[i]) : 0.f;
3317 ASSERT( yp[i] <= 1.00001 );
3330 ipLo =
gv.
bin[nd]->chrg[nz]->ipThresInf;
3339 gv.
bin[nd]->chrg[nz]->cs_pdt.reserve( len );
3342 double c1 = -CS_PDT*(double)Zg;
3343 double ThresInf =
gv.
bin[nd]->chrg[nz]->ThresInf;
3344 double cnv_GR_pH =
gv.
bin[nd]->cnv_GR_pH;
3348 double x = (anu[i] - ThresInf)*INV_DELTA_E;
3349 double cs = c1*x/
POW2(1.+(1./3.)*
POW2(x));
3352 gv.
bin[nd]->chrg[nz]->cs_pdt[i] =
MAX2(cs,0.)*cnv_GR_pH;
3360 gv.
bin[nd]->chrg[nz]->fac1.reserve( len );
3361 gv.
bin[nd]->chrg[nz]->fac2.reserve( len );
3362 gv.
bin[nd]->chrg[nz]->fac1.alloc( ipLo, nfill );
3363 gv.
bin[nd]->chrg[nz]->fac2.alloc( ipLo, nfill );
3367 gv.
bin[nd]->chrg[nz]->fac1.realloc( nfill );
3368 gv.
bin[nd]->chrg[nz]->fac2.realloc( nfill );
3373 for(
long i=
max(ipLo,ipStart); i < nfill; i++ )
3375 double cs1,cs2,cs_tot,cool1,cool2,ehat1,ehat2;
3378 PE_init(nd,nz,i,&cs1,&cs2,&cs_tot,&cool1,&cool2,&ehat1,&ehat2);
3380 gv.
bin[nd]->chrg[nz]->fac1[i] = cs_tot*anu[i]-cs1*cool1-cs2*cool2;
3381 gv.
bin[nd]->chrg[nz]->fac2[i] = cs1*ehat1 + cs2*ehat2;
3383 ASSERT(
gv.
bin[nd]->chrg[nz]->fac1[i] >= 0. &&
gv.
bin[nd]->chrg[nz]->fac2[i] >= 0. );
3395 gv.
bin[nd]->chrg[nz]->FracPop = -DBL_MAX;
3397 gv.
bin[nd]->chrg[nz]->RSum1 = -DBL_MAX;
3398 gv.
bin[nd]->chrg[nz]->RSum2 = -DBL_MAX;
3399 gv.
bin[nd]->chrg[nz]->ESum1a = -DBL_MAX;
3400 gv.
bin[nd]->chrg[nz]->ESum1b = -DBL_MAX;
3401 gv.
bin[nd]->chrg[nz]->ESum2 = -DBL_MAX;
3403 gv.
bin[nd]->chrg[nz]->tedust = 1.f;
3405 gv.
bin[nd]->chrg[nz]->hcon1 = -DBL_MAX;
3406 gv.
bin[nd]->chrg[nz]->hots1 = -DBL_MAX;
3407 gv.
bin[nd]->chrg[nz]->bolflux1 = -DBL_MAX;
3408 gv.
bin[nd]->chrg[nz]->pe1 = -DBL_MAX;
3410 gv.
bin[nd]->chrg[nz]->BolFlux = -DBL_MAX;
3411 gv.
bin[nd]->chrg[nz]->GrainHeat = -DBL_MAX;
3412 gv.
bin[nd]->chrg[nz]->GrainHeatColl = -DBL_MAX;
3413 gv.
bin[nd]->chrg[nz]->GasHeatPhotoEl = -DBL_MAX;
3414 gv.
bin[nd]->chrg[nz]->GasHeatTherm = -DBL_MAX;
3415 gv.
bin[nd]->chrg[nz]->GrainCoolTherm = -DBL_MAX;
3416 gv.
bin[nd]->chrg[nz]->ChemEnIon = -DBL_MAX;
3417 gv.
bin[nd]->chrg[nz]->ChemEnH2 = -DBL_MAX;
3419 gv.
bin[nd]->chrg[nz]->HeatingRate2 = -DBL_MAX;
3422 ASSERT(
gv.
bin[nd]->chrg[nz]->ipThresInf <=
gv.
bin[nd]->chrg[nz]->ipThresInfVal );
3442 double ThermExp =
gv.
bin[nd]->chrg[nz]->ThresInf*TE1RYD/
gv.
bin[nd]->tedust;
3445 # if defined( WD_TEST2 ) || defined( IGNORE_THERMIONIC )
3446 gv.
bin[nd]->chrg[nz]->ThermRate = 0.;
3473 long Zg =
gv.
bin[nd]->chrg[nz]->DustZ;
3479 else if( pcase ==
PE_SIL )
3483 fprintf(
ioQQQ,
" Yfunc: unknown type for PE effect: %d\n" , pcase );
3489 y2pa( Elo, Ehi, Zg, Ehp, y2pr.ptr0(), ilo, ihi );
3490 y2s( Elo, Ehi, Zg, y2pr.ptr0(), Ehs, y2sec.ptr0(), ilo, ihi );
3494 for(
long i=ilo; i < ihi; ++i )
3495 y01cap[i] =
min(y0[i],maxval[i]);
3499 for(
long i=ilo; i < ihi; ++i )
3500 y01cap[i] =
min(y0[i]*y1[i],maxval[i]);
3503 for(
long i=ilo; i < ihi; ++i )
3509 Yp[i] = y2pr[i]*y01cap[i];
3514 Ys[i] = y2sec[i]*f3*y01cap[i];
3556 long Zg =
gv.
bin[nd]->chrg[nz]->DustZ;
3558 y2pa( Elo, Ehiloc.
ptr0(), Zg, Ehp, y2pr.ptr0(), ilo, ilo+1 );
3560 if( y2pr[ilo] > 0.f )
3562 y2s( Elo, Ehiloc.
ptr0(), Zg, y2pr.ptr0(), Ehs, y2sec.
ptr0(), ilo, ilo+1 );
3564 for(
long i=ilo+1; i < ihi; ++i )
3574 else if( pcase ==
PE_SIL )
3578 fprintf(
ioQQQ,
" Yfunc: unknown type for PE effect: %d\n" , pcase );
3586 for(
long i=ilo; i < ihi; ++i )
3589 Yp[i] = y2pr[ilo]*y01cap;
3590 Ys[i] = y2sec[ilo]*f3*y01cap;
3595 memset( Yp+ilo, 0,
size_t((ihi-ilo)*
sizeof(Yp[ilo])) );
3596 memset( Ys+ilo, 0,
size_t((ihi-ilo)*
sizeof(Ys[ilo])) );
3597 memset( Ehp+ilo, 0,
size_t((ihi-ilo)*
sizeof(Ehp[ilo])) );
3598 memset( Ehs+ilo, 0,
size_t((ihi-ilo)*
sizeof(Ehs[ilo])) );
3616 y0b01( nd, nz, yzero, ilo, ihi );
3619 realnum Ethres_low = 20./EVRYD;
3620 realnum Ethres_high = 50./EVRYD;
3625 long ipr2a =
max(ilo,ip20);
3626 long ipr2b =
min(ip50,ihi);
3627 long ipr3a =
max(ilo,ip50);
3630 y0b01( nd, nz, yzero, ipr1a, ipr2b );
3632 for(
int i=ipr2a; i < ipr2b; i++ )
3634 vlog( arg.ptr0(), val.ptr0(), ipr2a, ipr2b );
3635 for(
int i=ipr2a; i < ipr2b; i++ )
3636 arg[i] =
gv.
bin[nd]->y0b06[i]/yzero[i];
3637 vlog( arg.ptr0(), val2.
ptr0(), ipr2a, ipr2b );
3638 for(
int i=ipr2a; i < ipr2b; i++ )
3640 arg[i] = val2[i]*val[i]*1.0913566679f;
3641 vexp( arg.ptr0(), val.ptr0(), ipr2a, ipr2b );
3642 for(
int i=ipr2a; i < ipr2b; i++ )
3644 for(
int i=ipr3a; i < ipr3b; i++ )
3645 yzero[i] =
gv.
bin[nd]->y0b06[i];
3647 for(
int i=ilo; i < ihi; i++ )
3648 ASSERT( yzero[i] > 0.f );
3666 for(
int i=ilo; i < ihi; i++ )
3671 yzero[i] =
realnum(xv/((1./9.e-3) + (3.7e-2/9.e-3)*xv));
3676 for(
int i=ilo; i < ihi; i++ )
3679 yzero[i] =
realnum(xv/(2.+10.*xv));
3683 fprintf(
ioQQQ,
" y0b01: unknown type for PE effect: %d\n" , pcase );
3707 yzero =
gv.
bin[nd]->sd[ns]->p[i]*leola*(1. - leola*log(1.+1./leola));
3710 double x = 1./leola;
3711 yzero =
gv.
bin[nd]->sd[ns]->p[i]*(((-1./5.*x+1./4.)*x-1./3.)*x+1./2.);
3725 double bf, beta =
gv.
bin[nd]->AvRadius*
gv.
bin[nd]->inv_att_len[i];
3727 bf =
pow2(beta) - 2.*beta + 2. - 2.*exp(-beta);
3729 bf = ((1./60.*beta - 1./12.)*beta + 1./3.)*
pow3(beta);
3733 af =
pow2(alpha) - 2.*alpha + 2. - 2.*exp(-alpha);
3735 af = ((1./60.*alpha - 1./12.)*alpha + 1./3.)*
pow3(alpha);
3737 double yone =
pow2(beta/alpha)*af/bf;
3757 for(
long i=ilo; i < ihi; ++i )
3761 double x = Elo/Ehi[i];
3762 Ehp[i] = 0.5f*
realnum(Ehi[i]*(1.-2.*x)/(1.-3.*x));
3764 y2pr[i] = ( abs(x) > 1e-4 ) ?
3766 ASSERT( Ehp[i] > 0.f && Ehp[i] <=
realnum(Ehi[i]) && y2pr[i] > 0.f && y2pr[i] <= 1.f );
3777 for(
long i=ilo; i < ihi; ++i )
3781 Ehp[i] = 0.5f*
realnum(Elo+Ehi[i]);
3813 memset( arg.data(), 0, size_t((ihi-ilo)*
sizeof(arg[ilo])) );
3814 memset( arg2.data(), 0, size_t((ihi-ilo)*
sizeof(arg2[ilo])) );
3815 memset( lgVecFun.
data(), 0, size_t((ihi-ilo)*
sizeof(lgVecFun[ilo])) );
3816 long iveclo = -1, ivechi = -1;
3819 double sR0 = (1. + yl*yl);
3820 double sqR0 = sqrt(sR0);
3824 double asinh_yl = asinh(yl);
3827 for(
long i=ilo; i < ihi; ++i )
3833 vsqrt(arg2.ptr0(), val2.ptr0(), ilo, ihi);
3835 for(
long i=ilo; i < ihi; ++i )
3837 if( !
gv.
lgWD01 && Ehi[i] > 0. && y2pr[i] > 0.f )
3844 double x2 = x*x, x3 = x2*x, x4 = x3*x, x5 = x4*x;
3845 double yh2 = yh*yh, yh3 = yh2*yh, yh4 = yh3*yh, yh5 = yh4*yh;
3846 double help1 = 2.*x-yh;
3847 double help2 = (6.*x3-15.*yh*x2+12.*yh2*x-3.*yh3)/4.;
3848 double help3 = (22.*x5-95.*yh*x4+164.*yh2*x3-141.*yh3*x2+60.*yh4*x-10.*yh5)/16.;
3849 N0[i] = yh*(help1 - help2 + help3)/x2;
3851 help1 = (3.*x-yh)/3.;
3852 help2 = (15.*x3-25.*yh*x2+15.*yh2*x-3.*yh3)/20.;
3853 help3 = (1155.*x5-3325.*yh*x4+4305.*yh2*x3-2961.*yh3*x2+1050.*yh4*x-150.*yh5)/
3855 E0[i] =
ETILDE*yh2*(help1 - help2 + help3)/x2;
3859 double sqRh = val2[i];
3860 alpha[i] = sqRh/(sqRh - 1.);
3861 if( yh/sqR0 < 0.01 )
3864 double z = yh*(yh - 2.*yl)/sR0;
3865 N0[i] = ((((7./256.*z-5./128.)*z+1./16.)*z-1./8.)*z+1./2.)*z/(sqRh-1.);
3867 double yl2 = yl*yl, yl3 = yl2*yl, yl4 = yl3*yl;
3868 double help1 = yl/2.;
3869 double help2 = (2.*yl2-1.)/3.;
3870 double help3 = (6.*yl3-9.*yl)/8.;
3871 double help4 = (8.*yl4-24.*yl2+3.)/10.;
3873 E0[i] = -alpha[i]*Ehi[i]*(((help4*h+help3)*h+help2)*h+help1)*h/sqR0;
3877 N0[i] = alpha[i]*(1./sqR0 - 1./sqRh);
3883 E0[i] = alpha[i]*
ETILDE*(asinh_yl - yh/sqRh);
3886 ASSERT( N0[i] > 0. && N0[i] <= 1. );
3891 for(
long i=ilo; i < ihi; ++i )
3893 if( !
gv.
lgWD01 && Ehi[i] > 0. && y2pr[i] > 0.f )
3896 E0[i] += alpha[i]*
ETILDE*val[i];
3897 Ehs[i] =
realnum(E0[i]/N0[i]);
3910 for(
long i=ilo; i < ihi; ++i )
3916 vsqrt(arg2.ptr0(), val2.ptr0(), ilo, ihi);
3918 for(
long i=ilo; i < ihi; ++i )
3920 if( !
gv.
lgWD01 && Ehi[i] > Elo && y2pr[i] > 0.f )
3927 double sqRh = val2[i];
3928 alpha[i] = sqRh/(sqRh - 1.);
3939 Ehs[i] =
realnum(Ehi[i] - (Ehi[i]-Elo)*((-37./840.*x2 + 1./10.)*x2 + 1./3.));
3945 for(
long i=ilo; i < ihi; ++i )
3947 if( !
gv.
lgWD01 && Ehi[i] > Elo && y2pr[i] > 0.f )
3969 long Zg =
gv.
bin[nd]->chrg[nz]->DustZ;
3971 double d[5], phi_s_up[
LIMELM+1], phi_s_dn[2];
3972 phi_s_up[0] =
gv.
bin[nd]->chrg[nz]->ThresSurf;
3973 for(
long i=1; i <=
LIMELM; i++ )
3975 phi_s_dn[0] =
gv.
bin[nd]->chrg[nz]->ThresSurfInc;
3979 for(
long nelem=0; nelem <
LIMELM; nelem++ )
3983 for(
long ion=0; ion <= nelem+1; ion++ )
3986 &
gv.
bin[nd]->chrg[nz]->RecomZ0[nelem][ion],
3987 &
gv.
bin[nd]->chrg[nz]->RecomEn[nelem][ion],
3988 &
gv.
bin[nd]->chrg[nz]->ChemEn[nelem][ion]);
3998 double *ThresInfVal,
4000 double *ThresSurfVal,
4003 bool lgUseTunnelCorr)
4011 double dZg = (double)Zg;
4018 double IP_v =
gv.
bin[nd]->DustWorkFcn + dstpot - 0.5*
one_elec(nd) +
4040 fprintf(
ioQQQ,
" GetPotValues detected unknown type for ionization pot: %d\n" , pcase );
4046 IP_v =
MAX2(IP,IP_v);
4051 double help = fabs(dZg+1);
4054 if( lgUseTunnelCorr )
4065 *ThresInf = IP - *Emin;
4066 *ThresInfVal = IP_v - *Emin;
4067 *ThresSurf = *ThresInf;
4068 *ThresSurfVal = *ThresInfVal;
4074 *ThresInfVal = IP_v;
4075 *ThresSurf = *ThresInf - dstpot;
4076 *ThresSurfVal = *ThresInfVal - dstpot;
4091 const double phi_s_up[],
4092 const double phi_s_dn[],
4106 long Zg =
gv.
bin[nd]->chrg[nz]->DustZ;
4107 double phi_s = phi_s_up[0];
4115 *ChemEn -= (
realnum)(phi_s - phi_s_up[0]);
4118 phi_s = phi_s_up[save-ion];
4123 else if( ion <= nelem &&
gv.
bin[nd]->chrg[nz]->DustZ >
gv.
bin[nd]->LowestZg &&
4129 long Zg =
gv.
bin[nd]->chrg[nz]->DustZ;
4130 double phi_s = phi_s_dn[0];
4138 *ChemEn += (
realnum)(phi_s - phi_s_dn[0]);
4144 phi_s = phi_s_dn[ion-
save];
4148 }
while( ion <= nelem && Zg >
gv.
bin[nd]->LowestZg &&
4171 # ifndef IGNORE_GRAIN_ION_COLLISIONS
4173 for(
long nz=0; nz <
gv.
bin[nd]->nChrg; nz++ )
4176 double fac1 = gptr->
FracPop*fac0;
4181 for(
long ion=0; ion <=
LIMELM; ion++ )
4188 double fac2 = eta*fac1;
4193 for(
long nelem=
MAX2(0,ion-1); nelem <
LIMELM; nelem++ )
4214 for(
long nelem=0; nelem <
LIMELM; nelem++ )
4220 for(
size_t nd=0; nd <
gv.
bin.size(); nd++ )
4228 gv.
bin[nd]->cnv_CM3_pH = 1./
gv.
bin[nd]->cnv_H_pCM3;
4230 gv.
bin[nd]->cnv_CM3_pGR =
gv.
bin[nd]->cnv_H_pGR/
gv.
bin[nd]->cnv_H_pCM3;
4231 gv.
bin[nd]->cnv_GR_pCM3 = 1./
gv.
bin[nd]->cnv_CM3_pGR;
4236 for(
long nelem=0; nelem <
LIMELM; nelem++ )
4256 for(
size_t nd=0; nd <
gv.
bin.size(); nd++ )
4258 double dstAbund =
gv.
bin[nd]->dstAbund;
4276 for(
long nz=0; nz <
gv.
bin[nd]->nChrg; nz++ )
4279 if( gptr->
DustZ <= -1 )
4281 double FracAbund = gptr->
FracPop*dstAbund;
4312 fprintf(
ioQQQ,
" GrainTemperature starts for grain %s\n",
gv.
bin[nd]->chDstLab );
4329 gv.
bin[nd]->GasHeatPhotoEl = 0.;
4331 gv.
bin[nd]->GrainCoolTherm = 0.;
4332 gv.
bin[nd]->thermionic = 0.;
4337 gv.
bin[nd]->BolFlux = 0.;
4341 for(
long nz=0; nz <
gv.
bin[nd]->nChrg; nz++ )
4346 bool lgReEvaluate1 = gptr->
hcon1 < 0.;
4347 bool lgReEvaluate2 = gptr->
hots1 < 0.;
4352 double hcon1, hots1, pe1, bolflux1, hla1;
4357 gptr->
hcon1 = hcon1;
4361 hcon1 = gptr->
hcon1;
4374 if( gptr->
DustZ <= -1 )
4380 gptr->
hots1 = hots1;
4386 hots1 = gptr->
hots1;
4409 ASSERT( hcon1 >= 0. && hots1 >= 0. && hla1 >= 0. && bolflux1 >= 0. && pe1 >= 0. );
4422 gptr->
FracPop*bolflux1*EN1RYD*
gv.
bin[nd]->cnv_H_pCM3 );
4434 double EhatThermionic = 2.*BOLTZMANN*
gv.
bin[nd]->tedust +
MAX2(gptr->
PotSurf*EN1RYD,0.);
4435 gv.
bin[nd]->GrainCoolTherm += rate * (EhatThermionic + gptr->
ThresSurf*EN1RYD);
4436 gv.
bin[nd]->thermionic += rate * (EhatThermionic - gptr->
PotSurf*EN1RYD);
4440 double norm = EN1RYD*
gv.
bin[nd]->cnv_H_pCM3;
4451 gv.
bin[nd]->BolFlux *= norm;
4462 gv.
bin[nd]->GasHeatPhotoEl *= norm;
4493 gv.
bin[nd]->GrainHeat = *hcon + *hots + dcheat -
gv.
bin[nd]->GrainCoolTherm;
4496 gv.
bin[nd]->GrainHeatColl = dcheat;
4502 if(
gv.
bin[nd]->GrainHeat > 0. )
4507 double y, x = log(
MAX2(DBL_MIN,
gv.
bin[nd]->GrainHeat*
gv.
bin[nd]->cnv_CM3_pH));
4514 gv.
bin[nd]->GrainHeat = -1.;
4515 gv.
bin[nd]->tedust = -1.;
4524 double y, x = log(
gv.
bin[nd]->tedust);
4526 gv.
bin[nd]->GrainHeat = exp(y)*
gv.
bin[nd]->cnv_H_pCM3;
4534 fprintf(
ioQQQ,
" >GrainTemperature finds %s Tdst %.5e hcon %.4e ",
4535 gv.
bin[nd]->chDstLab,
gv.
bin[nd]->tedust, *hcon);
4536 fprintf(
ioQQQ,
"hots %.4e dcheat %.4e GrainCoolTherm %.4e\n",
4537 *hots, dcheat,
gv.
bin[nd]->GrainCoolTherm );
4572 *cs1 =
gv.
bin[nd]->dstab1[i]*gptr->
yhat[i];
4575 *ehat1 = gptr->
ehat[i];
4596 if( gptr->
DustZ <= -1 )
4601 ASSERT( *ehat1 > 0. && *cool1 > 0. );
4611 if( gptr->
DustZ <= -1 && i >= ipLo2 )
4620 ASSERT( *ehat2 > 0. && *cool2 > 0. );
4629 *cs_tot =
gv.
bin[nd]->dstab1[i] + *cs2;
4643 double Accommodation,
4644 CollisionRateElectr,
4671 const double H2_FORMATION_GRAIN_HEATING[
H2_TOP] = { 0.20, 0.4, 1.72 };
4691 gv.
bin[nd]->ChemEn = 0.;
4694 for( nz=0; nz <
gv.
bin[nd]->nChrg; nz++ )
4702 double ChemEn1 = 0.;
4708 for( ion=0; ion <=
LIMELM; ion++ )
4717 CollisionRateIon = 0.;
4719 CoolPotentialGas = 0.;
4720 HeatRecombination = 0.;
4727 for( nelem=
MAX2(0,ion-1); nelem <
LIMELM; nelem++ )
4731 double CollisionRateOne;
4736 #if defined( IGNORE_GRAIN_ION_COLLISIONS )
4738 #elif defined( WD_TEST2 )
4739 Stick = ( ion == gptr->
RecomZ0[nelem][ion] ) ?
4742 Stick = ( ion == gptr->
RecomZ0[nelem][ion] ) ?
4750 CollisionRateIon += CollisionRateOne;
4759 if( ion >= gptr->
RecomZ0[nelem][ion] )
4761 CoolPotential += CollisionRateOne * (double)ion *
4763 CoolPotentialGas += CollisionRateOne *
4764 (
double)gptr->
RecomZ0[nelem][ion] *
4769 CoolPotential += CollisionRateOne * (double)ion *
4771 CoolPotentialGas += CollisionRateOne *
4772 (
double)gptr->
RecomZ0[nelem][ion] *
4779 HeatRecombination += CollisionRateOne *
4781 HeatChem += CollisionRateOne * gptr->
ChemEn[nelem][ion];
4790 HeatCollisions = CollisionRateIon * 2.*BOLTZMANN*
phycon.
te*xi;
4795 CoolPotential *= eta*EN1RYD;
4796 CoolPotentialGas *= eta*EN1RYD;
4798 HeatRecombination *= eta*EN1RYD;
4799 HeatChem *= eta*EN1RYD;
4801 CoolEmitted = CollisionRateIon * 2.*BOLTZMANN*
gv.
bin[nd]->tedust*eta;
4804 Heat1 += HeatCollisions - CoolPotential + HeatRecombination - CoolEmitted;
4809 Cool1 += HeatCollisions - CoolEmitted - CoolPotentialGas;
4811 ChemEn1 += HeatChem;
4815 HeatIons += gptr->
FracPop*Heat1;
4828 Stick = ( gptr->
DustZ <= -1 ) ?
gv.
bin[nd]->StickElecNeg :
gv.
bin[nd]->StickElecPos;
4831 ve = sqrt(8.*BOLTZMANN/PI/ELECTRON_MASS*
phycon.
te);
4834 CollisionRateElectr = Stick*
dense.
eden*ve;
4842 HeatCollisions = CollisionRateElectr*2.*BOLTZMANN*
phycon.
te*xi;
4845 CoolPotential = CollisionRateElectr * (double)ion*gptr->
PotSurfInc*eta*EN1RYD;
4847 HeatRecombination = CollisionRateElectr * gptr->
ThresSurfInc*eta*EN1RYD;
4853 HeatCollisions = 0.;
4855 HeatRecombination = 0.;
4860 HeatBounce = CollisionRateElectr * 2.*BOLTZMANN*
phycon.
te*xi;
4863 CoolBounce = CollisionRateElectr *
4865 CoolBounce =
MAX2(CoolBounce,0.);
4870 HeatElectrons = HeatCollisions-CoolPotential+HeatRecombination+HeatBounce-CoolBounce;
4871 Heat1 += HeatElectrons;
4873 CoolElectrons = HeatCollisions+HeatBounce-CoolBounce;
4874 Cool1 += CoolElectrons;
4895 gv.
bin[nd]->GrainCoolTherm*
gv.
bin[nd]->cnv_CM3_pH;
4901 HeatTot += gptr->
FracPop*Heat1;
4904 CoolTot += gptr->
FracPop*Cool1;
4917 Accommodation = 2.*
gv.
bin[nd]->atomWeight*WeightMol/
POW2(
gv.
bin[nd]->atomWeight+WeightMol);
4919 #ifndef IGNORE_GRAIN_ION_COLLISIONS
4922 sqrt(8.*BOLTZMANN/PI/ATOMIC_MASS_UNIT/WeightMol*
phycon.
te);
4929 H2_FORMATION_GRAIN_HEATING[
ipH2]*EN1EV;
4931 gv.
bin[nd]->ChemEnH2 /=
gv.
bin[nd]->IntArea/4.*
gv.
bin[nd]->cnv_H_pCM3;
4933 CollisionRateMol = 0.;
4934 gv.
bin[nd]->ChemEnH2 = 0.;
4939 Accommodation = 2.*
gv.
bin[nd]->atomWeight*WeightMol/
POW2(
gv.
bin[nd]->atomWeight+WeightMol);
4940 #ifndef IGNORE_GRAIN_ION_COLLISIONS
4942 sqrt(8.*BOLTZMANN/PI/ATOMIC_MASS_UNIT/WeightMol*
phycon.
te);
4944 CollisionRateMol = 0.;
4948 HeatCollisions = CollisionRateMol * 2.*BOLTZMANN*
phycon.
te;
4949 CoolEmitted = CollisionRateMol * 2.*BOLTZMANN*
gv.
bin[nd]->tedust;
4951 HeatMolecules = HeatCollisions - CoolEmitted +
gv.
bin[nd]->ChemEnH2;
4952 HeatTot += HeatMolecules;
4955 CoolMolecules = HeatCollisions - CoolEmitted;
4956 CoolTot += CoolMolecules;
4958 gv.
bin[nd]->RateUp = 0.;
4959 gv.
bin[nd]->RateDn = 0.;
4961 for( nz=0; nz <
gv.
bin[nd]->nChrg; nz++ )
4967 gv.
bin[nd]->RateUp +=
gv.
bin[nd]->chrg[nz]->FracPop*rate_up;
4968 gv.
bin[nd]->RateDn +=
gv.
bin[nd]->chrg[nz]->FracPop*rate_dn;
4971 HeatCor += (
gv.
bin[nd]->chrg[nz]->FracPop*rate_up*
gv.
bin[nd]->chrg[nz]->ThresSurf -
4972 gv.
bin[nd]->chrg[nz]->FracPop*rate_dn*
gv.
bin[nd]->chrg[nz]->ThresSurfInc +
4973 gv.
bin[nd]->chrg[nz]->FracPop*rate_up*
gv.
bin[nd]->chrg[nz]->PotSurf -
4974 gv.
bin[nd]->chrg[nz]->FracPop*rate_dn*
gv.
bin[nd]->chrg[nz]->PotSurfInc)*EN1RYD;
4982 fprintf(
ioQQQ,
" molecules heat/cool: %.4e %.4e heatcor: %.4e\n",
4983 HeatMolecules*
gv.
bin[nd]->IntArea/4.*
gv.
bin[nd]->cnv_H_pCM3,
4984 CoolMolecules*
gv.
bin[nd]->IntArea/4.*
gv.
bin[nd]->cnv_H_pCM3,
4985 HeatCor*
gv.
bin[nd]->IntArea/4.*
gv.
bin[nd]->cnv_H_pCM3 );
4992 gv.
bin[nd]->ChemEnH2 *=
gv.
bin[nd]->IntArea/4.*
gv.
bin[nd]->cnv_H_pCM3;
5009 gv.
bin[nd]->HeatingRate1 = (HeatMolecules+HeatIons+HeatCor)*
gv.
bin[nd]->IntArea/4.;
5044 for(
size_t nd=0; nd <
gv.
bin.size(); nd++ )
5051 dmomen += help[i]*(
gv.
bin[nd]->dstab1[i] +
gv.
bin[nd]->pure_sc1[i]*
gv.
bin[nd]->asym[i]);
5054 dmomen *= EN1RYD*4./
gv.
bin[nd]->IntArea;
5072 phi2lm =
POW2(psi)*alam;
5075 for( loop = 0; loop < 50 && fabs(corr-1.) > 0.001; loop++ )
5077 vdold =
gv.
bin[nd]->DustDftVel;
5081 g0 = 1.5045*si*sqrt(1.+0.4418*si*si);
5082 g2 = si/(1.329 +
POW3(si));
5090 g0 = 1.5045*si*sqrt(1.+0.4418*si*si);
5091 g2 = si/(1.329 +
POW3(si));
5092 fdrag += fac*
dense.
eden*(g0 + phi2lm*g2);
5096 g0 = 1.5045*si*sqrt(1.+0.4418*si*si);
5101 g0 = 1.5045*si*sqrt(1.+0.4418*si*si);
5102 g2 = si/(1.329 +
POW3(si));
5108 volmom = dmomen/SPEEDLIGHT;
5112 corr = sqrt(volmom/fdrag);
5119 gv.
bin[nd]->DustDftVel = 0.;
5124 fprintf(
ioQQQ,
" %2ld new drift velocity:%10.2e momentum absorbed:%10.2e\n",
5125 loop,
gv.
bin[nd]->DustDftVel, volmom );
5145 if( nd >=
gv.
bin.size() )
5147 fprintf(
ioQQQ,
"invalid parameter for GrnVryDpth\n" );
5166 return max(1.e-10,GrnVryDpth_v);
STATIC double GrainElecRecomb1(size_t, long, double *, double *)
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
static const double ETILDE
STATIC void UpdatePot(size_t, long, long, double[], double[])
double rate_h2_form_grains_used
const int FILENAME_PATH_LENGTH_2
flex_arr< double > cs_pdt
STATIC double PlanckIntegral(double, size_t, long)
NORETURN void TotalInsanity(void)
double widflx(size_t i) const
STATIC void GrainIonColl(size_t, long, long, long, const double[], const double[], long *, realnum *, realnum *)
void vfast_asinh(const double x[], double y[], long nlo, long nhi)
STATIC void GetPotValues(size_t, long, double *, double *, double *, double *, double *, double *, bool)
STATIC void GrainChargeTemp()
static const bool NO_TUNNEL
long int IonHigh[LIMELM+1]
AEInfo * AugerData[LIMELM]
double pot2chrg(double x, long nd)
vector< realnum > inv_att_len
STATIC void GrainChrgTransferRates(long)
STATIC void Yfunc(long, long, const realnum[], const realnum[], const realnum[], double, const double[], const double[], realnum[], realnum[], realnum[], realnum[], long, long)
STATIC void InitEmissivities()
flex_arr< realnum > yhat_primary
const char * strstr_s(const char *haystack, const char *needle)
pe_type which_pe[MAT_TOP]
STATIC void y0b01(size_t, long, realnum[], long, long)
sys_float sexp(sys_float x)
double fudge(long int ipnt)
realnum elmSumAbund[LIMELM]
molezone * findspecieslocal(const char buf[])
static const double THERMCONST
void vlog(const double x[], double y[], long nlo, long nhi)
double anu(size_t i) const
vector< vector< double > > AvNumber
STATIC void UpdateRecomZ0(size_t, long)
double elec_esc_length(double e, long nd)
void splint(const double xa[], const double ya[], const double y2a[], long int n, double x, double *y)
vector< realnum > GraphiteEmission
H2_type which_H2distr[MAT_TOP]
t_elementnames elementnames
t_iso_sp iso_sp[NISO][LIMELM]
void vexp(const double x[], double y[], long nlo, long nhi)
long int nflux_with_check
static long int nCalledGrainDrive
double xIonDense[LIMELM][LIMELM+1]
vector< vector< double > > Energy
static const long T_LOOP_MAX
const double * anuptr() const
vector< double > dstab1_x_anu
bool fp_equal(sys_float x, sys_float y, int n=3)
STATIC double y0psa(size_t, long, long, double)
size_t ipointC(double anu) const
long ipoint(double energy_ryd)
realnum dstAbundThresholdNear
void spline(const double x[], const double y[], long int n, double yp1, double ypn, double y2a[])
void y2pa(double, const double[], long, realnum[], realnum[], long, long)
long int IonLow[LIMELM+1]
double chrg2pot(double x, long nd)
long int nsShells[LIMELM][LIMELM]
const bool ENABLE_QUANTUM_HEATING
vector< realnum > GrainEmission
strg_type which_strg[MAT_TOP]
pot_type which_pot[MAT_TOP]
realnum dstAbundThresholdFar
STATIC void GrainTemperature(size_t, realnum *, double *, double *, double *)
enth_type which_enth[MAT_TOP]
static const long MAGIC_AUGER_DATA
STATIC void GrainCollHeating(size_t, realnum *, realnum *)
static const double TOLER
STATIC void GrainUpdateRadius1()
STATIC void UpdatePot2(size_t, long)
double rate_h2_form_grains_CT02
double powi(double, long int)
STATIC void InitBinAugerData(size_t, long, long)
STATIC void GrainScreen(long, size_t, long, double *, double *)
double rate_h2_form_grains_ELRD
STATIC void GetFracPop(size_t, long, double[], double[], long *)
sys_float safe_div(sys_float x, sys_float y, sys_float res_0by0)
static const bool INCL_TUNNEL
realnum GetAveVelocity(realnum massAMU)
double anu3(size_t i) const
realnum HeatCoolRelErrorAllowed
void ConvFail(const char chMode[], const char chDetail[])
double rate_h2_form_grains_HM79
double heating(long nelem, long ion)
static const long CT_LOOP_MAX
TransitionProxy trans(const long ipHi, const long ipLo)
char chElementSym[LIMELM][CHARS_ELEMENT_SYM]
STATIC void ReadAugerData()
void setHeating(long nelem, long ion, double heating)
realnum gas_phase[LIMELM]
STATIC double ThetaNu(double)
realnum AtomicWeight[LIMELM]
ial_type which_ial[MAT_TOP]
static const double STICK_ELEC
vector< flex_arr< realnum > > y01A
double reduce_ab(const double *a, const double *b, long ilo, long ihi)
STATIC void NewChargeData(long)
void y2s(double, const double[], long, const realnum[], realnum[], realnum[], long, long)
zmin_type which_zmin[MAT_TOP]
double reduce_ab_a(const double *a, const double *b, long ilo, long ihi, double *sum_a)
realnum ChemEn[LIMELM][LIMELM+1]
void vsqrt(const double x[], double y[], long nlo, long nhi)
void realloc(size_type end)
vector< double > pure_sc1
STATIC void GrainUpdateRadius2()
STATIC void y0b(size_t, long, realnum[], long, long)
#define DEBUG_ENTRY(funcname)
double powpq(double x, int p, int q)
STATIC double GrnVryDpth(size_t)
static const long BRACKET_MAX
realnum RecomEn[LIMELM][LIMELM+1]
realnum GrainHeatScaleFactor
long RecomZ0[LIMELM][LIMELM+1]
vector< double > IonThres
STATIC void GetNextLine(const char *, FILE *, char[])
void splint_safe(const double xa[], const double ya[], const double y2a[], long int n, double x, double *y, bool *lgOutOfBounds)
double reduce_abc(const double *a, const double *b, const double *c, long ilo, long ihi)
size_t ithreshC(double threshold) const
STATIC double GrnStdDpth(long)
realnum UV_Cont_rel2_Habing_TH85_depth
int fprintf(const Output &stream, const char *format,...)
static const double INV_ETILDE
sys_float SDIV(sys_float x)
double pow(double x, int i)
void setConvIonizFail(const char *reason, double oldval, double newval)
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
STATIC double y1psa(size_t, long, double)
STATIC double GrainElecEmis1(size_t, long, double *, double *, double *, double *)
vector< unsigned int > nData
char chElementName[LIMELM][CHARS_ELEMENT_NAME]
vector< realnum > SilicateEmission
STATIC void GrainCharge(size_t, double *)
vector< string > ReadRecord
STATIC void UpdatePot1(size_t, long, long, long)
STATIC void PE_init(size_t, long, long, double *, double *, double *, double *, double *, double *, double *)
void SetNChrgStates(long nChrg)
realnum GrainChTrRate[LIMELM][LIMELM+1][LIMELM+1]
bool lgPAHsInIonizedRegion
double phfit(long int nz, long int ne, long int is, double e)
long int ipHeavy[LIMELM][LIMELM]
static double HEAT_TOLER_BIN
static const double STICK_ION