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