#ifndef lint
static char *sccsid = "@(#)alipoly.c	1.1 (UKC) 9/12/87";
#endif  lint
/*
 *	Draw image for hidden line removal stuff to bomb out on.
 *
 *	Y
 *	^
 *	|  Z
 *	| /
 *	|/
 *	+----> X
 *
 *	phi is horizontal scanning, starting on positive X axis and going
 *	towards Z axis.
 *
 *	theta is elevation, horizontal = 0.0, positive towards y-axis.
 *
 *	Vectors (once clipped) are plotted in the order and direction
 *	in which they were supplied here, so you can reasonably try to
 *	optimise head motion by drawing continuously as far as possible.
 */

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

#define RES 5		/* stepsize in degrees, should be a factor of 180 */

#define PI 3.141592653589793

/* convert array index to angle */
#define PtoPHI(p) ((p*RES) * (PI/180))
#define TtoTHETA(t) ((t*RES - 90) * (PI/180))
/* defines to save typing */
#define MAXP (360/RES)
#define MAXT (180/RES)
#define LOWER (lower/RES)
#define UPPER (upper/RES)

static double x[MAXP+1][MAXT+1];
static double y[MAXP+1][MAXT+1];
static double z[MAXP+1][MAXT+1];

double maxp;	/* power in major lobe */
double prel(), maxpwr();
static triangle(), line();
extern char *optarg;
extern optind;
/*
 * Function in power.c
 * Given azimuth-elevation in radians, return
 * power flowing in that direction
 */
extern double power();

genpoly(argc,argv)
	char **argv;
{
	int p, t;		/* phi and theta loop counters */
	double sintilt, costilt;/* precomputed for tilting */
	double sinphi[MAXP+1], cosphi[MAXP+1];	/* precomputed */
	double tilt = 0.0;	/* see option decoding */
	double rotate = PI/180.0;	/* 1 degree frig */
	int lower = 0;
	int upper = 180;
	char c;

	while((c = getopt(argc,argv,"t:r:l:u:")) != EOF)
	{	switch(c)
		{
		case 't':	/* tilt object towards you in degrees */
			tilt = (double)atoi(optarg);
			tilt = tilt * PI/180.0;
			break;
		case 'r':	/* rotate object in degrees */
			rotate = (double)atoi(optarg);
			rotate = rotate + 1.0;	/* avoid nasty black lines */
			rotate = rotate * PI/180.0;
			break;
		case 'l':
			lower = atoi(optarg)+90;
			if(lower < 0)
			{	error("lower elevation limit too low");
				exit(1);
			}
			break;
		case 'u':
			upper = atoi(optarg)+90;
			if(upper > 180)
			{	error("upper elevation limit too high");
				exit(1);
			}
			break;
		case '?':
		default:
			error("unknown option %s", optarg);
			exit(1);
		}
	}
	if(optind != argc)
	{	error("bad arg %s",argv[optind]);
		exit(1);
	}
	if(lower >= upper)
	{	error("elevation bounds inside out");
		exit(1);
	}

	/* precompute sin( and cos(phi) */
	for (p = 0; p < MAXP; p++) {
		double phi;
		phi = PtoPHI(p);
		sinphi[p] = sin(phi);
		cosphi[p] = cos(phi);
	}
	/* wraparound */
	sinphi[MAXP] = sinphi[0]; cosphi[MAXP] = cosphi[0];

	maxp = maxpwr();

	for (t = LOWER; t <= UPPER; t++) {
		double theta, sintheta, costheta;
		theta = TtoTHETA(t);
		sintheta = sin(theta);
		costheta = cos(theta);
		for (p = 0; p < MAXP; p++) {
			double r; /* radius of object in this direction */
			double phi = PtoPHI(p) + rotate;
			/* image-generating function */
			r = prel(phi,theta);
			x[p][t] = r * cosphi[p] * costheta;
			z[p][t] = r * sinphi[p] * costheta;
			y[p][t] = r * sintheta;
		}
		/* wraparound */
		x[p][t] = x[0][t];
		y[p][t] = y[0][t];
		z[p][t] = z[0][t];
	}

	/* tilt-um */
	sintilt = sin(tilt);
	costilt = cos(tilt);
	for (t = 0; t <= MAXT; t++) {
		for (p = 0; p <= MAXP; p++) {
			double newy, newz; /* each depends on the old other */

			newz = z[p][t] * costilt - y[p][t] * sintilt;
			newy = y[p][t] * costilt + z[p][t] * sintilt;
			z[p][t] = newz;
			y[p][t] = newy;
		}
	}

	axes(-rotate,tilt);

	/* generate all obscuring quadrilateral polygons as triangle pairs */
	for (p = 0; p < MAXP; p++) {
		/* bottom triangle */
		triangle(p, 1, p, 0, p+1, 1);

		/* square tiles */
		for (t = 2; t < MAXT; t++) {
			triangle(p, t, p, t-1, p+1, t-1);
			triangle(p, t, p+1, t-1, p+1, t);
		}

		/* top triangle */
		triangle(p, MAXT, p, MAXT-1, p+1, MAXT-1);
	}
	
	/* draw lines of longitude(?) around y-axis */
	for (p = 0; p < MAXP; p += 1+5/RES) {
		for (t = 0; t < MAXT; t++) {
			line(p, t, p, t+1);
		}
	}

	/* lines of latitude, ring-shaped. */
	for (t = 1; t < MAXT; t += 1+5/RES) {
		for (p = 0; p < MAXP; p++) {
			line(p, t, p+1, t);
		}
	}
}

/*
 * Draw some sort of axis representation
 */
static
axes(R,T)
	double R,T;
{
	double x1,y1,z1;
	double x2,y2,z2;

	/*
	 * x-axis
	 */
	viewpoint(R,T, -1.0,0.0,0.0, &x1,&y1,&z1);
	viewpoint(R,T, 1.0,0.0,0.0, &x2,&y2,&z2);
	addline(x1,y1,z1, x2,y2,z2);

	/*
	 * z-axis
	 */
	viewpoint(R,T, 0.0,0.0,-1.0, &x1,&y1,&z1);
	viewpoint(R,T, 0.0,0.0,1.0, &x2,&y2,&z2);
	addline(x1,y1,z1, x2,y2,z2);

	/*
	 * +ve y-axis
	 */
	viewpoint(R,T, 0.0,1.0,0.0, &x2,&y2,&z2);
	addline(0.0,0.0,0.0, x2,y2,z2);
}

/*
 * Transform x,z,y for viewpoint rotation and tilt
 */
viewpoint(R,T, x,y,z, newx,newy,newz)
	double R,T;	/* rotation and tilt in radians */
	double x,y,z;
	double *newx,*newy,*newz;
{
	double z1;

	z1 = z*cos(R) + x*sin(R);
	*newx = x*cos(R) - z*sin(R);

	*newz = z1*cos(T) - y*sin(T);
	*newy = y*cos(T) + z1*sin(T);
}

static
triangle(p1, t1, p2, t2, p3, t3)
{
	double x1,y1,z1;
	double x2,y2,z2;
	double x3,y3,z3;

	x1 = x[p1][t1]; x2 = x[p2][t2]; x3 = x[p3][t3];
	y1 = y[p1][t1]; y2 = y[p2][t2]; y3 = y[p3][t3];
	z1 = z[p1][t1]; z2 = z[p2][t2]; z3 = z[p3][t3];

	if(x1 == x2 && x1 == x3 &&
		y1 == y2 && y1 == y3 &&
		z1 == z2 && z1 == z3)
			return;

	addtri(x1,y1,z1, x2,y2,z2, x3,y3,z3);
}

static
line(p1, t1, p2, t2)
{
	double x1,y1,z1;
	double x2,y2,z2;

	x1 = x[p1][t1]; x2 = x[p2][t2];
	y1 = y[p1][t1]; y2 = y[p2][t2];
	z1 = z[p1][t1]; z2 = z[p2][t2];

	if(x1 == x2 && y1 == y2 && z1 == z2) return;
	addline(x1,y1,z1, x2,y2,z2);
}

/*
 * alan's bits
 */

#define DBRANGE	40.0

/*
 * find the maximum power radiated in any direction, ie
 * the max 'gain' of the aerial. The result is used to
 * normalise subsequent calls to 'power', via 'prel'.
 */
double
maxpwr() 
{	double a,e;
	double p,mp;
	double incr = PI/12.0;

	for(mp = 0.0, a = 0.0; a < 2*PI; a += incr)
		for(e = -PI/2; e < PI/2; e += incr)
			if((p = power(a,e)) > mp)
				mp = p;
	return(mp);
}

/*
 * Return a normalised power in the range 0 to 1.0
 * args are azimuth-elevation (radians).
 */
double
prel(a,e)
	double a,e;
{	double power();
	double p;

	p = power(a,e)/maxp;
	if(p < 0.00001) return(0.0);
	p = (DBRANGE + 10.0*log10(p)) / DBRANGE;
	if(p < 0.0) p = 0.0;
	return(p);
}