/*SetCon derive intensity of incident continuum */ #include #include "cddefines.h" #include "physconst.h" #include "incidspec.h" #include "nhe1lvl.h" #include "nhe3lvl.h" #include "opacpoint.h" #include "hydrogenic.h" #include "ionfracs.h" #include "egray.h" #include "prhd.h" #include "showme.h" #include "noexec.h" #include "smbeta.h" #include "tepowers.h" #include "compton.h" #include "hextra.h" #include "recoil.h" #include "trace.h" #include "flxsav.h" #include "exptau.h" #include "e2tau.h" #include "i3727.h" #include "ircoil.h" #include "u.h" #include "heavy.h" #include "qs.h" #include "rfield.h" #include "nxray.h" #include "he.h" #include "phycon.h" #include "called.h" #include "occ1hi.h" #include "phe1lv.h" #include "hphoto.h" #include "habing.h" #include "tenden.h" #include "timesc.h" #include "flxfnt.h" #include "totlum.h" #include "neutrn.h" #include "nh.h" #include "occmax.h" #include "h.h" #include "con0.h" #include "nhe1.h" #include "plasnu.h" #include "bit32.h" #include "hionrad.h" #include "secondaries.h" #include "elecden.h" #include "elmton.h" #include "opac.h" #include "he3n.h" #include "forcte.h" #include "linelabl.h" #include "ffun.h" #include "ipoint.h" #include "cfit.h" #include "smeets.h" #include "ptrcer.h" #include "extin.h" #include "sumcon.h" #include "tfidle.h" #include "sumcontinuum.h" #include "velset.h" #include "pressure.h" #include "setcon.h" void SetCon(int *lgOK) { int lgCheckOK; long int i, ip, j, n; float EdenHeav, ex1ryd, factor, occ1, p, p1, p2, p3, p4, p5, p6, p7, pgn, phe, pheii, qgn, temp; float xIoniz; double rec, wanu[4], alf, bet, fntest, fsum, ecrit, tcompr, tcomp, r2ov1, r3ov2; static double aweigh[4], fweigh[4]; double amean, amean2, amean3, peak, wfun[4]; static int _aini = 1; if( _aini ){ /* Do 1 TIME INITIALIZATIONS! */ aweigh[0] = -0.4305682; aweigh[1] = -0.1699905; aweigh[2] = 0.1699905; aweigh[3] = 0.4305682; fweigh[0] = 0.1739274; fweigh[1] = 0.3260726; fweigh[2] = 0.3260726; fweigh[3] = 0.1739274; _aini = 0; } # ifdef DEBUG_FUN fputs( "<+>SetCon()\n", debug_fp ); # endif /* set continuum */ if( trace.lgTrace ) { fprintf( ioQQQ, " SetCon called.\n" ); } /* define factors to convert FLUX array into photon occupation array OCCNUM * by multiplication */ factor = (float)(EN1RYD/PI4/FR1RYD/HNU3C2); /*------------------------------------------------------------- */ lgCheckOK = TRUE; for( i=0; i < rfield.nupper; i++ ) { /* this was original anu array with no continuum averaging */ rfield.anu[i] = rfield.AnuOrg[i]; rfield.ContBoltz[i] = 0.; opac.tmn[i] = 1.; exptau.ExpmTau[i] = 1.; e2tau.e2TauAbs[i] = 1.; fsum = 0.; amean = 0.; amean2 = 0.; amean3 = 0.; for( j=0; j < 4; j++ ) { wanu[j] = rfield.anu[i] + rfield.widflx[i]*aweigh[j]; wfun[j] = fweigh[j]*ffun(wanu[j]); fsum += wfun[j]; amean += wanu[j]*wfun[j]; amean2 += wanu[j]*wanu[j]*wfun[j]; amean3 += wanu[j]*wanu[j]*wanu[j]*wfun[j]; } rfield.flux[i] = (float)(fsum*rfield.widflx[i]); if( rfield.flux[i] > 0. ) { rfield.anu[i] = (float)(amean2/amean); rfield.anu2[i] = (float)(amean3/amean); /* fix conversion factor for occupation number */ } else if( rfield.flux[i] == 0. ) { rfield.anu2[i] = rfield.anu[i]*rfield.anu[i]; } else { rfield.anu2[i] = rfield.anu[i]*rfield.anu[i]; fprintf( ioQQQ, " negative continuum returned at%6ld%10.2e%10.2e\n", i, rfield.anu[i], rfield.flux[i] ); lgCheckOK = FALSE; } rfield.anu3[i] = rfield.anu2[i]*rfield.anu[i]; rfield.refcon[i] = 0.; rfield.convoc[i] = factor/rfield.widflx[i]/rfield.anu2[i]; /* following are Compton exchange factors from Tarter */ alf = 1./(1. + rfield.anu[i]*(1.1792e-4 + 7.084e-10*rfield.anu[i])); bet = 1. - alf*rfield.anu[i]*(1.1792e-4 + 2.*7.084e-10*rfield.anu[i])/ 4.; Compton.csigh[i] = (float)(alf*rfield.anu2[i]*3.858e-25); Compton.csigc[i] = (float)(alf*bet*rfield.anu[i]*3.858e-25); } if( !lgCheckOK ) { ShowMe(); puts( "[Stop in setcon]" ); exit(1); } if( trace.lgTrace && trace.lgComBug ) { fprintf( ioQQQ, "\n\n Compton heating, cooling coefficients \n" ); for( i=0; i < rfield.nupper; i += 2 ) { fprintf( ioQQQ, "%6ld%10.2e%10.2e%10.2e", i, rfield.anu[i], Compton.csigh[i], Compton.csigc[i] ); } fprintf( ioQQQ, "\n" ); } /* option to check frequencies in real time, drive pointers command */ if( trace.lgPtrace ) ptrcer(); for( i=0; i < rfield.nupper; i++ ) { /* define array of LOG10( nu(ryd)) */ rfield.anulog[i] = (float)log10(rfield.anu[i]); } /* extinguish continuum if set on */ extin(&ex1ryd); /* now find peak of hydrogen ionizing continuum - for PDR calculations * this will remain equal to 1 since the loop will not execute */ prhd.ipeak = 1; peak = 0.; for( i=nh.ipHn[0][IP1S]-1; i < rfield.nupper; i++ ) { if( rfield.flux[i]*rfield.anu[i]/rfield.widflx[i] > peak ) { /* prhd.ipeak points to largest f_nu at H-ionizing energies * and is passed to other parts of code */ /* i+1 to keep ipeak on fortran version of energy array */ prhd.ipeak = i+1; peak = rfield.flux[i]*rfield.anu[i]/rfield.widflx[i]; } } /* say what type of cpu this is, if desired */ if( trace.lgTrace ) { fprintf( ioQQQ, " SetCon: The peak of the H-ion continuum is at%10.3e\n", rfield.anu[prhd.ipeak-1] ); if( bit32.lgBit32 ) { fprintf( ioQQQ, " SetCon: this is a 32-bit cpu.\n" ); } else { fprintf( ioQQQ, " SetCon: this is a long-word cpu.\n" ); } } /* find highest energy to consider in continuum flux array * peak is the peak product nu*flux */ peak = rfield.flux[prhd.ipeak-1]/rfield.widflx[prhd.ipeak-1]* rfield.anu2[prhd.ipeak-1]; if( peak > 1e38 && bit32.lgBit32 ) { fprintf( ioQQQ, " SetCon: this appears to be a 32-bit cpu.\n" ); fprintf( ioQQQ, " The continuum is too intense to compute with this cpu. Use a long-word cpu or a fainter continuum. (This is the nu*f_nu test)\n" ); fprintf( ioQQQ, " Sorry.\n" ); puts( "[Stop in setcon]" ); exit(1); } /* FluxFaint set in block data; normally 1e-10 */ /* this will be faiintest level of continuum we want to consider. * peak was set above, is peak of hydrogen ionizing radiation field, * and is zero if no H-ionizing radiation */ fntest = peak*flxfnt.FluxFaint; if( fntest > 0. ) { /* this test is not done in pdr conditions where NO H-ionizing radiation, * since fntest is zero*/ i = rfield.nupper; while( i > prhd.ipeak && rfield.flux[i-1]*rfield.anu2[i-1]/ rfield.widflx[i-1] < fntest ) { --i; } } else { /* when no H-ionizing radiation set to Lyman edge */ i = nh.ipHn[0][IP1S]; } /* * this line of code dates from 1979 and IOA Cambridge. removed july 19 95 * I think it was the last line of the original Cambridge source nflux = MAX( ineon(1)+4 , i ) */ /* >>chng 99 april 28, reinstate the flxfnt limit with nflux */ rfield.nflux = i; /* trim down nflux, was set to NCELL, the dimension of all vectors, in zero.c, * in CreatePoint was set to nupper, the number of cells needed to get up to the * high energy limit of the code */ while( rfield.flux[rfield.nflux-1] < SMALLFLOAT && rfield.nflux > 1 ) { --rfield.nflux; } if( rfield.nflux == 1 ) { fprintf( ioQQQ, " This incident continuum appears to have no radiation.\n" ); fprintf( ioQQQ, " Sorry.\n" ); puts( "[Stop in setcon]" ); exit(1); } /* check that continuum defined everywhere - look for zero's and comment if present */ con0.lgCon0 = FALSE; ip = 0; for( i=0; i < rfield.nflux; i++ ) { if( rfield.flux[i] == 0. ) { if( called.lgTalk && !con0.lgCon0 ) { fprintf( ioQQQ, " Setcon: continuum has zero intensity starting at %11.4e Ryd.\n", rfield.anu[i] ); con0.lgCon0 = TRUE; } ip += 1; } } if( con0.lgCon0 && called.lgTalk ) { fprintf( ioQQQ, "%6ld cells in the incident continuum have zero intensity. Problems???\n", ip ); } /*begin sanity check */ lgCheckOK = TRUE; for( i=1; i < rfield.nflux; i++ ) { if( rfield.flux[i] < 0. ) { fprintf( ioQQQ, " Continuum has negative intensity at%11.4e Ryd=%10.2e %4.4s %4.4s\n", rfield.anu[i], rfield.flux[i], LineLabl.chLineLabel[i] , LineLabl.chContLabel[i] ); lgCheckOK = FALSE; } else if( rfield.anu[i] <= rfield.anu[i-1] ) { fprintf( ioQQQ, " Continuum energies not in increasing order: energies follow\n" ); fprintf( ioQQQ, "%3ld %10.2e%3ld %10.2e%3ld %10.2e\n", i -1 , rfield.anu[i-1], i, rfield.anu[i], i +1, rfield.anu[i+1] ); lgCheckOK = FALSE; } } if( !lgCheckOK ) { ShowMe(); puts( "[Stop in setcon]" ); exit(1); } /*end sanity check */ /* turn off recoil ionization if high energy < 190R */ if( rfield.anu[rfield.nflux-1] <= 190 ) { recoil.lgCompRecoil = FALSE; } /* sum photons and energy, save mean */ /* sum from low energy to Balmer edge */ sumcon(1,nh.ipHn[0][IP2P]-1,&qs.qrad,&prhd.pradio,&p1); /* sum over Balmer continuum */ sumcon(nh.ipHn[0][2],nh.ipHn[0][IP1S]-1,&qs.qbal,&prhd.pbal,&p1); /* sum from Lyman edge to HeI edge */ sumcon(nh.ipHn[0][IP1S],nhe1Com.nhe1[0]-1,&prhd.q,&p,&p2); /* sum from HeI to HeII edges */ sumcon(nhe1Com.nhe1[0],nh.ipHn[1][IP1S]-1,&qs.qhe,&phe,&p3); /* sum from Lyman edge to carbon k-shell */ sumcon(nh.ipHn[1][IP1S],nxray.ipCKshell-1,&qs.qheii,&pheii,&p4); /* sum from c k-shell to gamma ray - where pairs start */ sumcon( nxray.ipCKshell , egray.ipEnerGammaRay-1 , &prhd.qx, &prhd.xpow , &p5); /* complete sum up to high energy limit */ sumcon(egray.ipEnerGammaRay,rfield.nflux,&prhd.qgam,&prhd.GammaLumin, &p6); /* find to estimate photoerosion timescale */ n = ipoint(7.35e5); sumcon(n,rfield.nflux,&qgn,&pgn,&p7); timesc.TimeErode = qgn; /* find Compton temp */ tcompr = (p1 + p2 + p3 + p4 + p5 + p6)/(prhd.pradio + prhd.pbal + p + phe + pheii + prhd.xpow + prhd.GammaLumin); tcomp = tcompr/(4.*6.272e-6); if( trace.lgTrace ) { fprintf( ioQQQ, " mean photon energy=%10.3eR =%10.3eK low, high nu=%12.4e%12.4e\n", tcompr, tcomp, rfield.anu[0] - rfield.widflx[0]/2., rfield.anu[rfield.nflux-1] + rfield.widflx[rfield.nflux-1]/2. ); } /* this is total power in ionizing radiation */ prhd.powion = p + phe + pheii + prhd.xpow + prhd.GammaLumin; /* this is the total photon luminosity */ totlum.TotalLumin = prhd.pradio + prhd.powion + prhd.pbal; totlum.totlsv = totlum.TotalLumin; /* total H-ionizing photon number, */ qs.qhtot = prhd.q + qs.qhe + qs.qheii + prhd.qx + prhd.qgam; /* ftotal photon number, all energies */ qs.qtot = qs.qhtot + qs.qbal + qs.qrad; if( prhd.powion <= 0. && called.lgTalk ) { HionRad.lgHionRad = TRUE; fprintf( ioQQQ, " >>>>There is no hydrogen-ionizing radiation.<<<<\n" ); fprintf( ioQQQ, " >>>>Was this intended?<<<<\n" ); fprintf( ioQQQ, " \n" ); /* >>chng 97 mar 18, add following sanity check - stop if no Paschen * ionizing radiaion since even metals will be totally neutral */ sumcon(nh.ipHn[0][3],nh.ipHn[0][2]-1,&factor,&temp,&p1); if( factor <= 0. ) { fprintf( ioQQQ, " >>>>There is no sodium-ionizing radiation.<<<<\n" ); fprintf( ioQQQ, " >>>>This code is not approriate for these conditions.<<<<\n" ); /* punch out before we crash */ *lgOK = FALSE; # ifdef DEBUG_FUN fputs( " <->SetCon()\n", debug_fp ); # endif return; } } else { HionRad.lgHionRad = FALSE; } /* option to add energy deposition due to fast neutrons, as frc of tot lum * efficiency default is unity */ if( neutrn.lgNeutrnHeatOn ) { neutrn.totneu = (float)(neutrn.effneu*totlum.TotalLumin*pow(10.,neutrn.frcneu)); } else { neutrn.totneu = (float)0.; } /* temp correspond to energy density, printed in STARTR */ tenden.TEnerDen = (float)pow(totlum.TotalLumin/SPEEDLIGHT/7.56464e-15,0.25); /* sanity check for single blackbody, that energy density temperature * is smaller than black body temperature */ if( IncidSpec.nspec==1 && strcmp( spnorm.chSpType[IncidSpec.nspec-1], "BLACK" )==0 ) { /* single black body, now confirm that TEnerDen is less than this temperature, * >>>chng 99 may 2, * in lte these are very very close, factor of 1.00001 allows for numerical * errors, and apparently slightly different atomic coef in different parts * of code. eventaully all mustuse physonst.h and agree exactly */ if( tenden.TEnerDen > 1.0001*IncidSpec.slope[IncidSpec.nspec-1] ) { fprintf( ioQQQ, "\n WARNING: The energy density temperature (%g) is greater than the" " black body temperature (%g). This is unphysical.\n\n", tenden.TEnerDen , IncidSpec.slope[IncidSpec.nspec-1]); } } /* simple estimate of Hbeta flux if all phot go to H */ smbeta.SimHBeta = (float)(qs.qhtot*4.75e-13); /* incident continuum nu*f_nu at Hbeta and Ly-alpha */ smbeta.cn4861 = (float)(ffun(0.1875)*HPLANCK*FR1RYD*0.1875*0.1875); smbeta.cn1216 = (float)(ffun(0.75)*HPLANCK*FR1RYD*0.75*0.75); smbeta.sv4861 = smbeta.cn4861; smbeta.sv1216 = smbeta.cn1216; /* flux density nu*Fnu = erg / s / cm2 * EX1RYD is optional extinction factor at 1 ryd */ prhd.fx1ryd = (float)(ffun(1.000)*HPLANCK*ex1ryd*FR1RYD); /* check for plasma frequency - then zero out incident continuum * for energies below this * this is critical electron density, shielding of incident continuum * if electron density is greater than this */ ecrit = POW2(rfield.anu[0]/2.729e-12); if( phycon.hden*1.2 > ecrit ) { plasnu.lgPlasNu = TRUE; plasnu.plsfrq = (float)(2.729e-12*sqrt(phycon.hden*1.2)); plasnu.plsfrqmax = plasnu.plsfrq; plasnu.ipPlasma = ipoint(plasnu.plsfrq); /* save max pointer too */ plasnu.ipPlasmax = plasnu.ipPlasma; /* now loop over all affected energies, setting incident continuum * to zero there, and counting all as reflected */ for( i=0; i < plasnu.ipPlasma; i++ ) { /* count as reflected continuum */ rfield.refcon[i] = rfield.flux[i]; /* set continuum to zero there */ rfield.flux[i] = 0.; } } else { plasnu.lgPlasNu = FALSE; plasnu.ipPlasma = 0; plasnu.plsfrqmax = 0.; plasnu.plsfrq = 0.; } if( plasnu.ipPlasma > 1 && called.lgTalk ) { fprintf( ioQQQ, " !The plasma frequency is%10.2e Ryd. The incident continuum is set to 0 below this.\n", plasnu.plsfrq ); } occmaxCom.occmax = 0.; occmaxCom.tbrmax = 0.; for( i=0; i < rfield.nupper; i++ ) { /* set up occupation number array */ rfield.occnum[i] = rfield.flux[i]*rfield.convoc[i]; if( rfield.occnum[i] > occmaxCom.occmax ) { occmaxCom.occmax = rfield.occnum[i]; occmaxCom.occmnu = rfield.anu[i]; } /* following product is continuum brightness temperature */ if( rfield.occnum[i]*TE1RYD*rfield.anu[i] > occmaxCom.tbrmax ) { occmaxCom.tbrmax = (float)(rfield.occnum[i]*TE1RYD*rfield.anu[i]); occmaxCom.tbrmnu = rfield.anu[i]; } /* save continuum for next iteration */ flxsav.FluxSave[i] = rfield.flux[i]; } /* if continuum brightness temp is large, where does it fall below * 1e4k??? */ if( occmaxCom.tbrmax > 1e4 ) { i = ipoint(occmaxCom.tbrmnu); while( i < rfield.nupper && (rfield.occnum[i-1]*TE1RYD* rfield.anu[i-1] > 1e4) ) { i += 1; } occmaxCom.tbr4nu = rfield.anu[i-1]; } else { occmaxCom.tbr4nu = 0.; } /* if continuum occ num is large, where does it fall below 1? */ if( occmaxCom.occmax > 1. ) { i = ipoint(occmaxCom.occmnu)-1; while( i < rfield.nupper && (rfield.occnum[i] > 1.) ) { ++i; } occmaxCom.occ1nu = rfield.anu[i]; } else { occmaxCom.occ1nu = 0.; } /* remember if incident radiation field is less than 10*Habing ISM */ if( totlum.TotalLumin < 1.8e-2 ) { habing.lgHabing = TRUE; } /* fix ionization parameter (per hydrogen) at inner edge */ u.uh = (float)(qs.qhtot/phycon.hden/SPEEDLIGHT); u.uheii = (float)((qs.qheii + prhd.qx)/phycon.hden/SPEEDLIGHT); /* guess first temperature and neutral h density */ if( forcte.ForceTemp > 0. ) { phycon.te = forcte.ForceTemp; } else { if( u.uh > 0. ) { phycon.te = (float)(20000.+log10(u.uh)*5000.); phycon.te = (float)MAX2(8000. , phycon.te ); } else { phycon.te = (float)1000.; } } /* this is an option to stop after printing header only */ if( noexec.lgNoExec ) { return; } /* following needed for tfidle */ phycon.eden = phycon.hden; h.hi = 0.; /* this must be zero sinze TotalPressure will do radiation pressure due to H */ hydro.hn[0][IP1S] = 0.; rec = TotalPressure(); tfidle(); velset(); /* estimate secondary ionization rate - probably 0, but possible extra * SetCsupra set with "set csupra" */ /* >>>chng 99 apr 29, added cosmic ray ionization since this is used to get * helium ionization fraction, and was zero in pdr, so He turned off at start, * and never turned back on */ /* coef on cryden is from highen.c */ Secondaries.csupra = Secondaries.SetCsupra + hextra.cryden*2e-9f; /********************************************************************* * * * esimate hydrogen's level of ionization * * * *********************************************************************/ /* first estimate level of hydrogen ionization */ rec = (-9.9765209 + 0.158607055*TePowers.telogn[0] + 0.30112749* TePowers.telogn[1] - 0.063969007*TePowers.telogn[2] + 0.0012691546* TePowers.telogn[3])/(1. + 0.035055422*TePowers.telogn[0] - 0.037621619*TePowers.telogn[1] + 0.0076175048*TePowers.telogn[2] - 0.00023432613*TePowers.telogn[3]); rec = pow(10.,rec)/phycon.te*phycon.eden; xIoniz = (float)(qs.qhtot*2e-18 + smeets(1)*phycon.eden); r2ov1 = xIoniz/rec; h.hi = (float)(phycon.hden/(1. + r2ov1)); h.hii = (float)phycon.hden - h.hi; if( h.hii > 1e-30 ) { hydro.hn[0][IP1S] = h.hi/h.hii; } else { hydro.hn[0][IP1S] = 0.; } /* save into master ionization array */ xIonFracs[1][0] = h.hi; xIonFracs[2][0] = h.hii; /* now save estimates of whether induced recombination is going * to be important -this is a code pacesetter since GammaBn is slower * than GammaK */ HPhoto.lgHInducImp = FALSE; for( i=IP1S; i <= nhlevel; i++ ) { if( rfield.occnum[nh.ipHn[0][i]-1] > 0.01 ) HPhoto.lgHInducImp = TRUE; } /********************************************************************* * * * esimate helium's level of ionization * * * *********************************************************************/ /* only if helium is turned on */ if( elmton.lgElmtOn[1] ) { /* next estimate level of helium singly ionized */ cfit(2,2,phycon.te,&xIoniz); xIoniz = (float)(xIoniz*phycon.eden + qs.qhe*1e-18); r2ov1 = xIoniz/rec; /* now estimate level of helium doubly ionized */ cfit(2,1,phycon.te,&xIoniz); xIoniz = (float)(xIoniz*phycon.eden + qs.qheii*1e-18); /* rough charge dependence */ rec *= 4.; r3ov2 = xIoniz/rec; /* now set level of helium ionization */ if( r2ov1 > 0. ) { xIonFracs[2][1] = (float)(he.ahe/(1./r2ov1 + 1. + r3ov2)); xIonFracs[1][1] = (float)(xIonFracs[2][1]/r2ov1); } else { /* no He ionizing radiation */ xIonFracs[2][1] = 0.; xIonFracs[1][1] = he.ahe; } xIonFracs[3][1] = (float)(xIonFracs[2][1]*r3ov2); he.hei1 = xIonFracs[1][1]; /* save into master ionization array */ he.hei = xIonFracs[1][1]; he.heii = xIonFracs[2][1]; he.heiii = xIonFracs[3][1]; if( xIonFracs[3][1] > 1e-30 ) { hydro.hn[1][IP1S] = xIonFracs[2][1]/xIonFracs[3][1]; } else { hydro.hn[1][IP1S] = 0.; } } else { /* case where helium is turned off */ xIonFracs[2][1] = 0.; xIonFracs[1][1] = 0.; xIonFracs[3][1] = 0.; he.hei1 = 0.; he.hei = 0.; he.heii = 0.; he.heiii = 0.; hydro.hn[1][IP1S] = 0.; } /* estimate electrons from heavies, assuming each at least * 1 times ionized */ EdenHeav = 0.; for( i=2; i < LIMELM; i++ ) { if( elmton.lgElmtOn[i] ) { EdenHeav += xIonFracs[0][i]; } } /* estimate of electron density */ phycon.eden = (float)(xIonFracs[2][0] + xIonFracs[2][1] + 2.*xIonFracs[3][1] + EdenHeav + ElecDen.EdenExtra); if( ElecDen.EdenSet > 0. ) { phycon.eden = ElecDen.EdenSet; } ElecDen.EdenHCorr = (float)phycon.eden; /* carbon is nearly always at least once ionized */ if( elmton.lgElmtOn[5] ) { ElecDen.EdenMet = xIonFracs[0][5]; ElecDen.EdenNotH = (float)phycon.eden - h.hii; } else { ElecDen.EdenMet = 0.; ElecDen.EdenNotH = 0.; } if( phycon.eden < 0. ) { fprintf( ioQQQ, " Negative electron density results in SetCon.\n" ); fprintf( ioQQQ, "%10.2e%10.2e%10.2e%10.2e%10.2e%10.2e\n", phycon.eden, xIonFracs[2][0], xIonFracs[2][1], xIonFracs[3][1], xIonFracs[0][5], ElecDen.EdenExtra ); ShowMe(); puts( "[Stop in setcon]" ); exit(1); } ElecDen.EdenTrue = (float)phycon.eden; if( trace.lgTrace && trace.lgNeBug ) { fprintf( ioQQQ, " EDEN set to%12.4e by SetCon.\n", phycon.eden ); } /* he triplets only one that needs eden */ he.hei3 = (float)(xIonFracs[2][1]*5.79e-2/phycon.te*(1. + 3110* pow(phycon.te/1e4,-0.51)/phycon.eden)); /* for later line opacities we need ratioal estimates */ if( xIonFracs[2][1] > 0. ) { /* if heii is 0 then leave at values set in block data */ phe1lv.he1n[0] = xIonFracs[1][1]/xIonFracs[2][1]; he3nCom.he3n[0] = he.hei3/xIonFracs[2][1]; } occ1 = (float)(prhd.fx1ryd/HNU3C2/PI4/FR1RYD); /* what is occupation number at 1 Ryd? */ if( occ1 > 1. ) { occ1hi.lgOcc1Hi = TRUE; } else { occ1hi.lgOcc1Hi = FALSE; } if( trace.lgTrace && trace.lgConBug ) { /* print some useful pointers to ionization edges */ fprintf( ioQQQ, " H2,1=%5ld%5ld NX=%5ld IRC=%5ld\n", nh.ipHn[0][2], nh.ipHn[0][IP1S], nxray.ipCKshell, ircoil.ipRecoilIon ); fprintf( ioQQQ, " HE3, 1=%5ld%5ld HE2=%5ld\n", he.nhei3, nhe1Com.nhe1[0], nh.ipHn[1][IP1S] ); fprintf( ioQQQ, " CARBON" ); for( i=0; i < 6; i++ ) { fprintf( ioQQQ, "%5ld", Heavy.ipHeavy[i][5] ); } fprintf( ioQQQ, "\n" ); fprintf( ioQQQ, " OXY" ); for( i=0; i < 8; i++ ) { fprintf( ioQQQ, "%5ld", Heavy.ipHeavy[i][7] ); } fprintf( ioQQQ, "%5ld%5ld%5ld\n", OpacPoint.ipo3exc[0], i3727.i2d, i3727.i2p ); fprintf( ioQQQ, "\n\n PHOTONS PER CELL (NOT RYD)\n" ); fprintf( ioQQQ, "\n\n nu, flux, wid, occ \n" ); fprintf( ioQQQ, " " ); for( i=0; i < rfield.nflux; i++ ) { fprintf( ioQQQ, "%4ld%10.2e%10.2e%10.2e%10.2e", i, rfield.anu[i], rfield.flux[i], rfield.widflx[i], rfield.occnum[i] ); } fprintf( ioQQQ, " \n" ); } /* zero out some continua related to the ots rates, * prototype and routine in SumContinuum. This is done here since summed cont will * be set to rfield */ ZeroOTSContinuum(); if( trace.lgTrace ) { fprintf( ioQQQ, " SetCon returns, nflux=%5ld anu(nflux)=%11.4e eden=%10.2e\n", rfield.nflux, rfield.anu[rfield.nflux-1], phycon.eden ); } *lgOK = TRUE; # ifdef DEBUG_FUN fputs( " <->SetCon()\n", debug_fp ); # endif return; }