/** imdemo.cpp
*
*   Create demonstration images
*
***/

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

#define SIZE 128

void image_func(IMAGE *ip, double (*f)(double x, double y), double domain);

double cone(double x, double y);
double cyl(double x, double y);
double gauss(double x, double y);
double exponential(double x, double y);
double airy(double x, double y);
double sinc2(double x, double y);
double sinc2xy(double x, double y);
double tri2(double x, double y);

const double pi = 4.0*atan(1.0);

#define min(x,y)  ((x)<(y)? (x): (y))

struct Choice
{
	double (*f)(double x, double y);
	char *name;
};

Choice choices[] = {
	{cone,"cone(r)"},
	{cyl,"cyl(r)"},
	{gauss,"gauss(r)"},
	{airy,"airy(r)"},
	{exponential,"exponential(r)"},
	{sinc2,"sinc(r)^2"},
	{sinc2xy,"(sinc(x)*sinc(y))^2"},
	{tri2, "tri(x)*tri(y)"}
};

int nchoices = sizeof(choices)/ sizeof(choices[0]);

int main(int argc, char *argv[])
{
	IMAGE *ip;
	char *filename;
	int i, size;
	int choice=0;
	double domain = 4.0;
	if (argc<2) {
		printf("Program usage:\n");
		printf("  imdemo [choice] [filename] [size] [domain]\n");
	}
	if (argc>1) choice = atoi(argv[1]);
	else {
		printf("\nEnter choice:\n");
		for (i=0; i<nchoices; i++)
			printf(" (%2d) %s\n",i+1,choices[i].name);
		scanf(" %d",&choice);
	}
	if (choice>0 && choice<=nchoices) {
		printf("\nGenerating image (%d) %s\n", choice, choices[choice-1].name);
	}
	else {
		printf("choice (%d) not in list\n",choice);
		return -1;
	}

	filename=(argc>2? argv[2]: "demo.fmg");

	size = SIZE;
	if (argc>3) sscanf(argv[3],"%d",&size);
	if ((size<2) || (size>1024)) size = SIZE;

	if (argc>4) sscanf(argv[4]," %lf",&domain);

	ip = make_image(filename,size,size,PIX_FLOAT);
	if (!ip) {
		printf("error creating image file\n");
		return -1;
	}

	printf("image file: %s\n",filename);
	printf("image size: %d x %d\n",size,size);
	printf("image domain: (%g x %g)\n", domain,domain);

	image_func(ip,choices[choice-1].f,domain);

	return 0;
}


void image_func(IMAGE *ip, double (*f)(double x, double y), double domain)
{
	int i, j;
	int mx, my;
	double x, y;
	double dx, dy;
	pixel *buf;
	float *fbuf;

	mx = ip->hlen/2;
	my = ip->vlen/2;
	dx = domain / ip->hlen;
	dy = domain / ip->vlen;

	buf = make_buffer(ip);
	fbuf = (float *) buf;

	for (j=0; j<ip->vlen;  j++) {
		y = dy*(my-j);
		for (i=0; i<ip->hlen; i++) {
			x = dx*(i-mx);
			fbuf[i] = (float) (*f)(x,y);
		}
		put_line(ip,j,buf,PIX_FLOAT);
	}
	free_buffer(buf);
}

/*---------------------------------------------------------------------
*
*    One dimensional functions
*/

double rect(double x)
{
	x = fabs(x);
	return (x>0.5? 0.0: (x<0.5? 1.0: 0.5) );
}

double sinc(double x)
{
	if (fabs(x)<=1e-8) return 1.0;
	x *= pi;
	return sin(x)/x;
}

double tri(double x)
{
	x = fabs(x);
	return (x>1.0? 0.0: 1.0-x);
}

/*-------------------------------------------------------------
*
*   Two dimensional functions
*/

/*
*  Create a conical intensity pattern
*     same as tri(r)
*/

double cone(double x, double y)
{
	double r;
	r = hypot(x,y);
	return  (r<1.0? 1.0-r: 0.0);
}

/*
*   Cylinder function
*      same as rect(r)
*/

double cyl(double x, double y)
{
	double r;
	r = hypot(x,y);
	return (r>0.5? 0.0: (r<0.5? 1.0: 0.5) );
}

double gauss(double x, double y)
{
	double rsq;
	rsq = x*x+y*y;
	return exp(-pi*rsq);
}

double tri2(double x, double y)
{
	return tri(x)*tri(y);
}

/**
*
*  Create a airy(r) intensity pattern
*
**/

double airy(double x, double y)
{
	double r, z;
	r = pi*hypot(x,y);
	z = (r>0.0? 2.0*j1(r)/r:1.0);
	return z*z;
}

double exponential(double x, double y)
{
	double r;
	r = hypot(x,y);
	return exp(-2.0*r);
}


/**
*
*  Create a sinc squared intensity pattern
*
**/

double sinc2(double x, double y)
{
	double r, z;
	r = hypot(x,y);
	z = sinc(r);
	return z*z;
}

double sinc2xy(double x, double y)
{
	double z;
	z = sinc(x) * sinc(y);
	return z*z;
}