/*mihals generate continuum from Mihalas stellar atmosphere */
#include "cddefines.h"
#ifdef NDATA
#	undef NDATA
#endif
#define	NDATA	47
#include "bounds.h"
#include "incidspec.h"
#include "mihalas.h"

void mihals(long int *nstar, 
  double par1, 
  double par2)
{
	long int i, 
	  mtemp;
	double alogg, 
	  fac1, 
	  fac2, 
	  temp;

	/* wavelength scale in inverse microns */
	static double wl[NDATA]={87.7939,75.0851,62.3770,62.3760,53.1367,
	  43.8975,43.8965,38.6940,38.6930,31.7719,24.8510,24.8500,22.3465,
	  19.8430,19.8420,17.3087,14.7749,13.5080,12.2411,11.6077,10.9750,
	  10.9740,10.2883,9.7549,8.2307,6.5541,4.8780,4.8770,3.8105,2.7440,
	  2.7430,2.4006,2.0577,1.7564,1.7554,1.5242,1.2199,1.2189,0.9526,
	  0.6864,0.6854,0.5335,0.4395,0.4385,0.2212,0.0033,1e-4};

	/* these puppy's are in MAGNITUDES!!!!! */
	static double t30g4[NDATA]={51.490,43.717,35.900,35.898,30.234,
	  24.636,22.972,19.877,20.031,15.728,11.647,11.199,9.782,8.417,
	  4.205,3.631,2.850,2.413,1.965,1.741,1.520,-1.931,-1.923,-1.909,
	  -1.828,-1.645,-1.317,-1.323,-0.993,-0.510,-0.748,-0.517,-0.243,
	  0.048,0.048,0.312,0.740,0.695,1.186,1.852,1.840,2.359,2.764,
	  2.761,4.202,13.107,20.699};

	static double t325g4[NDATA]={43.608,36.844,30.128,30.114,25.301,
	  20.573,17.983,15.386,13.113,10.637,7.128,6.695,5.527,4.367,1.771,
	  1.418,0.933,0.646,0.341,0.185,0.029,-2.016,-2.011,-1.999,-1.926,
	  -1.750,-1.421,-1.437,-1.088,-0.569,-0.764,-0.529,-0.248,0.049,
	  0.049,0.319,0.754,0.718,1.214,1.884,1.875,2.395,2.801,2.798,
	  4.238,13.143,20.735};

	static double t35g4[NDATA]={37.814,31.886,26.019,25.986,21.790,
	  17.670,14.982,12.661,4.193,4.020,2.826,2.747,2.082,1.283,0.696,
	  0.388,0.017,-0.192,-0.410,-0.521,-0.631,-2.056,-2.049,-2.038,
	  -1.964,-1.790,-1.462,-1.488,-1.138,-0.605,-0.761,-0.527,-0.248,
	  0.049,0.048,0.319,0.756,0.728,1.225,1.895,1.886,2.408,2.814,
	  2.812,4.249,13.184,20.776};

	static double t375g4[NDATA]={33.776,28.424,23.152,23.069,19.318,
	  15.642,13.211,11.053,2.397,1.888,1.150,1.143,0.731,0.221,0.097,
	  -0.171,-0.475,-0.637,-0.803,-0.886,-0.966,-2.102,-2.093,-2.079,
	  -2.000,-1.820,-1.487,-1.518,-1.166,-0.628,-0.756,-0.524,-0.246,
	  0.048,0.047,0.317,0.753,0.730,1.225,1.893,1.888,2.407,2.812,
	  2.811,4.249,13.207,20.799};

	static double t40g4[NDATA]={31.030,26.058,21.192,21.030,17.596,
	  14.242,11.981,9.539,1.543,0.939,0.289,0.289,0.002,-0.337,-0.355,
	  -0.587,-0.838,-0.969,-1.100,-1.164,-1.225,-2.158,-2.146,-2.130,
	  -2.040,-1.856,-1.514,-1.547,-1.191,-0.647,-0.758,-0.525,-0.247,
	  0.049,0.046,0.316,0.752,0.734,1.228,1.896,1.892,2.411,2.817,
	  2.816,4.256,13.215,20.807};

	static double t45g4[NDATA]={26.168,21.629,17.233,16.584,13.538,
	  10.627,3.420,2.306,0.525,-0.082,-0.662,-0.662,-0.864,-1.063,
	  -1.060,-1.242,-1.418,-1.503,-1.583,-1.620,-1.655,-2.257,-2.239,
	  -2.219,-2.123,-1.925,-1.570,-1.595,-1.231,-0.675,-0.765,-0.530,
	  -0.249,0.049,0.046,0.319,0.759,0.743,1.240,1.913,1.910,2.431,
	  2.839,2.838,4.281,13.235,20.827};

	static double t50g4[NDATA]={23.826,19.658,15.626,14.645,11.833,
	  9.190,0.621,-0.005,-0.199,-0.733,-1.241,-1.241,-1.412,-1.572,
	  -1.571,-1.717,-1.843,-1.897,-1.942,-1.961,-1.976,-2.330,-2.307,
	  -2.283,-2.178,-1.972,-1.610,-1.622,-1.254,-0.690,-0.768,-0.532,
	  -0.250,0.049,0.047,0.321,0.764,0.750,1.249,1.925,1.922,2.445,
	  2.853,2.853,4.296,13.256,20.848};

	static double t55g4[NDATA]={21.808,17.921,14.170,12.897,10.251,
	  7.788,-0.428,-0.817,-0.837,-1.281,-1.706,-1.706,-1.848,-1.979,
	  -1.979,-2.094,-2.184,-2.216,-2.237,-2.242,-2.242,-2.369,-2.342,
	  -2.315,-2.204,-1.994,-1.631,-1.530,-1.262,-0.695,-0.765,-0.531,
	  -0.250,0.049,0.048,0.322,0.766,0.752,1.252,1.929,1.927,2.451,
	  2.860,2.859,4.304,13.274,20.866};

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

	/* following atmospheres NLTE from Mihalas, NCAR-TN/STR-76 */

	/* 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 mihals]" );
		exit(1);
	}

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

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

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

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

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

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

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

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

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

	else if( mtemp < 300 || mtemp > 550 )
	{
		fprintf( ioQQQ, " This temp is not inside table.\n" );
		puts( "[Stop in mihals]" );
		exit(1);
	}

	else
	{
		/*  must interpolate on grid */
		if( mtemp > 300 && mtemp <= 325 )
		{
			fac1 = (temp - 30000.)/2500.;
			fac2 = 1. - fac1;
			for( i=0; i < NDATA; i++ )
			{
				IncidSpec.tslop[i][IncidSpec.nspec-1] = (float)(fac2*t30g4[NDATA-i-1] + 
				  fac1*t325g4[NDATA-i-1]);
			}
		}

		else if( mtemp > 325 && mtemp <= 350 )
		{
			fac1 = (temp - 32500.)/2500.;
			fac2 = 1. - fac1;
			for( i=0; i < NDATA; i++ )
			{
				IncidSpec.tslop[i][IncidSpec.nspec-1] = (float)(fac2*t325g4[NDATA-i-1] + 
				  fac1*t35g4[NDATA-i-1]);
			}
		}

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

		else if( mtemp > 375 && mtemp <= 400 )
		{
			fac1 = (temp - 37500.)/2500.;
			fac2 = 1. - fac1;
			for( i=0; i < NDATA; i++ )
			{
				IncidSpec.tslop[i][IncidSpec.nspec-1] = (float)(fac2*t375g4[NDATA-i-1] + 
				  fac1*t40g4[NDATA-i-1]);
			}
		}

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

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

		else
		{
			/* actually- else if( mtemp.gt.500.and.mtemp.le.550 ) THEN */
			fac1 = (temp - 50000.)/5000.;
			fac2 = 1. - fac1;
			for( i=0; i < NDATA; i++ )
			{
				IncidSpec.tslop[i][IncidSpec.nspec-1] = (float)(fac2*t50g4[NDATA-i-1] + 
				  fac1*t55g4[NDATA-i-1]);
			}
		}
	}

	/* convert from magnitudes to LOG10( f-nu ) */
	for( i=0; i < NDATA; i++ )
	{
		IncidSpec.tslop[i][IncidSpec.nspec-1] = 
			(float)(-IncidSpec.tslop[i][IncidSpec.nspec-1]/2.5);
	}
	*nstar = NDATA;

	/* 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 ;

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