00001
00002
00003
00004 #include "cddefines.h"
00005 #include "physconst.h"
00006 #include "rfield.h"
00007 #include "hmi.h"
00008 #include "trace.h"
00009 #include "conv.h"
00010 #include "ionbal.h"
00011 #include "thermal.h"
00012 #include "phycon.h"
00013 #include "doppvel.h"
00014 #include "taulines.h"
00015 #include "mole.h"
00016 #include "heavy.h"
00017 #include "thirdparty.h"
00018 #include "dense.h"
00019 #include "ipoint.h"
00020 #include "elementnames.h"
00021 #include "grainvar.h"
00022 #include "grains.h"
00023
00024
00025
00026
00027
00028
00029
00030 #define FREE_CHECK(PTR) { ASSERT( PTR != NULL ); free( PTR ); PTR = NULL; }
00031 #define FREE_SAFE(PTR) { if( PTR != NULL ) free( PTR ); PTR = NULL; }
00032
00033 #define POT2CHRG(X) ((X)*EVRYD/ELEM_CHARGE*gv.bin[nd]->Capacity - 1.)
00034 #define CHRG2POT(X) (((X)+1.)*ELEM_CHARGE/EVRYD/gv.bin[nd]->Capacity)
00035 #define ONE_ELEC (ELEM_CHARGE/EVRYD/gv.bin[nd]->Capacity)
00036
00037
00038 #define EPSP(Z_0,Z) ( ((Z) <= (Z_0)) ? 0. : (((Z) >= (Z_0)+1) ? 1. : (double)((Z)-(Z_0))) )
00039
00040 static const int INCL_TUNNEL = true;
00041 static const int NO_TUNNEL = false;
00042
00043 static const int ALL_STAGES = true;
00044
00045
00046
00047
00048 static long int nCalledGrainDrive;
00049
00050
00051
00052
00053
00054 static const long NTOP = NDEMS/5;
00055
00056
00057
00058
00059
00060 static const double THERMCONST = PI4*ELECTRON_MASS*POW2(BOLTZMANN)/POW3(HPLANCK);
00061
00062
00063 static const double STICK_ELEC = 0.5;
00064 static const double STICK_ION = 1.0;
00065
00066
00067 static const double MEAN_PATH = 1.e-7;
00068
00069
00070
00071 static const double TOLER = CONSERV_TOL/10.;
00072 static const long BRACKET_MAX = 50L;
00073
00074
00075
00076 static const int NCHU = NCHS/3;
00077
00078
00079
00080
00081 static const long CT_LOOP_MAX = 25L;
00082
00083
00084 static const long T_LOOP_MAX = 50L;
00085
00086
00087 static double HEAT_TOLER = DBL_MAX;
00088 static double HEAT_TOLER_BIN = DBL_MAX;
00089 static double CHRG_TOLER = DBL_MAX;
00090
00091
00092
00093
00094
00095 static bool lgGvInitialized = false;
00096
00097
00098 static const double AC0 = 3.e-9;
00099 static const double AC1G = 4.e-8;
00100 static const double AC2G = 7.e-8;
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111 STATIC void InitEmissivities(void);
00112
00113 STATIC double PlanckIntegral(double,long,long);
00114
00115 STATIC void NewChargeData(long);
00116
00117 STATIC double GrnStdDpth(long);
00118
00119 STATIC void GrainChargeTemp(void);
00120
00121 STATIC void GrainCharge(long,double*);
00122
00123 STATIC double GrainElecRecomb1(long,long,double*,double*);
00124
00125 STATIC double GrainElecEmis1(long,long,double*,double*,
00126 double*,double*,double*);
00127
00128
00129 STATIC void GrainScreen(long,long,long,double*,double*);
00130
00131 STATIC double ThetaNu(double);
00132
00133 STATIC void UpdatePot(long,long,long,double[],double[]);
00134
00135 STATIC void GetFracPop(long,long,double[],double[],long*);
00136
00137 STATIC void UpdatePot1(long,long,long,long);
00138
00139 STATIC void UpdatePot2(long,long);
00140
00141 STATIC long HighestIonStage(void);
00142
00143 STATIC void UpdateRecomZ0(long,long,bool);
00144
00145 STATIC void GetPotValues(long,long,double*,double*,double*,
00146 double*,double*,bool);
00147
00148
00149
00150 STATIC void GrainIonColl(long,long,long,long,const double[],const double[],long*,
00151 float*,float*);
00152
00153 STATIC void GrainChrgTransferRates(long);
00154
00155 STATIC void GrainUpdateRadius1(void);
00156
00157 STATIC void GrainUpdateRadius2(bool);
00158
00159 STATIC void GrainTemperature(long,float*,double*,double*,
00160 double*);
00161
00162 STATIC void GrainCollHeating(long,float*,float*);
00163
00164 STATIC double GrnVryDpth(long);
00165
00166
00167
00168
00169
00170
00171
00172
00173 void GrainZero(void)
00174 {
00175 long int nelem, ion, ion_to;
00176
00177 DEBUG_ENTRY( "GrainZero()" );
00178
00179 gv.lgAnyDustVary = false;
00180 gv.TotalEden = 0.;
00181 gv.dHeatdT = 0.;
00182 gv.lgQHeatAll = false;
00183
00184
00185 gv.lgGrainElectrons = true;
00186 gv.lgQHeatOn = true;
00187 gv.lgDHetOn = true;
00188 gv.lgDColOn = true;
00189 gv.GrainMetal = 1.;
00190 gv.dstAbundThresholdNear = 1.e-6f;
00191 gv.dstAbundThresholdFar = 1.e-3f;
00192 gv.lgBakes = false;
00193 gv.nChrgRequested = NCHRG_DEFAULT;
00194 gv.ReadPtr = 0L;
00195
00196 gv.lgReevaluate = true;
00197
00198 gv.lgNegGrnDrg = false;
00199
00200
00201 nCalledGrainDrive = 0;
00202
00203
00204
00205 gv.lgBakesPAH_heat = false;
00206
00207
00208
00209 gv.lgGrainPhysicsOn = true;
00210
00211
00212
00213 gv.GrainHeatScaleFactor = 1.f;
00214
00215
00216
00217
00218 gv.which_enth[MAT_CAR] = ENTH_CAR;
00219 gv.which_zmin[MAT_CAR] = ZMIN_CAR;
00220 gv.which_pot[MAT_CAR] = POT_CAR;
00221 gv.which_ial[MAT_CAR] = IAL_CAR;
00222 gv.which_pe[MAT_CAR] = PE_CAR;
00223 gv.which_strg[MAT_CAR] = STRG_CAR;
00224 gv.which_H2distr[MAT_CAR] = H2_CAR;
00225
00226 gv.which_enth[MAT_SIL] = ENTH_SIL;
00227 gv.which_zmin[MAT_SIL] = ZMIN_SIL;
00228 gv.which_pot[MAT_SIL] = POT_SIL;
00229 gv.which_ial[MAT_SIL] = IAL_SIL;
00230 gv.which_pe[MAT_SIL] = PE_SIL;
00231 gv.which_strg[MAT_SIL] = STRG_SIL;
00232 gv.which_H2distr[MAT_SIL] = H2_SIL;
00233
00234 gv.which_enth[MAT_PAH] = ENTH_PAH;
00235 gv.which_zmin[MAT_PAH] = ZMIN_CAR;
00236 gv.which_pot[MAT_PAH] = POT_CAR;
00237 gv.which_ial[MAT_PAH] = IAL_CAR;
00238 gv.which_pe[MAT_PAH] = PE_CAR;
00239 gv.which_strg[MAT_PAH] = STRG_CAR;
00240 gv.which_H2distr[MAT_PAH] = H2_CAR;
00241
00242 gv.which_enth[MAT_CAR2] = ENTH_CAR2;
00243 gv.which_zmin[MAT_CAR2] = ZMIN_CAR;
00244 gv.which_pot[MAT_CAR2] = POT_CAR;
00245 gv.which_ial[MAT_CAR2] = IAL_CAR;
00246 gv.which_pe[MAT_CAR2] = PE_CAR;
00247 gv.which_strg[MAT_CAR2] = STRG_CAR;
00248 gv.which_H2distr[MAT_CAR2] = H2_CAR;
00249
00250 gv.which_enth[MAT_SIL2] = ENTH_SIL2;
00251 gv.which_zmin[MAT_SIL2] = ZMIN_SIL;
00252 gv.which_pot[MAT_SIL2] = POT_SIL;
00253 gv.which_ial[MAT_SIL2] = IAL_SIL;
00254 gv.which_pe[MAT_SIL2] = PE_SIL;
00255 gv.which_strg[MAT_SIL2] = STRG_SIL;
00256 gv.which_H2distr[MAT_SIL2] = H2_SIL;
00257
00258 gv.which_enth[MAT_PAH2] = ENTH_PAH2;
00259 gv.which_zmin[MAT_PAH2] = ZMIN_CAR;
00260 gv.which_pot[MAT_PAH2] = POT_CAR;
00261 gv.which_ial[MAT_PAH2] = IAL_CAR;
00262 gv.which_pe[MAT_PAH2] = PE_CAR;
00263 gv.which_strg[MAT_PAH2] = STRG_CAR;
00264 gv.which_H2distr[MAT_PAH2] = H2_CAR;
00265
00266 for( nelem=0; nelem < LIMELM; nelem++ )
00267 {
00268 for( ion=0; ion <= nelem+1; ion++ )
00269 {
00270 for( ion_to=0; ion_to <= nelem+1; ion_to++ )
00271 {
00272 gv.GrainChTrRate[nelem][ion][ion_to] = 0.f;
00273 }
00274 }
00275 }
00276
00277
00278
00279
00280 strcpy( gv.chPAH_abundance_fcn , "H0" );
00281
00282
00283
00284 ReturnGrainBins();
00285
00286 DEBUG_EXIT( "GrainZero()" );
00287 return;
00288 }
00289
00290
00291
00292
00293 void GrainStartIter(void)
00294 {
00295 long nd,
00296 i;
00297
00298 DEBUG_ENTRY( "GrainStartIter()" );
00299
00300 if( gv.lgDustOn && gv.lgGrainPhysicsOn )
00301 {
00302 for( i=0; i < rfield.nupper; i++ )
00303 {
00304
00305 gv.SilicateEmission[i] = 0.;
00306 gv.GraphiteEmission[i] = 0.;
00307 }
00308
00309 for( nd=0; nd < gv.nBin; nd++ )
00310 {
00311
00312
00313 gv.bin[nd]->dstpotsav = gv.bin[nd]->dstpot;
00314 gv.bin[nd]->qtmin = ( gv.bin[nd]->qtmin_zone1 > 0. ) ?
00315 gv.bin[nd]->qtmin_zone1 : DBL_MAX;
00316 gv.bin[nd]->avdust = 0.;
00317 gv.bin[nd]->avdpot = 0.;
00318 gv.bin[nd]->avdft = 0.;
00319 gv.bin[nd]->avDGRatio = 0.;
00320 gv.bin[nd]->TeGrainMax = -1.f;
00321 gv.bin[nd]->lgEverQHeat = false;
00322 gv.bin[nd]->QHeatFailures = 0L;
00323 gv.bin[nd]->lgQHTooWide = false;
00324 gv.bin[nd]->lgPAHsInIonizedRegion = false;
00325 gv.bin[nd]->nChrgOrg = gv.bin[nd]->nChrg;
00326 }
00327 }
00328
00329 DEBUG_EXIT( "GrainStartIter()" );
00330 return;
00331 }
00332
00333
00334
00335
00336 void GrainRestartIter(void)
00337 {
00338 long nd;
00339
00340 DEBUG_ENTRY( "GrainRestartIter()" );
00341
00342 if( gv.lgDustOn && gv.lgGrainPhysicsOn )
00343 {
00344 for( nd=0; nd < gv.nBin; nd++ )
00345 {
00346
00347
00348 gv.bin[nd]->dstpot = gv.bin[nd]->dstpotsav;
00349 gv.bin[nd]->nChrg = gv.bin[nd]->nChrgOrg;
00350 }
00351 }
00352
00353 DEBUG_EXIT( "GrainRestartIter()" );
00354 return;
00355 }
00356
00357
00358
00359 void SetNChrgStates(long nChrg)
00360 {
00361 DEBUG_ENTRY( "SetNChrgStates()" );
00362
00363 ASSERT( nChrg >= 2 && nChrg <= NCHU );
00364 gv.nChrgRequested = nChrg;
00365
00366 DEBUG_EXIT( "SetNChrgStates()" );
00367 return;
00368 }
00369
00370
00371 long NewGrainBin(void)
00372 {
00373 long nd, nz;
00374
00375 DEBUG_ENTRY( "NewGrainBin()" );
00376
00377 ASSERT( lgGvInitialized );
00378
00379 if( gv.nBin >= NDUST )
00380 {
00381 fprintf( ioQQQ, " The code has run out of grain bins; increase NDUST and recompile.\n" );
00382 puts( "[Stop in NewGrainBin]" );
00383 cdEXIT(EXIT_FAILURE);
00384 }
00385 nd = gv.nBin;
00386
00387 ASSERT( gv.bin[nd] == NULL );
00388 gv.bin[nd] = (GrainBin *)MALLOC(sizeof(GrainBin));
00389
00390
00391 gv.bin[nd]->dstab1 = NULL;
00392 gv.bin[nd]->pure_sc1 = NULL;
00393 gv.bin[nd]->asym = NULL;
00394 gv.bin[nd]->inv_att_len = NULL;
00395 gv.bin[nd]->y1 = NULL;
00396 for( nz=0; nz < NCHS; nz++ )
00397 gv.bin[nd]->chrg[nz] = NULL;
00398
00399 gv.lgDustOn = true;
00400 gv.bin[nd]->lgQHeat = false;
00401 gv.bin[nd]->qnflux = LONG_MAX;
00402 gv.bin[nd]->lgDustFunc = false;
00403 gv.bin[nd]->DustDftVel = FLT_MAX;
00404 gv.bin[nd]->TeGrainMax = FLT_MAX;
00405
00406 gv.bin[nd]->nChrg = gv.nChrgRequested;
00407
00408
00409
00410 gv.bin[nd]->tedust = 1.f;
00411 gv.bin[nd]->GrainHeat = DBL_MAX/10.;
00412 gv.bin[nd]->GrainGasCool = DBL_MAX/10.;
00413
00414
00415 gv.bin[nd]->EnergyCheck = 0.;
00416 gv.bin[nd]->dstAbund = -FLT_MAX;
00417 gv.bin[nd]->dstfactor = 1.f;
00418 gv.bin[nd]->cnv_H_pGR = -DBL_MAX;
00419 gv.bin[nd]->cnv_H_pCM3 = -DBL_MAX;
00420 gv.bin[nd]->cnv_CM3_pGR = -DBL_MAX;
00421 gv.bin[nd]->cnv_CM3_pH = -DBL_MAX;
00422 gv.bin[nd]->cnv_GR_pH = -DBL_MAX;
00423 gv.bin[nd]->cnv_GR_pCM3 = -DBL_MAX;
00424
00425 gv.bin[nd]->rate_h2_form_grains_used = 0.;
00426
00427 gv.nBin++;
00428
00429 DEBUG_EXIT( "NewGrainBin()" );
00430 return nd;
00431 }
00432
00433
00434 void ReturnGrainBins(void)
00435 {
00436 long nd, nz;
00437
00438 DEBUG_ENTRY( "ReturnGrainBins()" );
00439
00440 if( lgGvInitialized )
00441 {
00442
00443 for( nd=0; nd < gv.nBin; nd++ )
00444 {
00445 ASSERT( gv.bin[nd] != NULL );
00446
00447 FREE_SAFE( gv.bin[nd]->dstab1 );
00448 FREE_SAFE( gv.bin[nd]->pure_sc1 );
00449 FREE_SAFE( gv.bin[nd]->asym );
00450 FREE_SAFE( gv.bin[nd]->y1 );
00451 FREE_SAFE( gv.bin[nd]->inv_att_len );
00452
00453 for( nz=0; nz < NCHS; nz++ )
00454 {
00455 if( gv.bin[nd]->chrg[nz] != NULL )
00456 {
00457 FREE_SAFE( gv.bin[nd]->chrg[nz]->yhat );
00458 FREE_SAFE( gv.bin[nd]->chrg[nz]->cs_pdt );
00459 FREE_SAFE( gv.bin[nd]->chrg[nz]->fac1 );
00460 FREE_SAFE( gv.bin[nd]->chrg[nz]->fac2 );
00461
00462 FREE_CHECK( gv.bin[nd]->chrg[nz] );
00463 }
00464 }
00465
00466 FREE_CHECK( gv.bin[nd] );
00467 }
00468
00469 FREE_SAFE( gv.anumin );
00470 FREE_SAFE( gv.anumax );
00471 FREE_SAFE( gv.dstab );
00472 FREE_SAFE( gv.dstsc );
00473 FREE_SAFE( gv.GrainEmission );
00474 FREE_SAFE( gv.GraphiteEmission );
00475 FREE_SAFE( gv.SilicateEmission );
00476 }
00477 else
00478 {
00479
00480
00481 for( nd=0; nd < NDUST; nd++ )
00482 {
00483 gv.bin[nd] = NULL;
00484 }
00485
00486 gv.anumin = NULL;
00487 gv.anumax = NULL;
00488 gv.dstab = NULL;
00489 gv.dstsc = NULL;
00490 gv.GrainEmission = NULL;
00491 gv.GraphiteEmission = NULL;
00492 gv.SilicateEmission = NULL;
00493
00494 lgGvInitialized = true;
00495 }
00496
00497 gv.lgDustOn = false;
00498 gv.nBin = 0;
00499
00500 DEBUG_EXIT( "ReturnGrainBins()" );
00501 return;
00502 }
00503
00504
00505
00506
00507
00508 void GrainsInit(void)
00509 {
00510 long int i,
00511 nelem,
00512 nd,
00513 nz;
00514
00515 DEBUG_ENTRY( "GrainsInit()" );
00516
00517 if( trace.lgTrace && trace.lgDustBug )
00518 {
00519 fprintf( ioQQQ, " GrainsInit called.\n" );
00520 }
00521
00522
00523
00524 ASSERT( gv.anumin == NULL );
00525 gv.anumin = (float*)MALLOC((size_t)((unsigned)rfield.nupper*sizeof(float)));
00526 ASSERT( gv.anumax == NULL );
00527 gv.anumax = (float*)MALLOC((size_t)((unsigned)rfield.nupper*sizeof(float)));
00528 ASSERT( gv.dstab == NULL );
00529 gv.dstab = (double*)MALLOC((size_t)((unsigned)rfield.nupper*sizeof(double)));
00530 ASSERT( gv.dstsc == NULL );
00531 gv.dstsc = (double*)MALLOC((size_t)((unsigned)rfield.nupper*sizeof(double)));
00532 ASSERT( gv.GrainEmission == NULL );
00533 gv.GrainEmission = (float*)MALLOC((size_t)((unsigned)rfield.nupper*sizeof(float)));
00534 ASSERT( gv.GraphiteEmission == NULL );
00535 gv.GraphiteEmission = (float*)MALLOC((size_t)((unsigned)rfield.nupper*sizeof(float)));
00536 ASSERT( gv.SilicateEmission == NULL );
00537 gv.SilicateEmission = (float*)MALLOC((size_t)((unsigned)rfield.nupper*sizeof(float)));
00538
00539
00540 ASSERT( gv.nBin >= 0 && gv.nBin < NDUST );
00541 for( nd=gv.nBin; nd < NDUST; nd++ )
00542 {
00543 ASSERT( gv.bin[nd] == NULL );
00544 }
00545
00546
00547 for( nelem=0; nelem < LIMELM; nelem++ )
00548 {
00549 gv.elmSumAbund[nelem] = 0.f;
00550 }
00551
00552 for( i=0; i < rfield.nupper; i++ )
00553 {
00554 gv.dstab[i] = 0.;
00555 gv.dstsc[i] = 0.;
00556
00557 gv.GrainEmission[i] = 0.;
00558 gv.SilicateEmission[i] = 0.;
00559 gv.GraphiteEmission[i] = 0.;
00560 }
00561
00562 if( !gv.lgDustOn )
00563 {
00564
00565 gv.GrainHeatInc = 0.;
00566 gv.GrainHeatDif = 0.;
00567 gv.GrainHeatLya = 0.;
00568 gv.GrainHeatCollSum = 0.;
00569 gv.GrainHeatSum = 0.;
00570 gv.GasCoolColl = 0.;
00571 thermal.heating[0][13] = 0.;
00572 thermal.heating[0][14] = 0.;
00573 thermal.heating[0][25] = 0.;
00574
00575 if( trace.lgTrace && trace.lgDustBug )
00576 {
00577 fprintf( ioQQQ, " GrainsInit exits.\n" );
00578 }
00579
00580 DEBUG_EXIT( "GrainsInit()" );
00581 return;
00582 }
00583
00584 HEAT_TOLER = conv.HeatCoolRelErrorAllowed / 3.;
00585 HEAT_TOLER_BIN = HEAT_TOLER / sqrt((double)gv.nBin);
00586 CHRG_TOLER = conv.EdenErrorAllowed / 3.;
00587
00588
00589 for( nd=0; nd < gv.nBin; nd++ )
00590 {
00591 double help,atoms,p_rad,ThresInf,dum[4];
00592 long low1,low2,low3,lowm;
00593
00594
00595 ASSERT( gv.bin[nd] != NULL );
00596
00597
00598 if( gv.lgQHeatAll )
00599 {
00600 gv.bin[nd]->lgQHeat = true;
00601 }
00602
00603
00604 if( !gv.lgQHeatOn )
00605 {
00606 gv.bin[nd]->lgQHeat = false;
00607 }
00608
00609
00610 if( thermal.ConstGrainTemp > 0. )
00611 {
00612 gv.bin[nd]->lgQHeat = false;
00613 }
00614
00615 #ifndef IGNORE_QUANTUM_HEATING
00616 gv.bin[nd]->lgQHTooWide = false;
00617 gv.bin[nd]->qtmin = DBL_MAX;
00618 #endif
00619
00620 if( gv.bin[nd]->lgDustFunc || gv.bin[nd]->matType == MAT_PAH || gv.bin[nd]->matType == MAT_PAH2 )
00621 gv.lgAnyDustVary = true;
00622
00623
00624
00625 gv.bin[nd]->dstAbund = -FLT_MAX;
00626
00627 gv.bin[nd]->qtmin_zone1 = -1.;
00628
00629 if( gv.bin[nd]->DustWorkFcn < rfield.anu[0] || gv.bin[nd]->DustWorkFcn > rfield.anu[rfield.nupper] )
00630 {
00631 fprintf( ioQQQ, " Grain work function for %s has insane value: %.4e\n",
00632 gv.bin[nd]->chDstLab,gv.bin[nd]->DustWorkFcn );
00633 puts( "[Stop in GrainsInit]" );
00634 cdEXIT(EXIT_FAILURE);
00635 }
00636
00637
00638 ASSERT( gv.bin[nd]->y1 == NULL );
00639 gv.bin[nd]->y1 = (float*)MALLOC((size_t)((unsigned)rfield.nupper*sizeof(float)));
00640 for( nz=0; nz < NCHS; nz++ )
00641 {
00642 ASSERT( gv.bin[nd]->chrg[nz] == NULL );
00643 gv.bin[nd]->chrg[nz] = (ChargeBin *)MALLOC(sizeof(ChargeBin));
00644 gv.bin[nd]->chrg[nz]->yhat = NULL;
00645 gv.bin[nd]->chrg[nz]->cs_pdt = NULL;
00646 gv.bin[nd]->chrg[nz]->fac1 = NULL;
00647 gv.bin[nd]->chrg[nz]->fac2 = NULL;
00648 gv.bin[nd]->chrg[nz]->nfill = 0;
00649
00650 gv.bin[nd]->chrg[nz]->DustZ = LONG_MIN;
00651 gv.bin[nd]->chrg[nz]->FracPop = -DBL_MAX;
00652 gv.bin[nd]->chrg[nz]->tedust = 1.f;
00653 }
00654
00655 for( i=0; i < rfield.nupper; i++ )
00656 {
00657 double alpha,beta,af,bf;
00658
00659 beta = gv.bin[nd]->AvRadius*gv.bin[nd]->inv_att_len[i];
00660 if( beta > 1.e-4 )
00661 {
00662 bf = POW2(beta) - 2.*beta + 2. - 2.*exp(-beta);
00663 }
00664 else
00665 {
00666 bf = POW3(beta)/3.;
00667 }
00668 alpha = beta + gv.bin[nd]->AvRadius/MEAN_PATH;
00669 if( alpha > 1.e-4 )
00670 {
00671 af = POW2(alpha) - 2.*alpha + 2. - 2.*exp(-alpha);
00672 }
00673 else
00674 {
00675 af = POW3(alpha)/3.;
00676 }
00677
00678
00679
00680 gv.bin[nd]->y1[i] = (float)(POW2(beta/alpha)*af/bf);
00681 }
00682
00683 ASSERT( gv.bin[nd]->nChrg >= 2 && gv.bin[nd]->nChrg <= NCHU );
00684
00685
00686
00687 if( gv.lgBakes )
00688 {
00689
00690
00691 help = ceil(POT2CHRG(gv.bin[nd]->BandGap-gv.bin[nd]->DustWorkFcn+ONE_ELEC/2.));
00692 low1 = nint(help);
00693 }
00694 else
00695 {
00696 zmin_type zcase = gv.which_zmin[gv.bin[nd]->matType];
00697
00698 switch( zcase )
00699 {
00700 case ZMIN_CAR:
00701 help = gv.bin[nd]->AvRadius*1.e7;
00702 help = ceil(-(1.2*POW2(help)+3.9*help+0.2)/1.44);
00703 low1 = nint(help);
00704
00705
00706 break;
00707 case ZMIN_SIL:
00708 help = gv.bin[nd]->AvRadius*1.e7;
00709 help = ceil(-(0.7*POW2(help)+2.5*help+0.8)/1.44);
00710 low1 = nint(help);
00711
00712
00713 break;
00714 case ZMIN_BAKES:
00715
00716 help = ceil(POT2CHRG(gv.bin[nd]->BandGap-gv.bin[nd]->DustWorkFcn+ONE_ELEC/2.));
00717 low1 = nint(help);
00718 break;
00719 default:
00720 fprintf( ioQQQ, " GrainsInit detected unknown Zmin type: %d\n" , zcase );
00721 puts( "[Stop in GrainsInit]" );
00722 cdEXIT(1);
00723 }
00724 }
00725
00726
00727 ASSERT( help > (double)(LONG_MIN+1) );
00728
00729
00730
00731
00732
00733
00734 low2 = low1;
00735 GetPotValues(nd,low2,&ThresInf,&dum[0],&dum[1],&dum[2],&dum[3],INCL_TUNNEL);
00736 if( ThresInf < 0. )
00737 {
00738 low3 = 0;
00739
00740
00741 while( low3-low2 > 1 )
00742 {
00743 lowm = (low2+low3)/2;
00744 GetPotValues(nd,lowm,&ThresInf,&dum[0],&dum[1],&dum[2],&dum[3],INCL_TUNNEL);
00745 if( ThresInf < 0. )
00746 low2 = lowm;
00747 else
00748 low3 = lowm;
00749 }
00750 low2 = low3;
00751 }
00752
00753
00754
00755
00756 gv.bin[nd]->LowestZg = MAX2(low1,low2);
00757 gv.bin[nd]->LowestPot = CHRG2POT(gv.bin[nd]->LowestZg);
00758
00759
00760
00761
00762
00763
00764 gv.bin[nd]->StickElecPos = STICK_ELEC*(1. - exp(-gv.bin[nd]->AvRadius/MEAN_PATH));
00765 atoms = gv.bin[nd]->AvVol*gv.bin[nd]->dustp[0]/ATOMIC_MASS_UNIT/gv.bin[nd]->atomWeight;
00766 p_rad = 1./(1.+exp(20.-atoms));
00767 gv.bin[nd]->StickElecNeg = gv.bin[nd]->StickElecPos*p_rad;
00768
00769
00770
00771
00772 gv.bin[nd]->dstAbund = (float)(gv.bin[nd]->dstfactor*gv.GrainMetal*GrnVryDpth(nd));
00773 ASSERT( gv.bin[nd]->dstAbund > 0.f );
00774
00775 gv.bin[nd]->cnv_H_pCM3 = dense.gas_phase[ipHYDROGEN]*gv.bin[nd]->dstAbund;
00776 gv.bin[nd]->cnv_CM3_pH = 1./gv.bin[nd]->cnv_H_pCM3;
00777
00778 gv.bin[nd]->cnv_CM3_pGR = gv.bin[nd]->cnv_H_pGR/gv.bin[nd]->cnv_H_pCM3;
00779 gv.bin[nd]->cnv_GR_pCM3 = 1./gv.bin[nd]->cnv_CM3_pGR;
00780 }
00781
00782
00783
00784
00785
00786
00787 for( nelem=0; nelem < LIMELM; nelem++ )
00788 {
00789 gv.elmSumAbund[nelem] = 0.f;
00790 }
00791
00792 for( nd=0; nd < gv.nBin; nd++ )
00793 {
00794 for( nelem=0; nelem < LIMELM; nelem++ )
00795 {
00796 gv.elmSumAbund[nelem] += gv.bin[nd]->elmAbund[nelem]*(float)gv.bin[nd]->cnv_H_pCM3;
00797 }
00798 }
00799
00800 gv.nfill = -1;
00801 gv.nzone = -1;
00802 gv.lgAnyNegCharge = false;
00803 gv.GrnRecomTe = -1.f;
00804
00805
00806
00807 for( i=0; i < rfield.nupper; i++ )
00808 {
00809
00810
00811 gv.dstab[i] = -DBL_MAX;
00812 gv.dstsc[i] = -DBL_MAX;
00813 }
00814
00815 for( i=0; i < rfield.nupper; i++ )
00816 {
00817 gv.anumin[i] = rfield.anu[i] - 0.5f*rfield.widflx[i];
00818 gv.anumax[i] = rfield.anu[i] + 0.5f*rfield.widflx[i];
00819 }
00820
00821 InitEmissivities();
00822 InitEnthalpy();
00823
00824 if( trace.lgDustBug && trace.lgTrace )
00825 {
00826 fprintf( ioQQQ, " There are %ld grain types turned on.\n", gv.nBin );
00827
00828 fprintf( ioQQQ, " grain depletion factors, dstfactor*GrainMetal=" );
00829 for( nd=0; nd < gv.nBin; nd++ )
00830 {
00831 fprintf( ioQQQ, "%10.2e", gv.bin[nd]->dstfactor*gv.GrainMetal );
00832 }
00833 fprintf( ioQQQ, "\n" );
00834
00835 fprintf( ioQQQ, " nChrg =" );
00836 for( nd=0; nd < gv.nBin; nd++ )
00837 {
00838 fprintf( ioQQQ, " %ld", gv.bin[nd]->nChrg );
00839 }
00840 fprintf( ioQQQ, "\n" );
00841
00842 fprintf( ioQQQ, " lowest charge (e) =" );
00843 for( nd=0; nd < gv.nBin; nd++ )
00844 {
00845 fprintf( ioQQQ, "%10.2e", POT2CHRG(gv.bin[nd]->LowestPot) );
00846 }
00847 fprintf( ioQQQ, "\n" );
00848
00849 fprintf( ioQQQ, " lgDustFunc flag for user requested custom depth dependence:" );
00850 for( nd=0; nd < gv.nBin; nd++ )
00851 {
00852 fprintf( ioQQQ, "%2c", TorF(gv.bin[nd]->lgDustFunc) );
00853 }
00854 fprintf( ioQQQ, "\n" );
00855
00856 fprintf( ioQQQ, " Quantum heating flag:" );
00857 for( nd=0; nd < gv.nBin; nd++ )
00858 {
00859 fprintf( ioQQQ, "%2c", TorF(gv.bin[nd]->lgQHeat) );
00860 }
00861 fprintf( ioQQQ, "\n" );
00862
00863
00864 fprintf( ioQQQ, " NU(Ryd), Abs cross sec per proton\n" );
00865
00866 fprintf( ioQQQ, " Ryd " );
00867 for( nd=0; nd < gv.nBin; nd++ )
00868 {
00869 fprintf( ioQQQ, " %-12.12s", gv.bin[nd]->chDstLab );
00870 }
00871 fprintf( ioQQQ, "\n" );
00872
00873 for( i=0; i < rfield.nupper; i += 40 )
00874 {
00875 fprintf( ioQQQ, "%10.2e", rfield.anu[i] );
00876 for( nd=0; nd < gv.nBin; nd++ )
00877 {
00878 fprintf( ioQQQ, " %10.2e ", gv.bin[nd]->dstab1[i] );
00879 }
00880 fprintf( ioQQQ, "\n" );
00881 }
00882
00883 fprintf( ioQQQ, " NU(Ryd), Sct cross sec per proton\n" );
00884
00885 fprintf( ioQQQ, " Ryd " );
00886 for( nd=0; nd < gv.nBin; nd++ )
00887 {
00888 fprintf( ioQQQ, " %-12.12s", gv.bin[nd]->chDstLab );
00889 }
00890 fprintf( ioQQQ, "\n" );
00891
00892 for( i=0; i < rfield.nupper; i += 40 )
00893 {
00894 fprintf( ioQQQ, "%10.2e", rfield.anu[i] );
00895 for( nd=0; nd < gv.nBin; nd++ )
00896 {
00897 fprintf( ioQQQ, " %10.2e ", gv.bin[nd]->pure_sc1[i] );
00898 }
00899 fprintf( ioQQQ, "\n" );
00900 }
00901
00902 fprintf( ioQQQ, " NU(Ryd), Q abs\n" );
00903
00904 fprintf( ioQQQ, " Ryd " );
00905 for( nd=0; nd < gv.nBin; nd++ )
00906 {
00907 fprintf( ioQQQ, " %-12.12s", gv.bin[nd]->chDstLab );
00908 }
00909 fprintf( ioQQQ, "\n" );
00910
00911 for( i=0; i < rfield.nupper; i += 40 )
00912 {
00913 fprintf( ioQQQ, "%10.2e", rfield.anu[i] );
00914 for( nd=0; nd < gv.nBin; nd++ )
00915 {
00916 fprintf( ioQQQ, " %10.2e ", gv.bin[nd]->dstab1[i]*4./gv.bin[nd]->IntArea );
00917 }
00918 fprintf( ioQQQ, "\n" );
00919 }
00920
00921 fprintf( ioQQQ, " NU(Ryd), Q sct\n" );
00922
00923 fprintf( ioQQQ, " Ryd " );
00924 for( nd=0; nd < gv.nBin; nd++ )
00925 {
00926 fprintf( ioQQQ, " %-12.12s", gv.bin[nd]->chDstLab );
00927 }
00928 fprintf( ioQQQ, "\n" );
00929
00930 for( i=0; i < rfield.nupper; i += 40 )
00931 {
00932 fprintf( ioQQQ, "%10.2e", rfield.anu[i] );
00933 for( nd=0; nd < gv.nBin; nd++ )
00934 {
00935 fprintf( ioQQQ, " %10.2e ", gv.bin[nd]->pure_sc1[i]*4./gv.bin[nd]->IntArea );
00936 }
00937 fprintf( ioQQQ, "\n" );
00938 }
00939
00940 fprintf( ioQQQ, " NU(Ryd), asymmetry factor\n" );
00941
00942 fprintf( ioQQQ, " Ryd " );
00943 for( nd=0; nd < gv.nBin; nd++ )
00944 {
00945 fprintf( ioQQQ, " %-12.12s", gv.bin[nd]->chDstLab );
00946 }
00947 fprintf( ioQQQ, "\n" );
00948
00949 for( i=0; i < rfield.nupper; i += 40 )
00950 {
00951 fprintf( ioQQQ, "%10.2e", rfield.anu[i] );
00952 for( nd=0; nd < gv.nBin; nd++ )
00953 {
00954 fprintf( ioQQQ, " %10.2e ", gv.bin[nd]->asym[i] );
00955 }
00956 fprintf( ioQQQ, "\n" );
00957 }
00958
00959 fprintf( ioQQQ, " NU(Ryd), yield enhancement\n" );
00960
00961 fprintf( ioQQQ, " Ryd " );
00962 for( nd=0; nd < gv.nBin; nd++ )
00963 {
00964 fprintf( ioQQQ, " %-12.12s", gv.bin[nd]->chDstLab );
00965 }
00966 fprintf( ioQQQ, "\n" );
00967
00968 for( i=0; i < rfield.nupper; i += 40 )
00969 {
00970 fprintf( ioQQQ, "%10.2e", rfield.anu[i] );
00971 for( nd=0; nd < gv.nBin; nd++ )
00972 {
00973 fprintf( ioQQQ, " %10.2e ", gv.bin[nd]->y1[i] );
00974 }
00975 fprintf( ioQQQ, "\n" );
00976 }
00977
00978 fprintf( ioQQQ, " GrainsInit exits.\n" );
00979 }
00980
00981 DEBUG_EXIT( "GrainsInit()" );
00982 return;
00983 }
00984
00985 STATIC void InitEmissivities(void)
00986 {
00987 double fac,
00988 fac2,
00989 mul,
00990 tdust;
00991 long int i,
00992 nd;
00993
00994 DEBUG_ENTRY( "InitEmissivities()" );
00995
00996 if( trace.lgTrace && trace.lgDustBug )
00997 {
00998 fprintf( ioQQQ, " InitEmissivities starts\n" );
00999 fprintf( ioQQQ, " ND Tdust Emis BB Check 4pi*a^2*<Q>\n" );
01000 }
01001
01002 ASSERT( NTOP >= 2 && NDEMS > 2*NTOP );
01003 fac = exp(log(GRAIN_TMID/GRAIN_TMIN)/(double)(NDEMS-NTOP));
01004 tdust = GRAIN_TMIN;
01005 for( i=0; i < NDEMS-NTOP; i++ )
01006 {
01007 gv.dsttmp[i] = log(tdust);
01008 for( nd=0; nd < gv.nBin; nd++ )
01009 {
01010 gv.bin[nd]->dstems[i] = log(PlanckIntegral(tdust,nd,i));
01011 }
01012 tdust *= fac;
01013 }
01014
01015
01016 fac2 = exp(log(GRAIN_TMAX/GRAIN_TMID/powi(fac,NTOP-1))/(double)((NTOP-1)*NTOP/2));
01017 for( i=NDEMS-NTOP; i < NDEMS; i++ )
01018 {
01019 gv.dsttmp[i] = log(tdust);
01020 for( nd=0; nd < gv.nBin; nd++ )
01021 {
01022 gv.bin[nd]->dstems[i] = log(PlanckIntegral(tdust,nd,i));
01023 }
01024 fac *= fac2;
01025 tdust *= fac;
01026 }
01027
01028
01029 mul = 1.;
01030 ASSERT( fabs(gv.dsttmp[0] - log(GRAIN_TMIN)) < 10.*mul*DBL_EPSILON );
01031 mul = sqrt((double)(NDEMS-NTOP));
01032 ASSERT( fabs(gv.dsttmp[NDEMS-NTOP] - log(GRAIN_TMID)) < 10.*mul*DBL_EPSILON );
01033 mul = (double)NTOP + sqrt((double)NDEMS);
01034 ASSERT( fabs(gv.dsttmp[NDEMS-1] - log(GRAIN_TMAX)) < 10.*mul*DBL_EPSILON );
01035
01036
01037 for( nd=0; nd < gv.nBin; nd++ )
01038 {
01039
01040 spline(gv.bin[nd]->dstems,gv.dsttmp,NDEMS,2e31,2e31,gv.bin[nd]->dstslp);
01041 spline(gv.dsttmp,gv.bin[nd]->dstems,NDEMS,2e31,2e31,gv.bin[nd]->dstslp2);
01042 }
01043
01044 # if 0
01045
01046 nd = nint(fudge(0));
01047 ASSERT( nd >= 0 && nd < gv.nBin );
01048 for( i=0; i < 2000; i++ )
01049 {
01050 double x,y,z;
01051 z = pow(10.,-40. + (double)i/50.);
01052 splint(gv.bin[nd]->dstems,gv.dsttmp,gv.bin[nd]->dstslp,NDEMS,log(z),&y);
01053 if( exp(y) > GRAIN_TMIN && exp(y) < GRAIN_TMAX )
01054 {
01055 x = PlanckIntegral(exp(y),nd,3);
01056 printf(" input %.6e temp %.6e output %.6e rel. diff. %.6e\n",z,exp(y),x,(x-z)/z);
01057 }
01058 }
01059 cdEXIT(EXIT_SUCCESS);
01060 # endif
01061
01062 DEBUG_EXIT( "InitEmissivities()" );
01063 return;
01064 }
01065
01066
01067
01068 STATIC double PlanckIntegral(double tdust,
01069 long int nd,
01070 long int ip)
01071 {
01072 long int i;
01073
01074 double arg,
01075 ExpM1,
01076 integral1 = 0.,
01077 integral2 = 0.,
01078 Planck1,
01079 Planck2,
01080 TDustRyg,
01081 x;
01082
01083 DEBUG_ENTRY( "PlanckIntegral()" );
01084
01085
01086
01087
01088
01089
01090
01091
01092
01093
01094
01095 TDustRyg = TE1RYD/tdust;
01096
01097 x = 0.999*log(DBL_MAX);
01098
01099 for( i=0; i < rfield.nupper; i++ )
01100 {
01101
01102 arg = TDustRyg*rfield.anu[i];
01103
01104
01105 if( arg < 1.e-5 )
01106 {
01107
01108 ExpM1 = arg*(1. + arg/2.);
01109 }
01110 else
01111 {
01112
01113 ExpM1 = exp(MIN2(x,arg)) - 1.;
01114 }
01115
01116 Planck1 = PI4*2.*HPLANCK/POW2(SPEEDLIGHT)*POW2(FR1RYD)*POW2(FR1RYD)*
01117 rfield.anu3[i]/ExpM1*rfield.widflx[i];
01118 Planck2 = Planck1*gv.bin[nd]->dstab1[i];
01119
01120
01121 if( i == 0 )
01122 {
01123 integral1 = Planck1/rfield.widflx[0]*rfield.anu[0]/3.;
01124 integral2 = Planck2/rfield.widflx[0]*rfield.anu[0]/5.;
01125 }
01126
01127 if( Planck1/integral1 < DBL_EPSILON && Planck2/integral2 < DBL_EPSILON )
01128 break;
01129
01130 integral1 += Planck1;
01131 integral2 += Planck2;
01132 }
01133
01134
01135 if( trace.lgDustBug && trace.lgTrace && ip%10 == 0 )
01136 {
01137 fprintf( ioQQQ, " %4ld %11.4e %11.4e %11.4e %11.4e\n", nd, tdust,
01138 integral2, integral1/4./5.67051e-5/powi(tdust,4), integral2*4./integral1 );
01139 }
01140
01141 ASSERT( integral2 > 0. );
01142
01143 DEBUG_EXIT( "PlanckIntegral()" );
01144 return integral2;
01145 }
01146
01147
01148
01149 STATIC void NewChargeData(long nd)
01150 {
01151 long nz;
01152
01153 DEBUG_ENTRY( "NewChargeData()" );
01154
01155 for( nz=0; nz < NCHS; nz++ )
01156 {
01157 gv.bin[nd]->chrg[nz]->RSum1 = -DBL_MAX;
01158 gv.bin[nd]->chrg[nz]->RSum2 = -DBL_MAX;
01159 gv.bin[nd]->chrg[nz]->ESum1a = -DBL_MAX;
01160 gv.bin[nd]->chrg[nz]->ESum1b = -DBL_MAX;
01161 gv.bin[nd]->chrg[nz]->ESum1c = -DBL_MAX;
01162 gv.bin[nd]->chrg[nz]->ESum2 = -DBL_MAX;
01163
01165 gv.bin[nd]->chrg[nz]->ThermRate = -DBL_MAX;
01166 gv.bin[nd]->chrg[nz]->HeatingRate2 = -DBL_MAX;
01167 gv.bin[nd]->chrg[nz]->GrainHeat = -DBL_MAX;
01168
01169 gv.bin[nd]->chrg[nz]->hots1 = -DBL_MAX;
01170 gv.bin[nd]->chrg[nz]->bolflux1 = -DBL_MAX;
01171 gv.bin[nd]->chrg[nz]->pe1 = -DBL_MAX;
01172 }
01173
01174 if( fabs(phycon.te-gv.GrnRecomTe) > 10.f*FLT_EPSILON*gv.GrnRecomTe )
01175 {
01176 for( nz=0; nz < NCHS; nz++ )
01177 {
01178 memset( gv.bin[nd]->chrg[nz]->eta, 0, (LIMELM+2)*sizeof(double) );
01179 memset( gv.bin[nd]->chrg[nz]->xi, 0, (LIMELM+2)*sizeof(double) );
01180 }
01181 }
01182
01183 if( nzone != gv.nzone )
01184 {
01185 for( nz=0; nz < NCHS; nz++ )
01186 {
01187 gv.bin[nd]->chrg[nz]->hcon1 = -DBL_MAX;
01188 }
01189 }
01190
01191 DEBUG_EXIT( "NewChargeData()" );
01192 return;
01193 }
01194
01195
01196
01197
01198 STATIC double GrnStdDpth(long int nd)
01199 {
01200 double GrnStdDpth_v;
01201
01202 DEBUG_ENTRY( "GrnStdDpth()" );
01203
01204
01205
01206
01207
01208
01209 if( gv.bin[nd]->matType == MAT_PAH || gv.bin[nd]->matType == MAT_PAH2 )
01210 {
01211
01212 if( strcmp( gv.chPAH_abundance_fcn , "H0" )==0 )
01213 {
01214
01215
01216
01217
01218 GrnStdDpth_v = dense.xIonDense[ipHYDROGEN][0]/dense.gas_phase[ipHYDROGEN];
01219 }
01220 else if( strcmp( gv.chPAH_abundance_fcn , "CON" )==0 )
01221 {
01222
01223 GrnStdDpth_v = 1.;
01224 }
01225 else
01226 TotalInsanity();
01227 }
01228 else
01229 {
01230
01231 GrnStdDpth_v = 1.;
01232 }
01233
01234 ASSERT( GrnStdDpth_v > 0. && GrnStdDpth_v <= 1.000001 );
01235
01236
01237 DEBUG_EXIT( "GrnStdDpth()" );
01238 return GrnStdDpth_v;
01239 }
01240
01241
01242
01243 void GrainDrive(void)
01244 {
01245 DEBUG_ENTRY( "GrainDrive()" );
01246
01247
01248 if( gv.lgDustOn && gv.lgGrainPhysicsOn )
01249 {
01250 static double tesave=-1.f;
01251 static long int nzonesave=-1;
01252
01253
01254
01255
01256
01257
01258
01259 if( gv.lgReevaluate || conv.lgSearch || nzonesave != nzone ||
01260
01261
01262 phycon.te !=tesave ||
01263
01264
01265 fabs(gv.TotalEden)/dense.eden > conv.EdenErrorAllowed/5. ||
01266
01267
01268 (fabs( gv.GasCoolColl ) + fabs( thermal.heating[0][13] ))/SDIV(thermal.ctot)>0.1 )
01269
01270
01271
01272 {
01273 nzonesave = nzone;
01274 tesave = phycon.te;
01275
01276 if( trace.nTrConvg >= 5 )
01277 {
01278 fprintf( ioQQQ, " grain5 calling GrainChargeTemp and GrainDrift\n");
01279 }
01280
01281
01282 GrainChargeTemp();
01283
01284
01285 }
01286 }
01287 else if( gv.lgDustOn && !gv.lgGrainPhysicsOn )
01288 {
01289
01290
01291 if( nCalledGrainDrive == 0 )
01292 {
01293 long nelem, ion, ion_to, nd;
01294
01295
01296
01297 gv.GasHeatPhotoEl = 0.;
01298 for( nd=0; nd < gv.nBin; nd++ )
01299 {
01300 long nz;
01301
01302
01303 gv.bin[nd]->lgPAHsInIonizedRegion = false;
01304
01305 for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
01306 {
01307 gv.bin[nd]->chrg[nz]->DustZ = nz;
01308 gv.bin[nd]->chrg[nz]->FracPop = ( nz == 0 ) ? 1. : 0.;
01309 gv.bin[nd]->chrg[nz]->nfill = 0;
01310 gv.bin[nd]->chrg[nz]->tedust = 100.f;
01311 }
01312
01313 gv.bin[nd]->AveDustZ = 0.;
01314 gv.bin[nd]->dstpot = CHRG2POT(0.);
01315
01316 gv.bin[nd]->tedust = 100.f;
01317 gv.bin[nd]->TeGrainMax = 100.;
01318
01319
01320 gv.bin[nd]->BolFlux = 0.;
01321 gv.bin[nd]->GrainCoolTherm = 0.;
01322 gv.bin[nd]->GasHeatPhotoEl = 0.;
01323 gv.bin[nd]->GrainHeat = 0.;
01324 gv.bin[nd]->GrainHeatColl = 0.;
01325 gv.bin[nd]->ChemEn = 0.;
01326 gv.bin[nd]->ChemEnH2 = 0.;
01327 gv.bin[nd]->thermionic = 0.;
01328
01329 gv.bin[nd]->lgUseQHeat = false;
01330 gv.bin[nd]->lgEverQHeat = false;
01331 gv.bin[nd]->QHeatFailures = 0;
01332
01333 gv.bin[nd]->DustDftVel = 0.;
01334
01335 gv.bin[nd]->avdust = gv.bin[nd]->tedust;
01336 gv.bin[nd]->avdft = 0.f;
01337 gv.bin[nd]->avdpot = (float)(gv.bin[nd]->dstpot*EVRYD);
01338 gv.bin[nd]->avDGRatio = -1.f;
01339
01340
01341
01342 if( 0 && gv.lgBakesPAH_heat )
01343 {
01344
01345
01346
01347
01348
01349 gv.bin[nd]->GasHeatPhotoEl = 1.e-24*hmi.UV_Cont_rel2_Habing_TH85_depth*
01350 dense.gas_phase[ipHYDROGEN]*(4.87e-2/(1.0+4e-3*pow((hmi.UV_Cont_rel2_Habing_TH85_depth*
01351 sqrt(phycon.te)/dense.eden),0.73)) + 3.65e-2*pow(phycon.te/1.e4,0.7)/
01352 (1.+2.e-4*(hmi.UV_Cont_rel2_Habing_TH85_depth*sqrt(phycon.te)/dense.eden)))/gv.nBin *
01353 gv.GrainHeatScaleFactor;
01354 gv.GasHeatPhotoEl += gv.bin[nd]->GasHeatPhotoEl;
01355 }
01356 }
01357
01358 gv.lgAnyNegCharge = false;
01359
01360 gv.TotalEden = 0.;
01361 gv.GrnElecDonateMax = 0.f;
01362 gv.GrnElecHoldMax = 0.f;
01363
01364 for( nelem=0; nelem < LIMELM; nelem++ )
01365 {
01366 for( ion=0; ion <= nelem+1; ion++ )
01367 {
01368 for( ion_to=0; ion_to <= nelem+1; ion_to++ )
01369 {
01370 gv.GrainChTrRate[nelem][ion][ion_to] = 0.f;
01371 }
01372 }
01373 }
01374
01375
01376 gv.GrainHeatInc = 0.;
01377 gv.GrainHeatDif = 0.;
01378 gv.GrainHeatLya = 0.;
01379 gv.GrainHeatCollSum = 0.;
01380 gv.GrainHeatSum = 0.;
01381 gv.GrainHeatChem = 0.;
01382 gv.GasCoolColl = 0.;
01383 gv.TotalDustHeat = 0.f;
01384 gv.dphmax = 0.f;
01385 gv.dclmax = 0.f;
01386
01387 thermal.heating[0][13] = 0.;
01388 thermal.heating[0][14] = 0.;
01389 thermal.heating[0][25] = 0.;
01390 }
01391
01392 if( nCalledGrainDrive == 0 || gv.lgAnyDustVary )
01393 {
01394 GrainUpdateRadius1();
01395 GrainUpdateRadius2(false);
01396 }
01397 }
01398
01399 ++nCalledGrainDrive;
01400
01401 DEBUG_EXIT( "GrainDrive()" );
01402 return;
01403 }
01404
01405
01406 STATIC void GrainChargeTemp(void)
01407 {
01408 bool lgAnyNegCharge = false;
01409 long int i,
01410 ion,
01411 ion_to,
01412 nelem,
01413 nd,
01414 nz;
01415 float dccool = FLT_MAX;
01416 double delta,
01417 GasHeatNet,
01418 hcon = DBL_MAX,
01419 hla = DBL_MAX,
01420 hots = DBL_MAX,
01421 oldtemp,
01422 oldTotalEden,
01423 ratio,
01424 ThermRatio;
01425
01426 static long int oldZone = -1;
01427 static double oldTe = -DBL_MAX,
01428 oldHeat = -DBL_MAX;
01429
01430 DEBUG_ENTRY( "GrainChargeTemp()" );
01431
01432 if( trace.lgTrace && trace.lgDustBug )
01433 {
01434 fprintf( ioQQQ, "\n GrainChargeTemp called lgSearch%2c\n\n", TorF(conv.lgSearch) );
01435 }
01436
01437 oldTotalEden = gv.TotalEden;
01438
01439
01440 gv.GrainHeatInc = 0.;
01441 gv.GrainHeatDif = 0.;
01442 gv.GrainHeatLya = 0.;
01443 gv.GrainHeatCollSum = 0.;
01444 gv.GrainHeatSum = 0.;
01445 gv.GrainHeatChem = 0.;
01446
01447 gv.GasCoolColl = 0.;
01448 gv.GasHeatPhotoEl = 0.;
01449 gv.GasHeatTherm = 0.;
01450
01451 gv.TotalEden = 0.;
01452
01453 for( nelem=0; nelem < LIMELM; nelem++ )
01454 {
01455 for( ion=0; ion <= nelem+1; ion++ )
01456 {
01457 for( ion_to=0; ion_to <= nelem+1; ion_to++ )
01458 {
01459 gv.GrainChTrRate[nelem][ion][ion_to] = 0.f;
01460 }
01461 }
01462 }
01463
01464 gv.HighestIon = HighestIonStage();
01465
01466
01467 GrainUpdateRadius1();
01468
01469 for( nd=0; nd < gv.nBin; nd++ )
01470 {
01471 double one;
01472 double ChTdBracketLo = 0., ChTdBracketHi = -DBL_MAX;
01473 long relax = ( conv.lgSearch ) ? 3 : 1;
01474
01475
01476 if( gv.bin[nd]->matType == MAT_PAH || gv.bin[nd]->matType == MAT_PAH2 )
01477 {
01478 if( dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN] > 0.50 )
01479 {
01480 gv.bin[nd]->lgPAHsInIonizedRegion = true;
01481 }
01482 }
01483
01484
01485
01486
01487 NewChargeData(nd);
01488
01489 if( trace.lgTrace && trace.lgDustBug )
01490 {
01491 fprintf( ioQQQ, " >>GrainChargeTemp starting grain %s\n",
01492 gv.bin[nd]->chDstLab );
01493 }
01494
01495 delta = 2.*TOLER;
01496
01497 for( i=0; i < relax*CT_LOOP_MAX && delta > TOLER; ++i )
01498 {
01499 char which[20];
01500 long j;
01501 double TdBracketLo = 0., TdBracketHi = -DBL_MAX;
01502 double ThresEst = 0.;
01503 oldtemp = gv.bin[nd]->tedust;
01504
01505
01506
01507
01508
01509 GrainCharge(nd,&ThermRatio);
01510
01511 ASSERT( gv.bin[nd]->GrainHeat > 0. );
01512 ASSERT( gv.bin[nd]->tedust >= GRAIN_TMIN && gv.bin[nd]->tedust <= GRAIN_TMAX );
01513
01514
01515
01516
01517
01518 gv.bin[nd]->lgTdustConverged = false;
01519 for( j=0; j < relax*T_LOOP_MAX; ++j )
01520 {
01521 double oldTemp2 = gv.bin[nd]->tedust;
01522 double oldHeat2 = gv.bin[nd]->GrainHeat;
01523 double oldCool = gv.bin[nd]->GrainGasCool;
01524
01525
01526 GrainTemperature(nd,&dccool,&hcon,&hots,&hla);
01527
01528 gv.bin[nd]->GrainGasCool = dccool;
01529
01530 if( trace.lgTrace && trace.lgDustBug )
01531 {
01532 fprintf( ioQQQ, " >>loop %ld BracketLo %.6e BracketHi %.6e",
01533 j, TdBracketLo, TdBracketHi );
01534 }
01535
01536
01537
01538
01539
01540 if( fabs(gv.bin[nd]->GrainHeat-oldHeat2) < HEAT_TOLER*gv.bin[nd]->GrainHeat &&
01541 fabs(gv.bin[nd]->GrainGasCool-oldCool) < HEAT_TOLER_BIN*thermal.ctot )
01542 {
01543 gv.bin[nd]->lgTdustConverged = true;
01544 if( trace.lgTrace && trace.lgDustBug )
01545 fprintf( ioQQQ, " converged\n" );
01546 break;
01547 }
01548
01549
01550 if( gv.bin[nd]->tedust < oldTemp2 )
01551 TdBracketHi = oldTemp2;
01552 else
01553 TdBracketLo = oldTemp2;
01554
01555
01556
01557
01558
01559
01562
01563 if( TdBracketHi > TdBracketLo )
01564 {
01565
01566
01567 if( ( j >= 2 && TdBracketLo > 0. ) ||
01568 gv.bin[nd]->tedust <= TdBracketLo ||
01569 gv.bin[nd]->tedust >= TdBracketHi )
01570 {
01571 gv.bin[nd]->tedust = (float)(0.5*(TdBracketLo + TdBracketHi));
01572 if( trace.lgTrace && trace.lgDustBug )
01573 fprintf( ioQQQ, " bisection\n" );
01574 }
01575 else
01576 {
01577 if( trace.lgTrace && trace.lgDustBug )
01578 fprintf( ioQQQ, " iteration\n" );
01579 }
01580 }
01581 else
01582 {
01583 if( trace.lgTrace && trace.lgDustBug )
01584 fprintf( ioQQQ, " iteration\n" );
01585 }
01586
01587 ASSERT( gv.bin[nd]->tedust >= GRAIN_TMIN && gv.bin[nd]->tedust <= GRAIN_TMAX );
01588 }
01589
01590 if( gv.bin[nd]->lgTdustConverged )
01591 {
01592
01593 if( gv.bin[nd]->tedust < oldtemp )
01594 ChTdBracketHi = oldtemp;
01595 else
01596 ChTdBracketLo = oldtemp;
01597 }
01598 else
01599 {
01600 bool lgBoundErr;
01601 double y, x = log(gv.bin[nd]->tedust);
01602
01603 splint_safe(gv.dsttmp,gv.bin[nd]->dstems,gv.bin[nd]->dstslp2,NDEMS,x,&y,&lgBoundErr);
01604 gv.bin[nd]->GrainHeat = exp(y)*gv.bin[nd]->cnv_H_pCM3;
01605
01606 fprintf( ioQQQ," PROBLEM temperature of grain species %s (Tg=%.3eK) not converged\n",
01607 gv.bin[nd]->chDstLab , gv.bin[nd]->tedust );
01608 ConvFail("grai","");
01609 if( lgAbort )
01610 {
01611 DEBUG_EXIT( "GrainChargeTemp()" );
01612 return;
01613 }
01614 }
01615
01616 ASSERT( gv.bin[nd]->GrainHeat > 0. );
01617 ASSERT( gv.bin[nd]->tedust >= GRAIN_TMIN && gv.bin[nd]->tedust <= GRAIN_TMAX );
01618
01619
01620
01621
01622
01623
01625 ratio = gv.bin[nd]->tedust/oldtemp;
01626 for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
01627 {
01628 ThresEst += gv.bin[nd]->chrg[nz]->FracPop*gv.bin[nd]->chrg[nz]->ThresInf;
01629 }
01630 delta = ThresEst*TE1RYD/gv.bin[nd]->tedust*(ratio - 1.);
01632 delta = ( delta < 0.9*log(DBL_MAX) ) ?
01633 ThermRatio*fabs(POW2(ratio)*exp(delta)-1.) : DBL_MAX;
01634
01635
01636 if( delta > TOLER )
01637 {
01638 if( trace.lgTrace && trace.lgDustBug )
01639 strncpy( which, "iteration", sizeof(which) );
01640
01641
01642
01643
01644
01645
01648
01649 if( ChTdBracketHi > ChTdBracketLo )
01650 {
01651
01652
01653 if( ( i >= 2 && ChTdBracketLo > 0. ) ||
01654 gv.bin[nd]->tedust <= ChTdBracketLo ||
01655 gv.bin[nd]->tedust >= ChTdBracketHi )
01656 {
01657 gv.bin[nd]->tedust = (float)(0.5*(ChTdBracketLo + ChTdBracketHi));
01658 if( trace.lgTrace && trace.lgDustBug )
01659 strncpy( which, "bisection", sizeof(which) );
01660 }
01661 }
01662 }
01663
01664 if( trace.lgTrace && trace.lgDustBug )
01665 {
01666 fprintf( ioQQQ, " >>GrainChargeTemp finds delta=%.4e, ", delta );
01667 fprintf( ioQQQ, " old/new temp=%.5e %.5e, ", oldtemp, gv.bin[nd]->tedust );
01668 if( delta > TOLER )
01669 fprintf( ioQQQ, "doing another %s\n", which );
01670 else
01671 fprintf( ioQQQ, "converged\n" );
01672 }
01673 }
01674 if( delta > TOLER )
01675 {
01676 fprintf( ioQQQ, " PROBLEM charge/temperature not converged for %s zone %.2f\n",
01677 gv.bin[nd]->chDstLab , fnzone );
01678 ConvFail("grai","");
01679 }
01680
01681 for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
01682 {
01683 if( gv.bin[nd]->chrg[nz]->DustZ <= -1 )
01684 lgAnyNegCharge = true;
01685 }
01686
01687
01688
01689
01690 if( ionbal.lgGrainIonRecom )
01691 GrainChrgTransferRates(nd);
01692
01693
01694
01695
01696
01697
01698 gv.GrainHeatInc += hcon;
01699
01700 gv.GrainHeatDif += hots;
01701
01702
01703 gv.GrainHeatLya += hla;
01704
01705
01706 gv.GrainHeatCollSum += gv.bin[nd]->GrainHeatColl;
01707
01708
01709
01710
01711 gv.GrainHeatSum += gv.bin[nd]->GrainHeat;
01712
01713
01714 gv.GrainHeatChem += gv.bin[nd]->ChemEn + gv.bin[nd]->ChemEnH2;
01715
01716
01717
01718 gv.GasCoolColl += dccool;
01719 gv.GasHeatPhotoEl += gv.bin[nd]->GasHeatPhotoEl;
01720 gv.GasHeatTherm += gv.bin[nd]->thermionic;
01721
01722
01723
01724 one = 0.;
01725 for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
01726 {
01727 one += gv.bin[nd]->chrg[nz]->FracPop*(double)gv.bin[nd]->chrg[nz]->DustZ*
01728 gv.bin[nd]->cnv_GR_pCM3;
01729 }
01730
01731 gv.TotalEden += one;
01732 {
01733
01734 enum {DEBUG_LOC=false};
01735
01736 if( DEBUG_LOC )
01737 {
01738 fprintf(ioQQQ," DEBUG grn chr nz\t%.2f\teden\t%.3e\tnd\t%li",
01739 fnzone,
01740 dense.eden,
01741 nd);
01742 fprintf(ioQQQ,"\tne\t%.2e\tAveDustZ\t%.2e\t%.2e\t%.2e\t%.2e",
01743 one,
01744 gv.bin[nd]->AveDustZ,
01745 gv.bin[nd]->chrg[0]->FracPop,(double)gv.bin[nd]->chrg[0]->DustZ,
01746 gv.bin[nd]->cnv_GR_pCM3);
01747 fprintf(ioQQQ,"\n");
01748 }
01749 }
01750
01751 if( trace.lgTrace && trace.lgDustBug )
01752 {
01753 fprintf(ioQQQ," %s Pot %.5e Thermal %.5e GasCoolColl %.5e" ,
01754 gv.bin[nd]->chDstLab, gv.bin[nd]->dstpot, gv.bin[nd]->GrainHeat, dccool );
01755 fprintf(ioQQQ," GasPEHeat %.5e GasThermHeat %.5e ChemHeat %.5e\n\n" ,
01756 gv.bin[nd]->GasHeatPhotoEl, gv.bin[nd]->thermionic, gv.bin[nd]->ChemEn );
01757 }
01758 }
01759
01760
01761 GasHeatNet = gv.GasHeatPhotoEl + gv.GasHeatTherm - gv.GasCoolColl;
01762
01763 if( fabs(phycon.te-gv.GrnRecomTe) > 10.f*FLT_EPSILON*gv.GrnRecomTe )
01764 {
01765 oldZone = gv.nzone;
01766 oldTe = gv.GrnRecomTe;
01767 oldHeat = gv.GasHeatNet;
01768 }
01769
01770
01771 if( nzone == oldZone && fabs(phycon.te-oldTe) > 10.f*FLT_EPSILON*phycon.te )
01772 {
01773 gv.dHeatdT = (GasHeatNet-oldHeat)/(phycon.te-oldTe);
01774 }
01775
01776
01777 if( nzone != gv.nzone || fabs(phycon.te-gv.GrnRecomTe) > 10.f*FLT_EPSILON*gv.GrnRecomTe ||
01778 fabs(gv.GasHeatNet-GasHeatNet) > HEAT_TOLER*thermal.ctot ||
01779 fabs(gv.TotalEden-oldTotalEden) > CHRG_TOLER*dense.eden )
01780 {
01781
01782
01783
01784 conv.lgConvIoniz = false;
01785 if( fabs(gv.TotalEden-oldTotalEden) > CHRG_TOLER*dense.eden )
01786 {
01787 strcpy( conv.chConvIoniz, "grn eden chg" );
01788 conv.BadConvIoniz[0] = oldTotalEden;
01789 conv.BadConvIoniz[1] = gv.TotalEden;
01790 }
01791 else if( fabs(gv.GasHeatNet-GasHeatNet) > HEAT_TOLER*thermal.ctot )
01792 {
01793 strcpy( conv.chConvIoniz, "grn het chg" );
01794 conv.BadConvIoniz[0] = gv.GasHeatNet;
01795 conv.BadConvIoniz[1] = GasHeatNet;
01796 }
01797 else if( fabs(phycon.te-gv.GrnRecomTe) > 10.f*FLT_EPSILON*gv.GrnRecomTe )
01798 {
01799 strcpy( conv.chConvIoniz, "grn ter chg" );
01800 conv.BadConvIoniz[0] = gv.GrnRecomTe;
01801 conv.BadConvIoniz[1] = phycon.te;
01802 }
01803 else if( nzone != gv.nzone )
01804 {
01805 strcpy( conv.chConvIoniz, "grn zon chg" );
01806 conv.BadConvIoniz[0] = gv.nzone;
01807 conv.BadConvIoniz[1] = nzone;
01808 }
01809 else
01810 TotalInsanity();
01811 }
01812
01813
01814
01815
01816
01817 gv.nzone = nzone;
01818 gv.GrnRecomTe = (float)phycon.te;
01819 gv.GasHeatNet = GasHeatNet;
01820
01821
01822
01823
01824 GrainUpdateRadius2(lgAnyNegCharge);
01825
01826 #ifdef WD_TEST2
01827 printf("wd test: proton fraction %.5e Total DustZ %.6f heating/cooling rate %.5e %.5e\n",
01828 dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN],
01829 gv.bin[0]->AveDustZ,gv.GasHeatPhotoEl/dense.gas_phase[ipHYDROGEN]/fudge(0),
01830 gv.GasCoolColl/dense.gas_phase[ipHYDROGEN]/fudge(0));
01831 #endif
01832
01833 if( trace.lgTrace )
01834 {
01835
01836 enum {DEBUG_LOC=false};
01837
01838 if( DEBUG_LOC )
01839 {
01840 fprintf( ioQQQ, " %.2f Grain surface charge transfer rates\n", fnzone );
01841
01842 for( nelem=0; nelem < LIMELM; nelem++ )
01843 {
01844 if( dense.lgElmtOn[nelem] )
01845 {
01846 printf( " %s:", elementnames.chElementSym[nelem] );
01847 for( ion=dense.IonLow[nelem]; ion <= dense.IonHigh[nelem]; ++ion )
01848 {
01849 for( ion_to=0; ion_to <= nelem+1; ion_to++ )
01850 {
01851 if( gv.GrainChTrRate[nelem][ion][ion_to] > 0.f )
01852 {
01853 printf( " %ld->%ld %.2e", ion, ion_to,
01854 gv.GrainChTrRate[nelem][ion][ion_to] );
01855 }
01856 }
01857 }
01858 printf( "\n" );
01859 }
01860 }
01861 }
01862
01863 fprintf( ioQQQ, " %.2f Grain contribution to electron density %.2e\n",
01864 fnzone , gv.TotalEden );
01865
01866 fprintf( ioQQQ, " Grain electons: " );
01867 for( nd=0; nd < gv.nBin; nd++ )
01868 {
01869 double sum = 0.;
01870 for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
01871 {
01872 sum += gv.bin[nd]->chrg[nz]->FracPop*(double)gv.bin[nd]->chrg[nz]->DustZ*
01873 gv.bin[nd]->cnv_GR_pCM3;
01874 }
01875 fprintf( ioQQQ, " %.2e", sum );
01876 }
01877 fprintf( ioQQQ, "\n" );
01878
01879 fprintf( ioQQQ, " Grain potentials:" );
01880 for( nd=0; nd < gv.nBin; nd++ )
01881 {
01882 fprintf( ioQQQ, " %.2e", gv.bin[nd]->dstpot );
01883 }
01884 fprintf( ioQQQ, "\n" );
01885
01886 fprintf( ioQQQ, " Grain temperatures:" );
01887 for( nd=0; nd < gv.nBin; nd++ )
01888 {
01889 fprintf( ioQQQ, " %.2e", gv.bin[nd]->tedust );
01890 }
01891 fprintf( ioQQQ, "\n" );
01892
01893 fprintf( ioQQQ, " GrainCollCool: %.6e\n", gv.GasCoolColl );
01894 }
01895
01896
01897
01898
01899
01900
01901
01902
01903 DEBUG_EXIT( "GrainChargeTemp()" );
01904 return;
01905 }
01906
01907
01908 STATIC void GrainCharge(long int nd,
01909 double *ThermRatio)
01910 {
01911 bool lgBigError,
01912 lgInitial;
01913 long backup,
01914 i,
01915 loopMax,
01916 newZlo,
01917 nz,
01918 power,
01919 stride,
01920 stride0,
01921 Zlo;
01922 double crate,
01923 csum1,
01924 csum1a,
01925 csum1b,
01926 csum1c,
01927 csum2,
01928 csum3,
01929 netloss0 = -DBL_MAX,
01930 netloss1 = -DBL_MAX,
01931 rate_up[NCHU],
01932 rate_dn[NCHU],
01933 step;
01934
01935 DEBUG_ENTRY( "GrainCharge()" );
01936
01937
01938 if( trace.lgTrace && trace.lgDustBug )
01939 {
01940 fprintf( ioQQQ, " Charge loop, search %c,", TorF(conv.lgSearch) );
01941 }
01942
01943 ASSERT( nd >= 0 && nd < gv.nBin );
01944
01945 for( nz=0; nz < NCHS; nz++ )
01946 {
01947 gv.bin[nd]->chrg[nz]->FracPop = -DBL_MAX;
01948 }
01949
01950
01951
01952
01953
01954
01955
01956
01957
01958
01959
01960
01961
01962
01963
01964
01965
01966
01967
01968
01969
01970
01971
01972
01973
01974
01975
01976
01977
01978
01979
01980
01981
01982
01983
01984
01985
01986
01987
01988
01989
01990
01991
01992
01993
01994
01995
01996
01997
01998 lgInitial = ( gv.bin[nd]->chrg[0]->DustZ == LONG_MIN );
01999
02000 backup = gv.bin[nd]->nChrg;
02001 gv.bin[nd]->nChrg = MAX2(gv.bin[nd]->nChrg,3);
02002
02003 stride0 = gv.bin[nd]->nChrg-1;
02004
02005
02006 if( lgInitial )
02007 {
02008 double xxx;
02009 step = MAX2((double)(-gv.bin[nd]->LowestZg),1.);
02010 power = (int)(log(step)/log((double)stride0));
02011 power = MAX2(power,0);
02012 xxx = powi((double)stride0,power);
02013 stride = nint(xxx);
02014 Zlo = gv.bin[nd]->LowestZg;
02015 }
02016 else
02017 {
02018
02019 stride = 1;
02020 Zlo = gv.bin[nd]->chrg[0]->DustZ;
02021 }
02022 UpdatePot( nd, Zlo, stride, rate_up, rate_dn );
02023
02024
02025 for( i=0; i < BRACKET_MAX; i++ )
02026 {
02027 netloss0 = rate_up[0] - rate_dn[0];
02028 netloss1 = rate_up[gv.bin[nd]->nChrg-1] - rate_dn[gv.bin[nd]->nChrg-1];
02029
02030 if( netloss0*netloss1 <= 0. )
02031 break;
02032
02033 if( netloss1 > 0. )
02034 Zlo += (gv.bin[nd]->nChrg-1)*stride;
02035
02036 if( i > 0 )
02037 stride *= stride0;
02038
02039 if( netloss1 < 0. )
02040 Zlo -= (gv.bin[nd]->nChrg-1)*stride;
02041
02042 Zlo = MAX2(Zlo,gv.bin[nd]->LowestZg);
02043 UpdatePot( nd, Zlo, stride, rate_up, rate_dn );
02044 }
02045
02046 if( netloss0*netloss1 > 0. ) {
02047 fprintf( ioQQQ, " insanity: could not bracket grain charge for %s\n", gv.bin[nd]->chDstLab );
02048 ShowMe();
02049 puts( "[Stop in GrainCharge]" );
02050 cdEXIT(EXIT_FAILURE);
02051 }
02052
02053
02054 while( stride > 1 )
02055 {
02056 stride /= stride0;
02057
02058 netloss0 = rate_up[0] - rate_dn[0];
02059 for( nz=0; nz < gv.bin[nd]->nChrg-1; nz++ )
02060 {
02061 netloss1 = rate_up[nz+1] - rate_dn[nz+1];
02062
02063 if( netloss0*netloss1 <= 0. )
02064 {
02065 Zlo = gv.bin[nd]->chrg[nz]->DustZ;
02066 break;
02067 }
02068
02069 netloss0 = netloss1;
02070 }
02071 UpdatePot( nd, Zlo, stride, rate_up, rate_dn );
02072 }
02073
02074 ASSERT( netloss0*netloss1 <= 0. );
02075
02076 gv.bin[nd]->nChrg = backup;
02077
02078
02079 loopMax = ( lgInitial ) ? 4*gv.bin[nd]->nChrg : 2*gv.bin[nd]->nChrg;
02080
02081 lgBigError = true;
02082 for( i=0; i < loopMax; i++ )
02083 {
02084 GetFracPop( nd, Zlo, rate_up, rate_dn, &newZlo );
02085
02086 if( newZlo == Zlo )
02087 {
02088 lgBigError = false;
02089 break;
02090 }
02091
02092 Zlo = newZlo;
02093 UpdatePot( nd, Zlo, 1, rate_up, rate_dn );
02094 }
02095
02096 if( lgBigError ) {
02097 fprintf( ioQQQ, " insanity: could not converge grain charge for %s\n", gv.bin[nd]->chDstLab );
02098 ShowMe();
02099 puts( "[Stop in GrainCharge]" );
02100 cdEXIT(EXIT_FAILURE);
02101 }
02102
02105 gv.bin[nd]->lgChrgConverged = true;
02106
02107 gv.bin[nd]->AveDustZ = 0.;
02108 crate = csum3 = 0.;
02109 for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
02110 {
02111 double d[5];
02112
02113 gv.bin[nd]->AveDustZ += gv.bin[nd]->chrg[nz]->FracPop*gv.bin[nd]->chrg[nz]->DustZ;
02114
02115 crate += gv.bin[nd]->chrg[nz]->FracPop*GrainElecEmis1(nd,nz,&d[0],&d[1],&d[2],&d[3],&d[4]);
02116 csum3 += gv.bin[nd]->chrg[nz]->FracPop*d[4];
02117 }
02118 gv.bin[nd]->dstpot = CHRG2POT(gv.bin[nd]->AveDustZ);
02119 *ThermRatio = ( crate > 0. ) ? csum3/crate : 0.;
02120
02121 ASSERT( *ThermRatio >= 0. );
02122
02123 if( trace.lgTrace && trace.lgDustBug )
02124 {
02125 double d[5];
02126
02127 fprintf( ioQQQ, "\n" );
02128
02129 crate = csum1a = csum1b = csum1c = csum2 = csum3 = 0.;
02130 for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
02131 {
02132 crate += gv.bin[nd]->chrg[nz]->FracPop*GrainElecEmis1(nd,nz,&d[0],&d[1],&d[2],&d[3],&d[4]);
02133 csum1a += gv.bin[nd]->chrg[nz]->FracPop*d[0];
02134 csum1b += gv.bin[nd]->chrg[nz]->FracPop*d[1];
02135 csum1c += gv.bin[nd]->chrg[nz]->FracPop*d[2];
02136 csum2 += gv.bin[nd]->chrg[nz]->FracPop*d[3];
02137 csum3 += gv.bin[nd]->chrg[nz]->FracPop*d[4];
02138 }
02139
02140 fprintf( ioQQQ, " ElecEm rate1a=%.4e, rate1b=%.4e, rate1c=%.4e, ", csum1a, csum1b, csum1c );
02141 fprintf( ioQQQ, "rate2=%.4e, rate3=%.4e, sum=%.4e\n", csum2, csum3, crate );
02142 if( crate > 0. )
02143 {
02144 fprintf( ioQQQ, " rate1a/sum=%.4e, rate1b/sum=%.4e, rate1c/sum=%.4e, ",
02145 csum1a/crate, csum1b/crate, csum1c/crate );
02146 fprintf( ioQQQ, "rate2/sum=%.4e, rate3/sum=%.4e\n", csum2/crate, csum3/crate );
02147 }
02148
02149 crate = csum1 = csum2 = 0.;
02150 for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
02151 {
02152 crate += gv.bin[nd]->chrg[nz]->FracPop*GrainElecRecomb1(nd,nz,&d[0],&d[1]);
02153 csum1 += gv.bin[nd]->chrg[nz]->FracPop*d[0];
02154 csum2 += gv.bin[nd]->chrg[nz]->FracPop*d[1];
02155 }
02156
02157 fprintf( ioQQQ, " ElecRc rate1=%.4e, rate2=%.4e, sum=%.4e\n", csum1, csum2, crate );
02158 if( crate > 0. )
02159 {
02160 fprintf( ioQQQ, " rate1/sum=%.4e, rate2/sum=%.4e\n", csum1/crate, csum2/crate );
02161 }
02162
02163 fprintf( ioQQQ, " Charging rates:" );
02164 for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
02165 {
02166 fprintf( ioQQQ, " Zg %ld up %.4e dn %.4e",
02167 gv.bin[nd]->chrg[nz]->DustZ, rate_up[nz], rate_dn[nz] );
02168 }
02169 fprintf( ioQQQ, "\n" );
02170
02171 fprintf( ioQQQ, " FracPop: fnzone %.2f nd %ld AveZg %.5e", fnzone, nd, gv.bin[nd]->AveDustZ );
02172 for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
02173 {
02174 fprintf( ioQQQ, " Zg %ld %.5f", Zlo+nz, gv.bin[nd]->chrg[nz]->FracPop );
02175 }
02176 fprintf( ioQQQ, "\n" );
02177
02178 fprintf( ioQQQ, " >Grain potential:%12.12s %.6fV",
02179 gv.bin[nd]->chDstLab, gv.bin[nd]->dstpot*EVRYD );
02180 for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
02181 {
02182 fprintf( ioQQQ, " Thres[%ld]: %.4f V ThresVal[%ld]: %.4f V",
02183 gv.bin[nd]->chrg[nz]->DustZ, gv.bin[nd]->chrg[nz]->ThresInf*EVRYD,
02184 gv.bin[nd]->chrg[nz]->DustZ, gv.bin[nd]->chrg[nz]->ThresInfVal*EVRYD );
02185 }
02186 fprintf( ioQQQ, "\n" );
02187 }
02188
02189 DEBUG_EXIT( "GrainCharge()" );
02190 return;
02191 }
02192
02193
02194
02195 STATIC double GrainElecRecomb1(long nd,
02196 long nz,
02197 double *sum1,
02198 double *sum2)
02199 {
02200 long ion,
02201 nelem;
02202 double eta,
02203 rate,
02204 Stick,
02205 ve,
02206 xi;
02207
02208 DEBUG_ENTRY( "GrainElecRecomb1()" );
02209
02210 ASSERT( nd >= 0 && nd < gv.nBin );
02211 ASSERT( nz >= 0 && nz < gv.bin[nd]->nChrg );
02212
02213
02214
02215 if( gv.bin[nd]->chrg[nz]->RSum1 >= 0. )
02216 {
02217 *sum1 = gv.bin[nd]->chrg[nz]->RSum1;
02218 *sum2 = gv.bin[nd]->chrg[nz]->RSum2;
02219 rate = *sum1 + *sum2;
02220
02221 DEBUG_EXIT( "GrainElecRecomb1()" );
02222 return rate;
02223 }
02224
02225
02226 ion = -1;
02227
02228
02229 ve = sqrt(8.*BOLTZMANN/PI/ELECTRON_MASS*phycon.te);
02230
02231 Stick = ( gv.bin[nd]->chrg[nz]->DustZ <= -1 ) ? gv.bin[nd]->StickElecNeg : gv.bin[nd]->StickElecPos;
02232
02233
02234 GrainScreen(ion,nd,nz,&eta,&xi);
02235
02236 *sum1 = EPSP(gv.bin[nd]->LowestZg,gv.bin[nd]->chrg[nz]->DustZ)*Stick*dense.eden*ve*eta;
02237
02238
02239 *sum2 = 0.;
02240
02241 #ifndef IGNORE_GRAIN_ION_COLLISIONS
02242 for( ion=0; ion <= LIMELM; ion++ )
02243 {
02244 double CollisionRateAll = 0.;
02245
02246 for( nelem=MAX2(ion-1,0); nelem < LIMELM; nelem++ )
02247 {
02248 if( dense.lgElmtOn[nelem] && dense.xIonDense[nelem][ion] > 0. &&
02249 gv.bin[nd]->chrg[nz]->RecomZ0[nelem][ion] > ion )
02250 {
02251
02252 CollisionRateAll += STICK_ION*dense.xIonDense[nelem][ion]*DoppVel.AveVel[nelem]*
02253 (double)(gv.bin[nd]->chrg[nz]->RecomZ0[nelem][ion]-ion);
02254 }
02255 }
02256
02257 if( CollisionRateAll > 0. )
02258 {
02259
02260
02261
02262 GrainScreen(ion,nd,nz,&eta,&xi);
02263 *sum2 += CollisionRateAll*eta;
02264 }
02265 }
02266 #endif
02267
02268 rate = *sum1 + *sum2;
02269
02270
02271 gv.bin[nd]->chrg[nz]->RSum1 = *sum1;
02272 gv.bin[nd]->chrg[nz]->RSum2 = *sum2;
02273
02274 ASSERT( *sum1 >= 0. && *sum2 >= 0. );
02275
02276 DEBUG_EXIT( "GrainElecRecomb1()" );
02277 return rate;
02278 }
02279
02280
02281
02282 STATIC double GrainElecEmis1(long nd,
02283 long nz,
02284 double *sum1a,
02285 double *sum1b,
02286 double *sum1c,
02287 double *sum2,
02288 double *sum3)
02289 {
02290 long int i,
02291 ion,
02292 ipLo,
02293 nelem;
02294 double eta,
02295 rate,
02296 xi;
02297
02298 DEBUG_ENTRY( "GrainElecEmis1()" );
02299
02300 ASSERT( nd >= 0 && nd < gv.nBin );
02301 ASSERT( nz >= 0 && nz < gv.bin[nd]->nChrg );
02302
02303
02304
02305 if( gv.bin[nd]->chrg[nz]->ESum1a >= 0. )
02306 {
02307 *sum1a = gv.bin[nd]->chrg[nz]->ESum1a;
02308 *sum1b = gv.bin[nd]->chrg[nz]->ESum1b;
02309 *sum1c = gv.bin[nd]->chrg[nz]->ESum1c;
02310 *sum2 = gv.bin[nd]->chrg[nz]->ESum2;
02311
02312 *sum3 = 4.*gv.bin[nd]->chrg[nz]->ThermRate;
02313 rate = *sum1a + *sum1b + *sum1c + *sum2 + *sum3;
02314
02315 DEBUG_EXIT( "GrainElecEmis1()" );
02316 return rate;
02317 }
02318
02319
02320 *sum1a = 0.;
02321 ipLo = gv.bin[nd]->chrg[nz]->ipThresInfVal;
02322 for( i=ipLo; i < rfield.nflux; i++ )
02323 {
02324 # ifdef WD_TEST2
02325 *sum1a += rfield.flux[i]*gv.bin[nd]->dstab1[i]*gv.bin[nd]->chrg[nz]->yhat[i-ipLo];
02326 # else
02327 *sum1a += rfield.SummedCon[i]*gv.bin[nd]->dstab1[i]*gv.bin[nd]->chrg[nz]->yhat[i-ipLo];
02328 # endif
02329 }
02330
02331 *sum1a /= gv.bin[nd]->IntArea/4.;
02332
02333
02334
02335
02336
02337
02338
02339
02340
02345 *sum1b = 0.;
02346
02347 # ifdef WD_TEST2
02348 *sum1b = 0.;
02349 # else
02350 *sum1b /= gv.bin[nd]->IntArea/4.;
02351 # endif
02352
02353 *sum1c = 0.;
02354 if( gv.bin[nd]->chrg[nz]->DustZ <= -1 )
02355 {
02356 ipLo = gv.bin[nd]->chrg[nz]->ipThresInf;
02357 for( i=ipLo; i < rfield.nflux; i++ )
02358 {
02359
02360 # ifdef WD_TEST2
02361 *sum1c += rfield.flux[i]*gv.bin[nd]->chrg[nz]->cs_pdt[i-ipLo];
02362 # else
02363 *sum1c += rfield.SummedCon[i]*gv.bin[nd]->chrg[nz]->cs_pdt[i-ipLo];
02364 # endif
02365 }
02366 *sum1c /= gv.bin[nd]->IntArea/4.;
02367 }
02368
02369
02370 *sum2 = 0.;
02371 # ifndef IGNORE_GRAIN_ION_COLLISIONS
02372 for( ion=0; ion <= LIMELM; ion++ )
02373 {
02374 double CollisionRateAll = 0.;
02375
02376 for( nelem=MAX2(ion-1,0); nelem < LIMELM; nelem++ )
02377 {
02378 if( dense.lgElmtOn[nelem] && dense.xIonDense[nelem][ion] > 0. &&
02379 ion > gv.bin[nd]->chrg[nz]->RecomZ0[nelem][ion] )
02380 {
02381
02382 CollisionRateAll += STICK_ION*dense.xIonDense[nelem][ion]*DoppVel.AveVel[nelem]*
02383 (double)(ion-gv.bin[nd]->chrg[nz]->RecomZ0[nelem][ion]);
02384 }
02385 }
02386
02387 if( CollisionRateAll > 0. )
02388 {
02389
02390
02391
02392 GrainScreen(ion,nd,nz,&eta,&xi);
02393 *sum2 += CollisionRateAll*eta;
02394 }
02395 }
02396 # endif
02397
02398
02399
02400
02401 *sum3 = 4.*gv.bin[nd]->chrg[nz]->ThermRate;
02402
02403 rate = *sum1a + *sum1b + *sum1c + *sum2 + *sum3;
02404
02405
02406 gv.bin[nd]->chrg[nz]->ESum1a = *sum1a;
02407 gv.bin[nd]->chrg[nz]->ESum1b = *sum1b;
02408 gv.bin[nd]->chrg[nz]->ESum1c = *sum1c;
02409 gv.bin[nd]->chrg[nz]->ESum2 = *sum2;
02410
02411 ASSERT( *sum1a >= 0. && *sum1b >= 0. && *sum1c >= 0. && *sum2 >= 0. && *sum3 >= 0. );
02412
02413 DEBUG_EXIT( "GrainElecEmis1()" );
02414 return rate;
02415 }
02416
02417
02418
02419
02420 STATIC void GrainScreen(long ion,
02421 long nd,
02422 long nz,
02423 double *eta,
02424 double *xi)
02425 {
02426
02427
02428
02429
02430 long ind = ion+1;
02431
02432 DEBUG_ENTRY( "GrainScreen()" );
02433
02434 ASSERT( ind >= 0 && ind < LIMELM+2 );
02435
02436 if( gv.bin[nd]->chrg[nz]->eta[ind] > 0. )
02437 {
02438 *eta = gv.bin[nd]->chrg[nz]->eta[ind];
02439 *xi = gv.bin[nd]->chrg[nz]->xi[ind];
02440
02441 DEBUG_EXIT( "GrainScreen()" );
02442 return;
02443 }
02444
02445
02446
02447 if( ion == 0 )
02448 {
02449 *eta = 1.;
02450 *xi = 1.;
02451 }
02452 else
02453 {
02454
02455
02456
02457
02458
02459
02460 double nu = (double)gv.bin[nd]->chrg[nz]->DustZ/(double)ion;
02461 double tau = gv.bin[nd]->Capacity*BOLTZMANN*phycon.te*1.e-7/POW2((double)ion*ELEM_CHARGE);
02462 if( nu < 0. )
02463 {
02464 *eta = (1. - nu/tau)*(1. + sqrt(2./(tau - 2.*nu)));
02465 *xi = (1. - nu/(2.*tau))*(1. + 1./sqrt(tau - nu));
02466 }
02467 else if( nu == 0. )
02468 {
02469 *eta = 1. + sqrt(PI/(2.*tau));
02470 *xi = 1. + 0.75*sqrt(PI/(2.*tau));
02471 }
02472 else
02473 {
02474 double theta_nu = ThetaNu(nu);
02475
02476 double xxx = 1. + 1./sqrt(4.*tau+3.*nu);
02477 *eta = POW2(xxx)*exp(-theta_nu/tau);
02478 # ifdef WD_TEST2
02479 *xi = (1. + nu/(2.*tau))*(1. + 1./sqrt(3./(2.*tau)+3.*nu))*exp(-theta_nu/tau);
02480 # else
02481
02482
02483
02484 xxx = 0.25*pow(nu/tau,0.75)/(pow(nu/tau,0.75) + pow((25.+3.*nu)/5.,0.75)) +
02485 (1. + 0.75*sqrt(PI/(2.*tau)))/(1. + sqrt(PI/(2.*tau)));
02486 *xi = (MIN2(xxx,1.) + theta_nu/(2.*tau))*(*eta);
02487 # endif
02488 }
02489
02490 ASSERT( *eta >= 0. && *xi >= 0. );
02491 }
02492
02493 gv.bin[nd]->chrg[nz]->eta[ind] = *eta;
02494 gv.bin[nd]->chrg[nz]->xi[ind] = *xi;
02495
02496 DEBUG_EXIT( "GrainScreen()" );
02497
02498 return;
02499 }
02500
02501
02502 STATIC double ThetaNu(double nu)
02503 {
02504 double theta_nu;
02505
02506 DEBUG_ENTRY( "ThetaNu()" );
02507
02508 if( nu > 0. )
02509 {
02510 double xxx;
02511 const double REL_TOLER = 10.*DBL_EPSILON;
02512
02513
02514 double xi_nu = 1. + 1./sqrt(3.*nu);
02515 double xi_nu2 = POW2(xi_nu);
02516 do
02517 {
02518 double old = xi_nu;
02519
02520
02521 double fnu = 2.*xi_nu2 - 1. - nu*xi_nu*POW2(xi_nu2 - 1.);
02522 double dfdxi = 4.*xi_nu - nu*((5.*xi_nu2 - 6.)*xi_nu2 + 1.);
02523 xi_nu -= fnu/dfdxi;
02524 xi_nu2 = POW2(xi_nu);
02525 xxx = fabs(old-xi_nu);
02526 } while( xxx > REL_TOLER*xi_nu );
02527
02528 theta_nu = nu/xi_nu - 1./(2.*xi_nu2*(xi_nu2-1.));
02529 }
02530 else
02531 {
02532 theta_nu = 0.;
02533 }
02534
02535 DEBUG_EXIT( "ThetaNu()" );
02536 return theta_nu;
02537 }
02538
02539
02540
02541 STATIC void UpdatePot(long nd,
02542 long Zlo,
02543 long stride,
02544 double rate_up[],
02545 double rate_dn[])
02546 {
02547 long nz,
02548 Zg;
02549 double BoltzFac,
02550 HighEnergy;
02551
02552 DEBUG_ENTRY( "UpdatePot()" );
02553
02554 ASSERT( nd >= 0 && nd < gv.nBin );
02555 ASSERT( Zlo >= gv.bin[nd]->LowestZg );
02556 ASSERT( stride >= 1 );
02557
02558
02559
02560
02561
02562
02563
02564
02565
02566
02567 if( trace.lgTrace && trace.lgDustBug )
02568 {
02569 fprintf( ioQQQ, " %ld/%ld", Zlo, stride );
02570 }
02571
02572 for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
02573 {
02574 long ind, zz;
02575 double d[5];
02576 ChargeBin *ptr;
02577
02578 Zg = Zlo + nz*stride;
02579
02580
02581 ind = NCHS-1;
02582 for( zz=0; zz < NCHS-1; zz++ )
02583 {
02584 if( gv.bin[nd]->chrg[zz]->DustZ == Zg )
02585 {
02586 ind = zz;
02587 break;
02588 }
02589 }
02590
02591
02592
02593 ptr = gv.bin[nd]->chrg[ind];
02594
02595
02596 for( zz=ind-1; zz >= nz; zz-- )
02597 gv.bin[nd]->chrg[zz+1] = gv.bin[nd]->chrg[zz];
02598
02599 gv.bin[nd]->chrg[nz] = ptr;
02600
02601 if( gv.bin[nd]->chrg[nz]->DustZ != Zg )
02602 UpdatePot1(nd,nz,Zg,0);
02603 else if( gv.bin[nd]->chrg[nz]->nfill < rfield.nflux )
02604 UpdatePot1(nd,nz,Zg,gv.bin[nd]->chrg[nz]->nfill);
02605
02606 UpdatePot2(nd,nz);
02607
02608 rate_up[nz] = GrainElecEmis1(nd,nz,&d[0],&d[1],&d[2],&d[3],&d[4]);
02609 rate_dn[nz] = GrainElecRecomb1(nd,nz,&d[0],&d[1]);
02610
02611
02612 ASSERT( gv.bin[nd]->chrg[nz]->DustZ == Zg );
02613 ASSERT( gv.bin[nd]->chrg[nz]->nfill >= rfield.nflux );
02614 ASSERT( rate_up[nz] >= 0. && rate_dn[nz] >= 0. );
02615 }
02616
02617
02618
02619
02620
02621
02622 BoltzFac = (-log(CONSERV_TOL) + 8.)*BOLTZMANN/EN1RYD;
02623 HighEnergy = 0.;
02624 for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
02625 {
02626
02627 HighEnergy = MAX2(HighEnergy,
02628 MAX2(gv.bin[nd]->chrg[nz]->ThresInfInc,0.) + BoltzFac*MAX2(phycon.te,gv.bin[nd]->tedust));
02629 }
02630 HighEnergy = MIN2(HighEnergy,rfield.anu[rfield.nupper-1]);
02631 gv.bin[nd]->qnflux2 = ipoint(HighEnergy);
02632 gv.bin[nd]->qnflux = MAX2(rfield.nflux,gv.bin[nd]->qnflux2);
02633
02634 ASSERT( gv.bin[nd]->qnflux <= rfield.nupper );
02635
02636 DEBUG_EXIT( "UpdatePot()" );
02637 return;
02638 }
02639
02640
02641
02642 STATIC void GetFracPop(long nd,
02643 long Zlo,
02644 double rate_up[],
02645 double rate_dn[],
02646 long *newZlo)
02647 {
02648 bool lgRedo;
02649 long i,
02650 nz;
02651 double netloss[2],
02652 pop[2][NCHU-1];
02653
02654
02655 DEBUG_ENTRY( "GetFracPop()" );
02656
02657 ASSERT( nd >= 0 && nd < gv.nBin );
02658 ASSERT( Zlo >= gv.bin[nd]->LowestZg );
02659
02660
02661
02662
02663
02664
02665
02666
02667 do
02668 {
02669 for( i=0; i < 2; i++ )
02670 {
02671 long j, k;
02672 double sum;
02673
02674 sum = pop[i][0] = 1.;
02675 for( j=1; j < gv.bin[nd]->nChrg-1; j++ )
02676 {
02677 nz = i + j;
02678 if( rate_dn[nz] > 10.*rate_up[nz-1]/sqrt(DBL_MAX) )
02679 {
02680 pop[i][j] = pop[i][j-1]*rate_up[nz-1]/rate_dn[nz];
02681 sum += pop[i][j];
02682 }
02683 else
02684 {
02685 for( k=0; k < j; k++ )
02686 {
02687 pop[i][k] = 0.;
02688 }
02689 pop[i][j] = 1.;
02690 sum = pop[i][j];
02691 }
02692
02693 if( pop[i][j] > sqrt(DBL_MAX) )
02694 {
02695 for( k=0; k <= j; k++ )
02696 {
02697 pop[i][k] /= DBL_MAX/10.;
02698 }
02699 sum /= DBL_MAX/10.;
02700 }
02701 }
02702 netloss[i] = 0.;
02703 for( j=0; j < gv.bin[nd]->nChrg-1; j++ )
02704 {
02705 nz = i + j;
02706 pop[i][j] /= sum;
02707 netloss[i] += pop[i][j]*(rate_up[nz] - rate_dn[nz]);
02708 }
02709 }
02710
02711
02712
02713 if( netloss[0]*netloss[1] > 0. )
02714 *newZlo = ( netloss[1] > 0. ) ? Zlo + 1 : Zlo - 1;
02715 else
02716 *newZlo = Zlo;
02717
02718
02719
02720
02721
02722
02723 if( gv.bin[nd]->nChrg > 2 &&
02724 ( *newZlo < gv.bin[nd]->LowestZg ||
02725 ( *newZlo == Zlo && pop[1][gv.bin[nd]->nChrg-2] < DBL_EPSILON ) ) )
02726 {
02727 gv.bin[nd]->nChrg--;
02728 lgRedo = true;
02729 }
02730 else
02731 {
02732 lgRedo = false;
02733 }
02734
02735 # if 0
02736 printf( " fnzone %.2f nd %ld Zlo %ld newZlo %ld netloss %.4e %.4e nChrg %ld lgRedo %d\n",
02737 fnzone, nd, Zlo, *newZlo, netloss[0], netloss[1], gv.bin[nd]->nChrg, lgRedo );
02738 # endif
02739 }
02740 while( lgRedo );
02741
02742 if( *newZlo < gv.bin[nd]->LowestZg )
02743 {
02744 fprintf( ioQQQ, " could not converge charge state populations for %s\n", gv.bin[nd]->chDstLab );
02745 ShowMe();
02746 puts( "[Stop in GetFracPop]" );
02747 cdEXIT(EXIT_FAILURE);
02748 }
02749
02750 if( *newZlo == Zlo )
02751 {
02752 double frac0, frac1;
02753 # ifndef NDEBUG
02754 double test1, test2, test3, x1, x2;
02755 # endif
02756
02757
02758 frac0 = netloss[1]/(netloss[1]-netloss[0]);
02759 frac1 = -netloss[0]/(netloss[1]-netloss[0]);
02760
02761 gv.bin[nd]->chrg[0]->FracPop = frac0*pop[0][0];
02762 gv.bin[nd]->chrg[gv.bin[nd]->nChrg-1]->FracPop = frac1*pop[1][gv.bin[nd]->nChrg-2];
02763 for( nz=1; nz < gv.bin[nd]->nChrg-1; nz++ )
02764 {
02765 gv.bin[nd]->chrg[nz]->FracPop = frac0*pop[0][nz] + frac1*pop[1][nz-1];
02766 }
02767
02768 # ifndef NDEBUG
02769 test1 = test2 = test3 = 0.;
02770 for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
02771 {
02772 ASSERT( gv.bin[nd]->chrg[nz]->FracPop >= 0. );
02773 test1 += gv.bin[nd]->chrg[nz]->FracPop;
02774 test2 += gv.bin[nd]->chrg[nz]->FracPop*rate_up[nz];
02775 test3 += gv.bin[nd]->chrg[nz]->FracPop*(rate_up[nz]-rate_dn[nz]);
02776 }
02777 x1 = fabs(test1-1.);
02778 x2 = fabs(test3/test2);
02779 ASSERT( MAX2(x1,x2) < 10.*sqrt((double)gv.bin[nd]->nChrg)*DBL_EPSILON );
02780 # endif
02781 }
02782
02783 DEBUG_EXIT( "GetFracPop()" );
02784 return;
02785 }
02786
02787
02788
02789
02790
02791
02792
02793
02794
02795
02796
02797 STATIC void UpdatePot1(long nd,
02798 long nz,
02799 long Zg,
02800 long ipStart)
02801 {
02802 long i,
02803 ipLo,
02804 ipHi,
02805 nfill;
02806 size_t alloc_size;
02807 pe_type pcase;
02808 double c1,
02809 cnv_GR_pH,
02810 d[2],
02811 DustWorkFcn,
02812 Elo,
02813 PotSurf,
02814 ThresInf,
02815 ThresInfVal,
02816 ThresSurfVal;
02817
02818 float *yone = gv.bin[nd]->y1;
02819 float *anu = rfield.anu;
02820 float *yhat;
02821 double *cs_pdt, *fac1, *fac2;
02822
02823
02824
02825
02826 const double INV_DELTA_E = EVRYD/3.;
02827 const double CS_PDT = 1.2e-17;
02828
02829 DEBUG_ENTRY( "UpdatePot1()" );
02830
02831
02832
02833
02834
02835
02836
02837 if( ipStart == 0 )
02838 {
02839 gv.bin[nd]->chrg[nz]->DustZ = Zg;
02840
02841
02842 memset( gv.bin[nd]->chrg[nz]->eta, 0, (LIMELM+2)*sizeof(double) );
02843 memset( gv.bin[nd]->chrg[nz]->xi, 0, (LIMELM+2)*sizeof(double) );
02844
02845 GetPotValues(nd,Zg,&gv.bin[nd]->chrg[nz]->ThresInf,&gv.bin[nd]->chrg[nz]->ThresInfVal,
02846 &gv.bin[nd]->chrg[nz]->ThresSurf,&gv.bin[nd]->chrg[nz]->ThresSurfVal,
02847 &gv.bin[nd]->chrg[nz]->PotSurf,INCL_TUNNEL);
02848
02849
02850
02851 GetPotValues(nd,Zg-1,&gv.bin[nd]->chrg[nz]->ThresInfInc,&d[0],&gv.bin[nd]->chrg[nz]->ThresSurfInc,
02852 &d[1],&gv.bin[nd]->chrg[nz]->PotSurfInc,NO_TUNNEL);
02853
02854 ipLo = 0;
02855
02856 ipHi = rfield.nupper-1;
02857
02858 while( ipHi-ipLo > 1 )
02859 {
02860 long ipMd = (ipLo+ipHi)/2;
02861 if( gv.anumin[ipMd] > (float)gv.bin[nd]->chrg[nz]->ThresInfVal )
02862 ipHi = ipMd;
02863 else
02864 ipLo = ipMd;
02865 }
02866 gv.bin[nd]->chrg[nz]->ipThresInfVal = ipLo;
02867 }
02868 else
02869 {
02870 ipLo = gv.bin[nd]->chrg[nz]->ipThresInfVal;
02871 }
02872
02873
02874 gv.bin[nd]->chrg[nz]->nfill = rfield.nflux;
02875 nfill = gv.bin[nd]->chrg[nz]->nfill;
02876
02877
02878 FREE_SAFE( gv.bin[nd]->chrg[nz]->yhat );
02879
02880 if( nfill > ipLo )
02881 {
02882 alloc_size = (size_t)((unsigned)(nfill-ipLo)*sizeof(float));
02883
02884 gv.bin[nd]->chrg[nz]->yhat = (float*)MALLOC(alloc_size);
02885
02886 yhat = gv.bin[nd]->chrg[nz]->yhat;
02887 ThresSurfVal = gv.bin[nd]->chrg[nz]->ThresSurfVal;
02888 DustWorkFcn = gv.bin[nd]->DustWorkFcn;
02889 pcase = gv.which_pe[gv.bin[nd]->matType];
02890 Elo = -gv.bin[nd]->chrg[nz]->PotSurf;
02891 ThresInfVal = gv.bin[nd]->chrg[nz]->ThresInfVal;
02892
02893
02894 for( i=ipLo; i < nfill; i++ )
02895 {
02896 double xv,yzero,y2,Ehi;
02897 xv = MAX2((anu[i] - ThresSurfVal)/DustWorkFcn,0.);
02898
02899 switch( pcase )
02900 {
02901 case PE_CAR:
02902
02903 xv = xv*POW2(xv)*POW2(xv);
02904 yzero = xv/((1./9.e-3) + (3.7e-2/9.e-3)*xv);
02905 break;
02906 case PE_SIL:
02907
02908 yzero = xv/(2.+10.*xv);
02909 break;
02910 default:
02911 fprintf( ioQQQ, " UpdatePot1 detected unknown type for PE effect: %d\n" , pcase );
02912 puts( "[Stop in UpdatePot1]" );
02913 cdEXIT(1);
02914 }
02915 if( Zg > -1 )
02916 {
02917 Ehi = anu[i] - ThresInfVal;
02918 y2 = ( Ehi > 0. ) ? POW2(Ehi) * (Ehi - 3.*Elo) / POW3(Ehi - Elo) : 0.;
02919 }
02920 else
02921 {
02922 y2 = 1.;
02923 }
02924 yhat[i-ipLo] = (float)(y2*MIN2(yzero*yone[i],1.));
02925 }
02926 }
02927
02928 if( ipStart == 0 )
02929 {
02930
02931
02932
02933 ipLo = 0;
02934
02935 ipHi = rfield.nupper-1;
02936
02937 while( ipHi-ipLo > 1 )
02938 {
02939 long ipMd = (ipLo+ipHi)/2;
02940 if( gv.anumin[ipMd] > (float)gv.bin[nd]->chrg[nz]->ThresInf )
02941 ipHi = ipMd;
02942 else
02943 ipLo = ipMd;
02944 }
02945 gv.bin[nd]->chrg[nz]->ipThresInf = ipLo;
02946 }
02947 else
02948 {
02949 ipLo = gv.bin[nd]->chrg[nz]->ipThresInf;
02950 }
02951
02952 if( Zg > -1 )
02953 ipLo = nfill;
02954
02955
02956 FREE_SAFE( gv.bin[nd]->chrg[nz]->cs_pdt );
02957
02958 if( nfill > ipLo )
02959 {
02960 alloc_size = (size_t)((unsigned)(nfill-ipLo)*sizeof(double));
02961
02962 gv.bin[nd]->chrg[nz]->cs_pdt = (double*)MALLOC(alloc_size);
02963
02964 cs_pdt = gv.bin[nd]->chrg[nz]->cs_pdt;
02965 c1 = -CS_PDT*(double)Zg;
02966 ThresInf = gv.bin[nd]->chrg[nz]->ThresInf;
02967 cnv_GR_pH = gv.bin[nd]->cnv_GR_pH;
02968
02969 for( i=ipLo; i < nfill; i++ )
02970 {
02971 double x = (anu[i] - ThresInf)*INV_DELTA_E;
02972 double cs = c1*x/POW2(1.+(1./3.)*POW2(x));
02973
02974
02975 cs_pdt[i-ipLo] = MAX2(cs,0.)*cnv_GR_pH;
02976 }
02977 }
02978
02979
02980 FREE_SAFE( gv.bin[nd]->chrg[nz]->fac1 );
02981 FREE_SAFE( gv.bin[nd]->chrg[nz]->fac2 );
02982
02983 ipLo = gv.bin[nd]->chrg[nz]->ipThresInf;
02984 if( nfill > ipLo )
02985 {
02986 alloc_size = (size_t)((unsigned)(nfill-ipLo)*sizeof(double));
02987
02988 gv.bin[nd]->chrg[nz]->fac1 = (double*)MALLOC(alloc_size);
02989 gv.bin[nd]->chrg[nz]->fac2 = (double*)MALLOC(alloc_size);
02990
02991 fac1 = gv.bin[nd]->chrg[nz]->fac1;
02992 fac2 = gv.bin[nd]->chrg[nz]->fac2;
02993 PotSurf = gv.bin[nd]->chrg[nz]->PotSurf;
02994
02995 for( i=ipLo; i < nfill; i++ )
02996 {
02997 double cs1,cs2,cs_tot,cool1,cool2,ehat1,ehat2;
02998
02999
03000 PE_init(nd,nz,i,&cs1,&cs2,&cs_tot,&cool1,&cool2,&ehat1,&ehat2);
03001
03002 fac1[i-ipLo] = MAX2(cs_tot*anu[i]-cs1*cool1-cs2*cool2,0.);
03003 fac2[i-ipLo] = cs1*(ehat1-PotSurf) + cs2*(ehat2-PotSurf);
03004 }
03005 }
03006
03007 if( ipStart == 0 )
03008 {
03009
03010
03011
03012 UpdateRecomZ0(nd,nz,ALL_STAGES);
03013 }
03014
03015
03016 gv.bin[nd]->chrg[nz]->FracPop = -DBL_MAX;
03017
03018 gv.bin[nd]->chrg[nz]->RSum1 = -DBL_MAX;
03019 gv.bin[nd]->chrg[nz]->RSum2 = -DBL_MAX;
03020 gv.bin[nd]->chrg[nz]->ESum1a = -DBL_MAX;
03021 gv.bin[nd]->chrg[nz]->ESum1b = -DBL_MAX;
03022 gv.bin[nd]->chrg[nz]->ESum1c = -DBL_MAX;
03023 gv.bin[nd]->chrg[nz]->ESum2 = -DBL_MAX;
03024
03025 gv.bin[nd]->chrg[nz]->tedust = 1.f;
03026
03027 gv.bin[nd]->chrg[nz]->hcon1 = -DBL_MAX;
03028 gv.bin[nd]->chrg[nz]->hots1 = -DBL_MAX;
03029 gv.bin[nd]->chrg[nz]->bolflux1 = -DBL_MAX;
03030 gv.bin[nd]->chrg[nz]->pe1 = -DBL_MAX;
03031
03032 gv.bin[nd]->chrg[nz]->BolFlux = -DBL_MAX;
03033 gv.bin[nd]->chrg[nz]->GrainHeat = -DBL_MAX;
03034 gv.bin[nd]->chrg[nz]->GrainHeatColl = -DBL_MAX;
03035 gv.bin[nd]->chrg[nz]->GasHeatPhotoEl = -DBL_MAX;
03036 gv.bin[nd]->chrg[nz]->GasHeatTherm = -DBL_MAX;
03037 gv.bin[nd]->chrg[nz]->GrainCoolTherm = -DBL_MAX;
03038 gv.bin[nd]->chrg[nz]->ChemEnIon = -DBL_MAX;
03039 gv.bin[nd]->chrg[nz]->ChemEnH2 = -DBL_MAX;
03040
03041 gv.bin[nd]->chrg[nz]->HeatingRate2 = -DBL_MAX;
03042
03043
03044 ASSERT( gv.bin[nd]->chrg[nz]->ipThresInf <= gv.bin[nd]->chrg[nz]->ipThresInfVal );
03045
03046 DEBUG_EXIT( "UpdatePot1()" );
03047 return;
03048 }
03049
03050
03051
03052
03053
03054
03055
03056
03057
03058
03059
03060 STATIC void UpdatePot2(long nd,
03061 long nz)
03062 {
03063 double ThermExp;
03064
03065 DEBUG_ENTRY( "UpdatePot2()" );
03066
03067
03068 ThermExp = gv.bin[nd]->chrg[nz]->ThresInf*TE1RYD/gv.bin[nd]->tedust;
03069
03070 gv.bin[nd]->chrg[nz]->ThermRate = THERMCONST*gv.bin[nd]->ThermEff*POW2(gv.bin[nd]->tedust)*exp(-ThermExp);
03071 # if defined( WD_TEST2 ) || defined( IGNORE_THERMIONIC )
03072 gv.bin[nd]->chrg[nz]->ThermRate = 0.;
03073 # endif
03074
03075 DEBUG_EXIT( "UpdatePot2()" );
03076 return;
03077 }
03078
03079
03080
03081 STATIC long HighestIonStage(void)
03082 {
03083 long high,
03084 ion,
03085 nelem;
03086
03087 DEBUG_ENTRY( "HighestIonStage()" );
03088
03089 high = 0;
03090 for( nelem=LIMELM-1; nelem >= 0; nelem-- )
03091 {
03092 if( dense.lgElmtOn[nelem] )
03093 {
03094 for( ion=nelem+1; ion >= 0; ion-- )
03095 {
03096 if( ion == high || dense.xIonDense[nelem][ion] > 0. )
03097 break;
03098 }
03099 high = MAX2(high,ion);
03100 }
03101 if( nelem <= high )
03102 break;
03103 }
03104
03105 DEBUG_EXIT( "HighestIonStage()" );
03106 return high;
03107 }
03108
03109
03110 STATIC void UpdateRecomZ0(long nd,
03111 long nz,
03112 bool lgAllIonStages)
03113 {
03114 long hi_ion,
03115 i,
03116 ion,
03117 nelem,
03118 Zg;
03119 double d[4],
03120 phi_s_up[LIMELM+1],
03121 phi_s_dn[2];
03122
03123 DEBUG_ENTRY( "UpdateRecomZ0()" );
03124
03125 Zg = gv.bin[nd]->chrg[nz]->DustZ;
03126
03127 hi_ion = ( lgAllIonStages ) ? LIMELM : gv.HighestIon;
03128
03129 phi_s_up[0] = gv.bin[nd]->chrg[nz]->ThresSurf;
03130 for( i=1; i <= LIMELM; i++ )
03131 {
03132 if( i <= hi_ion )
03133 GetPotValues(nd,Zg+i,&d[0],&d[1],&phi_s_up[i],&d[2],&d[3],INCL_TUNNEL);
03134 else
03135 phi_s_up[i] = -DBL_MAX;
03136 }
03137 phi_s_dn[0] = gv.bin[nd]->chrg[nz]->ThresSurfInc;
03138 GetPotValues(nd,Zg-2,&d[0],&d[1],&phi_s_dn[1],&d[2],&d[3],NO_TUNNEL);
03139
03140
03141 for( nelem=0; nelem < LIMELM; nelem++ )
03142 {
03143 if( dense.lgElmtOn[nelem] )
03144 {
03145 for( ion=0; ion <= nelem+1; ion++ )
03146 {
03147 if( lgAllIonStages || dense.xIonDense[nelem][ion] > 0. )
03148 {
03149 GrainIonColl(nd,nz,nelem,ion,phi_s_up,phi_s_dn,
03150 &gv.bin[nd]->chrg[nz]->RecomZ0[nelem][ion],
03151 &gv.bin[nd]->chrg[nz]->RecomEn[nelem][ion],
03152 &gv.bin[nd]->chrg[nz]->ChemEn[nelem][ion]);
03153 }
03154 else
03155 {
03156 gv.bin[nd]->chrg[nz]->RecomZ0[nelem][ion] = ion;
03157 gv.bin[nd]->chrg[nz]->RecomEn[nelem][ion] = 0.f;
03158 gv.bin[nd]->chrg[nz]->ChemEn[nelem][ion] = 0.f;
03159 }
03160 }
03161 }
03162 }
03163
03164 DEBUG_EXIT( "UpdateRecomZ0()" );
03165 return;
03166 }
03167
03168 STATIC void GetPotValues(long nd,
03169 long Zg,
03170 double *ThresInf,
03171 double *ThresInfVal,
03172 double *ThresSurf,
03173 double *ThresSurfVal,
03174 double *PotSurf,
03175 bool lgUseTunnelCorr)
03176 {
03177 double dstpot,
03178 dZg = (double)Zg,
03179 IP_v;
03180
03181 DEBUG_ENTRY( "GetPotValues()" );
03182
03183
03184
03185
03186
03187
03188 dstpot = CHRG2POT(dZg);
03189
03190
03191
03192
03193 IP_v = gv.bin[nd]->DustWorkFcn + dstpot - 0.5*ONE_ELEC + (dZg+2.)*AC0/gv.bin[nd]->AvRadius*ONE_ELEC;
03194
03195
03196
03197
03198
03199 if( Zg <= -1 )
03200 {
03201 pot_type pcase = gv.which_pot[gv.bin[nd]->matType];
03202 double Emin,IP;
03203
03204 IP = gv.bin[nd]->DustWorkFcn - gv.bin[nd]->BandGap + dstpot - 0.5*ONE_ELEC;
03205 switch( pcase )
03206 {
03207 case POT_CAR:
03208 IP -= AC1G/(gv.bin[nd]->AvRadius+AC2G)*ONE_ELEC;
03209 break;
03210 case POT_SIL:
03211
03212 break;
03213 default:
03214 fprintf( ioQQQ, " GetPotValues detected unknown type for ionization pot: %d\n" , pcase );
03215 puts( "[Stop in GetPotValues]" );
03216 cdEXIT(1);
03217 }
03218
03219
03220
03221 IP_v = MAX2(IP,IP_v);
03222
03223 if( Zg < -1 )
03224 {
03225
03226 double help = fabs(dZg+1);
03227
03228 Emin = -ThetaNu(help)*ONE_ELEC;
03229 if( lgUseTunnelCorr )
03230 {
03231
03232 Emin *= 1. - 2.124e-4/(pow(gv.bin[nd]->AvRadius,0.45f)*pow(help,0.26));
03233 }
03234 }
03235 else
03236 {
03237 Emin = 0.;
03238 }
03239
03240 *ThresInf = IP - Emin;
03241 *ThresInfVal = IP_v - Emin;
03242 *ThresSurf = *ThresInf;
03243 *ThresSurfVal = *ThresInfVal;
03244 *PotSurf = Emin;
03245 }
03246 else
03247 {
03248 *ThresInf = IP_v;
03249 *ThresInfVal = IP_v;
03250 *ThresSurf = *ThresInf - dstpot;
03251 *ThresSurfVal = *ThresInfVal - dstpot;
03252 *PotSurf = dstpot;
03253 }
03254
03255 DEBUG_EXIT( "GetPotValues()" );
03256 return;
03257 }
03258
03259
03260
03261
03262
03263 STATIC void GrainIonColl(long int nd,
03264 long int nz,
03265 long int nelem,
03266 long int ion,
03267 const double phi_s_up[],
03268 const double phi_s_dn[],
03269 long *Z0,
03270 float *ChEn,
03271 float *ChemEn)
03272 {
03273 long Zg;
03274 double d[4];
03275 double phi_s;
03276
03277 long save = ion;
03278
03279 DEBUG_ENTRY( "GrainIonColl()" );
03280 if( ion > 0 && rfield.anu[Heavy.ipHeavy[nelem][ion-1]-1] > (float)phi_s_up[0] )
03281 {
03282
03283 *ChEn = 0.f;
03284 *ChemEn = 0.f;
03285 Zg = gv.bin[nd]->chrg[nz]->DustZ;
03286 phi_s = phi_s_up[0];
03287 do
03288 {
03289 *ChEn += rfield.anu[Heavy.ipHeavy[nelem][ion-1]-1] - (float)phi_s;
03290 *ChemEn += rfield.anu[Heavy.ipHeavy[nelem][ion-1]-1];
03291
03292
03293
03294 *ChemEn -= (float)(phi_s - phi_s_up[0]);
03295 --ion;
03296 ++Zg;
03297 phi_s = phi_s_up[save-ion];
03298 } while( ion > 0 && rfield.anu[Heavy.ipHeavy[nelem][ion-1]-1] > (float)phi_s );
03299
03300 *Z0 = ion;
03301 }
03302 else if( ion <= nelem && gv.bin[nd]->chrg[nz]->DustZ > gv.bin[nd]->LowestZg &&
03303 rfield.anu[Heavy.ipHeavy[nelem][ion]-1] < (float)phi_s_dn[0] )
03304 {
03305
03306 *ChEn = 0.f;
03307 *ChemEn = 0.f;
03308 Zg = gv.bin[nd]->chrg[nz]->DustZ;
03309 phi_s = phi_s_dn[0];
03310 do
03311 {
03312 *ChEn += (float)phi_s - rfield.anu[Heavy.ipHeavy[nelem][ion]-1];
03313 *ChemEn -= rfield.anu[Heavy.ipHeavy[nelem][ion]-1];
03314
03315
03316
03317 *ChemEn += (float)(phi_s - phi_s_dn[0]);
03318 ++ion;
03319 --Zg;
03320
03321 if( ion-save < 2 )
03322 phi_s = phi_s_dn[ion-save];
03323 else
03324 GetPotValues(nd,Zg-1,&d[0],&d[1],&phi_s,&d[2],&d[3],NO_TUNNEL);
03325
03326 } while( ion <= nelem && Zg > gv.bin[nd]->LowestZg &&
03327 rfield.anu[Heavy.ipHeavy[nelem][ion]-1] < (float)phi_s );
03328 *Z0 = ion;
03329 }
03330 else
03331 {
03332
03333 *ChEn = 0.f;
03334 *ChemEn = 0.f;
03335 *Z0 = ion;
03336 }
03337
03338
03339 DEBUG_EXIT( "GrainIonColl()" );
03340 return;
03341 }
03342
03343
03344
03345 STATIC void GrainChrgTransferRates(long nd)
03346 {
03347 long nz;
03348 double fac0 = STICK_ION*gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3;
03349
03350 DEBUG_ENTRY( "GrainChrgTransferRates()" );
03351
03352 # ifndef IGNORE_GRAIN_ION_COLLISIONS
03353
03354 for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
03355 {
03356 long ion;
03357 ChargeBin *gptr = gv.bin[nd]->chrg[nz];
03358 double fac1 = gptr->FracPop*fac0;
03359
03360 if( fac1 == 0. )
03361 continue;
03362
03363 for( ion=0; ion <= LIMELM; ion++ )
03364 {
03365 long nelem;
03366 double eta, fac2, xi;
03367
03368
03369
03370 GrainScreen(ion,nd,nz,&eta,&xi);
03371
03372 fac2 = eta*fac1;
03373
03374 if( fac2 == 0. )
03375 continue;
03376
03377 for( nelem=MAX2(0,ion-1); nelem < LIMELM; nelem++ )
03378 {
03379 if( dense.lgElmtOn[nelem] && ion != gptr->RecomZ0[nelem][ion] )
03380 {
03381 gv.GrainChTrRate[nelem][ion][gptr->RecomZ0[nelem][ion]] +=
03382 (float)(fac2*DoppVel.AveVel[nelem]);
03383 }
03384 }
03385 }
03386 }
03387 # endif
03388
03389 DEBUG_EXIT( "GrainChrgTransferRates()" );
03390 return;
03391 }
03392
03393
03394
03395
03396 STATIC void GrainUpdateRadius1(void)
03397 {
03398 long nelem,
03399 nd;
03400
03401 DEBUG_ENTRY( "GrainUpdateRadius1()" );
03402
03403 for( nelem=0; nelem < LIMELM; nelem++ )
03404 {
03405 gv.elmSumAbund[nelem] = 0.f;
03406 }
03407
03408
03409 for( nd=0; nd < gv.nBin; nd++ )
03410 {
03411 gv.bin[nd]->dstAbund = (float)(gv.bin[nd]->dstfactor*gv.GrainMetal*GrnVryDpth(nd));
03412 ASSERT( gv.bin[nd]->dstAbund > 0.f );
03413
03414
03415 gv.bin[nd]->cnv_H_pCM3 = dense.gas_phase[ipHYDROGEN]*gv.bin[nd]->dstAbund;
03416 gv.bin[nd]->cnv_CM3_pH = 1./gv.bin[nd]->cnv_H_pCM3;
03417
03418 gv.bin[nd]->cnv_CM3_pGR = gv.bin[nd]->cnv_H_pGR/gv.bin[nd]->cnv_H_pCM3;
03419 gv.bin[nd]->cnv_GR_pCM3 = 1./gv.bin[nd]->cnv_CM3_pGR;
03420
03421
03422
03423
03424 for( nelem=0; nelem < LIMELM; nelem++ )
03425 {
03426 gv.elmSumAbund[nelem] += gv.bin[nd]->elmAbund[nelem]*(float)gv.bin[nd]->cnv_H_pCM3;
03427 }
03428 }
03429
03430 DEBUG_EXIT( "GrainUpdateRadius1()" );
03431 return;
03432 }
03433
03434
03435
03436
03437 STATIC void GrainUpdateRadius2(bool lgAnyNegCharge)
03438 {
03439 bool lgChangeCS_PDT;
03440 long i,
03441 nd,
03442 nz;
03443
03444 DEBUG_ENTRY( "GrainUpdateRadius2()" );
03445
03446
03447 lgChangeCS_PDT = gv.lgAnyNegCharge || lgAnyNegCharge;
03448
03449 if( rfield.nflux > gv.nfill || ( gv.lgAnyDustVary && nzone != gv.nzone ) || lgChangeCS_PDT )
03450 {
03451
03452 gv.nfill = rfield.nflux;
03453 gv.lgAnyNegCharge = lgAnyNegCharge;
03454
03455 for( i=0; i < rfield.nupper; i++ )
03456 {
03457 gv.dstab[i] = 0.;
03458 gv.dstsc[i] = 0.;
03459 }
03460
03461
03462 for( nd=0; nd < gv.nBin; nd++ )
03463 {
03464
03465 for( i=0; i < rfield.nflux; i++ )
03466 {
03467
03468
03469
03470
03471
03472
03473 gv.dstab[i] += gv.bin[nd]->dstab1[i]*gv.bin[nd]->dstAbund;
03474 gv.dstsc[i] += gv.bin[nd]->pure_sc1[i]*gv.bin[nd]->asym[i]*gv.bin[nd]->dstAbund;
03475
03476
03477
03478
03479
03480 for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
03481 {
03482 ChargeBin *gptr = gv.bin[nd]->chrg[nz];
03483 long ipLo = gptr->ipThresInf;
03484
03485
03486 if( gptr->DustZ <= -1 && i >= ipLo )
03487 gv.dstab[i] +=
03488 gptr->FracPop*gptr->cs_pdt[i-ipLo]*gv.bin[nd]->dstAbund;
03489 }
03490 }
03491 }
03492 for( i=0; i < rfield.nflux; i++ )
03493 {
03494
03495 ASSERT( gv.dstab[i] > 0. && gv.dstsc[i] > 0. );
03496 }
03497 }
03498
03499 DEBUG_EXIT( "GrainUpdateRadius2()" );
03500 return;
03501 }
03502
03503
03504
03505 STATIC void GrainTemperature(long int nd,
03506 float *dccool,
03507 double *hcon,
03508 double *hots,
03509 double *hla)
03510 {
03511 long int i,
03512 ipLya,
03513 nz;
03514 double EhatThermionic,
03515 norm,
03516 rate,
03517 x,
03518 y;
03519 float dcheat;
03520
03521 DEBUG_ENTRY( "GrainTemperature()" );
03522
03523
03524 ASSERT( nd >= 0 && nd < gv.nBin );
03525
03526 if( trace.lgTrace && trace.lgDustBug )
03527 {
03528 fprintf( ioQQQ, " GrainTemperature starts for grain %s\n", gv.bin[nd]->chDstLab );
03529 }
03530
03531
03532
03533
03534
03535 *hcon = 0.;
03536
03537 *hots = 0.;
03538
03539 *hla = 0.;
03540
03541 ipLya = EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].ipCont - 1;
03542
03543
03544
03545 gv.bin[nd]->GasHeatPhotoEl = 0.;
03546
03547 gv.bin[nd]->GrainCoolTherm = 0.;
03548 gv.bin[nd]->thermionic = 0.;
03549
03550 dcheat = 0.f;
03551 *dccool = 0.f;
03552
03553 gv.bin[nd]->BolFlux = 0.;
03554
03555
03556
03557 for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
03558 {
03559 ChargeBin *gptr = gv.bin[nd]->chrg[nz];
03560 long ipLo = gptr->ipThresInf;
03561
03562 double hcon1 = 0.;
03563 double hots1 = 0.;
03564 double hla1 = 0.;
03565 double bolflux1 = 0.;
03566 double pe1 = 0.;
03567
03568
03569 bool lgReEvaluate1 = gptr->hcon1 < 0.;
03570 bool lgReEvaluate2 = gptr->hots1 < 0.;
03571 bool lgReEvaluate = lgReEvaluate1 || lgReEvaluate2;
03572
03573
03574 if( lgReEvaluate )
03575 {
03576 long loopmax = MIN2(gptr->ipThresInf,rfield.nflux);
03577 for( i=0; i < loopmax; i++ )
03578 {
03579 double fac = gv.bin[nd]->dstab1[i]*rfield.anu[i];
03580
03581 if( lgReEvaluate1 )
03582 hcon1 += rfield.flux[i]*fac;
03583
03584 if( lgReEvaluate2 )
03585 {
03586 hots1 += rfield.SummedDif[i]*fac;
03587 # ifndef NDEBUG
03588 bolflux1 += rfield.SummedCon[i]*fac;
03589 # endif
03590 }
03591 }
03592 }
03593
03594
03595
03596
03597
03598
03599
03600 if( lgReEvaluate1 )
03601 {
03602 for( i=gptr->ipThresInf; i < rfield.nflux; i++ )
03603 {
03604 hcon1 += rfield.flux[i]*gptr->fac1[i-ipLo];
03605 }
03606
03607 gptr->hcon1 = hcon1;
03608 }
03609 else
03610 {
03611 hcon1 = gptr->hcon1;
03612 }
03613
03614 if( lgReEvaluate2 )
03615 {
03616 for( i=gptr->ipThresInf; i < rfield.nflux; i++ )
03617 {
03618
03619
03620 hots1 += rfield.SummedDif[i]*gptr->fac1[i-ipLo];
03621
03622 #ifdef WD_TEST2
03623 pe1 += rfield.flux[i]*gptr->fac2[i-ipLo];
03624 #else
03625 pe1 += rfield.SummedCon[i]*gptr->fac2[i-ipLo];
03626 #endif
03627 # ifndef NDEBUG
03628 bolflux1 += rfield.SummedCon[i]*gv.bin[nd]->dstab1[i]*rfield.anu[i];
03629 if( gptr->DustZ <= -1 )
03630 bolflux1 += rfield.SummedCon[i]*gptr->cs_pdt[i-ipLo]*rfield.anu[i];
03631 # endif
03632 }
03633 gptr->hots1 = hots1;
03634 gptr->bolflux1 = bolflux1;
03635 gptr->pe1 = pe1;
03636 }
03637 else
03638 {
03639 hots1 = gptr->hots1;
03640 bolflux1 = gptr->bolflux1;
03641 pe1 = gptr->pe1;
03642 }
03643
03644
03645
03646
03647
03648 if( ipLya < MIN2(gptr->ipThresInf,rfield.nflux) )
03649 {
03650 hla1 = rfield.otslin[ipLya]*gv.bin[nd]->dstab1[ipLya]*0.75;
03651 }
03652 else if( ipLya < rfield.nflux )
03653 {
03654
03655 hla1 = rfield.otslin[ipLya]*gptr->fac1[ipLya-ipLo];
03656 }
03657 else
03658 {
03659 hla1 = 0.;
03660 }
03661
03662 ASSERT( hcon1 >= 0. && hots1 >= 0. && hla1 >= 0. && bolflux1 >= 0. && pe1 >= 0. );
03663
03664 *hcon += gptr->FracPop*hcon1;
03665 *hots += gptr->FracPop*hots1;
03666 *hla += gptr->FracPop*hla1;
03667 gv.bin[nd]->BolFlux += gptr->FracPop*bolflux1;
03668 if( gv.lgDHetOn )
03669 gv.bin[nd]->GasHeatPhotoEl += gptr->FracPop*pe1;
03670
03671 # ifndef NDEBUG
03672 if( trace.lgTrace && trace.lgDustBug )
03673 {
03674 fprintf( ioQQQ, " Zg %ld bolflux: %.4e\n", gptr->DustZ,
03675 gptr->FracPop*bolflux1*EN1RYD*gv.bin[nd]->cnv_H_pCM3 );
03676 }
03677 # endif
03678
03679
03680
03681
03682
03683
03684
03685 rate = gptr->FracPop*gptr->ThermRate*gv.bin[nd]->IntArea*gv.bin[nd]->cnv_H_pCM3;
03686
03687 EhatThermionic = 2.*BOLTZMANN*gv.bin[nd]->tedust + MAX2(gptr->PotSurf*EN1RYD,0.);
03688 gv.bin[nd]->GrainCoolTherm += rate * (EhatThermionic + gptr->ThresSurf*EN1RYD);
03689 gv.bin[nd]->thermionic += rate * (EhatThermionic - gptr->PotSurf*EN1RYD);
03690 }
03691
03692
03693 norm = EN1RYD*gv.bin[nd]->cnv_H_pCM3;
03694
03695
03696 *hcon *= norm;
03697
03698
03699 *hots *= norm;
03700
03701
03702 *hla *= norm;
03703
03704 gv.bin[nd]->BolFlux *= norm;
03705
03706
03707
03708
03709
03710
03711
03712 GrainCollHeating(nd,&dcheat,dccool);
03713
03714
03715 gv.bin[nd]->GasHeatPhotoEl *= norm;
03716
03717 if( gv.lgBakesPAH_heat )
03718 {
03719
03720
03721
03722
03723
03724 gv.bin[nd]->GasHeatPhotoEl = 1.e-24*hmi.UV_Cont_rel2_Habing_TH85_depth*
03725 dense.gas_phase[ipHYDROGEN]*(4.87e-2/(1.0+4e-3*pow((hmi.UV_Cont_rel2_Habing_TH85_depth*
03726
03727 phycon.sqrte/dense.eden),0.73)) + 3.65e-2*pow(phycon.te/1.e4,0.7)/
03728 (1.+2.e-4*(hmi.UV_Cont_rel2_Habing_TH85_depth*phycon.sqrte/dense.eden)))/gv.nBin;
03729
03730 }
03731
03732
03733
03734 gv.bin[nd]->GasHeatPhotoEl *= gv.GrainHeatScaleFactor;
03735
03736
03737
03738
03739
03740
03741
03742
03743
03744
03745
03746
03747
03748
03749
03750
03751 gv.bin[nd]->GrainHeat = *hcon + *hots + dcheat - gv.bin[nd]->GrainCoolTherm;
03752
03753
03754 gv.bin[nd]->GrainHeatColl = dcheat;
03755
03756
03757
03758
03759
03760 if( gv.bin[nd]->GrainHeat > 0. )
03761 {
03762 bool lgOutOfBounds;
03763
03764
03765 x = log(MAX2(DBL_MIN,gv.bin[nd]->GrainHeat*gv.bin[nd]->cnv_CM3_pH));
03766
03767 splint_safe(gv.bin[nd]->dstems,gv.dsttmp,gv.bin[nd]->dstslp,NDEMS,x,&y,&lgOutOfBounds);
03768 gv.bin[nd]->tedust = (float)exp(y);
03769 }
03770 else
03771 {
03772 gv.bin[nd]->GrainHeat = -1.;
03773 gv.bin[nd]->tedust = -1.;
03774 }
03775
03776 if( thermal.ConstGrainTemp > 0. )
03777 {
03778 bool lgOutOfBounds;
03779
03780 gv.bin[nd]->tedust = thermal.ConstGrainTemp;
03781
03782 x = log(gv.bin[nd]->tedust);
03783 splint_safe(gv.dsttmp,gv.bin[nd]->dstems,gv.bin[nd]->dstslp2,NDEMS,x,&y,&lgOutOfBounds);
03784 gv.bin[nd]->GrainHeat = exp(y)*gv.bin[nd]->cnv_H_pCM3;
03785 }
03786
03787
03788 gv.bin[nd]->TeGrainMax = (float)MAX2(gv.bin[nd]->TeGrainMax,gv.bin[nd]->tedust);
03789
03790 if( trace.lgTrace && trace.lgDustBug )
03791 {
03792 fprintf( ioQQQ, " >GrainTemperature finds %s Tdst %.5e hcon %.4e ",
03793 gv.bin[nd]->chDstLab, gv.bin[nd]->tedust, *hcon);
03794 fprintf( ioQQQ, "hots %.4e dcheat %.4e GrainCoolTherm %.4e\n",
03795 *hots, dcheat, gv.bin[nd]->GrainCoolTherm );
03796 }
03797
03798 DEBUG_EXIT( "GrainTemperature()" );
03799 return;
03800 }
03801
03802
03803
03804 void PE_init(long nd,
03805 long nz,
03806 long i,
03807 double *cs1,
03808 double *cs2,
03809 double *cs_tot,
03810 double *cool1,
03811 double *cool2,
03812 double *ehat1,
03813 double *ehat2)
03814 {
03815 ChargeBin *gptr = gv.bin[nd]->chrg[nz];
03816 long ipLo1 = gptr->ipThresInfVal;
03817 long ipLo2 = gptr->ipThresInf;
03818
03819 DEBUG_ENTRY( "PE_init()" );
03820
03821
03822 ASSERT( nd >= 0 && nd < gv.nBin );
03823 ASSERT( nz >= 0 && nz < gv.bin[nd]->nChrg );
03824 ASSERT( i >= 0 && i < rfield.nflux );
03825
03826
03827 *cs1 = ( i >= ipLo1 ) ? gv.bin[nd]->dstab1[i]*gptr->yhat[i-ipLo1] : 0.;
03828 *cs2 = ( gptr->DustZ <= -1 && i >= ipLo2 ) ? gptr->cs_pdt[i-ipLo2] : 0.;
03829 *cs_tot = gv.bin[nd]->dstab1[i] + *cs2;
03830
03831
03832
03833 if( gptr->DustZ > -1 )
03834 {
03835
03836 double Elo = -gptr->PotSurf;
03837 double Ehi = rfield.anu[i] - gptr->ThresInfVal;
03838 *ehat1 = ( Ehi > 0. ) ?
03839 0.5*Ehi*(Ehi-2.*Elo)/(Ehi-3.*Elo) + gptr->PotSurf : 0.;
03840 *ehat2 = 0.;
03841 }
03842 else
03843 {
03844
03845 *ehat1 = 0.5*MAX2(rfield.anu[i] - gptr->ThresSurfVal,0.);
03846
03847 *ehat2 = MAX2(rfield.anu[i] - gptr->ThresInf,0.);
03848 }
03849
03850
03851 *cool1 = gptr->ThresSurfVal + *ehat1;
03852
03853
03854
03855
03856 if( gptr->DustZ <= -1 )
03857 *cool1 += gptr->ThresSurf - gptr->ThresSurfVal;
03858
03859 *cool2 = gptr->ThresSurf + *ehat2;
03860
03861
03862 ASSERT( *ehat1 > 0. || ( *ehat1 == 0. && *cs1 == 0. ) );
03863 ASSERT( *ehat2 > 0. || ( *ehat2 == 0. && *cs2 == 0. ) );
03864 ASSERT( *cool1 >= 0. && *cool2 >= 0. );
03865
03866 DEBUG_EXIT( "PE_init()" );
03867 return;
03868 }
03869
03870
03871
03872 STATIC void GrainCollHeating(long int nd,
03873 float *dcheat,
03874 float *dccool)
03875 {
03876 long int ion,
03877 nelem,
03878 nz;
03879 H2_type ipH2;
03880 double Accommodation,
03881 CollisionRateElectr,
03882 CollisionRateMol,
03883 CollisionRateIon,
03884 CoolTot,
03885 CoolBounce,
03886 CoolEmitted,
03887 CoolElectrons,
03888 CoolMolecules,
03889 CoolPotential,
03890 CoolPotentialGas,
03891 eta,
03892 HeatTot,
03893 HeatBounce,
03894 HeatCollisions,
03895 HeatElectrons,
03896 HeatIons,
03897 HeatMolecules,
03898 HeatRecombination,
03899 HeatChem,
03900 HeatCor,
03901 Stick,
03902 ve,
03903 WeightMol,
03904 xi;
03905
03906
03907
03908 const double H2_FORMATION_GRAIN_HEATING[H2_TOP] = { 0.20, 0.4, 1.72 };
03909
03910 DEBUG_ENTRY( "GrainCollHeating()" );
03911
03912
03913
03914
03915
03916
03917
03918
03919
03920
03921
03922
03923 HeatTot = 0.;
03924 CoolTot = 0.;
03925
03926 HeatIons = 0.;
03927
03928 gv.bin[nd]->ChemEn = 0.;
03929
03930
03931 for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
03932 {
03933 ChargeBin *gptr = gv.bin[nd]->chrg[nz];
03934
03935
03936
03937 double Heat1 = 0.;
03938 double Cool1 = 0.;
03939 double ChemEn1 = 0.;
03940
03941
03942
03943
03944
03945 for( ion=0; ion <= LIMELM; ion++ )
03946 {
03947
03948
03949
03950
03951
03952
03953
03954 CollisionRateIon = 0.;
03955 CoolPotential = 0.;
03956 CoolPotentialGas = 0.;
03957 HeatRecombination = 0.;
03958 HeatChem = 0.;
03959
03960
03961
03962 GrainScreen(ion,nd,nz,&eta,&xi);
03963
03964 for( nelem=MAX2(0,ion-1); nelem < LIMELM; nelem++ )
03965 {
03966 if( dense.lgElmtOn[nelem] && dense.xIonDense[nelem][ion] > 0. )
03967 {
03968 double CollisionRateOne;
03969
03970
03971
03972
03973 #if defined( IGNORE_GRAIN_ION_COLLISIONS )
03974 Stick = 0.;
03975 #elif defined( WD_TEST2 )
03976 Stick = ( ion == gptr->RecomZ0[nelem][ion] ) ?
03977 0. : STICK_ION;
03978 #else
03979 Stick = ( ion == gptr->RecomZ0[nelem][ion] ) ?
03980 gv.bin[nd]->AccomCoef[nelem] : STICK_ION;
03981 #endif
03982
03983
03984
03985 CollisionRateOne = Stick*dense.xIonDense[nelem][ion]*DoppVel.AveVel[nelem];
03986 CollisionRateIon += CollisionRateOne;
03987
03988
03989
03990
03991
03992
03993
03994
03995 if( ion >= gptr->RecomZ0[nelem][ion] )
03996 {
03997 CoolPotential += CollisionRateOne * (double)ion *
03998 gptr->PotSurf;
03999 CoolPotentialGas += CollisionRateOne *
04000 (double)gptr->RecomZ0[nelem][ion] *
04001 gptr->PotSurf;
04002 }
04003 else
04004 {
04005 CoolPotential += CollisionRateOne * (double)ion *
04006 gptr->PotSurfInc;
04007 CoolPotentialGas += CollisionRateOne *
04008 (double)gptr->RecomZ0[nelem][ion] *
04009 gptr->PotSurfInc;
04010 }
04011
04012
04013
04014
04015 HeatRecombination += CollisionRateOne *
04016 gptr->RecomEn[nelem][ion];
04017 HeatChem += CollisionRateOne * gptr->ChemEn[nelem][ion];
04018 }
04019 }
04020
04021
04022
04023
04024
04025
04026 HeatCollisions = CollisionRateIon * 2.*BOLTZMANN*phycon.te*xi;
04027
04028
04029
04030
04031 CoolPotential *= eta*EN1RYD;
04032 CoolPotentialGas *= eta*EN1RYD;
04033
04034 HeatRecombination *= eta*EN1RYD;
04035 HeatChem *= eta*EN1RYD;
04036
04037 CoolEmitted = CollisionRateIon * 2.*BOLTZMANN*gv.bin[nd]->tedust*eta;
04038
04039
04040 Heat1 += HeatCollisions - CoolPotential + HeatRecombination - CoolEmitted;
04041
04042
04043
04044
04045 Cool1 += HeatCollisions - CoolEmitted - CoolPotentialGas;
04046
04047 ChemEn1 += HeatChem;
04048 }
04049
04050
04051 HeatIons += gptr->FracPop*Heat1;
04052
04053 if( trace.lgTrace && trace.lgDustBug )
04054 {
04055 fprintf( ioQQQ, " Zg %ld ions heat/cool: %.4e %.4e\n", gptr->DustZ,
04056 gptr->FracPop*Heat1*gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3,
04057 gptr->FracPop*Cool1*gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3 );
04058 }
04059
04060
04061
04062
04063 ion = -1;
04064 Stick = ( gptr->DustZ <= -1 ) ? gv.bin[nd]->StickElecNeg : gv.bin[nd]->StickElecPos;
04065
04066
04067 ve = sqrt(8.*BOLTZMANN/PI/ELECTRON_MASS*phycon.te);
04068
04069
04070 CollisionRateElectr = Stick*dense.eden*ve;
04071
04072
04073
04074 GrainScreen(ion,nd,nz,&eta,&xi);
04075
04076 if( gptr->DustZ > gv.bin[nd]->LowestZg )
04077 {
04078 HeatCollisions = CollisionRateElectr*2.*BOLTZMANN*phycon.te*xi;
04079
04080
04081 CoolPotential = CollisionRateElectr * (double)ion*gptr->PotSurfInc*eta*EN1RYD;
04082
04083 HeatRecombination = CollisionRateElectr * gptr->ThresSurfInc*eta*EN1RYD;
04084 HeatBounce = 0.;
04085 CoolBounce = 0.;
04086 }
04087 else
04088 {
04089 HeatCollisions = 0.;
04090 CoolPotential = 0.;
04091 HeatRecombination = 0.;
04092
04093
04094
04095
04096 HeatBounce = CollisionRateElectr * 2.*BOLTZMANN*phycon.te*xi;
04097
04098
04099 CoolBounce = CollisionRateElectr *
04100 (-gptr->ThresSurfInc-gptr->PotSurfInc)*EN1RYD*eta;
04101 CoolBounce = MAX2(CoolBounce,0.);
04102 }
04103
04104
04105
04106 HeatElectrons = HeatCollisions-CoolPotential+HeatRecombination+HeatBounce-CoolBounce;
04107 Heat1 += HeatElectrons;
04108
04109 CoolElectrons = HeatCollisions+HeatBounce-CoolBounce;
04110 Cool1 += CoolElectrons;
04111
04112 if( trace.lgTrace && trace.lgDustBug )
04113 {
04114 fprintf( ioQQQ, " Zg %ld electrons heat/cool: %.4e %.4e\n", gptr->DustZ,
04115 gptr->FracPop*HeatElectrons*gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3,
04116 gptr->FracPop*CoolElectrons*gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3 );
04117 }
04118
04119
04120
04121
04122
04123
04124
04125
04126
04127
04128
04129
04130 gptr->HeatingRate2 = HeatElectrons*gv.bin[nd]->IntArea/4. -
04131 gv.bin[nd]->GrainCoolTherm*gv.bin[nd]->cnv_CM3_pH;
04132
04133
04134
04135
04136
04137 HeatTot += gptr->FracPop*Heat1;
04138
04139
04140 CoolTot += gptr->FracPop*Cool1;
04141
04142 gv.bin[nd]->ChemEn += gptr->FracPop*ChemEn1;
04143 }
04144
04145
04146
04147
04148
04149
04150
04151
04152 WeightMol = 2.*dense.AtomicWeight[ipHYDROGEN];
04153 Accommodation = 2.*gv.bin[nd]->atomWeight*WeightMol/POW2(gv.bin[nd]->atomWeight+WeightMol);
04154
04155 #ifndef IGNORE_GRAIN_ION_COLLISIONS
04156
04157 CollisionRateMol = Accommodation*hmi.H2_total*
04158 sqrt(8.*BOLTZMANN/PI/ATOMIC_MASS_UNIT/WeightMol*phycon.te);
04159
04160
04161 ipH2 = gv.which_H2distr[gv.bin[nd]->matType];
04162
04163
04164 gv.bin[nd]->ChemEnH2 = gv.bin[nd]->rate_h2_form_grains_used*dense.xIonDense[ipHYDROGEN][0]*
04165 H2_FORMATION_GRAIN_HEATING[ipH2]*EN1EV;
04166
04167 gv.bin[nd]->ChemEnH2 /= gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3;
04168 #else
04169 CollisionRateMol = 0.;
04170 gv.bin[nd]->ChemEnH2 = 0.;
04171 #endif
04172
04173
04174 WeightMol = dense.AtomicWeight[ipCARBON] + dense.AtomicWeight[ipOXYGEN];
04175 Accommodation = 2.*gv.bin[nd]->atomWeight*WeightMol/POW2(gv.bin[nd]->atomWeight+WeightMol);
04176 #ifndef IGNORE_GRAIN_ION_COLLISIONS
04177 CollisionRateMol += Accommodation*findspecies("CO")->hevmol*
04178 sqrt(8.*BOLTZMANN/PI/ATOMIC_MASS_UNIT/WeightMol*phycon.te);
04179 #else
04180 CollisionRateMol = 0.;
04181 #endif
04182
04183
04184 HeatCollisions = CollisionRateMol * 2.*BOLTZMANN*phycon.te;
04185 CoolEmitted = CollisionRateMol * 2.*BOLTZMANN*gv.bin[nd]->tedust;
04186
04187 HeatMolecules = HeatCollisions - CoolEmitted + gv.bin[nd]->ChemEnH2;
04188 HeatTot += HeatMolecules;
04189
04190
04191 CoolMolecules = HeatCollisions - CoolEmitted;
04192 CoolTot += CoolMolecules;
04193
04194 gv.bin[nd]->RateUp = 0.;
04195 gv.bin[nd]->RateDn = 0.;
04196 HeatCor = 0.;
04197 for( nz=0; nz < gv.bin[nd]->nChrg; nz++ )
04198 {
04199 double d[4],Auger;
04200 double rate_dn = GrainElecRecomb1(nd,nz,&d[0],&d[1]);
04201 double rate_up = GrainElecEmis1(nd,nz,&d[0],&Auger,&d[1],&d[2],&d[3]);
04202
04203 gv.bin[nd]->RateUp += gv.bin[nd]->chrg[nz]->FracPop*rate_up;
04204 gv.bin[nd]->RateDn += gv.bin[nd]->chrg[nz]->FracPop*rate_dn;
04205
04206
04207
04210 HeatCor += (gv.bin[nd]->chrg[nz]->FracPop*(rate_up-Auger)*gv.bin[nd]->chrg[nz]->ThresSurf -
04211 gv.bin[nd]->chrg[nz]->FracPop*rate_dn*gv.bin[nd]->chrg[nz]->ThresSurfInc +
04212 gv.bin[nd]->chrg[nz]->FracPop*(rate_up-Auger)*gv.bin[nd]->chrg[nz]->PotSurf -
04213 gv.bin[nd]->chrg[nz]->FracPop*rate_dn*gv.bin[nd]->chrg[nz]->PotSurfInc)*EN1RYD;
04214 }
04215
04216
04217 HeatTot += HeatCor;
04218
04219 if( trace.lgTrace && trace.lgDustBug )
04220 {
04221 fprintf( ioQQQ, " molecules heat/cool: %.4e %.4e heatcor: %.4e\n",
04222 HeatMolecules*gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3,
04223 CoolMolecules*gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3,
04224 HeatCor*gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3 );
04225 }
04226
04227 *dcheat = (float)(HeatTot*gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3);
04228 *dccool = ( gv.lgDColOn ) ? (float)(CoolTot*gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3) : 0.f;
04229
04230 gv.bin[nd]->ChemEn *= gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3;
04231 gv.bin[nd]->ChemEnH2 *= gv.bin[nd]->IntArea/4.*gv.bin[nd]->cnv_H_pCM3;
04232
04233
04234
04235
04236
04237
04238
04239
04240
04241
04242
04243
04244
04245
04246
04247
04248 gv.bin[nd]->HeatingRate1 = (HeatMolecules+HeatIons+HeatCor)*gv.bin[nd]->IntArea/4.;
04249
04250
04251
04252 DEBUG_EXIT( "GrainCollHeating()" );
04253 return;
04254 }
04255
04256
04257
04258 void GrainDrift(void)
04259 {
04260 long int i,
04261 loop,
04262 nd;
04263 float *help;
04264 double alam,
04265 corr,
04266 dmomen,
04267 fac,
04268 fdrag,
04269 g0,
04270 g2,
04271 phi2lm,
04272 psi,
04273 rdust,
04274 si,
04275 vdold,
04276 volmom;
04277
04278 DEBUG_ENTRY( "GrainDrift()" );
04279
04280
04281 help = (float*)MALLOC((size_t)((unsigned)rfield.nflux*sizeof(float)));
04282 for( i=0; i < rfield.nflux; i++ )
04283 {
04284 help[i] = (rfield.flux[i]+rfield.ConInterOut[i]+rfield.outlin[i]+rfield.outlin_noplot[i])*
04285 rfield.anu[i];
04286 }
04287
04288 for( nd=0; nd < gv.nBin; nd++ )
04289 {
04290
04291 dmomen = 0.;
04292 for( i=0; i < rfield.nflux; i++ )
04293 {
04294
04295 dmomen += help[i]*(gv.bin[nd]->dstab1[i] + gv.bin[nd]->pure_sc1[i]*gv.bin[nd]->asym[i]);
04296 }
04297 ASSERT( dmomen >= 0. );
04298 dmomen *= EN1RYD*4./gv.bin[nd]->IntArea;
04299
04300
04301 fac = 2*BOLTZMANN*phycon.te;
04302
04303
04304
04305 psi = gv.bin[nd]->dstpot*TE1RYD/phycon.te;
04306 if( psi > 0. )
04307 {
04308 rdust = 1.e-6;
04309 alam = log(20.702/rdust/psi*phycon.sqrte/dense.SqrtEden);
04310 }
04311 else
04312 {
04313 alam = 0.;
04314 }
04315
04316 phi2lm = POW2(psi)*alam;
04317 corr = 2.;
04318
04319 for( loop = 0; loop < 50 && fabs(corr-1.) > 0.001; loop++ )
04320 {
04321 vdold = gv.bin[nd]->DustDftVel;
04322
04323
04324 si = gv.bin[nd]->DustDftVel/phycon.sqrte*7.755e-5;
04325 g0 = 1.5045*si*sqrt(1.+0.4418*si*si);
04326 g2 = si/(1.329 + POW3(si));
04327
04328
04329
04330 fdrag = fac*dense.xIonDense[ipHYDROGEN][1]*(g0 + phi2lm*g2);
04331
04332
04333 si = gv.bin[nd]->DustDftVel/phycon.sqrte*1.816e-6;
04334 g0 = 1.5045*si*sqrt(1.+0.4418*si*si);
04335 g2 = si/(1.329 + POW3(si));
04336 fdrag += fac*dense.eden*(g0 + phi2lm*g2);
04337
04338
04339 si = gv.bin[nd]->DustDftVel/phycon.sqrte*7.755e-5;
04340 g0 = 1.5045*si*sqrt(1.+0.4418*si*si);
04341 fdrag += fac*(dense.xIonDense[ipHYDROGEN][0] + 1.1*dense.xIonDense[ipHELIUM][0])*g0;
04342
04343
04344 si = gv.bin[nd]->DustDftVel/phycon.sqrte*1.551e-4;
04345 g0 = 1.5045*si*sqrt(1.+0.4418*si*si);
04346 g2 = si/(1.329 + POW3(si));
04347 fdrag += fac*dense.xIonDense[ipHELIUM][1]*(g0 + phi2lm*g2);
04348
04349
04350
04351
04352 volmom = dmomen/SPEEDLIGHT;
04353
04354 if( fdrag > 0. )
04355 {
04356 corr = sqrt(volmom/fdrag);
04357 gv.bin[nd]->DustDftVel = (float)(vdold*corr);
04358 }
04359 else
04360 {
04361 corr = 1.;
04362 gv.lgNegGrnDrg = true;
04363 gv.bin[nd]->DustDftVel = 0.;
04364 }
04365
04366 if( trace.lgTrace && trace.lgDustBug )
04367 {
04368 fprintf( ioQQQ, " %2ld new drift velocity:%10.2e momentum absorbed:%10.2e\n",
04369 loop, gv.bin[nd]->DustDftVel, volmom );
04370 }
04371 }
04372 }
04373
04374 free( help );
04375
04376 DEBUG_EXIT( "GrainDrift()" );
04377 return;
04378 }
04379
04380
04381 STATIC double GrnVryDpth(
04382
04383
04384
04385
04386
04387 long int nd)
04388 {
04389 double GrnVryDpth_v;
04390
04391 DEBUG_ENTRY( "GrnVryDpth()" );
04392
04393
04394
04395
04396
04397
04398 if( gv.bin[nd]->lgDustFunc )
04399 {
04400
04401
04402
04403
04404
04405
04406
04407 GrnVryDpth_v = dense.xIonDense[ipHYDROGEN][0]/dense.gas_phase[ipHYDROGEN];
04408 }
04409 else
04410 {
04411
04412
04413
04414 GrnVryDpth_v = GrnStdDpth(nd);
04415 }
04416
04417 ASSERT( GrnVryDpth_v > 0. && GrnVryDpth_v <= 1.000001 );
04418
04419 DEBUG_EXIT( "GrnVryDpth()" );
04420 return GrnVryDpth_v;
04421 }
04422
04423