00001
00002
00003
00004
00005
00006 #include "cddefines.h"
00007 #include "atmdat.h"
00008 #include "conv.h"
00009 #include "dense.h"
00010 #include "opacity.h"
00011 #include "elementnames.h"
00012 #include "h2.h"
00013 #include "helike.h"
00014 #include "helike_cs.h"
00015 #include "hmi.h"
00016 #include "mole.h"
00017 #include "hydrogenic.h"
00018 #include "ionbal.h"
00019 #include "iso.h"
00020 #include "phycon.h"
00021 #include "rfield.h"
00022 #include "secondaries.h"
00023 #include "taulines.h"
00024 #include "thermal.h"
00025 #include "trace.h"
00026
00027 void iso_collapsed_update( void )
00028 {
00029 for( long nelem=ipHYDROGEN; nelem<=ipHELIUM; nelem++ )
00030 {
00031 if (dense.lgElmtOn[nelem])
00032 {
00033 for( long ipISO=ipH_LIKE; ipISO<MIN2(NISO,nelem+1); ipISO++ )
00034 {
00035 if ((dense.IonHigh[nelem] >= nelem - ipISO &&
00036 dense.IonLow[nelem] <= nelem - ipISO) || !conv.nTotalIoniz)
00037 {
00038 iso_collapsed_bnl_set( ipISO, nelem );
00039
00040 iso_collapsed_Aul_update( ipISO, nelem );
00041
00042 iso_collapsed_lifetimes_update( ipISO, nelem );
00043
00044 iso_cascade( ipISO, nelem );
00045 }
00046 }
00047 }
00048 }
00049 }
00050
00051 void iso_update_rates( void )
00052 {
00053 for( long nelem=ipHYDROGEN; nelem<LIMELM; nelem++ )
00054 {
00055 if (!dense.lgElmtOn[nelem])
00056 continue;
00057 for( long ipISO=ipH_LIKE; ipISO<MIN2(NISO,nelem+1); ipISO++ )
00058 {
00059 if ((dense.IonHigh[nelem] >= nelem - ipISO &&
00060 dense.IonLow[nelem] <= nelem - ipISO) || !conv.nTotalIoniz)
00061 {
00062
00063
00064 iso_collide( ipISO, nelem );
00065
00066
00067
00068 if( iso_ctrl.lgContinuumLoweringEnabled[ipISO] && !conv.nPres2Ioniz )
00069 iso_continuum_lower( ipISO, nelem );
00070
00071
00072 iso_radiative_recomb( ipISO , nelem );
00073
00074
00075 iso_photo( ipISO , nelem );
00076
00077
00078 if( iso_ctrl.lgRandErrGen[ipISO] && nzone==0 && !iso_sp[ipISO][nelem].lgErrGenDone )
00079 {
00080 iso_error_generation(ipISO, nelem );
00081 }
00082
00083 iso_radiative_recomb_effective( ipISO, nelem );
00084
00085 iso_ionize_recombine( ipISO , nelem );
00086
00087 ionbal.RateRecomTot[nelem][nelem-ipISO] = ionbal.RateRecomIso[nelem][ipISO];
00088 }
00089
00091 ASSERT( ipISO <= ipHE_LIKE );
00092
00093 t_iso_sp* sp = &iso_sp[ipISO][nelem];
00094 for( vector<two_photon>::iterator tnu = sp->TwoNu.begin(); tnu != sp->TwoNu.end(); ++tnu )
00095 {
00096 CalcTwoPhotonRates( *tnu, rfield.lgInducProcess && iso_ctrl.lgInd2nu_On );
00097 }
00098 }
00099 }
00100 }
00101
00102 void iso_solve(long ipISO, long nelem, double &maxerr)
00103 {
00104 DEBUG_ENTRY( "iso_solve()" );
00105
00106 maxerr = 0.;
00107
00108 if( dense.lgElmtOn[nelem] )
00109 {
00110
00111
00112 if( (dense.IonHigh[nelem] >= nelem - ipISO) &&
00113 (dense.IonLow[nelem] <= nelem - ipISO) )
00114 {
00115
00116 double renorm;
00117 iso_level( ipISO , nelem, renorm );
00118 if (fabs(renorm-1.0) > maxerr)
00119 maxerr = fabs(renorm-1.0);
00120
00121
00122 if( ipISO == ipH_LIKE )
00123 HydroLevel(nelem);
00124 }
00125 else
00126 {
00127
00128 iso_sp[ipISO][nelem].st[0].Pop() = 0.;
00129 for( long ipHi=1; ipHi < iso_sp[ipISO][nelem].numLevels_max; ipHi++ )
00130 {
00131 iso_sp[ipISO][nelem].st[ipHi].Pop() = 0.;
00132 for( long ipLo=0; ipLo < ipHi; ipLo++ )
00133 {
00134 if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() <= iso_ctrl.SmallA )
00135 continue;
00136
00137
00138 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().PopOpc() = 0.;
00139 }
00140 }
00141 }
00142
00143 ASSERT( (*iso_sp[ipISO][nelem].trans(iso_ctrl.nLyaLevel[ipISO],0).Lo()).Pop() == iso_sp[ipISO][nelem].st[0].Pop() );
00144 }
00145
00146 return;
00147 }
00148
00149 void IonHydro( void )
00150 {
00151 DEBUG_ENTRY( "IonHydro()" );
00152
00153
00154
00155
00156 {
00157
00158
00159
00160
00161 enum {DEBUG_LOC=false};
00162
00163 if(DEBUG_LOC )
00164 {
00165 fprintf(ioQQQ,"DEBUG \t%.2f\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
00166 fnzone,
00167 hmi.H2_total ,
00168 findspecieslocal("H2")->den,
00169 dense.xIonDense[ipHYDROGEN][0],
00170 dense.xIonDense[ipHYDROGEN][1],
00171 hmi.H2_Solomon_dissoc_rate_used_H2g,
00172 hmi.H2_Solomon_dissoc_rate_BD96_H2g,
00173 hmi.H2_Solomon_dissoc_rate_TH85_H2g);
00174 }
00175 }
00176 #if 0
00177
00178 if( dense.lgSetIoniz[ipHYDROGEN] )
00179 {
00180 realnum dense_ation = dense.xIonDense[ipHYDROGEN][0]+dense.xIonDense[ipHYDROGEN][1];
00181 dense.xIonDense[ipHYDROGEN][1] = dense.SetIoniz[ipHYDROGEN][1]*dense_ation;
00182 dense.xIonDense[ipHYDROGEN][0] = dense.SetIoniz[ipHYDROGEN][0]*dense_ation;
00183
00184
00185 iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop() = dense.xIonDense[ipHYDROGEN][0];
00186 }
00187 else
00188 {
00189
00190
00191
00192
00193
00194
00195 ion_solver( ipHYDROGEN , false );
00196 }
00197
00198 fixit();
00199
00200
00201 double renorm;
00202 iso_renorm( ipHYDROGEN, ipH_LIKE, renorm );
00203 #endif
00204 ion_solver( ipHYDROGEN , false );
00205
00206
00207
00208
00209
00210 if( iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH2p].Pop()/MAX2(SMALLDOUBLE,iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop()) > 0.1 &&
00211 iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop() > SMALLDOUBLE )
00212 {
00213 hydro.lgHiPop2 = true;
00214 hydro.pop2mx = (realnum)MAX2(iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH2p].Pop()/iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop(),
00215 hydro.pop2mx);
00216 }
00217
00218 double gamtot = iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].gamnc + secondaries.csupra[ipHYDROGEN][0];
00219
00220 double coltot = iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ColIoniz +
00221 iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Coll().ColUL( colliders ) / dense.EdenHCorr *
00222 4. * iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH2p].Boltzmann();
00223
00224
00225 if( iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].RateLevel2Cont > SMALLFLOAT )
00226 {
00227 hydro.H_ion_frac_photo =
00228 (realnum)(iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].gamnc/iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].RateLevel2Cont );
00229
00230
00231 hydro.H_ion_frac_collis =
00232 (realnum)(iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ColIoniz*dense.eden/iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].RateLevel2Cont);
00233
00234
00235
00236 secondaries.sec2total =
00237 (realnum)(secondaries.csupra[ipHYDROGEN][0] / iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].RateLevel2Cont);
00238
00239
00240 atmdat.HIonFrac = atmdat.HCharExcIonTotal / iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].RateLevel2Cont;
00241 }
00242 else
00243 {
00244 hydro.H_ion_frac_collis = 0.;
00245 hydro.H_ion_frac_photo = 0.;
00246 secondaries.sec2total = 0.;
00247 atmdat.HIonFrac = 0.;
00248 }
00249
00250 if( trace.lgTrace )
00251 {
00252 fprintf( ioQQQ, " Hydrogenic return %.2f ",fnzone);
00253 fprintf(ioQQQ,"H0:%.3e ", dense.xIonDense[ipHYDROGEN][0]);
00254 fprintf(ioQQQ,"H+:%.3e ", dense.xIonDense[ipHYDROGEN][1]);
00255 fprintf(ioQQQ,"H2:%.3e ", hmi.H2_total);
00256 fprintf(ioQQQ,"H-:%.3e ", findspecieslocal("H-")->den);
00257 fprintf(ioQQQ,"ne:%.3e ", dense.eden);
00258 fprintf( ioQQQ, " REC, COL, GAMT= ");
00259
00260 fprintf(ioQQQ,"%.2e ", iso_sp[ipH_LIKE][ipHYDROGEN].RadRec_effec );
00261 fprintf(ioQQQ,"%.2e ", coltot);
00262 fprintf(ioQQQ,"%.2e ", gamtot);
00263 fprintf( ioQQQ, " CSUP=");
00264 PrintE82( ioQQQ, secondaries.csupra[ipHYDROGEN][0]);
00265 fprintf( ioQQQ, "\n");
00266 }
00267
00268 return;
00269 }
00270
00271
00272 void iso_renorm( long nelem, long ipISO, double &renorm )
00273 {
00274 DEBUG_ENTRY( "iso_renorm()" );
00275
00276 const bool lgJustAssert = false, lgJustTest = true;
00277 double sum_atom_iso;
00278
00279 renorm = 1.0;
00280
00281 if (ipISO > nelem)
00282 return;
00283
00284
00285
00286
00287
00288 sum_atom_iso = 0.;
00289
00290 for( long level=0; level < iso_sp[ipISO][nelem].numLevels_local; level++ )
00291 {
00292
00293 sum_atom_iso += iso_sp[ipISO][nelem].st[level].Pop();
00294 }
00295
00296
00297
00298 if (sum_atom_iso <= SMALLFLOAT)
00299 {
00300
00301 if (dense.xIonDense[nelem][nelem-ipISO] > 2.0*SMALLFLOAT)
00302 sum_atom_iso = 0.5*dense.xIonDense[nelem][nelem-ipISO];
00303 else
00304 sum_atom_iso = 1.0;
00305 if ( !lgJustTest )
00306 iso_sp[ipISO][nelem].st[0].Pop() = sum_atom_iso;
00307 }
00308
00309
00310 renorm = dense.xIonDense[nelem][nelem-ipISO] / sum_atom_iso;
00311
00312 if (lgJustTest)
00313 return;
00314
00315 if (0)
00316 fprintf (ioQQQ, "Iso_Renorm %ld %ld %g %g %g[%g]\n",
00317 ipISO,nelem,renorm-1.,
00318 dense.xIonDense[nelem][nelem-ipISO], sum_atom_iso,
00319 iso_sp[ipISO][nelem].st[0].Pop());
00320 if (lgJustAssert)
00321 {
00322 ASSERT(fp_equal(renorm,1.0));
00323 return;
00324 }
00325 if (conv.lgConvIoniz() && dense.xIonDense[nelem][nelem-ipISO] > SMALLFLOAT &&
00326 fabs(renorm - 1.0) > conv.IonizErrorAllowed)
00327 {
00328 conv.setConvIonizFail( "Iso vs. ion",
00329 dense.xIonDense[nelem][nelem-ipISO],
00330 sum_atom_iso);
00331
00332 }
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342 for( long ipHi=0; ipHi < iso_sp[ipISO][nelem].numLevels_local; ipHi++ )
00343 {
00344 iso_sp[ipISO][nelem].st[ipHi].Pop() *= renorm;
00345 }
00346
00347
00348 for( long ipHi=0; ipHi < iso_sp[ipISO][nelem].numLevels_local; ipHi++ )
00349 {
00350 for( long ipLo=0; ipLo < ipHi; ipLo++ )
00351 {
00352 if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() <= iso_ctrl.SmallA )
00353 continue;
00354
00355
00356 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().PopOpc() *= renorm;
00357 }
00358 }
00359
00360 return;
00361 }
00362
00363 void iso_departure_coefficients( long ipISO, long nelem )
00364 {
00365 DEBUG_ENTRY( "iso_departure_coefficients()" );
00366
00367 for( long level=0; level < iso_sp[ipISO][nelem].numLevels_local; level++ )
00368 {
00369 double denom = dense.xIonDense[nelem][nelem+1-ipISO]*
00370 iso_sp[ipISO][nelem].fb[level].PopLTE*dense.eden;
00371
00372 if( iso_sp[ipISO][nelem].fb[level].PopLTE > 0. && denom > SMALLFLOAT )
00373 iso_sp[ipISO][nelem].fb[level].DepartCoef = safe_div(
00374 iso_sp[ipISO][nelem].st[level].Pop(), denom );
00375 else
00376 iso_sp[ipISO][nelem].fb[level].DepartCoef = 0.;
00377 }
00378
00379 for( long level=iso_sp[ipISO][nelem].numLevels_local; level < iso_sp[ipISO][nelem].numLevels_max; level++ )
00380 iso_sp[ipISO][nelem].fb[level].DepartCoef = 0.;
00381
00382 return;
00383 }
00384
00385
00386 void iso_prt_pops( long ipISO, long nelem, bool lgPrtDeparCoef )
00387 {
00388 long int in, il, is, i, ipLo, nResolved, ipFirstCollapsed=LONG_MIN;
00389 char chPrtType[2][12]={"populations","departure"};
00390
00391 char chSpin[3][9]= {"singlets", "doublets", "triplets"};
00392
00393 #define ITEM_TO_PRINT(A_) ( lgPrtDeparCoef ? iso_sp[ipISO][nelem].fb[A_].DepartCoef : iso_sp[ipISO][nelem].st[A_].Pop() )
00394
00395 DEBUG_ENTRY( "iso_prt_pops()" );
00396
00397 ASSERT( ipISO < NISO );
00398
00399 for( is = 1; is<=3; ++is)
00400 {
00401 if( ipISO == ipH_LIKE && is != 2 )
00402 continue;
00403 else if( ipISO == ipHE_LIKE && is != 1 && is != 3 )
00404 continue;
00405
00406 ipFirstCollapsed= iso_sp[ipISO][nelem].numLevels_local-iso_sp[ipISO][nelem].nCollapsed_local;
00407 nResolved = iso_sp[ipISO][nelem].st[ipFirstCollapsed-1].n();
00408 ASSERT( nResolved == iso_sp[ipISO][nelem].n_HighestResolved_local );
00409 ASSERT(nResolved > 0 );
00410
00411
00412 fprintf(ioQQQ," %s %s %s %s\n",
00413 iso_ctrl.chISO[ipISO],
00414 elementnames.chElementSym[nelem],
00415 chSpin[is-1],
00416 chPrtType[lgPrtDeparCoef]);
00417
00418
00419 fprintf(ioQQQ," n\\l=> ");
00420 for( i =0; i < nResolved; ++i)
00421 {
00422 fprintf(ioQQQ,"%2ld ",i);
00423 }
00424 fprintf(ioQQQ,"\n");
00425
00426
00427 for( in = 1; in <= nResolved; ++in)
00428 {
00429 if( is==3 && in==1 )
00430 continue;
00431
00432 fprintf(ioQQQ," %2ld ",in);
00433
00434 for( il = 0; il < in; ++il)
00435 {
00436 if( ipISO==ipHE_LIKE && (in==2) && (il==1) && (is==3) )
00437 {
00438 fprintf( ioQQQ, "%9.3e ", ITEM_TO_PRINT(ipHe2p3P0) );
00439 fprintf( ioQQQ, "%9.3e ", ITEM_TO_PRINT(ipHe2p3P1) );
00440 fprintf( ioQQQ, "%9.3e ", ITEM_TO_PRINT(ipHe2p3P2) );
00441 }
00442 else
00443 {
00444 ipLo = iso_sp[ipISO][nelem].QuantumNumbers2Index[in][il][is];
00445 fprintf( ioQQQ, "%9.3e ", ITEM_TO_PRINT(ipLo) );
00446 }
00447 }
00448 fprintf(ioQQQ,"\n");
00449 }
00450 }
00451
00452 for( il = ipFirstCollapsed; il < iso_sp[ipISO][nelem].numLevels_local; ++il)
00453 {
00454 in = iso_sp[ipISO][nelem].st[il].n();
00455
00456 fprintf(ioQQQ," %2ld ",in);
00457 fprintf( ioQQQ, "%9.3e ", ITEM_TO_PRINT(il) );
00458 fprintf(ioQQQ,"\n");
00459 }
00460 return;
00461 }
00462
00463
00464 void AGN_He1_CS( FILE *ioPun )
00465 {
00466
00467 long int i;
00468
00469
00470 const int NTE = 5;
00471 double TeList[NTE] = {6000.,10000.,15000.,20000.,25000.};
00472 double TempSave;
00473
00474 DEBUG_ENTRY( "AGN_He1_CS()" );
00475
00476
00477 fprintf(ioPun, "Te\t2 3s 33s\n");
00478
00479
00480 TempSave = phycon.te;
00481
00482 for( i=0; i<NTE; ++i )
00483 {
00484 TempChange(TeList[i] , false);
00485
00486 fprintf(ioPun , "%.0f\t",
00487 TeList[i] );
00488 fprintf(ioPun , "%.2f\t",
00489 HeCSInterp( 1 , ipHe3s3S , ipHe2s3S, ipELECTRON ) );
00490 fprintf(ioPun , "%.2f\t",
00491 HeCSInterp( 1 , ipHe3p3P , ipHe2s3S, ipELECTRON ) );
00492 fprintf(ioPun , "%.2f\t",
00493 HeCSInterp( 1 , ipHe3d3D , ipHe2s3S, ipELECTRON ) );
00494 fprintf(ioPun , "%.3f\t",
00495 HeCSInterp( 1 , ipHe3d1D , ipHe2s3S, ipELECTRON ) );
00496
00497 fprintf(ioPun , "%.1f\n",
00498 HeCSInterp( 1 , ipHe2p3P0 , ipHe2s3S, ipELECTRON ) +
00499 HeCSInterp( 1 , ipHe2p3P1 , ipHe2s3S, ipELECTRON ) +
00500 HeCSInterp( 1 , ipHe2p3P2 , ipHe2s3S, ipELECTRON ));
00501 }
00502
00503
00504 TempChange(TempSave , false);
00505 return;
00506 }