cloudy trunk

cont_createpointers.cpp

Go to the documentation of this file.
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 /*ContCreatePointers set up pointers for lines and continua called by cloudy after input read in 
00004  * and continuum mesh has been set */
00005 /*fiddle adjust energy bounds of certain cells so that match ionization edges exactly */
00006 /*ipShells assign continuum energy pointers to shells for all atoms,
00007  * called by ContCreatePointers */
00008 /*LimitSh sets upper energy limit to subshell integrations */
00009 /*ContBandsCreate - read set of continuum bands to enter total emission into line stack*/
00010 #include "cddefines.h"
00011 #include "physconst.h"
00012 #include "lines_service.h"
00013 #include "iso.h"
00014 #include "secondaries.h"
00015 #include "taulines.h"
00016 #include "elementnames.h"
00017 #include "ionbal.h"
00018 #include "rt.h"
00019 #include "opacity.h"
00020 #include "yield.h"
00021 #include "dense.h"
00022 #include "he.h"
00023 #include "fe.h"
00024 #include "rfield.h"
00025 #include "oxy.h"
00026 #include "atomfeii.h"
00027 #include "atoms.h"
00028 #include "trace.h"
00029 #include "hmi.h"
00030 #include "heavy.h"
00031 #include "predcont.h"
00032 #include "atmdat.h"
00033 #include "ipoint.h"
00034 #include "h2.h"
00035 #include "continuum.h"
00036 
00037 /*LimitSh sets upper energy limit to subshell integrations */
00038 STATIC long LimitSh(long int ion, 
00039   long int nshell, 
00040   long int nelem);
00041 
00042 STATIC void ipShells(
00043         /* nelem is the atomic number on the C scale, Li is 2 */
00044         long int nelem);
00045 
00046 /*ContBandsCreate - read set of continuum bands to enter total emission into line*/
00047 STATIC void ContBandsCreate(
00048         /* chFile is optional filename, if void then use default bands,
00049          * if not void then use file specified,
00050          * return value is 0 for success, 1 for failure */
00051          const char chFile[] );
00052 
00053 /* upper level for two-photon emission in H and He iso sequences */
00054 #define TwoS    (1+ipISO)
00055 
00056 /*fiddle adjust energy bounds of certain cells so that match ionization edges exactly */
00057 STATIC void fiddle(long int ipnt, 
00058   double exact);
00059 
00060 void ContCreatePointers(void)
00061 {
00062         char chLab[5];
00063         long int 
00064           i, 
00065           ion, 
00066           ipHi, 
00067           ipLo, 
00068           ipISO,
00069           iWL_Ang,
00070           j, 
00071           nelem,
00072           nshells;
00073         double energy,
00074                 xnew;
00075         /* counter to say whether pointers have ever been evaluated */
00076         static int nCalled = 0;
00077 
00078         DEBUG_ENTRY( "ContCreatePointers()" );
00079 
00080         /* create the hydrogen atom for this core load, routine creates space then zeros it out
00081          * on first call, on second and later calls it only zeros things out */
00082         iso_create();
00083 
00084         /* create internal static variables needed to do the H2 molecule */
00085         H2_Create();
00086 
00087         /* nCalled is local static variable defined 0 when defined. 
00088          * it is incremented below, so that space only allocated one time per coreload. */
00089         if( nCalled > 0 )
00090         {
00091                 if( trace.lgTrace )
00092                         fprintf( ioQQQ, " ContCreatePointers called, not evaluating.\n" );
00093                 return;
00094         }
00095         else
00096         {
00097                 if( trace.lgTrace )
00098                         fprintf( ioQQQ, " ContCreatePointers called first time.\n" );
00099                 ++nCalled;
00100         }
00101 
00102         for( i=0; i < rfield.nupper; i++ )
00103         {
00104                 /* this is array of labels for lines and continua, set to blanks at first */
00105                 strcpy( rfield.chContLabel[i], "    ");
00106                 strcpy( rfield.chLineLabel[i], "    ");
00107         }
00108 
00109         /* we will generate a set of array indices to ionization edges for
00110          * the first thirty elements.  First set all array indices to
00111          * totally bogus values so we will crash if misused */
00112         for( nelem=0; nelem<LIMELM; ++nelem )
00113         {
00114                 if( dense.lgElmtOn[nelem] )
00115                 {
00116                         for( ion=0; ion<LIMELM; ++ion )
00117                         {
00118                                 for( nshells=0; nshells<7; ++nshells )
00119                                 {
00120                                         for( j=0; j<3; ++j )
00121                                         {
00122                                                 opac.ipElement[nelem][ion][nshells][j] = INT_MIN;
00123                                         }
00124                                 }
00125                         }
00126                 }
00127         }
00128 
00129         /* pointer to excited state of O+2 */
00130         opac.ipo3exc[0] = ipContEnergy(3.855,"O3ex");
00131 
00132         /* main hydrogenic arrays - THIS OCCURS TWICE!! H and He here, then the
00133          * remaining hydrogenic species near the bottom.  This is so that H and HeII get
00134          * their labels stuffed into the arrays, and the rest of the hydrogenic series 
00135          * get whatever is left over after the level 1 lines.
00136          * to find second block, search on "ipZ=2" */
00137         /* NB note that no test for H or He being on exists here - we will always
00138          * define the H and He arrays even when He is off, since we need to
00139          * know where the he edges are for the bookkeeping that occurs in continuum
00140          * binning routines */
00141         /* this loop is over H, He-like only - DO NOT CHANGE TO NISO */
00142         for( ipISO=ipH_LIKE; ipISO<=ipHE_LIKE; ++ipISO )
00143         {
00144                 /* this will be over HI, HeII, then HeI only */
00145                 for( nelem=ipISO; nelem < 2; nelem++ )
00146                 {
00147                         /* generate label for this ion */
00148                         sprintf( chLab, "%2s%2ld",elementnames.chElementSym[nelem], nelem+1-ipISO );
00149 
00150                         /* array index for continuum edges for ground */
00151                         iso.ipIsoLevNIonCon[ipISO][nelem][0] = ipContEnergy(iso.xIsoLevNIonRyd[ipISO][nelem][0],chLab);
00152                         for( ipHi=1; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ )
00153                         {
00154                                 /* array index for continuum edges for excited levels */
00155                                 iso.ipIsoLevNIonCon[ipISO][nelem][ipHi] = ipContEnergy(iso.xIsoLevNIonRyd[ipISO][nelem][ipHi],chLab);
00156 
00157                                 /* define all line array indices */
00158                                 for( ipLo=0; ipLo < ipHi; ipLo++ )
00159                                 {
00160                                         if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul <= iso.SmallA )
00161                                                 continue;
00162 
00163                                         /* some lines have negative or zero energy */
00164                                         /* >>chng 03 apr 22, this check was if less than or equal to zero,
00165                                          * changed to lowest energy point so that ultra low energy transitions are
00166                                          * not considered.      */
00167                                         if( Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN * WAVNRYD < continuum.filbnd[0] )
00168                                                 continue;
00169 
00170                                         /* some energies are negative for inverted levels */
00171                                         Transitions[ipISO][nelem][ipHi][ipLo].ipCont = 
00172                                                 ipLineEnergy(Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN * WAVNRYD , chLab ,
00173                                                 iso.ipIsoLevNIonCon[ipISO][nelem][ipLo]);
00174                                         Transitions[ipISO][nelem][ipHi][ipLo].Emis->ipFine = 
00175                                                 ipFineCont(Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN * WAVNRYD );
00176                                         /* check that energy scales are the same, to within energy resolution of arrays */
00177 #                                       ifndef NDEBUG
00178                                         if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->ipFine > 0 )
00179                                         {
00180                                                 realnum anuCoarse = rfield.anu[Transitions[ipISO][nelem][ipHi][ipLo].ipCont-1];
00181                                                 realnum anuFine = rfield.fine_anu[Transitions[ipISO][nelem][ipHi][ipLo].Emis->ipFine];
00182                                                 realnum widCoarse = rfield.widflx[Transitions[ipISO][nelem][ipHi][ipLo].ipCont-1];
00183                                                 realnum widFine = anuFine - rfield.fine_anu[Transitions[ipISO][nelem][ipHi][ipLo].Emis->ipFine-1];
00184                                                 realnum width = MAX2( widFine , widCoarse );
00185                                                 /* NB - can only assert more than width of coarse cell 
00186                                                  * since close to ionization limit, coarse index is 
00187                                                  * kept below next ionization edge 
00188                                                  * >>chng 05 mar 02, pre factor below had been 1.5, chng to 2
00189                                                  * tripped near H grnd photo threshold */
00190                                                 ASSERT( fabs(anuCoarse - anuFine) / anuCoarse < 
00191                                                         2.*width/anuCoarse);
00192                                         }
00193 #                                       endif
00194                                 }
00195                         }/* ipHi loop */
00196                 }/* nelem loop */
00197         }/* ipISO */
00198 
00199         /* need to increase the cell for the HeII balmer continuum by one, so that
00200          * hydrogen Lyman continuum will pick it up. */
00201         nelem = ipHELIUM;
00202         /* copy label over to next higher cell */
00203         if( strcmp( rfield.chContLabel[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH2s]-1] , "He 2" ) == 0)
00204         {
00205                 strcpy( rfield.chContLabel[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH2s]], 
00206                                  rfield.chContLabel[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH2s]-1] );
00207                 /* set previous spot to blank */
00208                 strcpy( rfield.chContLabel[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH2s]-1] , "    ");
00209         }
00210         else if( strcmp( rfield.chContLabel[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH2s]-1] , "H  1" ) == 0)
00211         {
00212                 /* set previous spot to blank */
00213                 strcpy( rfield.chContLabel[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH2s]] , "He 2");
00214         }
00215         else
00216         {
00217                 fprintf(ioQQQ," insanity heii pointer fix contcreatepointers\n");
00218         }
00219         /* finally increment the two HeII pointers so that they are above the Lyman continuum */
00220         ++iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH2s];
00221         ++iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH2p];
00222 
00223         /* we will define an array of either 1 or 0 to show whether photooelectron
00224          * energy is large enough to produce secondary ionizations
00225          * below 100eV no secondary ionization - this is the threshold*/
00226         secondaries.ipSecIon = ipoint(7.353);
00227 
00228         /* this is highest energy where k-shell opacities are counted
00229          * can be adjusted with "set kshell" command */
00230         continuum.KshellLimit = ipoint(continuum.EnergyKshell);
00231 
00232         /* pointers for molecules
00233          * H2+ dissociation energy 2.647 eV but cs small below 0.638 Ryd */
00234         opac.ih2pnt[0] = ipContEnergy(0.502,"H2+ ");
00235         opac.ih2pnt[1] = ipoint(1.03);
00236 
00237         /* pointers for most prominent PAH features
00238          * energies given to ipContEnergy are only to put lave in the right place
00239          * wavelengths are rough observed values of blends
00240          * 7.6 microns */
00241         i = ipContEnergy(0.0117, "PAH " );
00242 
00243         /* feature near 6.2 microns */
00244         i = ipContEnergy(0.0147, "PAH " );
00245 
00246         /* 3.3 microns */
00247         i = ipContEnergy(0.028, "PAH " );
00248 
00249         /* 11.2 microns */
00250         i = ipContEnergy(0.0080, "PAH " );
00251 
00252         /* 12.3 microns */
00253         i = ipContEnergy(0.0077, "PAH " );
00254 
00255         /* 13.5 microns */
00256         i = ipContEnergy(0.0069, "PAH " );
00257 
00258 
00259         /* fix pointers for hydrogen and helium */
00260 
00261         /* pointer to Bowen 374A resonance line */
00262         he.ip374 = ipLineEnergy(1.92,"He 2",0);
00263         he.ip660 = ipLineEnergy(1.38,"He 2",0);
00264 
00265         /* set up series of continuum pointers to be used in routine lines to
00266          * enter continuum points into line array */
00267         for( i=0; i < NPREDCONT; i++ )
00268         {
00269                 /* EnrPredCont contains array of energy points, as set in zerologic.c */
00270                 ipPredCont[i] = ipoint(EnrPredCont[i]) - 1;
00271         }
00272 
00273         /* pointer to energy defining effect x-ray column */
00274         rt.ipxry = ipoint(73.5);
00275 
00276         /* pointer to Hminus edge at 0.754eV */
00277         hmi.iphmin = ipContEnergy(0.05544,"H-  ");
00278 
00279         /* pointer to threshold for H2 photoionization at 15.4 eV */
00280         opac.ipH2_photo_thresh = ipContEnergy( 15.4/EVRYD , "H2  ");
00281 
00282         hmi.iheh1 = ipoint(1.6);
00283         hmi.iheh2 = ipoint(2.3);
00284 
00285         /* pointer to carbon k-shell ionization */
00286         opac.ipCKshell = ipoint(20.6);
00287 
00288         /* pointer to threshold for pair production */
00289         opac.ippr = ipoint(7.51155e4) + 1;
00290 
00291         /* pointer to x-ray - gamma ray bound; 100 kev */
00292         rfield.ipEnerGammaRay = ipoint(rfield.EnerGammaRay);
00293 
00294         t_fe2ovr_la::Inst().init_pointers();
00295 
00296         /* make low energy edges of some cells exact,
00297          * index is on c scale
00298          * 0.99946 is correct energy of hydrogen in inf mass ryd */
00299         fiddle(iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1,0.99946);
00300         fiddle(iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH2p]-1,0.24986);
00301         /* confirm that labels are in correct location */
00302         ASSERT( strcmp( rfield.chContLabel[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s]-1], "H  1" ) ==0 );
00303         ASSERT( strcmp( rfield.chContLabel[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH2p]-1], "H  1" ) ==0 );
00304 
00305         fiddle(iso.ipIsoLevNIonCon[ipH_LIKE][ipHELIUM][ipH1s]-1,4.00000);
00306         ASSERT( strcmp( rfield.chContLabel[iso.ipIsoLevNIonCon[ipH_LIKE][ipHELIUM][ipH1s]-1], "He 2" ) ==0 );
00307 
00308         /* pointer to excited state of O+2 */
00309         fiddle(opac.ipo3exc[0]-1,3.855);
00310 
00311         /* now fix widflx array so that it is correct */
00312         for( i=1; i<rfield.nupper-1; ++i )
00313         {
00314                 rfield.widflx[i] = ((rfield.anu[i+1] - rfield.anu[i]) + (rfield.anu[i] - rfield.anu[i-1]))/2.f;
00315         }
00316 
00317         /* these are indices for centers of B and V filters,
00318          * taken from table on page 202 of Allen, AQ, 3rd ed */
00319         /* the B filter array offset */
00320         rfield.ipB_filter = ipoint( RYDLAM / WL_B_FILT );
00321         /* the V filter array offset */
00322         rfield.ipV_filter = ipoint( RYDLAM / WL_V_FILT );
00323 
00324         /* these are the lower and upper bounds for the G0 radiation field
00325          * used by Tielens & Hollenbach in their PDR work */
00326         rfield.ipG0_TH85_lo =  ipoint(  6.0 / EVRYD );
00327         rfield.ipG0_TH85_hi =  ipoint( 13.6 / EVRYD );
00328 
00329         /* this is the limits for Draine & Bertoldi Habing field */
00330         rfield.ipG0_DB96_lo =  ipoint(  RYDLAM / 1110. );
00331         rfield.ipG0_DB96_hi =  ipoint( RYDLAM / 912. );
00332 
00333         /* this is special form of G0 that could be used in some future version, for now,
00334          * use default, TH85 */
00335         rfield.ipG0_spec_lo = ipoint(  6.0 / EVRYD );
00336         rfield.ipG0_spec_hi = ipoint( RYDLAM / 912. );
00337 
00338         /* this index is to 1000A to obtain the extinction at 1000A */
00339         rfield.ip1000A = ipoint( RYDLAM / 1000. );
00340 
00341         /* now save current form of array */
00342         for( i=0; i < rfield.nupper; i++ )
00343         {
00344                 rfield.AnuOrg[i] = rfield.anu[i];
00345                 rfield.anusqr[i] = (realnum)sqrt(rfield.AnuOrg[i]);
00346         }
00347 
00348         /* following order of elements is in roughly decreasing abundance
00349          * when ipShells gets a cell for the valence IP it does so through
00350          * ipContEnergy, which makes sure that no two ions get the same cell
00351          * earliest elements have most precise ip mapping */
00352 
00353         /* set up shell pointers for hydrogen */
00354         nelem = ipHYDROGEN;
00355         ion = 0;
00356 
00357         /* the number of shells */
00358         Heavy.nsShells[nelem][0] = 1;
00359 
00360         /*pointer to ionization threshold in energy array*/
00361         Heavy.ipHeavy[nelem][ion] = iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH1s];
00362         opac.ipElement[nelem][ion][0][0] = iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH1s];
00363 
00364         /* upper limit to energy integration */
00365         opac.ipElement[nelem][ion][0][1] = rfield.nupper;
00366 
00367         if( dense.lgElmtOn[ipHELIUM] )
00368         {
00369                 /* set up shell pointers for helium */
00370                 nelem = ipHELIUM;
00371                 ion = 0;
00372 
00373                 /* the number of shells */
00374                 Heavy.nsShells[nelem][0] = 1;
00375 
00376                 /*pointer to ionization threshold in energy array*/
00377                 Heavy.ipHeavy[nelem][ion] = iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0];
00378                 opac.ipElement[nelem][ion][0][0] = iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0];
00379 
00380                 /* upper limit to energy integration */
00381                 opac.ipElement[nelem][ion][0][1] = rfield.nupper;
00382 
00383                 /* (hydrogenic) helium ion */
00384                 ion = 1;
00385                 /* the number of shells */
00386                 Heavy.nsShells[nelem][1] = 1;
00387 
00388                 /*pointer to ionization threshold in energy array*/
00389                 Heavy.ipHeavy[nelem][ion] = iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH1s];
00390                 opac.ipElement[nelem][ion][0][0] = iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH1s];
00391 
00392                 /* upper limit to energy integration */
00393                 opac.ipElement[nelem][ion][0][1] = rfield.nupper;
00394         }
00395 
00396         /* check that ionization potential of neutral carbon valence shell is
00397          * positive */
00398         ASSERT( t_ADfA::Inst().ph1(2,5,5,0) > 0. );
00399 
00400         /* now fill in all sub-shell ionization array indices for elements heavier than He,
00401          * this must be done after previous loop on iso.ipIsoLevNIonCon[ipH_LIKE] since hydrogenic species use
00402          * iso.ipIsoLevNIonCon[ipH_LIKE] rather than ipoint in getting array index within continuum array */
00403         for( i=NISO; i<LIMELM; ++i )
00404         {
00405                 if( dense.lgElmtOn[i])
00406                 {
00407                         /* i is the atomic number on the c scale, 2 for Li */
00408                         ipShells(i);
00409                 }
00410         }
00411 
00412         /* most of these are set in ipShells, but not h-like or he-like, so do these here */
00413         Heavy.Valence_IP_Ryd[ipHYDROGEN][0] = t_ADfA::Inst().ph1(0,0,ipHYDROGEN,0)/EVRYD* 0.9998787;
00414         Heavy.Valence_IP_Ryd[ipHELIUM][0] = t_ADfA::Inst().ph1(0,1,ipHELIUM,0)/EVRYD* 0.9998787;
00415         Heavy.Valence_IP_Ryd[ipHELIUM][1] = t_ADfA::Inst().ph1(0,0,ipHELIUM,0)/EVRYD* 0.9998787;
00416         for( nelem=2; nelem<LIMELM; ++nelem )
00417         {
00418                 Heavy.Valence_IP_Ryd[nelem][nelem-1] = t_ADfA::Inst().ph1(0,1,nelem,0)/EVRYD* 0.9998787;
00419                 Heavy.Valence_IP_Ryd[nelem][nelem] = t_ADfA::Inst().ph1(0,0,nelem,0)/EVRYD* 0.9998787;
00420                 if( dense.lgElmtOn[nelem])
00421                 {
00422                         /* now confirm that all are properly set */
00423                         for( j=0; j<=nelem; ++j )
00424                         {
00425                                 ASSERT( Heavy.Valence_IP_Ryd[nelem][j] > 0.05 );
00426                         }
00427                         for( j=0; j<nelem; ++j )
00428                         {
00429                                 ASSERT( Heavy.Valence_IP_Ryd[nelem][j] < Heavy.Valence_IP_Ryd[nelem][j+1]);
00430                         }
00431                 }
00432         }
00433 
00434         /* array indices to bound Compton electron recoil ionization of all elements */
00435         for( nelem=0; nelem<LIMELM; ++nelem )
00436         {
00437                 if( dense.lgElmtOn[nelem])
00438                 {
00439                         for( ion=0; ion<nelem+1; ++ion )
00440                         {
00441                                 /* this is the threshold energy to Compton ionize valence shell electrons */
00442                                 energy = sqrt( Heavy.Valence_IP_Ryd[nelem][ion] * EN1RYD * ELECTRON_MASS * SPEEDLIGHT * SPEEDLIGHT ) / EN1RYD;
00443                                 /* the array index for this energy */
00444                                 ionbal.ipCompRecoil[nelem][ion] = ipoint( energy );
00445                         }
00446                 }
00447         }
00448 
00449         /* oxygen pointers for excited states
00450          * IO3 is pointer to O++ exc state, is set above */
00451         oxy.i2d = ipoint(1.242);
00452         oxy.i2p = ipoint(1.367);
00453         opac.ipo1exc[0] = ipContEnergy(0.856,"O1ex");
00454         opac.ipo1exc[1] = ipoint(2.0);
00455 
00456         /* upper limit for excited state photoionization
00457          * do not use ipContEnergy since only upper limit */
00458         opac.ipo3exc[1] = ipoint(5.0);
00459 
00460         /* upper level of 4363 */
00461         opac.ipo3exc3[0] = ipContEnergy(3.646,"O3ex");
00462         opac.ipo3exc3[1] = ipoint(5.0);
00463 
00464         /* following are various pointers for OI - Ly-beta pump problem
00465          * these are delta energies for Boltzmann factors */
00470         atoms.ipoiex[0] = ipoint(0.7005);
00471         atoms.ipoiex[1] = ipoint(0.10791);
00472         atoms.ipoiex[2] = ipoint(0.06925);
00473         atoms.ipoiex[3] = ipoint(0.01151);
00474         atoms.ipoiex[4] = ipoint(0.01999);
00475 
00476         /* >>chng 97 jan 27, move nitrogen after oxygen so that O gets the
00477          * most accurate pointers
00478          * Nitrogen
00479          * in1(1) is thresh for photo from excited state */
00480         opac.in1[0] = ipContEnergy(0.893,"N1ex");
00481 
00482         /* upper limit */
00483         opac.in1[1] = ipoint(2.);
00484 
00485         if( (trace.lgTrace && trace.lgConBug) || (trace.lgTrace && trace.lgPointBug) )
00486         {
00487                 fprintf( ioQQQ, "   ContCreatePointers:%ld energy cells used. N(1R):%4ld N(1.8):%4ld  N(4Ryd):%4ld N(O3)%4ld  N(x-ray):%5ld N(rcoil)%5ld\n", 
00488                   rfield.nupper, 
00489                   iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][ipH1s], 
00490                   iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][ipH1s], 
00491                   iso.ipIsoLevNIonCon[ipH_LIKE][ipHELIUM][ipH1s], 
00492                   opac.ipo3exc[0], 
00493                   opac.ipCKshell, 
00494                   ionbal.ipCompRecoil[ipHYDROGEN][0] );
00495 
00496 
00497                 fprintf( ioQQQ, "   ContCreatePointers: ipEnerGammaRay: %5ld IPPRpari produc%5ld\n", 
00498                   rfield.ipEnerGammaRay, opac.ippr );
00499 
00500                 fprintf( ioQQQ, "   ContCreatePointers: H pointers;" );
00501                 for( i=0; i <= 6; i++ )
00502                 {
00503                         fprintf( ioQQQ, "%4ld%4ld", i, iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][i] );
00504                 }
00505                 fprintf( ioQQQ, "\n" );
00506 
00507                 fprintf( ioQQQ, "   ContCreatePointers: Oxy pnters;" );
00508 
00509                 for( i=1; i <= 8; i++ )
00510                 {
00511                         fprintf( ioQQQ, "%4ld%4ld", i, Heavy.ipHeavy[ipOXYGEN][i-1] );
00512                 }
00513                 fprintf( ioQQQ, "\n" );
00514         }
00515 
00516         /* Magnesium
00517          * following is energy for phot of MG+ from exc state producing 2798 */
00518         opac.ipmgex = ipoint(0.779);
00519 
00520         /* lower, upper edges of Ca+ excited term photoionizaiton */
00521         opac.ica2ex[0] = ipContEnergy(0.72,"Ca2x");
00522         opac.ica2ex[1] = ipoint(1.);
00523 
00524         /* set up factors and pointers for Fe continuum fluorescence */
00525         fe.ipfe10 = ipoint(2.605);
00526 
00527         /* following is WL(CM)**2/(8PI) * conv fac for RYD to NU *A21 */
00528         fe.pfe10 = (realnum)(2.00e-18/rfield.widflx[fe.ipfe10-1]);
00529 
00530         /* this is 353 pump, f=0.032 */
00531         fe.pfe11a = (realnum)(4.54e-19/rfield.widflx[fe.ipfe10-1]);
00532 
00533         /* this is 341.1 f=0.012 */
00534         fe.pfe11b = (realnum)(2.75e-19/rfield.widflx[fe.ipfe10-1]);
00535         fe.pfe14 = (realnum)(1.15e-18/rfield.widflx[fe.ipfe10-1]);
00536 
00537         /* set up energy pointers for line optical depth arrays
00538          * this also increments flux, sets other parameters for lines */
00539 
00540         /* level 1 heavy elements line array */
00541         for( i=1; i <= nLevel1; i++ )
00542         {
00543                 /* put null terminated line label into chLab using line array*/
00544                 chIonLbl(chLab, &TauLines[i]);
00545 
00546                 TauLines[i].ipCont = 
00547                         ipLineEnergy(TauLines[i].EnergyWN * WAVNRYD, chLab ,0);
00548                 TauLines[i].Emis->ipFine = 
00549                         ipFineCont(TauLines[i].EnergyWN * WAVNRYD );
00550                 /* for debugging pointer - index on f not c scale, 
00551                  * this will find all lines that entered a given cell 
00552                    if( TauLines[i].ipCont==603 )
00553                         fprintf(ioQQQ,"( level1 %s\n", chLab);*/
00554 
00555                 if( TauLines[i].Emis->gf > 0. )
00556                 {
00557                         /* the gf value was entered
00558                          * derive the A, call to function is gf, wl (A), g_up */
00559                         TauLines[i].Emis->Aul = 
00560                                 (realnum)(eina(TauLines[i].Emis->gf,
00561                           TauLines[i].EnergyWN,TauLines[i].Hi->g));
00562                 }
00563                 else if( TauLines[i].Emis->Aul > 0. )
00564                 {
00565                         /* the Einstein A value was entered
00566                          * derive the gf, call to function is A, wl (A), g_up */
00567                         TauLines[i].Emis->gf = 
00568                                 (realnum)(GetGF(TauLines[i].Emis->Aul,
00569                           TauLines[i].EnergyWN,TauLines[i].Hi->g));
00570                 }
00571                 else
00572                 {
00573                         fprintf( ioQQQ, " level 1 line does not have valid gf or A\n" );
00574                         fprintf( ioQQQ, " This is ContCreatePointers\n" );
00575                         cdEXIT(EXIT_FAILURE);
00576                 }
00577 
00578                 /* used to get damping constant */
00579                 TauLines[i].Emis->dampXvel = 
00580                         (realnum)(TauLines[i].Emis->Aul / TauLines[i].EnergyWN/PI4);
00581 
00582                 /* derive the abs coefficient, call to function is gf, wl (A), g_low */
00583                 TauLines[i].Emis->opacity = 
00584                         (realnum)(abscf(TauLines[i].Emis->gf,
00585                   TauLines[i].EnergyWN,TauLines[i].Lo->g));
00586                 /*fprintf(ioQQQ,"TauLinesss\t%li\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
00587                         i,TauLines[i].Emis->opacity , TauLines[i].Emis->gf , TauLines[i].Emis->Aul ,TauLines[i].EnergyWN,TauLines[i].Lo->g);*/
00588 
00589                 /* excitation energy of transition in degrees kelvin */
00590                 TauLines[i].EnergyK = 
00591                         (realnum)(T1CM)*TauLines[i].EnergyWN;
00592 
00593                 /* energy of photon in ergs */
00594                 TauLines[i].EnergyErg = 
00595                         (realnum)(ERG1CM)*TauLines[i].EnergyWN;
00596 
00597                 /*line wavelength must be gt 0 */
00598                 ASSERT( TauLines[i].WLAng > 0 );
00599         }
00600 
00601         /*Beginning of the atmolEmis*/
00602         for( i=0; i <linesAdded2; i++)
00603         {
00604                 atmolEmis[i].gf = (realnum)GetGF(atmolEmis[i].Aul,
00605                         atmolEmis[i].tran->EnergyWN,atmolEmis[i].tran->Hi->g);
00606                 atmolEmis[i].dampXvel = (realnum)(atmolEmis[i].Aul/
00607                                         atmolEmis[i].tran->EnergyWN/PI4);
00608                 atmolEmis[i].damp = -1000.0;
00609                 /*Put null terminated line label into chLab*/
00610                 strncpy(chLab,atmolEmis[i].tran->Hi->chLabel,4);
00611                 chLab[4]='\0';
00612 
00613                 atmolEmis[i].tran->ipCont = 
00614                         ipLineEnergy(atmolEmis[i].tran->EnergyWN * WAVNRYD, chLab ,0);
00615                 atmolEmis[i].ipFine = ipFineCont(atmolEmis[i].tran->EnergyWN * WAVNRYD );
00616                 /* derive the abs coefficient, call to function is gf, wl (A), g_low */
00617                 atmolEmis[i].opacity = 
00618                         (realnum)(abscf(atmolEmis[i].gf,atmolEmis[i].tran->EnergyWN,
00619                                 atmolEmis[i].tran->Lo->g));
00620                 /* excitation energy of in degrees kelvin */
00621                 atmolEmis[i].tran->EnergyK = 
00622                         (realnum)(T1CM)*atmolEmis[i].tran->EnergyWN;
00623                 /* energy of photon in ergs */
00624                 atmolEmis[i].tran->EnergyErg = 
00625                         (realnum)(ERG1CM)*atmolEmis[i].tran->EnergyWN ;                                                           
00626         }
00627         /*end of the atmolEmis*/
00628 
00629         /* set the ipCont struc element for the H2 molecule */
00630         H2_ContPoint();
00631 
00632         for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00633         {
00634                 /* do remaining part of the he iso sequence */
00635                 for( nelem=2; nelem < LIMELM; nelem++ )
00636                 {
00637                         if( dense.lgElmtOn[nelem])
00638                         {
00639                                 /* generate label for this ion */
00640                                 sprintf( chLab, "%2s%2ld",elementnames.chElementSym[nelem], nelem+1-ipISO );
00641                                 /* Lya itself is the only transition below n=3 - we explicitly do not
00642                                  * want to generate pointers for 2s-1s or 2p-2s */
00643                                 /* array index for continuum edges for levels in He-like ions  */
00644                                 iso.ipIsoLevNIonCon[ipISO][nelem][0] = 
00645                                         ipContEnergy(iso.xIsoLevNIonRyd[ipISO][nelem][0],chLab);
00646 
00647                                 for( ipHi=1; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ )
00648                                 {
00649                                         /* array index for continuum edges for levels in He-like ions  */
00650                                         iso.ipIsoLevNIonCon[ipISO][nelem][ipHi] = ipContEnergy(iso.xIsoLevNIonRyd[ipISO][nelem][ipHi],chLab);
00651 
00652                                         /* define all he-like line pointers */
00653                                         for( ipLo=0; ipLo < ipHi; ipLo++ )
00654                                         {
00655 
00656                                                 if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul <= iso.SmallA )
00657                                                         continue;
00658 
00659                                                 /* some lines have negative or zero energy */
00660                                                 /* >>chng 03 apr 22, this check was if less than or equal to zero,
00661                                                  * changed to lowest energy point so that ultra low energy transitions are
00662                                                  * not considered.      */
00663                                                 if( Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN * WAVNRYD < continuum.filbnd[0] )
00664                                                         continue;
00665 
00666                                                 Transitions[ipISO][nelem][ipHi][ipLo].ipCont = 
00667                                                         ipLineEnergy(Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN * WAVNRYD , chLab ,
00668                                                         iso.ipIsoLevNIonCon[ipISO][nelem][ipLo]);
00669                                                 Transitions[ipISO][nelem][ipHi][ipLo].Emis->ipFine = 
00670                                                         ipFineCont(Transitions[ipISO][nelem][ipHi][ipLo].EnergyWN * WAVNRYD );
00671                                         }
00672                                 }
00673                                 iso.ipIsoLevNIonCon[ipISO][nelem][0] = ipContEnergy(iso.xIsoLevNIonRyd[ipISO][nelem][0],chLab);
00674                         }
00675                 }
00676         }
00677         for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00678         {
00679                 /* this will be over HI, HeII, then HeI only */
00680                 for( nelem=ipISO; nelem < LIMELM; nelem++ )
00681                 {
00682                         if( dense.lgElmtOn[nelem])
00683                         {
00684                                 ipLo = 0;
00685                                 /* these are the extra Lyman lines */
00686                                 for( ipHi=2; ipHi < iso.nLyman[ipISO]; ipHi++ )
00687                                 {
00688                                         /* some energies are negative for inverted levels */
00689                                         /*lint -e772 chLab not initialized */
00690                                         ExtraLymanLines[ipISO][nelem][ipHi].ipCont = 
00691                                                 ipLineEnergy(ExtraLymanLines[ipISO][nelem][ipHi].EnergyWN * WAVNRYD , chLab ,
00692                                                 iso.ipIsoLevNIonCon[ipISO][nelem][ipLo]);
00693 
00694                                         ExtraLymanLines[ipISO][nelem][ipHi].Emis->ipFine = 
00695                                                 ipFineCont(ExtraLymanLines[ipISO][nelem][ipHi].EnergyWN * WAVNRYD );
00696                                         /*lint +e772 chLab not initialized */
00697                                 }
00698 
00699                                 if( iso.lgDielRecom[ipISO] )
00700                                 {
00701                                         for( ipHi=0; ipHi<iso.numLevels_max[ipISO][nelem]; ipHi++ )
00702                                         {
00703                                                 
00704                                                 SatelliteLines[ipISO][nelem][ipHi].ipCont = ipLineEnergy(
00705                                                         SatelliteLines[ipISO][nelem][ipHi].EnergyWN * WAVNRYD , chLab , 
00706                                                         iso.ipIsoLevNIonCon[ipISO][nelem][0]);
00707 
00708                                                 SatelliteLines[ipISO][nelem][ipHi].Emis->ipFine =  
00709                                                         ipFineCont(SatelliteLines[ipISO][nelem][ipHi].EnergyWN * WAVNRYD );
00710                                         }
00711                                 }
00712                         }
00713                 }
00714         }
00715 
00716         /* some lines do not exist, the flag indicating this is ipCont == -1 */
00717         /* for H-like sequence, these are 2p-2s (energies degenerate) and 2s-1s, two photon */
00718         ipISO = ipHYDROGEN;
00719         /* do remaining part of the he iso sequence */
00720         for( nelem=ipISO; nelem < LIMELM; nelem++ )
00721         {
00722                 if( dense.lgElmtOn[nelem])
00723                 {
00724                         /* for H-like sequence want 2p-2s (energies degenerate) and 2s-1s, two photon */
00725                         Transitions[ipISO][nelem][ipH2p][ipH2s].ipCont = -1;
00726                         //Transitions[ipISO][nelem][ipH2p][ipH2s].Emis->ipFine = -1;
00727                         Transitions[ipISO][nelem][ipH2s][ipH1s].ipCont = -1;
00728                         //Transitions[ipISO][nelem][ipH2s][ipH1s].Emis->ipFine = -1;
00729                 }
00730         }
00731 
00732         fixit(); /* is this redundant? */
00733         /* for He-like sequence the majority of the transitions are bogus - A set to special value in this case */
00734         ipISO = ipHELIUM;
00735         /* do remaining part of the he iso sequence */
00736         for( nelem=ipISO; nelem < LIMELM; nelem++ )
00737         {
00738                 if( dense.lgElmtOn[nelem])
00739                 {
00740                         /* this is the two photon transition in the singlets */
00741                         Transitions[ipISO][nelem][ipHe2s1S][ipHe1s1S].ipCont = -1;
00742                         Transitions[ipISO][nelem][ipHe2s1S][ipHe1s1S].Emis->ipFine = -1;
00743 
00744                         for( ipHi=1; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ )
00745                         {
00746                                 for( ipLo=0; ipLo < ipHi; ipLo++ )
00747                                 {
00748                                         if( Transitions[ipISO][nelem][ipHi][ipLo].ipCont <= 0 )
00749                                                 continue;
00750 
00751                                         if( fabs(Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul - iso.SmallA) < SMALLFLOAT )
00752                                         {
00753                                                 /* iso.SmallA is value assigned to bogus transitions */
00754                                                 Transitions[ipISO][nelem][ipHi][ipLo].ipCont = -1;
00755                                                 Transitions[ipISO][nelem][ipHi][ipLo].Emis->ipFine = -1;
00756                                         }
00757                                 }
00758                         }
00759                 }
00760         }
00761 
00762         /* co carbon monoxide line array */
00763         for( i=0; i < nCORotate; i++ )
00764         {
00765                 /* excitation energy of transition in degrees kelvin */
00766                 C12O16Rotate[i].EnergyK = 
00767                         (realnum)(T1CM)*C12O16Rotate[i].EnergyWN;
00768                 C13O16Rotate[i].EnergyK = 
00769                         (realnum)(T1CM)*C13O16Rotate[i].EnergyWN;
00770 
00771                 /* energy of photon in ergs */
00772                 C12O16Rotate[i].EnergyErg = 
00773                         (realnum)(ERG1CM)*C12O16Rotate[i].EnergyWN;
00774                 C13O16Rotate[i].EnergyErg = 
00775                         (realnum)(ERG1CM)*C13O16Rotate[i].EnergyWN;
00776 
00777                 /* put null terminated line label into chLab using line array*/
00778                 chIonLbl(chLab, &C12O16Rotate[i]);
00779                 chIonLbl(chLab, &C13O16Rotate[i]);
00780 
00781                 C12O16Rotate[i].ipCont = 
00782                         ipLineEnergy(C12O16Rotate[i].EnergyWN * WAVNRYD, "12CO" ,0);
00783                 C12O16Rotate[i].Emis->ipFine = 
00784                         ipFineCont(C12O16Rotate[i].EnergyWN * WAVNRYD );
00785                 C13O16Rotate[i].ipCont = 
00786                         ipLineEnergy(C13O16Rotate[i].EnergyWN * WAVNRYD, "13CO" ,0);
00787                 C13O16Rotate[i].Emis->ipFine = 
00788                         ipFineCont(C13O16Rotate[i].EnergyWN * WAVNRYD );
00789 
00790                 if( C12O16Rotate[i].Emis->gf > 0. )
00791                 {
00792                         /* the gf value was entered
00793                          * derive the A, call to function is gf, wl (A), g_up */
00794                         C12O16Rotate[i].Emis->Aul = 
00795                                 (realnum)(eina(C12O16Rotate[i].Emis->gf,
00796                           C12O16Rotate[i].EnergyWN,C12O16Rotate[i].Hi->g));
00797                 }
00798                 else if( C12O16Rotate[i].Emis->Aul > 0. )
00799                 {
00800                         /* the Einstein A value was entered
00801                          * derive the gf, call to function is A, wl (A), g_up */
00802                         C12O16Rotate[i].Emis->gf = 
00803                                 (realnum)(GetGF(C12O16Rotate[i].Emis->Aul,
00804                           C12O16Rotate[i].EnergyWN,C12O16Rotate[i].Hi->g));
00805                 }
00806                 else
00807                 {
00808                         fprintf( ioQQQ, " 12CO line does not have valid gf or A\n" );
00809                         fprintf( ioQQQ, " This is ContCreatePointers\n" );
00810                         cdEXIT(EXIT_FAILURE);
00811                 }
00812                 if( C13O16Rotate[i].Emis->gf > 0. )
00813                 {
00814                         /* the gf value was entered
00815                          * derive the A, call to function is gf, wl (A), g_up */
00816                         C13O16Rotate[i].Emis->Aul = 
00817                                 (realnum)(eina(C13O16Rotate[i].Emis->gf,
00818                           C13O16Rotate[i].EnergyWN,C13O16Rotate[i].Hi->g));
00819                 }
00820                 else if( C13O16Rotate[i].Emis->Aul > 0. )
00821                 {
00822                         /* the Einstein A value was entered
00823                          * derive the gf, call to function is A, wl (A), g_up */
00824                         C13O16Rotate[i].Emis->gf = 
00825                                 (realnum)(GetGF(C13O16Rotate[i].Emis->Aul,
00826                           C13O16Rotate[i].EnergyWN,C13O16Rotate[i].Hi->g));
00827                 }
00828                 else
00829                 {
00830                         fprintf( ioQQQ, " 13CO line does not have valid gf or A\n" );
00831                         fprintf( ioQQQ, " This is ContCreatePointers\n" );
00832                         cdEXIT(EXIT_FAILURE);
00833                 }
00834 
00835                 /*line wavelength must be gt 0 */
00836                 ASSERT( C12O16Rotate[i].WLAng > 0 );
00837                 ASSERT( C13O16Rotate[i].WLAng > 0 );
00838 
00839                 C12O16Rotate[i].Emis->dampXvel = 
00840                         (realnum)(C12O16Rotate[i].Emis->Aul/
00841                   C12O16Rotate[i].EnergyWN/PI4);
00842 
00843                 C13O16Rotate[i].Emis->dampXvel = 
00844                         (realnum)(C13O16Rotate[i].Emis->Aul/
00845                   C13O16Rotate[i].EnergyWN/PI4);
00846 
00847                 /* derive the abs coefficient, call to function is gf, wl (A), g_low */
00848                 C12O16Rotate[i].Emis->opacity = 
00849                         (realnum)(abscf(C12O16Rotate[i].Emis->gf,
00850                   C12O16Rotate[i].EnergyWN,C12O16Rotate[i].Lo->g));
00851                 C13O16Rotate[i].Emis->opacity = 
00852                         (realnum)(abscf(C13O16Rotate[i].Emis->gf,
00853                   C13O16Rotate[i].EnergyWN,C13O16Rotate[i].Lo->g));
00854         }
00855 
00856         /* inner shell transitions */
00857         for( i=0; i<nUTA; ++i )
00858         {
00859                 if( UTALines[i].Emis->Aul > 0. )
00860                 {
00861                         UTALines[i].Emis->dampXvel = 
00862                                 (realnum)(UTALines[i].Emis->Aul/ UTALines[i].EnergyWN/PI4);
00863 
00864                         /* derive the abs coefficient, call to function is gf, wl (A), g_low */
00865                         UTALines[i].Emis->opacity = 
00866                                 (realnum)(abscf( UTALines[i].Emis->gf, UTALines[i].EnergyWN, UTALines[i].Lo->g));
00867 
00868                         /* excitation energy of transition in degrees kelvin */
00869                         UTALines[i].EnergyK = 
00870                                 (realnum)(T1CM*UTALines[i].EnergyWN);
00871 
00872                         /* energy of photon in ergs */
00873                         UTALines[i].EnergyErg = 
00874                                 (realnum)(ERG1CM*UTALines[i].EnergyWN);
00875 
00876                         /* put null terminated line label into chLab using line array*/
00877                         chIonLbl(chLab, &UTALines[i]);
00878 
00879                         /* get pointer to energy in continuum mesh */
00880                         UTALines[i].ipCont = ipLineEnergy(UTALines[i].EnergyWN * WAVNRYD , chLab,0 );
00881                         UTALines[i].Emis->ipFine = ipFineCont(UTALines[i].EnergyWN * WAVNRYD  );
00882 
00883                         /* find heating per absorption,
00884                          * first find threshold for this shell in ergs */
00885                         /* ionization threshold in erg */
00886                         double thresh = Heavy.Valence_IP_Ryd[UTALines[i].Hi->nelem-1][UTALines[i].Hi->IonStg-1] *EN1RYD;
00887                         UTALines[i].Coll.heat = (UTALines[i].EnergyErg-thresh);
00888                         ASSERT( UTALines[i].Coll.heat> 0. );
00889                 }
00890         }
00891 
00892         /* level 2 heavy element lines */
00893         for( i=0; i < nWindLine; i++ )
00894         {
00895 
00896                 /* derive the A, call to function is gf, wl (A), g_up */
00897                 TauLine2[i].Emis->Aul = 
00898                         (realnum)(eina(TauLine2[i].Emis->gf,
00899                   TauLine2[i].EnergyWN,TauLine2[i].Hi->g));
00900 
00901                 /* coefficient needed for damping constant - units cm s-1 */
00902                 TauLine2[i].Emis->dampXvel = 
00903                         (realnum)(TauLine2[i].Emis->Aul/
00904                   TauLine2[i].EnergyWN/PI4);
00905 
00906                 /* derive the abs coefficient, call to function is gf, wl (A), g_low */
00907                 TauLine2[i].Emis->opacity = 
00908                         (realnum)(abscf(TauLine2[i].Emis->gf,
00909                   TauLine2[i].EnergyWN,TauLine2[i].Lo->g));
00910 
00911                 /* excitation energy of transition in degrees kelvin */
00912                 TauLine2[i].EnergyK = 
00913                         (realnum)(T1CM*TauLine2[i].EnergyWN);
00914 
00915                 /* energy of photon in ergs */
00916                 TauLine2[i].EnergyErg = 
00917                         (realnum)(ERG1CM*TauLine2[i].EnergyWN);
00918 
00919                 /* put null terminated line label into chLab using line array*/
00920                 chIonLbl(chLab, &TauLine2[i]);
00921 
00922                 /* get pointer to energy in continuum mesh */
00923                 TauLine2[i].ipCont = ipLineEnergy(TauLine2[i].EnergyWN * WAVNRYD , chLab,0 );
00924                 TauLine2[i].Emis->ipFine = ipFineCont(TauLine2[i].EnergyWN * WAVNRYD );
00925                 /*if( TauLine2[i].ipCont==751 )
00926                         fprintf(ioQQQ,"( atom_level2 %s\n", chLab);*/
00927         }
00928 
00929         /* hyperfine structure lines */
00930         for( i=0; i < nHFLines; i++ )
00931         {
00932                 /* derive the A, call to function is gf, wl (A), g_up 
00933                 HFLines[i].Emis->Aul = 
00934                         (realnum)(eina(HFLines[i].Emis->gf,
00935                   HFLines[i].EnergyWN,HFLines[i].Hi->g));*/
00936 
00937                 /* coefficient needed for damping constant */
00938                 HFLines[i].Emis->dampXvel = 
00939                         (realnum)(HFLines[i].Emis->Aul/
00940                         HFLines[i].EnergyWN/PI4);
00941                 HFLines[i].Emis->damp = 1e-20f;
00942 
00943                 /* derive the abs coefficient, call to function is gf, wl (A), g_low */
00944                 HFLines[i].Emis->opacity = 
00945                         (realnum)(abscf(HFLines[i].Emis->gf,
00946                         HFLines[i].EnergyWN,
00947                         HFLines[i].Lo->g));
00948                 /* gf from this and 21 cm do not agree, A for HFS is 10x larger than level1 dat */
00949                 /*fprintf(ioQQQ,"HFLinesss\t%li\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
00950                         i,HFLines[i].Emis->opacity , HFLines[i].Emis->gf , HFLines[i].Emis->Aul , HFLines[i].EnergyWN,HFLines[i].Lo->g);*/
00951 
00952                 /* excitation energy of transition in degrees kelvin */
00953                 HFLines[i].EnergyK = 
00954                         (realnum)(T1CM*HFLines[i].EnergyWN);
00955 
00956                 /* energy of photon in ergs */
00957                 HFLines[i].EnergyErg = 
00958                         (realnum)(ERG1CM*HFLines[i].EnergyWN);
00959 
00960                 /* put null terminated line label into chLab using line array*/
00961                 chIonLbl(chLab, &HFLines[i]);
00962 
00963                 /* get pointer to energy in continuum mesh */
00964                 HFLines[i].ipCont = ipLineEnergy(HFLines[i].EnergyWN * WAVNRYD , chLab,0 );
00965                 HFLines[i].Emis->ipFine = ipFineCont(HFLines[i].EnergyWN * WAVNRYD );
00966         }
00967 
00968         /* Verner's FeII lines - do first so following labels will over write this
00969          * only call if large FeII atom is turned on */
00970         FeIIPoint();
00971 
00972         /* the group of inner shell fluorescent lines */
00973         for( i=0; i < t_yield::Inst().nlines(); ++i )
00974         {
00975                 strcpy( chLab , elementnames.chElementSym[t_yield::Inst().nelem(i)] );
00976                 strcat( chLab , elementnames.chIonStage[t_yield::Inst().ion_emit(i)] );
00977 
00978                 t_yield::Inst().set_ipoint( i, ipLineEnergy( t_yield::Inst().energy(i) , chLab , 0 ) );
00979         }
00980 
00981         /* ================================================================================== */
00982         /*        two photon two-photon 2-nu 2 nu 2 photon 2-photon                           */
00983 
00984         /* for induced two photon emission we need mirror image set
00985          * of continuum indices for continuum between Lya and half Lya */
00986         /* first MALLOC LIMELM dimension of space */
00987         /* >>chng 02 jun 28, malloc this NISO instead of 2.     */
00988         iso.ipSym2nu.reserve( NISO );
00989 
00990         /* now loop over the two iso-sequences with two photon continua */
00991         for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00992         {
00993                 iso.ipSym2nu.reserve( ipISO, LIMELM );
00994 
00995                 /* set up two photon emission */
00996                 for( nelem=ipISO; nelem<LIMELM; ++nelem )
00997                 {
00998                         if( dense.lgElmtOn[nelem] )
00999                         {
01000                                 double E2nu = Transitions[ipISO][nelem][TwoS][0].EnergyWN * WAVNRYD;
01001 
01002                                 /* pointer to the Lya energy */
01003                                 iso.ipTwoPhoE[ipISO][nelem] = ipoint(E2nu);
01004                                 /* this half-energy is only used to get induced rates in two photon */
01005                                 iso.ipHalfTwoPhoE[ipISO][nelem] = ipoint(E2nu / 2.);
01006                                 while( rfield.anu[iso.ipTwoPhoE[ipISO][nelem]] > E2nu )
01007                                 {
01008                                         --iso.ipTwoPhoE[ipISO][nelem];
01009                                 }
01010                                 while( rfield.anu[iso.ipHalfTwoPhoE[ipISO][nelem]] > E2nu / 2. )
01011                                 {
01012                                         --iso.ipHalfTwoPhoE[ipISO][nelem];
01013                                 }
01014 
01015                                 iso.ipSym2nu.reserve( ipISO, nelem, iso.ipTwoPhoE[ipISO][nelem] );
01016                         }
01017                 }
01018         }
01019 
01020         iso.ipSym2nu.alloc();
01021         iso.As2nu.alloc( iso.ipSym2nu.clone() );
01022 
01023         /* now loop over the two iso-sequences with two photon continua */
01024         for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
01025         {
01026                 /* set up two photon emission */
01027                 for( nelem=ipISO; nelem<LIMELM; ++nelem )
01028                 {
01029                         if( dense.lgElmtOn[nelem] )
01030                         {
01031                                 double E2nu = Transitions[ipISO][nelem][TwoS][0].EnergyWN * WAVNRYD;
01032                                 double Aul = Transitions[ipISO][nelem][TwoS][0].Emis->Aul;
01033                                 double SumShapeFunction = 0., Renorm= 0.;
01034 
01035                                 /* >>chng 02 aug 14, change upper limit to full Lya energy */
01036                                 for( i=0; i < iso.ipTwoPhoE[ipISO][nelem]; i++ )
01037                                 {
01038                                         /* energy is symmetric energy, the other side of half E2nu */
01039                                         energy = E2nu - rfield.anu[i];
01040                                         /* this is needed since mirror image of cell next to two-nu energy
01041                                          * may be slightly negative */
01042                                         energy = MAX2( energy, rfield.anu[0] + rfield.widflx[0]/2. );
01043                                         /* find index for this symmetric energy */
01044                                         iso.ipSym2nu[ipISO][nelem][i] = ipoint(energy);
01045                                         while( rfield.anu[iso.ipSym2nu[ipISO][nelem][i]] > E2nu ||
01046                                                 iso.ipSym2nu[ipISO][nelem][i] >= iso.ipTwoPhoE[ipISO][nelem])
01047                                         {
01048                                                 --iso.ipSym2nu[ipISO][nelem][i];
01049                                         }
01050                                         ASSERT( iso.ipSym2nu[ipISO][nelem][i] >= 0 );
01051                                 }
01052 
01053                                 /* ipTwoPhoE is the cell holding the 2nu energy itself, and we do not want
01054                                  * to include that in the following sum */
01055                                 for( i=0; i<iso.ipTwoPhoE[ipISO][nelem]; i++ )
01056                                 {
01057                                         double ShapeFunction;
01058 
01059                                         ASSERT( rfield.anu[i]<=E2nu );
01060 
01061                                         ShapeFunction = atmdat_2phot_shapefunction( rfield.AnuOrg[i]/E2nu, ipISO, nelem )*rfield.widflx[i]/E2nu;
01062                                         SumShapeFunction += ShapeFunction;
01063 
01064                                         /* >>refer      HI      2nu     Spitzer, L., & Greenstein, J., 1951, ApJ, 114, 407 */
01065                                         /* As2nu will add up to the A, so its units are s-1     */ 
01066                                         iso.As2nu[ipISO][nelem][i] = (realnum)( Aul * ShapeFunction );
01067                                 }
01068 
01069                                 /* The spline function in twophoton.c causes a bit of an error in the integral of the
01070                                  * shape function.  So we renormalize the integral to 1.        */
01071                                 Renorm = 1./SumShapeFunction;
01072 
01073                                 for( i=0; i<iso.ipTwoPhoE[ipISO][nelem]; i++ )
01074                                 {
01075                                         iso.As2nu[ipISO][nelem][i] *= (realnum)Renorm;
01076                                 }
01077 
01078                                 /* The result should be VERY close to 1.        */
01079                                 ASSERT( fabs( SumShapeFunction*Renorm - 1. ) < 0.00001 );
01080                         }
01081                 }
01082         }
01083 
01084         {
01085                 /* this is an option to print out one of the two photon continua */
01086                 enum {DEBUG_LOC=false};
01087                 if( DEBUG_LOC )
01088                 {       
01089 #                       define  NCRS    21
01090                         double limit,ener[NCRS]={
01091                           0.,     0.03738,  0.07506,  0.1124,  0.1498,  0.1875,
01092                           0.225,  0.263,    0.300,    0.3373,  0.375,   0.4127,
01093                           0.4500, 0.487,    0.525,    0.5625,  0.6002,  0.6376,
01094                           0.6749, 0.7126,   0.75};
01095 
01096                         nelem = ipHYDROGEN;
01097                         ipISO = ipHYDROGEN;
01098 
01099                         limit = iso.ipTwoPhoE[ipISO][nelem];
01100 
01101                         for( i=0; i < NCRS; i++ )
01102                         {
01103                                 fprintf(ioQQQ,"%.3e\t%.3e\n", ener[i] , 
01104                                         atmdat_2phot_shapefunction( ener[i]/0.75, ipISO, nelem ) );
01105                         }
01106 
01107                         xnew = 0.;
01109                         for( i=0; i < limit; i++ )
01110                         {
01111                                 double fach = iso.As2nu[ipISO][nelem][i]/2.*rfield.anu2[i]/rfield.widflx[i]*EN1RYD;
01112                                 fprintf(ioQQQ,"%.3e\t%.3e\t%.3e\n", 
01113                                         rfield.anu[i] , 
01114                                         iso.As2nu[ipISO][nelem][i] / rfield.widflx[i] , 
01115                                         fach );
01116                                 xnew += iso.As2nu[ipISO][nelem][i];
01117                         }
01118                         fprintf(ioQQQ," sum is %.3e\n", xnew );
01119                         cdEXIT(EXIT_FAILURE);
01120                 }
01121         }
01122 
01123         {
01124                 /* this is an option to print out one of the two photon continua */
01125                 enum {DEBUG_LOC=false};
01126                 if( DEBUG_LOC )
01127                 {       
01128                         for( i=0; i<11; ++i )
01129                         {
01130                                 char chLsav[10];
01131                                 TauDummy.WLAng = (realnum)(PI * pow(10.,(double)i));
01132                                 strcpy( chLsav, chLineLbl(&TauDummy) );
01133                                 fprintf(ioQQQ,"%.2f\t%s\n", TauDummy.WLAng , chLsav );
01134                         }
01135                         cdEXIT(EXIT_FAILURE);
01136                 }
01137         }
01138 
01139         /* option to print out whole thing with "trace lines" command */
01140         if( trace.lgTrLine )
01141         {
01142                 fprintf( ioQQQ, "       WL(Ang)   E(RYD)   IP   gl  gu      gf       A        damp     abs K\n" );
01143                 for( i=1; i <= nLevel1; i++ )
01144                 {
01145                         strcpy( chLab, chLineLbl(&TauLines[i]) );
01146                         iWL_Ang = (long)TauLines[i].WLAng;
01147                         if( iWL_Ang > 1000000 )
01148                         {
01149                                 iWL_Ang /= 10000;
01150                         }
01151                         else if( iWL_Ang > 10000 )
01152                         {
01153                                 iWL_Ang /= 1000;
01154                         }
01155 
01156                         fprintf( ioQQQ, " %10.10s%5ld%10.3e %4li%4ld%4ld%10.2e%10.2e%10.2e%10.2e\n", 
01157                           chLab, iWL_Ang, RYDLAM/TauLines[i].WLAng, 
01158                           TauLines[i].ipCont, (long)(TauLines[i].Lo->g), 
01159                           (long)(TauLines[i].Hi->g), TauLines[i].Emis->gf, 
01160                           TauLines[i].Emis->Aul, TauLines[i].Emis->dampXvel, 
01161                           TauLines[i].Emis->opacity );
01162                 }
01163 
01164                 /*Atomic Or Molecular lines*/
01165                 for( i=0; i <linesAdded2; i++)
01166                 {
01167                         strcpy( chLab, chLineLbl(atmolEmis[i].tran) );
01168 
01169                         iWL_Ang = (long)atmolEmis[i].tran->WLAng;
01170 
01171                         if( iWL_Ang > 1000000 )
01172                         {
01173                                 iWL_Ang /= 10000;
01174                         }
01175                         else if( iWL_Ang > 10000 )
01176                         {
01177                                 iWL_Ang /= 1000;
01178                         }
01179                         fprintf( ioQQQ, " %10.10s%5ld%10.3e %4li%4ld%4ld%10.2e%10.2e%10.2e%10.2e\n", 
01180                           chLab, iWL_Ang, RYDLAM/atmolEmis[i].tran->WLAng, 
01181                           atmolEmis[i].tran->ipCont, (long)(atmolEmis[i].tran->Lo->g), 
01182                           (long)(atmolEmis[i].tran->Hi->g),atmolEmis[i].gf, 
01183                           atmolEmis[i].Aul,atmolEmis[i].dampXvel, 
01184                           atmolEmis[i].opacity);
01185                 }
01186 
01187                 /* C12O16 lines */
01188                 for( i=0; i < nCORotate; i++ )
01189                 {
01190                         strcpy( chLab, chLineLbl(&C12O16Rotate[i]) );
01191 
01192                         iWL_Ang = (long)C12O16Rotate[i].WLAng;
01193 
01194                         if( iWL_Ang > 1000000 )
01195                         {
01196                                 iWL_Ang /= 10000;
01197                         }
01198                         else if( iWL_Ang > 10000 )
01199                         {
01200                                 iWL_Ang /= 1000;
01201                         }
01202                         fprintf( ioQQQ, " %10.10s%5ld%10.3e %4li%4ld%4ld%10.2e%10.2e%10.2e%10.2e\n", 
01203                           chLab, iWL_Ang, RYDLAM/C12O16Rotate[i].WLAng, 
01204                           C12O16Rotate[i].ipCont, (long)(C12O16Rotate[i].Lo->g), 
01205                           (long)(C12O16Rotate[i].Hi->g), C12O16Rotate[i].Emis->gf, 
01206                           C12O16Rotate[i].Emis->Aul, C12O16Rotate[i].Emis->dampXvel, 
01207                           C12O16Rotate[i].Emis->opacity );
01208                 }
01209 
01210                 /* C13O16 lines */
01211                 for( i=0; i < nCORotate; i++ )
01212                 {
01213                         strcpy( chLab, chLineLbl(&C13O16Rotate[i]) );
01214 
01215                         iWL_Ang = (long)C13O16Rotate[i].WLAng;
01216 
01217                         if( iWL_Ang > 1000000 )
01218                         {
01219                                 iWL_Ang /= 10000;
01220                         }
01221                         else if( iWL_Ang > 10000 )
01222                         {
01223                                 iWL_Ang /= 1000;
01224                         }
01225                         fprintf( ioQQQ, " %10.10s%5ld%10.3e %4li%4ld%4ld%10.2e%10.2e%10.2e%10.2e\n", 
01226                           chLab, iWL_Ang, RYDLAM/C13O16Rotate[i].WLAng, 
01227                           C13O16Rotate[i].ipCont, (long)(C13O16Rotate[i].Lo->g), 
01228                           (long)(C13O16Rotate[i].Hi->g), C13O16Rotate[i].Emis->gf, 
01229                           C13O16Rotate[i].Emis->Aul, C13O16Rotate[i].Emis->dampXvel, 
01230                           C13O16Rotate[i].Emis->opacity );
01231                 }
01232 
01233                 for( i=0; i < nWindLine; i++ )
01234                 {
01235                         strcpy( chLab, chLineLbl(&TauLine2[i]) );
01236 
01237                         iWL_Ang = (long)TauLine2[i].WLAng;
01238 
01239                         if( iWL_Ang > 1000000 )
01240                         {
01241                                 iWL_Ang /= 10000;
01242                         }
01243                         else if( iWL_Ang > 10000 )
01244                         {
01245                                 iWL_Ang /= 1000;
01246                         }
01247                         fprintf( ioQQQ, " %10.10s%5ld%10.3e %4li%4ld%4ld%10.2e%10.2e%10.2e%10.2e\n", 
01248                           chLab, iWL_Ang, RYDLAM/TauLine2[i].WLAng, 
01249                           TauLine2[i].ipCont, (long)(TauLine2[i].Lo->g), 
01250                           (long)(TauLine2[i].Hi->g), TauLine2[i].Emis->gf, 
01251                           TauLine2[i].Emis->Aul, TauLine2[i].Emis->dampXvel, 
01252                           TauLine2[i].Emis->opacity );
01253                 }
01254                 for( i=0; i < nHFLines; i++ )
01255                 {
01256                         strcpy( chLab, chLineLbl(&HFLines[i]) );
01257 
01258                         iWL_Ang = (long)HFLines[i].WLAng;
01259 
01260                         if( iWL_Ang > 1000000 )
01261                         {
01262                                 iWL_Ang /= 10000;
01263                         }
01264                         else if( iWL_Ang > 10000 )
01265                         {
01266                                 iWL_Ang /= 1000;
01267                         }
01268                         fprintf( ioQQQ, " %10.10s%5ld%10.3e %4li%4ld%4ld%10.2e%10.2e%10.2e%10.2e\n", 
01269                           chLab, iWL_Ang, RYDLAM/HFLines[i].WLAng, 
01270                           HFLines[i].ipCont, (long)(HFLines[i].Lo->g), 
01271                           (long)(HFLines[i].Hi->g), HFLines[i].Emis->gf, 
01272                           HFLines[i].Emis->Aul, HFLines[i].Emis->dampXvel, 
01273                           HFLines[i].Emis->opacity );
01274                 }
01275         }
01276 
01277         /* this is an option to kill fine structure line optical depths */
01278         if( !rt.lgFstOn )
01279         {
01280                 for( i=1; i <= nLevel1; i++ )
01281                 {
01282                         if( TauLines[i].EnergyWN < 10000. )
01283                         {
01284                                 TauLines[i].Emis->opacity = 0.;
01285                         }
01286                 }
01287 
01288                 /*Atomic or Molecular Lines-Humeshkar Nemala*/
01289                 for( i=0; i <linesAdded2; i++)
01290                 {
01291                         if(atmolEmis[i].tran->EnergyWN < 10000. )
01292                         {
01293                                 atmolEmis[i].opacity = 0.;
01294                         }
01295                 }
01296         }
01297 
01298         /* read in continuum bands data set */
01299         ContBandsCreate( "" );
01300 
01301         /* we're done adding lines and states to the stacks.  
01302          * This flag is used to make sure we never add them again in this coreload. */
01303         lgLinesAdded = true;
01304         lgStatesAdded = true;
01305 
01306         return;
01307 }
01308 
01309 /*fiddle adjust energy bounds of cell with index ipnt so that lower energy
01310  * matches ionization edges exactly, called by createpoint */
01311 /* ipnt is index on c scale */
01312 STATIC void fiddle(long int ipnt, 
01313   double exact)
01314 {
01315         realnum Ehi, 
01316           Elo,
01317           OldEner;
01318 
01319         DEBUG_ENTRY( "fiddle()" );
01320 
01321         /* make low edge of cell exact for photo integrals */
01322         ASSERT( ipnt >= 0 );
01323         ASSERT( ipnt < rfield.nupper-1 );
01324 
01325         /* upper edge of higher cell*/
01326         Ehi = rfield.anu[ipnt] + rfield.widflx[ipnt]*0.5f;
01327         /* lower edge of lower cell */
01328         Elo = rfield.anu[ipnt-1] - rfield.widflx[ipnt-1]*0.5f;
01329 
01330         /* >>chng 02 nov 11, do nothing if already very close,
01331          * because VERY high res models had negative widths for some very close edges */
01332         if( fabs( Elo/exact - 1. ) < 0.001 ) 
01333                 return;
01334 
01335         ASSERT( Elo <= exact );
01336 
01337         OldEner = rfield.anu[ipnt];
01338 
01339         /* centroid of desired lower energy and upper edge */
01340         rfield.anu[ipnt] = (realnum)((Ehi + exact)/2.);
01341         /* same for lower */
01342         rfield.anu[ipnt-1] = (realnum)((exact + Elo)/2.);
01343 
01344         rfield.widflx[ipnt] = (realnum)(Ehi - exact);
01345         rfield.widflx[ipnt-1] = (realnum)(exact - Elo);
01346 
01347         /* bring upper cell down a bit too, since we dont want large change in widflx */
01348         rfield.anu[ipnt+1] -= (OldEner - rfield.anu[ipnt])/2.f;
01349 
01350         /* sanity check */
01351         ASSERT( rfield.widflx[ipnt-1] > 0. );
01352         ASSERT( rfield.widflx[ipnt] > 0. );
01353         return;
01354 }
01355 
01356 /*ipShells assign continuum energy pointers to shells for all atoms,
01357  * called by ContCreatePointers */
01358 STATIC void ipShells(
01359         /* nelem is the atomic number on the C scale, Li is 2 */
01360         long int nelem)
01361 {
01362         char chLab[5];
01363         long int 
01364           imax, 
01365           ion, 
01366           nelec, 
01367           ns, 
01368           nshell;
01369         /* following value cannot be used - will be set to proper threshold */
01370         double thresh=-DBL_MAX;
01371 
01372         DEBUG_ENTRY( "ipShells()" );
01373 
01374         ASSERT( nelem >= NISO);
01375         ASSERT( nelem < LIMELM );
01376 
01377         /* fills in pointers to valence shell ionization threshold
01378          * PH1(a,b,c,d)
01379          * a=1 => thresh, others fitting parameters
01380          * b atomic number
01381          * c number of electrons
01382          * d shell number 7-1 */
01383 
01384         /* threshold in Ryd
01385          * ion=0 for atom, up to nelem-1 for helium like, hydrogenic is elsewhere */
01386         for( ion=0; ion < nelem; ion++ )
01387         {
01388                 /* generate label for ionization stage */
01389                 /* this is short form of element name */
01390                 strcpy( chLab, elementnames.chElementSym[nelem] );
01391 
01392                 /* this is a number between 1 and 31 */
01393                 strcat( chLab, elementnames.chIonStage[ion] );
01394 
01395                 /* this is the iso sequence - must not redo sequence if done as iso */
01396                 long int ipISO = nelem-ion;
01397 
01398                 /* number of bound electrons */
01399                 nelec = ipISO+1;
01400 
01401                 /* nsShells(nelem,ion) is the number of shells for ion with nelec electrons,
01402                  * physical not c scale */
01403                 imax = Heavy.nsShells[nelem][ion];
01404 
01405                 /* loop on all inner shells, valence shell */
01406                 for( nshell=0; nshell < imax; nshell++ )
01407                 {
01408                         /* ionization potential of this shell in rydbergs */
01409                         thresh = (double)(t_ADfA::Inst().ph1(nshell,nelec-1,nelem,0)/EVRYD* 0.9998787);
01410                         if( thresh <= 0.1 )
01411                         {
01412                                 /* negative ip shell does not exist, set upper limit
01413                                  * to less than lower limit so this never looped upon
01414                                  * these are used as flags by LimitSh to check whether
01415                                  * this is a real shell - if 1 or 2 is changed - change LimitSh!! */
01416                                 opac.ipElement[nelem][ion][nshell][0] = 2;
01417                                 opac.ipElement[nelem][ion][nshell][1] = 1;
01418                         }
01419                         else
01420                         {
01421                                 /* this is lower bound to energy range for this shell */
01422                                 /* >>chng 02 may 27, change to version of ip with label, so that
01423                                  * inner shell edges will appear */
01424                                 /*opac.ipElement[nelem][ion][nshell][0] = ipoint(thresh);*/
01425                                 opac.ipElement[nelem][ion][nshell][0] = 
01426                                         ipContEnergy( thresh , chLab );
01427 
01428                                 /* this is upper bound to energy range for this shell 
01429                                  * LimitSh is an integer function, returns pointer
01430                                  * to threshold of next major shell.  For k-shell it
01431                                  * returns the values KshellLimit, default=7.35e4
01432                                  * >>chng 96 sep 26, had been below, result zero cross sec at 
01433                                  * many energies where opacity project did not produce state specific 
01434                                  * cross section */
01435                                 opac.ipElement[nelem][ion][nshell][1] = 
01436                                         LimitSh(ion+1,  nshell+1,nelem+1);
01437                                 ASSERT( opac.ipElement[nelem][ion][nshell][1] > 0);
01438                         }
01439                 }
01440 
01441                 ASSERT( imax > 0 && imax <= 7 );
01442 
01443                 /* this will be index pointing to valence edge */
01444                 /* [0] is pointer to threshold in energy array */
01445                 opac.ipElement[nelem][ion][imax-1][0] = 
01446                         ipContEnergy(thresh, chLab);
01447 
01448                 /* pointer to valence electron ionization potential */
01449                 Heavy.ipHeavy[nelem][ion] = opac.ipElement[nelem][ion][imax-1][0];
01450 
01451                 /* ionization potential of valence shell in Ryd 
01452                  * thresh was evaluated above, now has last value, the valence shell */
01453                 Heavy.Valence_IP_Ryd[nelem][ion] = thresh;
01454 
01455                 Heavy.xLyaHeavy[nelem][ion] = 0.;
01456                 if( ipISO >= NISO )
01457                 {
01458                         /* this is set of 3/4 of valence shell IP, this is important
01459                          * source of ots deep in cloud */
01460                         Heavy.ipLyHeavy[nelem][ion] = 
01461                                 ipLineEnergy(thresh*0.75,chLab , 0);
01462 
01463 
01464                         Heavy.ipBalHeavy[nelem][ion] = 
01465                                 ipLineEnergy(thresh*0.25,chLab , 0);
01466                 }
01467                 else
01468                 {
01469                         /* do not treat this simple way since done exactly with iso 
01470                          * sequences */
01471                         Heavy.ipLyHeavy[nelem][ion] = -1;
01472                         Heavy.ipBalHeavy[nelem][ion] = -1;
01473                 }
01474         }
01475 
01476         /* above loop did up to hydrogenic, now do hydrogenic - 
01477          * hydrogenic is special since arrays already set up */
01478         Heavy.nsShells[nelem][nelem] = 1;
01479 
01480         /* this is lower limit to range */
01481         /* hydrogenic photoionization set to special hydro array 
01482          * this is pointer to threshold energy */
01483         /* this statement is in ContCreatePointers but has not been done when this routine called */
01484         /*iso.ipIsoLevNIonCon[ipH_LIKE][ipZ][ipLo] = ipContEnergy(iso.xIsoLevNIonRyd[ipH_LIKE][ipZ][ipLo],chLab);*/
01485         /*opac.ipElement[nelem][nelem][0][0] = iso.ipIsoLevNIonCon[ipH_LIKE][nelem][ipH1s];*/
01486         opac.ipElement[nelem][nelem][0][0] = ipoint( iso.xIsoLevNIonRyd[ipH_LIKE][nelem][ipH1s] );
01487         ASSERT( opac.ipElement[nelem][nelem][0][0] > 0 );
01488 
01489         /* this is the high-energy limit */
01490         opac.ipElement[nelem][nelem][0][1] = continuum.KshellLimit;
01491 
01492         Heavy.ipHeavy[nelem][nelem] = opac.ipElement[nelem][nelem][0][0];
01493 
01494         /* this is for backwards computability with Cambridge code */
01495         if( trace.lgTrace && trace.lgPointBug )
01496         {
01497                 for( ion=0; ion < (nelem+1); ion++ )
01498                 {
01499                         fprintf( ioQQQ, "Ion:%3ld%3ld %2.2s%2.2s total shells:%3ld\n", 
01500                           nelem, ion+1, elementnames.chElementSym[nelem], elementnames.chIonStage[ion]
01501                           , Heavy.nsShells[nelem][ion] );
01502                         for( ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
01503                         {
01504                                 fprintf( ioQQQ, " shell%3ld %2.2s range eV%10.2e-%8.2e\n", 
01505                                   ns+1, Heavy.chShell[ns], rfield.anu[opac.ipElement[nelem][ion][ns][0]-1]*
01506                                   EVRYD, rfield.anu[opac.ipElement[nelem][ion][ns][1]-1]*EVRYD );
01507                         }
01508                 }
01509         }
01510         return;
01511 }
01512 
01513 /*LimitSh sets upper energy limit to subshell integrations */
01514 STATIC long LimitSh(long int ion, 
01515   long int nshell, 
01516   long int nelem)
01517 {
01518         long int LimitSh_v;
01519 
01520         DEBUG_ENTRY( "LimitSh()" );
01521 
01522         /* this routine returns the high-energy limit to the energy range
01523          * for photoionization of a given shell
01524          * */
01525         if( nshell == 1 )
01526         {
01527                 /* this limit is high-energy limit to code unless changed with set kshell */
01528                 LimitSh_v = continuum.KshellLimit;
01529 
01530         }
01531         else if( nshell == 2 )
01532         {
01533                 /* this is 2s shell, upper limit is 1s
01534                  * >>chng 96 oct 08, up to high-energy limit
01535                  * LimitSh = ipElement(nelem,ion , 1,1)-1 */
01536                 LimitSh_v = continuum.KshellLimit;
01537 
01538         }
01539         else if( nshell == 3 )
01540         {
01541                 /* this is 2p shell, upper limit is 1s
01542                  * >>chng 96 oct 08, up to high-energy limit
01543                  * LimitSh = ipElement(nelem,ion , 1,1)-1 */
01544                 LimitSh_v = continuum.KshellLimit;
01545 
01546         }
01547         else if( nshell == 4 )
01548         {
01549                 /* this is 3s shell, upper limit is 2p
01550                  * >>chng 96 oct 08, up to K-shell edge
01551                  * LimitSh = ipElement(nelem,ion , 3,1)-1 */
01552                 LimitSh_v = opac.ipElement[nelem-1][ion-1][0][0] - 1;
01553 
01554         }
01555         else if( nshell == 5 )
01556         {
01557                 /* this is 3p shell, upper limit is 2p
01558                  * >>chng 96 oct 08, up to K-shell edge
01559                  * LimitSh = ipElement(nelem,ion , 3,1)-1 */
01560                 LimitSh_v = opac.ipElement[nelem-1][ion-1][0][0] - 1;
01561 
01562         }
01563         else if( nshell == 6 )
01564         {
01565                 /* this is 3d shell, upper limit is 2p
01566                  * >>chng 96 oct 08, up to K-shell edge
01567                  * LimitSh = ipElement(nelem,ion , 3,1)-1 */
01568                 LimitSh_v = opac.ipElement[nelem-1][ion-1][0][0] - 1;
01569 
01570         }
01571         else if( nshell == 7 )
01572         {
01573                 /* this is 4s shell, upper limit is 3d */
01574                 if( opac.ipElement[nelem-1][ion-1][5][0] < 3 )
01575                 {
01576                         /* this is check for empty shell 6, 3d
01577                          * if so then set to 3p instead */
01578                         LimitSh_v = opac.ipElement[nelem-1][ion-1][4][0] - 
01579                           1;
01580                 }
01581                 else
01582                 {
01583                         LimitSh_v = opac.ipElement[nelem-1][ion-1][5][0] - 
01584                           1;
01585                 }
01586                 /* >>chng 96 sep 26, set upper limit down to 2s */
01587                 LimitSh_v = opac.ipElement[nelem-1][ion-1][2][0] - 1;
01588 
01589         }
01590         else
01591         {
01592                 fprintf( ioQQQ, " LimitSh cannot handle nshell as large as%4ld\n", 
01593                   nshell );
01594                 cdEXIT(EXIT_FAILURE);
01595         }
01596         return( LimitSh_v );
01597 }
01598 
01599 /*ContBandsCreate - read set of continuum bands to enter total emission into line*/
01600 STATIC void ContBandsCreate(
01601         /* chFile is optional filename, if void then use default bands,
01602          * if not void then use file specified,
01603          * return value is 0 for success, 1 for failure */
01604          const char chFile[] )
01605 {
01606         char chLine[FILENAME_PATH_LENGTH_2];
01607         const char* chFilename;
01608         FILE *ioDATA;
01609 
01610         bool lgEOL;
01611         long int i,k;
01612 
01613         /* keep track of whether we have been called - want to be
01614          * called a total of one time */
01615         static bool lgCalled=false;
01616 
01617         DEBUG_ENTRY( "ContBandsCreate()" );
01618 
01619         /* do nothing if second or later call*/
01620         if( lgCalled )
01621         {
01622                 /* success */
01623                 return;
01624         }
01625         lgCalled = true;
01626 
01627         /* use default filename if void string, else use file specified */
01628         chFilename = ( strlen(chFile) == 0 ) ? "bands_continuum.dat" : chFile;
01629 
01630         /* get continuum band data  */
01631         if( trace.lgTrace )
01632         {
01633                 fprintf( ioQQQ, " ContBandsCreate opening %s:", chFilename );
01634         }
01635 
01636         ioDATA = open_data( chFilename, "r" );
01637 
01638         /* now count how many bands are in the file */
01639         continuum.nContBand = 0;
01640 
01641         /* first line is a magic number and does not count as a band*/
01642         if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
01643         {
01644                 fprintf( ioQQQ, " ContBandsCreate could not read first line of %s.\n", chFilename );
01645                 cdEXIT(EXIT_FAILURE);
01646         }
01647         while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
01648         {
01649                 /* we want to count the lines that do not start with #
01650                  * since these contain data */
01651                 if( chLine[0] != '#')
01652                         ++continuum.nContBand;
01653         }
01654 
01655         /* now rewind the file so we can read it a second time*/
01656         if( fseek( ioDATA , 0 , SEEK_SET ) != 0 )
01657         {
01658                 fprintf( ioQQQ, " ContBandsCreate could not rewind %s.\n", chFilename );
01659                 cdEXIT(EXIT_FAILURE);
01660         }
01661 
01662         continuum.ContBandWavelength = (realnum *)MALLOC(sizeof(realnum)*(unsigned)(continuum.nContBand) );
01663         continuum.chContBandLabels = (char **)MALLOC(sizeof(char *)*(unsigned)(continuum.nContBand) );
01664         continuum.ipContBandLow = (long int *)MALLOC(sizeof(long int)*(unsigned)(continuum.nContBand) );
01665         continuum.ipContBandHi = (long int *)MALLOC(sizeof(long int)*(unsigned)(continuum.nContBand) );
01666 
01667         /* now make second dim, id wavelength, and lower and upper bounds */
01668         for( i=0; i<continuum.nContBand; ++i )
01669         {
01670                 /* array of labels, each 4 long plus 0 at [4] */
01671                 continuum.chContBandLabels[i] = (char *)MALLOC(sizeof(char)*(unsigned)(5) );
01672         }
01673 
01674         /* first line is a versioning magic number - now confirm that it is valid */
01675         if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
01676         {
01677                 fprintf( ioQQQ, " ContBandsCreate could not read first line of %s.\n", chFilename );
01678                 cdEXIT(EXIT_FAILURE);
01679         }
01680         /* bands_continuum magic number here <- this string is in band_continuum.dat
01681          * with comment to search for this to find magic number 
01682          * >>chng 06 aug 07, update bands and magic number */
01683         {
01684                 long int m1 , m2 , m3,
01685                         myr = 6,
01686                         mmo = 8,
01687                         mdy = 7;
01688                 i = 1;
01689                 m1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
01690                 m2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
01691                 m3 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
01692                 if( ( m1 != myr ) ||
01693                     ( m2 != mmo ) ||
01694                     ( m3 != mdy ) )
01695                 {
01696                         fprintf( ioQQQ, 
01697                                 " ContBandsCreate: the version of the data file %s I found (%li %li %li)is not the current version (%li %li %li).\n", 
01698                                 chFilename ,
01699                                 m1 , m2 , m3 ,
01700                                 myr , mmo , mdy );
01701                         fprintf( ioQQQ, 
01702                                 " ContBandsCreate: you need to update this file.\n");
01703                         cdEXIT(EXIT_FAILURE);
01704                 }
01705         }
01706 
01707         /* now read in data again, but save it this time */
01708         k = 0;
01709         while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
01710         {
01711                 /* we want to count the lines that do not start with #
01712                  * since these contain data */
01713                 if( chLine[0] != '#')
01714                 {
01715                         double xLow , xHi;
01716                         /* copy 4 char label plus termination */
01717                         strncpy( continuum.chContBandLabels[k] , chLine , 4 );
01718                         continuum.chContBandLabels[k][4] = 0;
01719 
01720                         /* now get central band wavelength 
01721                          * >>chng 06 aug 11 from 4 to 6, the first 4 char are labels and
01722                          * these can contain numbers, next comes a space, then the number */
01723                         i = 6;
01724                         continuum.ContBandWavelength[k] = (realnum)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
01725                         /* >>chng 06 feb 21, multiply by 1e4 to convert micron wavelength into Angstroms,
01726                          * which is assumed by the code.  before this correction the band centroid 
01727                          * wavelength was given in the output incorrectly listed as Angstroms.
01728                          * results were correct just label was wrong */
01729                         continuum.ContBandWavelength[k] *= 1e4f;
01730 
01731                         /* these are short and long wave limits, which are high and
01732                          * low energy limits - these are now wl in microns */
01733                         xHi = (realnum)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
01734                         xLow = (realnum)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
01735                         if( lgEOL )
01736                         {
01737                                 fprintf( ioQQQ, " There should have been 3 numbers on this band line.   Sorry.\n" );
01738                                 fprintf( ioQQQ, " string==%s==\n" ,chLine );
01739                                 cdEXIT(EXIT_FAILURE);
01740                         }
01741 
01742                         /* make sure bands bounds are in correct order, shorter - longer wavelength*/
01743                         if( xHi >= xLow )
01744                         {
01745                                 fprintf( ioQQQ, " ContBandWavelength band %li has improper bounds.\n" ,i);
01746                                 cdEXIT(EXIT_FAILURE);
01747                         }
01748 
01749                         /* get continuum index - RYDLAM is 911.6A = 1 Ryd so 1e4 converts 
01750                          * micron to Angstrom */
01751                         continuum.ipContBandHi[k] = ipoint( RYDLAM / (xHi*1e4) );
01752                         continuum.ipContBandLow[k] = ipoint( RYDLAM / (xLow*1e4) );
01753                         /*fprintf(ioQQQ,"DEBUG bands_continuum %s %.3e %li %li \n",
01754                                 continuum.chContBandLabels[k],
01755                                 continuum.ContBandWavelength[k],
01756                                 continuum.ipContBandHi[k] , 
01757                                 continuum.ipContBandLow[k]);*/
01758 
01759                         if( trace.lgTrace && trace.lgConBug )
01760                         {
01761                                 if( k==0 )
01762                                         fprintf( ioQQQ, "   ContCreatePointer trace bands\n");
01763                                 fprintf( ioQQQ, 
01764                                         "     band %ld label %s low wl= %.3e low ipnt= %li "
01765                                         " hi wl= %.3e hi ipnt= %li \n", 
01766                                         k, 
01767                                         continuum.chContBandLabels[k] ,
01768                                         xLow,
01769                                         continuum.ipContBandLow[k] ,
01770                                         xHi,
01771                                         continuum.ipContBandHi[k]  );
01772                         }
01773                         /*fprintf(ioQQQ,
01774                                 "DEBUG band data %s %f %f %f %f %f \n", 
01775                                 continuum.chContBandLabels[k],
01776                                 continuum.ContBandWavelength[k],
01777                                 xHi,
01778                                 rfield.anu[continuum.ipContBandHi[k]-1],
01779                                 xLow,
01780                                 rfield.anu[continuum.ipContBandLow[k]-1]);*/
01781                         ++k;
01782                 }
01783         }
01784         /* now validate this incoming data */
01785         for( i=0; i<continuum.nContBand; ++i )
01786         {
01787                 /* make sure all are positive */
01788                 if( continuum.ContBandWavelength[i] <=0. )
01789                 {
01790                         fprintf( ioQQQ, " ContBandWavelength band %li has non-positive entry.\n",i );
01791                         cdEXIT(EXIT_FAILURE);
01792                 }
01793         }
01794 
01795         fclose(ioDATA);
01796         return;
01797 }
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

Generated for cloudy by doxygen 1.7.3