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