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

/*
Written by Kenneth E. Moore
May 11, 2000

Modified January 24, 2005 to support rectangular periodic arrays
*/

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

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

/* 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
data[2]  = y
data[3]  = z
data[4]  = wavelength in µm
data[5]  = millimeters per unit length (1.0 for mm, 25.4 for inches, 10.0 for cm and 1000.0 for meters)
data[6]  = n (to be computed by the dll and returned)
data[7]  = n * dn/dx (to be computed by the dll and returned)
data[8]  = n * dn/dy (to be computed by the dll and returned)
data[9]  = n * dn/dz (to be computed by the dll and returned)
data[10] = parameter 1 from user input
data[11] = parameter 2 from user input
etc... up to data[maxdata] where maxdata = int(data[0])

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

/* this DLL models a GRIN media similar to the GRIN 1 surface type. */
/* this is n(r) = N + A*r^2 + B*r */

// Modified to be periodic with period dx and dy if both are greater than zero to support arrays

int __declspec(dllexport) APIENTRY UserGrinDefinition(double *data)
	{
   double r, r2, index, N, A, B, x, y, dx, dy;

   N = data[10];
   A = data[11];
   B = data[12];
	dx= data[13];
	dy= data[14];
   x = data[1];
   y = data[2];

	if (dx > 0 && dy > 0)
		{
		// adjust to center of array
		if (x > +0.5*dx) x -= dx*(double)ceil(+(x - 0.5*dx)/dx);
		if (x < -0.5*dx) x += dx*(double)ceil(-(x + 0.5*dx)/dx);
		if (y > +0.5*dy) y -= dy*(double)ceil(+(y - 0.5*dy)/dy);
		if (y < -0.5*dy) y += dy*(double)ceil(-(y + 0.5*dy)/dy);
		}

   r2 = x*x + y*y;
	r  = sqrt(r2);

   index = N + A*r2 + B*r;
	data[6] = index;

   if (index < 1.0)
		{
		return -1;
		}

	/* now the derivatives, note they are multiplied by index */
   data[7] = 2.0*A*x*index;
   data[8] = 2.0*A*y*index;
	if (B != 0.0 && r > 0.0)
		{
		data[7] += B*x*index/r;
		data[8] += B*y*index/r;
		}
   data[9] = 0.0;
	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,"N0");
	if (i == 2) strcpy(data,"Nr2");
	if (i == 3) strcpy(data,"Nr1");
	if (i == 4) strcpy(data,"Array dx");
	if (i == 5) strcpy(data,"Array dy");
	return 0;
	}
