00001
00002
00003
00004
00005 #include "cddefines.h"
00006 #include "phycon.h"
00007 #include "heavy.h"
00008 #include "mole.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 "save.h"
00017 #include "elementnames.h"
00018 #include "atmdat.h"
00019 #include "ionbal.h"
00020
00021
00022 void ion_recomb(
00023
00024 bool lgPrintIt,
00025
00026 long int nelem
00027 )
00028 {
00029 long int i,
00030 ion,
00031 limit;
00032
00033 DEBUG_ENTRY( "ion_recomb(false,)" );
00034
00035 ASSERT( nelem < LIMELM);
00036 ASSERT( nelem > 1 );
00037
00038
00039 ASSERT( dense.IonLow[nelem] >= 0 );
00040 ASSERT( dense.IonLow[nelem] <= nelem+1 );
00041
00042 atmdat.nsbig = MAX2(dense.IonHigh[nelem]+1,atmdat.nsbig);
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053 for( ion=0; ion <= dense.IonHigh[nelem]-1; ion++ )
00054 {
00055 ionbal.RateRecomTot[nelem][ion] = 0.;
00056 }
00057
00058 for( long ipISO=ipH_LIKE; ipISO<MIN2(NISO,nelem+1); ipISO++ )
00059 {
00060 if( ((dense.IonHigh[nelem] >= nelem - ipISO) &&
00061 (dense.IonLow[nelem] <= nelem - ipISO)) )
00062 {
00063 ionbal.RateRecomTot[nelem][nelem-ipISO] = ionbal.RateRecomIso[nelem][ipISO];
00064 }
00065 }
00066
00067 limit = MIN2(nelem-NISO,dense.IonHigh[nelem]-1);
00068 ASSERT( limit >= -1 );
00069
00070
00071 Heavy.xLyaHeavy[nelem][nelem] = 0.;
00072 Heavy.xLyaHeavy[nelem][nelem-1] = 0.;
00073
00074
00075 for( ion=dense.IonLow[nelem]; ion <= limit; ion++ )
00076 {
00077
00078
00079
00080
00081 ASSERT(ionbal.DR_Badnell_rate_coef[nelem][ion] >= 0);
00082 ASSERT(ionbal.RR_rate_coef_used[nelem][ion] >= 0);
00083
00084
00085 ionbal.RateRecomTot[nelem][ion] =
00086 dense.eden* (
00087 ionbal.RR_rate_coef_used[nelem][ion] +
00088 ionbal.DR_Badnell_rate_coef[nelem][ion] +
00089 ionbal.CotaRate[ion] );
00090
00091
00092
00093 # define FRAC_LINE 1.
00094
00095
00096 Heavy.xLyaHeavy[nelem][ion] = (realnum)(dense.eden*
00097 (ionbal.RR_rate_coef_used[nelem][ion]+ionbal.DR_Badnell_rate_coef[nelem][ion])*FRAC_LINE );
00098 }
00099
00100
00101 if( save.lgioRecom || lgPrintIt )
00102 {
00103
00104 FILE *ioOut;
00105 if( lgPrintIt )
00106 ioOut = ioQQQ;
00107 else
00108 ioOut = save.ioRecom;
00109
00110
00111 fprintf( ioOut,
00112 " %s recombination coefficients fnzone:%.2f \tte\t%.4e\tne\t%.4e\n",
00113 elementnames.chElementName[nelem] , fnzone , phycon.te , dense.eden );
00114
00115
00116
00117 limit = dense.IonHigh[nelem];
00118
00119 for( i=0; i<limit; ++i )
00120 fprintf( ioOut, "%10ld",i+1);
00121 fprintf( ioOut, "\n");
00122
00123 for( i=0; i < limit; i++ )
00124 {
00125 fprintf( ioOut, "%10.2e", ionbal.RR_rate_coef_used[nelem][i] );
00126 }
00127 fprintf( ioOut, " radiative used vs Z\n" );
00128
00129 for( i=0; i < limit; i++ )
00130 {
00131 fprintf( ioOut, "%10.2e", ionbal.RR_Verner_rate_coef[nelem][i] );
00132 }
00133 fprintf( ioOut, " old Verner vs Z\n" );
00134
00135 for( i=0; i < limit; i++ )
00136 {
00137 fprintf( ioOut, "%10.2e", ionbal.RR_Badnell_rate_coef[nelem][i] );
00138 }
00139 fprintf( ioOut, " new Badnell vs Z\n" );
00140
00141 for( i=0; i < limit; i++ )
00142 {
00143
00144
00145
00146 fprintf( ioOut, "%10.2e", ionbal.CX_recomb_rate_used[nelem][i]/SDIV(dense.xIonDense[ipHYDROGEN][0]) );
00147 }
00148 fprintf( ioOut, " CT/n(H0)\n" );
00149
00150 for( i=0; i < limit; i++ )
00151 {
00152 fprintf( ioOut, "%10.2e", ionbal.CotaRate[ion] );
00153 }
00154 fprintf( ioOut, " 3body vs Z /ne\n" );
00155
00156
00157 for( i=0; i < dense.IonHigh[nelem]; i++ )
00158 {
00159 fprintf( ioOut, "%10.2e", gv.GrainChTrRate[nelem][i+1][i]/dense.eden );
00160 }
00161 fprintf( ioOut, " Grain vs Z /ne\n" );
00162 fprintf( ioOut, " old Nussbaumer Storey DR vs Z\n" );
00163
00164 for( i=0; i < limit; i++ )
00165 {
00166 fprintf( ioOut, "%10.2e", ionbal.DR_Badnell_rate_coef[nelem][i] );
00167 }
00168 fprintf( ioOut, " new Badnell DR vs Z\n" );
00169
00170
00171 for( i=0; i < limit; i++ )
00172 {
00173 fprintf( ioOut, "%10.2e", ionbal.RateRecomTot[nelem][i] );
00174 }
00175 fprintf( ioOut,
00176 " total rec rate (with density) for %s\n",
00177 elementnames.chElementSym[nelem] );
00178 for( i=0; i < limit; i++ )
00179 {
00180 fprintf( ioOut, "%10.2e", ionbal.RateRecomTot[nelem][i]/dense.eden );
00181 }
00182 fprintf( ioOut,
00183 " total rec rate / ne for %s\n\n",
00184 elementnames.chElementSym[nelem] );
00185
00186
00187 if( dense.IonHigh[nelem] > 11 )
00188 {
00189 limit = MIN2(29,dense.IonHigh[nelem]);
00190 fprintf( ioOut, " R " );
00191 for( i=11; i < limit; i++ )
00192 {
00193 fprintf( ioOut, "%10.2e", dense.eden*ionbal.CotaRate[ion] );
00194 }
00195 fprintf( ioOut, "\n" );
00196
00197 fprintf( ioOut, " " );
00198 for( i=11; i < limit; i++ )
00199 {
00200 fprintf( ioOut, "%10.2e", ionbal.RateRecomTot[nelem][i] );
00201 }
00202 fprintf( ioOut, "\n\n" );
00203 }
00204 }
00205
00206
00207
00208 limit = MIN2(nelem-NISO,dense.IonHigh[nelem]-1);
00209 for( i=dense.IonLow[nelem]; i <= limit; i++ )
00210 {
00211 ASSERT( Heavy.xLyaHeavy[nelem][i] > 0. );
00212 ASSERT( ionbal.RateRecomTot[nelem][i] > 0. );
00213 }
00214 return;
00215 }
00216
00217
00218 void ion_recombAGN( FILE * io )
00219 {
00220 # define N1LIM 3
00221 # define N2LIM 4
00222 double te1[N1LIM]={ 5000., 10000., 20000.};
00223 double te2[N2LIM]={ 20000.,50000.,100000.,1e6};
00224
00225 double BreakEnergy = 100./13.0;
00226 long int nelem, ion , i;
00227
00228 char chString[100],
00229 chOutput[100];
00230
00231 double TempSave = phycon.te;
00232
00233 double EdenSave = dense.eden;
00234
00235 DEBUG_ENTRY( "ion_recomb(false,)" );
00236
00237 EdenChange( 1. );
00238
00239
00240
00241 fprintf(io,"X+i\\Te");
00242 for( i=0; i<N1LIM; ++i )
00243 {
00244 phycon.te = te1[i];
00245 fprintf(io,"\t%.0f K",phycon.te);
00246 }
00247 fprintf(io,"\n");
00248
00249
00250 for( nelem=ipLITHIUM; nelem<LIMELM; ++nelem )
00251 {
00252
00253 if( abund.lgAGN[nelem] )
00254 {
00255 for( ion=0; ion<=nelem; ++ion )
00256 {
00257 ASSERT( Heavy.Valence_IP_Ryd[nelem][ion] > 0.05 );
00258
00259 if( Heavy.Valence_IP_Ryd[nelem][ion] > BreakEnergy )
00260 break;
00261
00262
00263 sprintf(chOutput,"%s",
00264 elementnames.chElementSym[nelem]);
00265
00266 if( chOutput[1]==' ' )
00267 chOutput[1] = chOutput[2];
00268
00269 if( ion==0 )
00270 {
00271 sprintf(chString,"0 ");
00272 }
00273 else if( ion==1 )
00274 {
00275 sprintf(chString,"+ ");
00276 }
00277 else
00278 {
00279 sprintf(chString,"+%li ",ion);
00280 }
00281 strcat( chOutput , chString );
00282 fprintf(io,"%5s",chOutput );
00283
00284 for( i=0; i<N1LIM; ++i )
00285 {
00286 TempChange(te1[i] , false);
00287 dense.IonLow[nelem] = 0;
00288 dense.IonHigh[nelem] = nelem+1;
00289 if( ConvBase(0) )
00290 fprintf(ioQQQ,"PROBLEM ConvBase returned error.\n");
00291 fprintf(io,"\t%.2e",ionbal.RateRecomTot[nelem][ion]);
00292 }
00293 fprintf(io,"\n");
00294 }
00295 fprintf(io,"\n");
00296 }
00297 }
00298
00299
00300 fprintf(io,"X+i\\Te");
00301 for( i=0; i<N2LIM; ++i )
00302 {
00303 TempChange(te2[i] , false);
00304 fprintf(io,"\t%.0f K",phycon.te);
00305 }
00306 fprintf(io,"\n");
00307
00308
00309 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
00310 {
00311
00312 if( abund.lgAGN[nelem] )
00313 {
00314 for( ion=0; ion<=nelem; ++ion )
00315 {
00316 ASSERT( Heavy.Valence_IP_Ryd[nelem][ion] > 0.05 );
00317
00318 if( Heavy.Valence_IP_Ryd[nelem][ion] <= BreakEnergy )
00319 continue;
00320
00321
00322 fprintf(io,"%s",
00323 elementnames.chElementSym[nelem]);
00324
00325 if( ion==0 )
00326 {
00327 fprintf(io,"0 ");
00328 }
00329 else if( ion==1 )
00330 {
00331 fprintf(io,"+ ");
00332 }
00333 else
00334 {
00335 fprintf(io,"+%li",ion);
00336 }
00337
00338 for( i=0; i<N2LIM; ++i )
00339 {
00340 TempChange(te2[i] , false);
00341 dense.IonLow[nelem] = 0;
00342 dense.IonHigh[nelem] = nelem+1;
00343 if( ConvBase(0) )
00344 fprintf(ioQQQ,"PROBLEM ConvBase returned error.\n");
00345 fprintf(io,"\t%.2e",ionbal.RateRecomTot[nelem][ion]);
00346 }
00347 fprintf(io,"\n");
00348 }
00349 fprintf(io,"\n");
00350 }
00351 }
00352
00353 TempChange(TempSave , true);
00354 EdenChange( EdenSave );
00355 return;
00356 }