00001
00002
00003
00004 #include "cddefines.h"
00005 #include "hydrogenic.h"
00006 #include "iso.h"
00007 #include "opacity.h"
00008 #include "trace.h"
00009 #include "phycon.h"
00010 #include "ionbal.h"
00011 #include "punch.h"
00012 #include "elementnames.h"
00013 #include "atmdat.h"
00014 #include "rt.h"
00015
00016 void HydroRecom(
00017
00018 long int nelem)
00019 {
00020 long int n;
00021 bool lgOK;
00022 double extra,
00023 SumCaseB ,
00024 SumTopOff;
00025
00026 static double TeUsed[LIMELM]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
00027 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
00028
00029 DEBUG_ENTRY( "HydroRecom()" );
00030
00031
00032 ASSERT( nelem >= 0);
00033 ASSERT( nelem < LIMELM );
00034
00035
00036
00037 for( n=ipH1s; n < iso.numLevels_local[ipH_LIKE][nelem]; n++ )
00038 {
00039
00040 iso.RadRecomb[ipH_LIKE][nelem][n][ipRecEsc] =
00041 RT_recom_effic(iso.ipIsoLevNIonCon[ipH_LIKE][nelem][n]);
00042
00043
00044
00045 iso.RadRecomb[ipH_LIKE][nelem][n][ipRecEsc] =
00046 MAX2(iso.RadRecomb[ipH_LIKE][nelem][n][ipRecEsc], opac.otsmin);
00047 }
00048
00049
00050 for( n=iso.numLevels_local[ipH_LIKE][nelem]; n < iso.numLevels_max[ipH_LIKE][nelem]; n++ )
00051 {
00052 iso.RadRecomb[ipH_LIKE][nelem][n][ipRecEsc] = 0.;
00053 }
00054
00055 # if 0
00056 n = iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][3];
00057 if( nelem==0 ) fprintf(ioQQQ,"%i\thydrogrecom optical depths\t%.2e\t%.2e\t%.2e \n",
00058 nzone,
00059 opac.TauAbsGeo[0][n-1],
00060 opac.TauAbsGeo[1][n-1],
00061 opac.TauAbsGeo[1][n-1] - opac.TauAbsGeo[0][n-1]);
00062 # endif
00063
00064
00065 if( opac.lgCaseB )
00066 {
00067 iso.RadRecomb[ipH_LIKE][nelem][ipH1s][ipRecEsc] = 1e-10;
00068 }
00069
00070
00071
00072 for( n=ipH1s; n < iso.numLevels_local[ipH_LIKE][nelem]; n++ )
00073 {
00074
00075
00076
00077 iso.RadRecomb[ipH_LIKE][nelem][n][ipRecNetEsc] =
00078 MIN2(1.,iso.RadRecomb[ipH_LIKE][nelem][n][ipRecEsc]+
00079 (1.-iso.RadRecomb[ipH_LIKE][nelem][n][ipRecEsc])*iso.ConOpacRatio[ipH_LIKE][nelem][n]);
00080
00081 ASSERT( iso.RadRecomb[ipH_LIKE][nelem][n][ipRecNetEsc] >= 0. &&
00082 iso.RadRecomb[ipH_LIKE][nelem][n][ipRecNetEsc] <= 1. );
00083 }
00084
00085
00086 for( n=iso.numLevels_local[ipH_LIKE][nelem]; n < iso.numLevels_max[ipH_LIKE][nelem]; n++ )
00087 {
00088 iso.RadRecomb[ipH_LIKE][nelem][n][ipRecNetEsc] = 0.;
00089 }
00090
00091
00092
00093
00094 if( (trace.lgTrace && trace.lgIsoTraceFull[ipH_LIKE]) && (nelem == trace.ipIsoTrace[ipH_LIKE]) )
00095 {
00096
00097 fprintf( ioQQQ, " HydroRecom recomb effic Z%li ",nelem );
00098 for( n=ipH1s; n < iso.numLevels_max[ipH_LIKE][nelem]; n++ )
00099 {
00100 fprintf( ioQQQ,PrintEfmt("%10.3e", iso.RadRecomb[ipH_LIKE][nelem][n][ipRecEsc] ));
00101 }
00102 fprintf( ioQQQ, "\n" );
00103
00104
00105 fprintf( ioQQQ, " HydroRecom recomb net effic" );
00106
00107 for( n=ipH1s; n < iso.numLevels_local[ipH_LIKE][nelem]; n++ )
00108 {
00109 fprintf( ioQQQ,PrintEfmt("%10.3e", iso.RadRecomb[ipH_LIKE][nelem][n][ipRecNetEsc]) );
00110 }
00111 fprintf( ioQQQ, "\n" );
00112
00113
00114 fprintf( ioQQQ, " HydroRecom in optic dep%10.3e",
00115 opac.TauAbsGeo[0][iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH1s]-1] );
00116
00117 for( n=2; n < iso.numLevels_local[ipH_LIKE][nelem]; n++ )
00118 {
00119 fprintf( ioQQQ,PrintEfmt("%10.3e", opac.TauAbsGeo[0][iso.ipIsoLevNIonCon[ipH_LIKE][nelem][n]-1] ));
00120 }
00121 fprintf( ioQQQ, "\n" );
00122
00123
00124 fprintf( ioQQQ, " HydroRecom out op depth%10.3e",
00125 opac.TauAbsGeo[1][iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH1s]-1] );
00126
00127 for( n=2; n < iso.numLevels_local[ipH_LIKE][nelem]; n++ )
00128 {
00129 fprintf( ioQQQ,PrintEfmt("%10.3e", opac.TauAbsGeo[1][iso.ipIsoLevNIonCon[ipH_LIKE][nelem][n]-1] ));
00130 }
00131 fprintf( ioQQQ, "\n" );
00132 }
00133
00134 if( (fabs(1.-TeUsed[nelem]/phycon.te)> 0.001) || hydro.lgReevalRecom )
00135 {
00136 TeUsed[nelem] = phycon.te;
00137
00138
00139 iso.RadRecomb[ipH_LIKE][nelem][ipH1s][ipRecRad] = t_ADfA::Inst().H_rad_rec(nelem+1,ipH1s,phycon.te);
00140
00141
00142
00143 SumCaseB = 0.;
00144
00145 for( n=ipH2s; n < iso.numLevels_local[ipH_LIKE][nelem]; n++ )
00146 {
00147 iso.RadRecomb[ipH_LIKE][nelem][n][ipRecRad] = t_ADfA::Inst().H_rad_rec(nelem+1,n,phycon.te);
00148 ASSERT( iso.RadRecomb[ipH_LIKE][nelem][n][ipRecRad] > 0. );
00149 SumCaseB += iso.RadRecomb[ipH_LIKE][nelem][n][ipRecRad];
00150 }
00151
00152
00153 for( n=iso.numLevels_local[ipH_LIKE][nelem]; n < iso.numLevels_max[ipH_LIKE][nelem]; n++ )
00154 {
00155 iso.RadRecomb[ipH_LIKE][nelem][n][ipRecRad] = 0.;
00156 }
00157
00158
00159 iso.RadRec_caseB[ipH_LIKE][nelem] = t_ADfA::Inst().H_rad_rec(nelem+1,-1,phycon.te);
00160
00161
00162 ionbal.RR_rate_coef_used[nelem][nelem] = iso.RadRec_caseB[ipH_LIKE][nelem];
00163
00164
00165
00166
00167
00168
00169 iso.RadRec_caseB[ipH_LIKE][nelem] -= iso.RadRecomb[ipH_LIKE][nelem][ipH1s][ipRecRad];
00170
00171
00172 extra = iso.RadRec_caseB[ipH_LIKE][nelem] - SumCaseB;
00173
00174
00175 ASSERT( (iso.numLevels_max[ipH_LIKE][nelem] > iso.nTopOff[ipH_LIKE][nelem]) ||
00176 iso.lgLevelsLowered[ipH_LIKE][nelem] );
00177
00178
00179
00180
00181
00182
00183
00184 if( extra > 0. && !iso.lgLevelsLowered[ipH_LIKE][nelem] )
00185 {
00186 if( strcmp(hydro.chHTopType,"scal") == 0 )
00187 {
00188
00189
00190 SumTopOff = 0;
00191 for( n=iso.nTopOff[ipH_LIKE][nelem]; n < iso.numLevels_max[ipH_LIKE][nelem]; n++ )
00192 {
00193 SumTopOff += iso.RadRecomb[ipH_LIKE][nelem][n][ipRecRad];
00194 }
00195
00196
00197 if( SumTopOff > 1e-30 )
00198 {
00199 extra = (1. + extra/SumTopOff);
00200 }
00201 else
00202 {
00203 extra = 1.;
00204 }
00205
00206
00207 for( n=iso.nTopOff[ipH_LIKE][nelem]; n < iso.numLevels_max[ipH_LIKE][nelem]; n++ )
00208 {
00209 iso.RadRecomb[ipH_LIKE][nelem][n][ipRecRad] *= extra;
00210 }
00211 }
00212 else if( strcmp(hydro.chHTopType," add") == 0 )
00213 {
00214
00215
00216 extra /= (iso.numLevels_max[ipH_LIKE][nelem] - iso.nTopOff[ipH_LIKE][nelem]);
00217 ASSERT( extra >= 0. );
00218 for( n=iso.nTopOff[ipH_LIKE][nelem]; n < iso.numLevels_max[ipH_LIKE][nelem]; n++ )
00219 {
00220 iso.RadRecomb[ipH_LIKE][nelem][n][ipRecRad] += extra;
00221 }
00222 }
00223 else
00224 {
00225
00226 fprintf( ioQQQ, " HydroRecom has insane chHTopType=%4.4s\n",
00227 hydro.chHTopType );
00228 ShowMe();
00229 puts( "[Stop in hydrorecom]" );
00230 cdEXIT(EXIT_FAILURE);
00231 }
00232 }
00233
00234
00235 lgOK = true;
00236
00237 for( n=ipH1s; n < iso.numLevels_local[ipH_LIKE][nelem]; n++ )
00238 {
00239 if( iso.RadRecomb[ipH_LIKE][nelem][n][ipRecRad] <= 0. )
00240 {
00241 fprintf( ioQQQ,
00242 " HydroRecom non-positive recombination coefficient for Z=%3ld lev n=%3ld rec=%11.2e te=%11.2e\n",
00243 nelem, n, iso.RadRecomb[ipH_LIKE][nelem][n][ipRecRad] , phycon.te);
00244 lgOK = false;
00245 }
00246 }
00247
00248 if( !lgOK )
00249 {
00250 ShowMe();
00251 puts( "[Stop in HydroRecom]" );
00252 cdEXIT(EXIT_FAILURE);
00253 }
00254
00255
00256 if( (trace.lgTrace && trace.lgIsoTraceFull[ipH_LIKE]) && (nelem == trace.ipIsoTrace[ipH_LIKE]) )
00257 {
00258 fprintf( ioQQQ, " HydroRecom eval rec cof" );
00259 for( n=ipH1s; n < iso.numLevels_max[ipH_LIKE][nelem]; n++ )
00260 {
00261 fprintf( ioQQQ,PrintEfmt("%10.3e", iso.RadRecomb[ipH_LIKE][nelem][n][ipRecRad]) );
00262 }
00263 fprintf( ioQQQ, "\n" );
00264 }
00265 }
00266
00267
00268
00269 ASSERT( iso.RadRecomb[ipH_LIKE][nelem][ipH1s][ipRecRad] > 0. );
00270
00271 ASSERT( iso.RadRecomb[ipH_LIKE][nelem][iso.numLevels_local[ipH_LIKE][nelem]-1][ipRecRad] > 0. );
00272
00273
00274 iso.RadRec_effec[ipH_LIKE][nelem] = 0.;
00275
00276 for( n=ipH1s; n < iso.numLevels_local[ipH_LIKE][nelem]; n++ )
00277 {
00278 iso.RadRec_effec[ipH_LIKE][nelem] += iso.RadRecomb[ipH_LIKE][nelem][n][ipRecRad]*
00279 iso.RadRecomb[ipH_LIKE][nelem][n][ipRecNetEsc];
00280 }
00281
00282
00283 if( punch.lgioRecom )
00284 {
00285 if( nelem==0 )
00286 {
00287 fprintf( punch.ioRecom, "Hydrogenic species\n");
00288 }
00289
00290 fprintf( punch.ioRecom, "H-like %2li %s ",
00291 nelem+1 , elementnames.chElementSym[nelem] );
00292 fprintf( punch.ioRecom,PrintEfmt("%9.2e", iso.RadRec_caseB[ipH_LIKE][nelem] ));
00293 fprintf( punch.ioRecom, "\n" );
00294 }
00295
00296 if( trace.lgTrace && trace.lgHBug )
00297 {
00298 fprintf( ioQQQ, " HydroRecom Z=%3ld total rec coef", nelem );
00299 fprintf( ioQQQ,PrintEfmt("%10.3e", iso.RadRec_effec[ipH_LIKE][nelem] ));
00300 fprintf( ioQQQ, " case A=" );
00301 fprintf( ioQQQ,PrintEfmt("%10.3e",
00302 iso.RadRec_caseB[ipH_LIKE][nelem] + iso.RadRecomb[ipH_LIKE][nelem][ipH1s][ipRecRad] ) );
00303 fprintf( ioQQQ, " caseB=");
00304 fprintf( ioQQQ,PrintEfmt("%10.3e", iso.RadRec_caseB[ipH_LIKE][nelem] ));
00305 fprintf( ioQQQ, "\n" );
00306 }
00307
00308 DEBUG_EXIT( "HydroRecom()" );
00309 return;
00310 }
00311
00312
00313