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


	/*******************************************************************
	 * Calculate the far-field diffraction pattern of some specified
	 *   aperture.  
	 *
	 * MWR 11/19/2024
	 *
	 * This version employs Fresnel approximations, not Fraunhofer
	 *   approximations, and so should be more accurate in more
	 *   near-ish-field regimes -- such as planned starshields.
	 *
	 * MWR 11/20/2024
	 *
	 * Added a fourth stanza in output, which is (1 - normalized intensity)
	 *   at each point on the screen.
	 *
	 */


#define MYPI 3.14159265359

static int debug = 1;


	/* some physical constants for starshield */

	/* wavelength of light (in meters) */
#define LAMBDA 1e-6
	/* maximum size of starshield (in meters) */
#define SHIELD_SIZE 10.0
	/* distance from shield to observer (in meters) */
#define SHIELD_DISTANCE 1.0e7
	/* size of projection screen (in meters) */
#define SCREEN_SIZE 1.0e0
	/* speed of light (m/s) */
#define LIGHTSPEED  299792458.0


	/* number of elements in half of the aperture mask, and the screen */
#define HALF_ELEMS   50 
#define NUMPIX  (2*HALF_ELEMS + 1)


	/* we break mask and screen into this many pieces */
static double aperture[NUMPIX][NUMPIX];
static double screen[NUMPIX][NUMPIX];

	/* useful global variables  */
	/* wave number of light (rad/meter) */
static double k = (2.0*MYPI) / LAMBDA;	


	/* static functions in this source code file */
static int init_aperture (void);
static int init_screen (void);
static int compute_pattern (void);
static double aperture_pixel_size (void);


int 
main (
	int argc,
	char *argv[]
	)
{

	


	if (init_aperture() != 0) {
		printf("init_aperture returns with error.  Quitting. \n");
		exit(1);
	}

	if (init_screen() != 0) {
		printf("init_screen returns with error.  Quitting. \n");
		exit(1);
	}


	if (compute_pattern() != 0) {
		printf("compute_pattern returns with error.  Quitting \n");
		exit(1);
	}


	return(0);
}



	/************************************************************
	 * PROCEDURE: init_aperture
	 *
	 * DESCRIPTION: Set the elements of the obstructing screen.
	 *              Pixels which are transparent are set to 1,
	 *              pixels which are opaque are set to 0.
	 *
	 *              Actually, allow the ability to set some
	 *              pixels to partially-transparent values,
	 *              between 0 and 1.  This will allow us to 
	 *              simulate apodized apertures.
	 *
	 *              Add a version which makes petal-based apodization.
	 *
	 * RETURNS:  
	 *                 0             if all goes well
	 *                 1             if a problem occurs
	 */

static int
init_aperture
	(
	void
	)
{
	static int i, j;


#if 0
	/* start by setting all pixels to 1 = transparent */
	for (i = 0; i < NUMPIX; i++) {
		for (j = 0; j < NUMPIX; j++) {
			aperture[i][j] = 1;
		}
	}

	if (debug > 0) {
		int i, j;

		printf("# here comes aperture \n");
		for (i = 0; i < NUMPIX; i++) {
			for (j = 0; j < NUMPIX; j++) {
				printf(" %4.2f", aperture[i][j]);
			}
			printf("\n");
		}
		printf("\n");
		printf("\n");
	}

	return(0);
#endif


	/* start by setting all pixels to 0 = opaque */
	for (i = 0; i < NUMPIX; i++) {
		for (j = 0; j < NUMPIX; j++) {
			aperture[i][j] = 0;
		}
	}

	/* now set chosen pixels to 1 to mark transparent sections of mask */
#if 0
	aperture[HALF_ELEMS][HALF_ELEMS] = 1;
	aperture[HALF_ELEMS][HALF_ELEMS+12] = 1;
#endif

#if 0
	/*
	 * This section makes a circular aperture of some given physical size.
	 * No sophistication around edges: each pixel is either all-in,
	 * or all-out.
	 */
	double hole_radius = 5.0;		/* radius of hole, in meters */
	double hole_center_x = 0.0;	/* center of hole, rel to center of screen */
	double hole_center_y = 0.0;	/* center of hole, rel to center of screen */
	double pix_size = aperture_pixel_size();
	int screen_center_x = HALF_ELEMS;	/* center of screen, in pixel units */
	int screen_center_y = HALF_ELEMS;	/* center of screen, in pixel units */

	double hole_radius_sq = hole_radius*hole_radius;
	double x, y, dx, dy, dist_sq;

#undef APOD_C

	for (i = 0; i < NUMPIX; i++) {
		for (j = 0; j < NUMPIX; j++) {
			x = (i - screen_center_x)*pix_size;
			y = (j - screen_center_y)*pix_size;
			dx = (x - hole_center_x);
			dy = (y - hole_center_y);
			dist_sq = dx*dx + dy*dy;

#ifdef OPAQUE
			if (dist_sq <= hole_radius_sq) {
				aperture[i][j] = 1.0;
			}	
#endif
#ifdef APOD_B
				/* here we apodize from 4.0 -> 5.0 in linear manner */
			if (dist_sq <= (hole_radius*0.80)*(hole_radius*0.80)) {
				aperture[i][j] = 1.0;
			} 
			else if (dist_sq <= hole_radius_sq) {
				aperture[i][j] = 1.0 - (dist_sq / hole_radius_sq);
			}
#endif
#ifdef APOD_C
				/* here we apodize from 3.0 -> 5.0 in quadratic manner */
			if (dist_sq <= (hole_radius*0.60)*(hole_radius*0.60)) {
				aperture[i][j] = 1.0;
			} 
			else if (dist_sq <= hole_radius_sq) {
				aperture[i][j] = 1.0 - ((dist_sq / hole_radius_sq)*(dist_sq / hole_radius_sq));
			}
#endif
#endif

#if 1
			/* 
		 	 * This section makes an apodized mask with petal design.
			 *   See Hu's Ph.D. dissertation, eqn. 2.17, for explanation.
		 	 */
	double hole_radius = 5.0;		/* radius of hole, in meters */
	double hole_center_x = 0.0;	/* center of hole, rel to center of screen */
	double hole_center_y = 0.0;	/* center of hole, rel to center of screen */
	double pix_size = aperture_pixel_size();
	int screen_center_x = HALF_ELEMS;	/* center of screen, in pixel units */
	int screen_center_y = HALF_ELEMS;	/* center of screen, in pixel units */

	double hole_radius_sq = hole_radius*hole_radius;
	double x, y, dx, dy, dist_sq;

	int num_petals = 12;
	

	for (i = 0; i < NUMPIX; i++) {
		for (j = 0; j < NUMPIX; j++) {
			x = (i - screen_center_x)*pix_size;
			y = (j - screen_center_y)*pix_size;
			dx = (x - hole_center_x);
			dy = (y - hole_center_y);
			dist_sq = dx*dx + dy*dy;

			if (dist_sq < hole_radius_sq) {
	
				/* 
				 * this is the section where we check to see if pixel 
				 *   is inside a petal
				 */
				double r = sqrt(dist_sq);
				double theta = atan2(dy, dx);
				if (debug > 1) {
					printf("# i %4d j %4d   dx %5.2f  dy %5.2f  r %5.2f theta %5.1f \n",
							i, j, dx, dy, r, theta);
				}

				/* compute A(r) value */
				double A = 0.0;
				if (r < 0.7*hole_radius) {
					A = 1.0;
				}
				else {
					A = 2.0 - ((2.0*r)/hole_radius);
				}
				if (debug > 1) {
					printf("#   A  %6.3f \n", A);
				}

				/* 
				 * loop over all possible petals, see if this (r,theta)
				 *   combination falls inside a petal.
				 */
				int n;
				int inside_petal = 0;
	
				for (n = 0; n < num_petals; n++) {
					double min = (2.0*MYPI*n - MYPI*A)/((double) num_petals);
					double max = (2.0*MYPI*n + MYPI*A)/((double) num_petals);
					if ((theta >= min) && (theta <= max)) {
						inside_petal = 1;
					}
					if (((theta + MYPI) >= min) && ((theta + MYPI) <= max)) {
						inside_petal = 1;
					}
				}

				if (inside_petal == 1) {
					aperture[i][j] = 1.0;
				}


			}


#endif


			

		}
	}



	if (debug > 0) {
		int i, j;

		printf("# here comes aperture \n");
		for (i = 0; i < NUMPIX; i++) {
			for (j = 0; j < NUMPIX; j++) {
				printf(" %4.2f", aperture[i][j]);
			}
			printf("\n");
		}
		printf("\n");
		printf("\n");
	}


	return(0);
}
	



	/************************************************************
	 * PROCEDURE: init_screen
	 *
	 * DESCRIPTION: Set the elements of the projection screen
	 *              to 0.  As we compute the contributions from 
	 *              light passing through the screen, we'll 
	 *              add the contributions to this array, accumulating
	 *              the eventual image.
	 *
	 * RETURNS:  
	 *                 0             if all goes well
	 *                 1             if a problem occurs
	 */

static int
init_screen
	(
	void
	)
{
	static int i, j;

	/* start by setting all pixels to 0 = black */
	for (i = 0; i < NUMPIX; i++) {
		for (j = 0; j < NUMPIX; j++) {
			screen[i][j] = 0;
		}
	}

	if (debug > 0) {
		int i, j;

		printf("# here comes screen \n");
		for (i = 0; i < NUMPIX; i++) {
			for (j = 0; j < NUMPIX; j++) {
				printf(" %1.0f", screen[i][j]);
			}
			printf("\n");
		}
		printf("\n");
		printf("\n");
	}


	return(0);
}
	



	/************************************************************
	 * PROCEDURE: compute_pattern
	 *
	 * DESCRIPTION: Perform the calculation of the intensity
	 *              measured on the screen, for every location on the
	 *              screen.  Accumulate results in the 
	 *              "screen[][]" array.
	 *
	 * RETURNS:  
	 *                 0             if all goes well
	 *                 1             if a problem occurs
	 */

static int
compute_pattern
	(
	void
	)
{
	int i_x, i_y, i_u, i_v;
	double x, y, u, v;
	double delta_x, delta_y, delta_u, delta_v;
	double theta_x, theta_y;
	
	
	/* for convenience in calculations below */
	double R = SHIELD_DISTANCE;

	if (debug > 0) {
		printf("#  R is %10.4e  k is %10.4e \n", R, k);
	}

	/* 
	 * now compute the stepsize in u (or v), in meters, in aperture array 
	 */
	delta_u = aperture_pixel_size();
	delta_v = delta_u;
	double aperture_max = delta_u * HALF_ELEMS;

	/* 
	 * we will consider all (x,y) coords over a finite range,
	 * which is the linear size over which we expect interesting
	 * amounts of light to fall on the projection screen.
	 * 
	 * The exact value is somewhat arbitrary, but we need to 
	 * get the order of magnitude right.
	 */
	double screen_max = 5.0*(SCREEN_SIZE);
	/* 
	 * now compute the stepsize in x (or y), in radian, in screen array 
	 */
	delta_x = (screen_max) / (double) HALF_ELEMS;
	delta_y = delta_x;

	if (debug > 0) {
		printf("#  x will run %8.1f to %8.1f  steps %10.3f \n",
					0.0 - screen_max, screen_max, delta_x);
		printf("#  u will run %8.1f to %8.1f  steps %10.3f \n",
					0.0 - aperture_max, aperture_max, delta_u);
	}
	
	
	x = (0.0 - screen_max);
	for (i_x = 0; i_x < NUMPIX; i_x++) {

		theta_x = x/R;

		y = (0.0 - screen_max);
		for (i_y = 0; i_y < NUMPIX; i_y++) {

			theta_y = y/R;

			/* these will accumulate contributions from waves in aperture */
			double cos_sum = 0.0;
			double sin_sum = 0.0;

			u = (0.0 - aperture_max);
			for (i_u = 0; i_u < NUMPIX; i_u++) {

				v = (0.0 - aperture_max);
				for (i_v = 0; i_v < NUMPIX; i_v++) {

					if (debug > 1) {
						printf(" x %10.4e  y %10.4e  u %10.4e  v %10.4e \n",
								x, y, u, v);
					}


					/* compute the contribution of this (u,v) to (x,y) */
					/*   this could be optimized quite a bit */
#if 0
					double arg = (k*theta_x*u) + (k*theta_y*v);
#else
					/* this is to shut up compiler while we test code */
					theta_x = theta_x + 1;
					theta_y = theta_y + 1;

					/* this is the actual calculation we need */
					double arg = (k/(2.0*R)) * ( (2*x*u - u*u) + (2*y*v - v*v) );
#endif
					double cos_term = aperture[i_u][i_v] * cos(arg);
					double sin_term = aperture[i_u][i_v] * sin(arg);

					cos_sum += cos_term*delta_u*delta_v;
					sin_sum += sin_term*delta_u*delta_v;

					if (debug > 1) {
						printf("     arg %10.4e  cos_term %10.4e  sin_term %10.4e  cos_sum %10.4e  sin_sum %10.4e \n",
								arg, cos_term, sin_term, cos_sum, sin_sum);
					}

					v += delta_v;
				}

				u += delta_u;
			}

			/* 
			 * okay, we've computed all the waves which strike (x, y)
			 * on the screen.  Now we can compute the intensity here
			 */
			double psi = (cos_sum + sin_sum);
			screen[i_x][i_y] = psi*psi;


			y += delta_y;
		}

		x += delta_x;
	}


	/* find the maximum value in the pattern -- could be useful */
	double max_intensity = 0.0;
	for (i_x = 0; i_x < NUMPIX; i_x++) {
		for (i_y = 0; i_y < NUMPIX; i_y++) {
			if (screen[i_x][i_y] > max_intensity) {
				max_intensity = screen[i_x][i_y];
			}
		}
	}
	if (debug > 0) {
		printf("# max screen intensity is %10.4e \n", max_intensity);
	}




	if (debug > 0) {
		int i, j;

		printf("# here comes screen \n");
		for (i = 0; i < NUMPIX; i++) {
			for (j = 0; j < NUMPIX; j++) {
				printf(" %8.3e", screen[i][j]);
			}
			printf("\n");
		}
		printf("\n");
	}


	/* this ia an attempt to show inverse of intensity */
	if (debug > 0) {
		int i, j;

		printf("\n\n");
		printf("# here comes inverse, sort of \n");
		for (i = 0; i < NUMPIX; i++) {
			for (j = 0; j < NUMPIX; j++) {
				printf(" %8.3e", (1.0 - (screen[i][j] / max_intensity)));
			}
			printf("\n");
		}
		printf("\n");
	}



	return(0);
}
	


	/************************************************************
	 * PROCEDURE: aperture_pixel_size
	 *
	 * DESCRIPTION: Compute the physical size (in meters)
	 *              of one pixel in the aperture array; in other
	 *              words, the physical size (in meters) of one
	 *              pixel in the mask.
	 *
	 *              This routine effectively sets the PHYSICAL SIZE
	 *              (in meters) of the aperture mask array.
	 *
	 * RETURNS:     the pixel size (in meters) 
	 */

static double
aperture_pixel_size
	(
	void
	)
{
	/* 
	 * we will consider all (u,v) coords over a finite range,
	 * even through we're supposed to integrate to infinity.
	 * This is the maximum distance, in meters, from the center
	 * of the aperture screen for which we perform integration.
	 * 
	 * The exact value is somewhat arbitrary, but we need to 
	 * get the order of magnitude right.
	 */
	double aperture_max = 3.0*(SHIELD_SIZE);
	
	/* 
	 * now compute the stepsize in u (or v), in meters, in aperture array 
	 *    In other words, the linear size (in meters) of each pixel.
	 */
	double pixel_size = (aperture_max) / (double) HALF_ELEMS;

	return(pixel_size);
}