#include <windows.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

/*
Written by Kenneth E. Moore
Complete re-write November 24, 2004 to use b and x scatter model as used by ZEMAX. 
*/

int __declspec(dllexport) APIENTRY UserScatterDefinition(double *data);
int __declspec(dllexport) APIENTRY UserParamNames(char *data);
void CrossProduct(double x1, double y1, double z1, double x2, double y2, double z2, double *x3, double *y3, double *z3);
void Normalize(double *x, double *y, double *z);

BOOL WINAPI DllMain (HANDLE hInst, ULONG ul_reason_for_call, LPVOID lpReserved)
	{
   return TRUE;
   }

/* the data is stored as follows:
data[ 0]  = the total number of values in the passed data array
data[ 1]  = x position of specular ray
data[ 2]  = y position of specular ray
data[ 3]  = z position of specular ray
data[ 4]  = x cosine of specular ray, on output it is the scattered ray
data[ 5]  = y cosine of specular ray, on output it is the scattered ray
data[ 6]  = z cosine of specular ray, on output it is the scattered ray
data[ 7]  = x normal
data[ 8]  = y normal
data[ 9]  = z normal

data[10] = 0 initially
if the DLL scatters the ray return 1 in data[10].
if the DLL returns full polarization data return 2 in data[10].

data[11] = millimeters per unit length (1.0 for mm, 25.4 for inches, 10.0 for cm and 1000.0 for meters)
data[12] = relative energy (to be computed by the dll and returned)
data[16] = a random value to use as a seed
data[17] = wavelength in µm

data[20] = incident Ex real
data[21] = incident Ex imaginary
data[22] = incident Ey real
data[23] = incident Ey imaginary
data[24] = incident Ez real
data[25] = incident Ez imaginary

data 40-45 need to be computed if the DLL sets data[10] = 2
data[40] = output Ex real
data[41] = output Ex imaginary
data[42] = output Ey real
data[43] = output Ey imaginary
data[44] = output Ez real
data[45] = output Ez imaginary

data[51] = input parameter 1 from user
data[52] = input parameter 2 from user
etc... up to data[maxdata] where maxdata = int(data[0])

Return 0 if it works; else return -1.

*/

/* this DLL models a surface that is part specular, part Gaussian distribution 1, and part Gaussian distribution 2 */

int __declspec(dllexport) APIENTRY UserScatterDefinition(double *data)
	{
	int use_lamb;
	double px, py, pz, qx, qy, qz, nx, ny, nz, sx, sy, sz;
	double MAG, random_number;
	double frac1, frac2, fraclamb, sigma, fac;

	/* we need some way of randomizing the random number generator */
	srand((unsigned int) data[16]);

	/* return 100% transmission */
	data[12] = 1.0;

	/*
	The main thing is to make sure the scattered ray lies within a
	hemisphere centered on the normal vector. The angle between the
	normal and the scattered ray must be less than 90 degrees.
	We check this before returning, see below.
	*/

	sx = data[4]; // specular ray
	sy = data[5];
	sz = data[6];

	nx = data[7]; // normal vector
	ny = data[8];
	nz = data[9];

	frac1     = data[51];
	frac2     = data[53];
	fraclamb  = data[55];

	/* grab a random number to see if we scattered using distribution 1, 2, or none */
	/* note this code does not do any error checking... frac1 + frac2 should be no more than 1.0. */

	random_number = (double)rand()/(double)RAND_MAX;
	use_lamb = 0;

	if (random_number < frac1)
		{
		/* use distrib 1 */
		sigma = data[52];
		}
	else
		{
		if (random_number < frac1 + frac2)
			{
			/* use distrib 2 */
			sigma = data[54];
			}
		else
			{
			if (random_number < frac1 + frac2 + fraclamb)
				{
				/* use Lambertian */
				use_lamb = 1;
				}
			else
				{
				/* specular! */
				return 0;
				}
			}
		}

	/* return flag to indicate we scattered */
	data[10] = 1.0;

	/* find vectors p,q perpindicular to n, such that the specular ray lies along the p-n plane */
	px = sx;
	py = sy;
	pz = sz;

	MAG = px*nx + py*ny + pz*nz;

	if (MAG > 0.9999999)
		{
		/* specular is very close to the normal, any p will do */
		if (fabs(nz) < 0.9)
			{
			px = 0.0;
			py = 0.0;
			pz = 1.0;
			}
		else
			{
			px = 1.0;
			py = 0.0;
			pz = 0.0;
			}
		}

	/* this creates q normal to n and p */
	CrossProduct(nx, ny, nz, px, py, pz, &qx, &qy, &qz);
	Normalize(&qx, &qy, &qz);
	/* this creates p normal to both q and n */
	CrossProduct(nx, ny, nz, qx, qy, qz, &px, &py, &pz);
	//Normalize(&px, &py, &pz); // not needed since s and q are orthonormal already

	get_another_scattered_ray:;

	if (use_lamb)
		{
		double A, B, C;

		/* for Lambertian, the scattering is independent of the specular ray */

		MAG = 2.0;
		while (MAG >= 1.0)
			{
			random_number = (double)rand()/(double)RAND_MAX;
			B = 2.0*random_number - 1.0;
			random_number = (double)rand()/(double)RAND_MAX;
			C = 2.0*random_number - 1.0;
			MAG = B*B+C*C;
			}

		A = sqrt(1.0 - MAG);

		/* this is the scattered ray */
		data[4] = A*nx + B*px + C*qx; 
		data[5] = A*ny + B*py + C*qy; 
		data[6] = A*nz + B*pz + C*qz; 
		}
	else
		{
		double tx, ty, tz, bx, by, bo;
		/* use Gaussian distribution for x. For a description of this method see Numerical Recipes in C, Cambridge Press */

		/* the specular ray, when projected onto the surface, lies along p */
		bo = px*sx + py*sy + pz*sz;
		
		get_new_b:;

		MAG = 2.0;
		while (MAG >= 1.0)
			{
			random_number = (double)rand()/(double)RAND_MAX;
			bx = 2.0*random_number - 1.0;
			random_number = (double)rand()/(double)RAND_MAX;
			by = 2.0*random_number - 1.0;
			MAG = bx*bx+by*by;
			}

		fac = sqrt(-log(MAG)/MAG);
		bx *= fac;
		by *= fac;
		bx *= sigma;
		by *= sigma;
		/* t is the scattered ray vector as projected on to the surface */
		tx = bo*px + bx*px + by*qx;
		ty = bo*py + bx*py + by*qy;
		tz = bo*pz + bx*pz + by*qz;
		MAG = tx*tx+ty*ty+tz*tz;
		if (MAG > 1.0) goto get_new_b;

		MAG = sqrt(1.0 - MAG);
		data[4] = tx+MAG*nx;
		data[5] = ty+MAG*ny;
		data[6] = tz+MAG*nz;
		/* that's the scattered ray! */
		}
	
	/*
	NOW, make sure we didn't scatter in such a way that the ray is now more than 90 degrees from normal!
	*/

	if (nx*data[4] + ny*data[5] + nz*data[6] < 0.0)
		{
		goto get_another_scattered_ray;
		}

	return 0;
   }

int __declspec(dllexport) APIENTRY UserParamNames(char *data)
	{
	/* this function returns the name of the parameter requested */
	int i;
	i = (int) data[0];
	strcpy(data,"");
	if (i == 1) strcpy(data,"Frac Sig 1");
	if (i == 2) strcpy(data,"Sigma 1");
	if (i == 3) strcpy(data,"Frac Sig 2");
	if (i == 4) strcpy(data,"Sigma 2");
	if (i == 5) strcpy(data,"Frac Lamb");
	return 0;
	}

void CrossProduct(double x1, double y1, double z1, double x2, double y2, double z2, double *x3, double *y3, double *z3)
{
*x3 = y1*z2 - z1*y2;
*y3 = z1*x2 - x1*z2;
*z3 = x1*y2 - y1*x2;
}

void Normalize(double *x, double *y, double *z)
{
double temp;
temp = (*x)*(*x) +(*y)*(*y)+(*z)*(*z);
temp = sqrt(temp);
if (temp == 0) return;
temp = 1.0/temp;
*x *= temp;
*y *= temp;
*z *= temp;
}
