/*pgaunt called by punch gaunts command to output gaunt factors */
#include "cddefines.h"
#define	NENR	31
#define	NTE	21
#include "logte.h"
#include "gffsub.h"
#include "pgaunt.h"

void pgaunt(FILE* io)
{
	long int i, 
	  ier[NENR][NTE], 
	  if1, 
	  ite, 
	  j;

	float ener[NENR], 
	  g[NENR][NTE], 
	  gfac, 
	  ste[NTE], 
	  z;

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

	/* this routine is called from the PUNCH GAUNTS command
	 * it drives the gaunt factor routine to punch gaunts over full range */

	for( i=0; i < NTE; i++ )
	{
		ste[i] = (float)(0.5*i);
	}

	for( i=0; i < NENR; i++ )
	{
		ener[i] = (float)(0.5*i - 5.);
	}

	/* ener is log of energy */
	for( ite=0; ite < NTE; ite++ )
	{
		logte.alogte = ste[ite];
		/* nuclear charge */
		z = (float)log10(1.);

		for( j=0; j < NENR; j++ )
		{
			gffsub(z,&ener[j],(float*)&gfac,1,1,&if1);
			g[j][ite] = gfac;
			ier[j][ite] = if1;
		}
	}

	/* now punch out the results */
	fprintf( io, "      " );
	for( i=1; i <= NTE; i++ )
	{
		fprintf( io, "%6.1f", ste[i-1] );
	}
	fprintf( io, "\n" );

	for( j=0; j < NENR; j++ )
	{
		fprintf( io, "%6.2f", ener[j] );
		for( ite=0; ite < NTE; ite++ )
		{
			fprintf( io, "%6.2f", log10(g[j][ite]) );
		}
		fprintf( io, "\n" );
	}

	fprintf( io, "      " );
	for( i=0; i < NTE; i++ )
	{
		fprintf( io, "%6.1f", ste[i] );
	}
	fprintf( io, "\n" );

	for( j=0; j < NENR; j++ )
	{
		fprintf( io, "%6.3f", ener[j] );
		for( ite=0; ite < NTE; ite++ )
		{
			fprintf( io, "%6ld", ier[j][ite] );
		}
		fprintf( io, "\n" );
	}

	/* now do the same thing as triplets */
	fprintf( io, "      " );
	for( i=0; i < NTE; i++ )
	{
		fprintf( io, "%6.1f", ste[i] );
	}
	fprintf( io, "\n" );

	for( i=0; i < NTE; i++ )
	{
		for( j=0; j < NENR; j++ )
		{
			fprintf( io, "%6.2f%6.2f%6.2f\n", ste[i], ener[j], 
			  log10(g[j][i]) );
		}
	}


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