00001
00002
00003
00004
00005
00006 #include "cddefines.h"
00007 #include "physconst.h"
00008 #include "radius.h"
00009 #include "dense.h"
00010 #include "hyperfine.h"
00011 #include "magnetic.h"
00012 #include "hmi.h"
00013 #include "phycon.h"
00014 #include "geometry.h"
00015 #include "mean.h"
00016
00017 t_mean mean;
00018
00019 t_mean::t_mean()
00020 {
00021 DEBUG_ENTRY( "t_mean()" );
00022
00023
00024 mean.xIonMean.reserve(3);
00025 for( int j=0; j < 3; ++j )
00026 {
00027
00028 mean.xIonMean.reserve(j,LIMELM);
00029 for( int nelem=0; nelem < LIMELM; ++nelem )
00030 {
00031 int limit = max(3,nelem+2);
00032
00033 mean.xIonMean.reserve(j,nelem,limit);
00034 for( int ion=0; ion < limit; ++ion )
00035
00036 mean.xIonMean.reserve(j,nelem,ion,2);
00037 }
00038 }
00039 mean.xIonMean.alloc();
00040 mean.xIonEdenMean.alloc( mean.xIonMean.clone() );
00041 mean.TempIonMean.alloc( mean.xIonMean.clone() );
00042 mean.TempIonEdenMean.alloc( mean.xIonMean.clone() );
00043 mean.TempB_HarMean.alloc(3,2);
00044 mean.TempHarMean.alloc(3,2);
00045 mean.TempH_21cmSpinMean.alloc(3,2);
00046 mean.TempMean.alloc(3,2);
00047 mean.TempEdenMean.alloc(3,2);
00048 }
00049
00050
00051 void t_mean::MeanZero()
00052 {
00053 DEBUG_ENTRY( "t_mean::MeanZero()" );
00054
00055
00056
00057
00058 mean.xIonMean.zero();
00059 mean.xIonEdenMean.zero();
00060 mean.TempIonMean.zero();
00061 mean.TempIonEdenMean.zero();
00062 mean.TempB_HarMean.zero();
00063 mean.TempHarMean.zero();
00064 mean.TempH_21cmSpinMean.zero();
00065 mean.TempMean.zero();
00066 mean.TempEdenMean.zero();
00067
00068 return;
00069 }
00070
00071
00072 void t_mean::MeanInc()
00073 {
00074
00075
00076
00077 DEBUG_ENTRY( "t_mean::MeanInc()" );
00078
00079
00080 double Jac[3] = { radius.drad_x_fillfac, radius.darea_x_fillfac, radius.dVeffVol };
00081
00082 for( int d=0; d < 3; ++d )
00083 {
00084 double dc;
00085 for( int nelem=0; nelem < LIMELM; nelem++ )
00086 {
00087 int limit = max(3,nelem+2);
00088
00089
00090 double norm = dense.gas_phase[nelem]*Jac[d];
00091 for( int ion=0; ion < limit; ion++ )
00092 {
00093 if( nelem == ipHYDROGEN && ion == 2 )
00094
00095
00096
00097 dc = 2.*hmi.H2_total*Jac[d];
00098 else
00099 dc = dense.xIonDense[nelem][ion]*Jac[d];
00100 mean.xIonMean[d][nelem][ion][0] += dc;
00101 mean.xIonMean[d][nelem][ion][1] += norm;
00102 mean.TempIonMean[d][nelem][ion][0] += dc*phycon.te;
00103 mean.TempIonMean[d][nelem][ion][1] += dc;
00104
00105 dc *= dense.eden;
00106 mean.xIonEdenMean[d][nelem][ion][0] += dc;
00107 mean.xIonEdenMean[d][nelem][ion][1] += norm*dense.eden;
00108 mean.TempIonEdenMean[d][nelem][ion][0] += dc*phycon.te;
00109 mean.TempIonEdenMean[d][nelem][ion][1] += dc;
00110 }
00111 }
00112
00113
00114
00115
00116
00117
00118
00119
00120 dc = ( hyperfine.Tspin21cm > SMALLFLOAT ) ? dense.xIonDense[ipHYDROGEN][0]*Jac[d]/phycon.te : 0.;
00121
00122 mean.TempB_HarMean[d][0] += sqrt(fabs(magnetic.pressure)*PI8) * dc;
00123
00124 mean.TempB_HarMean[d][1] += dc;
00125
00126
00127
00128 dc = dense.xIonDense[ipHYDROGEN][0]*Jac[d];
00129
00130
00131
00132
00133 mean.TempHarMean[d][0] += dc;
00134 mean.TempHarMean[d][1] += dc/phycon.te;
00135
00136
00137 mean.TempH_21cmSpinMean[d][0] += dc;
00138 mean.TempH_21cmSpinMean[d][1] += dc / SDIV( hyperfine.Tspin21cm );
00139
00140 dc = Jac[d];
00141 mean.TempMean[d][0] += dc*phycon.te;
00142 mean.TempMean[d][1] += dc;
00143
00144 dc = Jac[d]*dense.eden;
00145 mean.TempEdenMean[d][0] += dc*phycon.te;
00146 mean.TempEdenMean[d][1] += dc;
00147 }
00148 return;
00149 }
00150
00151
00152 void t_mean::MeanIon(
00153
00154 char chType,
00155
00156 long int nelem,
00157
00158 long int dim,
00159
00160 long int *n,
00161
00162
00163 realnum arlog[],
00164
00165 bool lgDensity ) const
00166 {
00167 DEBUG_ENTRY( "t_mean::MeanIon()" );
00168
00169
00170 int limit = max( 3, nelem+2 );
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180 if( !dense.lgElmtOn[nelem] )
00181 {
00182
00183 for( int ion=0; ion < limit; ion++ )
00184 arlog[ion] = -30.f;
00185 *n = 0;
00186 return;
00187 }
00188
00189
00190 *n = limit;
00191 while( *n > 0 && mean.xIonMean[0][nelem][*n-1][0] <= 0. )
00192 {
00193 arlog[*n-1] = -30.f;
00194 --*n;
00195 }
00196
00197 double meanv, normv;
00198 for( int ion=0; ion < *n; ion++ )
00199 {
00200
00201 if( chType == 'i' )
00202 {
00203 if( lgDensity )
00204 {
00205 meanv = mean.xIonEdenMean[dim][nelem][ion][0];
00206 normv = mean.xIonEdenMean[dim][nelem][ion][1];
00207 }
00208 else
00209 {
00210 meanv = mean.xIonMean[dim][nelem][ion][0];
00211 normv = mean.xIonMean[dim][nelem][ion][1];
00212 }
00213 arlog[ion] = ( meanv > 0. ) ? (realnum)log10(max(1e-30,meanv/normv)) : -30.f;
00214 }
00215
00216 else if( chType == 't' )
00217 {
00218 if( lgDensity )
00219 {
00220 meanv = mean.TempIonEdenMean[dim][nelem][ion][0];
00221 normv = mean.TempIonEdenMean[dim][nelem][ion][1];
00222 }
00223 else
00224 {
00225 meanv = mean.TempIonMean[dim][nelem][ion][0];
00226 normv = mean.TempIonMean[dim][nelem][ion][1];
00227 }
00228 arlog[ion] = ( normv > SMALLFLOAT ) ? (realnum)log10(max(1e-30,meanv/normv)) : -30.f;
00229 }
00230 else
00231 {
00232 fprintf(ioQQQ," MeanIon called with insane job: %c \n",chType);
00233 TotalInsanity();
00234 }
00235 }
00236 return;
00237 }