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 #define RJRWLOOP 1
00022 #define RJRWMACRO 1
00023
00024 void OpacityAddTotal(void)
00025 {
00026 long int i,
00027 ion,
00028 ipHi,
00029 ipISO,
00030 ipop,
00031 limit,
00032 low,
00033 nelem,
00034 n;
00035 double DepartCoefInv ,
00036 fac,
00037 sum;
00038 float SaveOxygen1 ,
00039 SaveCarbon1;
00040
00041 DEBUG_ENTRY( "OpacityAddTotal()" );
00042
00043
00044
00045 OpacityZero();
00046
00047
00048
00049
00050 for( i=0; i < rfield.nflux; i++ )
00051 {
00052
00053 opac.opacity_sct[i] += opac.OpacStack[i-1+opac.iopcom]*
00054 dense.eden;
00055 }
00056
00057
00058
00059
00060 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
00061 {
00062 if( dense.lgElmtOn[nelem] )
00063 {
00064 for( ion=0; ion<nelem+1; ++ion )
00065 {
00066 float factor = dense.xIonDense[nelem][ion];
00067
00068
00069
00070 if( nelem==ipHYDROGEN )
00071 factor += hmi.H2_total*2.f;
00072 if( factor > 0. )
00073 {
00074
00075 long loop_min = ionbal.ipCompRecoil[nelem][ion]-1;
00076 long loop_max = rfield.nflux;
00077
00078
00079 factor *= ionbal.nCompRecoilElec[nelem-ion];
00080 for( i=loop_min; i < loop_max; i++ )
00081 {
00082
00083 opac.opacity_abs[i] += opac.OpacStack[i-1+opac.iopcom]*factor;
00084 }
00085 }
00086 }
00087 }
00088 }
00089
00090
00091
00093 sum = dense.gas_phase[ipHYDROGEN] + 4.*dense.gas_phase[ipHELIUM];
00094 OpacityAdd1Subshell(opac.ioppr,opac.ippr,rfield.nflux,(float)sum,'s');
00095
00096
00097
00098
00099
00100
00101 if(RJRWLOOP)
00102 {
00103 static double *TotBremsAllIons;
00104 static bool lgFirstTime=true;
00105 double BremsThisEner,bfac, bhfac,sumion[LIMELM+1];
00106 long int ion_lo , ion_hi;
00107
00108 if( lgFirstTime )
00109 {
00110
00111 TotBremsAllIons = (double *)MALLOC((unsigned long)rfield.nupper*sizeof(double));
00112 lgFirstTime = false;
00113 }
00114 bfac = (dense.eden/1e20)/phycon.sqrte/1e10;
00115
00116
00117 sumion[0] = 0.;
00118 for(ion=1; ion<=LIMELM; ++ion )
00119 {
00120 sumion[ion] = 0.;
00121 for( nelem=ipHELIUM; nelem < LIMELM; ++nelem )
00122 {
00123 if( dense.lgElmtOn[nelem] && ion<=nelem+1 )
00124 {
00125 sumion[ion] += dense.xIonDense[nelem][ion];
00126 }
00127 }
00128
00129 sumion[ion] *= POW2((double)ion)*bfac;
00130 }
00131
00132
00133
00134
00135 ion_lo = 1;
00136 while( sumion[ion_lo]==0 && ion_lo<LIMELM-1 )
00137 ++ion_lo;
00138 ion_hi = LIMELM;
00139 while( sumion[ion_hi]==0 && ion_hi>0 )
00140 --ion_hi;
00141
00142 bhfac = bfac*(dense.xIonDense[ipHYDROGEN][1]+hmi.Hmolec[ipMH2p]+hmi.Hmolec[ipMH3p]);
00143 for( i=0; i < rfield.nflux; i++ )
00144 {
00145
00146 nelem = ipHYDROGEN;
00147 ion = 1;
00148
00149
00150
00151 BremsThisEner = bhfac * rfield.gff[ion][i]*
00152 (1.+opac.OpacStack[i-1+opac.iphmra]*iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH1s]);
00153
00154
00155
00156
00157
00158
00159
00160
00161 for(ion=ion_lo; ion<=ion_hi; ++ion )
00162 {
00163 BremsThisEner += sumion[ion]*rfield.gff[ion][i];
00164 }
00165 TotBremsAllIons[i] = BremsThisEner;
00166
00167
00168 if( rfield.ContBoltz[i] < 0.995 )
00169 {
00170 TotBremsAllIons[i] *= opac.OpacStack[i-1+opac.ipBrems]*
00171 (1. - rfield.ContBoltz[i]);
00172 }
00173 else
00174 {
00175 TotBremsAllIons[i] *= opac.OpacStack[i-1+opac.ipBrems]*
00176 rfield.anu[i]*TE1RYD/phycon.te;
00177 }
00178 opac.FreeFreeOpacity[i] = TotBremsAllIons[i];
00179
00180
00181
00182 opac.opacity_abs[i] += opac.FreeFreeOpacity[i];
00183 }
00184 }
00185 else
00186 {
00187 for( i=0; i < rfield.nflux; i++ )
00188 {
00189 double TotBremsAllIons = 0., BremsThisIon;
00190
00191
00192 nelem = ipHYDROGEN;
00193 ion = 1;
00194 BremsThisIon = (dense.xIonDense[ipHYDROGEN][1]+hmi.Hmolec[ipMH2p]+hmi.Hmolec[ipMH3p]) *
00195 POW2( (double)ion )/1e10*rfield.gff[ion][i] *(dense.eden/1e20)/phycon.sqrte;
00196
00197
00198
00199 BremsThisIon *= (1.+opac.OpacStack[i-1+opac.iphmra]*iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH1s]);
00200 TotBremsAllIons = BremsThisIon;
00201
00202
00203 ASSERT( (dense.xIonDense[nelem][ion]==0.) || (TotBremsAllIons > 0.) );
00204
00205
00206
00207 for( nelem=ipHELIUM; nelem < LIMELM; nelem++ )
00208 {
00209
00210
00211
00212 for( ion=MAX2(1,dense.IonLow[nelem]); ion<=dense.IonHigh[nelem]; ++ion )
00213 {
00214
00215 BremsThisIon = POW2( (double)ion )*dense.xIonDense[nelem][ion]/1e10*rfield.gff[ion][i] *
00216 (dense.eden/1e20)/phycon.sqrte;
00217 TotBremsAllIons += BremsThisIon;
00218 }
00219 }
00220
00221 if( rfield.ContBoltz[i] < 0.995 )
00222 {
00223 TotBremsAllIons *= opac.OpacStack[i-1+opac.ipBrems]*
00224 (1. - rfield.ContBoltz[i]);
00225 }
00226 else
00227 {
00228 TotBremsAllIons *= opac.OpacStack[i-1+opac.ipBrems]*
00229 rfield.anu[i]*TE1RYD/phycon.te;
00230 }
00231 opac.FreeFreeOpacity[i] = TotBremsAllIons;
00232
00233
00234 ASSERT( (opac.FreeFreeOpacity[i] > 0.) || (dense.xIonDense[ipHYDROGEN][1] == 0.) );
00235 opac.opacity_abs[i] += opac.FreeFreeOpacity[i];
00236 }
00237 }
00238
00239
00240 if( hmi.hmidep > SMALLFLOAT )
00241 {
00242 DepartCoefInv = 1./hmi.hmidep;
00243 }
00244 else
00245 {
00246
00247
00248 DepartCoefInv = 1.;
00249 }
00250
00251 limit = iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0];
00252 if(limit > rfield.nflux)
00253 limit = rfield.nflux;
00254
00255 for( i=hmi.iphmin-1; i < limit; i++ )
00256 {
00257 #if RJRWMACRO
00258 double factor;
00259 factor = 1. - rfield.ContBoltz[i]*DepartCoefInv;
00260 if(factor > 0)
00261 #if 1
00262 opac.opacity_abs[i] += (opac.OpacStack[i-hmi.iphmin+opac.iphmop]*
00263 hmi.Hmolec[ipMHm]*factor);
00264 #else
00265 opac.opacity_abs[i] += (opac.OpacStack[i-hmi.iphmin+opac.iphmop]*
00266 hmi.hminus*factor);
00267 #endif
00268 #else
00269 opac.opacity_abs[i] += (opac.OpacStack[i-hmi.iphmin+opac.iphmop]*
00270 hmi.hminus*MAX2(0.,(1. - rfield.ContBoltz[i]*DepartCoefInv)));
00271 #endif
00272 }
00273
00274
00275 #if RJRWMACRO
00276 limit = opac.ih2pnt[1];
00277 if(limit > rfield.nflux)
00278 limit = rfield.nflux;
00279 #else
00280 limit = MIN2(rfield.nflux,opac.ih2pnt[1]);
00281 #endif
00282 for( i=opac.ih2pnt[0]-1; i < limit; i++ )
00283 {
00284 opac.opacity_abs[i] += hmi.Hmolec[ipMH2p]*opac.OpacStack[i-opac.ih2pnt[0]+
00285 opac.ih2pof];
00286 }
00287
00288
00289 if( dense.xIonDense[ipHYDROGEN][1] <= 0. )
00290 {
00291 fac = dense.xIonDense[ipHYDROGEN][0];
00292 }
00293 else
00294 {
00295 fac = dense.xIonDense[ipHYDROGEN][1]*iso.Pop2Ion[ipH_LIKE][ipHYDROGEN][ipH1s];
00296 }
00297
00298
00299 #if RJRWMACRO
00300 limit = iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s];
00301 if(limit > rfield.nflux)
00302 limit = rfield.nflux;
00303 #else
00304 limit = MIN2(rfield.nflux, iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]);
00305 #endif
00306 for( i=0; i < limit; i++ )
00307 {
00308 opac.opacity_sct[i] += (fac*opac.OpacStack[i-1+opac.ipRayScat]);
00309 }
00310
00311
00312 if( iso.DepartCoef[ipH_LIKE][ipHYDROGEN][ipH1s] > 1e-30 && !conv.lgSearch )
00313 {
00314 #if RJRWMACRO
00315 float factor;
00316 factor = (float)(rfield.ContBoltz[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1]/iso.DepartCoef[ipH_LIKE][ipHYDROGEN][ipH1s]);
00317 if(opac.stimax[0] < factor)
00318 opac.stimax[0] = factor;
00319 #else
00320 opac.stimax[0] =
00321 MAX2(opac.stimax[0],
00322 (float)(rfield.ContBoltz[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1]/iso.DepartCoef[ipH_LIKE][ipHYDROGEN][ipH1s]));
00323 #endif
00324 }
00325
00326 if( iso.DepartCoef[ipH_LIKE][ipHYDROGEN][ipH2p] > 1e-30 && !conv.lgSearch )
00327 {
00328 #if RJRWMACRO
00329 float factor;
00330 factor = (float)(rfield.ContBoltz[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH2p]-1]/iso.DepartCoef[ipH_LIKE][ipHYDROGEN][ipH2p]);
00331 if(opac.stimax[1] < factor)
00332 opac.stimax[1] = factor;
00333 #else
00334 opac.stimax[1] =
00335 MAX2(opac.stimax[1],
00336 (float)(rfield.ContBoltz[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH2p]-1]/iso.DepartCoef[ipH_LIKE][ipHYDROGEN][ipH2p]));
00337 #endif
00338 }
00339
00340 # if 0
00341
00342 if( !conv.lgSearch )
00343 {
00344 if( EmisLines[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].PopOpc < 0. )
00345 {
00346 hydro.lgHLyaMased = true;
00347 }
00348 }
00349 # endif
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359 low = opac.ipH2_photo_thresh;
00360 ipHi = rfield.nupper;
00361 ipop = opac.ipH2_photo_opac_offset;
00362
00363
00364
00365
00366 OpacityAdd1Subshell( ipop , low , ipHi , hmi.H2_total , 's' );
00367
00368
00369
00370
00371
00372
00373 SaveOxygen1 = dense.xIonDense[ipOXYGEN][0];
00374 SaveCarbon1 = dense.xIonDense[ipCARBON][0];
00375 dense.xIonDense[ipOXYGEN][0] += findspecies("CO")->hevmol + findspecies("H2Ogrn")->hevmol;
00376 dense.xIonDense[ipCARBON][0] += findspecies("CO")->hevmol;
00377
00378
00379
00380 for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
00381 {
00382
00383 if( dense.lgElmtOn[nelem] )
00384 {
00385 OpacityAdd1Element(nelem);
00386 }
00387 }
00388
00389
00390 dense.xIonDense[ipOXYGEN][0] = SaveOxygen1;
00391 dense.xIonDense[ipCARBON][0] = SaveCarbon1;
00392
00393
00394
00395
00396
00397
00398 OpacityAdd1Subshell(opac.in1[2],opac.in1[0],opac.in1[1],
00399 dense.xIonDense[ipNITROGEN][0]*atoms.p2nit , 'v' );
00400
00401
00402
00403 OpacityAdd1Subshell(opac.ipo1exc[2],opac.ipo1exc[0],opac.ipo1exc[1],
00404 dense.xIonDense[ipOXYGEN][0]*oxy.poiexc,'v');
00405
00406
00407 OpacityAdd1Subshell(opac.ipo3exc[2],opac.ipo3exc[0],opac.ipo3exc[1],
00408 dense.xIonDense[ipOXYGEN][2]*oxy.poiii2,'v');
00409
00410 OpacityAdd1Subshell(opac.ipo3exc3[2],opac.ipo3exc3[0],opac.ipo3exc3[1],
00411 dense.xIonDense[ipOXYGEN][2]*oxy.poiii3,'v');
00412
00413
00414
00415 OpacityAdd1Subshell(opac.ipOpMgEx,opac.ipmgex,iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s],
00416 dense.xIonDense[ipMAGNESIUM][1]* atoms.popmg2,'v');
00417
00418
00419
00420 OpacityAdd1Subshell(opac.ica2op,opac.ica2ex[0],opac.ica2ex[1],
00421 ca.popca2ex,'v');
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00434 {
00435 for( nelem=ipISO; nelem < LIMELM; nelem++ )
00436 {
00437
00438 if( dense.lgElmtOn[nelem] )
00439 {
00440
00441 if( nzone < 1 )
00442 {
00443
00444 for( n=0; n < iso.numLevels_local[ipISO][nelem]; n++ )
00445 {
00446 if(iso.ipIsoLevNIonCon[ipISO][nelem][n] < rfield.nflux )
00447 {
00448
00449
00450
00451
00452 if( opac.opacity_abs[iso.ipIsoLevNIonCon[ipISO][nelem][n]-1] > 0. )
00453 {
00454 #if RJRWLOOP
00455
00456
00457 long int ip = iso.ipIsoLevNIonCon[ipISO][nelem][n];
00458 double t2 = csphot(
00459 ip ,
00460 ip ,
00461 iso.ipOpac[ipISO][nelem][n] );
00462
00463 iso.ConOpacRatio[ipISO][nelem][n] = 1.f-(float)(
00464 (iso.Pop2Ion[ipISO][nelem][n]*dense.xIonDense[nelem][nelem+1-ipISO]*t2)/
00465 opac.opacity_abs[ip-1]);
00466 #else
00467
00468
00469 double t2 = csphot(
00470 iso.ipIsoLevNIonCon[ipISO][nelem][n] ,
00471 iso.ipIsoLevNIonCon[ipISO][nelem][n] ,
00472 iso.ipOpac[ipISO][nelem][n] );
00473
00474 iso.ConOpacRatio[ipISO][nelem][n] = (float)(
00475 (opac.opacity_abs[iso.ipIsoLevNIonCon[ipISO][nelem][n]-1] -
00476 iso.Pop2Ion[ipISO][nelem][n]*dense.xIonDense[nelem][nelem+1-ipISO]*t2)/
00477 opac.opacity_abs[iso.ipIsoLevNIonCon[ipISO][nelem][n]-1]);
00478 #endif
00479 }
00480 else
00481 {
00482 iso.ConOpacRatio[ipISO][nelem][n] = 0.;
00483 }
00484 }
00485 }
00486 }
00487
00488
00489 for( n=0; n < iso.numLevels_local[ipISO][nelem]; n++ )
00490 {
00491
00492 if(iso.ipIsoLevNIonCon[ipISO][nelem][n] < rfield.nflux )
00493 {
00494
00495
00496
00497
00498 if( opac.opacity_abs[iso.ipIsoLevNIonCon[ipISO][nelem][n]-1] > 0. )
00499 {
00500
00501 if( iso.DepartCoef[ipISO][nelem][n] > 1e-30 && (!conv.lgSearch ) )
00502 {
00503
00504 fac = 1./iso.DepartCoef[ipISO][nelem][n];
00505 }
00506 else if( conv.lgSearch )
00507 {
00508
00509
00510 fac = 0.;
00511 }
00512 else
00513 {
00514 fac = 1.;
00515 }
00516
00519
00520
00521
00522
00523
00524 if( opac.opacity_abs[iso.ipIsoLevNIonCon[ipISO][nelem][n]-1] > 0. )
00525 {
00526
00528 #if RJRWMACRO
00529 long int ip = iso.ipIsoLevNIonCon[ipISO][nelem][n];
00530
00531 double t2 = csphot(
00532 ip ,
00533 ip ,
00534 iso.ipOpac[ipISO][nelem][n] );
00535
00536 double opacity_this_species =
00537 iso.Pop2Ion[ipISO][nelem][n]*dense.xIonDense[nelem][nelem+1-ipISO]*t2*
00538 (1. - fac*rfield.ContBoltz[ip-1]);
00539
00540 double opacity_fraction = 1. - opacity_this_species / opac.opacity_abs[ip-1];
00541 if(opacity_fraction < 0)
00542 opacity_fraction = 0.;
00543 #else
00544 double t2 = csphot(
00545 iso.ipIsoLevNIonCon[ipISO][nelem][n] ,
00546 iso.ipIsoLevNIonCon[ipISO][nelem][n] ,
00547 iso.ipOpac[ipISO][nelem][n] );
00548
00549 int ip = iso.ipIsoLevNIonCon[ipISO][nelem][n]-1;
00550
00551 double opacity_this_species =
00552 iso.Pop2Ion[ipISO][nelem][n]*dense.xIonDense[nelem][nelem+1-ipISO]*t2*
00553 (1. - fac*rfield.ContBoltz[iso.ipIsoLevNIonCon[ipISO][nelem][n]-1]);
00554
00555 double opacity_fraction =
00556
00557 MAX2( 0. , opac.opacity_abs[ip] - opacity_this_species ) / opac.opacity_abs[ip];
00558 #endif
00559
00560 iso.ConOpacRatio[ipISO][nelem][n] = (float)(
00561 iso.ConOpacRatio[ipISO][nelem][n]* 0.75 + 0.25*opacity_fraction );
00562
00563
00564
00565
00566
00567
00568
00569
00570
00571
00572
00573
00574 #if RJRWMACRO
00575 if(iso.ConOpacRatio[ipISO][nelem][n] < 0.)
00576 iso.ConOpacRatio[ipISO][nelem][n] = 0.;
00577 #endif
00578 }
00579 else
00580 {
00581 iso.ConOpacRatio[ipISO][nelem][n] = 0.;
00582 }
00583 #if (!RJRWMACRO)
00584 iso.ConOpacRatio[ipISO][nelem][n] = (float)MAX2(0.,iso.ConOpacRatio[ipISO][nelem][n]);
00585 #endif
00586 }
00587 else
00588 {
00589 iso.ConOpacRatio[ipISO][nelem][n] = 0.;
00590 }
00591 }
00592 else
00593 {
00594 iso.ConOpacRatio[ipISO][nelem][n] = 0.;
00595 }
00596 }
00597 }
00598 }
00599 }
00600
00601
00602 if( gv.lgDustOn )
00603 {
00604
00605
00606 for( i=0; i < rfield.nflux; i++ )
00607 {
00608
00609 opac.opacity_sct[i] += gv.dstsc[i]*dense.gas_phase[ipHYDROGEN];
00610 opac.opacity_abs[i] += gv.dstab[i]*dense.gas_phase[ipHYDROGEN];
00611 }
00612 }
00613
00614
00615 for( i=0; i < rfield.nflux; i++ )
00616 {
00617
00618 opac.opacity_abs[i] += opac.OpacStatic[i];
00619
00620
00621 }
00622
00623
00624 for( i=0; i < rfield.nflux; i++ )
00625 {
00626 opac.albedo[i] = opac.opacity_sct[i]/
00627 (opac.opacity_sct[i] + opac.opacity_abs[i]);
00628 }
00629
00630
00631 if( conv.lgSearch )
00632 OpacityZeroOld();
00633
00634 if( trace.lgTrace )
00635 {
00636 fprintf( ioQQQ, " OpacityAddTotal returns; grd rec eff (opac) for Hn=1,4%10.2e%10.2e%10.2e%10.2e HeI,II;%10.2e\n",
00637 iso.ConOpacRatio[ipH_LIKE][ipHYDROGEN][ipH1s],
00638 iso.ConOpacRatio[ipH_LIKE][ipHYDROGEN][ipH2p],
00639 iso.ConOpacRatio[ipH_LIKE][ipHYDROGEN][3],
00640 iso.ConOpacRatio[ipH_LIKE][ipHYDROGEN][4],
00641 iso.ConOpacRatio[ipHE_LIKE][ipHELIUM][ipHe1s1S] );
00642 }
00643 {
00644
00645
00646 enum {DEBUG_LOC=false};
00647
00648 if( DEBUG_LOC && (nzone>=378) )
00649 {
00650 if( nzone > 380 )
00651 exit(1);
00652 for( i=0; i<rfield.nflux; ++i )
00653 {
00654 fprintf(ioQQQ,"rtotsbugggg\t%li\t%.3e\t%.3e\t%.3e\t%.3e\n",
00655 conv.nPres2Ioniz,
00656 rfield.anu[i],
00657 opac.opacity_abs[i],
00658 rfield.otscon[i],
00659 rfield.otslin[i]);
00660 }
00661 }
00662 }
00663
00664 DEBUG_EXIT( "OpacityAddTotal()" );
00665 return;
00666 }
00667