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