00001
00002
00003
00004 #include "cddefines.h"
00005 #include "physconst.h"
00006 #include "trace.h"
00007 #include "heavy.h"
00008 #include "radius.h"
00009 #include "magnetic.h"
00010 #include "hextra.h"
00011 #include "thermal.h"
00012 #include "dense.h"
00013 #include "doppvel.h"
00014 #include "ionbal.h"
00015 #include "rfield.h"
00016 #include "opacity.h"
00017 #include "gammas.h"
00018 #include "highen.h"
00019 #include "pressure.h"
00020
00021 void highen( void )
00022 {
00023 long int i,
00024 ion,
00025 nelem;
00026
00027 double CosRayDen,
00028 crsphi,
00029 heatin,
00030 sqthot;
00031
00032 double hin;
00033
00034 DEBUG_ENTRY( "highen()" );
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049 if( rfield.comoff > .0 && ionbal.lgCompRecoil )
00050 {
00051 for( nelem=0; nelem<LIMELM; ++nelem )
00052 {
00053 for( ion=0; ion<nelem+1; ++ion )
00054 {
00055 ionbal.CompRecoilIonRate[nelem][ion] = 0.;
00056 ionbal.CompRecoilHeatRate[nelem][ion] = 0.;
00057 if( dense.xIonDense[nelem][ion] > 0. )
00058 {
00059
00060
00061
00062
00063 for( i=ionbal.ipCompRecoil[nelem][ion]; i < rfield.nflux; ++i)
00064 {
00065 double recoil_energy;
00066 crsphi = opac.OpacStack[i-1+opac.iopcom] * rfield.SummedCon[i] *
00067
00068 ionbal.nCompRecoilElec[nelem-ion];
00069
00070
00071
00072
00073 ionbal.CompRecoilIonRate[nelem][ion] += crsphi;
00074
00075
00076
00077
00078
00079 recoil_energy = rfield.anu2[i] /
00080 ( EN1RYD * ELECTRON_MASS * SPEEDLIGHT * SPEEDLIGHT) * EN1RYD* EN1RYD -
00081 Heavy.Valence_IP_Ryd[nelem][ion];
00082
00083
00084 ionbal.CompRecoilHeatRate[nelem][ion] += crsphi*recoil_energy;
00085 }
00086
00087 ionbal.CompRecoilHeatRate[nelem][ion] *= EN1RYD;
00088
00089 ASSERT( ionbal.CompRecoilHeatRate[nelem][ion] >= 0.);
00090 ASSERT( ionbal.CompRecoilIonRate[nelem][ion] >= 0.);
00091 }
00092 }
00093 }
00094 }
00095 else
00096 {
00097 for( nelem=0; nelem<LIMELM; ++nelem )
00098 {
00099 for( ion=0; ion<nelem+1; ++ion )
00100 {
00101 ionbal.CompRecoilIonRate[nelem][ion] = 0.;
00102 ionbal.CompRecoilHeatRate[nelem][ion] = 0.;
00103 }
00104 }
00105 }
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119 if( hextra.cryden > 0. )
00120 {
00121 ASSERT( hextra.crtemp > 0. );
00122
00123
00124 if( hextra.lg_CR_B_equipartition )
00125 {
00126
00127
00128
00129 CosRayDen = hextra.background_density *
00130
00131
00132
00133 magnetic.energydensity /
00134 (CR_EDEN_GAL_BACK_EV_CMM3 * EN1EV );
00135 }
00136 else
00137 {
00138
00139 CosRayDen = hextra.cryden*pow(radius.Radius/radius.rinner,(double)hextra.crpowr);
00140 }
00141
00142
00143 hextra.cr_energydensity = CosRayDen/hextra.background_density *
00144 (CR_EDEN_GAL_BACK_EV_CMM3 * EN1EV );
00145
00146
00147 sqthot = sqrt(hextra.crtemp);
00148
00149
00150
00151
00152
00153 ionbal.CosRayHeatThermalElectrons = 5.5e-14/sqthot*CosRayDen;
00154
00155
00156 ionbal.CosRayIonRate = 1.22e-4/sqthot*
00157 log(2.72*pow(hextra.crtemp/1e8,0.097))*CosRayDen;
00158
00159
00160 if( hextra.crtemp > 2e9 )
00161 {
00162
00163
00164 ionbal.CosRayIonRate *= 3.;
00165
00166 }
00167 else
00168 {
00169
00170 ionbal.CosRayIonRate *= 10.;
00171 }
00172
00173
00174
00175
00176
00177
00178
00179
00180 ionbal.CosRayHeatNeutralParticles = ionbal.CosRayIonRate*2.574*EN1RYD;
00181
00182 if( trace.lgTrace )
00183 {
00184 fprintf( ioQQQ, " highen: cosmic ray density;%10.2e CRion rate;%10.2e CR heat rate=;%10.2e CRtemp;%10.2e\n",
00185 CosRayDen, ionbal.CosRayIonRate, ionbal.CosRayHeatNeutralParticles, hextra.crtemp );
00186 }
00187 }
00188 else
00189 {
00190 ionbal.CosRayIonRate = 0.;
00191 ionbal.CosRayHeatNeutralParticles = 0.;
00192 }
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210 if( (hextra.TurbHeat+DoppVel.DispScale) != 0. )
00211 {
00212
00213
00214
00215
00216 if( hextra.lgHextraDepth )
00217 {
00218
00219 ionbal.ExtraHeatRate =
00220 hextra.TurbHeat*sexp(radius.depth /hextra.turrad);
00221
00222
00223 if( hextra.turback != 0. )
00224 {
00225 ionbal.ExtraHeatRate +=
00226 hextra.TurbHeat*sexp((hextra.turback - radius.depth) /hextra.turrad);
00227 }
00228 }
00229 else if( hextra.lgHextraDensity )
00230 {
00231
00232 ionbal.ExtraHeatRate =
00233 hextra.TurbHeat*dense.gas_phase[ipHYDROGEN]/hextra.HextraScaleDensity;
00234 }
00235 else if( hextra.lgHextraSS )
00236 {
00237
00238 ionbal.ExtraHeatRate =
00239 hextra.HextraSSalpha*pressure.PresGasCurr*pow(hextra.HextraSS_M*GRAV_CONST*
00240 pow(hextra.HextraSSradius,-3.0),0.5);
00241 }
00242
00243 else if( DoppVel.DispScale > 0. )
00244 {
00245 double turb = DoppVel.TurbVel * sexp( radius.depth / DoppVel.DispScale );
00246
00247
00248 ionbal.ExtraHeatRate = 3.45e-28 / 2.82824 * turb * turb * turb *
00249 ( dense.gas_phase[ipHYDROGEN] / 1e10 ) / (DoppVel.DispScale/1e13);
00250 }
00251 else
00252 {
00253
00254 ionbal.ExtraHeatRate = hextra.TurbHeat;
00255 }
00256 }
00257
00258 else
00259 {
00260 ionbal.ExtraHeatRate = 0.;
00261 }
00262
00263
00264
00265
00266
00267
00268 if( hextra.lgNeutrnHeatOn )
00269 {
00270
00271
00272
00273
00274 ionbal.xNeutronHeatRate = hextra.totneu*hextra.CrsSecNeutron;
00275 }
00276 else
00277 {
00278 ionbal.xNeutronHeatRate = 0.;
00279 }
00280
00281
00282
00283
00284
00285
00286
00287 t_phoHeat photoHeat;
00288 ionbal.PairProducPhotoRate[0] = GammaK(opac.ippr,rfield.nflux,opac.ioppr,1.,&photoHeat);
00289 ionbal.PairProducPhotoRate[1] = photoHeat.HeatLowEnr;
00290 ionbal.PairProducPhotoRate[2] = photoHeat.HeatHiEnr;
00291
00292
00293
00294
00295
00296
00297 rfield.cmcool = 0.;
00298 rfield.cmheat = 0.;
00299 heatin = 0.;
00300
00301 if( rfield.comoff >0.01 )
00302 {
00303 for( i=0; i < rfield.nflux; i++ )
00304 {
00305
00306
00307
00308
00309 rfield.comup[i] = (double)(rfield.flux[0][i]+rfield.ConInterOut[i]+
00310 rfield.outlin[0][i]+ rfield.outlin_noplot[i])*rfield.csigc[i]*(dense.eden*4.e0*
00311 6.338e-6*1e-15);
00312
00313 rfield.cmcool += rfield.comup[i];
00314
00315
00316
00317
00318 rfield.comdn[i] = (double)(rfield.flux[0][i]+rfield.ConInterOut[i]+
00319 rfield.outlin[0][i]+ rfield.outlin_noplot[i])*rfield.csigh[i]*dense.eden*1e-15;
00320
00321
00322 hin = (double)(rfield.flux[0][i]+rfield.ConInterOut[i]+rfield.outlin[0][i]+rfield.outlin_noplot[i])*
00323 rfield.csigh[i]*rfield.OccNumbIncidCont[i]*dense.eden*1e-15;
00324 rfield.comdn[i] += hin;
00325 heatin += hin;
00326
00327
00328 rfield.cmheat += rfield.comdn[i];
00329 }
00330
00331
00332 if( rfield.cmheat > 0. )
00333 {
00334 rfield.cinrat = MAX2(rfield.cinrat,heatin/rfield.cmheat);
00335 }
00336
00337 if( trace.lgTrace && trace.lgComBug )
00338 {
00339 fprintf( ioQQQ, " HIGHEN: COOL num=%8.2e HEAT num=%8.2e\n",
00340 rfield.cmcool, rfield.cmheat );
00341 }
00342 }
00343
00344
00345
00346 thermal.heating[0][19] = rfield.cmheat;
00347
00348 if( trace.lgTrace && trace.lgComBug )
00349 {
00350 fprintf( ioQQQ,
00351 " HIGHEN finds heating fracs= frac(compt)=%10.2e "
00352 " f(pair)%10.2e totHeat=%10.2e\n",
00353 thermal.heating[0][19]/thermal.htot,
00354 thermal.heating[0][21]/thermal.htot,
00355 thermal.htot );
00356 }
00357 return;
00358 }