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 /*FeII_Colden maintain H2 column densities within X */ 00004 /*FeIILevelPops main FeII routine, called by CoolIron to evaluate iron cooling */ 00005 /*FeIICreate read in needed data from file 00006 * convert form of FeII data, as read in from file within routine FeIICreate 00007 * into physical form. called by atmdat_readin */ 00008 /*FeIIPunPop - punch level populations */ 00009 /*AssertFeIIDep called by assert FeII depart coef command */ 00010 /*FeIIPrint print FeII information */ 00011 /*FeIICollRatesBoltzmann evaluate collision strenths, 00012 * both interpolating on r-mat and creating g-bar 00013 * find Boltzmann factors, evaluate collisional rate coefficients */ 00014 /*FeIIPrint print output from large FeII atom, called by prtzone */ 00015 /*FeIISumBand sum up large FeII emission over certain bands, called in lineset4 */ 00016 /*FeII_RT_TauInc called once per zone in RT_tau_inc to increment large FeII atom line optical depths */ 00017 /*FeII_RT_tau_reset reset optical depths for large FeII atom, called by update after each iteration */ 00018 /*FeIIPoint called by ContCreatePointers to create pointers for lines in large FeII atom */ 00019 /*FeIIAccel called by rt_line_driving to compute radiative acceleration due to FeII lines */ 00020 /*FeIIAddLines save accumulated FeII intensities called by lineset4 */ 00021 /*FeIIPunchLines punch FeII lines at end of calculation, if punch verner set, called by punch_do*/ 00022 /*FeIIPunchOpticalDepth punch FeII line optical depth at end of calculation, called by punch_do*/ 00023 /*FeIIPunchLevels punch FeII levels and energies */ 00024 /*FeII_LineZero zero out storage for large FeII atom, called by tauout */ 00025 /*FeIIIntenZero zero out intensity of FeII atom */ 00026 /*FeIIZero initialize some variables, called by zero one time before commands parsed */ 00027 /*FeIIReset reset some variables, called by zero */ 00028 /*FeIIPunData punch line data */ 00029 /*FeIIPunDepart punch some departure coef for large atom, 00030 * set with punch FeII departure command*/ 00031 /*FeIIPun1Depart send the departure coef for physical level nPUN to unit ioPUN */ 00032 /*FeIIContCreate create FeII continuum bins to add lines into ncell cells between wavelengths lambda low and high, 00033 * returns number of cells used */ 00034 /*FeIIBandsCreate returns number of FeII bands */ 00035 /*FeII_RT_Out do outward rates for FeII, called by RT_diffuse */ 00036 /*FeII_OTS do ots rates for FeII, called by RT_OTS */ 00037 /*FeII_RT_Make called by RT_line_all, does large FeII atom radiative transfer */ 00038 /*FeIILyaPump find rate of Lya excitation of the FeII atom */ 00039 /*FeIIOvrLap handle overlapping FeII lines */ 00040 /*ParseAtomFeII parse the atom FeII command */ 00041 /*FeIIPunchLineStuff include FeII lines in punched optical depths, etc, called from PunchLineStuff */ 00042 #include "cddefines.h" 00043 #include "cddrive.h" 00044 #include "thermal.h" 00045 #include "physconst.h" 00046 #include "doppvel.h" 00047 #include "taulines.h" 00048 #include "dense.h" 00049 #include "rfield.h" 00050 #include "radius.h" 00051 #include "lines_service.h" 00052 #include "ipoint.h" 00053 #include "thirdparty.h" 00054 #include "hydrogenic.h" 00055 #include "lines.h" 00056 #include "rt.h" 00057 #include "trace.h" 00058 #include "punch.h" 00059 #include "phycon.h" 00060 #include "atomfeii.h" 00061 #include "iso.h" 00062 #include "pressure.h" 00063 00064 /* FeIIOvrLap handle overlapping FeII lines */ 00065 STATIC void FeIIOvrLap(void); 00066 00067 /* add FeII lines into ncell cells between wavelengths lambda low and high */ 00068 STATIC void FeIIContCreate(double xLamLow , double xLamHigh , long int ncell ); 00069 00070 /* read in the FeII Bands file, and sets nFeIIBands, the number of bands, 00071 * if argument is "" then use default file with bands, 00072 * if filename within "" then use it instead, 00073 * return value is 0 if success, 1 if failed */ 00074 STATIC int FeIIBandsCreate( const char chFile[] ); 00075 00076 /* this will be the smallest collision strength we will permit with the gbar 00077 * collision strengths, or for the data that have zeroes */ 00078 /* >>chng 99 jul 24, this was 1e-9 in the old fortran version and in baldwin et al. 96, 00079 * and has been changed several times, and affects results. this is the smallest 00080 * non-zero collision strength in the r-matrix calculations */ 00081 realnum CS2SMALL = (realnum)1e-5; 00082 /* routines used only within this file */ 00083 00084 /*FeIICollRatesBoltzmann evaluate collision strenths, 00085 * both interpolating on r-mat and creating g-bar 00086 * find Boltzmann factors, evaluate collisional rate coefficients */ 00087 STATIC void FeIICollRatesBoltzmann(void); 00088 /* find rate of Lya excitation of the FeII atom */ 00089 STATIC void FeIILyaPump(void); 00090 00091 /*extern realnum Fe2LevN[NFE2LEVN][NFE2LEVN][NTA];*/ 00092 /*extern realnum Fe2LevN[ipHi][ipLo].t[NTA];*/ 00093 /*realnum ***Fe2LevN;*/ 00094 /* >>chng 06 mar 01, boost to global namespace */ 00095 /*transition **Fe2LevN;*/ 00096 /* flag for the collision strength */ 00097 int **ncs1; 00098 00099 /* all following variables have file scope 00100 #define NFEIILINES 68635 */ 00101 00102 /* this is size of nnPradDat array */ 00103 #define NPRADDAT 159 00104 00105 /* band wavelength, lower and upper bounds, in vacuum Angstroms */ 00106 /* FeII_Bands[n][3], where n is the number of bands in fe2Bands.dat 00107 * these bands are defined in fe2Bands.dat and read in at startup 00108 * of calculation */ 00109 realnum **FeII_Bands; 00110 00111 /* continuum wavelengths, lower and upper bounds, in vacuum Angstroms, 00112 * third is integrated intensity */ 00113 /* FeII_Cont[n][3], where n is the number of cells needed 00114 * these bands are defined in FeIIContCreate */ 00115 realnum **FeII_Cont; 00116 00117 /* this is the number of bands read in from bands_Fe2.dat */ 00118 long int nFeIIBands; 00119 00120 /* number of bands in continuum array */ 00121 long int nFeIIConBins; 00122 00123 /* the dim of this vector this needs one extra since there are 20 numbers per line, 00124 * highest not ever used for anything */ 00125 /*long int nnPradDat[NPRADDAT+1];*/ 00126 static long int *nnPradDat; 00127 00128 /* malloced in feiidata */ 00129 /* realnum sPradDat[NPRADDAT][NPRADDAT][8];*/ 00130 /* realnum sPradDat[ipHi][ipLo].t[8];*/ 00131 static realnum ***sPradDat; 00132 00133 /* array used to integrate FeII line intensities over model, 00134 * Fe2SavN[upper][lower], 00135 *static realnum Fe2SavN[NFE2LEVN][NFE2LEVN];*/ 00136 static double **Fe2SavN; 00137 00138 /* save effective transition rates */ 00139 static double **Fe2A; 00140 00141 /* induced transition rates */ 00142 static double **Fe2LPump, **Fe2CPump; 00143 00144 /* energies read in from fe2energies.dat data file */ 00145 realnum *Fe2Energies; 00146 00147 /* collision rates */ 00148 static realnum **Fe2Coll; 00149 00150 /* Fe2DepCoef[NFE2LEVN]; 00151 realnum cli[NFEIILINES], 00152 cfe[NFEIILINES];*/ 00153 static double 00154 /* departure coefficients */ 00155 *Fe2DepCoef , 00156 /* level populations */ 00157 *Fe2LevelPop , 00158 /* column densities */ 00159 *Fe2ColDen , 00160 /* this will become array of Boltzmann factors */ 00161 *FeIIBoltzmann; 00162 /*FeIIBoltzmann[NFE2LEVN] ,*/ 00163 00164 static double EnerLyaProf1, 00165 EnerLyaProf4, 00166 PhotOccNumLyaCenter; 00167 static double 00168 /* the yVector - will become level populations after matrix inversion */ 00169 *yVector, 00170 /* this is used to call matrix routines */ 00171 /*xMatrix[NFE2LEVN][NFE2LEVN] ,*/ 00172 **xMatrix , 00173 /* this will become the very large 1-D array that 00174 * is passed to the matrix inversion routine*/ 00175 *amat; 00176 00177 00178 /*FeII_Colden maintain H2 column densities within X */ 00179 void FeII_Colden( const char *chLabel ) 00180 { 00181 long int n; 00182 00183 DEBUG_ENTRY( "FeII_Colden()" ); 00184 00185 /* >>chng 05 nov 29, FeII always on, always want to evaluate this */ 00186 /* nothing to do if no FeII 00187 if( !FeII. lgFeIION ) 00188 return;*/ 00189 00190 if( strcmp(chLabel,"ZERO") == 0 ) 00191 { 00192 /* zero out column density */ 00193 for( n=0; n < FeII.nFeIILevel; ++n ) 00194 { 00195 /* space for the rotation quantum number */ 00196 Fe2ColDen[n] = 0.; 00197 } 00198 } 00199 00200 else if( strcmp(chLabel,"ADD ") == 0 ) 00201 { 00202 /* add together column densities */ 00203 for( n=0; n < FeII.nFeIILevel; ++n ) 00204 { 00205 /* state-specific FeII column density */ 00206 Fe2ColDen[n] += Fe2LevelPop[n]*radius.drad_x_fillfac; 00207 } 00208 } 00209 00210 /* check for the print option, which we can't do, else we have a problem */ 00211 else if( strcmp(chLabel,"PRIN") != 0 ) 00212 { 00213 fprintf( ioQQQ, " FeII_Colden does not understand the label %s\n", 00214 chLabel ); 00215 cdEXIT(EXIT_FAILURE); 00216 } 00217 00218 return; 00219 } 00220 00221 /* 00222 *===================================================================== 00223 */ 00224 /* FeIICreate read in FeII data from files on disk. called by atmdat_readin 00225 * but only if FeII. lgFeIION is true, set with atom FeII verner command */ 00226 void FeIICreate(void) 00227 { 00228 FILE *ioDATA; 00229 00230 char chLine[FILENAME_PATH_LENGTH_2]; 00231 00232 long int i, 00233 ipHi , 00234 ipLo, 00235 lo, 00236 ihi, 00237 k, 00238 m1, 00239 m2, 00240 m3; 00241 00242 DEBUG_ENTRY( "FeIICreate()" ); 00243 00244 if( lgFeIIMalloc ) 00245 { 00246 /* we have already been called one time, just bail out */ 00247 00248 return; 00249 } 00250 00251 /* now set flag so never done again - this is set false when initi 00252 * when this is true it is no longer possible to change the number of levels 00253 * in the model atom with the atom FeII levels command */ 00254 lgFeIIMalloc = true; 00255 00256 /* remember how many levels this was, so that in future calculations 00257 * we can reset the atom to full value */ 00258 FeII.nFeIILevelAlloc = FeII.nFeIILevel; 00259 00260 /* set up array to save FeII transition probabilities */ 00261 Fe2A = (double **)MALLOC(sizeof(double *)*(unsigned long)FeII.nFeIILevelAlloc ); 00262 00263 /* second dimension, lower level, for line save array */ 00264 for( ipHi=0; ipHi < FeII.nFeIILevelAlloc; ++ipHi ) 00265 { 00266 Fe2A[ipHi]=(double*)MALLOC(sizeof(double )*(unsigned long)FeII.nFeIILevelAlloc); 00267 } 00268 00269 /* set up array to save FeII pumping rates */ 00270 Fe2CPump = (double **)MALLOC(sizeof(double *)*(unsigned long)FeII.nFeIILevelAlloc ); 00271 00272 /* set up array to save FeII pumping rates */ 00273 Fe2LPump = (double **)MALLOC(sizeof(double *)*(unsigned long)FeII.nFeIILevelAlloc ); 00274 00275 /* second dimension, lower level, for line save array */ 00276 for( ipHi=0; ipHi < FeII.nFeIILevelAlloc; ++ipHi ) 00277 { 00278 Fe2CPump[ipHi]=(double*)MALLOC(sizeof(double )*(unsigned long)FeII.nFeIILevelAlloc); 00279 00280 Fe2LPump[ipHi]=(double*)MALLOC(sizeof(double )*(unsigned long)FeII.nFeIILevelAlloc); 00281 } 00282 00283 /* set up array to save FeII collision rates */ 00284 Fe2Energies = (realnum *)MALLOC(sizeof(realnum)*(unsigned long)FeII.nFeIILevelAlloc ); 00285 00286 /* set up array to save FeII collision rates */ 00287 Fe2Coll = (realnum **)MALLOC(sizeof(realnum *)*(unsigned long)FeII.nFeIILevelAlloc ); 00288 00289 /* second dimension, lower level, for line save array */ 00290 for( ipHi=0; ipHi < FeII.nFeIILevelAlloc; ++ipHi ) 00291 { 00292 Fe2Coll[ipHi]=(realnum*)MALLOC(sizeof(realnum )*(unsigned long)FeII.nFeIILevelAlloc); 00293 } 00294 00295 /* MALLOC space for the 2D matrix array */ 00296 xMatrix = (double **)MALLOC(sizeof(double *)*(unsigned long)FeII.nFeIILevelAlloc ); 00297 00298 /* now do the second dimension */ 00299 for( i=0; i<FeII.nFeIILevelAlloc; ++i ) 00300 { 00301 xMatrix[i] = (double *)MALLOC(sizeof(double)*(unsigned long)FeII.nFeIILevelAlloc ); 00302 } 00303 /* MALLOC space for the 1-yVector array */ 00304 amat=(double*)MALLOC( (sizeof(double)*(unsigned long)(FeII.nFeIILevelAlloc*FeII.nFeIILevelAlloc) )); 00305 00306 /* MALLOC space for the 1-yVector array */ 00307 yVector=(double*)MALLOC( (sizeof(double)*(unsigned long)(FeII.nFeIILevelAlloc) )); 00308 00309 /* set up array to save FeII line intensities */ 00310 Fe2SavN = (double **)MALLOC(sizeof(double *)*(unsigned long)FeII.nFeIILevelAlloc ); 00311 00312 /* second dimension, lower level, for line save array */ 00313 for( ipHi=1; ipHi < FeII.nFeIILevelAlloc; ++ipHi ) 00314 { 00315 Fe2SavN[ipHi]=(double*)MALLOC(sizeof(double )*(unsigned long)ipHi); 00316 } 00317 00318 /* now MALLOC space for energy level table*/ 00319 nnPradDat = (long*)MALLOC( (NPRADDAT+1)*sizeof(long) ); 00320 00321 /*Fe2DepCoef[NFE2LEVN];*/ 00322 Fe2DepCoef = (double*)MALLOC( (unsigned long)FeII.nFeIILevelAlloc*sizeof(double) ); 00323 00324 /*Fe2LevelPop[NFE2LEVN];*/ 00325 Fe2LevelPop = (double*)MALLOC( (unsigned long)FeII.nFeIILevelAlloc*sizeof(double) ); 00326 00327 /*Fe2ColDen[NFE2LEVN];*/ 00328 Fe2ColDen = (double*)MALLOC( (unsigned long)FeII.nFeIILevelAlloc*sizeof(double) ); 00329 00330 /*FeIIBoltzmann[NFE2LEVN];*/ 00331 FeIIBoltzmann = (double*)MALLOC( (unsigned long)FeII.nFeIILevelAlloc*sizeof(double) ); 00332 00333 00334 /* MALLOC the realnum sPradDat[NPRADDAT][NPRADDAT][8] array */ 00335 /* MALLOC the realnum sPradDat[ipHi][ipLo].t[8] array */ 00336 sPradDat = ((realnum ***)MALLOC(NPRADDAT*sizeof(realnum **))); 00337 00338 /* >>chng 00 dec 06, changed lower limit of loop to 1, Tru64 alpha's will not allocate 0 bytes!, PvH */ 00339 sPradDat[0] = NULL; 00340 for(ipHi=1; ipHi < NPRADDAT; ipHi++) 00341 { 00342 /* >>chng 00 dec 06, changed sizeof(realnum) to sizeof(realnum*), PvH */ 00343 sPradDat[ipHi] = (realnum **)MALLOC((unsigned long)ipHi*sizeof(realnum *)); 00344 00345 /* now make space for the second dimension */ 00346 for( ipLo=0; ipLo< ipHi; ipLo++ ) 00347 { 00348 sPradDat[ipHi][ipLo] = (realnum *)MALLOC(8*sizeof(realnum )); 00349 } 00350 } 00351 00352 /* now set junk to make sure reset before used */ 00353 for(ipHi=0; ipHi < NPRADDAT; ipHi++) 00354 { 00355 for( ipLo=0; ipLo< ipHi; ipLo++ ) 00356 { 00357 for( k=0; k<8; ++k ) 00358 { 00359 sPradDat[ipHi][ipLo][k] = -FLT_MAX; 00360 } 00361 } 00362 } 00363 00364 /* now create the main emission line array and a helper array for the cs flag */ 00365 Fe2LevN=(transition**)MALLOC(sizeof(transition *)*(unsigned long)FeII.nFeIILevelAlloc); 00366 ncs1=(int**)MALLOC(sizeof(int *)*(unsigned long)FeII.nFeIILevelAlloc); 00367 00368 for( ipHi=1; ipHi < FeII.nFeIILevelAlloc; ++ipHi ) 00369 { 00370 Fe2LevN[ipHi]=(transition*)MALLOC(sizeof(transition)*(unsigned long)ipHi); 00371 00372 ncs1[ipHi]=(int*)MALLOC(sizeof(int)*(unsigned long)ipHi); 00373 } 00374 00375 #if 0 00376 /* now that its created, set to junk */ 00377 for( ipLo=0; ipLo < (FeII.nFeIILevelAlloc - 1); ipLo++ ) 00378 { 00379 for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevelAlloc; ipHi++ ) 00380 { 00381 TransitionJunk( &Fe2LevN[ipHi][ipLo] ); 00382 } 00383 } 00384 #endif 00385 00386 /* now assign state and Emis pointers */ 00387 Fe2LevN[1][0].Lo = AddState2Stack(); 00388 /* now that its created, set to junk */ 00389 for( ipHi=1; ipHi < FeII.nFeIILevelAlloc; ipHi++ ) 00390 { 00391 /* add this upper level */ 00392 Fe2LevN[ipHi][0].Hi = AddState2Stack(); 00393 for( ipLo=0; ipLo < ipHi; ipLo++ ) 00394 { 00395 if( ipLo == 0 ) 00396 { 00397 Fe2LevN[ipHi][ipLo].Lo = Fe2LevN[1][0].Lo; 00398 } 00399 else 00400 { 00401 /* lower state is same as a previous upper state. */ 00402 Fe2LevN[ipHi][ipLo].Lo = Fe2LevN[ipLo][0].Hi; 00403 } 00404 00405 Fe2LevN[ipHi][ipLo].Hi = Fe2LevN[ipHi][0].Hi; 00406 Fe2LevN[ipHi][ipLo].Emis = AddLine2Stack( true ); 00407 } 00408 } 00409 00410 if( trace.lgTrace ) 00411 { 00412 fprintf( ioQQQ," FeIICreate opening fe2nn.dat:"); 00413 } 00414 00415 ioDATA = open_data( "fe2nn.dat", "r" ); 00416 00417 ASSERT( ioDATA !=NULL ); 00418 /* read in the fe2nn.dat file - this gives the Zheng and Pradhan number of level 00419 * for every cloudy level. So nnPradDat[1] is the index in the cloudy level 00420 * counting for level 1 of Zheng & Pradan 00421 * note that the order of some levels is different, the nnPradDat file goes 00422 * 21 23 22 - also that many levels are missing, the file goes 95 99 94 93 116 00423 */ 00424 for( i=0; i < 8; i++ ) 00425 { 00426 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 00427 { 00428 fprintf( ioQQQ, " fe2nn.dat error reading data\n" ); 00429 cdEXIT(EXIT_FAILURE); 00430 } 00431 00432 /* get these integers from fe2nn.dat */ 00433 k = 20*i; 00434 /* NPRADDAT is size of nnPradDat array, 159+1, make sure we do not exceed it */ 00435 ASSERT( k+19 < NPRADDAT+1 ); 00436 sscanf( chLine , 00437 "%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld", 00438 &nnPradDat[k+0], &nnPradDat[k+1], &nnPradDat[k+2], &nnPradDat[k+3], &nnPradDat[k+4], 00439 &nnPradDat[k+5], &nnPradDat[k+6], &nnPradDat[k+7], &nnPradDat[k+8], &nnPradDat[k+9], 00440 &nnPradDat[k+10],&nnPradDat[k+11], &nnPradDat[k+12],&nnPradDat[k+13],&nnPradDat[k+14], 00441 &nnPradDat[k+15],&nnPradDat[k+16], &nnPradDat[k+17],&nnPradDat[k+18],&nnPradDat[k+19] 00442 ); 00443 # if !defined(NDEBUG) 00444 for( m1=0; m1<20; ++m1 ) 00445 { 00446 ASSERT( nnPradDat[k+m1] >= 0 && nnPradDat[k+m1] <= NFE2LEVN ); 00447 } 00448 # endif 00449 } 00450 fclose(ioDATA); 00451 00452 /* now get energies */ 00453 if( trace.lgTrace ) 00454 { 00455 fprintf( ioQQQ," FeIICreate opening fe2energies.dat:"); 00456 } 00457 00458 ioDATA = open_data( "fe2energies.dat", "r" ); 00459 00460 /* file now open, read the data */ 00461 for( ipHi=0; ipHi < FeII.nFeIILevelAlloc; ipHi++ ) 00462 { 00463 /* keep reading until have non-comment line, one that does not 00464 * start with # */ 00465 chLine[0] = '#'; 00466 while( chLine[0] == '#' ) 00467 { 00468 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 00469 { 00470 fprintf( ioQQQ, " fe2energies.dat error reading data\n" ); 00471 cdEXIT(EXIT_FAILURE); 00472 } 00473 } 00474 00475 /* first and last number on this line */ 00476 double help; 00477 sscanf( chLine, "%lf", &help ); 00478 Fe2Energies[ipHi] = (realnum)help; 00479 } 00480 fclose(ioDATA); 00481 00482 /* now get radiation data. 00483 * These are from the following data sources: 00484 >>refer fe2 as Nahar, S., 1995, A&A 293, 967 00485 >>refer fe2 as Quinet, P., LeDourneuf, M., & Zeipppen C.J., 1996, A&AS, 120, 361 00486 >>refer fe2 as Furh, J.R., Martin, G.A., & Wiese, W.L., 1988; J Phys Chem Ref Data 17, Suppl 4 00487 >>refer fe2 as Giridhar, S., & Arellano Ferro, A., 1995; Ref Mexicana Astron Astrofis 31, 23 00488 >>refer fe2 as Kurucz, R.L., 1995, SAO CD ROM 23 00489 */ 00490 00491 if( trace.lgTrace ) 00492 { 00493 fprintf( ioQQQ," FeIICreate opening fe2rad.dat:"); 00494 } 00495 00496 ioDATA = open_data( "fe2rad.dat", "r" ); 00497 00498 /* get the first line, this is a version number */ 00499 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 00500 { 00501 fprintf( ioQQQ, " fe2rad.dat error reading data\n" ); 00502 cdEXIT(EXIT_FAILURE); 00503 } 00504 /* scan off three integers */ 00505 sscanf( chLine ,"%ld%ld%ld",&lo, &ihi, &m1); 00506 if( lo!=3 || ihi!=1 || m1!=28 ) 00507 { 00508 fprintf( ioQQQ, " fe2rad.dat has the wrong magic numbers, expected 3 1 28\n" ); 00509 cdEXIT(EXIT_FAILURE); 00510 } 00511 00512 /* file now open, read the data */ 00513 for( ipHi=1; ipHi < FeII.nFeIILevelAlloc; ipHi++ ) 00514 { 00515 for( ipLo=0; ipLo < ipHi; ipLo++ ) 00516 { 00517 /* following double since used in sscanf */ 00518 double save[2]; 00519 /* keep reading until have non-comment line, one that does not 00520 * start with # */ 00521 chLine[0] = '#'; 00522 while( chLine[0] == '#' ) 00523 { 00524 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 00525 { 00526 fprintf( ioQQQ, " fe2nn.dat error reading data\n" ); 00527 cdEXIT(EXIT_FAILURE); 00528 } 00529 } 00530 00531 /* first and last number on this line */ 00532 sscanf( chLine , 00533 "%ld%ld%ld%ld%lf%lf%ld", 00534 &lo, &ihi, &m1, &m2 , 00535 &save[0], &save[1] , &m3); 00536 /* the pointeres ihi and lo within this array were 00537 * read in from the line above */ 00538 Fe2LevN[ihi-1][lo-1].Lo->g = (realnum)m1; 00539 Fe2LevN[ihi-1][lo-1].Hi->g = (realnum)m2; 00540 Fe2LevN[ihi-1][lo-1].Emis->Aul = (realnum)save[0]; 00541 /*>>chng 06 apr 10, option to use table of energies */ 00542 # define USE_OLD true 00543 if( USE_OLD ) 00544 Fe2LevN[ihi-1][lo-1].EnergyWN = (realnum)save[1]; 00545 else 00546 Fe2LevN[ihi-1][lo-1].EnergyWN = Fe2Energies[ihi-1]-Fe2Energies[lo-1]; 00547 if( Fe2LevN[ihi-1][lo-1].Emis->Aul == 1e3f ) 00548 { 00549 /*realnum xmicron; 00550 xmicron = 1e4f/ Fe2LevN[ihi-1][lo-1].EnergyWN;*/ 00551 /* these are made-up intercombination lines, set gf to 1e-5 00552 * then get better A */ 00553 Fe2LevN[ihi-1][lo-1].Emis->gf = 1e-5f; 00554 00555 Fe2LevN[ihi-1][lo-1].Emis->Aul = Fe2LevN[ihi-1][lo-1].Emis->gf*(realnum)TRANS_PROB_CONST* 00556 POW2(Fe2LevN[ihi-1][lo-1].EnergyWN)*Fe2LevN[ihi-1][lo-1].Lo->g/Fe2LevN[ihi-1][lo-1].Hi->g; 00557 00558 /*fprintf(ioQQQ," %e from 1e3 to %e\n",xmicron , Fe2LevN[ihi-1][lo-1].Emis->Aul );*/ 00559 } 00560 /* this is the last column of fe2rad.dat, and is 1, 2, or 3. 00561 * 1 signifies that transition is permitted, 00562 * 2 is semi-forbidden 00563 * 3 forbidden, within lowest 63 levels are forbidden, first permitted 00564 * transition is from 64 */ 00565 ncs1[ihi-1][lo-1] = (int)m3; 00566 /*if( m3==1 ) 00567 fprintf(ioQQQ,"DEBUG fe2 new old energ\t%.5e\t%.5e\t%.3e\n", 00568 Fe2Energies[ihi-1]-Fe2Energies[lo-1], 00569 Fe2LevN[ihi-1][lo-1].EnergyWN , 00570 (Fe2Energies[ihi-1]-Fe2Energies[lo-1]-Fe2LevN[ihi-1][lo-1].EnergyWN)/ 00571 Fe2LevN[ihi-1][lo-1].EnergyWN 00572 );*/ 00573 } 00574 } 00575 fclose(ioDATA); 00576 00577 /* now read collision data in fe2col.dat 00578 * These are from the following sources 00579 >>refer fe2 cs Zhang, H.L., & Pradhan, A.K., 1995, A&A 293, 953 00580 >>refer fe2 cs Bautista, M., (private communication), 00581 >>refer fe2 cs Mewe, R., 1972, A&AS 20, 215 (the g-bar approximation) 00582 */ 00583 00584 if( trace.lgTrace ) 00585 { 00586 fprintf( ioQQQ," FeIICreate opening fe2col.dat:"); 00587 } 00588 00589 ioDATA = open_data( "fe2col.dat", "r" ); 00590 00591 ASSERT( ioDATA !=NULL); 00592 for( ipHi=1; ipHi<NPRADDAT; ++ipHi ) 00593 { 00594 for( ipLo=0; ipLo<ipHi; ++ipLo ) 00595 { 00596 /* double since used in sscanf */ 00597 double save[8]; 00598 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 00599 { 00600 fprintf( ioQQQ, " fe2col.dat error reading data\n" ); 00601 cdEXIT(EXIT_FAILURE); 00602 } 00603 sscanf( chLine , 00604 "%ld%ld%lf%lf%lf%lf%lf%lf%lf%lf", 00605 &m1, &m2, 00606 &save[0], &save[1] , &save[2],&save[3], &save[4] , &save[5], 00607 &save[6], &save[7] 00608 ); 00609 for( k=0; k<8; ++k ) 00610 { 00611 /* the max is here because there are some zeroes in the data file. 00612 * this is unphysical but is part of their distribution. As a result 00613 * don't let the cs fall below 0.01 */ 00614 sPradDat[m2-1][m1-1][k] = max(CS2SMALL , (realnum)save[k] ); 00615 } 00616 } 00617 } 00618 00619 fclose( ioDATA ); 00620 00621 /*generate needed opacity data for the large FeII atom */ 00622 00623 /* this routine is called only one time when cloudy is init 00624 * for the very first time. It converts the FeII data stored 00625 * in the large FeII arrays into the array storage needed by cloudy 00626 * for its line optical depth arrays 00627 */ 00628 00629 /* convert large FeII line arrays into standard heavy el ar */ 00630 for( ipLo=0; ipLo < (FeII.nFeIILevelAlloc - 1); ipLo++ ) 00631 { 00632 for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevelAlloc; ipHi++ ) 00633 { 00634 /* pull information out of existing data arrays */ 00635 00636 ASSERT( Fe2LevN[ipHi][ipLo].EnergyWN > 0. ); 00637 ASSERT( Fe2LevN[ipHi][ipLo].Emis->Aul> 0. ); 00638 00639 /* now put into standard internal line format 00640 Fe2LevN[ipHi][ipLo].WLAng = (realnum)(1.e8/Fe2LevN[ipHi][ipLo].EnergyWN); */ 00641 /* >>chng 03 oct 28, above neglected index of refraction of air - 00642 * convert to below */ 00643 Fe2LevN[ipHi][ipLo].WLAng = 00644 (realnum)(1.0e8/ 00645 Fe2LevN[ipHi][ipLo].EnergyWN/ 00646 RefIndex( Fe2LevN[ipHi][ipLo].EnergyWN )); 00647 00648 /* generate gf from A */ 00649 Fe2LevN[ipHi][ipLo].Emis->gf = 00650 (realnum)(Fe2LevN[ipHi][ipLo].Emis->Aul*Fe2LevN[ipHi][ipLo].Hi->g/ 00651 TRANS_PROB_CONST/POW2(Fe2LevN[ipHi][ipLo].EnergyWN)); 00652 00653 Fe2LevN[ipHi][ipLo].Hi->IonStg = 2; 00654 Fe2LevN[ipHi][ipLo].Hi->nelem = 26; 00655 /* which redistribution function?? 00656 * For resonance line use K2 (-1), 00657 * for subordinate lines use CRD with wings */ 00658 /* >>chng 01 mar 09, all had been 1, inc redis with wings */ 00659 /* to reset this, so that the code works as it did pre 01 mar 29, 00660 * use command 00661 * atom FeII redistribution resonance pdr 00662 * atom FeII redistribution subordinate pdr */ 00663 if( ipLo == 0 ) 00664 { 00665 Fe2LevN[ipHi][ipLo].Emis->iRedisFun = FeII.ipRedisFcnResonance; 00666 } 00667 else 00668 { 00669 /* >>chng 01 feb 27, had been -1, crd with core only, 00670 * change to crd with wings as per discussion with Ivan Hubeny */ 00671 Fe2LevN[ipHi][ipLo].Emis->iRedisFun = FeII.ipRedisFcnSubordinate; 00672 } 00673 Fe2LevN[ipHi][ipLo].Emis->phots = 0.; 00674 Fe2LevN[ipHi][ipLo].Emis->ots = 0.; 00675 Fe2LevN[ipHi][ipLo].Emis->FracInwd = 1.; 00676 } 00677 } 00678 00679 /* finally get FeII bands, this sets */ 00680 if( FeIIBandsCreate("") ) 00681 { 00682 fprintf( ioQQQ," failed to open FeII bands file \n"); 00683 cdEXIT(EXIT_FAILURE); 00684 } 00685 00686 /* now establish the FeII continuum, these are set to 00687 * 1000, 7000, and 1000 in FeIIZero in this file, and 00688 * are reset with the atom FeII continuum command */ 00689 FeIIContCreate( FeII.fe2con_wl1 , FeII.fe2con_wl2 , FeII.nfe2con ); 00690 00691 /* this must be true */ 00692 ASSERT( FeII.nFeIILevelAlloc == FeII.nFeIILevel ); 00693 00694 return; 00695 } 00696 00697 /* 00698 *===================================================================== 00699 */ 00700 /*********************************************************************** 00701 *** version of FeIILevelPops.f with overlap in procces 05.28.97 ooooooo 00702 ******change in common block *te* sqrte 05.28.97 00703 *******texc is fixed 03.28.97 00704 *********this version has work on overlap 00705 *********this version has # of zones (ZoneCnt) 07.03.97 00706 *********taux - optical depth depends on iter correction 03.03.97 00707 *********calculate cooling (Fe2_large_cool) and output fecool from Cloudy 01.29.97god 00708 *********and cooling derivative (ddT_Fe2_large_cool) 00709 ************ this program for 371 level iron model 12/14/1995 00710 ************ this program for 371 level iron model 1/11/1996 00711 ************ this program for 371 level iron model 3/21/1996 00712 ************ this program without La 3/27/1996 00713 ************ this program for 371 level iron model 4/9/1996 00714 ************ includes:FeIICollRatesBoltzmann;cooling;overlapping in lines */ 00715 void FeIILevelPops( void ) 00716 { 00717 long int i, 00718 ipHi , 00719 ipLo , 00720 n; 00721 /* used in test for non-positive level populations */ 00722 bool lgPopNeg; 00723 00724 double 00725 EnergyWN, 00726 pop , 00727 sum; 00728 00729 int32 info; 00730 int32 ipiv[NFE2LEVN]; 00731 00732 DEBUG_ENTRY( "FeIILevelPops()" ); 00733 00734 if( trace.lgTrace ) 00735 { 00736 fprintf( ioQQQ," FeIILevelPops fe2 pops called\n"); 00737 } 00738 00739 /* FeII.lgSimulate was set true with simulate flag on atom FeII command, 00740 * for bebugging without actually calling the routine */ 00741 if( FeII.lgSimulate ) 00742 { 00743 00744 return; 00745 } 00746 00747 /* zero out some arrays */ 00748 for( n=0; n<FeII.nFeIILevel; ++n) 00749 { 00750 for( ipLo=0; ipLo<FeII.nFeIILevel; ++ipLo ) 00751 { 00752 Fe2CPump[ipLo][n] = 0.; 00753 Fe2LPump[ipLo][n] = 0.; 00754 Fe2A[ipLo][n] = 0.; 00755 xMatrix[ipLo][n] = 0.; 00756 } 00757 } 00758 00759 /* make the g-bar collision strengths and do linear interpolation on r-matrix data. 00760 * this also sets Boltzmann factors for all levels, 00761 * sets values of FeColl used below, but only if temp has changed */ 00762 FeIICollRatesBoltzmann(); 00763 00764 /* pumping and spontantous decays */ 00765 for( n=0; n<FeII.nFeIILevel; ++n) 00766 { 00767 for( ipHi=n+1; ipHi<FeII.nFeIILevel; ++ipHi ) 00768 { 00769 /* continuum pumping rate from n to upper ipHi */ 00770 Fe2CPump[n][ipHi] = Fe2LevN[ipHi][n].Emis->pump; 00771 00772 /* continuum pumping rate from ipHi to lower n */ 00773 Fe2CPump[ipHi][n] = Fe2LevN[ipHi][n].Emis->pump* 00774 Fe2LevN[ipHi][n].Hi->g/Fe2LevN[ipHi][n].Lo->g; 00775 00776 /* spontaneous decays */ 00777 Fe2A[ipHi][n] = Fe2LevN[ipHi][n].Emis->Aul*(Fe2LevN[ipHi][n].Emis->Pesc + Fe2LevN[ipHi][n].Emis->Pelec_esc + 00778 Fe2LevN[ipHi][n].Emis->Pdest); 00779 } 00780 } 00781 00782 /* now do pumping of atom by Lya */ 00783 FeIILyaPump(); 00784 00785 /* ************************************************************************** 00786 * 00787 * final rates into matrix 00788 * 00789 ***************************************************************************/ 00790 00791 /* fill in xMatrix with matrix elements */ 00792 for( n=0; n<FeII.nFeIILevel; ++n) 00793 { 00794 /* all processes leaving level n going down*/ 00795 for( ipLo=0; ipLo<n; ++ipLo ) 00796 { 00797 xMatrix[n][n] = xMatrix[n][n] + Fe2CPump[n][ipLo] + Fe2LPump[n][ipLo]+ Fe2A[n][ipLo] + 00798 Fe2Coll[n][ipLo]*dense.eden; 00799 } 00800 /* all processes leaving level n going up*/ 00801 for( ipHi=n+1; ipHi<FeII.nFeIILevel; ++ipHi ) 00802 { 00803 xMatrix[n][n] = xMatrix[n][n] + Fe2CPump[n][ipHi] + Fe2LPump[n][ipHi] + Fe2Coll[n][ipHi]*dense.eden; 00804 } 00805 /* all processes entering level n from below*/ 00806 for( ipLo=0; ipLo<n; ++ipLo ) 00807 { 00808 xMatrix[ipLo][n] = xMatrix[ipLo][n] - Fe2CPump[ipLo][n] - Fe2LPump[ipLo][n] - Fe2Coll[ipLo][n]*dense.eden; 00809 } 00810 /* all processes entering level n from above*/ 00811 for( ipHi=n+1; ipHi<FeII.nFeIILevel; ++ipHi ) 00812 { 00813 xMatrix[ipHi][n] = xMatrix[ipHi][n] - Fe2CPump[ipHi][n] - Fe2LPump[ipHi][n] - Fe2Coll[ipHi][n]*dense.eden - 00814 Fe2A[ipHi][n]; 00815 } 00816 00817 /* the top row of the matrix is the sum of level populations */ 00818 xMatrix[n][0] = 1.0; 00819 } 00820 00821 { 00822 /* option to print out entire matrix */ 00823 enum {DEBUG_LOC=false}; 00824 if( DEBUG_LOC ) 00825 { 00826 /* print the matrices */ 00827 for( n=0; n<FeII.nFeIILevel; ++n) 00828 { 00829 /*fprintf(ioQQQ,"\n");*/ 00830 /* now print the matrix*/ 00831 for( ipLo=0; ipLo<FeII.nFeIILevel; ++ipLo ) 00832 { 00833 fprintf(ioQQQ," %.1e", xMatrix[n][ipLo]); 00834 } 00835 fprintf(ioQQQ,"\n"); 00836 } 00837 } 00838 } 00839 00840 /* define the Y Vector. The oth element is the sum of all level populations 00841 * adding up to the total population. The remaining elements are the level 00842 * balance equations adding up to zero */ 00843 yVector[0] = 1.0; 00844 for( n=1; n < FeII.nFeIILevel; n++ ) 00845 { 00846 yVector[n] = 0.0; 00847 } 00848 00849 /* create the 1-yVector array that will save vector, 00850 * this is the macro trick */ 00851 # ifdef AMAT 00852 # undef AMAT 00853 # endif 00854 # define AMAT(I_,J_) (*(amat+(I_)*FeII.nFeIILevel+(J_))) 00855 00856 /* copy current contents of xMatrix array over to special amat array*/ 00857 for( ipHi=0; ipHi < FeII.nFeIILevel; ipHi++ ) 00858 { 00859 for( i=0; i < FeII.nFeIILevel; i++ ) 00860 { 00861 AMAT(i,ipHi) = xMatrix[i][ipHi]; 00862 } 00863 } 00864 00865 info = 0; 00866 00867 /* do the linear algebra to find the level populations */ 00868 getrf_wrapper(FeII.nFeIILevel, FeII.nFeIILevel, amat, FeII.nFeIILevel, ipiv, &info); 00869 getrs_wrapper('N', FeII.nFeIILevel, 1, amat, FeII.nFeIILevel, ipiv, yVector, FeII.nFeIILevel, &info); 00870 00871 if( info != 0 ) 00872 { 00873 fprintf( ioQQQ, "DISASTER FeIILevelPops: dgetrs finds singular or ill-conditioned matrix\n" ); 00874 cdEXIT(EXIT_FAILURE); 00875 } 00876 00877 /* yVector now contains the level populations */ 00878 00879 /* this better be false after this loop - if not then non-positive level pops */ 00880 lgPopNeg = false; 00881 /* copy all level pops over to Fe2LevelPop */ 00882 for( ipLo=0; ipLo < FeII.nFeIILevel; ipLo++ ) 00883 { 00884 if(yVector[ipLo] < 0. ) 00885 { 00886 lgPopNeg = true; 00887 fprintf(ioQQQ,"PROBLEM FeIILevelPops finds non-positive level population, level is %ld pop is %g\n", 00888 ipLo , yVector[ipLo] ); 00889 } 00890 /* this is now correct level population, cm^-3 */ 00891 Fe2LevelPop[ipLo] = yVector[ipLo] * dense.xIonDense[ipIRON][1]; 00892 } 00893 if( lgPopNeg ) 00894 { 00895 /* this is here to use the lgPopNeg value for something, if not here then 00896 * lint and some compilers will note that var is never used */ 00897 fprintf(ioQQQ , "PROBLEM FeIILevelPops exits with negative level populations.\n"); 00898 } 00899 00900 /* >>chng 06 jun 24, make sure remainder of populations up through max 00901 * limit are zero - this makes safe the case where the number 00902 * of levels actually computed has been trimmed down from largest 00903 * possible number of levels, for instance, in cool gas */ 00904 for( ipLo=FeII.nFeIILevel; ipLo < FeII.nFeIILevelAlloc; ++ipLo ) 00905 { 00906 Fe2LevelPop[ipLo] = 0.; 00907 } 00908 00909 /* now set line opacities, intensities, and level populations 00910 * >>chng 06 jun ipLo should go up to FeII.nFeIILevel-1 since this 00911 * is the largest lower level with non-zero population */ 00912 for( ipLo=0; ipLo < (FeII.nFeIILevel - 1); ipLo++ ) 00913 { 00914 /* >>chng 06 jun 24, upper level should go to limit 00915 * of all that were allocated */ 00916 /*for( ipHi=ipLo+1; ipHi < FeII.nFeIILevel; ipHi++ )*/ 00917 for( ipHi=ipLo+1; ipHi < FeII.nFeIILevelAlloc; ipHi++ ) 00918 { 00919 /* >>chng 06 jun 24, in all of these the product 00920 * yVector[ipHi]*dense.xIonDense[ipIRON][1] has been replaced 00921 * with Fe2LevelPop[ipLo] - should always have been this way, 00922 * and saves a mult */ 00923 Fe2LevN[ipHi][ipLo].Emis->PopOpc = (Fe2LevelPop[ipLo] - 00924 Fe2LevelPop[ipHi]*Fe2LevN[ipHi][ipLo].Lo->g/Fe2LevN[ipHi][ipLo].Hi->g); 00925 00926 Fe2LevN[ipHi][ipLo].Lo->Pop = Fe2LevelPop[ipLo]; 00927 00928 Fe2LevN[ipHi][ipLo].Hi->Pop = Fe2LevelPop[ipHi]; 00929 00930 Fe2LevN[ipHi][ipLo].Emis->phots = Fe2LevelPop[ipHi]* 00931 Fe2LevN[ipHi][ipLo].Emis->Aul*(Fe2LevN[ipHi][ipLo].Emis->Pesc + Fe2LevN[ipHi][ipLo].Emis->Pelec_esc ); 00932 00933 Fe2LevN[ipHi][ipLo].Emis->xIntensity = Fe2LevN[ipHi][ipLo].Emis->phots* 00934 Fe2LevN[ipHi][ipLo].EnergyErg; 00935 00936 /* ratio of collisional (new) to pumped excitations */ 00937 /* >>chng 02 mar 04, add test MAX2 to prevent div by zero */ 00938 Fe2LevN[ipHi][ipLo].Emis->ColOvTot = (realnum)(Fe2Coll[ipLo][ipHi]*dense.eden / 00939 SDIV( Fe2Coll[ipLo][ipHi]*dense.eden + Fe2CPump[ipLo][ipHi] + Fe2LPump[ipLo][ipHi] ) ); 00940 } 00941 } 00942 00943 /* only do this if large atom is on and more than ground terms computed - 00944 * do not if only lowest levels are computed */ 00945 if( FeII.lgFeIILargeOn ) 00946 { 00947 /* the hydrogen Lya destruction rate, then probability */ 00948 hydro.dstfe2lya = 0.; 00949 EnergyWN = 0.; 00950 /* count how many photons were removed-added */ 00951 for( ipLo=0; ipLo < (FeII.nFeIILevel - 1); ipLo++ ) 00952 { 00953 for( ipHi=ipLo+1; ipHi < FeII.nFeIILevel; ipHi++ ) 00954 { 00955 EnergyWN += Fe2LPump[ipLo][ipHi] + Fe2LPump[ipHi][ipLo]; 00956 hydro.dstfe2lya += (realnum)( 00957 Fe2LevN[ipHi][ipLo].Lo->Pop*Fe2LPump[ipLo][ipHi] - 00958 Fe2LevN[ipHi][ipLo].Hi->Pop*Fe2LPump[ipHi][ipLo]); 00959 } 00960 } 00961 /* the destruction prob comes from 00962 * dest rate = n(2p) * A(21) * PDest */ 00963 pop = StatesElem[ipH_LIKE][ipHYDROGEN][ipH2p].Pop*dense.xIonDense[ipHYDROGEN][1]; 00964 if( pop > SMALLFLOAT ) 00965 { 00966 hydro.dstfe2lya /= (realnum)(pop * Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->Aul); 00967 } 00968 else 00969 { 00970 hydro.dstfe2lya = 0.; 00971 } 00972 /* NB - do not update hydro.dstfe2lya here if large FeII not on since 00973 * done in FeII overlap */ 00974 } 00975 00976 { 00977 enum {DEBUG_LOC=false}; 00978 if( DEBUG_LOC) 00979 { 00980 /*lint -e644 EnergyWN not init */ 00981 fprintf(ioQQQ," sum all %.1e dest rate%.1e escR= %.1e\n", 00982 EnergyWN,hydro.dstfe2lya, 00983 Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->Pesc); 00984 /*lint +e644 EnergyWN not init */ 00985 } 00986 } 00987 00988 /* next two blocks determine departure coefficients for the atom */ 00989 00990 /* first sum up partition function for the model atom */ 00991 Fe2DepCoef[0] = 1.0; 00992 sum = 1.0; 00993 for( i=1; i < FeII.nFeIILevel; i++ ) 00994 { 00995 /* Boltzmann factor relative to ground times ratio of statistical weights */ 00996 Fe2DepCoef[i] = Fe2DepCoef[0]*FeIIBoltzmann[i]* 00997 Fe2LevN[i][0].Hi->g/ Fe2LevN[1][0].Lo->g; 00998 /* this sum is the partition function for the atom */ 00999 sum += Fe2DepCoef[i]; 01000 } 01001 01002 /* now renormalize departure coefficients with partition function */ 01003 for( i=0; i < FeII.nFeIILevel; i++ ) 01004 { 01005 /* divide by total partition function, Fe2DepCoef is now the fraction of the 01006 * population that is in this level in TE */ 01007 Fe2DepCoef[i] /= sum; 01008 01009 /* convert to true departure coefficient */ 01010 if( Fe2DepCoef[i]>SMALLFLOAT ) 01011 { 01012 Fe2DepCoef[i] = yVector[i]/Fe2DepCoef[i]; 01013 } 01014 else 01015 { 01016 Fe2DepCoef[i] = 0.; 01017 } 01018 } 01019 01020 /*calculate cooling, heating, and cooling derivative */ 01021 /* this is net cooling, cooling minus heating */ 01022 FeII.Fe2_large_cool = 0.0f; 01023 /* this is be heating, not heating minus cooling */ 01024 FeII.Fe2_large_heat = 0.f; 01025 /* derivative of cooling */ 01026 FeII.ddT_Fe2_large_cool = 0.0f; 01027 01028 /* compute heating and cooling due to model atom */ 01029 for( ipLo=0; ipLo < (FeII.nFeIILevel - 1); ipLo++ ) 01030 { 01031 for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevel; ipHi++ ) 01032 { 01033 double OneLine; 01034 01035 /* net cooling due to single line */ 01036 OneLine = (Fe2Coll[ipLo][ipHi]*Fe2LevelPop[ipLo] - Fe2Coll[ipHi][ipLo]*Fe2LevelPop[ipHi])* 01037 dense.eden*Fe2LevN[ipHi][ipLo].EnergyErg; 01038 01039 /* net cooling due to this line */ 01040 FeII.Fe2_large_cool += (realnum)MAX2(0., OneLine); 01041 01042 /* net heating due to line */ 01043 FeII.Fe2_large_heat += (realnum)MAX2(0., -OneLine); 01044 01045 /* derivative of FeII cooling */ 01046 if( OneLine >= 0. ) 01047 { 01048 /* net coolant, exponential dominates deriv */ 01049 FeII.ddT_Fe2_large_cool += (realnum)OneLine* 01050 (Fe2LevN[ipHi][0].EnergyK*thermal.tsq1 - thermal.halfte); 01051 } 01052 else 01053 { 01054 /* net heating, te^-1/2 dominates */ 01055 FeII.ddT_Fe2_large_cool -= (realnum)OneLine*thermal.halfte; 01056 } 01057 } 01058 } 01059 01060 return; 01061 } 01062 01063 /*FeIICollRatesBoltzmann evaluate collision strenths, 01064 * both interpolating on r-mat and creating g-bar 01065 * find Boltzmann factors, evaluate collisional rate coefficients */ 01066 /* >>05 dec 03, verified that this rountine only changes temperature dependent 01067 * quantities - nothing with temperature dependence is done */ 01068 STATIC void FeIICollRatesBoltzmann(void) 01069 { 01070 /* this will be used to only reevaluate cs when temp changes */ 01071 static double OldTemp = -1.; 01072 long int i, 01073 ipLo , 01074 ipHi, 01075 ipTrim; 01076 realnum ag, 01077 cg, 01078 dg, 01079 gb, 01080 y; 01081 realnum FracLowTe , FracHighTe; 01082 static realnum tt[8]={1e3f,3e3f,5e3f,7e3f,1e4f,12e3f,15e3f,2e4f}; 01083 long int ipTemp, 01084 ipTempp1 , 01085 ipLoFe2, 01086 ipHiFe2; 01087 static long int nFeII_old = -1; 01088 01089 DEBUG_ENTRY( "FeIICollRatesBoltzmann()" ); 01090 01091 /* return if temperature has not changed */ 01092 /* >>chng 06 feb 14, add test on number of levels changing */ 01093 if( fp_equal( phycon.te, OldTemp ) && FeII.nFeIILevel == nFeII_old ) 01094 { 01095 01096 return; 01097 } 01098 OldTemp = phycon.te; 01099 nFeII_old = FeII.nFeIILevel; 01100 01101 /* find g-bar collision strength for all levels 01102 * expression for g-bar taken from 01103 *>>refer Fe2 g-bar Mewe, R. 1972, A&A, 20, 215 */ 01104 for( ipLo=0; ipLo < (FeII.nFeIILevel - 1); ipLo++ ) 01105 { 01106 for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevel; ipHi++ ) 01107 { 01108 /* these coefficients are from page 221 or Mewe 1972 */ 01109 if( ncs1[ipHi][ipLo] == 3 ) 01110 { 01111 /************************forbidden tr**************************/ 01112 ag = 0.15f; 01113 cg = 0.f; 01114 dg = 0.f; 01115 } 01116 /************************allowed tr*******************************/ 01117 else if( ncs1[ipHi][ipLo] == 1 ) 01118 { 01119 ag = 0.6f; 01120 cg = 0.f; 01121 dg = 0.28f; 01122 } 01123 /************************semiforbidden****************************/ 01124 else if( ncs1[ipHi][ipLo] == 2 ) 01125 { 01126 ag = 0.0f; 01127 cg = 0.1f; 01128 dg = 0.0f; 01129 } 01130 else 01131 { 01132 /* this branch is impossible, since cs must be 1, 2, or 3 */ 01133 ag = 0.0f; 01134 cg = 0.1f; 01135 dg = 0.0f; 01136 fprintf(ioQQQ,">>>PROBLEM impossible ncs1 in FeIILevelPops\n"); 01137 } 01138 01139 /*y = 1.438826f*Fe2LevN[ipHi][ipLo].EnergyWN/ phycon.te;*/ 01140 y = Fe2LevN[ipHi][ipLo].EnergyWN/ (realnum)phycon.te_wn; 01141 01142 gb = (realnum)(ag + (-cg*POW2(y) + dg)*(log(1.0+1.0/y) - 0.4/ 01143 POW2(y + 1.0)) + cg*y); 01144 01145 Fe2LevN[ipHi][ipLo].Coll.cs = 23.861f*1e5f*gb* 01146 Fe2LevN[ipHi][ipLo].Emis->Aul*Fe2LevN[ipHi][ipLo].Hi->g/ 01147 POW3(Fe2LevN[ipHi][ipLo].EnergyWN); 01148 01149 /* confirm that collision strength is positive */ 01150 ASSERT( Fe2LevN[ipHi][ipLo].Coll.cs > 0.); 01151 01152 /* g-bar cs becomes unphysically small for forbidden transitions - 01153 * this sets a lower limit to the final cs - CS2SMALL is defined above */ 01154 Fe2LevN[ipHi][ipLo].Coll.cs = MAX2( CS2SMALL, Fe2LevN[ipHi][ipLo].Coll.cs); 01155 /* this was the logic used in the old fortran version, 01156 * and reproduces the results in Baldwin et al '96 01157 if( Fe2LevN[ipHi][ipLo].cs < 1e-10 ) 01158 { 01159 Fe2LevN[ipHi][ipLo].cs = 0.01f; 01160 } 01161 */ 01162 } 01163 } 01164 /* end loop setting Mewe 1972 g-bar approximation */ 01165 01166 /* we will interpolate on the set of listed collision strengths - 01167 * where in this set are we? */ 01168 if( phycon.te <= tt[0] ) 01169 { 01170 /* temperature is lower than lowest tabulated, use the 01171 * lowest tabulated point */ 01172 /* ipTemp usually points to the cell cooler than te, ipTemp+1 to one higher, 01173 * here both are lowest */ 01174 ipTemp = 0; 01175 ipTempp1 = 0; 01176 FracHighTe = 0.; 01177 } 01178 else if( phycon.te > tt[7] ) 01179 { 01180 /* temperature is higher than highest tabulated, use the 01181 * highest tabulated point */ 01182 /* ipTemp usually points to the cell cooler than te, ipTemp+1 to one higher, 01183 * here both are highest */ 01184 ipTemp = 7; 01185 ipTempp1 = 7; 01186 FracHighTe = 0.; 01187 } 01188 else 01189 { 01190 /* where in the array is the temperature we need? */ 01191 ipTemp = -1; 01192 for( i=0; i < 8-1; i++ ) 01193 { 01194 if( phycon.te <= tt[i+1] ) 01195 { 01196 ipTemp = i; 01197 break; 01198 } 01199 01200 } 01201 /* this cannot possibly happen */ 01202 if( ipTemp < 0 ) 01203 { 01204 fprintf( ioQQQ, " Insanity while looking for temperature in coll str array, te=%g.\n", 01205 phycon.te ); 01206 cdEXIT(EXIT_FAILURE); 01207 } 01208 /* ipTemp points to the cell cooler than te, ipTemp+1 to one higher, 01209 * do linear interpolation between */ 01210 ipTemp = i; 01211 ipTempp1 = i+1; 01212 FracHighTe = ((realnum)phycon.te - tt[ipTemp])/(tt[ipTempp1] - tt[ipTemp]); 01213 } 01214 01215 /* this is fraction of final linear interpolated collision strength that 01216 * is weighted by the lower bound cs */ 01217 FracLowTe = 1.f - FracHighTe; 01218 01219 for( ipHi=1; ipHi < NPRADDAT; ipHi++ ) 01220 { 01221 for( ipLo=0; ipLo < ipHi; ipLo++ ) 01222 { 01223 /* ipHiFe2 should point to upper level of this pair, and 01224 * ipLoFe2 should point to lower level */ 01225 ipHiFe2 = MAX2( nnPradDat[ipHi] , nnPradDat[ipLo] ); 01226 ipLoFe2 = MIN2( nnPradDat[ipHi] , nnPradDat[ipLo] ); 01227 ASSERT( ipHiFe2-1 < NFE2LEVN ); 01228 ASSERT( ipHiFe2-1 >= 0 ); 01229 ASSERT( ipLoFe2-1 < NFE2LEVN ); 01230 ASSERT( ipLoFe2-1 >= 0 ); 01231 01232 /* do linear interpolation for CS, these are dimensioned NPRADDAT = 159 */ 01233 if( ipHiFe2-1 < FeII.nFeIILevel ) 01234 { 01235 /*fprintf(ioQQQ,"DEBUG\t%.2e", Fe2LevN[ipHiFe2-1][ipLoFe2-1].cs);*/ 01236 /* do linear interpolation */ 01237 Fe2LevN[ipHiFe2-1][ipLoFe2-1].Coll.cs = 01238 sPradDat[ipHi][ipLo][ipTemp] * FracLowTe + 01239 sPradDat[ipHi][ipLo][ipTempp1] * FracHighTe; 01240 /*fprintf(ioQQQ,"\t%.2e\n", Fe2LevN[ipHiFe2-1][ipLoFe2-1].cs);*/ 01241 01242 /* confirm that this is positive */ 01243 ASSERT( Fe2LevN[ipHiFe2-1][ipLoFe2-1].Coll.cs > 0. ); 01244 } 01245 } 01246 } 01247 01248 /* create Boltzmann factors for all levels */ 01249 FeIIBoltzmann[0] = 1.0; 01250 for( ipHi=1; ipHi < FeII.nFeIILevel; ipHi++ ) 01251 { 01252 /*FeIIBoltzmann[ipHi] = sexp( 1.438826f*Fe2LevN[ipHi][0].EnergyWN/phycon.te );*/ 01253 /* >>chng 99 may 21, from above to following slight different number from phyconst.h 01254 * >>chng 05 dec 01, from sexp to dsexp to avoid all 0 Boltzmann factors in low temps 01255 * now that FeII is on all the time */ 01256 FeIIBoltzmann[ipHi] = dsexp( Fe2LevN[ipHi][0].EnergyWN/phycon.te_wn ); 01257 } 01258 01259 /* now possibly trim down atom if Boltzmann factors for upper levels are zero */ 01260 ipTrim = FeII.nFeIILevel - 1; 01261 while( FeIIBoltzmann[ipTrim] == 0. && ipTrim > 0 ) 01262 { 01263 --ipTrim; 01264 } 01265 /* ipTrim now is the highest level with finite Boltzmann factor - 01266 * use this as the number of levels 01267 *>>chng 05 dec 01, from <=1 to <=0 - func_map does 10K, and large FeII had never 01268 * been tested in that limit. 1 is ok. */ 01269 ASSERT( ipTrim > 0 ); 01270 /* trim FeII.nFeIILevel so that FeIIBoltzmann is positive 01271 * in nearly all cases this does nothing since ipTrim and FeII.nFeIILevel 01272 * are equal . . . */ 01273 if( ipTrim+1 < FeII.nFeIILevel ) 01274 { 01275 /* >>chng 06 jun 27, zero out collision data */ 01276 for( ipLo=0; ipLo<FeII.nFeIILevel; ++ipLo ) 01277 { 01278 for( ipHi=ipTrim; ipHi<FeII.nFeIILevel; ++ipHi ) 01279 { 01280 Fe2Coll[ipLo][ipHi] = 0.; 01281 Fe2Coll[ipHi][ipLo] = 0.; 01282 } 01283 } 01284 } 01285 FeII.nFeIILevel = ipTrim+1; 01286 /*fprintf(ioQQQ," levels reset to %li\n",FeII.nFeIILevel);*/ 01287 01288 /* debug code to print out the collision strengths for some levels */ 01289 { 01290 enum {DEBUG_LOC=false}; 01291 if( DEBUG_LOC) 01292 { 01293 for( ipLo=0; ipLo < 52; ipLo++ ) 01294 { 01295 fprintf(ioQQQ,"%e %e\n", 01296 Fe2LevN[51][ipLo].Coll.cs,Fe2LevN[52][ipLo].Coll.cs); 01297 } 01298 cdEXIT(EXIT_FAILURE); 01299 } 01300 } 01301 01302 /* collisional excitation and deexcitation */ 01303 for( ipLo=0; ipLo<FeII.nFeIILevel; ++ipLo) 01304 { 01305 for( ipHi=ipLo+1; ipHi<FeII.nFeIILevel; ++ipHi ) 01306 { 01307 /* collisional deexcitation rate coefficient from ipHi to lower ipLo 01308 * note that it needs eden to become rate 01309 * units cm3 s-1 */ 01310 Fe2Coll[ipHi][ipLo] = (realnum)(COLL_CONST/phycon.sqrte*Fe2LevN[ipHi][ipLo].Coll.cs/ 01311 Fe2LevN[ipHi][ipLo].Hi->g); 01312 /*fprintf(ioQQQ,"DEBUG fe2coll %li %li %.2e", ipHi , ipLo , Fe2Coll[ipHi][ipLo] );*/ 01313 01314 /* collisional excitation rate coefficient from ipLo to upper ipHi 01315 * units cm3 s-1 01316 * FeIIBoltzmann guaranteed to be > 0 by nFeIILevel trimming */ 01317 Fe2Coll[ipLo][ipHi] = (realnum)(Fe2Coll[ipHi][ipLo]*FeIIBoltzmann[ipHi]/FeIIBoltzmann[ipLo]* 01318 Fe2LevN[ipHi][ipLo].Hi->g/Fe2LevN[ipHi][ipLo].Lo->g); 01319 /*fprintf(ioQQQ,"DEBUG fe2coll %.2e\n", Fe2Coll[ipLo][ipHi] );*/ 01320 } 01321 } 01322 01323 return; 01324 } 01325 01326 /* 01327 *===================================================================== 01328 */ 01329 /*subroutine FeIIPrint PhotOccNum_at_nu raspechatki naselennostej v cloudy.out ili v 01330 * otdel'nom file unit=33 01331 *!!nado takzhe vklyuchit raspechatku iz perekrytiya linii */ 01332 /*FeIIPrint - print output from large FeII atom, called by prtzone */ 01333 void FeIIPrint(void) 01334 { 01335 01336 DEBUG_ENTRY( "FeIIPrint()" ); 01337 01338 return; 01339 } 01340 01341 /* 01342 *===================================================================== 01343 */ 01344 /*FeIISumBand, sum up large FeII emission over certain bands, called in lineset4 01345 * units are erg s-1 cm-3, same units as xIntensity in line structure */ 01346 /* >>chng 06 apr 11, from using erg as energy unit to angstroms, 01347 * this fixes bug in physical constant and also uses air wavelengths for wl > 2000A */ 01348 double FeIISumBand( 01349 /* lower and upper range to wavelength in Angstroms */ 01350 realnum wl1, 01351 realnum wl2) 01352 { 01353 long int ipHi, 01354 ipLo; 01355 double SumBandFe2_v; 01356 /* >>chng 00 jun 02, demoted next two to realnum, PvH 01357 realnum ahi, 01358 alo;*/ 01359 01360 DEBUG_ENTRY( "FeIISumBand()" ); 01361 /*sum up large FeII emission over certain bands */ 01362 01363 if( dense.xIonDense[ipIRON][1] == 0. ) 01364 { 01365 01366 return( 0. ); 01367 } 01368 else 01369 { 01370 01371 SumBandFe2_v = 0.; 01372 01373 /* convert to line energy in ergs */ 01374 /*ahi = 1.99e-8f/wl1; 01375 alo = 1.99e-8f/wl2;*/ 01376 /*>>chng 06 apr 10, from above to below, imprecise converstion 01377 * factor for Angstroms into ergs, caused shift in spectrum by 01378 * about 600 km/s. Bug caught by Yumihiko Tsuzuki */ 01379 # if 0 01380 /* convert to line energy in ergs */ 01381 ahi = 1.99e-8f/wl1; 01382 alo = 1.99e-8f/wl2; 01383 ASSERT( ahi > alo ); 01384 for( ipLo=0; ipLo < (FeII.nFeIILevel - 1); ipLo++ ) 01385 { 01386 for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevel; ipHi++ ) 01387 { 01388 if( Fe2LevN[ipHi][ipLo].EnergyErg >= alo && 01389 Fe2LevN[ipHi][ipLo].EnergyErg <= ahi ) 01390 SumBandFe2_v += Fe2LevN[ipHi][ipLo].Emis->xIntensity; 01391 } 01392 } 01393 # endif 01394 01395 /* >>chng 06 apr 10, from above, with line energy in ergs, 01396 * to below, line energy in wavenumbers */ 01397 ASSERT( wl2 > wl1 ); 01398 for( ipHi=1; ipHi < FeII.nFeIILevel; ++ipHi ) 01399 { 01400 for( ipLo=0; ipLo < ipHi; ++ipLo ) 01401 { 01402 if( Fe2LevN[ipHi][ipLo].WLAng >= wl1 && 01403 Fe2LevN[ipHi][ipLo].WLAng < wl2 ) 01404 SumBandFe2_v += Fe2LevN[ipHi][ipLo].Emis->xIntensity; 01405 } 01406 } 01407 01408 return( SumBandFe2_v ); 01409 } 01410 } 01411 01412 /* 01413 *===================================================================== 01414 */ 01415 /*FeII_RT_TauInc called once per zone in RT_tau_inc to increment large FeII atom line optical depths */ 01416 void FeII_RT_TauInc(void) 01417 { 01418 long int ipHi, 01419 ipLo; 01420 01421 DEBUG_ENTRY( "FeII_RT_TauInc()" ); 01422 01423 for( ipLo=0; ipLo < (FeII.nFeIILevel - 1); ipLo++ ) 01424 { 01425 /* >>chng 06 jun 24, for upper level go all the way to the 01426 * largest possible number of levels so that optical depths 01427 * of UV transitions are correct even for very cold gas where 01428 * the high level populations are not computed */ 01429 for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevelAlloc; ipHi++ ) 01430 { 01431 /* >>chng 04 aug 28, add test on bogus line */ 01432 if( Fe2LevN[ipHi][ipLo].ipCont > 0 ) 01433 RT_line_one_tauinc( &Fe2LevN[ipHi][ipLo] , 01434 -8 , -8 , ipHi , ipLo ); 01435 } 01436 } 01437 01438 return; 01439 } 01440 01441 /* 01442 *===================================================================== 01443 */ 01444 /*FeII_RT_tau_reset reset optical depths for large FeII atom, called by update after each iteration */ 01445 void FeII_RT_tau_reset(void) 01446 { 01447 long int ipHi, 01448 ipLo; 01449 01450 DEBUG_ENTRY( "FeII_RT_tau_reset()" ); 01451 01452 /*fprintf(ioQQQ,"DEBUG FeIITauAver1 %li %.2e %.2e nFeIILevel %li \n", 01453 iteration, 01454 Fe2LevN[9][0].Emis->TauIn , 01455 Fe2LevN[9][0].Emis->TauTot, 01456 FeII.nFeIILevel);*/ 01457 01458 /* called at end of iteration */ 01459 /* >>chng 05 dec 07, had been FeII.nFeIILevel but this may have been trimmed down 01460 * if previous iteration went very deep - not reset until FeIIReset called */ 01461 for( ipHi=1; ipHi < FeII.nFeIILevelAlloc; ipHi++ ) 01462 { 01463 for( ipLo=0; ipLo < ipHi; ipLo++ ) 01464 { 01465 RT_line_one_tau_reset( &Fe2LevN[ipHi][ipLo],0.5 ); 01466 } 01467 } 01468 /*fprintf(ioQQQ,"DEBUG FeIITauAver2 %li %.2e %.2e\n", 01469 iteration, 01470 Fe2LevN[9][0].Emis->TauIn , 01471 Fe2LevN[9][0].Emis->TauTot);*/ 01472 01473 /* call the line overlap routine */ 01474 FeIIOvrLap(); 01475 /*fprintf(ioQQQ,"DEBUG FeIITauAver3 %li %.2e %.2e\n", 01476 iteration, 01477 Fe2LevN[9][0].Emis->TauIn , 01478 Fe2LevN[9][0].Emis->TauTot);*/ 01479 01480 return; 01481 } 01482 01483 /* 01484 *===================================================================== 01485 */ 01486 /*FeIIPoint called by ContCreatePointers to create pointers for lines in large FeII atom */ 01487 void FeIIPoint(void) 01488 { 01489 long int ipHi, 01490 ip, 01491 ipLo; 01492 01493 DEBUG_ENTRY( "FeIIPoint()" ); 01494 01495 /* routine called when cloudy sets continuum array indices for Fe2 lines, 01496 * once per coreload */ 01497 for( ipLo=0; ipLo < FeII.nFeIILevel-1; ipLo++ ) 01498 { 01499 for( ipHi=ipLo+1; ipHi < FeII.nFeIILevel; ipHi++ ) 01500 { 01501 01502 /* >>chng 02 feb 11, set continuum index to negative value for fake transition */ 01503 if( fabs(Fe2LevN[ipHi][ipLo].Emis->Aul- 1e-5 ) > 1e-8 ) 01504 { 01505 ip = ipoint(Fe2LevN[ipHi][ipLo].EnergyWN * WAVNRYD); 01506 Fe2LevN[ipHi][ipLo].ipCont = ip; 01507 01508 /* do not over write other pointers with FeII since FeII is everywhere */ 01509 if( strcmp( rfield.chLineLabel[ip-1], " " ) == 0 ) 01510 strcpy( rfield.chLineLabel[ip-1], "FeII" ); 01511 01512 Fe2LevN[ipHi][ipLo].Emis->ipFine = ipFineCont( Fe2LevN[ipHi][ipLo].EnergyWN * WAVNRYD); 01513 } 01514 else 01515 { 01516 Fe2LevN[ipHi][ipLo].ipCont = -1; 01517 Fe2LevN[ipHi][ipLo].Emis->ipFine = -1; 01518 } 01519 01520 Fe2LevN[ipHi][ipLo].Emis->dampXvel = 01521 (realnum)(Fe2LevN[ipHi][ipLo].Emis->Aul/ 01522 Fe2LevN[ipHi][ipLo].EnergyWN/PI4); 01523 01524 /* derive the abs coef, call to function is gf, wl (A), g_low */ 01525 Fe2LevN[ipHi][ipLo].Emis->opacity = 01526 (realnum)abscf(Fe2LevN[ipHi][ipLo].Emis->gf, 01527 Fe2LevN[ipHi][ipLo].EnergyWN, 01528 Fe2LevN[ipHi][ipLo].Lo->g); 01529 01530 /* excitation energy of transition in degrees kelvin */ 01531 Fe2LevN[ipHi][ipLo].EnergyK = 01532 (realnum)(T1CM*Fe2LevN[ipHi][ipLo].EnergyWN); 01533 01534 /* energy of photon in ergs */ 01535 Fe2LevN[ipHi][ipLo].EnergyErg = 01536 (realnum)(ERG1CM*Fe2LevN[ipHi][ipLo].EnergyWN); 01537 } 01538 } 01539 01540 return; 01541 } 01542 01543 /* 01544 *===================================================================== 01545 */ 01546 /*FeIIAccel called by rt_line_driving to compute radiative acceleration due to FeII lines */ 01547 void FeIIAccel(double *fe2drive) 01548 { 01549 long int ipHi, 01550 ipLo; 01551 01552 DEBUG_ENTRY( "FeIIAccel()" ); 01553 /*compute acceleration due to large FeII atom */ 01554 01555 /* this routine computes the line driven radiative acceleration 01556 * due to large FeII atom*/ 01557 01558 *fe2drive = 0.; 01559 for( ipLo=0; ipLo < (FeII.nFeIILevel - 1); ipLo++ ) 01560 { 01561 for( ipHi=ipLo+1; ipHi < FeII.nFeIILevel; ipHi++ ) 01562 { 01563 *fe2drive += Fe2LevN[ipHi][ipLo].Emis->pump* 01564 Fe2LevN[ipHi][ipLo].EnergyErg*Fe2LevN[ipHi][ipLo].Emis->PopOpc; 01565 } 01566 } 01567 01568 return; 01569 } 01570 01571 /* 01572 *===================================================================== 01573 */ 01574 /*FeII_RT_Make called by RT_line_all, does large FeII atom radiative transfer */ 01575 void FeII_RT_Make( bool lgDoEsc , bool lgUpdateFineOpac ) 01576 { 01577 long int ipHi, 01578 ipLo; 01579 01580 DEBUG_ENTRY( "FeII_RT_Make()" ); 01581 01582 if( trace.lgTrace ) 01583 { 01584 fprintf( ioQQQ," FeII_RT_Make called\n"); 01585 } 01586 01587 /* this routine drives calls to make RT relations with large FeII atom */ 01588 for( ipLo=0; ipLo < (FeII.nFeIILevel - 1); ipLo++ ) 01589 { 01590 /* >>chng 06 jun 24, go up to number allocated to keep UV resonance lines */ 01591 /*for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevel; ipHi++ )*/ 01592 for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevelAlloc; ipHi++ ) 01593 { 01594 /* only evaluate real transitions */ 01595 if( Fe2LevN[ipHi][ipLo].ipCont > 0 ) 01596 { 01597 /*RT_line_one do rt for emission line structure - 01598 * calls RT_line_static or RT_line_wind */ 01599 RT_line_one( &Fe2LevN[ipHi][ipLo] , lgDoEsc , lgUpdateFineOpac,true); 01600 } 01601 } 01602 } 01603 01604 return; 01605 } 01606 01607 /* 01608 *===================================================================== 01609 */ 01610 /* called in LineSet4 to add FeII lines to save array */ 01611 void FeIIAddLines( void ) 01612 { 01613 long int ipHi, 01614 ipLo; 01615 01616 DEBUG_ENTRY( "FeIIAddLines()" ); 01617 01618 /* this routine puts the emission from the large FeII atom 01619 * into an array to save and integrate them*/ 01620 01621 /* add lines option called from lines, add intensities into storage array */ 01622 01623 /* routine is called three different ways, ipass < 0 before space allocated, 01624 * =0 when time to generate labels (and we zero out array here), and ipass>0 01625 * when time to save intensities */ 01626 if( LineSave.ipass == 0 ) 01627 { 01628 /* we were called by lines, and we want to zero out Fe2SavN */ 01629 for( ipLo=0; ipLo < (FeII.nFeIILevel - 1); ipLo++ ) 01630 { 01631 for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevel; ipHi++ ) 01632 { 01633 Fe2SavN[ipHi][ipLo] = 0.; 01634 } 01635 } 01636 } 01637 01638 /* this call during calculations, save intensities */ 01639 else if( LineSave.ipass == 1 ) 01640 { 01641 /* total emission from vol */ 01642 for( ipLo=0; ipLo < (FeII.nFeIILevel - 1); ipLo++ ) 01643 { 01644 for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevel; ipHi++ ) 01645 { 01646 Fe2SavN[ipHi][ipLo] += 01647 radius.dVeff*Fe2LevN[ipHi][ipLo].Emis->xIntensity; 01648 } 01649 } 01650 } 01651 01652 { 01653 enum {DEBUG_LOC=false}; 01654 if( DEBUG_LOC /*&& (iteration==2)*/ ) 01655 { 01656 fprintf(ioQQQ," 69-1\t%li\t%e\n", nzone , Fe2SavN[68][0] ); 01657 } 01658 } 01659 01660 return; 01661 } 01662 01663 /*FeIIPunchLevels punch level energies and stat weights */ 01664 void FeIIPunchLevels( 01665 /* file we will punch to */ 01666 FILE *ioPUN ) 01667 { 01668 01669 long int ipHi; 01670 01671 DEBUG_ENTRY( "FeIIPunchLevels()" ); 01672 01673 fprintf(ioPUN , "%.2f\t%li\n", 0., (long)Fe2LevN[1][0].Lo->g ); 01674 for( ipHi=1; ipHi < FeII.nFeIILevel; ++ipHi ) 01675 { 01676 fprintf(ioPUN , "%.2f\t%li\n", 01677 Fe2LevN[ipHi][0].EnergyWN , 01678 (long)Fe2LevN[ipHi][0].Hi->g); 01679 } 01680 01681 return; 01682 } 01683 01684 /*FeIIPunchColden punch level energies, stat weights, column density */ 01685 void FeIIPunchColden( 01686 /* file we will punch to */ 01687 FILE *ioPUN ) 01688 { 01689 01690 long int n; 01691 01692 DEBUG_ENTRY( "FeIIPunchColden()" ); 01693 01694 fprintf(ioPUN , "%.2f\t%li\t%.3e\n", 0., (long)Fe2LevN[1][0].Lo->g , Fe2ColDen[0]); 01695 for( n=1; n < FeII.nFeIILevel; ++n ) 01696 { 01697 fprintf(ioPUN , "%.2f\t%li\t%.3e\n", 01698 Fe2LevN[n][0].EnergyWN , 01699 (long)Fe2LevN[n][0].Hi->g, 01700 Fe2ColDen[n]); 01701 } 01702 01703 return; 01704 } 01705 01706 01707 /*FeIIPunchOpticalDepth punch FeII optical depths, 01708 * called by punch_do to punch them out, 01709 * punch turned on with punch lines command */ 01710 void FeIIPunchOpticalDepth( 01711 /* file we will punch to */ 01712 FILE *ioPUN ) 01713 { 01714 long int 01715 ipHi, 01716 ipLo; 01717 01718 DEBUG_ENTRY( "FeIIPunchOpticalDepth()" ); 01719 01720 /* this routine punches the optical depths predicted by the large FeII atom */ 01721 /*>>chng 06 may 19, must use total number of levels allocated not current 01722 * number since this decreases as gas grows colder with depth nFeIILevel->nFeIILevelAlloc*/ 01723 for( ipLo=0; ipLo < (FeII.nFeIILevelAlloc - 1); ipLo++ ) 01724 { 01725 for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevelAlloc; ipHi++ ) 01726 { 01727 /* fe2ener(1) and (2) are lower and upper limits to range of 01728 * wavelengths of interest. entered in ryd with 01729 * punch FeII verner command, where they are converted 01730 * to wavenumbers, */ 01731 fprintf( ioPUN, "%ld\t%ld\t%.2f\t%.2e\n", 01732 ipHi+1, 01733 ipLo+1, 01734 Fe2LevN[ipHi][ipLo].WLAng , 01735 Fe2LevN[ipHi][ipLo].Emis->TauIn ); 01736 } 01737 } 01738 01739 return; 01740 } 01741 01742 /*FeIIPunchLines punch accumulated FeII intensities, at end of calculation, 01743 * called by punch_do to punch them out, 01744 * punch turned on with punch lines command */ 01745 void FeIIPunchLines( 01746 /* file we will punch to */ 01747 FILE *ioPUN ) 01748 { 01749 long int MaseHi, 01750 MaseLow, 01751 ipHi, 01752 ipLo; 01753 double hbeta, absint , renorm; 01754 /* >>chng 00 jun 02, demoted next two to realnum, PvH */ 01755 realnum TauMase, 01756 thresh; 01757 01758 DEBUG_ENTRY( "FeIIPunchLines()" ); 01759 01760 /* this routine puts the emission from the large FeII atom 01761 * into a line array, and eventually will punch it out */ 01762 01763 /* get the normalization line */ 01764 if( LineSv[LineSave.ipNormWavL].sumlin[LineSave.lgLineEmergent] > 0. ) 01765 { 01766 renorm = LineSave.ScaleNormLine/LineSv[LineSave.ipNormWavL].sumlin[LineSave.lgLineEmergent]; 01767 } 01768 else 01769 { 01770 renorm = 1.; 01771 } 01772 01773 fprintf( ioPUN, " up low log I, I/I(LineSave), Tau\n" ); 01774 01775 /* first look for any masing lines */ 01776 MaseLow = -1; 01777 MaseHi = -1; 01778 TauMase = 0.f; 01779 for( ipLo=0; ipLo < (FeII.nFeIILevelAlloc - 1); ipLo++ ) 01780 { 01781 for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevelAlloc; ipHi++ ) 01782 { 01783 if( Fe2LevN[ipHi][ipLo].Emis->TauIn < TauMase ) 01784 { 01785 TauMase = Fe2LevN[ipHi][ipLo].Emis->TauIn; 01786 MaseLow = ipLo; 01787 MaseHi = ipHi; 01788 } 01789 } 01790 } 01791 01792 if( TauMase < 0.f ) 01793 { 01794 fprintf( ioPUN, " Most negative optical depth was %4ld%4ld%10.2e\n", 01795 MaseHi, MaseLow, TauMase ); 01796 } 01797 01798 /* now print actual line intensities, Hbeta first */ 01799 if( cdLine("TOTL", 4861 , &hbeta , &absint)<=0 ) 01800 { 01801 fprintf( ioQQQ, " FeIILevelPops could not find Hbeta with cdLine.\n" ); 01802 cdEXIT(EXIT_FAILURE); 01803 } 01804 01805 fprintf( ioPUN, "Hbeta=%7.3f %10.2e\n", 01806 absint , 01807 hbeta ); 01808 01809 if( renorm > SMALLFLOAT ) 01810 { 01811 /* this is threshold for faintest line, normally 0, set with 01812 * number on punch verner command */ 01813 thresh = FeII.fe2thresh/(realnum)renorm; 01814 } 01815 else 01816 { 01817 thresh = 0.f; 01818 } 01819 01820 /*>>chng 06 may 19, must use total number of levels allocated not current 01821 * number since this decreases as gas grows colder with depth nFeIILevelAlloc*/ 01822 for( ipLo=0; ipLo < (FeII.nFeIILevelAlloc - 1); ipLo++ ) 01823 { 01824 for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevelAlloc; ipHi++ ) 01825 { 01826 /* fe2ener(1) and (2) are lower and upper limits to range of 01827 * wavelengths of interest. entered in ryd with 01828 * punch FeII verner command, where they are converted 01829 * to wavenumbers, */ 01830 if( (Fe2SavN[ipHi][ipLo] > thresh && 01831 Fe2LevN[ipHi][ipLo].EnergyWN > FeII.fe2ener[0]) && 01832 Fe2LevN[ipHi][ipLo].EnergyWN < FeII.fe2ener[1] ) 01833 { 01834 if( FeII.lgShortFe2 ) 01835 { 01836 /* short option on punch FeII line command 01837 * does not include rel intensity or optical dep */ 01838 fprintf( ioPUN, "%ld\t%ld\t%.2f\t%.3f\n", 01839 ipHi+1, 01840 ipLo+1, 01841 Fe2LevN[ipHi][ipLo].WLAng , 01842 log10(MAX2(1e-37,Fe2SavN[ipHi][ipLo])) + radius.Conv2PrtInten ); 01843 } 01844 else 01845 { 01846 /* long printout does */ 01847 fprintf( ioPUN, "%ld\t%ld\t%.2f\t%.3f\t%.2e\t%.2e\n", 01848 ipHi+1, 01849 ipLo+1, 01850 Fe2LevN[ipHi][ipLo].WLAng , 01851 log10(MAX2(1e-37,Fe2SavN[ipHi][ipLo])) + radius.Conv2PrtInten, 01852 Fe2SavN[ipHi][ipLo]*renorm , 01853 Fe2LevN[ipHi][ipLo].Emis->TauIn ); 01854 } 01855 } 01856 } 01857 } 01858 01859 return; 01860 } 01861 01862 01863 /* 01864 *===================================================================== 01865 */ 01866 /*FeII_LineZero zero out storage for large FeII atom, called in tauout */ 01867 void FeII_LineZero(void) 01868 { 01869 long int ipHi, 01870 ipLo; 01871 01872 DEBUG_ENTRY( "FeII_LineZero()" ); 01873 01874 /* this routine is called in routine zero and it 01875 * zero's out various elements of the FeII arrays 01876 * it is called on every iteration 01877 * */ 01878 for( ipLo=0; ipLo < (FeII.nFeIILevel - 1); ipLo++ ) 01879 { 01880 for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevel; ipHi++ ) 01881 { 01882 /* >>chng 03 feb 14, use EmLineZero rather than explicit sets */ 01883 TransitionZero( &Fe2LevN[ipHi][ipLo] ); 01884 } 01885 } 01886 01887 return; 01888 } 01889 /* 01890 *===================================================================== 01891 */ 01892 /*FeIIIntenZero zero out storage for large FeII atom, called in tauout */ 01893 void FeIIIntenZero(void) 01894 { 01895 long int ipHi, 01896 ipLo; 01897 01898 DEBUG_ENTRY( "FeIIIntenZero()" ); 01899 01900 /* this routine is called by routine cool_iron and it 01901 * zero's out various elements of the FeII arrays 01902 * */ 01903 Fe2LevelPop[0] = 0.; 01904 for( ipHi=1; ipHi < FeII.nFeIILevel; ipHi++ ) 01905 { 01906 /* >>chng 05 dec 14, zero Fe2LevelPop */ 01907 Fe2LevelPop[ipHi] = 0.; 01908 for( ipLo=0; ipLo < ipHi; ipLo++ ) 01909 { 01910 01911 /* population of lower level with correction for stim emission */ 01912 Fe2LevN[ipHi][ipLo].Emis->PopOpc = 0.; 01913 01914 /* population of lower level */ 01915 Fe2LevN[ipHi][ipLo].Lo->Pop = 0.; 01916 01917 /* population of upper level */ 01918 Fe2LevN[ipHi][ipLo].Hi->Pop = 0.; 01919 01920 /* following two heat exchange excitation, deexcitation */ 01921 Fe2LevN[ipHi][ipLo].Coll.cool = 0.; 01922 Fe2LevN[ipHi][ipLo].Coll.heat = 0.; 01923 01924 /* intensity of line */ 01925 Fe2LevN[ipHi][ipLo].Emis->xIntensity = 0.; 01926 01927 Fe2LevN[ipHi][ipLo].Emis->phots = 0.; 01928 Fe2LevN[ipHi][ipLo].Emis->ots = 0.; 01929 Fe2LevN[ipHi][ipLo].Emis->ColOvTot = 0.; 01930 } 01931 } 01932 01933 return; 01934 } 01935 01936 01937 /* 01938 *===================================================================== 01939 */ 01940 /* fill in IR lines from lowest 16 levels */ 01941 void FeIIFillLow16(void) 01942 { 01943 01944 DEBUG_ENTRY( "FeIIFillLow16()" ); 01945 01946 /* this is total cooling due to 16 level atom that would have been computed there */ 01947 /*FeII.Fe2_16levl_cool = 0.;*/ 01948 01949 /* all transitions within 16 levels of ground term */ 01950 FeII.fe21308 = Fe2LevN[12][7].Emis->xIntensity; 01951 FeII.fe21207 = Fe2LevN[11][6].Emis->xIntensity; 01952 FeII.fe21106 = Fe2LevN[10][5].Emis->xIntensity; 01953 FeII.fe21006 = Fe2LevN[9][5].Emis->xIntensity; 01954 FeII.fe21204 = Fe2LevN[11][3].Emis->xIntensity; 01955 FeII.fe21103 = Fe2LevN[10][2].Emis->xIntensity; 01956 FeII.fe21104 = Fe2LevN[10][3].Emis->xIntensity; 01957 FeII.fe21001 = Fe2LevN[9][0].Emis->xIntensity; 01958 FeII.fe21002 = Fe2LevN[9][1].Emis->xIntensity; 01959 FeII.fe20201 = Fe2LevN[1][0].Emis->xIntensity; 01960 FeII.fe20302 = Fe2LevN[2][1].Emis->xIntensity; 01961 FeII.fe20706 = Fe2LevN[6][5].Emis->xIntensity; 01962 FeII.fe20807 = Fe2LevN[7][6].Emis->xIntensity; 01963 FeII.fe20908 = Fe2LevN[8][7].Emis->xIntensity; 01964 FeII.fe21007 = Fe2LevN[9][6].Emis->xIntensity; 01965 FeII.fe21107 = Fe2LevN[10][6].Emis->xIntensity; 01966 FeII.fe21108 = Fe2LevN[10][7].Emis->xIntensity; 01967 FeII.fe21110 = Fe2LevN[10][9].Emis->xIntensity; 01968 FeII.fe21208 = Fe2LevN[11][7].Emis->xIntensity; 01969 FeII.fe21209 = Fe2LevN[11][8].Emis->xIntensity; 01970 FeII.fe21211 = Fe2LevN[11][10].Emis->xIntensity; 01971 FeII.fe21406 = Fe2LevN[13][5].Emis->xIntensity; 01972 FeII.fe21507 = Fe2LevN[14][6].Emis->xIntensity; 01973 FeII.fe21508 = Fe2LevN[14][7].Emis->xIntensity; 01974 FeII.fe21609 = Fe2LevN[15][8].Emis->xIntensity; 01975 01976 /* these are unique to this large model atom */ 01977 if( FeII.nFeIILevel > 43 ) 01978 { 01979 /* NB do not forget to decrement by one when going from physical (in left variables) 01980 * to c scale (in array) */ 01981 FeII.fe25to6 = Fe2LevN[25-1][5].Emis->xIntensity; 01982 FeII.fe27to7 = Fe2LevN[27-1][6].Emis->xIntensity; 01983 FeII.fe28to8 = Fe2LevN[28-1][7].Emis->xIntensity; 01984 FeII.fe29to9 = Fe2LevN[29-1][8].Emis->xIntensity; 01985 FeII.fe32to6 = Fe2LevN[32-1][5].Emis->xIntensity; 01986 FeII.fe33to7 = Fe2LevN[33-1][6].Emis->xIntensity; 01987 FeII.fe37to7 = Fe2LevN[37-1][6].Emis->xIntensity; 01988 FeII.fe39to8 = Fe2LevN[39-1][7].Emis->xIntensity; 01989 FeII.fe40to9 = Fe2LevN[40-1][8].Emis->xIntensity; 01990 FeII.fe37to6 = Fe2LevN[37-1][5].Emis->xIntensity; 01991 FeII.fe39to7 = Fe2LevN[39-1][6].Emis->xIntensity; 01992 FeII.fe40to8 = Fe2LevN[40-1][7].Emis->xIntensity; 01993 FeII.fe41to9 = Fe2LevN[41-1][8].Emis->xIntensity; 01994 FeII.fe39to6 = Fe2LevN[39-1][5].Emis->xIntensity; 01995 FeII.fe40to7 = Fe2LevN[40-1][6].Emis->xIntensity; 01996 FeII.fe41to8 = Fe2LevN[41-1][7].Emis->xIntensity; 01997 01998 FeII.fe42to6 = Fe2LevN[42-1][5].Emis->xIntensity; 01999 FeII.fe43to7 = Fe2LevN[43-1][6].Emis->xIntensity; 02000 FeII.fe42to7 = Fe2LevN[42-1][6].Emis->xIntensity; 02001 FeII.fe36to2 = Fe2LevN[36-1][1].Emis->xIntensity; 02002 FeII.fe36to3 = Fe2LevN[36-1][2].Emis->xIntensity; 02003 /*lint -e778 const expression eval to 0 */ 02004 FeII.fe32to1 = Fe2LevN[32-1][0].Emis->xIntensity; 02005 /*lint +e778 const expression eval to 0 */ 02006 FeII.fe33to2 = Fe2LevN[33-1][1].Emis->xIntensity; 02007 FeII.fe36to5 = Fe2LevN[36-1][4].Emis->xIntensity; 02008 FeII.fe32to2 = Fe2LevN[32-1][1].Emis->xIntensity; 02009 FeII.fe33to3 = Fe2LevN[33-1][2].Emis->xIntensity; 02010 FeII.fe30to3 = Fe2LevN[30-1][2].Emis->xIntensity; 02011 FeII.fe33to6 = Fe2LevN[33-1][5].Emis->xIntensity; 02012 FeII.fe24to2 = Fe2LevN[24-1][1].Emis->xIntensity; 02013 FeII.fe32to7 = Fe2LevN[32-1][6].Emis->xIntensity; 02014 FeII.fe35to8 = Fe2LevN[35-1][7].Emis->xIntensity; 02015 FeII.fe34to8 = Fe2LevN[34-1][7].Emis->xIntensity; 02016 FeII.fe27to6 = Fe2LevN[27-1][5].Emis->xIntensity; 02017 FeII.fe28to7 = Fe2LevN[28-1][6].Emis->xIntensity; 02018 FeII.fe30to8 = Fe2LevN[30-1][7].Emis->xIntensity; 02019 FeII.fe24to6 = Fe2LevN[24-1][5].Emis->xIntensity; 02020 FeII.fe29to8 = Fe2LevN[29-1][7].Emis->xIntensity; 02021 FeII.fe24to7 = Fe2LevN[24-1][6].Emis->xIntensity; 02022 FeII.fe22to7 = Fe2LevN[22-1][6].Emis->xIntensity; 02023 FeII.fe38to11 =Fe2LevN[38-1][11-1].Emis->xIntensity; 02024 FeII.fe19to8 = Fe2LevN[19-1][7].Emis->xIntensity; 02025 FeII.fe17to6 = Fe2LevN[17-1][5].Emis->xIntensity; 02026 FeII.fe18to7 = Fe2LevN[18-1][6].Emis->xIntensity; 02027 FeII.fe18to8 = Fe2LevN[18-1][7].Emis->xIntensity; 02028 FeII.fe17to7 = Fe2LevN[17-1][6].Emis->xIntensity; 02029 } 02030 else 02031 { 02032 FeII.fe25to6 = 0.; 02033 FeII.fe27to7 = 0.; 02034 FeII.fe28to8 = 0.; 02035 FeII.fe29to9 = 0.; 02036 FeII.fe32to6 = 0.; 02037 FeII.fe33to7 = 0.; 02038 FeII.fe37to7 = 0.; 02039 FeII.fe39to8 = 0.; 02040 FeII.fe40to9 = 0.; 02041 FeII.fe37to6 = 0.; 02042 FeII.fe39to7 = 0.; 02043 FeII.fe40to8 = 0.; 02044 FeII.fe41to9 = 0.; 02045 FeII.fe39to6 = 0.; 02046 FeII.fe40to7 = 0.; 02047 FeII.fe41to8 = 0.; 02048 02049 FeII.fe42to6 = 0.; 02050 FeII.fe43to7 = 0.; 02051 FeII.fe42to7 = 0.; 02052 FeII.fe36to2 = 0.; 02053 FeII.fe36to3 = 0.; 02054 FeII.fe32to1 = 0.; 02055 FeII.fe33to2 = 0.; 02056 FeII.fe36to5 = 0.; 02057 FeII.fe32to2 = 0.; 02058 FeII.fe33to3 = 0.; 02059 FeII.fe30to3 = 0.; 02060 FeII.fe33to6 = 0.; 02061 FeII.fe24to2 = 0.; 02062 FeII.fe32to7 = 0.; 02063 FeII.fe35to8 = 0.; 02064 FeII.fe34to8 = 0.; 02065 FeII.fe27to6 = 0.; 02066 FeII.fe28to7 = 0.; 02067 FeII.fe30to8 = 0.; 02068 FeII.fe24to6 = 0.; 02069 FeII.fe29to8 = 0.; 02070 FeII.fe24to7 = 0.; 02071 FeII.fe22to7 = 0.; 02072 FeII.fe38to11 =0.; 02073 FeII.fe19to8 = 0.; 02074 FeII.fe17to6 = 0.; 02075 FeII.fe18to7 = 0.; 02076 FeII.fe18to8 = 0.; 02077 FeII.fe17to7 = 0.; 02078 } 02079 02080 if( FeII.nFeIILevel > 80 ) 02081 { 02082 FeII.fe80to28 = Fe2LevN[80-1][28-1].Emis->xIntensity; 02083 } 02084 else 02085 { 02086 FeII.fe80to28 =0.; 02087 } 02088 02089 return; 02090 } 02091 /* 02092 *===================================================================== 02093 * punch line data for FeII atom */ 02094 void FeIIPunData( 02095 /* io unit for punch */ 02096 FILE* ioPUN , 02097 /* punch all levels if true, only subset if false */ 02098 bool lgDoAll ) 02099 { 02100 long int ipLo , ipHi; 02101 02102 DEBUG_ENTRY( "FeIIPunData()" ); 02103 02104 if( lgDoAll ) 02105 { 02106 fprintf( ioQQQ, 02107 " FeIIPunData ALL option not implemented yet 1\n" ); 02108 cdEXIT(EXIT_FAILURE); 02109 } 02110 else 02111 { 02112 long int nSkip=0, limit=MIN2(64, FeII.nFeIILevel); 02113 02114 /* false, only punch subset in above init */ 02115 /* first 64 do all lines */ 02116 for( ipHi=1; ipHi<limit; ++ipHi ) 02117 { 02118 for( ipLo=0; ipLo<ipHi; ++ipLo ) 02119 { 02120 fprintf(ioPUN,"%4li%4li ",ipLo,ipHi ); 02121 Punch1LineData( &Fe2LevN[ipHi][ipLo] , ioPUN, false); 02122 } 02123 } 02124 fprintf( ioPUN , "\n"); 02125 02126 if( limit==64 ) 02127 { 02128 /* higher than 64 only do real transitions - the majority of those above 02129 * n = 64 have no radiative transitions but still exist to hold collision, 02130 * energy, and other line data - they will have Aul == 1e-5 */ 02131 for( ipHi=limit; ipHi<FeII.nFeIILevel; ++ipHi ) 02132 { 02133 /*>>chng 06 jun 23, ipLo had ranged from limit up to ipHi, 02134 * so never printed lines from ipHi to ipLo<64 */ 02135 for( ipLo=0; ipLo<ipHi; ++ipLo ) 02136 { 02137 /*fprintf(ioPUN , "%li %li\n", ipHi , ipLo ); 02138 if( ipHi==70 && ipLo==0 ) 02139 Punch1LineData( &Fe2LevN[ipHi][ipLo] , ioPUN , false); */ 02140 /* ncs1 of 3 and Aul = 1e-5 indicate transition that is not optically 02141 * allowed with fake cs */ 02142 if( ncs1[ipHi][ipLo] == 3 && 02143 (fabs(Fe2LevN[ipHi][ipLo].Emis->Aul-1e-5) < 1e-8 ) ) 02144 { 02145 ++nSkip; 02146 } 02147 else 02148 { 02149 /* add one to ipLo and ipHi so that printed number is on atomic, 02150 * not C, scale */ 02151 fprintf(ioPUN,"%4li%4li ",ipLo+1,ipHi+1 ); 02152 Punch1LineData( &Fe2LevN[ipHi][ipLo] , ioPUN , false); 02153 } 02154 } 02155 } 02156 fprintf( ioPUN , " %li lines skiped\n",nSkip); 02157 } 02158 } 02159 02160 return; 02161 } 02162 02163 /* 02164 *===================================================================== 02165 */ 02166 void FeIIPunDepart( 02167 /* io unit for punch */ 02168 FILE* ioPUN , 02169 /* punch all levels if true, only subset if false */ 02170 bool lgDoAll ) 02171 { 02172 /* numer of levels with dep coef that we will punch out */ 02173 # define NLEVDEP 11 02174 /* these are the levels on the physical, not c, scale (count from 1) */ 02175 const int LevDep[NLEVDEP]={1,10,25,45,64,124,206,249,295,347,371}; 02176 long int i; 02177 static bool lgFIRST=true; 02178 02179 DEBUG_ENTRY( "FeIIPunDepart()" ); 02180 02181 /* on first call only, print levels that we will punch later */ 02182 if( lgFIRST && !lgDoAll ) 02183 { 02184 /* but all the rest do */ 02185 for( i=0; i<NLEVDEP; ++i ) 02186 { 02187 fprintf( ioPUN , "%i\t", LevDep[i] ); 02188 } 02189 fprintf( ioPUN , "\n"); 02190 lgFIRST = false; 02191 } 02192 02193 if( lgDoAll ) 02194 { 02195 /* true, punch all levels, one per line */ 02196 for( i=1; i<=FeII.nFeIILevel; ++i ) 02197 { 02198 FeIIPun1Depart( ioPUN , i ); 02199 fprintf( ioPUN , "\n"); 02200 } 02201 } 02202 else 02203 { 02204 /* false, only punch subset in above init */ 02205 for( i=0; i<NLEVDEP; ++i ) 02206 { 02207 FeIIPun1Depart( ioPUN , LevDep[i] ); 02208 fprintf( ioPUN , "\t"); 02209 } 02210 fprintf( ioPUN , "\n"); 02211 } 02212 02213 return; 02214 } 02215 02216 /* 02217 *===================================================================== 02218 */ 02219 void FeIIPun1Depart( 02220 /* the io unit where the print should be directed */ 02221 FILE * ioPUN , 02222 /* the physical (not c) number of the level */ 02223 long int nPUN ) 02224 { 02225 02226 DEBUG_ENTRY( "FeIIPun1Depart()" ); 02227 02228 ASSERT( nPUN > 0 ); 02229 /* need real assert to keep lint happy */ 02230 assert( ioPUN != NULL ); 02231 02232 /* print the level departure coef on ioPUN if the level was computed, 02233 * print a zero if it was not */ 02234 if( nPUN <= FeII.nFeIILevel ) 02235 { 02236 fprintf( ioPUN , "%e ",Fe2DepCoef[nPUN-1] ); 02237 } 02238 else 02239 { 02240 fprintf( ioPUN , "%e ",0. ); 02241 } 02242 02243 return; 02244 } 02245 02246 /* 02247 *===================================================================== 02248 */ 02249 void FeIIReset(void) 02250 { 02251 02252 DEBUG_ENTRY( "FeIIReset()" ); 02253 02254 /* space has been allocated, so reset number of levels to whatever it was */ 02255 FeII.nFeIILevel = FeII.nFeIILevelAlloc; 02256 02257 return; 02258 } 02259 02260 /* 02261 *===================================================================== 02262 */ 02263 /*FeIIZero initialize some variables, called by zero one time before commands parsed*/ 02264 void FeIIZero(void) 02265 { 02266 02267 DEBUG_ENTRY( "FeIIZero()" ); 02268 02269 /* heating, cooling, and deriv wrt temperature */ 02270 FeII.Fe2_large_cool = 0.; 02271 FeII.ddT_Fe2_large_cool = 0.; 02272 FeII.Fe2_large_heat = 0.; 02273 02274 /* flag saying that lya pumping of FeII in large atom is turned on */ 02275 FeII.lgLyaPumpOn = true; 02276 02277 /*FeII. lgFeIION = false;*/ 02278 /* >>chng 05 nov 29, test effects of always including FeII ground term with large atom 02279 * but if ground term only is done, still also do simple UV approximation */ 02280 /*FeII. lgFeIION = true;*/ 02281 /* will not compute large FeII atom */ 02282 FeII.lgFeIILargeOn = false; 02283 02284 /* energy range of large FeII atom is zero to infinity */ 02285 FeII.fe2ener[0] = 0.; 02286 FeII.fe2ener[1] = 1e8; 02287 02288 /* pre mar 01, these had both been 1, ipPRD */ 02289 /* resonance lines, ipCRD is -1 */ 02290 FeII.ipRedisFcnResonance = ipCRD; 02291 /* subordinate lines, ipCRDW is 2 */ 02292 FeII.ipRedisFcnSubordinate = ipCRDW; 02293 02294 /* set zero for the threshold of weakest large FeII atom line to punch */ 02295 FeII.fe2thresh = 0.; 02296 02297 /* normally do not constantly reevaluate the atom, set true with 02298 * SLOW key on atom FeII command */ 02299 FeII.lgSlow = false; 02300 02301 /* option to print each call to FeIILevelPops, set with print option on atom FeII */ 02302 FeII.lgPrint = false; 02303 02304 /* option to only simulate calls to FeIILevelPops */ 02305 FeII.lgSimulate = false; 02306 02307 /* set number of levels for the atom */ 02308 /* number of levels for the large FeII atom, changed with the atom FeII levels command */ 02309 if( lgFeIIMalloc ) 02310 { 02311 /* space has been allocated, so reset number of levels to whatever it was */ 02312 FeII.nFeIILevel = FeII.nFeIILevelAlloc; 02313 } 02314 else 02315 { 02316 /* space not allocated yet, set to largest possible number of levels */ 02317 FeII.nFeIILevel = NFE2LEVN; 02318 /* >>chng 05 nov 29, test effects of always including FeII ground term with large atom 02319 * but if ground term only is done, still also do simple UV approximation 02320 * set this to only ground term, - will reset to NFE2LEVN when atom FeII parsed if levels not set */ 02321 FeII.nFeIILevel = 16; 02322 } 02323 02324 /* lower and upper wavelength bounds, Angstroms, for the FeII continuum */ 02325 FeII.fe2con_wl1 = 1000.; 02326 FeII.fe2con_wl2 = 7000.; 02327 FeII.nfe2con = 1000; 02328 02329 return; 02330 } 02331 02332 /*FeIIPunPop - punch level populations */ 02333 void FeIIPunPop( 02334 /* io unit for punch */ 02335 FILE* ioPUN , 02336 /* punch range of levels if true, only selected subset if false */ 02337 bool lgPunchRange , 02338 /* if ipPunchRange is true, this is lower bound of range on C scale */ 02339 long int ipRangeLo , 02340 /* if ipPunchRange is true, this is upper bound of range on C scale */ 02341 long int ipRangeHi , 02342 /* flag saying whether to punch density (cm-3, true) or relative population (flase) */ 02343 bool lgPunchDensity ) 02344 { 02345 /* numer of levels with dep coef that we will punch out */ 02346 # define NLEVPOP 11 02347 /* these are the levels on the physical, not c, scale (count from 1) */ 02348 const int LevPop[NLEVPOP]={1,10,25,45,64,124,206,249,295,347,371}; 02349 long int i; 02350 static bool lgFIRST=true; 02351 realnum denominator = 1.f; 02352 02353 DEBUG_ENTRY( "FeIIPunPop()" ); 02354 02355 /* this implements the relative option on punch FeII populations command */ 02356 if( !lgPunchDensity ) 02357 denominator = SDIV( dense.xIonDense[ipIRON][1] ); 02358 02359 /* on first call only, print levels that we will punch later, 02360 * but only if we will only punch selected levels*/ 02361 if( lgFIRST && !lgPunchRange ) 02362 { 02363 /* but all the rest do */ 02364 for( i=0; i<NLEVPOP; ++i ) 02365 { 02366 /* indices for particular levels */ 02367 fprintf( ioPUN , "%i\t", LevPop[i] ); 02368 } 02369 fprintf( ioPUN , "\n"); 02370 lgFIRST = false; 02371 } 02372 02373 if( lgPunchRange ) 02374 { 02375 /* confirm sane level indices */ 02376 ASSERT( ipRangeLo>=0 && ipRangeLo<ipRangeHi && ipRangeHi<=FeII.nFeIILevel ); 02377 02378 /* true, punch all levels across line, 02379 * both call with physical level so that list is physical */ 02380 for( i=ipRangeLo; i<ipRangeHi; ++i ) 02381 { 02382 /* routine takes levels on physics scale */ 02383 fprintf( ioPUN , "%.3e\t", Fe2LevelPop[i]/denominator ); 02384 } 02385 fprintf( ioPUN , "\n"); 02386 } 02387 else 02388 { 02389 /* false, only punch subset in above init, 02390 * both call with physical level so that list is physical */ 02391 for( i=0; i<NLEVPOP; ++i ) 02392 { 02393 fprintf( ioPUN , "%.3e\t", Fe2LevelPop[LevPop[i]-1]/denominator ); 02394 } 02395 fprintf( ioPUN , "\n"); 02396 } 02397 02398 return; 02399 } 02400 02401 /* 02402 *===================================================================== 02403 */ 02404 #if 0 02405 void FeIIPun1Pop( 02406 /* the io unit where the print should be directed */ 02407 FILE * ioPUN , 02408 /* the physical (not c) number of the level */ 02409 long int nPUN ) 02410 { 02411 DEBUG_ENTRY( "FeIIPun1Pop()" ); 02412 02413 ASSERT( nPUN > 0 ); 02414 /* need real assert to keep lint happy */ 02415 assert( ioPUN != NULL ); 02416 02417 /* print the level population on ioPUN if the level was computed, 02418 * print a zero if it was not, 02419 * note that nPUN is on physical scale, so test is <=*/ 02420 if( nPUN <= FeII.nFeIILevel ) 02421 { 02422 fprintf( ioPUN , "%e ",Fe2LevelPop[nPUN-1] ); 02423 } 02424 else 02425 { 02426 fprintf( ioPUN , "%e ",0. ); 02427 } 02428 02429 return; 02430 } 02431 #endif 02432 02433 /* 02434 *===================================================================== 02435 */ 02436 STATIC int FeIIBandsCreate( 02437 /* chFile is optional filename, if void then use default bands, 02438 * if not void then use file specified, 02439 * return value is 0 for success, 1 for failure */ 02440 const char chFile[] ) 02441 { 02442 02443 char chLine[FILENAME_PATH_LENGTH_2]; 02444 const char* chFile1; 02445 FILE *ioDATA; 02446 02447 bool lgEOL; 02448 long int i,k; 02449 02450 /* keep track of whether we have been called - want to be 02451 * called a total of one time */ 02452 static bool lgCalled=false; 02453 02454 DEBUG_ENTRY( "FeIIBandsCreate()" ); 02455 02456 /* return previous number of bands if this is second or later call*/ 02457 if( lgCalled ) 02458 { 02459 /* success */ 02460 return 0; 02461 } 02462 lgCalled = true; 02463 02464 /* use default filename if void string, else use file specified */ 02465 chFile1 = ( strlen(chFile) == 0 ) ? "bands_Fe2.dat" : chFile; 02466 02467 /* get FeII band data */ 02468 if( trace.lgTrace ) 02469 { 02470 fprintf( ioQQQ, " FeIICreate opening %s:", chFile1 ); 02471 } 02472 02473 ioDATA = open_data( chFile1, "r" ); 02474 02475 /* now count how many bands are in the file */ 02476 nFeIIBands = 0; 02477 02478 /* first line is a version number and does not count */ 02479 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 02480 { 02481 fprintf( ioQQQ, " FeIICreate could not read first line of %s.\n", chFile1 ); 02482 return 1; 02483 } 02484 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL ) 02485 { 02486 /* we want to count the lines that do not start with # 02487 * since these contain data */ 02488 if( chLine[0] != '#') 02489 ++nFeIIBands; 02490 } 02491 02492 /* now rewind the file so we can read it a second time*/ 02493 if( fseek( ioDATA , 0 , SEEK_SET ) != 0 ) 02494 { 02495 fprintf( ioQQQ, " FeIICreate could not rewind %s.\n", chFile1 ); 02496 return 1; 02497 } 02498 02499 FeII_Bands = (realnum **)MALLOC(sizeof(realnum *)*(unsigned)(nFeIIBands) ); 02500 02501 /* now make second dim, id wavelength, and lower and upper bounds */ 02502 for( i=0; i<nFeIIBands; ++i ) 02503 { 02504 FeII_Bands[i] = (realnum *)MALLOC(sizeof(realnum)*(unsigned)(3) ); 02505 } 02506 02507 /* first line is a version number - now confirm that it is valid */ 02508 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 02509 { 02510 fprintf( ioQQQ, " FeIICreate could not read first line of %s.\n", chFile1 ); 02511 return 1; 02512 } 02513 i = 1; 02514 if( ( (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL) != 99 ) || 02515 ( (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL) != 12 ) || 02516 ( (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL) != 1 ) ) 02517 { 02518 fprintf( ioQQQ, 02519 " FeIICreate: the version of %s is not the current version.\n", chFile1 ); 02520 return 1; 02521 } 02522 02523 /* now read in data again, but save it this time */ 02524 k = 0; 02525 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL ) 02526 { 02527 /* we want to count the lines that do not start with # 02528 * since these contain data */ 02529 if( chLine[0] != '#') 02530 { 02531 i = 1; 02532 FeII_Bands[k][0] = (realnum)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL); 02533 if( lgEOL ) 02534 { 02535 fprintf( ioQQQ, " There should have been a number on this band line 1. Sorry.\n" ); 02536 fprintf( ioQQQ, "string==%s==\n" ,chLine ); 02537 return 1; 02538 } 02539 FeII_Bands[k][1] = (realnum)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL); 02540 if( lgEOL ) 02541 { 02542 fprintf( ioQQQ, " There should have been a number on this band line 2. Sorry.\n" ); 02543 fprintf( ioQQQ, "string==%s==\n" ,chLine ); 02544 return 1; 02545 } 02546 FeII_Bands[k][2] = (realnum)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL); 02547 if( lgEOL ) 02548 { 02549 fprintf( ioQQQ, " There should have been a number on this band line 3. Sorry.\n" ); 02550 fprintf( ioQQQ, "string==%s==\n" ,chLine ); 02551 return 1; 02552 } 02553 /*fprintf(ioQQQ, 02554 " band data %f %f %f \n", FeII_Bands[k][0],FeII_Bands[k][1],FeII_Bands[k][2]);*/ 02555 ++k; 02556 } 02557 } 02558 /* now validate this incoming data */ 02559 for( i=0; i<nFeIIBands; ++i ) 02560 { 02561 /* make sure all are positive */ 02562 if( FeII_Bands[i][0] <=0. || FeII_Bands[i][1] <=0. || FeII_Bands[i][2] <=0. ) 02563 { 02564 fprintf( ioQQQ, " FeII band %li has none positive entry.\n",i ); 02565 return 1; 02566 } 02567 /* make sure bands bounds are in correct order, shorter - longer wavelength*/ 02568 if( FeII_Bands[i][1] >= FeII_Bands[i][2] ) 02569 { 02570 fprintf( ioQQQ, " FeII band %li has improper bounds.\n" ,i); 02571 return 1; 02572 } 02573 02574 } 02575 02576 fclose(ioDATA); 02577 02578 /* return success */ 02579 return 0; 02580 } 02581 /* 02582 *===================================================================== 02583 */ 02584 void AssertFeIIDep( double *pred , double *BigError , double *StdDev ) 02585 { 02586 long int n; 02587 double arg , 02588 error , 02589 sum2; 02590 02591 DEBUG_ENTRY( "AssertFeIIDep()" ); 02592 02593 if( FeII.lgSimulate ) 02594 { 02595 02596 *pred = 0.; 02597 *BigError = 0.; 02598 *StdDev = 0.; 02599 02600 return; 02601 } 02602 02603 /* sanity check */ 02604 ASSERT( FeII.nFeIILevel > 0 ); 02605 02606 /* find sum of deviation of departure coef from unity */ 02607 *BigError = 0.; 02608 *pred = 0.; 02609 sum2 = 0; 02610 for( n=0; n<FeII.nFeIILevel; ++n ) 02611 { 02612 /* get mean */ 02613 *pred += Fe2DepCoef[n]; 02614 02615 /* error is the largest deviation from unity for any single level*/ 02616 error = fabs( Fe2DepCoef[n] -1. ); 02617 /* remember biggest deviation */ 02618 *BigError = MAX2( *BigError , error ); 02619 02620 /* get standard deviation */ 02621 sum2 += POW2( Fe2DepCoef[n] ); 02622 } 02623 02624 /* now get standard deviation */ 02625 arg = sum2 - POW2( *pred ) / (double)FeII.nFeIILevel; 02626 ASSERT( (arg >= 0.) ); 02627 *StdDev = sqrt( arg / (double)(FeII.nFeIILevel - 1.) ); 02628 02629 /* this is average value, should be unity */ 02630 *pred /= (double)(FeII.nFeIILevel); 02631 02632 return; 02633 } 02634 02635 /* 02636 *===================================================================== 02637 */ 02638 /* do ots rates for FeII, called by RT_OTS */ 02639 void FeII_OTS( void ) 02640 { 02641 long int ipLo , 02642 ipHi; 02643 02644 DEBUG_ENTRY( "FeII_OTS()" ); 02645 02646 for( ipLo=0; ipLo < (FeII.nFeIILevel - 1); ipLo++ ) 02647 { 02648 for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevel; ipHi++ ) 02649 { 02650 /* >>chng 02 feb 11, skip bogus transitions */ 02651 if( Fe2LevN[ipHi][ipLo].ipCont < 1) 02652 continue; 02653 02654 /* ots rates, the destp prob was set in hydropesc */ 02655 Fe2LevN[ipHi][ipLo].Emis->ots = 02656 Fe2LevN[ipHi][ipLo].Emis->Aul* 02657 Fe2LevN[ipHi][ipLo].Hi->Pop* 02658 Fe2LevN[ipHi][ipLo].Emis->Pdest; 02659 02660 ASSERT( Fe2LevN[ipHi][ipLo].Emis->ots >= 0. ); 02661 02662 /* finally dump the ots rate into the stack */ 02663 RT_OTS_AddLine(Fe2LevN[ipHi][ipLo].Emis->ots, 02664 Fe2LevN[ipHi][ipLo].ipCont ); 02665 } 02666 } 02667 02668 return; 02669 } 02670 02671 /* 02672 *===================================================================== 02673 */ 02674 /*FeII_RT_Out - do outward rates for FeII, 02675 * called by RT_diffuse, which is called by cloudy */ 02676 void FeII_RT_Out(void) 02677 { 02678 long int ipLo , ipHi; 02679 02680 DEBUG_ENTRY( "FeIIRTOut()" ); 02681 02682 /* only do this if Fe+ exists */ 02683 if( dense.xIonDense[ipIRON][1] > 0. ) 02684 { 02685 /* outward line photons */ 02686 for( ipLo=0; ipLo < (FeII.nFeIILevel - 1); ipLo++ ) 02687 { 02688 for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevel; ipHi++ ) 02689 { 02690 /* >>chng 02 feb 11, skip bogus transitions */ 02691 if( Fe2LevN[ipHi][ipLo].ipCont < 1) 02692 continue; 02693 /* >>chng 03 sep 09, change to outline, ouutlin is 02694 * not exactly the same in the two - tmn missing in outline */ 02695 outline( &Fe2LevN[ipHi][ipLo] ); 02696 02697 } 02698 } 02699 } 02700 02701 return; 02702 } 02703 02704 /* 02705 *===================================================================== 02706 */ 02707 STATIC void FeIIContCreate( 02708 /* wavelength of low-lambda end */ 02709 double xLamLow , 02710 /* wavelength of high end */ 02711 double xLamHigh , 02712 /* number of cells to break this into */ 02713 long int ncell ) 02714 { 02715 02716 double step , wl1; 02717 long int i; 02718 02719 /* keep track of whether we have been called - want to be 02720 * called a total of one time */ 02721 static bool lgCalled=false; 02722 02723 DEBUG_ENTRY( "FeIIContCreate()" ); 02724 02725 /* return previous number of bands if this is second or later call*/ 02726 if( lgCalled ) 02727 { 02728 /* return value of number of bands, may be used by calling program*/ 02729 return; 02730 } 02731 lgCalled = true; 02732 02733 /* how many cells will be needed to go from xLamLow to xLamHigh in ncell steps */ 02734 nFeIIConBins = ncell; 02735 02736 FeII_Cont = (realnum **)MALLOC(sizeof(realnum *)*(unsigned)(nFeIIConBins) ); 02737 02738 /* now make second dim, id wavelength, and lower and upper bounds */ 02739 for( i=0; i<nFeIIConBins; ++i ) 02740 { 02741 FeII_Cont[i] = (realnum *)MALLOC(sizeof(realnum)*(unsigned)(3) ); 02742 } 02743 02744 step = log10( xLamHigh/xLamLow)/ncell; 02745 wl1 = log10( xLamLow); 02746 FeII_Cont[0][1] = (realnum)pow(10. ,wl1); 02747 for( i=1; i<ncell; ++i) 02748 { 02749 /* lower bound of cell in Angstroms */ 02750 FeII_Cont[i][1] = (realnum)pow(10. , ( wl1 + i*step ) ); 02751 } 02752 02753 for( i=0; i<(ncell-1); ++i) 02754 { 02755 /* upper bound of cell in Angstroms */ 02756 FeII_Cont[i][2] = FeII_Cont[i+1][1]; 02757 } 02758 FeII_Cont[ncell-1][2] = (realnum)(pow(10., log10(FeII_Cont[ncell-2][2]) + step) ); 02759 02760 for( i=0; i<ncell; ++i) 02761 { 02762 /* wavelength (A) as it will appear in the printout */ 02763 FeII_Cont[i][0] = ( FeII_Cont[i][1] + FeII_Cont[i][2] ) /2.f; 02764 } 02765 { 02766 enum {DEBUG_LOC=false}; 02767 if( DEBUG_LOC ) 02768 { 02769 FILE *ioDEBUG; 02770 ioDEBUG = fopen( "fe2cont.txt", "w" ); 02771 if( ioDEBUG==NULL ) 02772 { 02773 fprintf( ioQQQ," fe2con could not open fe2cont.txt for writing\n"); 02774 cdEXIT(EXIT_FAILURE); 02775 } 02776 for( i=0; i<ncell; ++i) 02777 { 02778 /* wavelength as it will appear in the printout */ 02779 fprintf(ioDEBUG,"%.1f\t%.1f\t%.1f\n", 02780 FeII_Cont[i][0] , FeII_Cont[i][1] , FeII_Cont[i][2] ); 02781 } 02782 fclose( ioDEBUG); 02783 } 02784 } 02785 02786 return; 02787 } 02788 02789 /* FeIIOvrLap handle overlapping FeII lines */ 02790 STATIC void FeIIOvrLap(void) 02791 { 02792 02793 DEBUG_ENTRY( "FeIIOvrLap()" ); 02794 02795 return; 02796 } 02797 02798 /*ParseAtomFeII parse the atom FeII command */ 02799 void ParseAtomFeII(char *chCard ) 02800 { 02801 long int i; 02802 bool lgEOL=false; 02803 02804 DEBUG_ENTRY( "ParseAtomFeII()" ); 02805 02806 /* turn on the large verner atom */ 02807 FeII.lgFeIILargeOn = true; 02808 /* >>chng 05 nov 29, reset number of levels when called, so increased above default of 16 */ 02809 if( lgFeIIMalloc ) 02810 { 02811 /* space has been allocated, so reset number of levels to whatever it was */ 02812 FeII.nFeIILevel = FeII.nFeIILevelAlloc; 02813 } 02814 else 02815 { 02816 /* space not allocated yet, set to largest possible number of levels */ 02817 FeII.nFeIILevel = NFE2LEVN; 02818 } 02819 02820 /* levels keyword is to adjust number of levels. But this only has effect 02821 * BEFORE space is allocated for the FeII arrays */ 02822 if( nMatch("LEVE",chCard) ) 02823 { 02824 /* do option only if space not yet allocated */ 02825 if( !lgFeIIMalloc ) 02826 { 02827 i = 5; 02828 /* number of levels for hydrogen at, 2s is this plus one */ 02829 FeII.nFeIILevel = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 02830 if( lgEOL ) 02831 { 02832 /* hit eol with no number - this is a problem */ 02833 NoNumb(chCard); 02834 } 02835 if( FeII.nFeIILevel <16 ) 02836 { 02837 fprintf( ioQQQ, " This would be too few levels, must have at least 16.\n" ); 02838 cdEXIT(EXIT_FAILURE); 02839 } 02840 else if( FeII.nFeIILevel > NFE2LEVN ) 02841 { 02842 fprintf( ioQQQ, " This would be too many levels.\n" ); 02843 cdEXIT(EXIT_FAILURE); 02844 } 02845 } 02846 } 02847 02848 /* slow keyword means do not try to avoid evaluating atom */ 02849 else if( nMatch("SLOW",chCard) ) 02850 { 02851 FeII.lgSlow = true; 02852 } 02853 02854 /* redistribution keyword changes form of redistribution function */ 02855 else if( nMatch("REDI",chCard) ) 02856 { 02857 int ipRedis=0; 02858 /* there are three functions, PRD_, CRD_, and CRDW, 02859 * representing partial redistribution, 02860 * complete redistribution (doppler core only, no wings) 02861 * and complete with wings */ 02862 /* partial redistribution */ 02863 if( nMatch(" PRD",chCard) ) 02864 { 02865 ipRedis = ipPRD; 02866 } 02867 /* complete redistribution */ 02868 else if( nMatch(" CRD",chCard) ) 02869 { 02870 ipRedis = ipCRD; 02871 } 02872 /* complete redistribution with wings */ 02873 else if( nMatch("CRDW",chCard) ) 02874 { 02875 ipRedis = ipCRDW; 02876 } 02877 02878 /* if not SHOW option (handled below) then we have a problem */ 02879 else if( !nMatch("SHOW",chCard) ) 02880 { 02881 fprintf(ioQQQ," There should have been a second keyword on this command.\n"); 02882 fprintf(ioQQQ," Options are _PRD, _CRD, CRDW (_ is space). Sorry.\n"); 02883 cdEXIT(EXIT_FAILURE); 02884 } 02885 02886 /* resonance lines */ 02887 if( nMatch("RESO",chCard) ) 02888 { 02889 FeII.ipRedisFcnResonance = ipRedis; 02890 } 02891 /* subordinate lines */ 02892 else if( nMatch("SUBO",chCard) ) 02893 { 02894 FeII.ipRedisFcnSubordinate = ipRedis; 02895 } 02896 /* the show option, say what we are assuming */ 02897 else if( nMatch("SHOW",chCard) ) 02898 { 02899 fprintf(ioQQQ," FeII resonance lines are "); 02900 if( FeII.ipRedisFcnResonance ==ipCRDW ) 02901 { 02902 fprintf(ioQQQ,"complete redistribution with wings\n"); 02903 } 02904 else if( FeII.ipRedisFcnResonance ==ipCRD ) 02905 { 02906 fprintf(ioQQQ,"complete redistribution with core only.\n"); 02907 } 02908 else if( FeII.ipRedisFcnResonance ==ipPRD ) 02909 { 02910 fprintf(ioQQQ,"partial redistribution.\n"); 02911 } 02912 else 02913 { 02914 fprintf(ioQQQ," PROBLEM Impossible value for ipRedisFcnResonance.\n"); 02915 TotalInsanity(); 02916 } 02917 02918 fprintf(ioQQQ," FeII subordinate lines are "); 02919 if( FeII.ipRedisFcnSubordinate ==ipCRDW ) 02920 { 02921 fprintf(ioQQQ,"complete redistribution with wings\n"); 02922 } 02923 else if( FeII.ipRedisFcnSubordinate ==ipCRD ) 02924 { 02925 fprintf(ioQQQ,"complete redistribution with core only.\n"); 02926 } 02927 else if( FeII.ipRedisFcnSubordinate ==ipPRD ) 02928 { 02929 fprintf(ioQQQ,"partial redistribution.\n"); 02930 } 02931 else 02932 { 02933 fprintf(ioQQQ," PROBLEM Impossible value for ipRedisFcnSubordinate.\n"); 02934 TotalInsanity(); 02935 } 02936 } 02937 else 02938 { 02939 fprintf(ioQQQ," here should have been a second keyword on this command.\n"); 02940 fprintf(ioQQQ," Options are RESONANCE, SUBORDINATE. Sorry.\n"); 02941 cdEXIT(EXIT_FAILURE); 02942 } 02943 } 02944 02945 /* trace keyword - print info for each call to FeIILevelPops */ 02946 else if( nMatch("TRAC",chCard) ) 02947 { 02948 FeII.lgPrint = true; 02949 } 02950 02951 /* only simulate the FeII atom, do not actually call it */ 02952 else if( nMatch("SIMU",chCard) ) 02953 { 02954 /* option to only simulate calls to FeIILevelPops */ 02955 FeII.lgSimulate = true; 02956 } 02957 02958 /* only simulate the FeII atom, do not actually call it */ 02959 else if( nMatch("CONT",chCard) ) 02960 { 02961 /* the continuum output with the punch FeII continuum command */ 02962 i=5; 02963 02964 /* number of levels for hydrogen at, 2s is this plus one */ 02965 FeII.fe2con_wl1 = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 02966 FeII.fe2con_wl2 = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 02967 FeII.nfe2con = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 02968 if( lgEOL ) 02969 { 02970 fprintf(ioQQQ," there are three numbers on the FeII continuum command, start and end wavelengths, and number of intervals.\n"); 02971 /* hit eol with no number - this is a problem */ 02972 NoNumb(chCard); 02973 } 02974 /* check that all are positive */ 02975 if( FeII.fe2con_wl1<=0. || FeII.fe2con_wl2<=0. || FeII.nfe2con<= 0 ) 02976 { 02977 fprintf(ioQQQ," there are three numbers on the FeII continuum command, start and end wavelengths, and number of intervals.\n"); 02978 fprintf(ioQQQ," all three must be greater than zero, sorry.\n"); 02979 cdEXIT(EXIT_FAILURE); 02980 } 02981 /* make sure that wl1 is less than wl2 */ 02982 if( FeII.fe2con_wl1>= FeII.fe2con_wl2 ) 02983 { 02984 fprintf(ioQQQ," there are three numbers on the FeII continuum command, start and end wavelengths, and number of intervals.\n"); 02985 fprintf(ioQQQ," the second wavelength must be greater than the first, sorry.\n"); 02986 cdEXIT(EXIT_FAILURE); 02987 } 02988 } 02989 /* note no fall-through error since routine can be called with no options, 02990 * to turn on the large atom */ 02991 02992 return; 02993 } 02994 02995 void PunFeII( FILE * io ) 02996 { 02997 long int n, ipHi; 02998 02999 DEBUG_ENTRY( "PunFeII()" ); 03000 03001 for( n=0; n<FeII.nFeIILevel-1; ++n) 03002 { 03003 for( ipHi=n+1; ipHi<FeII.nFeIILevel; ++ipHi ) 03004 { 03005 if( Fe2LevN[ipHi][n].ipCont > 0 ) 03006 fprintf(io,"%li\t%li\t%.2e\n", n , ipHi , 03007 Fe2LevN[ipHi][n].Emis->TauIn ); 03008 } 03009 } 03010 03011 return; 03012 } 03013 03014 /* include FeII lines in punched optical depths, etc, called from PunchLineStuff */ 03015 void FeIIPunchLineStuff( FILE * io , realnum xLimit , long index) 03016 { 03017 long int n, ipHi; 03018 03019 DEBUG_ENTRY( "FeIIPunchLineStuff()" ); 03020 03021 for( n=0; n<FeII.nFeIILevel-1; ++n) 03022 { 03023 for( ipHi=n+1; ipHi<FeII.nFeIILevel; ++ipHi ) 03024 { 03025 pun1Line( &Fe2LevN[ipHi][n] , io , xLimit , index , 1. ); 03026 } 03027 } 03028 03029 return; 03030 } 03031 03032 /* rad pre due to FeII lines called in PresTotCurrent*/ 03033 double FeIIRadPress(void) 03034 { 03035 03036 /* will be used to check on size of opacity, was capped at this value */ 03037 realnum smallfloat=SMALLFLOAT*10.f; 03038 double press, 03039 RadPres1; 03040 # undef DEBUGFE 03041 # ifdef DEBUGFE 03042 double RadMax; 03043 long ipLoMax , ipHiMax; 03044 # endif 03045 long int ipLo, ipHi; 03046 03047 DEBUG_ENTRY( "FeIIRadPress()" ); 03048 03049 press = 0.; 03050 # ifdef DEBUGFE 03051 RadMax = 0; 03052 ipLoMax = -1; 03053 ipHiMax = -1; 03054 # endif 03055 /* loop over all levels to find radiation pressure */ 03056 for( ipHi=1; ipHi<FeII.nFeIILevel; ++ipHi ) 03057 { 03058 for( ipLo=0; ipLo<ipHi; ++ipLo) 03059 { 03060 /* >>chng 04 aug 27, add test on ipCont for bogus line */ 03061 if( Fe2LevN[ipHi][ipLo].ipCont > 0 && 03062 Fe2LevN[ipHi][ipLo].Hi->Pop > 1e-30 ) 03063 { 03064 if( Fe2LevN[ipHi][ipLo].Hi->Pop > smallfloat && 03065 Fe2LevN[ipHi][ipLo].Emis->PopOpc > smallfloat ) 03066 { 03067 RadPres1 = PressureRadiationLine( &Fe2LevN[ipHi][ipLo], 1.0 ); 03068 03069 # ifdef DEBUGFE 03070 if( RadPres1 > RadMax ) 03071 { 03072 RadMax = RadPres1; 03073 ipLoMax = ipLo; 03074 ipHiMax = ipHi; 03075 } 03076 # endif 03077 press += RadPres1; 03078 } 03079 } 03080 } 03081 } 03082 03083 # ifdef DEBUGFE 03084 /* option to print radiation pressure */ 03085 if( iteration > 1 || nzone > 1558 ) 03086 { 03087 fprintf(ioQQQ,"DEBUG FeIIRadPress finds P(FeII) = %.2e %.2e %li %li width %.2e\n", 03088 press , 03089 RadMax , 03090 ipLoMax , 03091 ipHiMax , 03092 RT_LineWidth(&Fe2LevN[9][0]) ); 03093 DumpLine( &Fe2LevN[9][0] ); 03094 } 03095 # endif 03096 # undef DEBUGFE 03097 03098 return press; 03099 } 03100 03101 #if 0 03102 /* internal energy of FeII called in PresTotCurrent */ 03103 double FeII_InterEnergy(void) 03104 { 03105 double energy; 03106 long int n, ipHi; 03107 03108 DEBUG_ENTRY( "FeII_InterEnergy()" ); 03109 03110 broken(); /* the code below contains serious bugs! It is supposed to implement a loop 03111 * over quantum states in order to evaluate the potential energy stored in 03112 * excited states of Fe II. The code below however implements loops over all 03113 * combinations of upper and lower levels! */ 03114 03115 energy = 0.; 03116 for( n=0; n<FeII.nFeIILevel-1; ++n) 03117 { 03118 for( ipHi=n+1; ipHi<FeII.nFeIILevel; ++ipHi ) 03119 { 03120 if( Fe2LevN[ipHi][n].Hi->Pop > 1e-30 ) 03121 { 03122 energy += 03123 Fe2LevN[ipHi][n].Hi->Pop* Fe2LevN[ipHi][n].EnergyErg; 03124 } 03125 } 03126 } 03127 03128 return energy; 03129 } 03130 #endif 03131 03132 /* HP cc cannot compile this routine with any optimization */ 03133 #if defined(__HP_aCC) 03134 # pragma OPT_LEVEL 1 03135 #endif 03136 /*FeIILyaPump find rate of Lya excitation of the FeII atom */ 03137 STATIC void FeIILyaPump(void) 03138 { 03139 03140 long int ipLo , 03141 ipHi; 03142 double EnerLyaProf2, 03143 EnerLyaProf3, 03144 EnergyWN, 03145 Gup_ov_Glo, 03146 PhotOccNum_at_nu, 03147 PumpRate, 03148 de, 03149 FeIILineWidthHz, 03150 taux; 03151 03152 DEBUG_ENTRY( "FeIILyaPump()" ); 03153 03154 /* lgLyaPumpOn is false if no Lya pumping, with no FeII command */ 03155 /* get rates FeII atom is pumped */ 03156 03157 /* elsewhere in this file the dest prob hydro.dstfe2lya is defined from 03158 * quantites derived here, and the resulting populations */ 03159 if( FeII.lgLyaPumpOn ) 03160 { 03161 03162 /*************trapeze form La profile:de,EnerLyaProf1,EnerLyaProf2,EnerLyaProf3,EnerLyaProf4************************* 03163 * */ 03164 /* width of Lya in cm^-1 */ 03165 /* HLineWidth has units of cm/s, as was evaluated in PresTotCurrent */ 03166 /* the factor is 1/2 of E(Lya, cm^-1_/c */ 03167 de = 1.37194e-06*hydro.HLineWidth*2.0/3.0; 03168 /* 82259 is energy of Lya in wavenumbers, so these are the form of the trapezoid */ 03169 EnerLyaProf1 = 82259.0 - de*2.0; 03170 EnerLyaProf2 = 82259.0 - de; 03171 EnerLyaProf3 = 82259.0 + de; 03172 EnerLyaProf4 = 82259.0 + de*2.0; 03173 03174 /* find Lya photon occupation number */ 03175 if( Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Hi->Pop > SMALLFLOAT ) 03176 { 03177 /* This is the photon occupation number at the Lya line center */ 03178 PhotOccNumLyaCenter = 03179 MAX2(0.,1.0- Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->Pelec_esc - 03180 Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->Pesc)/ 03181 (StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop/StatesElem[ipH_LIKE][ipHYDROGEN][ipH2p].Pop*3. - 1.0); 03182 } 03183 else 03184 { 03185 /* lya excitation temperature not available */ 03186 PhotOccNumLyaCenter = 0.; 03187 } 03188 } 03189 else 03190 { 03191 PhotOccNumLyaCenter = 0.; 03192 de = 0.; 03193 EnerLyaProf1 = 0.; 03194 EnerLyaProf2 = 0.; 03195 EnerLyaProf3 = 0.; 03196 EnerLyaProf4 = 0.; 03197 } 03198 03199 /* is Lya pumping enabled? */ 03200 if( FeII.lgLyaPumpOn ) 03201 { 03202 for( ipLo=0; ipLo < (FeII.nFeIILevel - 1); ipLo++ ) 03203 { 03204 for( ipHi=ipLo + 1; ipHi < FeII.nFeIILevel; ipHi++ ) 03205 { 03206 /* on first iteration optical depth in line is inward only, on later 03207 * iterations is total optical depth */ 03208 if( iteration == 1 ) 03209 { 03210 taux = Fe2LevN[ipHi][ipLo].Emis->TauIn; 03211 } 03212 else 03213 { 03214 taux = Fe2LevN[ipHi][ipLo].Emis->TauTot; 03215 } 03216 03217 /* Gup_ov_Glo is ratio of g values */ 03218 Gup_ov_Glo = Fe2LevN[ipHi][ipLo].Hi->g/Fe2LevN[ipHi][ipLo].Lo->g; 03219 03220 /* the energy of the FeII line */ 03221 EnergyWN = Fe2LevN[ipHi][ipLo].EnergyWN; 03222 03223 if( EnergyWN >= EnerLyaProf1 && EnergyWN <= EnerLyaProf4 && taux > 1 ) 03224 { 03225 /* this branch, line is within the Lya profile */ 03226 03227 /* 03228 * Lya source function, at peak is PhotOccNumLyaCenter, 03229 * 03230 * Prof2 Prof3 03231 * ---------- 03232 * / \ 03233 * / \ 03234 * / \ 03235 * ====================== 03236 * Prof1 Prof4 03237 * 03238 */ 03239 03240 if( EnergyWN < EnerLyaProf2 ) 03241 { 03242 /* linear interpolation on edge of trapazoid */ 03243 PhotOccNum_at_nu = PhotOccNumLyaCenter*(EnergyWN - EnerLyaProf1)/ de; 03244 } 03245 else if( EnergyWN < EnerLyaProf3 ) 03246 { 03247 /* this is the central plateau */ 03248 PhotOccNum_at_nu = PhotOccNumLyaCenter; 03249 } 03250 else 03251 { 03252 /* linear interpolation on edge of trapazoid */ 03253 PhotOccNum_at_nu = PhotOccNumLyaCenter*(EnerLyaProf4 - EnergyWN)/de; 03254 } 03255 03256 /* at this point Lya source function at FeII line energy is defined, but 03257 * we need to multiply by line width in Hz, 03258 * >>refer fe2 pump rate Netzer, H., Elitzur, M., & Ferland, G.J., 1985, ApJ, 299, 752-762*/ 03259 03261 /* width of FeII line in Hz */ 03262 FeIILineWidthHz = 1.e8/(EnergyWN*299792.5)*sqrt(log(taux))*DoppVel.doppler[ipIRON]; 03263 03264 /* final Lya pumping rate, s^-1*/ 03265 PumpRate = FeIILineWidthHz * PhotOccNum_at_nu * Fe2LevN[ipHi][ipLo].Emis->Aul* 03266 powi(82259.0f/EnergyWN,3); 03267 /* above must be bogus, use just occ num times A */ 03268 PumpRate = Fe2LevN[ipHi][ipLo].Emis->Aul* PhotOccNum_at_nu; 03269 03270 /* Lya pumping rate from ipHi to lower n */ 03271 Fe2LPump[ipHi][ipLo] += PumpRate; 03272 03273 /* Lya pumping rate from n to upper ipHi */ 03274 Fe2LPump[ipLo][ipHi] += PumpRate*Gup_ov_Glo; 03275 } 03276 } 03277 } 03278 } 03279 03280 return; 03281 } 03282 03283 /* end work around bugs in HP compiler */ 03284 #if defined(__HP_aCC) 03285 #pragma OPTIMIZE OFF 03286 #pragma OPTIMIZE ON 03287 #endif