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