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

/*
Written by Kenneth E. Moore
August 28, 2000
Modified October 13, 2003 to use more realistic scatter distribution.
*/

int __declspec(dllexport) APIENTRY UserBulkDefinition(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);

#define PI 3.14159265358979323846

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
data[2]  = y position
data[3]  = z position
data[4]  = x cosine
data[5]  = y cosine
data[6]  = z cosine
data[7] = wavelength in microns
data[8] = index of refraction

data[9] = Unscattered path length
if data[9] is longer than the randomly generated scatter position, then the DLL should scatter the ray.
if data[9] is shorter than the randomly generated scatter position, then no scattering should occur!

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[13] = output phase added to ray
data[14] = mean path value from dialog box
data[15] = angle value from dialog box
data[16] = a random value to use as a seed
data[17] = the glue distance... do not scatter a distance shorter than this

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 35-37 need to be computed if the DLL sets data[10] = 1 or 2
data[35] = output x cosine
data[36] = output y cosine
data[37] = output z cosine

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 imaginary scattering material */

int __declspec(dllexport) APIENTRY UserBulkDefinition(double *data)
	{
	double mean, path, angl, random_number, t, theta, phi;
	double sx, sy, sz, nx, ny, nz, tx, ty, tz, rx, ry, rz;
	double cosphi, sinphi, costheta, sintheta;

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

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

	path = data[9];  // the unscattered path
	mean = data[14]; // the mean path
	angl = data[15]; // the maximum scatter angle

	if (mean <= 0.0) return 0; // we will not scatter

	/* Decide if the ray scatters at all. Use exponential probability. */
	get_new_random_number:;
	random_number = (double)rand()/(double)RAND_MAX;
	t = -mean*log(1-random_number);
	if (t <= data[17]) goto get_new_random_number; // the scatter path must be longer than the glue distance!

	/* here is the scatter/no scatter decision point */
	if (t > path) return 0; // we will not scatter

	/* ok, the ray will scatter at distance t. */
	data[9] = t;
	data[10] = 1;

	/* pick theta. */
	random_number = (double)rand()/(double)RAND_MAX;
	theta = angl*PI/180.0*random_number;

	/* we know have theta, measured from the specular. Need phi. */
	random_number = (double)rand()/(double)RAND_MAX;
	phi = 2.0*PI*random_number;

	/* now construct a ray that makes an angle of theta with respect to the specular cosines */
	rx = data[4]; // specular cosines
	ry = data[5];
	rz = data[6];

	/* find a vector n that is normal to r */
	if (fabs(rz) > 0.9)
		{
		sx = 1.0;
		sy = 0.0;
		sz = 0.0;
		}
	else
		{
		sx = 0.0;
		sy = 0.0;
		sz = 1.0;
		}

	CrossProduct(sx, sy, sz, rx, ry, rz, &nx, &ny, &nz);
	Normalize(&nx, &ny, &nz);
	/* find a vector t that is normal to r and n */
	CrossProduct(nx, ny, nz, rx, ry, rz, &tx, &ty, &tz);
	Normalize(&tx, &ty, &tz); // not really required...

	/* find a random vector s that lies in the plane of n and t */
	cosphi = cos(phi);
	sinphi = sin(phi);
	sx = cosphi*nx + sinphi*tx;
	sy = cosphi*ny + sinphi*ty;
	sz = cosphi*nz + sinphi*tz;

	/* the scattered vector is cos(theta)*r + sin(theta)*s */
	costheta = cos(theta);
	sintheta = sin(theta);
	sx = costheta*rx + sintheta*sx;
	sy = costheta*ry + sintheta*sy;
	sz = costheta*rz + sintheta*sz;
	
	/* scattered ray cosines to return to ZEMAX */
	data[35] = sx;
	data[36] = sy;
	data[37] = sz;

	/* attenuate the ray here if desired */
	data[12] = data[51];
	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,"Transmission");
	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;
}
