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