00001
00002
00003
00004 #include "cddefines.h"
00005 #include "physconst.h"
00006 #include "iso.h"
00007 #include "hydrogenic.h"
00008 #include "rfield.h"
00009 #include "colden.h"
00010 #include "geometry.h"
00011 #include "opacity.h"
00012 #include "thermal.h"
00013 #include "doppvel.h"
00014 #include "dense.h"
00015 #include "mole.h"
00016 #include "h2.h"
00017 #include "timesc.h"
00018 #include "hmi.h"
00019 #include "lines_service.h"
00020 #include "taulines.h"
00021 #include "trace.h"
00022 #include "wind.h"
00023 #include "phycon.h"
00024 #include "pressure.h"
00025 #include "grainvar.h"
00026 #include "molcol.h"
00027 #include "conv.h"
00028 #include "hyperfine.h"
00029 #include "mean.h"
00030 #include "struc.h"
00031 #include "radius.h"
00032 #include "gravity.h"
00033
00034 void radius_increment(void)
00035 {
00036
00037 long int nelem,
00038 i,
00039 ion,
00040 nzone_minus_1 ,
00041 mol;
00042 double
00043 avWeight,
00044 t;
00045
00046 double ajmass,
00047 Error,
00048 rjeans;
00049
00050 DEBUG_ENTRY( "radius_increment()" );
00051
00052
00053
00054
00055
00056
00057 nzone_minus_1 = MAX2( 0, nzone-1 );
00058 ASSERT(nzone_minus_1>=0 && nzone_minus_1 < struc.nzlim );
00059
00060 struc.heatstr[nzone_minus_1] = thermal.htot;
00061 struc.coolstr[nzone_minus_1] = thermal.ctot;
00062 struc.testr[nzone_minus_1] = (realnum)phycon.te;
00063
00064
00065 struc.DenParticles[nzone_minus_1] = dense.pden;
00066
00067 struc.hden[nzone_minus_1] = (realnum)dense.gas_phase[ipHYDROGEN];
00068
00069 struc.DenMass[nzone_minus_1] = dense.xMassDensity;
00070 struc.volstr[nzone_minus_1] = (realnum)radius.dVeffAper;
00071 struc.drad[nzone_minus_1] = (realnum)radius.drad;
00072 struc.drad_x_fillfac[nzone_minus_1] = (realnum)radius.drad_x_fillfac;
00073 struc.histr[nzone_minus_1] = dense.xIonDense[ipHYDROGEN][0];
00074 struc.hiistr[nzone_minus_1] = dense.xIonDense[ipHYDROGEN][1];
00075 struc.ednstr[nzone_minus_1] = (realnum)dense.eden;
00076 struc.o3str[nzone_minus_1] = dense.xIonDense[ipOXYGEN][2];
00077 struc.pressure[nzone_minus_1] = (realnum)pressure.PresTotlCurr;
00078 struc.windv[nzone_minus_1] = (realnum)wind.windv;
00079 struc.AccelTotalOutward[nzone_minus_1] = wind.AccelTotalOutward;
00080 struc.AccelGravity[nzone_minus_1] = wind.AccelGravity;
00081 struc.pres_radiation_lines_curr[nzone_minus_1] = (realnum)pressure.pres_radiation_lines_curr;
00082 struc.GasPressure[nzone_minus_1] = (realnum)pressure.PresGasCurr;
00083 struc.depth[nzone_minus_1] = (realnum)radius.depth;
00084
00085 struc.xLyman_depth[nzone_minus_1] = opac.TauAbsFace[iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon];
00086 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
00087 {
00088 struc.gas_phase[nelem][nzone_minus_1] = dense.gas_phase[nelem];
00089 for( ion=0; ion<nelem+2; ++ion )
00090 {
00091 struc.xIonDense[nelem][ion][nzone_minus_1] = dense.xIonDense[nelem][ion];
00092 }
00093 }
00094 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00095 {
00096 for( nelem=ipISO; nelem<LIMELM; ++nelem)
00097 {
00098 if( dense.lgElmtOn[nelem] )
00099 {
00100 for( long level=0; level < iso_sp[ipISO][nelem].numLevels_max; ++level )
00101 {
00102 struc.StatesElem[nelem][nelem-ipISO][level][nzone_minus_1] = (realnum)iso_sp[ipISO][nelem].st[level].Pop();
00103 }
00104 }
00105 }
00106 }
00107
00108
00109 for(mol=0;mol<mole_global.num_calc;mol++)
00110 {
00111 struc.molecules[mol][nzone_minus_1] = (realnum) mole.species[mol].den;
00112 }
00113 struc.H2_abund[nzone_minus_1] = hmi.H2_total;
00114
00115
00116
00117 if( lgAbort )
00118 {
00119 return;
00120 }
00121
00122 if( trace.lgTrace )
00123 {
00124 fprintf( ioQQQ,
00125 " radius_increment called; radius=%10.3e rinner=%10.3e DRAD=%10.3e drNext=%10.3e ROUTER=%10.3e DEPTH=%10.3e\n",
00126 radius.Radius, radius.rinner, radius.drad, radius.drNext,
00127 radius.StopThickness[iteration-1], radius.depth );
00128 }
00129
00130
00131 Error = fabs(dense.eden - dense.EdenTrue)/SDIV(dense.EdenTrue);
00132 if( Error > conv.BigEdenError )
00133 {
00134 conv.BigEdenError = (realnum)Error;
00135 dense.nzEdenBad = nzone;
00136 }
00137 conv.AverEdenError += (realnum)Error;
00138
00139 dense.EdenMax = MAX2( dense.eden , dense.EdenMax);
00140 dense.EdenMin = MIN2( dense.eden , dense.EdenMin);
00141
00142
00143 Error = fabs(thermal.ctot - thermal.htot) / thermal.ctot;
00144 conv.BigHeatCoolError = MAX2((realnum)Error , conv.BigHeatCoolError );
00145 conv.AverHeatCoolError += (realnum)Error;
00146
00147
00148 Error = fabs(pressure.PresTotlError);
00149 conv.BigPressError = MAX2((realnum)Error , conv.BigPressError );
00150 conv.AverPressError += (realnum)Error;
00151
00152
00153 dense.xMassTotal += (realnum)(dense.xMassDensity * radius.dVeffVol);
00154
00155
00156 timesc.time_therm_long = MAX2(timesc.time_therm_long,1.5*dense.pden*BOLTZMANN*phycon.te/
00157 thermal.ctot);
00158
00159
00160
00161 t = H21cm_H_atom( phycon.te )* dense.xIonDense[ipHYDROGEN][0] +
00162
00163
00164 H21cm_electron( phycon.te )*dense.eden;
00165
00166
00167 if( t > SMALLFLOAT )
00168 timesc.TimeH21cm = MAX2( 1./t, timesc.TimeH21cm );
00169
00170
00171 if( (double)dense.xIonDense[ipCARBON][0]*(double)dense.xIonDense[ipOXYGEN][0] > SMALLFLOAT )
00172 {
00173 int ipCO = findspecies("CO")->index;
00174
00175 if (ipCO != -1)
00176 timesc.BigCOMoleForm = MAX2( timesc.BigCOMoleForm, 1./SDIV(mole.species[ipCO].snk ));
00177 }
00178
00179
00180 timesc.time_H2_Dest_longest = MAX2(timesc.time_H2_Dest_longest, timesc.time_H2_Dest_here );
00181
00182
00183 timesc.time_H2_Form_longest = MAX2( timesc.time_H2_Form_longest , timesc.time_H2_Form_here );
00184
00185
00186
00187
00188 if( thermal.lgUnstable )
00189 thermal.nUnstable += 1;
00190
00191
00192 if( !rfield.lgUSphON && dense.xIonDense[ipHYDROGEN][0]/dense.gas_phase[ipHYDROGEN] > 0.49 )
00193 {
00194 rfield.rstrom = (realnum)radius.Radius;
00195 rfield.lgUSphON = true;
00196 }
00197
00198
00199 wind.AccelMax = (realnum)MAX2(wind.AccelMax,wind.AccelTotalOutward);
00200
00201
00202 wind.AccelAver += wind.AccelTotalOutward*(realnum)radius.drad_x_fillfac;
00203 wind.acldr += (realnum)radius.drad_x_fillfac;
00204
00205
00206 pressure.pinzon = (realnum)(wind.AccelTotalOutward*dense.xMassDensity*
00207 radius.drad_x_fillfac*geometry.DirectionalCosin);
00208
00209
00210 pressure.PresInteg += pressure.pinzon;
00211
00212
00213
00214 static realnum AccelElecScatZone1;
00215 if( nzone == 1 )
00216 AccelElecScatZone1 = wind.AccelElectron;
00217 pressure.pinzon_PresIntegElecThin = (realnum)(AccelElecScatZone1*dense.xMassDensity/radius.r1r0sq*
00218 radius.drad_x_fillfac*geometry.DirectionalCosin);
00219 pressure.PresIntegElecThin += pressure.pinzon_PresIntegElecThin;
00220
00221
00222 GravitationalPressure();
00223 pressure.IntegRhoGravity += pressure.RhoGravity;
00224
00225
00226 timesc.sound_speed_isothermal = sqrt(pressure.PresGasCurr/dense.xMassDensity);
00227
00228 timesc.sound_speed_adiabatic = sqrt(5.*pressure.PresGasCurr/(3.*dense.xMassDensity) );
00229 timesc.sound += radius.drad_x_fillfac / timesc.sound_speed_isothermal;
00230
00231
00232
00233
00234
00235
00236 if( iteration > 1 && nzone_minus_1 < struc.nzonePreviousIteration )
00237 {
00238
00239
00240 realnum TempChange = (realnum)
00241 fabs( (phycon.te-struc.testr[nzone_minus_1])/phycon.te);
00242
00243 struc.TempChangeMax = MAX2( struc.TempChangeMax , TempChange );
00244 }
00245 else
00246 {
00247
00248 struc.TempChangeMax = 0.;
00249 }
00250
00251
00252
00253 colden.dlnenp += dense.eden*(double)(dense.xIonDense[ipHYDROGEN][1])*radius.drad_x_fillfac;
00254 colden.dlnenHep += dense.eden*(double)(dense.xIonDense[ipHELIUM][1])*radius.drad_x_fillfac;
00255 colden.dlnenHepp += dense.eden*(double)(dense.xIonDense[ipHELIUM][2])*radius.drad_x_fillfac;
00256 colden.dlnenCp += dense.eden*(double)(dense.xIonDense[ipCARBON][1])*radius.drad_x_fillfac;
00257
00258
00259 for( unsigned i = 0; i < mole.species.size(); ++i )
00260 {
00261 if( mole.species[i].levels != NULL )
00262 {
00263 for( qList::iterator st = mole.species[i].levels->begin(); st != mole.species[i].levels->end(); ++st )
00264 {
00265 (*st).ColDen() += radius.drad_x_fillfac * (*st).Pop();
00266 }
00267 }
00268 }
00269
00270
00271 if( hyperfine.Tspin21cm > SMALLFLOAT )
00272 colden.H0_ov_Tspin += (double)(dense.xIonDense[ipHYDROGEN][0]) / hyperfine.Tspin21cm*radius.drad_x_fillfac;
00273
00274
00275 colden.OH_ov_Tspin += (double)(findspecieslocal("OH")->den) / phycon.te*radius.drad_x_fillfac;
00276
00277
00278 hydro.TexcLya = (realnum)TexcLine( iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s) );
00279
00280
00281 if( hydro.TexcLya > phycon.te )
00282 {
00283 hydro.nLyaHot += 1;
00284 if( hydro.TexcLya > hydro.TLyaMax )
00285 {
00286 hydro.TLyaMax = hydro.TexcLya;
00287 hydro.TeLyaMax = (realnum)phycon.te;
00288 hydro.nZTLaMax = nzone;
00289 }
00290 }
00291
00292
00293 colden.colden[ipCOL_HTOT] += (realnum)(dense.gas_phase[ipHYDROGEN]*radius.drad_x_fillfac);
00294 ASSERT( colden.colden[ipCOL_HTOT] >SMALLFLOAT );
00295 colden.colden[ipCOL_HMIN] += (realnum)(findspecieslocal("H-")->den*radius.drad_x_fillfac);
00296
00297
00298 colden.colden[ipCOL_H2g] += (realnum)(findspecieslocal("H2")->den*(realnum)radius.drad_x_fillfac);
00299 colden.colden[ipCOL_H2s] += (realnum)(findspecieslocal("H2*")->den*(realnum)radius.drad_x_fillfac);
00300
00301 colden.coldenH2_ov_vel += hmi.H2_total*(realnum)radius.drad_x_fillfac / GetDopplerWidth(2.f*dense.AtomicWeight[ipHYDROGEN]);
00302 colden.colden[ipCOL_HeHp] += (realnum)(findspecieslocal("HeH+")->den*(realnum)radius.drad_x_fillfac);
00303 colden.colden[ipCOL_H2p] += (realnum)(findspecieslocal("H2+")->den*(realnum)radius.drad_x_fillfac);
00304 colden.colden[ipCOL_H3p] += (realnum)(findspecieslocal("H3+")->den*(realnum)radius.drad_x_fillfac);
00305 colden.colden[ipCOL_Hp] += dense.xIonDense[ipHYDROGEN][1]*(realnum)radius.drad_x_fillfac;
00306 colden.colden[ipCOL_H0] += dense.xIonDense[ipHYDROGEN][0]*(realnum)radius.drad_x_fillfac;
00307
00308 colden.colden[ipCOL_elec] += (realnum)(dense.eden*radius.drad_x_fillfac);
00309
00310 colden.He123S += iso_sp[ipHE_LIKE][ipHELIUM].st[ipHe2s3S].Pop()*(realnum)radius.drad_x_fillfac;
00311
00312
00313 h2.ortho_colden += h2.ortho_density*radius.drad_x_fillfac;
00314 h2.para_colden += h2.para_density*radius.drad_x_fillfac;
00315 if( hmi.H2_total > SMALLFLOAT )
00316 ASSERT( fabs( h2.ortho_density + h2.para_density - hmi.H2_total ) / hmi.H2_total < 1e-4 );
00317
00318
00319
00320
00321
00322 colden.H0_21cm_upper += ((*HFLines[0].Hi()).Pop()*radius.drad_x_fillfac);
00323 colden.H0_21cm_lower += ((*HFLines[0].Lo()).Pop()*radius.drad_x_fillfac);
00324
00325
00326
00327
00328
00329 for( i=0; i < 5; i++ )
00330 {
00331
00332 colden.C2Colden[i] += colden.C2Pops[i]*(realnum)radius.drad_x_fillfac;
00333
00334 colden.Si2Colden[i] += colden.Si2Pops[i]*(realnum)radius.drad_x_fillfac;
00335 }
00336 for( i=0; i < 3; i++ )
00337 {
00338
00339 colden.C1Colden[i] += colden.C1Pops[i]*(realnum)radius.drad_x_fillfac;
00340
00341 colden.O1Colden[i] += colden.O1Pops[i]*(realnum)radius.drad_x_fillfac;
00342 }
00343 for( i=0; i < 4; i++ )
00344 {
00345
00346 colden.C3Colden[i] += colden.C3Pops[i]*(realnum)radius.drad_x_fillfac;
00347 }
00348
00349
00350 molcol("ADD ",ioQQQ);
00351
00352
00353 mean.MeanInc();
00354
00355
00356
00357
00358 avWeight = 0.;
00359 for( nelem=0; nelem < LIMELM; nelem++ )
00360 {
00361 avWeight += dense.gas_phase[nelem]*dense.AtomicWeight[nelem];
00362 }
00363 avWeight /= dense.gas_phase[ipHYDROGEN];
00364
00365
00366 rfield.opac_mag_B_point = 0.;
00367 rfield.opac_mag_V_point = 0.;
00368 rfield.opac_mag_B_extended = 0.;
00369 rfield.opac_mag_V_extended = 0.;
00370 for( size_t nd=0; nd < gv.bin.size(); nd++ )
00371 {
00372
00373
00374
00375
00376 rfield.opac_mag_B_point += (gv.bin[nd]->dstab1[rfield.ipB_filter-1] +
00377 gv.bin[nd]->pure_sc1[rfield.ipB_filter-1])*double(gv.bin[nd]->dstAbund)*
00378 double(dense.gas_phase[ipHYDROGEN]) * OPTDEP2EXTIN;
00379
00380 rfield.opac_mag_V_point += (gv.bin[nd]->dstab1[rfield.ipV_filter-1] +
00381 gv.bin[nd]->pure_sc1[rfield.ipV_filter-1])*double(gv.bin[nd]->dstAbund)*
00382 double(dense.gas_phase[ipHYDROGEN]) * OPTDEP2EXTIN;
00383
00384
00385
00386
00387 rfield.opac_mag_B_extended += (gv.bin[nd]->dstab1[rfield.ipB_filter-1] +
00388 gv.bin[nd]->pure_sc1[rfield.ipB_filter-1]*gv.bin[nd]->asym[rfield.ipB_filter-1])*
00389 double(gv.bin[nd]->dstAbund)*double(dense.gas_phase[ipHYDROGEN]) * OPTDEP2EXTIN;
00390
00391 rfield.opac_mag_V_extended += (gv.bin[nd]->dstab1[rfield.ipV_filter-1] +
00392 gv.bin[nd]->pure_sc1[rfield.ipV_filter-1]*gv.bin[nd]->asym[rfield.ipV_filter-1])*
00393 double(gv.bin[nd]->dstAbund)*double(dense.gas_phase[ipHYDROGEN]) * OPTDEP2EXTIN;
00394
00395 gv.bin[nd]->avdust += gv.bin[nd]->tedust*(realnum)radius.drad_x_fillfac;
00396 gv.bin[nd]->avdft += gv.bin[nd]->DustDftVel*(realnum)radius.drad_x_fillfac;
00397 gv.bin[nd]->avdpot += (realnum)(gv.bin[nd]->dstpot*EVRYD*radius.drad_x_fillfac);
00398 gv.bin[nd]->avDGRatio += (realnum)(gv.bin[nd]->dustp[1]*gv.bin[nd]->dustp[2]*
00399 gv.bin[nd]->dustp[3]*gv.bin[nd]->dustp[4]*gv.bin[nd]->dstAbund/avWeight*
00400 radius.drad_x_fillfac);
00401 }
00402
00403
00404
00405
00406
00407 rfield.extin_mag_B_point += rfield.opac_mag_B_point*radius.drad_x_fillfac;
00408 rfield.extin_mag_V_point += rfield.opac_mag_V_point*radius.drad_x_fillfac;
00409 rfield.extin_mag_B_extended += rfield.opac_mag_B_extended*radius.drad_x_fillfac;
00410 rfield.extin_mag_V_extended += rfield.opac_mag_V_extended*radius.drad_x_fillfac;
00411
00412
00413 colden.TotMassColl += dense.xMassDensity*(realnum)radius.drad_x_fillfac;
00414 colden.tmas += (realnum)phycon.te*dense.xMassDensity*(realnum)radius.drad_x_fillfac;
00415 colden.wmas += dense.wmole*dense.xMassDensity*(realnum)radius.drad_x_fillfac;
00416
00417
00418 rjeans = 7.79637 + (phycon.alogte - log10(dense.wmole) - log10(dense.xMassDensity*
00419 geometry.FillFac))/2.;
00420
00421
00422 ajmass = 3.*(rjeans - 0.30103) + log10(4.*PI/3.*dense.xMassDensity*
00423 geometry.FillFac);
00424
00425
00426 colden.rjnmin = MIN2(colden.rjnmin,(realnum)rjeans);
00427 colden.ajmmin = MIN2(colden.ajmmin,(realnum)ajmass);
00428
00429 if( trace.lgTrace )
00430 {
00431 fprintf( ioQQQ, " radius_increment returns\n" );
00432 }
00433 return;
00434 }