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 /*lines main routine to put emission line intensities into line stack, 00004 * calls lineset1, 2, 3, 4 */ 00005 /*Drive_cdLine do the drive cdLine command */ 00006 #include "cddefines.h" 00007 #include "physconst.h" 00008 #include "taulines.h" 00009 #include "thermal.h" 00010 #include "yield.h" 00011 #include "ionbal.h" 00012 #include "cddrive.h" 00013 #include "trace.h" 00014 #include "dense.h" 00015 #include "prt.h" 00016 #include "rt.h" 00017 #include "coolheavy.h" 00018 #include "rfield.h" 00019 #include "phycon.h" 00020 #include "elementnames.h" 00021 #include "iso.h" 00022 #include "hyperfine.h" 00023 #include "hydrogenic.h" 00024 #include "lines_service.h" 00025 #include "atmdat.h" 00026 #include "lines.h" 00027 #include "radius.h" 00028 STATIC void Drive_cdLine( void ); 00029 00030 void lines(void) 00031 { 00032 char chLabel[5]; 00033 static bool lgRecOn; 00034 long int i, 00035 ipnt, 00036 nelem; 00037 double BigstExtra, 00038 ExtraCool, 00039 f2, sum; 00040 00041 DEBUG_ENTRY( "lines()" ); 00042 00043 /* LineSave.ipass 00044 * -1 - space not yet allocated - just count number of lines entered into stack 00045 * 0 - first call with space allocated - must create labels and add in wavelengths 00046 * +1 - later calls - just add intensity 00047 */ 00048 00049 /* major routines used here: 00050 * 00051 * PutLine( tarray ) 00052 * this uses information in tarray to give various 00053 * contributions to lines, and their intensities 00054 * 00055 * PutExtra( extra ) 00056 * extra is some extra intensity to add to the line 00057 * it will go into the totl contribution put out by PutLine, 00058 * and this contribution should be indicated by independent 00059 * call to linadd 00060 * */ 00061 00062 if( trace.lgTrace ) 00063 { 00064 fprintf( ioQQQ, " lines called\n" ); 00065 } 00066 00067 /* the drive cdline command, which checks that all lines can be pulled out by cdLine */ 00068 if( trace.lgDrv_cdLine && LineSave.ipass > 0 ) 00069 Drive_cdLine(); 00070 00071 /* total luminosity radiated by this model - will be compared with energy in incident 00072 * continuum when calculation is complete */ 00073 thermal.power += thermal.htot*radius.dVeff; 00074 00075 /* remember the total free-free heating */ 00076 thermal.FreeFreeTotHeat += CoolHeavy.brems_heat_total*radius.dVeff; 00077 00078 /* total Compton heating - cooling */ 00079 rfield.comtot += rfield.cmheat*radius.dVeff; 00080 thermal.totcol += thermal.ctot*radius.dVeff; 00081 00082 /* up up induced recombination cooling */ 00083 for( nelem=0; nelem<LIMELM; ++nelem ) 00084 { 00085 hydro.cintot += iso.RecomInducCool_Rate[ipH_LIKE][nelem]*radius.dVeff; 00086 } 00087 00088 /* nsum is line pointer within large stack of line intensities */ 00089 LineSave.nsum = 0; 00090 LineSave.nComment = 0; 00091 00092 /* this is used by lindst to proportion inward and outward. should be 50% for 00093 * optically thin line. putline sets this to actual value for particular line 00094 * and calls lindst then rests to 50% */ 00095 rt.fracin = 0.5; 00096 00097 /* last arg in call to lindst and linadd is info on line 00098 * info is char variable indicating type of line this is 00099 * 'c' cooling 00100 * 'h' heating 00101 * 'i' information only 00102 * 'r' recombination line 00103 * 00104 * all components of lines are entered into the main line stack here 00105 * when printing, filters exist to not print Inwd component */ 00106 00107 /* initialize routine that can generate pointers for forbidden lines, 00108 * these are lines that are not transferred otherwise, 00109 * in following routines there will be pairs of calls, first to 00110 * PntForLine to get pointer, then lindst to add to stack */ 00111 PntForLine(0.,"FILL",&i); 00112 00113 /* evaluate rec coefficient for rec lines of C, N, O 00114 * some will be used in LineSet2 and then zeroed out, 00115 * others left alone and used below */ 00116 t_ADfA::Inst().rec_lines(phycon.te,LineSave.RecCoefCNO); 00117 00118 /* initialize ExtraInten with a zero */ 00119 PutExtra(0.); 00120 00121 /* put in something impossible in element 0 of line stack */ 00122 linadd(0.f,0,"zero",'i' , "null placeholder"); 00123 00124 /* this is compared with true volume in final. The number can't 00125 * actually be unity since this would overflow on a 32 bit machine */ 00126 /* integrate the volume as a sanity check */ 00127 linadd( 1.e-10 , 1 , "Unit" , 'i' , "unit integration placeholder"); 00128 00129 /* initial set of general properties */ 00130 lines_general(); 00131 00132 /* do all continua */ 00133 lines_continuum(); 00134 00135 /* information about grains */ 00136 lines_grains(); 00137 00138 /* update all satellite lines */ 00139 iso_satellite_update(); 00140 00141 /* do all hydrogenic ions */ 00142 lines_hydro(); 00143 00144 /* enter He-iso sequence lines */ 00145 lines_helium(); 00146 00147 #if 0 00148 /* This is Ryan's code for dumping lots of Helium lines according to 00149 * quantum number rather than wavelength, principally for comparison with Rob 00150 * Bauman's code. */ 00151 if( iteration > 1 ) 00152 { 00153 fprintf( ioQQQ,"ipHi\tipLo\tnu\tlu\tsu\tnl\tll\tsl\tWL\tintens\n" ); 00154 for( long ipHi=5; ipHi<= iso.numLevels_local[ipHE_LIKE][ipHELIUM] - iso.nCollapsed_local[ipHE_LIKE][ipHELIUM]; ipHi++ ) 00155 { 00156 for( long ipLo=0; ipLo<ipHi; ipLo++ ) 00157 { 00158 if( Transitions[ipHE_LIKE][ipHELIUM][ipHi][ipLo].ipCont > 0 ) 00159 { 00160 double relint, absint; 00161 00162 if( cdLine("He 1", 00163 Transitions[ipHE_LIKE][ipHELIUM][ipHi][ipLo].WLAng, 00164 &relint, &absint ) ) 00165 { 00166 //Transitions[ipHE_LIKE][ipHELIUM][ipHi][ipLo].Hi->chLabel 00167 00168 if( Transitions[ipHE_LIKE][ipHELIUM][ipHi][ipLo].WLAng < 1.1E4 && 00169 Transitions[ipHE_LIKE][ipHELIUM][ipHi][ipLo].WLAng > 3.59E3 && 00170 ipLo!=3 && ipLo!=4 && relint >= 0.0009 ) 00171 { 00172 fprintf( ioQQQ,"%li\t%li\t%li\t%li\t%li\t%li\t%li\t%li\t%e\t%e\n", 00173 ipHi, 00174 ipLo, 00175 StatesElem[ipHE_LIKE][ipHELIUM][ipHi].n, 00176 StatesElem[ipHE_LIKE][ipHELIUM][ipHi].l, 00177 StatesElem[ipHE_LIKE][ipHELIUM][ipHi].S, 00178 StatesElem[ipHE_LIKE][ipHELIUM][ipLo].n, 00179 StatesElem[ipHE_LIKE][ipHELIUM][ipLo].l, 00180 StatesElem[ipHE_LIKE][ipHELIUM][ipLo].S, 00181 Transitions[ipHE_LIKE][ipHELIUM][ipHi][ipLo].WLAng, 00182 relint ); 00183 } 00184 } 00185 } 00186 } 00187 } 00188 } 00189 #endif 00190 00191 /* do heavies, lithium through neon */ 00192 lines_lv1_li_ne(); 00193 00194 /* do heavies, sodium through argon */ 00195 lines_lv1_na_ar(); 00196 00197 /* do heavies, potassium through zinc */ 00198 lines_lv1_k_zn(); 00199 00200 /* add up line intensities for certain set of lines */ 00201 sum = PrtLineSum( " SUM" ); 00202 /* zero out the location that will receive this information, 00203 * remember that memory not allocated until ipass >= 0 */ 00204 if( LineSave.ipass > 0 ) 00205 LineSv[LineSave.nsum].sumlin[LineSave.lgLineEmergent] = 0.; 00206 /* optional sum of certain emission lines, set with "print sum" */ 00207 linadd(sum/radius.dVeff,0,"Stoy",'i' , 00208 "Stoy method energy sum "); 00209 00210 /* next come some recombination lines */ 00211 i = StuffComment( "recombination" ); 00212 linadd( 0., (realnum)i , "####", 'i' , 00213 "recombination lines"); 00214 00215 /*********************************************************************** 00216 * large number of C, N, and O recombination lines * 00217 *************************************************************************/ 00218 00219 if( LineSave.ipass <= 0 ) 00220 { 00221 /* initialize to true, may be set false if density is very high below */ 00222 lgRecOn = true; 00223 } 00224 00225 for( i=0; i < 471; i++ ) 00226 { 00227 /* generate label for the line */ 00228 strcpy( chLabel, elementnames.chElementSym[(long)(LineSave.RecCoefCNO[0][i])-1] ); 00229 strcat( chLabel, elementnames.chIonStage[(long)(LineSave.RecCoefCNO[0][i]- 00230 LineSave.RecCoefCNO[1][i]+1.01)-1] ); 00231 00232 /* number of rec per unit vol 00233 * >>chng 97 aug 28, do not predict rec intensities at high densities 00234 * blr.in and oldblr.in go nuts if this is added in outward beam 00235 * these recombination coefficients were computed in the low density limit - 00236 * they were not intended for high densities */ 00237 if( dense.eden < 1e8 && lgRecOn ) 00238 { 00239 f2 = LineSave.RecCoefCNO[3][i]*dense.eden* 00240 dense.xIonDense[(long)(LineSave.RecCoefCNO[0][i])-1][(long)(LineSave.RecCoefCNO[0][i]-LineSave.RecCoefCNO[1][i]+2)-1]; 00241 00242 /* convert to intensity */ 00243 f2 = f2*1.99e-8/LineSave.RecCoefCNO[2][i]; 00244 } 00245 else 00246 { 00247 lgRecOn = false; 00248 f2 = 0.; 00249 } 00250 /* finally stuff it into the stack */ 00251 PntForLine(LineSave.RecCoefCNO[2][i],chLabel,&ipnt); 00252 lindst(f2,LineSave.RecCoefCNO[2][i],chLabel,ipnt, 'r',true , 00253 "recombination line"); 00254 } 00255 00256 /* next come the atom_level2 lines */ 00257 i = StuffComment( "level2 lines" ); 00258 linadd( 0., (realnum)i , "####", 'i' , 00259 "level2 lines"); 00260 00261 /* add in all the other level 2 wind lines 00262 * Dima's 6k lines */ 00263 ExtraCool = 0.; 00264 BigstExtra = 0.; 00265 for( i=0; i < nWindLine; i++ ) 00266 { 00267 if( TauLine2[i].Hi->IonStg < TauLine2[i].Hi->nelem+1-NISO ) 00268 { 00269 PutLine(&TauLine2[i],"level 2 line"); 00270 if( TauLine2[i].Coll.cool > BigstExtra ) 00271 { 00272 BigstExtra = TauLine2[i].Coll.cool; 00273 thermal.ipMaxExtra = i+1; 00274 } 00275 ExtraCool += TauLine2[i].Coll.cool; 00276 } 00277 } 00278 /* keep track of how important this is */ 00279 thermal.GBarMax = MAX2(thermal.GBarMax,(realnum)(ExtraCool/thermal.ctot)); 00280 00281 /* next come the hyperfine structure lines */ 00282 i = StuffComment( "hyperfine structure" ); 00283 linadd( 0., (realnum)i , "####", 'i' , 00284 "hyperfine structure lines "); 00285 00286 /* this is total cooling due to all HF lines */ 00287 linadd( hyperfine.cooling_total, 0., "hfin", 'i' , 00288 "total cooling all hyperfine structure lines "); 00289 00290 /* remember largest local cooling for possible printout in comments */ 00291 hyperfine.cooling_max = (realnum)MAX2(hyperfine.cooling_max,hyperfine.cooling_total/thermal.ctot); 00292 00293 /* the hyperfine lines */ 00294 for( i=0; i < nHFLines; i++ ) 00295 { 00296 PutLine(&HFLines[i], 00297 "hyperfine structure line"); 00298 } 00299 00300 /* Databases: Atoms & Molecules: Added Oct 26 07*/ 00301 i = StuffComment( "database lines" ); 00302 linadd( 0., (realnum)i , "####", 'i' , 00303 "database lines "); 00304 00305 /* Atoms & Molecules */ 00306 for( i=0; i < linesAdded2; i++ ) 00307 { 00308 /* \todo 2 say which database in the comment */ 00309 PutLine(atmolEmis[i].tran, "lines from third-party databases", atmolEmis[i].tran->Hi->chLabel); 00310 } 00311 00312 /* next come the inner shell fluorescent lines */ 00313 i = StuffComment( "inner shell" ); 00314 linadd( 0., (realnum)i , "####", 'i' , 00315 "inner shell lines"); 00316 00317 /* the group of inner shell fluorescent lines */ 00318 for( i=0; i < t_yield::Inst().nlines(); ++i ) 00319 { 00320 double xInten = 00321 /* density of parent ion, cm-3 */ 00322 dense.xIonDense[t_yield::Inst().nelem(i)][t_yield::Inst().ion(i)] * 00323 /* photo rate per atom per second, s-1 */ 00324 ionbal.PhotoRate_Shell[t_yield::Inst().nelem(i)][t_yield::Inst().ion(i)][t_yield::Inst().nshell(i)][0]* 00325 /* fluor yield - dimensionless */ 00326 t_yield::Inst().yield(i) * 00327 /* photon energy - ryd, converted into ergs */ 00328 t_yield::Inst().energy(i) * EN1RYD; 00329 00330 /* create label if initializing line stack */ 00331 if( LineSave.ipass == 0 ) 00332 { 00333 /* only generate the line label if it is going to be used */ 00334 strcpy( chLabel , elementnames.chElementSym[t_yield::Inst().nelem(i)] ); 00335 strcat( chLabel , elementnames.chIonStage[t_yield::Inst().ion_emit(i)] ); 00336 # if 0 00337 /* only print yields for atoms */ 00338 if( t_yield::Inst().ion(i) == 0 && t_yield::Inst().nelem(i) == ipIRON ) 00339 fprintf(ioQQQ,"DEBUGyeild\t%s\t%.3f\t%.3e\n", 00340 /* line designation, energy in eV, yield */ 00341 chLabel , t_yield::Inst().energy(i)*EVRYD, t_yield::Inst().yield(i) ); 00342 # endif 00343 } 00344 00345 /* the group of inner shell fluorescent lines */ 00346 lindst( 00347 /* intensity of line */ 00348 xInten, 00349 /* wavelength of line in Angstroms */ 00350 (realnum)RYDLAM / t_yield::Inst().energy(i), 00351 /* label */ 00352 chLabel , 00353 /* continuum array offset for line as set in ipoint */ 00354 t_yield::Inst().ipoint(i), 00355 /* type of line - count as a recombination line */ 00356 'r', 00357 /* include line in continuum? */ 00358 true , 00359 "inner shell line"); 00360 } 00361 00362 /* >>chng 06 jan 03, confirm that number of lines never changes once we have 00363 * created the labels */ 00364 { 00365 static long nLineSave=-1 , ndLineSave=-1; 00366 if( LineSave.ipass == 0 ) 00367 { 00368 nLineSave = LineSave.nsum; 00369 ndLineSave = LineSave.nsum; 00370 } 00371 else if( LineSave.ipass > 0 ) 00372 { 00373 /* this can't happen */ 00374 if( nLineSave<= 0 || ndLineSave < 0 ) 00375 TotalInsanity(); 00376 00377 /* now make sure that we have the same number of lines as we had previously 00378 * created labels. This would not pass if we did not add in exactly the same 00379 * number of lines on each pass */ 00380 if( nLineSave != LineSave.nsum ) 00381 { 00382 fprintf( ioQQQ, "DISASTER number of lines in LineSave.nsum changed between pass 0 and 1 - this is impossible\n" ); 00383 fprintf( ioQQQ, "DISASTER LineSave.nsum is %li and nLineSave is %li\n", 00384 LineSave.nsum , 00385 nLineSave); 00386 ShowMe(); 00387 cdEXIT(EXIT_FAILURE); 00388 } 00389 if( ndLineSave != LineSave.nsum ) 00390 { 00391 fprintf( ioQQQ, "DISASTER number of lines in LineSave.nsum changed between pass 0 and 1 - this is impossible\n" ); 00392 fprintf( ioQQQ, "DISASTER LineSave.nsum is %li and ndLineSave is %li\n", 00393 LineSave.nsum , 00394 ndLineSave); 00395 ShowMe(); 00396 cdEXIT(EXIT_FAILURE); 00397 } 00398 } 00399 } 00400 00401 /* now do all molecules - do last since so many H2 lines */ 00402 lines_molecules(); 00403 00404 if( trace.lgTrace ) 00405 { 00406 fprintf( ioQQQ, " lines returns\n" ); 00407 } 00408 return; 00409 } 00410 00411 /*Drive_cdLine do the drive cdLine command */ 00412 STATIC void Drive_cdLine( void ) 00413 { 00414 long int j; 00415 bool lgMustPrintHeader = true; 00416 double absval , rel; 00417 00418 DEBUG_ENTRY( "Drive_cdLine()" ); 00419 00420 for( j=1; j < LineSave.nsum; j++ ) 00421 { 00422 if( cdLine( LineSv[j].chALab , LineSv[j].wavelength , &absval , &rel ) <= 0 ) 00423 { 00424 /* print header if first miss */ 00425 if( lgMustPrintHeader ) 00426 { 00427 fprintf(ioQQQ,"n\tlab\twl\n"); 00428 lgMustPrintHeader = false; 00429 } 00430 00431 fprintf(ioQQQ,"%li\t%s\t%f\n", j, LineSv[j].chALab , LineSv[j].wavelength ); 00432 } 00433 } 00434 fprintf( ioQQQ, " Thanks for checking on the cdLine routine!\n" ); 00435 cdEXIT(EXIT_FAILURE); 00436 }