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 "hydrogenic.h"
00017 #include "ionbal.h"
00018 #include "iso.h"
00019 #include "phycon.h"
00020 #include "secondaries.h"
00021 #include "taulines.h"
00022 #include "thermal.h"
00023 #include "trace.h"
00024
00025 valarray<double> PumpSave;
00026
00027 void iso_drive( void )
00028 {
00029 DEBUG_ENTRY( "iso_drive()" );
00030
00031 iso_charge_transfer_update();
00032
00033
00034 if( !conv.nTotalIoniz )
00035 {
00036 for( long ipISO=ipH_LIKE; ipISO<NISO; ipISO++ )
00037 for( long nelem=ipISO; nelem<LIMELM; nelem++ )
00038 {
00039 iso_update_rates( ipISO, nelem );
00040 }
00041 }
00042 for( long ipISO=ipH_LIKE; ipISO<NISO; ipISO++ )
00043 for( long nelem=ipISO; nelem<LIMELM; nelem++ )
00044 {
00045 iso_update_rates( ipISO, nelem );
00046 iso_solve( ipISO, nelem );
00047 }
00048
00049 return;
00050 }
00051
00052
00053 void iso_update_rates(long ipISO, long nelem )
00054 {
00055 DEBUG_ENTRY( "iso_update_rates()" );
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065 if( ipISO==ipH_LIKE && nelem==ipHYDROGEN && hydro.xLymanPumpingScaleFactor!=1.f )
00066 {
00067 if( PumpSave.size() != (unsigned int)iso.numLevels_local[ipISO][nelem] )
00068 PumpSave.resize( iso.numLevels_local[ipISO][nelem] );
00069 for( long ipHi=1; ipHi < iso.numLevels_local[ipH_LIKE][ipHYDROGEN]; ++ipHi )
00070 {
00071 PumpSave[ipHi] = Transitions[ipISO][ipHYDROGEN][ipHi][0].Emis->pump;
00072 Transitions[ipISO][ipHYDROGEN][ipHi][0].Emis->pump *= hydro.xLymanPumpingScaleFactor;
00073 }
00074 }
00075
00076
00077 if( dense.lgElmtOn[nelem] )
00078 {
00079
00080
00081 if( ((dense.IonHigh[nelem] >= nelem - ipISO) &&
00082 (dense.IonLow[nelem] <= nelem - ipISO)) || !conv.nTotalIoniz )
00083 {
00084
00085
00086 if( iso.lgContinuumLoweringEnabled[ipISO] )
00087 iso_continuum_lower( ipISO, nelem );
00088
00089
00090 iso_collide( ipISO, nelem );
00091
00092
00093 iso_photo( ipISO , nelem );
00094
00095
00096 iso_radiative_recomb( ipISO , nelem );
00097
00098
00099 if( iso.lgRandErrGen[ipISO] && nzone==0 && !iso.lgErrGenDone[ipISO][nelem] )
00100 {
00101 iso_error_generation(ipISO, nelem );
00102 }
00103
00104 if( opac.lgRedoStatic )
00105 {
00106 if( nelem<=ipHELIUM )
00107 {
00108 iso_collapsed_bnl_set( ipISO, nelem );
00109
00110
00111
00112 iso_collapsed_Aul_update( ipISO, nelem );
00113
00114 iso_collapsed_lifetimes_update( ipISO, nelem );
00115 }
00116
00117 iso_cascade( ipISO, nelem );
00118
00119 iso_radiative_recomb_effective( ipISO, nelem );
00120 }
00121
00122
00123
00124 iso_ionize_recombine( ipISO , nelem );
00125 }
00126 }
00127
00128 return;
00129 }
00130
00131 void iso_solve(long ipISO, long nelem)
00132 {
00133 DEBUG_ENTRY( "iso_solve()" );
00134
00135
00136 if( dense.lgElmtOn[nelem] )
00137 {
00138
00139 if( dense.lgSetIoniz[nelem] )
00140 {
00141 dense.xIonDense[nelem][nelem+1-ipISO] = dense.SetIoniz[nelem][nelem+1-ipISO]*dense.gas_phase[nelem];
00142 dense.xIonDense[nelem][nelem-ipISO] = dense.SetIoniz[nelem][nelem-ipISO]*dense.gas_phase[nelem];
00143 }
00144
00145
00146
00147 if( (dense.IonHigh[nelem] >= nelem - ipISO) &&
00148 (dense.IonLow[nelem] <= nelem - ipISO) )
00149 {
00150
00151 iso_level( ipISO , nelem );
00152
00153
00154
00155 if( ipISO == ipH_LIKE )
00156 HydroLevel(nelem);
00157 }
00158 else
00159 {
00160 iso.xIonSimple[ipISO][nelem] = 0.;
00161
00162
00163 StatesElemNEW[nelem][nelem-ipISO][0].Pop = 0.;
00164 for( long ipHi=1; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ )
00165 {
00166 StatesElemNEW[nelem][nelem-ipISO][ipHi].Pop = 0.;
00167 for( long ipLo=0; ipLo < ipHi; ipLo++ )
00168 {
00169 if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul <= iso.SmallA )
00170 continue;
00171
00172
00173 Transitions[ipISO][nelem][ipHi][ipLo].Emis->PopOpc = 0.;
00174 }
00175 }
00176 }
00177
00178
00179 iso_departure_coefficients( ipISO, nelem );
00180
00181 ASSERT( Transitions[ipISO][nelem][iso.nLyaLevel[ipISO]][0].Lo->Pop == StatesElemNEW[nelem][nelem-ipISO][0].Pop );
00182 }
00183
00184 return;
00185 }
00186
00187 void IonHydro( void )
00188 {
00189 double coltot,
00190 gamtot,
00191 sum,
00192 error;
00193 bool lgH_chem_conv;
00194 int loop_H_chem;
00195 double solomon_assumed;
00196
00197 DEBUG_ENTRY( "IonHydro()" );
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210 lgH_chem_conv = false;
00211 loop_H_chem = 0;
00212 while( loop_H_chem < 5 && !lgH_chem_conv )
00213 {
00214
00215
00216 if( hmi.H2_rate_destroy < SMALLFLOAT )
00217 solomon_assumed = 1.;
00218 else
00219 solomon_assumed = hmi.H2_Solomon_dissoc_rate_used_H2g/hmi.H2_rate_destroy;
00220
00221
00222 hmole();
00223
00224
00225
00226 H2_RTMake();
00227 H2_LevelPops();
00228
00229
00230
00231 lgH_chem_conv = true;
00232
00233 if( h2.lgH2ON && hmi.lgBigH2_evaluated && hmi.lgH2_Chemistry_BigH2 )
00234 {
00235
00236 if( fabs( solomon_assumed - hmi.H2_Solomon_dissoc_rate_BigH2_H2g/SDIV(hmi.H2_rate_destroy) ) >
00237 conv.EdenErrorAllowed/5.)
00238 {
00239 lgH_chem_conv = false;
00240 }
00241 }
00242 ++loop_H_chem;
00243 }
00244
00245 {
00246
00247
00248
00249
00250 enum {DEBUG_LOC=false};
00251
00252 if(DEBUG_LOC )
00253 {
00254 fprintf(ioQQQ,"DEBUG \t%.2f\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
00255 fnzone,
00256 hmi.H2_total ,
00257 hmi.Hmolec[ipMH2g],
00258 dense.xIonDense[ipHYDROGEN][0],
00259 dense.xIonDense[ipHYDROGEN][1],
00260 hmi.H2_Solomon_dissoc_rate_used_H2g,
00261 hmi.H2_Solomon_dissoc_rate_BD96_H2g,
00262 hmi.H2_Solomon_dissoc_rate_TH85_H2g);
00263 }
00264 }
00265
00266 if( dense.lgSetIoniz[ipHYDROGEN] )
00267 {
00268 dense.xIonDense[ipHYDROGEN][1] = dense.SetIoniz[ipHYDROGEN][1]*dense.gas_phase[ipHYDROGEN];
00269 dense.xIonDense[ipHYDROGEN][0] = dense.SetIoniz[ipHYDROGEN][0]*dense.gas_phase[ipHYDROGEN];
00270
00271
00272 StatesElemNEW[ipHYDROGEN][ipHYDROGEN-ipH_LIKE][ipH1s].Pop = dense.xIonDense[ipHYDROGEN][0];
00273 }
00274 else
00275 {
00276
00277
00278
00279
00280
00281
00282 ion_solver( ipHYDROGEN , false );
00283 }
00284
00285
00286
00287
00288
00289
00290 sum = dense.xIonDense[ipHYDROGEN][0] + dense.xIonDense[ipHYDROGEN][1];
00291 for(long mol=0;mol<N_H_MOLEC;mol++)
00292 {
00293 sum += hmi.Hmolec[mol]*hmi.nProton[mol];
00294 }
00295
00296 sum -= hmi.Hmolec[ipMH]+hmi.Hmolec[ipMHp];
00297
00298
00299
00300
00301
00302 error = ( dense.gas_phase[ipHYDROGEN] - sum )/dense.gas_phase[ipHYDROGEN];
00303
00304 {
00305
00306
00307
00308
00309 enum {DEBUG_LOC=false};
00310
00311 if(DEBUG_LOC && (fabs(error) > 1e-4) )
00312 fprintf(ioQQQ,"PROBLEM hydrogenic zone %li hden %.4e, sum %.4e (h-s)/h %.3e \n", nzone, dense.gas_phase[ipHYDROGEN] , sum ,
00313 error );
00314 }
00315
00316 fixit();
00317
00318
00319 HydroRenorm();
00320
00321
00322
00323
00324
00325 if( StatesElemNEW[ipHYDROGEN][ipHYDROGEN-ipH_LIKE][ipH2p].Pop/MAX2(SMALLDOUBLE,StatesElemNEW[ipHYDROGEN][ipHYDROGEN-ipH_LIKE][ipH1s].Pop) > 0.1 &&
00326 StatesElemNEW[ipHYDROGEN][ipHYDROGEN-ipH_LIKE][ipH1s].Pop > SMALLDOUBLE )
00327 {
00328 hydro.lgHiPop2 = true;
00329 hydro.pop2mx = (realnum)MAX2(StatesElemNEW[ipHYDROGEN][ipHYDROGEN-ipH_LIKE][ipH2p].Pop/StatesElemNEW[ipHYDROGEN][ipHYDROGEN-ipH_LIKE][ipH1s].Pop,
00330 hydro.pop2mx);
00331 }
00332
00333 gamtot = iso.gamnc[ipH_LIKE][ipHYDROGEN][ipH1s] + secondaries.csupra[ipHYDROGEN][0];
00334
00335 coltot = iso.ColIoniz[ipH_LIKE][ipHYDROGEN][ipH1s] +
00336 Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Coll.ColUL*4.*
00337 iso.Boltzmann[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s];
00338
00339
00340 if( iso.RateLevel2Cont[ipH_LIKE][ipHYDROGEN][ipH1s] > SMALLFLOAT )
00341 {
00342 hydro.H_ion_frac_photo =
00343 (realnum)(iso.gamnc[ipH_LIKE][ipHYDROGEN][ipH1s]/iso.RateLevel2Cont[ipH_LIKE][ipHYDROGEN][ipH1s] );
00344
00345
00346 hydro.H_ion_frac_collis =
00347 (realnum)(iso.ColIoniz[ipH_LIKE][ipHYDROGEN][ipH1s]*dense.eden/iso.RateLevel2Cont[ipH_LIKE][ipHYDROGEN][ipH1s]);
00348
00349
00350
00351 secondaries.sec2total =
00352 (realnum)(secondaries.csupra[ipHYDROGEN][0] / iso.RateLevel2Cont[ipH_LIKE][ipHYDROGEN][ipH1s]);
00353
00354
00355 atmdat.HIonFrac = atmdat.HCharExcIonTotal / iso.RateLevel2Cont[ipH_LIKE][ipHYDROGEN][ipH1s];
00356 }
00357 else
00358 {
00359 hydro.H_ion_frac_collis = 0.;
00360 hydro.H_ion_frac_photo = 0.;
00361 secondaries.sec2total = 0.;
00362 atmdat.HIonFrac = 0.;
00363 }
00364
00365 if( trace.lgTrace )
00366 {
00367 fprintf( ioQQQ, " Hydrogenic return %.2f ",fnzone);
00368 fprintf(ioQQQ,"H0:%.3e ", dense.xIonDense[ipHYDROGEN][0]);
00369 fprintf(ioQQQ,"H+:%.3e ", dense.xIonDense[ipHYDROGEN][1]);
00370 fprintf(ioQQQ,"H2:%.3e ", hmi.H2_total);
00371 fprintf(ioQQQ,"H-:%.3e ", hmi.Hmolec[ipMHm]);
00372 fprintf(ioQQQ,"ne:%.3e ", dense.eden);
00373 fprintf( ioQQQ, " REC, COL, GAMT= ");
00374
00375 fprintf(ioQQQ,"%.2e ", iso.RadRec_effec[ipH_LIKE][ipHYDROGEN] );
00376 fprintf(ioQQQ,"%.2e ", coltot);
00377 fprintf(ioQQQ,"%.2e ", gamtot);
00378 fprintf( ioQQQ, " CSUP=");
00379 PrintE82( ioQQQ, secondaries.csupra[ipHYDROGEN][0]);
00380 fprintf( ioQQQ, "\n");
00381 }
00382
00383
00384 if( hydro.xLymanPumpingScaleFactor!=1.f )
00385 {
00386 for( long ipHi=1; ipHi < iso.numLevels_local[ipH_LIKE][ipHYDROGEN]; ++ipHi )
00387 {
00388 Transitions[ipH_LIKE][ipHYDROGEN][ipHi][0].Emis->pump = PumpSave[ipHi];
00389 }
00390 }
00391
00392 return;
00393 }
00394
00395
00396 void HydroRenorm( void )
00397 {
00398 double sum_atom_iso , renorm;
00399
00400 DEBUG_ENTRY( "HydroRenorm()" );
00401
00402
00403
00404
00405
00406 sum_atom_iso = 0.;
00407
00408 for( long level=ipH1s; level < iso.numLevels_local[ipH_LIKE][ipHYDROGEN]; level++ )
00409 {
00410
00411 sum_atom_iso += StatesElemNEW[ipHYDROGEN][ipHYDROGEN-ipH_LIKE][level].Pop;
00412 }
00413
00414
00415 if( sum_atom_iso > SMALLFLOAT )
00416 {
00417 renorm = dense.xIonDense[ipHYDROGEN][0] / sum_atom_iso;
00418 }
00419 else
00420 {
00421 renorm = 0.;
00422 }
00423 ASSERT( renorm < BIGFLOAT );
00424
00425
00426 StatesElemNEW[ipHYDROGEN][ipHYDROGEN-ipH_LIKE][ipH1s].Pop *= renorm;
00427
00428
00429
00430
00431
00432 for( long ipHi=ipH2s; ipHi < iso.numLevels_local[ipH_LIKE][ipHYDROGEN]; ipHi++ )
00433 {
00434 StatesElemNEW[ipHYDROGEN][ipHYDROGEN-ipH_LIKE][ipHi].Pop *= renorm;
00435
00436 for( long ipLo=ipH1s; ipLo < ipHi; ipLo++ )
00437 {
00438 if( Transitions[ipH_LIKE][ipHYDROGEN][ipHi][ipLo].Emis->Aul <= iso.SmallA )
00439 continue;
00440
00441
00442 Transitions[ipH_LIKE][ipHYDROGEN][ipHi][ipLo].Emis->PopOpc *= renorm;
00443 }
00444 }
00445
00446 return;
00447 }
00448
00449 void iso_departure_coefficients( long ipISO, long nelem )
00450 {
00451 DEBUG_ENTRY( "iso_departure_coefficients()" );
00452
00453 for( long level=0; level < iso.numLevels_local[ipISO][nelem]; level++ )
00454 {
00455 double denom = dense.xIonDense[nelem][nelem+1-ipISO]*
00456 iso.PopLTE[ipISO][nelem][level]*dense.eden;
00457
00458 if( iso.PopLTE[ipISO][nelem][level] > 0. && denom > SMALLFLOAT )
00459 iso.DepartCoef[ipISO][nelem][level] = safe_div(
00460 StatesElemNEW[nelem][nelem-ipISO][level].Pop, denom );
00461 else
00462 iso.DepartCoef[ipISO][nelem][level] = 0.;
00463 }
00464
00465 for( long level=iso.numLevels_local[ipISO][nelem]; level < iso.numLevels_max[ipISO][nelem]; level++ )
00466 iso.DepartCoef[ipISO][nelem][level] = 0.;
00467
00468 return;
00469 }
00470
00471
00472 void iso_prt_pops( long ipISO, long nelem, bool lgPrtDeparCoef )
00473 {
00474 long int in, il, is, i, ipLo, nResolved, ipFirstCollapsed=LONG_MIN;
00475 char chPrtType[2][12]={"populations","departure"};
00476
00477 char chSpin[3][9]= {"singlets", "doublets", "triplets"};
00478
00479 #define ITEM_TO_PRINT(A_) ( lgPrtDeparCoef ? iso.DepartCoef[ipISO][nelem][A_] : StatesElemNEW[nelem][nelem-ipISO][A_].Pop )
00480
00481 DEBUG_ENTRY( "iso_prt_pops()" );
00482
00483 ASSERT( ipISO < NISO );
00484
00485 for( is = 1; is<=3; ++is)
00486 {
00487 if( ipISO == ipH_LIKE && is != 2 )
00488 continue;
00489 else if( ipISO == ipHE_LIKE && is != 1 && is != 3 )
00490 continue;
00491
00492 ipFirstCollapsed= iso.numLevels_local[ipISO][nelem]-iso.nCollapsed_local[ipISO][nelem];
00493 nResolved = StatesElemNEW[nelem][nelem-ipISO][ipFirstCollapsed-1].n;
00494 ASSERT( nResolved == iso.n_HighestResolved_local[ipISO][nelem] );
00495 ASSERT(nResolved > 0 );
00496
00497
00498 fprintf(ioQQQ," %s %s %s %s\n",
00499 iso.chISO[ipISO],
00500 elementnames.chElementSym[nelem],
00501 chSpin[is-1],
00502 chPrtType[lgPrtDeparCoef]);
00503
00504
00505 fprintf(ioQQQ," n\\l=> ");
00506 for( i =0; i < nResolved; ++i)
00507 {
00508 fprintf(ioQQQ,"%2ld ",i);
00509 }
00510 fprintf(ioQQQ,"\n");
00511
00512
00513 for( in = 1; in <= nResolved; ++in)
00514 {
00515 if( is==3 && in==1 )
00516 continue;
00517
00518 fprintf(ioQQQ," %2ld ",in);
00519
00520 for( il = 0; il < in; ++il)
00521 {
00522 if( ipISO==ipHE_LIKE && (in==2) && (il==1) && (is==3) )
00523 {
00524 fprintf( ioQQQ, "%9.3e ", ITEM_TO_PRINT(ipHe2p3P0) );
00525 fprintf( ioQQQ, "%9.3e ", ITEM_TO_PRINT(ipHe2p3P1) );
00526 fprintf( ioQQQ, "%9.3e ", ITEM_TO_PRINT(ipHe2p3P2) );
00527 }
00528 else
00529 {
00530 ipLo = iso.QuantumNumbers2Index[ipISO][nelem][in][il][is];
00531 fprintf( ioQQQ, "%9.3e ", ITEM_TO_PRINT(ipLo) );
00532 }
00533 }
00534 fprintf(ioQQQ,"\n");
00535 }
00536 }
00537
00538 for( il = ipFirstCollapsed; il < iso.numLevels_local[ipISO][nelem]; ++il)
00539 {
00540 in = StatesElemNEW[nelem][nelem-ipISO][il].n;
00541
00542 fprintf(ioQQQ," %2ld ",in);
00543 fprintf( ioQQQ, "%9.3e ", ITEM_TO_PRINT(il) );
00544 fprintf(ioQQQ,"\n");
00545 }
00546 return;
00547 }
00548
00549
00550 void AGN_He1_CS( FILE *ioPun )
00551 {
00552
00553 long int i;
00554
00555
00556 const int NTE = 5;
00557 double TeList[NTE] = {6000.,10000.,15000.,20000.,25000.};
00558 double TempSave;
00559
00560 DEBUG_ENTRY( "AGN_He1_CS()" );
00561
00562
00563 fprintf(ioPun, "Te\t2 3s 33s\n");
00564
00565
00566 TempSave = phycon.te;
00567
00568 for( i=0; i<NTE; ++i )
00569 {
00570 TempChange(TeList[i] , false);
00571
00572 fprintf(ioPun , "%.0f\t",
00573 TeList[i] );
00574 fprintf(ioPun , "%.2f\t",
00575 HeCSInterp( 1 , ipHe3s3S , ipHe2s3S, ipELECTRON ) );
00576 fprintf(ioPun , "%.2f\t",
00577 HeCSInterp( 1 , ipHe3p3P , ipHe2s3S, ipELECTRON ) );
00578 fprintf(ioPun , "%.2f\t",
00579 HeCSInterp( 1 , ipHe3d3D , ipHe2s3S, ipELECTRON ) );
00580 fprintf(ioPun , "%.3f\t",
00581 HeCSInterp( 1 , ipHe3d1D , ipHe2s3S, ipELECTRON ) );
00582
00583 fprintf(ioPun , "%.1f\n",
00584 HeCSInterp( 1 , ipHe2p3P0 , ipHe2s3S, ipELECTRON ) +
00585 HeCSInterp( 1 , ipHe2p3P1 , ipHe2s3S, ipELECTRON ) +
00586 HeCSInterp( 1 , ipHe2p3P2 , ipHe2s3S, ipELECTRON ));
00587 }
00588
00589
00590 TempChange(TempSave , false);
00591 return;
00592 }