cloudy: /home66/gary/public_html/cloudy/c13_branch/source/hydrolevel.cpp Source File

00001 /* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
00002  * others.  For conditions of distribution and use see copyright notice in license.txt */
00003 /*HydroLevel calls iso_level to solve for ionization balance 
00004  * and level populations of model hydrogen atom */
00005 /*PrtHydroTrace1 print trace info for hydrogen-like species */
00006 #include "cddefines.h"
00007 #include "taulines.h"
00008 #include "iso.h"
00009 #include "dense.h"
00010 #include "secondaries.h"
00011 #include "trace.h"
00012 #include "phycon.h"
00013 #include "ionbal.h"
00014 #include "hydrogenic.h"
00015 
00016 /*PrtHydroTrace1a print trace info for hydrogen-like species */
00017 STATIC void PrtHydroTrace1a(long nelem )
00018 {
00019         double colfrc, 
00020           phtfrc, 
00021           secfrc;
00022 
00023         DEBUG_ENTRY( "PrtHydroTrace1a()" );
00024 
00025         /* identify how atom is ionized for full trace */
00026         if( iso_sp[ipH_LIKE][nelem].xIonSimple > 0. )
00027         {
00028                 /* fraction of ionization due to photoionization */
00029                 phtfrc = iso_sp[ipH_LIKE][nelem].fb[ipH1s].gamnc/((dense.eden*(iso_sp[ipH_LIKE][nelem].RadRec_effec + 
00030                         ionbal.CotaRate[nelem]) )*
00031                         iso_sp[ipH_LIKE][nelem].xIonSimple);
00032 
00033                 /* fraction of ionization due to collisional ionization */
00034                 colfrc = (iso_sp[ipH_LIKE][nelem].fb[ipH1s].ColIoniz*dense.EdenHCorr )/
00035                         ((dense.eden*(iso_sp[ipH_LIKE][nelem].RadRec_effec + 
00036                         ionbal.CotaRate[0]) )*
00037                         iso_sp[ipH_LIKE][nelem].xIonSimple);
00038 
00039                 /* fraction of ionization due to secondary ionization */
00040                 secfrc = secondaries.csupra[nelem][nelem]/((dense.eden*(iso_sp[ipH_LIKE][nelem].RadRec_effec + 
00041                         ionbal.CotaRate[0]) )*
00042                         iso_sp[ipH_LIKE][nelem].xIonSimple);
00043         }
00044         else
00045         {
00046                 phtfrc = 0.;
00047                 colfrc = 0.;
00048                 secfrc = 0.;
00049         }
00050 
00051         fprintf( ioQQQ, "     HydroLevel Z=%2ld called, simple II/I=",nelem);
00052         PrintE93( ioQQQ, iso_sp[ipH_LIKE][nelem].xIonSimple);
00053         fprintf( ioQQQ," PhotFrc:");
00054         PrintE93( ioQQQ,phtfrc);
00055         fprintf(ioQQQ," ColFrc:");
00056         PrintE93( ioQQQ,colfrc);
00057         fprintf( ioQQQ," SecFrc");
00058         PrintE93( ioQQQ, secfrc);
00059         fprintf( ioQQQ,"  Te:");
00060         PrintE93( ioQQQ,phycon.te);
00061         fprintf( ioQQQ," eden:");
00062         PrintE93( ioQQQ,dense.eden);
00063         fprintf( ioQQQ,"\n"); 
00064         return;
00065 }
00066 
00067 /*PrtHydroTrace1 print trace info for hydrogen-like species */
00068 STATIC void PrtHydroTrace1(long nelem )
00069 {
00070         long int ipHi , ipLo , i;
00071 
00072         DEBUG_ENTRY( "PrtHydroTrace1()" );
00073 
00074         fprintf( ioQQQ, 
00075                 "       HydroLevel%3ld finds arrays, with optical depths defined? %li induced 2ph=%12.3e\n", 
00076                 nelem, iteration, iso_sp[ipH_LIKE][nelem].TwoNu[0].induc_up );
00077         /* 06 aug 28, from numLevels_max to _local. */
00078         for( ipHi=ipH2s; ipHi < iso_sp[ipH_LIKE][nelem].numLevels_local; ipHi++ )
00079         {
00080                 fprintf( ioQQQ, "up:%2ld", ipHi );
00081                 fprintf( ioQQQ, "lo" );
00082                 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
00083                 {
00084                         fprintf( ioQQQ, "%9ld", ipLo );
00085                 }
00086                 fprintf( ioQQQ, "\n" );
00087 
00088                 fprintf( ioQQQ, "%3ld", ipHi );
00089                 fprintf( ioQQQ, " A*esc" );
00090                 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
00091                 {
00092                         fprintf( ioQQQ,PrintEfmt("%9.2e",  iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().Aul()*
00093                                 iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().Pesc() ));
00094                 }
00095                 fprintf( ioQQQ, "\n" );
00096 
00097                 fprintf( ioQQQ, "%3ld", ipHi );
00098                 fprintf( ioQQQ, " A*ees" );
00099                 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
00100                 {
00101                         fprintf( ioQQQ,PrintEfmt("%9.2e",  iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().Aul()*
00102                                 iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().Pelec_esc() ));
00103                 }
00104                 fprintf( ioQQQ, "\n" );
00105 
00106                 fprintf( ioQQQ, "%3ld", ipHi );
00107                 fprintf( ioQQQ, " tauin" );
00108                 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
00109                 {
00110                         fprintf( ioQQQ,PrintEfmt("%9.2e",  iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().TauIn() ));
00111                 }
00112                 fprintf( ioQQQ, "\n" );
00113 
00114                 fprintf( ioQQQ, "%3ld", ipHi );
00115                 fprintf( ioQQQ, " t tot" );
00116                 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
00117                 {
00118                         fprintf( ioQQQ,PrintEfmt("%9.2e", iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().TauTot() ));
00119                 }
00120                 fprintf( ioQQQ, "\n" );
00121 
00122                 fprintf( ioQQQ, "%3ld", ipHi );
00123                 fprintf( ioQQQ, " Esc  " );
00124                 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
00125                 {
00126                         fprintf( ioQQQ,PrintEfmt("%9.2e", iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().Pesc() ));
00127                 }
00128                 fprintf( ioQQQ, "\n" );
00129 
00130                 fprintf( ioQQQ, "%3ld", ipHi );
00131                 fprintf( ioQQQ, " Eesc " );
00132                 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
00133                 {
00134                         fprintf( ioQQQ,PrintEfmt("%9.2e", iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().Pelec_esc() ));
00135                 }
00136                 fprintf( ioQQQ, "\n" );
00137 
00138                 fprintf( ioQQQ, "%3ld", ipHi );
00139                 fprintf( ioQQQ, " Dest " );
00140                 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
00141                 {
00142                         fprintf( ioQQQ,PrintEfmt("%9.2e",  iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().Pdest()) );
00143                 }
00144                 fprintf( ioQQQ, "\n" );
00145 
00146                 fprintf( ioQQQ, "%3ld", ipHi );
00147                 fprintf( ioQQQ, " A*dst" );
00148                 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
00149                 {
00150                         fprintf( ioQQQ,PrintEfmt("%9.2e",  iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().Aul()*
00151                                 iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().Pdest() ));
00152                 }
00153                 fprintf( ioQQQ, "\n" );
00154 
00155                 fprintf( ioQQQ, "%3ld", ipHi );
00156                 fprintf( ioQQQ, " StrkE" );
00157                 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
00158                 {
00159                         fprintf( ioQQQ,PrintEfmt("%9.2e",  iso_sp[ipH_LIKE][nelem].ex[ipHi][ipLo].pestrk_up ));
00160                 }
00161                 fprintf( ioQQQ, "\n" );
00162 
00163                 fprintf( ioQQQ, "%3ld", ipHi );
00164                 fprintf( ioQQQ, " B(ul)" );
00165                 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
00166                 {
00167                         fprintf( ioQQQ,PrintEfmt("%9.2e", iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().pump()*
00168                                 iso_sp[ipH_LIKE][nelem].st[ipLo].g()/iso_sp[ipH_LIKE][nelem].st[ipHi].g() ));
00169                 }
00170                 fprintf( ioQQQ, "\n" );
00171 
00172                 fprintf( ioQQQ, "%3ld", ipHi );
00173                 fprintf( ioQQQ, " tcont" );
00174                 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
00175                 {
00176                         fprintf( ioQQQ,PrintEfmt("%9.2e",  iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().TauCon() ));
00177                 }
00178                 fprintf( ioQQQ, "\n" );
00179 
00180                 fprintf( ioQQQ, "%3ld", ipHi );
00181                 fprintf( ioQQQ, " C(ul)" );
00182                 for( ipLo=ipH1s; ipLo < ipHi; ipLo++ )
00183                 {
00184                         fprintf( ioQQQ,PrintEfmt("%9.2e", iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Coll().ColUL( colliders ) ));
00185                 }
00186                 fprintf( ioQQQ, "\n" );
00187 
00188                 if( ipHi == 2 )
00189                 {
00190                         fprintf( ioQQQ, "    FeIIo");
00191                         fprintf( ioQQQ,PrintEfmt("%9.2e", 
00192                                 hydro.dstfe2lya* iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().Aul() ));
00193                         fprintf( ioQQQ, "\n");
00194                 }
00195         }
00196 
00197         fprintf( ioQQQ, "         " );
00198         /* 06 aug 28, from numLevels_max to _local. */
00199         for( i=1; i < iso_sp[ipH_LIKE][nelem].numLevels_local; i++ )
00200         {
00201                 fprintf( ioQQQ, "%9ld", i );
00202         }
00203         fprintf( ioQQQ, "\n" );
00204         return;
00205 }
00206 
00207 /*HydroLevel calls iso_level to solve for ionization balance 
00208  * and level populations of model hydrogen atom */
00209 void HydroLevel(long int nelem)
00210 {
00211         long int i; 
00212         int ipISO = ipH_LIKE;
00213 
00214         DEBUG_ENTRY( "HydroLevel()" );
00215 
00216         /* check that we were called with valid charge */
00217         ASSERT( nelem >= 0);
00218         ASSERT( nelem < LIMELM );
00219 
00220         /* option to print some rates */
00221         if( (trace.lgTrace && trace.lgIsoTraceFull[ipISO]) && (nelem == trace.ipIsoTrace[ipISO]) )
00222                 PrtHydroTrace1(nelem);
00223 
00224         if( trace.lgHBug && trace.lgTrace )
00225                 PrtHydroTrace1a(nelem);
00226 
00227         /* this is main trace h-like printout */
00228         if( (trace.lgIsoTraceFull[ipISO] && trace.lgTrace) && (nelem == trace.ipIsoTrace[ipISO]) )
00229         {
00230                 fprintf( ioQQQ, "       HLEV HGAMNC" );
00231                 PrintE93( ioQQQ, iso_sp[ipISO][nelem].fb[ipH1s].gamnc );
00232                 /* 06 aug 28, from numLevels_max to _local. */
00233                 for( i=ipH2s; i < iso_sp[ipISO][nelem].numLevels_local; i++ )
00234                 {
00235                         fprintf(ioQQQ,PrintEfmt("%9.2e", iso_sp[ipISO][nelem].fb[i].gamnc ));
00236                 }
00237                 fprintf( ioQQQ, "\n" );
00238 
00239                 fprintf( ioQQQ, "       HLEV TOTCAP" );
00240                 /* 06 aug 28, from numLevels_max to _local. */
00241                 for( i=1; i < iso_sp[ipISO][nelem].numLevels_local; i++ )
00242                 {
00243                         fprintf(ioQQQ,PrintEfmt("%9.2e", iso_sp[ipISO][nelem].fb[i].RateCont2Level/dense.eden ));
00244                 }
00245                 fprintf( ioQQQ," tot");
00246                 fprintf( ioQQQ,PrintEfmt("%10.2e", ionbal.RateRecomTot[nelem][nelem-ipISO]/dense.eden ) );
00247                 fprintf( ioQQQ, "\n" );
00248 
00249                 fprintf( ioQQQ, "       HLEV IND Rc" );
00250                 /* 06 aug 28, from numLevels_max to _local. */
00251                 for( i=ipH1s; i < iso_sp[ipISO][nelem].numLevels_local; i++ )
00252                 {
00253                         fprintf(ioQQQ,PrintEfmt("%9.2e", iso_sp[ipISO][nelem].fb[i].RecomInducRate*iso_sp[ipISO][nelem].fb[i].PopLTE ));
00254                 }
00255                 fprintf( ioQQQ, "\n" );
00256 
00257                 /* incuded recombination rate coefficients */
00258                 fprintf( ioQQQ, "       IND Rc LTE " );
00259                 /* 06 aug 28, from numLevels_max to _local. */
00260                 for( i=ipH1s; i < iso_sp[ipISO][nelem].numLevels_local; i++ )
00261                 {
00262                         fprintf(ioQQQ,PrintEfmt("%9.2e",
00263                                 iso_sp[ipISO][nelem].fb[i].gamnc*iso_sp[ipISO][nelem].fb[i].PopLTE ));
00264                 }
00265                 fprintf( ioQQQ, "\n" );
00266 
00267                 /* LTE level populations */
00268                 fprintf( ioQQQ, "       HLEV   HLTE" );
00269                 /* 06 aug 28, from numLevels_max to _local. */
00270                 for( i=ipH1s; i < iso_sp[ipISO][nelem].numLevels_local; i++ )
00271                 {
00272                         fprintf(ioQQQ,PrintEfmt("%9.2e", iso_sp[ipISO][nelem].fb[i].PopLTE ));
00273                 }
00274                 fprintf( ioQQQ, "\n" );
00275 
00276                 /* fraction of total ionization due to collisions from given level */
00277                 fprintf( ioQQQ, "       HLEVfr cion" );
00278                 /* 06 aug 28, from numLevels_max to _local. */
00279                 for( i=ipH1s; i < iso_sp[ipISO][nelem].numLevels_local; i++ )
00280                 {
00281                         fprintf(ioQQQ,PrintEfmt("%9.2e", 
00282                                 iso_sp[ipISO][nelem].fb[i].ColIoniz*dense.EdenHCorr/MAX2(1e-30,iso_sp[ipISO][nelem].fb[i].RateLevel2Cont) ) );
00283                 }
00284                 fprintf( ioQQQ, "\n" );
00285 
00286                 /* fraction of total ionization due to photoionization from given level */
00287                 if( ionbal.RateRecomTot[nelem][nelem]> 0. )
00288                 {
00289                         fprintf( ioQQQ, "       HLEVfrPhIon" );
00290                         /* 06 aug 28, from numLevels_max to _local. */
00291                         for( i=ipH1s; i < iso_sp[ipISO][nelem].numLevels_local; i++ )
00292                         {
00293                                 fprintf(ioQQQ,PrintEfmt("%9.2e", 
00294                                         iso_sp[ipISO][nelem].fb[i].gamnc/MAX2(1e-30,iso_sp[ipISO][nelem].fb[i].RateLevel2Cont) ) );
00295                         }
00296                         fprintf( ioQQQ, "\n" );
00297                 }
00298 
00299                 fprintf( ioQQQ, "       HLEV     HN" );
00300                 /* 06 aug 28, from numLevels_max to _local. */
00301                 for( i=ipH1s; i < iso_sp[ipISO][nelem].numLevels_local; i++ )
00302                 {
00303                         fprintf(ioQQQ,PrintEfmt("%9.2e", iso_sp[ipISO][nelem].st[i].Pop() ));
00304                 }
00305                 fprintf( ioQQQ, "\n" );
00306 
00307                 fprintf( ioQQQ, "       HLEV   b(n)" );
00308                 /* 06 aug 28, from numLevels_max to _local. */
00309                 for( i=ipH1s; i < iso_sp[ipISO][nelem].numLevels_local; i++ )
00310                 {
00311                         fprintf(ioQQQ,PrintEfmt("%9.2e", iso_sp[ipISO][nelem].fb[i].DepartCoef ));
00312                 }
00313                 fprintf( ioQQQ, "\n" );
00314 
00315                 fprintf( ioQQQ, "       HLEV   X12tot");
00316                 fprintf(ioQQQ,PrintEfmt("%9.2e", secondaries.x12tot ));
00317                 fprintf( ioQQQ," Grn dest:");
00318                 fprintf(ioQQQ,PrintEfmt("%9.2e",
00319                   ionbal.RateIoniz[nelem][nelem][nelem+1] ));
00320                 fprintf(ioQQQ, "\n"); 
00321         }
00322 
00323         if( trace.lgTrace )
00324         {
00325                 /* iso.RecomTotal[nelem] is gross rec coef, computed here while filling in matrix
00326                  * elements, all physical processes included. 
00327                  * RadRec_effec is total effective radiative only */
00328                 fprintf( ioQQQ, "       HydroLevel Z:%2ld return %s te=",
00329                         nelem,
00330                         iso_sp[ipISO][nelem].chTypeAtomUsed );
00331                 PrintE93( ioQQQ,phycon.te);
00332                 fprintf( ioQQQ," density=%.4e", dense.xIonDense[nelem][nelem-ipISO] );
00333 
00334                 fprintf( ioQQQ," simple=%.4e",iso_sp[ipISO][nelem].xIonSimple);
00335 
00336                 fprintf( ioQQQ," b1=%.2e",iso_sp[ipISO][nelem].fb[ipH1s].DepartCoef);
00337 
00338                 fprintf( ioQQQ," ion rate=%.4e",ionbal.RateIonizTot(nelem,nelem-ipISO) );
00339 
00340                 fprintf( ioQQQ," TotRec=%.4e",ionbal.RateRecomTot[nelem][nelem-ipISO]);
00341 
00342                 fprintf( ioQQQ," RadRec=%.4e",iso_sp[ipISO][nelem].RadRec_effec);
00343                 fprintf( ioQQQ, "\n");
00344         }
00345         return;
00346 }