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


	/*******************************************************************
	 * Calculate the far-field diffraction pattern of some specified
	 *   aperture.  
	 *
	 * MWR 11/19/2024
	 */


#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.
	 *
	 * RETURNS:  
	 *                 0             if all goes well
	 *                 1             if a problem occurs
	 */

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

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

	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) {
				aperture[i][j] = 1.0;
			}	
		}
	}



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

		printf("# here comes aperture \n");
		for (i = 0; i < NUMPIX; i++) {
			for (j = 0; j < NUMPIX; j++) {
				printf(" %1.0f", 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 */
					double arg = (k*theta_x*u) + (k*theta_y*v);
					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");
	}


	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);
}