00001
00002
00003
00004
00005
00014 #include "cddefines.h"
00015 #include "lines_service.h"
00016 #include "iso.h"
00017 #include "geometry.h"
00018 #include "h2.h"
00019 #include "mole.h"
00020 #include "hyperfine.h"
00021 #include "opacity.h"
00022 #include "dense.h"
00023 #include "heavy.h"
00024 #include "grainvar.h"
00025 #include "elementnames.h"
00026 #include "rfield.h"
00027 #include "dynamics.h"
00028 #include "thermal.h"
00029 #include "hmi.h"
00030 #include "coolheavy.h"
00031 #include "timesc.h"
00032 #include "doppvel.h"
00033 #include "stopcalc.h"
00034 #include "colden.h"
00035 #include "phycon.h"
00036 #include "rt.h"
00037 #include "trace.h"
00038 #include "wind.h"
00039 #include "save.h"
00040 #include "taulines.h"
00041 #include "pressure.h"
00042 #include "iterations.h"
00043 #include "struc.h"
00044 #include "radius.h"
00045 #include "dark_matter.h"
00046
00047
00048 STATIC void ContRate(double *freqm,
00049 double *opacm);
00050
00051
00052 STATIC void GrainRateDr(double *grfreqm,
00053 double *gropacm);
00054
00055
00056
00057 int radius_next(void)
00058 {
00059 static double OHIIoHI,
00060 OldHeat = 0.,
00061 OldTe = 0.,
00062 OlddTdStep = 0.,
00063 OldWindVelocity = 0.,
00064 Old_H2_heat_cool;
00065
00066 const double BigRadius = 1e30;
00067
00068 DEBUG_ENTRY( "radius_next()" );
00069
00070
00071
00072
00073
00074
00075
00076
00077 if( trace.lgTrace )
00078 fprintf( ioQQQ, " radius_next called\n" );
00079
00080
00081 bool lgFirstCall = ( nzone == 0 );
00082
00083 multimap<double,string> drChoice;
00084
00085
00086
00087 if( !lgFirstCall && OHIIoHI > 1e-15 &&
00088 (dense.xIonDense[ipHYDROGEN][1] > 0. && dense.xIonDense[ipHYDROGEN][0] > 0.) )
00089 {
00090 double atomic_frac = dense.xIonDense[ipHYDROGEN][1]/dense.xIonDense[ipHYDROGEN][0];
00091 double drHion;
00092
00093
00094 double x = 1. - atomic_frac/OHIIoHI;
00095 if( atomic_frac > 0.05 && atomic_frac < 0.9 )
00096 {
00097
00098
00099
00100 drHion = radius.drad*MAX2( 0.2 , 0.05/MAX2(1e-10,x) );
00101 }
00102 else if( x > 0. )
00103 {
00104
00105
00106 drHion = radius.drad*MAX2( 0.2 , 0.2/MAX2(1e-10,x) );
00107 }
00108 else
00109 {
00110 drHion = BigRadius;
00111 }
00112 double SaveOHIIoHI = OHIIoHI;
00113 OHIIoHI = dense.xIonDense[ipHYDROGEN][1]/dense.xIonDense[ipHYDROGEN][0];
00114 ostringstream oss;
00115 oss << "change in H ionization old=" << scientific << setprecision(3);
00116 oss << SaveOHIIoHI << " new=" << OHIIoHI;
00117 drChoice.insert( pair<const double,string>( drHion, oss.str() ) );
00118 }
00119 else
00120 {
00121 if( dense.xIonDense[ipHYDROGEN][1] > 0. && dense.xIonDense[ipHYDROGEN][0] > 0. )
00122 OHIIoHI = dense.xIonDense[ipHYDROGEN][1]/dense.xIonDense[ipHYDROGEN][0];
00123 else
00124 OHIIoHI = 0.;
00125 }
00126
00127 if( rt.dTauMase < -0.01 )
00128 {
00129
00130
00131 double drHMase = radius.drad*MAX2(0.1,-0.2/rt.dTauMase);
00132 ostringstream oss;
00133
00134 oss << "H maser dTauMase=" << scientific << setprecision(2) << rt.dTauMase << " ";
00135 oss << rt.mas_species << " " << rt.mas_ion << " " << rt.mas_hi << " " << rt.mas_lo;
00136 drChoice.insert( pair<const double,string>( drHMase, oss.str() ) );
00137 }
00138
00139
00140
00141
00142 double change_heavy_frac_big = -1.;
00143 double frac_change_heavy_frac_big = -1.;
00144 const realnum CHANGE_ION_HEAV = 0.2f;
00145 const realnum CHANGE_ION_HHE = 0.15f;
00146 if( nzone > 4 )
00147 {
00148 long ichange_heavy_nelem = -1L;
00149 long ichange_heavy_ion = -1L;
00150 double dr_change_heavy = BigRadius;
00151 for( int nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
00152 {
00153 if( dense.lgElmtOn[nelem] )
00154 {
00155 realnum change;
00156
00157
00158 realnum frac_limit;
00159 if( nelem<=ipHELIUM || nelem==ipCARBON )
00160 {
00161 frac_limit = 1e-4f;
00162 change = CHANGE_ION_HHE;
00163 }
00164 else
00165 {
00166
00167
00168
00169 frac_limit = struc.dr_ionfrac_limit/2.f;
00170 change = CHANGE_ION_HEAV;
00171 }
00172
00173 for( int ion=0; ion<=nelem+1; ++ion )
00174 {
00175 realnum abundnew = dense.xIonDense[nelem][ion]/SDIV( dense.gas_phase[nelem]);
00176 realnum abundold = struc.xIonDense[nelem][ion][nzone-3]/SDIV( struc.gas_phase[nelem][nzone-3]);
00177 if( abundnew > frac_limit && abundold > frac_limit )
00178 {
00179
00180 realnum abundolder = struc.xIonDense[nelem][ion][nzone-4]/SDIV( struc.gas_phase[nelem][nzone-4]);
00181 realnum abundoldest = struc.xIonDense[nelem][ion][nzone-5]/SDIV( struc.gas_phase[nelem][nzone-5]);
00182
00183
00184 double change_heavy_frac = fabs(abundnew-abundold)/MIN2(abundold,abundnew);
00185
00186 if( (change_heavy_frac > change) && (change_heavy_frac > change_heavy_frac_big) &&
00187
00188
00189 ((abundnew-abundold)>0.) &&
00190 ((abundold-abundolder)>0.) &&
00191 ((abundolder-abundoldest)>0.) )
00192 {
00193 ichange_heavy_nelem = nelem;
00194 ichange_heavy_ion = ion;
00195 change_heavy_frac_big = change_heavy_frac;
00196 frac_change_heavy_frac_big = abundnew;
00197
00198
00199
00200
00201
00202
00203
00204 dr_change_heavy = radius.drad * MAX2(0.25, change / change_heavy_frac );
00205 }
00206 }
00207 }
00208 }
00209 }
00210
00211
00212
00213 ostringstream oss;
00214 if(ichange_heavy_nelem>=0)
00215 {
00216 oss << "change in ionization, element " << elementnames.chElementName[ichange_heavy_nelem];
00217 oss << " ion (C scale) " << ichange_heavy_ion << " rel change " << scientific << setprecision(2);
00218 oss << change_heavy_frac_big << " ion frac " << frac_change_heavy_frac_big;
00219 }
00220 else
00221 {
00222 oss << "none";
00223 dr_change_heavy = BigRadius;
00224 }
00225 drChoice.insert( pair<const double,string>( dr_change_heavy, oss.str() ) );
00226 }
00227
00228
00229
00230 if( !co.lgUMISTrates )
00231 {
00232
00233
00234
00235 double drLeiden_hack = MAX2( 0.02 , 0.05*rfield.extin_mag_V_point ) /
00236 SDIV( geometry.FillFac * rfield.opac_mag_V_point );
00237 drChoice.insert( pair<const double,string>( drLeiden_hack, "Leiden hack" ) );
00238 }
00239
00240
00241 static double extin_mag_V_point_old = -1.;
00242 if( nzone>1 )
00243 {
00244 const double extin_mag_V_limit = 20.;
00245 double dExtin = (rfield.extin_mag_V_point - extin_mag_V_point_old)/extin_mag_V_limit;
00246 if( dExtin > SMALLFLOAT )
00247 {
00248 double drExtintion = radius.drad / dExtin;
00249 ostringstream oss;
00250 oss << "change in V mag extinction; current=" << scientific << setprecision(3) <<
00251 rfield.extin_mag_V_point;
00252 oss << fixed << " delta=" << dExtin;
00253 drChoice.insert( pair<const double,string>( drExtintion, oss.str() ) );
00254 }
00255
00256 }
00257 extin_mag_V_point_old = rfield.extin_mag_V_point;
00258
00259
00260
00261 if( nzone > 1 && !thermal.lgTemperatureConstant &&
00262
00263 !( strcmp(dense.chDenseLaw,"SINE") == 0 && dense.lgDenFlucOn ) )
00264 {
00265 double dHdStep = fabs(thermal.htot-OldHeat)/thermal.htot;
00266 if( dHdStep > 0. )
00267 {
00268 double drdHdStep;
00269 if( dense.gas_phase[ipHYDROGEN] >= 1e13 )
00270 {
00271 drdHdStep = radius.drad*MAX2(0.8,0.075/dHdStep);
00272 }
00273 else if( dense.gas_phase[ipHYDROGEN] >= 1e11 )
00274 {
00275 drdHdStep = radius.drad*MAX2(0.8,0.100/dHdStep);
00276 }
00277 else
00278 {
00279
00280
00281
00282
00283 drdHdStep = radius.drad*MAX2(0.8,0.15/dHdStep);
00284 }
00285 ostringstream oss;
00286 oss << "change in heating; current=" << scientific << setprecision(3) << thermal.htot;
00287 oss << fixed << " delta=" << dHdStep;
00288 drChoice.insert( pair<const double,string>( drdHdStep, oss.str() ) );
00289 }
00290 }
00291 OldHeat = thermal.htot;
00292
00293
00294
00295 if( strcmp(dense.chDenseLaw,"CPRE") == 0 && pressure.lgContRadPresOn )
00296 {
00297 if( nzone > 2 && pressure.pinzon > 0. )
00298 {
00299
00300
00301
00302
00303
00304 double drConPres = 0.05*pressure.PresTotlCurr/
00305 (wind.AccelTotalOutward*dense.xMassDensity*geometry.FillFac);
00306 drChoice.insert( pair<const double,string>( drConPres, "change in con accel" ) );
00307 }
00308 }
00309
00310
00311
00312 if( strcmp(dense.chDenseLaw,"CPRE") == 0 && (dark.lgNFW_Set || pressure.gravity_symmetry>=0) )
00313 {
00314 if( fabs( pressure.RhoGravity ) > 0. )
00315 {
00316 double drGravPres = 0.05 * pressure.PresTotlCurr * radius.drad / fabs( pressure.RhoGravity );
00317 ASSERT( drGravPres > 0. );
00318 drChoice.insert( pair<const double,string>( drGravPres, "change in gravitational pressure" ) );
00319 }
00320 }
00321
00322
00323
00324 double dTdStep = FLT_MAX;
00325 if( nzone > 1 )
00326 {
00327
00328 dTdStep = (phycon.te-OldTe)/phycon.te;
00329
00330
00331
00332
00333
00334 double x = 0.03;
00335
00336
00337
00338 if( dTdStep*OlddTdStep > 0. )
00339 {
00340
00341
00342
00343
00344
00345 double absdt = fabs(dTdStep);
00346 double drdTdStep = radius.drad*MAX2(0.5,x/absdt);
00347 ostringstream oss;
00348 oss << "change in temperature; current=" << scientific << setprecision(3) << phycon.te;
00349 oss << fixed << " dT/T=" << dTdStep;
00350 drChoice.insert( pair<const double,string>( drdTdStep, oss.str() ) );
00351 }
00352 }
00353 OlddTdStep = dTdStep;
00354 OldTe = phycon.te;
00355
00356
00357
00358
00359
00360
00361 if( !thermal.lgTemperatureConstant && !dynamics.lgRecom )
00362 {
00363 double freqm = 0., opacm = 0.;
00364
00365 ContRate(&freqm,&opacm);
00366
00367
00368
00369
00370
00371
00372 double drInter = 0.3/MAX2(1e-30,opacm*geometry.FillFac*geometry.DirectionalCosin);
00373 ostringstream oss;
00374 oss << "cont inter nu=" << scientific << setprecision(2) << freqm << " opac=" << opacm;
00375 drChoice.insert( pair<const double,string>( drInter, oss.str() ) );
00376 }
00377
00378
00379
00380
00381
00382
00383 if( !wind.lgStatic() && !pressure.lgSonicPoint && !pressure.lgStrongDLimbo )
00384 {
00385 double v = fabs(wind.windv);
00386
00387 double dVelRelative = fabs(wind.windv-OldWindVelocity)/
00388 MAX2(v,0.1*timesc.sound_speed_isothermal);
00389
00390 const double WIND_CHNG_VELOCITY_RELATIVE = 0.01;
00391
00392
00393
00394
00395 double winddr;
00396 if( dVelRelative > WIND_CHNG_VELOCITY_RELATIVE && nzone>1 )
00397 {
00398
00399 double factor = WIND_CHNG_VELOCITY_RELATIVE / dVelRelative;
00400
00401 factor = MAX2(0.8 , factor );
00402 winddr = radius.drad * factor;
00403 }
00404 else
00405 {
00406 winddr = BigRadius;
00407 }
00408
00409
00410 if( dynamics.lgAdvection && iteration > dynamics.n_initial_relax)
00411 {
00412 winddr = MIN2( winddr , dynamics.dRad );
00413
00414
00415 dVelRelative = dynamics.dRad;
00416 }
00417 ostringstream oss;
00418 oss << "Wind, dVelRelative=" << scientific << setprecision(3);
00419 oss << dVelRelative << " windv=" << wind.windv;
00420 drChoice.insert( pair<const double,string>( winddr, oss.str() ) );
00421 }
00422
00423 OldWindVelocity = wind.windv;
00424
00425 const double DNGLOB = 0.10;
00426
00427
00428 if( strcmp(dense.chDenseLaw,"GLOB") == 0 )
00429 {
00430 if( radius.glbdst < 0. )
00431 {
00432 fprintf( ioQQQ, " Globule distance is negative, internal overflow has occured, sorry.\n" );
00433 fprintf( ioQQQ, " This is routine radius_next, GLBDST is%10.2e\n",
00434 radius.glbdst );
00435 cdEXIT(EXIT_FAILURE);
00436 }
00437 double factor = radius.glbden*pow(radius.glbrad/radius.glbdst,radius.glbpow);
00438 double fac2 = radius.glbden*pow(radius.glbrad/(radius.glbdst - (realnum)radius.drad),radius.glbpow);
00439
00440 double GlobDr = ( fac2/factor > 1. + DNGLOB ) ? radius.drad*DNGLOB/(fac2/factor - 1.) : BigRadius;
00441 GlobDr = MIN2(GlobDr,radius.glbdst/20.);
00442 ostringstream oss;
00443 oss << "GLOB law, HDEN=" << scientific << setprecision(2) << dense.gas_phase[ipHYDROGEN];
00444 drChoice.insert( pair<const double,string>( GlobDr, oss.str() ) );
00445 }
00446
00447
00448 double hdnew = 0.;
00449 if( strncmp( dense.chDenseLaw , "DLW" , 3) == 0 )
00450 {
00451
00452 if( strcmp(dense.chDenseLaw,"DLW1") == 0 )
00453 {
00454 hdnew = dense_fabden(radius.Radius+radius.drad,radius.depth+
00455 radius.drad);
00456 }
00457 else if( strcmp(dense.chDenseLaw,"DLW2") == 0 )
00458 {
00459 hdnew = dense_tabden(radius.Radius+radius.drad,radius.depth+
00460 radius.drad);
00461 }
00462 else if( strcmp(dense.chDenseLaw,"DLW3") == 0 )
00463 {
00464 hdnew = dense_parametric_wind(radius.Radius+radius.drad);
00465 }
00466 else
00467 {
00468 fprintf( ioQQQ, " dlw insanity in radius_next\n" );
00469 cdEXIT(EXIT_FAILURE);
00470 }
00471 double drTab = fabs(hdnew-dense.gas_phase[ipHYDROGEN])/MAX2(hdnew,dense.gas_phase[ipHYDROGEN]);
00472 drTab = radius.drad*MAX2(0.2,0.10/MAX2(0.01,drTab));
00473 ostringstream oss;
00474 oss << "spec den law, new old den " << scientific << setprecision(2);
00475 oss << hdnew << " " << dense.gas_phase[ipHYDROGEN];
00476 drChoice.insert( pair<const double,string>( drTab, oss.str() ) );
00477 }
00478
00479
00480
00481 if( strcmp(dense.chDenseLaw,"DLW1") == 0 )
00482 {
00483 double dnew = fabs(dense_fabden(radius.Radius+radius.drad,radius.depth+
00484 radius.drad)/dense.gas_phase[ipHYDROGEN]-1.);
00485
00486 double SpecDr;
00487 if( dnew == 0. )
00488 {
00489 SpecDr = radius.drad*3.;
00490 }
00491 else if( dnew/DNGLOB > 1.0 )
00492 {
00493 SpecDr = radius.drad/(dnew/DNGLOB);
00494 }
00495 else
00496 {
00497 SpecDr = BigRadius;
00498 }
00499 ostringstream oss;
00500 oss << "special law, HDEN=" << scientific << setprecision(2) << dense.gas_phase[ipHYDROGEN];
00501 drChoice.insert( pair<const double,string>( SpecDr, oss.str() ) );
00502 }
00503
00504
00505
00506
00507
00508 if( thermal.heating[0][13]/thermal.htot > 0.2 )
00509 {
00510 double grfreqm = 0., gropacm = 0.;
00511
00512
00513 GrainRateDr(&grfreqm,&gropacm);
00514 double DrGrainHeat = 1.0/MAX2(1e-30,gropacm*geometry.FillFac*geometry.DirectionalCosin);
00515 ostringstream oss;
00516 oss << "grain heating nu=" << scientific << setprecision(2) << grfreqm << " opac=" << gropacm;
00517 drChoice.insert( pair<const double,string>( DrGrainHeat, oss.str() ) );
00518 }
00519
00520
00521
00522
00523 if( thermal.heating[0][22]/thermal.htot > 0.2 )
00524 {
00525 long level = -1, ipStrong = -1;
00526 double Strong = 0.;
00527 FndLineHt(&level,&ipStrong,&Strong);
00528
00529 if( Strong/thermal.htot > 0.1 )
00530 {
00531 const transition *t = NULL;
00532 if( level == 1 )
00533
00534 t = &TauLines[ipStrong];
00535 else if( level == 2 )
00536
00537 t = &TauLine2[ipStrong];
00538 else if( level == 3 )
00539
00540
00541 t = &HFLines[ipStrong];
00542 else if( level == 4 )
00543
00544 t = dBaseLines[ipStrong].tran;
00545 else
00546 {
00547
00548 fprintf( ioQQQ, " PROBLEM radius_next Strong line heating set, but not level.\n" );
00549 TotalInsanity();
00550 }
00551
00552 ASSERT( t != NULL );
00553
00554 double TauInwd = t->Emis->TauIn;
00555 double DopplerWidth = t->Emis->dampXvel/t->Emis->damp;
00556 double TauDTau = t->Emis->PopOpc*t->Emis->opacity/DopplerWidth;
00557
00558 fixit();
00559
00560
00561
00562
00563
00564 double drLineHeat = ( TauDTau > 0. ) ? MAX2(1.,TauInwd)*0.4/TauDTau : BigRadius;
00565 ostringstream oss;
00566 oss << "level " << level << " line heating, " << chLineLbl(t) << " TauIn " << scientific;
00567 oss << setprecision(2) << TauInwd << " " << t->Emis->pump << " " << t->Emis->Pesc;
00568 drChoice.insert( pair<const double,string>( drLineHeat, oss.str() ) );
00569 }
00570 }
00571
00572
00573 if( lgFirstCall )
00574 Old_H2_heat_cool = 0.;
00575 else if( !thermal.lgTemperatureConstant )
00576 {
00577
00578
00579 double H2_heat_cool = (fabs(hmi.HeatH2Dexc_used)+fabs(hmi.HeatH2Dish_used)) / thermal.htot;
00580
00581 double dH2_heat_cool = fabs( H2_heat_cool - Old_H2_heat_cool );
00582 if( H2_heat_cool > 0.1 && Old_H2_heat_cool>0. && dH2_heat_cool>SMALLFLOAT )
00583 {
00584
00585
00586 double drH2_heat_cool = radius.drad*MAX2(0.5,0.05/dH2_heat_cool);
00587 ostringstream oss;
00588 oss << "change in H2 heating/cooling, d(c,h)/H " << scientific << setprecision(2);
00589 oss << dH2_heat_cool;
00590 drChoice.insert( pair<const double,string>( drH2_heat_cool, oss.str() ) );
00591 }
00592 }
00593 Old_H2_heat_cool = (fabs(hmi.HeatH2Dexc_used)+fabs(hmi.HeatH2Dish_used)) / thermal.htot;
00594
00595
00596
00597
00598 int mole_dr_change = -1;
00599 if( nzone >= 4 )
00600 {
00601
00602 double Old_H2_abund;
00603
00604 int i;
00605 Old_H2_abund = struc.H2_molec[ipMH2g][nzone-3] + struc.H2_molec[ipMH2s][nzone-3];
00606
00607
00608
00609
00610
00611
00612
00613
00614
00615
00616 if( 2.*hmi.H2_total/dense.gas_phase[ipHYDROGEN] > 1e-8 )
00617 {
00618 double fac = 0.2;
00619
00620 double dH2_abund = 2.*fabs( hmi.H2_total - Old_H2_abund ) / hmi.H2_total;
00621
00622
00623
00624
00625
00626
00627
00628
00629
00630
00631 dH2_abund = SDIV(dH2_abund);
00632 double drH2_abund = radius.drad*MAX2(0.2,fac/dH2_abund );
00633 ostringstream oss;
00634 oss << "change in H2 abundance, d(H2)/H " << scientific << setprecision(2) << dH2_abund;
00635 drChoice.insert( pair<const double,string>( drH2_abund, oss.str() ) );
00636 }
00637
00638
00639
00640 double dr_mole_abund = BigRadius;
00641
00642
00643 double dCO_abund = 0.;
00644 for( i=0; i < mole.num_comole_calc; ++i )
00645 {
00646 realnum abund, abund_limit;
00647 if( !dense.lgElmtOn[COmole[i]->nelem_hevmol] || COmole[i]->n_nuclei <= 1 )
00648 continue;
00649
00650
00651
00652
00653 abund = COmole[i]->hevmol/dense.gas_phase[COmole[i]->nelem_hevmol];
00654
00655
00656
00657
00658 if( COmole[i]->lgGas_Phase )
00659 {
00660 abund_limit = 1e-3f;
00661 }
00662 else
00663 {
00664
00665
00666 abund_limit = 1e-5f;
00667 }
00668
00669 if( abund > abund_limit )
00670 {
00671 double drmole_one, relative_change, relative_change_denom;
00672
00673
00674
00675 if( struc.CO_molec[i][nzone-3] > SMALLFLOAT )
00676 {
00677 relative_change_denom = MIN2( COmole[i]->hevmol , struc.CO_molec[i][nzone-3] );
00678 }
00679 else
00680 {
00681 relative_change_denom = COmole[i]->hevmol;
00682 }
00683
00684 relative_change =
00685 fabs( COmole[i]->hevmol - struc.CO_molec[i][nzone-3] ) / relative_change_denom;
00686
00687
00688
00689
00690
00691
00692
00693
00694
00695 relative_change = SDIV(relative_change);
00696
00697
00698 if( relative_change > 1. )
00699 relative_change = 1./relative_change;
00700
00701
00702 drmole_one = radius.drad*MAX2(0.6,0.05/relative_change );
00703
00704 if( drmole_one < dr_mole_abund )
00705 {
00706
00707 dr_mole_abund = drmole_one;
00708 mole_dr_change = i;
00709 dCO_abund = relative_change;
00710 }
00711 }
00712 }
00713 if( mole_dr_change >= 0 )
00714 {
00715 ostringstream oss;
00716 oss << "change in heavy elem mole abundance, d(mole)/elem=" << scientific << setprecision(2);
00717 oss << dCO_abund << " mole=" << mole_dr_change << "=" << COmole[mole_dr_change]->label;
00718 drChoice.insert( pair<const double,string>( dr_mole_abund, oss.str() ) );
00719 }
00720 }
00721
00722
00723 drChoice.insert( pair<const double,string>( H2_DR(), "change in big H2 Solomon rate line opt depth" ) );
00724
00725
00726
00727
00728
00729
00730
00731
00732 double drmax = ( nzone < 5 ) ? 4.*radius.drad : 1.3*radius.drad;
00733 drChoice.insert( pair<const double,string>( drmax, "DRMAX" ) );
00734
00735
00736
00737 double recom_dr_last_iter = BigRadius;
00738 if( dynamics.lgTimeDependentStatic && dynamics.lgRecom )
00739 {
00740
00741 static long int nzone_recom = -1;
00742 if( nzone <= 1 )
00743 {
00744
00745 nzone_recom = 0;
00746 }
00747 else if( radius.depth < struc.depth_last[struc.nzonePreviousIteration-1] &&
00748 nzone_recom < struc.nzonePreviousIteration )
00749 {
00750 ASSERT( nzone_recom>=0 && nzone_recom<struc.nzonePreviousIteration );
00751
00752
00753
00754 while( struc.depth_last[nzone_recom] < radius.depth &&
00755 nzone_recom < struc.nzonePreviousIteration-1 )
00756 ++nzone_recom;
00757 recom_dr_last_iter = struc.drad_last[nzone_recom]*3.;
00758 drChoice.insert( pair<const double,string>( recom_dr_last_iter,
00759 "previous iteration recom logic" ) );
00760 }
00761 }
00762
00763
00764
00765
00766
00767 if( nzone > 2 )
00768 {
00769
00770 double Efrac_old = struc.ednstr[nzone-3]/struc.gas_phase[ipHYDROGEN][nzone-3];
00771 double Efrac_new = dense.eden/dense.gas_phase[ipHYDROGEN];
00772 double dEfrac = fabs(Efrac_old-Efrac_new) / Efrac_new;
00773 double drEfrac;
00774 if( dEfrac > SMALLFLOAT )
00775 {
00776 double fac = 0.04;
00777
00778
00779
00780 if( dense.eden_from_metals > 0.5 )
00781 {
00782 fac = 0.04;
00783 }
00784
00785
00786
00787 else if( iso.RecomCollisFrac[ipH_LIKE][ipHYDROGEN] > 0.8 &&
00788 dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN] > 0.1 &&
00789 dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN] < 0.8 )
00790
00791 {
00792 fac = 0.02;
00793 }
00794
00795 drEfrac = radius.drad*MAX2(0.1,fac/dEfrac);
00796 }
00797 else
00798 {
00799 drEfrac = BigRadius;
00800 }
00801 ostringstream oss;
00802 oss << "change in elec den, rel chng: " << scientific << setprecision(3) << dEfrac;
00803 oss << " cur=" << Efrac_new << " old=" << Efrac_old;
00804 drChoice.insert( pair<const double,string>( drEfrac, oss.str() ) );
00805 }
00806
00807
00808
00809 if( nzone > 20 )
00810 {
00811
00812 drChoice.insert( pair<const double,string>( radius.depth/10., "relative depth" ) );
00813 }
00814
00815
00816
00817 double thickness_total = BigRadius;
00818 double DepthToGo = BigRadius;
00819 if( StopCalc.HColStop < 5e29 )
00820 {
00821 double coleff = SDIV( dense.gas_phase[ipHYDROGEN]*geometry.FillFac);
00822 DepthToGo = MIN2(DepthToGo ,
00823 (StopCalc.HColStop-colden.colden[ipCOL_HTOT]) / coleff );
00824
00825 thickness_total = MIN2(thickness_total , StopCalc.HColStop / coleff );
00826 }
00827
00828 if( StopCalc.colpls < 5e29 )
00829 {
00830 double coleff = (double)SDIV( (dense.xIonDense[ipHYDROGEN][1])*geometry.FillFac);
00831 DepthToGo = MIN2(DepthToGo ,
00832 (StopCalc.colpls-colden.colden[ipCOL_Hp]) / coleff );
00833 thickness_total = MIN2(thickness_total , StopCalc.colpls / coleff );
00834 }
00835
00836 if( StopCalc.col_h2 < 5e29 )
00837 {
00838
00839 double coleff = (double)SDIV( hmi.H2_total*geometry.FillFac);
00840 DepthToGo = MIN2(DepthToGo ,
00841 (StopCalc.col_h2-colden.colden[ipCOL_H2g]-colden.colden[ipCOL_H2s]) / coleff );
00842 thickness_total = MIN2(thickness_total , StopCalc.col_h2 / coleff );
00843 }
00844
00845 if( StopCalc.col_h2_nut < 5e29 )
00846 {
00847
00848 double coleff = (double)SDIV( (2*hmi.H2_total+dense.xIonDense[ipHYDROGEN][1])*geometry.FillFac);
00849 DepthToGo = MIN2(DepthToGo ,
00850 (StopCalc.col_h2_nut-(2*(colden.colden[ipCOL_H2g]+colden.colden[ipCOL_H2s])+dense.xIonDense[ipHYDROGEN][1])) / coleff );
00851 thickness_total = MIN2(thickness_total , StopCalc.col_h2_nut / coleff );
00852 }
00853
00854 if( StopCalc.col_H0_ov_Tspin < 5e29 )
00855 {
00856
00857 double coleff = (double)SDIV( dense.xIonDense[ipHYDROGEN][0] / hyperfine.Tspin21cm*geometry.FillFac );
00858 DepthToGo = MIN2(DepthToGo ,
00859 (StopCalc.col_H0_ov_Tspin - colden.H0_ov_Tspin) / coleff );
00860 thickness_total = MIN2(thickness_total , StopCalc.col_H0_ov_Tspin / coleff );
00861 }
00862
00863 if( StopCalc.col_monoxco < 5e29 )
00864 {
00865
00866 double coleff = (double)SDIV( (findspecies("CO")->hevmol)*geometry.FillFac);
00867 DepthToGo = MIN2(DepthToGo ,
00868 (StopCalc.col_monoxco-findspecies("CO")->hevcol) / coleff );
00869 thickness_total = MIN2(thickness_total , StopCalc.col_monoxco / coleff );
00870 }
00871
00872 if( StopCalc.colnut < 5e29 )
00873 {
00874 double coleff = (double)SDIV( (dense.xIonDense[ipHYDROGEN][0])*geometry.FillFac);
00875 DepthToGo = MIN2(DepthToGo ,
00876 (StopCalc.colnut - colden.colden[ipCOL_H0]) / coleff );
00877 thickness_total = MIN2(thickness_total , StopCalc.colnut / coleff );
00878 }
00879
00880
00881 if( radius.StopThickness[iteration-1] < 5e29 )
00882 {
00883 thickness_total = MIN2(thickness_total , radius.StopThickness[iteration-1] );
00884 DepthToGo = MIN2(DepthToGo ,
00885 radius.StopThickness[iteration-1] - radius.depth );
00886 }
00887
00888
00889 if( StopCalc.iptnu != rfield.nupper )
00890 {
00891
00892 double dt = SDIV(opac.opacity_abs[StopCalc.iptnu-1]*geometry.FillFac);
00893 DepthToGo = MIN2(DepthToGo ,
00894 (StopCalc.tauend-opac.TauAbsGeo[0][StopCalc.iptnu-1])/dt);
00895 }
00896
00897 if( DepthToGo <= 0. )
00898 TotalInsanity();
00899
00900
00901 if( DepthToGo < BigRadius )
00902 {
00903 double depth_min = MIN2( DepthToGo , thickness_total/100. );
00904 DepthToGo = MAX2( DepthToGo / 10. , depth_min );
00905
00906
00907
00908
00909
00910 double drThickness = MIN2( thickness_total/10. , DepthToGo );
00911 drChoice.insert( pair<const double,string>( drThickness, "depth to go" ) );
00912 }
00913
00914
00915
00916
00917 double AV_to_go = BigRadius;
00918 if( rfield.opac_mag_V_extended > SMALLFLOAT && rfield.opac_mag_V_point > SMALLFLOAT )
00919 {
00920 double SAFETY = 1. + 10.*DBL_EPSILON;
00921
00922
00923
00924
00925 double ave = log10(SAFETY*StopCalc.AV_extended - rfield.extin_mag_V_extended) -
00926 log10(rfield.opac_mag_V_extended);
00927 double avp = log10(SAFETY*StopCalc.AV_point - rfield.extin_mag_V_point) -
00928 log10(rfield.opac_mag_V_point);
00929 AV_to_go = MIN2( ave , avp );
00930 if( AV_to_go < 37. )
00931 {
00932 AV_to_go = pow(10., AV_to_go)/geometry.FillFac;
00933
00934
00935 AV_to_go *= 1.0001;
00936 }
00937 else
00938 AV_to_go = BigRadius;
00939
00940 }
00941 drChoice.insert( pair<const double,string>( AV_to_go, "A_V to go" ) );
00942
00943
00944
00945 double drSphere = radius.Radius*0.04;
00946 drChoice.insert( pair<const double,string>( drSphere, "sphericity" ) );
00947
00948
00949 double dRTaue = radius.drChange/(dense.eden*6.65e-25*geometry.FillFac);
00950
00951
00952 if( thermal.lgTemperatureConstant )
00953 dRTaue *= 3.;
00954 drChoice.insert( pair<const double,string>( dRTaue, "optical depth to electron scattering" ) );
00955
00956
00957 if( dense.flong != 0. )
00958 {
00959 double drFluc = 0.628/2./dense.flong;
00960
00961
00962 drFluc /= 2.;
00963 drChoice.insert( pair<const double,string>( drFluc, "density fluctuations" ) );
00964 }
00965
00966
00967
00968
00969
00970 double drStart = ( nzone <= 1 ) ? radius.drad : BigRadius;
00971 drChoice.insert( pair<const double,string>( drStart, "capped to old DR in first zone" ) );
00972
00973
00974 double rfacmax = radius.lgSdrmaxRel ? radius.Radius : 1.;
00975 drChoice.insert( pair<const double,string>( rfacmax*radius.sdrmax, "sdrmax" ) );
00976
00977
00978
00979
00980
00981
00982
00983
00984 if( nzone >= 5 && !dynamics.lgTimeDependentStatic )
00985 {
00986
00987
00988 if( phycon.te > 200. && phycon.te < 3000. &&
00989
00990
00991 (struc.testr[nzone-3] - phycon.te)/phycon.te > 0.02 &&
00992 (struc.testr[nzone-4] - phycon.te)/phycon.te > 0.02 &&
00993 (struc.testr[nzone-5] - phycon.te)/phycon.te > 0.02 )
00994 {
00995
00996
00997 double drThermalFront = radius.drad * 0.91;
00998 drChoice.clear();
00999 drChoice.insert( pair<const double,string>( drThermalFront, "thermal front logic" ) );
01000 }
01001 }
01002
01003 ASSERT( drChoice.size() > 0 );
01004
01005
01006 if( drChoice.begin()->first <= 0. )
01007 {
01008 multimap<double,string>::const_iterator ptr = drChoice.begin();
01009 fprintf( ioQQQ, " radius_next finds insane drNext: %.2e\n", ptr->first );
01010 fprintf( ioQQQ, " all drs follow:\n" );
01011 while( ptr != drChoice.end() )
01012 {
01013 fprintf( ioQQQ, " %.2e %s\n", ptr->first, ptr->second.c_str() );
01014 ++ptr;
01015 }
01016 cdEXIT(EXIT_FAILURE);
01017 }
01018
01019
01020 double rfacmin = radius.lgSdrminRel ? radius.Radius : 1.;
01021 if( drChoice.begin()->first < rfacmin*radius.sdrmin )
01022 {
01023 drChoice.clear();
01024 drChoice.insert( pair<const double,string>( rfacmin*radius.sdrmin, "sdrmin" ) );
01025 }
01026
01027
01028
01029
01030
01031
01032
01033
01034
01035
01036
01037
01038
01039
01040 if( strcmp(dense.chDenseLaw,"DLW1") != 0 &&
01041 strcmp(dense.chDenseLaw ,"GLOB") != 0 &&
01042 dense.flong == 0. &&
01043
01044 drChoice.begin()->first != DepthToGo &&
01045
01046 drChoice.begin()->first != AV_to_go )
01047 {
01048
01049
01050
01051
01052
01053
01054
01055
01056 if( drChoice.begin()->first*radius.Radius/radius.rinner*dense.gas_phase[ipHYDROGEN] < radius.drMinimum )
01057 {
01058 radius.drNext = radius.drMinimum/dense.gas_phase[ipHYDROGEN];
01059
01060
01061 radius.lgDrMinUsed = true;
01062
01063 lgAbort = true;
01064
01065
01066 --nzone;
01067 fprintf( ioQQQ,
01068 "\n DISASTER PROBLEM radius_next finds dr too small and aborts. "
01069 "This is zone %ld iteration %ld. The thickness was %.2e\n",
01070 nzone,
01071 iteration,
01072 radius.drNext);
01073 fprintf( ioQQQ,
01074 " If this simulation seems reasonable you can override this limit "
01075 "with the command SET DRMIN %.2e\n\n",
01076 radius.drNext*2);
01077
01078 lgAbort = true;
01079 return 1;
01080 }
01081 }
01082
01083
01084 const double Z = 1.0001;
01085
01086
01087
01088
01089
01090 double drOuterRadius = (radius.StopThickness[iteration-1]-radius.depth)*Z;
01091 drChoice.insert( pair<const double,string>( drOuterRadius, "outer radius" ) );
01092
01093
01094 radius.drNext = drChoice.begin()->first;
01095
01096
01097 if( drChoice.begin()->second.find( "H maser" ) != string::npos )
01098 {
01099 rt.lgMaserSetDR = true;
01100 }
01101
01102
01103
01104
01105
01106 bool lgDoPun = ( save.lgDROn && !( save.lgDRPLst && !iterations.lgLastIt ) );
01107 if( (trace.lgTrace && trace.lgDrBug) || lgDoPun )
01108 {
01109 fprintf( save.ipDRout , "%ld\t%.5e\t%.3e\t%.3e\t%s\n", nzone, radius.depth,
01110 radius.drNext, radius.Depth2Go, drChoice.begin()->second.c_str() );
01111 }
01112
01113 ASSERT( radius.drNext > 0. );
01114
01115 if( trace.lgTrace )
01116 {
01117 fprintf( ioQQQ, " radius_next chooses next drad drNext=%.4e; this drad was %.4e\n",
01118 radius.drNext, radius.drad );
01119 }
01120
01121 return 0;
01122 }
01123
01124
01125 STATIC void ContRate(double *freqm,
01126 double *opacm)
01127 {
01128 long int i,
01129 ipHi,
01130 iplow,
01131 limit;
01132 double FreqH,
01133 Freq_nonIonizing,
01134 Opac_Hion,
01135 Opac_nonIonizing,
01136 Rate_max_Hion,
01137 Rate_max_nonIonizing;
01138
01139 DEBUG_ENTRY( "ContRate()" );
01140
01141
01142
01143
01144
01145
01146 Rate_max_nonIonizing = 0.;
01147 Freq_nonIonizing = 0.;
01148 Opac_nonIonizing = 0.;
01149
01150
01151 *opacm = -1.;
01152 *freqm = -1.;
01153
01154
01155
01156 if( dense.lgElmtOn[ipCARBON] )
01157 {
01158
01159 ipHi = Heavy.ipHeavy[ipCARBON][0] - 1;
01160 }
01161 else
01162 {
01163
01164 ipHi = iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH2s]-1;
01165 }
01166
01167
01168 for( i=rfield.ipPlasma; i < ipHi; i++ )
01169 {
01170
01171
01172 if( rfield.anu[i]*rfield.flux[0][i]/rfield.widflx[i]*(opac.opacity_abs[i] -
01173 gv.dstab[i]*dense.gas_phase[ipHYDROGEN]) > Rate_max_nonIonizing )
01174 {
01175 Rate_max_nonIonizing = rfield.anu[i]*rfield.flux[0][i]/rfield.widflx[i]*
01176 (opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN]);
01177 Freq_nonIonizing = rfield.anu[i];
01178 Opac_nonIonizing = opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN];
01179 }
01180 }
01181
01182
01183
01184
01185 if( CoolHeavy.brems_heat_total/thermal.htot > 0.05 )
01186 {
01187
01188
01189 iplow = MAX2(1 , rfield.ipEnergyBremsThin );
01190 }
01191 else
01192 {
01193
01194
01195 iplow = ipHi;
01196 }
01197
01198 iplow = MAX2( rfield.ipPlasma , iplow );
01199
01200
01201 limit = MIN2(rfield.nflux,iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1);
01202 for( i=iplow; i < limit; i++ )
01203 {
01204 if( rfield.anu[i]*rfield.flux[0][i]/rfield.widflx[i]*(opac.opacity_abs[i] -
01205 gv.dstab[i]*dense.gas_phase[ipHYDROGEN]) > Rate_max_nonIonizing )
01206 {
01207 Rate_max_nonIonizing = rfield.anu[i]*rfield.flux[0][i]/rfield.widflx[i]*
01208 (opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN]);
01209 Freq_nonIonizing = rfield.anu[i];
01210 Opac_nonIonizing = opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN];
01211 }
01212 }
01213
01214
01215 Rate_max_Hion = 0.;
01216 Opac_Hion = 0.;
01217 FreqH = 0.;
01218
01219
01220 for( i=iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1; i < rfield.nflux; i++ )
01221 {
01222 if( rfield.anu[i]*rfield.flux[0][i]/rfield.widflx[i]*(opac.opacity_abs[i] -
01223 gv.dstab[i]*dense.gas_phase[ipHYDROGEN]) > Rate_max_Hion )
01224 {
01225
01226 Rate_max_Hion = rfield.anu[i]*rfield.flux[0][i]/rfield.widflx[i]*
01227 (opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN]);
01228 FreqH = rfield.anu[i];
01229 Opac_Hion = opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN];
01230 }
01231 }
01232
01233
01234 if( Rate_max_nonIonizing < 1e-30 && Opac_Hion > SMALLFLOAT )
01235 {
01236
01237 *opacm = Opac_Hion;
01238 *freqm = FreqH;
01239 }
01240
01241
01242 else if( Opac_Hion > Opac_nonIonizing && Rate_max_Hion/Rate_max_nonIonizing > 1e-10 && Opac_Hion > SMALLFLOAT )
01243 {
01244
01245 *opacm = Opac_Hion;
01246 *freqm = FreqH;
01247 }
01248 else
01249 {
01250
01251 *opacm = Opac_nonIonizing;
01252 *freqm = Freq_nonIonizing;
01253 }
01254
01255 # if 0
01256
01257
01258
01259
01260
01261 if( gv.lgDustOn() && gv.lgGrainPhysicsOn )
01262 {
01263 double fluxGrainPeak = -1.;
01264 long int ipGrainPeak = -1;
01265 long int ipGrainPeak2 = -1;
01266
01267 i = 0;
01268 while( rfield.anu[i] < 0.0912 && i < (rfield.nflux-2) )
01269 {
01270
01271
01272 if( gv.GrainEmission[i]*rfield.anu[i]*opac.opacity_abs[i] > fluxGrainPeak )
01273 {
01274 ipGrainPeak = i;
01275 fluxGrainPeak = gv.GrainEmission[i]*rfield.anu[i]*opac.opacity_abs[i];
01276 }
01277 ++i;
01278 }
01279 ASSERT( fluxGrainPeak>=0. && ipGrainPeak >=0 );
01280
01281
01282
01283
01284
01285
01286
01287 i = ipGrainPeak;
01288
01289
01290
01291
01292 if( fluxGrainPeak > 0. )
01293 {
01294 while( rfield.anu[i] < 0.0912 && i < (rfield.nflux-2) )
01295 {
01296
01297 if( gv.GrainEmission[i]*rfield.anu[i]*opac.opacity_abs[i] > 0.05*fluxGrainPeak )
01298 {
01299
01300 ipGrainPeak2 = i;
01301 }
01302 ++i;
01303 }
01304 ASSERT( ipGrainPeak2>=0 );
01305
01306
01307 if( opac.opacity_abs[ipGrainPeak2] > *opacm )
01308 {
01309
01310 *opacm = opac.opacity_abs[ipGrainPeak2]*10.;
01311 *freqm = rfield.anu[ipGrainPeak2];
01312
01313
01314 }
01315 }
01316 }
01317 # endif
01318
01319 {
01320
01321 enum {DEBUG_LOC=false};
01322 if( DEBUG_LOC )
01323 {
01324 fprintf(ioQQQ,"conratedebug \t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
01325 Rate_max_nonIonizing,Freq_nonIonizing,Opac_nonIonizing,
01326 Rate_max_Hion,FreqH ,Opac_Hion,*freqm,*opacm
01327 );
01328
01329 }
01330 }
01331
01332
01333
01334
01335
01336
01337 ASSERT( *opacm >= 0. && *freqm >= 0. );
01338 return;
01339 }
01340
01341
01342 STATIC void GrainRateDr(double *grfreqm,
01343 double *gropacm)
01344 {
01345 long int i,
01346 iplow;
01347 double xMax;
01348
01349 DEBUG_ENTRY( "GrainRateDr()" );
01350
01351
01352
01353
01354
01355
01356
01357
01358 if( CoolHeavy.brems_heat_total/thermal.htot > 0.05 )
01359 {
01360
01361
01362 iplow = MAX2(1 , rfield.ipEnergyBremsThin );
01363 }
01364 else
01365 {
01366
01367
01368 if( dense.lgElmtOn[ipCARBON] )
01369 {
01370
01371 iplow = Heavy.ipHeavy[ipCARBON][0];
01372 }
01373 else
01374 {
01375
01376 iplow = iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH2s];
01377 }
01378 }
01379
01380 xMax = -1.;
01381
01382 for( i=iplow-1; i < Heavy.ipHeavy[ipHYDROGEN][0]; i++ )
01383 {
01384 if( rfield.anu[i]*rfield.flux[0][i]/rfield.widflx[i]*opac.opacity_abs[i] > xMax )
01385 {
01386 xMax = rfield.anu[i]*rfield.flux[0][i]/rfield.widflx[i]*
01387 opac.opacity_abs[i];
01388 *grfreqm = rfield.anu[i];
01389 *gropacm = opac.opacity_abs[i];
01390 }
01391 }
01392
01393
01394 if( dense.lgElmtOn[ipHELIUM] )
01395 {
01396 for( i=Heavy.ipHeavy[ipHYDROGEN][0]-1; i < Heavy.ipHeavy[ipHELIUM][1]; i++ )
01397 {
01398 if( rfield.anu[i]*rfield.flux[0][i]/rfield.widflx[i]*opac.opacity_abs[i] > xMax )
01399 {
01400 xMax = rfield.anu[i]*rfield.flux[0][i]/rfield.widflx[i]*
01401 opac.opacity_abs[i];
01402 *grfreqm = rfield.anu[i];
01403 *gropacm = opac.opacity_abs[i];
01404 }
01405 }
01406 }
01407
01408
01409
01410 if( xMax <= 0. )
01411 {
01412 *gropacm = 0.;
01413 *grfreqm = 0.;
01414 }
01415 return;
01416 }