00001
00002
00003
00004
00005
00006 #include "cddefines.h"
00007 #include "hmi.h"
00008 #include "trace.h"
00009 #include "grainvar.h"
00010 #include "rfield.h"
00011 #include "mole.h"
00012 #include "dense.h"
00013 #include "conv.h"
00014
00015
00016
00017
00018 int eden_sum(void)
00019 {
00020 long int i,
00021 ion,
00022 nelem;
00023
00024 double sum_all_ions ,
00025 sum_metals ,
00026 hmole_eden,
00027 eden_ions[LIMELM];
00028
00029 DEBUG_ENTRY( "eden_sum()" );
00030
00031
00032 dense.EdenTrue = dense.EdenExtra;
00033
00034
00035 sum_all_ions = 0.;
00036 sum_metals = 0.;
00037 for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
00038 {
00039 if( nelem==ipLITHIUM )
00040 sum_metals = 0.;
00041 eden_ions[nelem] = dense.xIonDense[nelem][1];
00042
00043
00044
00045
00046 for( ion=2; ion <= nelem+1; ion++ )
00047 {
00048
00049 eden_ions[nelem] += ion*dense.xIonDense[nelem][ion];
00050 }
00051 sum_all_ions += eden_ions[nelem];
00052 sum_metals += eden_ions[nelem];
00053 }
00054 dense.EdenTrue += sum_all_ions;
00055 ASSERT( dense.EdenTrue >= 0. );
00056
00057
00058 co.comole_eden = 0.;
00059 for( i=0; i < mole.num_comole_calc; i++ )
00060 {
00061 if(COmole[i]->n_nuclei != 1)
00062 {
00063 co.comole_eden += COmole[i]->hevmol*COmole[i]->nElec;
00064 }
00065 }
00066
00067
00068 dense.EdenTrue += co.comole_eden;
00069 ASSERT( dense.EdenTrue >= 0. );
00070
00071
00072 hmole_eden = 0.;
00073 for(i=0;i<N_H_MOLEC;i++)
00074 {
00075
00076 hmole_eden += hmi.Hmolec[i]*hmi.nElectron[i];
00077 }
00078
00079
00080
00081
00082
00083 if( (-hmole_eden) > dense.EdenTrue/4. && conv.lgSearch )
00084 {
00085
00086 dense.EdenTrue *= 0.9;
00087 }
00088 else if( (-hmole_eden) > dense.EdenTrue )
00089 {
00090
00091
00092
00093 fprintf(ioQQQ," PROBLEM sum eden from hmole too neg, set to limt. EdenTrue:%.3e hmole_eden:%.3e \n",
00094 dense.EdenTrue,
00095 hmole_eden);
00096 dense.EdenTrue = dense.EdenTrue/2.;
00097 }
00098 else
00099 {
00100
00101 dense.EdenTrue += hmole_eden;
00102 }
00103 ASSERT(dense.EdenTrue >= 0. );
00104
00105
00106
00107 if( dense.EdenSet > 0. )
00108 {
00109 dense.EdenTrue = dense.EdenSet;
00110 dense.eden_from_metals = 1.;
00111 }
00112 else
00113 {
00114
00115
00116 dense.eden_from_metals = sum_metals / SDIV(dense.EdenTrue);
00117
00118 }
00119
00120
00121
00122
00123
00124
00125
00126
00127 if( dense.EdenTrue+gv.TotalEden*gv.lgGrainElectrons >= 0. )
00128 {
00129
00130
00131 dense.EdenTrue += gv.TotalEden*gv.lgGrainElectrons;
00132 }
00133 else
00134 {
00135
00136 if( !conv.lgSearch )
00137 fprintf(ioQQQ,
00138 " PROBLEM eden grain neg limt %.2f eden %.4e new eden bef grn %.4e"
00139 "grain eden %.4e set edentrue to %.4e Search?%c\n",
00140 fnzone,
00141 dense.eden,
00142 dense.EdenTrue,
00143 gv.TotalEden,
00144 dense.eden*0.9,
00145 TorF(conv.lgSearch));
00146
00147
00148 dense.EdenTrue += gv.TotalEden*gv.lgGrainElectrons;
00149
00150
00151
00152 }
00153
00154
00155
00156
00157 if( trace.lgTrace || trace.lgESOURCE )
00158 {
00159 fprintf( ioQQQ,
00160 " eden_sum zn:%.2f current:%.4e new true:%.4e ions:%.4e comole:%.4e hmole:%.4e grain:%.4e extra:%.4e sum:%.4e LaOTS:%.4e\n",
00161 fnzone ,
00162 dense.eden ,
00163 dense.EdenTrue ,
00164 sum_all_ions ,
00165 co.comole_eden ,
00166 hmole_eden ,
00167 gv.TotalEden*gv.lgGrainElectrons,
00168 dense.EdenExtra ,
00169 sum_all_ions + co.comole_eden + hmole_eden + gv.TotalEden*gv.lgGrainElectrons + dense.EdenExtra,
00170 rfield.otslin[2182] );
00171
00172 if( 1 || trace.lgNeBug )
00173 {
00174
00175 for(nelem=ipHYDROGEN; nelem < LIMELM; nelem++)
00176 {
00177 if( nelem==0 )
00178 {
00179 fprintf( ioQQQ, " eden_sum H -Ne:" );
00180 }
00181 else if( nelem==10 )
00182 {
00183 fprintf( ioQQQ, "\n eden_sum Na-Ca:" );
00184 }
00185 else if( nelem==20 )
00186 {
00187 fprintf( ioQQQ, "\n eden_sum Sc-Zn:" );
00188 }
00189 fprintf( ioQQQ, " %.4e", eden_ions[nelem] );
00190 if( nelem==29 )
00191 {
00192 fprintf( ioQQQ, "\n" );
00193 }
00194
00195
00196 }
00197
00198 }
00199 }
00200
00201
00202
00203 if( 0 && dense.EdenTrue < 0. )
00204 {
00205 fprintf( ioQQQ, "eden_sum finds non-positive electron density.\n" );
00206 fprintf( ioQQQ,
00207 " eden_sum: EdenTrue to%12.4e fm%12.4e Ne(H):%10.2e Ne(He):%10.2e Ne(C):\n",
00208 dense.EdenTrue,
00209 dense.eden,
00210 dense.xIonDense[ipHYDROGEN][1],
00211 dense.xIonDense[ipHELIUM][1] + 2.*dense.xIonDense[ipHELIUM][2] );
00212 ShowMe();
00213 puts( "[Stop in eden_sum]" );
00214 cdEXIT(EXIT_FAILURE);
00215 }
00216 else if( dense.EdenTrue == 0. )
00217 {
00218 fprintf( ioQQQ, "\nDISASTER PROBLEM eden_sum finds an electron density of zero, this is unphysical.\n" );
00219 fprintf( ioQQQ, "There appears to be no source of ionization.\n");
00220 fprintf( ioQQQ, "Consider adding background cosmic rays with COSMIC RAY BACKGROUND and BACKGROUND commands.\n\n");
00221 lgAbort = true;
00222
00223 DEBUG_EXIT( "eden_sum()" );
00224
00225 return 1;
00226 }
00227
00228 {
00229
00230 enum {DEBUG_LOC=false};
00231
00232 if( DEBUG_LOC )
00233 {
00234 fprintf(ioQQQ,"esumdebugg\t%li\t%.2e\t%.2e\n",
00235 nzone,
00236 dense.eden, dense.EdenTrue);
00237 }
00238 }
00239
00240
00241 dense.eden = MAX2( SMALLFLOAT , dense.eden );
00242
00243 DEBUG_EXIT( "eden_sum()" );
00244
00245 return 0;
00246 }
00247