/*ReadTable called by TABLE READ to read in continuum from PUNCH TRANSMITTED CONTINUUM */
#include "cddefines.h"
#include "physconst.h"
#include "rfield.h"
#include "opac.h"
#include "cnread.h"
#include "ffmtread.h"
#include "readtable.h"

void ReadTable(FILE * io)
{
	char chLine[140];
	int lgEOL;
	long int i, 
	  nPoints;
	double Differ, 
	  EnerLast;

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

	/* read in first two lines of header */
	if( NULL==fgets( chLine , sizeof(chLine) , io ) )
	{
		fprintf( ioQQQ, " error 1 reading input continuum file.\n" );
		puts( "[Stop in readtable]" );
		exit(1);
	}

	/* first number on this line is version number, confirm that
	 * this is greater than 90.01 due to bug in punch tran con in 90.00  */
	i = 1;
	/* first number should be version number of code */
	Differ = FFmtRead(chLine,&i,76,&lgEOL);
	if( Differ < 91.0 )
	{
		fprintf( ioQQQ, "%80.80s\n", chLine );
		fprintf( ioQQQ, " This transmitted continuum file was produced by an older version of Cloudy.  \n" );
		fprintf( ioQQQ, " It must be recreated with version 90.04 or later.\n" );
		fprintf( ioQQQ, " The version I read in was%6.2f\n", 
		  Differ );
		fprintf( ioQQQ, " Sorry.\n" );
		puts( "[Stop in readtable]" );
		exit(1);
	}
	
	if( NULL==fgets( chLine , sizeof(chLine) , io ) )
	{
		fprintf( ioQQQ, " error 2 reading input continuum file.\n" );
		puts( "[Stop in readtable]" );
		exit(1);
	}

	/* now read in the file of numbers */
	i = 0;
	/* keep reading until we hit eol or run out of room in the array */
	while( (NULL!=fgets(chLine,sizeof(chLine),io)) && (i<NCELL) )
	{
		sscanf( chLine, "%f%f ", &opac.tmn[i], &cnread.ConTabRead[i] );
		++i;
	}
	/* put pointer at last good value */
	nPoints = i - 1;

	/* check that sane number of values entered */
	if( nPoints < 1 )
	{
		fprintf( ioQQQ, " ReadTable, fort file for TABLE READ has too few points, =%5ld\n", 
		  nPoints );
		puts( "[Stop in readtable]" );
		exit(1);
	}

	/* check on units of energy scale, convert to Rydbergs if necessary
	 * tmn is scale read in from punch file, anu is correct scale
	 * EnerLast is energy of last point in rydbergs */
	EnerLast = opac.tmn[nPoints];
	if( fabs(opac.tmn[0]-rfield.anu[0])/rfield.anu[0] > 1e-3 )
	{
		/* first energy does not appear to have been in Rydbergs */
		if( opac.tmn[0] < 0. )
		{
			fprintf( ioQQQ, " ReadTable:  the first energy in table read file is not positive.  Something is wrong.\n" );
			puts( "[Stop in readtable]" );
			exit(1);
		}
		else if( fabs(opac.tmn[0]-RYDLAM/rfield.anu[0]*1e-4)/opac.tmn[0] < 
		  1e-3 )
		{
			/* wavelength in microns */
			EnerLast = RYDLAM/opac.tmn[nPoints]/1e4;
		}
		else if( fabs(opac.tmn[0]-RYDLAM/rfield.anu[0])/opac.tmn[0] < 
		  1e-3 )
		{
			/* wavelength in Angstroms */
			EnerLast = RYDLAM/opac.tmn[nPoints];
		}
		else if( fabs(opac.tmn[0]-rfield.anu[0]*EVRYD*1e-3)/opac.tmn[0] < 
		  1e-3 )
		{
			/* wavelength in keV */
			EnerLast = opac.tmn[nPoints]/EVRYD/1e-3;
		}
		else if( fabs(opac.tmn[0]-rfield.anu[0]*EVRYD)/opac.tmn[0] < 
		  1e-3 )
		{
			/* wavelength in eV */
			EnerLast = opac.tmn[nPoints]/EVRYD;
		}
		else
		{
			fprintf( ioQQQ, " ReadTable:  the energy scale in the table read file makes no sense to me.\n" );
			puts( "[Stop in readtable]" );
			exit(1);
		}
	}

	/* now check now the energies of the highest points agree */
	Differ = fabs(EnerLast-rfield.anu[nPoints])/rfield.anu[nPoints];
	if( Differ > 0.001 )
	{
		fprintf( ioQQQ, " ReadTable: The energy mesh of the file read in by the TABLE READ command does not agree with this version of Cloudy.\n" );
		fprintf( ioQQQ, " ReadTable: Was the file generated by an older version of the code?\n" );
		fprintf( ioQQQ, " ReadTable: Number of points read in=%5ld\n", 
		  nPoints );
		fprintf( ioQQQ, " ReadTable: input, internal energies=%12.4e%12.4e\n", 
		  opac.tmn[nPoints], rfield.anu[nPoints] );
		puts( "[Stop in readtable]" );
		exit(1);
	}

	/* make sure rest of array has valid zeros */
	for( i=nPoints + 1; i < NCELL; i++ )
	{
		cnread.ConTabRead[i] = 0.;
	}

	fclose(io);


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