cloudy trunk
|
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 }