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