00001
00002
00003
00004
00005 #include "cddefines.h"
00006 #include "phycon.h"
00007 #include "heavy.h"
00008 #include "hmi.h"
00009 #include "grainvar.h"
00010 #include "dense.h"
00011 #include "conv.h"
00012 #include "thermal.h"
00013 #include "iso.h"
00014 #include "abund.h"
00015 #include "save.h"
00016 #include "elementnames.h"
00017 #include "atmdat.h"
00018 #include "ionbal.h"
00019
00020
00021 void ion_recomb(
00022
00023 bool lgPrintIt,
00024 const double *dicoef,
00025 const double *dite,
00026 const double ditcrt[],
00027 const double aa[],
00028 const double bb[],
00029 const double cc[],
00030 const double dd[],
00031 const double ff[],
00032
00033 long int nelem
00034 )
00035 {
00036 #define DICOEF(I_,J_) (*(dicoef+(I_)*(nelem+1)+(J_)))
00037 #define DITE(I_,J_) (*(dite+(I_)*(nelem+1)+(J_)))
00038 long int i,
00039 ion,
00040 limit;
00041 double
00042 fac2,
00043 factor,
00044 DielRecomRateCoef_HiT[LIMELM],
00045 DielRecomRateCoef_LowT[LIMELM],
00046 ChargeTransfer[LIMELM];
00047
00048
00049 static double RecNoise[LIMELM][LIMELM];
00050 static bool lgNoiseNeedEval=true;
00051
00052 DEBUG_ENTRY( "ion_recomb(false,)" );
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063 ASSERT( nelem < LIMELM);
00064 ASSERT( nelem > 1 );
00065
00066
00067 ASSERT( dense.IonLow[nelem] >= 0 );
00068 ASSERT( dense.IonLow[nelem] <= nelem+1 );
00069
00070 atmdat.nsbig = MAX2(dense.IonHigh[nelem]+1,atmdat.nsbig);
00071 fac2 = 1e-14*phycon.sqrte;
00072
00073
00074
00075
00076 if( lgNoiseNeedEval )
00077 {
00078 int n;
00079 if( ionbal.guess_noise !=0. )
00080 {
00081 for( n=ipHYDROGEN; n<LIMELM; ++n )
00082 {
00083 for( ion=0; ion<=n; ++ion )
00084 {
00085
00086
00087 RecNoise[n][ion] = pow(10., RandGauss( 0. , ionbal.guess_noise ) );
00088 }
00089 }
00090 }
00091 else
00092 {
00093 for( n=ipHYDROGEN; n<LIMELM; ++n )
00094 {
00095 for( ion=0; ion<=n; ++ion )
00096 {
00097 RecNoise[n][ion] = 1.;
00098 }
00099 }
00100 }
00101 lgNoiseNeedEval = false;
00102 }
00103
00104
00105
00106
00107
00108
00109 limit = MIN2(nelem-NISO,dense.IonHigh[nelem]-1);
00110 ASSERT( limit >= -1 );
00111
00112
00113
00114
00115 for( ion=0; ion <= limit; ion++ )
00116 {
00117 ionbal.RateRecomTot[nelem][ion] = 0.;
00118 ChargeTransfer[ion] = 0.;
00119 DielRecomRateCoef_LowT[ion] = 0.;
00120 DielRecomRateCoef_HiT[ion] = 0.;
00121 ionbal.RR_rate_coef_used[nelem][ion] = 0.;
00122 ionbal.DR_rate_coef_used[nelem][ion] = 0.;
00123 }
00124 for( ion=limit+1; ion < LIMELM; ion++ )
00125 {
00126
00127
00128
00129 ChargeTransfer[ion] = -FLT_MAX;
00130 DielRecomRateCoef_LowT[ion] = -FLT_MAX;
00131 DielRecomRateCoef_HiT[ion] = -FLT_MAX;
00132 }
00133
00134 DielRecomRateCoef_HiT[nelem] = 0.;
00135 DielRecomRateCoef_HiT[nelem-1] = 0.;
00136
00137 DielRecomRateCoef_LowT[nelem] = 0.;
00138 DielRecomRateCoef_LowT[nelem-1] = 0.;
00139
00140
00141 Heavy.xLyaHeavy[nelem][nelem] = 0.;
00142 Heavy.xLyaHeavy[nelem][nelem-1] = 0.;
00143
00144
00145 for( ion=dense.IonLow[nelem]; ion <= limit; ion++ )
00146 {
00147
00148
00149
00150
00151 long n_bnd_elec_after_recomb = nelem+1 - ion;
00152
00153 #if 0
00154
00155
00156 ChargeTransfer[ion] =
00157
00158 atmdat.HeCharExcRecTo[nelem][ion]*
00159
00160 StatesElemNEW[ipHELIUM][ipHELIUM-ipHE_LIKE][ipHe1s1S].Pop +
00161
00162 atmdat.HCharExcRecTo[nelem][ion]*
00163
00164 StatesElemNEW[ipHYDROGEN][ipHYDROGEN-ipH_LIKE][ipH1s].Pop;
00165 #else
00166 ChargeTransfer[ion] =
00167
00168 atmdat.HeCharExcRecTo[nelem][ion]*
00169
00170 dense.xIonDense[ipHELIUM][0] +
00171
00172 atmdat.HCharExcRecTo[nelem][ion]*
00173
00174 dense.xIonDense[ipHYDROGEN][0];
00175 #endif
00176
00177
00178
00179 if( ion==0 && nelem>ipHELIUM && atmdat.lgCTOn )
00180 ChargeTransfer[ion] += hmi.hmin_ct_firstions * hmi.Hmolec[ipMHm];
00181
00182
00183 if( ionbal.lgRR_recom_Badnell_use && ionbal.lgRR_Badnell_rate_coef_exist[nelem][ion] )
00184 {
00185 ionbal.RR_rate_coef_used[nelem][ion] = ionbal.RR_Badnell_rate_coef[nelem][ion];
00186 }
00187 else
00188 {
00189 ionbal.RR_rate_coef_used[nelem][ion] = ionbal.RR_Verner_rate_coef[nelem][ion];
00190 }
00191
00192
00193 DielRecomRateCoef_HiT[ion] = 0.;
00194
00195 if( nelem==ipIRON )
00196 {
00197
00198 DielRecomRateCoef_HiT[ion] = atmdat_dielrec_fe(ion+1,phycon.te);
00199 }
00200 else if( phycon.te > (ditcrt[ion]*0.1) )
00201 {
00202 DielRecomRateCoef_HiT[ion] = ionbal.DielSupprs[0][ion]/phycon.te32*
00203 DICOEF(0,ion)*exp(-DITE(0,ion)/phycon.te)*
00204 (1. + DICOEF(1,ion)*
00205 sexp(DITE(1,ion)/phycon.te));
00206 }
00207
00208
00209
00210
00211 DielRecomRateCoef_LowT[ion] = 0.;
00212 if( ((n_bnd_elec_after_recomb-1) != 2) &&
00213 ((n_bnd_elec_after_recomb-1) != 10) &&
00214 ((n_bnd_elec_after_recomb-1) != 18) )
00215 {
00216
00217 if( ff[ion] != 0. && nelem != ipIRON )
00218 {
00219 double t4m1 = 1e4/phycon.te;
00220 double tefac = ff[ion]*t4m1;
00221
00222
00223
00224 if( tefac > -5. )
00225 {
00226 factor = (((aa[ion]*t4m1+bb[ion])*t4m1+cc[ion])*t4m1+dd[ion])* sexp(tefac);
00227 DielRecomRateCoef_LowT[ion] = ionbal.DielSupprs[1][ion]*fac2*
00228 MAX2(0.,factor );
00229 }
00230 else
00231 {
00232 DielRecomRateCoef_LowT[ion] = 0.;
00233 }
00234 }
00235 else if( ionbal.lg_guess_coef )
00236 {
00237
00238 DielRecomRateCoef_LowT[ion] = ionbal.DR_Badnell_rate_coef_mean_ion[ion];
00239
00240
00241
00242 DielRecomRateCoef_LowT[ion] *= RecNoise[nelem][ion];
00243 }
00244 }
00245
00246
00247 ionbal.DR_old_rate_coef[nelem][ion] = DielRecomRateCoef_HiT[ion] + DielRecomRateCoef_LowT[ion];
00248
00249
00250 if( ionbal.lgDR_recom_Badnell_use && ionbal.lgDR_Badnell_rate_coef_exist[nelem][ion] )
00251 {
00252 ionbal.DR_rate_coef_used[nelem][ion] = ionbal.DR_Badnell_rate_coef[nelem][ion];
00253 }
00254 else
00255 {
00256 ionbal.DR_rate_coef_used[nelem][ion] = ionbal.DR_old_rate_coef[nelem][ion];
00257 }
00258
00259
00260 ionbal.RateRecomTot[nelem][ion] =
00261 dense.eden* (
00262 ionbal.RR_rate_coef_used[nelem][ion] +
00263 ionbal.DR_rate_coef_used[nelem][ion] +
00264 ionbal.CotaRate[ion] ) +
00265 ChargeTransfer[ion];
00266
00267
00268
00269 # define FRAC_LINE 1.
00270
00271
00272 Heavy.xLyaHeavy[nelem][ion] = (realnum)(dense.eden*
00273 (ionbal.RR_rate_coef_used[nelem][ion]+ionbal.DR_rate_coef_used[nelem][ion])*FRAC_LINE );
00274 }
00275
00276
00277 if( save.lgioRecom || lgPrintIt )
00278 {
00279
00280 FILE *ioOut;
00281 if( lgPrintIt )
00282 ioOut = ioQQQ;
00283 else
00284 ioOut = save.ioRecom;
00285
00286
00287 fprintf( ioOut,
00288 " %s recombination coefficients fnzone:%.2f \tte\t%.4e\tne\t%.4e\n",
00289 elementnames.chElementName[nelem] , fnzone , phycon.te , dense.eden );
00290
00291
00292
00293 limit = dense.IonHigh[nelem];
00294
00295 for( i=0; i<limit; ++i )
00296 fprintf( ioOut, "%10ld",i+1);
00297 fprintf( ioOut, "\n");
00298
00299 for( i=0; i < limit; i++ )
00300 {
00301 fprintf( ioOut, "%10.2e", ionbal.RR_rate_coef_used[nelem][i] );
00302 }
00303 fprintf( ioOut, " radiative used vs Z\n" );
00304
00305 for( i=0; i < limit; i++ )
00306 {
00307 fprintf( ioOut, "%10.2e", ionbal.RR_Verner_rate_coef[nelem][i] );
00308 }
00309 fprintf( ioOut, " old Verner vs Z\n" );
00310
00311 for( i=0; i < limit; i++ )
00312 {
00313 fprintf( ioOut, "%10.2e", ionbal.RR_Badnell_rate_coef[nelem][i] );
00314 }
00315 fprintf( ioOut, " new Badnell vs Z\n" );
00316
00317 for( i=0; i < limit; i++ )
00318 {
00319
00320
00321
00322 fprintf( ioOut, "%10.2e", ChargeTransfer[i]/SDIV(dense.xIonDense[ipHYDROGEN][0]) );
00323 }
00324 fprintf( ioOut, " CT/n(H0)\n" );
00325
00326 for( i=0; i < limit; i++ )
00327 {
00328 fprintf( ioOut, "%10.2e", ionbal.CotaRate[ion] );
00329 }
00330 fprintf( ioOut, " 3body vs Z /ne\n" );
00331
00332
00333 for( i=0; i < dense.IonHigh[nelem]; i++ )
00334 {
00335 fprintf( ioOut, "%10.2e", gv.GrainChTrRate[nelem][i+1][i]/dense.eden );
00336 }
00337 fprintf( ioOut, " Grain vs Z /ne\n" );
00338
00339 for( i=0; i < limit; i++ )
00340 {
00341 fprintf( ioOut, "%10.2e", DielRecomRateCoef_HiT[i] );
00342 }
00343 fprintf( ioOut, " Burgess vs Z\n" );
00344
00345 for( i=0; i < limit; i++ )
00346 {
00347 fprintf( ioOut, "%10.2e", ionbal.DR_old_rate_coef[nelem][i] );
00348 }
00349 fprintf( ioOut, " old Nussbaumer Storey DR vs Z\n" );
00350
00351 for( i=0; i < limit; i++ )
00352 {
00353 fprintf( ioOut, "%10.2e", ionbal.DR_Badnell_rate_coef[nelem][i] );
00354 }
00355 fprintf( ioOut, " new Badnell DR vs Z\n" );
00356
00357 for( i=0; i < limit; i++ )
00358 {
00359 fprintf( ioOut, "%10.2e", DielRecomRateCoef_LowT[i] );
00360 }
00361 fprintf( ioOut, " low T DR used vs Z\n" );
00362
00363
00364 for( i=0; i < limit; i++ )
00365 {
00366 fprintf( ioOut, "%10.2e", ionbal.RateRecomTot[nelem][i] );
00367 }
00368 fprintf( ioOut,
00369 " total rec rate (with density) for %s\n",
00370 elementnames.chElementSym[nelem] );
00371 for( i=0; i < limit; i++ )
00372 {
00373 fprintf( ioOut, "%10.2e", ionbal.RateRecomTot[nelem][i]/dense.eden );
00374 }
00375 fprintf( ioOut,
00376 " total rec rate / ne for %s\n\n",
00377 elementnames.chElementSym[nelem] );
00378
00379
00380 if( dense.IonHigh[nelem] > 11 )
00381 {
00382 limit = MIN2(29,dense.IonHigh[nelem]);
00383 fprintf( ioOut, " R " );
00384 for( i=11; i < limit; i++ )
00385 {
00386 fprintf( ioOut, "%10.2e", dense.eden*ionbal.CotaRate[ion] );
00387 }
00388 fprintf( ioOut, "\n" );
00389
00390 fprintf( ioOut, " B " );
00391 for( i=11; i < limit; i++ )
00392 {
00393 fprintf( ioOut, "%10.2e", DielRecomRateCoef_HiT[i] );
00394 }
00395 fprintf( ioOut, "\n" );
00396
00397 fprintf( ioOut, " NS" );
00398 for( i=11; i < limit; i++ )
00399 {
00400 fprintf( ioOut, "%10.2e", DielRecomRateCoef_LowT[i] );
00401 }
00402 fprintf( ioOut, "\n" );
00403
00404 fprintf( ioOut, " " );
00405 for( i=11; i < limit; i++ )
00406 {
00407 fprintf( ioOut, "%10.2e", ionbal.RateRecomTot[nelem][i] );
00408 }
00409 fprintf( ioOut, "\n\n" );
00410 }
00411 }
00412
00413
00414
00415 limit = MIN2(nelem-NISO,dense.IonHigh[nelem]-1);
00416 for( i=dense.IonLow[nelem]; i <= limit; i++ )
00417 {
00418 ASSERT( Heavy.xLyaHeavy[nelem][i] > 0. );
00419 ASSERT( ionbal.RateRecomTot[nelem][i] > 0. );
00420 }
00421 return;
00422 # undef DITE
00423 # undef DICOEF
00424 }
00425
00426
00427 void ion_recombAGN( FILE * io )
00428 {
00429 # define N1LIM 3
00430 # define N2LIM 4
00431 double te1[N1LIM]={ 5000., 10000., 20000.};
00432 double te2[N2LIM]={ 20000.,50000.,100000.,1e6};
00433
00434 double BreakEnergy = 100./13.0;
00435 long int nelem, ion , i;
00436
00437 char chString[100],
00438 chOutput[100];
00439
00440 double TempSave = phycon.te;
00441
00442 double EdenSave = dense.eden;
00443
00444 DEBUG_ENTRY( "ion_recomb(false,)" );
00445
00446 dense.eden = 1.;
00447
00448
00449
00450 fprintf(io,"X+i\\Te");
00451 for( i=0; i<N1LIM; ++i )
00452 {
00453 phycon.te = te1[i];
00454 fprintf(io,"\t%.0f K",phycon.te);
00455 }
00456 fprintf(io,"\n");
00457
00458
00459 for( nelem=ipLITHIUM; nelem<LIMELM; ++nelem )
00460 {
00461
00462 if( abund.lgAGN[nelem] )
00463 {
00464 for( ion=0; ion<=nelem; ++ion )
00465 {
00466 ASSERT( Heavy.Valence_IP_Ryd[nelem][ion] > 0.05 );
00467
00468 if( Heavy.Valence_IP_Ryd[nelem][ion] > BreakEnergy )
00469 break;
00470
00471
00472 sprintf(chOutput,"%s",
00473 elementnames.chElementSym[nelem]);
00474
00475 if( chOutput[1]==' ' )
00476 chOutput[1] = chOutput[2];
00477
00478 if( ion==0 )
00479 {
00480 sprintf(chString,"0 ");
00481 }
00482 else if( ion==1 )
00483 {
00484 sprintf(chString,"+ ");
00485 }
00486 else
00487 {
00488 sprintf(chString,"+%li ",ion);
00489 }
00490 strcat( chOutput , chString );
00491 fprintf(io,"%5s",chOutput );
00492
00493 for( i=0; i<N1LIM; ++i )
00494 {
00495 TempChange(te1[i] , false);
00496 dense.IonLow[nelem] = 0;
00497 dense.IonHigh[nelem] = nelem+1;
00498 if( ConvBase(0) )
00499 fprintf(ioQQQ,"PROBLEM ConvBase returned error.\n");
00500 fprintf(io,"\t%.2e",ionbal.RateRecomTot[nelem][ion]);
00501 }
00502 fprintf(io,"\n");
00503 }
00504 fprintf(io,"\n");
00505 }
00506 }
00507
00508
00509 fprintf(io,"X+i\\Te");
00510 for( i=0; i<N2LIM; ++i )
00511 {
00512 TempChange(te2[i] , false);
00513 fprintf(io,"\t%.0f K",phycon.te);
00514 }
00515 fprintf(io,"\n");
00516
00517
00518 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
00519 {
00520
00521 if( abund.lgAGN[nelem] )
00522 {
00523 for( ion=0; ion<=nelem; ++ion )
00524 {
00525 ASSERT( Heavy.Valence_IP_Ryd[nelem][ion] > 0.05 );
00526
00527 if( Heavy.Valence_IP_Ryd[nelem][ion] <= BreakEnergy )
00528 continue;
00529
00530
00531 fprintf(io,"%s",
00532 elementnames.chElementSym[nelem]);
00533
00534 if( ion==0 )
00535 {
00536 fprintf(io,"0 ");
00537 }
00538 else if( ion==1 )
00539 {
00540 fprintf(io,"+ ");
00541 }
00542 else
00543 {
00544 fprintf(io,"+%li",ion);
00545 }
00546
00547 for( i=0; i<N2LIM; ++i )
00548 {
00549 TempChange(te2[i] , false);
00550 dense.IonLow[nelem] = 0;
00551 dense.IonHigh[nelem] = nelem+1;
00552 if( ConvBase(0) )
00553 fprintf(ioQQQ,"PROBLEM ConvBase returned error.\n");
00554 fprintf(io,"\t%.2e",ionbal.RateRecomTot[nelem][ion]);
00555 }
00556 fprintf(io,"\n");
00557 }
00558 fprintf(io,"\n");
00559 }
00560 }
00561
00562 TempChange(TempSave , true);
00563 dense.eden = EdenSave;
00564 return;
00565 }