/*SumHeat evaluate heating and secondary ionization for current conditions */ /*ZeroHeat is called by ionize */ #include "cddefines.h" #include "physconst.h" #include "ionrange.h" #include "coolants.h" #include "heat.h" #include "collionrate.h" #include "nsshells.h" #include "trace.h" #include "secondaries.h" #include "showme.h" #include "coolheavy.h" #include "heating.h" #include "hmi.h" #include "hydrogenic.h" #include "elmton.h" #include "tcool.h" #include "hcolrec.h" #include "photrate.h" #include "ionfracs.h" #include "phycon.h" #include "zonecnt.h" #include "numderiv.h" #include "holod.h" #include "hpheat.h" #include "hchargtran.h" #include "sumheat.h" /* this is the faintest relative heating we will print */ #define FAINT_HEAT 0.02 void SumHeat( void ) { /* use to dim some vectors */ # ifdef NDIM # undef NDIM # endif # define NDIM 40 long int i, ipnt, ipsave[NDIM], ipZ, j, jpsave[NDIM], limit , nion; double HeatingLo , HeatingHi , secmet , smetla; long ns; static long int nhit=0, nzSave=0; double PhotoHeat, OtherHeat , deriv, oldfac, save[NDIM]; static double oldheat=0., oldtemp=0.; # ifdef DEBUG_FUN fputs( "<+>SumHeat()\n", debug_fp ); # endif /* routine to sum up total heating, and print agents if needed * it also computes the true deriv, dH/dT */ double xe; # ifdef DEBUG_FUN fputs( "<+>nockon()\n", debug_fp ); # endif /******************************************************************* * * reevaluate the secondary ionization and excitation rates * *******************************************************************/ /* phycon.ElecFrac is electron fraction */ /* XE = EDEN/(total neutrals, roughly H) * analytic fits to Shull and Van Steenberg (1985; Ap.J. 298, 268). * lgSecOFF turns off secondary ionizations, sets heating effic to 100% */ if( Secondaries.lgSecOFF || phycon.ElecFrac > 0.95 ) { Secondaries.heatef = 1.; Secondaries.efionz = 0.; Secondaries.exctef = 0.; } else { /* at very low ionization - as per * >>>refer Xu and McCray 1991, Ap.J. 375, 190. * everything goes to asymptote that is not present in Shull and * Van Steenberg - do this by not letting phycon.ElecFrac fall below 1e-4 */ xe = MAX2(phycon.ElecFrac,1e-4); /* this is heating efficiency for high-energy photoejections. Fully * ionized, this is unity. It is a scale factor that multiplies the * high energy heating rate */ Secondaries.heatef = (float)(0.9971*(1. - pow(1. - pow(xe,0.2663),1.3163))); /* secondary ionizations - this is the number of secondary ionizations * produced for each ryd of heat input. it only multiplies the high-energy * heating rate. It is zero for a fully ionized gas */ Secondaries.sec2prim = (float)(0.3908*pow(1. - pow(xe,0.4092),1.7592)); /* by dividing by the energy of one ryd, we can use this factor to multiply * the heating rate in ergs */ Secondaries.efionz = (float)(Secondaries.sec2prim/EN1RYD); /* secondary excitation of L-alpha, actually all Lyman lines */ Secondaries.exctef = (float)(0.4766*pow(1. - pow(xe,0.2735),1.5221)); Secondaries.exctef /= (float)EN1RYD; } /******************************************************************* * * get total heating from all species * *******************************************************************/ /* get total heating */ PhotoHeat = 0.; /* sum up various contributors to heating and cooling from the heavy elements */ CoolHeavy.colmet = 0.; secmet = 0.; smetla = 0.; /* this loop includes hydrogenic, which is special case due to presence * of substantial excited state populations */ for( ipZ=0; ipZ>chng 99 may 8, had not included rec cool to ground */ /* correct for recombination cooling * this is surprisingly important in lte.in - with off * about 1000K too hot, with on about 3000K too cool. * also run lte.in with he/h=1, */ /*HeatingCom.heating[ipZ][ipZ] -= (float)( HydroRecCool(1,ipZ)*hydro.hrec[ipZ][0][ipRecNetEsc]* phycon.eden*xIonFracs[ipZ+2][ipZ]);*/ /* add to total heating */ PhotoHeat += HeatingCom.heating[ipZ][ipZ]; } /* do not include hdyrogenic in following loop */ limit = MIN2( IonRange.IonHigh[ipZ]-1 , ipZ ); /* >>chng 99 may 8, this had included hydrogenic, but without correction * for fraction in excited states, so overestimated heating in dense * models near lte */ for( nion=IonRange.IonLow[ipZ]-1; nion < limit; nion++ ) /*for( nion=IonRange.IonLow[ipZ]-1; nion < (IonRange.IonHigh[ipZ]-1); nion++ )*/ { /* this will be heating for a single stage of ionization */ HeatingLo = 0.; HeatingHi = 0.; for( ns=0; ns < nsShellsCom.nsShells[nion][ipZ]; ns++ ) { /* heating by various sub-shells */ HeatingLo += PhotRate.PhotoRate[1][ns][nion][ipZ]; HeatingHi += PhotRate.PhotoRate[2][ns][nion][ipZ]; } /* total heating, all shells, for this stage of ionization */ HeatingCom.heating[nion][ipZ] = xIonFracs[nion+1][ipZ]* (HeatingLo + HeatingHi*Secondaries.heatef); /* add to total heating */ PhotoHeat += HeatingCom.heating[nion][ipZ]; /* HeavIonCool is cooling due to collisional ionization of heavy elements * CollidRate[1][nion][nelem] is cooling, erg/s/atom, eval in CollidIonize */ CoolHeavy.colmet += xIonFracs[nion+1][ipZ]*CollIonRate.CollidRate[1][nion][ipZ]; /* secondary ionization rate, */ secmet += HeatingHi*Secondaries.efionz* xIonFracs[nion+1][ipZ]; /* LA excit rate, =0 if ionized */ smetla += HeatingHi*Secondaries.exctef* xIonFracs[nion+1][ipZ]; } } } /* convert secmet to proper final units */ secmet /= phycon.CollidDensity; smetla /= phycon.CollidDensity; /* find heating due to charge transfer */ HeatingCom.heating[24][0] = SumCTHeat(); /* heating due to pair production */ HeatingCom.heating[21][0] = PhotRate.PhotoRate[2][0][10][0]*Secondaries.heatef*(phycon.hden + 4.F*xIonFracs[0][1]); /* last term above is number of nuclei in helium */ /* this is heating due to fast neutrons, assumed to secondary ionize */ HeatingCom.heating[20][0] = PhotRate.PhotoRate[2][0][11][0]*Secondaries.heatef*phycon.hden; /* turbulent heating, assumed to be a non-ionizing heat agent, added here */ HeatingCom.heating[20][0] += PhotRate.PhotoRate[1][0][11][0]; /* cosmic ray heating, counts as non-ionizing heat since already result of cascade, * this was per secondary ionization, so must multiply by number of sec that * occur per primary */ HeatingCom.heating[20][0] += PhotRate.PhotoRate[1][0][12][0]* (xIonFracs[1][0] + xIonFracs[1][1] + hmi.htwo)*Secondaries.sec2prim; /* bound compton heating */ HeatingCom.heating[18][0] = PhotRate.PhotoRate[2][0][13][0]*Secondaries.heatef*xIonFracs[1][0] + PhotRate.PhotoRate[2][0][14][0]*Secondaries.heatef*xIonFracs[1][1]; /* now sum up all heating agents not included in sum over normal bound levels above */ OtherHeat = 0.; for( ipZ=0; ipZ>>chng 99 may 08, following loop had started at ipZ+3, and so missed [1][0], * which is excited states of hydrogenic species. increase this limit */ for( i=ipZ+1; i < LIMELM; i++ ) { OtherHeat += HeatingCom.heating[i][ipZ]; } } heat.htot = OtherHeat + PhotoHeat; /* check whether Lya heating is important * if so then only very small temp changes can be made */ if( HeatingCom.heating[10][0]/heat.htot > 0.15 ) { HColRec.lgHCoolng = TRUE; } else { HColRec.lgHCoolng = FALSE; } /* add on line heating to this array, heatl was evaluated in sumcool * not zero, because of roundoff error */ if( tcool.heatl/tcool.ctot < -1e-15 ) { fprintf( ioQQQ, " SumHeat gets negative line heating,%10.2e%10.2e this is insane.\n", tcool.heatl, tcool.ctot ); fprintf( ioQQQ, " this is zone%4ld\n", ZoneCnt.nzone ); ShowMe(); puts( "[Stop in sumheat]" ); exit(1); } /******************************************************************* * * secondary ionization and excitation rates * *******************************************************************/ /* secondary rates due to hydrogen */ Secondaries.HyIon = (float)( PhotRate.PhotoRate[2][0][0][0]*Secondaries.efionz*xIonFracs[1][0]/phycon.CollidDensity); Secondaries.HyExc =(float)( PhotRate.PhotoRate[2][0][0][0]*Secondaries.exctef*1.3333F*xIonFracs[1][0]/phycon.CollidDensity); /* secondary ioinzation due to helium */ Secondaries.seche =(float)( (PhotRate.PhotoRate[2][0][0][1]*Secondaries.efionz*xIonFracs[1][1] + PhotRate.PhotoRate[2][0][1][1]*Secondaries.efionz*xIonFracs[2][1])/phycon.CollidDensity); /* total secondary excitation of Lya due to helium, in all following, fac of 1.33 is * asymptotic limit for all hydrogen lines above Lya */ Secondaries.helax =(float)( (PhotRate.PhotoRate[2][0][0][1]*Secondaries.exctef*1.3333F*xIonFracs[1][1] + PhotRate.PhotoRate[2][0][1][1]*Secondaries.exctef*1.3333F*xIonFracs[2][1])/ phycon.CollidDensity); /* the terms Secondaries.seccmp & Secondaries.scmpla contain energies added in highen. * The only significant one is usually bound compton heating */ /* add on heating due to pair production, ionizatoin then excitation */ Secondaries.seccmp = (float)( PhotRate.PhotoRate[2][0][10][0]*Secondaries.efionz* (phycon.hden + 4.F*xIonFracs[0][1])/phycon.CollidDensity); Secondaries.scmpla = (float)( PhotRate.PhotoRate[2][0][10][0]*Secondaries.exctef*1.3333F* (phycon.hden + 4.F*xIonFracs[0][1])/phycon.CollidDensity); /* this is heating due to fast neutrons, assumed to secondary ionize */ Secondaries.seccmp += (float)( PhotRate.PhotoRate[2][0][11][0]*Secondaries.efionz*phycon.hden/phycon.CollidDensity); Secondaries.scmpla += (float)( PhotRate.PhotoRate[2][0][11][0]*Secondaries.exctef*1.3333F*phycon.hden/ phycon.CollidDensity); /* cosmic ray ionization */ /* >>>chng 99 apr 29, term in PhotoRate was not present */ Secondaries.seccmp += (float)( /* this term is cosmic ray ionization, set in highen, did not multiply by * collider density in highen, so do not divide by it here */ /* >>chng 99 jun 29, added efionz to cosmic ray rate, also multiply * by number of secondaries per primary*/ PhotRate.PhotoRate[0][0][12][0]*Secondaries.sec2prim + /* this is the cosmic ray heating rate */ PhotRate.PhotoRate[2][0][12][0]*Secondaries.efionz); /* cosmic ray Lya excitation rate */ Secondaries.scmpla += (float)( PhotRate.PhotoRate[2][0][12][0]*Secondaries.exctef*1.3333f* /* multiply by atomic H and He densities */ (xIonFracs[1][0] + xIonFracs[1][1])/phycon.CollidDensity); /* bound compton heating */ Secondaries.seccmp += (float)((PhotRate.PhotoRate[2][0][13][0]*Secondaries.efionz*xIonFracs[1][0] + PhotRate.PhotoRate[2][0][14][0]*Secondaries.efionz*xIonFracs[1][1])/ phycon.CollidDensity); Secondaries.scmpla += (float)((PhotRate.PhotoRate[2][0][13][0]*Secondaries.exctef*1.3333F*xIonFracs[1][0] + PhotRate.PhotoRate[2][0][14][0]*Secondaries.exctef*1.3333F*xIonFracs[1][1])/ phycon.CollidDensity); /* find total suprathermal collisional ionization rate * CSUPRA is col ionize rate from all processes (inv sec) * X12tot is LA excit rate, inv sec * SECCMP evaluated in HIGHEN, =ioniz rate for cosmic rays, sec bound */ /* option to force secondary ionization rate, normally false */ if( Secondaries.lgCSetOn ) { Secondaries.csupra = Secondaries.SetCsupra; Secondaries.x12tot = Secondaries.SetCsupra; } else { /* CollidDensity is total neutral particle density */ Secondaries.csupra = Secondaries.seccmp + Secondaries.HyIon + Secondaries.seche + (float)secmet ; /* now do same for Ly-alpha */ Secondaries.x12tot = Secondaries.scmpla + Secondaries.HyExc + Secondaries.helax + (float)smetla; } if( trace.lgTrace && Secondaries.csupra > 0. ) { fprintf( ioQQQ, " SumHeat retrn CSUPRA %9.2e SECCMP %6.3f SecHI %6.3f SECHE %6.3f SECMET %6.3f efrac= %9.2e \n", Secondaries.csupra, Secondaries.seccmp/Secondaries.csupra, Secondaries.HyIon/Secondaries.csupra, Secondaries.seche/Secondaries.csupra, secmet/Secondaries.csupra , phycon.ElecFrac ); } /******************************************************************* * * get derivative of total heating * *******************************************************************/ /* now get derivative of heating due to photoionization, * >>chng 96 jan 14 *>>>>>NB heating going down with iincreasing temp is negative */ heat.dHTotDT = -0.7*PhotoHeat/phycon.te; /* oldfac was factor used in old implimentation * any term using it should be rethought */ oldfac = -0.7/phycon.te; /* carbon monoxide */ heat.dHTotDT += HeatingCom.heating[9][0]*oldfac; /* ly alpha collisoinal heating */ heat.dHTotDT += HeatingCom.heating[10][0]*oldfac; /* line heating */ heat.dHTotDT += HeatingCom.heating[22][0]*oldfac; /* free free heating, propto T^-0.5 * >>chng 96 oct 30, from heating(1,12) to FreeFreeHeat, * do cooling separately assume FreeFreeHeat goes as T^-3/2 * dHTotDT = dHTotDT + heating(1,12) * (-0.5/te) */ heat.dHTotDT += hydro.FreeFreeHeat*(-1.5/phycon.te); /* grain heating by photo * >>chng 97 july 2, change deriv in following, better fit to what came out */ heat.dHTotDT += HeatingCom.heating[13][0]*0.1/phycon.te; /* grain heating by collsional */ heat.dHTotDT += HeatingCom.heating[14][0]*oldfac; /* helium triplets heating */ heat.dHTotDT += HeatingCom.heating[2][1]*oldfac; /* hydrogen molecule dissociation */ heat.dHTotDT += HeatingCom.heating[17][0]*oldfac; /* >>chng 96 nov 15, had been oldfac below, wrong sign * H- heating, goes up with temp since rad assoc does */ heat.dHTotDT += HeatingCom.heating[15][0]*0.7/phycon.te; /* H2+ heating */ heat.dHTotDT += HeatingCom.heating[16][0]*oldfac; /* Balmer continuum and all other excited states * - T goes up so does pop and heating * but this all screwed up by change in eden */ heat.dHTotDT += HeatingCom.heating[1][0]*oldfac; /* all three body heating, opposite of coll ion cooling */ heat.dHTotDT += HeatingCom.heating[3][0]*oldfac; /* bound electron recoil heating */ heat.dHTotDT += HeatingCom.heating[18][0]*oldfac; /* compton heating */ heat.dHTotDT += HeatingCom.heating[19][0]*oldfac; /* extra heating source, usually zero */ heat.dHTotDT += HeatingCom.heating[20][0]*oldfac; /* pair production */ heat.dHTotDT += HeatingCom.heating[21][0]*oldfac; /* deriv of heating due to collisions of H lines, heating(1,24) */ for( ipZ=0; ipZ 0. ) { heat.dHTotDT += holod.dfcool; } /* possibility of getting emperical heating derivative * normally false, set true with 'set numerical derivatives' command */ if( NumDeriv.lgNumDeriv ) { if( ((ZoneCnt.nzone > 2 && ZoneCnt.nzone == nzSave) && oldtemp != phycon.te) && nhit > 4 ) { /* hnit is number of trys on this zone - use to stop numerical problems * do not evaluate numerical deriv until well into solution */ deriv = (oldheat - heat.htot)/(oldtemp - phycon.te); heat.dHTotDT = deriv; } else { deriv = heat.dHTotDT; } /* this happens when new zone starts */ if( ZoneCnt.nzone != nzSave ) { nhit = 0; } nzSave = ZoneCnt.nzone; nhit += 1; oldheat = heat.htot; oldtemp = phycon.te; } if( trace.lgTrace && trace.lgHeatBug ) { ipnt = 0; /* this loops through the 2D array that contains all agents counted in htot */ for( i=0; i < LIMELM; i++ ) { for( j=0; j < LIMELM; j++ ) { if( HeatingCom.heating[j][i]/heat.htot > FAINT_HEAT ) { ipsave[ipnt] = i; jpsave[ipnt] = j; save[ipnt] = HeatingCom.heating[j][i]/heat.htot; /* increment ipnt, but do not let it get too big */ ipnt = MIN2(NDIM-1,ipnt+1); } } } /* now check for possible line heating not counted in 1,23 * there should not be any significant heating source heree * since they would not be counted in deriv correctly */ for( i=0; i < Coolants.ncltot; i++ ) { if( Coolants.heatnt[i]/heat.htot > FAINT_HEAT ) { ipsave[ipnt] = -1; jpsave[ipnt] = i; save[ipnt] = Coolants.heatnt[i]/heat.htot; ipnt = MIN2(NDIM-1,ipnt+1); } } fprintf( ioQQQ, " SumHeat HTOT%11.3e Te:%10.3e dH/dT%10.2e other %10.2e photo %10.2e\n", heat.htot, phycon.te, heat.dHTotDT , OtherHeat , PhotoHeat); fprintf( ioQQQ, " " ); for( i=0; i < ipnt; i++ ) { /*lint -e644 these three are initialized above */ fprintf( ioQQQ, "%3ld%3ld%6.3f", jpsave[i], ipsave[i], save[i] ); /*lint +e644 these three are initialized above */ /* throw a cr every n numbers */ if( !(i%8) && i>0 && i!=(ipnt-1) ) { fprintf( ioQQQ, "\n " ); } } fprintf( ioQQQ, "\n" ); } # ifdef DEBUG_FUN fputs( " <->SumHeat()\n", debug_fp ); # endif return; } /* =============================================================================*/ /* ZeroHeat is called by ionize */ void ZeroHeat( void ) { long int i , j; for( i=0; i < LIMELM; i++ ) { for( j=0; j < LIMELM; j++ ) { HeatingCom.heating[j][i] = 0.; } } # ifdef DEBUG_FUN fputs( " <->ZeroHeat()\n", debug_fp ); # endif return; }