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.lgComptonOn )
00025 {
00026 return;
00027 }
00028
00029 if( rfield.lgComptonOn )
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 realnum tauzone = rfield.fine_opac_zone[i]*(realnum)radius.drad_x_fillfac;
00159 rfield.fine_opt_depth[i] += tauzone;
00160 }
00161 rfield.trans_coef_total_stale = true;
00162 }
00163
00164
00165
00166 double DilutionHere = POW2((radius.Radius - radius.drad*radius.dRadSign)/
00167 radius.Radius);
00168
00169 rfield.EnergyIncidCont = 0.;
00170 rfield.EnergyDiffCont = 0.;
00171
00172
00173
00174
00175
00176
00177
00178 for( long i=0; i < rfield.nflux; i++ )
00179 {
00180 double dTau_abs = opac.opacity_abs[i]*radius.drad_x_fillfac;
00181 double dTau_sct = opac.opacity_sct[i]*radius.drad_x_fillfac;
00182
00183
00184 opac.TauAbsGeo[0][i] += (realnum)(dTau_abs);
00185 opac.TauScatGeo[0][i] += (realnum)(dTau_sct);
00186
00187
00188 opac.TauAbsFace[i] += (realnum)(dTau_abs);
00189 opac.TauScatFace[i] += (realnum)(dTau_sct);
00190
00191
00192 opac.TauTotalGeo[0][i] = opac.TauAbsGeo[0][i] + opac.TauScatGeo[0][i];
00193
00194
00195
00196
00197 opac.ExpZone[i] = sexp(dTau_abs*geometry.DirectionalCosin);
00198
00199
00200 opac.ExpmTau[i] *= (realnum)opac.ExpZone[i];
00201
00202
00203 opac.E2TauAbsFace[i] = (realnum)e2(opac.TauAbsFace[i]);
00204 ASSERT( opac.E2TauAbsFace[i] <= 1. && opac.E2TauAbsFace[i] >= 0. );
00205
00206
00207 if( iteration > 1 )
00208 {
00209
00210 realnum tau = MAX2(SMALLFLOAT , opac.TauAbsTotal[i] - opac.TauAbsFace[i] );
00211 opac.E2TauAbsOut[i] = (realnum)e2( tau );
00212 ASSERT( opac.E2TauAbsOut[i]>=0. && opac.E2TauAbsOut[i]<=1. );
00213 }
00214
00215
00216 double AttenuationDilutionFactor = opac.ExpZone[i]*DilutionHere;
00217 ASSERT( AttenuationDilutionFactor <= 1.0 );
00218
00219
00220 rfield.flux_beam_const[i] *= (realnum)AttenuationDilutionFactor;
00221 rfield.flux_beam_time[i] *= (realnum)AttenuationDilutionFactor;
00222 rfield.flux_isotropic[i] *= (realnum)AttenuationDilutionFactor;
00223 rfield.flux[0][i] = rfield.flux_beam_const[i] + rfield.flux_beam_time[i] +
00224 rfield.flux_isotropic[i];
00225
00226
00227 rfield.SummedCon[i] = rfield.flux[0][i] + rfield.SummedDif[i];
00228
00229
00230 rfield.outlin[0][i] *= (realnum)AttenuationDilutionFactor;
00231 rfield.outlin_noplot[i] *= (realnum)AttenuationDilutionFactor;
00232
00233
00234 TestCode();
00235
00236 rfield.ConInterOut[i] += rfield.DiffuseEscape[i]*(realnum)radius.dVolOutwrd;
00237 rfield.ConInterOut[i] *= (realnum)AttenuationDilutionFactor;
00238
00239
00240 rfield.ConEmitOut[0][i] *= (realnum)AttenuationDilutionFactor;
00241 rfield.ConEmitOut[0][i] += rfield.ConEmitLocal[nzone][i]*(realnum)radius.dVolOutwrd*opac.tmn[i];
00242
00243
00244 rfield.OccNumbIncidCont[i] = rfield.flux[0][i]*rfield.convoc[i];
00245
00246
00247 rfield.OccNumbDiffCont[i] = rfield.ConEmitLocal[nzone][i]*rfield.convoc[i];
00248
00249
00250 rfield.OccNumbContEmitOut[i] = rfield.ConEmitOut[0][i]*rfield.convoc[i];
00251
00252
00253 rfield.EnergyIncidCont += rfield.flux[0][i]*rfield.anu[i];
00254 rfield.EnergyDiffCont += (rfield.outlin[0][i] + rfield.outlin_noplot[i] +
00255 rfield.ConInterOut[i])* rfield.anu[i];
00256 }
00257
00258
00259 rfield.EnergyIncidCont *= (realnum)EN1RYD;
00260 rfield.EnergyDiffCont *= (realnum)EN1RYD;
00261
00262
00263
00264
00265
00266
00267 if( rfield.flux[0][iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].ipIsoLevNIonCon-1]>SMALLFLOAT &&
00268 (rfield.flux[0][iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].ipIsoLevNIonCon-1]/
00269 SDIV(rfield.flux_total_incident[0][iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].ipIsoLevNIonCon-1]) ) > SMALLFLOAT &&
00270 !opac.lgScatON &&
00271 radius.depth/radius.Radius < 1e-4 )
00272 {
00273
00274
00275
00276
00277
00278 double tau_effec = -log((double)rfield.flux[0][iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].ipIsoLevNIonCon-1]/
00279 (double)opac.tmn[iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].ipIsoLevNIonCon-1]/
00280 (double)rfield.flux_total_incident[0][iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].ipIsoLevNIonCon-1]);
00281
00282
00283 double tau_true = opac.TauAbsFace[iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].ipIsoLevNIonCon-1]*geometry.DirectionalCosin;
00284
00285
00286
00287 if( fabs( tau_effec - tau_true ) / MAX2(tau_effec , tau_true) > 0.01 &&
00288
00289
00290 fabs(tau_effec-tau_true)>MAX2(0.001,1.-opac.tmn[iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].ipIsoLevNIonCon-1]) )
00291 {
00292
00293
00294 fprintf( ioQQQ,
00295 " PROBLEM radius_increment Lyman continuum insanity zone %li, effective tau=%g, atomic tau=%g simple tau=%g\n",
00296 nzone,
00297 tau_effec ,
00298 tau_true ,
00299 6.327e-18*(dense.xIonDense[ipHYDROGEN][0]*radius.drad_x_fillfac+colden.colden[ipCOL_H0]) );
00300 TotalInsanity();
00301 }
00302 }
00303
00304
00305 if( opac.lgScatON )
00306 {
00307 for( long i=0; i < rfield.nflux; i++ )
00308 {
00309
00310
00311 double AttenuationScatteringFactor = 1./(1. + radius.drad_x_fillfac*opac.opacity_sct[i]);
00312 ASSERT( AttenuationScatteringFactor <= 1.0 );
00313 rfield.flux_beam_const[i] *= (realnum)AttenuationScatteringFactor;
00314 rfield.flux_beam_time[i] *= (realnum)AttenuationScatteringFactor;
00315 rfield.flux_isotropic[i] *= (realnum)AttenuationScatteringFactor;
00316 rfield.flux[0][i] = rfield.flux_beam_const[i] + rfield.flux_beam_time[i] +
00317 rfield.flux_isotropic[i];
00318
00319 rfield.ConInterOut[i] *= (realnum)AttenuationScatteringFactor;
00320 rfield.ConEmitOut[0][i] *= (realnum)AttenuationScatteringFactor;
00321 rfield.outlin[0][i] *= (realnum)AttenuationScatteringFactor;
00322 rfield.outlin_noplot[i] *= (realnum)AttenuationScatteringFactor;
00323 }
00324 }
00325
00326
00327
00328 realnum Dilution = (realnum)POW2( radius.rinner / (radius.Radius-radius.drad/2.) );
00329
00330
00331
00332
00333
00334
00335
00336 rfield.ConEmitLocal[nzone][rfield.nflux] = 1.e-10f * Dilution;
00337 rfield.DiffuseEscape[rfield.nflux] = 1.e-10f * Dilution;
00338
00339 rfield.ConInterOut[rfield.nflux] +=
00340 rfield.DiffuseEscape[rfield.nflux]*(realnum)radius.dVolOutwrd;
00341
00342
00343 ASSERT( opac.opacity_abs[rfield.nflux] == 0. );
00344
00345
00346 cmshft();
00347
00348
00349 if( hextra.lgNeutrnHeatOn )
00350 {
00351
00352 hextra.totneu *= (realnum)sexp(hextra.CrsSecNeutron*
00353 dense.gas_phase[ipHYDROGEN]*radius.drad_x_fillfac*geometry.DirectionalCosin);
00354
00355 hextra.totneu *= (realnum)DilutionHere;
00356 }
00357
00358
00359
00360
00361
00362 if( !geometry.lgSphere )
00363 {
00364 double Reflec_Diffuse_Cont;
00365
00366
00367 for( long i=rfield.ipPlasma-1; i < rfield.nflux; i++ )
00368 {
00369 if( opac.TauAbsGeo[0][i] < 30. )
00370 {
00371
00372
00373 Reflec_Diffuse_Cont = rfield.ConEmitLocal[nzone][i]/2.*
00374 radius.drad_x_fillfac * opac.E2TauAbsFace[i]*radius.r1r0sq;
00375
00376
00377 rfield.ConEmitReflec[0][i] += (realnum)(Reflec_Diffuse_Cont);
00378
00379
00380 rfield.ConRefIncid[0][i] += (realnum)(rfield.flux[0][i]*opac.opacity_sct[i]*
00381 radius.drad_x_fillfac*opac.E2TauAbsFace[i]*radius.r1r0sq);
00382
00383
00384 rfield.reflin[0][i] += (realnum)(rfield.outlin[0][i]*opac.opacity_sct[i]*
00385 radius.drad_x_fillfac*opac.E2TauAbsFace[i]*radius.r1r0sq);
00386 }
00387 }
00388 }
00389 return;
00390 }