tools.hpp

//////////////////////////
// tools.hpp
//////////////////////////
//
// Collection of analysis tools for gevolution
//
// Author: Julian Adamek (Université de Genève & Observatoire de Paris & Queen Mary University of London)
//
// Last modified: April 2019
//
//////////////////////////

#ifndef TOOLS_HEADER
#define TOOLS_HEADER

#ifndef Cplx
#define Cplx Imag
#endif

#define KTYPE_GRID 0
#define KTYPE_LINEAR 1

using namespace std;
using namespace LATfield2;

#ifdef FFT3D
//////////////////////////
// extractCrossSpectrum
//////////////////////////
// Description:
//   generates the cross spectrum for two Fourier images
//
// Arguments:
//   fld1FT     reference to the first Fourier image for which the cross spectrum should be extracted
//   fld2FT     reference to the second Fourier image for which the cross spectrum should be extracted
//   kbin       allocated array that will contain the central k-value for the bins
//   power      allocated array that will contain the average power in each bin
//   kscatter   allocated array that will contain the k-scatter for each bin
//   pscatter   allocated array that will contain the scatter in power for each bin
//   occupation allocated array that will count the number of grid points contributing to each bin
//   numbin     number of bins (minimum size of all arrays)
//   ktype      flag indicating which definition of momentum to be used
//                  0: grid momentum
//                  1: linear (default)
//   comp1      for component-wise cross spectra, the component for the first field (ignored if negative)
//   comp2      for component-wise cross spectra, the component for the second field (ignored if negative)
//
// Returns:
//
//////////////////////////

void extractCrossSpectrum(Field<Cplx> &fld1FT, Field<Cplx> &fld2FT, Real *kbin, Real *power, Real *kscatter, Real *pscatter, int *occupation, const int numbins, const bool deconvolve = true, const int ktype = KTYPE_LINEAR, const int comp1 = -1, const int comp2 = -1)
{
	int i, weight;
	const int linesize = fld1FT.lattice().size(1);
	Real *typek2;
	Real *sinc;
	Real k2max, k2, s;
	rKSite k(fld1FT.lattice());
	Cplx p;

	typek2 = (Real *)malloc(linesize * sizeof(Real));
	sinc = (Real *)malloc(linesize * sizeof(Real));

	if (ktype == KTYPE_GRID)
	{
		for (i = 0; i < linesize; i++)
		{
			typek2[i] = 2. * (Real)linesize * sin(M_PI * (Real)i / (Real)linesize);
			typek2[i] *= typek2[i];
		}
	}
	else
	{
		for (i = 0; i <= linesize / 2; i++)
		{
			typek2[i] = 2. * M_PI * (Real)i;
			typek2[i] *= typek2[i];
		}
		for (; i < linesize; i++)
		{
			typek2[i] = 2. * M_PI * (Real)(linesize - i);
			typek2[i] *= typek2[i];
		}
	}

	sinc[0] = 1.;
	if (deconvolve)
	{
		for (i = 1; i <= linesize / 2; i++)
		{
			sinc[i] = sin(M_PI * (float)i / (float)linesize) * (float)linesize / (M_PI * (float)i);
		}
	}
	else
	{
		for (i = 1; i <= linesize / 2; i++)
		{
			sinc[i] = 1.;
		}
	}
	for (; i < linesize; i++)
	{
		sinc[i] = sinc[linesize - i];
	}

	k2max = 3. * typek2[linesize / 2];

	for (i = 0; i < numbins; i++)
	{
		kbin[i] = 0.;
		power[i] = 0.;
		kscatter[i] = 0.;
		pscatter[i] = 0.;
		occupation[i] = 0;
	}

	for (k.first(); k.test(); k.next())
	{
		if (k.coord(0) == 0 && k.coord(1) == 0 && k.coord(2) == 0)
			continue;
		else if (k.coord(0) == 0)
			weight = 1;
		else if ((k.coord(0) == linesize / 2) && (linesize % 2 == 0))
			weight = 1;
		else
			weight = 2;

		k2 = typek2[k.coord(0)] + typek2[k.coord(1)] + typek2[k.coord(2)];
		s = sinc[k.coord(0)] * sinc[k.coord(1)] * sinc[k.coord(2)];
		s *= s;

		if (comp1 >= 0 && comp2 >= 0 && comp1 < fld1FT.components() && comp2 < fld2FT.components())
		{
			p = fld1FT(k, comp1) * fld2FT(k, comp2).conj();
		}
		else if (fld1FT.symmetry() == LATfield2::symmetric)
		{
			p = fld1FT(k, 0, 1) * fld2FT(k, 0, 1).conj();
			p += fld1FT(k, 0, 2) * fld2FT(k, 0, 2).conj();
			p += fld1FT(k, 1, 2) * fld2FT(k, 1, 2).conj();
			p *= 2.;
			p += fld1FT(k, 0, 0) * fld2FT(k, 0, 0).conj();
			p += fld1FT(k, 1, 1) * fld2FT(k, 1, 1).conj();
			p += fld1FT(k, 2, 2) * fld2FT(k, 2, 2).conj();
		}
		else
		{
			p = Cplx(0., 0.);
			for (i = 0; i < fld1FT.components(); i++)
				p += fld1FT(k, i) * fld2FT(k, i).conj();
		}

		i = (int)floor((double)((Real)numbins * sqrt(k2 / k2max)));
		if (i < numbins)
		{
			kbin[i] += weight * sqrt(k2);
			kscatter[i] += weight * k2;
			power[i] += weight * p.real() * k2 * sqrt(k2) / s;
			pscatter[i] += weight * p.real() * p.real() * k2 * k2 * k2 / s / s;
			occupation[i] += weight;
		}
	}

	free(typek2);
	free(sinc);

	if (parallel.isRoot())
	{
#ifdef SINGLE
		MPI_Reduce(MPI_IN_PLACE, (void *)kbin, numbins, MPI_FLOAT, MPI_SUM, 0, parallel.lat_world_comm());
		MPI_Reduce(MPI_IN_PLACE, (void *)kscatter, numbins, MPI_FLOAT, MPI_SUM, 0, parallel.lat_world_comm());
		MPI_Reduce(MPI_IN_PLACE, (void *)power, numbins, MPI_FLOAT, MPI_SUM, 0, parallel.lat_world_comm());
		MPI_Reduce(MPI_IN_PLACE, (void *)pscatter, numbins, MPI_FLOAT, MPI_SUM, 0, parallel.lat_world_comm());
#else
		MPI_Reduce(MPI_IN_PLACE, (void *)kbin, numbins, MPI_DOUBLE, MPI_SUM, 0, parallel.lat_world_comm());
		MPI_Reduce(MPI_IN_PLACE, (void *)kscatter, numbins, MPI_DOUBLE, MPI_SUM, 0, parallel.lat_world_comm());
		MPI_Reduce(MPI_IN_PLACE, (void *)power, numbins, MPI_DOUBLE, MPI_SUM, 0, parallel.lat_world_comm());
		MPI_Reduce(MPI_IN_PLACE, (void *)pscatter, numbins, MPI_DOUBLE, MPI_SUM, 0, parallel.lat_world_comm());
#endif
		MPI_Reduce(MPI_IN_PLACE, (void *)occupation, numbins, MPI_INT, MPI_SUM, 0, parallel.lat_world_comm());

		for (i = 0; i < numbins; i++)
		{
			if (occupation[i] > 0)
			{
				kscatter[i] = sqrt(kscatter[i] * occupation[i] - kbin[i] * kbin[i]) / occupation[i];
				if (!isfinite(kscatter[i]))
					kscatter[i] = 0.;
				kbin[i] = kbin[i] / occupation[i];
				power[i] /= occupation[i];
				pscatter[i] = sqrt(pscatter[i] / occupation[i] - power[i] * power[i]);
				if (!isfinite(pscatter[i]))
					pscatter[i] = 0.;
			}
		}
	}
	else
	{
#ifdef SINGLE
		MPI_Reduce((void *)kbin, NULL, numbins, MPI_FLOAT, MPI_SUM, 0, parallel.lat_world_comm());
		MPI_Reduce((void *)kscatter, NULL, numbins, MPI_FLOAT, MPI_SUM, 0, parallel.lat_world_comm());
		MPI_Reduce((void *)power, NULL, numbins, MPI_FLOAT, MPI_SUM, 0, parallel.lat_world_comm());
		MPI_Reduce((void *)pscatter, NULL, numbins, MPI_FLOAT, MPI_SUM, 0, parallel.lat_world_comm());
#else
		MPI_Reduce((void *)kbin, NULL, numbins, MPI_DOUBLE, MPI_SUM, 0, parallel.lat_world_comm());
		MPI_Reduce((void *)kscatter, NULL, numbins, MPI_DOUBLE, MPI_SUM, 0, parallel.lat_world_comm());
		MPI_Reduce((void *)power, NULL, numbins, MPI_DOUBLE, MPI_SUM, 0, parallel.lat_world_comm());
		MPI_Reduce((void *)pscatter, NULL, numbins, MPI_DOUBLE, MPI_SUM, 0, parallel.lat_world_comm());
#endif
		MPI_Reduce((void *)occupation, NULL, numbins, MPI_INT, MPI_SUM, 0, parallel.lat_world_comm());
	}
}

//////////////////////////
// extractPowerSpectrum
//////////////////////////
// Description:
//   generates the power spectrum for a Fourier image
//
// Arguments:
//   fldFT      reference to the Fourier image for which the power spectrum should be extracted
//   kbin       allocated array that will contain the central k-value for the bins
//   power      allocated array that will contain the average power in each bin
//   kscatter   allocated array that will contain the k-scatter for each bin
//   pscatter   allocated array that will contain the scatter in power for each bin
//   occupation allocated array that will count the number of grid points contributing to each bin
//   numbin     number of bins (minimum size of all arrays)
//   ktype      flag indicating which definition of momentum to be used
//                  0: grid momentum
//                  1: linear (default)
//
// Returns:
//
//////////////////////////

void extractPowerSpectrum(Field<Cplx> &fldFT, Real *kbin, Real *power, Real *kscatter, Real *pscatter, int *occupation, const int numbins, const bool deconvolve = true, const int ktype = KTYPE_LINEAR)
{
	extractCrossSpectrum(fldFT, fldFT, kbin, power, kscatter, pscatter, occupation, numbins, deconvolve, ktype);
}
#endif

//////////////////////////
// writePowerSpectrum
//////////////////////////
// Description:
//   writes power spectra as tabulated data into ASCII file
//
// Arguments:
//   kbin           array containing the central values of k for each bin
//   power          array containing the central values of P(k) for each bin
//   kscatter       array containing the statistical error on k for each bin
//   pscatter       array containing the statistical error on P(k) for each bin
//   occupation     array containing the number of k-modes contributing to each bin
//   numbins        total number of bins (length of the arrays)
//   rescalek       unit conversion factor for k
//   rescalep       unit conversion factor for P(k)
//   filename       output file name
//   description    descriptive header
//   a              scale factor for this spectrum
//   z_target       target redshift for this output (used only if EXACT_OUTPUT_REDSHIFTS is defined)
//
// Returns:
//
//////////////////////////

void writePowerSpectrum(Real *kbin, Real *power, Real *kscatter, Real *pscatter, int *occupation, const int numbins, const Real rescalek, const Real rescalep, const char *filename, const char *description, double a, const double z_target = -1)
{
	if (parallel.isRoot())
	{
#ifdef EXACT_OUTPUT_REDSHIFTS
		Real *power2 = (Real *)malloc(numbins * sizeof(Real));

		for (int i = 0; i < numbins; i++)
			power2[i] = power[i] / rescalep;

		if (1. / a < z_target + 1.)
		{
			FILE *infile = fopen(filename, "r");
			double weight = 1.;
			int count = 0;
			if (infile != NULL)
			{
				fscanf(infile, "%*[^\n]\n");
				if (fscanf(infile, "# redshift z=%lf\n", &weight) != 1)
				{
					COUT << " error parsing power spectrum file header for interpolation (EXACT_OUTPUT_REDSHIFTS)" << endl;
					weight = 1.;
				}
				else
				{
					weight = (weight - z_target) / (1. + weight - 1. / a);
					fscanf(infile, "%*[^\n]\n");
					for (int i = 0; i < numbins; i++)
					{
						if (occupation[i] > 0)
						{
#ifdef SINGLE
							if (fscanf(infile, " %*e %e %*e %*e %*d \n", power2 + i) != 1)
#else
							if (fscanf(infile, " %*e %le %*e %*e %*d \n", power2 + i) != 1)
#endif
							{
								COUT << " error parsing power spectrum file data " << i << " for interpolation (EXACT_OUTPUT_REDSHIFTS)" << endl;
								break;
							}
							else
								count++;
						}
					}
				}
				fclose(infile);

				for (int i = 0; i < numbins; i++)
					power2[i] = (1. - weight) * power2[i] + weight * power[i] / rescalep;

				a = 1. / (z_target + 1.);
			}
		}
#endif // EXACT_OUTPUT_REDSHIFTS
		FILE *outfile = fopen(filename, "w");
		if (outfile == NULL)
		{
			COUT << " error opening file for power spectrum output!" << endl;
		}
		else
		{
			fprintf(outfile, "# %s\n", description);
			fprintf(outfile, "# redshift z=%f\n", (1. / a) - 1.);
			fprintf(outfile, "# k              Pk             sigma(k)       sigma(Pk)      count\n");
			for (int i = 0; i < numbins; i++)
			{
				if (occupation[i] > 0)
#ifdef EXACT_OUTPUT_REDSHIFTS
					fprintf(outfile, "  %e   %e   %e   %e   %d\n", kbin[i] / rescalek, power2[i], kscatter[i] / rescalek, pscatter[i] / rescalep / sqrt(occupation[i]), occupation[i]);
#else
					fprintf(outfile, "  %e   %e   %e   %e   %d\n", kbin[i] / rescalek, power[i] / rescalep, kscatter[i] / rescalek, pscatter[i] / rescalep / sqrt(occupation[i]), occupation[i]);
#endif
			}
			fclose(outfile);
		}
#ifdef EXACT_OUTPUT_REDSHIFTS
		free(power2);
#endif
	}
}

//////////////////////////
// computeVectorDiagnostics
//////////////////////////
// Description:
//   computes some diagnostics for the spin-1 perturbation
//
// Arguments:
//   Bi         reference to the real-space vector field to analyze
//   mdivB      will contain the maximum value of the divergence of Bi
//   mcurlB     will contain the maximum value of the curl of Bi
//
// Returns:
//
//////////////////////////

void computeVectorDiagnostics(Field<Real> &Bi, Real &mdivB, Real &mcurlB)
{
	Real b1, b2, b3, b4;
	const Real linesize = (Real)Bi.lattice().sizeLocal(0);
	Site x(Bi.lattice());

	mdivB = 0.;
	mcurlB = 0.;

	for (x.first(); x.test(); x.next())
	{
		b1 = fabs((Bi(x, 0) - Bi(x - 0, 0)) + (Bi(x, 1) - Bi(x - 1, 1)) + (Bi(x, 2) - Bi(x - 2, 2))) * linesize;
		if (b1 > mdivB)
			mdivB = b1;
		b1 = 0.5 * (Bi(x, 0) + Bi(x + 0, 1) - Bi(x + 1, 0) - Bi(x, 1) + Bi(x + 2, 0) + Bi(x + 0 + 2, 1) - Bi(x + 1 + 2, 0) - Bi(x + 2, 1)) * linesize;
		b2 = 0.5 * (Bi(x, 0) + Bi(x + 0, 2) - Bi(x + 2, 0) - Bi(x, 2) + Bi(x + 1, 0) + Bi(x + 0 + 1, 2) - Bi(x + 2 + 1, 0) - Bi(x + 1, 2)) * linesize;
		b3 = 0.5 * (Bi(x, 2) + Bi(x + 2, 1) - Bi(x + 1, 2) - Bi(x, 1) + Bi(x + 0, 2) + Bi(x + 2 + 0, 1) - Bi(x + 1 + 0, 2) - Bi(x + 0, 1)) * linesize;
		b4 = sqrt(b1 * b1 + b2 * b2 + b3 * b3);
		if (b4 > mcurlB)
			mcurlB = b4;
	}

	parallel.max<Real>(mdivB);
	parallel.max<Real>(mcurlB);
}

//////////////////////////
// computeTensorDiagnostics
//////////////////////////
// Description:
//   computes some diagnostics for the spin-2 perturbation
//
// Arguments:
//   hij        reference to the real-space tensor field to analyze
//   mdivh      will contain the maximum value of the divergence of hij
//   mtraceh    will contain the maximum value of the trace of hij
//   mnormh     will contain the maximum value of the norm of hij
//
// Returns:
//
//////////////////////////

void computeTensorDiagnostics(Field<Real> &hij, Real &mdivh, Real &mtraceh, Real &mnormh)
{
	Real d1, d2, d3;
	const Real linesize = (Real)hij.lattice().sizeLocal(0);
	Site x(hij.lattice());

	mdivh = 0.;
	mtraceh = 0.;
	mnormh = 0.;

	for (x.first(); x.test(); x.next())
	{
		d1 = (hij(x + 0, 0, 0) - hij(x, 0, 0) + hij(x, 0, 1) - hij(x - 1, 0, 1) + hij(x, 0, 2) - hij(x - 2, 0, 2)) * linesize;
		d2 = (hij(x + 1, 1, 1) - hij(x, 1, 1) + hij(x, 0, 1) - hij(x - 0, 0, 1) + hij(x, 1, 2) - hij(x - 2, 1, 2)) * linesize;
		d3 = (hij(x + 2, 2, 2) - hij(x, 2, 2) + hij(x, 0, 2) - hij(x - 0, 0, 2) + hij(x, 1, 2) - hij(x - 1, 1, 2)) * linesize;
		d1 = sqrt(d1 * d1 + d2 * d2 + d3 * d3);
		if (d1 > mdivh)
			mdivh = d1;
		d1 = fabs(hij(x, 0, 0) + hij(x, 1, 1) + hij(x, 2, 2));
		if (d1 > mtraceh)
			mtraceh = d1;
		d1 = sqrt(hij(x, 0, 0) * hij(x, 0, 0) + 2. * hij(x, 0, 1) * hij(x, 0, 1) + 2. * hij(x, 0, 2) * hij(x, 0, 2) + hij(x, 1, 1) * hij(x, 1, 1) + 2. * hij(x, 1, 2) * hij(x, 1, 2) + hij(x, 2, 2) * hij(x, 2, 2));
		if (d1 > mnormh)
			mnormh = d1;
	}

	parallel.max<Real>(mdivh);
	parallel.max<Real>(mtraceh);
	parallel.max<Real>(mnormh);
}

//////////////////////////
// findIntersectingLightcones
//////////////////////////
// Description:
//   determines periodic copies of light cone vertex for which the present
//   look-back interval may overlap with a given spatial domain
//
// Arguments:
//   lightcone  reference to structure describing light cone geometry
//   outer      outer (far) limit of look-back interval
//   inner      inner (close) limit of look-back interval
//   domain     array of domain boundaries
//   vertex     will contain array of relevant vertex locations
//
// Returns:
//   number of vertices found
//
//////////////////////////

int findIntersectingLightcones(lightcone_geometry &lightcone, double outer, double inner, double *domain, double vertex[MAX_INTERSECTS][3])
{
	int range = (int)ceil(outer) + 1;
	int u, v, w, n = 0;
	double corner[8][3];
	double rdom, dist;

	corner[0][0] = domain[0];
	corner[0][1] = domain[1];
	corner[0][2] = domain[2];

	corner[1][0] = domain[3];
	corner[1][1] = domain[1];
	corner[1][2] = domain[2];

	corner[2][0] = domain[0];
	corner[2][1] = domain[4];
	corner[2][2] = domain[2];

	corner[3][0] = domain[3];
	corner[3][1] = domain[4];
	corner[3][2] = domain[2];

	corner[4][0] = domain[0];
	corner[4][1] = domain[1];
	corner[4][2] = domain[5];

	corner[5][0] = domain[3];
	corner[5][1] = domain[1];
	corner[5][2] = domain[5];

	corner[6][0] = domain[0];
	corner[6][1] = domain[4];
	corner[6][2] = domain[5];

	corner[7][0] = domain[3];
	corner[7][1] = domain[4];
	corner[7][2] = domain[5];

	for (u = -range; u <= range; u++)
	{
		for (v = -range; v <= range; v++)
		{
			for (w = -range; w <= range; w++)
			{
				if (n >= MAX_INTERSECTS)
				{
					cout << COLORTEXT_YELLOW << " /!\\ warning" << COLORTEXT_RESET << ": maximum number of lightcone intersects exceeds MAX_INTERSECTS = " << MAX_INTERSECTS << " for domain (" << domain[0] << ", " << domain[1] << ", " << domain[2] << ") - (" << domain[3] << ", " << domain[4] << ", " << domain[5] << "); some data may be missing in output!" << endl;
					return MAX_INTERSECTS;
				}
				vertex[n][0] = lightcone.vertex[0] + u;
				vertex[n][1] = lightcone.vertex[1] + v;
				vertex[n][2] = lightcone.vertex[2] + w;

				// first, check if domain lies outside outer sphere
				if (vertex[n][0] < domain[0])
				{
					if (vertex[n][1] < domain[1])
					{
						if (vertex[n][2] < domain[2])
						{
							if (sqrt((vertex[n][0] - corner[0][0]) * (vertex[n][0] - corner[0][0]) + (vertex[n][1] - corner[0][1]) * (vertex[n][1] - corner[0][1]) + (vertex[n][2] - corner[0][2]) * (vertex[n][2] - corner[0][2])) > outer)
								continue;
						}
						else if (vertex[n][2] > domain[5])
						{
							if (sqrt((vertex[n][0] - corner[4][0]) * (vertex[n][0] - corner[4][0]) + (vertex[n][1] - corner[4][1]) * (vertex[n][1] - corner[4][1]) + (vertex[n][2] - corner[4][2]) * (vertex[n][2] - corner[4][2])) > outer)
								continue;
						}
						else if (sqrt((vertex[n][0] - domain[0]) * (vertex[n][0] - domain[0]) + (vertex[n][1] - domain[1]) * (vertex[n][1] - domain[1])) > outer)
							continue;
					}
					else if (vertex[n][1] > domain[4])
					{
						if (vertex[n][2] < domain[2])
						{
							if (sqrt((vertex[n][0] - corner[2][0]) * (vertex[n][0] - corner[2][0]) + (vertex[n][1] - corner[2][1]) * (vertex[n][1] - corner[2][1]) + (vertex[n][2] - corner[2][2]) * (vertex[n][2] - corner[2][2])) > outer)
								continue;
						}
						else if (vertex[n][2] > domain[5])
						{
							if (sqrt((vertex[n][0] - corner[6][0]) * (vertex[n][0] - corner[6][0]) + (vertex[n][1] - corner[6][1]) * (vertex[n][1] - corner[6][1]) + (vertex[n][2] - corner[6][2]) * (vertex[n][2] - corner[6][2])) > outer)
								continue;
						}
						else if (sqrt((vertex[n][0] - domain[0]) * (vertex[n][0] - domain[0]) + (vertex[n][1] - domain[4]) * (vertex[n][1] - domain[4])) > outer)
							continue;
					}
					else
					{
						if (vertex[n][2] < domain[2])
						{
							if (sqrt((vertex[n][0] - domain[0]) * (vertex[n][0] - domain[0]) + (vertex[n][2] - domain[2]) * (vertex[n][2] - domain[2])) > outer)
								continue;
						}
						else if (vertex[n][2] > domain[5])
						{
							if (sqrt((vertex[n][0] - domain[0]) * (vertex[n][0] - domain[0]) + (vertex[n][2] - domain[5]) * (vertex[n][2] - domain[5])) > outer)
								continue;
						}
						else if (domain[0] - vertex[n][0] > outer)
							continue;
					}
				}
				else if (vertex[n][0] > domain[3])
				{
					if (vertex[n][1] < domain[1])
					{
						if (vertex[n][2] < domain[2])
						{
							if (sqrt((vertex[n][0] - corner[1][0]) * (vertex[n][0] - corner[1][0]) + (vertex[n][1] - corner[1][1]) * (vertex[n][1] - corner[1][1]) + (vertex[n][2] - corner[1][2]) * (vertex[n][2] - corner[1][2])) > outer)
								continue;
						}
						else if (vertex[n][2] > domain[5])
						{
							if (sqrt((vertex[n][0] - corner[5][0]) * (vertex[n][0] - corner[5][0]) + (vertex[n][1] - corner[5][1]) * (vertex[n][1] - corner[5][1]) + (vertex[n][2] - corner[5][2]) * (vertex[n][2] - corner[5][2])) > outer)
								continue;
						}
						else if (sqrt((vertex[n][0] - domain[3]) * (vertex[n][0] - domain[3]) + (vertex[n][1] - domain[1]) * (vertex[n][1] - domain[1])) > outer)
							continue;
					}
					else if (vertex[n][1] > domain[4])
					{
						if (vertex[n][2] < domain[2])
						{
							if (sqrt((vertex[n][0] - corner[3][0]) * (vertex[n][0] - corner[3][0]) + (vertex[n][1] - corner[3][1]) * (vertex[n][1] - corner[3][1]) + (vertex[n][2] - corner[3][2]) * (vertex[n][2] - corner[3][2])) > outer)
								continue;
						}
						else if (vertex[n][2] > domain[5])
						{
							if (sqrt((vertex[n][0] - corner[7][0]) * (vertex[n][0] - corner[7][0]) + (vertex[n][1] - corner[7][1]) * (vertex[n][1] - corner[7][1]) + (vertex[n][2] - corner[7][2]) * (vertex[n][2] - corner[7][2])) > outer)
								continue;
						}
						else if (sqrt((vertex[n][0] - domain[3]) * (vertex[n][0] - domain[3]) + (vertex[n][1] - domain[4]) * (vertex[n][1] - domain[4])) > outer)
							continue;
					}
					else
					{
						if (vertex[n][2] < domain[2])
						{
							if (sqrt((vertex[n][0] - domain[3]) * (vertex[n][0] - domain[3]) + (vertex[n][2] - domain[2]) * (vertex[n][2] - domain[2])) > outer)
								continue;
						}
						else if (vertex[n][2] > domain[5])
						{
							if (sqrt((vertex[n][0] - domain[3]) * (vertex[n][0] - domain[3]) + (vertex[n][2] - domain[5]) * (vertex[n][2] - domain[5])) > outer)
								continue;
						}
						else if (vertex[n][0] - domain[3] > outer)
							continue;
					}
				}
				else
				{
					if (vertex[n][1] < domain[1])
					{
						if (vertex[n][2] < domain[2])
						{
							if (sqrt((vertex[n][1] - domain[1]) * (vertex[n][1] - domain[1]) + (vertex[n][2] - domain[2]) * (vertex[n][2] - domain[2])) > outer)
								continue;
						}
						else if (vertex[n][2] > domain[5])
						{
							if (sqrt((vertex[n][1] - domain[1]) * (vertex[n][1] - domain[1]) + (vertex[n][2] - domain[5]) * (vertex[n][2] - domain[5])) > outer)
								continue;
						}
						else if (domain[1] - vertex[n][1] > outer)
							continue;
					}
					else if (vertex[n][1] > domain[4])
					{
						if (vertex[n][2] < domain[2])
						{
							if (sqrt((vertex[n][1] - domain[4]) * (vertex[n][1] - domain[4]) + (vertex[n][2] - domain[2]) * (vertex[n][2] - domain[2])) > outer)
								continue;
						}
						else if (vertex[n][2] > domain[5])
						{
							if (sqrt((vertex[n][1] - domain[4]) * (vertex[n][1] - domain[4]) + (vertex[n][2] - domain[5]) * (vertex[n][2] - domain[5])) > outer)
								continue;
						}
						else if (vertex[n][1] - domain[4] > outer)
							continue;
					}
					else if (vertex[n][2] - domain[5] > outer || domain[2] - vertex[n][2] > outer)
						continue;
				}

				if (sqrt((corner[0][0] - vertex[n][0]) * (corner[0][0] - vertex[n][0]) + (corner[0][1] - vertex[n][1]) * (corner[0][1] - vertex[n][1]) + (corner[0][2] - vertex[n][2]) * (corner[0][2] - vertex[n][2])) < inner && sqrt((corner[1][0] - vertex[n][0]) * (corner[1][0] - vertex[n][0]) + (corner[1][1] - vertex[n][1]) * (corner[1][1] - vertex[n][1]) + (corner[1][2] - vertex[n][2]) * (corner[1][2] - vertex[n][2])) < inner && sqrt((corner[2][0] - vertex[n][0]) * (corner[2][0] - vertex[n][0]) + (corner[2][1] - vertex[n][1]) * (corner[2][1] - vertex[n][1]) + (corner[2][2] - vertex[n][2]) * (corner[2][2] - vertex[n][2])) < inner && sqrt((corner[3][0] - vertex[n][0]) * (corner[3][0] - vertex[n][0]) + (corner[3][1] - vertex[n][1]) * (corner[3][1] - vertex[n][1]) + (corner[3][2] - vertex[n][2]) * (corner[3][2] - vertex[n][2])) < inner && sqrt((corner[4][0] - vertex[n][0]) * (corner[4][0] - vertex[n][0]) + (corner[4][1] - vertex[n][1]) * (corner[4][1] - vertex[n][1]) + (corner[4][2] - vertex[n][2]) * (corner[4][2] - vertex[n][2])) < inner && sqrt((corner[5][0] - vertex[n][0]) * (corner[5][0] - vertex[n][0]) + (corner[5][1] - vertex[n][1]) * (corner[5][1] - vertex[n][1]) + (corner[5][2] - vertex[n][2]) * (corner[5][2] - vertex[n][2])) < inner && sqrt((corner[6][0] - vertex[n][0]) * (corner[6][0] - vertex[n][0]) + (corner[6][1] - vertex[n][1]) * (corner[6][1] - vertex[n][1]) + (corner[6][2] - vertex[n][2]) * (corner[6][2] - vertex[n][2])) < inner && sqrt((corner[7][0] - vertex[n][0]) * (corner[7][0] - vertex[n][0]) + (corner[7][1] - vertex[n][1]) * (corner[7][1] - vertex[n][1]) + (corner[7][2] - vertex[n][2]) * (corner[7][2] - vertex[n][2])) < inner)
					continue; // domain lies within inner sphere

				rdom = 0.5 * sqrt((domain[3] - domain[0]) * (domain[3] - domain[0]) + (domain[4] - domain[1]) * (domain[4] - domain[1]) + (domain[5] - domain[2]) * (domain[5] - domain[2]));
				dist = sqrt((0.5 * domain[0] + 0.5 * domain[3] - vertex[n][0]) * (0.5 * domain[0] + 0.5 * domain[3] - vertex[n][0]) + (0.5 * domain[1] + 0.5 * domain[4] - vertex[n][1]) * (0.5 * domain[1] + 0.5 * domain[4] - vertex[n][1]) + (0.5 * domain[2] + 0.5 * domain[5] - vertex[n][2]) * (0.5 * domain[2] + 0.5 * domain[5] - vertex[n][2]));

				if (dist <= rdom) // vertex lies within domain enclosing sphere
				{
					n++;
					continue;
				}

				if (((0.5 * domain[0] + 0.5 * domain[3] - vertex[n][0]) * lightcone.direction[0] + (0.5 * domain[1] + 0.5 * domain[4] - vertex[n][1]) * lightcone.direction[1] + (0.5 * domain[2] + 0.5 * domain[5] - vertex[n][2]) * lightcone.direction[2]) / dist >= lightcone.opening) // center of domain lies within opening
				{
					n++;
					continue;
				}

				if (dist > outer && acos(((0.5 * domain[0] + 0.5 * domain[3] - vertex[n][0]) * lightcone.direction[0] + (0.5 * domain[1] + 0.5 * domain[4] - vertex[n][1]) * lightcone.direction[1] + (0.5 * domain[2] + 0.5 * domain[5] - vertex[n][2]) * lightcone.direction[2]) / dist) - acos(lightcone.opening) <= acos((outer * outer + dist * dist - rdom * rdom) / (2. * outer * dist))) // enclosing sphere within opening
				{
					n++;
					continue;
				}

				if (dist <= outer && acos(((0.5 * domain[0] + 0.5 * domain[3] - vertex[n][0]) * lightcone.direction[0] + (0.5 * domain[1] + 0.5 * domain[4] - vertex[n][1]) * lightcone.direction[1] + (0.5 * domain[2] + 0.5 * domain[5] - vertex[n][2]) * lightcone.direction[2]) / dist) - acos(lightcone.opening) <= asin(rdom / dist)) // enclosing sphere within opening
				{
					n++;
				}
			}
		}
	}

	return n;
}

//////////////////////////
// hourMinSec
//////////////////////////
// Description:
//   generates formatted output for cpu-time: hh..h:mm:ss.s
//
// Arguments:
//   seconds    number of seconds
//
// Returns:
//   formatted string
//
//////////////////////////

string hourMinSec(double seconds)
{
	string output;
	char ptr[20];
	int h, m, s, f;

	h = (int)floor(seconds / 3600.);
	seconds -= 3600. * h;
	m = (int)floor(seconds / 60.);
	seconds -= 60. * m;
	s = (int)floor(seconds);
	seconds -= s;
	f = (int)floor(10. * seconds);
	snprintf(ptr,20, "%d:%02d:%02d.%d", h, m, s, f);

	output.reserve(20);
	output.assign(ptr);

	return output;
}

#endif