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

double pintr(double penlo, 
  double penhi)
{
	int lgMore;
	long int i, 
	  ip, 
	  j;
	double en1, 
	  en2, 
	  fsum, 
	  pintr_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( "<+>pintr()\n", debug_fp );
#	endif

	/* computes log of luminosity in radiation over some intergal
	 * answer is in Ryd per sec */

	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 - 20; 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)*wanu;
				fsum += wfun;
			}
			sum += fsum*rfield.widflx[i];
		}

	}
	else
	{
		lgMore = TRUE;
		en1 = penlo/SKIP;
		while( lgMore )
		{
			en1 *= SKIP;
			en2 = en1*SKIP;
			if( en2 >= penhi )
			{
				en2 = penhi;
				lgMore = FALSE;
			}
			sum += qg32(en1,en2,pow1);
		}
	}

	if( sum > 0. )
	{
		pintr_v = log10(sum);
	}
	else
	{
		pintr_v = -38.;
	}


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