190 if( strcmp(chLabel,
"ZERO") == 0 )
200 else if( strcmp(chLabel,
"ADD ") == 0 )
211 else if( strcmp(chLabel,
"PRIN") != 0 )
213 fprintf(
ioQQQ,
" FeII_Colden does not understand the label %s\n",
315 Fe2SavN[ipHi]=(
double*)
MALLOC(
sizeof(
double )*(
unsigned long)ipHi);
340 for(ipHi=1; ipHi <
NPRADDAT; ipHi++)
346 for( ipLo=0; ipLo< ipHi; ipLo++ )
353 for(ipHi=0; ipHi <
NPRADDAT; ipHi++)
355 for( ipLo=0; ipLo< ipHi; ipLo++ )
372 ncs1[ipHi]=(
int*)
MALLOC(
sizeof(
int)*(
unsigned long)ipHi);
393 for( ipLo=0; ipLo < ipHi; ipLo++ )
412 fprintf(
ioQQQ,
" FeIICreate opening fe2nn.dat:");
424 for( i=0; i < 8; i++ )
428 fprintf(
ioQQQ,
" fe2nn.dat error reading data\n" );
435 ASSERT( k+19 < NPRADDAT+1 );
437 "%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld%4ld",
443 # if !defined(NDEBUG)
444 for( m1=0; m1<20; ++m1 )
455 fprintf(
ioQQQ,
" FeIICreate opening fe2energies.dat:");
458 ioDATA =
open_data(
"fe2energies.dat",
"r" );
466 while( chLine[0] ==
'#' )
470 fprintf(
ioQQQ,
" fe2energies.dat error reading data\n" );
477 sscanf( chLine,
"%lf", &help );
493 fprintf(
ioQQQ,
" FeIICreate opening fe2rad.dat:");
501 fprintf(
ioQQQ,
" fe2rad.dat error reading data\n" );
505 sscanf( chLine ,
"%ld%ld%ld",&lo, &ihi, &m1);
506 if( lo!=3 || ihi!=1 || m1!=28 )
508 fprintf(
ioQQQ,
" fe2rad.dat has the wrong magic numbers, expected 3 1 28\n" );
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 ,
535 &save[0], &save[1] , &m3);
542 # define USE_OLD true
547 if(
Fe2LevN[ihi-1][lo-1].Emis->Aul == 1e3f )
565 ncs1[ihi-1][lo-1] = (int)m3;
586 fprintf(
ioQQQ,
" FeIICreate opening fe2col.dat:");
592 for( ipHi=1; ipHi<
NPRADDAT; ++ipHi )
594 for( ipLo=0; ipLo<ipHi; ++ipLo )
600 fprintf(
ioQQQ,
" fe2col.dat error reading data\n" );
604 "%ld%ld%lf%lf%lf%lf%lf%lf%lf%lf",
606 &save[0], &save[1] , &save[2],&save[3], &save[4] , &save[5],
682 fprintf(
ioQQQ,
" failed to open FeII bands file \n");
736 fprintf(
ioQQQ,
" FeIILevelPops fe2 pops called\n");
795 for( ipLo=0; ipLo<n; ++ipLo )
806 for( ipLo=0; ipLo<n; ++ipLo )
823 enum {DEBUG_LOC=
false};
854 # define AMAT(I_,J_) (*(amat+(I_)*FeII.nFeIILevel+(J_)))
873 fprintf(
ioQQQ,
"DISASTER FeIILevelPops: dgetrs finds singular or ill-conditioned matrix\n" );
887 fprintf(
ioQQQ,
"PROBLEM FeIILevelPops finds non-positive level population, level is %ld pop is %g\n",
897 fprintf(
ioQQQ ,
"PROBLEM FeIILevelPops exits with negative level populations.\n");
957 Fe2LevN[ipHi][ipLo].Lo->Pop*Fe2LPump[ipLo][ipHi] -
977 enum {DEBUG_LOC=
false};
981 fprintf(
ioQQQ,
" sum all %.1e dest rate%.1e escR= %.1e\n",
1071 static double OldTemp = -1.;
1081 realnum FracLowTe , FracHighTe;
1082 static realnum tt[8]={1e3f,3e3f,5e3f,7e3f,1e4f,12e3f,15e3f,2e4f};
1087 static long int nFeII_old = -1;
1109 if(
ncs1[ipHi][ipLo] == 3 )
1117 else if(
ncs1[ipHi][ipLo] == 1 )
1124 else if(
ncs1[ipHi][ipLo] == 2 )
1136 fprintf(
ioQQQ,
">>>PROBLEM impossible ncs1 in FeIILevelPops\n");
1142 gb = (
realnum)(ag + (-cg*
POW2(y) + dg)*(log(1.0+1.0/y) - 0.4/
1143 POW2(y + 1.0)) + cg*y);
1192 for( i=0; i < 8-1; i++ )
1204 fprintf(
ioQQQ,
" Insanity while looking for temperature in coll str array, te=%g.\n",
1212 FracHighTe = ((
realnum)
phycon.
te - tt[ipTemp])/(tt[ipTempp1] - tt[ipTemp]);
1217 FracLowTe = 1.f - FracHighTe;
1219 for( ipHi=1; ipHi <
NPRADDAT; ipHi++ )
1221 for( ipLo=0; ipLo < ipHi; ipLo++ )
1228 ASSERT( ipHiFe2-1 >= 0 );
1230 ASSERT( ipLoFe2-1 >= 0 );
1238 sPradDat[ipHi][ipLo][ipTemp] * FracLowTe +
1239 sPradDat[ipHi][ipLo][ipTempp1] * FracHighTe;
1290 enum {DEBUG_LOC=
false};
1293 for( ipLo=0; ipLo < 52; ipLo++ )
1295 fprintf(
ioQQQ,
"%e %e\n",
1355 double SumBandFe2_v;
1388 if(
Fe2LevN[ipHi][ipLo].EnergyErg >= alo &&
1389 Fe2LevN[ipHi][ipLo].EnergyErg <= ahi )
1400 for( ipLo=0; ipLo < ipHi; ++ipLo )
1402 if(
Fe2LevN[ipHi][ipLo].WLAng >= wl1 &&
1403 Fe2LevN[ipHi][ipLo].WLAng < wl2 )
1408 return( SumBandFe2_v );
1432 if(
Fe2LevN[ipHi][ipLo].ipCont > 0 )
1434 -8 , -8 , ipHi , ipLo );
1463 for( ipLo=0; ipLo < ipHi; ipLo++ )
1503 if( fabs(
Fe2LevN[ipHi][ipLo].Emis->Aul- 1e-5 ) > 1e-8 )
1584 fprintf(
ioQQQ,
" FeII_RT_Make called\n");
1595 if(
Fe2LevN[ipHi][ipLo].ipCont > 0 )
1653 enum {DEBUG_LOC=
false};
1673 fprintf(ioPUN ,
"%.2f\t%li\n", 0., (
long)
Fe2LevN[1][0].Lo->g );
1676 fprintf(ioPUN ,
"%.2f\t%li\n",
1677 Fe2LevN[ipHi][0].EnergyWN ,
1678 (
long)Fe2LevN[ipHi][0].Hi->g);
1694 fprintf(ioPUN ,
"%.2f\t%li\t%.3e\n", 0., (
long)
Fe2LevN[1][0].Lo->g ,
Fe2ColDen[0]);
1697 fprintf(ioPUN ,
"%.2f\t%li\t%.3e\n",
1698 Fe2LevN[n][0].EnergyWN ,
1699 (
long)Fe2LevN[n][0].Hi->g,
1731 fprintf( ioPUN,
"%ld\t%ld\t%.2f\t%.2e\n",
1735 Fe2LevN[ipHi][ipLo].Emis->TauIn );
1753 double hbeta, absint , renorm;
1773 fprintf( ioPUN,
" up low log I, I/I(LineSave), Tau\n" );
1783 if(
Fe2LevN[ipHi][ipLo].Emis->TauIn < TauMase )
1794 fprintf( ioPUN,
" Most negative optical depth was %4ld%4ld%10.2e\n",
1795 MaseHi, MaseLow, TauMase );
1799 if(
cdLine(
"TOTL", 4861 , &hbeta , &absint)<=0 )
1801 fprintf(
ioQQQ,
" FeIILevelPops could not find Hbeta with cdLine.\n" );
1805 fprintf( ioPUN,
"Hbeta=%7.3f %10.2e\n",
1830 if( (
Fe2SavN[ipHi][ipLo] > thresh &&
1838 fprintf( ioPUN,
"%ld\t%ld\t%.2f\t%.3f\n",
1847 fprintf( ioPUN,
"%ld\t%ld\t%.2f\t%.3f\t%.2e\t%.2e\n",
1908 for( ipLo=0; ipLo < ipHi; ipLo++ )
2100 long int ipLo , ipHi;
2107 " FeIIPunData ALL option not implemented yet 1\n" );
2116 for( ipHi=1; ipHi<limit; ++ipHi )
2118 for( ipLo=0; ipLo<ipHi; ++ipLo )
2120 fprintf(ioPUN,
"%4li%4li ",ipLo,ipHi );
2124 fprintf( ioPUN ,
"\n");
2135 for( ipLo=0; ipLo<ipHi; ++ipLo )
2142 if(
ncs1[ipHi][ipLo] == 3 &&
2143 (fabs(
Fe2LevN[ipHi][ipLo].Emis->Aul-1e-5) < 1e-8 ) )
2151 fprintf(ioPUN,
"%4li%4li ",ipLo+1,ipHi+1 );
2156 fprintf( ioPUN ,
" %li lines skiped\n",nSkip);
2175 const int LevDep[
NLEVDEP]={1,10,25,45,64,124,206,249,295,347,371};
2177 static bool lgFIRST=
true;
2182 if( lgFIRST && !lgDoAll )
2187 fprintf( ioPUN ,
"%i\t", LevDep[i] );
2189 fprintf( ioPUN ,
"\n");
2199 fprintf( ioPUN ,
"\n");
2208 fprintf( ioPUN ,
"\t");
2210 fprintf( ioPUN ,
"\n");
2230 assert( ioPUN != NULL );
2240 fprintf( ioPUN ,
"%e ",0. );
2339 long int ipRangeLo ,
2341 long int ipRangeHi ,
2343 bool lgPunchDensity )
2348 const int LevPop[
NLEVPOP]={1,10,25,45,64,124,206,249,295,347,371};
2350 static bool lgFIRST=
true;
2356 if( !lgPunchDensity )
2361 if( lgFIRST && !lgPunchRange )
2367 fprintf( ioPUN ,
"%i\t", LevPop[i] );
2369 fprintf( ioPUN ,
"\n");
2380 for( i=ipRangeLo; i<ipRangeHi; ++i )
2383 fprintf( ioPUN ,
"%.3e\t",
Fe2LevelPop[i]/denominator );
2385 fprintf( ioPUN ,
"\n");
2393 fprintf( ioPUN ,
"%.3e\t",
Fe2LevelPop[LevPop[i]-1]/denominator );
2395 fprintf( ioPUN ,
"\n");
2415 assert( ioPUN != NULL );
2426 fprintf( ioPUN ,
"%e ",0. );
2440 const char chFile[] )
2444 const char* chFile1;
2465 chFile1 = ( strlen(chFile) == 0 ) ?
"bands_Fe2.dat" : chFile;
2470 fprintf(
ioQQQ,
" FeIICreate opening %s:", chFile1 );
2479 if(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioDATA ) == NULL )
2481 fprintf(
ioQQQ,
" FeIICreate could not read first line of %s.\n", chFile1 );
2484 while(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioDATA ) != NULL )
2488 if( chLine[0] !=
'#')
2493 if( fseek( ioDATA , 0 , SEEK_SET ) != 0 )
2495 fprintf(
ioQQQ,
" FeIICreate could not rewind %s.\n", chFile1 );
2508 if(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioDATA ) == NULL )
2510 fprintf(
ioQQQ,
" FeIICreate could not read first line of %s.\n", chFile1 );
2519 " FeIICreate: the version of %s is not the current version.\n", chFile1 );
2525 while(
read_whole_line( chLine , (
int)
sizeof(chLine) , ioDATA ) != NULL )
2529 if( chLine[0] !=
'#')
2535 fprintf(
ioQQQ,
" There should have been a number on this band line 1. Sorry.\n" );
2536 fprintf(
ioQQQ,
"string==%s==\n" ,chLine );
2542 fprintf(
ioQQQ,
" There should have been a number on this band line 2. Sorry.\n" );
2543 fprintf(
ioQQQ,
"string==%s==\n" ,chLine );
2549 fprintf(
ioQQQ,
" There should have been a number on this band line 3. Sorry.\n" );
2550 fprintf(
ioQQQ,
"string==%s==\n" ,chLine );
2564 fprintf(
ioQQQ,
" FeII band %li has none positive entry.\n",i );
2570 fprintf(
ioQQQ,
" FeII band %li has improper bounds.\n" ,i);
2618 *BigError =
MAX2( *BigError , error );
2651 if(
Fe2LevN[ipHi][ipLo].ipCont < 1)
2678 long int ipLo , ipHi;
2691 if(
Fe2LevN[ipHi][ipLo].ipCont < 1)
2744 step = log10( xLamHigh/xLamLow)/ncell;
2745 wl1 = log10( xLamLow);
2747 for( i=1; i<ncell; ++i)
2753 for( i=0; i<(ncell-1); ++i)
2760 for( i=0; i<ncell; ++i)
2766 enum {DEBUG_LOC=
false};
2770 ioDEBUG = fopen(
"fe2cont.txt",
"w" );
2773 fprintf(
ioQQQ,
" fe2con could not open fe2cont.txt for writing\n");
2776 for( i=0; i<ncell; ++i)
2779 fprintf(ioDEBUG,
"%.1f\t%.1f\t%.1f\n",
2822 if(
nMatch(
"LEVE",chCard) )
2837 fprintf(
ioQQQ,
" This would be too few levels, must have at least 16.\n" );
2842 fprintf(
ioQQQ,
" This would be too many levels.\n" );
2849 else if(
nMatch(
"SLOW",chCard) )
2855 else if(
nMatch(
"REDI",chCard) )
2863 if(
nMatch(
" PRD",chCard) )
2868 else if(
nMatch(
" CRD",chCard) )
2873 else if(
nMatch(
"CRDW",chCard) )
2879 else if( !
nMatch(
"SHOW",chCard) )
2881 fprintf(
ioQQQ,
" There should have been a second keyword on this command.\n");
2882 fprintf(
ioQQQ,
" Options are _PRD, _CRD, CRDW (_ is space). Sorry.\n");
2887 if(
nMatch(
"RESO",chCard) )
2892 else if(
nMatch(
"SUBO",chCard) )
2897 else if(
nMatch(
"SHOW",chCard) )
2899 fprintf(
ioQQQ,
" FeII resonance lines are ");
2902 fprintf(
ioQQQ,
"complete redistribution with wings\n");
2906 fprintf(
ioQQQ,
"complete redistribution with core only.\n");
2910 fprintf(
ioQQQ,
"partial redistribution.\n");
2914 fprintf(
ioQQQ,
" PROBLEM Impossible value for ipRedisFcnResonance.\n");
2918 fprintf(
ioQQQ,
" FeII subordinate lines are ");
2921 fprintf(
ioQQQ,
"complete redistribution with wings\n");
2925 fprintf(
ioQQQ,
"complete redistribution with core only.\n");
2929 fprintf(
ioQQQ,
"partial redistribution.\n");
2933 fprintf(
ioQQQ,
" PROBLEM Impossible value for ipRedisFcnSubordinate.\n");
2939 fprintf(
ioQQQ,
" here should have been a second keyword on this command.\n");
2940 fprintf(
ioQQQ,
" Options are RESONANCE, SUBORDINATE. Sorry.\n");
2946 else if(
nMatch(
"TRAC",chCard) )
2952 else if(
nMatch(
"SIMU",chCard) )
2959 else if(
nMatch(
"CONT",chCard) )
2970 fprintf(
ioQQQ,
" there are three numbers on the FeII continuum command, start and end wavelengths, and number of intervals.\n");
2977 fprintf(
ioQQQ,
" there are three numbers on the FeII continuum command, start and end wavelengths, and number of intervals.\n");
2978 fprintf(
ioQQQ,
" all three must be greater than zero, sorry.\n");
2984 fprintf(
ioQQQ,
" there are three numbers on the FeII continuum command, start and end wavelengths, and number of intervals.\n");
2985 fprintf(
ioQQQ,
" the second wavelength must be greater than the first, sorry.\n");
3005 if(
Fe2LevN[ipHi][n].ipCont > 0 )
3006 fprintf(io,
"%li\t%li\t%.2e\n", n , ipHi ,
3007 Fe2LevN[ipHi][n].Emis->TauIn );
3043 long ipLoMax , ipHiMax;
3045 long int ipLo, ipHi;
3058 for( ipLo=0; ipLo<ipHi; ++ipLo)
3061 if(
Fe2LevN[ipHi][ipLo].ipCont > 0 &&
3062 Fe2LevN[ipHi][ipLo].Hi->Pop > 1e-30 )
3064 if(
Fe2LevN[ipHi][ipLo].Hi->Pop > smallfloat &&
3070 if( RadPres1 > RadMax )
3087 fprintf(
ioQQQ,
"DEBUG FeIIRadPress finds P(FeII) = %.2e %.2e %li %li width %.2e\n",
3120 if(
Fe2LevN[ipHi][n].Hi->Pop > 1e-30 )
3133 #if defined(__HP_aCC)
3134 # pragma OPT_LEVEL 1
3142 double EnerLyaProf2,
3170 EnerLyaProf2 = 82259.0 - de;
3171 EnerLyaProf3 = 82259.0 + de;
3223 if( EnergyWN >=
EnerLyaProf1 && EnergyWN <= EnerLyaProf4 && taux > 1 )
3240 if( EnergyWN < EnerLyaProf2 )
3245 else if( EnergyWN < EnerLyaProf3 )
3265 PumpRate = FeIILineWidthHz * PhotOccNum_at_nu *
Fe2LevN[ipHi][ipLo].
Emis->
Aul*
3266 powi(82259.0f/EnergyWN,3);
3274 Fe2LPump[ipLo][ipHi] += PumpRate*Gup_ov_Glo;
3284 #if defined(__HP_aCC)
3285 #pragma OPTIMIZE OFF