00001
00002
00003
00004
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 "colden.h"
00012 #include "geometry.h"
00013 #include "opacity.h"
00014 #include "thermal.h"
00015 #include "doppvel.h"
00016 #include "dense.h"
00017 #include "mole.h"
00018 #include "h2.h"
00019 #include "timesc.h"
00020 #include "neutrn.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 "rt.h"
00035 #include "radius.h"
00036
00037 static void cmshft(void);
00038
00039 #if !defined(NDEBUG)
00040
00041
00042 static void pnegopc(void);
00043 #endif
00044
00045 void radius_increment(void)
00046 {
00047 # if !defined(NDEBUG)
00048 bool lgFlxNeg;
00049 # endif
00050
00051 long int nelem,
00052 i,
00053 j ,
00054 ion,
00055 nzone_minus_1 ,
00056 mol;
00057 double
00058 avWeight,
00059 DilutionHere ,
00060 escatc,
00061 fmol,
00062 opfac,
00063 relec,
00064 rforce,
00065 t;
00066
00067 double ajmass,
00068 Error,
00069 rjeans;
00070
00071 float dradfac;
00072
00073 DEBUG_ENTRY( "radius_increment()" );
00074
00075
00076
00077 if( lgAbort )
00078 {
00079 DEBUG_EXIT( "radius_increment()" );
00080 return;
00081 }
00082
00083
00084 # if !defined(NDEBUG)
00085 if( !dynamics.lgAdvection )
00086 {
00087 long int ipISO;
00088 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00089 {
00090 for( nelem=ipISO; nelem<LIMELM; ++nelem )
00091 {
00092
00093
00094
00095 if( dense.lgElmtOn[nelem] && dense.xIonDense[nelem][nelem-ipISO]/dense.eden>conv.EdenErrorAllowed / 10. &&
00096
00097 !(iso.chTypeAtomUsed[ipISO][nelem][0]=='L' && iso.xIonSimple[ipISO][nelem]<1e-10) )
00098 {
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110 # define ERR_CHK 1.001
00111 double err_chk = ERR_CHK;
00112
00113
00114
00115 if( iso.chTypeAtomUsed[ipISO][nelem][0]=='L' )
00116 err_chk *= 10.;
00117
00118 if( iso.Pop2Ion[ipISO][nelem][0] >
00119 dense.xIonDense[nelem][nelem-ipISO]/(double)SDIV(dense.xIonDense[nelem][nelem+1-ipISO])*err_chk &&
00120
00121 iso.pop_ion_ov_neut[ipISO][nelem] > 1e-15 )
00122 {
00123 fprintf(ioQQQ,
00124 " DISASTER radius_increment finds inconsistent populations, Pop2Ion zn %.2f ipISO %li nelem %li %.4e %.4e solver:%s\n",
00125 fnzone,
00126 ipISO , nelem ,
00127 iso.Pop2Ion[ipISO][nelem][0]*dense.xIonDense[nelem][nelem+1-ipISO] ,
00128 dense.xIonDense[nelem][nelem-ipISO] ,
00129 iso.chTypeAtomUsed[ipISO][nelem]);
00130 fprintf(ioQQQ," level solver %s found pop_ion_ov_neut of %.5e",
00131 iso.chTypeAtomUsed[ipISO][nelem] ,
00132 iso.pop_ion_ov_neut[ipISO][nelem] );
00133 fprintf(ioQQQ," ion_solver found pop_ion_ov_neut of %.5e\n",
00134 dense.xIonDense[nelem][nelem+1-ipISO]/SDIV( dense.xIonDense[nelem][nelem-ipISO] ) );
00135 fprintf(ioQQQ,
00136 "simple %.3e Test shows Pop2Ion (%.6e) > xIonDense[nelem]/[nelem+1] (%.6e) \n",
00137 iso.xIonSimple[ipISO][nelem],
00138 iso.Pop2Ion[ipISO][nelem][0] ,
00139 dense.xIonDense[nelem][nelem-ipISO]/(double)SDIV(dense.xIonDense[nelem][nelem+1-ipISO]) );
00140 }
00141
00142
00143
00144
00145 ASSERT( !conv.lgSearch );
00146 ASSERT(
00147 iso.Pop2Ion[ipISO][nelem][0] <=
00148 dense.xIonDense[nelem][nelem-ipISO]/(double)SDIV(dense.xIonDense[nelem][nelem+1-ipISO])*err_chk ||
00149
00150 iso.pop_ion_ov_neut[ipISO][nelem] < 1e-15);
00151
00152
00153 if(
00154 EmisLines[ipISO][nelem][iso.nLyaLevel[ipISO]][0].PopLo*dense.xIonDense[nelem][nelem+1-ipISO] >
00155 dense.xIonDense[nelem][nelem-ipISO]*err_chk &&
00156
00157 iso.pop_ion_ov_neut[ipISO][nelem] > 1e-15)
00158 {
00159 fprintf(ioQQQ," DISASTER PopLo fnz %.2f ipISO %li nelem %li pop_ion_ov_neut %.2e 2pops %.6e %.6e solver:%s\n",
00160 fnzone,
00161 ipISO , nelem ,
00162 iso.pop_ion_ov_neut[ipISO][nelem] ,
00163 EmisLines[ipISO][nelem][iso.nLyaLevel[ipISO]][0].PopLo*dense.xIonDense[nelem][nelem+1-ipISO],
00164 dense.xIonDense[nelem][nelem-ipISO]*err_chk ,
00165 iso.chTypeAtomUsed[ipISO][nelem]);
00166 }
00167 ASSERT(
00168
00169 (EmisLines[ipISO][nelem][iso.nLyaLevel[ipISO]][0].PopLo*dense.xIonDense[nelem][nelem+1-ipISO] <=
00170 dense.xIonDense[nelem][nelem-ipISO]*err_chk) ||
00171
00172 iso.pop_ion_ov_neut[ipISO][nelem] < 1e-15);
00173 # undef ERR_CHK
00174 }
00175 }
00176 }
00177 }
00178 # endif
00179
00180 if( trace.lgTrace )
00181 {
00182 fprintf( ioQQQ,
00183 " radius_increment called; radius=%10.3e rinner=%10.3e DRAD=%10.3e drNext=%10.3e ROUTER=%10.3e DEPTH=%10.3e\n",
00184 radius.Radius, radius.rinner, radius.drad, radius.drNext,
00185 radius.router[iteration-1], radius.depth );
00186 }
00187
00188
00189 Error = fabs(dense.eden - dense.EdenTrue)/SDIV(dense.EdenTrue);
00190 if( Error > conv.BigEdenError )
00191 {
00192 conv.BigEdenError = (float)Error;
00193 dense.nzEdenBad = nzone;
00194 }
00195 conv.AverEdenError += (float)Error;
00196
00197
00198 Error = fabs(thermal.ctot - thermal.htot) / thermal.ctot;
00199 conv.BigHeatCoolError = MAX2((float)Error , conv.BigHeatCoolError );
00200 conv.AverHeatCoolError += (float)Error;
00201
00202
00203 Error = fabs(pressure.PresTotlCurr - pressure.PresTotlCorrect) / pressure.PresTotlCorrect;
00204 conv.BigPressError = MAX2((float)Error , conv.BigPressError );
00205 conv.AverPressError += (float)Error;
00206
00207
00208 dense.xMassTotal += dense.xMassDensity * (float)(radius.drad_x_fillfac*radius.r1r0sq);
00209
00210
00211 timesc.time_therm_long = MAX2(timesc.time_therm_long,1.5*dense.pden*BOLTZMANN*phycon.te/
00212 thermal.ctot);
00213
00214
00215
00216 ASSERT( N_H_MOLEC == 8);
00217 fmol = (hmi.Hmolec[ipMHm] + 2.*(hmi.H2_total + hmi.Hmolec[ipMH2p]))/dense.gas_phase[ipHYDROGEN];
00218
00219
00220 hmi.BiggestH2 = MAX2(hmi.BiggestH2,(float)fmol);
00221
00222
00223 for( i=0; i<mole.num_comole_calc; ++i )
00224 {
00225 if( dense.lgElmtOn[COmole[i]->nelem_hevmol] )
00226 {
00227 float frac = COmole[i]->hevmol / SDIV(dense.gas_phase[COmole[i]->nelem_hevmol]);
00228
00229 frac *= (float) COmole[i]->nElem[COmole[i]->nelem_hevmol];
00230 #if 0
00231
00232 if( COmole[i]->nelem_hevmol==ipOXYGEN && (i==ipO2 || i==ipO2P || i==ipNO2 || i==ipNO2P) )
00233 frac *= 2.f;
00234 else if( COmole[i]->nelem_hevmol==ipNITROGEN && ( i==ipN2O || i==ipN2 ) )
00235 frac *= 2.f;
00236 else if( COmole[i]->nelem_hevmol==ipSULPHUR && (i==ipS2||i==ipS2P ) )
00237 frac *= 2.f;
00238 else if( COmole[i]->nelem_hevmol==ipCARBON && (i==ipC2||i==ipC2P ) )
00239 frac *= 2.f;
00240 #endif
00241 COmole[i]->xMoleFracMax = MAX2(frac,COmole[i]->xMoleFracMax);
00242 }
00243 }
00244
00245
00246
00247 t = H21cm_H_atom( phycon.te )* dense.xIonDense[ipHYDROGEN][0] +
00248
00249
00250 H21cm_electron( phycon.te )*dense.eden;
00251
00252 if( t > SMALLFLOAT )
00253 {
00254 timesc.TimeH21cm = MAX2( 1./t, timesc.TimeH21cm );
00255 }
00256
00257
00258 if( dense.xIonDense[ipCARBON][0]*dense.xIonDense[ipOXYGEN][0] > SMALLFLOAT )
00259 {
00260 double timemin;
00261 double product = SDIV(dense.xIonDense[ipCARBON][0]*dense.xIonDense[ipOXYGEN][0]);
00262 struct molecule *co_molecule;
00263 co_molecule = findspecies("CO");
00264 timemin = MIN2(timesc.AgeCOMoleDest[co_molecule->index] ,
00265 timesc.AgeCOMoleDest[co_molecule->index] * co_molecule->hevmol/ product );
00266
00267 timesc.BigCOMoleForm = MAX2( timesc.BigCOMoleForm , timemin );
00268 }
00269
00270
00271 timesc.time_H2_Dest_longest = MAX2(timesc.time_H2_Dest_longest, timesc.time_H2_Dest_here );
00272
00273
00274 timesc.time_H2_Form_longest = MAX2( timesc.time_H2_Form_longest , timesc.time_H2_Form_here );
00275
00276
00277
00278 if( thermal.lgUnstable )
00279 {
00280 thermal.nUnstable += 1;
00281 }
00282
00283
00284 if( !rfield.lgUSphON && dense.xIonDense[ipHYDROGEN][0]/dense.gas_phase[ipHYDROGEN] > 0.49 )
00285 {
00286 rfield.rstrom = (float)radius.Radius;
00287 rfield.lgUSphON = true;
00288 }
00289
00290
00291
00292 DilutionHere = POW2((radius.Radius - radius.drad*radius.dRadSign)/
00293 radius.Radius);
00294
00295 relec = 0.;
00296 rforce = 0.;
00297
00298 if( trace.lgTrace && trace.lgConBug )
00299 {
00300 fprintf( ioQQQ, " Energy, flux, OTS:\n" );
00301 for( i=0; i < rfield.nflux; i++ )
00302 {
00303 fprintf( ioQQQ, "%6ld%10.2e%10.2e%10.2e", i, rfield.anu[i],
00304 rfield.flux[i] + rfield.outlin[i] + rfield.ConInterOut[i],
00305 rfield.otscon[i] + rfield.otslin[i] + rfield.outlin_noplot[i]);
00306 }
00307 fprintf( ioQQQ, "\n" );
00308 }
00309
00310
00311
00312
00313
00314
00315
00316 # if !defined(NDEBUG)
00317 lgFlxNeg = false;
00318 for( i=0; i < rfield.nflux; i++ )
00319 {
00320 if( rfield.flux[i] < 0. )
00321 {
00322 fprintf( ioQQQ, " radius_increment finds negative intensity in flux.\n" );
00323 fprintf( ioQQQ, " Intensity, frequency, pointer=%11.3e%11.3e%6ld\n",
00324 rfield.flux[i], rfield.anu[i], i );
00325 lgFlxNeg = true;
00326 }
00327 if( rfield.otscon[i] < 0. )
00328 {
00329 fprintf( ioQQQ, " radius_increment finds negative intensity in otscon.\n" );
00330 fprintf( ioQQQ, " Intensity, frequency, pointer=%11.3e%11.3e%6ld\n",
00331 rfield.otscon[i], rfield.anu[i], i );
00332 lgFlxNeg = true;
00333 }
00334 if( opac.tmn[i] < 0. )
00335 {
00336 fprintf( ioQQQ, " radius_increment finds negative tmn.\n" );
00337 fprintf( ioQQQ, " value, frequency, pointer=%11.3e%11.3e%6ld %4.4s\n",
00338 opac.tmn[i], rfield.anu[i], i, rfield.chLineLabel[i] );
00339 lgFlxNeg = true;
00340 }
00341 if( rfield.otslin[i] < 0. )
00342 {
00343 fprintf( ioQQQ, " radius_increment finds negative intensity in otslin.\n" );
00344 fprintf( ioQQQ, " Intensity, frequency, pointer=%11.3e%11.3e%6ld %4.4s\n",
00345 rfield.otslin[i], rfield.anu[i], i, rfield.chLineLabel[i] );
00346 lgFlxNeg = true;
00347 }
00348 if( rfield.outlin[i] < 0. )
00349 {
00350 fprintf( ioQQQ, " radius_increment finds negative intensity in outlin.\n" );
00351 fprintf( ioQQQ, " Intensity, frequency, pointer=%11.3e%11.3e%6ld %4.4s\n",
00352 rfield.outlin[i], rfield.anu[i], i, rfield.chLineLabel[i] );
00353 lgFlxNeg = true;
00354 }
00355 if( rfield.ConInterOut[i] < 0. )
00356 {
00357 fprintf( ioQQQ, " radius_increment finds negative intensity in ConInterOut.\n" );
00358 fprintf( ioQQQ, " Intensity, frequency, pointer=%11.3e%11.3e%6ld %4.4s\n",
00359 rfield.ConInterOut[i], rfield.anu[i], i, rfield.chContLabel[i] );
00360 lgFlxNeg = true;
00361 }
00362 if( opac.opacity_abs[i] < 0. )
00363 {
00364 opac.lgOpacNeg = true;
00365
00366
00367 pnegopc();
00368 }
00369 }
00370 if( lgFlxNeg )
00371 {
00372 fprintf( ioQQQ, " Insanity has occurred, this is zone%4ld\n",
00373 nzone );
00374 ShowMe();
00375 puts( "[Stop in radius_increment]" );
00376 cdEXIT(EXIT_FAILURE);
00377 }
00378
00379 # endif
00380
00381
00382
00383
00384
00385
00386 if( rfield.lgOpacityFine )
00387 {
00388 dradfac = (float)radius.drad_x_fillfac;
00389
00390 for( i=0; i<rfield.nfine; ++i )
00391 {
00392 rfield.fine_opt_depth[i] +=
00393 rfield.fine_opac_zone[i]*dradfac;
00394 }
00395
00396
00397
00398
00399 if( opac.lgScatON )
00400 {
00401
00402 for( i=0; i < rfield.nflux-1; i++ )
00403 {
00404
00405
00406
00407 if( rfield.ipnt_coarse_2_fine[i] && rfield.ipnt_coarse_2_fine[i+1] )
00408 {
00409 rfield.trans_coef_zone[i] = 0.;
00410 for( j=rfield.ipnt_coarse_2_fine[i]; j<rfield.ipnt_coarse_2_fine[i+1]; ++j )
00411 {
00412
00413 rfield.trans_coef_zone[i] += 1.f / (1.f + rfield.fine_opac_zone[j]*dradfac );
00414 }
00415
00416
00417
00418 if( rfield.ipnt_coarse_2_fine[i+1]>rfield.ipnt_coarse_2_fine[i] )
00419 {
00420 rfield.trans_coef_zone[i] /= (rfield.ipnt_coarse_2_fine[i+1]-rfield.ipnt_coarse_2_fine[i]);
00421 }
00422 else
00423 {
00424
00425 rfield.trans_coef_zone[i] =
00426 1.f / (1.f + rfield.fine_opac_zone[rfield.ipnt_coarse_2_fine[i]]*dradfac );
00427 }
00428 }
00429 else
00430 {
00431
00432 rfield.trans_coef_zone[i] = 1.f;
00433 }
00434
00435
00436 rfield.trans_coef_total[i] *= rfield.trans_coef_zone[i];
00437 }
00438 }
00439 }
00440
00441
00442 escatc = 6.65e-25*dense.eden;
00443
00444
00445
00446 for( i=0; i < rfield.nflux; i++ )
00447 {
00448
00449 opac.TauAbsGeo[0][i] += (float)(opac.opacity_abs[i]*radius.drad_x_fillfac);
00450 opac.TauScatGeo[0][i] += (float)(opac.opacity_sct[i]*radius.drad_x_fillfac);
00451
00452
00453 opac.TauAbsFace[i] += (float)(opac.opacity_abs[i]*radius.drad_x_fillfac);
00454 opac.TauScatFace[i] += (float)(opac.opacity_sct[i]*radius.drad_x_fillfac);
00455
00456
00457 opac.TauTotalGeo[0][i] = opac.TauAbsGeo[0][i] + opac.TauScatGeo[0][i];
00458
00459
00460
00461
00462
00463
00464 rforce += (rfield.flux[i] + rfield.outlin[i] + rfield.outlin_noplot[i]+ rfield.ConInterOut[i])*
00465 rfield.anu[i]*(opac.opacity_abs[i] +
00466 opac.opacity_sct[i]);
00467
00468 relec += ((rfield.flux[i] + rfield.outlin[i] + rfield.outlin_noplot[i]+ rfield.ConInterOut[i])*
00469 escatc)*rfield.anu[i];
00470
00471
00472
00473
00474 opac.ExpZone[i] = sexp(opac.opacity_abs[i]*radius.drad_x_fillfac*geometry.DirectionalCosin);
00475
00476
00477 opac.ExpmTau[i] *= (float)opac.ExpZone[i];
00478
00479
00480 opac.E2TauAbsFace[i] = (float)e2(opac.TauAbsFace[i]);
00481 ASSERT( opac.E2TauAbsFace[i] <= 1. && opac.E2TauAbsFace[i] >= 0. );
00482
00483
00484 if( iteration > 1 )
00485 {
00486
00487 float tau = MAX2(SMALLFLOAT , opac.TauAbsTotal[i] - opac.TauAbsFace[i] );
00488 opac.E2TauAbsOut[i] = (float)e2( tau );
00489 }
00490
00491
00492 opfac = opac.ExpZone[i]*DilutionHere;
00493 ASSERT( opfac <= 1.0 );
00494
00495 rfield.flux[i] *= (float)opfac;
00496
00497
00498
00499 rfield.SummedCon[i] = rfield.flux[i] + rfield.SummedDif[i];
00500
00501
00502 rfield.ConInterOut[i] *= (float)opfac;
00503 rfield.outlin[i] *= (float)opfac;
00504 rfield.outlin_noplot[i] *= (float)opfac;
00505
00506
00507
00508
00509
00510
00511 rfield.ConEmitOut[i] *= (float)opfac;
00512 rfield.ConEmitOut[i] += rfield.ConEmitLocal[nzone][i]*(float)radius.dVolOutwrd*opac.tmn[i];
00513
00514
00515 rfield.OccNumbIncidCont[i] = rfield.flux[i]*rfield.convoc[i];
00516
00517
00518
00519 rfield.OccNumbDiffCont[i] = rfield.ConEmitLocal[nzone][i]*rfield.convoc[i];
00520
00521
00522 rfield.OccNumbContEmitOut[i] = rfield.ConEmitOut[i]*rfield.convoc[i];
00523
00524 }
00525
00526
00527
00528
00529
00530
00531
00532
00533
00534 if( rfield.flux[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]>SMALLFLOAT &&
00535 (rfield.flux[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]/SDIV(rfield.FluxSave[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]) ) > SMALLFLOAT &&
00536 !opac.lgScatON &&
00537 radius.depth/radius.Radius < 1e-4 )
00538 {
00539
00540
00541
00542
00543
00544 double tau_effec = -log((double)rfield.flux[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]/
00545 (double)opac.tmn[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]/
00546 (double)rfield.FluxSave[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]);
00547
00548
00549
00550
00551 double tau_true = opac.TauAbsFace[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]*geometry.DirectionalCosin;
00552
00553
00554
00555 if( fabs( tau_effec - tau_true ) / t > 0.01 &&
00556
00557
00558 fabs(tau_effec-tau_true)>MAX2(0.001,1.-opac.tmn[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]) )
00559 {
00560
00561
00562 fprintf( ioQQQ,
00563 " PROBLEM radius_increment Lyman continuum insanity zone %li, effective tau=%g, atomic tau=%g simple tau=%g\n",
00564 nzone,
00565 tau_effec ,
00566 tau_true ,
00567 6.327e-18*(dense.xIonDense[ipHYDROGEN][0]*radius.drad_x_fillfac+colden.colden[ipCOL_H0]) );
00568 TotalInsanity();
00569 }
00570 }
00571
00572
00573
00574
00575 if( opac.lgScatON )
00576 {
00577 for( i=0; i < rfield.nflux; i++ )
00578 {
00579
00580
00581
00582
00583
00584 opfac = 1./(1. + radius.drad_x_fillfac*opac.opacity_sct[i]);
00585 ASSERT( opfac <= 1.0 );
00586
00587
00588 rfield.flux[i] *= (float)opfac;
00589 rfield.ConInterOut[i] *= (float)opfac;
00590 rfield.ConEmitOut[i] *= (float)opfac;
00591 rfield.outlin[i] *= (float)opfac;
00592 rfield.outlin_noplot[i] *= (float)opfac;
00593 }
00594 }
00595
00596
00597 cmshft();
00598
00599 relec *= EN1RYD;
00600
00601
00602 t = 1./SPEEDLIGHT/dense.xMassDensity;
00603
00604 wind.AccelLine = (float)(RT_line_driving()*t)*pressure.lgLineRadPresOn;
00605
00606 wind.AccelCont = (float)(rforce*EN1RYD*t)*pressure.lgContRadPresOn;
00607
00608 wind.AccelPres = 0.;
00609
00610 wind.AccelTot = wind.AccelCont + wind.AccelLine + wind.AccelPres;
00611
00612 wind.AccelMax = (float)MAX2(wind.AccelMax,wind.AccelTot);
00613
00614
00615
00616
00617
00618 if( relec*t > SMALLFLOAT )
00619 {
00620 wind.fmul = (float)((wind.AccelLine + wind.AccelCont)/(relec*t));
00621 }
00622 else
00623 {
00624 wind.fmul = 0.;
00625 }
00626
00627
00628 pressure.pinzon = (float)(wind.AccelTot*dense.xMassDensity*radius.drad_x_fillfac*geometry.DirectionalCosin);
00629
00630
00631 pressure.PresInteg += pressure.pinzon;
00632
00633
00634 timesc.sound_speed_isothermal = sqrt(pressure.PresGasCurr/dense.xMassDensity);
00635
00636 timesc.sound_speed_adiabatic = sqrt(5.*pressure.PresGasCurr/(3.*dense.xMassDensity) );
00637 timesc.sound += radius.drad_x_fillfac / timesc.sound_speed_isothermal;
00638
00639
00640 if( neutrn.lgNeutrnHeatOn )
00641 {
00642
00643 neutrn.totneu *= (float)sexp(neutrn.CrsSecNeutron*dense.gas_phase[ipHYDROGEN]*radius.drad_x_fillfac*geometry.DirectionalCosin);
00644
00645 neutrn.totneu *= (float)DilutionHere;
00646 }
00647
00648
00649
00650
00651
00652 if( !geometry.lgSphere )
00653 {
00654 double Reflec_Diffuse_Cont;
00655
00656
00657
00658 for( i=rfield.ipPlasma-1; i < rfield.nflux; i++ )
00659 {
00660 if( opac.TauAbsGeo[0][i] < 30. )
00661 {
00662
00663
00664
00665
00666
00667
00668
00669
00670
00671
00672
00673 Reflec_Diffuse_Cont = rfield.ConEmitLocal[nzone][i]/2.*
00674 radius.drad_x_fillfac * opac.E2TauAbsFace[i]*radius.r1r0sq;
00675
00676
00677
00678 rfield.ConEmitReflec[i] += (float)(Reflec_Diffuse_Cont);
00679
00680
00681 rfield.ConRefIncid[i] += (float)(rfield.flux[i]*opac.opacity_sct[i]*
00682 radius.drad_x_fillfac*opac.E2TauAbsFace[i]*radius.r1r0sq);
00683
00684
00685
00686 rfield.reflin[i] += (float)(rfield.outlin[i]*opac.opacity_sct[i]*
00687 radius.drad_x_fillfac*opac.E2TauAbsFace[i]*radius.r1r0sq);
00688 }
00689 }
00690 }
00691
00692
00693
00694
00695
00696
00697 aver("zone",1.,1.," ");
00698 aver("doit",phycon.te,1.," Te ");
00699 aver("doit",phycon.te,dense.eden," Te(Ne) ");
00700 aver("doit",phycon.te,dense.eden*dense.xIonDense[ipHYDROGEN][1]," Te(NeNp) ");
00701 aver("doit",phycon.te,dense.eden*dense.xIonDense[ipHELIUM][1] ," Te(NeHe+)");
00702 aver("doit",phycon.te,dense.eden*dense.xIonDense[ipHELIUM][2] ,"Te(NeHe2+)");
00703 aver("doit",phycon.te,dense.eden*dense.xIonDense[ipOXYGEN][1] ," Te(NeO+) " );
00704 aver("doit",phycon.te,dense.eden*dense.xIonDense[ipOXYGEN][2] ," Te(NeO2+)");
00705
00706 aver("doit",phycon.te,hmi.H2_total," Te(H2) ");
00707 aver("doit",dense.gas_phase[ipHYDROGEN],1.," N(H) ");
00708 aver("doit",dense.eden,dense.xIonDense[ipOXYGEN][2]," Ne(O2+) ");
00709 aver("doit",dense.eden,dense.xIonDense[ipHYDROGEN][1]," Ne(Np) ");
00710
00711
00712
00713
00714
00715 nzone_minus_1 = MAX2( 0, nzone-1 );
00716
00717 struc.nzone = nzone_minus_1+1;
00718 ASSERT(nzone_minus_1>=0 && nzone_minus_1 < struc.nzlim );
00719 struc.testr[nzone_minus_1] = (float)phycon.te;
00720
00721 struc.DenParticles[nzone_minus_1] = dense.pden;
00722
00723 struc.hden[nzone_minus_1] = (float)dense.gas_phase[ipHYDROGEN];
00724
00725 struc.DenMass[nzone_minus_1] = dense.xMassDensity;
00726 struc.heatstr[nzone_minus_1] = thermal.htot;
00727 struc.coolstr[nzone_minus_1] = thermal.ctot;
00728 struc.volstr[nzone_minus_1] = (float)radius.dVeff;
00729 struc.drad[nzone_minus_1] = (float)radius.drad;
00730 struc.drad_x_fillfac[nzone_minus_1] = (float)radius.drad_x_fillfac;
00731 struc.histr[nzone_minus_1] = dense.xIonDense[ipHYDROGEN][0];
00732 struc.hiistr[nzone_minus_1] = dense.xIonDense[ipHYDROGEN][1];
00733 struc.ednstr[nzone_minus_1] = (float)dense.eden;
00734 struc.o3str[nzone_minus_1] = dense.xIonDense[ipOXYGEN][2];
00735 struc.pressure[nzone_minus_1] = (float)pressure.PresTotlCurr;
00736 struc.pres_radiation_lines_curr[nzone_minus_1] = (float)pressure.pres_radiation_lines_curr;
00737 struc.GasPressure[nzone_minus_1] = (float)pressure.PresGasCurr;
00738 struc.depth[nzone_minus_1] = (float)radius.depth;
00739
00740 struc.xLyman_depth[nzone_minus_1] = opac.TauAbsFace[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]];
00741 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
00742 {
00743 struc.gas_phase[nelem][nzone_minus_1] = dense.gas_phase[nelem];
00744 for( ion=0; ion<nelem+2; ++ion )
00745 {
00746 struc.xIonDense[nelem][ion][nzone_minus_1] = dense.xIonDense[nelem][ion];
00747 }
00748 }
00749
00750
00751 for(mol=0;mol<N_H_MOLEC;mol++)
00752 {
00753 struc.H2_molec[mol][nzone_minus_1] = hmi.Hmolec[mol];
00754 }
00755
00756
00757 for(mol=0;mol<mole.num_comole_calc;mol++)
00758 {
00759 struc.CO_molec[mol][nzone_minus_1] = COmole[mol]->hevmol;
00760 }
00761
00762 colden.dlnenp += dense.eden*(double)(dense.xIonDense[ipHYDROGEN][1])*radius.drad_x_fillfac;
00763 colden.dlnenHep += dense.eden*(double)(dense.xIonDense[ipHELIUM][1])*radius.drad_x_fillfac;
00764 colden.dlnenHepp += dense.eden*(double)(dense.xIonDense[ipHELIUM][2])*radius.drad_x_fillfac;
00765
00766
00767 colden.H0_ov_Tspin += (double)(dense.xIonDense[ipHYDROGEN][0]) / hyperfine.Tspin21cm*radius.drad_x_fillfac;
00768
00769
00770 colden.OH_ov_Tspin += (double)(findspecies("OH")->hevmol) / phycon.te*radius.drad_x_fillfac;
00771
00772
00773 hydro.TexcLya = (float)TexcLine( &EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s] );
00774
00775
00776 if( hydro.TexcLya > phycon.te )
00777 {
00778 hydro.nLyaHot += 1;
00779 if( hydro.TexcLya > hydro.TLyaMax )
00780 {
00781 hydro.TLyaMax = hydro.TexcLya;
00782 hydro.TeLyaMax = (float)phycon.te;
00783 hydro.nZTLaMax = nzone;
00784 }
00785 }
00786
00787
00788 colden.colden[ipCOL_HTOT] += (float)(dense.gas_phase[ipHYDROGEN]*radius.drad_x_fillfac);
00789 ASSERT( colden.colden[ipCOL_HTOT] >SMALLFLOAT );
00790 colden.colden[ipCOL_HMIN] += hmi.Hmolec[ipMHm]*(float)radius.drad_x_fillfac;
00791
00792
00793 colden.colden[ipCOL_H2g] += hmi.Hmolec[ipMH2g]*(float)radius.drad_x_fillfac;
00794 colden.colden[ipCOL_H2s] += hmi.Hmolec[ipMH2s]*(float)radius.drad_x_fillfac;
00795
00796 colden.coldenH2_ov_vel += hmi.H2_total*(float)radius.drad_x_fillfac / DoppVel.doppler[LIMELM+2];
00797 colden.colden[ipCOL_HeHp] += hmi.Hmolec[ipMHeHp]*(float)radius.drad_x_fillfac;
00798 colden.colden[ipCOL_H2p] += hmi.Hmolec[ipMH2p]*(float)radius.drad_x_fillfac;
00799 colden.colden[ipCOL_H3p] += hmi.Hmolec[ipMH3p]*(float)radius.drad_x_fillfac;
00800 colden.colden[ipCOL_Hp] += dense.xIonDense[ipHYDROGEN][1]*(float)radius.drad_x_fillfac;
00801 colden.colden[ipCOL_H0] += dense.xIonDense[ipHYDROGEN][0]*(float)radius.drad_x_fillfac;
00802
00803 colden.colden[ipCOL_elec] += (float)(dense.eden*radius.drad_x_fillfac);
00804
00805 colden.He123S += dense.xIonDense[ipHELIUM][1]*
00806 iso.Pop2Ion[ipHE_LIKE][ipHELIUM][ipHe2s3S]*(float)radius.drad_x_fillfac;
00807
00808
00809 h2.ortho_colden += h2.ortho_density*radius.drad_x_fillfac;
00810 h2.para_colden += h2.para_density*radius.drad_x_fillfac;
00811 ASSERT( fabs( h2.ortho_density + h2.para_density - hmi.H2_total ) / SDIV( hmi.H2_total ) < 1e-4 );
00812
00813
00814
00815
00816
00817 colden.H0_21cm_upper += (HFLines[0].PopHi*radius.drad_x_fillfac);
00818 colden.H0_21cm_lower +=(HFLines[0].PopLo*radius.drad_x_fillfac);
00819
00820
00821
00822
00823
00824 for( i=0; i < 5; i++ )
00825 {
00826
00827 colden.C2Colden[i] += colden.C2Pops[i]*(float)radius.drad_x_fillfac;
00828
00829 colden.Si2Colden[i] += colden.Si2Pops[i]*(float)radius.drad_x_fillfac;
00830 }
00831 for( i=0; i < 3; i++ )
00832 {
00833
00834 colden.C1Colden[i] += colden.C1Pops[i]*(float)radius.drad_x_fillfac;
00835
00836 colden.O1Colden[i] += colden.O1Pops[i]*(float)radius.drad_x_fillfac;
00837 }
00838 for( i=0; i < 4; i++ )
00839 {
00840
00841 colden.C3Colden[i] += colden.C3Pops[i]*(float)radius.drad_x_fillfac;
00842 }
00843
00844
00845 molcol("ADD ",ioQQQ);
00846
00847
00848 MeanInc();
00849
00850
00851
00852
00853 avWeight = 0.;
00854 for( nelem=0; nelem < LIMELM; nelem++ )
00855 {
00856 avWeight += dense.gas_phase[nelem]*dense.AtomicWeight[nelem];
00857 }
00858 avWeight /= dense.gas_phase[ipHYDROGEN];
00859
00860
00861 rfield.opac_mag_B_point = 0.;
00862 rfield.opac_mag_V_point = 0.;
00863 rfield.opac_mag_B_extended = 0.;
00864 rfield.opac_mag_V_extended = 0.;
00865 if( gv.lgGrainPhysicsOn )
00866 {
00867 long int nd;
00868 for( nd=0; nd < gv.nBin; nd++ )
00869 {
00870 gv.bin[nd]->avdust += gv.bin[nd]->tedust*(float)radius.drad_x_fillfac;
00871 gv.bin[nd]->avdft += gv.bin[nd]->DustDftVel*(float)radius.drad_x_fillfac;
00872 gv.bin[nd]->avdpot += (float)(gv.bin[nd]->dstpot*EVRYD*radius.drad_x_fillfac);
00873 gv.bin[nd]->avDGRatio += (float)(gv.bin[nd]->dustp[1]*gv.bin[nd]->dustp[2]*
00874 gv.bin[nd]->dustp[3]*gv.bin[nd]->dustp[4]*gv.bin[nd]->dstAbund/avWeight*
00875 radius.drad_x_fillfac);
00876
00877
00878
00879
00880
00881 rfield.extin_mag_B_point += (gv.bin[nd]->dstab1[rfield.ipB_filter-1] +
00882 gv.bin[nd]->pure_sc1[rfield.ipB_filter-1])*gv.bin[nd]->dstAbund*
00883 radius.drad_x_fillfac*dense.gas_phase[ipHYDROGEN] * OPTDEP2EXTIN;
00884 rfield.opac_mag_B_point += (gv.bin[nd]->dstab1[rfield.ipB_filter-1] +
00885 gv.bin[nd]->pure_sc1[rfield.ipB_filter-1])*gv.bin[nd]->dstAbund*
00886 dense.gas_phase[ipHYDROGEN] * OPTDEP2EXTIN;
00887
00888 rfield.extin_mag_V_point += (gv.bin[nd]->dstab1[rfield.ipV_filter-1] +
00889 gv.bin[nd]->pure_sc1[rfield.ipV_filter-1])*gv.bin[nd]->dstAbund*
00890 radius.drad_x_fillfac*dense.gas_phase[ipHYDROGEN] * OPTDEP2EXTIN;
00891 rfield.opac_mag_V_point += (gv.bin[nd]->dstab1[rfield.ipV_filter-1] +
00892 gv.bin[nd]->pure_sc1[rfield.ipV_filter-1])*gv.bin[nd]->dstAbund*
00893 dense.gas_phase[ipHYDROGEN] * OPTDEP2EXTIN;
00894
00895
00896
00897
00898 rfield.extin_mag_B_extended += (gv.bin[nd]->dstab1[rfield.ipB_filter-1] +
00899 gv.bin[nd]->pure_sc1[rfield.ipB_filter]*gv.bin[nd]->asym[rfield.ipB_filter-1] )*gv.bin[nd]->dstAbund*
00900 radius.drad_x_fillfac*dense.gas_phase[ipHYDROGEN] * OPTDEP2EXTIN;
00901 rfield.opac_mag_B_extended += (gv.bin[nd]->dstab1[rfield.ipB_filter-1] +
00902 gv.bin[nd]->pure_sc1[rfield.ipB_filter]*gv.bin[nd]->asym[rfield.ipB_filter-1] )*gv.bin[nd]->dstAbund*
00903 dense.gas_phase[ipHYDROGEN] * OPTDEP2EXTIN;
00904
00905 rfield.extin_mag_V_extended += (gv.bin[nd]->dstab1[rfield.ipV_filter-1] +
00906 gv.bin[nd]->pure_sc1[rfield.ipV_filter]*gv.bin[nd]->asym[rfield.ipV_filter-1] )*gv.bin[nd]->dstAbund*
00907 radius.drad_x_fillfac*dense.gas_phase[ipHYDROGEN] * OPTDEP2EXTIN;
00908 rfield.opac_mag_V_extended += (gv.bin[nd]->dstab1[rfield.ipV_filter-1] +
00909 gv.bin[nd]->pure_sc1[rfield.ipV_filter]*gv.bin[nd]->asym[rfield.ipV_filter-1] )*gv.bin[nd]->dstAbund*
00910 dense.gas_phase[ipHYDROGEN] * OPTDEP2EXTIN;
00911 }
00912 }
00913
00914
00915 colden.TotMassColl += dense.xMassDensity*(float)radius.drad_x_fillfac;
00916 colden.tmas += (float)phycon.te*dense.xMassDensity*(float)radius.drad_x_fillfac;
00917 colden.wmas += dense.wmole*dense.xMassDensity*(float)radius.drad_x_fillfac;
00918
00919
00920 rjeans = 7.79637 + (phycon.alogte - log10(dense.wmole) - log10(dense.xMassDensity*
00921 geometry.FillFac))/2.;
00922
00923
00924 ajmass = 3.*(rjeans - 0.30103) + log10(4.*PI/3.*dense.xMassDensity*
00925 geometry.FillFac);
00926
00927
00928 colden.rjnmin = MIN2(colden.rjnmin,(float)rjeans);
00929 colden.ajmmin = MIN2(colden.ajmmin,(float)ajmass);
00930
00931 if( trace.lgTrace )
00932 {
00933 fprintf( ioQQQ, " radius_increment returns\n" );
00934 }
00935
00936 DEBUG_EXIT( "radius_increment()" );
00937 return;
00938 }
00939
00940 #if !defined(NDEBUG)
00941
00942 static void pnegopc(void)
00943 {
00944 long int i;
00945 FILE *ioFile;
00946
00947 DEBUG_ENTRY( "pnegopc()" );
00948
00949 if( opac.lgNegOpacIO )
00950 {
00951
00952 if( NULL==(ioFile=fopen("negopc.txt","w")) )
00953 {
00954 fprintf( ioQQQ,"pnegopc could not open negopc.txt for writing\n");
00955 puts( "[Stop in opac0]" );
00956 cdEXIT(EXIT_FAILURE);
00957 }
00958 for( i=0; i < rfield.nflux; i++ )
00959 {
00960 fprintf( ioFile, "%10.2e %10.2e \n", rfield.anu[i],
00961 opac.opacity_abs[i] );
00962 }
00963 fclose( ioFile);
00964 }
00965
00966
00967 DEBUG_EXIT( "pnegopc()" );
00968 return;
00969 }
00970 #endif
00971
00972
00973
00974
00975
00976
00977
00978
00979
00980 static void cmshft(void)
00981 {
00982 long int i;
00983
00984 DEBUG_ENTRY( "cmshft()" );
00985
00986
00987 if( rfield.comoff == 0. )
00988 {
00989 DEBUG_EXIT( "cmshft()" );
00990 return;
00991 }
00992
00993 if( rfield.comoff != 0. )
00994 {
00995 DEBUG_EXIT( "cmshft()" );
00996 return;
00997 }
00998
00999
01000 for( i=0; i < rfield.nflux; i++ )
01001 {
01002 continue;
01003
01004
01005 }
01006
01007 DEBUG_EXIT( "cmshft()" );
01008 return;
01009 }
01010