/*
outputCross.c

output cross section
*/

#include "oth.h"
#include "complex.h"
#include "oth_prototype.h"


static void _outputCross(FILE *fp, const double bcs[], const double fcs[])
{
	fprintf(fp, "%s\n", "=== cross section ===");
	fprintf(fp, "  %s\n", "frequency[Hz] backward[m*m]    forward[m*m]");
	for (int ifreq = 0; ifreq < NFrequency; ifreq++) {
		fprintf(fp, "%13.5e%16.5e%16.5e\n", Frequency[ifreq], bcs[ifreq], fcs[ifreq]);
	}
	fflush(fp);
}


void outputCross(FILE *fp)
{
	// setup
	setup_farfield();

	// alloc
	double *ffctr = (double *)malloc(NFrequency * sizeof(double));
	double *bcs   = (double *)malloc(NFrequency * sizeof(double));
	double *fcs   = (double *)malloc(NFrequency * sizeof(double));

	// far field factor
	farfactor(ffctr);

	// calculation
	for (int ifreq = 0; ifreq < NFrequency; ifreq++) {
		d_complex_t etheta, ephi;
		// BCS
		farfield(ifreq,       Planewave.theta,       Planewave.phi, ffctr[ifreq], &etheta, &ephi);
		bcs[ifreq] = d_norm(etheta) + d_norm(ephi);
		// FCS
		farfield(ifreq, 180 - Planewave.theta, 180 + Planewave.phi, ffctr[ifreq], &etheta, &ephi);
		fcs[ifreq] = d_norm(etheta) + d_norm(ephi);
	}

	// output
	_outputCross(stdout, bcs, fcs);
	_outputCross(fp,     bcs, fcs);

	// free
	free(ffctr);
	free(bcs);
	free(fcs);
}
