/*pcontn print information about continuum if requested with PRINT CONTINUUM command */
#include "cddefines.h"
#define	NFLUXSV	360
#define	NXBD	9
#include "rfield.h"
#include "radius.h"
#include "nxray.h"
#include "prnline.h"
#include "nh.h"
#include "flxsav.h"
#include "pcontn.h"

void pcontn(void)
{
	long int i, 
	  i1, 
	  j, 
	  ninc, 
	  nline;
	float fluxsv[NFLUXSV], 
	  xbdsav[NXBD];

	/* energies for the x-ray bands */
	static double xraybd[NXBD + 1]={
		7.3498,
		36.8,
		73.5,
		110.3,
		147.1,
		183.8,
		220.6,
		367.6,
		551.5,
		735.3};

#	ifdef DEBUG_FUN
	fputs( "<+>pcontn()\n", debug_fp );
#	endif
	/*print information about continuum if requested with PRINT CONTINUUM command */
	/* this is number of ranges for adding x-ray bands */

	/* this stuff was always printed before version 84, and is now an option
	 * to turn on information with PRINT CONTINUUM command */

	/* the default - not turned on, just return */
	if( !PrnLine.lgPrtCont )
	{ 
#	ifdef DEBUG_FUN
		fputs( " <->pcontn()\n", debug_fp );
#	endif
		return;
	}

	/* derive x-ray fluxes in various bands for later printout
	 *
	 * >>chng 97 june 8, get rid of statement labels */
	if( xraybd[0] < rfield.anu[rfield.nflux-1] )
	{
		i1 = nxray.ipCKshell - 10;
		/* get out to lower bound of first range */
		while( rfield.anu[i1-1] < xraybd[0] && i1 < rfield.nflux )
		{
			i1 += 1;
		}
		/* now add up over that range */
		for( i=0; i < NXBD; i++ )
		{
			xbdsav[i] = 0.;
			while( rfield.anu[i1-1] < xraybd[i+1] && i1 < rfield.nflux )
			{
				xbdsav[i] += rfield.flux[i1-1] + rfield.ConOutNoInter[i1-1] + 
				  rfield.outlin[i1-1] + rfield.outcon[i1-1];
				i1 += 1;
			}
			xbdsav[i] *= (float)radius.r1r0sq;
		}
	}
	else
	{
		/* continuum not defined at x-ray energies, but zero it for safety */
		for( i=0; i < NXBD; i++ )
		{
			xbdsav[i] = 0.;
		}
	}

	/* now print x-ray flux (photons s^-1, flux or lum) in various bands */
	if( xbdsav[0] > 0 )
	{
		fprintf( ioQQQ, "\n 0.1-0.5kev:" );
		for(i=0; i < NXBD; i++)
			fprintf( ioQQQ, "%8.2e", xbdsav[i] );
		fprintf( ioQQQ, " 0.5-1.0kev:\n" );
	}

	/* renorm fLUX array before printing it */
	if( nh.ipHn[0][IP1S] - nh.ipHn[0][2] + 1 > NFLUXSV )
	{
		fprintf( ioQQQ, " PCONTN: not enough cells in FluxSave, need%4ld\n", 
		  nh.ipHn[0][IP1S] - nh.ipHn[0][2] + 1 );
		puts( "[Stop in pcontn]" );
		exit(1);
	}

	for( i=nh.ipHn[0][2]-1; i < nh.ipHn[0][IP1S]; i++ )
	{
		if( flxsav.FluxSave[i] > 0. )
		{
			fluxsv[i-nh.ipHn[0][2]+1] = (float)((rfield.flux[i] + rfield.outlin[i] + 
			  rfield.outcon[i] + rfield.ConOutNoInter[i] )*radius.r1r0sq/
			  flxsav.FluxSave[i]);
		}
		else
		{
			fluxsv[i-nh.ipHn[0][2]+1] = 0.;
		}
	}
	/* renormalize whole continuum */
	for( i=0; i < rfield.nflux; i++ )
	{
		rfield.flux[i] = (float)(((rfield.flux[i] + rfield.ConOutNoInter[i] + 
		  rfield.outcon[i]/rfield.anu[i])/rfield.widflx[i] + rfield.outlin[i])*
		  radius.r1r0sq);
	}

	fprintf( ioQQQ, 
		"\n\n                                                        Normalised continuum\n" );

	for( i=nh.ipHn[0][2]; i <= nh.ipHn[0][IP1S]; i += 3 )
	{
		fprintf( ioQQQ, "%7.3f%6.3f", rfield.anu[i-1], fluxsv[i-nh.ipHn[0][2]] );
	}
	fprintf( ioQQQ, "\n" );

	fprintf( ioQQQ, 
		"\n                                                  Emergent continuum - phot/ryd/cm2/s (r in)\n" );

	nline = ((rfield.nflux - nh.ipHn[0][2] - 1)/4 + 7)/7;
	ninc = nline*4;
	for( j=0; j < nline; j++ )
	{
		i1 = j*4 + nh.ipHn[0][2];
		fprintf( ioQQQ, " " );

		for( i=i1; i<rfield.nflux; i = i + ninc)
		{
			fprintf( ioQQQ, "%6.3f%10.2e", rfield.anu[i], 
			  rfield.flux[i] + rfield.outlin[i] + rfield.outcon[i] + 
			  rfield.ConOutNoInter[i]);
		}
		fprintf( ioQQQ, "\n" );
	}


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