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

/*
Written by Kenneth E. Moore
Complete re-write August 11, 2004 to fix the model for r
*/

int __declspec(dllexport) APIENTRY UserSourceDefinition(double *data);
int __declspec(dllexport) APIENTRY UserParamNames(char *data);

#define PI 3.14159265358979323846

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

/*
For the call to UserSourceDefinition, the data is stored as follows:

ndata = the number of elements in the array data
data[0]  = the total number of values in the passed data array
data[1]  = x (to be computed by the dll and returned)
data[2]  = y (to be computed by the dll and returned)
data[3]  = z (to be computed by the dll and returned)
data[4]  = l (to be computed by the dll and returned)
data[5]  = m (to be computed by the dll and returned)
data[6]  = n (to be computed by the dll and returned)
data[7]  = relative intensity (to be computed by the dll and returned)

data[20] = wavelength in µm
data[21] = millimeters per unit length (1.0 for mm, 25.4 for inches, 10.0 for cm and 1000.0 for meters)
data[22] = a random number seed
data[30] = parameter 1 from user input
data[31] = parameter 2 from user input
etc... up to data[maxdata] where maxdata = int(data[0])

The DLL must compute data 1 through 7 given the other data,
and return 0 if it works; else return -1
*/

/*

This DLL models a flat disk source which has an na which
depends upon the radial aperture.

The ray starting point is chosen from a distribution of the form

I(r) = A + B r^2 + C r^4

where 0 <= r <= R.

Once the ray starting position is determined, the ray angles fall
within a cone of angle given by

NA(r) = D + E r^2 + F r^4

within this cone, the ray distribution is uniform.

*/

int __declspec(dllexport) APIENTRY UserSourceDefinition(double *data)
	{
   double R, A, B, C, D, E, F, P;
	double R2, R4;
	double r2, r4;
	double x, y, z, l, m, n, lmr, r, p;
	double costhetam, random_number, random_angle;
	double na, theta_m, del_r, G, Gn;

	srand((DWORD)data[22]);

   R = data[30];
	A = data[31];
   B = data[32];
   C = data[33];
   D = data[34];
   E = data[35];
   F = data[36];

	/* integrate the distribution over a circle */
	R2 = R*R;
	R4 = R2*R2;
	P = (A/2.0 + B*R2/4.0 + C*R4/6.0)*R2;

	/* normalize to 1 */
	if (P <= 0.0)
		{
		r = 0.0;
		goto jump_here_if_r_is_zero;
		}
	else
		{
		A /= P;
		B /= P;
		C /= P;
		}

	/* pick a random number between 0 and 1 */
	p = (double)rand()/(double)RAND_MAX;

	/* find r that yields this p. We can't solve this directly, so use Newton's method */
	r = R*p; // first guess
	del_r = 1.0;
	while(fabs(del_r)>(1e-12))
		{
		r2 = r*r;
		r4 = r2*r2;
		G = (A/2.0 + B*r2/4.0 + C*r4/6.0)*r2 - p;
		Gn = (A + B*r2 + C*r4)*r;
		if (Gn == 0.0) del_r = 0.0;
		else del_r = -G/Gn;
		r += del_r;
		}

	jump_here_if_r_is_zero:;

	/* okay, we have r... now pick an angle */
	random_angle = 2.0*PI * (double)rand()/(double)RAND_MAX;
	x = r*cos(random_angle);
	y = r*sin(random_angle);
	z = 0.0;

	/* at this aperture, what is the NA? */
	r2 = r*r;
	r4 = r2*r2;
	na = D + E*r2 + F*r4;
	if (na < 0.0) na = 0.0;
	if (na > 1.0) na = 1.0;

	/* now choose the ray cosines */
	theta_m = asin(na);
	costhetam = cos(theta_m);

	random_number = (double)rand()/(double)RAND_MAX;
	n = 1.0 - random_number*(1.0 - costhetam);

	random_number = (double)rand()/(double)RAND_MAX;
	random_angle = 2.0*PI*random_number;

	lmr = sqrt(1.0 - n*n);
	l = lmr * cos(random_angle);
	m = lmr * sin(random_angle);

	data[1] = x;
	data[2] = y;
	data[3] = z;
	data[4] = l;
	data[5] = m;
	data[6] = n;
	data[7] = 1.0;
	return 0;
   }

/*

For the call to UserSourceDefinition, the data is stored as follows:
data[0] is an integer indicating the number of the parameter whose name should be returned by the function.
  
*/
int __declspec(dllexport) APIENTRY UserParamNames(char *data)
	{
	/* this function returns the name of the parameter requested */
	unsigned int i;
	i = (unsigned int) (unsigned char) data[0];
	strcpy(data,"");
	if (i == 1) strcpy(data,"R");
	if (i == 2) strcpy(data,"A");
	if (i == 3) strcpy(data,"B");
	if (i == 4) strcpy(data,"C");
	if (i == 5) strcpy(data,"D");
	if (i == 6) strcpy(data,"E");
	if (i == 7) strcpy(data,"F");
	return 0;
	}
