cloudy  trunk
lines_service.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 /*PutLine enter local line intensity into the intensity stack for eventual printout */
00004 /*PutExtra enter and 'extra' intensity source for some line */
00005 /*linadd enter lines into the line storage array, called once per zone */
00006 /*DumpLine print various information about an emission line vector, 
00007  * used in debugging, print to std out, ioQQQ */
00008 /*TexcLine derive excitation temperature of line from contents of line array */
00009 /*GetGF convert Einstein A into oscillator strength */
00010 /*emit_frac returns fraction of populations the produce emission */
00011 /*abscf convert gf into absorption coefficient */
00012 /*chIonLbl use information in line array to generate a null terminated ion label in "Fe 2" */
00013 /*chLineLbl use information in line transfer arrays to generate a line label */
00014 /*RefIndex calculates the index of refraction of air using the line energy in wavenumbers,
00015  * used to convert vacuum wavelengths to air wavelengths. */
00016 /*eina convert a gf into an Einstein A */
00017 /*WavlenErrorGet - find difference between two wavelengths */
00018 /*PutCS enter a collision strength into an individual line vector */
00019 /*lindst add local line intensity to line luminosity stack */
00020 /*PntForLine generate pointer for forbidden line */
00021 /*TransitionZero zeros out TransitionZero */
00022 /*EmLineJunk set all elements of transition struc to dangerous values */
00023 /*EmLineZero set all elements of transition struc to zero */
00024 /*LineConvRate2CS convert down coll rate back into electron cs in case other parts of code need this for reference */
00025 /*lgTauGood returns true is we have not overrun optical depth scale */
00026 /*OccupationNumberLine - derive the photon occupation number at line center for any line */
00027 /*outline - adds line photons to reflin and outlin */
00028 /*MakeCS compute collision strength by g-bar approximations */
00029 /*gbar1 compute g-bar collision strength using Mewe approximations */
00030 /*gbar0 compute g-bar gaunt factor for neutrals */
00031 /*totlin sum total intensity of cooling, recombination, or intensity lines */
00032 /*FndLineHt search through line heat arrays to find the strongest heat source */
00033 /*ConvRate2CS convert down coll rate back into electron cs in case other parts of code need this for reference */
00034 #include "cddefines.h"
00035 #include "physconst.h"
00036 #include "phycon.h"
00037 #include "lines.h"
00038 #include "radius.h"
00039 #include "geometry.h"
00040 #include "elementnames.h"
00041 #include "rt.h"
00042 #include "dense.h"
00043 #include "rfield.h"
00044 #include "opacity.h"
00045 #include "ipoint.h"
00046 #include "iso.h"
00047 #include "taulines.h"
00048 #include "hydrogenic.h"
00049 #include "lines_service.h"
00050 
00051 /*outline - adds line photons to reflin and outlin */
00052 void outline( transition *t )
00053 {
00054         long int ip = t->ipCont-1;
00055         double xInWrd = t->Emis->phots*t->Emis->FracInwd;
00056 
00057         DEBUG_ENTRY( "outline()" );
00058 
00059         ASSERT( t->Emis->phots >= 0. );
00060         ASSERT( t->Emis->FracInwd >= 0. );
00061         ASSERT( radius.BeamInIn >= 0. );
00062         ASSERT( radius.BeamInOut >= 0. );
00063 
00064         /* the reflected part */
00065         rfield.reflin[ip] += (realnum)(xInWrd*radius.BeamInIn);
00066 
00067         /* inward beam that goes out since sphere set */
00068         rfield.outlin[ip] += (realnum)(xInWrd*radius.BeamInOut*opac.tmn[ip]*
00069                 t->Emis->ColOvTot);
00070 
00071         /* outward part */
00073         rfield.outlin[ip] += (realnum)(t->Emis->phots*
00074                 (1. - t->Emis->FracInwd)*radius.BeamOutOut* t->Emis->ColOvTot);
00075         return;
00076 }
00077 
00078 /*emit_frac returns fraction of populations the produce emission */
00079 double emit_frac( transition *t )
00080 {
00081         DEBUG_ENTRY( "emit_frac()" );
00082 
00083         ASSERT( t->ipCont > 0 );
00084 
00085         /* collisional deexcitation and destruction by background opacities
00086          * are loss of photons without net emission */
00087         double deexcit_loss = t->Coll.cs * dense.cdsqte + t->Emis->Aul*t->Emis->Pdest;
00088         /* this is what is observed */
00089         double rad_deexcit = t->Emis->Aul*(t->Emis->Pelec_esc + t->Emis->Pesc);
00090         return rad_deexcit/(deexcit_loss + rad_deexcit);
00091 }
00092 
00093 /*DumpLine print various information about an emission line vector, 
00094  * used in debugging, print to std out, ioQQQ */
00095 void DumpLine(transition * t)
00096 {
00097         char chLbl[11];
00098 
00099         DEBUG_ENTRY( "DumpLine()" );
00100 
00101         ASSERT( t->ipCont > 0 );
00102 
00103         /* routine to print contents of line arrays */
00104         strcpy( chLbl, chLineLbl(t) );
00105 
00106         fprintf( ioQQQ, 
00107                 "   %10.10s Te%.2e eden%.1e CS%.2e Aul%.1e Tex%.2e cool%.1e het%.1e conopc%.1e albdo%.2e\n", 
00108           chLbl, 
00109           phycon.te, 
00110           dense.eden, 
00111           t->Coll.cs, 
00112           t->Emis->Aul, 
00113           TexcLine(t), 
00114           t->Coll.cool, 
00115           t->Coll.heat ,
00116           opac.opacity_abs[t->ipCont-1],
00117           opac.albedo[t->ipCont-1]);
00118 
00119         fprintf( ioQQQ, 
00120                 "  Tin%.1e Tout%.1e Esc%.1e eEsc%.1e DesP%.1e Pump%.1e OTS%.1e PopL,U %.1e %.1e PopOpc%.1e\n", 
00121           t->Emis->TauIn, 
00122           t->Emis->TauTot, 
00123           t->Emis->Pesc, 
00124           t->Emis->Pelec_esc, 
00125           t->Emis->Pdest, 
00126           t->Emis->pump, 
00127           t->Emis->ots, 
00128           t->Lo->Pop, 
00129           t->Hi->Pop ,
00130           t->Emis->PopOpc );
00131         return;
00132 }
00133 
00134 
00135 /*OccupationNumberLine - derive the photon occupation number at line center for any line */
00136 double OccupationNumberLine( transition * t )
00137 {
00138         double OccupationNumberLine_v;
00139 
00140         DEBUG_ENTRY( "OccupationNumberLine()" );
00141 
00142         ASSERT( t->ipCont > 0 );
00143 
00144         /* routine to evaluate line photon occupation number */
00145         if( t->Lo->Pop > SMALLFLOAT )
00146         {
00147                 /* the lower population with correction for stimulated emission */
00148                 double PopLo_corr = t->Lo->Pop / t->Lo->g - t->Hi->Pop / t->Hi->g;
00149                 OccupationNumberLine_v = ( t->Hi->Pop / t->Hi->g )/SDIV(PopLo_corr ) *
00150                         (1. - t->Emis->Pesc);
00151         }
00152         else
00153         {
00154                 OccupationNumberLine_v = 0.;
00155         }
00156         return( OccupationNumberLine_v );
00157 }
00158 
00159 /*TexcLine derive excitation temperature of line from contents of line array */
00160 double TexcLine(transition * t)
00161 {
00162         double TexcLine_v;
00163 
00164         DEBUG_ENTRY( "TexcLine()" );
00165 
00166         /* routine to evaluate line excitation temp using contents of line array
00167          * */
00168         if( t->Hi->Pop * t->Lo->Pop > 0. )
00169         {
00170                 TexcLine_v = ( t->Hi->Pop / t->Hi->g )/( t->Lo->Pop / t->Lo->g );
00171                 TexcLine_v = log(TexcLine_v);
00172                 /* protect against infinite temp limit */
00173                 if( fabs(TexcLine_v) > SMALLFLOAT )
00174                 {
00175                         TexcLine_v = - t->EnergyK / TexcLine_v;
00176                 }
00177         }
00178         else
00179         {
00180                 TexcLine_v = 0.;
00181         }
00182         return( TexcLine_v );
00183 }
00184 
00185 /*eina convert a gf into an Einstein A */
00186 double eina(double gf,
00187           double enercm, 
00188           double gup)
00189 {
00190         double eina_v;
00191 
00192         DEBUG_ENTRY( "eina()" );
00193 
00194         /* derive the transition prob, given the following
00195          * call to function is gf, energy in cm^-1, g_up
00196          * gf is product of g and oscillator strength
00197          * eina = ( gf / 1.499e-8 ) / (wl/1e4)**2 / gup  */
00198         eina_v = (gf/gup)*TRANS_PROB_CONST*POW2(enercm);
00199         return( eina_v );
00200 }
00201 
00202 /*GetGF convert Einstein A into oscillator strength */
00203 double GetGF(double trans_prob, 
00204           double enercm, 
00205           double gup)
00206 {
00207         double GetGF_v;
00208 
00209         DEBUG_ENTRY( "GetGF()" );
00210 
00211         ASSERT( enercm > 0. );
00212         ASSERT( trans_prob > 0. );
00213         ASSERT( gup > 0.);
00214 
00215         /* derive the transition prob, given the following
00216          * call to function is gf, energy in cm^-1, g_up
00217          * gf is product of g and oscillator strength
00218          * trans_prob = ( GetGF/gup) / 1.499e-8 / ( 1e4/enercm )**2 */
00219         GetGF_v = trans_prob*gup/TRANS_PROB_CONST/POW2(enercm);
00220         return( GetGF_v );
00221 }
00222 
00223 /*abscf convert gf into absorption coefficient */
00224 double abscf(double gf, 
00225           double enercm, 
00226           double gl)
00227 {
00228         double abscf_v;
00229 
00230         DEBUG_ENTRY( "abscf()" );
00231 
00232         ASSERT(gl > 0. && enercm > 0. && gf > 0. );
00233 
00234         /* derive line absorption coefficient, given the following:
00235          * gf, enercm, g_low
00236          * gf is product of g and oscillator strength */
00237         abscf_v = 1.4974e-6*(gf/gl)*(1e4/enercm);
00238         return( abscf_v );
00239 }
00240 
00241 /*chIonLbl use information in line array to generate a null terminated ion label in "Fe 2" */
00242 void chIonLbl(char *chIonLbl_v , transition * t )
00243 {
00244         /*static char chIonLbl_v[5];*/
00245 
00246         DEBUG_ENTRY( "chIonLbl()" );
00247 
00248         /* function to use information within the line array
00249          * to generate an ion label, giving element and
00250          * ionization stage
00251          * */
00252         if( t->Hi->nelem < 0 )
00253         {
00254                 /* this line is to be ignored */
00255                 strcpy( chIonLbl_v, "Dumy" );
00256         }
00257         else if( t->Hi->nelem-1 >= LIMELM )
00258         {
00259                 /* this is one of the molecules, either 12CO or 13CO */
00260 
00261                 /* >>chng 02 may 15, from to chElementNameShort to go mole right */
00262                 strcpy( chIonLbl_v , elementnames.chElementNameShort[t->Hi->nelem-1] );
00263 
00264                 /* chIonStage is four char null terminated string, starting with "_1__" 
00265                 strcat( chIonLbl_v, "CO");*/
00266         }
00267 
00268         else
00269         {
00270                 ASSERT( t->Hi->nelem >= 1 );
00271                 /* ElmntSym.chElementSym is null terminated, 2 ch + null, string giving very
00272                  * short form of element name */
00273                 strcpy( chIonLbl_v , elementnames.chElementSym[t->Hi->nelem -1] );
00274 
00275                 /* chIonStage is four char null terminated string, starting with "_1__" */
00276                 strcat( chIonLbl_v, elementnames.chIonStage[t->Hi->IonStg-1]);
00277         }
00278         /* chIonLbl is four char null terminated string */
00279         return/*( chIonLbl_v )*/;
00280 }
00281 
00282 /*chLineLbl use information in line transfer arrays to generate a line label */
00283 /* this label is null terminated */
00284 /* ContCreatePointers has test this with full range of wavelengths */
00285 char* chLineLbl(transition * t )
00286 {
00287         static char chLineLbl_v[11];
00288 
00289         DEBUG_ENTRY( "chLineLbl()" );
00290 
00291 
00292         /* function to use information within the line array
00293          * to generate a line label, giving element and
00294          * ionization stage
00295          * rhs are set in large block data */
00296 
00297         /* NB this function is profoundly slow due to sprintf statement
00298          * also - it cannot be evaluated within a write statement itself*/
00299 
00300         if( t->WLAng > (realnum)INT_MAX )
00301         {
00302                 sprintf( chLineLbl_v, "%2.2s%2.2s%5i%c", 
00303                         elementnames.chElementSym[t->Hi->nelem -1], 
00304                         elementnames.chIonStage[t->Hi->IonStg-1], 
00305                    (int)(t->WLAng/1e8), 'c' );
00306         }
00307         else if( t->WLAng > 9999999. )
00308         {
00309                 /* wavelength is very large, convert to centimeters */
00310                 sprintf( chLineLbl_v, "%2.2s%2.2s%5.2f%c", 
00311                         elementnames.chElementSym[t->Hi->nelem -1], 
00312                         elementnames.chIonStage[t->Hi->IonStg-1], 
00313                         t->WLAng/1e8, 'c' );
00314         }
00315         else if( t->WLAng > 999999. )
00316         {
00317                 /* wavelength is very large, convert to microns */
00318                 sprintf( chLineLbl_v, "%2.2s%2.2s%5i%c", 
00319                         elementnames.chElementSym[t->Hi->nelem -1], 
00320                         elementnames.chIonStage[t->Hi->IonStg-1], 
00321                         (int)(t->WLAng/1e4), 'm' );
00322         }
00323         else if( t->WLAng > 99999. )
00324         {
00325                 /* wavelength is very large, convert to microns */
00326                 sprintf( chLineLbl_v, "%2.2s%2.2s%5.1f%c", 
00327                         elementnames.chElementSym[t->Hi->nelem -1], 
00328                         elementnames.chIonStage[t->Hi->IonStg-1], 
00329                         t->WLAng/1e4, 'm' );
00330         }
00331         else if( t->WLAng > 9999. )
00332         {
00333                 sprintf( chLineLbl_v, "%2.2s%2.2s%5.2f%c", 
00334                         elementnames.chElementSym[ t->Hi->nelem -1], 
00335                         elementnames.chIonStage[t->Hi->IonStg-1], 
00336                    t->WLAng/1e4, 'm' );
00337         }
00338         else if( t->WLAng >= 100. )
00339         {
00340                 sprintf( chLineLbl_v, "%2.2s%2.2s%5i%c", 
00341                         elementnames.chElementSym[ t->Hi->nelem -1], 
00342                         elementnames.chIonStage[t->Hi->IonStg-1], 
00343                    (int)t->WLAng, 'A' );
00344         }
00345         /* the following two formats should be changed for the
00346          * wavelength to get more precision */
00347         else if( t->WLAng >= 10. )
00348         {
00349                 sprintf( chLineLbl_v, "%2.2s%2.2s%5.1f%c", 
00350                         elementnames.chElementSym[ t->Hi->nelem -1], 
00351                         elementnames.chIonStage[t->Hi->IonStg-1], 
00352                    t->WLAng, 'A' );
00353         }
00354         else
00355         {
00356                 sprintf( chLineLbl_v, "%2.2s%2.2s%5.2f%c", 
00357                         elementnames.chElementSym[ t->Hi->nelem -1], 
00358                         elementnames.chIonStage[t->Hi->IonStg-1], 
00359                    t->WLAng, 'A' );
00360         }
00361         /* make sure that string ends with null character - this should
00362          * be redundant */
00363         ASSERT( chLineLbl_v[10]==0 );
00364         return( chLineLbl_v );
00365 }
00366 
00367 /*RefIndex calculates the index of refraction of air using the line energy in wavenumbers,
00368  * used to convert vacuum wavelengths to air wavelengths. */
00369 double RefIndex(double EnergyWN )
00370 {
00371         double RefIndex_v, 
00372           WaveMic, 
00373           xl, 
00374           xn;
00375 
00376         DEBUG_ENTRY( "RefIndex()" );
00377 
00378         ASSERT( EnergyWN > 0. );
00379 
00380         /* the wavelength in microns */
00381         WaveMic = 1.e4/EnergyWN;
00382 
00383         /* only do index of refraction if longward of 2000A */
00384         if( WaveMic > 0.2 )
00385         {
00386                 /* longward of 2000A
00387                  * xl is 1/WaveMic^2 */
00388                 xl = 1.0/WaveMic/WaveMic;
00389                 /* use a formula from 
00390                  *>>refer       air     index refraction        Allen, C.W. 1973, Astrophysical quantities, 
00391                  *>>refercon    3rd Edition (AQ), p.124 */
00392                 xn = 255.4/(41. - xl);
00393                 xn += 29498.1/(146. - xl);
00394                 xn += 64.328;
00395                 RefIndex_v = xn/1.e6 + 1.;
00396         }
00397         else
00398         {
00399                 RefIndex_v = 1.;
00400         }
00401         ASSERT( RefIndex_v > 0. );
00402         return( RefIndex_v );
00403 }
00404 
00405 /*PutCS enter a collision strength into an individual line vector */
00406 void PutCS(double cs, 
00407   transition * t)
00408 {
00409 
00410         DEBUG_ENTRY( "PutCS()" );
00411 
00412         /* collision strength must not be negative, had been test for being positive,
00413          * but called with zero - did not check why 98 jul 5 */
00414         ASSERT( cs >= 0. );
00415 
00416         t->Coll.cs = (realnum)cs;
00417         return;
00418 }
00419 
00420 /*WavlenErrorGet - given the real wavelength in A for a line
00421  * routine will find the error expected between the real 
00422  * wavelength and the wavelength printed in the output, with 4 sig figs,
00423  * function returns difference between exact and 4 sig fig wl, so 
00424  * we have found correct line is fabs(d wl) < return */
00425 realnum WavlenErrorGet( realnum wavelength )
00426 {
00427         double a;
00428         realnum errorwave;
00429 
00430         DEBUG_ENTRY( "WavlenErrorGet()" );
00431 
00432         ASSERT( LineSave.sig_figs <= 6 );
00433 
00434         if( wavelength > 0. )
00435         {
00436                 /* normal case, positive (non zero) wavelength */
00437                 a = log10( wavelength+FLT_EPSILON);
00438                 a = floor(a);
00439         }
00440         else
00441         {
00442                 /* might be called with wl of zero, this is that case */
00443                 /* errorwave = 1e-4f; */
00444                 a = 0.;
00445         }
00446 
00447         errorwave = 5.f * (realnum)pow( 10., a - (double)LineSave.sig_figs );
00448         return errorwave;
00449 }
00450 
00451 /*linadd enter lines into the line storage array, called once per zone for each line*/
00452 void linadd(
00453   double xInten,        /* xInten - local emissivity per unit vol, no fill fac */
00454   realnum wavelength,   /* realnum wavelength */
00455   const char *chLab,/* string label for ion */
00456   char chInfo,          /* character type of entry for line - given below */
00457                                         /* 'c' cooling, 'h' heating, 'i' info only, 'r' recom line */
00458   const char *chComment )
00459 {
00460         DEBUG_ENTRY( "linadd()" );
00461 
00462         /* main routine to actually enter lines into the line storage array
00463          * called at top level within routine lines
00464          * called series of times in routine PutLine for lines transferred
00465          */
00466 
00467         /* three values, -1 is just counting, 0 if init, 1 for calculation */
00468         if( LineSave.ipass > 0 )
00469         {
00470                 /* not first pass, sum lines only
00471                  * total emission from vol */
00472                 LineSv[LineSave.nsum].sumlin[0] += xInten*radius.dVeff;
00473                 /* local emissivity in line */
00474                 LineSv[LineSave.nsum].emslin[0] = xInten;
00475                 /* no need to increment or set [1] version since this is called with no continuum
00476                  * index, don't know what to do */
00477                 if( wavelength > 0 )
00478                 {
00479                         /* only put informational lines, like "Q(H) 4861", in this stack
00480                          * many continua have a wavelength of zero and are proper intensities,
00481                          * it would be wrong to predict their transferred intensity */
00482                         LineSv[LineSave.nsum].emslin[1] = LineSv[LineSave.nsum].emslin[0];
00483                         LineSv[LineSave.nsum].sumlin[1] = LineSv[LineSave.nsum].sumlin[0]; 
00484                 }
00485         }
00486 
00487         else if( LineSave.ipass == 0 )
00488         {
00489                 /* first call to stuff lines in array, confirm that label is one of
00490                  * the four correct ones */
00491                 /* >>chng 04 apr 24, last two had been != so this test never really happened; PvH */
00492                 ASSERT( (chInfo == 'c') || (chInfo == 'h') || (chInfo == 'i') || (chInfo == 'r' ) );
00493 
00494                 /* then save it into array */
00495                 LineSv[LineSave.nsum].chSumTyp = (char)chInfo;
00496 
00497                 /* number of lines ok, set parameters for first pass */
00498                 LineSv[LineSave.nsum].sumlin[0] = 0.;
00499                 LineSv[LineSave.nsum].emslin[0] = 0.;
00500                 LineSv[LineSave.nsum].wavelength = wavelength;
00501                 LineSv[LineSave.nsum].sumlin[1] = 0.;
00502                 LineSv[LineSave.nsum].emslin[1] = 0.;
00503                 LineSv[LineSave.nsum].chComment = chComment;
00504                 /* check that null is correct, string overruns have 
00505                  * been a problem in the past */
00506                 ASSERT( strlen( chLab )<5 );
00507                 strcpy( LineSv[LineSave.nsum].chALab, chLab );
00508         }
00509 
00510         /* increment the line counter */
00511         ++LineSave.nsum;
00512 
00513         /* routine can be called with negative LineSave.ipass, in this case
00514          * we are just counting number of lines for current setup */
00515         return;
00516 }
00517 
00518 
00519 /* this is energy in Ryd as set by call to PntForLine */
00520 static double EnergyRyd;
00521 
00522 /*emergent_line find absorption due to continuous opacity */
00523 double emergent_line( 
00524         /* emissivity [erg cm-3 s-1] in inward direction */
00525         double emissivity_in , 
00526         /* emissivity [erg cm-3 s-1] in outward direction */
00527         double emissivity_out , 
00528         /* array index for continuum frequency */
00529         long int ipCont )
00530 {
00531 
00532         double emergent_in , emergent_out;
00533         long int i = ipCont-1;
00534 
00535         DEBUG_ENTRY( "emergent_line()" );
00536 
00537         ASSERT( i >= 0 && i < rfield.nupper-1 );
00538 
00539         /* do nothing if first iteration since we do not know the outward-looking
00540          * optical depths.  In version C07.02.00 we assumed an infinite optical
00541          * depth in the outward direction, which would be appropriate for a 
00542          * HII region on the surface of a molecular cloud.  This converged onto
00543          * the correct solution in later iterations, but on the first iteration
00544          * this underestimated total emission if the infinite cloud were not
00545          * present.  With C07.02.01 we make no assuptions about what is in the
00546          * outward direction and simply use the local emission. 
00547          * Behavior is unchanged on later iterations */
00548         if( iteration == 1 )
00549         {
00550                 /* first iteration - do not know outer optical depths so assume very large 
00551                  * optical depths */
00552                 emergent_in = emissivity_in*opac.E2TauAbsFace[i];
00553                 emergent_out = emissivity_out;
00554         }
00555         else
00556         {
00557                 if( geometry.lgSphere )
00558                 {
00559                         /* second or later iteration in closed or spherical geometry */
00560                         /* inwardly directed emission must get to central hole then across entire
00561                          * far side of shell */
00562                         emergent_in = emissivity_in  * opac.E2TauAbsFace[i] *opac.E2TauAbsTotal[i];
00563 
00564                         /* E2 is outwardly directed emission to get to outer edge of cloud */
00565                         emergent_out = emissivity_out * opac.E2TauAbsOut[i];
00566                 }
00567                 else
00568                 {
00569                         /* open geometry in second or later iteration, outer optical depths are known 
00570                          * this is light emitted into the outer direction and backscattered
00571                          * into the inner */
00572                         double reflected = emissivity_out * opac.albedo[i] * (1.-opac.E2TauAbsOut[i]);
00573                         /* E2 is to get to central hole */
00574                         emergent_in = (emissivity_in + reflected) * opac.E2TauAbsFace[i];
00575                         /* E2 is to get to outer edge */
00576                         emergent_out = (emissivity_out - reflected) * opac.E2TauAbsOut[i];
00577                 }
00578         }
00579         /* return the net emissivity times a vol element */
00580         return( (emergent_in + emergent_out)  );
00581 }
00582 
00583 /*lindst add line with destruction and outward */
00584 void lindst(
00585   double xInten, 
00586   realnum wavelength, 
00587   const char *chLab, 
00588   long int ipnt, 
00589   char chInfo, 
00590   bool lgOutToo,
00591   const char *chComment )
00592 {
00593         double saveemis;
00594         DEBUG_ENTRY( "lindst()" );
00595 
00596         /* >>chng 06 feb 08, add test on xInten positive, no need to evaluate
00597          * for majority of zero */
00598         if( LineSave.ipass > 0 && xInten > 0. )
00599         {
00600                 /* LineSave.ipass > 0, integration across simulation, sum lines only 
00601                  * emissivity, emission per unit vol, for this zone */
00602                 LineSv[LineSave.nsum].emslin[0] = xInten;
00603                 /* integrated intensity or luminosity, the emissivity times the volume */
00604                 LineSv[LineSave.nsum].sumlin[0] += xInten*radius.dVeff;
00605 
00606                 /* add line to outward beam 
00607                  * there are lots of lines that are sums of other lines, or
00608                  * just for info of some sort.  These have flag lgOutToo false.
00609                  * Note that the EnergyRyd variable only has a rational
00610                  * value if PntForLine was called just before this routine - in all
00611                  * cases where this did not happen the flag is false. */
00612                 if( lgOutToo )
00613                 {
00614                         rfield.outlin[ipnt-1] += (realnum)(xInten/(rfield.anu[ipnt-1]*EN1RYD)*
00615                           radius.dVolOutwrd*opac.ExpZone[ipnt-1]);
00616                         rfield.reflin[ipnt-1] += (realnum)(xInten/(rfield.anu[ipnt-1]*EN1RYD)*
00617                           radius.dVolReflec);
00618                 }
00619                 /* main production pass, sum lines only */
00620                 if( ipnt <= rfield.nflux )
00621                 {
00622                         /* emergent_line accounts for destruction by absorption outside
00623                          * the line-forming region */
00624                         saveemis = emergent_line( 
00625                                 xInten*rt.fracin , xInten*(1.-rt.fracin) , ipnt );
00626                         LineSv[LineSave.nsum].emslin[1] = saveemis;
00627                         LineSv[LineSave.nsum].sumlin[1] += saveemis*radius.dVeff; 
00628                 }
00629         }
00630         else if( LineSave.ipass == 0 )
00631         {
00632                 ASSERT( ipnt > 0 );
00633 
00634                 /* first call to stuff lines in array, confirm that label is one of
00635                  * the four correct ones */
00636                 /* >>chng 04 apr 24, last two had been != so this test never really happened; PvH */
00637                 ASSERT( (chInfo == 'c') || (chInfo == 'h') || (chInfo == 'i') || (chInfo == 'r' ) );
00638                 LineSv[LineSave.nsum].chSumTyp = (char)chInfo;
00639                 /* number of lines ok, set parameters for first pass */
00640                 LineSv[LineSave.nsum].sumlin[0] = 0.;
00641                 LineSv[LineSave.nsum].emslin[0] = 0.;
00642                 LineSv[LineSave.nsum].wavelength = wavelength;
00643                 LineSv[LineSave.nsum].emslin[1] = 0.;
00644                 LineSv[LineSave.nsum].sumlin[1] = 0.;
00645                 LineSv[LineSave.nsum].chComment = chComment;
00646 
00647                 /* check that null is correct, string overruns have 
00648                 * been a problem in the past */
00649                 ASSERT( strlen( chLab )<5 );
00650                 strcpy( LineSv[LineSave.nsum].chALab, chLab );
00651         }
00652 
00653         /* increment the line counter */
00654         ++LineSave.nsum;
00655 
00656         EnergyRyd = 0.;
00657         return;
00658 }
00659 
00660 /*PntForLine generate pointer for forbidden line */
00661 void PntForLine(
00662   /* wavelength of transition in Angstroms */
00663   double wavelength, 
00664   /* label for this line */
00665   const char *chLabel,
00666   /* this is array index on the f, not c scale,
00667    * for the continuum cell holding the line */
00668   long int *ipnt)
00669 {
00670         /* 
00671          * maximum number of forbidden lines - this is a good bet since
00672          * new lines do not go into this group, and lines are slowly 
00673          * moving to level 1 
00674          */
00675 #       define  MAXFORLIN       1000
00676         static long int ipForLin[MAXFORLIN]={0};
00677 
00678         /* number of forbidden lines entered into continuum array */
00679         static long int nForLin;
00680 
00681         DEBUG_ENTRY( "PntForLine()" );
00682 
00683         /* must be 0 or greater */
00684         ASSERT( wavelength >= 0. );
00685 
00686         if( wavelength == 0. )
00687         {
00688                 /* zero is special flag to initialize */
00689                 nForLin = 0;
00690                 /* ipLineEnergy will only put in line label if nothing already there */
00691                 EnergyRyd = 0.;
00692         }
00693         else
00694         {
00695 
00696                 if( LineSave.ipass > 0 )
00697                 {
00698                         /* not first pass, sum lines only */
00699                         *ipnt = ipForLin[nForLin];
00700                 }
00701                 else if( LineSave.ipass == 0 )
00702                 {
00703                         /* check if number of lines in arrays exceeded */
00704                         if( nForLin >= MAXFORLIN )
00705                         {
00706                                 fprintf( ioQQQ, "PROBLEM %5ld lines is too many for PntForLine.\n", 
00707                                   nForLin );
00708                                 fprintf( ioQQQ, " Increase the value of maxForLine everywhere in the code.\n" );
00709                                 cdEXIT(EXIT_FAILURE);
00710                         }
00711 
00712                         /* ipLineEnergy will only put in line label if nothing already there */
00713                         EnergyRyd = RYDLAM/wavelength;
00714                         ipForLin[nForLin] = ipLineEnergy(EnergyRyd,chLabel , 0);
00715                         *ipnt = ipForLin[nForLin];
00716                 }
00717                 else
00718                 {
00719                         /* this is case where we are only counting lines */
00720                         *ipnt = 0;
00721                 }
00722                 ++nForLin;
00723         }
00724         return;
00725 }
00726 
00727 static realnum ExtraInten;
00728 
00729 /*PutLine enter local line intensity into the intensity stack for eventual printout */
00730 void PutLine(transition * t, const char *chComment, const char *chLabelTemp)
00731 {
00732         char chLabel[5];
00733         realnum wl;
00734         double xIntensity,
00735                 other,
00736                 xIntensity_in;
00737 
00738         DEBUG_ENTRY( "PutLine()" );
00739 
00740         /* routine to use line array data to generate input
00741          * for emission line array */
00742         ASSERT( t->ipCont > 0. );
00743 
00744         strcpy( chLabel, chLabelTemp );
00745         chLabel[4] = '\0';
00746 
00747         /* if ipass=0 then we must generate label info since first pass
00748          * gt.0 then only need line intensity data */
00749         if( LineSave.ipass == 0 )
00750         {
00751                 /* save wavelength */
00752                 wl = t->WLAng;
00753 
00754                 xIntensity = 0.;
00755         }
00756         else
00757         {
00758                 /* both the counting and integrating modes comes here */
00759                 /*>>chng 06 mar 10, these are not actually used so set to impossible values */
00760                 /*strcpy( chLabel, "    " );*/
00761                 chLabel[0] = '\n';
00762                 wl = 0.;
00763                 /* total line intensity or luminosity 
00764                  * these may not be defined in initial calls so define here */
00765                 xIntensity = t->Emis->xIntensity + ExtraInten;
00766         }
00767         /* initial counting case, where ipass == -1, just ignored above, call linadd below */
00768 
00769         /* ExtraInten is option that allows extra intensity (i.e., recomb)
00770          * to be added to this line  with Call PutExtra( exta )
00771          * in main lines this extra
00772          * contribution must be identified explicitly */
00773         ExtraInten = 0.;
00774         /*linadd(xIntensity,wl,chLabel,'i');*/
00775         /*lindst add line with destruction and outward */
00776         rt.fracin = t->Emis->FracInwd;
00777         lindst(xIntensity, 
00778           wl, 
00779           chLabel, 
00780           t->ipCont, 
00781           /* this is information only - has been counted in cooling already */
00782           'i', 
00783           /* do not add to outward beam, also done separately */
00784           false,
00785           chComment);
00786         rt.fracin = 0.5;
00787 
00788         /* inward part of line - do not move this away from previous lines
00789          * since xIntensity is used here */
00790         xIntensity_in = xIntensity*t->Emis->FracInwd;
00791         linadd(xIntensity_in,wl,"Inwd",'i',chComment);
00792 
00793         /* cooling part of line */
00794         other = t->Coll.cool;
00795         linadd(other,wl,"Coll",'i',chComment);
00796 
00797         /* fluorescent excited part of line */
00798         other = t->Emis->PopOpc * t->Emis->pump * (1.-t->Emis->ColOvTot) * t->EnergyErg;
00799         linadd(other,wl,"Pump",'i',chComment);
00800 
00801         /* heating part of line */
00802         other = t->Coll.heat;
00803         linadd(other,wl,"Heat",'i',chComment);
00804         return;
00805 }
00806 
00807 /*PutLine enter local line intensity into the intensity stack for eventual printout */
00808 void PutLine(transition * t,
00809   const char *chComment)
00810 {
00811         char chLabel[5];
00812         realnum wl;
00813         double xIntensity,
00814                 other,
00815                 xIntensity_in;
00816 
00817         DEBUG_ENTRY( "PutLine()" );
00818 
00819         /* routine to use line array data to generate input
00820          * for emission line array */
00821         ASSERT( t->ipCont > 0. );
00822 
00823         /* if ipass=0 then we must generate label info since first pass
00824          * gt.0 then only need line intensity data */
00825         if( LineSave.ipass == 0 )
00826         {
00827                 /* these variables not used by linadd if ipass ne 0 */
00828                 /* chIonLbl is function that generates a null terminated 4 char string, of form "C  2" */
00829                 chIonLbl(chLabel, t);
00830 
00831                 /* save wavelength */
00832                 wl = t->WLAng;
00833 
00834                 xIntensity = 0.;
00835         }
00836         else
00837         {
00838                 /* both the counting and integrating modes comes here */
00839                 /*>>chng 06 mar 10, these are not actually used so set to impossible values */
00840                 /*strcpy( chLabel, "    " );*/
00841                 chLabel[0] = '\n';
00842                 wl = 0.;
00843                 /* total line intensity or luminosity 
00844                  * these may not be defined in initial calls so define here */
00845                 xIntensity = t->Emis->xIntensity + ExtraInten;
00846         }
00847         /* initial counting case, where ipass == -1, just ignored above, call linadd below */
00848 
00849         /* ExtraInten is option that allows extra intensity (i.e., recomb)
00850          * to be added to this line  with Call PutExtra( exta )
00851          * in main lines this extra
00852          * contribution must be identified explicitly */
00853         ExtraInten = 0.;
00854         /*linadd(xIntensity,wl,chLabel,'i');*/
00855         /*lindst add line with destruction and outward */
00856         rt.fracin = t->Emis->FracInwd;
00857         lindst(xIntensity, 
00858           wl, 
00859           chLabel, 
00860           t->ipCont, 
00861           /* this is information only - has been counted in cooling already */
00862           'i', 
00863           /* do not add to outward beam, also done separately */
00864           false,
00865           chComment);
00866         rt.fracin = 0.5;
00867 
00868         /* inward part of line - do not move this away from previous lines
00869          * since xIntensity is used here */
00870         xIntensity_in = xIntensity*t->Emis->FracInwd;
00871         linadd(xIntensity_in,wl,"Inwd",'i',chComment);
00872 
00873         /* cooling part of line */
00874         other = t->Coll.cool;
00875         linadd(other,wl,"Coll",'i',chComment);
00876 
00877         /* fluorescent excited part of line */
00878         other = t->Emis->PopOpc * t->Emis->pump * (1.-t->Emis->ColOvTot) * t->EnergyErg;
00879         linadd(other,wl,"Pump",'i',chComment);
00880 
00881         /* heating part of line */
00882         other = t->Coll.heat;
00883         linadd(other,wl,"Heat",'i',chComment);
00884         return;
00885 }
00886 
00887 /* ==================================================================== */
00888 /*PutExtra enter and 'extra' intensity source for some line */
00889 void PutExtra(double Extra)
00890 {
00891 
00892         DEBUG_ENTRY( "PutExtra()" );
00893 
00894         ExtraInten = (realnum)Extra;
00895         return;
00896 }
00897 
00898 void TransitionJunk( transition *t )
00899 {
00900 
00901         DEBUG_ENTRY( "TransitionJunk()" );
00902 
00903                 /* wavelength, usually in A, used for printout */
00904         t->WLAng = -FLT_MAX;
00905 
00906                 /* transition energy in degrees kelvin*/
00907         t->EnergyK = -FLT_MAX;
00908                 /* transition energy in ergs */
00909         t->EnergyErg = -FLT_MAX;
00910                 /* transition energy in wavenumbers */
00911         t->EnergyWN = -FLT_MAX;
00912 
00913         /* array offset for radiative transition within continuum array 
00914          * is negative if transition is non-radiative. */
00915         t->ipCont = -10000;
00916 
00917         CollisionJunk( &t->Coll );
00918 
00919         /* set these equal to NULL first. That will cause the code to crash if
00920          * the variables are ever used before being deliberately set. */
00921         t->Emis = &DummyEmis;
00922         
00923         t->Lo = NULL;
00924         t->Hi = NULL;
00925         return;
00926 }
00927 
00928 
00929 /*EmLineJunk set all elements of transition struc to dangerous values */
00930 void EmLineJunk( emission *t )
00931 {
00932 
00933         DEBUG_ENTRY( "EmLineJunk()" );
00934 
00935         /* optical depth in continuum to ill face */
00936         t->TauCon = -FLT_MAX;
00937 
00938         /* inward and total line optical depths */
00939         t->TauIn = -FLT_MAX;
00940         t->TauTot = -FLT_MAX;
00941 
00942         /* type of redistribution function, */
00943         t->iRedisFun = INT_MIN;
00944 
00945         /* array offset for line within fine opacity */
00946         t->ipFine = -10000;
00947 
00948         /* inward fraction */
00949         t->FracInwd = -FLT_MAX;
00950 
00951         /* continuum pumping rate */
00952         t->pump = -FLT_MAX;
00953 
00954         /* line intensity */
00955         t->xIntensity = -FLT_MAX;
00956 
00957         /* number of photons emitted per sec in the line */
00958         t->phots = -FLT_MAX;
00959 
00960         /* gf value */
00961         t->gf = -FLT_MAX;
00962 
00963         /* escape and destruction probs */
00964         t->Pesc = -FLT_MAX;
00965         t->Pdest = -FLT_MAX;
00966         t->Pelec_esc = -FLT_MAX;
00967 
00968         /* damping constant, and number related to it */
00969         t->dampXvel = -FLT_MAX;
00970         t->damp = -FLT_MAX;
00971 
00972         /* ratio of collisional to radiative excitation*/
00973         t->ColOvTot = -FLT_MAX;
00974 
00975         /* line opacity */
00976         t->opacity = -FLT_MAX;
00977 
00978         t->PopOpc = -FLT_MAX;
00979 
00980         /* transition prob, Einstein A upper to lower */
00981         t->Aul = -FLT_MAX;
00982 
00983         /* ots rate */
00984         t->ots = -FLT_MAX;
00985         return;
00986 }
00987 
00988 /*CollisionJunk set all elements of transition struc to dangerous values */
00989 void CollisionJunk( collision * t )
00990 {
00991 
00992         DEBUG_ENTRY( "CollisionJunk()" );
00993 
00995         t->ColUL = -FLT_MAX;
00996 
00997         /* Coll->cooling and Coll->heating due to collisional excitation */
00998         t->cool = -FLT_MAX;
00999         t->heat = -FLT_MAX;
01000 
01001         /* collision strengths for transition */
01002         t->cs = -FLT_MAX;
01003         for( long i=0; i<ipNCOLLIDER; i++ )
01004                 t->csi[i] = -FLT_MAX;
01005         return;
01006 }
01007 
01008 /*StateJunk set all elements of transition struc to dangerous values */
01009 void StateJunk( quantumState * t )
01010 {
01011 
01012         DEBUG_ENTRY( "StateJunk()" );
01013 
01014         /* t->chLabel = ''; */
01015 
01017         t->g = -FLT_MAX;
01018 
01020         t->Pop = -FLT_MAX;
01021 
01023         t->IonStg = -10000;
01024 
01026         t->nelem = -10000;
01027         return;
01028 }
01029 
01030 /*TransitionZero zeros out TransitionZero */
01031 void TransitionZero( transition *t )
01032 {
01033 
01034         DEBUG_ENTRY( "TransitionZero()" );
01035 
01036         CollisionZero( &t->Coll );
01037 
01038         StateZero( t->Lo );
01039         StateZero( t->Hi );
01040         EmLineZero( t->Emis );
01041 
01042         return;
01043 }
01044 
01045 /*EmLineZero zeros out the emission line structure */
01046 void EmLineZero( emission * t )
01047 {
01048 
01049         DEBUG_ENTRY( "EmLineZero()" );
01050 
01051         /* total optical depth in all overlapping lines to illuminated face,
01052          * used for pumping */
01053         t->TauCon = opac.taumin;
01054 
01055         /* inward and total line optical depths */
01056         /* >>chng 03 feb 14, from 0 to opac.taumin */
01057         t->TauIn = opac.taumin;
01058 
01059         /* total optical depths */
01060         t->TauTot = 1e20f;
01061 
01062         /* inward fraction */
01063         /* >>chng 03 feb 14, from 0 to 1 */
01064         t->FracInwd = 1.;
01065 
01066         /* continuum pumping rate */
01067         t->pump = 0.;
01068 
01069         /* line intensity */
01070         t->xIntensity = 0.;
01071 
01072         /* number of photons emitted per sec in the line */
01073         t->phots = 0.;
01074 
01075         /* escape and destruction probs */
01076         /* >>chng 03 feb 14, change from 0 to 1 */
01077         t->Pesc = 1.;
01078         t->Pdest = 0.;
01079         t->Pelec_esc = 0.;
01080 
01081         /* ratio of collisional to radiative excitation*/
01082         t->ColOvTot = 0.;
01083 
01084         /* pop that enters net opacity */
01085         t->PopOpc = 0.;
01086 
01087         /* ots rate */
01088         t->ots = 0.;
01089         return;
01090 }
01091 
01092 /*CollisionZero zeros out the structure */
01093 void CollisionZero( collision * t )
01094 {
01095 
01096         DEBUG_ENTRY( "CollisionZero()" );
01097 
01098         /* Coll->cooling and Coll->heating due to collisional excitation */
01099         t->cool = 0.;
01100         t->heat = 0.;
01101         return;
01102 }
01103 
01104 /*StateZero zeros out the structure */
01105 void StateZero( quantumState * t )
01106 {
01107 
01108         DEBUG_ENTRY( "StateZero()" );
01109 
01111         t->Pop = 0.;
01112         return;
01113 }
01114 
01115 /*LineConvRate2CS convert down coll rate back into electron cs in case other parts of code need this for reference */
01116 void LineConvRate2CS( transition* t , realnum rate )
01117 {
01118 
01119         DEBUG_ENTRY( "LineConvRate2CS()" );
01120 
01121         /* return is collision strength, convert from collision rate from 
01122          * upper to lower, this assumes pure electron collisions, but that will
01123          * also be assumed by anything that uses cs, for self-consistency */
01124         t->Coll.cs = rate * t->Hi->g / (realnum)dense.cdsqte;
01125 
01126         /* change assert to non-negative - there can be cases (Iin H2) where cs has
01127          * underflowed to 0 on some platforms */
01128         ASSERT( t->Coll.cs >= 0. );
01129         return;
01130 }
01131 
01132 /*ConvRate2CS convert down coll rate back into electron cs in case other parts of code need this for reference */
01133 double ConvRate2CS( realnum gHi , realnum rate )
01134 {
01135 
01136         double cs;
01137 
01138         DEBUG_ENTRY( "ConvRate2CS()" );
01139 
01140         /* return is collision strength, convert from collision rate from 
01141          * upper to lower, this assumes pure electron collisions, but that will
01142          * also be assumed by anything that uses cs, for self-consistency */
01143         cs = rate * gHi / dense.cdsqte;
01144 
01145         /* change assert to non-negative - there can be cases (Iin H2) where cs has
01146          * underflowed to 0 on some platforms */
01147         ASSERT( cs >= 0. );
01148         return cs;
01149 }
01150 
01151 /* returns true is we have not overrun optical depth scale */
01152 bool lgTauGood( transition * t)
01153 {
01154         bool lgGoodTau;
01155 
01156         DEBUG_ENTRY( "lgTauGood()" );
01157 
01158         if( (t->Emis->TauTot*0.9 - t->Emis->TauIn < 0. && t->Emis->TauIn > 0.) && 
01159             ! fp_equal( t->Emis->TauTot, opac.taumin ) )
01160         {
01161                 /* do not do anything if we have overrun the scale, */
01162                 lgGoodTau = false;
01163         }
01164         else
01165         {
01166                 lgGoodTau = true;
01167         }
01168         return lgGoodTau;
01169 }
01170 
01171 /*gbar0 compute g-bar gaunt factor for neutrals */
01172 STATIC void gbar0(double ex, 
01173           realnum *g)
01174 {
01175         double a, 
01176           b, 
01177           c, 
01178           d, 
01179           y;
01180 
01181         DEBUG_ENTRY( "gbar0()" );
01182 
01183         /* written by Dima Verner
01184          *
01185          * Calculation of the effective Gaunt-factor by use of 
01186          * >>refer      gbar    cs      Van Regemorter, H., 1962, ApJ 136, 906
01187          * fits for neutrals
01188          *  Input parameters: 
01189          * ex - energy ryd - now K
01190          * t  - temperature in K
01191          *  Output parameter:
01192          * g  - effective Gaunt factor
01193          * */
01194 
01195         /* y = ex*157813.7/te */
01196         y = ex/phycon.te;
01197         if( y < 0.01 )
01198         {
01199                 *g = (realnum)(0.29*(log(1.0+1.0/y) - 0.4/POW2(y + 1.0))/exp(y));
01200         }
01201         else
01202         {
01203                 if( y > 10.0 )
01204                 {
01205                         *g = (realnum)(0.066/sqrt(y));
01206                 }
01207                 else
01208                 {
01209                         a = 1.5819068e-02;
01210                         b = 1.3018207e00;
01211                         c = 2.6896230e-03;
01212                         d = 2.5486007e00;
01213                         d = log(y/c)/d;
01214                         *g = (realnum)(a + b*exp(-0.5*POW2(d)));
01215                 }
01216         }
01217         return;
01218 }
01219 
01220 /*gbar1 compute g-bar collision strength using Mewe approximations */
01221 STATIC void gbar1(double ex, 
01222           realnum *g)
01223 {
01224         double y;
01225 
01226         DEBUG_ENTRY( "gbar1()" );
01227 
01228         /*      *written by Dima Verner
01229          *
01230          * Calculation of the effective Gaunt-factor by use of 
01231          * >>refer      gbar    cs      Mewe,R., 1972, A&A 20, 215
01232          * fits for permitted transitions in ions MgII, CaII, FeII (delta n = 0)
01233          * Input parameters: 
01234          * ex - excitation energy in Ryd - now K
01235          * te  - temperature in K
01236          * Output parameter:
01237          * g  - effective Gaunt factor
01238          */
01239 
01240         /* y = ex*157813.7/te */
01241         y = ex/phycon.te;
01242         *g = (realnum)(0.6 + 0.28*(log(1.0+1.0/y) - 0.4/POW2(y + 1.0)));
01243         return;
01244 }
01245 
01246 /*MakeCS compute collision strength by g-bar approximations */
01247 void MakeCS(transition * t)
01248 {
01249         long int ion;
01250         double Abun, 
01251           cs;
01252         realnum
01253           gbar;
01254 
01255         DEBUG_ENTRY( "MakeCS()" );
01256 
01257         /* routine to get cs from various approximations */
01258 
01259         /* check if abundance greater than 0 */
01260         ion = t->Hi->IonStg;
01261 
01262         Abun = dense.xIonDense[ t->Hi->nelem -1 ][ ion-1 ];
01263         if( Abun <= 0. )
01264         {
01265                 gbar = 1.;
01266         }
01267         else
01268         {
01269                 /* check if neutral or ion */
01270                 if( ion == 1 )
01271                 {
01272                         /* neutral - compute gbar for eventual cs */
01273                         gbar0(t->EnergyK,&gbar);
01274                 }
01275                 else
01276                 {
01277                         /* ion - compute gbar for eventual cs */
01278                         gbar1(t->EnergyK,&gbar);
01279                 }
01280         }
01281 
01282         /* above was g-bar, convert to cs */
01283         cs = gbar*(14.5104/WAVNRYD)*t->Emis->gf/t->EnergyWN;
01284 
01285         /* stuff the cs in place */
01286         t->Coll.cs = (realnum)cs;
01287         return;
01288 }
01289 
01290 /*totlin sum total intensity of cooling, recombination, or intensity lines */
01291 double totlin(
01292         /* chInfor is 1 char, 
01293         'i' information, 
01294         'r' recombination or 
01295         'c' collision */
01296         int chInfo)
01297 {
01298         long int i;
01299         double totlin_v;
01300 
01301         DEBUG_ENTRY( "totlin()" );
01302 
01303         /* routine goes through set of entered line
01304          * intensities and picks out those which have
01305          * types agreeing with chInfo.  Valid types are
01306          * 'c', 'r', and 'i'
01307          *begin sanity check */
01308         if( (chInfo != 'i' && chInfo != 'r') && chInfo != 'c' )
01309         {
01310                 fprintf( ioQQQ, " TOTLIN does not understand chInfo=%c\n", 
01311                   chInfo );
01312                 cdEXIT(EXIT_FAILURE);
01313         }
01314         /*end sanity check */
01315 
01316         /* now find sum of lines of type chInfo */
01317         totlin_v = 0.;
01318         for( i=0; i < LineSave.nsum; i++ )
01319         {
01320                 if( LineSv[i].chSumTyp == chInfo )
01321                 {
01322                         totlin_v += LineSv[i].sumlin[0];
01323                 }
01324         }
01325         return( totlin_v );
01326 }
01327 
01328 
01329 /*FndLineHt search through line heat arrays to find the strongest heat source */
01330 void FndLineHt(long int *level, 
01331   /* this is the index of the strongest line in the array on the c scale */
01332   long int *ipStrong, 
01333   double *Strong)
01334 {
01335         long int i; 
01336 
01337         DEBUG_ENTRY( "FndLineHt()" );
01338 
01339         *Strong = 0.;
01340         *level = 0;
01341 
01342         /* do the level 1 lines, 0 is dummy line, <=nLevel1 is correct for c scale */
01343         for( i=1; i <= nLevel1; i++ )
01344         {
01345                 /* check if a line was the major heat agent */
01346                 if( TauLines[i].Coll.heat > *Strong )
01347                 {
01348                         *ipStrong = i;
01349                         *level = 1;
01350                         *Strong = TauLines[i].Coll.heat;
01351                 }
01352         }
01353 
01354         /* now do the level 2 lines */
01355         for( i=0; i < nWindLine; i++ )
01356         {
01357                 if( TauLine2[i].Hi->IonStg < TauLine2[i].Hi->nelem+1-NISO )
01358                 {
01359                         /* check if a line was the major heat agent */
01360                         if( TauLine2[i].Coll.heat > *Strong )
01361                         {
01362                                 *ipStrong = i;
01363                                 *level = 2;
01364                                 *Strong = TauLine2[i].Coll.heat;
01365                         }
01366                 }
01367         }
01368 
01369         /* now do co carbon monoxide lines */
01370         for( i=0; i < nCORotate; i++ )
01371         {
01372                 /* check if a line was the major heat agent */
01373                 if( C12O16Rotate[i].Coll.heat > *Strong )
01374                 {
01375                         *ipStrong = i;
01376                         *level = 3;
01377                         *Strong = C12O16Rotate[i].Coll.heat;
01378                 }
01379         }
01380         for( i=0; i < nCORotate; i++ )
01381         {
01382                 /* check if a line was the major heat agent */
01383                 if( C13O16Rotate[i].Coll.heat > *Strong )
01384                 {
01385                         *ipStrong = i;
01386                         *level = 4;
01387                         *Strong = C13O16Rotate[i].Coll.heat;
01388                 }
01389         }
01390 
01391         /* now do the hyperfine structure lines */
01392         for( i=0; i < nHFLines; i++ )
01393         {
01394                 /* check if a line was the major heat agent */
01395                 if( HFLines[i].Coll.heat > *Strong )
01396                 {
01397                         *ipStrong = i;
01398                         *level = 5;
01399                         *Strong = HFLines[i].Coll.heat;
01400                 }
01401         }
01402 
01403         /*Chianti & Leiden Atoms & Molecules */
01404         for( i=0; i <linesAdded2; i++)
01405         {
01406                 /* check if a line was the major heat agent */
01407                 if( atmolEmis[i].tran->Coll.heat > *Strong )
01408                 {
01409                         *ipStrong = i;
01410                         *level = 6;
01411                         *Strong = atmolEmis[i].tran->Coll.heat;
01412                 }
01413         }
01414         return;
01415 }
01416 
01417 quantumState *AddState2Stack( void )
01418 {
01419         DEBUG_ENTRY( "AddState2Stack()" );
01420 
01421         ASSERT( !lgStatesAdded );
01422 
01423         currentState = new quantumState;
01424 
01425         StateJunk( currentState );
01426 
01427         if( statesAdded == 0 )
01428         {
01429                 GenericStates = currentState;
01430                 GenericStates->next = NULL;
01431                 lastState = GenericStates;
01432         }
01433         else
01434         {
01435                 StateZero( currentState );
01436                 lastState->next = currentState;
01437                 lastState = lastState->next;
01438         }
01439 
01440         statesAdded++;
01441 
01442         return currentState;
01443 }
01444 
01445 emission *AddLine2Stack( bool lgRadiativeTrans )
01446 {
01447         DEBUG_ENTRY( "AddLine2Stack()" );
01448 
01449         if( !lgRadiativeTrans )
01450         {
01451                 return &DummyEmis;
01452         }
01453         else
01454         {
01455                 ASSERT( lgLinesAdded == false );
01456 
01457                 currentLine = new emission;
01458 
01459                 EmLineJunk( currentLine );
01460 
01461                 if( linesAdded == 0 )
01462                 {
01463                         GenericLines = currentLine;
01464                         GenericLines->next = NULL;
01465                         lastLine = GenericLines;
01466                 }
01467                 else
01468                 {
01469                         /* 
01470                         \todo 2 Does doing EmLineZero here defeat the purpose of EmLineJunk? 
01471                         * maybe we should pass full set of Emis components, fill everything in 
01472                         * here, and THEN use EmLineZero?  */
01473                         EmLineZero( currentLine );
01474 
01475                         lastLine->next = currentLine;
01476                         lastLine = lastLine->next;
01477                 }
01478 
01479                 linesAdded++;
01480                 return currentLine;
01481         }
01482 }
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

Generated for cloudy by doxygen 1.7.6.1