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