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