/*kurucz79 obtain interpolated Kurucz stellar atmosphere */
#include "cddefines.h"
#define	NDATA	56
#include "incidspec.h"
#include "bounds.h"
#include "kurucz79.h"

void kurucz79(long int *nstar, 
  double par1, 
  double par2)
{
	long int i, 
	  mtemp;

	double alogg, 
	  fac1, 
	  fac2, 
	  sc30, 
	  sc35, 
	  sc40, 
	  sc45, 
	  sc50, 
	  temp;
	static int lgFirst = TRUE;
	/* wavelengths in nanometers */
	static double wl[NDATA]={22.9,23.4,24.3,25.2,25.6,26.0,26.6,27.5,
	  28.3,29.3,30.1,30.9,32.5,34.2,35.1,35.6,36.5,37.5,38.4,40.0,
	  41.5,42.5,44.5,47.2,49.4,50.6,51.6,54.0,56.4,58.7,61.2,64.0,
	  67.1,70.5,73.6,77.0,81.0,85.0,89.0,92.0,94.0,118.7,143.3,194.0,
	  256.2,331.2,406.2,471.2,592.5,775.0,1035.,1335.,5000.,6500.,
	  1.e4,9e6};

	/* following atmospheres LTE from Kurucz 1979, Ap.J. Sup 40, 1. and
	 * Kurucz (1989) private communication, newer opacities,
	 * log f_lam */
	static double t30g4[NDATA]={2.67,3.75,4.16,4.53,4.32,4.81,4.92,
	  5.19,5.42,5.65,5.53,5.83,6.30,6.63,6.77,6.82,6.97,7.09,7.24,
	  7.44,7.61,7.70,7.91,8.13,8.31,9.55,9.74,9.74,9.81,9.83,9.90,
	  9.94,9.98,10.02,10.06,10.12,10.14,10.13,10.19,10.72,11.24,11.16,
	  10.96,10.59,10.25,9.87,9.63,9.40,9.03,8.59,8.12,7.69,5.43,4.97,
	  4.22,-7.60};
	static double t35g45[NDATA]={4.96,5.32,5.55,5.86,6.00,6.19,6.53,
	  6.78,7.07,7.23,7.59,7.84,8.32,8.60,8.8,9.00,9.30,9.30,9.40,9.60,
	  9.70,9.73,9.83,9.92,10.00,10.34,10.92,10.78,10.97,10.97,10.92,
	  10.99,11.00,10.88,11.02,10.99,10.98,10.95,10.97,11.30,11.57,
	  11.36,11.19,10.81,10.43,10.04,9.76,9.53,9.16,8.72,8.24,7.80,
	  5.52,5.07,4.32,-7.50};
	static double t40g45[NDATA]={6.93,7.39,7.60,7.90,8.02,8.13,8.85,
	  9.00,9.23,9.35,9.53,10.06,10.35,10.49,10.55,10.80,11.10,11.08,
	  10.95,11.33,11.32,11.28,11.33,11.27,11.28,11.28,11.53,11.38,
	  11.51,11.53,11.45,11.54,11.52,11.41,11.52,11.43,11.45,11.42,
	  11.40,11.43,11.68,11.51,11.31,10.93,10.54,10.15,9.85,9.62,9.24,
	  8.79,8.31,7.88,5.60,5.14,4.39,-7.43};
	static double t45g45[NDATA]={8.32,9.34,9.40,10.14,10.10,10.51,11.46,
	  11.40,11.40,11.52,10.68,11.55,11.89,11.95,11.85,11.84,11.84,
	  11.77,11.48,11.72,11.74,11.73,11.80,11.77,11.80,11.79,11.83,
	  11.80,11.81,11.80,11.79,11.77,11.76,11.73,11.73,11.68,11.67,
	  11.63,11.61,11.64,11.79,11.62,11.37,10.99,10.60,10.20,9.90,9.66,
	  9.28,8.84,8.36,7.92,5.64,5.18,4.43,-7.39};
	static double t50g45[NDATA]={9.36,11.12,11.11,11.65,11.54,11.68,
	  12.12,12.07,12.09,12.17,11.65,12.04,12.22,12.27,12.16,12.13,
	  12.13,12.08,11.77,12.00,12.03,12.03,12.10,12.07,12.11,12.09,
	  12.10,12.06,12.06,12.04,12.03,11.98,11.98,11.93,11.93,11.86,
	  11.85,11.81,11.78,11.78,11.88,11.70,11.45,11.06,10.66,10.26,
	  9.95,9.71,9.33,8.88,8.40,7.97,5.68,5.22,4.47,-7.35};

#	ifdef DEBUG_FUN
	fputs( "<+>kurucz79()\n", debug_fp );
#	endif

	if( lgFirst )
	{
		sc30 = t30g4[38];
		sc35 = t35g45[38];
		sc40 = t40g45[38];
		sc45 = t45g45[38];
		sc50 = t50g45[38];

		for( i=0; i < NDATA; i++ )
		{
			t30g4[i] -= sc30;
			t35g45[i] -= sc35;
			t40g45[i] -= sc40;
			t45g45[i] -= sc45;
			t50g45[i] -= sc50;
		}
		lgFirst = FALSE;
	}

	/* PAR1 and 2 are LOG(g) and temp */
	if( par1 > 10. )
	{
		temp = par1;
		mtemp = (long)((par1 + 10.)/100.);
		alogg = par2;
	}
	else
	{
		temp = par2;
		mtemp = (long)((par2 + 10.)/100.);
		alogg = par1;
	}

	if( alogg != 4. )
	{
		fprintf( ioQQQ, " only LOG(g)=4 in table at present.\n" );
		puts( "[Stop in kurucz79]" );
		exit(1);
	}

	for( i=0; i < NDATA; i++ )
	{
		/* convert nanometers to Rydbergs */
		IncidSpec.tnu[NDATA-i-1][IncidSpec.nspec-1] = (float)(91.16/wl[i]);
	}

	if( mtemp == 300 )
	{
		for( i=0; i < NDATA; i++ )
		{
			IncidSpec.tslop[NDATA-i-1][IncidSpec.nspec-1] = (float)t30g4[i];
		}
	}

	else if( mtemp == 350 )
	{
		for( i=0; i < NDATA; i++ )
		{
			IncidSpec.tslop[NDATA-i-1][IncidSpec.nspec-1] = (float)t35g45[i];
		}
	}

	else if( mtemp == 400 )
	{
		for( i=0; i < NDATA; i++ )
		{
			IncidSpec.tslop[NDATA-i-1][IncidSpec.nspec-1] = (float)t40g45[i];
		}
	}

	else if( mtemp == 450 )
	{
		for( i=0; i < NDATA; i++ )
		{
			IncidSpec.tslop[NDATA-i-1][IncidSpec.nspec-1] = (float)t45g45[i];
		}
	}

	else if( mtemp == 500 )
	{
		for( i=0; i < NDATA; i++ )
		{
			IncidSpec.tslop[NDATA-i-1][IncidSpec.nspec-1] = (float)t50g45[i];
		}
	}

	else if( mtemp < 300 || mtemp > 500 )
	{
		fprintf( ioQQQ, " This temp is not inside Kurucz 79 table.\n" );
		puts( "[Stop in kurucz79]" );
		exit(1);
	}
	else
	{
		/* must interpolate on grid */
		if( mtemp > 300 && mtemp <= 350 )
		{
			fac1 = (temp - 30000.)/5000.;
			fac2 = 1. - fac1;

			for( i=0; i < NDATA; i++ )
			{
				IncidSpec.tslop[NDATA-i-1][IncidSpec.nspec-1] = (float)(fac2*
				  t30g4[i] + fac1*t35g45[i]);
			}
		}

		else if( mtemp > 350 && mtemp <= 400 )
		{
			fac1 = (temp - 35000.)/5000.;
			fac2 = 1. - fac1;
			for( i=0; i < NDATA; i++ )
			{
				IncidSpec.tslop[NDATA-i-1][IncidSpec.nspec-1] = (float)(fac2*
				  t35g45[i] + fac1*t40g45[i]);
			}
		}

		else if( mtemp > 400 && mtemp <= 450 )
		{
			fac1 = (temp - 40000.)/5000.;
			fac2 = 1. - fac1;
			for( i=0; i < NDATA; i++ )
			{
				IncidSpec.tslop[NDATA-i-1][IncidSpec.nspec-1] = (float)(fac2*
				  t40g45[i] + fac1*t45g45[i]);
			}
		}

		else if( mtemp > 450 && mtemp <= 500 )
		{
			fac1 = (temp - 45000.)/5000.;
			fac2 = 1. - fac1;
			for( i=0; i < NDATA; i++ )
			{
				IncidSpec.tslop[NDATA-i-1][IncidSpec.nspec-1] = (float)(fac2*
				  t45g45[i] + fac1*t50g45[i]);
			}
		}
	}

	/* now convert to f-nu (above is f-lam) */
	for( i=0; i < NDATA; i++ )
	{
		IncidSpec.tslop[i][IncidSpec.nspec-1] += 
			(float)(-2.*log10(IncidSpec.tnu[i][IncidSpec.nspec-1]));
	}

	/* we now have the final product loaded in tnu and tslop, but need to
	 * lower the lower limit to the current low energy limit of the code,
	 * assume all stars are on Rayleigh-Jeans tail of planck function,
	 * so fnu propto nu^2 */

	/* this is the log of the low-energy f_nu */
	IncidSpec.tslop[0][IncidSpec.nspec-1] += 
		(float)(2.*log10( bounds.emm / IncidSpec.tnu[0][IncidSpec.nspec-1] ) );

	/* now reset lowest energy to current bounds of code */
	IncidSpec.tnu[0][IncidSpec.nspec-1] = bounds.emm ;

	*nstar = NDATA;

#	ifdef DEBUG_FUN
	fputs( " <->kurucz79()\n", debug_fp );
#	endif
	return;
}