/*--------------------------------------------------------------------------*/
/* fresnel.c - some simple Fresnel zone calculations
/*--------------------------------------------------------------------------*/
/*
Created 7 May by ecsd using formulas from section 3.7 (pp 90-99) of
"Wireless Communications, Principles and Practice" by Theodore S. Rappaport,
Prentice Hall.

rev 1.1 did fresnel zone radii only
rev 1.2 did single obstructions

$Header: /u0/ecsd/src/RCS/fresnel.c,v 1.4 1998/05/11 08:41:41 ecsd Exp ecsd $
$Log: fresnel.c,v $
Revision 1.4  1998/05/11 08:41:41  ecsd
first steps toward total path analysis, add frequency override (for
testing.) has calculation bugs somewhere.

Revision 1.3  1998/05/08 05:01:21  ecsd
added comments

*/
/*--------------------------------------------------------------------------*/

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

FILE *LOGfp;

#define PI 3.141592653589793	/* yup */
#define RAD 57.29577951			/* degrees per radian */

const double c = 2.99792458e8;	/* speed of light */

#define M2F 0.3048				/* multiply meters by this to get feet */
#define F2M 3.28084				/* multiply feet by this to get meters */
#define feet					/* empty label for documentation */
#define meters					/* empty label for documentation */

#define BUFLEN	256

/*--------------------------------------------------------------------------*/

#define NZONES	5				/* number of fresnel zones to look at */

double frequency = 2.4000e9;	/* lowpoint of 2.4GHz band, worst case */
double wavelength;				/* will be calculated from frequency */

/*--------------------------------------------------------------------------*/
int mprintf( const char *fmt, ... )
{
int n;
va_list args;
va_start(args,fmt);
vfprintf(LOGfp,fmt,args);
n = vfprintf(stderr,fmt,args);
va_end(args);
return(n);
}
/*--------------------------------------------------------------------------*/

typedef struct obstruction OBSTRUCTION;
struct obstruction
	{
	OBSTRUCTION *succ, *pred;
	double height, r1, r2;
	char *name;
	};

OBSTRUCTION *OBSchain;

/*--------------------------------------------------------------------------*/

/*--------------------------------------------------------------------------*/
/* Estimate Loss from Fresnel-Kirchhoff parameter
/*--------------------------------------------------------------------------*/

double FKDloss( double v )
{
if( v <= -1.0 )
	return( 0 );
else if( v <= 0.0 )
	return( 20.0 * log10( 0.5 - 0.62 * v ) );
else if( v <= 1.0 )
	return( 20.0 * log10( 0.5 * exp(-0.95 * v ) ) );
else if( v <= 2.4 )
	return( 20.0 * log10( 0.4 - sqrt(0.1184 - (0.38-0.1*v)*(0.38-0.1*v)) ) );
else
	return( 20.0 * log10( 0.225 / v ) );
}

/*--------------------------------------------------------------------------*/
/* FresnelZoneRadius(n,d) - radius of n-th Fresnel zone
/* On a straight line of sight path (heights ignored.)
/* Total path length is D, distance from source is d. 
/* MKS units. No sanity check on d.
/*--------------------------------------------------------------------------*/

double FresnelZoneRadius( int n, double D, double d )
{
return( sqrt( n * wavelength * d * (D - d) / D ) );
}

/*--------------------------------------------------------------------------*/
/* Print first k Fresnel Zone radii. Arguments in meters, output in feet.
/*--------------------------------------------------------------------------*/

void print_fresnel_radii_header( int tab, int k )
{
int n;
if(tab) mprintf("\t");
mprintf("Dist(ft)");
for( n = 1; n < (k+1); n++ ) mprintf("  zone %d",n);
mprintf("\n");
}

void print_fresnel_radii( int tab, int k, double D, double d )
{
int n;
if(tab) mprintf("\t");
mprintf("%8.2lf",M2F*d);
for( n = 1; n < (k+1); n++ )
	mprintf("  %6.2lf", M2F * FresnelZoneRadius(n,D,d) );
mprintf("\n");
}

/*--------------------------------------------------------------------------*/

/*--------------------------------------------------------------------------*/

int compare_by_distance( const void *Va, const void *Vb )
{
OBSTRUCTION *OBSa = (OBSTRUCTION *)Va, *OBSb = (OBSTRUCTION *)Vb;
if(OBSa->height < OBSb->height) return(-1);
if(OBSa->height > OBSb->height) return(1);
return(0);
}

int compare_by_height( const void *Va, const void *Vb )
{
OBSTRUCTION *OBSa = (OBSTRUCTION *)Va, *OBSb = (OBSTRUCTION *)Vb;
if(OBSa->r1 < OBSb->r1) return(-1);
if(OBSa->r1 > OBSb->r1) return(1);
return(0);
}

/*--------------------------------------------------------------------------*/

/*--------------------------------------------------------------------------*/

int main(void)
{
int i,k,n,nPeaks;
char buf[BUFLEN],*cp;
OBSTRUCTION *OBSp, *OBSq;
OBSTRUCTION **OBSdistance, **OBSheight;

double B;	/* distance between near and far antennas, horizontal */
double D;	/* distance between near and far antennas, spatial */
double H1;	/* height of first antenna */
double H2;	/* height of second antenna */

double v;	/* Fresnel-Kirchoff diffraction parameter */
double h,d,x,y,dx;

LOGfp = fopen("f.log","w");
if(!LOGfp) exit(1);

wavelength = c / frequency;
mprintf("Welcome to analysis of path loss due to diffraction.\n\n");
mprintf("Defaults [Frequency = %.4leHz, Wavelength = %.4lfcm]\n\n"
	,frequency,wavelength*100.0
	);

fprintf(stderr,"Enter alternative frequency or 0 to take default: ");
fgets(buf,BUFLEN,stdin);
sscanf(buf,"%lf%n",&x,&i);
cp = buf + i; *(cp + strlen(cp) - 1) = '\0';
if(x)
	{
	if(*cp)
		if(!strcasecmp(cp,"thz")) x *= 1e12;
		else if(!strcasecmp(cp,"ghz")) x *= 1e9;
		else if(!strcasecmp(cp,"mhz")) x *= 1e6;
		else if(!strcasecmp(cp,"khz")) x *= 1e3;
		else if(strcasecmp(cp,"hz"))
			{ mprintf("Bad unit \"%s\", forget it.\n",cp); exit(0); }
	frequency = x;
	wavelength = c / frequency;
	mprintf("Using [Frequency = %.4leHz, Wavelength = %.4lfcm]\n\n"
		,frequency,wavelength*100.0
		);
	}

fprintf(stderr,"Enter horizontal distance (ft) between antennas: ");
scanf("%lf",&B);
if(B < 0) { mprintf("Goodbye.\n"); exit(0); }
if(B < 100)
	{ mprintf("I'd rather look at 100 or more feet... goodbye.\n"); exit(0); }
B *= F2M;

fprintf(stderr,"Enter height (ft) above sea level of first antenna: ");
scanf("%lf",&H1);
H1 *= F2M;

fprintf(stderr,"Enter height (ft) above sea level of second antenna: ");
scanf("%lf",&H2);
H2 *= F2M;

D = H1 - H2; D = sqrt( B * B + D * D );

mprintf("Spatial line of sight distance is %.2lf feet...\n",M2F*D);

/*--------------------------------------------------*/
/* Print table of Fresnel Zone radii
/*--------------------------------------------------*/

mprintf("Fresnel Zone Radii. Minimum clearance = 0.55 * r(z1).\n\n" );

/* Establish not-too-fine graduation of path length */
/*if(D > 3048 meters) dx = 1000.0 feet; else */
if(D > 304.8 meters) dx = 100.0 feet;
else if(D > 30.48 meters) dx = 50.0 feet;
else dx = 10.0 feet;
dx *= F2M;

print_fresnel_radii_header(0,NZONES);
for( d = dx; d < (M2F*D/2.0); d += dx) print_fresnel_radii(0,NZONES,D,d);
print_fresnel_radii(0,NZONES,D,D/2.0);	/* point of maximum radii */

#if 0
/*--------------------------------------------------*/
/* Obstruction analysis
/*--------------------------------------------------*/

mprintf("\nPath Analysis ...\n");

mprintf("Enter a series of peaks...\n");
for( nPeaks = 0; ; )
	{
	mprintf("\tEnter height (ft) above sea level: ");
	scanf("%lf",&h); h *= F2M;
	mprintf("\tEnter distance (ft) from first antenna (0 to end data entry): ");
	scanf("%lf",&d); if(!d) break; d *= F2M;
	/*mprintf("\tYou can even give this peak a name! Wow: ");*/
	bzero(buf,BUFLEN); /*fgets(buf,BUFLEN,stdin);*/
	i = strlen(buf); *(buf + i - 1) = '\0';
	OBSp = (OBSTRUCTION *)malloc(sizeof(OBSTRUCTION) + i);
	OBSp->height = h; OBSp->r1 = d; OBSp->r2 = B - d;
	OBSp->name = (char *)OBSp + sizeof(OBSTRUCTION);
	strcpy(OBSp->name,buf);
	if(!OBSchain) OBSchain = OBSp->succ = OBSp->pred = OBSp;
	else
		{
		OBSp->succ = OBSchain;
		OBSq = OBSp->pred = OBSchain->pred;
		OBSchain->pred = OBSq->succ = OBSp;
		}
	nPeaks++;
	}

/* Sort pointers by both distance and height */
k = n * sizeof(OBSTRUCTION *);
OBSdistance = (OBSTRUCTION **)malloc(k);
OBSheight = (OBSTRUCTION **)malloc(k);
for( OBSp = OBSchain, i = nPeaks; --i >= 0; )
	{ OBSdistance[i] = OBSp; OBSp = OBSp->succ; if(OBSp == OBSchain) break; }
bcopy(OBSdistance,OBSheight,k);
qsort(OBSdistance,nPeaks,sizeof(OBSTRUCTION *),compare_by_distance);
qsort(OBSheight,nPeaks,sizeof(OBSTRUCTION *),compare_by_height);

for( i = 0; i < nPeaks; i++ )
	{
	OBSp = OBSdistance[i];
	if(*(OBSp->name)) mprintf("Peak %s:\n",OBSp->name);
		else mprintf("Peak %d:\n",i);
	mprintf("\tr1 = %.1lfft, height = %.1lfft\n",M2F*OBSp->height,M2F*OBSp->r1);
	x = H1 - OBSp->height; x = sqrt( x * x + OBSp->r1 * OBSp->r1 );
	mprintf("\tSpatial path distance to peak is %.1lfft.\n",M2F*x);
	print_fresnel_radii_header(1,NZONES);
	print_fresnel_radii(1,NZONES,D,x);
	mprintf("\tMinimum clearance for zone 1 is %.1lf feet.\n"
		, 0.55 * M2F * FresnelZoneRadius(1,D,x)
		);
	/* Diffraction loss assuming this is the only obstruction */
	y = atan2( OBSp->height - H1, D - x ) + atan2(OBSp->height,x);
	v = y * sqrt( 2.0 * x * (D - x) / (wavelength * D) );
	mprintf("\tFresnel-Kirchoff diffraction parameter is %.3lf.\n",v);
	mprintf("\tDiffraction loss is roughly %.1lf dB.\n",FKDloss(v));
	}

/* Picquenard Analysis */
mprintf("Picquenard analysis coming soon ...\n");

#endif

}

/*--------------------------------------------------------------------------*/
/* end - fresnel.c
/*--------------------------------------------------------------------------*/