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 "conv.h"
00027 #include "rfield.h"
00028 #include "dynamics.h"
00029 #include "thermal.h"
00030 #include "hmi.h"
00031 #include "coolheavy.h"
00032 #include "timesc.h"
00033 #include "doppvel.h"
00034 #include "stopcalc.h"
00035 #include "colden.h"
00036 #include "phycon.h"
00037 #include "rt.h"
00038 #include "trace.h"
00039 #include "wind.h"
00040 #include "punch.h"
00041 #include "taulines.h"
00042 #include "pressure.h"
00043 #include "iterations.h"
00044 #include "struc.h"
00045 #include "radius.h"
00046
00047 #if 0
00048
00049 static void ChkRate(
00050
00051 long int nelem,
00052
00053 double *dDestRate,
00054
00055 double *DestRateOld,
00056 double *DestRateNew,
00057
00058 long int *istage);
00059 #endif
00060
00061
00062 static void ContRate(double *freqm,
00063 double *opacm);
00064
00065
00066 static void GrainRateDr(double *grfreqm,
00067 double *gropacm);
00068
00069
00070
00071
00072 int radius_next(void)
00073 {
00074 char chLbl[11];
00075 bool lgDoPun;
00076
00077 long int level , ipStrong;
00078
00079 double thickness_total , drThickness , DepthToGo , AV_to_go;
00080 int mole_dr_change;
00081
00082 double DrGrainHeat,
00083 GlobDr,
00084 SaveOHIIoHI,
00085 SpecDr,
00086 Strong,
00087 TauDTau,
00088 TauInwd,
00089 drSolomon_BigH2 ,
00090 dEfrac,
00091 dHdStep,
00092 dRTaue,
00093 dTdStep,
00094 dnew,
00095 drConPres,
00096 drH2_heat_cool = 0. ,
00097 dH2_heat_cool = 0.,
00098 drH2_abund = 0. ,
00099 dr_mole_abund = 0.,
00100 dH2_abund=0.,
00101 dCO_abund=0.,
00102 drDepth,
00103 drDest,
00104 drEfrac,
00105 drFail,
00106 drFluc,
00107 drHMase,
00108 drHe1ovHe2,
00109 drHion,
00110 drInter,
00111 drLeiden_hack ,
00112 drLineHeat,
00113 drSphere,
00114 drTab,
00115 drdHdStep,
00116 drdTdStep,
00117 drThermalFront ,
00118 drmax,
00119 dVeldRad,
00120 fac2,
00121 factor,
00122 freqm,
00123 grfreqm=0.,
00124 gropacm=0.,
00125 hdnew,
00126 opacm,
00127 recom_dr_last_iter ,
00128 OldDR ,
00129 winddr,
00130 x;
00131
00132 double change_heavy_frac=-1. , change_heavy_frac_big , dr_change_heavy ,
00133 frac_change_heavy_frac_big, Efrac_old , Efrac_new;
00134 long int ichange_heavy_nelem=-1 , nelem , ion , ichange_heavy_ion=-1;
00135
00136 static double OHIIoHI,
00137 OldHeat = 0.,
00138 OldTe = 0.,
00139 OlddTdStep = 0.,
00140 OldH2Abund=0.,
00141 OldWindVelocity=0.,
00142 Old_He_atom_ov_ion = 0,
00143 Old_H2_heat_cool;
00144 static long int iteration_last=-1;
00145
00146 static double BigRadius = 1e30;
00147 bool lgFirstCall;
00148
00149 DEBUG_ENTRY( "radius_next()" );
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173 if( (iteration != iteration_last) && (nzone==0) )
00174 {
00175
00176 iteration_last = iteration;
00177 lgFirstCall = true;
00178 }
00179 else
00180 lgFirstCall = false;
00181
00182 if( trace.lgTrace )
00183 {
00184 fprintf( ioQQQ, " radius_next called\n" );
00185 }
00186
00187
00188 OldDR = radius.drad;
00189
00190
00191
00192
00193 if( lgFirstCall )
00194 {
00195 if( dense.xIonDense[ipHYDROGEN][1] > 0. && dense.xIonDense[ipHYDROGEN][0] > 0. )
00196 {
00197 OHIIoHI = dense.xIonDense[ipHYDROGEN][1]/dense.xIonDense[ipHYDROGEN][0];
00198 }
00199 else
00200 {
00201 OHIIoHI = 0.;
00202 }
00203 SaveOHIIoHI = OHIIoHI;
00204 drHion = BigRadius;
00205
00206
00207 }
00208
00209 else if( (dense.xIonDense[ipHYDROGEN][1] > 0. && dense.xIonDense[ipHYDROGEN][0] > 0.) && OHIIoHI > 1e-15 )
00210 {
00211 double atomic_frac = (dense.xIonDense[ipHYDROGEN][1]/dense.xIonDense[ipHYDROGEN][0]);
00212
00213
00214 x = 1. - atomic_frac /OHIIoHI;
00215 if( atomic_frac > 0.05 && atomic_frac < 0.9 )
00216 {
00217
00218
00219
00220 drHion = radius.drad*MAX2( 0.2 , 0.05/MAX2(1e-10,x) );
00221 }
00222 else if( x > 0. )
00223 {
00224
00225
00226 drHion = radius.drad*MAX2( 0.2 , 0.2/MAX2(1e-10,x) );
00227 }
00228 else
00229 {
00230 drHion = BigRadius;
00231 }
00232 SaveOHIIoHI = OHIIoHI;
00233 OHIIoHI = dense.xIonDense[ipHYDROGEN][1]/dense.xIonDense[ipHYDROGEN][0];
00234 }
00235
00236 else
00237 {
00238 SaveOHIIoHI = OHIIoHI;
00239 if( dense.xIonDense[ipHYDROGEN][1] > 0. && dense.xIonDense[ipHYDROGEN][0] > 0. )
00240 {
00241 OHIIoHI = dense.xIonDense[ipHYDROGEN][1]/dense.xIonDense[ipHYDROGEN][0];
00242 }
00243 else
00244 {
00245 OHIIoHI = 0.;
00246 }
00247 drHion = BigRadius;
00248 }
00249
00250
00251
00252 if( rt.dTauMase < -0.01 )
00253 {
00254
00255
00256 drHMase = radius.drad*MAX2(0.1,-0.2/rt.dTauMase);
00257 }
00258 else
00259 {
00260 drHMase = BigRadius;
00261 }
00262
00263
00264 Old_He_atom_ov_ion = 0.;
00265 drHe1ovHe2 = BigRadius;
00266
00267
00268
00269
00270 dr_change_heavy = BigRadius;
00271 change_heavy_frac_big = -1.;
00272 frac_change_heavy_frac_big = -1.;
00273
00274 # define CHANGE_ION_HEAV 0.2f
00275 # define CHANGE_ION_HHE 0.15f
00276 if( nzone > 4 )
00277 {
00278
00279 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
00280 {
00281
00282 if( dense.lgElmtOn[nelem] )
00283 {
00284 float change;
00285
00286
00287 float frac_limit;
00288 if( nelem<=ipHELIUM )
00289 {
00290 frac_limit = 1e-4f;
00291 change = CHANGE_ION_HHE;
00292 }
00293 else
00294 {
00295
00296
00297
00298 frac_limit = struc.dr_ionfrac_limit/2.f;
00299 change = CHANGE_ION_HEAV;
00300 }
00301
00302 for( ion=0; ion<=nelem+1; ++ion )
00303 {
00304 float abundnew = dense.xIonDense[nelem][ion]/SDIV( dense.gas_phase[nelem]);
00305 float abundold = struc.xIonDense[nelem][ion][nzone-3]/SDIV( struc.gas_phase[nelem][nzone-3]);
00306 if( abundnew > frac_limit && abundold > frac_limit )
00307 {
00308
00309
00310
00311 float abundolder = struc.xIonDense[nelem][ion][nzone-4]/SDIV( struc.gas_phase[nelem][nzone-4]);
00312 float abundoldest = struc.xIonDense[nelem][ion][nzone-5]/SDIV( struc.gas_phase[nelem][nzone-5]);
00313
00314
00315
00316 change_heavy_frac = fabs(abundnew-abundold)/MIN2(abundold,abundnew);
00317
00318 if( (change_heavy_frac > change) && (change_heavy_frac > change_heavy_frac_big) &&
00319
00320
00321 ((abundnew-abundold)>0.) &&
00322 ((abundold-abundolder)>0.) &&
00323 ((abundolder-abundoldest)>0.) )
00324 {
00325 ichange_heavy_nelem = nelem;
00326 ichange_heavy_ion = ion;
00327 change_heavy_frac_big = change_heavy_frac;
00328 frac_change_heavy_frac_big = abundnew;
00329
00330
00331
00332
00333
00334 dr_change_heavy = radius.drad * MAX2(0.25, change / change_heavy_frac );
00335 }
00336 }
00337 }
00338 }
00339 }
00340 }
00341
00342
00343
00344 if(!co.lgUMISTrates)
00345 {
00346
00347
00348
00349 drLeiden_hack = MAX2( 0.02 , 0.05*rfield.extin_mag_V_point) / SDIV( geometry.FillFac * rfield.opac_mag_V_point );
00350 }
00351 else
00352 {
00353 drLeiden_hack = BigRadius;
00354 }
00355
00356
00357
00358
00359
00360 if( nzone <= 1 || thermal.lgTSetOn )
00361 {
00362 drdHdStep = BigRadius;
00363 dHdStep = FLT_MAX;
00364 }
00365 else
00366 {
00367 dHdStep = fabs(thermal.htot-OldHeat)/thermal.htot;
00368 if( dHdStep > 0. )
00369 {
00370 if( dense.gas_phase[ipHYDROGEN] >= 1e13 )
00371 {
00372
00373 drdHdStep = radius.drad*MAX2(0.8,0.075/dHdStep);
00374 }
00375 else if( dense.gas_phase[ipHYDROGEN] >= 1e11 )
00376 {
00377
00378 drdHdStep = radius.drad*MAX2(0.8,0.100/dHdStep);
00379 }
00380 else
00381 {
00382
00383
00384
00385
00386 drdHdStep = radius.drad*MAX2(0.8,0.15/dHdStep);
00387 }
00388 }
00389 else
00390 {
00391 drdHdStep = BigRadius;
00392 }
00393 }
00394 OldHeat = thermal.htot;
00395
00396
00397
00398 if( strcmp(dense.chDenseLaw,"CPRE") == 0 && pressure.lgContRadPresOn )
00399 {
00400 if( nzone > 2 && pressure.pinzon > 0. )
00401 {
00402
00403
00404
00405
00406
00407 drConPres = 0.05*pressure.PresTotlCurr/(wind.AccelTot*
00408 dense.xMassDensity*geometry.FillFac);
00409 }
00410 else
00411 {
00412 drConPres = BigRadius;
00413 }
00414 }
00415 else
00416 {
00417 drConPres = BigRadius;
00418 }
00419
00420
00421
00422
00423 if( nzone <= 1 )
00424 {
00425 drdTdStep = BigRadius;
00426 dTdStep = FLT_MAX;
00427 OlddTdStep = dTdStep;
00428 }
00429 else
00430 {
00431
00432 dTdStep = (phycon.te-OldTe)/phycon.te;
00433
00434
00435
00436
00437 x = conv.HeatCoolRelErrorAllowed*2.;
00438 x = MAX2( 0.01 , x );
00439 x = MIN2( 0.05 , x );
00440
00441 x = 0.03;
00442
00443
00444
00445 if( dTdStep*OlddTdStep > 0. )
00446 {
00447
00448
00449
00450
00451
00452 double absdt = fabs(dTdStep);
00453 drdTdStep = radius.drad*MAX2(0.5,x/absdt);
00454 }
00455 else
00456 {
00457 drdTdStep = BigRadius;
00458 }
00459 }
00460 OlddTdStep = dTdStep;
00461 OldTe = phycon.te;
00462
00463
00464
00465
00466
00467
00468 if( !thermal.lgTSetOn && !dynamics.lgRecom )
00469 {
00470
00471
00472
00473 ContRate(&freqm,&opacm);
00474
00475
00476
00477
00478
00479
00480 drInter = 0.3/MAX2(1e-30,opacm*geometry.FillFac*geometry.DirectionalCosin);
00481 }
00482 else
00483 {
00484 drInter = BigRadius;
00485 freqm = 0.;
00486 opacm = 1.;
00487 }
00488
00489 # if 0
00490
00491
00492 ChkRate(ipCARBON,&dDRCarbon,&DestOldCarb, &DestNewCarb,&icarstag);
00493 ChkRate(ipNITROGEN,&dDRNitrogen,&DestOldNit, &DestNewNit,&initstag);
00494 ChkRate(ipOXYGEN,&dDROxygen,&DestOldOxy, &DestNewOxy,&ioxystag);
00495 ChkRate(ipIRON,&dDRIron,&DestOldIron, &DestNewIron,&iironstag);
00496
00497 dDestRate = MAX4(dDROxygen,dDRIron,dDRCarbon,dDRNitrogen);
00498
00499 if( dDRCarbon == dDestRate )
00500 {
00501 dDestRate = dDRCarbon;
00502 DestOld = DestOldCarb;
00503 DestNew = DestNewCarb;
00504 istage = icarstag;
00505 strcpy( chDestAtom, "Carbon " );
00506 }
00507
00508 else if( dDRNitrogen == dDestRate )
00509 {
00510 dDestRate = dDRNitrogen;
00511 DestOld = DestOldNit;
00512 DestNew = DestNewNit;
00513 istage = initstag;
00514 strcpy( chDestAtom, "Nitrogen" );
00515 }
00516
00517 else if( dDROxygen == dDestRate )
00518 {
00519 dDestRate = dDROxygen;
00520 DestOld = DestOldOxy;
00521 DestNew = DestNewOxy;
00522 strcpy( chDestAtom, "Oxygen " );
00523 istage = ioxystag;
00524 }
00525
00526 else if( dDRIron == dDestRate )
00527 {
00528 dDestRate = dDRIron;
00529 DestOld = DestOldIron;
00530 DestNew = DestNewIron;
00531 istage = iironstag;
00532 strcpy( chDestAtom, "Iron " );
00533 }
00534
00535 else
00536 {
00537 fprintf( ioQQQ, " insanity following ChkRate\n" );
00538 ShowMe();
00539 puts( "[Stop in radius_next]" );
00540 cdEXIT(EXIT_FAILURE);
00541 }
00542
00543
00544 if( dDestRate > 0. )
00545 {
00546
00547
00548
00549
00550
00551
00552
00553
00554
00555
00556
00557
00558 drDest = radius.drad*MAX2(0.5,0.20/dDestRate);
00559
00560 }
00561 else
00562 {
00563 drDest = BigRadius;
00564 }
00565 # endif
00566 drDest = BigRadius;
00567
00568
00569
00570
00571
00572
00573
00574 if( wind.windv!=0. && !pressure.lgSonicPoint && !pressure.lgStrongDLimbo )
00575 {
00576 double v = fabs(wind.windv);
00577
00578 dVeldRad = fabs(wind.windv-OldWindVelocity)/
00579 MAX2(v,0.1*timesc.sound_speed_isothermal)/radius.drad;
00580
00581 if( 1.1*dVeldRad*radius.drad > 0.03 )
00582 {
00583
00584 winddr = 0.03/dVeldRad;
00585 }
00586 else
00587 {
00588 winddr = 1.1*radius.drad;
00589 }
00590
00591
00592
00593 if( dynamics.lgAdvection )
00594 {
00595
00596 winddr = MIN2( winddr , dynamics.dRad / 20. );
00597
00598
00599 dVeldRad = dynamics.dRad;
00600 }
00601 }
00602 else
00603 {
00604 winddr = BigRadius;
00605 dVeldRad = 0.;
00606 }
00607
00608 OldWindVelocity = wind.windv;
00609
00610
00611
00612 if( strcmp(dense.chDenseLaw,"GLOB") == 0 )
00613 {
00614 # define DNGLOB 0.10
00615 if( radius.glbdst < 0. )
00616 {
00617 fprintf( ioQQQ, " Globule distance is negative, internal overflow has occured, sorry.\n" );
00618 fprintf( ioQQQ, " This is routine radius_next, GLBDST is%10.2e\n",
00619 radius.glbdst );
00620 puts( "[Stop in radius_next]" );
00621 cdEXIT(EXIT_FAILURE);
00622 }
00623 factor = radius.glbden*pow(radius.glbrad/radius.glbdst,radius.glbpow);
00624 fac2 = radius.glbden*pow(radius.glbrad/(radius.glbdst - (float)radius.drad),radius.glbpow);
00625 if( fac2/factor > 1. + DNGLOB )
00626 {
00627
00628 GlobDr = radius.drad*DNGLOB/(fac2/factor - 1.);
00629 }
00630 else
00631 {
00632 GlobDr = BigRadius;
00633 }
00634
00635 GlobDr = MIN2(GlobDr,radius.glbdst/20.);
00636 }
00637 else
00638 {
00639 GlobDr = BigRadius;
00640 }
00641
00642
00643 hdnew = 0.;
00644 if( strncmp( dense.chDenseLaw , "DLW" , 3) == 0 )
00645 {
00646
00647 if( strcmp(dense.chDenseLaw,"DLW1") == 0 )
00648 {
00649 hdnew = dense_fabden(radius.Radius+radius.drad,radius.depth+
00650 radius.drad);
00651 }
00652 else if( strcmp(dense.chDenseLaw,"DLW2") == 0 )
00653 {
00654 hdnew = dense_tabden(radius.Radius+radius.drad,radius.depth+
00655 radius.drad);
00656 }
00657 else
00658 {
00659 fprintf( ioQQQ, " dlw insanity in radius_next\n" );
00660 puts( "[Stop in radius_next]" );
00661 cdEXIT(EXIT_FAILURE);
00662 }
00663 drTab = fabs(hdnew-dense.gas_phase[ipHYDROGEN])/MAX2(hdnew,dense.gas_phase[ipHYDROGEN]);
00664 drTab = radius.drad*MAX2(0.2,0.10/MAX2(0.01,drTab));
00665 }
00666 else
00667 {
00668 drTab = BigRadius;
00669 }
00670
00671
00672
00673 if( strcmp(dense.chDenseLaw,"DLW1") == 0 )
00674 {
00675 dnew = fabs(dense_fabden(radius.Radius+radius.drad,radius.depth+
00676 radius.drad)/dense.gas_phase[ipHYDROGEN]-1.);
00677
00678 if( dnew == 0. )
00679 {
00680 SpecDr = radius.drad*3.;
00681 }
00682 else if( dnew/DNGLOB > 1.0 )
00683 {
00684 SpecDr = radius.drad/(dnew/DNGLOB);
00685 }
00686 else
00687 {
00688 SpecDr = BigRadius;
00689 }
00690 }
00691 else
00692 {
00693 SpecDr = BigRadius;
00694 }
00695
00696
00697
00698
00699
00700 if( thermal.heating[0][13]/thermal.htot > 0.2 )
00701 {
00702
00703
00704 GrainRateDr(&grfreqm,&gropacm);
00705 DrGrainHeat = 1.0/MAX2(1e-30,gropacm*geometry.FillFac*geometry.DirectionalCosin);
00706 }
00707 else
00708 {
00709 gropacm = 0.;
00710 grfreqm = 0.;
00711 DrGrainHeat = BigRadius;
00712 }
00713
00714
00715
00716
00717 if( thermal.heating[0][22]/thermal.htot > 0.2 )
00718 {
00719 FndLineHt(&level,&ipStrong,&Strong);
00720 if( Strong/thermal.htot > 0.1 )
00721 {
00722 if( level == 1 )
00723 {
00724
00725
00726
00727 TauInwd = TauLines[ipStrong].TauIn;
00728 TauDTau = TauLines[ipStrong].PopOpc * TauLines[ipStrong].opacity /
00729 DoppVel.doppler[TauLines[ipStrong].nelem-1];
00730 }
00731 else if( level == 2 )
00732 {
00733
00734
00735
00736
00737 TauInwd = TauLine2[ipStrong].TauIn;
00738 TauDTau = TauLine2[ipStrong].PopOpc * TauLine2[ipStrong].opacity /
00739 DoppVel.doppler[TauLine2[ipStrong].nelem-1];
00740 }
00741 else if( level == 3 )
00742 {
00743
00744
00745
00746
00747 TauInwd = C12O16Rotate[ipStrong].TauIn;
00748 TauDTau = C12O16Rotate[ipStrong].PopOpc * C12O16Rotate[ipStrong].opacity /
00749 DoppVel.doppler[C12O16Rotate[ipStrong].nelem-1];
00750 }
00751 else if( level == 4 )
00752 {
00753
00754
00755
00756
00757 TauInwd = C13O16Rotate[ipStrong].TauIn;
00758 TauDTau = C13O16Rotate[ipStrong].PopOpc * C13O16Rotate[ipStrong].opacity /
00759 DoppVel.doppler[C13O16Rotate[ipStrong].nelem-1];
00760 }
00761 else if( level == 5 )
00762 {
00763
00764
00765 TauInwd = HFLines[ipStrong].TauIn;
00766 TauDTau = HFLines[ipStrong].PopOpc * HFLines[ipStrong].opacity /
00767 DoppVel.doppler[HFLines[ipStrong].nelem-1];
00768 }
00769 else
00770 {
00771
00772 fprintf( ioQQQ, " PROBLEM radius_next Strong line heating set, but not level.\n" );
00773 TotalInsanity();
00774 }
00775
00776
00777
00778 if( TauDTau > 0. )
00779 {
00780 drLineHeat = MAX2(1.,TauInwd)*0.4/TauDTau;
00781 }
00782 else
00783 {
00784 drLineHeat = BigRadius;
00785 }
00786 }
00787 else
00788 {
00789 TauInwd = 0.;
00790 drLineHeat = BigRadius;
00791 ipStrong = 0;
00792 Strong = 0.;
00793 }
00794 }
00795 else
00796 {
00797 TauInwd = 0.;
00798 drLineHeat = BigRadius;
00799 ipStrong = 0;
00800 level = 0;
00801 Strong = 0.;
00802 }
00803
00804
00805
00806 drH2_heat_cool = BigRadius;
00807 if( lgFirstCall )
00808 {
00809 Old_H2_heat_cool = 0.;
00810 }
00811 else if( !thermal.lgTSetOn )
00812 {
00813
00814
00815 double H2_heat_cool = (fabs(hmi.HeatH2Dexc_used)+fabs(hmi.HeatH2Dish_used)) / thermal.htot;
00816 if( H2_heat_cool > 0.1 )
00817 {
00818 dH2_heat_cool = fabs( H2_heat_cool - Old_H2_heat_cool );
00819 dH2_heat_cool = SDIV(dH2_heat_cool);
00820
00821
00822
00823 drH2_heat_cool = radius.drad*MAX2(0.5,0.05/dH2_heat_cool);
00824 }
00825 else
00826 {
00827 drH2_heat_cool = BigRadius;
00828 }
00829 }
00830 Old_H2_heat_cool = (fabs(hmi.HeatH2Dexc_used)+fabs(hmi.HeatH2Dish_used)) / thermal.htot;
00831
00832
00833
00834
00835 drH2_abund = BigRadius;
00836 dr_mole_abund = BigRadius;
00837 mole_dr_change = -1;
00838 if( nzone>=4 )
00839 {
00840
00841 double Old_H2_abund;
00842
00843 int i;
00844 Old_H2_abund = struc.H2_molec[ipMH2g][nzone-3] + struc.H2_molec[ipMH2s][nzone-3];
00845
00846
00847
00848
00849
00850
00851
00852
00853
00854
00855 if( 2.*hmi.H2_total/dense.gas_phase[ipHYDROGEN] > 1e-8 )
00856 {
00857 double fac = 0.2;
00858
00859 dH2_abund = 2.*fabs( hmi.H2_total - Old_H2_abund ) / hmi.H2_total;
00860
00861
00862
00863
00864
00865
00866
00867
00868
00869
00870
00871 dH2_abund = SDIV(dH2_abund);
00872 drH2_abund = radius.drad*MAX2(0.2,fac/dH2_abund );
00873 }
00874 else
00875 {
00876 drH2_abund = BigRadius;
00877 }
00878
00879
00880
00881 dr_mole_abund = BigRadius;
00882
00883
00884 for(i=0; i<mole.num_comole_calc; ++i)
00885 {
00886 float abund,
00887 abund_limit;
00888 if( !dense.lgElmtOn[COmole[i]->nelem_hevmol] || COmole[i]->n_nuclei <= 1 )
00889 continue;
00890
00891
00892
00893
00894 abund = COmole[i]->hevmol/dense.gas_phase[COmole[i]->nelem_hevmol];
00895
00896
00897
00898
00899 if( COmole[i]->lgGas_Phase )
00900 {
00901 abund_limit = 1e-3f;
00902 }
00903 else
00904 {
00905
00906
00907 abund_limit = 1e-5f;
00908 }
00909
00910 if( abund > abund_limit )
00911 {
00912 double drmole_one, relative_change, relative_change_denom;
00913
00914
00915
00916 if( struc.CO_molec[i][nzone-3]>SMALLFLOAT )
00917 {
00918 relative_change_denom = MIN2( COmole[i]->hevmol , struc.CO_molec[i][nzone-3] );
00919 }
00920 else
00921 {
00922 relative_change_denom = COmole[i]->hevmol;
00923 }
00924
00925 relative_change =
00926
00927
00928 fabs( COmole[i]->hevmol - struc.CO_molec[i][nzone-3] ) / relative_change_denom;
00929
00930
00931
00932
00933
00934
00935
00936
00937
00938 relative_change = SDIV(relative_change);
00939
00940
00941 if( relative_change > 1. )
00942 relative_change = 1./relative_change;
00943
00944
00945
00946
00947 drmole_one = radius.drad*MAX2(0.6,0.05/relative_change );
00948
00949 if( drmole_one < dr_mole_abund )
00950 {
00951
00952 dr_mole_abund = drmole_one;
00953 mole_dr_change = i;
00954 dCO_abund = relative_change;
00955 }
00956 }
00957 }
00958 }
00959
00960
00961 drSolomon_BigH2 = H2_DR();
00962
00963
00964
00965
00966
00967 if( nzone < 5 )
00968 {
00969
00970
00971
00972 drmax = 4.*radius.drad;
00973 }
00974 else
00975 {
00976 drmax = 1.3*radius.drad;
00977 }
00978
00979
00980
00981 # if 0
00982
00983 dr_ne_oscil = BigRadius;
00984 dr_te_oscil = BigRadius;
00985 if( nzone >= 11 )
00986 {
00987
00988
00989
00990 float errorHC = POW2(conv.HeatCoolRelErrorAllowed);
00991 float errorNe = (float)POW2(conv.EdenErrorAllowed );
00992 limit = nzone -2;
00993 ASSERT( limit < struc.nzlim );
00994 for( k=nzone - 10; k < limit; k++ )
00995 {
00996
00997
00998 if( (struc.testr[k-1] - struc.testr[k])/struc.testr[k]*
00999 (struc.testr[k] - struc.testr[k+1])/struc.testr[k] <
01000 -errorHC )
01001 {
01002 dr_te_oscil = radius.drad;
01003 }
01004
01005 if( (struc.ednstr[k-1] - struc.ednstr[k])/struc.ednstr[k]*
01006 (struc.ednstr[k] - struc.ednstr[k+1])/struc.ednstr[k] <
01007 -errorNe )
01008 {
01009
01010
01011 dr_ne_oscil = radius.drad;
01012 }
01013 }
01014 }
01015
01016
01017 dr_ne_oscil = BigRadius;
01018 dr_te_oscil = BigRadius;
01019 # endif
01020
01021
01022
01023
01024
01025 if( !conv.lgConvTemp )
01026 {
01027 drFail = radius.drad/2.;
01028 }
01029 else
01030 {
01031 drFail = BigRadius;
01032 }
01033
01034
01035
01036
01037
01038 recom_dr_last_iter = BigRadius;
01039 if( dynamics.lgStatic && dynamics.lgRecom )
01040 {
01041 static long int nzone_recom = -1;
01042 if( nzone==0 )
01043 {
01044
01045 nzone_recom = 0;
01046 }
01047 else if( radius.depth < struc.depth_last[struc.nzone_last-1] )
01048 {
01049
01050
01051
01052 while( struc.depth_last[nzone_recom] < radius.depth &&
01053 nzone_recom < struc.nzone_last-1 )
01054 ++nzone_recom;
01055 recom_dr_last_iter = struc.drad_last[nzone_recom]*3.;
01056 }
01057 else
01058 {
01059
01060 recom_dr_last_iter = BigRadius;
01061 }
01062 }
01063
01064
01065
01066
01067
01068 if( nzone > 2 )
01069 {
01070
01071 Efrac_old = struc.ednstr[nzone-3]/struc.gas_phase[ipHYDROGEN][nzone-3];
01072 Efrac_new = dense.eden/dense.gas_phase[ipHYDROGEN];
01073 dEfrac = fabs(Efrac_old-Efrac_new) / Efrac_new;
01074
01075 if( dEfrac > SMALLFLOAT )
01076 {
01077 double fac = 0.04;
01078
01079
01080
01081 if( dense.eden_from_metals > 0.5 )
01082 {
01083
01084
01085 fac = 0.04;
01086 }
01087
01088
01089
01090 else if( iso.RecomCollisFrac[ipH_LIKE][ipHYDROGEN] > 0.8 &&
01091 dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN]>0.1 &&
01092 dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN]<0.8 )
01093
01094 {
01095 fac = 0.02;
01096 }
01097
01098 drEfrac = radius.drad*MAX2(0.1,fac/dEfrac);
01099 }
01100 else
01101 {
01102 drEfrac = BigRadius;
01103 }
01104 }
01105 else
01106 {
01107 dEfrac = 0.;
01108 drEfrac = BigRadius;
01109 Efrac_old = 0.;
01110 Efrac_new = 0.;
01111 }
01112
01113
01114
01115
01116 if( nzone > 20 )
01117 {
01118
01119
01120 drDepth = radius.depth/10.;
01121 }
01122 else
01123 {
01124 drDepth = BigRadius;
01125 }
01126
01127
01128
01129 thickness_total = BigRadius;
01130 DepthToGo = BigRadius;
01131 if( StopCalc.HColStop < 5e29 )
01132 {
01133 double coleff = SDIV( dense.gas_phase[ipHYDROGEN]*geometry.FillFac);
01134 DepthToGo = MIN2(DepthToGo ,
01135 (StopCalc.HColStop-colden.colden[ipCOL_HTOT]) / coleff );
01136
01137 thickness_total = MIN2(thickness_total , StopCalc.HColStop / coleff );
01138 }
01139
01140 if( StopCalc.colpls < 5e29 )
01141 {
01142 double coleff = (double)SDIV( (dense.xIonDense[ipHYDROGEN][1])*geometry.FillFac);
01143 DepthToGo = MIN2(DepthToGo ,
01144 (StopCalc.colpls-colden.colden[ipCOL_Hp]) / coleff );
01145 thickness_total = MIN2(thickness_total , StopCalc.colpls / coleff );
01146 }
01147
01148 if( StopCalc.col_h2 < 5e29 )
01149 {
01150
01151 double coleff = (double)SDIV( hmi.H2_total*geometry.FillFac);
01152 DepthToGo = MIN2(DepthToGo ,
01153 (StopCalc.col_h2-colden.colden[ipCOL_H2g]-colden.colden[ipCOL_H2s]) / coleff );
01154 thickness_total = MIN2(thickness_total , StopCalc.col_h2 / coleff );
01155 }
01156
01157 if( StopCalc.col_h2_nut < 5e29 )
01158 {
01159
01160 double coleff = (double)SDIV( (2*hmi.H2_total+dense.xIonDense[ipHYDROGEN][1])*geometry.FillFac);
01161 DepthToGo = MIN2(DepthToGo ,
01162 (StopCalc.col_h2_nut-(2*(colden.colden[ipCOL_H2g]+colden.colden[ipCOL_H2s])+dense.xIonDense[ipHYDROGEN][1])) / coleff );
01163 thickness_total = MIN2(thickness_total , StopCalc.col_h2_nut / coleff );
01164 }
01165
01166 if( StopCalc.col_H0_ov_Tspin < 5e29 )
01167 {
01168
01169 double coleff = (double)SDIV( dense.xIonDense[ipHYDROGEN][0] / hyperfine.Tspin21cm*geometry.FillFac );
01170 DepthToGo = MIN2(DepthToGo ,
01171 (StopCalc.col_H0_ov_Tspin - colden.H0_ov_Tspin) / coleff );
01172 thickness_total = MIN2(thickness_total , StopCalc.col_H0_ov_Tspin / coleff );
01173 }
01174
01175 if( StopCalc.col_monoxco < 5e29 )
01176 {
01177
01178 double coleff = (double)SDIV( (findspecies("CO")->hevmol)*geometry.FillFac);
01179 DepthToGo = MIN2(DepthToGo ,
01180 (StopCalc.col_monoxco-findspecies("CO")->hevcol) / coleff );
01181 thickness_total = MIN2(thickness_total , StopCalc.col_monoxco / coleff );
01182 }
01183
01184 if( StopCalc.colnut < 5e29 )
01185 {
01186 double coleff = (double)SDIV( (dense.xIonDense[ipHYDROGEN][0])*geometry.FillFac);
01187 DepthToGo = MIN2(DepthToGo ,
01188 (StopCalc.colnut - colden.colden[ipCOL_H0]) / coleff );
01189 thickness_total = MIN2(thickness_total , StopCalc.colnut / coleff );
01190 }
01191
01192
01193 if( radius.router[iteration-1] < 5e29 )
01194 {
01195 thickness_total = MIN2(thickness_total , radius.router[iteration-1] );
01196 DepthToGo = MIN2(DepthToGo ,
01197 radius.router[iteration-1] - radius.depth );
01198 }
01199
01200
01201 if( StopCalc.iptnu != rfield.nupper )
01202 {
01203
01204 double dt = SDIV(opac.opacity_abs[StopCalc.iptnu-1]*geometry.FillFac);
01205 DepthToGo = MIN2(DepthToGo ,
01206 (StopCalc.tauend-opac.TauAbsGeo[0][StopCalc.iptnu-1])/dt);
01207 }
01208
01209
01210
01211 AV_to_go = BigRadius;
01212 if( rfield.opac_mag_V_extended > SMALLFLOAT && rfield.opac_mag_V_point>SMALLFLOAT )
01213 {
01214
01215
01216 double ave = log10(StopCalc.AV_extended - rfield.extin_mag_V_extended) -
01217 log10(rfield.opac_mag_V_extended);
01218 double avp = log10(StopCalc.AV_point - rfield.extin_mag_V_point) -
01219 log10(rfield.opac_mag_V_point);
01220 AV_to_go = MIN2( ave , avp );
01221 if( AV_to_go < 37. )
01222 {
01223 AV_to_go = pow(10., AV_to_go );
01224
01225
01226 AV_to_go *= 1.0001;
01227 }
01228 else
01229 AV_to_go = BigRadius;
01230
01231 }
01232
01233
01234
01235 if( DepthToGo <= 0. )
01236 {
01237 TotalInsanity();
01238 }
01239 else if( DepthToGo < BigRadius )
01240 {
01241
01242
01243
01244 drThickness = MIN2( thickness_total/10. , DepthToGo );
01245 }
01246 else
01247 {
01248 drThickness = BigRadius;
01249 }
01250
01251
01252
01253
01254
01255
01256
01257 drSphere = radius.Radius*0.04;
01258
01259
01260
01261 dRTaue = radius.drChange/(dense.eden*6.65e-25*geometry.FillFac);
01262
01263
01264 if( thermal.lgTSetOn )
01265 dRTaue *= 3.;
01266
01267
01268 if( dense.flong != 0. )
01269 {
01270 drFluc = 0.628/2./dense.flong;
01271
01272
01273 drFluc /= 2.;
01274 }
01275 else
01276 {
01277 drFluc = BigRadius;
01278 }
01279
01280
01281
01282
01283 if( strcmp(dense.chDenseLaw,"SINE") == 0 && dense.lgDenFlucOn )
01284 {
01285 drdHdStep = BigRadius;
01286 }
01287
01288
01289
01290
01291
01292 radius.drNext = MIN2( MIN4( drmax, drInter, drLineHeat, winddr ),
01293 MIN4( drFluc, GlobDr, DrGrainHeat, dr_change_heavy ) );
01294 radius.drNext = MIN2( MIN3( radius.drNext, SpecDr, drFail ),
01295 MIN3( drSphere, radius.sdrmax, dRTaue ) );
01296 radius.drNext = MIN3( MIN3( radius.drNext, drDest, drdTdStep ),
01297 MIN3( drdHdStep, drConPres, drTab ),
01298 MIN3( drSolomon_BigH2, drLeiden_hack, recom_dr_last_iter ) );
01299 radius.drNext = MIN3( MIN4( radius.drNext, drHion, drDepth, dr_mole_abund ),
01300 MIN4( AV_to_go, drEfrac, drHMase, drThickness ),
01301 MIN3( drHe1ovHe2, drH2_heat_cool, drH2_abund ) );
01302
01303
01304
01305
01306
01307
01308
01309
01310 if( nzone <= 1 && radius.drNext > OldDR )
01311 {
01312 radius.drNext = OldDR;
01313 }
01314
01315
01316 if( radius.drNext < radius.sdrmin )
01317 {
01318 radius.drNext = radius.sdrmin;
01319 }
01320
01323
01324
01325 if( !conv.lgConvPres || !conv.lgConvTemp )
01326 {
01327 radius.drNext = radius.drad;
01328 }
01329
01330
01331
01332
01333
01334
01335
01336 drThermalFront = BigRadius;
01337
01338 if( nzone >=5 && !dynamics.lgStatic )
01339 {
01340
01341
01342 if( phycon.te > 200. && phycon.te < 3000. &&
01343
01344
01345 (struc.testr[nzone-3] - phycon.te)/phycon.te > 0.02 &&
01346 (struc.testr[nzone-4] - phycon.te)/phycon.te > 0.02 &&
01347 (struc.testr[nzone-5] - phycon.te)/phycon.te > 0.02 )
01348 {
01349
01350
01351 drThermalFront = radius.drad * 0.91;
01352 radius.drNext = drThermalFront;
01353 }
01354 }
01355
01356
01357 if( radius.drNext <= 0. )
01358 {
01359 fprintf( ioQQQ, " radius_next finds insane drNext:%10.2e\n",
01360 radius.drNext );
01361 fprintf( ioQQQ, " all drs follow:%10.2e%10.2e%10.2e%10.2e\n all drs follow:%10.2e%10.2e%10.2e%10.2e%10.2e\n all drs follow:%10.2e%10.2e%10.2e%10.2e\n",
01362 drmax, drInter, drLineHeat, winddr, drFluc, GlobDr, SpecDr,
01363 drFail, drSphere, radius.sdrmax, dRTaue,
01364 OldH2Abund, drDepth );
01365 puts( "[Stop in radius_next]" );
01366 cdEXIT(EXIT_FAILURE);
01367 }
01368
01369
01370 if( radius.drNext == drHMase )
01371 {
01372 rt.lgMaserSetDR = true;
01373 }
01374
01375
01376
01377
01378 if( punch.lgDROn )
01379 {
01380 if( punch.lgDRPLst )
01381 {
01382
01383 if( iterations.lgLastIt )
01384 {
01385 lgDoPun = true;
01386 }
01387 else
01388 {
01389 lgDoPun = false;
01390 }
01391 }
01392 else
01393 {
01394 lgDoPun = true;
01395 }
01396 }
01397 else
01398 {
01399 lgDoPun = false;
01400 }
01401 if( (trace.lgTrace && trace.lgDrBug) || lgDoPun )
01402 {
01403 if( !conv.lgConvTemp && nzone > 0 )
01404 {
01405 fprintf( punch.ipDRout, "#>>>> A temperature failure occured.\n" );
01406 }
01407 if( !conv.lgConvPres && nzone > 0 )
01408 {
01409 fprintf( punch.ipDRout, "#>>>> A pressure failure occured.\n" );
01410 }
01411
01412
01413 fprintf( punch.ipDRout , "%ld\t%.5e\t%.3e\t%.3e\t", nzone, radius.depth, radius.drNext, radius.Depth2Go );
01414
01415
01416 if( radius.drNext == drLineHeat )
01417 {
01418 if( level == 1 )
01419 {
01420 strcpy( chLbl, chLineLbl(&TauLines[ipStrong]) );
01421 fprintf( punch.ipDRout, "level 1 line heating,%10.10s TauIn%10.2e%10.2e%10.2e\n",
01422 chLbl, TauInwd, TauLines[ipStrong].pump,
01423 TauLines[ipStrong].Pesc );
01424 }
01425 else if( level == 2 )
01426 {
01427 strcpy( chLbl, chLineLbl(&TauLine2[ipStrong]));
01428 fprintf( punch.ipDRout, "level 2 line heating,%10.10s TauIn%10.2e%10.2e%10.2e\n",
01429 chLbl, TauInwd, TauLine2[ipStrong].pump,
01430 TauLine2[ipStrong].Pesc );
01431 }
01432 else
01433 {
01434 fprintf( ioQQQ, " insanity pr line heat\n" );
01435 puts( "[Stop in radius_next]" );
01436 cdEXIT(EXIT_FAILURE);
01437 }
01438 }
01439
01440 else if( radius.drNext == drDepth )
01441 {
01442 fprintf( punch.ipDRout, "relative depth\n");
01443 }
01444
01445 else if( radius.drNext == drThermalFront )
01446 {
01447 fprintf( punch.ipDRout, "thermal front logic\n");
01448 }
01449
01450 else if( radius.drNext == dr_change_heavy )
01451 {
01452 fprintf( punch.ipDRout,
01453 "change in ionization, element %s ion %li rel change %.2e ion frac %.2e\n",
01454 elementnames.chElementName[ichange_heavy_nelem],
01455 ichange_heavy_ion ,
01456 change_heavy_frac_big ,
01457 frac_change_heavy_frac_big);
01458 }
01459
01460 else if( radius.drNext == drThickness )
01461 {
01462 fprintf( punch.ipDRout, "depth to go\n");
01463 }
01464
01465 else if( radius.drNext == AV_to_go )
01466 {
01467 fprintf( punch.ipDRout, "A_V to go\n");
01468 }
01469
01470 else if( radius.drNext == drTab )
01471 {
01472 fprintf( punch.ipDRout, "spec den law, new old den%10.2e%10.2e\n",
01473 hdnew, dense.gas_phase[ipHYDROGEN] );
01474 }
01475
01476 else if( radius.drNext == drHMase )
01477 {
01478 fprintf( punch.ipDRout, "H maser dTauMase=%10.2e %li %li %li %li\n",
01479 rt.dTauMase,
01480 rt.mas_species,
01481 rt.mas_ion,
01482 rt.mas_hi,
01483 rt.mas_lo );
01484 }
01485
01486 else if( radius.drNext == drHe1ovHe2 )
01487 {
01488
01489 fprintf( punch.ipDRout, "change in He0/He+ ionization, ratio %.2e\n",
01490 Old_He_atom_ov_ion );
01491 }
01492
01493 else if( radius.drNext == drH2_heat_cool )
01494 {
01495
01496 fprintf( punch.ipDRout, "change in H2 heating/cooling, d(c,h)/H %.2e\n",
01497 dH2_heat_cool );
01498 }
01499
01500 else if( radius.drNext == drH2_abund )
01501 {
01502
01503 fprintf( punch.ipDRout, "change in H2 abundance, d(H2)/H %.2e\n",
01504 dH2_abund );
01505 }
01506
01507 else if( radius.drNext == dr_mole_abund )
01508 {
01509
01510 fprintf( punch.ipDRout, "change in heav ele mole abundance, d(mole)/elem %.2e mole=%i=%s\n",
01511 dCO_abund , mole_dr_change , COmole[mole_dr_change]->label);
01512 }
01513
01514 else if( radius.drNext == drSolomon_BigH2 )
01515 {
01516
01517 fprintf( punch.ipDRout, "change in big H2 Solomon rate line opt depth\n");
01518 }
01519
01520 else if( radius.drNext == recom_dr_last_iter )
01521 {
01522
01523 fprintf( punch.ipDRout, "previous iteration recom logic\n");
01524 }
01525
01526 else if( radius.drNext == drHion )
01527 {
01528 fprintf( punch.ipDRout, "change in H ionization fm to%11.3e%11.3e\n",
01529 SaveOHIIoHI, OHIIoHI );
01530 }
01531
01532 else if( radius.drNext == drEfrac )
01533 {
01534 fprintf( punch.ipDRout,
01535 "change in elec den, rel chng:%11.3e, cur %g old=%g\n",
01536 dEfrac , Efrac_old , Efrac_new );
01537 }
01538
01539 else if( radius.drNext == drdHdStep )
01540 {
01541 fprintf( punch.ipDRout,
01542 "change in heating; current %10.3e delta=%10.3e\n",
01543 thermal.htot, dHdStep );
01544 }
01545
01546 else if( radius.drNext == drLeiden_hack )
01547 {
01548 fprintf( punch.ipDRout,
01549 "Leiden hack\n" );
01550 }
01551
01552 else if( radius.drNext == drConPres )
01553 {
01554 fprintf( punch.ipDRout, "change in con accel\n" );
01555 }
01556
01557 else if( radius.drNext == drdTdStep )
01558 {
01559 fprintf( punch.ipDRout,
01560 "change in temperature; current= %10.3e, dT/T= %.3f\n",
01561 phycon.te, dTdStep );
01562 }
01563
01564 else if( radius.drNext == radius.sdrmin )
01565 {
01566 fprintf( punch.ipDRout, "sdrmin\n" );
01567 }
01568
01569 else if( radius.drNext == radius.sdrmax )
01570 {
01571 fprintf( punch.ipDRout, "sdrmax\n" );
01572 }
01573
01574 else if( radius.drNext == drSphere )
01575 {
01576 fprintf( punch.ipDRout, "sphericity\n" );
01577 }
01578
01579 else if( radius.drNext == dRTaue )
01580 {
01581 fprintf( punch.ipDRout,
01582 "optical depth to electron scattering\n" );
01583 }
01584
01585 else if( radius.drNext == drFail )
01586 {
01587 fprintf( punch.ipDRout,
01588 "temperature failure.\n" );
01589 }
01590
01591 else if( radius.drNext == drmax )
01592 {
01593 fprintf( punch.ipDRout,
01594 "DRMAX; nu opc dr=%10.2e%10.2e%10.2e\n",
01595 freqm, opacm, radius.drChange/
01596 SDIV(opacm) );
01597 }
01598
01599 else if( radius.drNext == drInter )
01600 {
01601 fprintf( punch.ipDRout,
01602 "cont inter nu=%10.2e opac=%10.2e\n",
01603 freqm, opacm );
01604 }
01605
01606 else if( radius.drNext == DrGrainHeat )
01607 {
01608
01609 fprintf( punch.ipDRout,
01610 "grain heating nu=%10.2e opac=%10.2e\n",
01611 grfreqm, gropacm );
01612 }
01613
01614 else if( radius.drNext == winddr )
01615 {
01616 fprintf( punch.ipDRout,
01617 "Wind, dVeldRad=%10.3e\n",
01618 dVeldRad );
01619 }
01620
01621 else if( radius.drNext == drFluc )
01622 {
01623 fprintf( punch.ipDRout,
01624 "density fluctuations\n" );
01625 }
01626
01627 else if( radius.drNext == GlobDr )
01628 {
01629 fprintf( punch.ipDRout,
01630 "GLOB law new dr=%10.2e HDEN=%10.2e\n",
01631 GlobDr,
01632 dense.gas_phase[ipHYDROGEN] );
01633 }
01634
01635 else if( radius.drNext == SpecDr )
01636 {
01637 fprintf( punch.ipDRout,
01638 "special law new dr=%10.2e HDEN=%10.2e\n",
01639 SpecDr,
01640 dense.gas_phase[ipHYDROGEN] );
01641 }
01642
01643 else if( radius.drNext == OldDR )
01644 {
01645 fprintf( punch.ipDRout, "old DR.\n" );
01646 }
01647
01648 else
01649 {
01650 fprintf( punch.ipDRout,
01651 " %4ld radius_next keys from insanity %10.2e\n",
01652 nzone, radius.drNext );
01653
01654 fprintf( ioQQQ,
01655 " %4ld radius_next keys from insanity %10.2e\n",
01656 nzone, radius.drNext );
01657 puts( "[Stop in radius_next]" );
01658 cdEXIT(EXIT_FAILURE);
01659 }
01660
01661
01662 }
01663
01664
01665
01666
01667
01668
01669
01670
01671
01672
01673
01674
01675
01676
01677 if( ((strcmp(dense.chDenseLaw,"DLW1") != 0 &&
01678 strcmp(dense.chDenseLaw ,"GLOB") != 0) )&&
01679
01680
01681
01682 (dense.flong == 0.) &&
01683
01684 radius.drNext != DepthToGo )
01685 {
01686
01687
01688
01689
01690
01691 if( radius.drNext * dense.gas_phase[ipHYDROGEN] < radius.drMinimum )
01692 {
01693 radius.drNext = radius.drMinimum/ dense.gas_phase[ipHYDROGEN];
01694
01695
01696 radius.lgDrMinUsed = true;
01697
01698 lgAbort = true;
01699
01700
01701 --nzone;
01702 fprintf( ioQQQ,
01703 "\n DISASTER PROBLEM radius_next finds dr too small and aborts. This is zone %ld iteration %ld\n\n",
01704 nzone,
01705 iteration);
01706 return(1);
01707 }
01708 }
01709
01710
01711 # define Z 1.0001
01712
01713
01714
01715
01716
01717 radius.drNext = MIN2(radius.drNext,(radius.router[iteration-1]-
01718 radius.depth)*Z);
01719
01720
01721 if( radius.drNext < 0. )
01722 {
01723 radius.lgDrNeg = true;
01724 }
01725 else
01726 {
01727 radius.lgDrNeg = false;
01728 }
01729
01730 ASSERT( radius.drNext > 0. );
01731
01732 if( trace.lgTrace )
01733 {
01734 fprintf( ioQQQ, " radius_next chooses next drad drNext=%.4e; this drad was%12.4e\n",
01735 radius.drNext, radius.drad );
01736 }
01737
01738 DEBUG_EXIT( "radius_next()" );
01739
01740 return( 0 );
01741 }
01742
01743
01744 static void ContRate(double *freqm,
01745 double *opacm)
01746 {
01747 long int i,
01748 ipHi,
01749 iplow,
01750 limit;
01751 double FreqH,
01752 Freq_nonIonizing,
01753 Opac_Hion,
01754 Opac_nonIonizing,
01755 Rate_max_Hion,
01756 Rate_max_nonIonizing;
01757
01758 DEBUG_ENTRY( "ContRate()" );
01759
01760
01761
01762
01763
01764
01765 Rate_max_nonIonizing = 0.;
01766 Freq_nonIonizing = 0.;
01767 Opac_nonIonizing = 0.;
01768
01769
01770 *opacm = -1.;
01771 *freqm = -1.;
01772
01773
01774
01775 if( dense.lgElmtOn[ipCARBON] )
01776 {
01777
01778 ipHi = Heavy.ipHeavy[ipCARBON][0] - 1;
01779 }
01780 else
01781 {
01782
01783 ipHi = iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH2s]-1;
01784 }
01785
01786 for( i=1; i < ipHi; i++ )
01787 {
01788
01789
01790 if( rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*(opac.opacity_abs[i] -
01791 gv.dstab[i]*dense.gas_phase[ipHYDROGEN]) > Rate_max_nonIonizing )
01792 {
01793 Rate_max_nonIonizing = rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*
01794 (opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN]);
01795 Freq_nonIonizing = rfield.anu[i];
01796 Opac_nonIonizing = opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN];
01797 }
01798 }
01799
01800
01801
01802
01803 if( CoolHeavy.brems_heat_total/thermal.htot > 0.05 )
01804 {
01805
01806
01807 iplow = MAX2(1 , rfield.ipEnergyBremsThin );
01808 }
01809 else
01810 {
01811
01812
01813 iplow = ipHi;
01814 }
01815
01816
01817 limit = MIN2(rfield.nflux,iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1);
01818 for( i=iplow; i < limit; i++ )
01819 {
01820 if( rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*(opac.opacity_abs[i] -
01821 gv.dstab[i]*dense.gas_phase[ipHYDROGEN]) > Rate_max_nonIonizing )
01822 {
01823 Rate_max_nonIonizing = rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*
01824 (opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN]);
01825 Freq_nonIonizing = rfield.anu[i];
01826 Opac_nonIonizing = opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN];
01827 }
01828 }
01829
01830
01831 Rate_max_Hion = 0.;
01832 Opac_Hion = 0.;
01833 FreqH = 0.;
01834
01835
01836 for( i=iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1; i < rfield.nflux; i++ )
01837 {
01838 if( rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*(opac.opacity_abs[i] -
01839 gv.dstab[i]*dense.gas_phase[ipHYDROGEN]) > Rate_max_Hion )
01840 {
01841
01842 Rate_max_Hion = rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*
01843 (opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN]);
01844 FreqH = rfield.anu[i];
01845 Opac_Hion = opac.opacity_abs[i] - gv.dstab[i]*dense.gas_phase[ipHYDROGEN];
01846 }
01847 }
01848
01849
01850 if( Rate_max_nonIonizing < 1e-30 && Opac_Hion > SMALLFLOAT )
01851 {
01852
01853 *opacm = Opac_Hion;
01854 *freqm = FreqH;
01855 }
01856
01857
01858 else if( Opac_Hion > Opac_nonIonizing && Rate_max_Hion/Rate_max_nonIonizing > 1e-10 && Opac_Hion > SMALLFLOAT )
01859 {
01860
01861 *opacm = Opac_Hion;
01862 *freqm = FreqH;
01863 }
01864 else
01865 {
01866
01867 *opacm = Opac_nonIonizing;
01868 *freqm = Freq_nonIonizing;
01869 }
01870
01871 # if 0
01872
01873
01874
01875
01876
01877 if( gv.lgDustOn && gv.lgGrainPhysicsOn )
01878 {
01879 double fluxGrainPeak = -1.;
01880 long int ipGrainPeak = -1;
01881 long int ipGrainPeak2 = -1;
01882
01883 i = 0;
01884 while( rfield.anu[i] < 0.0912 && i < (rfield.nflux-2) )
01885 {
01886
01887
01888 if( gv.GrainEmission[i]*rfield.anu[i]*opac.opacity_abs[i] > fluxGrainPeak )
01889 {
01890 ipGrainPeak = i;
01891 fluxGrainPeak = gv.GrainEmission[i]*rfield.anu[i]*opac.opacity_abs[i];
01892 }
01893 ++i;
01894 }
01895 ASSERT( fluxGrainPeak>=0. && ipGrainPeak >=0 );
01896
01897
01898
01899
01900
01901
01902
01903 i = ipGrainPeak;
01904
01905
01906
01907
01908 if( fluxGrainPeak > 0. )
01909 {
01910 while( rfield.anu[i] < 0.0912 && i < (rfield.nflux-2) )
01911 {
01912
01913 if( gv.GrainEmission[i]*rfield.anu[i]*opac.opacity_abs[i] > 0.05*fluxGrainPeak )
01914 {
01915
01916 ipGrainPeak2 = i;
01917 }
01918 ++i;
01919 }
01920 ASSERT( ipGrainPeak2>=0 );
01921
01922
01923 if( opac.opacity_abs[ipGrainPeak2] > *opacm )
01924 {
01925
01926 *opacm = opac.opacity_abs[ipGrainPeak2]*10.;
01927 *freqm = rfield.anu[ipGrainPeak2];
01928
01929
01930 }
01931 }
01932 }
01933 # endif
01934
01935 {
01936
01937
01938 enum {DEBUG_LOC=false};
01939
01940 if( DEBUG_LOC )
01941 {
01942 fprintf(ioQQQ,"conratedebug \t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
01943 Rate_max_nonIonizing,Freq_nonIonizing,Opac_nonIonizing,
01944 Rate_max_Hion,FreqH ,Opac_Hion,*freqm,*opacm
01945 );
01946
01947 }
01948 }
01949
01950
01951
01952
01953
01954
01955 ASSERT( *opacm >= 0. && *freqm >= 0. );
01956
01957 DEBUG_EXIT( "ContRate()" );
01958 return;
01959 }
01960
01961
01962 static void GrainRateDr(double *grfreqm,
01963 double *gropacm)
01964 {
01965 long int i,
01966 iplow;
01967 double xMax;
01968
01969 DEBUG_ENTRY( "GrainRateDr()" );
01970
01971
01972
01973
01974
01975
01976
01977
01978 if( CoolHeavy.brems_heat_total/thermal.htot > 0.05 )
01979 {
01980
01981
01982 iplow = MAX2(1 , rfield.ipEnergyBremsThin );
01983 }
01984 else
01985 {
01986
01987
01988 if( dense.lgElmtOn[ipCARBON] )
01989 {
01990
01991 iplow = Heavy.ipHeavy[ipCARBON][0];
01992 }
01993 else
01994 {
01995
01996 iplow = iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH2s];
01997 }
01998
01999 }
02000
02001 xMax = -1.;
02002
02003 for( i=iplow-1; i < Heavy.ipHeavy[ipHYDROGEN][0]; i++ )
02004 {
02005 if( rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*opac.opacity_abs[i] > xMax )
02006 {
02007 xMax = rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*
02008 opac.opacity_abs[i];
02009 *grfreqm = rfield.anu[i];
02010 *gropacm = opac.opacity_abs[i];
02011 }
02012 }
02013
02014
02015 if( dense.lgElmtOn[ipHELIUM] )
02016 {
02017 for( i=Heavy.ipHeavy[ipHYDROGEN][0]-1; i < Heavy.ipHeavy[ipHELIUM][1]; i++ )
02018 {
02019 if( rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*opac.opacity_abs[i] > xMax )
02020 {
02021 xMax = rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*
02022 opac.opacity_abs[i];
02023 *grfreqm = rfield.anu[i];
02024 *gropacm = opac.opacity_abs[i];
02025 }
02026 }
02027 }
02028
02029
02030
02031 if( xMax <= 0. )
02032 {
02033 *gropacm = 0.;
02034 *grfreqm = 0.;
02035 }
02036
02037 DEBUG_EXIT( "GrainRateDr()" );
02038 return;
02039 }
02040 #if 0
02041
02042 static void ChkRate(
02043
02044 long int nelem,
02045
02046 double *dDestRate,
02047
02048 double *DestRateOld,
02049 double *DestRateNew,
02050
02051 long int *istage)
02052 {
02053 long int i;
02054
02055 double average,
02056 dDest;
02057
02058 static double OldDest[LIMELM][LIMELM];
02059
02060 DEBUG_ENTRY( "ChkRate()" );
02061
02062
02063 if( !dense.lgElmtOn[nelem] )
02064 {
02065 *dDestRate = 1e-3;
02066 *DestRateOld = 1e-3;
02067 *DestRateNew = 1e-3;
02068 *istage = 0;
02069
02070 DEBUG_EXIT( "ChkRate()" );
02071 return;
02072 }
02073
02074
02075
02076 *istage = 1;
02077 *dDestRate = 0.;
02078 *DestRateOld = 0.;
02079 *DestRateNew = 0.;
02080 *dDestRate = 0.;
02081
02082 if( nzone <= 1 )
02083 {
02084 for( i=0; i < nelem+1; i++ )
02085 {
02086 OldDest[nelem][i] = ionbal.RateIonizTot[nelem][i];
02087 }
02088 }
02089
02090 else if( dense.xIonDense[nelem][0]/dense.gas_phase[nelem] < 0.9 )
02091 {
02092
02093 for( i=0; i < (nelem); i++ )
02094 {
02095
02096
02097
02098 if( ((dense.xIonDense[nelem][i]/dense.gas_phase[nelem] > 0.01 &&
02099 dense.xIonDense[nelem][i]/dense.gas_phase[nelem] < 0.9) &&
02100 dense.xIonDense[nelem][i+1]/dense.gas_phase[nelem] > .05) &&
02101 OldDest[nelem][i] > 0. )
02102 {
02103
02104
02105
02106 if( ionbal.RateIonizTot[nelem][i] <= 0. )
02107 {
02108 fprintf( ioQQQ, " ChkRate gets insane destruction rate for ion%4ld%4ld%10.2e\n",
02109 nelem+1, i, ionbal.RateIonizTot[nelem][i] );
02110 puts( "[Stop in chkrate]" );
02111 cdEXIT(EXIT_FAILURE);
02112 }
02113
02114
02115
02116
02117
02118 average = (OldDest[nelem][i] + ionbal.RateIonizTot[nelem][i])* 0.5;
02119
02120 dDest = (OldDest[nelem][i] - ionbal.RateIonizTot[nelem][i])/ average;
02121
02122
02123 if( *dDestRate < dDest )
02124 {
02125
02126 *dDestRate = dDest;
02127 *istage = i+1;
02128 *DestRateOld = OldDest[nelem][i];
02129 *DestRateNew = ionbal.RateIonizTot[nelem][i];
02130 }
02131
02132 }
02133 OldDest[nelem][i] = ionbal.RateIonizTot[nelem][i];
02134 }
02135 }
02136
02137 DEBUG_EXIT( "ChkRate()" );
02138 return;
02139 }
02140 #endif
02141