00001
00002
00003
00004
00005 #include "cddefines.h"
00006 #include "elementnames.h"
00007 #include "radius.h"
00008 #include "heavy.h"
00009 #include "conv.h"
00010 #include "rfield.h"
00011 #include "phycon.h"
00012 #include "mole.h"
00013 #include "thermal.h"
00014 #include "iso.h"
00015 #include "dense.h"
00016 #include "conv.h"
00017 #include "struc.h"
00018 #include "ionbal.h"
00019
00020 void ion_trim(
00021
00022 long int nelem )
00023 {
00024
00025
00026 bool lgHi_Down = false;
00027 bool lgHi_Up = false;
00028 bool lgHi_Up_enable;
00029
00030 bool lgLo_Up = false;
00031 bool lgLo_Down = false;
00032 long int ion_lo_save = dense.IonLow[nelem],
00033 ion_hi_save = dense.IonHigh[nelem];
00034 long int ion;
00035 realnum trimhi , trimlo;
00036
00037
00038
00039
00040
00041
00042
00043 DEBUG_ENTRY( "ion_trim()" );
00044
00045
00046
00047 ASSERT( conv.nTotalIoniz>2 );
00048
00049
00050 ASSERT( nelem >= ipHELIUM && nelem < LIMELM );
00051 ASSERT( dense.IonLow[nelem] >= 0 );
00052 ASSERT( dense.IonHigh[nelem] <= nelem+1 );
00053
00054
00055 ASSERT( dense.IonLow[nelem] < dense.IonHigh[nelem] ||
00056 (dense.IonLow[nelem]==0 && dense.IonHigh[nelem]==0 ) ||
00057 dense.lgSetIoniz[nelem] );
00058
00059
00060
00061
00062 if( conv.lgSearch )
00063 {
00064 trimhi = (realnum)ionbal.trimhi * 1e-4f;
00065 trimlo = (realnum)ionbal.trimlo * 1e-4f;
00066 }
00067 else
00068 {
00069 trimhi = (realnum)ionbal.trimhi;
00070 trimlo = (realnum)ionbal.trimlo;
00071 }
00072
00073
00074
00075 if( nelem == ipHELIUM )
00076 {
00077
00078 trimlo = SMALLFLOAT;
00079
00080
00081
00082 if( dense.IonHigh[ipHELIUM] == 1 )
00083 {
00084 trimhi = MIN2( trimhi , 1e-20f );
00085 }
00086 else if( dense.IonHigh[ipHELIUM] == 2 )
00087 {
00088 if( conv.lgSearch )
00089 {
00090
00091
00092
00093
00094
00095
00096 trimhi = MIN2( trimhi , 1e-17f );
00097 }
00098 else
00099 {
00100
00101 trimhi = MIN2( trimhi , 1e-12f );
00102 }
00103 }
00104 }
00105
00106
00107
00108 if( mole.lgElem_in_chemistry[nelem] )
00109 {
00110 trimlo = SMALLFLOAT;
00111 if(dense.IonHigh[nelem] ==2)
00112 {
00113 trimhi = MIN2(trimhi, 1e-20f);
00114 }
00115 }
00116
00117
00118
00119
00120
00121
00122
00123 if( conv.lgSearch )
00124 {
00125
00126 while(
00127 (dense.xIonDense[nelem][dense.IonHigh[nelem]]/dense.gas_phase[nelem] < dense.density_low_limit ||
00128 dense.xIonDense[nelem][dense.IonHigh[nelem]] < dense.density_low_limit ) &&
00129
00130 dense.IonHigh[nelem] > dense.IonLow[nelem] )
00131 {
00132
00133
00134 dense.xIonDense[nelem][dense.IonHigh[nelem]] = 0.;
00135 thermal.heating[nelem][dense.IonHigh[nelem]-1] = 0.;
00136 if( dense.IonHigh[nelem] == nelem+1-NISO )
00137 {
00138 long int ipISO = nelem - dense.IonHigh[nelem];
00139 ASSERT( ipISO>=0 && ipISO<NISO );
00140 StatesElem[ipISO][nelem][0].Pop = 0.;
00141 }
00142
00143
00144 --dense.IonHigh[nelem];
00145 ASSERT( dense.IonHigh[nelem] >= 0);
00146
00147 lgHi_Down = true;
00148 }
00149
00150
00151 while(
00152 (dense.xIonDense[nelem][dense.IonLow[nelem]]/dense.gas_phase[nelem] < dense.density_low_limit ||
00153 dense.xIonDense[nelem][dense.IonLow[nelem]] < dense.density_low_limit ) &&
00154 dense.IonLow[nelem] < dense.IonHigh[nelem] - 1 )
00155 {
00156
00157
00158 dense.xIonDense[nelem][dense.IonLow[nelem]] = 0.;
00159
00160
00161 thermal.heating[nelem][dense.IonLow[nelem]] = 0.;
00162 if( dense.IonLow[nelem] == nelem+1-NISO )
00163 {
00164 long int ipISO = nelem - dense.IonLow[nelem];
00165 ASSERT( ipISO>=0 && ipISO<NISO );
00166 StatesElem[ipISO][nelem][0].Pop = 0.;
00167 }
00168
00169
00170 ++dense.IonLow[nelem];
00171 lgLo_Up = true;
00172 }
00173 }
00174
00175
00176
00177 ASSERT( dense.IonLow[nelem] < dense.IonHigh[nelem] ||
00178 (dense.IonLow[nelem]==0 && dense.IonHigh[nelem]==0 ) ||
00179 dense.lgSetIoniz[nelem] );
00180
00181
00182
00183
00184 if( dense.IonHigh[nelem] > 1 &&
00185 (dense.IonLow[nelem]==dense.IonHigh[nelem]-1) &&
00186 (dense.xIonDense[nelem][dense.IonLow[nelem]] < dense.density_low_limit) )
00187 {
00188 fprintf(ioQQQ,
00189 "PROBLEM DISASTER the density of ion %li of element %s is too small to be computed on this cpu.\n",
00190 dense.IonLow[nelem]+1,
00191 elementnames.chElementName[nelem]);
00192 fprintf(ioQQQ,
00193 "Turn off the element with the command ELEMENT XXX OFF or do not consider "
00194 "gas with low density, the current hydrogen density is %.2e cm-3.\n",
00195 dense.gas_phase[ipHYDROGEN]);
00196 fprintf(ioQQQ,
00197 "The outer radius of the cloud is %.2e cm - should a stopping "
00198 "radius be set?\n",
00199 radius.Radius );
00200 fprintf(ioQQQ,
00201 "The abort flag is being set - the calculation is stopping.\n");
00202 lgAbort = true;
00203 return;
00204 }
00205
00206
00207 if(
00208 dense.IonHigh[nelem] > dense.IonLow[nelem] &&
00209 ( (dense.xIonDense[nelem][dense.IonHigh[nelem]]/dense.gas_phase[nelem] <=
00210 trimhi ) ||
00211 (dense.xIonDense[nelem][dense.IonHigh[nelem]] <= dense.density_low_limit ) )
00212 )
00213 {
00214
00215
00216 if( dense.IonHigh[nelem]>1 ||
00217 (dense.IonHigh[nelem]==1&&dense.xIonDense[nelem][1]<100.*dense.density_low_limit) )
00218 {
00219 dense.xIonDense[nelem][dense.IonHigh[nelem]] = 0.;
00220 thermal.heating[nelem][dense.IonHigh[nelem]-1] = 0.;
00221 if( dense.IonHigh[nelem] == nelem+1-NISO )
00222 {
00223 long int ipISO = nelem - dense.IonHigh[nelem];
00224 ASSERT( ipISO>=0 && ipISO<NISO );
00225 StatesElem[ipISO][nelem][0].Pop = 0.;
00226 }
00227 --dense.IonHigh[nelem];
00228 lgHi_Down = true;
00229 }
00230 }
00231
00232
00233 lgHi_Up_enable = true;
00234
00235 if( !ionbal.lgTrimhiOn )
00236 lgHi_Up_enable = false;
00237
00238 if( dense.IonHigh[nelem]>=nelem+1 )
00239 lgHi_Up_enable = false;
00240
00241 if( lgHi_Down )
00242 lgHi_Up_enable = false;
00243
00244 if( dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN]<0.9 )
00245 lgHi_Up_enable = false;
00246
00247 if( lgHi_Up_enable )
00248 {
00249 realnum abundold = 0;
00250 if( nzone>2 )
00251 abundold = struc.xIonDense[nelem][dense.IonHigh[nelem]][nzone-3]/
00252 SDIV( struc.gas_phase[nelem][nzone-3]);
00253 realnum abundnew = dense.xIonDense[nelem][dense.IonHigh[nelem]]/dense.gas_phase[nelem];
00254
00255
00256 if( (Heavy.Valence_IP_Ryd[nelem][dense.IonHigh[nelem]] < rfield.anu[rfield.nflux-1] ||
00257 Heavy.Valence_IP_Ryd[nelem][dense.IonHigh[nelem]] < phycon.te_ryd*10.) &&
00258 dense.xIonDense[nelem][dense.IonHigh[nelem]]/dense.gas_phase[nelem] > 1e-4f &&
00259
00260 abundnew > abundold*1.01 )
00261 {
00262
00263
00264 ++dense.IonHigh[nelem];
00265 lgHi_Up = true;
00266
00267 dense.xIonDense[nelem][dense.IonHigh[nelem]] = (realnum)(dense.density_low_limit*10.);
00268 }
00269 }
00270
00271
00272
00273 ASSERT( dense.IonLow[nelem] < dense.IonHigh[nelem] ||
00274 (dense.IonLow[nelem]==0 && dense.IonHigh[nelem]==0 ) ||
00275 dense.lgSetIoniz[nelem] );
00276
00277
00278 if( dense.xIonDense[nelem][dense.IonLow[nelem]]/dense.gas_phase[nelem] > 1e-3f &&
00279 dense.IonLow[nelem] > 0 )
00280 {
00281
00282 --dense.IonLow[nelem];
00283 lgLo_Down = true;
00284
00285 dense.xIonDense[nelem][dense.IonLow[nelem]] = (realnum)dense.density_low_limit;
00286 }
00287
00288
00289 else if( nzone < 10 &&
00290 (dense.xIonDense[nelem][dense.IonLow[nelem]]/dense.gas_phase[nelem] <= (realnum)trimlo) &&
00291 (dense.IonLow[nelem] < (dense.IonHigh[nelem] - 2) ) )
00292 {
00293
00294 dense.xIonDense[nelem][dense.IonLow[nelem]] = 0.;
00295
00296 thermal.heating[nelem][dense.IonLow[nelem]] = 0.;
00297 if( dense.IonLow[nelem] == nelem+1-NISO )
00298 {
00299 long int ipISO = nelem - dense.IonLow[nelem];
00300 ASSERT( ipISO>=0 && ipISO<NISO );
00301 StatesElem[ipISO][nelem][0].Pop = 0.;
00302 }
00303 ++dense.IonLow[nelem];
00304 lgLo_Up = true;
00305 }
00306
00307 else if( (dense.xIonDense[nelem][dense.IonLow[nelem]] < dense.density_low_limit) &&
00308
00309
00310
00311 (dense.IonLow[nelem] < dense.IonHigh[nelem]-1) )
00312 {
00313 while(dense.xIonDense[nelem][dense.IonLow[nelem]] < dense.density_low_limit &&
00314
00315
00316 (dense.IonLow[nelem] < dense.IonHigh[nelem]-1) )
00317 {
00318
00319 dense.xIonDense[nelem][dense.IonLow[nelem]] = 0.;
00320
00321 thermal.heating[nelem][dense.IonLow[nelem]] = 0.;
00322 if( dense.IonLow[nelem] == nelem+1-NISO )
00323 {
00324 long int ipISO = nelem - dense.IonLow[nelem];
00325 ASSERT( ipISO>=0 && ipISO<NISO );
00326 StatesElem[ipISO][nelem][0].Pop = 0.;
00327 }
00328 ++dense.IonLow[nelem];
00329 lgLo_Up = true;
00330 }
00331 }
00332
00333
00334
00335 ASSERT( dense.IonLow[nelem] < dense.IonHigh[nelem] ||
00336 (dense.IonLow[nelem]==0 && dense.IonHigh[nelem]==0 ) ||
00337 dense.lgSetIoniz[nelem] );
00338
00339
00340
00341
00342
00343
00344 ASSERT( dense.IonLow[nelem] <= 1 ||
00345 dense.xIonDense[nelem][dense.IonLow[nelem]-1] == 0. );
00346
00347 ASSERT( (dense.IonLow[nelem]==0 && dense.IonHigh[nelem]==0 ) || lgLo_Up ||
00348 dense.xIonDense[nelem][dense.IonLow[nelem]] >= dense.density_low_limit ||
00349 dense.xIonDense[nelem][dense.IonLow[nelem]]/dense.gas_phase[nelem] >= dense.density_low_limit ||
00350
00351
00352 (dense.IonLow[nelem]==dense.IonHigh[nelem]-1 ));
00353
00354 {
00355
00356 enum {DEBUG_LOC=false};
00357 if( DEBUG_LOC && nelem == ipHELIUM )
00358 {
00359 if( lgHi_Down ||lgHi_Up ||lgLo_Up ||lgLo_Down )
00360 {
00361 fprintf(ioQQQ,"DEBUG TrimZone\t%li\t",nzone );
00362 if( lgHi_Down )
00363 {
00364 fprintf(ioQQQ,"high dn %li to %li",
00365 ion_hi_save ,
00366 dense.IonHigh[nelem] );
00367 }
00368 if( lgHi_Up )
00369 {
00370 fprintf(ioQQQ,"high up %li to %li",
00371 ion_hi_save ,
00372 dense.IonHigh[nelem] );
00373 }
00374 if( lgLo_Up )
00375 {
00376 fprintf(ioQQQ,"low up %li to %li",
00377 ion_lo_save ,
00378 dense.IonLow[nelem] );
00379 }
00380 if( lgLo_Down )
00381 {
00382 fprintf(ioQQQ,"low dn %li to %li",
00383 ion_lo_save ,
00384 dense.IonLow[nelem] );
00385 }
00386 fprintf(ioQQQ,"\n" );
00387 }
00388 }
00389 }
00390
00391
00392
00393 if(dense.IonHigh[nelem] < nelem+1 )
00394 ASSERT( dense.xIonDense[nelem][dense.IonHigh[nelem]+1] == 0. );
00395
00396
00397 if( lgHi_Down || lgHi_Up || lgLo_Up || lgLo_Down )
00398 {
00399 conv.lgIonStageTrimed = true;
00400 {
00401
00402 enum {DEBUG_LOC=false};
00403 if( DEBUG_LOC && nelem==ipHELIUM )
00404 {
00405 fprintf(ioQQQ,"DEBUG ion_trim zone\t%.2f \t%li\t", fnzone, nelem);
00406 if( lgHi_Down )
00407 fprintf(ioQQQ,"\tHi_Down");
00408 if( lgHi_Up )
00409 fprintf(ioQQQ,"\tHi_Up");
00410 if( lgLo_Up )
00411 fprintf(ioQQQ,"\tLo_Up");
00412 if( lgLo_Down )
00413 fprintf(ioQQQ,"\tLo_Down");
00414 fprintf(ioQQQ,"\n");
00415 }
00416 }
00417 }
00418
00419
00420
00421 for( ion=0; ion<dense.IonLow[nelem]; ++ion )
00422 {
00423 ASSERT( dense.xIonDense[nelem][ion] == 0. );
00424 }
00425 for( ion=dense.IonHigh[nelem]+1; ion<nelem+1; ++ion )
00426 {
00427 ASSERT( dense.xIonDense[nelem][ion] == 0. );
00428 }
00429
00430 for( ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion )
00431 {
00432
00433
00434
00435 dense.xIonDense[nelem][ion] = MAX2(dense.xIonDense[nelem][ion], SMALLFLOAT);
00436 }
00437 return;
00438 }