00001
00002
00003
00004
00005 #include "cddefines.h"
00006 #include "physconst.h"
00007 #include "iso.h"
00008 #include "grainvar.h"
00009 #include "ca.h"
00010 #include "rfield.h"
00011 #include "oxy.h"
00012 #include "mole.h"
00013 #include "dense.h"
00014 #include "atoms.h"
00015 #include "conv.h"
00016 #include "ionbal.h"
00017 #include "trace.h"
00018 #include "hmi.h"
00019 #include "phycon.h"
00020 #include "opacity.h"
00021 #include "taulines.h"
00022
00023 void OpacityAddTotal(void)
00024 {
00025 long int i,
00026 ion,
00027 ipHi,
00028 ipISO,
00029 ipop,
00030 limit,
00031 low,
00032 nelem,
00033 n;
00034 double DepartCoefInv ,
00035 fac,
00036 sum;
00037 realnum SaveOxygen1 ,
00038 SaveCarbon1;
00039
00040 DEBUG_ENTRY( "OpacityAddTotal()" );
00041
00042
00043
00044 OpacityZero();
00045
00046
00047 for( i=0; i < rfield.nflux; i++ )
00048 {
00049
00050 opac.opacity_sct[i] += opac.OpacStack[i-1+opac.iopcom]*
00051 dense.eden;
00052 }
00053
00054
00055 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
00056 {
00057 if( dense.lgElmtOn[nelem] )
00058 {
00059 for( ion=0; ion<nelem+1; ++ion )
00060 {
00061 realnum factor = dense.xIonDense[nelem][ion];
00062
00063
00064
00065 if( nelem==ipHYDROGEN )
00066 factor += hmi.H2_total*2.f;
00067 if( factor > 0. )
00068 {
00069
00070 long loop_min = ionbal.ipCompRecoil[nelem][ion]-1;
00071 long loop_max = rfield.nflux;
00072
00073
00074 factor *= ionbal.nCompRecoilElec[nelem-ion];
00075 for( i=loop_min; i < loop_max; i++ )
00076 {
00077
00078 opac.opacity_abs[i] += opac.OpacStack[i-1+opac.iopcom]*factor;
00079 }
00080 }
00081 }
00082 }
00083 }
00084
00085
00086
00088 sum = dense.gas_phase[ipHYDROGEN] + 4.*dense.gas_phase[ipHELIUM];
00089 OpacityAdd1Subshell(opac.ioppr,opac.ippr,rfield.nflux,(realnum)sum,'s');
00090
00091
00092
00093
00094
00095
00096 static double *TotBremsAllIons;
00097 static bool lgFirstTime=true;
00098 double BremsThisEner,bfac, sumion[LIMELM+1];
00099 long int ion_lo , ion_hi;
00100
00101 if( lgFirstTime )
00102 {
00103
00104 TotBremsAllIons = (double *)MALLOC((unsigned long)rfield.nupper*sizeof(double));
00105 lgFirstTime = false;
00106 }
00107 bfac = (dense.eden/1e20)/phycon.sqrte/1e10;
00108
00109
00110 sumion[0] = 0.;
00111 for(ion=1; ion<=LIMELM; ++ion )
00112 {
00113 sumion[ion] = 0.;
00114 for( nelem=ipHELIUM; nelem < LIMELM; ++nelem )
00115 {
00116 if( dense.lgElmtOn[nelem] && ion<=nelem+1 )
00117 {
00118 sumion[ion] += dense.xIonDense[nelem][ion];
00119 }
00120 }
00121
00122 sumion[ion] *= POW2((double)ion)*bfac;
00123 }
00124
00125
00126 sumion[1] += hmi.Hmolec[ipMHeHp] * POW2((realnum)1.f)*bfac;
00127
00128
00129
00130
00131
00132 ion_lo = 1;
00133 while( sumion[ion_lo]==0 && ion_lo<LIMELM-1 )
00134 ++ion_lo;
00135 ion_hi = LIMELM;
00136 while( sumion[ion_hi]==0 && ion_hi>0 )
00137 --ion_hi;
00138
00139 for( i=0; i < rfield.nflux; i++ )
00140 {
00141
00142 nelem = ipHYDROGEN;
00143 ion = 1;
00144
00145 BremsThisEner = bfac * rfield.gff[ion][i] * (dense.xIonDense[ipHYDROGEN][1]+hmi.Hmolec[ipMH2p]+hmi.Hmolec[ipMH3p]);
00146
00147
00148 BremsThisEner += bfac * rfield.gff[ion][i] * opac.OpacStack[i-1+opac.iphmra] * StatesElemNEW[ipHYDROGEN][ipHYDROGEN-ipH_LIKE][ipH1s].Pop;
00149
00150
00151
00152
00153
00154
00155
00156
00157 for(ion=ion_lo; ion<=ion_hi; ++ion )
00158 {
00159 BremsThisEner += sumion[ion]*rfield.gff[ion][i];
00160 }
00161 ASSERT( !isnan( BremsThisEner ) );
00162 TotBremsAllIons[i] = BremsThisEner;
00163
00164
00165 if( rfield.ContBoltz[i] < 0.995 )
00166 {
00167 TotBremsAllIons[i] *= opac.OpacStack[i-1+opac.ipBrems]*
00168 (1. - rfield.ContBoltz[i]);
00169 }
00170 else
00171 {
00172 TotBremsAllIons[i] *= opac.OpacStack[i-1+opac.ipBrems]*
00173 rfield.anu[i]*TE1RYD/phycon.te;
00174 }
00175 opac.FreeFreeOpacity[i] = TotBremsAllIons[i];
00176
00177
00178
00179 opac.opacity_abs[i] += opac.FreeFreeOpacity[i];
00180 }
00181
00182
00183
00184 if( hmi.hmidep > SMALLFLOAT )
00185 {
00186 DepartCoefInv = 1./hmi.hmidep;
00187 }
00188 else
00189 {
00190
00191
00192 DepartCoefInv = 1.;
00193 }
00194
00195 limit = iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0];
00196 if(limit > rfield.nflux)
00197 limit = rfield.nflux;
00198
00199 for( i=hmi.iphmin-1; i < limit; i++ )
00200 {
00201 double factor;
00202 factor = 1. - rfield.ContBoltz[i]*DepartCoefInv;
00203 if(factor > 0)
00204 opac.opacity_abs[i] += opac.OpacStack[i-hmi.iphmin+opac.iphmop]*
00205 hmi.Hmolec[ipMHm]*factor;
00206 }
00207
00208
00209 limit = opac.ih2pnt[1];
00210 if(limit > rfield.nflux)
00211 limit = rfield.nflux;
00212 for( i=opac.ih2pnt[0]-1; i < limit; i++ )
00213 {
00214 opac.opacity_abs[i] += hmi.Hmolec[ipMH2p]*opac.OpacStack[i-opac.ih2pnt[0]+opac.ih2pof];
00215 ASSERT( !isnan( opac.opacity_abs[i] ) );
00216 }
00217
00218
00219 if( dense.xIonDense[ipHYDROGEN][1] <= 0. )
00220 {
00221 fac = dense.xIonDense[ipHYDROGEN][0];
00222 }
00223 else
00224 {
00225 fac = StatesElemNEW[ipHYDROGEN][ipHYDROGEN-ipH_LIKE][ipH1s].Pop;
00226 }
00227
00228
00229 limit = iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s];
00230 if(limit > rfield.nflux)
00231 limit = rfield.nflux;
00232 for( i=0; i < limit; i++ )
00233 {
00234 opac.opacity_sct[i] += (fac*opac.OpacStack[i-1+opac.ipRayScat]);
00235 }
00236
00237
00238 if( iso.DepartCoef[ipH_LIKE][ipHYDROGEN][ipH1s] > 1e-30 && !conv.lgSearch )
00239 {
00240 realnum factor;
00241 factor = (realnum)(rfield.ContBoltz[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1]/iso.DepartCoef[ipH_LIKE][ipHYDROGEN][ipH1s]);
00242 if(opac.stimax[0] < factor)
00243 opac.stimax[0] = factor;
00244 }
00245
00246 if( iso.DepartCoef[ipH_LIKE][ipHYDROGEN][ipH2p] > 1e-30 && !conv.lgSearch )
00247 {
00248 realnum factor;
00249 factor = (realnum)(rfield.ContBoltz[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH2p]-1]/iso.DepartCoef[ipH_LIKE][ipHYDROGEN][ipH2p]);
00250 if(opac.stimax[1] < factor)
00251 opac.stimax[1] = factor;
00252 }
00253
00254 # if 0
00255
00256 if( !conv.lgSearch )
00257 {
00258 if( Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].PopOpc < 0. )
00259 {
00260 hydro.lgHLyaMased = true;
00261 }
00262 }
00263 # endif
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273 low = opac.ipH2_photo_thresh;
00274 ipHi = rfield.nupper;
00275 ipop = opac.ipH2_photo_opac_offset;
00276
00277
00278
00279
00280 OpacityAdd1Subshell( ipop , low , ipHi , hmi.H2_total , 's' );
00281
00282
00283
00284
00285
00286
00287 SaveOxygen1 = dense.xIonDense[ipOXYGEN][0];
00288 SaveCarbon1 = dense.xIonDense[ipCARBON][0];
00289 dense.xIonDense[ipOXYGEN][0] += findspecies("CO")->hevmol + findspecies("H2Ogrn")->hevmol;
00290 dense.xIonDense[ipCARBON][0] += findspecies("CO")->hevmol;
00291
00292
00293
00294 for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
00295 {
00296
00297 if( dense.lgElmtOn[nelem] )
00298 {
00299 OpacityAdd1Element(nelem);
00300 }
00301 }
00302
00303
00304 dense.xIonDense[ipOXYGEN][0] = SaveOxygen1;
00305 dense.xIonDense[ipCARBON][0] = SaveCarbon1;
00306
00307
00308
00309
00310
00311 OpacityAdd1Subshell(opac.in1[2],opac.in1[0],opac.in1[1],
00312 dense.xIonDense[ipNITROGEN][0]*atoms.p2nit , 'v' );
00313
00314
00315
00316 OpacityAdd1Subshell(opac.ipo1exc[2],opac.ipo1exc[0],opac.ipo1exc[1],
00317 dense.xIonDense[ipOXYGEN][0]*oxy.poiexc,'v');
00318
00319
00320 OpacityAdd1Subshell(opac.ipo3exc[2],opac.ipo3exc[0],opac.ipo3exc[1],
00321 dense.xIonDense[ipOXYGEN][2]*oxy.poiii2,'v');
00322
00323 OpacityAdd1Subshell(opac.ipo3exc3[2],opac.ipo3exc3[0],opac.ipo3exc3[1],
00324 dense.xIonDense[ipOXYGEN][2]*oxy.poiii3,'v');
00325
00326
00327
00328 OpacityAdd1Subshell(opac.ipOpMgEx,opac.ipmgex,iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s],
00329 dense.xIonDense[ipMAGNESIUM][1]* atoms.popmg2,'v');
00330
00331
00332
00333 OpacityAdd1Subshell(opac.ica2op,opac.ica2ex[0],opac.ica2ex[1],
00334 ca.popca2ex,'v');
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00347 {
00348 for( nelem=ipISO; nelem < LIMELM; nelem++ )
00349 {
00350
00351 if( dense.lgElmtOn[nelem] )
00352 {
00353
00354 if( nzone < 1 )
00355 {
00356
00357 for( n=0; n < iso.numLevels_local[ipISO][nelem]; n++ )
00358 {
00359 if(iso.ipIsoLevNIonCon[ipISO][nelem][n] < rfield.nflux )
00360 {
00361
00362
00363
00364
00365 if( opac.opacity_abs[iso.ipIsoLevNIonCon[ipISO][nelem][n]-1] > 0. )
00366 {
00367
00368
00369 long int ip = iso.ipIsoLevNIonCon[ipISO][nelem][n];
00370 double t2 = csphot(
00371 ip ,
00372 ip ,
00373 iso.ipOpac[ipISO][nelem][n] );
00374
00375 iso.ConOpacRatio[ipISO][nelem][n] = 1.-
00376 (StatesElemNEW[nelem][nelem-ipISO][n].Pop*t2)/
00377 opac.opacity_abs[ip-1];
00378 }
00379 else
00380 {
00381 iso.ConOpacRatio[ipISO][nelem][n] = 0.;
00382 }
00383 }
00384 }
00385 }
00386
00387
00388 for( n=0; n < iso.numLevels_local[ipISO][nelem]; n++ )
00389 {
00390
00391 if(iso.ipIsoLevNIonCon[ipISO][nelem][n] < rfield.nflux )
00392 {
00393
00394
00395
00396
00397 if( opac.opacity_abs[iso.ipIsoLevNIonCon[ipISO][nelem][n]-1] > 0. )
00398 {
00399
00400 if( iso.DepartCoef[ipISO][nelem][n] > 1e-30 && (!conv.lgSearch ) )
00401 {
00402
00403 fac = 1./iso.DepartCoef[ipISO][nelem][n];
00404 }
00405 else if( conv.lgSearch )
00406 {
00407
00408
00409 fac = 0.;
00410 }
00411 else
00412 {
00413 fac = 1.;
00414 }
00415
00418
00419
00420
00421
00422
00423 if( opac.opacity_abs[iso.ipIsoLevNIonCon[ipISO][nelem][n]-1] > 0. )
00424 {
00425
00426 long int ip = iso.ipIsoLevNIonCon[ipISO][nelem][n];
00427
00428 double t2 = csphot(
00429 ip ,
00430 ip ,
00431 iso.ipOpac[ipISO][nelem][n] );
00432
00433 double opacity_this_species =
00434 StatesElemNEW[nelem][nelem-ipISO][n].Pop*t2*
00435 (1. - fac*rfield.ContBoltz[ip-1]);
00436
00437 double opacity_fraction = 1. - opacity_this_species / opac.opacity_abs[ip-1];
00438 if(opacity_fraction < 0)
00439 opacity_fraction = 0.;
00440
00441
00442 iso.ConOpacRatio[ipISO][nelem][n] =
00443 iso.ConOpacRatio[ipISO][nelem][n]* 0.75 + 0.25*opacity_fraction;
00444
00445 if(iso.ConOpacRatio[ipISO][nelem][n] < 0.)
00446 iso.ConOpacRatio[ipISO][nelem][n] = 0.;
00447 }
00448 else
00449 {
00450 iso.ConOpacRatio[ipISO][nelem][n] = 0.;
00451 }
00452 }
00453 else
00454 {
00455 iso.ConOpacRatio[ipISO][nelem][n] = 0.;
00456 }
00457 }
00458 else
00459 {
00460 iso.ConOpacRatio[ipISO][nelem][n] = 0.;
00461 }
00462 }
00463 }
00464 }
00465 }
00466
00467
00468 if( gv.lgDustOn() )
00469 {
00470
00471
00472 for( i=0; i < rfield.nflux; i++ )
00473 {
00474
00475 opac.opacity_sct[i] += gv.dstsc[i]*dense.gas_phase[ipHYDROGEN];
00476 opac.opacity_abs[i] += gv.dstab[i]*dense.gas_phase[ipHYDROGEN];
00477 }
00478 }
00479
00480
00481 for( i=0; i < rfield.nflux; i++ )
00482 {
00483
00484 opac.opacity_abs[i] += opac.OpacStatic[i];
00485
00486
00487 }
00488
00489
00490 for( i=0; i < rfield.nflux; i++ )
00491 {
00492 opac.albedo[i] = opac.opacity_sct[i]/
00493 (opac.opacity_sct[i] + opac.opacity_abs[i]);
00494 }
00495
00496
00497 if( conv.lgSearch )
00498 OpacityZeroOld();
00499
00500 if( trace.lgTrace )
00501 {
00502 fprintf( ioQQQ, " OpacityAddTotal returns; grd rec eff (opac) for Hn=1,4%10.2e%10.2e%10.2e%10.2e%10.2e%10.2e%10.2e%10.2e%10.2e%10.2e HeI,II:%10.2e%10.2e\n",
00503 iso.ConOpacRatio[ipH_LIKE][ipHYDROGEN][ipH1s],
00504 iso.ConOpacRatio[ipH_LIKE][ipHYDROGEN][ipH2s],
00505 iso.ConOpacRatio[ipH_LIKE][ipHYDROGEN][ipH2p],
00506 iso.ConOpacRatio[ipH_LIKE][ipHYDROGEN][ipH3s],
00507 iso.ConOpacRatio[ipH_LIKE][ipHYDROGEN][ipH3p],
00508 iso.ConOpacRatio[ipH_LIKE][ipHYDROGEN][ipH3d],
00509 iso.ConOpacRatio[ipH_LIKE][ipHYDROGEN][ipH4s],
00510 iso.ConOpacRatio[ipH_LIKE][ipHYDROGEN][ipH4p],
00511 iso.ConOpacRatio[ipH_LIKE][ipHYDROGEN][ipH4d],
00512 iso.ConOpacRatio[ipH_LIKE][ipHYDROGEN][ipH4f],
00513 iso.ConOpacRatio[ipHE_LIKE][ipHELIUM][ipHe1s1S],
00514 iso.ConOpacRatio[ipH_LIKE][ipHELIUM][ipH1s] );
00515 }
00516
00517 {
00518
00519 enum {DEBUG_LOC=false};
00520 if( DEBUG_LOC && (nzone>=378) )
00521 {
00522 if( nzone > 380 )
00523 cdEXIT( EXIT_FAILURE );
00524 for( i=0; i<rfield.nflux; ++i )
00525 {
00526 fprintf(ioQQQ,"rtotsbugggg\t%li\t%.3e\t%.3e\t%.3e\t%.3e\n",
00527 conv.nPres2Ioniz,
00528 rfield.anu[i],
00529 opac.opacity_abs[i],
00530 rfield.otscon[i],
00531 rfield.otslin[i]);
00532 }
00533 }
00534 }
00535 return;
00536 }