cloudy
trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /*OpacityAddTotal derive total opacity for this position, 00004 * called by ConvBase */ 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 00022 void OpacityAddTotal(void) 00023 { 00024 long int i, 00025 ion, 00026 ipHi, 00027 ipISO, 00028 ipop, 00029 limit, 00030 low, 00031 nelem, 00032 n; 00033 double DepartCoefInv , 00034 fac, 00035 sum; 00036 realnum SaveOxygen1 , 00037 SaveCarbon1; 00038 00039 DEBUG_ENTRY( "OpacityAddTotal()" ); 00040 00041 /* OpacityZero will zero out scattering and absorption opacities, 00042 * and set OldOpacSave to opac to save it */ 00043 OpacityZero(); 00044 00045 /* free electron scattering opacity, Compton recoil energy loss */ 00046 /* >>chng 02 mar 26, make this loop only one over electrons alone */ 00047 /*for( i=0; i < (ionbal.ipCompRecoil[ipHYDROGEN][0] - 1); i++ )*/ 00048 for( i=0; i < rfield.nflux; i++ ) 00049 { 00050 /* scattering part of total opacity */ 00051 opac.opacity_sct[i] += opac.OpacStack[i-1+opac.iopcom]* 00052 dense.eden; 00053 } 00054 00055 /* opacity due to compton bound recoil ionization */ 00056 /* >>chng 01 dec 19, rewrite as loop over all elements */ 00057 /* >>chng 02 mar 26, add in ions */ 00058 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem ) 00059 { 00060 if( dense.lgElmtOn[nelem] ) 00061 { 00062 for( ion=0; ion<nelem+1; ++ion ) 00063 { 00064 realnum factor = dense.xIonDense[nelem][ion]; 00065 /*>>chng 05 nov 26, add molecular hydrogen assuming same 00066 * as two free hydrogen atoms - hence mult density by two 00067 *>>KEYWORD H2 bound Compton ionization */ 00068 if( nelem==ipHYDROGEN ) 00069 factor += hmi.H2_total*2.f; 00070 if( factor > 0. ) 00071 { 00072 // loop_min and loop_max are needed to work around a bug in icc 10.0 00073 long loop_min = ionbal.ipCompRecoil[nelem][ion]-1; 00074 long loop_max = rfield.nflux; 00075 /* ionbal.nCompRecoilElec number of electrons in valence shell 00076 * that can compton recoil ionize */ 00077 factor *= ionbal.nCompRecoilElec[nelem-ion]; 00078 for( i=loop_min; i < loop_max; i++ ) 00079 { 00080 /* add in bound hydrogen electron scattering, treated as absorption */ 00081 opac.opacity_abs[i] += opac.OpacStack[i-1+opac.iopcom]*factor; 00082 } 00083 } 00084 } 00085 } 00086 } 00087 00088 /* opacity due to pair production - does not matter what form these 00089 * elements are in */ 00091 sum = dense.gas_phase[ipHYDROGEN] + 4.*dense.gas_phase[ipHELIUM]; 00092 OpacityAdd1Subshell(opac.ioppr,opac.ippr,rfield.nflux,(realnum)sum,'s'); 00093 00094 /* free-free free free brems emission for all ions */ 00095 00096 /* >>chng 02 jul 21, use full set of ions and gaunt factor */ 00097 /* ipEnergyBremsThin is index to energy where gas goes optically thin to brems, 00098 * so this loop is over optically thick frequencies */ 00099 static double *TotBremsAllIons; 00100 static bool lgFirstTime=true; 00101 double BremsThisEner,bfac, bhfac,sumion[LIMELM+1]; 00102 long int ion_lo , ion_hi; 00103 00104 if( lgFirstTime ) 00105 { 00106 /* rfield.nupper will not change in one coreload, so just malloc this once */ 00107 TotBremsAllIons = (double *)MALLOC((unsigned long)rfield.nupper*sizeof(double)); 00108 lgFirstTime = false; 00109 } 00110 bfac = (dense.eden/1e20)/phycon.sqrte/1e10; 00111 /* gaunt factors depend only on photon energy and ion charge, so do 00112 * sum of ions here before entering into loop over photon energy */ 00113 sumion[0] = 0.; 00114 for(ion=1; ion<=LIMELM; ++ion ) 00115 { 00116 sumion[ion] = 0.; 00117 for( nelem=ipHELIUM; nelem < LIMELM; ++nelem ) 00118 { 00119 if( dense.lgElmtOn[nelem] && ion<=nelem+1 ) 00120 { 00121 sumion[ion] += dense.xIonDense[nelem][ion]; 00122 } 00123 } 00124 /* now include the charge, density, and temperature */ 00125 sumion[ion] *= POW2((double)ion)*bfac; 00126 } 00127 /* now find lowest and highest ion we need to consider - following loop is over 00128 * full continuum and eats time 00129 * >>chng 05 oct 20, following had tests on ion being within bounds, bug, 00130 * changed to ion_lo and ion_hi */ 00131 ion_lo = 1; 00132 while( sumion[ion_lo]==0 && ion_lo<LIMELM-1 ) 00133 ++ion_lo; 00134 ion_hi = LIMELM; 00135 while( sumion[ion_hi]==0 && ion_hi>0 ) 00136 --ion_hi; 00137 00138 bhfac = bfac*(dense.xIonDense[ipHYDROGEN][1]+hmi.Hmolec[ipMH2p]+hmi.Hmolec[ipMH3p]); 00139 for( i=0; i < rfield.nflux; i++ ) 00140 { 00141 /* do hydrogen first, before main loop since want to add on H- brems */ 00142 nelem = ipHYDROGEN; 00143 ion = 1;/* >>chng 02 nov 07, add charged ions as H+ */ 00144 00145 /* for case of hydrogen, add H- brems - OpacStack contains the ratio 00146 * of the H- to H brems cross section - multiply by this and the fraction in ground */ 00147 BremsThisEner = bhfac * rfield.gff[ion][i]* 00148 (1.+opac.OpacStack[i-1+opac.iphmra]*StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop); 00149 /*TotBremsAllIons[i] = BremsThisIon;*/ 00150 /* there is at least one standard test (grainlte) which has ZERO ionization - 00151 * this assert will pass that test (the ==0) and also the usual case */ 00152 /* ASSERT( (dense.xIonDense[nelem][ion]==0.) || (TotBremsAllIons[i] > 0.) );*/ 00153 00154 /* chng 02 may 16, by Ryan...do all brems for all ions in one fell swoop, 00155 * using gaunt factors from rfield.gff. */ 00156 /* >>chng 05 jul 11 reorganize loop for speed */ 00157 for(ion=ion_lo; ion<=ion_hi; ++ion ) 00158 { 00159 BremsThisEner += sumion[ion]*rfield.gff[ion][i]; 00160 } 00161 TotBremsAllIons[i] = BremsThisEner; 00162 00163 /* >>>chng 05 jul 12, bring these two loops together */ 00164 if( rfield.ContBoltz[i] < 0.995 ) 00165 { 00166 TotBremsAllIons[i] *= opac.OpacStack[i-1+opac.ipBrems]* 00167 (1. - rfield.ContBoltz[i]); 00168 } 00169 else 00170 { 00171 TotBremsAllIons[i] *= opac.OpacStack[i-1+opac.ipBrems]* 00172 rfield.anu[i]*TE1RYD/phycon.te; 00173 } 00174 opac.FreeFreeOpacity[i] = TotBremsAllIons[i]; 00175 /* >>chng 02 jul 23, from >0 to >=0, some models have no ionization, 00176 * like grainlte.in */ 00177 /*ASSERT( (opac.FreeFreeOpacity[i] > 0.) || (dense.xIonDense[ipHYDROGEN][1] == 0.) );*/ 00178 opac.opacity_abs[i] += opac.FreeFreeOpacity[i]; 00179 } 00180 00181 00182 /* H minus absorption, with correction for stimulated emission */ 00183 if( hmi.hmidep > SMALLFLOAT ) 00184 { 00185 DepartCoefInv = 1./hmi.hmidep; 00186 } 00187 else 00188 { 00189 /* the hmidep departure coef can become vastly small in totally 00190 * neutral gas (no electrons) */ 00191 DepartCoefInv = 1.; 00192 } 00193 00194 limit = iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0]; 00195 if(limit > rfield.nflux) 00196 limit = rfield.nflux; 00197 00198 for( i=hmi.iphmin-1; i < limit; i++ ) 00199 { 00200 double factor; 00201 factor = 1. - rfield.ContBoltz[i]*DepartCoefInv; 00202 if(factor > 0) 00203 opac.opacity_abs[i] += opac.OpacStack[i-hmi.iphmin+opac.iphmop]* 00204 hmi.Hmolec[ipMHm]*factor; 00205 } 00206 00207 /* H2P h2plus bound free opacity */ 00208 limit = opac.ih2pnt[1]; 00209 if(limit > rfield.nflux) 00210 limit = rfield.nflux; 00211 for( i=opac.ih2pnt[0]-1; i < limit; i++ ) 00212 { 00213 opac.opacity_abs[i] += hmi.Hmolec[ipMH2p]*opac.OpacStack[i-opac.ih2pnt[0]+opac.ih2pof]; 00214 } 00215 00216 /* get total population of hydrogen ground to do Rayleigh scattering */ 00217 if( dense.xIonDense[ipHYDROGEN][1] <= 0. ) 00218 { 00219 fac = dense.xIonDense[ipHYDROGEN][0]; 00220 } 00221 else 00222 { 00223 fac = dense.xIonDense[ipHYDROGEN][1]*StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop; 00224 } 00225 00226 /* Ly a damp wing opac (Rayleigh scattering) */ 00227 limit = iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]; 00228 if(limit > rfield.nflux) 00229 limit = rfield.nflux; 00230 for( i=0; i < limit; i++ ) 00231 { 00232 opac.opacity_sct[i] += (fac*opac.OpacStack[i-1+opac.ipRayScat]); 00233 } 00234 00235 /* remember largest correction for stimulated emission */ 00236 if( iso.DepartCoef[ipH_LIKE][ipHYDROGEN][ipH1s] > 1e-30 && !conv.lgSearch ) 00237 { 00238 realnum factor; 00239 factor = (realnum)(rfield.ContBoltz[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1]/iso.DepartCoef[ipH_LIKE][ipHYDROGEN][ipH1s]); 00240 if(opac.stimax[0] < factor) 00241 opac.stimax[0] = factor; 00242 } 00243 00244 if( iso.DepartCoef[ipH_LIKE][ipHYDROGEN][ipH2p] > 1e-30 && !conv.lgSearch ) 00245 { 00246 realnum factor; 00247 factor = (realnum)(rfield.ContBoltz[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH2p]-1]/iso.DepartCoef[ipH_LIKE][ipHYDROGEN][ipH2p]); 00248 if(opac.stimax[1] < factor) 00249 opac.stimax[1] = factor; 00250 } 00251 00252 # if 0 00253 /* check whether hydrogen or Helium singlets mased, if not in search mode */ 00254 if( !conv.lgSearch ) 00255 { 00256 if( Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].PopOpc < 0. ) 00257 { 00258 hydro.lgHLyaMased = true; 00259 } 00260 } 00261 # endif 00262 00263 /* >>chng 05 nov 25, use Yan et al. H2 photo cs 00264 * following reference gives cross section for all energies 00265 * >>refer H2 photo cs Yan, M., Sadeghpour, H.R., & Dalgarno, A., 1998, ApJ, 496, 1044 00266 * Wilms, J., Allen, A., & McCray, R. 2000, ApJ, 542, 914 */ 00267 /* >>chng 02 jan 16, approximate inclusion of H_2 photoelectric opacity 00268 * include H_2 in total photoelectric opacity as twice H0 cs */ 00269 /* set lower and upper limits to this range */ 00270 /*>>KEYWORD H2 photoionization opacity */ 00271 low = opac.ipH2_photo_thresh; 00272 ipHi = rfield.nupper; 00273 ipop = opac.ipH2_photo_opac_offset; 00274 /* OpacityAdd1Subshell just returns for static opacities if opac.lgRedoStatic not set*/ 00275 /* >>chng 05 nov 27, change on nov 25 had left 2*density from H0, so 00276 * twice the H2 density was used - * also changed to static opacity 00277 * this assumes that all v,J levels contribute the same opacity */ 00278 OpacityAdd1Subshell( ipop , low , ipHi , hmi.H2_total , 's' ); 00279 00280 /*>>KEYWORD CO photoionization opacity */ 00281 /* include photoionization of CO - assume C and O in CO each have 00282 * same photo cs as atom - this should only be significant in highly 00283 * shielded regions where only very hard photons penetrate 00284 * also H2O condensed onto grain surfaces - very important deep in cloud */ 00285 SaveOxygen1 = dense.xIonDense[ipOXYGEN][0]; 00286 SaveCarbon1 = dense.xIonDense[ipCARBON][0]; 00287 dense.xIonDense[ipOXYGEN][0] += findspecies("CO")->hevmol + findspecies("H2Ogrn")->hevmol; 00288 dense.xIonDense[ipCARBON][0] += findspecies("CO")->hevmol; 00289 00290 /* following loop adds standard opacities for first 30 elements 00291 * most heavy element opacity added here */ 00292 for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ ) 00293 { 00294 /* this element may be turned off */ 00295 if( dense.lgElmtOn[nelem] ) 00296 { 00297 OpacityAdd1Element(nelem); 00298 } 00299 } 00300 00301 /* now reset the abundances */ 00302 dense.xIonDense[ipOXYGEN][0] = SaveOxygen1; 00303 dense.xIonDense[ipCARBON][0] = SaveCarbon1; 00304 00305 /* following are opacities due to specific excited levels */ 00306 00307 /* nitrogen opacity 00308 * excited level of N+ */ 00309 OpacityAdd1Subshell(opac.in1[2],opac.in1[0],opac.in1[1], 00310 dense.xIonDense[ipNITROGEN][0]*atoms.p2nit , 'v' ); 00311 00312 /* oxygen opacity 00313 * excited level of Oo */ 00314 OpacityAdd1Subshell(opac.ipo1exc[2],opac.ipo1exc[0],opac.ipo1exc[1], 00315 dense.xIonDense[ipOXYGEN][0]*oxy.poiexc,'v'); 00316 00317 /* O2+ excited states */ 00318 OpacityAdd1Subshell(opac.ipo3exc[2],opac.ipo3exc[0],opac.ipo3exc[1], 00319 dense.xIonDense[ipOXYGEN][2]*oxy.poiii2,'v'); 00320 00321 OpacityAdd1Subshell(opac.ipo3exc3[2],opac.ipo3exc3[0],opac.ipo3exc3[1], 00322 dense.xIonDense[ipOXYGEN][2]*oxy.poiii3,'v'); 00323 00324 /* magnesium opacity 00325 * excited level of Mg+ */ 00326 OpacityAdd1Subshell(opac.ipOpMgEx,opac.ipmgex,iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s], 00327 dense.xIonDense[ipMAGNESIUM][1]* atoms.popmg2,'v'); 00328 00329 /* calcium opacity 00330 * photoionization of excited levels of Ca+ */ 00331 OpacityAdd1Subshell(opac.ica2op,opac.ica2ex[0],opac.ica2ex[1], 00332 ca.popca2ex,'v'); 00333 00334 /******************************************************************* 00335 * 00336 * complete evaluation of total opacity by adding in the static part and grains 00337 * 00338 *******************************************************************/ 00339 00340 /* this loop defines the variable iso.ConOpacRatio[ipH_LIKE][nelem][n], 00341 * the ratio of not H to Hydrogen opacity. for grain free environments 00342 * at low densities this is nearly zero. The correction includes 00343 * stimulated emission correction */ 00344 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO ) 00345 { 00346 for( nelem=ipISO; nelem < LIMELM; nelem++ ) 00347 { 00348 /* this element may be turned off */ 00349 if( dense.lgElmtOn[nelem] ) 00350 { 00351 /* this branch is for startup only */ 00352 if( nzone < 1 ) 00353 { 00354 /* >>chng 06 aug 17, should go to numLevels_local instead of _max */ 00355 for( n=0; n < iso.numLevels_local[ipISO][nelem]; n++ ) 00356 { 00357 if(iso.ipIsoLevNIonCon[ipISO][nelem][n] < rfield.nflux ) 00358 { 00359 /*>>chng 04 dec 12, had tested against 1e-30, set to zero if not 00360 * greater than this - caused oscillations as opacity fell below 00361 * and around this value - change to greater than 0 */ 00362 /*if( opac.opacity_abs[iso.ipIsoLevNIonCon[ipISO][nelem][n]-1] > 1e-30 )*/ 00363 if( opac.opacity_abs[iso.ipIsoLevNIonCon[ipISO][nelem][n]-1] > 0. ) 00364 { 00365 /* >>chng 02 may 06, use general form of threshold cs */ 00366 /*double t1 = atmdat_sth(n)/(POW2(nelem+1.-ipISO));*/ 00367 long int ip = iso.ipIsoLevNIonCon[ipISO][nelem][n]; 00368 double t2 = csphot( 00369 ip , 00370 ip , 00371 iso.ipOpac[ipISO][nelem][n] ); 00372 00373 iso.ConOpacRatio[ipISO][nelem][n] = 1.f-(realnum)( 00374 (StatesElem[ipISO][nelem][n].Pop*dense.xIonDense[nelem][nelem+1-ipISO]*t2)/ 00375 opac.opacity_abs[ip-1]); 00376 } 00377 else 00378 { 00379 iso.ConOpacRatio[ipISO][nelem][n] = 0.; 00380 } 00381 } 00382 } 00383 } 00384 /* end branch for startup only, start branch for all zones including startup */ 00385 /* >>chng 06 aug 17, should go to numLevels_local instead of _max */ 00386 for( n=0; n < iso.numLevels_local[ipISO][nelem]; n++ ) 00387 { 00388 /* ratios of other to total opacity for continua of all atoms done with iso model */ 00389 if(iso.ipIsoLevNIonCon[ipISO][nelem][n] < rfield.nflux ) 00390 { 00391 /*>>chng 04 dec 12, had tested against 1e-30, set to zero if not 00392 * greater than this - caused oscillations as opacity fell below 00393 * and around this value - change to greater than 0 */ 00394 /*if( opac.opacity_abs[iso.ipIsoLevNIonCon[ipISO][nelem][n]-1] > 1e-30 )*/ 00395 if( opac.opacity_abs[iso.ipIsoLevNIonCon[ipISO][nelem][n]-1] > 0. ) 00396 { 00397 /* first get departure coef */ 00398 if( iso.DepartCoef[ipISO][nelem][n] > 1e-30 && (!conv.lgSearch ) ) 00399 { 00400 /* this is the usual case, use inverse of departure coef */ 00401 fac = 1./iso.DepartCoef[ipISO][nelem][n]; 00402 } 00403 else if( conv.lgSearch ) 00404 { 00405 /* do not make correction for stim emission during search 00406 * for initial temperature solution, since trys are very non-equil */ 00407 fac = 0.; 00408 } 00409 else 00410 { 00411 fac = 1.; 00412 } 00413 00416 /* now get opaicty ratio with correction for stimulated emission */ 00417 /*>>chng 04 dec 12, had tested against 1e-30, set to zero if not 00418 * greater than this - caused oscillations as opacity fell below 00419 * and around this value - change to greater than 0 */ 00420 /*if( opac.opacity_abs[iso.ipIsoLevNIonCon[ipISO][nelem][n]-1] > 1e-30 )*/ 00421 if( opac.opacity_abs[iso.ipIsoLevNIonCon[ipISO][nelem][n]-1] > 0. ) 00422 { 00423 /* >>chng 02 may 06, use general form of threshold cs */ 00424 long int ip = iso.ipIsoLevNIonCon[ipISO][nelem][n]; 00425 00426 double t2 = csphot( 00427 ip , 00428 ip , 00429 iso.ipOpac[ipISO][nelem][n] ); 00430 00431 double opacity_this_species = 00432 StatesElem[ipISO][nelem][n].Pop*dense.xIonDense[nelem][nelem+1-ipISO]*t2* 00433 (1. - fac*rfield.ContBoltz[ip-1]); 00434 00435 double opacity_fraction = 1. - opacity_this_species / opac.opacity_abs[ip-1]; 00436 if(opacity_fraction < 0) 00437 opacity_fraction = 0.; 00438 00439 /* use mean of old and new ratios */ 00440 iso.ConOpacRatio[ipISO][nelem][n] = (realnum)( 00441 iso.ConOpacRatio[ipISO][nelem][n]* 0.75 + 0.25*opacity_fraction ); 00442 00443 if(iso.ConOpacRatio[ipISO][nelem][n] < 0.) 00444 iso.ConOpacRatio[ipISO][nelem][n] = 0.; 00445 } 00446 else 00447 { 00448 iso.ConOpacRatio[ipISO][nelem][n] = 0.; 00449 } 00450 } 00451 else 00452 { 00453 iso.ConOpacRatio[ipISO][nelem][n] = 0.; 00454 } 00455 } 00456 else 00457 { 00458 iso.ConOpacRatio[ipISO][nelem][n] = 0.; 00459 } 00460 } 00461 } 00462 } 00463 } 00464 00465 /* add dust grain opacity if dust present */ 00466 if( gv.lgDustOn ) 00467 { 00468 /* generate current grain opacities since may be function of depth */ 00469 /* >>chng 01 may 11, removed code to update grain opacities, already done by GrainChargeTemp */ 00470 for( i=0; i < rfield.nflux; i++ ) 00471 { 00472 /* units cm-1 */ 00473 opac.opacity_sct[i] += gv.dstsc[i]*dense.gas_phase[ipHYDROGEN]; 00474 opac.opacity_abs[i] += gv.dstab[i]*dense.gas_phase[ipHYDROGEN]; 00475 } 00476 } 00477 00478 /* check that opacity is sane */ 00479 for( i=0; i < rfield.nflux; i++ ) 00480 { 00481 /* OpacStatic was zeroed in OpacityZero, incremented in opacityadd1subshell */ 00482 opac.opacity_abs[i] += opac.OpacStatic[i]; 00483 /* make sure that opacity is positive */ 00484 /*ASSERT( opac.opacity_abs[i] > 0. );*/ 00485 } 00486 00487 /* compute gas albedo here */ 00488 for( i=0; i < rfield.nflux; i++ ) 00489 { 00490 opac.albedo[i] = opac.opacity_sct[i]/ 00491 (opac.opacity_sct[i] + opac.opacity_abs[i]); 00492 } 00493 00494 /* during search phase set old opacity array to current value */ 00495 if( conv.lgSearch ) 00496 OpacityZeroOld(); 00497 00498 if( trace.lgTrace ) 00499 { 00500 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", 00501 iso.ConOpacRatio[ipH_LIKE][ipHYDROGEN][ipH1s], 00502 iso.ConOpacRatio[ipH_LIKE][ipHYDROGEN][ipH2s], 00503 iso.ConOpacRatio[ipH_LIKE][ipHYDROGEN][ipH2p], 00504 iso.ConOpacRatio[ipH_LIKE][ipHYDROGEN][ipH3s], 00505 iso.ConOpacRatio[ipH_LIKE][ipHYDROGEN][ipH3p], 00506 iso.ConOpacRatio[ipH_LIKE][ipHYDROGEN][ipH3d], 00507 iso.ConOpacRatio[ipH_LIKE][ipHYDROGEN][ipH4s], 00508 iso.ConOpacRatio[ipH_LIKE][ipHYDROGEN][ipH4p], 00509 iso.ConOpacRatio[ipH_LIKE][ipHYDROGEN][ipH4d], 00510 iso.ConOpacRatio[ipH_LIKE][ipHYDROGEN][ipH4f], 00511 iso.ConOpacRatio[ipHE_LIKE][ipHELIUM][ipHe1s1S], 00512 iso.ConOpacRatio[ipH_LIKE][ipHELIUM][ipH1s] ); 00513 } 00514 { 00515 /* following should be set true to print strongest ots contributors */ 00516 /*@-redef@*/ 00517 enum {DEBUG_LOC=false}; 00518 /*@+redef@*/ 00519 if( DEBUG_LOC && (nzone>=378) ) 00520 { 00521 if( nzone > 380 ) 00522 cdEXIT( EXIT_FAILURE ); 00523 for( i=0; i<rfield.nflux; ++i ) 00524 { 00525 fprintf(ioQQQ,"rtotsbugggg\t%li\t%.3e\t%.3e\t%.3e\t%.3e\n", 00526 conv.nPres2Ioniz, 00527 rfield.anu[i], 00528 opac.opacity_abs[i], 00529 rfield.otscon[i], 00530 rfield.otslin[i]); 00531 } 00532 } 00533 } 00534 return; 00535 }