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