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