183 if( strcmp(chLabel,
"ZERO") == 0 )
186 for( n=0; n <
FeII.nFeIILevel_malloc; ++n )
193 else if( strcmp(chLabel,
"ADD ") == 0 )
196 for( n=0; n <
FeII.nFeIILevel_local; ++n )
204 else if( strcmp(chLabel,
"PRIN") != 0 )
206 fprintf(
ioQQQ,
" FeII_Colden does not understand the label %s\n",
251 FeII.nFeIILevel_malloc =
FeII.nFeIILevel_local;
254 Fe2A = (
double **)
MALLOC(
sizeof(
double *)*(
unsigned long)
FeII.nFeIILevel_malloc );
257 for( ipHi=0; ipHi <
FeII.nFeIILevel_malloc; ++ipHi )
259 Fe2A[ipHi]=(
double*)
MALLOC(
sizeof(
double )*(
unsigned long)
FeII.nFeIILevel_malloc);
263 Fe2CPump = (
double **)
MALLOC(
sizeof(
double *)*(
unsigned long)
FeII.nFeIILevel_malloc );
266 Fe2LPump = (
double **)
MALLOC(
sizeof(
double *)*(
unsigned long)
FeII.nFeIILevel_malloc );
269 for( ipHi=0; ipHi <
FeII.nFeIILevel_malloc; ++ipHi )
271 Fe2CPump[ipHi]=(
double*)
MALLOC(
sizeof(
double )*(
unsigned long)
FeII.nFeIILevel_malloc);
273 Fe2LPump[ipHi]=(
double*)
MALLOC(
sizeof(
double )*(
unsigned long)
FeII.nFeIILevel_malloc);
283 for( ipHi=0; ipHi <
FeII.nFeIILevel_malloc; ++ipHi )
289 xMatrix = (
double **)
MALLOC(
sizeof(
double *)*(
unsigned long)
FeII.nFeIILevel_malloc );
292 for( i=0; i<
FeII.nFeIILevel_malloc; ++i )
294 xMatrix[i] = (
double *)
MALLOC(
sizeof(
double)*(
unsigned long)
FeII.nFeIILevel_malloc );
297 amat=(
double*)
MALLOC( (
sizeof(
double)*(
unsigned long)(
FeII.nFeIILevel_malloc*
FeII.nFeIILevel_malloc) ));
300 yVector=(
double*)
MALLOC( (
sizeof(
double)*(
unsigned long)(
FeII.nFeIILevel_malloc) ));
303 Fe2SavN = (
double **)
MALLOC(
sizeof(
double *)*(
unsigned long)
FeII.nFeIILevel_malloc );
306 for( ipHi=1; ipHi <
FeII.nFeIILevel_malloc; ++ipHi )
308 Fe2SavN[ipHi]=(
double*)
MALLOC(
sizeof(
double )*(
unsigned long)ipHi);
333 for(ipHi=1; ipHi <
NPRADDAT; ipHi++)
339 for( ipLo=0; ipLo< ipHi; ipLo++ )
346 for(ipHi=0; ipHi <
NPRADDAT; ipHi++)
348 for( ipLo=0; ipLo< ipHi; ipLo++ )
358 for( ipHi=0; ipHi <
FeII.nFeIILevel_malloc; ipHi++ )
360 for( ipLo=0; ipLo <
FeII.nFeIILevel_malloc; ipLo++ )
362 FeII.FeIIAul[ipHi][ipLo] = -2.;
363 for(
int k=0; k<8; k++ )
365 FeII.FeIIColl[ipHi][ipLo][k] = -2.;
368 FeII.FeIINRGs[ipHi] = -2.;
369 FeII.FeIISTWT[ipHi] = -2.;
374 ncs1=(
int**)
MALLOC(
sizeof(
int *)*(
unsigned long)
FeII.nFeIILevel_malloc);
376 for( ipHi=1; ipHi <
FeII.nFeIILevel_malloc; ++ipHi )
379 ncs1[ipHi]=(
int*)
MALLOC(
sizeof(
int)*(
unsigned long)ipHi);
388 for( ipHi=1; ipHi <
FeII.nFeIILevel_malloc; ipHi++ )
390 for( ipLo=0; ipLo < ipHi; ipLo++ )
401 for( ipHi=1; ipHi <
FeII.nFeIILevel_malloc; ipHi++ )
404 for( ipLo=0; ipLo < ipHi; ipLo++ )
416 fprintf(
ioQQQ,
" FeIICreate opening fe2nn.dat:");
428 for( i=0; i < 8; i++ )
432 fprintf(
ioQQQ,
" fe2nn.dat error reading data\n" );
441 "%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld",
448 for( m1=0; m1<20; ++m1 )
459 fprintf(
ioQQQ,
" FeIICreate opening fe2energies.dat:");
462 ioDATA =
open_data(
"fe2energies.dat",
"r" );
465 for( ipHi=0; ipHi <
FeII.nFeIILevel_malloc; ipHi++ )
470 while( chLine[0] ==
'#' )
474 fprintf(
ioQQQ,
" fe2energies.dat error reading data\n" );
481 sscanf( chLine,
"%lf", &help );
490 fprintf(
ioQQQ,
" FeIICreate opening fe2rad.dat:");
497 fprintf(
ioQQQ,
" fe2rad.dat error reading data\n" );
501 sscanf( chLine ,
"%ld%ld%ld",&lo, &ihi, &m1);
502 const int nYrFe2Rad=8, nMoFe2Rad=8, nDyFe2Rad=24;
503 if( lo!=nYrFe2Rad || ihi!=nMoFe2Rad || m1!=nDyFe2Rad )
505 fprintf(
ioQQQ,
"DISASTER fe2rad.dat has the wrong magic numbers, expected "
506 "%2i %2i %2i and got %2li %2li %2li\n" ,
507 nYrFe2Rad, nMoFe2Rad, nDyFe2Rad,
513 for( ipHi=1; ipHi <
FeII.nFeIILevel_malloc; ipHi++ )
515 for( ipLo=0; ipLo < ipHi; ipLo++ )
522 while( chLine[0] ==
'#' )
526 fprintf(
ioQQQ,
" fe2nn.dat error reading data\n" );
533 "%ld%ld%ld%ld%lf%lf%ld",
534 &lo, &ihi, &m1, &m2 ,
540 FeII.FeIISTWT[lo-1] = (*tr.
Lo()).g();
542 FeII.FeIISTWT[ihi-1] = (*tr.
Hi()).g();
569 ncs1[ihi-1][lo-1] = (int)m3;
584 fprintf(
ioQQQ,
" FeIICreate opening fe2col.dat:");
590 for( ipHi=1; ipHi<
NPRADDAT; ++ipHi )
592 for( ipLo=0; ipLo<ipHi; ++ipLo )
598 fprintf(
ioQQQ,
" fe2col.dat error reading data\n" );
602 "%ld%ld%lf%lf%lf%lf%lf%lf%lf%lf",
616 FeII.FeIIColl[ipHiFe2-1][ipLoFe2-1][k] =
save[k];
632 for( ipLo=0; ipLo < (
FeII.nFeIILevel_malloc - 1); ipLo++ )
634 for( ipHi=ipLo + 1; ipHi <
FeII.nFeIILevel_malloc; ipHi++ )
655 (*tr.
Hi()).IonStg() = 2;
656 (*tr.
Hi()).nelem() = 26;
684 fprintf(
ioQQQ,
" failed to open FeII bands file \n");
739 fprintf(
ioQQQ,
" FeIILevelPops fe2 pops called\n");
744 if(
FeII.lgSimulate )
749 for( n=0; n<
FeII.nFeIILevel_local; ++n)
751 for( ipLo=0; ipLo<
FeII.nFeIILevel_local; ++ipLo )
766 for( n=0; n<
FeII.nFeIILevel_local; ++n)
768 for( ipHi=n+1; ipHi<
FeII.nFeIILevel_local; ++ipHi )
776 (*tr.
Hi()).g()/(*tr.
Lo()).g();
794 for( n=0; n<
FeII.nFeIILevel_local; ++n)
797 for( ipLo=0; ipLo<n; ++ipLo )
803 for( ipHi=n+1; ipHi<
FeII.nFeIILevel_local; ++ipHi )
808 for( ipLo=0; ipLo<n; ++ipLo )
813 for( ipHi=n+1; ipHi<
FeII.nFeIILevel_local; ++ipHi )
825 enum {DEBUG_LOC=
false};
829 for( n=0; n<
FeII.nFeIILevel_local; ++n)
833 for( ipLo=0; ipLo<
FeII.nFeIILevel_local; ++ipLo )
846 for( n=1; n <
FeII.nFeIILevel_local; n++ )
856# define AMAT(I_,J_) (*(amat+(I_)*FeII.nFeIILevel_local+(J_)))
859 for( ipHi=0; ipHi <
FeII.nFeIILevel_local; ipHi++ )
861 for( i=0; i <
FeII.nFeIILevel_local; i++ )
875 fprintf(
ioQQQ,
"DISASTER FeIILevelPops: dgetrs finds singular or ill-conditioned matrix\n" );
884 for( ipLo=0; ipLo <
FeII.nFeIILevel_local; ipLo++ )
889 fprintf(
ioQQQ,
"PROBLEM FeIILevelPops finds non-positive level population, level is %ld pop is %g\n",
899 fprintf(
ioQQQ ,
"PROBLEM FeIILevelPops exits with negative level populations.\n");
906 for( ipLo=
FeII.nFeIILevel_local; ipLo <
FeII.nFeIILevel_malloc; ++ipLo )
914 for( ipLo=0; ipLo < (
FeII.nFeIILevel_local - 1); ipLo++ )
918 for( ipHi=ipLo+1; ipHi <
FeII.nFeIILevel_malloc; ipHi++ )
947 if(
FeII.lgFeIILargeOn )
950 hydro.dstfe2lya = 0.;
953 for( ipLo=0; ipLo < (
FeII.nFeIILevel_local - 1); ipLo++ )
955 for( ipHi=ipLo+1; ipHi <
FeII.nFeIILevel_local; ipHi++ )
973 hydro.dstfe2lya = 0.;
980 enum {DEBUG_LOC=
false};
984 fprintf(
ioQQQ,
" sum all %.1e dest rate%.1e escR= %.1e\n",
985 EnergyWN,
hydro.dstfe2lya,
996 for( i=1; i <
FeII.nFeIILevel_local; i++ )
1006 for( i=0; i <
FeII.nFeIILevel_local; i++ )
1025 FeII.Fe2_large_cool = 0.0f;
1027 FeII.Fe2_large_heat = 0.f;
1029 FeII.ddT_Fe2_large_cool = 0.0f;
1032 for( ipLo=0; ipLo < (
FeII.nFeIILevel_local - 1); ipLo++ )
1034 for( ipHi=ipLo + 1; ipHi <
FeII.nFeIILevel_local; ipHi++ )
1075 static double OldTemp = -1.;
1085 realnum FracLowTe , FracHighTe;
1086 static realnum tt[8]={1e3f,3e3f,5e3f,7e3f,1e4f,12e3f,15e3f,2e4f};
1091 static long int nFeII_old = -1;
1103 nFeII_old =
FeII.nFeIILevel_local;
1108 for( ipLo=0; ipLo < (
FeII.nFeIILevel_local - 1); ipLo++ )
1110 for( ipHi=ipLo + 1; ipHi <
FeII.nFeIILevel_local; ipHi++ )
1113 if(
ncs1[ipHi][ipLo] == 3 )
1121 else if(
ncs1[ipHi][ipLo] == 1 )
1128 else if(
ncs1[ipHi][ipLo] == 2 )
1140 fprintf(
ioQQQ,
">>>PROBLEM impossible ncs1 in FeIILevelPops\n");
1147 gb = (
realnum)(ag + (-cg*
POW2(y) + dg)*(log(1.0+1.0/y) - 0.4/
1148 POW2(y + 1.0)) + cg*y);
1183 else if(
phycon.te > tt[7] )
1197 for( i=0; i < 8-1; i++ )
1199 if(
phycon.te <= tt[i+1] )
1209 fprintf(
ioQQQ,
" Insanity while looking for temperature in coll str array, te=%g.\n",
1217 FracHighTe = ((
realnum)
phycon.te - tt[ipTemp])/(tt[ipTempp1] - tt[ipTemp]);
1222 FracLowTe = 1.f - FracHighTe;
1224 for( ipHi=1; ipHi <
NPRADDAT; ipHi++ )
1226 for( ipLo=0; ipLo < ipHi; ipLo++ )
1233 ASSERT( ipHiFe2-1 >= 0 );
1235 ASSERT( ipLoFe2-1 >= 0 );
1238 if( ipHiFe2-1 <
FeII.nFeIILevel_local )
1244 sPradDat[ipHi][ipLo][ipTemp] * FracLowTe +
1245 sPradDat[ipHi][ipLo][ipTempp1] * FracHighTe;
1256 for( ipHi=1; ipHi <
FeII.nFeIILevel_local; ipHi++ )
1266 ipTrim =
FeII.nFeIILevel_local - 1;
1279 if( ipTrim+1 <
FeII.nFeIILevel_local )
1282 for( ipLo=0; ipLo<
FeII.nFeIILevel_local; ++ipLo )
1284 for( ipHi=ipTrim; ipHi<
FeII.nFeIILevel_local; ++ipHi )
1291 FeII.nFeIILevel_local = ipTrim+1;
1296 enum {DEBUG_LOC=
false};
1299 for( ipLo=0; ipLo < 52; ipLo++ )
1301 fprintf(
ioQQQ,
"%e %e\n",
1309 for( ipLo=0; ipLo<
FeII.nFeIILevel_local; ++ipLo)
1311 for( ipHi=ipLo+1; ipHi<
FeII.nFeIILevel_local; ++ipHi )
1325 (*tr.
Hi()).g()/(*tr.
Lo()).g());
1359 double *SumBandInward )
1363 double SumBandFe2_v;
1369 *SumBandInward = 0.;
1374 for( ipHi=1; ipHi <
FeII.nFeIILevel_local; ++ipHi )
1376 for( ipLo=0; ipLo < ipHi; ++ipLo )
1379 if( tr.
WLAng() >= wl1 &&
1389 return( SumBandFe2_v );
1403 for( ipLo=0; ipLo < (
FeII.nFeIILevel_local - 1); ipLo++ )
1409 for( ipHi=ipLo + 1; ipHi <
FeII.nFeIILevel_malloc; ipHi++ )
1441 for( ipHi=1; ipHi <
FeII.nFeIILevel_malloc; ipHi++ )
1443 for( ipLo=0; ipLo < ipHi; ipLo++ )
1466 for( ipLo=0; ipLo <
FeII.nFeIILevel_malloc-1; ipLo++ )
1468 for( ipHi=ipLo+1; ipHi <
FeII.nFeIILevel_malloc; ipHi++ )
1473 if( fabs(tr.
Emis().
Aul()- 1e-5 ) > 1e-8 )
1479 if( strcmp(
rfield.chLineLabel[ip-1],
" " ) == 0 )
1480 strcpy(
rfield.chLineLabel[ip-1],
"FeII" );
1521 for( ipLo=0; ipLo < (
FeII.nFeIILevel_local - 1); ipLo++ )
1523 for( ipHi=ipLo+1; ipHi <
FeII.nFeIILevel_local; ipHi++ )
1547 fprintf(
ioQQQ,
" FeII_RT_Make called\n");
1551 for( ipLo=0; ipLo < (
FeII.nFeIILevel_local - 1); ipLo++ )
1554 for( ipHi=ipLo + 1; ipHi <
FeII.nFeIILevel_malloc; ipHi++ )
1590 for(
long ipLo=0; ipLo < (
FeII.nFeIILevel_malloc - 1); ipLo++ )
1592 for(
long ipHi=ipLo + 1; ipHi <
FeII.nFeIILevel_malloc; ipHi++ )
1603 for(
long ipLo=0; ipLo < (
FeII.nFeIILevel_local - 1); ipLo++ )
1605 for(
long ipHi=ipLo + 1; ipHi <
FeII.nFeIILevel_local; ipHi++ )
1614 enum {DEBUG_LOC=
false};
1634 fprintf(ioPUN ,
"%.2f\t%li\n", 0., (
long)(*
Fe2LevN[
ipFe2LevN[1][0]].Lo()).
g() );
1635 for( ipHi=1; ipHi <
FeII.nFeIILevel_malloc; ++ipHi )
1638 fprintf(ioPUN ,
"%.2f\t%li\n",
1640 (
long)(*tr.
Hi()).g());
1657 for( n=1; n <
FeII.nFeIILevel_malloc; ++n )
1660 fprintf(ioPUN ,
"%.2f\t%li\t%.3e\n",
1662 (
long)(*tr.
Hi()).g(),
1686 for( ipLo=0; ipLo < (
FeII.nFeIILevel_malloc - 1); ipLo++ )
1688 for( ipHi=ipLo + 1; ipHi <
FeII.nFeIILevel_malloc; ipHi++ )
1695 fprintf( ioPUN,
"%ld\t%ld\t%.2f\t%.2e\n",
1717 double hbeta, absint , renorm;
1733 fprintf( ioPUN,
" up low log I, I/I(LineSave), Tau\n" );
1739 for( ipLo=0; ipLo < (
FeII.nFeIILevel_malloc - 1); ipLo++ )
1741 for( ipHi=ipLo + 1; ipHi <
FeII.nFeIILevel_malloc; ipHi++ )
1755 fprintf( ioPUN,
" Most negative optical depth was %4ld%4ld%10.2e\n",
1756 MaseHi, MaseLow, TauMase );
1760 if(
cdLine(
"TOTL", 4861.36f , &hbeta , &absint)<=0 )
1762 fprintf(
ioQQQ,
" FeIILevelPops could not find Hbeta with cdLine.\n" );
1766 fprintf( ioPUN,
"Hbeta=%7.3f %10.2e\n",
1779 for( ipLo=0; ipLo < (
FeII.nFeIILevel_malloc - 1); ipLo++ )
1781 for( ipHi=ipLo + 1; ipHi <
FeII.nFeIILevel_malloc; ipHi++ )
1788 if( (
Fe2SavN[ipHi][ipLo] > thresh &&
1792 if(
FeII.lgShortFe2 )
1796 fprintf( ioPUN,
"%ld\t%ld\t%.2f\t%.3f\n",
1805 fprintf( ioPUN,
"%ld\t%ld\t%.2f\t%.3f\t%.2e\t%.2e\n",
1836 for( ipLo=0; ipLo < (
FeII.nFeIILevel_malloc - 1); ipLo++ )
1838 for( ipHi=ipLo + 1; ipHi <
FeII.nFeIILevel_malloc; ipHi++ )
1862 for( ipHi=1; ipHi <
FeII.nFeIILevel_malloc; ipHi++ )
1866 for( ipLo=0; ipLo < ipHi; ipLo++ )
1874 (*tr.
Lo()).Pop() = 0.;
1877 (*tr.
Hi()).Pop() = 0.;
1960 long int ipLo , ipHi;
1967 " FeIIPunData ALL option not implemented yet 1\n" );
1972 long int nSkip=0, limit=
MIN2(64,
FeII.nFeIILevel_local);
1977 bool lgPrint =
true;
1978 for( ipHi=1; ipHi<limit; ++ipHi )
1980 for( ipLo=0; ipLo<ipHi; ++ipLo )
1986 fprintf( ioPUN ,
"\n");
1993 for( ipHi=limit; ipHi<
FeII.nFeIILevel_local; ++ipHi )
1997 for( ipLo=0; ipLo<ipHi; ++ipLo )
2005 if(
ncs1[ipHi][ipLo] == 3 &&
2006 (fabs(tr.
Emis().
Aul()-1e-5) < 1e-8 ) )
2019 fprintf( ioPUN ,
" %li lines skipped\n",nSkip);
2038 const int LevDep[
NLEVDEP]={1,10,25,45,64,124,206,249,295,347,371};
2040 static bool lgFIRST=
true;
2045 if( lgFIRST && !lgDoAll )
2050 fprintf( ioPUN ,
"%i\t", LevDep[i] );
2052 fprintf( ioPUN ,
"\n");
2059 for( i=1; i<=
FeII.nFeIILevel_local; ++i )
2062 fprintf( ioPUN ,
"\n");
2071 fprintf( ioPUN ,
"\t");
2073 fprintf( ioPUN ,
"\n");
2093 assert( ioPUN != NULL );
2097 if( nPUN <=
FeII.nFeIILevel_local )
2103 fprintf( ioPUN ,
"%e ",0. );
2118 FeII.nFeIILevel_local =
FeII.nFeIILevel_malloc;
2133 FeII.Fe2_large_cool = 0.;
2134 FeII.ddT_Fe2_large_cool = 0.;
2135 FeII.Fe2_large_heat = 0.;
2138 FeII.lgLyaPumpOn =
true;
2145 FeII.lgFeIILargeOn =
false;
2148 FeII.fe2ener[0] = 0.;
2149 FeII.fe2ener[1] = 1e8;
2158 FeII.fe2thresh = 0.;
2162 FeII.lgSlow =
false;
2165 FeII.lgPrint =
false;
2168 FeII.lgSimulate =
false;
2175 FeII.nFeIILevel_local =
FeII.nFeIILevel_malloc;
2182 FeII.nFeIILevel_local = 16;
2195 long int ipRangeLo ,
2197 long int ipRangeHi ,
2199 bool lgPunchDensity )
2204 const int LevPop[
NLEVPOP]={1,10,25,45,64,124,206,249,295,347,371};
2206 static bool lgFIRST=
true;
2212 if( !lgPunchDensity )
2217 if( lgFIRST && !lgPunchRange )
2223 fprintf( ioPUN ,
"%i\t", LevPop[i] );
2225 fprintf( ioPUN ,
"\n");
2232 ASSERT( ipRangeLo>=0 && ipRangeLo<ipRangeHi );
2236 long nHigh =
MIN2(ipRangeHi,
FeII.nFeIILevel_local );
2237 for( i=ipRangeLo; i<nHigh; ++i )
2240 fprintf( ioPUN ,
"%.3e\t",
Fe2LevelPop[i]/denominator );
2242 fprintf( ioPUN ,
"\n");
2250 fprintf( ioPUN ,
"%.3e\t",
Fe2LevelPop[LevPop[i]-1]/denominator );
2252 fprintf( ioPUN ,
"\n");
2272 assert( ioPUN != NULL );
2277 if( nPUN <=
FeII.nFeIILevel_local )
2283 fprintf( ioPUN ,
"%e ",0. );
2297 const char chFile[] )
2301 const char* chFile1;
2322 chFile1 = ( strlen(chFile) == 0 ) ?
"FeII_bands.ini" : chFile;
2327 fprintf(
ioQQQ,
" FeIICreate opening %s:", chFile1 );
2336 if(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioDATA ) == NULL )
2338 fprintf(
ioQQQ,
" FeIICreate could not read first line of %s.\n", chFile1 );
2341 while(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioDATA ) != NULL )
2345 if( chLine[0] !=
'#')
2350 if( fseek( ioDATA , 0 , SEEK_SET ) != 0 )
2352 fprintf(
ioQQQ,
" FeIICreate could not rewind %s.\n", chFile1 );
2365 if(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioDATA ) == NULL )
2367 fprintf(
ioQQQ,
" FeIICreate could not read first line of %s.\n", chFile1 );
2372 const long int iyr = 9, imo=6 , idy = 11;
2373 long iyrread, imoread , idyread;
2374 iyrread = (long)
FFmtRead(chLine,&i,
sizeof(chLine),&lgEOL);
2375 imoread = (long)
FFmtRead(chLine,&i,
sizeof(chLine),&lgEOL);
2376 idyread = (long)
FFmtRead(chLine,&i,
sizeof(chLine),&lgEOL);
2378 if(( iyrread != iyr ) ||
2379 ( imoread != imo ) ||
2380 ( idyread != idy ) )
2383 " PROBLEM FeIIBandsCreate: the version of %s is not the current version.\n", chFile1 );
2385 " FeIIBandsCreate: I expected the magic numbers %li %li %li but found %li %li %li.\n",
2386 iyr, imo , idy ,iyrread, imoread , idyread );
2393 while(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioDATA ) != NULL )
2397 if( chLine[0] !=
'#')
2403 fprintf(
ioQQQ,
" There should have been a number on this band line 1. Sorry.\n" );
2404 fprintf(
ioQQQ,
"string==%s==\n" ,chLine );
2410 fprintf(
ioQQQ,
" There should have been a number on this band line 2. Sorry.\n" );
2411 fprintf(
ioQQQ,
"string==%s==\n" ,chLine );
2417 fprintf(
ioQQQ,
" There should have been a number on this band line 3. Sorry.\n" );
2418 fprintf(
ioQQQ,
"string==%s==\n" ,chLine );
2432 fprintf(
ioQQQ,
" FeII band %li has none positive entry.\n",i );
2438 fprintf(
ioQQQ,
" FeII band %li has improper bounds.\n" ,i);
2478 for( n=0; n<
FeII.nFeIILevel_local; ++n )
2486 *BigError =
MAX2( *BigError , error );
2493 arg = sum2 -
POW2( *pred ) / (double)
FeII.nFeIILevel_local;
2495 *StdDev = sqrt( arg / (
double)(
FeII.nFeIILevel_local - 1.) );
2498 *pred /= (double)(
FeII.nFeIILevel_local);
2511 for( ipLo=0; ipLo < (
FeII.nFeIILevel_local - 1); ipLo++ )
2513 for( ipHi=ipLo + 1; ipHi <
FeII.nFeIILevel_local; ipHi++ )
2544 long int ipLo , ipHi;
2552 for( ipLo=0; ipLo < (
FeII.nFeIILevel_local - 1); ipLo++ )
2554 for( ipHi=ipLo + 1; ipHi <
FeII.nFeIILevel_local; ipHi++ )
2609 step = log10( xLamHigh/xLamLow)/ncell;
2610 wl1 = log10( xLamLow);
2624 FeII.lgFeIILargeOn =
true;
2629 FeII.nFeIILevel_local =
FeII.nFeIILevel_malloc;
2646 if(
FeII.nFeIILevel_local <16 )
2648 fprintf(
ioQQQ,
" This would be too few levels, must have at least 16.\n" );
2653 fprintf(
ioQQQ,
" This would be too many levels.\n" );
2660 else if( p.
nMatch(
"SLOW") )
2666 else if( p.
nMatch(
"REDI") )
2679 else if( p.
nMatch(
" CRD") )
2684 else if( p.
nMatch(
"CRDW") )
2690 else if( !p.
nMatch(
"SHOW") )
2692 fprintf(
ioQQQ,
" There should have been a second keyword on this command.\n");
2693 fprintf(
ioQQQ,
" Options are _PRD, _CRD, CRDW (_ is space). Sorry.\n");
2700 FeII.ipRedisFcnResonance = ipRedis;
2703 else if( p.
nMatch(
"SUBO") )
2705 FeII.ipRedisFcnSubordinate = ipRedis;
2708 else if( p.
nMatch(
"SHOW") )
2710 fprintf(
ioQQQ,
" FeII resonance lines are ");
2713 fprintf(
ioQQQ,
"complete redistribution with wings\n");
2715 else if(
FeII.ipRedisFcnResonance ==
ipCRD )
2717 fprintf(
ioQQQ,
"complete redistribution with core only.\n");
2719 else if(
FeII.ipRedisFcnResonance ==
ipPRD )
2721 fprintf(
ioQQQ,
"partial redistribution.\n");
2725 fprintf(
ioQQQ,
" PROBLEM Impossible value for ipRedisFcnResonance.\n");
2729 fprintf(
ioQQQ,
" FeII subordinate lines are ");
2732 fprintf(
ioQQQ,
"complete redistribution with wings\n");
2734 else if(
FeII.ipRedisFcnSubordinate ==
ipCRD )
2736 fprintf(
ioQQQ,
"complete redistribution with core only.\n");
2738 else if(
FeII.ipRedisFcnSubordinate ==
ipPRD )
2740 fprintf(
ioQQQ,
"partial redistribution.\n");
2744 fprintf(
ioQQQ,
" PROBLEM Impossible value for ipRedisFcnSubordinate.\n");
2750 fprintf(
ioQQQ,
" here should have been a second keyword on this command.\n");
2751 fprintf(
ioQQQ,
" Options are RESONANCE, SUBORDINATE. Sorry.\n");
2757 else if( p.
nMatch(
"TRAC") )
2759 FeII.lgPrint =
true;
2763 else if( p.
nMatch(
"SIMU") )
2766 FeII.lgSimulate =
true;
2770 else if( p.
nMatch(
"CONT") )
2779 if(
FeII.feconwlLo<=0. ||
FeII.feconwlHi<=0. ||
FeII.nfe2con<= 0 )
2781 fprintf(
ioQQQ,
" there are three numbers on the FeII continuum command, start and end wavelengths, and number of intervals.\n");
2782 fprintf(
ioQQQ,
" all three must be greater than zero, sorry.\n");
2786 if(
FeII.feconwlLo>=
FeII.feconwlHi )
2788 fprintf(
ioQQQ,
" there are three numbers on the FeII continuum command, start and end wavelengths, and number of intervals.\n");
2789 fprintf(
ioQQQ,
" the second wavelength must be greater than the first, sorry.\n");
2805 for( n=0; n<
FeII.nFeIILevel_local-1; ++n)
2807 for( ipHi=n+1; ipHi<
FeII.nFeIILevel_local; ++ipHi )
2811 fprintf(io,
"%li\t%li\t%.2e\n", n , ipHi ,
2826 for( n=0; n<
FeII.nFeIILevel_local-1; ++n)
2828 for( ipHi=n+1; ipHi<
FeII.nFeIILevel_local; ++ipHi )
2848 long ipLoMax , ipHiMax;
2850 long int ipLo, ipHi;
2864 for( ipHi=1; ipHi<
FeII.nFeIILevel_local; ++ipHi )
2866 for( ipLo=0; ipLo<ipHi; ++ipLo)
2871 (*tr.
Hi()).Pop() > 1e-30 )
2873 if( (*tr.
Hi()).Pop() > smallfloat &&
2879 if( RadPres1 > RadMax )
2896 fprintf(
ioQQQ,
"DEBUG FeIIRadPress finds P(FeII) = %.2e %.2e %li %li width %.2e\n",
2922 for(
long ipHi=1; ipHi<
FeII.nFeIILevel_local; ++ipHi )
2925 if( (*(*tr).Hi()).Pop() > 1e-30 )
2928 (*(*tr).Hi()).Pop() * (*tr).EnergyErg;
2937#if defined(__HP_aCC)
2946 double EnerLyaProf2,
2963 if(
FeII.lgLyaPumpOn )
2971 de = 1.37194e-06*
hydro.HLineWidth*2.0/3.0;
2974 EnerLyaProf2 = 82259.0 - de;
2975 EnerLyaProf3 = 82259.0 + de;
3004 if(
FeII.lgLyaPumpOn )
3006 for( ipLo=0; ipLo < (
FeII.nFeIILevel_local - 1); ipLo++ )
3008 for( ipHi=ipLo + 1; ipHi <
FeII.nFeIILevel_local; ipHi++ )
3023 Gup_ov_Glo = (*tr.
Hi()).g()/(*tr.
Lo()).g();
3028 if( EnergyWN >=
EnerLyaProf1 && EnergyWN <= EnerLyaProf4 && taux > 1 )
3045 if( EnergyWN < EnerLyaProf2 )
3050 else if( EnergyWN < EnerLyaProf3 )
3070 PumpRate = FeIILineWidthHz * PhotOccNum_at_nu * tr.
Emis().
Aul()*
3071 powi(82259.0f/EnergyWN,3);
3073 PumpRate = tr.
Emis().
Aul()* PhotOccNum_at_nu;
3079 Fe2LPump[ipLo][ipHi] += PumpRate*Gup_ov_Glo;
3089#if defined(__HP_aCC)
static double * Fe2DepCoef
void FeII_RT_TauInc(void)
void FeIIPunPop(FILE *ioPUN, bool lgPunchRange, long int ipRangeLo, long int ipRangeHi, bool lgPunchDensity)
static double * Fe2LevelPop
double FeIIRadPress(void)
STATIC void FeIILyaPump(void)
STATIC void FeIICollRatesBoltzmann(void)
void FeIIPunDepart(FILE *ioPUN, bool lgDoAll)
void FeII_RT_tau_reset(void)
STATIC bool lgFeIIEverCalled
void AssertFeIIDep(double *pred, double *BigError, double *StdDev)
static double * Fe2ColDen
void FeIIAccel(double *fe2drive)
STATIC int FeIIBandsCreate(const char chFile[])
void FeIIPun1Depart(FILE *ioPUN, long int nPUN)
void FeIIPunchLineStuff(FILE *io, realnum xLimit, long index)
void FeIIPunchColden(FILE *ioPUN)
void FeIIPunchLevels(FILE *ioPUN)
static double EnerLyaProf1
static double * FeIIBoltzmann
void FeIIPunchOpticalDepth(FILE *ioPUN)
static double PhotOccNumLyaCenter
static double ** Fe2LPump
double FeIISumBand(realnum wl1, realnum wl2, double *SumBandInward)
static double EnerLyaProf4
void ParseAtomFeII(Parser &p)
STATIC void FeIIContCreate(double xLamLow, double xLamHigh, long int ncell)
static double ** Fe2CPump
void FeIIPunData(FILE *ioPUN, bool lgDoAll)
static realnum *** sPradDat
static long int * nnPradDat
void FeIISaveLines(FILE *ioPUN)
static realnum ** Fe2Coll
void FeII_Colden(const char *chLabel)
double FeII_InterEnergy(void)
const int FILENAME_PATH_LENGTH_2
double FFmtRead(const char *chCard, long int *ipnt, long int last, bool *lgEOL)
bool fp_equal(sys_float x, sys_float y, int n=3)
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
NORETURN void TotalInsanity(void)
sys_float SDIV(sys_float x)
#define DEBUG_ENTRY(funcname)
double powi(double, long int)
long int cdLine(const char *chLabel, realnum wavelength, double *relint, double *absint)
realnum & col_str() const
double & ColOvTot() const
long int & ipFine() const
double & xIntensity() const
realnum & Pelec_esc() const
realnum & opacity() const
realnum & FracInwd() const
realnum & dampXvel() const
bool nMatch(const char *chKey) const
double getNumberCheck(const char *chDesc)
TransitionProxy::iterator iterator
void setHi(int ipHi) const
CollisionProxy Coll() const
void AddLine2Stack() const
realnum EnergyErg() const
realnum & EnergyWN() const
qList::iterator Lo() const
qList::iterator Hi() const
void setLo(int ipLo) const
EmissionList::reference Emis() const
void outline_resonance() const
long ipFineCont(double energy_ryd)
long ipoint(double energy_ryd)
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
realnum GetDopplerWidth(realnum massAMU)
t_iso_sp iso_sp[NISO][LIMELM]
double RefIndex(double EnergyWN)
double abscf(double gf, double enercm, double gl)
UNUSED const double TRANS_PROB_CONST
UNUSED const double COLL_CONST
double PressureRadiationLine(const TransitionProxy &t, realnum DopplerWidth)
void RT_line_one(const TransitionProxy &t, bool lgShield_this_zone, realnum pestrk, realnum DopplerWidth)
void RT_line_one_tau_reset(const TransitionProxy &t)
void RT_line_one_tauinc(const TransitionProxy &t, long int mas_species, long int mas_ion, long int mas_hi, long int mas_lo, realnum DopplerWidth)
double RT_LineWidth(const TransitionProxy &t, realnum DopplerWidth)
void RT_OTS_AddLine(double ots, long int ip)
void Save1Line(const TransitionProxy &t, FILE *io, realnum xLimit, long index, realnum DopplerWidth)
void Save1LineData(const TransitionProxy &t, FILE *io, bool lgCS_2, bool &lgPrint)
TransitionList Fe2LevN("Fe2LevN", &Fe2LevNStates)
multi_arr< int, 2 > ipFe2LevN
TransitionList TauLines("TauLines", &AnonStates)
vector< TransitionList > AllTransitions
void getrs_wrapper(char trans, long N, long nrhs, double *A, long lda, int32 *ipiv, double *B, long ldb, int32 *info)
void getrf_wrapper(long M, long N, double *A, long lda, int32 *ipiv, int32 *info)
void DumpLine(const TransitionProxy &t)