124static const double tlim[5]={0.,50.,150.,500.,DBL_MAX};
125static const double ppower[4]={2.00,1.30,0.68,0.00};
137 const vector<double>&,vector<double>&,
138 vector<double>&,vector<double>&,
142STATIC double TryDoubleStep(vector<double>&,vector<double>&,vector<double>&,vector<double>&,
143 vector<double>&,
const vector<double>&,
const vector<double>&,
double,
144 double*,
long,
size_t,
bool*);
148STATIC void ScanProbDistr(
const vector<double>&,
const vector<double>&,
long,
double,
long,
double,
149 long*,
long*,
long*,
long*,
QH_Code*);
152 vector<double>&,vector<double>&,vector<double>&,vector<double>&,
QH_Code*);
155 vector<double>&,
double*,
180 bool lgNoTdustFailures;
187 const double LIM1 = pow(2.e-6,1./2.);
188 const double LIM2 = pow(6.e-6,1./3.);
195 x = log(0.999*DBL_MAX);
198 for( i=0; i <
rfield.nflux; i++ )
201 gv.GrainEmission[i] = 0.;
202 gv.SilicateEmission[i] = 0.;
203 gv.GraphiteEmission[i] = 0.;
206 vector<double> qtemp(
NQGRID);
207 vector<double> qprob(
NQGRID);
209 for(
size_t nd=0; nd <
gv.bin.size(); nd++ )
217 lgLocalQHeat =
gv.bin[nd]->lgQHeat;
222 gv.dstAbundThresholdNear :
gv.dstAbundThresholdFar;
224 if( lgLocalQHeat &&
gv.bin[nd]->dstAbund >= threshold )
226 qheat(qtemp,qprob,&qnbin,nd);
228 if(
gv.bin[nd]->lgUseQHeat )
236 gv.bin[nd]->lgUseQHeat =
false;
245 if( lgLocalQHeat &&
gv.bin[nd]->lgUseQHeat )
250 for( j=qnbin-1; j >= 0 && xx > flux*DBL_EPSILON; j-- )
260 bfac = (1. + 0.5*f)*f;
263 xx = qprob[j]*
gv.bin[nd]->dstab1[i]*
rfield.anu2[i]*
282 bfac = (1. + 0.5*f)*f;
285 flux =
gv.bin[nd]->dstab1[i]*
rfield.anu2[i]*
rfield.widflx[i]/bfac;
290 flux *= factor*
gv.bin[nd]->cnv_H_pCM3;
366 lgNoTdustFailures =
true;
367 for(
size_t nd=0; nd <
gv.bin.size(); nd++ )
369 if( !
gv.bin[nd]->lgTdustConverged )
371 lgNoTdustFailures =
false;
378 for( i=0; i <
rfield.nflux; i++ )
383 for(
size_t nd=0; nd <
gv.bin.size(); nd++ )
385 if(
gv.bin[nd]->tedust <
gv.bin[nd]->Tsublimat )
395 ASSERT( fabs(BolFlux-
gv.GrainHeatSum) < Comparison1 );
398 for(
size_t nd=0; nd <
gv.bin.size(); nd++ )
400 if(
gv.bin[nd]->lgChrgConverged )
402 double ave = 0.5*(
gv.bin[nd]->RateDn+
gv.bin[nd]->RateUp);
409 if( lgNoTdustFailures &&
gv.lgDHetOn &&
gv.lgDColOn &&
thermal.ConstGrainTemp == 0. )
414 for(
size_t nd=0; nd <
gv.bin.size(); nd++ )
416 Comparison1 +=
gv.bin[nd]->BolFlux;
419 Comparison1 +=
MAX2(
gv.GasCoolColl,0.);
421 Comparison1 +=
gv.GrainHeatChem;
429 fabs(Comparison1-Comparison2)/Comparison2 <
CONSERV_TOL );
451 vector<double>& qprob,
495 fprintf(
ioQQQ,
"\n >>>>qheat called for grain %s\n",
gv.bin[nd]->chDstLab );
500 vector<double> phiTilde(
rfield.nupper);
501 vector<double> Phi(
rfield.nupper);
502 vector<double> PhiDrv(
rfield.nupper);
503 vector<double> dPdlnT(
NQGRID);
507 check +=
gv.bin[nd]->GrainHeatColl-
gv.bin[nd]->GrainCoolTherm;
517 for( i=
gv.bin[nd]->qnflux-1; i >= 0; i-- )
524 double fac = ( i <
gv.bin[nd]->qnflux-1 ) ? 1./
rfield.widflx[i] : 0.;
526 y = phiTilde[i]*
gv.bin[nd]->cnv_H_pGR*fac;
528 PhiDrv[i] = -0.5*(y + oldy);
536 integral += phiTilde[i]*
gv.bin[nd]->cnv_H_pCM3*
rfield.anu[i]*
EN1RYD;
540 c0 += Phi[i]*
rfield.widflx[i];
544 lgNegRate = lgNegRate || ( phiTilde[i] < 0. );
557 sprintf(fnam,
"Lambda_%2.2ld.asc",nd);
558 file = fopen(fnam,
"w");
559 for( i=0; i <
NDEMS; ++i )
560 fprintf(file,
"%e %e %e\n",
562 ufunct(exp(
gv.dsttmp[i]),nd,&lgBoundErr),
563 exp(
gv.bin[nd]->dstems[i])*
gv.bin[nd]->cnv_H_pGR/
EN1RYD);
566 sprintf(fnam,
"Phi_%2.2ld.asc",nd);
567 file = fopen(fnam,
"w");
568 for( i=0; i <
gv.bin[nd]->qnflux; ++i )
569 fprintf(file,
"%e %e\n",
rfield.anu[i],Phi[i]);
575 Tmax =
gv.bin[nd]->tedust;
577 Umax =
ufunct(Tmax,nd,&lgBoundErr);
583 deriv = y*c0/(
uderiv(Tmax,nd)*Tmax);
585 fwhm = sqrt(8.*
LN_TWO*c1/deriv);
588 FwhmRatio = fwhm/Umax;
607 ErrorCode = ErrorStart;
612 for(
int nz=0; nz <
gv.bin[nd]->nChrg; nz++ )
613 Rate2 +=
gv.bin[nd]->chrg[nz]->FracPop*
gv.bin[nd]->chrg[nz]->HeatingRate2;
615 fprintf(
ioQQQ,
" grain heating: %.4e, integral %.4e, total rate %.4e lgNegRate %c\n",
616 gv.bin[nd]->GrainHeat,integral,Phi[0],
TorF(lgNegRate));
617 fprintf(
ioQQQ,
" av grain temp %.4e av grain enthalpy (Ryd) %.4e\n",
618 gv.bin[nd]->tedust,Umax);
619 fprintf(
ioQQQ,
" fwhm^2/(4ln2*c2/c0): %.4e fwhm (Ryd) %.4e fwhm/Umax %.4e\n",
620 NumEvents,fwhm,FwhmRatio );
621 fprintf(
ioQQQ,
" HeatingRate1 %.4e HeatingRate2 %.4e lgQHTooWide %c\n",
622 gv.bin[nd]->HeatingRate1*
gv.bin[nd]->cnv_H_pCM3, Rate2*
gv.bin[nd]->cnv_H_pCM3,
623 TorF(
gv.bin[nd]->lgQHTooWide) );
628 maxBracket =
gv.bin[nd]->tedust;
647 for( i=0; i <
LOOP_MAX && !lgOK && !
gv.bin[nd]->lgQHTooWide; i++ )
649 if(
gv.bin[nd]->qtmin >=
gv.bin[nd]->tedust )
654 double MinEnth = exp(
gv.bin[nd]->DustEnth[0]);
655 Ulo =
MAX2(Ulo,MinEnth);
659 if(
gv.bin[nd]->qtmin <= minBracket ||
gv.bin[nd]->qtmin >= maxBracket )
660 gv.bin[nd]->qtmin = sqrt(minBracket*maxBracket);
664 ASSERT( minBracket <=
gv.bin[nd]->qtmin &&
gv.bin[nd]->qtmin <= maxBracket );
666 ErrorCode = ErrorStart;
671 GetProbDistr_LowLimit(nd,rel_tol,Umax,fwhm,Phi,PhiDrv,qtemp,qprob,dPdlnT,qnbin,
672 &new_tmin,&nWideFail,&ErrorCode);
688 GetProbDistr_HighLimit(nd,
STRICT,Umax,fwhm,qtemp,qprob,dPdlnT,&tol,qnbin,&dummy,
701 if( new_tmin < minBracket || new_tmin > maxBracket )
706 if( new_tmin <= minBracket )
707 new_tmin = sqrt(
gv.bin[nd]->qtmin*minBracket);
708 if( new_tmin >= maxBracket )
709 new_tmin = sqrt(
gv.bin[nd]->qtmin*maxBracket);
723 double help1 =
gv.bin[nd]->qtmin*sqrt(DefFac);
724 double help2 = sqrt(
gv.bin[nd]->qtmin*maxBracket);
725 minBracket =
gv.bin[nd]->qtmin;
726 new_tmin =
MIN2(help1,help2);
732 double help = sqrt(
gv.bin[nd]->qtmin*minBracket);
733 maxBracket =
gv.bin[nd]->qtmin;
738 new_tmin = sqrt(minBracket*maxBracket);
743 GetProbDistr_HighLimit(nd,
STRICT,Umax,fwhm,qtemp,qprob,dPdlnT,&tol,qnbin,&new_tmin,
747 gv.bin[nd]->qtmin = new_tmin;
759 fprintf(
ioQQQ,
" GetProbDistr returns code %d\n", ErrorCode );
762 fprintf(
ioQQQ,
" >>>>qheat trying another iteration, qtmin bracket %.4e %.4e",
763 minBracket,maxBracket );
764 fprintf(
ioQQQ,
" nWideFail %ld\n", nWideFail );
770 gv.bin[nd]->lgQHTooWide =
true;
774 if(
gv.bin[nd]->lgQHTooWide && !lgDelta )
792 if( !lgOK && !lgDelta )
794 double Umax2 = Umax*sqrt(tol);
795 double fwhm2 = fwhm*sqrt(tol);
802 GetProbDistr_HighLimit(nd,
RELAX,Umax2,fwhm2,qtemp,qprob,dPdlnT,&tol,qnbin,&dummy,
808 gv.bin[nd]->qtmin = dummy;
809 ErrorCode = ErrorCode2;
821 gv.bin[nd]->qtmin_zone1 =
gv.bin[nd]->qtmin;
823 gv.bin[nd]->lgUseQHeat = ( lgOK && !lgDelta );
824 gv.bin[nd]->lgEverQHeat = (
gv.bin[nd]->lgEverQHeat ||
gv.bin[nd]->lgUseQHeat );
829 fprintf(
ioQQQ,
" >>>>qheat converged with code: %d\n", ErrorCode );
834 ++
gv.bin[nd]->QHeatFailures;
835 fprintf(
ioQQQ,
" PROBLEM qheat did not converge grain %s in zone %ld, error code = %d\n",
836 gv.bin[nd]->chDstLab,
nzone,ErrorCode );
839 if(
gv.QHSaveFile != NULL && (
iterations.lgLastIt || !
gv.lgQHPunLast ) )
841 fprintf(
gv.QHSaveFile,
"\nDust Temperature Distribution: grain %s zone %ld\n",
844 fprintf(
gv.QHSaveFile,
"Equilibrium temperature: %.2f\n",
gv.bin[nd]->tedust );
846 if(
gv.bin[nd]->lgUseQHeat )
849 fprintf(
gv.QHSaveFile,
"Number of bins: %ld\n", *qnbin );
850 fprintf(
gv.QHSaveFile,
" Tgrain dP/dlnT\n" );
851 for( i=0; i < *qnbin; i++ )
853 fprintf(
gv.QHSaveFile,
"%.4e %.4e\n", qtemp[i],dPdlnT[i] );
858 fprintf(
gv.QHSaveFile,
"**** quantum heating was not used\n" );
867 vector<double>& phiTilde,
875 enum {DEBUG_LOC=
false};
888 for( i=0; i <
gv.bin[nd]->qnflux; i++ )
902 double NegHeatingRate = 0.;
904 for( nz=0; nz <
gv.bin[nd]->nChrg; nz++ )
912 check1 +=
rfield.SummedCon[i]*
gv.bin[nd]->dstab1[i]*
rfield.anu[i];
913 phiTilde[i] += gptr->
FracPop*
rfield.SummedCon[i]*
gv.bin[nd]->dstab1[i];
921 double cs1 = ( i >= ipLo2 ) ?
gv.bin[nd]->dstab1[i]*gptr->
yhat_primary[i] : 0.;
938 double ehat = gptr->
ehat[i];
939 double cool1,
sign = 1.;
942 if( gptr->
DustZ <= -1 )
963 double corr = xx/
rfield.anu[ipLo];
977 double integral = 0.;
978 for( i=0; i <
gv.bin[nd]->qnflux; i++ )
980 integral += phiTilde[i]*
gv.bin[nd]->cnv_H_pCM3*
rfield.anu[i]*
EN1RYD;
982 dprintf(
ioQQQ,
" integral test 1: integral %.6e %.6e\n", integral, sum );
1000 double Sum,ESum,DSum,Delta,E_av2,Corr;
1019 double ylo = -exp(-E0/fac);
1023 double Ehi =
gv.anumax[
gv.bin[nd]->qnflux-1]-Einf;
1024 double yhi = ((E0-Ehi)/fac-1.)*exp(-Ehi/fac);
1032 vector<double> RateArr(
gv.bin[nd]->qnflux);
1033 Sum = ESum = DSum = 0.;
1036 for( i=0; i <
gv.bin[nd]->qnflux2; i++ )
1038 Ehi =
gv.anumax[i] - Einf;
1042 yhi = ((E0-Ehi)/fac-1.)*exp(-Ehi/fac);
1044 RateArr[i] = rate*
MAX2(yhi-ylo,0.);
1046 ESum +=
rfield.anu[i]*RateArr[i];
1048 DSum +=
rfield.widflx[i]*RateArr[i];
1059 ASSERT( fabs(E_av-E_av2) <= Delta );
1062 for( i=0; i <
gv.bin[nd]->qnflux2; i++ )
1064 phiTilde[i] += RateArr[i]*Corr;
1071 double integral = 0.;
1072 for( i=0; i <
gv.bin[nd]->qnflux; i++ )
1074 integral += phiTilde[i]*
gv.bin[nd]->cnv_H_pCM3*
rfield.anu[i]*
EN1RYD;
1076 dprintf(
ioQQQ,
" integral test 2: integral %.6e %.6e\n", integral, sum );
1103 if(
gv.bin[nd]->HeatingRate1*
gv.bin[nd]->cnv_H_pCM3 > 0.05*
CONSERV_TOL*
gv.bin[nd]->GrainHeat )
1107 const double LIM2 = pow(3.e-6,1./3.);
1108 const double LIM3 = pow(8.e-6,1./4.);
1117 double rate =
gv.bin[nd]->HeatingRate1/E_av;
1123 double Ehi =
gv.anumax[
gv.bin[nd]->qnflux-1];
1124 double yhi = -(Ehi/fac+1.)*exp(-Ehi/fac);
1128 for( i=0; i <
gv.bin[nd]->qnflux2; i++ )
1133 double x =
gv.anumax[i]/fac;
1137 yhi = -(x+1.)*exp(-x);
1139 yhi = -(((1./3.)*x - 0.5)*x*x + 1.);
1141 yhi = -(1. - 0.5*x*x);
1144 phiTilde[i] += rate*
MAX2(yhi-ylo,0.);
1148 sum +=
gv.bin[nd]->HeatingRate1*
gv.bin[nd]->cnv_H_pCM3;
1152 double integral = 0.;
1153 for( i=0; i <
gv.bin[nd]->qnflux; i++ )
1155 integral += phiTilde[i]*
gv.bin[nd]->cnv_H_pCM3*
rfield.anu[i]*
EN1RYD;
1157 dprintf(
ioQQQ,
" integral test 3: integral %.6e %.6e\n", integral, sum );
1162 NegHeatingRate +=
gv.bin[nd]->HeatingRate1*
gv.bin[nd]->cnv_H_pCM3;
1168 if( NegHeatingRate < 0. )
1170 double scale_fac = (sum+NegHeatingRate)/sum;
1171 for( i=0; i <
gv.bin[nd]->qnflux; i++ )
1172 phiTilde[i] *= scale_fac;
1174 sum += NegHeatingRate;
1178 double integral = 0.;
1179 for( i=0; i <
gv.bin[nd]->qnflux; i++ )
1181 integral += phiTilde[i]*
gv.bin[nd]->cnv_H_pCM3*
rfield.anu[i]*
EN1RYD;
1183 dprintf(
ioQQQ,
" integral test 4: integral %.6e %.6e\n", integral, sum );
1208 const vector<double>& Phi,
1209 const vector<double>& PhiDrv,
1210 vector<double>& qtemp,
1211 vector<double>& qprob,
1212 vector<double>& dPdlnT,
1240 vector<double> delu(
NQGRID);
1241 vector<double> Lambda(
NQGRID);
1242 vector<double> p(
NQGRID);
1243 vector<double> u1(
NQGRID);
1253 fprintf(
ioQQQ,
" GetProbDistr_LowLimit called for grain %s\n",
gv.bin[nd]->chDstLab );
1254 fprintf(
ioQQQ,
" got qtmin1 %.4e\n",
gv.bin[nd]->qtmin);
1257 qtmin1 =
gv.bin[nd]->qtmin;
1260 u1[0] =
ufunct(qtemp[0],nd,&lgBoundErr);
1267 Lambda[0] = exp(y)*
gv.bin[nd]->cnv_H_pGR/
EN1RYD;
1271 delu[0] = 2.*Lambda[0]/Phi[0];
1273 dPdlnT[0] = p[0]*qtemp[0]*
uderiv(qtemp[0],nd);
1274 RadCooling = 0.5*p[0]*Lambda[0]*delu[0];
1275 NextStep = 0.01*Lambda[0]/Phi[0];
1277 if( NextStep < 0.05*DBL_EPSILON*u1[0] )
1282 else if( NextStep < 5.*DBL_EPSILON*u1[0] )
1293 lgAllNegSlope =
true;
1302 dlnpdlnU = u1[0]*Phi[0]/Lambda[0] - dlnLdlnT*u1[0]/(qtemp[0]*
uderiv(qtemp[0],nd));
1306 (void)
TryDoubleStep(u1,delu,p,qtemp,Lambda,Phi,PhiDrv,NextStep,&dCool,k,nd,&lgBoundErr);
1307 dPdlnT[2] = p[2]*qtemp[2]*
uderiv(qtemp[2],nd);
1309 if( dPdlnT[2] < dPdlnT[0] )
1321 double RelError =
TryDoubleStep(u1,delu,p,qtemp,Lambda,Phi,PhiDrv,NextStep,&dCool,k,nd,&lgBoundErr);
1338 NextStep *= sqrt(0.9*rel_tol*
QHEAT_TOL/RelError);
1341 if( NextStep < 5.*DBL_EPSILON*u1[k] )
1356 NextStep =
MIN2(NextStep,Lambda[k]/Phi[0]);
1359 dPdlnT[k-1] = p[k-1]*qtemp[k-1]*
uderiv(qtemp[k-1],nd);
1360 dPdlnT[k] = p[k]*qtemp[k]*
uderiv(qtemp[k],nd);
1362 lgAllNegSlope = lgAllNegSlope && ( dPdlnT[k] < dPdlnT[k-2] );
1364 if( dPdlnT[k-1] > maxVal )
1366 maxVal = dPdlnT[k-1];
1369 if( dPdlnT[k] > maxVal )
1375 RadCooling += dCool;
1385 if( p[k] > sqrt(DBL_MAX/100.) )
1393 for( j=0; j <= k; j++ )
1395 p[j] /= DBL_MAX/100.;
1396 dPdlnT[j] /= DBL_MAX/100.;
1398 RadCooling /= DBL_MAX/100.;
1404 ASSERT( p[k] > 0. && dPdlnT[k] > 0. && RadCooling > 0. );
1407 if( k > 0 && k%
NQTEST == 0 )
1409 double wid, avStep, factor;
1422 avStep = log(u1[k]/u1[0])/(double)k;
1426 factor = 1.1 + 3.9*(1.0 - sqrt((
double)k/(
double)
NQGRID));
1427 if( wid/avStep > factor*(
double)
NQGRID )
1443 fac = RadCooling*
gv.bin[nd]->cnv_GR_pCM3*
EN1RYD/
gv.bin[nd]->GrainHeat;
1464 ScanProbDistr(u1,dPdlnT,nbin,maxVal,nmax,qtmin1,&nstart,&nstart2,&nend,nWideFail,ErrorCode);
1471 for( j=0; j < nbin; j++ )
1480 for( j=nstart; j < nbin; j++ )
1484 *new_tmin = qtemp[j];
1498 if( dPdlnT[nstart2] < 0.999*dPdlnT[nstart2+
NSTARTUP] )
1503 double T1 = qtemp[nstart2];
1504 double T2 = qtemp[nstart2+
NSTARTUP];
1505 double delta_y = log(dPdlnT[nstart2+
NSTARTUP]/dPdlnT[nstart2]);
1506 double c2 = delta_y/(1./
POW3(T1)-1./
POW3(T2));
1508 *new_tmin = pow(c2/help,1./3.);
1514 double delta_x = log(qtemp[nstart2+
NSTARTUP]/qtemp[nstart2]);
1515 double delta_y = log(dPdlnT[nstart2+
NSTARTUP]/dPdlnT[nstart2]);
1517 *new_tmin = qtemp[nstart2]*exp(delta_x);
1518 if( *new_tmin < qtmin1 )
1520 *new_tmin = sqrt( *new_tmin*qtmin1 );
1531 ASSERT( *new_tmin <
gv.bin[nd]->tedust );
1567 nbin =
RebinQHeatResults(nd,nstart,nend,p,qtemp,qprob,dPdlnT,u1,delu,Lambda,ErrorCode);
1578 for( j=0; j < nbin; j++ )
1594 " zone %ld %s nbin %ld nok %ld nbad %ld total prob %.4e rel_tol %.3e new_tmin %.4e\n",
1595 nzone,
gv.bin[nd]->chDstLab,nbin,nok,nbad,sum,rel_tol,*new_tmin );
1606 vector<double>& delu,
1608 vector<double>& qtemp,
1609 vector<double>& Lambda,
1610 const vector<double>& Phi,
1611 const vector<double>& PhiDrv,
1641 ASSERT( index >= 0 && index <
NQGRID-2 && nd <
gv.bin.size() && step > 0. );
1646 uhi =
rfield.anu[
gv.bin[nd]->qnflux-1];
1650 for( i=0; i <= index; i++ )
1651 p_max =
MAX2(p_max,p[i]);
1656 for( i=1; i <= 2; i++ )
1660 long ipHi =
gv.bin[nd]->qnflux-1;
1663 u1[k] = u1[k-1] + delu[k];
1666 *lgBoundFail = *lgBoundFail || lgErr;
1667 Lambda[k] = exp(y)*
gv.bin[nd]->cnv_H_pGR/
EN1RYD;
1670 trap1 = trap2 = trap12 = 0.;
1673 for( j=jlo; j < k; j++ )
1675 umin = u1[k] - u1[j];
1685 else if( umin > ulo )
1707 while( ipHi-ipLo > 1 )
1709 long ipMd = (ipLo+ipHi)/2;
1715 bval_jk = Phi[ipLo] + (umin -
rfield.anu[ipLo])*PhiDrv[ipLo];
1733 trap2 = p[j]*bval_jk;
1735 sum += (trap1 + trap2)*delu[j];
1744 p[k] = (sum + trap1*delu[k])/(2.*Lambda[k] - Phi[0]*delu[k]);
1752 p2k = (sum2 + trap12*step)/(2.*Lambda[k] - Phi[0]*step);
1758 RelErrPk = fabs(p2k-p[k])/p[k];
1762 *cooling =
log_integral(u1[k-2],p[k-2]*Lambda[k-2],u1[k-1],p[k-1]*Lambda[k-1]);
1763 *cooling +=
log_integral(u1[k-1],p[k-1]*Lambda[k-1],u1[k],p[k]*Lambda[k]);
1766 cooling2 =
log_integral(u1[k-2],p[k-2]*Lambda[k-2],u1[k],p2k*Lambda[k]);
1769 RelErrCool = ( index > 0 ) ? fabs(cooling2-(*cooling))/(*cooling) : 0.;
1773 return MAX2(RelErrPk,RelErrCool)/3.;
1790 ASSERT( xx1 > 0. && yy1 > 0. && xx2 > 0. && yy2 > 0. );
1793 eps = log((xx2/xx1)*(yy2/yy1));
1794 if( fabs(eps) < 1.e-4 )
1796 result = xx1*yy1*xx*((eps/6. + 0.5)*eps + 1.);
1800 result = (xx2*yy2 - xx1*yy1)*xx/eps;
1808 const vector<double>& dPdlnT,
1824 const double MIN_FAC_LO = 1.e4;
1825 const double MIN_FAC_HI = 1.e4;
1830 ASSERT( nbin > 0 && nmax >= 0 && nmax < nbin && maxVal > 0. );
1840 for( i=nmax; i >= 0; --i )
1842 if( dPdlnT[i] < minVal )
1850 for( i=nmax; i > *nstart; --i )
1852 double deriv = log(dPdlnT[i]/dPdlnT[i-1])/log(u1[i]/u1[i-1]);
1853 if( deriv > deriv_max )
1862 lgSetLo = ( nmax >= *nend || maxVal/dPdlnT[*nend] < MIN_FAC_HI );
1868 lgSetHi = ( nmax <= *nstart || maxVal/dPdlnT[*nstart] < MIN_FAC_LO );
1870 lgSetHi = ( nmax <= *nstart2 || maxVal/dPdlnT[*nstart2] < MIN_FAC_LO );
1872 if( lgSetLo && lgSetHi )
1896 fprintf(
ioQQQ,
" ScanProbDistr nstart %ld nstart2 %ld nend %ld nmax %ld maxVal %.3e",
1897 *nstart,*nstart2,*nend,nmax,maxVal );
1898 fprintf(
ioQQQ,
" dPdlnT[nstart] %.3e dPdlnT[nstart2] %.3e dPdlnT[nend] %.3e code %d\n",
1899 dPdlnT[*nstart],dPdlnT[*nstart2],dPdlnT[*nend],*ErrorCode );
1917 vector<double>& qtemp,
1918 vector<double>& qprob,
1919 vector<double>& dPdlnT,
1921 vector<double>& delu,
1922 vector<double>& Lambda,
1944 ASSERT( nstart >= 0 && nstart < nend && nend <
NQGRID );
1957 vector<double> new_delu(
NQGRID);
1958 vector<double> new_dPdlnT(
NQGRID);
1959 vector<double> new_Lambda(
NQGRID);
1960 vector<double> new_p(
NQGRID);
1961 vector<double> new_qprob(
NQGRID);
1962 vector<double> new_qtemp(
NQGRID);
1963 vector<double> new_u1(
NQGRID);
1972 help = pow(qtemp[nend]/qtemp[i],1./(1.5*
NQMIN));
1991 double p1,p2,L1,L2,frac,slope;
1992 if( qtemp[i] <= T1 && T1 <= qtemp[i+1] )
1996 double xrlog = log(qtemp[i+1]/qtemp[i]);
1999 frac = log(T1/qtemp[i]);
2000 slope = log(p[i+1]/p[i])/xrlog;
2001 p1 = p[i]*exp(frac*slope);
2002 slope = log(Lambda[i+1]/Lambda[i])/xrlog;
2003 L1 = Lambda[i]*exp(frac*slope);
2008 p1 = sqrt(p[i]*p[i+1]);
2009 L1 = sqrt(Lambda[i]*Lambda[i+1]);
2019 if( qtemp[i] <= T2 && T2 <= qtemp[i+1] )
2022 help =
ufunct(T2,nd,&lgBoundErr);
2023 uu2 =
MIN2(help,u1[i+1]);
2025 double xrlog = log(qtemp[i+1]/qtemp[i]);
2028 frac = log(T2/qtemp[i]);
2029 slope = log(p[i+1]/p[i])/xrlog;
2030 p2 = p[i]*exp(frac*slope);
2031 slope = log(Lambda[i+1]/Lambda[i])/xrlog;
2032 L2 = Lambda[i]*exp(frac*slope);
2037 p2 = sqrt(p[i]*p[i+1]);
2038 L2 = sqrt(Lambda[i]*Lambda[i+1]);
2064 }
while( i < nend && ! lgDone );
2075 new_qprob[newnbin] = s0;
2076 new_Lambda[newnbin] = s1/s0;
2077 xx = log(new_Lambda[newnbin]*
EN1RYD*
gv.bin[nd]->cnv_GR_pH);
2080 new_qtemp[newnbin] = exp(y);
2081 new_u1[newnbin] =
ufunct(new_qtemp[newnbin],nd,&lgBoundErr);
2083 new_delu[newnbin] = wid;
2084 new_p[newnbin] = new_qprob[newnbin]/new_delu[newnbin];
2085 new_dPdlnT[newnbin] = new_p[newnbin]*new_qtemp[newnbin]*
uderiv(new_qtemp[newnbin],nd);
2088 RadCooling += new_qprob[newnbin]*new_Lambda[newnbin];
2096 if( newnbin <
NQMIN )
2102 fac = RadCooling*
EN1RYD*
gv.bin[nd]->cnv_GR_pCM3/
gv.bin[nd]->GrainHeat;
2106 fprintf(
ioQQQ,
" RebinQHeatResults found tol1 %.4e tol2 %.4e\n",
2107 fabs(u1[nend]/Ucheck-1.), fabs(fac-1.) );
2112 ASSERT( fabs(u1[nend]/Ucheck-1.) < 10.*sqrt((
double)(nend-nstart+newnbin))*DBL_EPSILON );
2117 for( i=0; i < newnbin; i++ )
2120 p[i] = new_p[i]/fac;
2121 qtemp[i] = new_qtemp[i];
2122 qprob[i] = new_qprob[i]/fac;
2123 dPdlnT[i] = new_dPdlnT[i]/fac;
2125 delu[i] = new_delu[i];
2126 Lambda[i] = new_Lambda[i];
2129 ASSERT( qtemp[i] > 0. && qprob[i] > 0. );
2142 vector<double>& qtemp,
2143 vector<double>& qprob,
2144 vector<double>& dPdlnT,
2186 fprintf(
ioQQQ,
" GetProbDistr_HighLimit called for grain %s\n",
gv.bin[nd]->chDstLab );
2194 help1 = Umax*exp(-fac1);
2195 help2 = exp(
gv.bin[nd]->DustEnth[0]);
2196 Ulo =
MAX2(help1,help2);
2202 if( fac2 > log(DBL_MAX/10.) )
2207 Uhi = Umax*exp(fac2);
2213 uu1 =
ufunct(T1,nd,&lgErr);
2214 lgBoundErr = lgBoundErr || lgErr;
2215 help1 = log(uu1/Umax);
2216 p1 = c1*exp(-c2*
POW2(help1));
2218 lgBoundErr = lgBoundErr || lgErr;
2219 L1 = exp(y)*
gv.bin[nd]->cnv_H_pGR/
EN1RYD;
2222 if( uu1*p1*L1 < 1.e5*DBL_MIN )
2229 help1 = pow(Thi/Tlo,1./(1.2*
NQMIN));
2240 uu2 =
ufunct(T2,nd,&lgErr);
2241 lgBoundErr = lgBoundErr || lgErr;
2242 help1 = log(uu2/Umax);
2243 p2 = c1*exp(-c2*
POW2(help1));
2245 lgBoundErr = lgBoundErr || lgErr;
2246 L2 = exp(y)*
gv.bin[nd]->cnv_H_pGR/
EN1RYD;
2253 Lambda[nbin] = s1/s0;
2254 xx = log(Lambda[nbin]*
EN1RYD*
gv.bin[nd]->cnv_GR_pH);
2256 lgBoundErr = lgBoundErr || lgErr;
2257 qtemp[nbin] = exp(y);
2259 p[nbin] = qprob[nbin]/delu[nbin];
2260 dPdlnT[nbin] = p[nbin]*qtemp[nbin]*
uderiv(qtemp[nbin],nd);
2263 RadCooling += qprob[nbin]*Lambda[nbin];
2272 }
while( T2 < Thi && nbin <
NQGRID );
2274 fac = RadCooling*
EN1RYD*
gv.bin[nd]->cnv_GR_pCM3/
gv.bin[nd]->GrainHeat;
2276 for( i=0; i < nbin; ++i )
2284 *new_tmin = qtemp[0];
2302 fprintf(
ioQQQ,
" GetProbDistr_HighLimit found tol1 %.4e tol2 %.4e\n",
2303 fabs(sum-1.), fabs(sum/fac-1.) );
2304 fprintf(
ioQQQ,
" zone %ld %s nbin %ld total prob %.4e new_tmin %.4e\n",
2305 nzone,
gv.bin[nd]->chDstLab,nbin,sum/fac,*new_tmin );
2321 hok[3] = {1275., 1670., 4359.},
2333 fprintf(
ioQQQ,
" uderiv called with non-positive temperature: %.6e\n" , temp );
2338 ecase =
gv.which_enth[
gv.bin[nd]->matType];
2342 numer = (4.15e-22/
EN1RYD)*pow(temp,3.3);
2343 dnumer = (3.3*4.15e-22/
EN1RYD)*pow(temp,2.3);
2344 denom = 1. + 6.51e-03*temp + 1.5e-06*temp*temp + 8.3e-07*pow(temp,2.3);
2345 ddenom = 6.51e-03 + 2.*1.5e-06*temp + 2.3*8.3e-07*pow(temp,1.3);
2348 deriv = (dnumer*denom - numer*ddenom)/
POW2(denom);
2359 for( j = 0; j < 4; j++ )
2361 if( temp >
tlim[j] && temp <=
tlim[j+1] )
2377 x = log10(
MIN2(temp,2000.));
2378 deriv = pow(10.,-21.26+3.1688*x-0.401894*
POW2(x))/
EN1RYD;
2390 else if( N_C <= 100. )
2392 N_H = 2.5*sqrt(N_C);
2399 for( i=0; i < 3; i++ )
2401 double help1 = hok[i]/temp;
2404 double help2 = exp(help1);
2405 double help3 = ( help1 < 1.e-7 ) ? help1*(1.+0.5*help1) : help2-1.;
2412 fprintf(
ioQQQ,
" uderiv called with unknown type for enthalpy function: %d\n", ecase );
2422 fprintf(
ioQQQ,
" uderiv finds non-positive derivative: %.6e, what's up?\n" , deriv );
2441 fprintf(
ioQQQ,
" ufunct called with non-positive temperature: %.6e\n" , temp );
2465 if( enthalpy <= 0. )
2467 fprintf(
ioQQQ,
" inv_ufunct called with non-positive enthalpy: %.6e\n" , enthalpy );
2491 for(
size_t nd=0; nd <
gv.bin.size(); nd++ )
2494 C_V2 =
uderiv(tdust2,nd);
2496 gv.bin[nd]->DustEnth[0] = C_V2*tdust2/4.;
2500 for(
long i=1; i <
NDEMS; i++ )
2503 tdust2 = exp(
gv.dsttmp[i]);
2504 C_V2 =
uderiv(tdust2,nd);
2505 tmid = sqrt(tdust1*tdust2);
2507 for(
long j=1; j < 4; j++ )
2509 if( tdust1 <
tlim[j] &&
tlim[j] < tdust2 )
2516 gv.bin[nd]->DustEnth[i] =
gv.bin[nd]->DustEnth[i-1] +
2525 for(
size_t nd=0; nd <
gv.bin.size(); nd++ )
2527 for(
long i=0; i <
NDEMS; i++ )
2529 gv.bin[nd]->DustEnth[i] = log(
gv.bin[nd]->DustEnth[i]);
2534 for(
size_t nd=0; nd <
gv.bin.size(); nd++ )
2555 ASSERT( n == 2 || n == 3 );
2562 res = 6.*1.202056903159594*
POW2(x);
2566 res = 24.*1.082323233711138*
POW3(x);
2575 nn = 4*
MAX2(4,2*(
long)(0.05/x));
2576 vector<double> xx(nn);
2577 vector<double> rr(nn);
2578 vector<double> aa(nn);
2579 vector<double> ww(nn);
2584 for( i=0; i < nn; i++ )
2586 double help1 = rr[i]/x;
2589 double help2 = exp(help1);
2590 double help3 = ( help1 < 1.e-7 ) ? help1*(1.+0.5*help1) : help2-1.;
2591 res += ww[i]*
powi(rr[i],n+1)*help2/
POW2(help3);
2596 return (
double)n*res;
int dprintf(FILE *fp, const char *format,...)
NORETURN void TotalInsanity(void)
#define DEBUG_ENTRY(funcname)
double powi(double, long int)
flex_arr< realnum > yhat_primary
void gauss_legendre(long int, vector< double > &, vector< double > &)
void gauss_init(long int, double, double, const vector< double > &, const vector< double > &, vector< double > &, vector< double > &)
static const double PROB_CUTOFF_HI
static const long WIDE_FAIL_MAX
static const long NSTARTUP
static const double MAX_EVENTS
static const double FWHM_RATIO2
static const double QHEAT_TOL
void GrainMakeDiffuse(void)
static const long LOOP_MAX
STATIC double ufunct(double, size_t, bool *)
static const double QT_RATIO
static const double SAFETY
STATIC void GetProbDistr_LowLimit(size_t, double, double, double, const vector< double > &, const vector< double > &, vector< double > &, vector< double > &, vector< double > &, long *, double *, long *, QH_Code *)
static const double DEN_SIL
STATIC void ScanProbDistr(const vector< double > &, const vector< double > &, long, double, long, double, long *, long *, long *, long *, QH_Code *)
STATIC double inv_ufunct(double, size_t, bool *)
double no_atoms(size_t nd)
static const double PROB_TOL
static const double DEF_FAC
static const double MW_SIL
STATIC double log_integral(double, double, double, double)
static const double ppower[4]
STATIC double DebyeDeriv(double, long)
STATIC long RebinQHeatResults(size_t, long, long, vector< double > &, vector< double > &, vector< double > &, vector< double > &, vector< double > &, vector< double > &, vector< double > &, QH_Code *)
static const double STRICT
static const double PROB_CUTOFF_LO
static const double FWHM_RATIO
STATIC double uderiv(double, size_t)
STATIC double TryDoubleStep(vector< double > &, vector< double > &, vector< double > &, vector< double > &, vector< double > &, const vector< double > &, const vector< double > &, double, double *, long, size_t, bool *)
STATIC void GetProbDistr_HighLimit(long, double, double, double, vector< double > &, vector< double > &, vector< double > &, double *, long *, double *, QH_Code *)
static const double RELAX
static const double tlim[5]
static const long MAX_LOOP
static const double cval[4]
void qheat(vector< double > &qtemp, vector< double > &qprob, long int *qnbin, size_t nd)
STATIC void qheat_init(size_t, vector< double > &, double *)
UNUSED const double FR1RYD
UNUSED const double SPEEDLIGHT
UNUSED const double BOLTZMANN
UNUSED const double LN_TWO
UNUSED const double EN1RYD
UNUSED const double ATOMIC_MASS_UNIT
UNUSED const double TE1RYD
void spldrv_safe(const double xa[], const double ya[], const double y2a[], long int n, double x, double *y, bool *lgOutOfBounds)
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[])