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 ASSERT( dense.IonLow[nelem] <= dense.IonHigh[nelem] );
00056
00057
00058
00059
00060 if( conv.lgSearch )
00061 {
00062 trimhi = (realnum)ionbal.trimhi * 1e-4f;
00063 trimlo = (realnum)ionbal.trimlo * 1e-4f;
00064 }
00065 else
00066 {
00067 trimhi = (realnum)ionbal.trimhi;
00068 trimlo = (realnum)ionbal.trimlo;
00069 }
00070
00071
00072
00073 if( nelem == ipHELIUM )
00074 {
00075
00076 trimlo = SMALLFLOAT;
00077
00078
00079
00080 if( dense.IonHigh[ipHELIUM] == 1 )
00081 {
00082 trimhi = MIN2( trimhi , 1e-20f );
00083 }
00084 else if( dense.IonHigh[ipHELIUM] == 2 )
00085 {
00086 if( conv.lgSearch )
00087 {
00088
00089
00090
00091
00092
00093
00094 trimhi = MIN2( trimhi , 1e-17f );
00095 }
00096 else
00097 {
00098
00099 trimhi = MIN2( trimhi , 1e-12f );
00100 }
00101 }
00102 }
00103
00104
00105
00106 if( !mole_global.lgNoMole )
00107 {
00108 trimlo = SMALLFLOAT;
00109 if(dense.IonHigh[nelem] ==2)
00110 {
00111 trimhi = MIN2(trimhi, 1e-20f);
00112 }
00113 }
00114
00115
00116
00117
00118
00119
00120
00121 if( conv.lgSearch )
00122 {
00123
00124 while(
00125 (dense.xIonDense[nelem][dense.IonHigh[nelem]]/dense.gas_phase[nelem] < dense.density_low_limit ||
00126 dense.xIonDense[nelem][dense.IonHigh[nelem]] < dense.density_low_limit ) &&
00127
00128 dense.IonHigh[nelem] > dense.IonLow[nelem] )
00129 {
00130
00131
00132 dense.xIonDense[nelem][dense.IonHigh[nelem]-1] +=
00133 dense.xIonDense[nelem][dense.IonHigh[nelem]];
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 for( long level = 0; level < iso_sp[ipISO][nelem].numLevels_max; level++ )
00141 iso_sp[ipISO][nelem].st[level].Pop() = 0.;
00142 }
00143
00144
00145 --dense.IonHigh[nelem];
00146 ASSERT( dense.IonHigh[nelem] >= 0);
00147
00148 lgHi_Down = true;
00149 }
00150
00151
00152 while(
00153 (dense.xIonDense[nelem][dense.IonLow[nelem]]/dense.gas_phase[nelem] < dense.density_low_limit ||
00154 dense.xIonDense[nelem][dense.IonLow[nelem]] < dense.density_low_limit ) &&
00155 dense.IonLow[nelem] < dense.IonHigh[nelem] - 1 )
00156 {
00157
00158
00159 dense.xIonDense[nelem][dense.IonLow[nelem]+1] +=
00160 dense.xIonDense[nelem][dense.IonLow[nelem]];
00161 dense.xIonDense[nelem][dense.IonLow[nelem]] = 0.;
00162
00163
00164 thermal.heating[nelem][dense.IonLow[nelem]] = 0.;
00165 if( dense.IonLow[nelem] == nelem+1-NISO )
00166 {
00167 long int ipISO = nelem - dense.IonLow[nelem];
00168 ASSERT( ipISO>=0 && ipISO<NISO );
00169 for( long level = 0; level < iso_sp[ipISO][nelem].numLevels_max; level++ )
00170 iso_sp[ipISO][nelem].st[level].Pop() = 0.;
00171 }
00172
00173
00174 ++dense.IonLow[nelem];
00175 lgLo_Up = true;
00176 }
00177 }
00178
00179
00180 ASSERT( dense.IonLow[nelem] <= dense.IonHigh[nelem] );
00181
00182
00183
00184
00185 if( dense.IonHigh[nelem] > 1 &&
00186 (dense.IonLow[nelem]==dense.IonHigh[nelem]-1) &&
00187 (dense.xIonDense[nelem][dense.IonLow[nelem]] < dense.density_low_limit) )
00188 {
00189 # if 0
00190 fprintf(ioQQQ,
00191 "PROBLEM DISASTER the density of ion %li of element %s is too small to be computed on this cpu.\n",
00192 dense.IonLow[nelem]+1,
00193 elementnames.chElementName[nelem]);
00194 fprintf(ioQQQ,
00195 "Turn off the element with the command ELEMENT XXX OFF or do not consider "
00196 "gas with low density, the current hydrogen density is %.2e cm-3.\n",
00197 dense.gas_phase[ipHYDROGEN]);
00198 fprintf(ioQQQ,
00199 "The outer radius of the cloud is %.2e cm - should a stopping "
00200 "radius be set?\n",
00201 radius.Radius );
00202 fprintf(ioQQQ,
00203 "The abort flag is being set - the calculation is stopping.\n");
00204 # endif
00205
00206 dense.xIonDense[nelem][dense.IonLow[nelem]] = (realnum)dense.density_low_limit;
00207 return;
00208 }
00209
00210
00211 if(
00212 dense.IonHigh[nelem] > dense.IonLow[nelem] &&
00213 ( (dense.xIonDense[nelem][dense.IonHigh[nelem]]/dense.gas_phase[nelem] <=
00214 trimhi ) ||
00215 (dense.xIonDense[nelem][dense.IonHigh[nelem]] <= dense.density_low_limit ) )
00216 )
00217 {
00218
00219
00220 if( dense.IonHigh[nelem]>1 ||
00221 (dense.IonHigh[nelem]==1&&dense.xIonDense[nelem][1]<100.*dense.density_low_limit) )
00222 {
00223 dense.xIonDense[nelem][dense.IonHigh[nelem]-1] +=
00224 dense.xIonDense[nelem][dense.IonHigh[nelem]];
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_sp[ipISO][nelem].numLevels_max; level++ )
00232 iso_sp[ipISO][nelem].st[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 dense.xIonDense[nelem][dense.IonHigh[nelem]-1] -=
00276 dense.xIonDense[nelem][dense.IonHigh[nelem]];
00277 long int ipISO = nelem - dense.IonHigh[nelem];
00278 if (ipISO >= 0 && ipISO < NISO)
00279 iso_sp[ipISO][nelem].st[0].Pop() = dense.xIonDense[nelem][dense.IonHigh[nelem]];
00280 }
00281 }
00282
00283
00284 ASSERT( dense.IonLow[nelem] <= dense.IonHigh[nelem] );
00285
00286
00287 if( dense.xIonDense[nelem][dense.IonLow[nelem]]/dense.gas_phase[nelem] > 1e-3f &&
00288 dense.IonLow[nelem] > 0 )
00289 {
00290
00291 --dense.IonLow[nelem];
00292 lgLo_Down = true;
00293
00294 dense.xIonDense[nelem][dense.IonLow[nelem]] = (realnum)dense.density_low_limit;
00295 dense.xIonDense[nelem][dense.IonLow[nelem]+1] -=
00296 dense.xIonDense[nelem][dense.IonLow[nelem]];
00297 long int ipISO = nelem - dense.IonLow[nelem];
00298 if (ipISO >= 0 && ipISO < NISO)
00299 iso_sp[ipISO][nelem].st[0].Pop() = dense.xIonDense[nelem][dense.IonLow[nelem]];
00300 }
00301
00302
00303 else if( nzone < 10 &&
00304 (dense.xIonDense[nelem][dense.IonLow[nelem]]/dense.gas_phase[nelem] <= (realnum)trimlo) &&
00305 (dense.IonLow[nelem] < (dense.IonHigh[nelem] - 2) ) )
00306 {
00307
00308 dense.xIonDense[nelem][dense.IonLow[nelem]+1] +=
00309 dense.xIonDense[nelem][dense.IonLow[nelem]];
00310 dense.xIonDense[nelem][dense.IonLow[nelem]] = 0.;
00311
00312 thermal.heating[nelem][dense.IonLow[nelem]] = 0.;
00313 if( dense.IonLow[nelem] == nelem+1-NISO )
00314 {
00315 long int ipISO = nelem - dense.IonLow[nelem];
00316 ASSERT( ipISO>=0 && ipISO<NISO );
00317 for( long level = 0; level < iso_sp[ipISO][nelem].numLevels_max; level++ )
00318 iso_sp[ipISO][nelem].st[level].Pop() = 0.;
00319 }
00320 ++dense.IonLow[nelem];
00321 lgLo_Up = true;
00322 }
00323
00324 else if( (dense.xIonDense[nelem][dense.IonLow[nelem]] < dense.density_low_limit) &&
00325
00326
00327
00328 (dense.IonLow[nelem] < dense.IonHigh[nelem]-1) )
00329 {
00330 while(dense.xIonDense[nelem][dense.IonLow[nelem]] < dense.density_low_limit &&
00331
00332
00333 (dense.IonLow[nelem] < dense.IonHigh[nelem]-1) )
00334 {
00335
00336 dense.xIonDense[nelem][dense.IonLow[nelem]+1] +=
00337 dense.xIonDense[nelem][dense.IonLow[nelem]];
00338 dense.xIonDense[nelem][dense.IonLow[nelem]] = 0.;
00339
00340 thermal.heating[nelem][dense.IonLow[nelem]] = 0.;
00341 if( dense.IonLow[nelem] == nelem+1-NISO )
00342 {
00343 long int ipISO = nelem - dense.IonLow[nelem];
00344 ASSERT( ipISO>=0 && ipISO<NISO );
00345 for( long level = 0; level < iso_sp[ipISO][nelem].numLevels_max; level++ )
00346 iso_sp[ipISO][nelem].st[level].Pop() = 0.;
00347 }
00348 ++dense.IonLow[nelem];
00349 lgLo_Up = true;
00350 }
00351 }
00352
00353
00354 ASSERT( dense.IonLow[nelem] <= dense.IonHigh[nelem] );
00355
00356
00357
00358
00359
00360
00361 ASSERT( dense.IonLow[nelem] <= 1 ||
00362 dense.xIonDense[nelem][dense.IonLow[nelem]-1] == 0. );
00363
00364 ASSERT( (dense.IonLow[nelem]==0 && dense.IonHigh[nelem]==0 ) || lgLo_Up ||
00365 dense.xIonDense[nelem][dense.IonLow[nelem]] >= dense.density_low_limit ||
00366 dense.xIonDense[nelem][dense.IonLow[nelem]]/dense.gas_phase[nelem] >= dense.density_low_limit ||
00367
00368
00369 (dense.IonLow[nelem]==dense.IonHigh[nelem]-1 ));
00370
00371 {
00372
00373 enum {DEBUG_LOC=false};
00374 if( DEBUG_LOC && nelem == ipHELIUM )
00375 {
00376 if( lgHi_Down ||lgHi_Up ||lgLo_Up ||lgLo_Down )
00377 {
00378 fprintf(ioQQQ,"DEBUG TrimZone\t%li\t",nzone );
00379 if( lgHi_Down )
00380 {
00381 fprintf(ioQQQ,"high dn %li to %li",
00382 ion_hi_save ,
00383 dense.IonHigh[nelem] );
00384 }
00385 if( lgHi_Up )
00386 {
00387 fprintf(ioQQQ,"high up %li to %li",
00388 ion_hi_save ,
00389 dense.IonHigh[nelem] );
00390 }
00391 if( lgLo_Up )
00392 {
00393 fprintf(ioQQQ,"low up %li to %li",
00394 ion_lo_save ,
00395 dense.IonLow[nelem] );
00396 }
00397 if( lgLo_Down )
00398 {
00399 fprintf(ioQQQ,"low dn %li to %li",
00400 ion_lo_save ,
00401 dense.IonLow[nelem] );
00402 }
00403 fprintf(ioQQQ,"\n" );
00404 }
00405 }
00406 }
00407
00408
00409
00410 if(dense.IonHigh[nelem] < nelem+1 )
00411 ASSERT( dense.xIonDense[nelem][dense.IonHigh[nelem]+1] == 0. );
00412
00413
00414 if( lgHi_Down || lgHi_Up || lgLo_Up || lgLo_Down )
00415 {
00416 conv.lgIonStageTrimed = true;
00417 {
00418
00419 enum {DEBUG_LOC=false};
00420 if( DEBUG_LOC && nelem==ipHELIUM )
00421 {
00422 fprintf(ioQQQ,"DEBUG ion_trim zone\t%.2f \t%li\t", fnzone, nelem);
00423 if( lgHi_Down )
00424 fprintf(ioQQQ,"\tHi_Down");
00425 if( lgHi_Up )
00426 fprintf(ioQQQ,"\tHi_Up");
00427 if( lgLo_Up )
00428 fprintf(ioQQQ,"\tLo_Up");
00429 if( lgLo_Down )
00430 fprintf(ioQQQ,"\tLo_Down");
00431 fprintf(ioQQQ,"\n");
00432 }
00433 }
00434 }
00435
00436
00437
00438 for( ion=0; ion<dense.IonLow[nelem]; ++ion )
00439 {
00440 ASSERT( dense.xIonDense[nelem][ion] == 0. );
00441 }
00442 for( ion=dense.IonHigh[nelem]+1; ion<nelem+1; ++ion )
00443 {
00444 ASSERT( dense.xIonDense[nelem][ion] == 0. );
00445 }
00446
00447 for( ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion )
00448 {
00449
00450
00451
00452 dense.xIonDense[nelem][ion] = MAX2(dense.xIonDense[nelem][ion], SMALLFLOAT);
00453 }
00454 return;
00455 }