37 return ((0.075*
pow2(x) - 1./6.)*
pow2(x) + 1.)*x;
38 else if( abs(x) <= 1./sqrt(DBL_EPSILON) )
41 return -log(sqrt(1. +
pow2(x)) - x);
43 return log(sqrt(1. +
pow2(x)) + x);
55#define FREE_CHECK(PTR) { ASSERT( PTR != NULL ); free( PTR ); PTR = NULL; }
56#define FREE_SAFE(PTR) { if( PTR != NULL ) free( PTR ); PTR = NULL; }
98static const double AC0 = 3.e-9;
99static const double AC1G = 4.e-8;
100static const double AC2G = 7.e-8;
137 if( e <=
gv.bin[nd]->le_thres )
140 return 3.e-6*
gv.bin[nd]->eec*sqrt(
pow3(e*
EVRYD*1.e-3));
179inline void Yfunc(
long,
long,
double,
double,
double,
double,
double,
double*,
double*,
190inline double y2pa(
double,
double,
long,
double*);
192inline double y2s(
double,
double,
long,
double*);
199 double*,
double*,
double*,
bool);
216 double*,
double*,
double*,
double*);
280 for(
unsigned int ns=0; ns <
sd.size(); ns++ )
284 for(
int nz=0; nz <
NCHS; nz++ )
362 for(
int nz=0; nz <
NCHS; nz++ )
368 for(
size_t nd=0; nd <
bin.size(); nd++ )
372 for(
int nelem=0; nelem <
LIMELM; nelem++ )
392 for(
int nelem=0; nelem <
LIMELM; nelem++ )
480 for(
int nelem=0; nelem <
LIMELM; nelem++ )
482 for(
int ion=0; ion <= nelem+1; ion++ )
484 for(
int ion_to=0; ion_to <= nelem+1; ion_to++ )
517 if(
gv.lgDustOn() &&
gv.lgGrainPhysicsOn )
519 gv.lgNegGrnDrg =
false;
520 gv.TotalDustHeat = 0.;
521 gv.GrnElecDonateMax = 0.;
522 gv.GrnElecHoldMax = 0.;
526 for(
size_t nd=0; nd <
gv.bin.size(); nd++ )
530 gv.bin[nd]->dstpotsav =
gv.bin[nd]->dstpot;
531 gv.bin[nd]->qtmin = (
gv.bin[nd]->qtmin_zone1 > 0. ) ?
532 gv.bin[nd]->qtmin_zone1 : DBL_MAX;
533 gv.bin[nd]->avdust = 0.;
534 gv.bin[nd]->avdpot = 0.;
535 gv.bin[nd]->avdft = 0.;
536 gv.bin[nd]->avDGRatio = 0.;
537 gv.bin[nd]->TeGrainMax = -1.f;
538 gv.bin[nd]->lgEverQHeat =
false;
539 gv.bin[nd]->QHeatFailures = 0L;
540 gv.bin[nd]->lgQHTooWide =
false;
541 gv.bin[nd]->lgPAHsInIonizedRegion =
false;
542 gv.bin[nd]->nChrgOrg =
gv.bin[nd]->nChrg;
555 if(
gv.lgDustOn() &&
gv.lgGrainPhysicsOn )
557 for(
size_t nd=0; nd <
gv.bin.size(); nd++ )
561 gv.bin[nd]->dstpot =
gv.bin[nd]->dstpotsav;
562 gv.bin[nd]->nChrg =
gv.bin[nd]->nChrgOrg;
575 gv.nChrgRequested = nChrg;
593 fprintf(
ioQQQ,
" GrainsInit called.\n" );
600 gv.GrainEmission.resize(
rfield.nupper );
601 gv.GraphiteEmission.resize(
rfield.nupper );
602 gv.SilicateEmission.resize(
rfield.nupper );
605 for( nelem=0; nelem <
LIMELM; nelem++ )
607 gv.elmSumAbund[nelem] = 0.f;
610 for( i=0; i <
rfield.nupper; i++ )
615 gv.GrainEmission[i] = 0.;
616 gv.SilicateEmission[i] = 0.;
617 gv.GraphiteEmission[i] = 0.;
623 gv.GrainHeatInc = 0.;
624 gv.GrainHeatDif = 0.;
625 gv.GrainHeatLya = 0.;
626 gv.GrainHeatCollSum = 0.;
627 gv.GrainHeatSum = 0.;
635 fprintf(
ioQQQ,
" GrainsInit exits.\n" );
650 for( i=1; i <
rfield.nupper; i++ )
652 gv.anumax[
rfield.nupper-1] = FLT_MAX;
656 for(
size_t nd=0; nd <
gv.bin.size(); nd++ )
658 double help,
atoms,p_rad,ThresInf,ThresInfVal,Emin,d[5];
659 long low1,low2,low3,lowm;
667 fprintf(
ioQQQ,
" Grain work function for %s has insane value: %.4e\n",
668 gv.bin[nd]->chDstLab,
gv.bin[nd]->DustWorkFcn );
675 gv.bin[nd]->lgQHeat =
true;
681 gv.bin[nd]->lgQHeat =
false;
685 if(
thermal.ConstGrainTemp > 0. )
687 gv.bin[nd]->lgQHeat =
false;
690#ifndef IGNORE_QUANTUM_HEATING
691 gv.bin[nd]->lgQHTooWide =
false;
692 gv.bin[nd]->qtmin = DBL_MAX;
697 gv.lgAnyDustVary =
true;
701 gv.bin[nd]->dstAbund = -FLT_MAX;
703 gv.bin[nd]->GrnDpth = 1.f;
705 gv.bin[nd]->qtmin_zone1 = -1.;
708 gv.bin[nd]->le_thres =
gv.lgWD01 ? FLT_MAX :
709 (
realnum)(pow(pow((
double)
gv.bin[nd]->dustp[0],0.85)/30.,2./3.)*1.e3/
EVRYD);
711 for(
long nz=0; nz <
NCHS; nz++ )
713 ASSERT(
gv.bin[nd]->chrg[nz] == NULL );
724 help =
gv.bin[nd]->AvRadius*1.e7;
725 help = ceil(-(1.2*
POW2(help)+3.9*help+0.2)/1.44);
729 help =
gv.bin[nd]->AvRadius*1.e7;
730 help = ceil(-(0.7*
POW2(help)+2.5*help+0.8)/1.44);
733 fprintf(
ioQQQ,
" GrainsInit detected unknown Zmin type: %d\n" , zcase );
738 ASSERT( help > (
double)(LONG_MIN+1) );
753 while( low3-low2 > 1 )
755 lowm = (low2+low3)/2;
768 gv.bin[nd]->LowestZg =
MAX2(low1,low2);
769 gv.bin[nd]->LowestPot =
chrg2pot(
gv.bin[nd]->LowestZg,nd);
773 ASSERT(
gv.bin[nd]->sd.size() == 0 );
777 gv.bin[nd]->sd[ns]->nelem = -1;
778 gv.bin[nd]->sd[ns]->ns = -1;
779 gv.bin[nd]->sd[ns]->ionPot =
gv.bin[nd]->DustWorkFcn;
784 if(
gv.bin[nd]->elmAbund[nelem] > 0. )
786 if(
gv.AugerData[nelem] == NULL )
788 fprintf(
ioQQQ,
" Grain Auger data are missing for element %s\n",
790 fprintf(
ioQQQ,
" Please include the NO GRAIN X-RAY TREATMENT command "
791 "to disable the Auger treatment in grains.\n" );
795 for(
unsigned int j=0; j <
gv.AugerData[nelem]->nSubShell; j++ )
801 gv.bin[nd]->sd[ns]->nelem = nelem;
802 gv.bin[nd]->sd[ns]->ns = j;
803 gv.bin[nd]->sd[ns]->ionPot =
gv.AugerData[nelem]->IonThres[j];
810 for( ns=0; ns <
gv.bin[nd]->sd.size(); ns++ )
813 double Ethres = ( ns == 0 ) ? ThresInfVal :
gv.bin[nd]->sd[ns]->ionPot;
821 long len =
rfield.nflux + 10 - ipLo;
837 for(
long n=0; n < sptr->
nData; n++ )
839 sptr->
AvNr[n] =
gv.AugerData[sptr->
nelem]->AvNumber[sptr->
ns][n];
840 sptr->
Ener[n] =
gv.AugerData[sptr->
nelem]->Energy[sptr->
ns][n];
842 sptr->
y01A[n].reserve( len );
848 gv.bin[nd]->y0b06.resize(
rfield.nupper );
862 p_rad = 1./(1.+exp(20.-
atoms));
863 gv.bin[nd]->StickElecNeg =
gv.bin[nd]->StickElecPos*p_rad;
869 gv.bin[nd]->dstAbund = (
realnum)(
gv.bin[nd]->dstfactor*
gv.GrainMetal*
gv.bin[nd]->GrnDpth);
870 ASSERT(
gv.bin[nd]->dstAbund > 0.f );
873 gv.bin[nd]->cnv_CM3_pH = 1./
gv.bin[nd]->cnv_H_pCM3;
875 gv.bin[nd]->cnv_CM3_pGR =
gv.bin[nd]->cnv_H_pGR/
gv.bin[nd]->cnv_H_pCM3;
876 gv.bin[nd]->cnv_GR_pCM3 = 1./
gv.bin[nd]->cnv_CM3_pGR;
884 for( nelem=0; nelem <
LIMELM; nelem++ )
886 gv.elmSumAbund[nelem] = 0.f;
889 for(
size_t nd=0; nd <
gv.bin.size(); nd++ )
891 for( nelem=0; nelem <
LIMELM; nelem++ )
893 gv.elmSumAbund[nelem] +=
gv.bin[nd]->elmAbund[nelem]*(
realnum)
gv.bin[nd]->cnv_H_pCM3;
902 for( i=0; i <
rfield.nupper; i++ )
906 gv.dstab[i] = -DBL_MAX;
907 gv.dstsc[i] = -DBL_MAX;
915 fprintf(
ioQQQ,
" There are %ld grain types turned on.\n", (
unsigned long)
gv.bin.size() );
917 fprintf(
ioQQQ,
" grain depletion factors, dstfactor*GrainMetal=" );
918 for(
size_t nd=0; nd <
gv.bin.size(); nd++ )
920 fprintf(
ioQQQ,
"%10.2e",
gv.bin[nd]->dstfactor*
gv.GrainMetal );
922 fprintf(
ioQQQ,
"\n" );
924 fprintf(
ioQQQ,
" nChrg =" );
925 for(
size_t nd=0; nd <
gv.bin.size(); nd++ )
927 fprintf(
ioQQQ,
" %ld",
gv.bin[nd]->nChrg );
929 fprintf(
ioQQQ,
"\n" );
931 fprintf(
ioQQQ,
" lowest charge (e) =" );
932 for(
size_t nd=0; nd <
gv.bin.size(); nd++ )
936 fprintf(
ioQQQ,
"\n" );
938 fprintf(
ioQQQ,
" nDustFunc flag for user requested custom depth dependence:" );
939 for(
size_t nd=0; nd <
gv.bin.size(); nd++ )
941 fprintf(
ioQQQ,
"%2i",
gv.bin[nd]->nDustFunc );
943 fprintf(
ioQQQ,
"\n" );
945 fprintf(
ioQQQ,
" Quantum heating flag:" );
946 for(
size_t nd=0; nd <
gv.bin.size(); nd++ )
948 fprintf(
ioQQQ,
"%2c",
TorF(
gv.bin[nd]->lgQHeat) );
950 fprintf(
ioQQQ,
"\n" );
953 fprintf(
ioQQQ,
" NU(Ryd), Abs cross sec per proton\n" );
955 fprintf(
ioQQQ,
" Ryd " );
956 for(
size_t nd=0; nd <
gv.bin.size(); nd++ )
958 fprintf(
ioQQQ,
" %-12.12s",
gv.bin[nd]->chDstLab );
960 fprintf(
ioQQQ,
"\n" );
962 for( i=0; i <
rfield.nupper; i += 40 )
965 for(
size_t nd=0; nd <
gv.bin.size(); nd++ )
967 fprintf(
ioQQQ,
" %10.2e ",
gv.bin[nd]->dstab1[i] );
969 fprintf(
ioQQQ,
"\n" );
972 fprintf(
ioQQQ,
" NU(Ryd), Sct cross sec per proton\n" );
974 fprintf(
ioQQQ,
" Ryd " );
975 for(
size_t nd=0; nd <
gv.bin.size(); nd++ )
977 fprintf(
ioQQQ,
" %-12.12s",
gv.bin[nd]->chDstLab );
979 fprintf(
ioQQQ,
"\n" );
981 for( i=0; i <
rfield.nupper; i += 40 )
984 for(
size_t nd=0; nd <
gv.bin.size(); nd++ )
986 fprintf(
ioQQQ,
" %10.2e ",
gv.bin[nd]->pure_sc1[i] );
988 fprintf(
ioQQQ,
"\n" );
991 fprintf(
ioQQQ,
" NU(Ryd), Q abs\n" );
993 fprintf(
ioQQQ,
" Ryd " );
994 for(
size_t nd=0; nd <
gv.bin.size(); nd++ )
996 fprintf(
ioQQQ,
" %-12.12s",
gv.bin[nd]->chDstLab );
998 fprintf(
ioQQQ,
"\n" );
1000 for( i=0; i <
rfield.nupper; i += 40 )
1003 for(
size_t nd=0; nd <
gv.bin.size(); nd++ )
1005 fprintf(
ioQQQ,
" %10.2e ",
gv.bin[nd]->dstab1[i]*4./
gv.bin[nd]->IntArea );
1007 fprintf(
ioQQQ,
"\n" );
1010 fprintf(
ioQQQ,
" NU(Ryd), Q sct\n" );
1012 fprintf(
ioQQQ,
" Ryd " );
1013 for(
size_t nd=0; nd <
gv.bin.size(); nd++ )
1015 fprintf(
ioQQQ,
" %-12.12s",
gv.bin[nd]->chDstLab );
1017 fprintf(
ioQQQ,
"\n" );
1019 for( i=0; i <
rfield.nupper; i += 40 )
1022 for(
size_t nd=0; nd <
gv.bin.size(); nd++ )
1024 fprintf(
ioQQQ,
" %10.2e ",
gv.bin[nd]->pure_sc1[i]*4./
gv.bin[nd]->IntArea );
1026 fprintf(
ioQQQ,
"\n" );
1029 fprintf(
ioQQQ,
" NU(Ryd), asymmetry factor\n" );
1031 fprintf(
ioQQQ,
" Ryd " );
1032 for(
size_t nd=0; nd <
gv.bin.size(); nd++ )
1034 fprintf(
ioQQQ,
" %-12.12s",
gv.bin[nd]->chDstLab );
1036 fprintf(
ioQQQ,
"\n" );
1038 for( i=0; i <
rfield.nupper; i += 40 )
1041 for(
size_t nd=0; nd <
gv.bin.size(); nd++ )
1043 fprintf(
ioQQQ,
" %10.2e ",
gv.bin[nd]->asym[i] );
1045 fprintf(
ioQQQ,
"\n" );
1048 fprintf(
ioQQQ,
" GrainsInit exits.\n" );
1062 static const char chFile[] =
"auger_spectrum.dat";
1066 sscanf( chString,
"%ld", &version );
1069 fprintf(
ioQQQ,
" File %s has wrong version number\n", chFile );
1070 fprintf(
ioQQQ,
" please update your installation...\n" );
1082 res = sscanf( chString,
"%ld", &nelem );
1093 res = sscanf( chString,
"%u", &ptr->
nSubShell );
1106 res = sscanf( chString,
"%u", &ss );
1107 ASSERT( res == 1 && ns == ss );
1110 res = sscanf( chString,
"%le", &ptr->
IonThres[ns] );
1115 res = sscanf( chString,
"%u", &ptr->
nData[ns] );
1121 for(
unsigned int n=0; n < ptr->
nData[ns]; n++ )
1124 res = sscanf(chString,
"%le %le",&ptr->
Energy[ns][n],&ptr->
AvNumber[ns][n]);
1131 ASSERT(
gv.AugerData[nelem] == NULL );
1132 gv.AugerData[nelem] = ptr;
1145 long i, ipLo, nelem;
1152 double norm =
gv.bin[nd]->cnv_H_pGR/
gv.bin[nd]->AvVol;
1155 for( ns=0; ns <
gv.bin[nd]->sd.size(); ns++ )
1157 ipLo =
max(
gv.bin[nd]->sd[ns]->ipLo, ipBegin );
1159 gv.bin[nd]->sd[ns]->p.realloc( ipEnd );
1163 for( i=ipLo; i < ipEnd; i++ )
1174 gv.bin[nd]->sd[ns]->p[i] = 0.;
1178 if(
gv.bin[nd]->elmAbund[nelem] == 0. )
1181 long nshmax =
Heavy.nsShells[nelem][0];
1183 for( nsh =
gv.AugerData[nelem]->nSubShell; nsh < nshmax; nsh++ )
1187 gv.bin[nd]->sd[ns]->p[i] +=
1188 (
realnum)(norm*
gv.bin[nd]->elmAbund[nelem]*cs*1e-18);
1192 temp[i] +=
gv.bin[nd]->sd[ns]->p[i];
1197 nelem =
gv.bin[nd]->sd[ns]->nelem;
1199 nsh =
gv.bin[nd]->sd[ns]->ns;
1201 gv.bin[nd]->sd[ns]->p[i] =
1202 (
realnum)(norm*
gv.bin[nd]->elmAbund[nelem]*cs*1e-18);
1203 temp[i] +=
gv.bin[nd]->sd[ns]->p[i];
1208 for( i=ipBegin; i < ipEnd && !
gv.lgWD01; i++ )
1212 gv.bin[nd]->inv_att_len[i] = temp[i];
1215 for( ns=0; ns <
gv.bin[nd]->sd.size(); ns++ )
1217 ipLo =
max(
gv.bin[nd]->sd[ns]->ipLo, ipBegin );
1219 for( i=ipLo; i < ipEnd; i++ )
1222 gv.bin[nd]->sd[ns]->p[i] /= temp[i];
1224 gv.bin[nd]->sd[ns]->p[i] = ( ns == 0 ) ? 1.f : 0.f;
1230 for( ns=0; ns <
gv.bin[nd]->sd.size(); ns++ )
1235 ipLo =
max( sptr->
ipLo, ipBegin );
1240 for( i=ipLo; i < ipEnd; i++ )
1242 double elec_en,yzero,yone;
1245 yzero =
y0psa( nd, ns, i, elec_en );
1249 yone =
y1psa( nd, i, elec_en );
1266 for( n=0; n < sptr->
nData; n++ )
1268 sptr->
y01A[n].realloc( ipEnd );
1270 for( i=ipLo; i < ipEnd; i++ )
1272 double yzero = sptr->
AvNr[n] *
y0psa( nd, ns, i, sptr->
Ener[n] );
1276 double yone =
y1psa( nd, i, sptr->
Ener[n] );
1298 fprintf(
ioQQQ,
" Could not read from %s\n", chFile );
1300 fprintf(
ioQQQ,
" EOF reached\n");
1304 while( chLine[0] ==
'#' );
1325 fprintf(
ioQQQ,
" InitEmissivities starts\n" );
1326 fprintf(
ioQQQ,
" ND Tdust Emis BB Check 4pi*a^2*<Q>\n" );
1334 gv.dsttmp[i] = log(tdust);
1335 for(
size_t nd=0; nd <
gv.bin.size(); nd++ )
1346 gv.dsttmp[i] = log(tdust);
1347 for(
size_t nd=0; nd <
gv.bin.size(); nd++ )
1360 mul = (double)
NTOP + sqrt((
double)
NDEMS);
1364 for(
size_t nd=0; nd <
gv.bin.size(); nd++ )
1374 ASSERT( nd >= 0 && nd <
gv.bin.size() );
1375 for( i=0; i < 2000; i++ )
1378 z = pow(10.,-40. + (
double)i/50.);
1383 printf(
" input %.6e temp %.6e output %.6e rel. diff. %.6e\n",z,exp(y),x,(x-z)/z);
1422 x = 0.999*log(DBL_MAX);
1424 for( i=0; i <
rfield.nupper; i++ )
1427 arg = TDustRyg*
rfield.anu[i];
1433 ExpM1 = arg*(1. + arg/2.);
1438 ExpM1 = exp(
MIN2(x,arg)) - 1.;
1443 Planck2 = Planck1*
gv.bin[nd]->dstab1[i];
1452 if( Planck1/integral1 < DBL_EPSILON && Planck2/integral2 < DBL_EPSILON )
1455 integral1 += Planck1;
1456 integral2 += Planck2;
1460 if(
trace.lgDustBug &&
trace.lgTrace && ip%10 == 0 )
1462 fprintf(
ioQQQ,
" %4ld %11.4e %11.4e %11.4e %11.4e\n", (
unsigned long)nd, tdust,
1463 integral2, integral1/4./5.67051e-5/
powi(tdust,4), integral2*4./integral1 );
1466 ASSERT( integral2 > 0. );
1478 for( nz=0; nz <
NCHS; nz++ )
1480 gv.bin[nd]->chrg[nz]->RSum1 = -DBL_MAX;
1481 gv.bin[nd]->chrg[nz]->RSum2 = -DBL_MAX;
1482 gv.bin[nd]->chrg[nz]->ESum1a = -DBL_MAX;
1483 gv.bin[nd]->chrg[nz]->ESum1b = -DBL_MAX;
1484 gv.bin[nd]->chrg[nz]->ESum2 = -DBL_MAX;
1487 gv.bin[nd]->chrg[nz]->ThermRate = -DBL_MAX;
1488 gv.bin[nd]->chrg[nz]->HeatingRate2 = -DBL_MAX;
1489 gv.bin[nd]->chrg[nz]->GrainHeat = -DBL_MAX;
1491 gv.bin[nd]->chrg[nz]->hots1 = -DBL_MAX;
1492 gv.bin[nd]->chrg[nz]->bolflux1 = -DBL_MAX;
1493 gv.bin[nd]->chrg[nz]->pe1 = -DBL_MAX;
1498 for( nz=0; nz <
NCHS; nz++ )
1500 memset(
gv.bin[nd]->chrg[nz]->eta, 0, (
LIMELM+2)*
sizeof(
double) );
1501 memset(
gv.bin[nd]->chrg[nz]->xi, 0, (
LIMELM+2)*
sizeof(
double) );
1507 for( nz=0; nz <
NCHS; nz++ )
1509 gv.bin[nd]->chrg[nz]->hcon1 = -DBL_MAX;
1520 double GrnStdDpth_v;
1534 if(
gv.chPAH_abundance ==
"H" )
1543 else if(
gv.chPAH_abundance ==
"H,H2" )
1552 else if(
gv.chPAH_abundance ==
"CON" )
1559 fprintf(
ioQQQ,
"Invalid argument to SET PAH: %s\n",
gv.chPAH_abundance.c_str());
1577 GrnStdDpth_v =
sexp(
pow3(
gv.bin[nd]->tedust /
gv.bin[nd]->Tsublimat ) );
1584 GrnStdDpth_v =
max(1.e-10,GrnStdDpth_v);
1586 return GrnStdDpth_v;
1596 if(
gv.lgDustOn() &&
gv.lgGrainPhysicsOn )
1598 static double tesave = -1.;
1599 static long int nzonesave = -1;
1607 if(
gv.lgReevaluate ||
conv.lgSearch || nzonesave !=
nzone ||
1612 fabs(
gv.TotalEden)/
dense.eden >
conv.EdenErrorAllowed/5. ||
1620 if(
trace.nTrConvg >= 5 )
1622 fprintf(
ioQQQ,
" grain5 calling GrainChargeTemp\n");
1631 else if(
gv.lgDustOn() && !
gv.lgGrainPhysicsOn )
1637 long nelem, ion, ion_to;
1641 gv.GasHeatPhotoEl = 0.;
1642 for(
size_t nd=0; nd <
gv.bin.size(); nd++ )
1647 gv.bin[nd]->lgPAHsInIonizedRegion =
false;
1649 for( nz=0; nz <
gv.bin[nd]->nChrg; nz++ )
1651 gv.bin[nd]->chrg[nz]->DustZ = nz;
1652 gv.bin[nd]->chrg[nz]->FracPop = ( nz == 0 ) ? 1. : 0.;
1653 gv.bin[nd]->chrg[nz]->nfill = 0;
1654 gv.bin[nd]->chrg[nz]->tedust = 100.f;
1657 gv.bin[nd]->AveDustZ = 0.;
1660 gv.bin[nd]->tedust = 100.f;
1661 gv.bin[nd]->TeGrainMax = 100.;
1664 gv.bin[nd]->BolFlux = 0.;
1665 gv.bin[nd]->GrainCoolTherm = 0.;
1666 gv.bin[nd]->GasHeatPhotoEl = 0.;
1667 gv.bin[nd]->GrainHeat = 0.;
1668 gv.bin[nd]->GrainHeatColl = 0.;
1669 gv.bin[nd]->ChemEn = 0.;
1670 gv.bin[nd]->ChemEnH2 = 0.;
1671 gv.bin[nd]->thermionic = 0.;
1673 gv.bin[nd]->lgUseQHeat =
false;
1674 gv.bin[nd]->lgEverQHeat =
false;
1675 gv.bin[nd]->QHeatFailures = 0;
1677 gv.bin[nd]->DustDftVel = 0.;
1679 gv.bin[nd]->avdust =
gv.bin[nd]->tedust;
1680 gv.bin[nd]->avdft = 0.f;
1682 gv.bin[nd]->avDGRatio = -1.f;
1686 if( 0 &&
gv.lgBakesPAH_heat )
1693 gv.bin[nd]->GasHeatPhotoEl = 1.e-24*
hmi.UV_Cont_rel2_Habing_TH85_depth*
1694 dense.gas_phase[
ipHYDROGEN]*(4.87e-2/(1.0+4e-3*pow((
hmi.UV_Cont_rel2_Habing_TH85_depth*
1696 (1.+2.e-4*(
hmi.UV_Cont_rel2_Habing_TH85_depth*sqrt(
phycon.te)/
dense.eden)))/
gv.bin.size() *
1697 gv.GrainHeatScaleFactor;
1698 gv.GasHeatPhotoEl +=
gv.bin[nd]->GasHeatPhotoEl;
1703 gv.GrnElecDonateMax = 0.f;
1704 gv.GrnElecHoldMax = 0.f;
1706 for( nelem=0; nelem <
LIMELM; nelem++ )
1708 for( ion=0; ion <= nelem+1; ion++ )
1710 for( ion_to=0; ion_to <= nelem+1; ion_to++ )
1712 gv.GrainChTrRate[nelem][ion][ion_to] = 0.f;
1718 gv.GrainHeatInc = 0.;
1719 gv.GrainHeatDif = 0.;
1720 gv.GrainHeatLya = 0.;
1721 gv.GrainHeatCollSum = 0.;
1722 gv.GrainHeatSum = 0.;
1723 gv.GrainHeatChem = 0.;
1724 gv.GasCoolColl = 0.;
1725 gv.TotalDustHeat = 0.f;
1764 static long int oldZone = -1;
1765 static double oldTe = -DBL_MAX,
1772 fprintf(
ioQQQ,
"\n GrainChargeTemp called lgSearch%2c\n\n",
TorF(
conv.lgSearch) );
1775 oldTotalEden =
gv.TotalEden;
1778 gv.GrainHeatInc = 0.;
1779 gv.GrainHeatDif = 0.;
1780 gv.GrainHeatLya = 0.;
1781 gv.GrainHeatCollSum = 0.;
1782 gv.GrainHeatSum = 0.;
1783 gv.GrainHeatChem = 0.;
1785 gv.GasCoolColl = 0.;
1786 gv.GasHeatPhotoEl = 0.;
1787 gv.GasHeatTherm = 0.;
1791 for( nelem=0; nelem <
LIMELM; nelem++ )
1793 for( ion=0; ion <= nelem+1; ion++ )
1795 for( ion_to=0; ion_to <= nelem+1; ion_to++ )
1797 gv.GrainChTrRate[nelem][ion][ion_to] = 0.f;
1807 for(
size_t nd=0; nd <
gv.bin.size(); nd++ )
1810 double ChTdBracketLo = 0., ChTdBracketHi = -DBL_MAX;
1811 long relax = (
conv.lgSearch ) ? 3 : 1;
1818 gv.bin[nd]->lgPAHsInIonizedRegion =
true;
1829 fprintf(
ioQQQ,
" >>GrainChargeTemp starting grain %s\n",
1830 gv.bin[nd]->chDstLab );
1835 for( i=0; i < relax*CT_LOOP_MAX && delta >
TOLER; ++i )
1839 double TdBracketLo = 0., TdBracketHi = -DBL_MAX;
1840 double ThresEst = 0.;
1841 oldtemp =
gv.bin[nd]->tedust;
1849 ASSERT(
gv.bin[nd]->GrainHeat > 0. );
1856 gv.bin[nd]->lgTdustConverged =
false;
1859 double oldTemp2 =
gv.bin[nd]->tedust;
1860 double oldHeat2 =
gv.bin[nd]->GrainHeat;
1861 double oldCool =
gv.bin[nd]->GrainGasCool;
1866 gv.bin[nd]->GrainGasCool = dccool;
1870 fprintf(
ioQQQ,
" >>loop %ld BracketLo %.6e BracketHi %.6e",
1871 j, TdBracketLo, TdBracketHi );
1878 if( fabs(
gv.bin[nd]->GrainHeat-oldHeat2) <
HEAT_TOLER*
gv.bin[nd]->GrainHeat &&
1881 gv.bin[nd]->lgTdustConverged =
true;
1883 fprintf(
ioQQQ,
" converged\n" );
1888 if(
gv.bin[nd]->tedust < oldTemp2 )
1889 TdBracketHi = oldTemp2;
1891 TdBracketLo = oldTemp2;
1901 if( TdBracketHi > TdBracketLo )
1905 if( ( j >= 2 && TdBracketLo > 0. ) ||
1906 gv.bin[nd]->tedust <= TdBracketLo ||
1907 gv.bin[nd]->tedust >= TdBracketHi )
1909 gv.bin[nd]->tedust = (
realnum)(0.5*(TdBracketLo + TdBracketHi));
1911 fprintf(
ioQQQ,
" bisection\n" );
1916 fprintf(
ioQQQ,
" iteration\n" );
1922 fprintf(
ioQQQ,
" iteration\n" );
1928 if(
gv.bin[nd]->lgTdustConverged )
1931 if(
gv.bin[nd]->tedust < oldtemp )
1932 ChTdBracketHi = oldtemp;
1934 ChTdBracketLo = oldtemp;
1939 double y, x = log(
gv.bin[nd]->tedust);
1942 gv.bin[nd]->GrainHeat = exp(y)*
gv.bin[nd]->cnv_H_pCM3;
1944 fprintf(
ioQQQ,
" PROBLEM temperature of grain species %s (Tg=%.3eK) not converged\n",
1945 gv.bin[nd]->chDstLab ,
gv.bin[nd]->tedust );
1953 ASSERT(
gv.bin[nd]->GrainHeat > 0. );
1962 ratio =
gv.bin[nd]->tedust/oldtemp;
1963 for( nz=0; nz <
gv.bin[nd]->nChrg; nz++ )
1965 ThresEst +=
gv.bin[nd]->chrg[nz]->FracPop*
gv.bin[nd]->chrg[nz]->ThresInf;
1967 delta = ThresEst*
TE1RYD/
gv.bin[nd]->tedust*(ratio - 1.);
1969 delta = ( delta < 0.9*log(DBL_MAX) ) ?
1970 ThermRatio*fabs(
POW2(ratio)*exp(delta)-1.) : DBL_MAX;
1976 which =
"iteration";
1986 if( ChTdBracketHi > ChTdBracketLo )
1990 if( ( i >= 2 && ChTdBracketLo > 0. ) ||
1991 gv.bin[nd]->tedust <= ChTdBracketLo ||
1992 gv.bin[nd]->tedust >= ChTdBracketHi )
1994 gv.bin[nd]->tedust = (
realnum)(0.5*(ChTdBracketLo + ChTdBracketHi));
1996 which =
"bisection";
2003 fprintf(
ioQQQ,
" >>GrainChargeTemp finds delta=%.4e, ", delta );
2004 fprintf(
ioQQQ,
" old/new temp=%.5e %.5e, ", oldtemp,
gv.bin[nd]->tedust );
2006 fprintf(
ioQQQ,
"doing another %s\n", which.c_str() );
2008 fprintf(
ioQQQ,
"converged\n" );
2013 fprintf(
ioQQQ,
" PROBLEM charge/temperature not converged for %s zone %.2f\n",
2021 if(
ionbal.lgGrainIonRecom )
2029 gv.GrainHeatInc += hcon;
2031 gv.GrainHeatDif += hots;
2034 gv.GrainHeatLya += hla;
2037 gv.GrainHeatCollSum +=
gv.bin[nd]->GrainHeatColl;
2042 gv.GrainHeatSum +=
gv.bin[nd]->GrainHeat;
2045 gv.GrainHeatChem +=
gv.bin[nd]->ChemEn +
gv.bin[nd]->ChemEnH2;
2049 gv.GasCoolColl += dccool;
2050 gv.GasHeatPhotoEl +=
gv.bin[nd]->GasHeatPhotoEl;
2051 gv.GasHeatTherm +=
gv.bin[nd]->thermionic;
2056 for( nz=0; nz <
gv.bin[nd]->nChrg; nz++ )
2058 one +=
gv.bin[nd]->chrg[nz]->FracPop*(double)
gv.bin[nd]->chrg[nz]->DustZ*
2059 gv.bin[nd]->cnv_GR_pCM3;
2062 gv.TotalEden += one;
2065 enum {DEBUG_LOC=
false};
2069 fprintf(
ioQQQ,
" DEBUG grn chr nz\t%.2f\teden\t%.3e\tnd\t%li",
2073 fprintf(
ioQQQ,
"\tne\t%.2e\tAveDustZ\t%.2e\t%.2e\t%.2e\t%.2e",
2075 gv.bin[nd]->AveDustZ,
2076 gv.bin[nd]->chrg[0]->FracPop,(
double)
gv.bin[nd]->chrg[0]->DustZ,
2077 gv.bin[nd]->cnv_GR_pCM3);
2078 fprintf(
ioQQQ,
"\n");
2084 fprintf(
ioQQQ,
" %s Pot %.5e Thermal %.5e GasCoolColl %.5e" ,
2085 gv.bin[nd]->chDstLab,
gv.bin[nd]->dstpot,
gv.bin[nd]->GrainHeat, dccool );
2086 fprintf(
ioQQQ,
" GasPEHeat %.5e GasThermHeat %.5e ChemHeat %.5e\n\n" ,
2087 gv.bin[nd]->GasHeatPhotoEl,
gv.bin[nd]->thermionic,
gv.bin[nd]->ChemEn );
2092 GasHeatNet =
gv.GasHeatPhotoEl +
gv.GasHeatTherm -
gv.GasCoolColl;
2097 oldTe =
gv.GrnRecomTe;
2098 oldHeat =
gv.GasHeatNet;
2104 gv.dHeatdT = (GasHeatNet-oldHeat)/(
phycon.te-oldTe);
2117 conv.setConvIonizFail(
"grn eden chg" , oldTotalEden,
gv.TotalEden);
2121 conv.setConvIonizFail(
"grn het chg" ,
gv.GasHeatNet, GasHeatNet);
2125 conv.setConvIonizFail(
"grn ter chg" ,
gv.GrnRecomTe,
phycon.te);
2129 conv.setConvIonizFail(
"grn zon chg" ,
gv.nzone,
nzone );
2146 gv.GasHeatNet = GasHeatNet;
2149 printf(
"wd test: proton fraction %.5e Total DustZ %.6f heating/cooling rate %.5e %.5e\n",
2158 enum {DEBUG_LOC=
false};
2162 fprintf(
ioQQQ,
" %.2f Grain surface charge transfer rates\n",
fnzone );
2164 for( nelem=0; nelem <
LIMELM; nelem++ )
2166 if(
dense.lgElmtOn[nelem] )
2169 for( ion=
dense.IonLow[nelem]; ion <=
dense.IonHigh[nelem]; ++ion )
2171 for( ion_to=0; ion_to <= nelem+1; ion_to++ )
2173 if(
gv.GrainChTrRate[nelem][ion][ion_to] > 0.f )
2175 printf(
" %ld->%ld %.2e", ion, ion_to,
2176 gv.GrainChTrRate[nelem][ion][ion_to] );
2185 fprintf(
ioQQQ,
" %.2f Grain contribution to electron density %.2e\n",
2188 fprintf(
ioQQQ,
" Grain electons: " );
2189 for(
size_t nd=0; nd <
gv.bin.size(); nd++ )
2192 for( nz=0; nz <
gv.bin[nd]->nChrg; nz++ )
2194 sum +=
gv.bin[nd]->chrg[nz]->FracPop*(double)
gv.bin[nd]->chrg[nz]->DustZ*
2195 gv.bin[nd]->cnv_GR_pCM3;
2197 fprintf(
ioQQQ,
" %.2e", sum );
2199 fprintf(
ioQQQ,
"\n" );
2201 fprintf(
ioQQQ,
" Grain potentials:" );
2202 for(
size_t nd=0; nd <
gv.bin.size(); nd++ )
2204 fprintf(
ioQQQ,
" %.2e",
gv.bin[nd]->dstpot );
2206 fprintf(
ioQQQ,
"\n" );
2208 fprintf(
ioQQQ,
" Grain temperatures:" );
2209 for(
size_t nd=0; nd <
gv.bin.size(); nd++ )
2211 fprintf(
ioQQQ,
" %.2e",
gv.bin[nd]->tedust );
2213 fprintf(
ioQQQ,
"\n" );
2215 fprintf(
ioQQQ,
" GrainCollCool: %.6e\n",
gv.GasCoolColl );
2248 netloss0 = -DBL_MAX,
2249 netloss1 = -DBL_MAX,
2259 fprintf(
ioQQQ,
" Charge loop, search %c,",
TorF(
conv.lgSearch) );
2264 for( nz=0; nz <
NCHS; nz++ )
2266 gv.bin[nd]->chrg[nz]->FracPop = -DBL_MAX;
2317 lgInitial = (
gv.bin[nd]->chrg[0]->DustZ == LONG_MIN );
2319 backup =
gv.bin[nd]->nChrg;
2320 gv.bin[nd]->nChrg =
MAX2(
gv.bin[nd]->nChrg,3);
2322 stride0 =
gv.bin[nd]->nChrg-1;
2328 step =
MAX2((
double)(-
gv.bin[nd]->LowestZg),1.);
2329 power = (int)(log(step)/log((
double)stride0));
2330 power =
MAX2(power,0);
2331 xxx =
powi((
double)stride0,power);
2333 Zlo =
gv.bin[nd]->LowestZg;
2339 Zlo =
gv.bin[nd]->chrg[0]->DustZ;
2341 UpdatePot( nd, Zlo, stride, rate_up, rate_dn );
2346 netloss0 = rate_up[0] - rate_dn[0];
2347 netloss1 = rate_up[
gv.bin[nd]->nChrg-1] - rate_dn[
gv.bin[nd]->nChrg-1];
2349 if( netloss0*netloss1 <= 0. )
2353 Zlo += (
gv.bin[nd]->nChrg-1)*stride;
2359 Zlo -= (
gv.bin[nd]->nChrg-1)*stride;
2361 Zlo =
MAX2(Zlo,
gv.bin[nd]->LowestZg);
2362 UpdatePot( nd, Zlo, stride, rate_up, rate_dn );
2365 if( netloss0*netloss1 > 0. ) {
2366 fprintf(
ioQQQ,
" insanity: could not bracket grain charge for %s\n",
gv.bin[nd]->chDstLab );
2376 netloss0 = rate_up[0] - rate_dn[0];
2377 for( nz=0; nz <
gv.bin[nd]->nChrg-1; nz++ )
2379 netloss1 = rate_up[nz+1] - rate_dn[nz+1];
2381 if( netloss0*netloss1 <= 0. )
2383 Zlo =
gv.bin[nd]->chrg[nz]->DustZ;
2387 netloss0 = netloss1;
2389 UpdatePot( nd, Zlo, stride, rate_up, rate_dn );
2392 ASSERT( netloss0*netloss1 <= 0. );
2394 gv.bin[nd]->nChrg = backup;
2397 loopMax = ( lgInitial ) ? 4*
gv.bin[nd]->nChrg : 2*
gv.bin[nd]->nChrg;
2400 for( i=0; i < loopMax; i++ )
2402 GetFracPop( nd, Zlo, rate_up, rate_dn, &newZlo );
2411 UpdatePot( nd, Zlo, 1, rate_up, rate_dn );
2415 fprintf(
ioQQQ,
" insanity: could not converge grain charge for %s\n",
gv.bin[nd]->chDstLab );
2422 gv.bin[nd]->lgChrgConverged =
true;
2424 gv.bin[nd]->AveDustZ = 0.;
2426 for( nz=0; nz <
gv.bin[nd]->nChrg; nz++ )
2430 gv.bin[nd]->AveDustZ +=
gv.bin[nd]->chrg[nz]->FracPop*
gv.bin[nd]->chrg[nz]->DustZ;
2432 crate +=
gv.bin[nd]->chrg[nz]->FracPop*
GrainElecEmis1(nd,nz,&d[0],&d[1],&d[2],&d[3]);
2433 csum3 +=
gv.bin[nd]->chrg[nz]->FracPop*d[3];
2435 gv.bin[nd]->dstpot =
chrg2pot(
gv.bin[nd]->AveDustZ,nd);
2436 *ThermRatio = ( crate > 0. ) ? csum3/crate : 0.;
2438 ASSERT( *ThermRatio >= 0. );
2444 fprintf(
ioQQQ,
"\n" );
2446 crate = csum1a = csum1b = csum2 = csum3 = 0.;
2447 for( nz=0; nz <
gv.bin[nd]->nChrg; nz++ )
2449 crate +=
gv.bin[nd]->chrg[nz]->FracPop*
2451 csum1a +=
gv.bin[nd]->chrg[nz]->FracPop*d[0];
2452 csum1b +=
gv.bin[nd]->chrg[nz]->FracPop*d[1];
2453 csum2 +=
gv.bin[nd]->chrg[nz]->FracPop*d[2];
2454 csum3 +=
gv.bin[nd]->chrg[nz]->FracPop*d[3];
2457 fprintf(
ioQQQ,
" ElecEm rate1a=%.4e, rate1b=%.4e, ", csum1a, csum1b );
2458 fprintf(
ioQQQ,
"rate2=%.4e, rate3=%.4e, sum=%.4e\n", csum2, csum3, crate );
2461 fprintf(
ioQQQ,
" rate1a/sum=%.4e, rate1b/sum=%.4e, ", csum1a/crate, csum1b/crate );
2462 fprintf(
ioQQQ,
"rate2/sum=%.4e, rate3/sum=%.4e\n", csum2/crate, csum3/crate );
2465 crate = csum1 = csum2 = 0.;
2466 for( nz=0; nz <
gv.bin[nd]->nChrg; nz++ )
2469 csum1 +=
gv.bin[nd]->chrg[nz]->FracPop*d[0];
2470 csum2 +=
gv.bin[nd]->chrg[nz]->FracPop*d[1];
2473 fprintf(
ioQQQ,
" ElecRc rate1=%.4e, rate2=%.4e, sum=%.4e\n", csum1, csum2, crate );
2476 fprintf(
ioQQQ,
" rate1/sum=%.4e, rate2/sum=%.4e\n", csum1/crate, csum2/crate );
2479 fprintf(
ioQQQ,
" Charging rates:" );
2480 for( nz=0; nz <
gv.bin[nd]->nChrg; nz++ )
2482 fprintf(
ioQQQ,
" Zg %ld up %.4e dn %.4e",
2483 gv.bin[nd]->chrg[nz]->DustZ, rate_up[nz], rate_dn[nz] );
2485 fprintf(
ioQQQ,
"\n" );
2487 fprintf(
ioQQQ,
" FracPop: fnzone %.2f nd %ld AveZg %.5e",
2488 fnzone, (
unsigned long)nd,
gv.bin[nd]->AveDustZ );
2489 for( nz=0; nz <
gv.bin[nd]->nChrg; nz++ )
2491 fprintf(
ioQQQ,
" Zg %ld %.5f", Zlo+nz,
gv.bin[nd]->chrg[nz]->FracPop );
2493 fprintf(
ioQQQ,
"\n" );
2495 fprintf(
ioQQQ,
" >Grain potential:%12.12s %.6fV",
2496 gv.bin[nd]->chDstLab,
gv.bin[nd]->dstpot*
EVRYD );
2497 for( nz=0; nz <
gv.bin[nd]->nChrg; nz++ )
2499 fprintf(
ioQQQ,
" Thres[%ld]: %.4f V ThresVal[%ld]: %.4f V",
2500 gv.bin[nd]->chrg[nz]->DustZ,
gv.bin[nd]->chrg[nz]->ThresInf*
EVRYD,
2501 gv.bin[nd]->chrg[nz]->DustZ,
gv.bin[nd]->chrg[nz]->ThresInfVal*
EVRYD );
2503 fprintf(
ioQQQ,
"\n" );
2526 ASSERT( nz >= 0 && nz <
gv.bin[nd]->nChrg );
2530 if(
gv.bin[nd]->chrg[nz]->RSum1 >= 0. )
2532 *sum1 =
gv.bin[nd]->chrg[nz]->RSum1;
2533 *sum2 =
gv.bin[nd]->chrg[nz]->RSum2;
2534 rate = *sum1 + *sum2;
2544 Stick = (
gv.bin[nd]->chrg[nz]->DustZ <= -1 ) ?
gv.bin[nd]->StickElecNeg :
gv.bin[nd]->StickElecPos;
2549 *sum1 = (
gv.bin[nd]->chrg[nz]->DustZ >
gv.bin[nd]->LowestZg ) ? Stick*
dense.eden*ve*eta : 0.;
2554#ifndef IGNORE_GRAIN_ION_COLLISIONS
2555 for( ion=0; ion <=
LIMELM; ion++ )
2557 double CollisionRateAll = 0.;
2559 for( nelem=
MAX2(ion-1,0); nelem <
LIMELM; nelem++ )
2561 if(
dense.lgElmtOn[nelem] &&
dense.xIonDense[nelem][ion] > 0. &&
2562 gv.bin[nd]->chrg[nz]->RecomZ0[nelem][ion] > ion )
2566 (double)(
gv.bin[nd]->chrg[nz]->RecomZ0[nelem][ion]-ion);
2570 if( CollisionRateAll > 0. )
2576 *sum2 += CollisionRateAll*eta;
2581 rate = *sum1 + *sum2;
2584 gv.bin[nd]->chrg[nz]->RSum1 = *sum1;
2585 gv.bin[nd]->chrg[nz]->RSum2 = *sum2;
2587 ASSERT( *sum1 >= 0. && *sum2 >= 0. );
2611 ASSERT( nz >= 0 && nz <
gv.bin[nd]->nChrg );
2615 if(
gv.bin[nd]->chrg[nz]->ESum1a >= 0. )
2617 *sum1a =
gv.bin[nd]->chrg[nz]->ESum1a;
2618 *sum1b =
gv.bin[nd]->chrg[nz]->ESum1b;
2619 *sum2 =
gv.bin[nd]->chrg[nz]->ESum2;
2621 *sum3 = 4.*
gv.bin[nd]->chrg[nz]->ThermRate;
2622 rate = *sum1a + *sum1b + *sum2 + *sum3;
2638 ipLo =
gv.bin[nd]->chrg[nz]->ipThresInfVal;
2639 for( i=ipLo; i <
rfield.nflux; i++ )
2642 *sum1a +=
rfield.flux[0][i]*
gv.bin[nd]->dstab1[i]*
gv.bin[nd]->chrg[nz]->yhat[i];
2644 *sum1a +=
rfield.SummedCon[i]*
gv.bin[nd]->dstab1[i]*
gv.bin[nd]->chrg[nz]->yhat[i];
2648 *sum1a /=
gv.bin[nd]->IntArea/4.;
2651 if(
gv.bin[nd]->chrg[nz]->DustZ <= -1 )
2653 ipLo =
gv.bin[nd]->chrg[nz]->ipThresInf;
2654 for( i=ipLo; i <
rfield.nflux; i++ )
2658 *sum1b +=
rfield.flux[0][i]*
gv.bin[nd]->chrg[nz]->cs_pdt[i];
2660 *sum1b +=
rfield.SummedCon[i]*
gv.bin[nd]->chrg[nz]->cs_pdt[i];
2663 *sum1b /=
gv.bin[nd]->IntArea/4.;
2668# ifndef IGNORE_GRAIN_ION_COLLISIONS
2669 for( ion=0; ion <=
LIMELM; ion++ )
2671 double CollisionRateAll = 0.;
2673 for( nelem=
MAX2(ion-1,0); nelem <
LIMELM; nelem++ )
2675 if(
dense.lgElmtOn[nelem] &&
dense.xIonDense[nelem][ion] > 0. &&
2676 ion >
gv.bin[nd]->chrg[nz]->RecomZ0[nelem][ion] )
2680 (double)(ion-
gv.bin[nd]->chrg[nz]->RecomZ0[nelem][ion]);
2684 if( CollisionRateAll > 0. )
2690 *sum2 += CollisionRateAll*eta;
2698 *sum3 = 4.*
gv.bin[nd]->chrg[nz]->ThermRate;
2702 rate = *sum1a + *sum1b + *sum2 + *sum3;
2705 gv.bin[nd]->chrg[nz]->ESum1a = *sum1a;
2706 gv.bin[nd]->chrg[nz]->ESum1b = *sum1b;
2707 gv.bin[nd]->chrg[nz]->ESum2 = *sum2;
2709 ASSERT( *sum1a >= 0. && *sum1b >= 0. && *sum2 >= 0. && *sum3 >= 0. );
2732 if(
gv.bin[nd]->chrg[nz]->eta[ind] > 0. )
2734 *eta =
gv.bin[nd]->chrg[nz]->eta[ind];
2735 *xi =
gv.bin[nd]->chrg[nz]->xi[ind];
2754 double nu = (double)
gv.bin[nd]->chrg[nz]->DustZ/(
double)ion;
2758 *eta = (1. - nu/tau)*(1. + sqrt(2./(tau - 2.*nu)));
2759 *xi = (1. - nu/(2.*tau))*(1. + 1./sqrt(tau - nu));
2763 *eta = 1. + sqrt(
PI/(2.*tau));
2764 *xi = 1. + 0.75*sqrt(
PI/(2.*tau));
2768 double theta_nu =
ThetaNu(nu);
2770 double xxx = 1. + 1./sqrt(4.*tau+3.*nu);
2771 *eta =
POW2(xxx)*exp(-theta_nu/tau);
2773 *xi = (1. + nu/(2.*tau))*(1. + 1./sqrt(3./(2.*tau)+3.*nu))*exp(-theta_nu/tau);
2778 xxx = 0.25*pow(nu/tau,0.75)/(pow(nu/tau,0.75) + pow((25.+3.*nu)/5.,0.75)) +
2779 (1. + 0.75*sqrt(
PI/(2.*tau)))/(1. + sqrt(
PI/(2.*tau)));
2780 *xi = (
MIN2(xxx,1.) + theta_nu/(2.*tau))*(*eta);
2784 ASSERT( *eta >= 0. && *xi >= 0. );
2787 gv.bin[nd]->chrg[nz]->eta[ind] = *eta;
2788 gv.bin[nd]->chrg[nz]->xi[ind] = *xi;
2803 const double REL_TOLER = 10.*DBL_EPSILON;
2806 double xi_nu = 1. + 1./sqrt(3.*nu);
2807 double xi_nu2 =
POW2(xi_nu);
2813 double fnu = 2.*xi_nu2 - 1. - nu*xi_nu*
POW2(xi_nu2 - 1.);
2814 double dfdxi = 4.*xi_nu - nu*((5.*xi_nu2 - 6.)*xi_nu2 + 1.);
2816 xi_nu2 =
POW2(xi_nu);
2817 xxx = fabs(old-xi_nu);
2818 }
while( xxx > REL_TOLER*xi_nu );
2820 theta_nu = nu/xi_nu - 1./(2.*xi_nu2*(xi_nu2-1.));
2845 ASSERT( Zlo >=
gv.bin[nd]->LowestZg );
2859 fprintf(
ioQQQ,
" %ld/%ld", Zlo, stride );
2862 if(
gv.bin[nd]->nfill <
rfield.nflux )
2868 for( nz=0; nz <
gv.bin[nd]->nChrg; nz++ )
2874 Zg = Zlo + nz*stride;
2878 for( zz=0; zz <
NCHS-1; zz++ )
2880 if(
gv.bin[nd]->chrg[zz]->DustZ == Zg )
2889 ptr =
gv.bin[nd]->chrg[ind];
2892 for( zz=ind-1; zz >= nz; zz-- )
2893 gv.bin[nd]->chrg[zz+1] =
gv.bin[nd]->chrg[zz];
2895 gv.bin[nd]->chrg[nz] = ptr;
2897 if(
gv.bin[nd]->chrg[nz]->DustZ != Zg )
2899 else if(
gv.bin[nd]->chrg[nz]->nfill <
rfield.nflux )
2908 ASSERT(
gv.bin[nd]->chrg[nz]->DustZ == Zg );
2910 ASSERT( rate_up[nz] >= 0. && rate_dn[nz] >= 0. );
2920 for( nz=0; nz <
gv.bin[nd]->nChrg; nz++ )
2923 HighEnergy =
MAX2(HighEnergy,
2924 MAX2(
gv.bin[nd]->chrg[nz]->ThresInfInc,0.) + BoltzFac*
MAX2(
phycon.te,
gv.bin[nd]->tedust));
2927 gv.bin[nd]->qnflux2 =
ipoint(HighEnergy);
2952 ASSERT( Zlo >=
gv.bin[nd]->LowestZg );
2963 for( i=0; i < 2; i++ )
2968 sum = pop[i][0] = 1.;
2969 for( j=1; j <
gv.bin[nd]->nChrg-1; j++ )
2972 if( rate_dn[nz] > 10.*rate_up[nz-1]/sqrt(DBL_MAX) )
2974 pop[i][j] = pop[i][j-1]*rate_up[nz-1]/rate_dn[nz];
2979 for( k=0; k < j; k++ )
2987 if( pop[i][j] > sqrt(DBL_MAX) )
2989 for( k=0; k <= j; k++ )
2991 pop[i][k] /= DBL_MAX/10.;
2997 for( j=0; j <
gv.bin[nd]->nChrg-1; j++ )
3001 netloss[i] += pop[i][j]*(rate_up[nz] - rate_dn[nz]);
3007 if( netloss[0]*netloss[1] > 0. )
3008 *newZlo = ( netloss[1] > 0. ) ? Zlo + 1 : Zlo - 1;
3017 if(
gv.bin[nd]->nChrg > 2 &&
3018 ( *newZlo <
gv.bin[nd]->LowestZg ||
3019 ( *newZlo == Zlo && pop[1][
gv.bin[nd]->nChrg-2] < DBL_EPSILON ) ) )
3021 gv.bin[nd]->nChrg--;
3030 printf(
" fnzone %.2f nd %ld Zlo %ld newZlo %ld netloss %.4e %.4e nChrg %ld lgRedo %d\n",
3031 fnzone, nd, Zlo, *newZlo, netloss[0], netloss[1],
gv.bin[nd]->nChrg, lgRedo );
3036 if( *newZlo <
gv.bin[nd]->LowestZg )
3038 fprintf(
ioQQQ,
" could not converge charge state populations for %s\n",
gv.bin[nd]->chDstLab );
3043 if( *newZlo == Zlo )
3045 double frac0, frac1;
3047 double test1, test2, test3,
x1,
x2;
3051 frac0 = netloss[1]/(netloss[1]-netloss[0]);
3052 frac1 = -netloss[0]/(netloss[1]-netloss[0]);
3054 gv.bin[nd]->chrg[0]->FracPop = frac0*pop[0][0];
3055 gv.bin[nd]->chrg[
gv.bin[nd]->nChrg-1]->FracPop = frac1*pop[1][
gv.bin[nd]->nChrg-2];
3056 for( nz=1; nz <
gv.bin[nd]->nChrg-1; nz++ )
3058 gv.bin[nd]->chrg[nz]->FracPop = frac0*pop[0][nz] + frac1*pop[1][nz-1];
3062 test1 = test2 = test3 = 0.;
3063 for( nz=0; nz <
gv.bin[nd]->nChrg; nz++ )
3065 ASSERT(
gv.bin[nd]->chrg[nz]->FracPop >= 0. );
3066 test1 +=
gv.bin[nd]->chrg[nz]->FracPop;
3067 test2 +=
gv.bin[nd]->chrg[nz]->FracPop*rate_up[nz];
3068 test3 +=
gv.bin[nd]->chrg[nz]->FracPop*(rate_up[nz]-rate_dn[nz]);
3070 x1 = fabs(test1-1.);
3105 double *anu =
rfield.anu;
3110 const double INV_DELTA_E =
EVRYD/3.;
3111 const double CS_PDT = 1.2e-17;
3123 gv.bin[nd]->chrg[nz]->DustZ = Zg;
3126 memset(
gv.bin[nd]->chrg[nz]->eta, 0, (
LIMELM+2)*
sizeof(
double) );
3127 memset(
gv.bin[nd]->chrg[nz]->xi, 0, (
LIMELM+2)*
sizeof(
double) );
3129 GetPotValues(nd,Zg,&
gv.bin[nd]->chrg[nz]->ThresInf,&
gv.bin[nd]->chrg[nz]->ThresInfVal,
3130 &
gv.bin[nd]->chrg[nz]->ThresSurf,&
gv.bin[nd]->chrg[nz]->ThresSurfVal,
3131 &
gv.bin[nd]->chrg[nz]->PotSurf,&
gv.bin[nd]->chrg[nz]->Emin,
INCL_TUNNEL);
3135 GetPotValues(nd,Zg-1,&
gv.bin[nd]->chrg[nz]->ThresInfInc,&d[0],&
gv.bin[nd]->chrg[nz]->ThresSurfInc,
3136 &d[1],&
gv.bin[nd]->chrg[nz]->PotSurfInc,&
gv.bin[nd]->chrg[nz]->EminInc,
NO_TUNNEL);
3138 gv.bin[nd]->chrg[nz]->ipThresInfVal =
3142 ipLo =
gv.bin[nd]->chrg[nz]->ipThresInfVal;
3145 gv.bin[nd]->chrg[nz]->nfill =
rfield.nflux;
3146 nfill =
gv.bin[nd]->chrg[nz]->nfill;
3149 long len = nfill + 10 - ipLo;
3152 gv.bin[nd]->chrg[nz]->yhat.reserve( len );
3153 gv.bin[nd]->chrg[nz]->yhat_primary.reserve( len );
3154 gv.bin[nd]->chrg[nz]->ehat.reserve( len );
3155 gv.bin[nd]->chrg[nz]->yhat.alloc( ipLo, nfill );
3156 gv.bin[nd]->chrg[nz]->yhat_primary.alloc( ipLo, nfill );
3157 gv.bin[nd]->chrg[nz]->ehat.alloc( ipLo, nfill );
3161 gv.bin[nd]->chrg[nz]->yhat.realloc( nfill );
3162 gv.bin[nd]->chrg[nz]->yhat_primary.realloc( nfill );
3163 gv.bin[nd]->chrg[nz]->ehat.realloc( nfill );
3170 DustWorkFcn =
gv.bin[nd]->DustWorkFcn;
3171 Elo = -
gv.bin[nd]->chrg[nz]->PotSurf;
3172 ThresInfVal =
gv.bin[nd]->chrg[nz]->ThresInfVal;
3173 Emin =
gv.bin[nd]->chrg[nz]->Emin;
3178 ASSERT(
gv.bin[nd]->sd[0]->ipLo <= ipLo );
3180 for( i=
max(ipLo,ipStart); i < nfill; i++ )
3184 double Yp1,Ys1,Ehp,Ehs,yp,ya,ys,eyp,eya,eys;
3185 double yzero,yone,Ehi,Ehi_band,Wcorr,Eel;
3189 eyp = eya = eys = 0.;
3192 Ehi = Ehi_band = anu[i] - ThresInfVal - Emin;
3193 Wcorr = ThresInfVal + Emin - GrainPot;
3194 Eel = anu[i] - DustWorkFcn;
3195 yzero =
y0b( nd, nz, i );
3196 yone = sptr->
y01[i];
3197 Yfunc(nd,nz,yzero*yone,sptr->
p[i],Elo,Ehi,Eel,&Yp1,&Ys1,&Ehp,&Ehs);
3204 unsigned int nsmax =
gv.bin[nd]->sd.size();
3205 for( ns=1; ns < nsmax; ns++ )
3207 sptr =
gv.bin[nd]->sd[ns];
3209 if( i < sptr->ipLo )
3212 Ehi = Ehi_band + Wcorr - sptr->
ionPot;
3213 Eel = anu[i] - sptr->
ionPot;
3214 Yfunc(nd,nz,sptr->
y01[i],sptr->
p[i],Elo,Ehi,Eel,&Yp1,&Ys1,&Ehp,&Ehs);
3221 long nmax = sptr->
nData;
3222 for( n=0; n < nmax; n++ )
3224 double max = sptr->
AvNr[n]*sptr->
p[i];
3225 Ehi = sptr->
Ener[n] - GrainPot;
3226 Eel = sptr->
Ener[n];
3227 Yfunc(nd,nz,sptr->
y01A[n][i],
max,Elo,Ehi,Eel,&Yp1,&Ys1,&Ehp,&Ehs);
3235 gv.bin[nd]->chrg[nz]->yhat[i] = (
realnum)(yp + ya + ys);
3236 gv.bin[nd]->chrg[nz]->yhat_primary[i] =
min((
realnum)yp,1.f);
3237 gv.bin[nd]->chrg[nz]->ehat[i] = (
gv.bin[nd]->chrg[nz]->yhat[i] > 0. ) ?
3238 (
realnum)((eyp + eya + eys)/
gv.bin[nd]->chrg[nz]->yhat[i]) : 0.f;
3241 ASSERT(
gv.bin[nd]->chrg[nz]->ehat[i] >= 0.f );
3250 gv.bin[nd]->chrg[nz]->ipThresInf =
3254 ipLo =
gv.bin[nd]->chrg[nz]->ipThresInf;
3256 len = nfill + 10 - ipLo;
3263 gv.bin[nd]->chrg[nz]->cs_pdt.reserve( len );
3264 gv.bin[nd]->chrg[nz]->cs_pdt.alloc( ipLo, nfill );
3268 gv.bin[nd]->chrg[nz]->cs_pdt.realloc( nfill );
3273 c1 = -CS_PDT*(double)Zg;
3274 ThresInf =
gv.bin[nd]->chrg[nz]->ThresInf;
3275 cnv_GR_pH =
gv.bin[nd]->cnv_GR_pH;
3277 for( i=
max(ipLo,ipStart); i < nfill; i++ )
3279 double x = (anu[i] - ThresInf)*INV_DELTA_E;
3280 double cs = c1*x/
POW2(1.+(1./3.)*
POW2(x));
3283 gv.bin[nd]->chrg[nz]->cs_pdt[i] =
MAX2(cs,0.)*cnv_GR_pH;
3291 gv.bin[nd]->chrg[nz]->fac1.reserve( len );
3292 gv.bin[nd]->chrg[nz]->fac2.reserve( len );
3293 gv.bin[nd]->chrg[nz]->fac1.alloc( ipLo, nfill );
3294 gv.bin[nd]->chrg[nz]->fac2.alloc( ipLo, nfill );
3298 gv.bin[nd]->chrg[nz]->fac1.realloc( nfill );
3299 gv.bin[nd]->chrg[nz]->fac2.realloc( nfill );
3304 for( i=
max(ipLo,ipStart); i < nfill; i++ )
3306 double cs1,cs2,cs_tot,cool1,cool2,ehat1,ehat2;
3309 PE_init(nd,nz,i,&cs1,&cs2,&cs_tot,&cool1,&cool2,&ehat1,&ehat2);
3311 gv.bin[nd]->chrg[nz]->fac1[i] = cs_tot*anu[i]-cs1*cool1-cs2*cool2;
3312 gv.bin[nd]->chrg[nz]->fac2[i] = cs1*ehat1 + cs2*ehat2;
3314 ASSERT(
gv.bin[nd]->chrg[nz]->fac1[i] >= 0. &&
gv.bin[nd]->chrg[nz]->fac2[i] >= 0. );
3326 gv.bin[nd]->chrg[nz]->FracPop = -DBL_MAX;
3328 gv.bin[nd]->chrg[nz]->RSum1 = -DBL_MAX;
3329 gv.bin[nd]->chrg[nz]->RSum2 = -DBL_MAX;
3330 gv.bin[nd]->chrg[nz]->ESum1a = -DBL_MAX;
3331 gv.bin[nd]->chrg[nz]->ESum1b = -DBL_MAX;
3332 gv.bin[nd]->chrg[nz]->ESum2 = -DBL_MAX;
3334 gv.bin[nd]->chrg[nz]->tedust = 1.f;
3336 gv.bin[nd]->chrg[nz]->hcon1 = -DBL_MAX;
3337 gv.bin[nd]->chrg[nz]->hots1 = -DBL_MAX;
3338 gv.bin[nd]->chrg[nz]->bolflux1 = -DBL_MAX;
3339 gv.bin[nd]->chrg[nz]->pe1 = -DBL_MAX;
3341 gv.bin[nd]->chrg[nz]->BolFlux = -DBL_MAX;
3342 gv.bin[nd]->chrg[nz]->GrainHeat = -DBL_MAX;
3343 gv.bin[nd]->chrg[nz]->GrainHeatColl = -DBL_MAX;
3344 gv.bin[nd]->chrg[nz]->GasHeatPhotoEl = -DBL_MAX;
3345 gv.bin[nd]->chrg[nz]->GasHeatTherm = -DBL_MAX;
3346 gv.bin[nd]->chrg[nz]->GrainCoolTherm = -DBL_MAX;
3347 gv.bin[nd]->chrg[nz]->ChemEnIon = -DBL_MAX;
3348 gv.bin[nd]->chrg[nz]->ChemEnH2 = -DBL_MAX;
3350 gv.bin[nd]->chrg[nz]->HeatingRate2 = -DBL_MAX;
3353 ASSERT(
gv.bin[nd]->chrg[nz]->ipThresInf <=
gv.bin[nd]->chrg[nz]->ipThresInfVal );
3375 ThermExp =
gv.bin[nd]->chrg[nz]->ThresInf*
TE1RYD/
gv.bin[nd]->tedust;
3377 gv.bin[nd]->chrg[nz]->ThermRate =
THERMCONST*
gv.bin[nd]->ThermEff*
POW2(
gv.bin[nd]->tedust)*exp(-ThermExp);
3378# if defined( WD_TEST2 ) || defined( IGNORE_THERMIONIC )
3379 gv.bin[nd]->chrg[nz]->ThermRate = 0.;
3400 long Zg =
gv.bin[nd]->chrg[nz]->DustZ;
3405 y2pr =
y2pa( Elo, Ehi, Zg, Ehp );
3412 *Yp = y2pr*
min(y01,maxval);
3414 y2sec =
y2s( Elo, Ehi, Zg, Ehs );
3417 else if( pcase ==
PE_SIL )
3421 fprintf(
ioQQQ,
" Yfunc: unknown type for PE effect: %d\n" , pcase );
3427 *Ys = y2sec*f3*
min(y01,maxval);
3450 yzero =
y0b01( nd, nz, i );
3453 double Eph =
rfield.anu[i];
3455 if( Eph <= 20./
EVRYD )
3456 yzero =
y0b01( nd, nz, i );
3457 else if( Eph < 50./
EVRYD )
3459 double y0a =
y0b01( nd, nz, i );
3460 double y0b =
gv.bin[nd]->y0b06[i];
3462 double frac = log(Eph*(
EVRYD/20.))*1.0913566679372915;
3464 yzero = y0a * exp(log(
y0b/y0a)*frac);
3467 yzero =
gv.bin[nd]->y0b06[i];
3485 xv =
MAX2((
rfield.anu[i] -
gv.bin[nd]->chrg[nz]->ThresSurfVal)/
gv.bin[nd]->DustWorkFcn,0.);
3492 yzero = xv/((1./9.e-3) + (3.7e-2/9.e-3)*xv);
3496 yzero = xv/(2.+10.*xv);
3499 fprintf(
ioQQQ,
" y0b01: unknown type for PE effect: %d\n" , pcase );
3514 double yzero, leola;
3518 ASSERT( i >=
gv.bin[nd]->sd[ns]->ipLo );
3527 yzero =
gv.bin[nd]->sd[ns]->p[i]*leola*(1. - leola*log(1.+1./leola));
3530 double x = 1./leola;
3531 yzero =
gv.bin[nd]->sd[ns]->p[i]*(((-1./5.*x+1./4.)*x-1./3.)*x+1./2.);
3544 double alpha, beta, af, bf, yone;
3548 beta =
gv.bin[nd]->AvRadius*
gv.bin[nd]->inv_att_len[i];
3550 bf =
pow2(beta) - 2.*beta + 2. - 2.*exp(-beta);
3552 bf = ((1./60.*beta - 1./12.)*beta + 1./3.)*
pow3(beta);
3556 af =
pow2(alpha) - 2.*alpha + 2. - 2.*exp(-alpha);
3558 af = ((1./60.*alpha - 1./12.)*alpha + 1./3.)*
pow3(alpha);
3560 yone =
pow2(beta/alpha)*af/bf;
3582 *Ehp = 0.5*Ehi*(1.-2.*x)/(1.-3.*x);
3584 ytwo = ( abs(x) > 1e-4 ) ? (1.-3.*x)/
pow3(1.-x) : 1. - (3. + 8.*x)*x*x;
3585 ASSERT( *Ehp > 0. && *Ehp <= Ehi && ytwo > 0. && ytwo <= 1. );
3597 *Ehp = 0.5*(Elo+Ehi);
3599 ASSERT( *Ehp >= Elo && *Ehp <= Ehi );
3623 if( !
gv.lgWD01 && Ehi > 0. )
3632 double x2 = x*x, x3 =
x2*x, x4 = x3*x, x5 = x4*x;
3633 double yh2 = yh*yh, yh3 = yh2*yh, yh4 = yh3*yh, yh5 = yh4*yh;
3634 double help1 = 2.*x-yh;
3635 double help2 = (6.*x3-15.*yh*
x2+12.*yh2*x-3.*yh3)/4.;
3636 double help3 = (22.*x5-95.*yh*x4+164.*yh2*x3-141.*yh3*
x2+60.*yh4*x-10.*yh5)/16.;
3637 N0 = yh*(help1 - help2 + help3)/
x2;
3639 help1 = (3.*x-yh)/3.;
3640 help2 = (15.*x3-25.*yh*
x2+15.*yh2*x-3.*yh3)/20.;
3641 help3 = (1155.*x5-3325.*yh*x4+4305.*yh2*x3-2961.*yh3*
x2+1050.*yh4*x-150.*yh5)/1680.;
3642 E0 =
ETILDE*yh2*(help1 - help2 + help3)/
x2;
3646 double sR0 = (1. + yl*yl);
3647 double sqR0 = sqrt(sR0);
3648 double sqRh = sqrt(1. + x*x);
3649 double alpha = sqRh/(sqRh - 1.);
3650 if( yh/sqR0 < 0.01 )
3653 double z = yh*(yh - 2.*yl)/sR0;
3654 N0 = ((((7./256.*z-5./128.)*z+1./16.)*z-1./8.)*z+1./2.)*z/(sqRh-1.);
3656 double yl2 = yl*yl, yl3 = yl2*yl, yl4 = yl3*yl;
3657 double help1 = yl/2.;
3658 double help2 = (2.*yl2-1.)/3.;
3659 double help3 = (6.*yl3-9.*yl)/8.;
3660 double help4 = (8.*yl4-24.*yl2+3.)/10.;
3662 E0 = -alpha*Ehi*(((help4*h + help3)*h + help2)*h + help1)*h/sqR0;
3666 N0 = alpha*(1./sqR0 - 1./sqRh);
3667 E0 = alpha*
ETILDE*(
ASINH(x*sqR0 + yl*sqRh) - yh/sqRh);
3670 ASSERT( N0 > 0. && N0 <= 1. );
3674 ASSERT( *Ehs > 0. && *Ehs <= Ehi );
3686 if( !
gv.lgWD01 && Ehi > Elo )
3694 double sqRh = sqrt(1. +
x2);
3695 double alpha = sqRh/(sqRh - 1.);
3701 *Ehs = Ehi - (Ehi-Elo)*((-37./840.*
x2 + 1./10.)*
x2 + 1./3.);
3704 ASSERT( *Ehs >= Elo && *Ehs <= Ehi );
3728 for( nelem=
LIMELM-1; nelem >= 0; nelem-- )
3730 if(
dense.lgElmtOn[nelem] )
3732 for( ion=nelem+1; ion >= 0; ion-- )
3734 if( ion == high ||
dense.xIonDense[nelem][ion] > 0. )
3737 high =
MAX2(high,ion);
3748 bool lgAllIonStages)
3761 Zg =
gv.bin[nd]->chrg[nz]->DustZ;
3763 hi_ion = ( lgAllIonStages ) ?
LIMELM :
gv.HighestIon;
3765 phi_s_up[0] =
gv.bin[nd]->chrg[nz]->ThresSurf;
3766 for( i=1; i <=
LIMELM; i++ )
3771 phi_s_up[i] = -DBL_MAX;
3773 phi_s_dn[0] =
gv.bin[nd]->chrg[nz]->ThresSurfInc;
3777 for( nelem=0; nelem <
LIMELM; nelem++ )
3779 if(
dense.lgElmtOn[nelem] )
3781 for( ion=0; ion <= nelem+1; ion++ )
3783 if( lgAllIonStages ||
dense.xIonDense[nelem][ion] > 0. )
3786 &
gv.bin[nd]->chrg[nz]->RecomZ0[nelem][ion],
3787 &
gv.bin[nd]->chrg[nz]->RecomEn[nelem][ion],
3788 &
gv.bin[nd]->chrg[nz]->ChemEn[nelem][ion]);
3792 gv.bin[nd]->chrg[nz]->RecomZ0[nelem][ion] = ion;
3793 gv.bin[nd]->chrg[nz]->RecomEn[nelem][ion] = 0.f;
3794 gv.bin[nd]->chrg[nz]->ChemEn[nelem][ion] = 0.f;
3805 double *ThresInfVal,
3807 double *ThresSurfVal,
3810 bool lgUseTunnelCorr)
3839 IP =
gv.bin[nd]->DustWorkFcn -
gv.bin[nd]->BandGap + dstpot - 0.5*
one_elec(nd);
3849 fprintf(
ioQQQ,
" GetPotValues detected unknown type for ionization pot: %d\n" , pcase );
3855 IP_v =
MAX2(IP,IP_v);
3860 double help = fabs(dZg+1);
3863 if( lgUseTunnelCorr )
3866 *Emin *= 1. - 2.124e-4/(pow(
gv.bin[nd]->AvRadius,(
realnum)0.45)*pow(help,0.26));
3874 *ThresInf = IP - *Emin;
3875 *ThresInfVal = IP_v - *Emin;
3876 *ThresSurf = *ThresInf;
3877 *ThresSurfVal = *ThresInfVal;
3883 *ThresInfVal = IP_v;
3884 *ThresSurf = *ThresInf - dstpot;
3885 *ThresSurfVal = *ThresInfVal - dstpot;
3900 const double phi_s_up[],
3901 const double phi_s_dn[],
3918 Zg =
gv.bin[nd]->chrg[nz]->DustZ;
3919 phi_s = phi_s_up[0];
3923 *ChemEn +=
rfield.anu[
Heavy.ipHeavy[nelem][ion-1]-1];
3927 *ChemEn -= (
realnum)(phi_s - phi_s_up[0]);
3930 phi_s = phi_s_up[
save-ion];
3935 else if( ion <= nelem &&
gv.bin[nd]->chrg[nz]->DustZ >
gv.bin[nd]->LowestZg &&
3941 Zg =
gv.bin[nd]->chrg[nz]->DustZ;
3942 phi_s = phi_s_dn[0];
3946 *ChemEn -=
rfield.anu[
Heavy.ipHeavy[nelem][ion]-1];
3950 *ChemEn += (
realnum)(phi_s - phi_s_dn[0]);
3955 phi_s = phi_s_dn[ion-
save];
3959 }
while( ion <= nelem && Zg >
gv.bin[nd]->LowestZg &&
3979 double fac0 =
STICK_ION*
gv.bin[nd]->IntArea/4.*
gv.bin[nd]->cnv_H_pCM3;
3983# ifndef IGNORE_GRAIN_ION_COLLISIONS
3985 for( nz=0; nz <
gv.bin[nd]->nChrg; nz++ )
3989 double fac1 = gptr->
FracPop*fac0;
3994 for( ion=0; ion <=
LIMELM; ion++ )
3997 double eta, fac2, xi;
4008 for( nelem=
MAX2(0,ion-1); nelem <
LIMELM; nelem++ )
4010 if(
dense.lgElmtOn[nelem] && ion != gptr->
RecomZ0[nelem][ion] )
4012 gv.GrainChTrRate[nelem][ion][gptr->
RecomZ0[nelem][ion]] +=
4031 for( nelem=0; nelem <
LIMELM; nelem++ )
4033 gv.elmSumAbund[nelem] = 0.f;
4037 for(
size_t nd=0; nd <
gv.bin.size(); nd++ )
4040 gv.bin[nd]->dstAbund = (
realnum)(
gv.bin[nd]->dstfactor*
gv.GrainMetal*
gv.bin[nd]->GrnDpth);
4041 ASSERT(
gv.bin[nd]->dstAbund > 0.f );
4045 gv.bin[nd]->cnv_CM3_pH = 1./
gv.bin[nd]->cnv_H_pCM3;
4047 gv.bin[nd]->cnv_CM3_pGR =
gv.bin[nd]->cnv_H_pGR/
gv.bin[nd]->cnv_H_pCM3;
4048 gv.bin[nd]->cnv_GR_pCM3 = 1./
gv.bin[nd]->cnv_CM3_pGR;
4053 for( nelem=0; nelem <
LIMELM; nelem++ )
4055 gv.elmSumAbund[nelem] +=
gv.bin[nd]->elmAbund[nelem]*(
realnum)
gv.bin[nd]->cnv_H_pCM3;
4068 for(
long i=0; i <
rfield.nupper; i++ )
4076 for(
size_t nd=0; nd <
gv.bin.size(); nd++ )
4078 realnum dstAbund =
gv.bin[nd]->dstAbund;
4081 for(
long i=0; i <
rfield.nflux; i++ )
4089 gv.dstab[i] +=
gv.bin[nd]->dstab1[i]*dstAbund;
4090 gv.dstsc[i] +=
gv.bin[nd]->pure_sc1[i]*
gv.bin[nd]->asym[i]*dstAbund;
4093 for(
long nz=0; nz <
gv.bin[nd]->nChrg; nz++ )
4096 if( gptr->
DustZ <= -1 )
4098 double FracPop = gptr->
FracPop;
4101 gv.dstab[i] += FracPop*gptr->
cs_pdt[i]*dstAbund;
4106 for(
long i=0; i <
rfield.nflux; i++ )
4109 ASSERT(
gv.dstab[i] > 0. &&
gv.dstsc[i] > 0. );
4125 double EhatThermionic,
4139 fprintf(
ioQQQ,
" GrainTemperature starts for grain %s\n",
gv.bin[nd]->chDstLab );
4156 gv.bin[nd]->GasHeatPhotoEl = 0.;
4158 gv.bin[nd]->GrainCoolTherm = 0.;
4159 gv.bin[nd]->thermionic = 0.;
4164 gv.bin[nd]->BolFlux = 0.;
4168 for( nz=0; nz <
gv.bin[nd]->nChrg; nz++ )
4175 double bolflux1 = 0.;
4179 bool lgReEvaluate1 = gptr->
hcon1 < 0.;
4180 bool lgReEvaluate2 = gptr->
hots1 < 0.;
4181 bool lgReEvaluate = lgReEvaluate1 || lgReEvaluate2;
4187 for( i=0; i < loopmax; i++ )
4189 double fac =
gv.bin[nd]->dstab1[i]*
rfield.anu[i];
4192 hcon1 +=
rfield.flux[0][i]*fac;
4196 hots1 +=
rfield.SummedDif[i]*fac;
4198 bolflux1 +=
rfield.SummedCon[i]*fac;
4217 gptr->
hcon1 = hcon1;
4221 hcon1 = gptr->
hcon1;
4238 bolflux1 +=
rfield.SummedCon[i]*
gv.bin[nd]->dstab1[i]*
rfield.anu[i];
4239 if( gptr->
DustZ <= -1 )
4243 gptr->
hots1 = hots1;
4249 hots1 = gptr->
hots1;
4260 hla1 =
rfield.otslin[ipLya]*
gv.bin[nd]->dstab1[ipLya]*0.75;
4262 else if( ipLya <
rfield.nflux )
4265 hla1 =
rfield.otslin[ipLya]*gptr->
fac1[ipLya];
4272 ASSERT( hcon1 >= 0. && hots1 >= 0. && hla1 >= 0. && bolflux1 >= 0. && pe1 >= 0. );
4277 gv.bin[nd]->BolFlux += gptr->
FracPop*bolflux1;
4279 gv.bin[nd]->GasHeatPhotoEl += gptr->
FracPop*pe1;
4284 fprintf(
ioQQQ,
" Zg %ld bolflux: %.4e\n", gptr->
DustZ,
4299 gv.bin[nd]->thermionic += rate * (EhatThermionic - gptr->
PotSurf*
EN1RYD);
4303 norm =
EN1RYD*
gv.bin[nd]->cnv_H_pCM3;
4314 gv.bin[nd]->BolFlux *= norm;
4325 gv.bin[nd]->GasHeatPhotoEl *= norm;
4327 if(
gv.lgBakesPAH_heat )
4334 gv.bin[nd]->GasHeatPhotoEl = 1.e-24*
hmi.UV_Cont_rel2_Habing_TH85_depth*
4335 dense.gas_phase[
ipHYDROGEN]*(4.87e-2/(1.0+4e-3*pow((
hmi.UV_Cont_rel2_Habing_TH85_depth*
4338 (1.+2.e-4*(
hmi.UV_Cont_rel2_Habing_TH85_depth*
phycon.sqrte/
dense.eden)))/
gv.bin.size();
4344 gv.bin[nd]->GasHeatPhotoEl *=
gv.GrainHeatScaleFactor;
4361 gv.bin[nd]->GrainHeat = *hcon + *hots + dcheat -
gv.bin[nd]->GrainCoolTherm;
4364 gv.bin[nd]->GrainHeatColl = dcheat;
4370 if(
gv.bin[nd]->GrainHeat > 0. )
4375 x = log(
MAX2(DBL_MIN,
gv.bin[nd]->GrainHeat*
gv.bin[nd]->cnv_CM3_pH));
4382 gv.bin[nd]->GrainHeat = -1.;
4383 gv.bin[nd]->tedust = -1.;
4386 if(
thermal.ConstGrainTemp > 0. )
4390 gv.bin[nd]->tedust =
thermal.ConstGrainTemp;
4392 x = log(
gv.bin[nd]->tedust);
4394 gv.bin[nd]->GrainHeat = exp(y)*
gv.bin[nd]->cnv_H_pCM3;
4398 gv.bin[nd]->TeGrainMax = (
realnum)
MAX2(
gv.bin[nd]->TeGrainMax,
gv.bin[nd]->tedust);
4402 fprintf(
ioQQQ,
" >GrainTemperature finds %s Tdst %.5e hcon %.4e ",
4403 gv.bin[nd]->chDstLab,
gv.bin[nd]->tedust, *hcon);
4404 fprintf(
ioQQQ,
"hots %.4e dcheat %.4e GrainCoolTherm %.4e\n",
4405 *hots, dcheat,
gv.bin[nd]->GrainCoolTherm );
4431 ASSERT( nz >= 0 && nz <
gv.bin[nd]->nChrg );
4440 *cs1 =
gv.bin[nd]->dstab1[i]*gptr->
yhat[i];
4443 *ehat1 = gptr->
ehat[i];
4464 if( gptr->
DustZ <= -1 )
4469 ASSERT( *ehat1 > 0. && *cool1 > 0. );
4479 if( gptr->
DustZ <= -1 && i >= ipLo2 )
4488 ASSERT( *ehat2 > 0. && *cool2 > 0. );
4497 *cs_tot =
gv.bin[nd]->dstab1[i] + *cs2;
4511 double Accommodation,
4512 CollisionRateElectr,
4539 const double H2_FORMATION_GRAIN_HEATING[
H2_TOP] = { 0.20, 0.4, 1.72 };
4559 gv.bin[nd]->ChemEn = 0.;
4562 for( nz=0; nz <
gv.bin[nd]->nChrg; nz++ )
4570 double ChemEn1 = 0.;
4576 for( ion=0; ion <=
LIMELM; ion++ )
4585 CollisionRateIon = 0.;
4587 CoolPotentialGas = 0.;
4588 HeatRecombination = 0.;
4595 for( nelem=
MAX2(0,ion-1); nelem <
LIMELM; nelem++ )
4597 if(
dense.lgElmtOn[nelem] &&
dense.xIonDense[nelem][ion] > 0. )
4599 double CollisionRateOne;
4604#if defined( IGNORE_GRAIN_ION_COLLISIONS )
4606#elif defined( WD_TEST2 )
4607 Stick = ( ion == gptr->
RecomZ0[nelem][ion] ) ?
4610 Stick = ( ion == gptr->
RecomZ0[nelem][ion] ) ?
4617 CollisionRateIon += CollisionRateOne;
4626 if( ion >= gptr->
RecomZ0[nelem][ion] )
4628 CoolPotential += CollisionRateOne * (double)ion *
4630 CoolPotentialGas += CollisionRateOne *
4631 (double)gptr->
RecomZ0[nelem][ion] *
4636 CoolPotential += CollisionRateOne * (double)ion *
4638 CoolPotentialGas += CollisionRateOne *
4639 (double)gptr->
RecomZ0[nelem][ion] *
4646 HeatRecombination += CollisionRateOne *
4648 HeatChem += CollisionRateOne * gptr->
ChemEn[nelem][ion];
4662 CoolPotential *= eta*
EN1RYD;
4663 CoolPotentialGas *= eta*
EN1RYD;
4665 HeatRecombination *= eta*
EN1RYD;
4668 CoolEmitted = CollisionRateIon * 2.*
BOLTZMANN*
gv.bin[nd]->tedust*eta;
4671 Heat1 += HeatCollisions - CoolPotential + HeatRecombination - CoolEmitted;
4676 Cool1 += HeatCollisions - CoolEmitted - CoolPotentialGas;
4678 ChemEn1 += HeatChem;
4682 HeatIons += gptr->
FracPop*Heat1;
4686 fprintf(
ioQQQ,
" Zg %ld ions heat/cool: %.4e %.4e\n", gptr->
DustZ,
4687 gptr->
FracPop*Heat1*
gv.bin[nd]->IntArea/4.*
gv.bin[nd]->cnv_H_pCM3,
4688 gptr->
FracPop*Cool1*
gv.bin[nd]->IntArea/4.*
gv.bin[nd]->cnv_H_pCM3 );
4695 Stick = ( gptr->
DustZ <= -1 ) ?
gv.bin[nd]->StickElecNeg :
gv.bin[nd]->StickElecPos;
4701 CollisionRateElectr = Stick*
dense.eden*ve;
4707 if( gptr->
DustZ >
gv.bin[nd]->LowestZg )
4712 CoolPotential = CollisionRateElectr * (double)ion*gptr->
PotSurfInc*eta*
EN1RYD;
4720 HeatCollisions = 0.;
4722 HeatRecombination = 0.;
4730 CoolBounce = CollisionRateElectr *
4732 CoolBounce =
MAX2(CoolBounce,0.);
4737 HeatElectrons = HeatCollisions-CoolPotential+HeatRecombination+HeatBounce-CoolBounce;
4738 Heat1 += HeatElectrons;
4740 CoolElectrons = HeatCollisions+HeatBounce-CoolBounce;
4741 Cool1 += CoolElectrons;
4745 fprintf(
ioQQQ,
" Zg %ld electrons heat/cool: %.4e %.4e\n", gptr->
DustZ,
4746 gptr->
FracPop*HeatElectrons*
gv.bin[nd]->IntArea/4.*
gv.bin[nd]->cnv_H_pCM3,
4747 gptr->
FracPop*CoolElectrons*
gv.bin[nd]->IntArea/4.*
gv.bin[nd]->cnv_H_pCM3 );
4762 gv.bin[nd]->GrainCoolTherm*
gv.bin[nd]->cnv_CM3_pH;
4768 HeatTot += gptr->
FracPop*Heat1;
4771 CoolTot += gptr->
FracPop*Cool1;
4773 gv.bin[nd]->ChemEn += gptr->
FracPop*ChemEn1;
4784 Accommodation = 2.*
gv.bin[nd]->atomWeight*WeightMol/
POW2(
gv.bin[nd]->atomWeight+WeightMol);
4786#ifndef IGNORE_GRAIN_ION_COLLISIONS
4788 CollisionRateMol = Accommodation*
hmi.H2_total*
4792 ipH2 =
gv.which_H2distr[
gv.bin[nd]->matType];
4795 gv.bin[nd]->ChemEnH2 =
gv.bin[nd]->rate_h2_form_grains_used*
dense.xIonDense[
ipHYDROGEN][0]*
4798 gv.bin[nd]->ChemEnH2 /=
gv.bin[nd]->IntArea/4.*
gv.bin[nd]->cnv_H_pCM3;
4800 CollisionRateMol = 0.;
4801 gv.bin[nd]->ChemEnH2 = 0.;
4806 Accommodation = 2.*
gv.bin[nd]->atomWeight*WeightMol/
POW2(
gv.bin[nd]->atomWeight+WeightMol);
4807#ifndef IGNORE_GRAIN_ION_COLLISIONS
4811 CollisionRateMol = 0.;
4816 CoolEmitted = CollisionRateMol * 2.*
BOLTZMANN*
gv.bin[nd]->tedust;
4818 HeatMolecules = HeatCollisions - CoolEmitted +
gv.bin[nd]->ChemEnH2;
4819 HeatTot += HeatMolecules;
4822 CoolMolecules = HeatCollisions - CoolEmitted;
4823 CoolTot += CoolMolecules;
4825 gv.bin[nd]->RateUp = 0.;
4826 gv.bin[nd]->RateDn = 0.;
4828 for( nz=0; nz <
gv.bin[nd]->nChrg; nz++ )
4834 gv.bin[nd]->RateUp +=
gv.bin[nd]->chrg[nz]->FracPop*rate_up;
4835 gv.bin[nd]->RateDn +=
gv.bin[nd]->chrg[nz]->FracPop*rate_dn;
4838 HeatCor += (
gv.bin[nd]->chrg[nz]->FracPop*rate_up*
gv.bin[nd]->chrg[nz]->ThresSurf -
4839 gv.bin[nd]->chrg[nz]->FracPop*rate_dn*
gv.bin[nd]->chrg[nz]->ThresSurfInc +
4840 gv.bin[nd]->chrg[nz]->FracPop*rate_up*
gv.bin[nd]->chrg[nz]->PotSurf -
4841 gv.bin[nd]->chrg[nz]->FracPop*rate_dn*
gv.bin[nd]->chrg[nz]->PotSurfInc)*
EN1RYD;
4849 fprintf(
ioQQQ,
" molecules heat/cool: %.4e %.4e heatcor: %.4e\n",
4850 HeatMolecules*
gv.bin[nd]->IntArea/4.*
gv.bin[nd]->cnv_H_pCM3,
4851 CoolMolecules*
gv.bin[nd]->IntArea/4.*
gv.bin[nd]->cnv_H_pCM3,
4852 HeatCor*
gv.bin[nd]->IntArea/4.*
gv.bin[nd]->cnv_H_pCM3 );
4855 *dcheat = (
realnum)(HeatTot*
gv.bin[nd]->IntArea/4.*
gv.bin[nd]->cnv_H_pCM3);
4856 *dccool = (
gv.lgDColOn ) ? (
realnum)(CoolTot*
gv.bin[nd]->IntArea/4.*
gv.bin[nd]->cnv_H_pCM3) : 0.f;
4858 gv.bin[nd]->ChemEn *=
gv.bin[nd]->IntArea/4.*
gv.bin[nd]->cnv_H_pCM3;
4859 gv.bin[nd]->ChemEnH2 *=
gv.bin[nd]->IntArea/4.*
gv.bin[nd]->cnv_H_pCM3;
4876 gv.bin[nd]->HeatingRate1 = (HeatMolecules+HeatIons+HeatCor)*
gv.bin[nd]->IntArea/4.;
4904 vector<realnum> help(
rfield.nflux );
4905 for( i=0; i <
rfield.nflux; i++ )
4911 for(
size_t nd=0; nd <
gv.bin.size(); nd++ )
4915 for( i=0; i <
rfield.nflux; i++ )
4918 dmomen += help[i]*(
gv.bin[nd]->dstab1[i] +
gv.bin[nd]->pure_sc1[i]*
gv.bin[nd]->asym[i]);
4921 dmomen *=
EN1RYD*4./
gv.bin[nd]->IntArea;
4932 alam = log(20.702/rdust/psi*
phycon.sqrte/
dense.SqrtEden);
4939 phi2lm =
POW2(psi)*alam;
4942 for( loop = 0; loop < 50 && fabs(corr-1.) > 0.001; loop++ )
4944 vdold =
gv.bin[nd]->DustDftVel;
4947 si =
gv.bin[nd]->DustDftVel/
phycon.sqrte*7.755e-5;
4948 g0 = 1.5045*si*sqrt(1.+0.4418*si*si);
4949 g2 = si/(1.329 +
POW3(si));
4956 si =
gv.bin[nd]->DustDftVel/
phycon.sqrte*1.816e-6;
4957 g0 = 1.5045*si*sqrt(1.+0.4418*si*si);
4958 g2 = si/(1.329 +
POW3(si));
4959 fdrag += fac*
dense.eden*(g0 + phi2lm*g2);
4962 si =
gv.bin[nd]->DustDftVel/
phycon.sqrte*7.755e-5;
4963 g0 = 1.5045*si*sqrt(1.+0.4418*si*si);
4967 si =
gv.bin[nd]->DustDftVel/
phycon.sqrte*1.551e-4;
4968 g0 = 1.5045*si*sqrt(1.+0.4418*si*si);
4969 g2 = si/(1.329 +
POW3(si));
4979 corr = sqrt(volmom/fdrag);
4980 gv.bin[nd]->DustDftVel = (
realnum)(vdold*corr);
4985 gv.lgNegGrnDrg =
true;
4986 gv.bin[nd]->DustDftVel = 0.;
4991 fprintf(
ioQQQ,
" %2ld new drift velocity:%10.2e momentum absorbed:%10.2e\n",
4992 loop,
gv.bin[nd]->DustDftVel, volmom );
5029 return max(1.e-10,GrnVryDpth_v);
const int FILENAME_PATH_LENGTH_2
sys_float sexp(sys_float x)
double fudge(long int ipnt)
sys_float safe_div(sys_float x, sys_float y, sys_float res_0by0)
const char * strstr_s(const char *haystack, const char *needle)
bool fp_equal(sys_float x, sys_float y, int n=3)
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
NORETURN void TotalInsanity(void)
sys_float SDIV(sys_float x)
#define DEBUG_ENTRY(funcname)
double powi(double, long int)
vector< unsigned int > nData
vector< vector< double > > Energy
vector< double > IonThres
vector< vector< double > > AvNumber
flex_arr< double > cs_pdt
long RecomZ0[LIMELM][LIMELM+1]
flex_arr< realnum > yhat_primary
realnum RecomEn[LIMELM][LIMELM+1]
realnum ChemEn[LIMELM][LIMELM+1]
vector< double > pure_sc1
double rate_h2_form_grains_used
bool lgPAHsInIonizedRegion
double rate_h2_form_grains_CT02
vector< realnum > inv_att_len
double rate_h2_form_grains_HM79
ial_type which_ial[MAT_TOP]
vector< realnum > GraphiteEmission
strg_type which_strg[MAT_TOP]
vector< realnum > GrainEmission
H2_type which_H2distr[MAT_TOP]
pot_type which_pot[MAT_TOP]
zmin_type which_zmin[MAT_TOP]
enth_type which_enth[MAT_TOP]
realnum dstAbundThresholdNear
vector< string > ReadRecord
realnum GrainHeatScaleFactor
pe_type which_pe[MAT_TOP]
AEInfo * AugerData[LIMELM]
realnum GrainChTrRate[LIMELM][LIMELM+1][LIMELM+1]
vector< realnum > SilicateEmission
realnum dstAbundThresholdFar
vector< flex_arr< realnum > > y01A
void reserve(size_type size)
void realloc(size_type end)
void alloc(size_type begin, size_type end)
double phfit(long int nz, long int ne, long int is, double e)
long ipoint(double energy_ryd)
void ConvFail(const char chMode[], const char chDetail[])
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
realnum GetAveVelocity(realnum massAMU)
t_elementnames elementnames
STATIC void NewChargeData(long)
double y2s(double, double, long, double *)
STATIC double ThetaNu(double)
void Yfunc(long, long, double, double, double, double, double, double *, double *, double *, double *)
static const long T_LOOP_MAX
STATIC void InitBinAugerData(size_t, long, long)
STATIC double GrnVryDpth(size_t)
double elec_esc_length(double e, long nd)
static const double STICK_ELEC
static const bool INCL_TUNNEL
void GrainRestartIter(void)
STATIC void UpdatePot2(size_t, long)
static double HEAT_TOLER_BIN
STATIC void PE_init(size_t, long, long, double *, double *, double *, double *, double *, double *, double *)
STATIC void GrainUpdateRadius1(void)
STATIC double GrainElecRecomb1(size_t, long, double *, double *)
STATIC void GetFracPop(size_t, long, double[], double[], long *)
STATIC double y1psa(size_t, long, double)
static const long MAGIC_AUGER_DATA
STATIC double y0psa(size_t, long, long, double)
STATIC void GrainCollHeating(size_t, realnum *, realnum *)
STATIC void GrainChargeTemp(void)
static const double STICK_ION
STATIC double GrnStdDpth(long)
STATIC void UpdatePot(size_t, long, long, double[], double[])
STATIC void ReadAugerData()
STATIC void GrainTemperature(size_t, realnum *, double *, double *, double *)
STATIC double PlanckIntegral(double, size_t, long)
STATIC void GrainCharge(size_t, double *)
void SetNChrgStates(long nChrg)
static long int nCalledGrainDrive
static const long CT_LOOP_MAX
STATIC void GetNextLine(const char *, FILE *, char[])
STATIC void UpdatePot1(size_t, long, long, long)
static const double TOLER
double y2pa(double, double, long, double *)
STATIC void GrainIonColl(size_t, long, long, long, const double[], const double[], long *, realnum *, realnum *)
static const double THERMCONST
double chrg2pot(double x, long nd)
static const double ETILDE
STATIC long HighestIonStage(void)
STATIC void UpdateRecomZ0(size_t, long, bool)
STATIC double y0b01(size_t, long, long)
STATIC double y0b(size_t, long, long)
static const bool ALL_STAGES
STATIC void GrainChrgTransferRates(long)
STATIC double GrainElecEmis1(size_t, long, double *, double *, double *, double *)
double pot2chrg(double x, long nd)
void GrainStartIter(void)
STATIC void GrainScreen(long, size_t, long, double *, double *)
STATIC void InitEmissivities(void)
static const long BRACKET_MAX
STATIC void GrainUpdateRadius2()
static const bool NO_TUNNEL
STATIC void GetPotValues(size_t, long, double *, double *, double *, double *, double *, double *, bool)
t_iso_sp iso_sp[NISO][LIMELM]
molezone * findspecieslocal(const char buf[])
UNUSED const double FR1RYD
UNUSED const double HPLANCK
UNUSED const double SPEEDLIGHT
UNUSED const double BOLTZMANN
UNUSED const double ELECTRON_MASS
UNUSED const double LN_TWO
UNUSED const double EN1RYD
UNUSED const double EVRYD
UNUSED const double EN1EV
UNUSED const double SQRT2
UNUSED const double ATOMIC_MASS_UNIT
UNUSED const double ELEM_CHARGE
UNUSED const double TE1RYD
long hunt_bisect(const T x[], long n, T xval)
void splint_safe(const double xa[], const double ya[], const double y2a[], long int n, double x, double *y, bool *lgOutOfBounds)
void spline(const double x[], const double y[], long int n, double yp1, double ypn, double y2a[])
void splint(const double xa[], const double ya[], const double y2a[], long int n, double x, double *y)