cloudy: /home66/gary/public_html/cloudy/c08_branch/source/prt_lines

00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
00002  * others.  For conditions of distribution and use see copyright notice in license.txt */
00003 /*lines_helium put He-like iso sequence into line intensity stack */
00004 /*TempInterp interpolates on a grid of values to produce predicted value at current Te.*/
00005 #include "cddefines.h"
00006 #include "dense.h"
00007 #include "helike.h"
00008 #include "iso.h"
00009 #include "atmdat.h"
00010 #include "lines.h"
00011 #include "lines_service.h"
00012 #include "phycon.h"
00013 #include "physconst.h"
00014 #include "taulines.h"
00015 #include "thirdparty.h"
00016 #include "trace.h"
00017 
00018 #define NUMTEMPS        22
00019 
00020 typedef struct 
00021 {
00022         /* index for upper and lower levels of line */
00023         long int ipHi;
00024         long int ipLo;
00025 
00026         char label[5];
00027 
00028 } stdLines;
00029 
00030 STATIC void GetStandardHeLines(void);
00031 STATIC double TempInterp2( double* TempArray , double* ValueArray, long NumElements );
00032 STATIC void DoSatelliteLines( long nelem );
00033 
00034 static bool lgFirstRun = true;
00035 static double CaABTemps[NUMTEMPS];
00036 static long NumLines;
00037 static double ***CaABIntensity;
00038 static stdLines **CaABLines;
00039 
00040 void lines_helium(void)
00041 {
00042         long ipISO = ipHE_LIKE;
00043         long int i, nelem, ipHi, ipLo;
00044         char chLabel[5]="    ";
00045 
00046         long int j;
00047 
00048         double 
00049           sum,
00050           Pop2_3S,
00051           photons_3889_plus_7065 = 0.;
00052 
00053         DEBUG_ENTRY( "lines_helium()" );
00054 
00055         if( trace.lgTrace )
00056                 fprintf( ioQQQ, "   prt_lines_helium called\n" );
00057 
00058         i = StuffComment( "He-like iso-sequence" );
00059         linadd( 0., (realnum)i , "####", 'i',
00060                 " start He-like iso sequence");
00061 
00062         /* read in Case A and B lines from data file    */
00063         if( lgFirstRun )
00064         {
00065                 GetStandardHeLines();
00066                 lgFirstRun = false;
00067         }
00068 
00069         /* store labels for all case b HeI lines in case we assert case b 
00070          * ipass == -1 only counting number of lines, =0, malloc then set wl */
00071         static bool lgMustMalloc=true;
00072         if( LineSave.ipass == 0 && atmdat.nCaseBHeI>0 && lgMustMalloc )
00073         {
00074                 /* second time through - on ipass=-1 we counted number of lines
00075                  * atmdat.nCaseBHeI, now create space but only if there are He I lines 
00076                  * this is not done if He is turned off */
00077                 atmdat.CaseBWlHeI = (realnum*)MALLOC( sizeof(realnum)*atmdat.nCaseBHeI);
00078                 lgMustMalloc=false;
00079         }
00080         atmdat.nCaseBHeI = 0;
00081 
00082         /* this is the main printout, where line intensities are entered into the stack */
00083         for( nelem=ipISO; nelem < LIMELM; nelem++ )
00084         {
00085                 if( dense.lgElmtOn[nelem] )
00086                 {
00087                         if( nelem == ipHELIUM )
00088                         {
00089                                 double *qTotEff;
00090 
00091                                 /* >>chng 06 aug 17, all of these from _max to _local */
00092                                 /* >>chng 06 dec 21, mistake - change back to _max */
00093                                 qTotEff = (double*)MALLOC(sizeof(double)*(unsigned)(iso.numLevels_max[ipHE_LIKE][nelem]) );
00094 
00095                                 qTotEff[0] = 0.;
00096                                 qTotEff[1] = 0.;
00097 
00098                                 for( i=ipHe2s3S+1; i<iso.numLevels_max[ipHE_LIKE][nelem]-iso.nCollapsed_max[ipHE_LIKE][nelem]; i++ )
00099                                 {
00100                                         qTotEff[i] = 0.;
00101                                         for( j = i; j<iso.numLevels_max[ipHE_LIKE][nelem]-iso.nCollapsed_max[ipHE_LIKE][nelem]; j++ )
00102                                         {
00103                                                 /*if( StatesElem[ipHE_LIKE][nelem][i].S == 3 )
00104                                                 {*/
00105                                                         qTotEff[i] += 
00106                                                                 Transitions[ipHE_LIKE][nelem][j][ipHe2s3S].Coll.ColUL*dense.EdenHCorr*
00107                                                                 iso.Boltzmann[ipHE_LIKE][nelem][j][ipHe2s3S] *
00108                                                                 (double)Transitions[ipHE_LIKE][nelem][j][ipHe2s3S].Hi->g / 
00109                                                                 (double)Transitions[ipHE_LIKE][nelem][j][ipHe2s3S].Lo->g*
00110                                                                 iso.CascadeProb[ipISO][nelem][j][i];
00111                                                 /*}*/
00112                                         }
00113                                 }
00114                                 
00115                                 /* get simple 2^3S pop, assume recombinations in are just 0.75 * case B */
00116                                 Pop2_3S = dense.eden*(0.75*iso.RadRec_caseB[ipHE_LIKE][nelem])/
00117                                         ( Transitions[ipHE_LIKE][nelem][ipHe2s3S][ipHe1s1S].Emis->Aul+ dense.eden*iso.qTot2S[ipISO][nelem]);
00118 
00119                                 for( i=0; i< NumLines; i++ )
00120                                 {
00121                                         ipHi = CaABLines[nelem][i].ipHi;
00122                                         ipLo = CaABLines[nelem][i].ipLo;
00123 
00124                                         /* >>chng 06 aug 17, from _max to _local */
00125                                         /* >>chng 06 dec 21, mistake - change back to _max */
00126                                         if( ipHi <= iso.n_HighestResolved_max[ipHE_LIKE][nelem]*(iso.n_HighestResolved_max[ipHE_LIKE][nelem]+1))
00127                                         {
00128                                                 double intens = TempInterp2( CaABTemps , CaABIntensity[nelem][i], NUMTEMPS )*
00129                                                         dense.xIonDense[nelem][nelem+1-ipISO]*dense.eden;
00130 
00131                                                 linadd( intens,
00132                                                         Transitions[ipHE_LIKE][nelem][ipHi][ipLo].WLAng,
00133                                                         CaABLines[nelem][i].label,'i',
00134                                                         "Case B intensity ");
00135 
00136                                                 if( nMatch("Ca B",CaABLines[nelem][i].label) )
00137                                                 {
00138                                                         /* all lines to/from 2^3Pj are stored as lines to/from 2^3P1, so make sure this loop never tries to 
00139                                                          * explicitly consider 2^3P0 or 2^3P2 */
00140                                                         ASSERT( ipLo!=ipHe2p3P0 && ipLo!=ipHe2p3P2 );
00141                                                         ASSERT( ipHi!=ipHe2p3P0 && ipHi!=ipHe2p3P2 );
00142 
00143                                                         double totBranch = iso.BranchRatio[ipISO][nelem][ipHi][ipLo];
00144                                                         if( ipLo==4 )
00145                                                                 totBranch += iso.BranchRatio[ipISO][nelem][ipHi][3] + iso.BranchRatio[ipISO][nelem][ipHi][5];
00146 
00147                                                         if( LineSave.ipass < 0 )
00148                                                                 ++atmdat.nCaseBHeI;
00149                                                         else if( LineSave.ipass == 0 )
00150                                                         {
00151                                                                 /* save wavelengths */
00152                                                                 atmdat.CaseBWlHeI[atmdat.nCaseBHeI] = 
00153                                                                         Transitions[ipHE_LIKE][nelem][ipHi][ipLo].WLAng;
00154                                                                 ++atmdat.nCaseBHeI;
00155                                                         }
00156 
00157                                                         if( ipHi==4 )
00158                                                         {
00159                                                                 linadd( intens + 
00160                                                                         Pop2_3S*dense.xIonDense[nelem][nelem+1-ipISO]*
00161                                                                         (
00162                                                                         qTotEff[ipHe2p3P0]*iso.BranchRatio[ipISO][nelem][ipHe2p3P0][ipLo]+
00163                                                                         qTotEff[ipHe2p3P1]*iso.BranchRatio[ipISO][nelem][ipHe2p3P1][ipLo]+
00164                                                                         qTotEff[ipHe2p3P2]*iso.BranchRatio[ipISO][nelem][ipHe2p3P2][ipLo]
00165                                                                         )*
00166                                                                         Transitions[ipHE_LIKE][nelem][ipHi][ipLo].EnergyErg,
00167                                                                         Transitions[ipHE_LIKE][nelem][ipHi][ipLo].WLAng,
00168                                                                         "+Col",'i',
00169                                                                         "Case B intensity with collisions included");
00170 
00171                                                         }
00172                                                         else
00173                                                         {
00174                                                                 /* chng 05 dec 14, branching ratio was missing here! 
00175                                                                  * not a big effect because lines with biggest collision
00176                                                                  * enhancements tend to be dominant decay route from upper level. */
00177                                                                 linadd( intens + 
00178                                                                         Pop2_3S*qTotEff[ipHi]*dense.xIonDense[nelem][nelem+1-ipISO]*totBranch*
00179                                                                         Transitions[ipHE_LIKE][nelem][ipHi][ipLo].EnergyErg,
00180                                                                         Transitions[ipHE_LIKE][nelem][ipHi][ipLo].WLAng,
00181                                                                         "+Col",'i',
00182                                                                         "Case B intensity with collisions included");
00183                                                         }
00184                                                 }
00185                                         }
00186                                         /* Make sure to at least do 4471 */
00187                                         else if( ipLo==ipHe2p3P1 && ipHi==ipHe4d3D && nMatch("Ca B",CaABLines[nelem][i].label) )
00188                                         {
00189                                                 double intens = TempInterp2( CaABTemps , CaABIntensity[nelem][i], NUMTEMPS )*
00190                                                         dense.xIonDense[nelem][nelem+1-ipISO]*dense.eden;
00191 
00192                                                 linadd( intens, 4471, CaABLines[nelem][i].label, 'i',
00193                                                         "Case B intensity ");
00194                                         }
00195 
00196                                 }
00197                                 free( qTotEff );
00198                         }
00199 
00200                         /* NB NB - low and high must be in this order so that all balmer, paschen,
00201                          * etc series line up correctly in final printout */
00202                         /* >>chng 01 jun 13, bring 23P lines back together */
00203                         /* two photon is special, not a line and no ipCont array index, add here */
00204                         Transitions[ipHE_LIKE][nelem][ipHe2s1S][ipHe1s1S].Emis->phots = 
00205                                 Transitions[ipHE_LIKE][nelem][ipHe2s1S][ipHe1s1S].Emis->Aul*
00206                                 StatesElem[ipHE_LIKE][nelem][ipHe2s1S].Pop*
00207                                 Transitions[ipHE_LIKE][nelem][ipHe2s1S][ipHe1s1S].Emis->Pesc*
00208                                 dense.xIonDense[nelem][nelem+1-ipISO];
00209 
00210                         Transitions[ipHE_LIKE][nelem][ipHe2s1S][ipHe1s1S].Emis->xIntensity = 
00211                                 Transitions[ipHE_LIKE][nelem][ipHe2s1S][ipHe1s1S].Emis->phots*
00212                                 Transitions[ipHE_LIKE][nelem][ipHe2s1S][ipHe1s1S].EnergyErg;
00213 
00214                         if( LineSave.ipass == 0 )
00215                         {
00216                                 /* chIonLbl is function that generates a null terminated 4 char string, of form "C  2" 
00217                                  * the result, chLable, is only used when ipass == 0, can be undefined otherwise */
00218                                 /* total two photon emission */
00219                                 chIonLbl(chLabel, &Transitions[ipHE_LIKE][nelem][ipHe2s1S][ipHe1s1S]);
00220                         }
00221 
00222                         linadd( Transitions[ipHE_LIKE][nelem][ipHe2s1S][ipHe1s1S].Emis->xIntensity , 0,chLabel,'r',
00223                                 " two photon continuum ");
00224 
00225                         linadd(
00226                                 StatesElem[ipHE_LIKE][nelem][ipHe2s1S].Pop*
00227                                 dense.xIonDense[nelem][nelem+1-ipISO]*
00228                                 iso.TwoNu_induc_dn[ipHE_LIKE][nelem]*
00229                                 Transitions[ipHE_LIKE][nelem][ipHe2s1S][ipHe1s1S].EnergyErg,
00230                                 22, chLabel ,'i',
00231                                 " induced two photon emission ");
00232 
00233                         /* here we will create an entry for the three lines 
00234                          * coming from 2 3P to 1 1S - the entry called TOTL will
00235                          * appear before the lines of the multiplet */
00236                         sum = 0.;
00237                         for( i=ipHe2p3P0; i <= ipHe2p3P2; i++ )
00238                         {
00239                                 if( Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].ipCont <= 0 )
00240                                         continue;
00241 
00242                                 sum += 
00243                                         Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Emis->Aul*
00244                                         StatesElem[ipHE_LIKE][nelem][i].Pop*
00245                                         Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].Emis->Pesc*
00246                                         dense.xIonDense[nelem][nelem+1-ipISO]*
00247                                         Transitions[ipHE_LIKE][nelem][i][ipHe1s1S].EnergyErg;
00248                         }
00249 
00250                         linadd(sum,Transitions[ipHE_LIKE][nelem][ipHe2p3P1][ipHe1s1S].WLAng,"TOTL",'i' ,
00251                                 " total emission in He-like intercombination lines from 2p3P to ground ");
00252 
00253                         /* now do real permitted lines */
00254                         for( ipLo=0; ipLo < ipHe2p3P0; ipLo++ )
00255                         {
00256                                 for( ipHi=ipLo+1; ipHi < iso.numPrintLevels[ipHE_LIKE][nelem]; ipHi++ )
00257                                 {
00258                                         /* >>chng 01 may 30, do not add fake he-like lines (majority) to line stack */
00259                                         /* >>chng 01 dec 11, use variable for smallest A */
00260                                         if( Transitions[ipHE_LIKE][nelem][ipHi][ipLo].ipCont < 1 ) 
00261                                                 continue;
00262 
00263                                         /* recombine fine-structure lines since the energies are
00264                                          * not resolved anyway. */
00265                                         if( iso.lgFSM[ipISO] && ( abs(StatesElem[ipHE_LIKE][nelem][ipHi].l -
00266                                                 StatesElem[ipHE_LIKE][nelem][ipLo].l)==1 )
00267                                                 && (StatesElem[ipHE_LIKE][nelem][ipLo].l>1) 
00268                                                 && (StatesElem[ipHE_LIKE][nelem][ipHi].l>1) 
00269                                                 && ( StatesElem[ipHE_LIKE][nelem][ipHi].n ==
00270                                                 StatesElem[ipHE_LIKE][nelem][ipLo].n ) )
00271                                         {
00272                                                 /* skip if both singlets. */
00273                                                 if( (StatesElem[ipHE_LIKE][nelem][ipHi].S==1) 
00274                                                         && (StatesElem[ipHE_LIKE][nelem][ipLo].S==1) )
00275                                                 {
00276                                                         continue;
00277                                                 }
00278                                                 else if( (StatesElem[ipHE_LIKE][nelem][ipHi].S==3) 
00279                                                         && (StatesElem[ipHE_LIKE][nelem][ipLo].S==3) )
00280                                                 {
00281 
00282                                                         /* singlet to singlet*/
00283                                                         Transitions[ipHE_LIKE][nelem][ipHi+1][ipLo+1].Emis->phots = 
00284                                                                 ( Transitions[ipHE_LIKE][nelem][ipHi][ipLo+1].Emis->Aul*
00285                                                                 StatesElem[ipHE_LIKE][nelem][ipHi].Pop*
00286                                                                 Transitions[ipHE_LIKE][nelem][ipHi][ipLo+1].Emis->Pesc +
00287                                                                 Transitions[ipHE_LIKE][nelem][ipHi+1][ipLo+1].Emis->Aul*
00288                                                                 StatesElem[ipHE_LIKE][nelem][ipHi+1].Pop*
00289                                                                 Transitions[ipHE_LIKE][nelem][ipHi+1][ipLo+1].Emis->Pesc )*
00290                                                                 dense.xIonDense[nelem][nelem+1-ipISO];
00291 
00292                                                         Transitions[ipHE_LIKE][nelem][ipHi+1][ipLo+1].Emis->xIntensity = 
00293                                                                 Transitions[ipHE_LIKE][nelem][ipHi+1][ipLo+1].Emis->phots *
00294                                                                 Transitions[ipHE_LIKE][nelem][ipHi+1][ipLo+1].EnergyErg;
00295 
00296                                                         PutLine(&Transitions[ipHE_LIKE][nelem][ipHi+1][ipLo+1],
00297                                                                 " ");
00298 
00299                                                         /* triplet to triplet */
00300                                                         Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->phots = 
00301                                                                 ( Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->Aul*
00302                                                                 StatesElem[ipHE_LIKE][nelem][ipHi].Pop*
00303                                                                 Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->Pesc +
00304                                                                 Transitions[ipHE_LIKE][nelem][ipHi+1][ipLo].Emis->Aul*
00305                                                                 StatesElem[ipHE_LIKE][nelem][ipHi+1].Pop*
00306                                                                 Transitions[ipHE_LIKE][nelem][ipHi+1][ipLo].Emis->Pesc )*
00307                                                                 dense.xIonDense[nelem][nelem+1-ipISO];
00308 
00309                                                         Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->xIntensity = 
00310                                                                 Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->phots *
00311                                                                 Transitions[ipHE_LIKE][nelem][ipHi][ipLo].EnergyErg;
00312 
00313                                                         PutLine(&Transitions[ipHE_LIKE][nelem][ipHi][ipLo],
00314                                                                 " ");
00315                                                 }
00316                                         }
00317 
00318                                         else if( ipLo==ipHe2s3S && ipHi == ipHe2p3P0 )
00319                                         {
00320                                                 /* here we will create an entry for the three lines 
00321                                                  * coming from 2 3P to 2 3S - the entry called TOTL will
00322                                                  * appear before the lines of the multiplet 
00323                                                  * for He I this is 10830 */
00324 
00325                                                 realnum av_wl = 0.;
00326                                                 sum = 0.;
00327                                                 for( i=ipHe2p3P0; i <= ipHe2p3P2; i++ )
00328                                                 {
00329                                                         sum += 
00330                                                                 Transitions[ipHE_LIKE][nelem][i][ipLo].Emis->Aul*
00331                                                                 StatesElem[ipHE_LIKE][nelem][i].Pop*
00332                                                                 Transitions[ipHE_LIKE][nelem][i][ipLo].Emis->Pesc*
00333                                                                 dense.xIonDense[nelem][nelem+1-ipISO]*
00334                                                                 Transitions[ipHE_LIKE][nelem][i][ipLo].EnergyErg;
00335                                                         av_wl += Transitions[ipHE_LIKE][nelem][i][ipLo].WLAng;
00336                                                 }
00337                                                 av_wl /= 3.;
00338 #                                               if 0
00339                                                 {
00340 #                                               include "elementnames.h"
00341 #                                               include "prt.h"
00342                                                 fprintf(ioQQQ,"DEBUG 2P - 2S multiplet wl %s ",
00343                                                         elementnames.chElementSym[nelem] );
00344                                                 prt_wl( ioQQQ , av_wl );
00345                                                 fprintf(ioQQQ,"\n" );
00346                                                 }
00347 #                                               endif
00348 
00349                                                 linadd(sum,av_wl,"TOTL",'i',
00350                                                         "total emission in He-like lines, use average of three line wavelengths " );
00351 
00352                                                 /* also add this with the regular label, so it is correctly picked up by assert case-b command */
00353                                                 linadd(sum,av_wl,chLabel,'i',
00354                                                         "total emission in He-like lines, use average of three line wavelengths " );
00355 
00356                                                 /*>>chng 05 sep 8, added the following so that the component
00357                                                  * from ipHe2p3P0 is printed, in addition to the total. */
00358 
00359                                                 /* find number of photons escaping cloud */
00360                                                 Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->phots = 
00361                                                         Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->Aul*
00362                                                         StatesElem[ipHE_LIKE][nelem][ipHi].Pop*
00363                                                         Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->Pesc*
00364                                                         dense.xIonDense[nelem][nelem+1-ipISO];
00365 
00366                                                 /* now find line intensity */
00367                                                 Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->xIntensity = 
00368                                                         Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->phots*
00369                                                         Transitions[ipHE_LIKE][nelem][ipHi][ipLo].EnergyErg;
00370 
00371                                                 if( iso.lgRandErrGen[ipISO] )
00372                                                 {
00373                                                         Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->phots *=
00374                                                                 iso.ErrorFactor[ipISO][nelem][ipHi][ipLo][IPRAD];
00375                                                         Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->xIntensity *= 
00376                                                                 iso.ErrorFactor[ipISO][nelem][ipHi][ipLo][IPRAD];
00377                                                 }
00378                                                 PutLine(&Transitions[ipHE_LIKE][nelem][ipHi][ipLo],
00379                                                         " ");
00380                                         }
00381 
00382                                         else
00383                                         {
00384 
00385                                                 /* find number of photons escaping cloud */
00386                                                 Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->phots = 
00387                                                         Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->Aul*
00388                                                         StatesElem[ipHE_LIKE][nelem][ipHi].Pop*
00389                                                         Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->Pesc*
00390                                                         dense.xIonDense[nelem][nelem+1-ipISO];
00391 
00392                                                 /* now find line intensity */
00393                                                 /* >>chng 01 jan 15, put cast double to force double evaluation */
00394                                                 Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->xIntensity = 
00395                                                         Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->phots*
00396                                                         Transitions[ipHE_LIKE][nelem][ipHi][ipLo].EnergyErg;
00397 
00398                                                 if( iso.lgRandErrGen[ipISO] )
00399                                                 {
00400                                                         Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->phots *=
00401                                                                 iso.ErrorFactor[ipISO][nelem][ipHi][ipLo][IPRAD];
00402                                                         Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->xIntensity *= 
00403                                                                 iso.ErrorFactor[ipISO][nelem][ipHi][ipLo][IPRAD];
00404                                                 }
00405 
00406                                                 /* 
00407                                                 fprintf(ioQQQ,"1 loop %li %li %.1f\n", ipLo,ipHi, 
00408                                                         Transitions[ipHE_LIKE][nelem][ipHi][ipLo].WLAng ); */
00409                                                 PutLine(&Transitions[ipHE_LIKE][nelem][ipHi][ipLo],
00410                                                         "total intensity of He-like line");
00411                                                 {
00412                                                         /* option to print particulars of some line when called
00413                                                          * a prettier print statement is near where chSpin is defined below*/
00414                                                         enum {DEBUG_LOC=false};
00415                                                         if( DEBUG_LOC )
00416                                                         {
00417                                                                 if( nelem==1 && ipLo==0 && ipHi==1 )
00418                                                                 {
00419                                                                         fprintf(ioQQQ,"he1 626 %.2e %.2e \n", 
00420                                                                                 Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->TauIn,
00421                                                                                 Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->TauTot
00422                                                                                 );
00423                                                                 }
00424                                                         }
00425                                                 }
00426                                         }
00427                                 }
00428                         }
00429 
00430                         /* this sum will bring together the three lines going to J levels within 23P */
00431                         for( ipHi=ipHe2p3P2+1; ipHi < iso.numPrintLevels[ipHE_LIKE][nelem]; ipHi++ )
00432                         {
00433                                 double sumcool , sumheat ,
00434                                         save , savecool , saveheat;
00435 
00436                                 sum = 0;
00437                                 sumcool = 0.;
00438                                 sumheat = 0.;
00439                                 for( ipLo=ipHe2p3P0; ipLo <= ipHe2p3P2; ++ipLo )
00440                                 {
00441                                         if( Transitions[ipHE_LIKE][nelem][ipHi][ipLo].ipCont <= 0 )
00442                                                 continue;
00443 
00444                                         /* find number of photons escaping cloud */
00445                                         Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->phots = 
00446                                                 Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->Aul*
00447                                                 StatesElem[ipHE_LIKE][nelem][ipHi].Pop*
00448                                                 Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->Pesc*
00449                                                 dense.xIonDense[nelem][nelem+1-ipISO];
00450 
00451                                         /* now find line intensity */
00452                                         /* >>chng 01 jan 15, put cast double to force double evaluation */
00453                                         Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->xIntensity = 
00454                                                 Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->phots*
00455                                                 Transitions[ipHE_LIKE][nelem][ipHi][ipLo].EnergyErg;
00456 
00457                                         if( iso.lgRandErrGen[ipISO] )
00458                                         {
00459                                                 Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->phots *=
00460                                                         iso.ErrorFactor[ipISO][nelem][ipHi][ipLo][IPRAD];
00461                                                 Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->xIntensity *= 
00462                                                         iso.ErrorFactor[ipISO][nelem][ipHi][ipLo][IPRAD];
00463                                         }
00464 
00465                                         sumcool += Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Coll.cool;
00466                                         sumheat += Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Coll.heat;
00467                                         sum += Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->xIntensity;
00468                                 }
00469 
00470                                 /* skip non-radiative lines */
00471                                 if( Transitions[ipHE_LIKE][nelem][ipHi][ipHe2p3P2].ipCont < 1 ) 
00472                                         continue;
00473 
00474                                 /* this will enter .xIntensity into the line stack */
00475                                 save = Transitions[ipHE_LIKE][nelem][ipHi][ipHe2p3P2].Emis->xIntensity;
00476                                 savecool = Transitions[ipHE_LIKE][nelem][ipHi][ipHe2p3P2].Coll.cool;
00477                                 saveheat = Transitions[ipHE_LIKE][nelem][ipHi][ipHe2p3P2].Coll.heat;
00478 
00479                                 Transitions[ipHE_LIKE][nelem][ipHi][ipHe2p3P2].Emis->xIntensity = sum;
00480                                 Transitions[ipHE_LIKE][nelem][ipHi][ipHe2p3P2].Coll.cool = sumcool;
00481                                 Transitions[ipHE_LIKE][nelem][ipHi][ipHe2p3P2].Coll.heat = sumheat;
00482 
00483                                 /*fprintf(ioQQQ,"2 loop %li %li %.1f\n", ipHe2p3P2,ipHi, 
00484                                         Transitions[ipHE_LIKE][nelem][ipHi][ipHe2p3P2].WLAng );*/
00485                                 PutLine(&Transitions[ipHE_LIKE][nelem][ipHi][ipHe2p3P2],
00486                                         "predicted line, all processes included");
00487 
00488                                 Transitions[ipHE_LIKE][nelem][ipHi][ipHe2p3P2].Emis->xIntensity = save;
00489                                 Transitions[ipHE_LIKE][nelem][ipHi][ipHe2p3P2].Coll.cool = savecool;
00490                                 Transitions[ipHE_LIKE][nelem][ipHi][ipHe2p3P2].Coll.heat = saveheat;
00491                         }
00492                         for( ipLo=ipHe2p3P2+1; ipLo < iso.numPrintLevels[ipHE_LIKE][nelem]-1; ipLo++ )
00493                         {
00494                                 for( ipHi=ipLo+1; ipHi < iso.numPrintLevels[ipHE_LIKE][nelem]; ipHi++ )
00495                                 {
00496                                         /* skip non-radiative lines */
00497                                         if( Transitions[ipHE_LIKE][nelem][ipHi][ipLo].ipCont < 1 ) 
00498                                                 continue;
00499 
00500                                         /* find number of photons escaping cloud */
00501                                         Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->phots = 
00502                                                 Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->Aul*
00503                                                 StatesElem[ipHE_LIKE][nelem][ipHi].Pop*
00504                                                 Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->Pesc*
00505                                                 dense.xIonDense[nelem][nelem+1-ipISO];
00506 
00507                                         /* now find line intensity */
00508                                         /* >>chng 01 jan 15, put cast double to force double evaluation */
00509                                         Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->xIntensity = 
00510                                                 Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->phots*
00511                                                 Transitions[ipHE_LIKE][nelem][ipHi][ipLo].EnergyErg;
00512 
00513                                         if( iso.lgRandErrGen[ipISO] )
00514                                         {
00515                                                 Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->phots *=
00516                                                         iso.ErrorFactor[ipISO][nelem][ipHi][ipLo][IPRAD];
00517                                                 Transitions[ipHE_LIKE][nelem][ipHi][ipLo].Emis->xIntensity *= 
00518                                                         iso.ErrorFactor[ipISO][nelem][ipHi][ipLo][IPRAD];
00519                                         }
00520 
00521                                         /* this will enter .xIntensity into the line stack */
00522                                         PutLine(&Transitions[ipHE_LIKE][nelem][ipHi][ipLo],
00523                                                 "predicted line, all processes included");
00524                                 }
00525                         }
00526 
00527                         /* Now put the satellite lines in */
00528                         if( iso.lgDielRecom[ipISO] )
00529                                 DoSatelliteLines(nelem);
00530                 }
00531         }
00532 
00533         if( iso.n_HighestResolved_max[ipHE_LIKE][ipHELIUM] >= 4 &&
00534                 ( iso.n_HighestResolved_max[ipH_LIKE][ipHYDROGEN] + iso.nCollapsed_max[ipH_LIKE][ipHYDROGEN] ) >=8 )
00535         {
00536                 /* For info only, add the total photon flux in 3889 and 7065,
00537                 * and 3188, 4713, and 5876. */
00538                 photons_3889_plus_7065 =
00539                         /* to 2p3P2 */
00540                         Transitions[ipHE_LIKE][ipHELIUM][ipHe3s3S][ipHe2p3P2].Emis->xIntensity/
00541                         Transitions[ipHE_LIKE][ipHELIUM][ipHe3s3S][ipHe2p3P2].EnergyErg +
00542                         Transitions[ipHE_LIKE][ipHELIUM][ipHe3d3D][ipHe2p3P2].Emis->xIntensity/
00543                         Transitions[ipHE_LIKE][ipHELIUM][ipHe3d3D][ipHe2p3P2].EnergyErg +
00544                         Transitions[ipHE_LIKE][ipHELIUM][ipHe4s3S][ipHe2p3P2].Emis->xIntensity/
00545                         Transitions[ipHE_LIKE][ipHELIUM][ipHe4s3S][ipHe2p3P2].EnergyErg +
00546                         /* to 2p3P1 */
00547                         Transitions[ipHE_LIKE][ipHELIUM][ipHe3s3S][ipHe2p3P1].Emis->xIntensity/
00548                         Transitions[ipHE_LIKE][ipHELIUM][ipHe3s3S][ipHe2p3P1].EnergyErg +
00549                         Transitions[ipHE_LIKE][ipHELIUM][ipHe3d3D][ipHe2p3P1].Emis->xIntensity/
00550                         Transitions[ipHE_LIKE][ipHELIUM][ipHe3d3D][ipHe2p3P1].EnergyErg +
00551                         Transitions[ipHE_LIKE][ipHELIUM][ipHe4s3S][ipHe2p3P1].Emis->xIntensity/
00552                         Transitions[ipHE_LIKE][ipHELIUM][ipHe4s3S][ipHe2p3P1].EnergyErg +
00553                         /* to 2p3P0 */
00554                         Transitions[ipHE_LIKE][ipHELIUM][ipHe3s3S][ipHe2p3P0].Emis->xIntensity/
00555                         Transitions[ipHE_LIKE][ipHELIUM][ipHe3s3S][ipHe2p3P0].EnergyErg +
00556                         Transitions[ipHE_LIKE][ipHELIUM][ipHe3d3D][ipHe2p3P0].Emis->xIntensity/
00557                         Transitions[ipHE_LIKE][ipHELIUM][ipHe3d3D][ipHe2p3P0].EnergyErg +
00558                         Transitions[ipHE_LIKE][ipHELIUM][ipHe4s3S][ipHe2p3P0].Emis->xIntensity/
00559                         Transitions[ipHE_LIKE][ipHELIUM][ipHe4s3S][ipHe2p3P0].EnergyErg +
00560                         /* to 2s3S */
00561                         Transitions[ipHE_LIKE][ipHELIUM][ipHe3p3P][ipHe2s3S].Emis->xIntensity/
00562                         Transitions[ipHE_LIKE][ipHELIUM][ipHe3p3P][ipHe2s3S].EnergyErg +
00563                         Transitions[ipHE_LIKE][ipHELIUM][ipHe4p3P][ipHe2s3S].Emis->xIntensity/
00564                         Transitions[ipHE_LIKE][ipHELIUM][ipHe4p3P][ipHe2s3S].EnergyErg;
00565 
00566                 long upperIndexofH8 = iso.QuantumNumbers2Index[ipH_LIKE][ipHYDROGEN][8][1][2];
00567 
00568                 /* Add in photon flux of H8 3889 */
00569                 photons_3889_plus_7065 += 
00570                         Transitions[ipH_LIKE][ipHYDROGEN][upperIndexofH8][1].Emis->xIntensity/
00571                         Transitions[ipH_LIKE][ipHYDROGEN][upperIndexofH8][1].EnergyErg;
00572 
00573                 /* now multiply by ergs of normalization line, so that relative flux of
00574                 * this line will be ratio of photon fluxes. */
00575                 photons_3889_plus_7065 *= (ERG1CM*1.e8)/LineSave.WavLNorm;
00576                 linadd( photons_3889_plus_7065, 3889., "Pho+", 'i',
00577                         "photon sum given in Porter et al. 2007 (astro-ph/0611579)");
00578         }
00579 
00580         /* ====================================================
00581          * end helium
00582          * ====================================================*/
00583 
00584         if( trace.lgTrace )
00585         {
00586                 fprintf( ioQQQ, "   lines_helium returns\n" );
00587         }
00588         return;
00589 }
00590 
00591 
00592 STATIC void GetStandardHeLines(void)
00593 {
00594         FILE *ioDATA;
00595         bool lgEOL, lgHIT;
00596         long i, i1, i2, j, nelem;
00597 
00598 #       define chLine_LENGTH 1000
00599         char chLine[chLine_LENGTH];
00600 
00601         DEBUG_ENTRY( "GetStandardHeLines()" );
00602 
00603         if( trace.lgTrace )
00604                 fprintf( ioQQQ," lines_helium opening he1_case_ab.dat\n");
00605 
00606         ioDATA = open_data( "he1_case_ab.dat", "r" );
00607 
00608         /* check that magic number is ok */
00609         if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00610         {
00611                 fprintf( ioQQQ, " lines_helium could not read first line of he1_case_ab.dat.\n");
00612                 cdEXIT(EXIT_FAILURE);
00613         }
00614         i = 1;
00615         /* first number is magic number, second is number of lines in file      */
00616         i1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00617         i2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00618         NumLines = i2;
00619 
00620         /* the following is to check the numbers that appear at the start of he1_case_ab.dat */
00621         if( i1 !=CASEABMAGIC )
00622         {
00623                 fprintf( ioQQQ, 
00624                         " lines_helium: the version of he1_case_ab.dat is not the current version.\n" );
00625                 fprintf( ioQQQ, 
00626                         " lines_helium: I expected to find the number %i and got %li instead.\n" ,
00627                         CASEABMAGIC, i1 );
00628                 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
00629                 cdEXIT(EXIT_FAILURE);
00630         }
00631 
00632         /* get the array of temperatures */
00633         lgHIT = false;
00634         while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
00635         {
00636                 /* only look at lines without '#' in first col */
00637                 if( chLine[0] != '#')
00638                 {
00639                         lgHIT = true;
00640                         j = 0;
00641                         lgEOL = false;
00642                         i = 1;
00643                         while( !lgEOL && j < NUMTEMPS)
00644                         {
00645                                 CaABTemps[j] = FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00646                                 ++j;
00647                         }
00648                         break;
00649                 }
00650         }
00651 
00652         if( !lgHIT )
00653         {
00654                 fprintf( ioQQQ, " lines_helium could not find line of temperatures in he1_case_ab.dat.\n");
00655                 cdEXIT(EXIT_FAILURE);
00656         }
00657 
00658         /* create space for array of values, if not already done */
00659         CaABIntensity = (double ***)MALLOC(sizeof(double **)*(unsigned)LIMELM );
00660         /* create space for array of values, if not already done */
00661         CaABLines = (stdLines **)MALLOC(sizeof(stdLines *)*(unsigned)LIMELM );
00662 
00663         for( nelem=ipHELIUM; nelem<LIMELM; ++nelem )
00664         {
00668                 if( nelem != ipHELIUM )
00669                         continue;
00670 
00671                 /* only grab core for elements that are turned on */
00672                 if( nelem == ipHELIUM || dense.lgElmtOn[nelem])
00673                 {
00674                         /* create space for array of values, if not already done */
00675                         CaABIntensity[nelem] = (double **)MALLOC(sizeof(double *)*(unsigned)(i2) );
00676                         CaABLines[nelem] = (stdLines *)MALLOC(sizeof(stdLines )*(unsigned)(i2) );
00677 
00678                         /* avoid allocating 0 bytes, some OS return NULL pointer, PvH 
00679                         CaABIntensity[nelem][0] = NULL;*/
00680                         for( j = 0; j < i2; ++j )
00681                         {
00682                                 CaABIntensity[nelem][j] = (double *)MALLOC(sizeof(double)*(unsigned)NUMTEMPS );
00683 
00684                                 CaABLines[nelem][j].ipHi = -1;
00685                                 CaABLines[nelem][j].ipLo = -1;
00686                                 strcpy( CaABLines[nelem][j].label , "    " );
00687 
00688                                 for( i=0; i<NUMTEMPS; ++i )
00689                                 {
00690                                         CaABIntensity[nelem][j][i] = 0.;
00691                                 }
00692                         }
00693                 }
00694         }
00695 
00696         /* now read in the case A and B line data */
00697         lgHIT = false;
00698         nelem = ipHELIUM;
00699         while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
00700         {
00701                 static long line = 0;
00702                 char *chTemp;
00703 
00704                 /* only look at lines without '#' in first col */
00705                 if( (chLine[0] == ' ') || (chLine[0]=='\n') )
00706                         break;
00707                 if( chLine[0] != '#')
00708                 {
00709                         lgHIT = true;
00710 
00711                         /* get lower and upper level index */
00712                         j = 1;
00713                         /* the first number is the wavelength, which is not used
00714                          * for anything, but must skip over it. */
00715                         i1 = (long)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
00716                         CaABLines[nelem][line].ipLo = (long)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
00717                         CaABLines[nelem][line].ipHi = (long)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
00718 
00719                         ASSERT( CaABLines[nelem][line].ipLo < CaABLines[nelem][line].ipHi );
00720                         /*if( CaABLines[nelem][line].ipHi >= iso.numLevels_max[ipHE_LIKE][nelem] )
00721                                 continue;*/
00722 
00723                         chTemp = chLine;
00724                         /* skip over 4 tabs to start of data */
00725                         for( i=0; i<3; ++i )
00726                         {
00727                                 if( (chTemp = strchr( chTemp, '\t' )) == NULL )
00728                                 {
00729                                         fprintf( ioQQQ, " lines_helium no init case A and B\n" );
00730                                         cdEXIT(EXIT_FAILURE);
00731                                 }
00732                                 ++chTemp;
00733                         }
00734 
00735                         strncpy( CaABLines[nelem][line].label, chTemp , 4 );
00736                         CaABLines[nelem][line].label[4] = 0;
00737 
00738                         for( i=0; i<NUMTEMPS; ++i )
00739                         {
00740                                 double b;
00741                                 if( (chTemp = strchr( chTemp, '\t' )) == NULL )
00742                                 {
00743                                         fprintf( ioQQQ, " lines_helium could not scan case A and B lines, current indices: %li %li\n",
00744                                                 CaABLines[nelem][line].ipHi,
00745                                                 CaABLines[nelem][line].ipLo );
00746                                         cdEXIT(EXIT_FAILURE);
00747                                 }
00748                                 ++chTemp;
00749                                 sscanf( chTemp , "%le" , &b );
00750                                 CaABIntensity[nelem][line][i] = pow(10.,b);
00751                         }
00752                         line++;
00753                 }
00754         }
00755 
00756         /* close the data file */
00757         fclose( ioDATA );
00758         return;
00759 }
00760 
00762 STATIC double TempInterp2( double* TempArray , double* ValueArray, long NumElements )
00763 {
00764         long int ipTe=-1;
00765         double rate = 0.;
00766         long i0;
00767 
00768         DEBUG_ENTRY( "TempInterp2()" );
00769 
00770         ASSERT( fabs( 1. - (double)phycon.alogte/log10(phycon.te) ) < 0.0001 );
00771 
00772         /* te totally unknown */
00773         if( phycon.alogte <= TempArray[0] )
00774         {
00775                 return ValueArray[0];
00776         }
00777         else if( phycon.alogte >= TempArray[NumElements-1] )
00778         {
00779                 return ValueArray[NumElements-1];
00780         }
00781 
00782         /* now search for temperature */
00783         ipTe = hunt_bisect( TempArray, NumElements, phycon.alogte );                    
00784 
00785         ASSERT( (ipTe >=0) && (ipTe < NumElements-1)  );
00786 
00787         /* Do a four-point interpolation */
00788         const int ORDER = 3; /* order of the fitting polynomial */
00789         i0 = max(min(ipTe-ORDER/2,NumElements-ORDER-1),0);
00790         rate = lagrange( &TempArray[i0], &ValueArray[i0], ORDER+1, phycon.alogte );
00791 
00792         return rate;
00793 }
00794 
00796 /* For double-ionization discussions, see Lindsay, Rejoub, & Stebbings 2002     */
00797 /* Also read Itza-Ortiz, Godunov, Wang, and McGuire 2001.       */
00798 STATIC void DoSatelliteLines( long nelem )
00799 {
00800         long ipISO = ipHE_LIKE;
00801         
00802         DEBUG_ENTRY( "DoSatelliteLines()" );
00803 
00804         ASSERT( iso.lgDielRecom[ipISO] && dense.lgElmtOn[nelem] );
00805 
00806         for( long i=0; i<iso.numPrintLevels[ipISO][nelem]; i++ )
00807         {
00808                 double dr_rate = iso.DielecRecomb[ipISO][nelem][i];
00809 
00810                 SatelliteLines[ipISO][nelem][i].Emis->phots = 
00811                         dr_rate * dense.eden * dense.xIonDense[nelem][nelem+1-ipISO];
00812                 
00813                 SatelliteLines[ipISO][nelem][i].Emis->xIntensity = 
00814                         SatelliteLines[ipISO][nelem][i].Emis->phots * ERG1CM * SatelliteLines[ipISO][nelem][i].EnergyWN;
00815 
00816                 PutLine( &SatelliteLines[ipISO][nelem][i], "iso satellite line", "Sate" );
00817         
00818                 //linadd( SatelliteLines[ipISO][nelem][i].Emis->xIntensity ,
00819                 //      SatelliteLines[ipISO][nelem][i].WLAng, "Sate",'i', "iso satellite line" );
00820         }
00821 
00822         return;
00823 }
/home66/gary/public_html/cloudy/c08_branch/source/prt_lines_helium.cpp
