/*pop371(void) - main feii routine, called by CoolIron to evaluate iron cooling */ /*FeIIData(void) - read in needed data from file convert form of feii data, as read in from file within routine FeIIData into physical form. called by CreateData */ /*FeIIPrint(void) print feii information */ /*FeIICollStrength(void) - make up collision data for feii */ /*FeIIPrint - print output from large feii atom, called by prtzone */ /*FeIISumBand, sum up large FeII emission over certain bands, called in lineset4 */ /*FeIITauInc called once per zone in tauinc to increment large FeII atom line optical depths */ /*FeIITauAver reset optical depths for large FeII atom, called by update after each iteration */ /*FeIIPoint called by CreatePoint to create pointers for lines in large FeII atom */ /*FeIIAccel called by forlin to compute radiative acceleration due to FeII lines */ /*RTFeIIMake called by RTMake, does large FeII atom radiative transfer */ /*FeIIAddLines save accumulated FeII intensities called by lineset4 */ /*FeIIPunchLines punch feii lines at end of calculation, if punch verner set, called by dopunch*/ /*FeIIEmitOut add large FeII emission to outward beam - called once per zone in RTDiffuse */ /*FeIITauInit zero out storage for large FeII atom, called in tauout */ /*FeIIZero initialize some variables, called by zero */ /*FeIIReset reset some variables, called by zero */ /*FeIIPunDepart(FILE* ioPUN ) punch some departure coef for large atom, set with punch feii departure command*/ /* send the departure coef for physical level nPUN to unit ioPUN */ /*FeIIPun1Depart( the io unit where the print should be directed FILE * ioPUN , the physical (not c) number of the level nPUN );*/ /*long int cdGetFeIIBands( void) returns number of feii bands */ #include #include #include "cddefines.h" #include "cddrive.h" #include "opac.h" #include "fe2cool.h" #include "physconst.h" #include "norm.h" #include "taulines.h" #include "rfield.h" #include "radius.h" #include "itercnt.h" #include "linelabl.h" #include "abscf.h" #include "ipoint.h" #include "ffmtread.h" #include "linpack.h" #include "tav.h" #include "fe2ovrlap.h" #include "hydrogenic.h" #include "linetauinc.h" #include "taumin.h" #include "sexp.h" #include "linesave.h" #include "rt.h" #include "wind.h" #include "fourpi.h" #include "path.h" #include "trace.h" #include "powi.h" #include "phycon.h" #include "ionfracs.h" #include "tepowers.h" #include "lyaext.h" #include "hwidth.h" #include "turb.h" #include "holod.h" #include "pop371.h" /* this will be the smallest collision strength we will permit with the gbar * collision strengths, or for the data that have zeroes */ /* >>chng 99 jul 24, this was 1e-9 in the old fortran version and in baldwin et al. 96, * and has been changed several times, and affects results. this is the smallest * non-zero collision strength in the r-matrix calculations */ #define CS2SMALL 1e-5f /* routines used only within this file */ /* evaluate collision strenths, both interpolating on r-mat and creating g-bar */ void FeIICollStrength(void); /*extern float Fe2LevN[NFE2LEVN][NFE2LEVN][NTA];*/ /*extern float Fe2LevN[ipHi][ipLo].t[NTA];*/ /*float ***Fe2LevN;*/ EmLine **Fe2LevN; /* all following variables have file scope */ #define NFEIILINES 68635 /* this is size of nnPradDat array */ #define NPRADDAT 159 /* band wavelength, lower and upper bounds, in vacuum Angstroms */ /* FeII_Bands[n][3], where n is the number of bands in fe2bands.dat * these bands are defined in fe2bands.dat and read in at startup * of calculation */ float **FeII_Bands; /* number of bands in fe2bands.dat */ long int FeIInBands; /* the dim of this vector this needs one extra since there are 20 numbers per line, * highest not ever used for anything */ /*long int nnPradDat[NPRADDAT+1];*/ static long int *nnPradDat; /* malloced in feiidata */ /* float sPradDat[NPRADDAT][NPRADDAT][8];*/ /* float sPradDat[ipHi][ipLo].t[8];*/ static float ***sPradDat; /* array used to integrate feii line intensities over model, * Fe2SavN[upper][lower], *static float Fe2SavN[NFE2LEVN][NFE2LEVN];*/ static float **Fe2SavN; /* Fe2DepCoef[NFE2LEVN]; float cli[NFEIILINES], cfe[NFEIILINES];*/ static float *cli , *cfe ; static double /* departure coefficients */ *Fe2DepCoef , /* level populations */ *Fe2LevelPop , /* this will become array of Boltzmann factors */ *FeIIBoltzmann ; /*FeIIBoltzmann[NFE2LEVN] ,*/ static long int *kli1, *kli2; static long int *kfeu, *kfed; static long int codeline, nFeIIOtherOverlap=0, nFeIISelfOverlap=0, codefe2; static double EnerLyaProf1, EnerLyaProf4, clya, hfelya, efelya, /* ask dima, why are these dimensioned 1000? */ tfelya[1000], bfelya[1000], cfelya[1000]; static long int iprov, ifelya, kfelya[1000], jfelya[1000]; /* *===================================================================== */ /* FeIIData read in feii data from files on disk. called by CreateData * but only if FeII.lgFeIION is true, set with atom feii verner command */ void FeIIData(void) { FILE *ioDATA; char chLine[132] , chFilename[134] ; long int i, ipHi , ipLo, lo, ihi, k, m1, m2, m3; # ifdef DEBUG_FUN fputs( "<+>FeIIData()\n", debug_fp ); # endif if( lgFeIIMalloc ) { /* we have already been called one time, just bail out */ # ifdef DEBUG_FUN fputs( " <->FeIIData()\n", debug_fp ); # endif return; } /* now set flag so never done again - this is set FALSE in cdDefines * when this is true it is no longer possible to change the number of levels * in the model atom with the atom feii levels command */ lgFeIIMalloc = TRUE; /* remember how many levels this was, so that in future calculations * we can reset the atom to full value */ nFeIILevelAlloc = nFeIILevel; /* set up array to save FeII line intensities */ if( (Fe2SavN = (float **)malloc(sizeof(float *)*(unsigned long)nFeIILevel )) ==NULL) { fprintf( ioQQQ, " FeIIData could not malloc Fe2SavN 1\n" ); puts( "[Stop in FeIIData]" ); exit(1); } /* second dimension, lower level, for line save array */ for( ipHi=1; ipHi < nFeIILevel; ++ipHi ) { if( (Fe2SavN[ipHi]=(float*)malloc(sizeof(float )*(unsigned long)ipHi)) ==NULL) { fprintf( ioQQQ, " FeIIData could not malloc Fe2SavN 2\n" ); puts( "[Stop in FeIIData]" ); exit(1); } } /* now malloc space for energy level table*/ if( (nnPradDat = (long*)malloc( (NPRADDAT+1)*sizeof(long) )) == NULL ) { fprintf(ioQQQ," FeIIData could not malloc nnPradDat\n"); puts( "[Stop in FeIIData]" ); exit(1); } if( ( (kli1 = (long*)malloc( nLevel1*sizeof(long) )) == NULL) || ( (kli2 = (long*)malloc( nLevel1*sizeof(long) )) == NULL ) ) { fprintf(ioQQQ," FeIIData could not malloc kli1 or kli2\n"); puts( "[Stop in FeIIData]" ); exit(1); } if( ( (kfeu = (long*)malloc( NFEIILINES*sizeof(long) )) == NULL) || ( (kfed = (long*)malloc( NFEIILINES*sizeof(long) )) == NULL ) ) { fprintf(ioQQQ," FeIIData could not malloc kfeu or kfed\n"); puts( "[Stop in FeIIData]" ); exit(1); } /*Fe2DepCoef[NFE2LEVN];*/ if( ( Fe2DepCoef = (double*)malloc( (unsigned long)nFeIILevel*sizeof(double) )) == NULL ) { fprintf(ioQQQ," FeIIData could not malloc Fe2DepCoef\n"); puts( "[Stop in FeIIData]" ); exit(1); } /*Fe2LevelPop[NFE2LEVN];*/ if( ( Fe2LevelPop = (double*)malloc( (unsigned long)nFeIILevel*sizeof(double) )) == NULL ) { fprintf(ioQQQ," FeIIData could not malloc Fe2LevelPop\n"); puts( "[Stop in FeIIData]" ); exit(1); } /*FeIIBoltzmann[NFE2LEVN];*/ if( ( FeIIBoltzmann = (double*)malloc( (unsigned long)nFeIILevel*sizeof(double) )) == NULL ) { fprintf(ioQQQ," FeIIData could not malloc FeIIBoltzmann\n"); puts( "[Stop in FeIIData]" ); exit(1); } /*float cli[NFEIILINES], cfe[NFEIILINES];*/ if( ( (cli = (float*)malloc( (unsigned long)nFeIILevel*sizeof(float) )) == NULL) || ( (cfe = (float*)malloc( (unsigned long)nFeIILevel*sizeof(float) )) == NULL ) ) { fprintf(ioQQQ," FeIIData could not malloc cli or cfe\n"); puts( "[Stop in FeIIData]" ); exit(1); } /* malloc the float sPradDat[NPRADDAT][NPRADDAT][8] array */ /* malloc the float sPradDat[ipHi][ipLo].t[8] array */ if( (sPradDat = ((float ***)malloc(NPRADDAT*sizeof(float *)))) == NULL ) { fprintf(ioQQQ," FeIIData could not malloc 1 sPradDat\n"); puts( "[Stop in FeIIData]" ); exit(1); } for (ipHi=0; ipHi < NPRADDAT; ipHi++) { if( (sPradDat[ipHi] = (float **)malloc((unsigned long)ipHi*sizeof(float ))) == NULL ) { fprintf(ioQQQ," FeIIData could not malloc 2 sPradDat\n"); puts( "[Stop in FeIIData]" ); exit(1); } /* now make space for the second dimension */ for( ipLo=0; ipLo< ipHi; ipLo++ ) { if( (sPradDat[ipHi][ipLo] = (float *)malloc(8*sizeof(float ))) == NULL ) { fprintf(ioQQQ," FeIIData could not malloc 3 sPradDat\n"); puts( "[Stop in FeIIData]" ); exit(1); } } } /* now set junk to make sure reset before used */ for (ipHi=0; ipHi < NPRADDAT; ipHi++) { for( ipLo=0; ipLo< ipHi; ipLo++ ) { for( k=0; k<8; ++k ) { sPradDat[ipHi][ipLo][k] = -FLT_MAX; } } } /* now create the main emission line arrays */ /*if( (Fe2LevN=(float***)malloc(sizeof(float **)*(unsigned long)nFeIILevel)) ==NULL)*/ if( (Fe2LevN=(EmLine**)malloc(sizeof(EmLine *)*(unsigned long)nFeIILevel)) ==NULL) { fprintf( ioQQQ, " FeIIData could not malloc Fe2LevN 1\n" ); puts( "[Stop in FeIIData]" ); exit(1); } for( ipHi=1; ipHi < nFeIILevel; ++ipHi ) { /*if( (Fe2LevN[ipHi]=(float**)malloc(sizeof(float *)*(unsigned long)ipHi))==NULL)*/ if( (Fe2LevN[ipHi]=(EmLine*)malloc(sizeof(EmLine)*(unsigned long)ipHi))==NULL) { fprintf( ioQQQ, " FeIIData could not malloc Fe2LevN 2\n" ); puts( "[Stop in FeIIData]" ); exit(1); } } /* check on path if path set */ /* path was parsed in getset */ if( lgDataPathSet == TRUE ) { /*path set, so look only there */ strcpy( chFilename , chDataPath ); strcat( chFilename , "fe2nn.dat" ); } else { /* path not set, check local space only */ strcpy( chFilename , "fe2nn.dat" ); } if( trace.lgTrace ) { fprintf( ioQQQ," FeIIData opening fe2nn.dat:"); } if( ( ioDATA = fopen( chFilename , "r" ) ) == NULL ) { fprintf( ioQQQ, " FeIIData could not open fe2nn.dat\n" ); if( lgDataPathSet == TRUE ) { fprintf( ioQQQ, " FeIIData could not open fe2nn.dat, even tried path.\n"); fprintf( ioQQQ, " path is *%s*\n",chDataPath ); fprintf( ioQQQ, " final path is *%s*\n",chFilename ); } puts( "[Stop in FeIIData]" ); exit(-1); } for( i=0; i < 8; i++ ) { if( fgets( chLine , sizeof(chLine) , ioDATA ) == NULL ) { fprintf( ioQQQ, " fe2nn.dat error reading data\n" ); puts( "[Stop in FeIIData]" ); exit(-1); } /* get these integers from fe2nn.dat */ k = 20*i; /* NPRADDAT is size of nnPradDat array, 159+1, make sure we do not exceed it */ assert( k+19 < NPRADDAT+1 ); sscanf( chLine , "%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld", &nnPradDat[k+0], &nnPradDat[k+1], &nnPradDat[k+2], &nnPradDat[k+3], &nnPradDat[k+4], &nnPradDat[k+5], &nnPradDat[k+6], &nnPradDat[k+7], &nnPradDat[k+8], &nnPradDat[k+9], &nnPradDat[k+10],&nnPradDat[k+11], &nnPradDat[k+12],&nnPradDat[k+13],&nnPradDat[k+14], &nnPradDat[k+15],&nnPradDat[k+16], &nnPradDat[k+17],&nnPradDat[k+18],&nnPradDat[k+19] ); # if !defined(NDEBUG) for( m1=0; m1<20; ++m1 ) { assert( nnPradDat[k+m1] >= 0 && nnPradDat[k+m1] <= NFE2LEVN ); } # endif } fclose(ioDATA); /* now get radiation data */ if( lgDataPathSet == TRUE ) { /*path set, so look only there */ strcpy( chFilename , chDataPath ); strcat( chFilename , "fe2rad.dat" ); } else { /* path not set, check local space only */ strcpy( chFilename , "fe2rad.dat" ); } if( trace.lgTrace ) { fprintf( ioQQQ," FeIIData opening fe2rad.dat:"); } if( ( ioDATA = fopen( chFilename , "r" ) ) == NULL ) { fprintf( ioQQQ, " FeIIData could not open fe2rad.dat\n" ); if( lgDataPathSet == TRUE ) { fprintf( ioQQQ, " FeIIData could not open fe2rad.dat, even tried path.\n"); fprintf( ioQQQ, " path is *%s*\n",chDataPath ); fprintf( ioQQQ, " final path is *%s*\n",chFilename ); } puts( "[Stop in FeIIData]" ); exit(-1); } /* file now open, read the data */ for( ipHi=1; ipHi < nFeIILevel; ipHi++ ) { for( ipLo=0; ipLo < ipHi; ipLo++ ) { float save[2]; if( fgets( chLine , sizeof(chLine) , ioDATA ) == NULL ) { fprintf( ioQQQ, " fe2nn.dat error reading data\n" ); puts( "[Stop in FeIIData]" ); exit(-1); } /* first and last number on this line */ sscanf( chLine , "%ld%ld%ld%ld%f%f%ld", &lo, &ihi, &m1, &m2 , &save[0], &save[1] , &m3); /* the pointeres ihi and lo within this array were * read in from the line above */ Fe2LevN[ihi-1][lo-1].gLo = (float)m1; Fe2LevN[ihi-1][lo-1].gHi = (float)m2; Fe2LevN[ihi-1][lo-1].Aul = save[0]; Fe2LevN[ihi-1][lo-1].EnergyWN = save[1]; Fe2LevN[ihi-1][lo-1].cs1 = (float)m3; } } fclose(ioDATA); /* now read collision data in fe2col.dat */ if( lgDataPathSet == TRUE ) { /*path set, so look only there */ strcpy( chFilename , chDataPath ); strcat( chFilename , "fe2col.dat" ); } else { /* path not set, check local space only */ strcpy( chFilename , "fe2col.dat" ); } if( trace.lgTrace ) { fprintf( ioQQQ," FeIIData opening fe2col.dat:"); } if( ( ioDATA = fopen( chFilename , "r" ) ) == NULL ) { fprintf( ioQQQ, " FeIIData could not open fe2col.dat\n" ); if( lgDataPathSet == TRUE ) { fprintf( ioQQQ, " FeIIData could not open fe2col.dat, even tried path.\n"); fprintf( ioQQQ, " path is *%s*\n",chDataPath ); fprintf( ioQQQ, " final path is *%s*\n",chFilename ); } puts( "[Stop in FeIIData]" ); exit(-1); } for( ipHi=1; ipHi 0. ); assert( Fe2LevN[ipHi][ipLo].Aul> 0. ); /* now put into standard internal line format */ Fe2LevN[ipHi][ipLo].WLAng = (float)(1.e8/Fe2LevN[ipHi][ipLo].EnergyWN); /* generate gf from A */ Fe2LevN[ipHi][ipLo].gf = (float)(Fe2LevN[ipHi][ipLo].Aul*Fe2LevN[ipHi][ipLo].gHi* 1.4992e-8*POW2(1e4/Fe2LevN[ipHi][ipLo].EnergyWN)); Fe2LevN[ipHi][ipLo].IonStg = 2; Fe2LevN[ipHi][ipLo].nelem = 26; Fe2LevN[ipHi][ipLo].iRedisFun = 1 ; } } /* finally get FeII bands, this sets */ i = cdGetFeIIBands(); # ifdef DEBUG_FUN fputs( " <->FeIIData()\n", debug_fp ); # endif return; } /* *===================================================================== */ /*********************************************************************** *** version of pop371.f with overlap in procces 05.28.97 ooooooo ******change in common block *te* sqrte 05.28.97 *******texc is fixed 03.28.97 *********this version has work on overlap *********this version has # of zones (ZoneCnt) 07.03.97 *********taux - optical depth depends on iter correction 03.03.97 *********calculate cooling (fcool) and output fecool from Cloudy 01.29.97god *********and cooling derivative (dfcool) ************ this program for 371 level iron model 12/14/1995 ************ this program for 371 level iron model 1/11/1996 ************ this program for 371 level iron model 3/21/1996 ************ this program without La 3/27/1996 ************ this program for 371 level iron model 4/9/1996 ************ includes:FeIICollStrength;cooling;overlapping in lines */ void pop371( void ) { long int i, ipHi , ipLo , info, kj, n; /* used in test for non-positive level populations */ int lgPopNeg ; double EnerLyaProf2, EnerLyaProf3, aeff, blya, clfe, de, defe, dlya, ee, gg, pp, sfe, sla, sum , sumfe, taux, tex, tex1; long int ipiv[NFE2LEVN]; double /* the yVector - will become level populations after matrix inversion */ yVector[NFE2LEVN], /*b[NFE2LEVN][NFE2LEVN], */ **b, /* this is used to call matrix routines */ /*xMatrix[NFE2LEVN][NFE2LEVN] ,*/ **xMatrix , /* gary change - this had been a 2-D array, but only 1 was ever used */ texc[NFE2LEVN], /* this will become the very large 1-D array that * is passed to the matrix inversion routine*/ *amat; # ifdef DEBUG_FUN fputs( "<+>pop371()\n", debug_fp ); # endif if( trace.lgTrace ) { fprintf( ioQQQ," pop371 fe2 pops called"); } /* ========================================================================== * malloc space for the 2D arrays */ if( (b = (double **)malloc(sizeof(double *)*(unsigned long)nFeIILevel ) ) ==NULL ) { fprintf( ioQQQ, " could not malloc ba array in 1D\n" ); puts( "[Stop in pop371]" ); exit(1); } if( (xMatrix = (double **)malloc(sizeof(double *)*(unsigned long)nFeIILevel ))==NULL ) { fprintf( ioQQQ, " could not malloc xMatrix arrays in 1D\n" ); puts( "[Stop in pop371]" ); exit(1); } /* now do the second dimension */ for( i=0; i SMALLFLOAT ) { clya = (1.0 - HydroLines[0][IP2P][0].Pesc)/ (exp(118356./LyaExT.TexcLya) - 1.0); } else { /* lya excitation temperature not available */ clya = 0.; } } else { clya = 0.f; de = 0.; EnerLyaProf1 = 0.; EnerLyaProf2 = 0.; EnerLyaProf3 = 0.; EnerLyaProf4 = 0.; } /* make the g-bar collision strengths and do linear interpolation on r-matrix data. * this also sets Boltzmann factors for all levels */ FeIICollStrength(); /*******To include overlapping with Fe-lines: parameter nFeIISelfOverlap>0 *******with other-lines: parameter nFeIIOtherOverlap >0;****************** *******they come first from Fe2OvrLap at the end*********** *******of first iterration********************************** */ sumfe = 0.0f; ifelya = 0; for( ipLo=0; ipLo < (nFeIILevel - 1); ipLo++ ) { for( ipHi=ipLo + 1; ipHi < nFeIILevel; ipHi++ ) { /* on first iteration optical depth in line is inward only, on later * iterations is total optical depth */ if( IterCnt.iter == 1 ) { taux = Fe2LevN[ipHi][ipLo].TauIn; } else { taux = Fe2LevN[ipHi][ipLo].TauTot; } /* gg is ratio of g values */ gg = Fe2LevN[ipHi][ipLo].gHi/Fe2LevN[ipHi][ipLo].gLo; /* collisional deexcitation */ b[ipHi][ipLo] = 8.629e-06/TePowers.sqrte*Fe2LevN[ipHi][ipLo].cs/ Fe2LevN[ipHi][ipLo].gHi*phycon.eden; /* this is collisional excitation to balance above collisional deexcitation */ b[ipLo][ipHi] = b[ipHi][ipLo]*gg*FeIIBoltzmann[ipHi]/FeIIBoltzmann[ipLo]; /* continuum pumping upward followed by downward, * pumping rate was established in */ b[ipLo][ipHi] += Fe2LevN[ipHi][ipLo].pump; b[ipHi][ipLo] += Fe2LevN[ipHi][ipLo].pump/ gg; /* if no La-pumping do go to 30 */ if( FeII.lgLyaPumpOn ) { ee = Fe2LevN[ipHi][ipLo].EnergyWN; /********line shows that this model includes lya after 2nd iteration * * proverka peresechneiya s Lya * esli peresekaetsya, iprov=1 */ iprov = -1; if( (ee >= EnerLyaProf1 && ee <= EnerLyaProf4) && IterCnt.iter > 1 ) { /* end Lya-pumping */ if( ee < EnerLyaProf2 ) { dlya = clya*(ee - EnerLyaProf1)/ de; iprov = 1; } else if( ee < EnerLyaProf3 ) { dlya = clya; iprov = 1; } else { iprov = 1; dlya = clya*(EnerLyaProf4 - ee)/de; } /* zapolnenie massiva opticheski tolstyh Fe2 linij, peresekayuschixsya * s Lya */ if( Fe2LevN[ipHi][ipLo].TauTot >= 0.5f && iprov > 0 ) { /* ifelya - nomer takoj linii po poryadku */ ifelya += 1; /* kfelya(ifelya) - verxnij uroven' takoj linii */ kfelya[ifelya-1] = ipHi+1; /* jfelya(ifelya) - nizhnij uroven' takoj linii */ jfelya[ifelya-1] = ipLo+1; /* tfelya(ifelya) - opticheskaya tolscha takoj linii */ tfelya[ifelya-1] = Fe2LevN[ipHi][ipLo].TauTot; /* hfelya - opticheskaya tolscha Lya */ hfelya = HydroLines[0][IP2P][0].TauIn; /* gjf change * hfelya=htau(2,1) * efelya - elektronnaya temperatura */ efelya = phycon.te; bfelya[ifelya-1] = b[ipHi][ipLo]; cfelya[ifelya-1] = (double)(dlya* Fe2LevN[ipHi][ipLo].Aul* powi(82259.0/ee,3)); } /*********to calculate Lya-distraction x :sfe,defe,clfe,sla,sumfe *********sfe-integral for 1 iron line *********defe-halfwidth iron line *********clfe-intensity of iron line in ee-point *********sla-integral of Lya *********sumfe- integral for all iron lines inside La square * */ if( (taux > 1 && Fe2LevN[ipHi][ipLo].PopHi > 0) && Fe2LevN[ipHi][ipLo].PopLo > 0 ) { pp = Fe2LevN[ipHi][ipLo].PopHi/ Fe2LevN[ipHi][ipLo].PopLo; clfe = (1.0f - Fe2LevN[ipHi][ipLo].Pesc)/ (gg/pp - 1.0f); defe = (double)(1.e8/(ee*299792.5)*pow(log(taux),0.5)* pow(2.97756*1e-04*phycon.te + powi(turb.TurbVel/ 1.e5,2),0.5) ); sfe = 2*defe*clfe; if( clfe >= dlya ) sfe = 2*defe*dlya; sumfe += sfe; sla = clya*(EnerLyaProf4 - EnerLyaProf1 + EnerLyaProf3 - EnerLyaProf2)/2; } blya = dlya*Fe2LevN[ipHi][ipLo].Aul* powi(82259.0f/ee,3); b[ipHi][ipLo] += blya; b[ipLo][ipHi] += blya*gg; } } /* beginning of Line-pumping */ /*fixit(); this code is dangerous since kli1 is not safely dimned*/ if( codeline != 2 ) { if( nFeIIOtherOverlap >= 1 ) { for( kj=0; kj < nFeIIOtherOverlap; kj++ ) { if( kli1[kj] == ipHi+1 && kli2[kj] == ipLo+1 ) { b[ipHi][ipLo] += cli[kj]; b[ipLo][ipHi] += cli[kj]*gg; } } } } /*******************Fe-pumping if no go to 222 * */ if( codefe2 != 2 ) { if( nFeIISelfOverlap >= 1 ) { for( kj=0; kj < nFeIISelfOverlap; kj++ ) { if( kfeu[kj] == ipHi+1 && kfed[kj] == ipLo+1 ) { b[ipHi][ipLo] += cfe[kj]; b[ipLo][ipHi] += cfe[kj]*gg; } } } } } } if( sla > 1e-20f && IterCnt.iter > 1 ) { hydro.dstfe2lya = (float)(sumfe/sla); } for( ipLo=0; ipLo < (nFeIILevel - 1); ipLo++ ) { xMatrix[ipLo][nFeIILevel-1] = -b[ipLo][nFeIILevel-1]; } xMatrix[nFeIILevel-1][nFeIILevel-1] = 0.0f; for( ipLo=0; ipLo < (nFeIILevel - 1); ipLo++ ) { aeff = Fe2LevN[nFeIILevel-1][ipLo].Aul*(Fe2LevN[nFeIILevel-1][ipLo].Pesc + Fe2LevN[nFeIILevel-1][ipLo].Pdest); xMatrix[nFeIILevel-1][nFeIILevel-1] += aeff + b[nFeIILevel-1][ipLo]; } for( n=1; n < (nFeIILevel - 1); n++ ) { for( ipLo=0; ipLo < n; ipLo++ ) { xMatrix[ipLo][n] = -b[ipLo][n]; } for( ipHi=n + 1; ipHi < nFeIILevel; ipHi++ ) { aeff = Fe2LevN[ipHi][n].Aul*(Fe2LevN[ipHi][n].Pesc + Fe2LevN[ipHi][n].Pdest); xMatrix[ipHi][n] = -aeff - b[ipHi][n]; } xMatrix[n][n] = 0.0f; for( ipLo=0; ipLo < n; ipLo++ ) { aeff = Fe2LevN[n][ipLo].Aul*(Fe2LevN[n][ipLo].Pesc + Fe2LevN[n][ipLo].Pdest); xMatrix[n][n] += aeff + b[n][ipLo]; } for( i=n + 1; i < nFeIILevel; i++ ) { xMatrix[n][n] += b[n][i]; } } /* the top row of the matrix is the sum of level populations */ for( i=0; i < nFeIILevel; i++ ) { xMatrix[i][0] = 1.0; } /* define the Y Vector. The oth element is the sum of all level populations * adding up to the total population. The remaining elements are the level * balance equations adding up to zero */ yVector[0] = 1.0; for( n=1; n < nFeIILevel; n++ ) { yVector[n] = 0.0; } /* create the 1-yVector array that will save vector, * this is the macro trick */ # ifdef AMAT # undef AMAT # endif # define AMAT(I_,J_) (*(amat+(I_)*nFeIILevel+(J_))) /* malloc space for the 1-yVector array */ if( (amat=(double*)malloc( (sizeof(double)*(unsigned long)(nFeIILevel*nFeIILevel) ))) == NULL ) { fprintf( ioQQQ, " pop371 malloc amat error\n" ); puts( "[Stop in pop371]" ); exit(1); } /* copy current contents of xMatrix array over to special amat array*/ for( ipHi=0; ipHi < nFeIILevel; ipHi++ ) { for( i=0; i < nFeIILevel; i++ ) { /*amat[i][ipHi] = z[i][ipHi];*/ AMAT(i,ipHi) = xMatrix[i][ipHi]; } } /*void DGETRF(long,long,double*,long,long[],long*);*/ /*DGETRF(NFE2LEVN,NFE2LEVN,(double*)xMatrix,NFE2LEVN, ipiv,&info);*/ DGETRF(nFeIILevel,nFeIILevel,amat,nFeIILevel, ipiv,&info); /*void DGETRS(long int TRANS, long int N, long int NRHS, double *A, long int LDA, long int IPIV[], double *B, long int LDB, long int *INFO);*/ /*DGETRS('N',1,NFE2LEVN,(double*)xMatrix,NFE2LEVN, ipiv,yVector,NFE2LEVN,&info);*/ /*DGETRS('N',1,NFE2LEVN,amat,NFE2LEVN, ipiv,yVector,NFE2LEVN,&info);*/ /* gary change, order of 2 and 3 parameter switched?? */ DGETRS('N',nFeIILevel,1,amat,nFeIILevel, ipiv,yVector,nFeIILevel,&info); /* free up the space we allocated above */ free( amat ); /* yVector now contains the level populations */ /* this better be false after this loop - if not then non-positive level pops */ lgPopNeg = FALSE; /* copy all level pops over to Fe2LevelPop */ for( ipLo=0; ipLo < nFeIILevel; ipLo++ ) { if(yVector[ipLo] <= 0. ) { lgPopNeg = TRUE; fprintf(ioQQQ," pop371 finds non-positive level population, leve is %ld pop is %g\n", ipLo , yVector[ipLo] ); } Fe2LevelPop[ipLo] = yVector[ipLo] * xIonFracs[2][25]; } if( lgPopNeg ) { fprintf(ioQQQ," TROUBLE pop371 finds non-positive level population\n" ); } /* now set line opacities */ for( ipLo=0; ipLo < (nFeIILevel - 1); ipLo++ ) { for( ipHi=ipLo+1; ipHi < nFeIILevel; ipHi++ ) { Fe2LevN[ipHi][ipLo].PopOpc = (float)((yVector[ipLo] - yVector[ipHi]*Fe2LevN[ipHi][ipLo].gLo/Fe2LevN[ipHi][ipLo].gHi)* xIonFracs[2][25]); Fe2LevN[ipHi][ipLo].PopLo = (float)yVector[ipLo]*xIonFracs[2][25]; Fe2LevN[ipHi][ipLo].PopHi = (float)yVector[ipHi]*xIonFracs[2][25]; Fe2LevN[ipHi][ipLo].phots = (float)(yVector[ipHi]* Fe2LevN[ipHi][ipLo].Aul*Fe2LevN[ipHi][ipLo].Pesc* xIonFracs[2][25]); Fe2LevN[ipHi][ipLo].xIntensity = Fe2LevN[ipHi][ipLo].phots* Fe2LevN[ipHi][ipLo].EnergyErg; } } /*********calculate cooling (fcool) and output fecool from Cloudy *********and cooling derivative (dfcool) * * this is net cooling, cooling minus heating */ holod.fcool = 0.0f; /* this is be heating, not heating minus cooling */ holod.feheat = 0.f; holod.dfcool = 0.0f; /* next two blocks determine departure coefficients for the atom */ /* first sum up partition function for the model atom */ Fe2DepCoef[0] = 1.0; sum = 1.0; for( i=1; i < nFeIILevel; i++ ) { /* boltzmann factor relative to ground times ratio of statistical weights */ Fe2DepCoef[i] = Fe2DepCoef[0]*FeIIBoltzmann[i]* Fe2LevN[i][0].gHi/ Fe2LevN[1][0].gLo; /* this sum is the partition function for the atom */ sum += Fe2DepCoef[i]; } /* now fill in departure coefficients */ for( i=0; i < nFeIILevel; i++ ) { /* divide by total partition function, Fe2DepCoef is now the fraction of the * population that is in this level in TE */ Fe2DepCoef[i] /= sum; /* convert to true departure coefficient */ if( Fe2DepCoef[i]>SMALLFLOAT ) { Fe2DepCoef[i] = yVector[i]/Fe2DepCoef[i]; } else { Fe2DepCoef[i] = 0.; } } /* compute heating and cooling due to model atom */ for( ipLo=0; ipLo < (nFeIILevel - 1); ipLo++ ) { for( ipHi=ipLo + 1; ipHi < nFeIILevel; ipHi++ ) { b[ipHi][ipLo] = 8.629e-06/TePowers.sqrte*Fe2LevN[ipHi][ipLo].cs/ Fe2LevN[ipHi][ipLo].gHi*phycon.eden; /* ratio of statistical weights */ gg = Fe2LevN[ipHi][ipLo].gHi/Fe2LevN[ipHi][ipLo].gLo; b[ipLo][ipHi] = b[ipHi][ipLo]*gg*FeIIBoltzmann[ipHi]/FeIIBoltzmann[ipLo]; holod.fcool += (float)(1.98648e-16f*(b[ipLo][ipHi]*yVector[ipLo] - b[ipHi][ipLo]*yVector[ipHi])*xIonFracs[2][25]* Fe2LevN[ipHi][ipLo].EnergyWN); holod.feheat += (float)(1.98648e-16f*(b[ipHi][ipLo]*yVector[ipHi])* xIonFracs[2][25]*Fe2LevN[ipHi][ipLo].EnergyWN); } } /* find excitation temperature of each level */ for( ipHi=1; ipHi < nFeIILevel; ipHi++ ) { tex = (yVector[ipHi]/Fe2LevN[ipHi][0].gHi)/ (yVector[0]/Fe2LevN[ipHi][0].gLo); tex1 = fabs(tex-1.0); if( tex1 <= 1.e-06 ) tex = phycon.te; /*texc[ipHi] = -1.438826f*Fe2LevN[ipHi][0].EnergyWN/ */ texc[ipHi] = -T1CM*Fe2LevN[ipHi][0].EnergyWN/ log(tex); } for( ipLo=0; ipLo < (nFeIILevel - 1); ipLo++ ) { for( ipHi=ipLo + 1; ipHi < nFeIILevel; ipHi++ ) { b[ipHi][ipLo] = 8.629e-06/TePowers.sqrte*Fe2LevN[ipHi][ipLo].cs/ Fe2LevN[ipHi][ipLo].gHi*phycon.eden; gg = Fe2LevN[ipHi][ipLo].gHi/Fe2LevN[ipHi][ipLo].gLo; b[ipLo][ipHi] = b[ipHi][ipLo]*gg*FeIIBoltzmann[ipHi]/FeIIBoltzmann[ipLo]; holod.dfcool += (float)(1.98648e-16f*(b[ipLo][ipHi]*yVector[ipLo] - b[ipHi][ipLo]*yVector[ipHi])*xIonFracs[2][25]* Fe2LevN[ipHi][ipLo].EnergyWN*(texc[ipHi]/ POW2(phycon.te) - 0.5f/phycon.te)); } } /* now free up the 2D arrays */ for( i=0; ipop371()\n", debug_fp ); # endif return; } /* *===================================================================== */ /*FeIICollStrength calculates collisional strenths for all lines that have no data */ void FeIICollStrength(void) { long int i, ipLo , ipHi; float ag, cg, dg, gb, y; float FracLowTe , FracHighTe; static float tt[8]={1e3f,3e3f,5e3f,7e3f,1e4f,12e3f,15e3f,2e4f}; long int ipTemp, ipTempp1 , ipLoFe2, ipHiFe2; # ifdef DEBUG_FUN fputs( "<+>FeIICollStrength()\n", debug_fp ); # endif for( ipLo=0; ipLo < (nFeIILevel - 1); ipLo++ ) { for( ipHi=ipLo + 1; ipHi < nFeIILevel; ipHi++ ) { if( Fe2LevN[ipHi][ipLo].cs1 == 3.0f ) { /************************forbidden tr**************************/ ag = 0.15f; cg = 0.f; dg = 0.f; } /************************allowed tr*******************************/ else if( Fe2LevN[ipHi][ipLo].cs1 == 1.0f ) { ag = 0.6f; cg = 0.f; dg = 0.28f; } /************************semiforbidden****************************/ else if( Fe2LevN[ipHi][ipLo].cs1 == 2.0f ) { ag = 0.0f; cg = 0.1f; dg = 0.0f; } else { /* this branch is impossible, since cs must be 1, 2, or 3 */ ag = 0.0f; cg = 0.1f; dg = 0.0f; fprintf(ioQQQ,">>>impossible iolncs1 in pop371\n"); } /*y = 1.438826f*Fe2LevN[ipHi][ipLo].EnergyWN/ phycon.te;*/ y = (float)(T1CM)*Fe2LevN[ipHi][ipLo].EnergyWN/ phycon.te; gb = (float)(ag + (-cg*POW2(y) + dg)*(log(1.0+1.0/y) - 0.4/ POW2(y + 1.0)) + cg*y); Fe2LevN[ipHi][ipLo].cs = 23.861f*1e5f*gb* Fe2LevN[ipHi][ipLo].Aul*Fe2LevN[ipHi][ipLo].gHi/ POW3(Fe2LevN[ipHi][ipLo].EnergyWN); /* confirm that collision strength is positive */ assert( Fe2LevN[ipHi][ipLo].cs > 0.); /* g-bar cs becomes unphysically small for forbidden transitions - * this sets a lower limit to the final cs - CS2SMALL is defined above */ Fe2LevN[ipHi][ipLo].cs = MAX2( CS2SMALL, Fe2LevN[ipHi][ipLo].cs); /* this was the logic used in the old fortran version, * and reproduces the results in Baldwin et al '96 if( Fe2LevN[ipHi][ipLo].cs < 1e-10 ) { Fe2LevN[ipHi][ipLo].cs = 0.01f; } */ } } /* we will interpolate on the set of listed collision strengths - * where in this set are we? */ if( phycon.te <= tt[0] ) { /* temperature is lower than lowest tabulated, use the * lowest tabulated point */ /* ipTemp usually points to the cell cooler than te, ipTemp+1 to one higher, * here both are lowest */ ipTemp = 0; ipTempp1 = 0; FracHighTe = 0.; } else if( phycon.te > tt[7] ) { /* temperature is higher than highest tabulated, use the * highest tabulated point */ /* ipTemp usually points to the cell cooler than te, ipTemp+1 to one higher, * here both are highest */ ipTemp = 7; ipTempp1 = 7; FracHighTe = 0.; } else { /* where in the array is the temperature we need? */ ipTemp = -1; for( i=0; i < 8-1; i++ ) { if( phycon.te <= tt[i+1] ) { ipTemp = i; break; } } /* this cannot possibly happen */ if( ipTemp < 0 ) { fprintf( ioQQQ, " Insanity while looking for temperature in coll str array, te=%g.\n", phycon.te ); puts( "[Stop in pop371]" ); exit(1); } /* ipTemp points to the cell cooler than te, ipTemp+1 to one higher, * do linear interpolation between */ ipTemp = i; ipTempp1 = i+1; FracHighTe = (phycon.te - tt[ipTemp])/(tt[ipTempp1] - tt[ipTemp]); } /* this is fraction of final linear interpolated collision strength that * is weighted by the lower bound cs */ FracLowTe = 1.f - FracHighTe; for( ipHi=1; ipHi < NPRADDAT; ipHi++ ) { for( ipLo=0; ipLo < ipHi; ipLo++ ) { /* ipHiFe2 should point to upper level of this pair, and * ipLoFe2 should point to lower level */ ipHiFe2 = MAX2( nnPradDat[ipHi] , nnPradDat[ipLo] ); ipLoFe2 = MIN2( nnPradDat[ipHi] , nnPradDat[ipLo] ); assert( ipHiFe2-1 < NFE2LEVN ); assert( ipHiFe2-1 >= 0 ); assert( ipLoFe2-1 < NFE2LEVN ); assert( ipLoFe2-1 >= 0 ); /* do linear interpolation for CS, these are dimensioned NPRADDAT = 159 */ if( ipHiFe2-1 < nFeIILevel ) { /* do linear interpolation */ Fe2LevN[ipHiFe2-1][ipLoFe2-1].cs = sPradDat[ipHi][ipLo][ipTemp] * FracLowTe + sPradDat[ipHi][ipLo][ipTempp1] * FracHighTe ; /* confirm that this is positive */ assert( Fe2LevN[ipHiFe2-1][ipLoFe2-1].cs > 0. ); } } } /* create boltzmann factors for all levels */ FeIIBoltzmann[0] = 1.0; for( ipHi=1; ipHi < nFeIILevel; ipHi++ ) { /*FeIIBoltzmann[ipHi] = (float)sexp( 1.438826f*Fe2LevN[ipHi][0].EnergyWN/phycon.te );*/ /* >>chng 99 may 21, from above to following slight different number from phyconst.h */ FeIIBoltzmann[ipHi] = sexp( T1CM*Fe2LevN[ipHi][0].EnergyWN/phycon.te ); } /* now possibly trim down atom if Boltzmann factors for upper levels are zero */ ipHi = nFeIILevel - 1; while( FeIIBoltzmann[ipHi] == 0. && ipHi > 0 ) { --ipHi; } /* ipHi now is the highest level with finite Boltzmann factor - * use this as the number of levels */ if( ipHi <= 1 ) { /* this is basically a sanity check and can't possibly happen */ fprintf( ioQQQ, " pop371 trimmed boltzmann factors too low.\n" ); puts( "[Stop in pop371]" ); exit(1); } else { /* is nearly all cases this does nothing since ipHi and nFeIILevel * are equal . . . */ nFeIILevel = ipHi+1; /*fprintf(ioQQQ," levels reset to %li\n",nFeIILevel);*/ } /* debug code to print out the collision strengths for some levels */ { enum {DEBUG=FALSE}; if( DEBUG) { for( ipLo=0; ipLo < 52; ipLo++ ) { fprintf(ioQQQ,"%e %e\n", Fe2LevN[51][ipLo].cs,Fe2LevN[52][ipLo].cs); } exit(1); } } # ifdef DEBUG_FUN fputs( " <->FeIICollStrength()\n", debug_fp ); # endif return; } /* *===================================================================== */ /*subroutine FeIIPrint dlya raspechatki naselennostej v cloudy.out ili v * otdel'nom file unit=33 *!!nado takzhe vklyuchit raspechatku iz perekrytiya linii */ /*FeIIPrint - print output from large feii atom, called by prtzone */ void FeIIPrint(void) { long int i, iPrtFeII; # ifdef DEBUG_FUN fputs( "<+>FeIIPrint()\n", debug_fp ); # endif if( FeII.lgPrint ) { iPrtFeII = 1; if( (FeII.lgLyaPumpOn && codeline == 1) && codefe2 == 1 ) { fprintf( ioQQQ, "lya-,line-, self-overlap are incuded\n" ); } if( (FeII.lgLyaPumpOn && codeline == 2) && codefe2 == 2 ) { fprintf( ioQQQ, "lya- included, line-, self-overlap are not incuded\n" ); } if( (!FeII.lgLyaPumpOn && codeline == 2) && codefe2 == 2 ) { fprintf( ioQQQ, "lya-,line-, self-overlap are not incuded\n" ); } if( iPrtFeII == 2 ) { fprintf( ioQQQ, "subroutine FeIIPrint is not active now\n" ); } if( ifelya > 0 ) { fprintf( ioQQQ, "Lya opt.width =%10.3e Te =%10.3e\n", hfelya, efelya ); fprintf( ioQQQ, "EnerLyaProf1 =%10.3e EnerLyaProf4 =%10.3e clya =%10.3e\n", EnerLyaProf1, EnerLyaProf4, clya ); for( i=0; i < ifelya; i++ ) { fprintf( ioQQQ, "Upper = %3ld Lower = %3ld Opt. width = %10.3e\n", kfelya[i], jfelya[i], tfelya[i] ); fprintf( ioQQQ, "A =%10.3e Blya =%10.3e\n", bfelya[i], cfelya[i] ); } } fprintf( ioQQQ, "Lya width =%10.3e\n", hwidth.HLineWidth[0] ); } # ifdef DEBUG_FUN fputs( " <->FeIIPrint()\n", debug_fp ); # endif return; } /* *===================================================================== */ /*FeIISumBand, sum up large FeII emission over certain bands, called in lineset4 */ double FeIISumBand(float wl1, float wl2) { long int ipHi, ipLo; double SumBandFe2_v, ahi, alo; # ifdef DEBUG_FUN fputs( "<+>FeIISumBand()\n", debug_fp ); # endif /*sum up large FeII emission over certain bands */ /* convert to line energy in ergs */ ahi = 1.99e-8/wl1; alo = 1.99e-8/wl2; SumBandFe2_v = 0.; for( ipLo=0; ipLo < (nFeIILevel - 1); ipLo++ ) { for( ipHi=ipLo + 1; ipHi < nFeIILevel; ipHi++ ) { if( Fe2LevN[ipHi][ipLo].EnergyErg >= alo && Fe2LevN[ipHi][ipLo].EnergyErg <= ahi ) SumBandFe2_v += Fe2LevN[ipHi][ipLo].xIntensity; } } # ifdef DEBUG_FUN fputs( " <->FeIISumBand()\n", debug_fp ); # endif return( SumBandFe2_v ); } /* *===================================================================== */ /*FeIITauInc called once per zone in tauinc to increment large FeII atom line optical depths */ void FeIITauInc(void) { long int ipHi, ipLo; # ifdef DEBUG_FUN fputs( "<+>FeIITauInc()\n", debug_fp ); # endif /*fixit(); prevent maser runaway by heavy-handed MAX */ for( ipLo=0; ipLo < (nFeIILevel - 1); ipLo++ ) { for( ipHi=ipLo + 1; ipHi < nFeIILevel; ipHi++ ) { /* needs to be fixed - this should not be necessary */ Fe2LevN[ipHi][ipLo].PopOpc = MAX2( tauminCom.taumin , Fe2LevN[ipHi][ipLo].PopOpc ); LineTauInc( &Fe2LevN[ipHi][ipLo] ); } } # ifdef DEBUG_FUN fputs( " <->FeIITauInc()\n", debug_fp ); # endif return; } /* *===================================================================== */ /*FeIITauAver reset optical depths for large FeII atom, called by update after each iteration */ void FeIITauAver(void) { long int ipHi, ipLo; # ifdef DEBUG_FUN fputs( "<+>FeIITauAver()\n", debug_fp ); # endif /* called at end of iteration */ for( ipLo=0; ipLo < (nFeIILevel - 1); ipLo++ ) { for( ipHi=ipLo + 1; ipHi < nFeIILevel; ipHi++ ) { tav( &Fe2LevN[ipHi][ipLo] ,0.5); } } /* call Katya's overlap routine */ Fe2OvrLap(); # ifdef DEBUG_FUN fputs( " <->FeIITauAver()\n", debug_fp ); # endif return; } /* *===================================================================== */ /*FeIIPoint called by CreatePoint to create pointers for lines in large FeII atom */ void FeIIPoint(void) { long int ipHi, ip, ipLo; # ifdef DEBUG_FUN fputs( "<+>FeIIPoint()\n", debug_fp ); # endif /* routine called when cloudy sets pointers for Fe2 lines, once per coreload */ for( ipLo=0; ipLo < (nFeIILevel - 1); ipLo++ ) { for( ipHi=ipLo + 1; ipHi < nFeIILevel; ipHi++ ) { Fe2LevN[ipHi][ipLo].EnergyRyd = (float)(Fe2LevN[ipHi][ipLo].EnergyWN*WAVNRYD); ip = ipoint(Fe2LevN[ipHi][ipLo].EnergyRyd); Fe2LevN[ipHi][ipLo].ipCont = ip; /* do not over write other pointers with feii since feii is everywhere */ if( strcmp(LineLabl.chLineLabel[ip-1]," ") == 0 ) strcpy( LineLabl.chLineLabel[ip-1], "FeII" ); Fe2LevN[ipHi][ipLo].damprel = (float)(Fe2LevN[ipHi][ipLo].Aul* Fe2LevN[ipHi][ipLo].WLAng*1e-8/PI4); /* derive the abs coef, call to function is gf, wl (A), g_low */ Fe2LevN[ipHi][ipLo].opacity = (float)abscf(Fe2LevN[ipHi][ipLo].gf, Fe2LevN[ipHi][ipLo].EnergyWN, Fe2LevN[ipHi][ipLo].gLo); /* excitation energy of transition in degrees kelvin */ Fe2LevN[ipHi][ipLo].EnergyK = (float)(T1CM*Fe2LevN[ipHi][ipLo].EnergyWN); /* energy of photon in ergs */ Fe2LevN[ipHi][ipLo].EnergyErg = (float)(ERG1CM*Fe2LevN[ipHi][ipLo].EnergyWN); } } # ifdef DEBUG_FUN fputs( " <->FeIIPoint()\n", debug_fp ); # endif return; } /* *===================================================================== */ /*FeIIAccel called by forlin to compute radiative acceleration due to FeII lines */ void FeIIAccel(double *fe2drive) { long int ipHi, ipLo; # ifdef DEBUG_FUN fputs( "<+>FeIIAccel()\n", debug_fp ); # endif /*compute acceleration due to Katya's FeII atom */ /* this routine computes the line driven radiative acceleration * due to katya's Fe2 atom*/ *fe2drive = 0.; for( ipLo=0; ipLo < (nFeIILevel - 1); ipLo++ ) { for( ipHi=ipLo+1; ipHi < nFeIILevel; ipHi++ ) { *fe2drive += Fe2LevN[ipHi][ipLo].pump* Fe2LevN[ipHi][ipLo].EnergyErg*Fe2LevN[ipHi][ipLo].PopOpc; } } # ifdef DEBUG_FUN fputs( " <->FeIIAccel()\n", debug_fp ); # endif return; } /* *===================================================================== */ /*RTFeIIMake called by RTMake, does large FeII atom radiative transfer */ void RTFeIIMake( int lgDoEsc ) { long int ipHi, ipLo; # ifdef DEBUG_FUN fputs( "<+>RTFeIIMake()\n", debug_fp ); # endif /* this routine drives calls to make RT relations with Katya's FeII atom */ if( wind.windv == 0. ) { /* static soluiton */ for( ipLo=0; ipLo < (nFeIILevel - 1); ipLo++ ) { for( ipHi=ipLo + 1; ipHi < nFeIILevel; ipHi++ ) { RTMakeStat( &Fe2LevN[ipHi][ipLo] , lgDoEsc ); } } } else { /* windy model */ for( ipLo=0; ipLo < (nFeIILevel - 1); ipLo++ ) { for( ipHi=ipLo + 1; ipHi < nFeIILevel; ipHi++ ) { RTMakeWind(&Fe2LevN[ipHi][ipLo] , lgDoEsc ); } } } # ifdef DEBUG_FUN fputs( " <->RTFeIIMake()\n", debug_fp ); # endif return; } /* *===================================================================== */ /* called in LineSet4 to add FeII lines to save array */ void FeIIAddLines( void ) { long int ipHi, ipLo; # ifdef DEBUG_FUN fputs( "<+>FeIIAddLines()\n", debug_fp ); # endif /* this routine puts the emission from the large FeII atom * into an array to save and integrate them*/ /* add lines option called from lines, add intensities into storage array */ /* routine is called three different ways, ipass < 0 before space allocated, * =0 when time to generate labels (and we zero out array here), and ipass>0 * when time to save intensities */ if( LineSave.ipass == 0 ) { /* we were called by lines, and we want to zero out Fe2SavN */ for( ipLo=0; ipLo < (nFeIILevel - 1); ipLo++ ) { for( ipHi=ipLo + 1; ipHi < nFeIILevel; ipHi++ ) { Fe2SavN[ipHi][ipLo] = 0.; } } } /* this call during calculations, save intensities */ else if( LineSave.ipass == 1 ) { /* total emission from vol */ for( ipLo=0; ipLo < (nFeIILevel - 1); ipLo++ ) { for( ipHi=ipLo + 1; ipHi < nFeIILevel; ipHi++ ) { Fe2SavN[ipHi][ipLo] += (float)radius.dVeff*Fe2LevN[ipHi][ipLo].xIntensity; } } } # ifdef DEBUG_FUN fputs( " <->FeIIAddLines()\n", debug_fp ); # endif return; } /* *===================================================================== */ /*FeIIPunchLines punch accumulated FeII intensities, at end of calculation, called by dopunch to punch them out, punch turned on with punch verner command */ void FeIIPunchLines( /* file we will punch to */ FILE *ioPUN ) { long int MaseHi, MaseLow, ipHi, ipLo; double hbeta, absint , renorm; double TauMase, thresh; # ifdef DEBUG_FUN fputs( "<+>FeIIPunchLines()\n", debug_fp ); # endif /* this routine puts the emission from the large FeII atom * into a line array, and eventually will punch it out */ /* get the normalization line */ if( LineSv[norm.ipNormWavL].sumlin > 0. ) { renorm = norm.ScaleNormLine/LineSv[norm.ipNormWavL].sumlin; } else { renorm = 1.; } fprintf( ioPUN, " up low log I, I/I(norm), Tau\n" ); /* first look for any masing lines */ MaseLow = -1; MaseHi = -1; TauMase = 0.; for( ipLo=0; ipLo < (nFeIILevel - 1); ipLo++ ) { for( ipHi=ipLo + 1; ipHi < nFeIILevel; ipHi++ ) { if( Fe2LevN[ipHi][ipLo].TauIn < TauMase ) { TauMase = Fe2LevN[ipHi][ipLo].TauIn; MaseLow = ipLo; MaseHi = ipHi; } } } if( TauMase < 0. ) { fprintf( ioPUN, " Most negative optical depth was %4ld%4ld%10.2e\n", MaseHi, MaseLow, TauMase ); } /* now print actual line intensities, Hbeta first */ if( !cdLine("TOTL", 4861 , &hbeta , &absint) ) { fprintf( ioQQQ, " pop371 could not find Hbeta with cdLine.\n" ); puts( "[Stop in pop371]" ); exit(1); } fprintf( ioPUN, "Hbeta=%7.3f %10.2e\n", absint + fourpi.pirsq, hbeta ); if( renorm > SMALLFLOAT ) { /* this is threshold for faintest line, normally 0, set with * number on punch verner command */ thresh = FeII.fe2thresh/renorm; } else { thresh = 0.; } for( ipLo=0; ipLo < (nFeIILevel - 1); ipLo++ ) { for( ipHi=ipLo + 1; ipHi < nFeIILevel; ipHi++ ) { /* fe2ener(1) and (2) are lower and upper limits to range of * wavelengths of interest. entered in ryd with * punch feii verner command, where they are converted * to wavenumbers, */ if( (Fe2SavN[ipHi][ipLo] > thresh && Fe2LevN[ipHi][ipLo].EnergyWN > FeII.fe2ener[0]) && Fe2LevN[ipHi][ipLo].EnergyWN < FeII.fe2ener[1] ) { if( FeII.lgShortFe2 ) { /* short printout does not include rel intensity or optical dep */ fprintf( ioPUN, "%3ld%3ld%7.3f\n", ipHi+1, ipLo+1, log10(MAX2(1e-37,Fe2SavN[ipHi][ipLo])) + fourpi.pirsq ); } else { /* long printout does */ fprintf( ioPUN, "%3ld%3ld %8.3f %10.2e %10.2e\n", ipHi+1, ipLo+1, log10(MAX2(1e-37,Fe2SavN[ipHi][ipLo])) + fourpi.pirsq, Fe2SavN[ipHi][ipLo]*renorm , Fe2LevN[ipHi][ipLo].TauIn ); } } } } # ifdef DEBUG_FUN fputs( " <->FeIIPunchLines()\n", debug_fp ); # endif return; } /* *===================================================================== */ /*FeIIEmitOut add large FeII emission to outward beam - called once per zone in RTDiffuse */ void FeIIEmitOut(double VolFac, double ref) { long int ipHi, ilo, ip; # ifdef DEBUG_FUN fputs( "<+>FeIIEmitOut()\n", debug_fp ); # endif for( ilo=0; ilo < (nFeIILevel - 1); ilo++ ) { for( ipHi=ilo + 1; ipHi < nFeIILevel; ipHi++ ) { /* pointer to line energy in continuum array */ ip = Fe2LevN[ipHi][ilo].ipCont; /* 666 error! make sure verner's put this in */ rfield.outlin[ip-1] += (float)(Fe2LevN[ipHi][ilo].phots* VolFac*opac.tmn[ip-1]*Fe2LevN[ipHi][ilo].ColOvTot); rfield.reflin[ip-1] += (float)(Fe2LevN[ipHi][ilo].phots* ref); } } # ifdef DEBUG_FUN fputs( " <->FeIIEmitOut()\n", debug_fp ); # endif return; } /* *===================================================================== */ /*FeIITauInit zero out storage for large FeII atom, called in tauout */ void FeIITauInit(void) { long int ipHi, ipLo; # ifdef DEBUG_FUN fputs( "<+>FeIITauInit()\n", debug_fp ); # endif /* this routine is called in routine zero and it * zero's out various elements of the FeII arrays * it is called on every iteration * */ for( ipLo=0; ipLo < (nFeIILevel - 1); ipLo++ ) { for( ipHi=ipLo + 1; ipHi < nFeIILevel; ipHi++ ) { /* inward optical depth */ Fe2LevN[ipHi][ipLo].TauIn = tauminCom.taumin; Fe2LevN[ipHi][ipLo].TauCon = tauminCom.taumin; Fe2LevN[ipHi][ipLo].ColOvTot = tauminCom.taumin; /* outward optical depth */ Fe2LevN[ipHi][ipLo].TauTot = 1e20f; /* escape probability */ Fe2LevN[ipHi][ipLo].Pesc = 1.; /* inward part of line */ Fe2LevN[ipHi][ipLo].FracInwd = 1.; /* destruction probability */ Fe2LevN[ipHi][ipLo].Pdest = 0.; /* line pumping rate */ Fe2LevN[ipHi][ipLo].pump = 0.; /* population of lower level with correction for stim emission */ Fe2LevN[ipHi][ipLo].PopOpc = 0.; /* population of lower level */ Fe2LevN[ipHi][ipLo].PopLo = 0.; /* following two heat exchange excitation, deexcitation */ Fe2LevN[ipHi][ipLo].cool = 0.; Fe2LevN[ipHi][ipLo].heat = 0.; /* intensity of line */ Fe2LevN[ipHi][ipLo].xIntensity = 0.; /* number of photons emitted in line */ Fe2LevN[ipHi][ipLo].phots = 0.; /* opacity in line */ Fe2LevN[ipHi][ipLo].dTau = 0.; } } # ifdef DEBUG_FUN fputs( " <->FeIITauInit()\n", debug_fp ); # endif return; } /* *===================================================================== */ /* fill in IR lines from lowest 16 levels - these are predicted in Fe2Lev16 */ void FeIIFillLow16(void) { # ifdef DEBUG_FUN fputs( "<+>FeIIFillLow16()\n", debug_fp ); # endif /* this is total cooling due to 16 level atom that would have been computed there */ fe2cool.Fe2L16Tot = 0.; /* all transitions within 16 levels of ground term */ fe2cool.fe21308 = Fe2LevN[12][7].xIntensity; fe2cool.fe21207 = Fe2LevN[11][6].xIntensity; fe2cool.fe21106 = Fe2LevN[10][5].xIntensity; fe2cool.fe21006 = Fe2LevN[9][5].xIntensity; fe2cool.fe21204 = Fe2LevN[11][3].xIntensity; fe2cool.fe21103 = Fe2LevN[10][2].xIntensity; fe2cool.fe21104 = Fe2LevN[10][3].xIntensity; fe2cool.fe21001 = Fe2LevN[9][0].xIntensity; fe2cool.fe21002 = Fe2LevN[9][1].xIntensity; fe2cool.fe20201 = Fe2LevN[1][0].xIntensity; fe2cool.fe20302 = Fe2LevN[2][1].xIntensity; fe2cool.fe20706 = Fe2LevN[6][5].xIntensity; fe2cool.fe20807 = Fe2LevN[7][6].xIntensity; fe2cool.fe20908 = Fe2LevN[8][7].xIntensity; fe2cool.fe21007 = Fe2LevN[9][6].xIntensity; fe2cool.fe21107 = Fe2LevN[10][6].xIntensity; fe2cool.fe21108 = Fe2LevN[10][7].xIntensity; fe2cool.fe21110 = Fe2LevN[10][9].xIntensity; fe2cool.fe21208 = Fe2LevN[11][7].xIntensity; fe2cool.fe21209 = Fe2LevN[11][8].xIntensity; fe2cool.fe21211 = Fe2LevN[11][10].xIntensity; fe2cool.fe21406 = Fe2LevN[13][5].xIntensity; fe2cool.fe21507 = Fe2LevN[14][6].xIntensity; fe2cool.fe21508 = Fe2LevN[14][7].xIntensity; fe2cool.fe21609 = Fe2LevN[15][8].xIntensity; /* these are unique to this large model atom */ if( nFeIILevel > 43 ) { /* NB do not forget to decrement by one when going from physical (in left variables) * to c scale (in array) */ fe2cool.fe25to6 = Fe2LevN[25-1][6-1].xIntensity; fe2cool.fe27to7 = Fe2LevN[27-1][7-1].xIntensity; fe2cool.fe28to8 = Fe2LevN[28-1][8-1].xIntensity; fe2cool.fe29to9 = Fe2LevN[29-1][9-1].xIntensity; fe2cool.fe32to6 = Fe2LevN[32-1][6-1].xIntensity; fe2cool.fe33to7 = Fe2LevN[33-1][7-1].xIntensity; fe2cool.fe37to7 = Fe2LevN[37-1][7-1].xIntensity; fe2cool.fe39to8 = Fe2LevN[39-1][8-1].xIntensity; fe2cool.fe40to9 = Fe2LevN[40-1][9-1].xIntensity; fe2cool.fe37to6 = Fe2LevN[37-1][6-1].xIntensity; fe2cool.fe39to7 = Fe2LevN[39-1][7-1].xIntensity; fe2cool.fe40to8 = Fe2LevN[40-1][8-1].xIntensity; fe2cool.fe41to9 = Fe2LevN[41-1][9-1].xIntensity; fe2cool.fe39to6 = Fe2LevN[39-1][6-1].xIntensity; fe2cool.fe40to7 = Fe2LevN[40-1][7-1].xIntensity; fe2cool.fe41to8 = Fe2LevN[41-1][8-1].xIntensity; fe2cool.fe42to6 = Fe2LevN[42-1][6-1].xIntensity; fe2cool.fe43to7 = Fe2LevN[43-1][7-1].xIntensity; fe2cool.fe42to7 = Fe2LevN[42-1][7-1].xIntensity; fe2cool.fe36to2 = Fe2LevN[36-1][2-1].xIntensity; fe2cool.fe36to3 = Fe2LevN[36-1][3-1].xIntensity; fe2cool.fe32to1 = Fe2LevN[32-1][1-1].xIntensity; fe2cool.fe33to2 = Fe2LevN[33-1][2-1].xIntensity; fe2cool.fe36to5 = Fe2LevN[36-1][5-1].xIntensity; fe2cool.fe32to2 = Fe2LevN[32-1][2-1].xIntensity; fe2cool.fe33to3 = Fe2LevN[33-1][3-1].xIntensity; fe2cool.fe30to3 = Fe2LevN[30-1][3-1].xIntensity; fe2cool.fe33to6 = Fe2LevN[33-1][6-1].xIntensity; fe2cool.fe24to2 = Fe2LevN[24-1][2-1].xIntensity; fe2cool.fe32to7 = Fe2LevN[32-1][7-1].xIntensity; fe2cool.fe35to8 = Fe2LevN[35-1][8-1].xIntensity; fe2cool.fe34to8 = Fe2LevN[34-1][8-1].xIntensity; fe2cool.fe27to6 = Fe2LevN[27-1][6-1].xIntensity; fe2cool.fe28to7 = Fe2LevN[28-1][7-1].xIntensity; fe2cool.fe30to8 = Fe2LevN[30-1][8-1].xIntensity; fe2cool.fe24to6 = Fe2LevN[24-1][6-1].xIntensity; fe2cool.fe29to8 = Fe2LevN[29-1][8-1].xIntensity; fe2cool.fe24to7 = Fe2LevN[24-1][7-1].xIntensity; fe2cool.fe22to7 = Fe2LevN[22-1][7-1].xIntensity; fe2cool.fe38to11 =Fe2LevN[38-1][11-1].xIntensity; fe2cool.fe19to8 = Fe2LevN[19-1][8-1].xIntensity; fe2cool.fe17to6 = Fe2LevN[17-1][6-1].xIntensity; fe2cool.fe18to7 = Fe2LevN[18-1][7-1].xIntensity; fe2cool.fe18to8 = Fe2LevN[18-1][8-1].xIntensity; fe2cool.fe17to7 = Fe2LevN[17-1][7-1].xIntensity; } else { fe2cool.fe25to6 = 0.; fe2cool.fe27to7 = 0.; fe2cool.fe28to8 = 0.; fe2cool.fe29to9 = 0.; fe2cool.fe32to6 = 0.; fe2cool.fe33to7 = 0.; fe2cool.fe37to7 = 0.; fe2cool.fe39to8 = 0.; fe2cool.fe40to9 = 0.; fe2cool.fe37to6 = 0.; fe2cool.fe39to7 = 0.; fe2cool.fe40to8 = 0.; fe2cool.fe41to9 = 0.; fe2cool.fe39to6 = 0.; fe2cool.fe40to7 = 0.; fe2cool.fe41to8 = 0.; fe2cool.fe42to6 = 0.; fe2cool.fe43to7 = 0.; fe2cool.fe42to7 = 0.; fe2cool.fe36to2 = 0.; fe2cool.fe36to3 = 0.; fe2cool.fe32to1 = 0.; fe2cool.fe33to2 = 0.; fe2cool.fe36to5 = 0.; fe2cool.fe32to2 = 0.; fe2cool.fe33to3 = 0.; fe2cool.fe30to3 = 0.; fe2cool.fe33to6 = 0.; fe2cool.fe24to2 = 0.; fe2cool.fe32to7 = 0.; fe2cool.fe35to8 = 0.; fe2cool.fe34to8 = 0.; fe2cool.fe27to6 = 0.; fe2cool.fe28to7 = 0.; fe2cool.fe30to8 = 0.; fe2cool.fe24to6 = 0.; fe2cool.fe29to8 = 0.; fe2cool.fe24to7 = 0.; fe2cool.fe22to7 = 0.; fe2cool.fe38to11 =0.; fe2cool.fe19to8 = 0.; fe2cool.fe17to6 = 0.; fe2cool.fe18to7 = 0.; fe2cool.fe18to8 = 0.; fe2cool.fe17to7 = 0.; } if( nFeIILevel > 80 ) { fe2cool.fe80to28 = Fe2LevN[80-1][28-1].xIntensity; } else { fe2cool.fe80to28 =0.; } # ifdef DEBUG_FUN fputs( " <->FeIIFillLow16()\n", debug_fp ); # endif return; } /* *===================================================================== */ void FeIIPunDepart( /* io unit for punch */ FILE* ioPUN , /* punch all levels if true, only subset if false */ int lgDoAll ) { /* numer of levels with dep coef that we will punch out */ # define NLEVDEP 11 /* these are the levels on the physical, not c, scale (count from 1) */ const int LevDep[NLEVDEP]={1,10,25,45,64,124,206,249,295,347,371}; long int i; static int lgFIRST=TRUE; # ifdef DEBUG_FUN fputs( "<+>FeIIPunDepart()\n", debug_fp ); # endif /* on first call only, print levels that we will punch later */ if( lgFIRST && !lgDoAll ) { /* but all the rest do */ for( i=0; iFeIIPunDepart()\n", debug_fp ); # endif return; } /* *===================================================================== */ void FeIIPun1Depart( /* the io unit where the print should be directed */ FILE * ioPUN , /* the physical (not c) number of the level */ long int nPUN ) { # ifdef DEBUG_FUN fputs( "<+>FeIIPun1Depart()\n", debug_fp ); # endif assert( nPUN > 0 ); assert( ioPUN != NULL ); /* print the level departure coef on ioPUN if the level was computed, * print a zero if it was not */ if( nPUN < nFeIILevel ) { fprintf( ioPUN , "%e ",Fe2DepCoef[nPUN-1] ); } else { fprintf( ioPUN , "%e ",0. ); } # ifdef DEBUG_FUN fputs( " <->FeIIPun1Depart()\n", debug_fp ); # endif return; } /* *===================================================================== */ void FeIIReset(void) { # ifdef DEBUG_FUN fputs( "<+>FeIIReset()\n", debug_fp ); # endif /* space has been allocated, so reset number of levels to whatever it was */ nFeIILevel = nFeIILevelAlloc; # ifdef DEBUG_FUN fputs( " <->FeIIReset()\n", debug_fp ); # endif return; } /* *===================================================================== */ void FeIIZero(void) { # ifdef DEBUG_FUN fputs( "<+>FeIIZero()\n", debug_fp ); # endif /* flag saying that lya pumping of feii in large atom is turned on */ FeII.lgLyaPumpOn = TRUE; /* will not compute large feii atom */ FeII.lgFeIION = FALSE; /* energy range of large FeII atom is zero to infinity */ FeII.fe2ener[0] = 0.; FeII.fe2ener[1] = 1e8; /* set zero for the threshold of weakest large FeII atom line to punch */ FeII.fe2thresh = 0.; /* normally do not constantly reevaluate the atom, set TRUE with * SLOW key on atom feii command */ FeII.lgSlow = FALSE; /* option to print each call to pop371, set with print option on atom feii */ FeII.lgPrint = FALSE; /* option to only simulate calls to pop371 */ FeII.lgSimulate = FALSE; /* set number of levels for the atom */ /* number of levels for the large FeII atom, changed with the atom feii levels command */ if( lgFeIIMalloc ) { /* space has been allocated, so reset number of levels to whatever it was */ nFeIILevel = nFeIILevelAlloc; } else { /* space not allocated yet, set to largest possible number of levels */ nFeIILevel = NFE2LEVN; } # ifdef DEBUG_FUN fputs( " <->FeIIZero()\n", debug_fp ); # endif return; } /* *===================================================================== */ void FeIIPunPop( /* io unit for punch */ FILE* ioPUN , /* punch all levels if true, only subset if false */ int lgDoAll ) { /* numer of levels with dep coef that we will punch out */ # define NLEVPOP 11 /* these are the levels on the physical, not c, scale (count from 1) */ const int LevPop[NLEVPOP]={1,10,25,45,64,124,206,249,295,347,371}; long int i; static int lgFIRST=TRUE; # ifdef DEBUG_FUN fputs( "<+>FeIIPunPop()\n", debug_fp ); # endif /* on first call only, print levels that we will punch later, * but only if we will only punch selected levels*/ if( lgFIRST && !lgDoAll ) { /* but all the rest do */ for( i=0; iFeIIPunPop()\n", debug_fp ); # endif return; } /* *===================================================================== */ void FeIIPun1Pop( /* the io unit where the print should be directed */ FILE * ioPUN , /* the physical (not c) number of the level */ long int nPUN ) { # ifdef DEBUG_FUN fputs( "<+>FeIIPun1Pop()\n", debug_fp ); # endif assert( nPUN > 0 ); assert( ioPUN != NULL ); /* print the level departure coef on ioPUN if the level was computed, * print a zero if it was not */ if( nPUN < nFeIILevel ) { fprintf( ioPUN , "%e ",Fe2LevelPop[nPUN-1] ); } else { fprintf( ioPUN , "%e ",0. ); } # ifdef DEBUG_FUN fputs( " <->FeIIPun1Pop()\n", debug_fp ); # endif return; } /* *===================================================================== */ long int cdGetFeIIBands(void) { char chLine[132] , chFilename[134] ; FILE *ioDATA; int lgEOL; long int i,k; /* keep track of whether we have been called - want to be * called a total of one time */ static int lgCalled=FALSE; # ifdef DEBUG_FUN fputs( "<+>cdGetFeIIBands()\n", debug_fp ); # endif /* return previous number of bands if this is second or later call*/ if( lgCalled ) { # ifdef DEBUG_FUN fputs( " <->cdGetFeIIBands()\n", debug_fp ); # endif /* return value of number of bands, may be used by calling program*/ return FeIInBands; } lgCalled = TRUE; /* get FeII band data */ /* check on path if path set */ /* path was parsed in getset */ if( lgDataPathSet == TRUE ) { /*path set, so look only there */ strcpy( chFilename , chDataPath ); strcat( chFilename , "fe2bands.dat" ); } else { /* path not set, check local space only */ strcpy( chFilename , "fe2bands.dat" ); } if( trace.lgTrace ) { fprintf( ioQQQ," FeIIData opening fe2bands.dat:"); } if( ( ioDATA = fopen( chFilename , "r" ) ) == NULL ) { fprintf( ioQQQ, " FeIIData could not open fe2bands.dat\n" ); if( lgDataPathSet == TRUE ) { fprintf( ioQQQ, " FeIIData could not open fe2bands.dat, even tried path.\n"); fprintf( ioQQQ, " path is *%s*\n",chDataPath ); fprintf( ioQQQ, " final path is *%s*\n",chFilename ); } puts( "[Stop in FeIIData]" ); exit(-1); } /* now count how many bands are in the file */ FeIInBands = 0; /* first line is a version number and does not count */ if( fgets( chLine , sizeof(chLine) , ioDATA ) == NULL ) { fprintf( ioQQQ, " FeIIData could not read first line of fe2bands.dat.\n"); puts( "[Stop in FeIIData]" ); exit(-1); } while( fgets( chLine , sizeof(chLine) , ioDATA ) != NULL ) { /* we want to count the lines that do not start with # * since these contain data */ if( chLine[0] != '#') ++FeIInBands; } /* now rewind the file so we can read it a second time*/ if( fseek( ioDATA , 0 , SEEK_SET ) != 0 ) { fprintf( ioQQQ, " FeIIData could not rewind fe2bands.dat.\n"); puts( "[Stop in FeIIData]" ); exit(-1); } if( (FeII_Bands = (float **)malloc(sizeof(float *)*(unsigned)(FeIInBands) ) )==NULL ) { fprintf( ioQQQ, " FeIIData could not malloc FeII_Bands 1\n" ); puts( "[Stop in FeIIData]" ); exit(1); } /* now make second dim, id wavelength, and lower and upper bounds */ for( i=0; i= FeII_Bands[i][2] ) { fprintf( ioQQQ, " FeII band %li has improper bounds.\n" ,i); puts( "[Stop in FeIIData]" ); exit(1); } } fclose(ioDATA); # ifdef DEBUG_FUN fputs( " <->cdGetFeIIBands()\n", debug_fp ); # endif /* return value of number of bands, may be used by calling program*/ return FeIInBands; }