/*qintr integrates Q for any continuum between two limits, used for normalization */
#include <string.h>
#include "cddefines.h"
#include "qg32.h"
#define	SKIP	2.
#include "rfield.h"
#include "incidspec.h"
#include "bounds.h"
#include "ipoint.h"
#include "ffun1.h"
#include "qintr.h"

double qintr(float *qenlo, 
  float *qenhi)
{
	int lgMore;
	long int i, 
	  ip, 
	  ipHi, 
	  ipLo, 
	  j;
	double en1, 
	  en2, 
	  fsum, 
	  hold, 
	  qintr_v, 
	  sum, 
	  wanu, 
	  wfun;
	static double aweigh[4], 
	  fweigh[4];
	static int _aini = 1;
	if( _aini ){ /* Do 1 TIME INITIALIZATIONS! */
		aweigh[0] = -0.4305682;
		aweigh[1] = -0.1699905;
		aweigh[2] = 0.1699905;
		aweigh[3] = 0.4305682;
		fweigh[0] = 0.1739274;
		fweigh[1] = 0.3260726;
		fweigh[2] = 0.3260726;
		fweigh[3] = 0.1739274;
		_aini = 0;
	}

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

	/* returns LOG of number of photons over energy interval */
	sum = 0.;
	if( strcmp(spnorm.chSpType[IncidSpec.ipspec],"LASER") == 0 )
	{
		/*  laser is special since delta function */
		ip = ipoint(IncidSpec.slope[IncidSpec.ipspec]);
		for( i=ip - 21; i < (ip + 20); i++ )
		{
			fsum = 0.;
			for( j=0; j < 4; j++ )
			{
				wanu = rfield.anu[i] + rfield.widflx[i]*aweigh[j];
				wfun = fweigh[j]*ffun1(wanu);
				fsum += wfun;
			}
			sum += fsum*rfield.widflx[i];
		}
	}

	else if( strcmp(spnorm.chSpType[IncidSpec.ipspec],"INTER") == 0 )
	{
		/* interpolate is special since might be very chaotic */
		ipLo = ipoint(*qenlo);
		ipHi = ipoint(*qenhi);
		/* there are possible edge effects due to finite size of cells 
		 * do start and end as special cases */
		/* this is the lowest energy cell, but not below continuum array*/
		hold = MAX2(bounds.emm,rfield.anu[ipLo]);
		sum = ffun1(hold)*rfield.widflx[ipLo];
		if( ipLo != ipHi )
		{
			/* this is highest energy point, but not above limit */
			hold = MIN2(bounds.egamry,rfield.anu[ipHi-1]);
			sum += ffun1(hold)*rfield.widflx[ipHi-1];
		}
		for( i=ipLo; i < (ipHi - 1); i++ )
		{
			sum += ffun1(rfield.anu[i])*rfield.widflx[i];
		}
	}
	else
	{
		lgMore = TRUE;
		en1 = *qenlo/SKIP;
		while( lgMore )
		{
			en1 *= SKIP;
			en2 = en1*SKIP;
			if( en2 >= *qenhi )
			{
				en2 = *qenhi;
				lgMore = FALSE;
			}
			sum += qg32(en1,en2,ffun1);
		}
	}

	if( sum <= 0. )
	{
		fprintf( ioQQQ, " Photon number sum in QINTR is%10.3e\n", 
		  sum );
		fprintf( ioQQQ, " This source has no ionizing radiation, and the number of ionizing photons was specified.\n" );
		fprintf( ioQQQ, " This was continuum source number%3ld\n", 
		  IncidSpec.ipspec );
		fprintf( ioQQQ, " Sorry, but I cannot go on.  ANU and FLUX arrays follow.  Enjoy.\n" );
		for( i=1; i <= NCELL; i++ )
		{
			fprintf( ioQQQ, "%10.2e%10.2e", rfield.anu[i-1], 
			  rfield.flux[i-1] );
		}
		fprintf( ioQQQ, "\n" );
		puts( "[Stop in qintr]" );
		exit(1);
	}
	else
	{
		qintr_v = log10(sum);
	}


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