/*DrvCaseBHS allow user to query hydrogen A's, asks for up, low level, gives A, drive hyas */
#include "cddefines.h"
#include "ffmtread.h"
#include "drvcasebhs.h"
#include "createdata.h"
#include "hsrate.h"  

void DrvCaseBHS(void)
{
	char chCard[133];
	int lgEOL,
		lgHit;
	long int i, 
	  n1, 
	  ipZ ,
	  n2;
	double Temperature,
		Density;

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

	/* this routine is entered with the command DRIVE CASEB, and
	 * interpolates on the Hummer & Storey case b data set */

	/* read in some external data files, but only if this is first call */
	fprintf(ioQQQ," I will get needed H data files. This will take a second.\n");
	CreateData();

	/* first get the charge, only H and He at present */
	lgHit = FALSE;
	ipZ = 0;
	while( !lgHit )
	{
		fprintf( ioQQQ, " Enter atomic number of species, either 1(H) or 2(He).\n" );
		if( fgets( chCard , sizeof(chCard) , stdin ) == NULL )
		{
			fprintf( ioQQQ, " error getting species \n" );
			puts( "[Stop in DrvCaseBHS]" );
			exit(1);
		}

		i = 1;
		ipZ = (long int)FFmtRead(chCard,&i,76,&lgEOL);
		if( lgEOL || ipZ< 1 || ipZ > 2 )
		{
			fprintf( ioQQQ, " must be either 1 or 2!\n" );
		}
		else
		{
			lgHit = TRUE;
		}
	}

	fprintf(ioQQQ," In the following temperatures <10 are log, >=10 linear.\n");
	fprintf(ioQQQ," The density is always a log.\n");
	fprintf(ioQQQ," The order of the quantum numbers do not matter.\n");
	fprintf(ioQQQ," The smallest must not be smaller than 2,\n");
	fprintf(ioQQQ," and the largest must not be larger than 25.\n");
	fprintf(ioQQQ," Units of emissivity are erg cm^3 s^-1\n\n");
	fprintf(ioQQQ," The limits of the HS tables are 2 <= n <= 25.\n");

	lgHit = TRUE;
	/* this is always true */
	while( lgHit )
	{
		fprintf( ioQQQ, " Enter 4 numbers, temperature, density, 2 quantum numbers, null line stop.\n" );
		if( fgets( chCard , sizeof(chCard) , stdin ) == NULL )
		{
			fprintf( ioQQQ, " Thanks for interpolating on the Hummer & Storey data set!\n" );
			puts( "[Stop in DrvCaseBHS]" );
			exit(1);
		}

		i = 1;
		Temperature = FFmtRead(chCard,&i,76,&lgEOL);
		if( lgEOL )
		{
			fprintf( ioQQQ, " error getting temperature!\n" );
			break;
		}

		/* log if less than 10 */
		if( Temperature < 10. )
		{
			Temperature = pow(10., Temperature );
		}
		fprintf(ioQQQ," Temperature is %g\n", Temperature );

		Density = FFmtRead(chCard,&i,76,&lgEOL);
		if( lgEOL )
		{
			fprintf( ioQQQ, " error getting density!\n" );
			break;
		}
		Density = pow(10., Density );
		fprintf(ioQQQ," Density is %g\n", Density );

		/* these quantum numbers can be in any order */
		n1 = (long)FFmtRead(chCard,&i,76,&lgEOL);
		if( lgEOL )
		{
			fprintf( ioQQQ, " error getting quantum number!\n" );
			break;
		}

		n2 = (long)FFmtRead(chCard,&i,76,&lgEOL);
		if( lgEOL )
		{
			fprintf( ioQQQ, " error getting quantum number!\n" );
			break;
		}

		if( MAX2( n1 , n2 ) > 25 )
		{
			fprintf( ioQQQ," The limits of the HS tables are 2 <= n <= 25.  Sorry.\n");
			break;
		}

		fprintf( ioQQQ, 
			" 4pJ(%ld,%ld)/n_e n_p=%11.3e\n", 
			 n1, n2, 
			 HSRate(n1,n2, ipZ,Temperature , Density) );
		/* this is check that we were in bounds */

	}

	fprintf( ioQQQ, " Thanks for interpolating on the Hummer & Storey data set!\n" );
	puts( "[Stop in DrvCaseBHS]" );
	exit(1);

}