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