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

/*
Written by Kenneth E. Moore
July 1, 2003
*/

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 microns

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 an ideal Lambertian scattering surface */

int __declspec(dllexport) APIENTRY UserScatterDefinition(double *data)
	{
	double px, py, pz, qx, qy, qz, nx, ny, nz;
	double A, B, C, MAG, random_number;

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

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

	/* return flag to indicate we scattered */
	data[10] = 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.
	*/

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

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

	/*
	Note the Lambertian distribution does not depend upon s.
	*/

	/* find vectors p,q perpindicular to n */
	/* first pick any vector other than n */
	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 */
	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 n and q are orthonormal already

	/* randomly choose new direction cosines */
	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);

	data[4] = A*nx + B*px + C*qx; 
	data[5] = A*ny + B*py + C*qy; 
	data[6] = A*nz + B*pz + C*qz; 
	/* that's the 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,"Fraction");
	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;
}
