00001
00002
00003
00004 #include "cddefines.h"
00005 #include "physconst.h"
00006 #include "rfield.h"
00007 #include "phycon.h"
00008 #include "dense.h"
00009 #include "hmi.h"
00010 #include "thermal.h"
00011 #include "trace.h"
00012 #include "thirdparty.h"
00013 #include "iterations.h"
00014 #include "grainvar.h"
00015 #include "grains.h"
00016
00017 inline double no_atoms(size_t nd)
00018 {
00019 return gv.bin[nd]->AvVol*gv.bin[nd]->dustp[0]/ATOMIC_MASS_UNIT/gv.bin[nd]->atomWeight;
00020 }
00021
00022
00023
00024
00025
00026
00027 typedef enum {
00028
00029
00030 QH_OK, QH_ANALYTIC, QH_ANALYTIC_RELAX, QH_DELTA,
00031
00032
00033 QH_NEGRATE_FAIL, QH_LOOP_FAIL, QH_ARRAY_FAIL, QH_THIGH_FAIL,
00034
00035
00036 QH_RETRY, QH_CONV_FAIL, QH_BOUND_FAIL, QH_DELTA_FAIL,
00037
00038
00039
00040 QH_NO_REBIN, QH_LOW_FAIL, QH_HIGH_FAIL, QH_STEP_FAIL,
00041
00042
00043 QH_FATAL, QH_WIDE_FAIL, QH_NBIN_FAIL, QH_REBIN_FAIL
00044 } QH_Code;
00045
00046
00047
00048
00049
00050 static const long NQMIN = 10L;
00051
00052
00053 static const double PROB_CUTOFF_LO = 1.e-15;
00054 static const double PROB_CUTOFF_HI = 1.e-20;
00055 static const double SAFETY = 1.e+8;
00056
00057
00058 static const long NSTARTUP = 5L;
00059
00060
00061
00062 static const double MAX_EVENTS = 150.;
00063
00064
00065
00066 static const long LOOP_MAX = 20L;
00067
00068
00069 static const double DEF_FAC = 3.;
00070
00071
00072 static const double PROB_TOL = 0.02;
00073
00074
00075
00076 static const long NQTEST = 500L;
00077
00078
00079
00080 static const double FWHM_RATIO = 0.1;
00081
00082
00083 static const double FWHM_RATIO2 = 0.007;
00084
00085
00086 static const long MAX_LOOP = 2*NQGRID;
00087
00088
00089 static const double QHEAT_TOL = 5.e-3;
00090
00091
00092 static const long WIDE_FAIL_MAX = 3;
00093
00094
00095 static const double STRICT = 1.;
00096 static const double RELAX = 15.;
00097
00098
00099
00100
00101
00102
00103 static const double QT_RATIO = 1.03;
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118 static const double DEN_SIL = 3.30;
00119
00120
00121 static const double MW_SIL = 24.6051;
00122
00123
00124 static const double tlim[5]={0.,50.,150.,500.,DBL_MAX};
00125 static const double ppower[4]={2.00,1.30,0.68,0.00};
00126 static const double cval[4]={
00127 1.40e3/DEN_SIL*ATOMIC_MASS_UNIT*MW_SIL/EN1RYD,
00128 2.20e4/DEN_SIL*ATOMIC_MASS_UNIT*MW_SIL/EN1RYD,
00129 4.80e5/DEN_SIL*ATOMIC_MASS_UNIT*MW_SIL/EN1RYD,
00130 3.41e7/DEN_SIL*ATOMIC_MASS_UNIT*MW_SIL/EN1RYD};
00131
00132
00133
00134 STATIC void qheat_init(size_t,vector<double>&,double*);
00135
00136 STATIC void GetProbDistr_LowLimit(size_t,double,double,double,const vector<double>&,
00137 const vector<double>&,vector<double>&,
00138 vector<double>&,vector<double>&,
00139 long*,double*,long*,QH_Code*);
00140
00141
00142 STATIC double TryDoubleStep(vector<double>&,vector<double>&,vector<double>&,vector<double>&,
00143 vector<double>&,const vector<double>&,const vector<double>&,double,
00144 double*,long,size_t,bool*);
00145
00146 STATIC double log_integral(double,double,double,double);
00147
00148 STATIC void ScanProbDistr(const vector<double>&,const vector<double>&,long,double,long,double,
00149 long*,long*,long*,long*,QH_Code*);
00150
00151 STATIC long RebinQHeatResults(size_t,long,long,vector<double>&,vector<double>&,vector<double>&,
00152 vector<double>&,vector<double>&,vector<double>&,vector<double>&,QH_Code*);
00153
00154 STATIC void GetProbDistr_HighLimit(long,double,double,double,vector<double>&,vector<double>&,
00155 vector<double>&,double*,
00156 long*,double*,QH_Code*);
00157
00158 STATIC double uderiv(double,size_t);
00159
00160 STATIC double ufunct(double,size_t,bool*);
00161
00162 STATIC double inv_ufunct(double,size_t,bool*);
00163
00164 STATIC double DebyeDeriv(double,long);
00165
00166
00167
00168
00169
00170 void GrainMakeDiffuse(void)
00171 {
00172 long i, j;
00173 double bfac,
00174 f,
00175 factor,
00176 flux,
00177 x;
00178
00179 # ifndef NDEBUG
00180 bool lgNoTdustFailures;
00181 double BolFlux,
00182 Comparison1,
00183 Comparison2;
00184 # endif
00185
00186
00187 const double LIM1 = pow(2.e-6,1./2.);
00188 const double LIM2 = pow(6.e-6,1./3.);
00189
00190 DEBUG_ENTRY( "GrainMakeDiffuse()" );
00191
00192 factor = 2.*PI4*POW2(FR1RYD/SPEEDLIGHT)*FR1RYD;
00193
00194
00195 x = log(0.999*DBL_MAX);
00196
00197
00198 for( i=0; i < rfield.nflux; i++ )
00199 {
00200
00201 gv.GrainEmission[i] = 0.;
00202 gv.SilicateEmission[i] = 0.;
00203 gv.GraphiteEmission[i] = 0.;
00204 }
00205
00206 vector<double> qtemp(NQGRID);
00207 vector<double> qprob(NQGRID);
00208
00209 for( size_t nd=0; nd < gv.bin.size(); nd++ )
00210 {
00211 strg_type scase;
00212 bool lgLocalQHeat;
00213 long qnbin=-200;
00214 realnum threshold;
00215 vector<realnum> Invalid;
00216 vector<realnum>& grn = Invalid;
00217 double xx;
00218
00219
00220 scase = gv.which_strg[gv.bin[nd]->matType];
00221 switch( scase )
00222 {
00223 case STRG_SIL:
00224 grn = gv.SilicateEmission;
00225 break;
00226 case STRG_CAR:
00227 grn = gv.GraphiteEmission;
00228 break;
00229 default:
00230 fprintf( ioQQQ, " GrainMakeDiffuse called with unknown storage class: %d\n", scase );
00231 cdEXIT(EXIT_FAILURE);
00232 }
00233
00234
00235 lgLocalQHeat = gv.bin[nd]->lgQHeat;
00236
00237
00238
00239 threshold = ( dense.xIonDense[ipHYDROGEN][0]+dense.xIonDense[ipHYDROGEN][1] > hmi.H2_total ) ?
00240 gv.dstAbundThresholdNear : gv.dstAbundThresholdFar;
00241
00242 if( lgLocalQHeat && gv.bin[nd]->dstAbund >= threshold )
00243 {
00244 qheat(qtemp,qprob,&qnbin,nd);
00245
00246 if( gv.bin[nd]->lgUseQHeat )
00247 {
00248 ASSERT( qnbin > 0 );
00249 }
00250 }
00251 else
00252 {
00253
00254 gv.bin[nd]->lgUseQHeat = false;
00255 }
00256
00257 flux = 1.;
00258
00259
00260 for( i=0; i < rfield.nflux && (realnum)flux > 0.f; i++ )
00261 {
00262 flux = 0.;
00263 if( lgLocalQHeat && gv.bin[nd]->lgUseQHeat )
00264 {
00265 xx = 1.;
00266
00267
00268 for( j=qnbin-1; j >= 0 && xx > flux*DBL_EPSILON; j-- )
00269 {
00270 f = TE1RYD/qtemp[j]*rfield.anu[i];
00271 if( f < x )
00272 {
00273
00274
00275 if( f > LIM2 )
00276 bfac = exp(f) - 1.;
00277 else if( f > LIM1 )
00278 bfac = (1. + 0.5*f)*f;
00279 else
00280 bfac = f;
00281 xx = qprob[j]*gv.bin[nd]->dstab1[i]*rfield.anu2[i]*
00282 rfield.widflx[i]/bfac;
00283 flux += xx;
00284 }
00285 else
00286 {
00287 xx = 0.;
00288 }
00289 }
00290 }
00291 else
00292 {
00293 f = TE1RYD/gv.bin[nd]->tedust*rfield.anu[i];
00294 if( f < x )
00295 {
00296
00297 if( f > LIM2 )
00298 bfac = exp(f) - 1.;
00299 else if( f > LIM1 )
00300 bfac = (1. + 0.5*f)*f;
00301 else
00302 bfac = f;
00303 flux = gv.bin[nd]->dstab1[i]*rfield.anu2[i]*rfield.widflx[i]/bfac;
00304 }
00305 }
00306
00307
00308 flux *= factor*gv.bin[nd]->cnv_H_pCM3;
00309
00310
00311
00312 gv.GrainEmission[i] += (realnum)flux;
00313
00314 grn[i] += (realnum)flux;
00315 }
00316 }
00317
00318 # ifndef NDEBUG
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373 lgNoTdustFailures = true;
00374 for( size_t nd=0; nd < gv.bin.size(); nd++ )
00375 {
00376 if( !gv.bin[nd]->lgTdustConverged )
00377 {
00378 lgNoTdustFailures = false;
00379 break;
00380 }
00381 }
00382
00383
00384 BolFlux = 0.;
00385 for( i=0; i < rfield.nflux; i++ )
00386 {
00387 BolFlux += gv.GrainEmission[i]*rfield.anu[i]*EN1RYD;
00388 }
00389 Comparison1 = 0.;
00390 for( size_t nd=0; nd < gv.bin.size(); nd++ )
00391 {
00392 if( gv.bin[nd]->tedust < gv.bin[nd]->Tsublimat )
00393 Comparison1 += CONSERV_TOL*gv.bin[nd]->GrainHeat;
00394 else
00395
00396
00397 Comparison1 += 10.*CONSERV_TOL*gv.bin[nd]->GrainHeat;
00398 }
00399
00400
00401
00402 ASSERT( fabs(BolFlux-gv.GrainHeatSum) < Comparison1 );
00403
00404
00405 for( size_t nd=0; nd < gv.bin.size(); nd++ )
00406 {
00407 if( gv.bin[nd]->lgChrgConverged )
00408 {
00409 double ave = 0.5*(gv.bin[nd]->RateDn+gv.bin[nd]->RateUp);
00410
00411
00412 ASSERT( lgAbort || fabs(gv.bin[nd]->RateDn-gv.bin[nd]->RateUp) < CONSERV_TOL*ave );
00413 }
00414 }
00415
00416 if( lgNoTdustFailures && gv.lgDHetOn && gv.lgDColOn && thermal.ConstGrainTemp == 0. )
00417 {
00418
00419
00420 Comparison1 = 0.;
00421 for( size_t nd=0; nd < gv.bin.size(); nd++ )
00422 {
00423 Comparison1 += gv.bin[nd]->BolFlux;
00424 }
00425
00426 Comparison1 += MAX2(gv.GasCoolColl,0.);
00427
00428 Comparison1 += gv.GrainHeatChem;
00429
00430
00431 Comparison2 = gv.GrainHeatSum+thermal.heating[0][13]+thermal.heating[0][14]+thermal.heating[0][25];
00432
00433
00434
00435 ASSERT( lgAbort || gv.GrainHeatScaleFactor != 1.f || gv.lgBakesPAH_heat ||
00436 fabs(Comparison1-Comparison2)/Comparison2 < CONSERV_TOL );
00437 }
00438 # endif
00439 return;
00440 }
00441
00442
00443
00444
00445
00446
00447
00448
00449
00450
00451
00452
00453
00454
00455
00456
00457 void qheat( vector<double>& qtemp,
00458 vector<double>& qprob,
00459 long int *qnbin,
00460 size_t nd)
00461 {
00462 bool lgBoundErr,
00463 lgDelta,
00464 lgNegRate,
00465 lgOK,
00466 lgTried;
00467 long int i,
00468 nWideFail;
00469 QH_Code ErrorCode,
00470 ErrorCode2,
00471 ErrorStart;
00472 double c0,
00473 c1,
00474 c2,
00475 check,
00476 DefFac,
00477 deriv,
00478 fwhm,
00479 FwhmRatio,
00480 integral,
00481 minBracket,
00482 maxBracket,
00483 new_tmin,
00484 NumEvents,
00485 oldy,
00486 rel_tol,
00487 Tmax,
00488 tol,
00489 Umax,
00490 xx,
00491 y;
00492
00493
00494 DEBUG_ENTRY( "qheat()" );
00495
00496
00497 ASSERT( gv.bin[nd]->lgQHeat );
00498 ASSERT( nd < gv.bin.size() );
00499
00500 if( trace.lgTrace && trace.lgDustBug )
00501 {
00502 fprintf( ioQQQ, "\n >>>>qheat called for grain %s\n", gv.bin[nd]->chDstLab );
00503 }
00504
00505
00506
00507 vector<double> phiTilde(rfield.nupper);
00508 vector<double> Phi(rfield.nupper);
00509 vector<double> PhiDrv(rfield.nupper);
00510 vector<double> dPdlnT(NQGRID);
00511
00512 qheat_init( nd, phiTilde, &check );
00513
00514 check += gv.bin[nd]->GrainHeatColl-gv.bin[nd]->GrainCoolTherm;
00515
00516 xx = integral = 0.;
00517 c0 = c1 = c2 = 0.;
00518 lgNegRate = false;
00519 oldy = 0.;
00520 tol = 1.;
00521
00522
00523
00524 for( i=gv.bin[nd]->qnflux-1; i >= 0; i-- )
00525 {
00526
00527
00528
00529
00530
00531 double fac = ( i < gv.bin[nd]->qnflux-1 ) ? 1./rfield.widflx[i] : 0.;
00532
00533 y = phiTilde[i]*gv.bin[nd]->cnv_H_pGR*fac;
00534
00535 PhiDrv[i] = -0.5*(y + oldy);
00536
00537 xx -= PhiDrv[i]*(rfield.anu[i+1]-rfield.anu[i]);
00538
00539 Phi[i] = xx;
00540
00541 # ifndef NDEBUG
00542
00543 integral += phiTilde[i]*gv.bin[nd]->cnv_H_pCM3*rfield.anu[i]*EN1RYD;
00544 # endif
00545
00546
00547 c0 += Phi[i]*rfield.widflx[i];
00548 c1 += Phi[i]*rfield.anu[i]*rfield.widflx[i];
00549 c2 += Phi[i]*POW2(rfield.anu[i])*rfield.widflx[i];
00550
00551 lgNegRate = lgNegRate || ( phiTilde[i] < 0. );
00552
00553 oldy = y;
00554 }
00555
00556
00557 ASSERT( fabs(check-integral)/check <= CONSERV_TOL );
00558
00559 # if 0
00560 {
00561 char fnam[50];
00562 FILE *file;
00563
00564 sprintf(fnam,"Lambda_%2.2ld.asc",nd);
00565 file = fopen(fnam,"w");
00566 for( i=0; i < NDEMS; ++i )
00567 fprintf(file,"%e %e %e\n",
00568 exp(gv.dsttmp[i]),
00569 ufunct(exp(gv.dsttmp[i]),nd,&lgBoundErr),
00570 exp(gv.bin[nd]->dstems[i])*gv.bin[nd]->cnv_H_pGR/EN1RYD);
00571 fclose(file);
00572
00573 sprintf(fnam,"Phi_%2.2ld.asc",nd);
00574 file = fopen(fnam,"w");
00575 for( i=0; i < gv.bin[nd]->qnflux; ++i )
00576 fprintf(file,"%e %e\n", rfield.anu[i],Phi[i]);
00577 fclose(file);
00578 }
00579 # endif
00580
00581
00582 Tmax = gv.bin[nd]->tedust;
00583
00584 Umax = ufunct(Tmax,nd,&lgBoundErr);
00585 ASSERT( !lgBoundErr );
00586
00587 spldrv_safe(gv.dsttmp,gv.bin[nd]->dstems,gv.bin[nd]->dstslp2,NDEMS,log(Tmax),&y,&lgBoundErr);
00588 ASSERT( !lgBoundErr );
00589
00590 deriv = y*c0/(uderiv(Tmax,nd)*Tmax);
00591
00592 fwhm = sqrt(8.*LN_TWO*c1/deriv);
00593
00594 NumEvents = POW2(fwhm)*c0/(4.*LN_TWO*c2);
00595 FwhmRatio = fwhm/Umax;
00596
00597
00598 lgDelta = ( FwhmRatio < FWHM_RATIO2 );
00599
00600 lgOK = lgDelta;
00601
00602 ErrorStart = QH_OK;
00603
00604 if( lgDelta )
00605 {
00606
00607 lgNegRate = false;
00608 ErrorStart = MAX2(ErrorStart,QH_DELTA);
00609 }
00610
00611 if( lgNegRate )
00612 ErrorStart = MAX2(ErrorStart,QH_NEGRATE_FAIL);
00613
00614 ErrorCode = ErrorStart;
00615
00616 if( trace.lgTrace && trace.lgDustBug )
00617 {
00618 double Rate2 = 0.;
00619 for( int nz=0; nz < gv.bin[nd]->nChrg; nz++ )
00620 Rate2 += gv.bin[nd]->chrg[nz]->FracPop*gv.bin[nd]->chrg[nz]->HeatingRate2;
00621
00622 fprintf( ioQQQ, " grain heating: %.4e, integral %.4e, total rate %.4e lgNegRate %c\n",
00623 gv.bin[nd]->GrainHeat,integral,Phi[0],TorF(lgNegRate));
00624 fprintf( ioQQQ, " av grain temp %.4e av grain enthalpy (Ryd) %.4e\n",
00625 gv.bin[nd]->tedust,Umax);
00626 fprintf( ioQQQ, " fwhm^2/(4ln2*c2/c0): %.4e fwhm (Ryd) %.4e fwhm/Umax %.4e\n",
00627 NumEvents,fwhm,FwhmRatio );
00628 fprintf( ioQQQ, " HeatingRate1 %.4e HeatingRate2 %.4e lgQHTooWide %c\n",
00629 gv.bin[nd]->HeatingRate1*gv.bin[nd]->cnv_H_pCM3, Rate2*gv.bin[nd]->cnv_H_pCM3,
00630 TorF(gv.bin[nd]->lgQHTooWide) );
00631 }
00632
00633
00634 minBracket = GRAIN_TMIN;
00635 maxBracket = gv.bin[nd]->tedust;
00636
00637
00638 lgTried = false;
00639
00640 nWideFail = 0;
00641
00642 DefFac = DEF_FAC;
00643
00644 rel_tol = 1.;
00645
00646
00647
00648
00649
00650
00651
00652
00653
00654 for( i=0; i < LOOP_MAX && !lgOK && !gv.bin[nd]->lgQHTooWide; i++ )
00655 {
00656 if( gv.bin[nd]->qtmin >= gv.bin[nd]->tedust )
00657 {
00658
00659
00660 double Ulo = Umax*exp(-sqrt(-log(PROB_CUTOFF_LO)/(4.*LN_TWO))*fwhm/Umax);
00661 double MinEnth = exp(gv.bin[nd]->DustEnth[0]);
00662 Ulo = MAX2(Ulo,MinEnth);
00663 gv.bin[nd]->qtmin = inv_ufunct(Ulo,nd,&lgBoundErr);
00664 ASSERT( !lgBoundErr );
00665
00666 if( gv.bin[nd]->qtmin <= minBracket || gv.bin[nd]->qtmin >= maxBracket )
00667 gv.bin[nd]->qtmin = sqrt(minBracket*maxBracket);
00668 }
00669 gv.bin[nd]->qtmin = MAX2(gv.bin[nd]->qtmin,GRAIN_TMIN);
00670
00671 ASSERT( minBracket <= gv.bin[nd]->qtmin && gv.bin[nd]->qtmin <= maxBracket );
00672
00673 ErrorCode = ErrorStart;
00674
00675
00676 if( FwhmRatio >= FWHM_RATIO && NumEvents <= MAX_EVENTS )
00677 {
00678 GetProbDistr_LowLimit(nd,rel_tol,Umax,fwhm,Phi,PhiDrv,qtemp,qprob,dPdlnT,qnbin,
00679 &new_tmin,&nWideFail,&ErrorCode);
00680
00681
00682 if( ErrorCode == QH_DELTA_FAIL && fwhm < Umax && !lgTried )
00683 {
00684 double dummy;
00685
00686
00687
00688
00689
00690
00691
00692 ErrorCode = MAX2(ErrorStart,QH_ANALYTIC);
00693
00694
00695 GetProbDistr_HighLimit(nd,STRICT,Umax,fwhm,qtemp,qprob,dPdlnT,&tol,qnbin,&dummy,
00696 &ErrorCode);
00697
00698 if( ErrorCode >= QH_RETRY )
00699 {
00700 ErrorCode = QH_DELTA_FAIL;
00701 lgTried = true;
00702 }
00703 }
00704
00705
00706 if( ErrorCode < QH_NO_REBIN )
00707 {
00708 if( new_tmin < minBracket || new_tmin > maxBracket )
00709 ++nWideFail;
00710
00711 if( nWideFail < WIDE_FAIL_MAX )
00712 {
00713 if( new_tmin <= minBracket )
00714 new_tmin = sqrt(gv.bin[nd]->qtmin*minBracket);
00715 if( new_tmin >= maxBracket )
00716 new_tmin = sqrt(gv.bin[nd]->qtmin*maxBracket);
00717 }
00718 else
00719 {
00720 ErrorCode = MAX2(ErrorCode,QH_WIDE_FAIL);
00721 }
00722
00723 if( ErrorCode == QH_CONV_FAIL )
00724 {
00725 rel_tol *= 0.9;
00726 }
00727 }
00728 else if( ErrorCode == QH_LOW_FAIL )
00729 {
00730 double help1 = gv.bin[nd]->qtmin*sqrt(DefFac);
00731 double help2 = sqrt(gv.bin[nd]->qtmin*maxBracket);
00732 minBracket = gv.bin[nd]->qtmin;
00733 new_tmin = MIN2(help1,help2);
00734
00735 DefFac += 1.5;
00736 }
00737 else if( ErrorCode == QH_HIGH_FAIL )
00738 {
00739 double help = sqrt(gv.bin[nd]->qtmin*minBracket);
00740 maxBracket = gv.bin[nd]->qtmin;
00741 new_tmin = MAX2(gv.bin[nd]->qtmin/DEF_FAC,help);
00742 }
00743 else
00744 {
00745 new_tmin = sqrt(minBracket*maxBracket);
00746 }
00747 }
00748 else
00749 {
00750 GetProbDistr_HighLimit(nd,STRICT,Umax,fwhm,qtemp,qprob,dPdlnT,&tol,qnbin,&new_tmin,
00751 &ErrorCode);
00752 }
00753
00754 gv.bin[nd]->qtmin = new_tmin;
00755
00756 lgOK = ( ErrorCode < QH_RETRY );
00757
00758 if( ErrorCode >= QH_FATAL )
00759 break;
00760
00761 if( ErrorCode != QH_LOW_FAIL )
00762 DefFac = DEF_FAC;
00763
00764 if( trace.lgTrace && trace.lgDustBug )
00765 {
00766 fprintf( ioQQQ, " GetProbDistr returns code %d\n", ErrorCode );
00767 if( !lgOK )
00768 {
00769 fprintf( ioQQQ, " >>>>qheat trying another iteration, qtmin bracket %.4e %.4e",
00770 minBracket,maxBracket );
00771 fprintf( ioQQQ, " nWideFail %ld\n", nWideFail );
00772 }
00773 }
00774 }
00775
00776 if( ErrorCode == QH_WIDE_FAIL )
00777 gv.bin[nd]->lgQHTooWide = true;
00778
00779
00780
00781 if( gv.bin[nd]->lgQHTooWide && !lgDelta )
00782 ErrorCode = MAX2(ErrorCode,QH_WIDE_FAIL);
00783
00784
00785
00786
00787
00788
00789
00790
00791
00792
00793
00794
00795
00796
00797
00798
00799 if( !lgOK && !lgDelta )
00800 {
00801 double Umax2 = Umax*sqrt(tol);
00802 double fwhm2 = fwhm*sqrt(tol);
00803
00804 for( i=0; i < LOOP_MAX; ++i )
00805 {
00806 double dummy;
00807
00808 ErrorCode2 = MAX2(ErrorStart,QH_ANALYTIC);
00809 GetProbDistr_HighLimit(nd,RELAX,Umax2,fwhm2,qtemp,qprob,dPdlnT,&tol,qnbin,&dummy,
00810 &ErrorCode2);
00811
00812 lgOK = ( ErrorCode2 < QH_RETRY );
00813 if( lgOK )
00814 {
00815 gv.bin[nd]->qtmin = dummy;
00816 ErrorCode = ErrorCode2;
00817 break;
00818 }
00819 else
00820 {
00821 Umax2 *= sqrt(tol);
00822 fwhm2 *= sqrt(tol);
00823 }
00824 }
00825 }
00826
00827 if( nzone == 1 )
00828 gv.bin[nd]->qtmin_zone1 = gv.bin[nd]->qtmin;
00829
00830 gv.bin[nd]->lgUseQHeat = ( lgOK && !lgDelta );
00831 gv.bin[nd]->lgEverQHeat = ( gv.bin[nd]->lgEverQHeat || gv.bin[nd]->lgUseQHeat );
00832
00833 if( lgOK )
00834 {
00835 if( trace.lgTrace && trace.lgDustBug )
00836 fprintf( ioQQQ, " >>>>qheat converged with code: %d\n", ErrorCode );
00837 }
00838 else
00839 {
00840 *qnbin = 0;
00841 ++gv.bin[nd]->QHeatFailures;
00842 fprintf( ioQQQ, " PROBLEM qheat did not converge grain %s in zone %ld, error code = %d\n",
00843 gv.bin[nd]->chDstLab,nzone,ErrorCode );
00844 }
00845
00846 if( gv.QHSaveFile != NULL && ( iterations.lgLastIt || !gv.lgQHPunLast ) )
00847 {
00848 fprintf( gv.QHSaveFile, "\nDust Temperature Distribution: grain %s zone %ld\n",
00849 gv.bin[nd]->chDstLab,nzone );
00850
00851 fprintf( gv.QHSaveFile, "Equilibrium temperature: %.2f\n", gv.bin[nd]->tedust );
00852
00853 if( gv.bin[nd]->lgUseQHeat )
00854 {
00855
00856 fprintf( gv.QHSaveFile, "Number of bins: %ld\n", *qnbin );
00857 fprintf( gv.QHSaveFile, " Tgrain dP/dlnT\n" );
00858 for( i=0; i < *qnbin; i++ )
00859 {
00860 fprintf( gv.QHSaveFile, "%.4e %.4e\n", qtemp[i],dPdlnT[i] );
00861 }
00862 }
00863 else
00864 {
00865 fprintf( gv.QHSaveFile, "**** quantum heating was not used\n" );
00866 }
00867 }
00868 return;
00869 }
00870
00871
00872
00873 STATIC void qheat_init(size_t nd,
00874 vector<double>& phiTilde,
00875 double *check)
00876 {
00877 long i,
00878 nz;
00879 double sum = 0.;
00880
00881
00882 enum {DEBUG_LOC=false};
00883
00884
00885 DEBUG_ENTRY( "qheat_init()" );
00886
00887
00888 ASSERT( gv.bin[nd]->lgQHeat );
00889 ASSERT( nd < gv.bin.size() );
00890
00891 *check = 0.;
00892
00893
00894
00895 for( i=0; i < gv.bin[nd]->qnflux; i++ )
00896 {
00897 phiTilde[i] = 0.;
00898 }
00899
00900
00901
00902
00903
00904
00905
00906
00907
00908
00909 double NegHeatingRate = 0.;
00910
00911 for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
00912 {
00913 double check1 = 0.;
00914 ChargeBin *gptr = gv.bin[nd]->chrg[nz];
00915
00916
00917 for( i=0; i < MIN2(gptr->ipThresInf,rfield.nflux); i++ )
00918 {
00919 check1 += rfield.SummedCon[i]*gv.bin[nd]->dstab1[i]*rfield.anu[i];
00920 phiTilde[i] += gptr->FracPop*rfield.SummedCon[i]*gv.bin[nd]->dstab1[i];
00921 }
00922
00923
00924
00925 for( i=gptr->ipThresInf; i < rfield.nflux; i++ )
00926 {
00927 long ipLo2 = gptr->ipThresInfVal;
00928 double cs1 = ( i >= ipLo2 ) ? gv.bin[nd]->dstab1[i]*gptr->yhat_primary[i] : 0.;
00929
00930 check1 += rfield.SummedCon[i]*gptr->fac1[i];
00931
00932 phiTilde[i] += gptr->FracPop*rfield.SummedCon[i]*MAX2(gv.bin[nd]->dstab1[i]-cs1,0.);
00933
00934
00935
00936
00937 if( cs1 > 0. )
00938 {
00939
00940
00941
00942
00943 double ratio = ( gv.lgWD01 ) ? 1. : gptr->yhat[i]/gptr->yhat_primary[i];
00944
00945 double ehat = gptr->ehat[i];
00946 double cool1, sign = 1.;
00947 realnum xx;
00948
00949 if( gptr->DustZ <= -1 )
00950 cool1 = gptr->ThresSurf + gptr->PotSurf + ehat;
00951 else
00952 cool1 = gptr->ThresSurfVal + gptr->PotSurf + ehat;
00953
00954
00955
00956
00957
00958
00959
00960 xx = rfield.anu[i] - (realnum)(ratio*cool1);
00961 if( xx < 0.f )
00962 {
00963 xx = -xx;
00964 sign = -1.;
00965 }
00966 long ipLo = hunt_bisect( &gv.anumin[0], i+1, xx );
00967
00968
00969
00970 double corr = xx/rfield.anu[ipLo];
00971 phiTilde[ipLo] += sign*corr*gptr->FracPop*rfield.SummedCon[i]*cs1;
00972 }
00973
00974
00975
00976 }
00977
00978 *check += gptr->FracPop*check1*EN1RYD*gv.bin[nd]->cnv_H_pCM3;
00979
00980 sum += gptr->FracPop*check1*EN1RYD*gv.bin[nd]->cnv_H_pCM3;
00981
00982 if( DEBUG_LOC )
00983 {
00984 double integral = 0.;
00985 for( i=0; i < gv.bin[nd]->qnflux; i++ )
00986 {
00987 integral += phiTilde[i]*gv.bin[nd]->cnv_H_pCM3*rfield.anu[i]*EN1RYD;
00988 }
00989 dprintf( ioQQQ, " integral test 1: integral %.6e %.6e\n", integral, sum );
00990 }
00991
00992
00993
00994
00995
00996
00997
00998
00999
01000
01001
01002
01003
01004
01005 if( gptr->HeatingRate2*gv.bin[nd]->cnv_H_pCM3 > 0.05*CONSERV_TOL*gv.bin[nd]->GrainHeat )
01006 {
01007 double Sum,ESum,DSum,Delta,E_av2,Corr;
01008 double fac = BOLTZMANN/EN1RYD*phycon.te;
01009
01010
01011
01012 double E0 = -(MIN2(gptr->PotSurfInc,0.) + MIN2(gptr->ThresInfInc,0.));
01013
01014
01015
01016 double Einf = gptr->ThresInfInc + MIN2(gptr->PotSurfInc,0.);
01017
01018
01019
01020
01021
01022 double E_av = MAX2(gptr->ThresInfInc,0.)*EN1RYD + 2.*BOLTZMANN*phycon.te;
01023
01024 double rate = gptr->HeatingRate2/E_av;
01025
01026 double ylo = -exp(-E0/fac);
01027
01028
01029
01030 double Ehi = gv.anumax[gv.bin[nd]->qnflux-1]-Einf;
01031 double yhi = ((E0-Ehi)/fac-1.)*exp(-Ehi/fac);
01032
01033 rate /= yhi-ylo;
01034
01035
01036 rate *= gptr->FracPop;
01037
01038
01039 vector<double> RateArr(gv.bin[nd]->qnflux);
01040 Sum = ESum = DSum = 0.;
01041
01042
01043 for( i=0; i < gv.bin[nd]->qnflux2; i++ )
01044 {
01045 Ehi = gv.anumax[i] - Einf;
01046 if( Ehi >= E0 )
01047 {
01048
01049 yhi = ((E0-Ehi)/fac-1.)*exp(-Ehi/fac);
01050
01051 RateArr[i] = rate*MAX2(yhi-ylo,0.);
01052 Sum += RateArr[i];
01053 ESum += rfield.anu[i]*RateArr[i];
01054 # ifndef NDEBUG
01055 DSum += rfield.widflx[i]*RateArr[i];
01056 # endif
01057 ylo = yhi;
01058 }
01059 else
01060 {
01061 RateArr[i] = 0.;
01062 }
01063 }
01064 E_av2 = ESum/Sum*EN1RYD;
01065 Delta = DSum/Sum*EN1RYD;
01066 ASSERT( fabs(E_av-E_av2) <= Delta );
01067 Corr = E_av/E_av2;
01068
01069 for( i=0; i < gv.bin[nd]->qnflux2; i++ )
01070 {
01071 phiTilde[i] += RateArr[i]*Corr;
01072 }
01073
01074 sum += gptr->FracPop*gptr->HeatingRate2*gv.bin[nd]->cnv_H_pCM3;
01075
01076 if( DEBUG_LOC )
01077 {
01078 double integral = 0.;
01079 for( i=0; i < gv.bin[nd]->qnflux; i++ )
01080 {
01081 integral += phiTilde[i]*gv.bin[nd]->cnv_H_pCM3*rfield.anu[i]*EN1RYD;
01082 }
01083 dprintf( ioQQQ, " integral test 2: integral %.6e %.6e\n", integral, sum );
01084 }
01085 }
01086 else
01087 {
01088 NegHeatingRate += gptr->FracPop*gptr->HeatingRate2*gv.bin[nd]->cnv_H_pCM3;
01089 }
01090 }
01091
01092
01093
01094
01095
01096
01097
01098
01099
01100
01101
01102
01103
01104
01105
01106
01107
01108
01109
01110 if( gv.bin[nd]->HeatingRate1*gv.bin[nd]->cnv_H_pCM3 > 0.05*CONSERV_TOL*gv.bin[nd]->GrainHeat )
01111 {
01112
01113
01114 const double LIM2 = pow(3.e-6,1./3.);
01115 const double LIM3 = pow(8.e-6,1./4.);
01116
01117
01118
01119
01120 double fac = BOLTZMANN/EN1RYD*MAX2(phycon.te,gv.bin[nd]->tedust);
01121
01122 double E_av = 2.*BOLTZMANN*MAX2(phycon.te,gv.bin[nd]->tedust);
01123
01124 double rate = gv.bin[nd]->HeatingRate1/E_av;
01125
01126 double ylo = -1.;
01127
01128
01129
01130 double Ehi = gv.anumax[gv.bin[nd]->qnflux-1];
01131 double yhi = -(Ehi/fac+1.)*exp(-Ehi/fac);
01132
01133 rate /= yhi-ylo;
01134
01135 for( i=0; i < gv.bin[nd]->qnflux2; i++ )
01136 {
01137
01138
01139
01140 double x = gv.anumax[i]/fac;
01141
01142
01143 if( x > LIM3 )
01144 yhi = -(x+1.)*exp(-x);
01145 else if( x > LIM2 )
01146 yhi = -(((1./3.)*x - 0.5)*x*x + 1.);
01147 else
01148 yhi = -(1. - 0.5*x*x);
01149
01150
01151 phiTilde[i] += rate*MAX2(yhi-ylo,0.);
01152 ylo = yhi;
01153 }
01154
01155 sum += gv.bin[nd]->HeatingRate1*gv.bin[nd]->cnv_H_pCM3;
01156
01157 if( DEBUG_LOC )
01158 {
01159 double integral = 0.;
01160 for( i=0; i < gv.bin[nd]->qnflux; i++ )
01161 {
01162 integral += phiTilde[i]*gv.bin[nd]->cnv_H_pCM3*rfield.anu[i]*EN1RYD;
01163 }
01164 dprintf( ioQQQ, " integral test 3: integral %.6e %.6e\n", integral, sum );
01165 }
01166 }
01167 else
01168 {
01169 NegHeatingRate += gv.bin[nd]->HeatingRate1*gv.bin[nd]->cnv_H_pCM3;
01170 }
01171
01172
01173
01174
01175 if( NegHeatingRate < 0. )
01176 {
01177 double scale_fac = (sum+NegHeatingRate)/sum;
01178 for( i=0; i < gv.bin[nd]->qnflux; i++ )
01179 phiTilde[i] *= scale_fac;
01180
01181 sum += NegHeatingRate;
01182
01183 if( DEBUG_LOC )
01184 {
01185 double integral = 0.;
01186 for( i=0; i < gv.bin[nd]->qnflux; i++ )
01187 {
01188 integral += phiTilde[i]*gv.bin[nd]->cnv_H_pCM3*rfield.anu[i]*EN1RYD;
01189 }
01190 dprintf( ioQQQ, " integral test 4: integral %.6e %.6e\n", integral, sum );
01191 }
01192 }
01193
01194 return;
01195 }
01196
01197
01198
01199
01200
01201
01202
01203
01204
01205
01206
01207
01208
01209
01210
01211 STATIC void GetProbDistr_LowLimit(size_t nd,
01212 double rel_tol,
01213 double Umax,
01214 double fwhm,
01215 const vector<double>& Phi,
01216 const vector<double>& PhiDrv,
01217 vector<double>& qtemp,
01218 vector<double>& qprob,
01219 vector<double>& dPdlnT,
01220 long int *qnbin,
01221 double *new_tmin,
01222 long *nWideFail,
01223 QH_Code *ErrorCode)
01224 {
01225 bool lgAllNegSlope,
01226 lgBoundErr;
01227 long int j,
01228 k,
01229 l,
01230 nbad,
01231 nbin,
01232 nend=0,
01233 nmax,
01234 nok,
01235 nstart=0,
01236 nstart2=0;
01237 double dCool,
01238 dlnLdlnT,
01239 dlnpdlnU,
01240 fac = 0.,
01241 maxVal,
01242 NextStep,
01243 qtmin1,
01244 RadCooling,
01245 sum,
01246 y;
01247 vector<double> delu(NQGRID);
01248 vector<double> Lambda(NQGRID);
01249 vector<double> p(NQGRID);
01250 vector<double> u1(NQGRID);
01251
01252
01253 DEBUG_ENTRY( "GetProbDistr_LowLimit()" );
01254
01255
01256 ASSERT( nd < gv.bin.size() );
01257
01258 if( trace.lgTrace && trace.lgDustBug )
01259 {
01260 fprintf( ioQQQ, " GetProbDistr_LowLimit called for grain %s\n", gv.bin[nd]->chDstLab );
01261 fprintf( ioQQQ, " got qtmin1 %.4e\n", gv.bin[nd]->qtmin);
01262 }
01263
01264 qtmin1 = gv.bin[nd]->qtmin;
01265 qtemp[0] = qtmin1;
01266
01267 u1[0] = ufunct(qtemp[0],nd,&lgBoundErr);
01268 ASSERT( !lgBoundErr );
01269
01270
01271 splint_safe(gv.dsttmp,gv.bin[nd]->dstems,gv.bin[nd]->dstslp2,NDEMS,log(qtemp[0]),&y,&lgBoundErr);
01272 ASSERT( !lgBoundErr );
01273
01274 Lambda[0] = exp(y)*gv.bin[nd]->cnv_H_pGR/EN1RYD;
01275
01276
01277
01278 delu[0] = 2.*Lambda[0]/Phi[0];
01279 p[0] = PROB_CUTOFF_LO;
01280 dPdlnT[0] = p[0]*qtemp[0]*uderiv(qtemp[0],nd);
01281 RadCooling = 0.5*p[0]*Lambda[0]*delu[0];
01282 NextStep = 0.01*Lambda[0]/Phi[0];
01283
01284 if( NextStep < 0.05*DBL_EPSILON*u1[0] )
01285 {
01286 *ErrorCode = MAX2(*ErrorCode,QH_STEP_FAIL);
01287 return;
01288 }
01289 else if( NextStep < 5.*DBL_EPSILON*u1[0] )
01290 {
01291
01292 NextStep *= 100.;
01293 }
01294
01295 nbad = 0;
01296 k = 0;
01297
01298 *qnbin = 0;
01299 *new_tmin = qtmin1;
01300 lgAllNegSlope = true;
01301 maxVal = dPdlnT[0];
01302 nmax = 0;
01303
01304
01305
01306
01307 spldrv_safe(gv.dsttmp,gv.bin[nd]->dstems,gv.bin[nd]->dstslp2,NDEMS,log(qtemp[0]),&dlnLdlnT,&lgBoundErr);
01308 ASSERT( !lgBoundErr );
01309 dlnpdlnU = u1[0]*Phi[0]/Lambda[0] - dlnLdlnT*u1[0]/(qtemp[0]*uderiv(qtemp[0],nd));
01310 if( dlnpdlnU < 0. )
01311 {
01312
01313 (void)TryDoubleStep(u1,delu,p,qtemp,Lambda,Phi,PhiDrv,NextStep,&dCool,k,nd,&lgBoundErr);
01314 dPdlnT[2] = p[2]*qtemp[2]*uderiv(qtemp[2],nd);
01315
01316 if( dPdlnT[2] < dPdlnT[0] )
01317 {
01318
01319
01320 *ErrorCode = MAX2(*ErrorCode,QH_HIGH_FAIL);
01321 return;
01322 }
01323 }
01324
01325
01326 for( l=0; l < MAX_LOOP; ++l )
01327 {
01328 double RelError = TryDoubleStep(u1,delu,p,qtemp,Lambda,Phi,PhiDrv,NextStep,&dCool,k,nd,&lgBoundErr);
01329
01330
01331
01332 if( lgBoundErr )
01333 {
01334 nbad += 2;
01335 *ErrorCode = MAX2(*ErrorCode,QH_THIGH_FAIL);
01336 break;
01337 }
01338
01339
01340 if( RelError > rel_tol*QHEAT_TOL )
01341 {
01342 nbad += 2;
01343
01344
01345 NextStep *= sqrt(0.9*rel_tol*QHEAT_TOL/RelError);
01346
01347
01348 if( NextStep < 5.*DBL_EPSILON*u1[k] )
01349 {
01350 *ErrorCode = MAX2(*ErrorCode,QH_STEP_FAIL);
01351 break;
01352 }
01353
01354 continue;
01355 }
01356 else
01357 {
01358
01359 k += 2;
01360
01361
01362 NextStep *= MIN2(pow(0.9*rel_tol*QHEAT_TOL/MAX2(RelError,1.e-50),1./3.),4.);
01363 NextStep = MIN2(NextStep,Lambda[k]/Phi[0]);
01364 }
01365
01366 dPdlnT[k-1] = p[k-1]*qtemp[k-1]*uderiv(qtemp[k-1],nd);
01367 dPdlnT[k] = p[k]*qtemp[k]*uderiv(qtemp[k],nd);
01368
01369 lgAllNegSlope = lgAllNegSlope && ( dPdlnT[k] < dPdlnT[k-2] );
01370
01371 if( dPdlnT[k-1] > maxVal )
01372 {
01373 maxVal = dPdlnT[k-1];
01374 nmax = k-1;
01375 }
01376 if( dPdlnT[k] > maxVal )
01377 {
01378 maxVal = dPdlnT[k];
01379 nmax = k;
01380 }
01381
01382 RadCooling += dCool;
01383
01384
01385
01386
01387
01388
01389
01390
01391
01392 if( p[k] > sqrt(DBL_MAX/100.) )
01393 {
01394 *ErrorCode = MAX2(*ErrorCode,QH_LOW_FAIL);
01395 break;
01396
01397 #if 0
01398
01399
01400 for( j=0; j <= k; j++ )
01401 {
01402 p[j] /= DBL_MAX/100.;
01403 dPdlnT[j] /= DBL_MAX/100.;
01404 }
01405 RadCooling /= DBL_MAX/100.;
01406 #endif
01407 }
01408
01409
01410
01411 ASSERT( p[k] > 0. && dPdlnT[k] > 0. && RadCooling > 0. );
01412
01413
01414 if( k > 0 && k%NQTEST == 0 )
01415 {
01416 double wid, avStep, factor;
01417
01418
01419 if( lgAllNegSlope )
01420 {
01421 *ErrorCode = MAX2(*ErrorCode,QH_HIGH_FAIL);
01422 break;
01423 }
01424
01425
01426
01427 wid = (sqrt(-log(PROB_CUTOFF_LO)/(4.*LN_TWO)) +
01428 sqrt(-log(PROB_CUTOFF_HI)/(4.*LN_TWO)))*fwhm/Umax;
01429 avStep = log(u1[k]/u1[0])/(double)k;
01430
01431
01432
01433 factor = 1.1 + 3.9*(1.0 - sqrt((double)k/(double)NQGRID));
01434 if( wid/avStep > factor*(double)NQGRID )
01435 {
01436 *ErrorCode = MAX2(*ErrorCode,QH_ARRAY_FAIL);
01437 break;
01438 }
01439 }
01440
01441
01442
01443 if( k >= NQGRID-2 )
01444 {
01445 *ErrorCode = MAX2(*ErrorCode,QH_ARRAY_FAIL);
01446 break;
01447 }
01448
01449
01450 fac = RadCooling*gv.bin[nd]->cnv_GR_pCM3*EN1RYD/gv.bin[nd]->GrainHeat;
01451
01452
01453 if( dPdlnT[k] < dPdlnT[k-1] && dPdlnT[k]/fac < PROB_CUTOFF_HI )
01454 {
01455 break;
01456 }
01457 }
01458
01459 if( l == MAX_LOOP )
01460 *ErrorCode = MAX2(*ErrorCode,QH_LOOP_FAIL);
01461
01462 nok = k;
01463 nbin = k+1;
01464
01465
01466 if( *ErrorCode < QH_NO_REBIN && nbin < NQMIN )
01467 *ErrorCode = MAX2(*ErrorCode,QH_NBIN_FAIL);
01468
01469
01470 if( *ErrorCode < QH_NO_REBIN )
01471 ScanProbDistr(u1,dPdlnT,nbin,maxVal,nmax,qtmin1,&nstart,&nstart2,&nend,nWideFail,ErrorCode);
01472
01473 if( *ErrorCode >= QH_NO_REBIN )
01474 {
01475 return;
01476 }
01477
01478 for( j=0; j < nbin; j++ )
01479 {
01480 p[j] /= fac;
01481 dPdlnT[j] /= fac;
01482 }
01483 RadCooling /= fac;
01484
01485
01486 *new_tmin = 0.;
01487 for( j=nstart; j < nbin; j++ )
01488 {
01489 if( dPdlnT[j] < PROB_CUTOFF_LO )
01490 {
01491 *new_tmin = qtemp[j];
01492 }
01493 else
01494 {
01495 if( j == nstart )
01496 {
01497
01498
01499
01500
01501 if( dPdlnT[j] > SAFETY*PROB_CUTOFF_LO && qtmin1 > 1.2*GRAIN_TMIN )
01502 *ErrorCode = MAX2(*ErrorCode,QH_BOUND_FAIL);
01503
01504
01505 if( dPdlnT[nstart2] < 0.999*dPdlnT[nstart2+NSTARTUP] )
01506 {
01507
01508
01509
01510 double T1 = qtemp[nstart2];
01511 double T2 = qtemp[nstart2+NSTARTUP];
01512 double delta_y = log(dPdlnT[nstart2+NSTARTUP]/dPdlnT[nstart2]);
01513 double c2 = delta_y/(1./POW3(T1)-1./POW3(T2));
01514 double help = c2/POW3(T1) + log(dPdlnT[nstart2]/PROB_CUTOFF_LO);
01515 *new_tmin = pow(c2/help,1./3.);
01516 }
01517
01518
01519 if( dPdlnT[j] > SAFETY*PROB_CUTOFF_LO && *new_tmin >= qtmin1 )
01520 {
01521 double delta_x = log(qtemp[nstart2+NSTARTUP]/qtemp[nstart2]);
01522 double delta_y = log(dPdlnT[nstart2+NSTARTUP]/dPdlnT[nstart2]);
01523 delta_x *= log(PROB_CUTOFF_LO/dPdlnT[nstart2])/delta_y;
01524 *new_tmin = qtemp[nstart2]*exp(delta_x);
01525 if( *new_tmin < qtmin1 )
01526
01527 *new_tmin = sqrt( *new_tmin*qtmin1 );
01528 else
01529
01530 *new_tmin = qtmin1/DEF_FAC;
01531 }
01532 }
01533 break;
01534 }
01535 }
01536 *new_tmin = MAX3(*new_tmin,qtmin1/DEF_FAC,GRAIN_TMIN);
01537
01538 ASSERT( *new_tmin < gv.bin[nd]->tedust );
01539
01540
01541 if( dPdlnT[nbin-1] > SAFETY*PROB_CUTOFF_HI )
01542 {
01543 if( *ErrorCode == QH_ARRAY_FAIL || *ErrorCode == QH_LOOP_FAIL )
01544 {
01545 ++(*nWideFail);
01546
01547 if( *nWideFail < WIDE_FAIL_MAX )
01548 {
01549
01550
01551 *ErrorCode = MAX2(*ErrorCode,QH_DELTA_FAIL);
01552 }
01553 else
01554 {
01555 *ErrorCode = MAX2(*ErrorCode,QH_WIDE_FAIL);
01556 }
01557 }
01558 else
01559 {
01560 *ErrorCode = MAX2(*ErrorCode,QH_BOUND_FAIL);
01561 }
01562 }
01563
01564
01565
01566
01567
01568
01569
01570
01571
01572
01573
01574 nbin = RebinQHeatResults(nd,nstart,nend,p,qtemp,qprob,dPdlnT,u1,delu,Lambda,ErrorCode);
01575
01576
01577 if( nbin == 0 )
01578 {
01579 return;
01580 }
01581
01582 *qnbin = nbin;
01583
01584 sum = 0.;
01585 for( j=0; j < nbin; j++ )
01586 {
01587 sum += qprob[j];
01588 }
01589
01590
01591
01592
01593
01594
01595 if( fabs(sum-1.) > PROB_TOL && qtmin1 > 1.2*GRAIN_TMIN )
01596 *ErrorCode = MAX2(*ErrorCode,QH_CONV_FAIL);
01597
01598 if( trace.lgTrace && trace.lgDustBug )
01599 {
01600 fprintf( ioQQQ,
01601 " zone %ld %s nbin %ld nok %ld nbad %ld total prob %.4e rel_tol %.3e new_tmin %.4e\n",
01602 nzone,gv.bin[nd]->chDstLab,nbin,nok,nbad,sum,rel_tol,*new_tmin );
01603 }
01604 return;
01605 }
01606
01607
01608
01609
01610
01611
01612 STATIC double TryDoubleStep(vector<double>& u1,
01613 vector<double>& delu,
01614 vector<double>& p,
01615 vector<double>& qtemp,
01616 vector<double>& Lambda,
01617 const vector<double>& Phi,
01618 const vector<double>& PhiDrv,
01619 double step,
01620 double *cooling,
01621 long index,
01622 size_t nd,
01623 bool *lgBoundFail)
01624 {
01625 long i,
01626 j,
01627 jlo,
01628 k=-1000;
01629 double bval_jk,
01630 cooling2,
01631 p2k,
01632 p_max,
01633 RelErrCool,
01634 RelErrPk,
01635 sum,
01636 sum2 = -DBL_MAX,
01637 trap1,
01638 trap12 = -DBL_MAX,
01639 trap2,
01640 uhi,
01641 ulo,
01642 umin,
01643 y;
01644
01645 DEBUG_ENTRY( "TryDoubleStep()" );
01646
01647
01648 ASSERT( index >= 0 && index < NQGRID-2 && nd < gv.bin.size() && step > 0. );
01649
01650 ulo = rfield.anu[0];
01651
01652
01653 uhi = rfield.anu[gv.bin[nd]->qnflux-1];
01654
01655
01656 p_max = 0.;
01657 for( i=0; i <= index; i++ )
01658 p_max = MAX2(p_max,p[i]);
01659 jlo = 0;
01660 while( p[jlo] < PROB_CUTOFF_LO*p_max )
01661 jlo++;
01662
01663 for( i=1; i <= 2; i++ )
01664 {
01665 bool lgErr;
01666 long ipLo = 0;
01667 long ipHi = gv.bin[nd]->qnflux-1;
01668 k = index + i;
01669 delu[k] = step/2.;
01670 u1[k] = u1[k-1] + delu[k];
01671 qtemp[k] = inv_ufunct(u1[k],nd,lgBoundFail);
01672 splint_safe(gv.dsttmp,gv.bin[nd]->dstems,gv.bin[nd]->dstslp2,NDEMS,log(qtemp[k]),&y,&lgErr);
01673 *lgBoundFail = *lgBoundFail || lgErr;
01674 Lambda[k] = exp(y)*gv.bin[nd]->cnv_H_pGR/EN1RYD;
01675
01676 sum = sum2 = 0.;
01677 trap1 = trap2 = trap12 = 0.;
01678
01679
01680 for( j=jlo; j < k; j++ )
01681 {
01682 umin = u1[k] - u1[j];
01683
01684 if( umin >= uhi )
01685 {
01686
01687
01688
01689
01690 continue;
01691 }
01692 else if( umin > ulo )
01693 {
01694
01695
01696
01697
01698
01699
01700 long ipStep = 1;
01701
01702 while( rfield.anu[ipLo] > (realnum)umin )
01703 {
01704 ipHi = ipLo;
01705 ipLo -= ipStep;
01706 if( ipLo <= 0 )
01707 {
01708 ipLo = 0;
01709 break;
01710 }
01711 ipStep *= 2;
01712 }
01713
01714 while( ipHi-ipLo > 1 )
01715 {
01716 long ipMd = (ipLo+ipHi)/2;
01717 if( rfield.anu[ipMd] > (realnum)umin )
01718 ipHi = ipMd;
01719 else
01720 ipLo = ipMd;
01721 }
01722 bval_jk = Phi[ipLo] + (umin - rfield.anu[ipLo])*PhiDrv[ipLo];
01723 }
01724 else
01725 {
01726 bval_jk = Phi[0];
01727 }
01728
01729
01730 trap12 = trap1;
01731 sum2 = sum;
01732
01733
01734
01735
01736
01737
01738
01739
01740 trap2 = p[j]*bval_jk;
01741
01742 sum += (trap1 + trap2)*delu[j];
01743 trap1 = trap2;
01744 }
01745
01746
01747
01748
01749
01750
01751 p[k] = (sum + trap1*delu[k])/(2.*Lambda[k] - Phi[0]*delu[k]);
01752
01753
01754 if( p[k] <= 0. )
01755 return 3.*QHEAT_TOL;
01756 }
01757
01758
01759 p2k = (sum2 + trap12*step)/(2.*Lambda[k] - Phi[0]*step);
01760
01761
01762 if( p2k <= 0. )
01763 return 3.*QHEAT_TOL;
01764
01765 RelErrPk = fabs(p2k-p[k])/p[k];
01766
01767
01768
01769 *cooling = log_integral(u1[k-2],p[k-2]*Lambda[k-2],u1[k-1],p[k-1]*Lambda[k-1]);
01770 *cooling += log_integral(u1[k-1],p[k-1]*Lambda[k-1],u1[k],p[k]*Lambda[k]);
01771
01772
01773 cooling2 = log_integral(u1[k-2],p[k-2]*Lambda[k-2],u1[k],p2k*Lambda[k]);
01774
01775
01776 RelErrCool = ( index > 0 ) ? fabs(cooling2-(*cooling))/(*cooling) : 0.;
01777
01778
01779
01780 return MAX2(RelErrPk,RelErrCool)/3.;
01781 }
01782
01783
01784
01785 STATIC double log_integral(double xx1,
01786 double yy1,
01787 double xx2,
01788 double yy2)
01789 {
01790 double eps,
01791 result,
01792 xx;
01793
01794 DEBUG_ENTRY( "log_integral()" );
01795
01796
01797 ASSERT( xx1 > 0. && yy1 > 0. && xx2 > 0. && yy2 > 0. );
01798
01799 xx = log(xx2/xx1);
01800 eps = log((xx2/xx1)*(yy2/yy1));
01801 if( fabs(eps) < 1.e-4 )
01802 {
01803 result = xx1*yy1*xx*((eps/6. + 0.5)*eps + 1.);
01804 }
01805 else
01806 {
01807 result = (xx2*yy2 - xx1*yy1)*xx/eps;
01808 }
01809 return result;
01810 }
01811
01812
01813
01814 STATIC void ScanProbDistr(const vector<double>& u1,
01815 const vector<double>& dPdlnT,
01816 long nbin,
01817 double maxVal,
01818 long nmax,
01819 double qtmin1,
01820 long *nstart,
01821 long *nstart2,
01822 long *nend,
01823 long *nWideFail,
01824 QH_Code *ErrorCode)
01825 {
01826 bool lgSetLo,
01827 lgSetHi;
01828 long i;
01829 double deriv_max,
01830 minVal;
01831 const double MIN_FAC_LO = 1.e4;
01832 const double MIN_FAC_HI = 1.e4;
01833
01834 DEBUG_ENTRY( "ScanProbDistr()" );
01835
01836
01837 ASSERT( nbin > 0 && nmax >= 0 && nmax < nbin && maxVal > 0. );
01838
01839
01840
01841
01842
01843
01844
01845 minVal = maxVal;
01846 *nstart = nmax;
01847 for( i=nmax; i >= 0; --i )
01848 {
01849 if( dPdlnT[i] < minVal )
01850 {
01851 *nstart = i;
01852 minVal = dPdlnT[i];
01853 }
01854 }
01855 deriv_max = 0.;
01856 *nstart2 = nmax;
01857 for( i=nmax; i > *nstart; --i )
01858 {
01859 double deriv = log(dPdlnT[i]/dPdlnT[i-1])/log(u1[i]/u1[i-1]);
01860 if( deriv > deriv_max )
01861 {
01862 *nstart2 = i-1;
01863 deriv_max = deriv;
01864 }
01865 }
01866 *nend = nbin-1;
01867
01868
01869 lgSetLo = ( nmax >= *nend || maxVal/dPdlnT[*nend] < MIN_FAC_HI );
01870
01871
01872
01873 if( lgSetLo )
01874
01875 lgSetHi = ( nmax <= *nstart || maxVal/dPdlnT[*nstart] < MIN_FAC_LO );
01876 else
01877 lgSetHi = ( nmax <= *nstart2 || maxVal/dPdlnT[*nstart2] < MIN_FAC_LO );
01878
01879 if( lgSetLo && lgSetHi )
01880 {
01881 ++(*nWideFail);
01882
01883 if( *nWideFail >= WIDE_FAIL_MAX )
01884 *ErrorCode = MAX2(*ErrorCode,QH_WIDE_FAIL);
01885 }
01886
01887 if( lgSetLo )
01888 *ErrorCode = MAX2(*ErrorCode,QH_LOW_FAIL);
01889
01890
01891
01892
01893
01894 if( lgSetHi && qtmin1 > 1.2*GRAIN_TMIN )
01895 *ErrorCode = MAX2(*ErrorCode,QH_HIGH_FAIL);
01896
01897
01898 if( *ErrorCode < QH_NO_REBIN && (*nend - *nstart) < NQMIN )
01899 *ErrorCode = MAX2(*ErrorCode,QH_NBIN_FAIL);
01900
01901 if( trace.lgTrace && trace.lgDustBug )
01902 {
01903 fprintf( ioQQQ, " ScanProbDistr nstart %ld nstart2 %ld nend %ld nmax %ld maxVal %.3e",
01904 *nstart,*nstart2,*nend,nmax,maxVal );
01905 fprintf( ioQQQ, " dPdlnT[nstart] %.3e dPdlnT[nstart2] %.3e dPdlnT[nend] %.3e code %d\n",
01906 dPdlnT[*nstart],dPdlnT[*nstart2],dPdlnT[*nend],*ErrorCode );
01907 }
01908
01909 if( *ErrorCode >= QH_NO_REBIN )
01910 {
01911 *nstart = -1;
01912 *nstart2 = -1;
01913 *nend = -2;
01914 }
01915 return;
01916 }
01917
01918
01919
01920 STATIC long RebinQHeatResults(size_t nd,
01921 long nstart,
01922 long nend,
01923 vector<double>& p,
01924 vector<double>& qtemp,
01925 vector<double>& qprob,
01926 vector<double>& dPdlnT,
01927 vector<double>& u1,
01928 vector<double>& delu,
01929 vector<double>& Lambda,
01930 QH_Code *ErrorCode)
01931 {
01932 long i,
01933 newnbin;
01934 double fac,
01935 help,
01936 mul_fac,
01937 PP1,
01938 PP2,
01939 RadCooling,
01940 T1,
01941 T2,
01942 Ucheck,
01943 uu1,
01944 uu2;
01945
01946 DEBUG_ENTRY( "RebinQHeatResults()" );
01947
01948
01949 ASSERT( nd < gv.bin.size() );
01950
01951 ASSERT( nstart >= 0 && nstart < nend && nend < NQGRID );
01952
01953
01954 for( i=nstart; i <= nend && dPdlnT[i] < PROB_CUTOFF_LO; i++ ) {}
01955
01956
01957 if( i >= NQGRID )
01958 {
01959 *ErrorCode = MAX2(*ErrorCode,QH_REBIN_FAIL);
01960 return 0;
01961 }
01962
01963
01964 vector<double> new_delu(NQGRID);
01965 vector<double> new_dPdlnT(NQGRID);
01966 vector<double> new_Lambda(NQGRID);
01967 vector<double> new_p(NQGRID);
01968 vector<double> new_qprob(NQGRID);
01969 vector<double> new_qtemp(NQGRID);
01970 vector<double> new_u1(NQGRID);
01971
01972 newnbin = 0;
01973
01974 T1 = qtemp[i];
01975 PP1 = p[i];
01976 uu1 = u1[i];
01977
01978
01979 help = pow(qtemp[nend]/qtemp[i],1./(1.5*NQMIN));
01980 mul_fac = MIN2(QT_RATIO,help);
01981
01982 Ucheck = u1[i];
01983 RadCooling = 0.;
01984
01985 while( i < nend )
01986 {
01987 bool lgBoundErr;
01988 bool lgDone= false;
01989 double s0 = 0.;
01990 double s1 = 0.;
01991 double wid = 0.;
01992 double xx,y;
01993
01994 T2 = T1*mul_fac;
01995
01996 do
01997 {
01998 double p1,p2,L1,L2,frac,slope;
01999 if( qtemp[i] <= T1 && T1 <= qtemp[i+1] )
02000 {
02001
02002
02003 frac = log(T1/qtemp[i]);
02004 slope = log(p[i+1]/p[i])/log(qtemp[i+1]/qtemp[i]);
02005 p1 = p[i]*exp(frac*slope);
02006 slope = log(Lambda[i+1]/Lambda[i])/log(qtemp[i+1]/qtemp[i]);
02007 L1 = Lambda[i]*exp(frac*slope);
02008 }
02009 else
02010 {
02011
02012
02013 p1 = p[i];
02014 L1 = Lambda[i];
02015 }
02016 if( qtemp[i] <= T2 && T2 <= qtemp[i+1] )
02017 {
02018
02019 help = ufunct(T2,nd,&lgBoundErr);
02020 uu2 = MIN2(help,u1[i+1]);
02021 ASSERT( !lgBoundErr );
02022 frac = log(T2/qtemp[i]);
02023 slope = log(p[i+1]/p[i])/log(qtemp[i+1]/qtemp[i]);
02024 p2 = p[i]*exp(frac*slope);
02025 slope = log(Lambda[i+1]/Lambda[i])/log(qtemp[i+1]/qtemp[i]);
02026 L2 = Lambda[i]*exp(frac*slope);
02027 lgDone = true;
02028 }
02029 else
02030 {
02031 uu2 = u1[i+1];
02032 p2 = p[i+1];
02033 L2 = Lambda[i+1];
02034
02035
02036 if( MAX2(p2,PP1)/MIN2(p2,PP1) > sqrt(SAFETY) )
02037 {
02038 lgDone = true;
02039 T2 = qtemp[i+1];
02040 }
02041 ++i;
02042 }
02043 PP2 = p2;
02044 wid += uu2 - uu1;
02045
02046 ASSERT( wid >= 0. );
02047 s0 += log_integral(uu1,p1,uu2,p2);
02048 s1 += log_integral(uu1,p1*L1,uu2,p2*L2);
02049 uu1 = uu2;
02050
02051 } while( i < nend && ! lgDone );
02052
02053
02054
02055
02056 if( s0 <= 0. )
02057 {
02058 ASSERT( wid == 0. );
02059 break;
02060 }
02061
02062 new_qprob[newnbin] = s0;
02063 new_Lambda[newnbin] = s1/s0;
02064 xx = log(new_Lambda[newnbin]*EN1RYD*gv.bin[nd]->cnv_GR_pH);
02065 splint_safe(gv.bin[nd]->dstems,gv.dsttmp,gv.bin[nd]->dstslp,NDEMS,xx,&y,&lgBoundErr);
02066 ASSERT( !lgBoundErr );
02067 new_qtemp[newnbin] = exp(y);
02068 new_u1[newnbin] = ufunct(new_qtemp[newnbin],nd,&lgBoundErr);
02069 ASSERT( !lgBoundErr );
02070 new_delu[newnbin] = wid;
02071 new_p[newnbin] = new_qprob[newnbin]/new_delu[newnbin];
02072 new_dPdlnT[newnbin] = new_p[newnbin]*new_qtemp[newnbin]*uderiv(new_qtemp[newnbin],nd);
02073
02074 Ucheck += wid;
02075 RadCooling += new_qprob[newnbin]*new_Lambda[newnbin];
02076
02077 T1 = T2;
02078 PP1 = PP2;
02079 ++newnbin;
02080 }
02081
02082
02083 if( newnbin < NQMIN )
02084 {
02085 *ErrorCode = MAX2(*ErrorCode,QH_REBIN_FAIL);
02086 return 0;
02087 }
02088
02089 fac = RadCooling*EN1RYD*gv.bin[nd]->cnv_GR_pCM3/gv.bin[nd]->GrainHeat;
02090
02091 if( trace.lgTrace && trace.lgDustBug )
02092 {
02093 fprintf( ioQQQ, " RebinQHeatResults found tol1 %.4e tol2 %.4e\n",
02094 fabs(u1[nend]/Ucheck-1.), fabs(fac-1.) );
02095 }
02096
02097
02098
02099 ASSERT( fabs(u1[nend]/Ucheck-1.) < 10.*sqrt((double)(nend-nstart+newnbin))*DBL_EPSILON );
02100
02101 if( fabs(fac-1.) > CONSERV_TOL )
02102 *ErrorCode = MAX2(*ErrorCode,QH_CONV_FAIL);
02103
02104 for( i=0; i < newnbin; i++ )
02105 {
02106
02107 p[i] = new_p[i]/fac;
02108 qtemp[i] = new_qtemp[i];
02109 qprob[i] = new_qprob[i]/fac;
02110 dPdlnT[i] = new_dPdlnT[i]/fac;
02111 u1[i] = new_u1[i];
02112 delu[i] = new_delu[i];
02113 Lambda[i] = new_Lambda[i];
02114
02115
02116 ASSERT( qtemp[i] > 0. && qprob[i] > 0. );
02117
02118
02119 }
02120 return newnbin;
02121 }
02122
02123
02124
02125 STATIC void GetProbDistr_HighLimit(long nd,
02126 double TolFac,
02127 double Umax,
02128 double fwhm,
02129 vector<double>& qtemp,
02130 vector<double>& qprob,
02131 vector<double>& dPdlnT,
02132 double *tol,
02133 long *qnbin,
02134 double *new_tmin,
02135 QH_Code *ErrorCode)
02136 {
02137 bool lgBoundErr,
02138 lgErr;
02139 long i,
02140 nbin;
02141 double c1,
02142 c2,
02143 delu[NQGRID],
02144 fac,
02145 fac1,
02146 fac2,
02147 help1,
02148 help2,
02149 L1,
02150 L2,
02151 Lambda[NQGRID],
02152 mul_fac,
02153 p[NQGRID],
02154 p1,
02155 p2,
02156 RadCooling,
02157 sum,
02158 T1,
02159 T2,
02160 Tlo,
02161 Thi,
02162 Ulo,
02163 Uhi,
02164 uu1,
02165 uu2,
02166 xx,
02167 y;
02168
02169 DEBUG_ENTRY( "GetProbDistr_HighLimit()" );
02170
02171 if( trace.lgTrace && trace.lgDustBug )
02172 {
02173 fprintf( ioQQQ, " GetProbDistr_HighLimit called for grain %s\n", gv.bin[nd]->chDstLab );
02174 }
02175
02176 c1 = sqrt(4.*LN_TWO/PI)/fwhm*exp(-POW2(fwhm/Umax)/(16.*LN_TWO));
02177 c2 = 4.*LN_TWO*POW2(Umax/fwhm);
02178
02179 fac1 = fwhm/Umax*sqrt(-log(PROB_CUTOFF_LO)/(4.*LN_TWO));
02180
02181 help1 = Umax*exp(-fac1);
02182 help2 = exp(gv.bin[nd]->DustEnth[0]);
02183 Ulo = MAX2(help1,help2);
02184
02185 Tlo = inv_ufunct(Ulo,nd,&lgBoundErr);
02186
02187 fac2 = fwhm/Umax*sqrt(-log(PROB_CUTOFF_HI)/(4.*LN_TWO));
02188
02189 if( fac2 > log(DBL_MAX/10.) )
02190 {
02191 *ErrorCode = MAX2(*ErrorCode,QH_WIDE_FAIL);
02192 return;
02193 }
02194 Uhi = Umax*exp(fac2);
02195 Thi = inv_ufunct(Uhi,nd,&lgBoundErr);
02196
02197 nbin = 0;
02198
02199 T1 = Tlo;
02200 uu1 = ufunct(T1,nd,&lgErr);
02201 lgBoundErr = lgBoundErr || lgErr;
02202 help1 = log(uu1/Umax);
02203 p1 = c1*exp(-c2*POW2(help1));
02204 splint_safe(gv.dsttmp,gv.bin[nd]->dstems,gv.bin[nd]->dstslp2,NDEMS,log(T1),&y,&lgErr);
02205 lgBoundErr = lgBoundErr || lgErr;
02206 L1 = exp(y)*gv.bin[nd]->cnv_H_pGR/EN1RYD;
02207
02208
02209 if( uu1*p1*L1 < 1.e5*DBL_MIN )
02210 {
02211 *ErrorCode = MAX2(*ErrorCode,QH_WIDE_FAIL);
02212 return;
02213 }
02214
02215
02216 help1 = pow(Thi/Tlo,1./(1.2*NQMIN));
02217 mul_fac = MIN2(QT_RATIO,help1);
02218
02219 sum = 0.;
02220 RadCooling = 0.;
02221
02222 do
02223 {
02224 double s0,s1,wid;
02225
02226 T2 = T1*mul_fac;
02227 uu2 = ufunct(T2,nd,&lgErr);
02228 lgBoundErr = lgBoundErr || lgErr;
02229 help1 = log(uu2/Umax);
02230 p2 = c1*exp(-c2*POW2(help1));
02231 splint_safe(gv.dsttmp,gv.bin[nd]->dstems,gv.bin[nd]->dstslp2,NDEMS,log(T2),&y,&lgErr);
02232 lgBoundErr = lgBoundErr || lgErr;
02233 L2 = exp(y)*gv.bin[nd]->cnv_H_pGR/EN1RYD;
02234
02235 wid = uu2 - uu1;
02236 s0 = log_integral(uu1,p1,uu2,p2);
02237 s1 = log_integral(uu1,p1*L1,uu2,p2*L2);
02238
02239 qprob[nbin] = s0;
02240 Lambda[nbin] = s1/s0;
02241 xx = log(Lambda[nbin]*EN1RYD*gv.bin[nd]->cnv_GR_pH);
02242 splint_safe(gv.bin[nd]->dstems,gv.dsttmp,gv.bin[nd]->dstslp,NDEMS,xx,&y,&lgErr);
02243 lgBoundErr = lgBoundErr || lgErr;
02244 qtemp[nbin] = exp(y);
02245 delu[nbin] = wid;
02246 p[nbin] = qprob[nbin]/delu[nbin];
02247 dPdlnT[nbin] = p[nbin]*qtemp[nbin]*uderiv(qtemp[nbin],nd);
02248
02249 sum += qprob[nbin];
02250 RadCooling += qprob[nbin]*Lambda[nbin];
02251
02252 T1 = T2;
02253 uu1 = uu2;
02254 p1 = p2;
02255 L1 = L2;
02256
02257 ++nbin;
02258
02259 } while( T2 < Thi && nbin < NQGRID );
02260
02261 fac = RadCooling*EN1RYD*gv.bin[nd]->cnv_GR_pCM3/gv.bin[nd]->GrainHeat;
02262
02263 for( i=0; i < nbin; ++i )
02264 {
02265 qprob[i] /= fac;
02266 dPdlnT[i] /= fac;
02267 }
02268
02269 *tol = sum/fac;
02270 *qnbin = nbin;
02271 *new_tmin = qtemp[0];
02272 *ErrorCode = MAX2(*ErrorCode,QH_ANALYTIC);
02273
02274
02275 if( TolFac > STRICT )
02276 *ErrorCode = MAX2(*ErrorCode,QH_ANALYTIC_RELAX);
02277
02278 if( lgBoundErr )
02279 *ErrorCode = MAX2(*ErrorCode,QH_THIGH_FAIL);
02280
02281 if( fabs(sum/fac-1.) > PROB_TOL )
02282 *ErrorCode = MAX2(*ErrorCode,QH_CONV_FAIL);
02283
02284 if( dPdlnT[0] > SAFETY*PROB_CUTOFF_LO || dPdlnT[nbin-1] > SAFETY*PROB_CUTOFF_HI )
02285 *ErrorCode = MAX2(*ErrorCode,QH_BOUND_FAIL);
02286
02287 if( trace.lgTrace && trace.lgDustBug )
02288 {
02289 fprintf( ioQQQ, " GetProbDistr_HighLimit found tol1 %.4e tol2 %.4e\n",
02290 fabs(sum-1.), fabs(sum/fac-1.) );
02291 fprintf( ioQQQ, " zone %ld %s nbin %ld total prob %.4e new_tmin %.4e\n",
02292 nzone,gv.bin[nd]->chDstLab,nbin,sum/fac,*new_tmin );
02293 }
02294 return;
02295 }
02296
02297
02298
02299 STATIC double uderiv(double temp,
02300 size_t nd)
02301 {
02302 enth_type ecase;
02303 long int i,
02304 j;
02305 double N_C,
02306 N_H;
02307 double deriv = 0.,
02308 hok[3] = {1275., 1670., 4359.},
02309 numer,
02310 dnumer,
02311 denom,
02312 ddenom,
02313 x;
02314
02315
02316 DEBUG_ENTRY( "uderiv()" );
02317
02318 if( temp <= 0. )
02319 {
02320 fprintf( ioQQQ, " uderiv called with non-positive temperature: %.6e\n" , temp );
02321 cdEXIT(EXIT_FAILURE);
02322 }
02323 ASSERT( nd < gv.bin.size() );
02324
02325 ecase = gv.which_enth[gv.bin[nd]->matType];
02326 switch( ecase )
02327 {
02328 case ENTH_CAR:
02329 numer = (4.15e-22/EN1RYD)*pow(temp,3.3);
02330 dnumer = (3.3*4.15e-22/EN1RYD)*pow(temp,2.3);
02331 denom = 1. + 6.51e-03*temp + 1.5e-06*temp*temp + 8.3e-07*pow(temp,2.3);
02332 ddenom = 6.51e-03 + 2.*1.5e-06*temp + 2.3*8.3e-07*pow(temp,1.3);
02333
02334
02335 deriv = (dnumer*denom - numer*ddenom)/POW2(denom);
02336 break;
02337 case ENTH_CAR2:
02338
02339
02340 deriv = (DebyeDeriv(temp/863.,2) + 2.*DebyeDeriv(temp/2504.,2))*BOLTZMANN/EN1RYD;
02341 break;
02342 case ENTH_SIL:
02343
02344
02345
02346 for( j = 0; j < 4; j++ )
02347 {
02348 if( temp > tlim[j] && temp <= tlim[j+1] )
02349 {
02350 deriv = cval[j]*pow(temp,ppower[j]);
02351 break;
02352 }
02353 }
02354 break;
02355 case ENTH_SIL2:
02356
02357
02358 deriv = (2.*DebyeDeriv(temp/500.,2) + DebyeDeriv(temp/1500.,3))*BOLTZMANN/EN1RYD;
02359 break;
02360 case ENTH_PAH:
02361
02362
02363
02364 x = log10(MIN2(temp,2000.));
02365 deriv = pow(10.,-21.26+3.1688*x-0.401894*POW2(x))/EN1RYD;
02366 break;
02367 case ENTH_PAH2:
02368
02369
02370
02371
02372 N_C = no_atoms(nd);
02373 if( N_C <= 25. )
02374 {
02375 N_H = 0.5*N_C;
02376 }
02377 else if( N_C <= 100. )
02378 {
02379 N_H = 2.5*sqrt(N_C);
02380 }
02381 else
02382 {
02383 N_H = 0.25*N_C;
02384 }
02385 deriv = 0.;
02386 for( i=0; i < 3; i++ )
02387 {
02388 double help1 = hok[i]/temp;
02389 if( help1 < 300. )
02390 {
02391 double help2 = exp(help1);
02392 double help3 = ( help1 < 1.e-7 ) ? help1*(1.+0.5*help1) : help2-1.;
02393 deriv += N_H/(N_C-2.)*POW2(help1)*help2/POW2(help3)*BOLTZMANN/EN1RYD;
02394 }
02395 }
02396 deriv += (DebyeDeriv(temp/863.,2) + 2.*DebyeDeriv(temp/2504.,2))*BOLTZMANN/EN1RYD;
02397 break;
02398 default:
02399 fprintf( ioQQQ, " uderiv called with unknown type for enthalpy function: %d\n", ecase );
02400 cdEXIT(EXIT_FAILURE);
02401 }
02402
02403
02404
02405 deriv *= MAX2(no_atoms(nd)-2.,1.);
02406
02407 if( deriv <= 0. )
02408 {
02409 fprintf( ioQQQ, " uderiv finds non-positive derivative: %.6e, what's up?\n" , deriv );
02410 cdEXIT(EXIT_FAILURE);
02411 }
02412 return( deriv );
02413 }
02414
02415
02416
02417 STATIC double ufunct(double temp,
02418 size_t nd,
02419 bool *lgBoundErr)
02420 {
02421 double enthalpy,
02422 y;
02423
02424 DEBUG_ENTRY( "ufunct()" );
02425
02426 if( temp <= 0. )
02427 {
02428 fprintf( ioQQQ, " ufunct called with non-positive temperature: %.6e\n" , temp );
02429 cdEXIT(EXIT_FAILURE);
02430 }
02431 ASSERT( nd < gv.bin.size() );
02432
02433
02434 splint_safe(gv.dsttmp,gv.bin[nd]->DustEnth,gv.bin[nd]->EnthSlp,NDEMS,log(temp),&y,lgBoundErr);
02435 enthalpy = exp(y);
02436
02437 ASSERT( enthalpy > 0. );
02438 return( enthalpy );
02439 }
02440
02441
02442
02443 STATIC double inv_ufunct(double enthalpy,
02444 size_t nd,
02445 bool *lgBoundErr)
02446 {
02447 double temp,
02448 y;
02449
02450 DEBUG_ENTRY( "inv_ufunct()" );
02451
02452 if( enthalpy <= 0. )
02453 {
02454 fprintf( ioQQQ, " inv_ufunct called with non-positive enthalpy: %.6e\n" , enthalpy );
02455 cdEXIT(EXIT_FAILURE);
02456 }
02457 ASSERT( nd < gv.bin.size() );
02458
02459
02460 splint_safe(gv.bin[nd]->DustEnth,gv.dsttmp,gv.bin[nd]->EnthSlp2,NDEMS,log(enthalpy),&y,lgBoundErr);
02461 temp = exp(y);
02462
02463 ASSERT( temp > 0. );
02464 return( temp );
02465 }
02466
02467
02468
02469 void InitEnthalpy(void)
02470 {
02471 double C_V1,
02472 C_V2,
02473 tdust1,
02474 tdust2;
02475
02476 DEBUG_ENTRY( "InitEnthalpy()" );
02477
02478 for( size_t nd=0; nd < gv.bin.size(); nd++ )
02479 {
02480 tdust2 = GRAIN_TMIN;
02481 C_V2 = uderiv(tdust2,nd);
02482
02483 gv.bin[nd]->DustEnth[0] = C_V2*tdust2/4.;
02484 tdust1 = tdust2;
02485 C_V1 = C_V2;
02486
02487 for( long i=1; i < NDEMS; i++ )
02488 {
02489 double tmid,Cmid;
02490 tdust2 = exp(gv.dsttmp[i]);
02491 C_V2 = uderiv(tdust2,nd);
02492 tmid = sqrt(tdust1*tdust2);
02493
02494 for( long j=1; j < 4; j++ )
02495 {
02496 if( tdust1 < tlim[j] && tlim[j] < tdust2 )
02497 {
02498 tmid = tlim[j];
02499 break;
02500 }
02501 }
02502 Cmid = uderiv(tmid,nd);
02503 gv.bin[nd]->DustEnth[i] = gv.bin[nd]->DustEnth[i-1] +
02504 log_integral(tdust1,C_V1,tmid,Cmid) +
02505 log_integral(tmid,Cmid,tdust2,C_V2);
02506 tdust1 = tdust2;
02507 C_V1 = C_V2;
02508 }
02509 }
02510
02511
02512 for( size_t nd=0; nd < gv.bin.size(); nd++ )
02513 {
02514 for( long i=0; i < NDEMS; i++ )
02515 {
02516 gv.bin[nd]->DustEnth[i] = log(gv.bin[nd]->DustEnth[i]);
02517 }
02518 }
02519
02520
02521 for( size_t nd=0; nd < gv.bin.size(); nd++ )
02522 {
02523
02524 spline(gv.dsttmp,gv.bin[nd]->DustEnth,NDEMS,2e31,2e31,gv.bin[nd]->EnthSlp);
02525 spline(gv.bin[nd]->DustEnth,gv.dsttmp,NDEMS,2e31,2e31,gv.bin[nd]->EnthSlp2);
02526 }
02527 return;
02528 }
02529
02530
02531
02532 STATIC double DebyeDeriv(double x,
02533 long n)
02534 {
02535 long i,
02536 nn;
02537 double res;
02538
02539 DEBUG_ENTRY( "DebyeDeriv()" );
02540
02541 ASSERT( x > 0. );
02542 ASSERT( n == 2 || n == 3 );
02543
02544 if( x < 0.001 )
02545 {
02546
02547 if( n == 2 )
02548 {
02549 res = 6.*1.202056903159594*POW2(x);
02550 }
02551 else if( n == 3 )
02552 {
02553 res = 24.*1.082323233711138*POW3(x);
02554 }
02555 else
02556
02557
02558 TotalInsanity();
02559 }
02560 else
02561 {
02562 nn = 4*MAX2(4,2*(long)(0.05/x));
02563 vector<double> xx(nn);
02564 vector<double> rr(nn);
02565 vector<double> aa(nn);
02566 vector<double> ww(nn);
02567 gauss_legendre(nn,xx,aa);
02568 gauss_init(nn,0.,1.,xx,aa,rr,ww);
02569
02570 res = 0.;
02571 for( i=0; i < nn; i++ )
02572 {
02573 double help1 = rr[i]/x;
02574 if( help1 < 300. )
02575 {
02576 double help2 = exp(help1);
02577 double help3 = ( help1 < 1.e-7 ) ? help1*(1.+0.5*help1) : help2-1.;
02578 res += ww[i]*powi(rr[i],n+1)*help2/POW2(help3);
02579 }
02580 }
02581 res /= POW2(x);
02582 }
02583 return (double)n*res;
02584 }