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

/*
Written by Kenneth E. Moore
October 13, 2003
*/

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 a bulk scatter distribution function where the probability of scattering an angle theta from the present path is given by:

p(theta) = c0 + c1*theta + c2*theta^2 + c3*theta^3 + c4*theta^4

The integral of this function from 0 to some theta_max is normalized to 1. Theta is in degrees.

The other angle of rotaion about the specular axis is random, uniform between 0 and 2 pi

*/

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

	/* return 100% transmission if we do not scatter */
	data[12] = 1.0;

	/* data from dialog box */
	mean = data[14]; // mean free path
	angl = data[15]; // maximum scatter angle, must be between 0 and 180
	tran = data[51]; // transmission, must be between 0 and 1
	c0   = data[52]; // coefficient on theta^0
	c1   = data[53]; // coefficient on theta^1
	c2   = data[54]; // coefficient on theta^2
	c3   = data[55]; // coefficient on theta^3
	c4   = data[56]; // coefficient on theta^4

	/* data from the ray trace */
	path = data[9]; //  the unscattered path length

	if (path <= 0.0) return 0; // we will not scatter
	if (mean <= 0.0) return 0; // we will not scatter
	if (angl <= 0.0 || angl > 180.0) return 0; // we will not scatter
	if (tran <= 0.0 || tran > 1.0) return 0; // we will not scatter

	/* randomly seed the random number generator */
	srand((unsigned int) data[16]);

	/* compute the integral of the polynomial, need this for normalization */
	theta = angl;
	F  = theta*(c0 + 0.5*theta*(c1 + (1.0/3.0)*theta*(c2 + 0.25*theta*(c3 + 0.20*theta*(c4)))));
	if (F <= 0.0) return 0; // we will not scatter
	
	/* normalize */
	c0 /= F;
	c1 /= F;
	c2 /= F;
	c3 /= F;
	c4 /= F;

	/* 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

	/* to get here means we will not make it to path, but instead will scatter at t */

	/* find the scatter angle theta */
	random_number = (double)rand()/(double)RAND_MAX;

	/* solve for theta using Newton's method */
	delta = 1.0;
	theta =  0.5*angl;
	while (fabs(delta) > 1E-6)
		{
		F  = theta*(c0 + 0.5*theta*(c1 + (1.0/3.0)*theta*(c2 + 0.25*theta*(c3 + 0.20*theta*(c4))))) - random_number;
		Fp = c0 + theta*(c1 + theta*(c2 + theta*(c3 + theta*(c4))));
		if (Fp == 0.0) return 0; // we will not scatter
		delta = -F/Fp;
		theta += delta;
		if (theta < 0 || theta > angl) return 0; // we will not scatter
		}
	theta = PI*theta/180.0; // convert to radians

	/* 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;

	/* return scatter distance, did scatter flag, and transmission set by user */
	data[9 ] = t;
	data[10] = 1;
	data[12] = tran;
	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");
	if (i == 2) strcpy(data,"C0");
	if (i == 3) strcpy(data,"C1");
	if (i == 4) strcpy(data,"C2");
	if (i == 5) strcpy(data,"C3");
	if (i == 6) strcpy(data,"C4");
	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;
}
