/*twopht hydrogen two photon emission, including induced emission */
#include "cddefines.h"
#include "taulines.h"
#include "ip2phtcom.h"
#include "rfield.h"
#include "rin2ph.h"
#include "em2nu.h"
#include "trace.h"
#include "max2pht.h"
#include "fluors.h"
#include "twopht.h"

void twopht(void)
{
	long int i, 
	  limit;
	double sum;

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

	rin2ph.ri2s1s = 0.;
	if( !fluors.lgFluor )
	{ 
#	ifdef DEBUG_FUN
		fputs( " <->twopht()\n", debug_fp );
#	endif
		return;
	}

	for( i=0; i < ip2phtCom.ipHalfLya; i++ )
	{
		rin2ph.ri2s1s += rfield.SummedOcc[i]*em2nu.Hy2nu[i]*
			rfield.occnum[ip2phtCom.ip2pht[i]-1];
	}

	/* integral is symmetric */
	rin2ph.ri2s1s *= 2.;

	/* remember maximum value */
	Max2Pht.xMax2Pht = (float)MAX2(Max2Pht.xMax2Pht,rin2ph.ri2s1s);

	rin2ph.ri1s2s = rin2ph.ri2s1s*1.f;

	if( trace.lgTrace && trace.lgBug2nu )
	{
		limit = HydroLines[0][IP2P][IP1S].ipCont;
		for( i=1; i <= limit; i += 10 )
		{
			fprintf( ioQQQ, "%12.4e%12.4e", rfield.anu[i-1], 
			  rfield.occnum[i-1] );
		}
		fprintf( ioQQQ, "\n" );

		/* sum two-nu A's to check OK */
		sum = 0.;
		for( i=0; i < limit; i++ )
		{
			sum += em2nu.Hy2nu[i];
		}
		fprintf( ioQQQ, "       TWOPHT finds induced two photon rate =%12.3e total A=%12.3e\n", 
		  rin2ph.ri2s1s, sum );
	}


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