00001
00002
00003
00004
00005 #include "cddefines.h"
00006 #include "rfield.h"
00007 #include "opacity.h"
00008 #include "dense.h"
00009 #include "colden.h"
00010 #include "opacity.h"
00011 #include "geometry.h"
00012 #include "trace.h"
00013 #include "radius.h"
00014 #include "iso.h"
00015 #include "hextra.h"
00016
00017
00018
00019 STATIC void cmshft(void)
00020 {
00021 DEBUG_ENTRY( "cmshft()" );
00022
00023
00024 if( rfield.comoff == 0. )
00025 {
00026 return;
00027 }
00028
00029 if( rfield.comoff != 0. )
00030 {
00031 return;
00032 }
00033
00034
00035 for( long i=0; i < rfield.nflux; i++ )
00036 {
00037 continue;
00038 }
00039 return;
00040 }
00041
00042
00043 #if !defined(NDEBUG)
00044
00045 STATIC void pnegopc(void)
00046 {
00047 FILE *ioFile;
00048
00049 DEBUG_ENTRY( "pnegopc()" );
00050
00051 if( opac.lgNegOpacIO )
00052 {
00053
00054 ioFile = open_data( "negopc.txt", "w", AS_LOCAL_ONLY );
00055 for( long i=0; i < rfield.nflux; i++ )
00056 {
00057 fprintf( ioFile, "%10.2e %10.2e \n", rfield.anu[i],
00058 opac.opacity_abs[i] );
00059 }
00060 fclose( ioFile);
00061 }
00062 return;
00063 }
00064 #endif
00065
00066
00067 void RT_continuum(void)
00068 {
00069
00070 DEBUG_ENTRY( "RT_continuum()" );
00071
00072 if( trace.lgTrace && trace.lgConBug )
00073 {
00074 fprintf( ioQQQ, " Energy, flux, OTS:\n" );
00075 for( long i=0; i < rfield.nflux; i++ )
00076 {
00077 fprintf( ioQQQ, "%6ld%10.2e%10.2e%10.2e", i, rfield.anu[i],
00078 rfield.flux[0][i] + rfield.outlin[0][i] + rfield.ConInterOut[i],
00079 rfield.otscon[i] + rfield.otslin[i] + rfield.outlin_noplot[i]);
00080 }
00081 fprintf( ioQQQ, "\n" );
00082 }
00083
00084
00085 # if !defined(NDEBUG)
00086 bool lgFlxNeg = false;
00087 for( long i=0; i < rfield.nflux; i++ )
00088 {
00089 if( rfield.flux[0][i] < 0. )
00090 {
00091 fprintf( ioQQQ, " radius_increment finds negative intensity in flux.\n" );
00092 fprintf( ioQQQ, " Intensity, frequency, pointer=%11.3e%11.3e%6ld\n",
00093 rfield.flux[0][i], rfield.anu[i], i );
00094 lgFlxNeg = true;
00095 }
00096 if( rfield.otscon[i] < 0. )
00097 {
00098 fprintf( ioQQQ, " radius_increment finds negative intensity in otscon.\n" );
00099 fprintf( ioQQQ, " Intensity, frequency, pointer=%11.3e%11.3e%6ld\n",
00100 rfield.otscon[i], rfield.anu[i], i );
00101 lgFlxNeg = true;
00102 }
00103 if( opac.tmn[i] < 0. )
00104 {
00105 fprintf( ioQQQ, " radius_increment finds negative tmn.\n" );
00106 fprintf( ioQQQ, " value, frequency, pointer=%11.3e%11.3e%6ld %4.4s\n",
00107 opac.tmn[i], rfield.anu[i], i, rfield.chLineLabel[i] );
00108 lgFlxNeg = true;
00109 }
00110 if( rfield.otslin[i] < 0. )
00111 {
00112 fprintf( ioQQQ, " radius_increment finds negative intensity in otslin.\n" );
00113 fprintf( ioQQQ, " Intensity, frequency, pointer=%11.3e%11.3e%6ld %4.4s\n",
00114 rfield.otslin[i], rfield.anu[i], i, rfield.chLineLabel[i] );
00115 lgFlxNeg = true;
00116 }
00117 if( rfield.outlin[0][i] < 0. )
00118 {
00119 fprintf( ioQQQ, " radius_increment finds negative intensity in outlin.\n" );
00120 fprintf( ioQQQ, " Intensity, frequency, pointer=%11.3e%11.3e%6ld %4.4s\n",
00121 rfield.outlin[0][i], rfield.anu[i], i, rfield.chLineLabel[i] );
00122 lgFlxNeg = true;
00123 }
00124 if( rfield.ConInterOut[i] < 0. )
00125 {
00126 fprintf( ioQQQ, " radius_increment finds negative intensity in ConInterOut.\n" );
00127 fprintf( ioQQQ, " Intensity, frequency, pointer=%11.3e%11.3e%6ld %4.4s\n",
00128 rfield.ConInterOut[i], rfield.anu[i], i, rfield.chContLabel[i] );
00129 lgFlxNeg = true;
00130 }
00131 if( opac.opacity_abs[i] < 0. )
00132 {
00133 opac.lgOpacNeg = true;
00134
00135
00136 pnegopc();
00137 }
00138 }
00139 if( lgFlxNeg )
00140 {
00141 fprintf( ioQQQ, " Insanity has occurred, this is zone%4ld\n",
00142 nzone );
00143 ShowMe();
00144 cdEXIT(EXIT_FAILURE);
00145 }
00146
00147 # endif
00148
00149
00150
00151
00152
00153 if( rfield.lgOpacityFine )
00154 {
00155
00156 for( long i=0; i<rfield.nfine; ++i )
00157 {
00158 rfield.fine_opt_depth[i] +=
00159 rfield.fine_opac_zone[i]*(realnum)radius.drad_x_fillfac;
00160 }
00161
00162
00163
00164 if( opac.lgScatON )
00165 {
00166
00167 for( long i=0; i < rfield.nflux-1; i++ )
00168 {
00169
00170
00171
00172 if( rfield.ipnt_coarse_2_fine[i] && rfield.ipnt_coarse_2_fine[i+1] )
00173 {
00174 rfield.trans_coef_zone[i] = 0.;
00175 for( long j=rfield.ipnt_coarse_2_fine[i]; j<rfield.ipnt_coarse_2_fine[i+1]; ++j )
00176 rfield.trans_coef_zone[i] += 1.f /
00177 (1.f + rfield.fine_opac_zone[j]*(realnum)radius.drad_x_fillfac );
00178
00179
00180
00181
00182 if( rfield.ipnt_coarse_2_fine[i+1]>rfield.ipnt_coarse_2_fine[i] )
00183 {
00184 rfield.trans_coef_zone[i] /= (rfield.ipnt_coarse_2_fine[i+1]-rfield.ipnt_coarse_2_fine[i]);
00185 }
00186 else
00187 {
00188
00189
00190 rfield.trans_coef_zone[i] = 1.f /
00191 (1.f + rfield.fine_opac_zone[rfield.ipnt_coarse_2_fine[i]]*
00192 (realnum)radius.drad_x_fillfac );
00193 }
00194 }
00195 else
00196 {
00197
00198
00199 rfield.trans_coef_zone[i] = 1.f;
00200 }
00201
00202
00203
00204 rfield.trans_coef_total[i] *= rfield.trans_coef_zone[i];
00205 }
00206 }
00207 }
00208
00209
00210
00211 double DilutionHere = POW2((radius.Radius - radius.drad*radius.dRadSign)/
00212 radius.Radius);
00213
00214 rfield.EnergyIncidCont = 0.;
00215 rfield.EnergyDiffCont = 0.;
00216
00217
00218
00219
00220
00221
00222
00223 for( long i=0; i < rfield.nflux; i++ )
00224 {
00225 double dTau_abs = opac.opacity_abs[i]*radius.drad_x_fillfac;
00226 double dTau_sct = opac.opacity_sct[i]*radius.drad_x_fillfac;
00227
00228
00229 opac.TauAbsGeo[0][i] += (realnum)(dTau_abs);
00230 opac.TauScatGeo[0][i] += (realnum)(dTau_sct);
00231
00232
00233 opac.TauAbsFace[i] += (realnum)(dTau_abs);
00234 opac.TauScatFace[i] += (realnum)(dTau_sct);
00235
00236
00237 opac.TauTotalGeo[0][i] = opac.TauAbsGeo[0][i] + opac.TauScatGeo[0][i];
00238
00239
00240
00241
00242 opac.ExpZone[i] = sexp(dTau_abs*geometry.DirectionalCosin);
00243
00244
00245 opac.ExpmTau[i] *= (realnum)opac.ExpZone[i];
00246
00247
00248 opac.E2TauAbsFace[i] = (realnum)e2(opac.TauAbsFace[i]);
00249 ASSERT( opac.E2TauAbsFace[i] <= 1. && opac.E2TauAbsFace[i] >= 0. );
00250
00251
00252 if( iteration > 1 )
00253 {
00254
00255 realnum tau = MAX2(SMALLFLOAT , opac.TauAbsTotal[i] - opac.TauAbsFace[i] );
00256 opac.E2TauAbsOut[i] = (realnum)e2( tau );
00257 ASSERT( opac.E2TauAbsOut[i]>=0. && opac.E2TauAbsOut[i]<=1. );
00258 }
00259
00260
00261 double AttenuationDilutionFactor = opac.ExpZone[i]*DilutionHere;
00262 ASSERT( AttenuationDilutionFactor <= 1.0 );
00263
00264
00265 rfield.flux_beam_const[i] *= (realnum)AttenuationDilutionFactor;
00266 rfield.flux_beam_time[i] *= (realnum)AttenuationDilutionFactor;
00267 rfield.flux_isotropic[i] *= (realnum)AttenuationDilutionFactor;
00268 rfield.flux[0][i] = rfield.flux_beam_const[i] + rfield.flux_beam_time[i] +
00269 rfield.flux_isotropic[i];
00270
00271
00272 rfield.SummedCon[i] = rfield.flux[0][i] + rfield.SummedDif[i];
00273
00274
00275 rfield.outlin[0][i] *= (realnum)AttenuationDilutionFactor;
00276 rfield.outlin_noplot[i] *= (realnum)AttenuationDilutionFactor;
00277
00278
00279 TestCode();
00280
00281 rfield.ConInterOut[i] += rfield.DiffuseEscape[i]*(realnum)radius.dVolOutwrd;
00282 rfield.ConInterOut[i] *= (realnum)AttenuationDilutionFactor;
00283
00284
00285 rfield.ConEmitOut[0][i] *= (realnum)AttenuationDilutionFactor;
00286 rfield.ConEmitOut[0][i] += rfield.ConEmitLocal[nzone][i]*(realnum)radius.dVolOutwrd*opac.tmn[i];
00287
00288
00289 rfield.OccNumbIncidCont[i] = rfield.flux[0][i]*rfield.convoc[i];
00290
00291
00292 rfield.OccNumbDiffCont[i] = rfield.ConEmitLocal[nzone][i]*rfield.convoc[i];
00293
00294
00295 rfield.OccNumbContEmitOut[i] = rfield.ConEmitOut[0][i]*rfield.convoc[i];
00296
00297
00298 rfield.EnergyIncidCont += rfield.flux[0][i]*rfield.anu[i];
00299 rfield.EnergyDiffCont += (rfield.outlin[0][i] + rfield.outlin_noplot[i] +
00300 rfield.ConInterOut[i])* rfield.anu[i];
00301 }
00302
00303
00304 rfield.EnergyIncidCont *= (realnum)EN1RYD;
00305 rfield.EnergyDiffCont *= (realnum)EN1RYD;
00306
00307
00308
00309
00310
00311
00312 if( rfield.flux[0][iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]>SMALLFLOAT &&
00313 (rfield.flux[0][iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]/
00314 SDIV(rfield.flux_total_incident[0][iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]) ) > SMALLFLOAT &&
00315 !opac.lgScatON &&
00316 radius.depth/radius.Radius < 1e-4 )
00317 {
00318
00319
00320
00321
00322
00323 double tau_effec = -log((double)rfield.flux[0][iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]/
00324 (double)opac.tmn[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]/
00325 (double)rfield.flux_total_incident[0][iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]);
00326
00327
00328 double tau_true = opac.TauAbsFace[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]*geometry.DirectionalCosin;
00329
00330
00331
00332 if( fabs( tau_effec - tau_true ) / MAX2(tau_effec , tau_true) > 0.01 &&
00333
00334
00335 fabs(tau_effec-tau_true)>MAX2(0.001,1.-opac.tmn[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]-1]) )
00336 {
00337
00338
00339 fprintf( ioQQQ,
00340 " PROBLEM radius_increment Lyman continuum insanity zone %li, effective tau=%g, atomic tau=%g simple tau=%g\n",
00341 nzone,
00342 tau_effec ,
00343 tau_true ,
00344 6.327e-18*(dense.xIonDense[ipHYDROGEN][0]*radius.drad_x_fillfac+colden.colden[ipCOL_H0]) );
00345 TotalInsanity();
00346 }
00347 }
00348
00349
00350 if( opac.lgScatON )
00351 {
00352 for( long i=0; i < rfield.nflux; i++ )
00353 {
00354
00355
00356 double AttenuationScatteringFactor = 1./(1. + radius.drad_x_fillfac*opac.opacity_sct[i]);
00357 ASSERT( AttenuationScatteringFactor <= 1.0 );
00358 rfield.flux_beam_const[i] *= (realnum)AttenuationScatteringFactor;
00359 rfield.flux_beam_time[i] *= (realnum)AttenuationScatteringFactor;
00360 rfield.flux_isotropic[i] *= (realnum)AttenuationScatteringFactor;
00361 rfield.flux[0][i] = rfield.flux_beam_const[i] + rfield.flux_beam_time[i] +
00362 rfield.flux_isotropic[i];
00363
00364 rfield.ConInterOut[i] *= (realnum)AttenuationScatteringFactor;
00365 rfield.ConEmitOut[0][i] *= (realnum)AttenuationScatteringFactor;
00366 rfield.outlin[0][i] *= (realnum)AttenuationScatteringFactor;
00367 rfield.outlin_noplot[i] *= (realnum)AttenuationScatteringFactor;
00368 }
00369 }
00370
00371
00372
00373 realnum Dilution = (realnum)POW2( radius.rinner / (radius.Radius-radius.drad/2.) );
00374
00375
00376
00377
00378
00379
00380
00381 rfield.ConEmitLocal[nzone][rfield.nflux] = 1.e-10f * Dilution;
00382 rfield.DiffuseEscape[rfield.nflux] = 1.e-10f * Dilution;
00383
00384 rfield.ConInterOut[rfield.nflux] +=
00385 rfield.DiffuseEscape[rfield.nflux]*(realnum)radius.dVolOutwrd;
00386
00387
00388 ASSERT( opac.opacity_abs[rfield.nflux] == 0. );
00389
00390
00391 cmshft();
00392
00393
00394 if( hextra.lgNeutrnHeatOn )
00395 {
00396
00397 hextra.totneu *= (realnum)sexp(hextra.CrsSecNeutron*
00398 dense.gas_phase[ipHYDROGEN]*radius.drad_x_fillfac*geometry.DirectionalCosin);
00399
00400 hextra.totneu *= (realnum)DilutionHere;
00401 }
00402
00403
00404
00405
00406
00407 if( !geometry.lgSphere )
00408 {
00409 double Reflec_Diffuse_Cont;
00410
00411
00412 for( long i=rfield.ipPlasma-1; i < rfield.nflux; i++ )
00413 {
00414 if( opac.TauAbsGeo[0][i] < 30. )
00415 {
00416
00417
00418 Reflec_Diffuse_Cont = rfield.ConEmitLocal[nzone][i]/2.*
00419 radius.drad_x_fillfac * opac.E2TauAbsFace[i]*radius.r1r0sq;
00420
00421
00422 rfield.ConEmitReflec[0][i] += (realnum)(Reflec_Diffuse_Cont);
00423
00424
00425 rfield.ConRefIncid[0][i] += (realnum)(rfield.flux[0][i]*opac.opacity_sct[i]*
00426 radius.drad_x_fillfac*opac.E2TauAbsFace[i]*radius.r1r0sq);
00427
00428
00429 rfield.reflin[0][i] += (realnum)(rfield.outlin[0][i]*opac.opacity_sct[i]*
00430 radius.drad_x_fillfac*opac.E2TauAbsFace[i]*radius.r1r0sq);
00431 }
00432 }
00433 }
00434 return;
00435 }