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( mole_global.lgLeidenHack )
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.+0.01*extin_mag_V_point_old;
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;
00526 TransitionProxy t = FndLineHt(&level);
00527
00528 if( t.Coll().heat()/thermal.htot > 0.1 )
00529 {
00530 ASSERT( t.associated() );
00531
00532 double TauInwd = t.Emis().TauIn();
00533 double DopplerWidth = t.Emis().dampXvel()/t.Emis().damp();
00534 double TauDTau = t.Emis().PopOpc()*t.Emis().opacity()/DopplerWidth;
00535
00536 fixit();
00537
00538
00539
00540
00541
00542 double drLineHeat = ( TauDTau > 0. ) ? MAX2(1.,TauInwd)*0.4/TauDTau : BigRadius;
00543 ostringstream oss;
00544 oss << "level " << level << " line heating, " << chLineLbl(t) << " TauIn " << scientific;
00545 oss << setprecision(2) << TauInwd << " " << t.Emis().pump() << " " << t.Emis().Pesc();
00546 drChoice.insert( pair<const double,string>( drLineHeat, oss.str() ) );
00547 }
00548 }
00549
00550
00551 if( lgFirstCall )
00552 Old_H2_heat_cool = 0.;
00553 else if( !thermal.lgTemperatureConstant )
00554 {
00555
00556
00557 double H2_heat_cool = (fabs(hmi.HeatH2Dexc_used)+fabs(hmi.HeatH2Dish_used)) / thermal.htot;
00558
00559 double dH2_heat_cool = fabs( H2_heat_cool - Old_H2_heat_cool );
00560 if( H2_heat_cool > 0.1 && Old_H2_heat_cool>0. && dH2_heat_cool>SMALLFLOAT )
00561 {
00562
00563
00564 double drH2_heat_cool = radius.drad*MAX2(0.5,0.05/dH2_heat_cool);
00565 ostringstream oss;
00566 oss << "change in H2 heating/cooling, d(c,h)/H " << scientific << setprecision(2);
00567 oss << dH2_heat_cool;
00568 drChoice.insert( pair<const double,string>( drH2_heat_cool, oss.str() ) );
00569 }
00570 }
00571 Old_H2_heat_cool = (fabs(hmi.HeatH2Dexc_used)+fabs(hmi.HeatH2Dish_used)) / thermal.htot;
00572
00573
00574
00575
00576 int mole_dr_change = -1;
00577 if( nzone >= 4 )
00578 {
00579
00580 double Old_H2_abund;
00581
00582 int i;
00583 Old_H2_abund = struc.H2_abund[nzone-3];
00584
00585
00586
00587
00588
00589
00590
00591
00592
00593
00594 if( 2.*hmi.H2_total/dense.gas_phase[ipHYDROGEN] > 1e-8 )
00595 {
00596 double fac = 0.2;
00597
00598 double dH2_abund = 2.*fabs( hmi.H2_total - Old_H2_abund ) / hmi.H2_total;
00599
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609 dH2_abund = SDIV(dH2_abund);
00610 double drH2_abund = radius.drad*MAX2(0.2,fac/dH2_abund );
00611 ostringstream oss;
00612 oss << "change in H2 abundance, d(H2)/H " << scientific << setprecision(2) << dH2_abund;
00613 drChoice.insert( pair<const double,string>( drH2_abund, oss.str() ) );
00614 }
00615
00616
00617
00618 double max_change = 0.0;
00619
00620
00621 for( i=0; i < mole_global.num_calc; ++i )
00622 {
00623 realnum abund,
00624 abund1,
00625 abund_limit;
00626 if( !mole_global.list[i]->isMonatomic() && mole_global.list[i]->parentLabel.empty() )
00627 {
00628
00629
00630
00631
00632 abund = 0.;
00633 for( molecule::nAtomsMap::iterator atom=mole_global.list[i]->nAtom.begin();
00634 atom != mole_global.list[i]->nAtom.end(); ++atom)
00635 {
00636 long int nelem = atom->first->el->Z-1;
00637 if (dense.lgElmtOn[nelem])
00638 {
00639 abund1 = (realnum)mole.species[i].den*atom->second/SDIV(dense.gas_phase[nelem]);
00640 if (abund1 > abund)
00641 {
00642 abund = abund1;
00643 }
00644 }
00645 }
00646
00647
00648
00649
00650 if( mole_global.list[i]->lgGas_Phase )
00651 {
00652 abund_limit = 1e-3f;
00653 }
00654 else
00655 {
00656
00657
00658 abund_limit = 1e-5f;
00659 }
00660
00661 if( abund > abund_limit )
00662 {
00663
00664 double relative_change =
00665 2.0 * fabs( mole.species[i].den - struc.molecules[i][nzone-3] ) /
00666 SDIV ( mole.species[i].den + struc.molecules[i][nzone-3] ) ;
00667 if (relative_change > max_change)
00668 {
00669 max_change = relative_change;
00670 mole_dr_change = i;
00671 }
00672 }
00673 }
00674 }
00675 if( mole_dr_change >= 0 )
00676 {
00677 double dr_mole_abund = radius.drad * MAX2(0.6, 0.05/SDIV(max_change));
00678 ostringstream oss;
00679 oss << "change in molecular abundance, old=" << scientific << setprecision(2);
00680 oss << struc.molecules[mole_dr_change][nzone-3] << " new=" << mole.species[mole_dr_change].den;
00681 oss << " mole=" << mole_dr_change << "=" << mole_global.list[mole_dr_change]->label;
00682 drChoice.insert( pair<const double,string>( dr_mole_abund, oss.str() ) );
00683 }
00684 }
00685
00686
00687 for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
00688 drChoice.insert( pair<const double,string>( (*diatom)->H2_DR(), "change in big H2 Solomon rate line opt depth" ) );
00689
00690
00691
00692
00693
00694
00695
00696
00697 double drmax = ( nzone < 5 ) ? 4.*radius.drad : 1.3*radius.drad;
00698 drChoice.insert( pair<const double,string>( drmax, "DRMAX" ) );
00699
00700
00701
00702 double recom_dr_last_iter = BigRadius;
00703 if( dynamics.lgTimeDependentStatic && dynamics.lgRecom )
00704 {
00705
00706 static long int nzone_recom = -1;
00707 if( nzone <= 1 )
00708 {
00709
00710 nzone_recom = 0;
00711 }
00712 else if( radius.depth < struc.depth_last[struc.nzonePreviousIteration-1] &&
00713 nzone_recom < struc.nzonePreviousIteration )
00714 {
00715 ASSERT( nzone_recom>=0 && nzone_recom<struc.nzonePreviousIteration );
00716
00717
00718
00719 while( struc.depth_last[nzone_recom] < radius.depth &&
00720 nzone_recom < struc.nzonePreviousIteration-1 )
00721 ++nzone_recom;
00722 recom_dr_last_iter = struc.drad_last[nzone_recom]*3.;
00723 drChoice.insert( pair<const double,string>( recom_dr_last_iter,
00724 "previous iteration recom logic" ) );
00725 }
00726 }
00727
00728
00729
00730
00731
00732 if( nzone > 2 )
00733 {
00734
00735 double Efrac_old = struc.ednstr[nzone-3]/struc.gas_phase[ipHYDROGEN][nzone-3];
00736 double Efrac_new = dense.eden/dense.gas_phase[ipHYDROGEN];
00737 double dEfrac = fabs(Efrac_old-Efrac_new) / Efrac_new;
00738 double drEfrac;
00739 if( dEfrac > SMALLFLOAT )
00740 {
00741 double fac = 0.04;
00742
00743
00744
00745 if( dense.eden_from_metals > 0.5 )
00746 {
00747 fac = 0.04;
00748 }
00749
00750
00751
00752 else if( iso_sp[ipH_LIKE][ipHYDROGEN].RecomCollisFrac > 0.8 &&
00753 dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN] > 0.1 &&
00754 dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN] < 0.8 )
00755
00756 {
00757 fac = 0.02;
00758 }
00759
00760 drEfrac = radius.drad*MAX2(0.1,fac/dEfrac);
00761 }
00762 else
00763 {
00764 drEfrac = BigRadius;
00765 }
00766 ostringstream oss;
00767 oss << "change in elec den, rel chng: " << scientific << setprecision(3) << dEfrac;
00768 oss << " cur=" << Efrac_new << " old=" << Efrac_old;
00769 drChoice.insert( pair<const double,string>( drEfrac, oss.str() ) );
00770 }
00771
00772
00773
00774 if( nzone > 20 )
00775 {
00776
00777 drChoice.insert( pair<const double,string>( radius.depth/10., "relative depth" ) );
00778 }
00779
00780
00781
00782 double thickness_total = BigRadius;
00783 double DepthToGo = BigRadius;
00784 if( StopCalc.HColStop < 5e29 )
00785 {
00786 double coleff = SDIV( dense.gas_phase[ipHYDROGEN]*geometry.FillFac);
00787 DepthToGo = MIN2(DepthToGo ,
00788 (StopCalc.HColStop-colden.colden[ipCOL_HTOT]) / coleff );
00789
00790 thickness_total = MIN2(thickness_total , StopCalc.HColStop / coleff );
00791 }
00792
00793 if( StopCalc.colpls < 5e29 )
00794 {
00795 double coleff = (double)SDIV( (dense.xIonDense[ipHYDROGEN][1])*geometry.FillFac);
00796 DepthToGo = MIN2(DepthToGo ,
00797 (StopCalc.colpls-colden.colden[ipCOL_Hp]) / coleff );
00798 thickness_total = MIN2(thickness_total , StopCalc.colpls / coleff );
00799 }
00800
00801 if( StopCalc.col_h2 < 5e29 )
00802 {
00803
00804 double coleff = (double)SDIV( hmi.H2_total*geometry.FillFac);
00805 DepthToGo = MIN2(DepthToGo ,
00806 (StopCalc.col_h2-colden.colden[ipCOL_H2g]-colden.colden[ipCOL_H2s]) / coleff );
00807 thickness_total = MIN2(thickness_total , StopCalc.col_h2 / coleff );
00808 }
00809
00810 if( StopCalc.col_h2_nut < 5e29 )
00811 {
00812
00813 double coleff = (double)SDIV( (2*hmi.H2_total+dense.xIonDense[ipHYDROGEN][1])*geometry.FillFac);
00814 DepthToGo = MIN2(DepthToGo ,
00815 (StopCalc.col_h2_nut-(2*(colden.colden[ipCOL_H2g]+colden.colden[ipCOL_H2s])+dense.xIonDense[ipHYDROGEN][1])) / coleff );
00816 thickness_total = MIN2(thickness_total , StopCalc.col_h2_nut / coleff );
00817 }
00818
00819 if( StopCalc.col_H0_ov_Tspin < 5e29 )
00820 {
00821
00822 double coleff = (double)SDIV( dense.xIonDense[ipHYDROGEN][0] / hyperfine.Tspin21cm*geometry.FillFac );
00823 DepthToGo = MIN2(DepthToGo ,
00824 (StopCalc.col_H0_ov_Tspin - colden.H0_ov_Tspin) / coleff );
00825 thickness_total = MIN2(thickness_total , StopCalc.col_H0_ov_Tspin / coleff );
00826 }
00827
00828 if( StopCalc.col_monoxco < 5e29 )
00829 {
00830
00831 double coleff = (double)SDIV( (findspecieslocal("CO")->den)*geometry.FillFac);
00832 DepthToGo = MIN2(DepthToGo ,
00833 (StopCalc.col_monoxco-findspecieslocal("CO")->column) / coleff );
00834 thickness_total = MIN2(thickness_total , StopCalc.col_monoxco / coleff );
00835 }
00836
00837 if( StopCalc.colnut < 5e29 )
00838 {
00839 double coleff = (double)SDIV( (dense.xIonDense[ipHYDROGEN][0])*geometry.FillFac);
00840 DepthToGo = MIN2(DepthToGo ,
00841 (StopCalc.colnut - colden.colden[ipCOL_H0]) / coleff );
00842 thickness_total = MIN2(thickness_total , StopCalc.colnut / coleff );
00843 }
00844
00845
00846 if( radius.StopThickness[iteration-1] < 5e29 )
00847 {
00848 thickness_total = MIN2(thickness_total , radius.StopThickness[iteration-1] );
00849 DepthToGo = MIN2(DepthToGo ,
00850 radius.StopThickness[iteration-1] - radius.depth );
00851 }
00852
00853
00854 if( StopCalc.iptnu != rfield.nupper )
00855 {
00856
00857 double dt = SDIV(opac.opacity_abs[StopCalc.iptnu-1]*geometry.FillFac);
00858 DepthToGo = MIN2(DepthToGo ,
00859 (StopCalc.tauend-opac.TauAbsGeo[0][StopCalc.iptnu-1])/dt);
00860 }
00861
00862 if( DepthToGo <= 0. )
00863 TotalInsanity();
00864
00865
00866 if( DepthToGo < BigRadius )
00867 {
00868 double depth_min = MIN2( DepthToGo , thickness_total/100. );
00869 DepthToGo = MAX2( DepthToGo / 10. , depth_min );
00870
00871
00872
00873
00874
00875 double drThickness = MIN2( thickness_total/10. , DepthToGo );
00876 drChoice.insert( pair<const double,string>( drThickness, "depth to go" ) );
00877 }
00878
00879
00880
00881
00882 double AV_to_go = BigRadius;
00883 if( rfield.opac_mag_V_extended > SMALLFLOAT && rfield.opac_mag_V_point > SMALLFLOAT )
00884 {
00885 double SAFETY = 1. + 10.*DBL_EPSILON;
00886
00887
00888
00889
00890 double ave = log10(SAFETY*StopCalc.AV_extended - rfield.extin_mag_V_extended) -
00891 log10(rfield.opac_mag_V_extended);
00892 double avp = log10(SAFETY*StopCalc.AV_point - rfield.extin_mag_V_point) -
00893 log10(rfield.opac_mag_V_point);
00894 AV_to_go = MIN2( ave , avp );
00895 if( AV_to_go < 37. )
00896 {
00897 AV_to_go = pow(10., AV_to_go)/geometry.FillFac;
00898
00899
00900 AV_to_go *= 1.0001;
00901 }
00902 else
00903 AV_to_go = BigRadius;
00904
00905 }
00906 drChoice.insert( pair<const double,string>( AV_to_go, "A_V to go" ) );
00907
00908
00909
00910 double drSphere = radius.Radius*0.04;
00911 drChoice.insert( pair<const double,string>( drSphere, "sphericity" ) );
00912
00913
00914 double dRTaue = radius.drChange/(dense.eden*6.65e-25*geometry.FillFac);
00915
00916
00917 if( thermal.lgTemperatureConstant )
00918 dRTaue *= 3.;
00919 drChoice.insert( pair<const double,string>( dRTaue, "optical depth to electron scattering" ) );
00920
00921
00922 if( dense.flong != 0. )
00923 {
00924 double drFluc = 0.628/2./dense.flong;
00925
00926
00927 drFluc /= 2.;
00928 drChoice.insert( pair<const double,string>( drFluc, "density fluctuations" ) );
00929 }
00930
00931
00932
00933
00934
00935 double drStart = ( nzone <= 1 ) ? radius.drad : BigRadius;
00936 drChoice.insert( pair<const double,string>( drStart, "capped to old DR in first zone" ) );
00937
00938
00939 double rfacmax = radius.lgSdrmaxRel ? radius.Radius : 1.;
00940 drChoice.insert( pair<const double,string>( rfacmax*radius.sdrmax, "sdrmax" ) );
00941
00942
00943
00944
00945
00946
00947
00948
00949 if( nzone >= 5 && !dynamics.lgTimeDependentStatic )
00950 {
00951
00952
00953 if( phycon.te > 200. && phycon.te < 3000. &&
00954
00955
00956 (struc.testr[nzone-3] - phycon.te)/phycon.te > 0.02 &&
00957 (struc.testr[nzone-4] - phycon.te)/phycon.te > 0.02 &&
00958 (struc.testr[nzone-5] - phycon.te)/phycon.te > 0.02 )
00959 {
00960
00961
00962 double drThermalFront = radius.drad * 0.91;
00963 drChoice.clear();
00964 drChoice.insert( pair<const double,string>( drThermalFront, "thermal front logic" ) );
00965 }
00966 }
00967
00968 ASSERT( drChoice.size() > 0 );
00969
00970
00971 if( drChoice.begin()->first <= 0. )
00972 {
00973 multimap<double,string>::const_iterator ptr = drChoice.begin();
00974 fprintf( ioQQQ, " radius_next finds insane drNext: %.2e\n", ptr->first );
00975 fprintf( ioQQQ, " all drs follow:\n" );
00976 while( ptr != drChoice.end() )
00977 {
00978 fprintf( ioQQQ, " %.2e %s\n", ptr->first, ptr->second.c_str() );
00979 ++ptr;
00980 }
00981 cdEXIT(EXIT_FAILURE);
00982 }
00983
00984
00985 if( !radius.lgFixed && drChoice.begin()->first < radius.depth*radius.sdrmin_rel_depth )
00986 {
00987 drChoice.clear();
00988 drChoice.insert( pair<const double,string>( radius.depth*radius.sdrmin_rel_depth, "sdrmin_rel_depth" ) );
00989 }
00990
00991
00992 double rfacmin = radius.lgSdrminRel ? radius.Radius : 1.;
00993 if( drChoice.begin()->first < rfacmin*radius.sdrmin )
00994 {
00995 drChoice.clear();
00996 drChoice.insert( pair<const double,string>( rfacmin*radius.sdrmin, "sdrmin" ) );
00997 }
00998
00999
01000
01001
01002
01003
01004
01005
01006
01007
01008
01009
01010
01011
01012 if( strcmp(dense.chDenseLaw,"DLW1") != 0 &&
01013 strcmp(dense.chDenseLaw ,"GLOB") != 0 &&
01014 dense.flong == 0. &&
01015
01016 drChoice.begin()->first != DepthToGo &&
01017
01018 drChoice.begin()->first != AV_to_go )
01019 {
01020
01021
01022
01023
01024
01025
01026
01027
01028 if( drChoice.begin()->first*radius.Radius/radius.rinner*dense.gas_phase[ipHYDROGEN] < radius.drMinimum )
01029 {
01030 radius.drNext = radius.drMinimum/dense.gas_phase[ipHYDROGEN];
01031
01032
01033 radius.lgDrMinUsed = true;
01034
01035 lgAbort = true;
01036
01037
01038 --nzone;
01039 fprintf( ioQQQ,
01040 "\n DISASTER PROBLEM radius_next finds dr too small and aborts. "
01041 "This is zone %ld iteration %ld. The thickness was %.2e\n",
01042 nzone,
01043 iteration,
01044 radius.drNext);
01045 fprintf( ioQQQ,
01046 " If this simulation seems reasonable you can override this limit "
01047 "with the command SET DRMIN %.2e\n\n",
01048 radius.drNext*2);
01049
01050 lgAbort = true;
01051 return 1;
01052 }
01053 }
01054
01055
01056 const double Z = 1.0001;
01057
01058
01059
01060
01061
01062 double drOuterRadius = (radius.StopThickness[iteration-1]-radius.depth)*Z;
01063 drChoice.insert( pair<const double,string>( drOuterRadius, "outer radius" ) );
01064
01065
01066 radius.drNext = drChoice.begin()->first;
01067
01068
01069 if( drChoice.begin()->second.find( "H maser" ) != string::npos )
01070 {
01071 rt.lgMaserSetDR = true;
01072 }
01073
01074
01075
01076
01077
01078 bool lgDoPun = ( save.lgDROn && !( save.lgDRPLst && !iterations.lgLastIt ) );
01079 if( (trace.lgTrace && trace.lgDrBug) || lgDoPun )
01080 {
01081 fprintf( save.ipDRout , "%ld\t%.5e\t%.3e\t%.3e\t%s\n", nzone, radius.depth,
01082 radius.drNext, radius.Depth2Go, drChoice.begin()->second.c_str() );
01083 }
01084
01085 ASSERT( radius.drNext > 0. );
01086
01087 if( trace.lgTrace )
01088 {
01089 fprintf( ioQQQ, " radius_next chooses next drad drNext=%.4e; this drad was %.4e\n",
01090 radius.drNext, radius.drad );
01091 }
01092
01093 return 0;
01094 }
01095
01096
01097 STATIC void ContRate(double *freqm,
01098 double *opacm)
01099 {
01100 long int i,
01101 ipHi,
01102 iplow,
01103 limit;
01104 double FreqH,
01105 Freq_nonIonizing,
01106 Opac_Hion,
01107 Opac_nonIonizing,
01108 Rate_max_Hion,
01109 Rate_max_nonIonizing;
01110
01111 DEBUG_ENTRY( "ContRate()" );
01112
01113
01114
01115
01116
01117
01118 Rate_max_nonIonizing = 0.;
01119 Freq_nonIonizing = 0.;
01120 Opac_nonIonizing = 0.;
01121
01122
01123 *opacm = -1.;
01124 *freqm = -1.;
01125
01126
01127
01128 if( dense.lgElmtOn[ipCARBON] )
01129 {
01130
01131 ipHi = Heavy.ipHeavy[ipCARBON][0] - 1;
01132 }
01133 else
01134 {
01135
01136 ipHi = iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH2s].ipIsoLevNIonCon-1;
01137 }
01138
01139
01140 for( i=rfield.ipPlasma; i < ipHi; i++ )
01141 {
01142
01143
01144 if( rfield.anu[i]*rfield.flux[0][i]/rfield.widflx[i]*(opac.opacity_abs[i] -
01145 gv.dstab[i]*dense.gas_phase[ipHYDROGEN]) > Rate_max_nonIonizing )
01146 {
01147 Rate_max_nonIonizing = rfield.anu[i]*rfield.flux[0][i]/rfield.widflx[i]*
01148 (opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN]);
01149 Freq_nonIonizing = rfield.anu[i];
01150 Opac_nonIonizing = opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN];
01151 }
01152 }
01153
01154
01155
01156
01157 if( CoolHeavy.brems_heat_total/thermal.htot > 0.05 )
01158 {
01159
01160
01161 iplow = MAX2(1 , rfield.ipEnergyBremsThin );
01162 }
01163 else
01164 {
01165
01166
01167 iplow = ipHi;
01168 }
01169
01170 iplow = MAX2( rfield.ipPlasma , iplow );
01171
01172
01173 limit = MIN2(rfield.nflux,iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon-1);
01174 for( i=iplow; i < limit; i++ )
01175 {
01176 if( rfield.anu[i]*rfield.flux[0][i]/rfield.widflx[i]*(opac.opacity_abs[i] -
01177 gv.dstab[i]*dense.gas_phase[ipHYDROGEN]) > Rate_max_nonIonizing )
01178 {
01179 Rate_max_nonIonizing = rfield.anu[i]*rfield.flux[0][i]/rfield.widflx[i]*
01180 (opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN]);
01181 Freq_nonIonizing = rfield.anu[i];
01182 Opac_nonIonizing = opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN];
01183 }
01184 }
01185
01186
01187 Rate_max_Hion = 0.;
01188 Opac_Hion = 0.;
01189 FreqH = 0.;
01190
01191
01192 for( i=iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon-1; i < rfield.nflux; i++ )
01193 {
01194 if( rfield.anu[i]*rfield.flux[0][i]/rfield.widflx[i]*(opac.opacity_abs[i] -
01195 gv.dstab[i]*dense.gas_phase[ipHYDROGEN]) > Rate_max_Hion )
01196 {
01197
01198 Rate_max_Hion = rfield.anu[i]*rfield.flux[0][i]/rfield.widflx[i]*
01199 (opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN]);
01200 FreqH = rfield.anu[i];
01201 Opac_Hion = opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN];
01202 }
01203 }
01204
01205
01206 if( Rate_max_nonIonizing < 1e-30 && Opac_Hion > SMALLFLOAT )
01207 {
01208
01209 *opacm = Opac_Hion;
01210 *freqm = FreqH;
01211 }
01212
01213
01214 else if( Opac_Hion > Opac_nonIonizing && Rate_max_Hion/SDIV(Rate_max_nonIonizing) > 1e-10 && Opac_Hion > SMALLFLOAT )
01215 {
01216
01217 *opacm = Opac_Hion;
01218 *freqm = FreqH;
01219 }
01220 else
01221 {
01222
01223 *opacm = Opac_nonIonizing;
01224 *freqm = Freq_nonIonizing;
01225 }
01226
01227 # if 0
01228
01229
01230
01231
01232
01233 if( gv.lgDustOn() && gv.lgGrainPhysicsOn )
01234 {
01235 double fluxGrainPeak = -1.;
01236 long int ipGrainPeak = -1;
01237 long int ipGrainPeak2 = -1;
01238
01239 i = 0;
01240 while( rfield.anu[i] < 0.0912 && i < (rfield.nflux-2) )
01241 {
01242
01243
01244 if( gv.GrainEmission[i]*rfield.anu[i]*opac.opacity_abs[i] > fluxGrainPeak )
01245 {
01246 ipGrainPeak = i;
01247 fluxGrainPeak = gv.GrainEmission[i]*rfield.anu[i]*opac.opacity_abs[i];
01248 }
01249 ++i;
01250 }
01251 ASSERT( fluxGrainPeak>=0. && ipGrainPeak >=0 );
01252
01253
01254
01255
01256
01257
01258
01259 i = ipGrainPeak;
01260
01261
01262
01263
01264 if( fluxGrainPeak > 0. )
01265 {
01266 while( rfield.anu[i] < 0.0912 && i < (rfield.nflux-2) )
01267 {
01268
01269 if( gv.GrainEmission[i]*rfield.anu[i]*opac.opacity_abs[i] > 0.05*fluxGrainPeak )
01270 {
01271
01272 ipGrainPeak2 = i;
01273 }
01274 ++i;
01275 }
01276 ASSERT( ipGrainPeak2>=0 );
01277
01278
01279 if( opac.opacity_abs[ipGrainPeak2] > *opacm )
01280 {
01281
01282 *opacm = opac.opacity_abs[ipGrainPeak2]*10.;
01283 *freqm = rfield.anu[ipGrainPeak2];
01284
01285
01286 }
01287 }
01288 }
01289 # endif
01290
01291 {
01292
01293 enum {DEBUG_LOC=false};
01294 if( DEBUG_LOC )
01295 {
01296 fprintf(ioQQQ,"conratedebug \t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
01297 Rate_max_nonIonizing,Freq_nonIonizing,Opac_nonIonizing,
01298 Rate_max_Hion,FreqH ,Opac_Hion,*freqm,*opacm
01299 );
01300
01301 }
01302 }
01303
01304
01305
01306
01307
01308
01309 ASSERT( *opacm >= 0. && *freqm >= 0. );
01310 return;
01311 }
01312
01313
01314 STATIC void GrainRateDr(double *grfreqm,
01315 double *gropacm)
01316 {
01317 long int i,
01318 iplow;
01319 double xMax;
01320
01321 DEBUG_ENTRY( "GrainRateDr()" );
01322
01323
01324
01325
01326
01327
01328
01329
01330 if( CoolHeavy.brems_heat_total/thermal.htot > 0.05 )
01331 {
01332
01333
01334 iplow = MAX2(1 , rfield.ipEnergyBremsThin );
01335 }
01336 else
01337 {
01338
01339
01340 if( dense.lgElmtOn[ipCARBON] )
01341 {
01342
01343 iplow = Heavy.ipHeavy[ipCARBON][0];
01344 }
01345 else
01346 {
01347
01348 iplow = iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH2s].ipIsoLevNIonCon;
01349 }
01350 }
01351
01352 xMax = -1.;
01353
01354 for( i=iplow-1; i < Heavy.ipHeavy[ipHYDROGEN][0]; i++ )
01355 {
01356 if( rfield.anu[i]*rfield.flux[0][i]/rfield.widflx[i]*opac.opacity_abs[i] > xMax )
01357 {
01358 xMax = rfield.anu[i]*rfield.flux[0][i]/rfield.widflx[i]*
01359 opac.opacity_abs[i];
01360 *grfreqm = rfield.anu[i];
01361 *gropacm = opac.opacity_abs[i];
01362 }
01363 }
01364
01365
01366 if( dense.lgElmtOn[ipHELIUM] )
01367 {
01368 for( i=Heavy.ipHeavy[ipHYDROGEN][0]-1; i < Heavy.ipHeavy[ipHELIUM][1]; i++ )
01369 {
01370 if( rfield.anu[i]*rfield.flux[0][i]/rfield.widflx[i]*opac.opacity_abs[i] > xMax )
01371 {
01372 xMax = rfield.anu[i]*rfield.flux[0][i]/rfield.widflx[i]*
01373 opac.opacity_abs[i];
01374 *grfreqm = rfield.anu[i];
01375 *gropacm = opac.opacity_abs[i];
01376 }
01377 }
01378 }
01379
01380
01381
01382 if( xMax <= 0. )
01383 {
01384 *gropacm = 0.;
01385 *grfreqm = 0.;
01386 }
01387 return;
01388 }