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