cloudy trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /*radius_increment do work associated with geometry increments of this zone, called before RT_tau_inc */ 00004 /*pnegopc punch negative opacities on io unit, iff 'set negopc' command was given */ 00005 #include "cddefines.h" 00006 #include "physconst.h" 00007 #include "iso.h" 00008 #include "hydrogenic.h" 00009 #include "dynamics.h" 00010 #include "rfield.h" 00011 #include "hextra.h" 00012 #include "colden.h" 00013 #include "geometry.h" 00014 #include "opacity.h" 00015 #include "thermal.h" 00016 #include "doppvel.h" 00017 #include "dense.h" 00018 #include "mole.h" 00019 #include "h2.h" 00020 #include "timesc.h" 00021 #include "hmi.h" 00022 #include "lines_service.h" 00023 #include "taulines.h" 00024 #include "trace.h" 00025 #include "wind.h" 00026 #include "phycon.h" 00027 #include "pressure.h" 00028 #include "grainvar.h" 00029 #include "molcol.h" 00030 #include "conv.h" 00031 #include "hyperfine.h" 00032 #include "mean.h" 00033 #include "struc.h" 00034 #include "radius.h" 00035 /*cmshft - so Compton scattering shift of spectrum */ 00036 STATIC void cmshft(void); 00037 00038 #if !defined(NDEBUG) 00039 /*pnegopc punch negative opacities on io unit, iff 'set negopc' command was 00040 * given. This function is only used in debug mode */ 00041 STATIC void pnegopc(void); 00042 #endif 00043 00044 void radius_increment(void) 00045 { 00046 # if !defined(NDEBUG) 00047 bool lgFlxNeg; 00048 # endif 00049 00050 long int nelem, 00051 i, 00052 j , 00053 ion, 00054 nzone_minus_1 , 00055 mol; 00056 double 00057 avWeight, 00058 DilutionHere , 00059 escatc, 00060 fmol, 00061 opfac, 00062 relec, 00063 rforce, 00064 t; 00065 00066 double ajmass, 00067 Error, 00068 rjeans; 00069 00070 realnum dradfac; 00071 00072 DEBUG_ENTRY( "radius_increment()" ); 00073 00074 /* when this sub is called radius is the outer edge of zone */ 00075 00076 if( lgAbort ) 00077 { 00078 return; 00079 } 00080 00081 /* following block only set of asserts to check that iso populations are sane */ 00082 # if !defined(NDEBUG) 00083 if( !dynamics.lgAdvection ) 00084 { 00085 long int ipISO; 00086 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO ) 00087 { 00088 for( nelem=ipISO; nelem<LIMELM; ++nelem ) 00089 { 00090 /* >>chng 05 feb 05, rm div by gas_phase for element, to eden, 00091 * see explain for this date below */ 00092 if( dense.lgElmtOn[nelem] && 00093 dense.xIonDense[nelem][nelem-ipISO]/dense.eden>conv.EdenErrorAllowed / 10. && 00094 !(iso.chTypeAtomUsed[ipISO][nelem][0]=='L' && 00095 iso.xIonSimple[ipISO][nelem]<1e-10) ) 00096 { 00097 /* we will check that ground pops are within this factor of the xIonDense value 00098 * these are only converged down to some extent - probably this should be related 00099 * to the ionization convergence error */ 00100 /* >>chng 05 feb 05, for very highly ionized sims, where N or O is He-like, 00101 * the error can approach 1%, and depends on the convergence criteria. 00102 * the code does not try to converge these quantities in an absolute sense, 00103 * only the quantities relative to the electron or hydrogen density. so it is not 00104 * safe to do the assert for small values 00105 * change is to only do this branch if abundance is large enough to be detected by 00106 * the ionization convergence solvers */ 00107 # define ERR_CHK 1.001 00108 double err_chk = ERR_CHK; 00109 /* >>chng 05 sep 02, when low-T solver used solution is approximate, 00110 * and it must not matter (lot-T solver should not be used if it 00111 * does matter - so use more liberal expectation */ 00112 if( iso.chTypeAtomUsed[ipISO][nelem][0]=='L' ) 00113 err_chk *= 10.; 00114 /* make sure that populations are valid, first check Pop2Ion */ 00115 if( StatesElem[ipISO][nelem][0].Pop > 00116 dense.xIonDense[nelem][nelem-ipISO]/(double)SDIV(dense.xIonDense[nelem][nelem+1-ipISO])*err_chk && 00117 /*>>chng 06 jul 22, only do check if pops are within tolerance of double */ 00118 iso.pop_ion_ov_neut[ipISO][nelem] > 1e-15 ) 00119 { 00120 fprintf(ioQQQ, 00121 " DISASTER radius_increment finds inconsistent populations, Pop2Ion zn %.2f ipISO %li nelem %li %.4e %.4e solver:%s\n", 00122 fnzone, 00123 ipISO , nelem , 00124 StatesElem[ipISO][nelem][0].Pop*dense.xIonDense[nelem][nelem+1-ipISO] , 00125 dense.xIonDense[nelem][nelem-ipISO] , 00126 iso.chTypeAtomUsed[ipISO][nelem]); 00127 fprintf(ioQQQ," level solver %s found pop_ion_ov_neut of %.5e", 00128 iso.chTypeAtomUsed[ipISO][nelem] , 00129 iso.pop_ion_ov_neut[ipISO][nelem] ); 00130 fprintf(ioQQQ," ion_solver found pop_ion_ov_neut of %.5e\n", 00131 dense.xIonDense[nelem][nelem+1-ipISO]/SDIV( dense.xIonDense[nelem][nelem-ipISO] ) ); 00132 fprintf(ioQQQ, 00133 "simple %.3e Test shows Pop2Ion (%.6e) > xIonDense[nelem]/[nelem+1] (%.6e) \n", 00134 iso.xIonSimple[ipISO][nelem], 00135 StatesElem[ipISO][nelem][0].Pop , 00136 dense.xIonDense[nelem][nelem-ipISO]/(double)SDIV(dense.xIonDense[nelem][nelem+1-ipISO]) ); 00137 } 00138 00139 /* >>chng 06 jul 24, add assert that this is not search 00140 * phase to confirm that removing other asserts on this 00141 * are ok - we cannot be in search phase in this routine 00142 * if ok, rm this and all ref to lgSearch, perhaps conv.h header */ 00143 ASSERT( !conv.lgSearch ); 00144 ASSERT( StatesElem[ipISO][nelem][0].Pop <= 00145 dense.xIonDense[nelem][nelem-ipISO]/ 00146 (double)SDIV(dense.xIonDense[nelem][nelem+1-ipISO])*err_chk || 00147 /*>>chng 06 jul 22, only do check if pops are within tolerance of double */ 00148 iso.pop_ion_ov_neut[ipISO][nelem] < 1e-15); 00149 00150 /* make sure that populations are valid, first check Pop2Ion */ 00151 if( StatesElem[ipISO][nelem][0].Pop* 00152 dense.xIonDense[nelem][nelem+1-ipISO] > 00153 dense.xIonDense[nelem][nelem-ipISO]*err_chk && 00154 /*>>chng 06 jul 22, only do check if pops are within tolerance of double */ 00155 iso.pop_ion_ov_neut[ipISO][nelem] > 1e-15) 00156 { 00157 fprintf(ioQQQ," DISASTER PopLo fnz %.2f ipISO %li nelem %li pop_ion_ov_neut %.2e 2pops %.6e %.6e solver:%s\n", 00158 fnzone, 00159 ipISO , nelem , 00160 iso.pop_ion_ov_neut[ipISO][nelem] , 00161 StatesElem[ipISO][nelem][0].Pop*dense.xIonDense[nelem][nelem+1-ipISO], 00162 dense.xIonDense[nelem][nelem-ipISO]*err_chk , 00163 iso.chTypeAtomUsed[ipISO][nelem]); 00164 } 00165 ASSERT( 00166 /*(conv.lgSearch || !conv.nPres2Ioniz) || */ 00167 (StatesElem[ipISO][nelem][0].Pop*dense.xIonDense[nelem][nelem+1-ipISO]<= 00168 dense.xIonDense[nelem][nelem-ipISO]*err_chk) || 00169 /*>>chng 06 jul 22, only do check if pops are within tolerance of double */ 00170 iso.pop_ion_ov_neut[ipISO][nelem] < 1e-15); 00171 # undef ERR_CHK 00172 } 00173 } 00174 } 00175 } 00176 # endif 00177 00178 if( trace.lgTrace ) 00179 { 00180 fprintf( ioQQQ, 00181 " radius_increment called; radius=%10.3e rinner=%10.3e DRAD=%10.3e drNext=%10.3e ROUTER=%10.3e DEPTH=%10.3e\n", 00182 radius.Radius, radius.rinner, radius.drad, radius.drNext, 00183 radius.router[iteration-1], radius.depth ); 00184 } 00185 00186 /* remember mean and largest errors on electron density */ 00187 Error = fabs(dense.eden - dense.EdenTrue)/SDIV(dense.EdenTrue); 00188 if( Error > conv.BigEdenError ) 00189 { 00190 conv.BigEdenError = (realnum)Error; 00191 dense.nzEdenBad = nzone; 00192 } 00193 conv.AverEdenError += (realnum)Error; 00194 00195 /* remember mean and largest errors between heating and cooling */ 00196 Error = fabs(thermal.ctot - thermal.htot) / thermal.ctot; 00197 conv.BigHeatCoolError = MAX2((realnum)Error , conv.BigHeatCoolError ); 00198 conv.AverHeatCoolError += (realnum)Error; 00199 00200 /* remember mean and largest pressure errors */ 00201 Error = fabs(pressure.PresTotlCurr - pressure.PresTotlCorrect) / pressure.PresTotlCorrect; 00202 conv.BigPressError = MAX2((realnum)Error , conv.BigPressError ); 00203 conv.AverPressError += (realnum)Error; 00204 00205 /* integrate total mass over model */ 00206 dense.xMassTotal += dense.xMassDensity * (realnum)(radius.drad_x_fillfac*radius.r1r0sq); 00207 00208 /* check cooling time for this zone, remember longest */ 00209 timesc.time_therm_long = MAX2(timesc.time_therm_long,1.5*dense.pden*BOLTZMANN*phycon.te/ 00210 thermal.ctot); 00211 00212 /* fmol is fraction of hydrogen in form of any molecule or ion */ 00213 /* if N_H_MOLEC != 8, we probably need to change this line... */ 00214 ASSERT( N_H_MOLEC == 8); 00215 fmol = (hmi.Hmolec[ipMHm] + 2.*(hmi.H2_total + hmi.Hmolec[ipMH2p]))/dense.gas_phase[ipHYDROGEN]; 00216 00217 /* remember the largest H2 fraction that occurs in the model */ 00218 hmi.BiggestH2 = MAX2(hmi.BiggestH2,(realnum)fmol); 00219 00220 /* largest fraction of atoms in molecules */ 00221 for( i=0; i<mole.num_comole_calc; ++i ) 00222 { 00223 if( dense.lgElmtOn[COmole[i]->nelem_hevmol] ) 00224 { 00225 realnum frac = COmole[i]->hevmol / SDIV(dense.gas_phase[COmole[i]->nelem_hevmol]); 00226 frac *= (realnum) COmole[i]->nElem[COmole[i]->nelem_hevmol]; 00227 COmole[i]->xMoleFracMax = MAX2(frac,COmole[i]->xMoleFracMax); 00228 } 00229 } 00230 00231 /* H 21 cm equilibrium timescale, H21cm returns H (not e) collisional 00232 * deexcitation rate (not cs) */ 00233 t = H21cm_H_atom( phycon.te )* dense.xIonDense[ipHYDROGEN][0] + 00234 /* >>chng 02 feb 14, add electron term as per discussion in */ 00235 /* >>refer H1 21cm Liszt, H., 2001, A&A, 371, 698 */ 00236 H21cm_electron( phycon.te )*dense.eden; 00237 00238 /* only update time scale if t is significant */ 00239 if( t > SMALLFLOAT ) 00240 timesc.TimeH21cm = MAX2( 1./t, timesc.TimeH21cm ); 00241 00242 /* remember longest CO timescale */ 00243 if( dense.xIonDense[ipCARBON][0]*dense.xIonDense[ipOXYGEN][0] > SMALLFLOAT ) 00244 { 00245 double timemin; 00246 double product = SDIV(dense.xIonDense[ipCARBON][0]*dense.xIonDense[ipOXYGEN][0]); 00247 struct molecule *co_molecule; 00248 co_molecule = findspecies("CO"); 00249 timemin = MIN2(timesc.AgeCOMoleDest[co_molecule->index] , 00250 timesc.AgeCOMoleDest[co_molecule->index] * co_molecule->hevmol/ product ); 00251 /* this is rate CO is destroyed, equal to formation rate in equilibrium */ 00252 timesc.BigCOMoleForm = MAX2( timesc.BigCOMoleForm , timemin ); 00253 } 00254 00255 /* remember longest H2 destruction timescale timescale */ 00256 timesc.time_H2_Dest_longest = MAX2(timesc.time_H2_Dest_longest, timesc.time_H2_Dest_here ); 00257 00258 /* remember longest H2 formation timescale timescale */ 00259 timesc.time_H2_Form_longest = MAX2( timesc.time_H2_Form_longest , timesc.time_H2_Form_here ); 00260 00261 /* increment counter if this zone possibly thermally unstable 00262 * this flag was set in conv_temp_eden_ioniz.cpp, 00263 * derivative of heating and cooling negative */ 00264 if( thermal.lgUnstable ) 00265 thermal.nUnstable += 1; 00266 00267 /* remember Stromgren radius - where hydrogen ionization falls below half */ 00268 if( !rfield.lgUSphON && dense.xIonDense[ipHYDROGEN][0]/dense.gas_phase[ipHYDROGEN] > 0.49 ) 00269 { 00270 rfield.rstrom = (realnum)radius.Radius; 00271 rfield.lgUSphON = true; 00272 } 00273 00274 /* ratio of inner to outer radii, at this point 00275 * radius is the outer radius of this zone */ 00276 DilutionHere = POW2((radius.Radius - radius.drad*radius.dRadSign)/ 00277 radius.Radius); 00278 00279 if( trace.lgTrace && trace.lgConBug ) 00280 { 00281 fprintf( ioQQQ, " Energy, flux, OTS:\n" ); 00282 for( i=0; i < rfield.nflux; i++ ) 00283 { 00284 fprintf( ioQQQ, "%6ld%10.2e%10.2e%10.2e", i, rfield.anu[i], 00285 rfield.flux[i] + rfield.outlin[i] + rfield.ConInterOut[i], 00286 rfield.otscon[i] + rfield.otslin[i] + rfield.outlin_noplot[i]); 00287 } 00288 fprintf( ioQQQ, "\n" ); 00289 } 00290 00291 /* diffuse continua factor 00292 * if lgSphere then all diffuse radiation is outward (COVRT=1) 00293 * lgSphere not set then COVRT=0.0 */ 00294 00295 /* begin sanity check */ 00296 /* only do this test in debug mode */ 00297 # if !defined(NDEBUG) 00298 lgFlxNeg = false; 00299 for( i=0; i < rfield.nflux; i++ ) 00300 { 00301 if( rfield.flux[i] < 0. ) 00302 { 00303 fprintf( ioQQQ, " radius_increment finds negative intensity in flux.\n" ); 00304 fprintf( ioQQQ, " Intensity, frequency, pointer=%11.3e%11.3e%6ld\n", 00305 rfield.flux[i], rfield.anu[i], i ); 00306 lgFlxNeg = true; 00307 } 00308 if( rfield.otscon[i] < 0. ) 00309 { 00310 fprintf( ioQQQ, " radius_increment finds negative intensity in otscon.\n" ); 00311 fprintf( ioQQQ, " Intensity, frequency, pointer=%11.3e%11.3e%6ld\n", 00312 rfield.otscon[i], rfield.anu[i], i ); 00313 lgFlxNeg = true; 00314 } 00315 if( opac.tmn[i] < 0. ) 00316 { 00317 fprintf( ioQQQ, " radius_increment finds negative tmn.\n" ); 00318 fprintf( ioQQQ, " value, frequency, pointer=%11.3e%11.3e%6ld %4.4s\n", 00319 opac.tmn[i], rfield.anu[i], i, rfield.chLineLabel[i] ); 00320 lgFlxNeg = true; 00321 } 00322 if( rfield.otslin[i] < 0. ) 00323 { 00324 fprintf( ioQQQ, " radius_increment finds negative intensity in otslin.\n" ); 00325 fprintf( ioQQQ, " Intensity, frequency, pointer=%11.3e%11.3e%6ld %4.4s\n", 00326 rfield.otslin[i], rfield.anu[i], i, rfield.chLineLabel[i] ); 00327 lgFlxNeg = true; 00328 } 00329 if( rfield.outlin[i] < 0. ) 00330 { 00331 fprintf( ioQQQ, " radius_increment finds negative intensity in outlin.\n" ); 00332 fprintf( ioQQQ, " Intensity, frequency, pointer=%11.3e%11.3e%6ld %4.4s\n", 00333 rfield.outlin[i], rfield.anu[i], i, rfield.chLineLabel[i] ); 00334 lgFlxNeg = true; 00335 } 00336 if( rfield.ConInterOut[i] < 0. ) 00337 { 00338 fprintf( ioQQQ, " radius_increment finds negative intensity in ConInterOut.\n" ); 00339 fprintf( ioQQQ, " Intensity, frequency, pointer=%11.3e%11.3e%6ld %4.4s\n", 00340 rfield.ConInterOut[i], rfield.anu[i], i, rfield.chContLabel[i] ); 00341 lgFlxNeg = true; 00342 } 00343 if( opac.opacity_abs[i] < 0. ) 00344 { 00345 opac.lgOpacNeg = true; 00346 /* this sub will punch negative opacities on io unit, 00347 * iff 'set negopc' command was given */ 00348 pnegopc(); 00349 } 00350 } 00351 if( lgFlxNeg ) 00352 { 00353 fprintf( ioQQQ, " Insanity has occurred, this is zone%4ld\n", 00354 nzone ); 00355 ShowMe(); 00356 cdEXIT(EXIT_FAILURE); 00357 } 00358 /*end sanity check*/ 00359 # endif 00360 00361 /* rfield.lgOpacityFine flag set false with no fine opacities command 00362 * tests show that always evaluating this changes fast run of 00363 * pn_paris from 26.7 sec to 35.1 sec 00364 * but must always update fine opacities since used for transmission */ 00365 if( rfield.lgOpacityFine ) 00366 { 00367 dradfac = (realnum)radius.drad_x_fillfac; 00368 /* increment the fine opacity array */ 00369 for( i=0; i<rfield.nfine; ++i ) 00370 { 00371 rfield.fine_opt_depth[i] += 00372 rfield.fine_opac_zone[i]*dradfac; 00373 } 00374 00375 /* evaluate fine opacity transmission coefficient for transmission 00376 * through this zone. This is assumed to be a scattering opacity, 00377 * only do this if including scattering */ 00378 if( opac.lgScatON ) 00379 { 00380 /* sum over coarse continuum */ 00381 for( i=0; i < rfield.nflux-1; i++ ) 00382 { 00383 /* find transmission coefficient if lower and upper bounds of coarse continuum is within 00384 * boundaries of fine continuum 00385 * unity is default */ 00386 if( rfield.ipnt_coarse_2_fine[i] && rfield.ipnt_coarse_2_fine[i+1] ) 00387 { 00388 rfield.trans_coef_zone[i] = 0.; 00389 for( j=rfield.ipnt_coarse_2_fine[i]; j<rfield.ipnt_coarse_2_fine[i+1]; ++j ) 00390 { 00391 /* find transmission coefficient for this zone */ 00392 rfield.trans_coef_zone[i] += 1.f / (1.f + rfield.fine_opac_zone[j]*dradfac ); 00393 } 00394 /* first branch is normal case, where fine continuum is finer than 00395 * coarse continuum. But, when end temp is very high, fine continuum is 00396 * very coarse, so may be just one cell, and following will not pass */ 00397 if( rfield.ipnt_coarse_2_fine[i+1]>rfield.ipnt_coarse_2_fine[i] ) 00398 { 00399 rfield.trans_coef_zone[i] /= (rfield.ipnt_coarse_2_fine[i+1]-rfield.ipnt_coarse_2_fine[i]); 00400 } 00401 else 00402 { 00403 /* in case where fine is coarser than coarse, just use firs cell */ 00404 rfield.trans_coef_zone[i] = 00405 1.f / (1.f + rfield.fine_opac_zone[rfield.ipnt_coarse_2_fine[i]]*dradfac ); 00406 } 00407 } 00408 else 00409 { 00410 /* this branch, not within find continuum energy range - transmission unity */ 00411 rfield.trans_coef_zone[i] = 1.f; 00412 } 00413 00414 /* the total is just the current total times the transmission coefficient in this zone. */ 00415 rfield.trans_coef_total[i] *= rfield.trans_coef_zone[i]; 00416 } 00417 } 00418 } 00419 00420 /* electron scattering opacity */ 00421 escatc = 6.65e-25*dense.eden; 00422 00423 /* this loop should not be to <= nflux since we not want to clobber the 00424 * continuum unit integration */ 00425 relec = 0.; 00426 rforce = 0.; 00427 for( i=0; i < rfield.nflux; i++ ) 00428 { 00429 /* sum total continuous optical depths */ 00430 opac.TauAbsGeo[0][i] += (realnum)(opac.opacity_abs[i]*radius.drad_x_fillfac); 00431 opac.TauScatGeo[0][i] += (realnum)(opac.opacity_sct[i]*radius.drad_x_fillfac); 00432 00433 /* following only optical depth to illuminated face */ 00434 opac.TauAbsFace[i] += (realnum)(opac.opacity_abs[i]*radius.drad_x_fillfac); 00435 opac.TauScatFace[i] += (realnum)(opac.opacity_sct[i]*radius.drad_x_fillfac); 00436 00437 /* these are total in inward direction, large if spherical */ 00438 opac.TauTotalGeo[0][i] = opac.TauAbsGeo[0][i] + opac.TauScatGeo[0][i]; 00439 00440 /* TMN is array of scale factors which account for attenuation 00441 * of continuum across the zone (necessary to conserve energy 00442 * at the 1% - 5% level.) sphere effects in 00443 * drNext was set by NEXTDR and will be next dr */ 00444 /* compute both total and Thomson scat rad acceleration */ 00445 rforce += (rfield.flux[i] + rfield.outlin[i] + rfield.outlin_noplot[i]+ 00446 rfield.ConInterOut[i])* rfield.anu[i]*(opac.opacity_abs[i] + 00447 opac.opacity_sct[i]); 00448 00449 relec += ((rfield.flux[i] + rfield.outlin[i] + rfield.outlin_noplot[i]+ 00450 rfield.ConInterOut[i])* escatc)*rfield.anu[i]; 00451 00452 /* attenuation of flux by optical depths IN THIS ZONE 00453 * AngleIllumRadian is 1/COS(theta), is usually 1, reset with illuminate command, 00454 * option for illumination of slab at an angle */ 00455 opac.ExpZone[i] = sexp(opac.opacity_abs[i]*radius.drad_x_fillfac* 00456 geometry.DirectionalCosin); 00457 00458 /* e(-tau) in inward direction, up to illuminated face */ 00459 opac.ExpmTau[i] *= (realnum)opac.ExpZone[i]; 00460 00461 /* e2(tau) in inward direction, up to illuminated face */ 00462 opac.E2TauAbsFace[i] = (realnum)e2(opac.TauAbsFace[i]); 00463 ASSERT( opac.E2TauAbsFace[i] <= 1. && opac.E2TauAbsFace[i] >= 0. ); 00464 00465 /* on second and later iterations define outward E2 */ 00466 if( iteration > 1 ) 00467 { 00468 /* e2 from current position to outer edge of shell */ 00469 realnum tau = MAX2(SMALLFLOAT , opac.TauAbsTotal[i] - opac.TauAbsFace[i] ); 00470 opac.E2TauAbsOut[i] = (realnum)e2( tau ); 00471 ASSERT( opac.E2TauAbsOut[i]>=0. && opac.E2TauAbsOut[i]<=1. ); 00472 } 00473 00474 /* DilutionHere is square of ratio of inner to outer radius */ 00475 opfac = opac.ExpZone[i]*DilutionHere; 00476 ASSERT( opfac <= 1.0 ); 00477 00478 /* >>chng 07 may 22, break continuum into three parts */ 00479 rfield.flux_beam_const[i] *= (realnum)opfac; 00480 rfield.flux_beam_time[i] *= (realnum)opfac; 00481 rfield.flux_isotropic[i] *= (realnum)opfac; 00482 rfield.flux[i] = rfield.flux_beam_const[i] + rfield.flux_beam_time[i] + 00483 rfield.flux_isotropic[i]; 00484 00485 /* update SummedCon here since flux changed */ 00486 rfield.SummedCon[i] = rfield.flux[i] + rfield.SummedDif[i]; 00487 00488 /* outward lines and continua */ 00489 rfield.ConInterOut[i] *= (realnum)opfac; 00490 rfield.outlin[i] *= (realnum)opfac; 00491 rfield.outlin_noplot[i] *= (realnum)opfac; 00492 00493 rfield.ConEmitOut[i] *= (realnum)opfac; 00494 rfield.ConEmitOut[i] += rfield.ConEmitLocal[nzone][i]*(realnum)radius.dVolOutwrd*opac.tmn[i]; 00495 00496 /* set occupation numbers, first attenuated incident continuum */ 00497 rfield.OccNumbIncidCont[i] = rfield.flux[i]*rfield.convoc[i]; 00498 00499 /* >>chng 00 oct 03, add diffuse continua */ 00500 /* local diffuse continua */ 00501 rfield.OccNumbDiffCont[i] = rfield.ConEmitLocal[nzone][i]*rfield.convoc[i]; 00502 00503 /* >>chng 05 feb 16, occupation number of outward diffuse continuum */ 00504 rfield.OccNumbContEmitOut[i] = rfield.ConEmitOut[i]*rfield.convoc[i]; 00505 } 00506 00507 /* begin sanity check, compare total Lyman continuum optical depth 00508 * with amount of extinction there */ 00509 00510 /* this is amount continuum attenuated to illuminated face, 00511 * but only do test if flux positive, not counting scattering opacity, 00512 * and correction for spherical dilution not important */ 00513 /* >>chng 02 dec 05, add test for small float, protect against models 00514 * where we have gone below smallfloat, and so float is ragged */ 00515 if( rfield.flux[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]>SMALLFLOAT && 00516 (rfield.flux[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]/ 00517 SDIV(rfield.flux_total_incident[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]) ) > SMALLFLOAT && 00518 !opac.lgScatON && 00519 radius.depth/radius.Radius < 1e-4 ) 00520 { 00521 /* ratio of current to incident continuum, converted to optical depth */ 00522 /* >>chng 99 apr 23, this crashed on alpha due to underflow to zero in argy 00523 * to log. defended two ways - above checks that ratio of fluxes is large enough, 00524 * and here convert to double. 00525 * error found by Peter van Hoof */ 00526 double tau_effec = -log((double)rfield.flux[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]/ 00527 (double)opac.tmn[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]/ 00528 (double)rfield.flux_total_incident[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]); 00529 00530 /* this is computed absorption optical depth to illuminated face */ 00531 /* >>chng 06 jul 11, add geometry factor for illumination at an angle - bug fix - should have 00532 * always been there */ 00533 double tau_true = opac.TauAbsFace[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]*geometry.DirectionalCosin; 00534 00535 /* first test is relative error, second is to absolute error and comes 00536 * in for very small tau, where differences are in the round-off */ 00537 if( fabs( tau_effec - tau_true ) / t > 0.01 && 00538 /* for very small inner optical depths, the tmn correction is major, 00539 * and this test is not precise */ 00540 fabs(tau_effec-tau_true)>MAX2(0.001,1.-opac.tmn[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]) ) 00541 { 00542 /* in print below must add extra HI col den since this will later be 00543 * incremented in RT_tau_inc */ 00544 fprintf( ioQQQ, 00545 " PROBLEM radius_increment Lyman continuum insanity zone %li, effective tau=%g, atomic tau=%g simple tau=%g\n", 00546 nzone, 00547 tau_effec , 00548 tau_true , 00549 6.327e-18*(dense.xIonDense[ipHYDROGEN][0]*radius.drad_x_fillfac+colden.colden[ipCOL_H0]) ); 00550 TotalInsanity(); 00551 } 00552 } 00553 /* end sanity check */ 00554 00555 /* do scattering opacity, intensity goes as 1/(1+tau) 00556 * scattering opacity is turned off when sphere is set */ 00557 if( opac.lgScatON ) 00558 { 00559 for( i=0; i < rfield.nflux; i++ ) 00560 { 00561 /* assume half forward scattering 00562 * opfac = 1. / ( 1. + dReff*0.5 * scatop(i) ) 00563 * >>chng 97 apr 25, remove .5 to get agreement with 00564 * Lightman and White equation 11 in small epsilon limit, 00565 * >>refer continuum RT Lightman, A.P., & White, T.R. 1988, ApJ, 335, 57 */ 00566 opfac = 1./(1. + radius.drad_x_fillfac*opac.opacity_sct[i]); 00567 ASSERT( opfac <= 1.0 ); 00568 rfield.flux_beam_const[i] *= (realnum)opfac; 00569 rfield.flux_beam_time[i] *= (realnum)opfac; 00570 rfield.flux_isotropic[i] *= (realnum)opfac; 00571 rfield.flux[i] = rfield.flux_beam_const[i] + rfield.flux_beam_time[i] + 00572 rfield.flux_isotropic[i]; 00573 00574 rfield.ConInterOut[i] *= (realnum)opfac; 00575 rfield.ConEmitOut[i] *= (realnum)opfac; 00576 rfield.outlin[i] *= (realnum)opfac; 00577 rfield.outlin_noplot[i] *= (realnum)opfac; 00578 } 00579 } 00580 00581 /* now do slight reshuffling of energy due to Compton scattering */ 00582 cmshft(); 00583 00584 /* remember the largest value */ 00585 wind.AccelMax = (realnum)MAX2(wind.AccelMax,wind.AccelTot); 00586 00587 /* force multiplier; relec can be zero for very low densities - so use double 00588 * form of safe_div - fmul is of order unity - wind.AccelLine and wind.AccelCont 00589 * are both floats to will underflow long before relec will - fmul is only used 00590 * in output, not in any physics */ 00591 relec *= (EN1RYD/SPEEDLIGHT/dense.xMassDensity); 00592 if( relec > SMALLFLOAT ) 00593 wind.fmul = (realnum)( (wind.AccelLine + wind.AccelCont) / relec); 00594 else 00595 wind.fmul = 0.; 00596 00597 /* keep track of average acceleration */ 00598 wind.AccelAver += wind.AccelTot*(realnum)radius.drad_x_fillfac; 00599 wind.acldr += (realnum)radius.drad_x_fillfac; 00600 00601 /* following is integral of radiative force */ 00602 pressure.pinzon = (realnum)(wind.AccelTot*dense.xMassDensity*radius.drad_x_fillfac*geometry.DirectionalCosin); 00603 /*fprintf(ioQQQ," debuggg pinzon %.2f %.2e %.2e %.2e\n", 00604 fnzone,pressure.pinzon,dense.xMassDensity,wind.AccelTot);*/ 00605 pressure.PresInteg += pressure.pinzon; 00606 00607 /* sound is sound travel time, sqrt term is sound speed */ 00608 timesc.sound_speed_isothermal = sqrt(pressure.PresGasCurr/dense.xMassDensity); 00609 /* adiabatic sound speed assuming mono-atomic gas - gamma is 5/3*/ 00610 timesc.sound_speed_adiabatic = sqrt(5.*pressure.PresGasCurr/(3.*dense.xMassDensity) ); 00611 timesc.sound += radius.drad_x_fillfac / timesc.sound_speed_isothermal; 00612 00613 /* attenuate neutrons if they are present */ 00614 if( hextra.lgNeutrnHeatOn ) 00615 { 00616 /* correct for optical depth effects */ 00617 hextra.totneu *= (realnum)sexp(hextra.CrsSecNeutron*dense.gas_phase[ipHYDROGEN]*radius.drad_x_fillfac*geometry.DirectionalCosin); 00618 /* correct for spherical effects */ 00619 hextra.totneu *= (realnum)DilutionHere; 00620 } 00621 00622 /* following radiation factors are extinguished by 1/r**2 and e- opacity 00623 * also bound electrons */ 00624 00625 /* do all emergent spectrum from illuminated face if model is NOT spherical */ 00626 if( !geometry.lgSphere ) 00627 { 00628 double Reflec_Diffuse_Cont; 00629 00630 /* >>chng 01 jul 14, from lower limit of 0 to plasma frequency - 00631 * never should have added diffuse emission from below the plasma frequency */ 00632 for( i=rfield.ipPlasma-1; i < rfield.nflux; i++ ) 00633 { 00634 if( opac.TauAbsGeo[0][i] < 30. ) 00635 { 00636 /* ConEmitLocal is diffuse emission per unit vol, fill factor 00637 * the 1/2 comes from isotropic emission */ 00638 Reflec_Diffuse_Cont = rfield.ConEmitLocal[nzone][i]/2.* 00639 radius.drad_x_fillfac * opac.E2TauAbsFace[i]*radius.r1r0sq; 00640 00641 /* ConEmitReflec - reflected diffuse continuum */ 00642 rfield.ConEmitReflec[i] += (realnum)(Reflec_Diffuse_Cont); 00643 00644 /* the reflected part of the incident continuum */ 00645 rfield.ConRefIncid[i] += (realnum)(rfield.flux[i]*opac.opacity_sct[i]* 00646 radius.drad_x_fillfac*opac.E2TauAbsFace[i]*radius.r1r0sq); 00647 00648 /* reflected line emission */ 00649 rfield.reflin[i] += (realnum)(rfield.outlin[i]*opac.opacity_sct[i]* 00650 radius.drad_x_fillfac*opac.E2TauAbsFace[i]*radius.r1r0sq); 00651 } 00652 } 00653 } 00654 00655 /* following is general method to find means weighted by various functions 00656 * called in IterStart to initialize to zero, called here to put in numbers 00657 * results will be weighted by radius and volume 00658 * this is the only place things must be entered to create averages 00659 * code is in mean.c */ 00660 aver("zone",1.,1.," "); 00661 aver("doit",phycon.te,1.," Te "); 00662 aver("doit",phycon.te,dense.eden," Te(Ne) "); 00663 aver("doit",phycon.te,dense.eden*dense.xIonDense[ipHYDROGEN][1]," Te(NeNp) "); 00664 aver("doit",phycon.te,dense.eden*dense.xIonDense[ipHELIUM][1] ," Te(NeHe+)"); 00665 aver("doit",phycon.te,dense.eden*dense.xIonDense[ipHELIUM][2] ,"Te(NeHe2+)"); 00666 aver("doit",phycon.te,dense.eden*dense.xIonDense[ipOXYGEN][1] ," Te(NeO+) " ); 00667 aver("doit",phycon.te,dense.eden*dense.xIonDense[ipOXYGEN][2] ," Te(NeO2+)"); 00668 /*aver("doit",phycon.te,hmi.Hmolec[ipMH2g]," Te(H2) ");*/ 00669 aver("doit",phycon.te,hmi.H2_total," Te(H2) "); 00670 aver("doit",dense.gas_phase[ipHYDROGEN],1.," N(H) "); 00671 aver("doit",dense.eden,dense.xIonDense[ipOXYGEN][2]," Ne(O2+) "); 00672 aver("doit",dense.eden,dense.xIonDense[ipHYDROGEN][1]," Ne(Np) "); 00673 00674 /* save information about structure of model, now used to get t^2 */ 00675 /* max because if program aborts during search phase, will get to here 00676 * with nzone = -1 */ 00677 nzone_minus_1 = MAX2( 0, nzone-1 ); 00678 /* this is number of struc.xx zones with valid data */ 00679 struc.nzone = nzone_minus_1+1; 00680 ASSERT(nzone_minus_1>=0 && nzone_minus_1 < struc.nzlim ); 00681 struc.testr[nzone_minus_1] = (realnum)phycon.te; 00682 /* number of particles per unit vol */ 00683 struc.DenParticles[nzone_minus_1] = dense.pden; 00684 /* >>chng 02 May 2001 rjrw: add hden for dilution */ 00685 struc.hden[nzone_minus_1] = (realnum)dense.gas_phase[ipHYDROGEN]; 00686 /* total grams per unit vol */ 00687 struc.DenMass[nzone_minus_1] = dense.xMassDensity; 00688 struc.heatstr[nzone_minus_1] = thermal.htot; 00689 struc.coolstr[nzone_minus_1] = thermal.ctot; 00690 struc.volstr[nzone_minus_1] = (realnum)radius.dVeff; 00691 struc.drad[nzone_minus_1] = (realnum)radius.drad; 00692 struc.drad_x_fillfac[nzone_minus_1] = (realnum)radius.drad_x_fillfac; 00693 struc.histr[nzone_minus_1] = dense.xIonDense[ipHYDROGEN][0]; 00694 struc.hiistr[nzone_minus_1] = dense.xIonDense[ipHYDROGEN][1]; 00695 struc.ednstr[nzone_minus_1] = (realnum)dense.eden; 00696 struc.o3str[nzone_minus_1] = dense.xIonDense[ipOXYGEN][2]; 00697 struc.pressure[nzone_minus_1] = (realnum)pressure.PresTotlCurr; 00698 struc.windv[nzone_minus_1] = (realnum)wind.windv; 00699 struc.AccelTot[nzone_minus_1] = wind.AccelTot; 00700 struc.AccelGravity[nzone_minus_1] = wind.AccelGravity; 00701 struc.pres_radiation_lines_curr[nzone_minus_1] = (realnum)pressure.pres_radiation_lines_curr; 00702 struc.GasPressure[nzone_minus_1] = (realnum)pressure.PresGasCurr; 00703 struc.depth[nzone_minus_1] = (realnum)radius.depth; 00704 /* save absorption optical depth from illuminated face to current position */ 00705 struc.xLyman_depth[nzone_minus_1] = opac.TauAbsFace[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]]; 00706 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem) 00707 { 00708 struc.gas_phase[nelem][nzone_minus_1] = dense.gas_phase[nelem]; 00709 for( ion=0; ion<nelem+2; ++ion ) 00710 { 00711 struc.xIonDense[nelem][ion][nzone_minus_1] = dense.xIonDense[nelem][ion]; 00712 } 00713 } 00714 00715 /* the hydrogen molecules */ 00716 for(mol=0;mol<N_H_MOLEC;mol++) 00717 { 00718 struc.H2_molec[mol][nzone_minus_1] = hmi.Hmolec[mol]; 00719 } 00720 00721 /* the heavy element molecules */ 00722 for(mol=0;mol<mole.num_comole_calc;mol++) 00723 { 00724 struc.CO_molec[mol][nzone_minus_1] = COmole[mol]->hevmol; 00725 } 00726 00727 colden.dlnenp += dense.eden*(double)(dense.xIonDense[ipHYDROGEN][1])*radius.drad_x_fillfac; 00728 colden.dlnenHep += dense.eden*(double)(dense.xIonDense[ipHELIUM][1])*radius.drad_x_fillfac; 00729 colden.dlnenHepp += dense.eden*(double)(dense.xIonDense[ipHELIUM][2])*radius.drad_x_fillfac; 00730 colden.dlnenCp += dense.eden*(double)(dense.xIonDense[ipCARBON][1])*radius.drad_x_fillfac; 00731 00732 /* >>chng 05 jan 09, add integral of n(H0) / Tspin */ 00733 colden.H0_ov_Tspin += (double)(dense.xIonDense[ipHYDROGEN][0]) / hyperfine.Tspin21cm*radius.drad_x_fillfac; 00734 00735 /* >>chng 05 Mar 07, add integral of n(OH) / Tspin */ 00736 colden.OH_ov_Tspin += (double)(findspecies("OH")->hevmol) / phycon.te*radius.drad_x_fillfac; 00737 00738 /* this is Lya excitation temperature */ 00739 hydro.TexcLya = (realnum)TexcLine( &Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s] ); 00740 00741 /* count number of times Lya excitation temp hotter than gas */ 00742 if( hydro.TexcLya > phycon.te ) 00743 { 00744 hydro.nLyaHot += 1; 00745 if( hydro.TexcLya > hydro.TLyaMax ) 00746 { 00747 hydro.TLyaMax = hydro.TexcLya; 00748 hydro.TeLyaMax = (realnum)phycon.te; 00749 hydro.nZTLaMax = nzone; 00750 } 00751 } 00752 00753 /* column densities in various species */ 00754 colden.colden[ipCOL_HTOT] += (realnum)(dense.gas_phase[ipHYDROGEN]*radius.drad_x_fillfac); 00755 ASSERT( colden.colden[ipCOL_HTOT] >SMALLFLOAT ); 00756 colden.colden[ipCOL_HMIN] += hmi.Hmolec[ipMHm]*(realnum)radius.drad_x_fillfac; 00757 /* >>chng 02 sep 20, from htwo to H2_total */ 00758 /* >>chng 05 mar 14, rather than H2_total, give H2g and H2s */ 00759 colden.colden[ipCOL_H2g] += hmi.Hmolec[ipMH2g]*(realnum)radius.drad_x_fillfac; 00760 colden.colden[ipCOL_H2s] += hmi.Hmolec[ipMH2s]*(realnum)radius.drad_x_fillfac; 00761 /* this is a special form of column density - should be proportional to total shielding */ 00762 colden.coldenH2_ov_vel += hmi.H2_total*(realnum)radius.drad_x_fillfac / DoppVel.doppler[LIMELM+2]; 00763 colden.colden[ipCOL_HeHp] += hmi.Hmolec[ipMHeHp]*(realnum)radius.drad_x_fillfac; 00764 colden.colden[ipCOL_H2p] += hmi.Hmolec[ipMH2p]*(realnum)radius.drad_x_fillfac; 00765 colden.colden[ipCOL_H3p] += hmi.Hmolec[ipMH3p]*(realnum)radius.drad_x_fillfac; 00766 colden.colden[ipCOL_Hp] += dense.xIonDense[ipHYDROGEN][1]*(realnum)radius.drad_x_fillfac; 00767 colden.colden[ipCOL_H0] += dense.xIonDense[ipHYDROGEN][0]*(realnum)radius.drad_x_fillfac; 00768 00769 colden.colden[ipCOL_elec] += (realnum)(dense.eden*radius.drad_x_fillfac); 00770 /* He I t2S column density*/ 00771 colden.He123S += dense.xIonDense[ipHELIUM][1]* 00772 StatesElem[ipHE_LIKE][ipHELIUM][ipHe2s3S].Pop*(realnum)radius.drad_x_fillfac; 00773 00774 /* the ortho and para column densities */ 00775 h2.ortho_colden += h2.ortho_density*radius.drad_x_fillfac; 00776 h2.para_colden += h2.para_density*radius.drad_x_fillfac; 00777 ASSERT( fabs( h2.ortho_density + h2.para_density - hmi.H2_total ) / SDIV( hmi.H2_total ) < 1e-4 ); 00778 /*fprintf(ioQQQ,"DEBUG ortho para\t%.3e\t%.3e\ttot\t%.3e\t or pa colden\t%.3e\t%.3e\n", 00779 h2.ortho_density, h2.para_density,hmi.H2_total, 00780 h2.ortho_colden , h2.para_colden);*/ 00781 00782 /*>>chng 27mar, GS, Column density of F=0 and F=1 levels of H0*/ 00783 colden.H0_21cm_upper += (HFLines[0].Hi->Pop*radius.drad_x_fillfac); 00784 colden.H0_21cm_lower +=(HFLines[0].Lo->Pop*radius.drad_x_fillfac); 00785 /*fprintf(ioQQQ,"DEBUG pophi-poplo\t%.3e\t%.3e\radius\t%.3e\t col_hi\t%.3e\t%.3e\n", 00786 HFLines[0].PopHi, HFLines[0].PopLo, radius.drad_x_fillfac, 00787 HFLines[0].PopHi*radius.drad_x_fillfac,colden.H0_21cm_upper );*/ 00788 00789 /* the CII and SiII atoms are resolved */ 00790 for( i=0; i < 5; i++ ) 00791 { 00792 /* pops and column density for C II atom */ 00793 colden.C2Colden[i] += colden.C2Pops[i]*(realnum)radius.drad_x_fillfac; 00794 /* pops and column density for SiII atom */ 00795 colden.Si2Colden[i] += colden.Si2Pops[i]*(realnum)radius.drad_x_fillfac; 00796 } 00797 for( i=0; i < 3; i++ ) 00798 { 00799 /* pops and column density for CI atom */ 00800 colden.C1Colden[i] += colden.C1Pops[i]*(realnum)radius.drad_x_fillfac; 00801 /* pops and column density for OI atom */ 00802 colden.O1Colden[i] += colden.O1Pops[i]*(realnum)radius.drad_x_fillfac; 00803 } 00804 for( i=0; i < 4; i++ ) 00805 { 00806 /* pops and column density for CIII atom */ 00807 colden.C3Colden[i] += colden.C3Pops[i]*(realnum)radius.drad_x_fillfac; 00808 } 00809 00810 /* now add total molecular column densities */ 00811 molcol("ADD ",ioQQQ); 00812 00813 /* increment forming the mean ionization and temperature */ 00814 MeanInc(); 00815 00816 /*-----------------------------------------------------------------------*/ 00817 00818 /* calculate average atomic weight per hydrogen of the plasma */ 00819 avWeight = 0.; 00820 for( nelem=0; nelem < LIMELM; nelem++ ) 00821 { 00822 avWeight += dense.gas_phase[nelem]*dense.AtomicWeight[nelem]; 00823 } 00824 avWeight /= dense.gas_phase[ipHYDROGEN]; 00825 00826 /* compute some average grain properties */ 00827 rfield.opac_mag_B_point = 0.; 00828 rfield.opac_mag_V_point = 0.; 00829 rfield.opac_mag_B_extended = 0.; 00830 rfield.opac_mag_V_extended = 0.; 00831 long int nd; 00832 for( nd=0; nd < gv.nBin; nd++ ) 00833 { 00834 00835 /* this is total extinction in magnitudes at V and B, for a point source 00836 * total absorption and scattering, 00837 * does not discount forward scattering to be similar to stellar extinction 00838 * measurements made within ism */ 00839 rfield.extin_mag_B_point += (gv.bin[nd]->dstab1[rfield.ipB_filter-1] + 00840 gv.bin[nd]->pure_sc1[rfield.ipB_filter-1])*gv.bin[nd]->dstAbund* 00841 radius.drad_x_fillfac*dense.gas_phase[ipHYDROGEN] * OPTDEP2EXTIN; 00842 rfield.opac_mag_B_point += (gv.bin[nd]->dstab1[rfield.ipB_filter-1] + 00843 gv.bin[nd]->pure_sc1[rfield.ipB_filter-1])*gv.bin[nd]->dstAbund* 00844 dense.gas_phase[ipHYDROGEN] * OPTDEP2EXTIN; 00845 00846 rfield.extin_mag_V_point += (gv.bin[nd]->dstab1[rfield.ipV_filter-1] + 00847 gv.bin[nd]->pure_sc1[rfield.ipV_filter-1])*gv.bin[nd]->dstAbund* 00848 radius.drad_x_fillfac*dense.gas_phase[ipHYDROGEN] * OPTDEP2EXTIN; 00849 00850 rfield.opac_mag_V_point += (gv.bin[nd]->dstab1[rfield.ipV_filter-1] + 00851 gv.bin[nd]->pure_sc1[rfield.ipV_filter-1])*gv.bin[nd]->dstAbund* 00852 dense.gas_phase[ipHYDROGEN] * OPTDEP2EXTIN; 00853 00854 /* this is total extinction in magnitudes at V and B, for an extended source 00855 * total absorption and scattering, 00856 * DOES discount forward scattering to apply for extended source like Orion */ 00857 rfield.extin_mag_B_extended += (gv.bin[nd]->dstab1[rfield.ipB_filter-1] + 00858 gv.bin[nd]->pure_sc1[rfield.ipB_filter]*gv.bin[nd]->asym[rfield.ipB_filter-1] )*gv.bin[nd]->dstAbund* 00859 radius.drad_x_fillfac*dense.gas_phase[ipHYDROGEN] * OPTDEP2EXTIN; 00860 rfield.opac_mag_B_extended += (gv.bin[nd]->dstab1[rfield.ipB_filter-1] + 00861 gv.bin[nd]->pure_sc1[rfield.ipB_filter]*gv.bin[nd]->asym[rfield.ipB_filter-1] )*gv.bin[nd]->dstAbund* 00862 dense.gas_phase[ipHYDROGEN] * OPTDEP2EXTIN; 00863 00864 rfield.extin_mag_V_extended += (gv.bin[nd]->dstab1[rfield.ipV_filter-1] + 00865 gv.bin[nd]->pure_sc1[rfield.ipV_filter]*gv.bin[nd]->asym[rfield.ipV_filter-1] )*gv.bin[nd]->dstAbund* 00866 radius.drad_x_fillfac*dense.gas_phase[ipHYDROGEN] * OPTDEP2EXTIN; 00867 rfield.opac_mag_V_extended += (gv.bin[nd]->dstab1[rfield.ipV_filter-1] + 00868 gv.bin[nd]->pure_sc1[rfield.ipV_filter]*gv.bin[nd]->asym[rfield.ipV_filter-1] )*gv.bin[nd]->dstAbund* 00869 dense.gas_phase[ipHYDROGEN] * OPTDEP2EXTIN; 00870 00871 gv.bin[nd]->avdust += gv.bin[nd]->tedust*(realnum)radius.drad_x_fillfac; 00872 gv.bin[nd]->avdft += gv.bin[nd]->DustDftVel*(realnum)radius.drad_x_fillfac; 00873 gv.bin[nd]->avdpot += (realnum)(gv.bin[nd]->dstpot*EVRYD*radius.drad_x_fillfac); 00874 gv.bin[nd]->avDGRatio += (realnum)(gv.bin[nd]->dustp[1]*gv.bin[nd]->dustp[2]* 00875 gv.bin[nd]->dustp[3]*gv.bin[nd]->dustp[4]*gv.bin[nd]->dstAbund/avWeight* 00876 radius.drad_x_fillfac); 00877 } 00878 00879 /* there are some quantities needed to calculation the Jeans mass and radius */ 00880 colden.TotMassColl += dense.xMassDensity*(realnum)radius.drad_x_fillfac; 00881 colden.tmas += (realnum)phycon.te*dense.xMassDensity*(realnum)radius.drad_x_fillfac; 00882 colden.wmas += dense.wmole*dense.xMassDensity*(realnum)radius.drad_x_fillfac; 00883 00884 /* now find minimum Jeans length and mass; length in cm */ 00885 rjeans = 7.79637 + (phycon.alogte - log10(dense.wmole) - log10(dense.xMassDensity* 00886 geometry.FillFac))/2.; 00887 00888 /* minimum Jeans mass in gm */ 00889 ajmass = 3.*(rjeans - 0.30103) + log10(4.*PI/3.*dense.xMassDensity* 00890 geometry.FillFac); 00891 00892 /* now remember smallest */ 00893 colden.rjnmin = MIN2(colden.rjnmin,(realnum)rjeans); 00894 colden.ajmmin = MIN2(colden.ajmmin,(realnum)ajmass); 00895 00896 if( trace.lgTrace ) 00897 { 00898 fprintf( ioQQQ, " radius_increment returns\n" ); 00899 } 00900 return; 00901 } 00902 00903 #if !defined(NDEBUG) 00904 /*pnegopc punch negative opacities on io unit, iff 'set negopc' command was given */ 00905 STATIC void pnegopc(void) 00906 { 00907 long int i; 00908 FILE *ioFile; 00909 00910 DEBUG_ENTRY( "pnegopc()" ); 00911 00912 if( opac.lgNegOpacIO ) 00913 { 00914 /* option to punch negative opacities */ 00915 ioFile = open_data( "negopc.txt", "w", AS_LOCAL_ONLY ); 00916 for( i=0; i < rfield.nflux; i++ ) 00917 { 00918 fprintf( ioFile, "%10.2e %10.2e \n", rfield.anu[i], 00919 opac.opacity_abs[i] ); 00920 } 00921 fclose( ioFile); 00922 } 00923 return; 00924 } 00925 #endif 00926 /* Note - when this file is compiled in fast optimized mode, 00927 * a good compiler will complain that the routine pnegopc is a 00928 * static, local routine, but that it has not been used. It is 00929 * indeed used , but only when NDEBUG is not set, so it is only 00930 * used when the code is compiled in debug mode. So this is 00931 * not a problem. */ 00932 00933 00934 /*cmshft - so Compton scattering shift of spectrum */ 00935 STATIC void cmshft(void) 00936 { 00937 long int i; 00938 00939 DEBUG_ENTRY( "cmshft()" ); 00940 00941 /* first check whether Compton scattering is in as heat/cool */ 00942 if( rfield.comoff == 0. ) 00943 { 00944 return; 00945 } 00946 00947 if( rfield.comoff != 0. ) 00948 { 00949 return; 00950 } 00951 00952 /* do reshuffle */ 00953 for( i=0; i < rfield.nflux; i++ ) 00954 { 00955 continue; 00956 /* watch this space for some really great code!!!! 00957 * COMUP needs factor of TE to be Compton cooling */ 00958 } 00959 return; 00960 }