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

/*
This simple program generates a 4 column file
suitable for importing to the grid surface types.
This generates data for a spherical surface on a grid

Kenneth Moore 11-16-94

Some bugs fixed 3-24-87

Added zeroes for unitflag and x,y decenters to header 5-03-99

*/
void GetSag(double x, double y, double *z);

/* these numbers must be odd */
#define NUM_X_PIXELS 25
#define NUM_Y_PIXELS 25

/* define radius and aperture of lens */
#define RADIUS 100
#define DIAMETER 21

void main(void)
{
int i, j, ic, jc;
double delx, dely, x, y, z, dx, dy, dxdy;
double tem1, tem2, tem3, tem4, epsilon;
char disp[120];
FILE *out;

/* STEP 1. First open the file */

out = fopen("PGRID.DAT", "wt");

/* STEP 2. Write four numbers for the header line */
delx = (double) DIAMETER / (NUM_X_PIXELS - 1);
dely = (double) DIAMETER / (NUM_Y_PIXELS - 1);

sprintf(disp, "%i %i %.9E %.9E 0 0 0\n", NUM_X_PIXELS, NUM_Y_PIXELS, delx, dely);
fputs(disp, out);

/* STEP 3. Write the rows and columns */

ic = (NUM_X_PIXELS+1)/2;
jc = (NUM_Y_PIXELS+1)/2;

epsilon = .0001 * min(delx, dely);

for (j = 1; j <= NUM_Y_PIXELS; j++)
	{
	for (i = 1;  i <= NUM_X_PIXELS; i++)
		{
		x =  (i - ic) * delx;
		y = -(j - jc) * dely;

		/* compute the sag */
		GetSag(x, y, &z);

		/* compute dz/dx using finite difference. Could use analytical formula instead */
		GetSag(x+epsilon, y, &tem1);
		GetSag(x-epsilon, y, &tem2);
		dx = (tem1 - tem2)/(2.0*epsilon);

		/* compute dz/dy */
		GetSag(x, y+epsilon, &tem1);
		GetSag(x, y-epsilon, &tem2);
		dy = (tem1 - tem2)/(2.0*epsilon);

		/* compute d2z/dxdy */
		GetSag(x+epsilon, y+epsilon, &tem1);
		GetSag(x-epsilon, y+epsilon, &tem2);
		tem3 = (tem1 - tem2)/(2.0*epsilon);
		GetSag(x+epsilon, y-epsilon, &tem1);
		GetSag(x-epsilon, y-epsilon, &tem2);
		tem4 = (tem1 - tem2)/(2.0*epsilon);
		dxdy = (tem3 - tem4)/(2.0*epsilon);

		sprintf(disp, "%.9E %.9E %.9E %.9E\n", z, dx, dy, dxdy);
		fputs(disp, out);
		}
	}

/* STEP 4. Close the file */
fclose(out);
}

void GetSag(double x, double y, double *z)
{
double c, r;

/* this could be any function, but it happens to be a sphere */

c = 1.0 / RADIUS;
r  = sqrt(x*x+y*y);
*z = c*r*r/(1.0+sqrt(1-(c*c*r*r)));
}
