/*--------------------------------------------------------------------------*/
/* fpath.c - full-fledged path analysis
/*--------------------------------------------------------------------------*/
/*
$Header: /usr/home/ecsd/RCS/fpath.c,v 1.10 1998/06/09 02:01:58 ecsd Exp ecsd $
$Log: fpath.c,v $
Revision 1.10  1998/06/09 02:01:58  ecsd
moving to reading in data, not assumed ordered.

Revision 1.9  1998/06/08 01:59:03  ecsd
seem to have worked out zone calculations and diffraction loss.
added in freespace loss - needs to be checked, seems too high.
check that algorithm ends properly, and two other flagged tasks.
appears to work for diffraction ...
next fix input of problem for more testing.

Revision 1.8  1998/06/05 21:07:03  ecsd
starting to debug

Revision 1.7  1998/06/05 01:31:48  ecsd
first working version, numbers are bad, loss should be gain, etc.
cleaned up, some intended code excised.

Revision 1.6  1998/06/04 02:22:36  ecsd
add FresnelZoneCutFromContext(s,z,t) and more documentation and credits/citations

Revision 1.5  1998/06/03 03:43:13  ecsd
replace inexact FKDloss() with exactFKloss() based on direct
computation of the complex curve. In particular, do not clip
gain to 0 when v > 1.0. Two source books assign different signs
to v, so exactFKloss() computes against -v.

Revision 1.4  1998/06/02 02:30:52  ecsd
add fresnel integral function, verified accuracy to e-4

Revision 1.3  1998/05/20 05:24:41  ecsd
add fresnel.c, more methodology
better approach to uniform problem

Revision 1.2  1998/05/19 05:06:11  ecsd
in transition from a.i to a.ij ... picquenard for more than 2?

Revision 1.1  1998/05/19 05:01:26  ecsd
Initial revision
*/
/*--------------------------------------------------------------------------*/

/*--------------------------------------------------------------------------*/
/*
Original 'fresnel.c' 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
rev 1.4 May 11th did first steps toward total path analysis,
	add frequency override (for testing.) has calculation bugs somewhere.

Approximation in FKDloss() given by William C. Y. Lee in
"Mobile Communications Engineering, Theory and Applications", 2ed.,
McGraw Hill, 1997. exactFKloss() also from WCY Lee, translating formula
4.14, section 4.2, p.141, using the Fresnel integral directly rather than
his approximation p.142, since we have no difficulty with the integral.
*/
/*--------------------------------------------------------------------------*/

/*--------------------------------------------------------------------------*/

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

/*--------------------------------------------------------------------------*/
/* Physics
/*--------------------------------------------------------------------------*/
#define PI		3.141592653589793
#define PI2		1.5707963267948965
#define RAD		57.29577951		/* degrees per radian */

#define M2F		.3048000		/* M2F * meters = feet */
#define F2M		3.280840		/* F2M * feet = meters */
#define feet					/* shows unit of a value */
#define meters					/* shows unit of a value */
#define dB						/* shows unit of a value */

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

#define FREQUENCY_DEFAULT	2.4425e9	/* midpoint of 2.4GHz band */

double frequency = FREQUENCY_DEFAULT;	/* midpoint of 2.4GHz band */
double wavelength;				/* calculate as (c/frequency) on change f */

/*--------------------------------------------------------------------------*/

/****************************************************************************/
/* fresnel integrals. code adapted from:
/* Cephes Math Library Release 2.1:  January, 1989
/* Copyright 1984, 1987, 1989 by Stephen L. Moshier
/* -ecsd: remove non IBMPC data (want to rewrite hex-coded constants as doubles)
/*--------------------------------------------------------------------------*/

/* S(x) for small x */
static short sn[24] = {
	0xb5c3,0x6d65,0x5fa3,0xc0a7, 0x05b6,0x172c,0xa1d0,0x4125,
	0xf5ed,0x24f6,0x0746,0xc18e, 0x7afe,0x60b7,0xfda8,0x41e2,
	0xf087,0x347b,0xa0ba,0xc224, 0x2457,0xe6e5,0x82cf,0x4252,
	};
static short sd[24] = {
	0x1072,0x3287,0x9605,0x4071, 0xdaa9,0xfe9c,0x4218,0x40e6,
	0x17b4,0xb8d0,0xbc2f,0x4153, 0xea9e,0xb5e5,0xfe51,0x41b8,
	0x22e9,0xe6b2,0xe664,0x4214, 0x42b5,0x73de,0xad3b,0x4261,
	};

/* C(x) for small x */
static short cn[24] = {
	0x62d3,0x2cfb,0xc80c,0xbe6a, 0x0f74,0x6210,0xee92,0x3ee3,
	0x8786,0x0ed6,0x2442,0xbf45, 0xcc13,0x1059,0x566a,0x3f93,
	0x9690,0x378a,0x4eac,0xbfca, 0x0000,0x0000,0x0000,0x3ff0,
	};
static short cd[28] = {
	0xc6b3,0x6a7f,0x9768,0x3d91, 0x75b9,0xdb03,0x7449,0x3e0f,
	0x5191,0x02a4,0xc708,0x3e80, 0x1d32,0xfa2a,0xa3ee,0x3ee9,
	0xd99d,0x3fdb,0x718f,0x3f4c, 0xfe74,0x6030,0x1a07,0x3fa5,
	0x0000,0x0000,0x0000,0x3ff0,
	};

/* Auxiliary function f(x) */
static short fn[40] = {
	0x2391,0xd1b0,0xfa91,0x3fda, 0x0e5e,0xd2b9,0x5b30,0x3fc2,
	0x7dad,0x7b15,0x98e5,0x3f87, 0x0efc,0xc495,0x9c70,0x3f36,
	0x7795,0xc535,0x7203,0x3ed3, 0x87a1,0x484e,0x67b5,0x3e60,
	0xa7bb,0x4998,0x1f0a,0x3ddc, 0xfdf7,0x1459,0x3557,0x3d48,
	0xa9e4,0xaf8f,0x5a2d,0x3ca3, 0x05f6,0x0e0b,0x36ef,0x3be6,
	};
static short fd[40] = {
	0x7a65,0xeb21,0x0cfe,0x3fe8, 0x6d7d,0xc1d4,0xec6e,0x3fbd,
	0x56f3,0xa6ed,0x615e,0x3f7a, 0xbd60,0x6017,0x704a,0x3f24,
	0x64f5,0x9232,0xf9b1,0x3ebe, 0x5916,0x9552,0x33b4,0x3e48,
	0xa061,0x33d3,0xcc85,0x3dc3, 0x0790,0x6b9f,0x926c,0x3d30,
	0xf48d,0x2352,0x0e5d,0x3c8a, 0x6141,0x12b9,0x9e94,0x3bcd,
	};

/* Auxiliary function g(x) */
static short gn[44] = {
	0xcaca,0xb420,0x2463,0x3fe0, 0x7481,0x67f8,0x3aaa,0x3fc9,
	0x3214,0x54ba,0x3718,0x3f93, 0xb33b,0x48d1,0x6a79,0x3f46,
	0xabf3,0xf9fa,0x2577,0x3ee8, 0x6091,0x4dea,0x621c,0x3e7a,
	0xeb09,0xf200,0x9a94,0x3dfe, 0x89cb,0x8766,0x0bf5,0x3d73,
	0xa964,0x3df8,0xc7a0,0x3cd8, 0x58a6,0xf173,0xdb24,0x3c2e,
	0xbe2b,0x624f,0x409d,0x3b6c,
	};
static short gd[44] = {
	0x55c1,0x23bc,0x996d,0x3ff7, 0xc34f,0xefa1,0x9dad,0x3fd5,
	0xaa09,0xe637,0xf811,0x3f99, 0xae28,0x0fa0,0xb206,0x3f4a,
	0x067b,0x2c0a,0xbf86,0x3eea, 0x7428,0xab1c,0x0071,0x3e7c,
	0xf1c6,0x8e3c,0xa861,0x3dff, 0xef2b,0x9c3d,0x6643,0x3d73,
	0x1936,0x37c8,0x00dc,0x3cd9, 0x8364,0x84ff,0xf5a1,0x3c2e,
	0xbe2b,0x624f,0x409d,0x3b6c,
	};

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

static double polevl( double x, short *coef, int n )
{
int i = n; double *p = (double *)coef; double ans = *p++;
do ans = ans * x  +  *p++; while( --i );
return(ans);
}

static double p1evl( double x, short *coef, int n )
{
double *p = (double *)coef; int i = n - 1; double ans = x + *p++;
do ans = ans * x  + *p++; while( --i );
return(ans);
}

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

void fresnelCS( double X, double *Cp, double *Sp )
{
double f, g, C, S, c, s, t, u, x, x2;

x = fabs(X);
x2 = x * x;
if( x > 36974.0 ) { C = 0.5; S = 0.5; }
else if( x2 < 2.5625 )
	{
	t = x2 * x2;
	S = x * x2 * polevl( t, sn, 5) / p1evl( t, sd, 6 );
	C = x * polevl( t, cn, 5) / polevl(t, cd, 6 );
	}
else
	{ /* Auxiliary functions for large argument  */
	u = PI * x2;
	t = 1.0/u;
	u = 1.0/(u * u);
	f = 1.0 - u * polevl( u, fn, 9) / p1evl(u, fd, 10);
	g = t * polevl( u, gn, 10) / p1evl(u, gd, 11);
	t = PI2 * x2;
	c = cos(t);
	s = sin(t);
	t = PI * x;
	C = 0.5  +  (f * s  -  g * c) / t;
	S = 0.5  -  (f * c  +  g * s) / t;
	}

if( X < 0.0 ) { C = -C; S = -S; }
*Cp = C;
*Sp = S;
return;
}
/****************************************************************************/

/*--------------------------------------------------------------------------*/
/*
/*            X
/*          / | \    %
/*         /  |  \          A
/*        /b  |   \       / | \
/*       *----------     -------
/*       |    |*            |
/*       |h0  |h1     *     |hr          %
/*       |    |           D | *
/*       |    |             |         *
/*       |    |             |                 *
/*       ---------------------------------------------*hn
/*       |<-->|                                       |
/*       d0   d1                                      dn
/*
/*--------------------------------------------------------------------------*/

typedef struct loc
	{
	double d;			/* distance coordinate						*/
	double h;			/* height coordinate						*/
	char *name;
	double *Beta;		/* forward elevation angle to point + i		*/
	double *Ipsilon;
	double FresnelZoneCut;
	double HeightAbovePath;
	int index;
	} L;

/* array of pointers to points sorted by height over/under source->target path
*/
L **Lh;

/* Initializer */
#define pt(d,h,name) { (d), (h), name }

/* pt(distance,height) - fyi */
L Ldata[] = { pt(0,100,"S")
	, pt(200, 70,"Santa Monica")
	, pt(300,240,"My big Nose")
	, pt(400,70,"Tokyo")
	, pt(1000, 70,"T")
	};
int Lcount = ( sizeof(Ldata) / sizeof(L) );

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

int L_compare_by_distance( const void *Va, const void *Vb )
{
L *La = (L *)Va, *Lb = (L *)Vb;
if(La->h < Lb->h) return(-1);
if(La->h > Lb->h) return(1);
return(0);
}

int L_compare_by_height( const void *Va, const void *Vb )
{
L *La = (L *)Va, *Lb = (L *)Vb;
if(La->d < Lb->d) return(-1);
if(La->d > Lb->d) return(1);
return(0);
}

/*--------------------------------------------------------------------------*/

/*--------------------------------------------------------------------------*/
/* Free Path Loss from source to target.
/* P = txpower * (GG = ant1gain * ant2gain)
/* p/P = GG * (wavelength/4*pi*d)^2
/* loss = 10*log10(p/P) = 20 ( log10(w/4pi) - log10(d) )
/*--------------------------------------------------------------------------*/

double FreePathLength( L *Ls, L *Lt )
{
double h = Ls->h - Lt->h, d = Lt->d - Ls->d;
return( sqrt(h*h + d*d) );
}

double PathLoss( double d, double wavelength )
{
double x = 4.0 * PI * d;
double z = 20.0 * log10(x);
double y = 20.0 * log10(wavelength);
printf("\n\t***pathloss(%.2lf) = %.2lf, wavelength/(d4pi) %.2lf/%.2lf\n"
		"\t*** = 20log10(w)] %+.3lf - 20log10(d4pi) %.3lf"
	, d, y - z, wavelength, x, y, z
	);
return( y - z );
}

/*--------------------------------------------------------------------------*/
/* Calculate loss from WCY Lee (v negative of Rappaport)
/* 4.13) Loss = 20 log10(F)
/* 4.14) F = (S + 0.5) / ( sqrt(2) * sin(Dphi+PI/4) )
/* ....) Dphi = atan2( (S+0.5) / (C+0.5) ) - PI/4
/* 4.15,6) C = C(v), S = S(v),
/* using -v with v = h * sqrt( (2/wavelength) * (1/r1 + 1/r2) )
/* recast F as spherical radius (length of line from 0 to pt on Cornu spiral)
/*--------------------------------------------------------------------------*/

double exactFKloss( double v )
{
double C, S, c, s, F;
fresnelCS(-v,&C,&S);
c = C + 0.5; s = S + 0.5; F = sqrt( (c*c + s*s) / 2.0 );
return( 20.0 * log10(F) );
}

/*--------------------------------------------------------------------------*/
/* Estimate Loss from Fresnel-Kirchhoff parameter (v from Rappaport)
/* This clips gain at 0 going positive
/*--------------------------------------------------------------------------*/

double FKDloss( double v )
{
double x;
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 )
	{ x = 0.38 - 0.1 * v; return( 20.0 * log10(0.4 - sqrt(0.1184 - x*x)) ); }
else
	return( 20.0 * log10( 0.225 / v ) );
}

/*--------------------------------------------------------------------------*/

/*--------------------------------------------------------------------------*/
/* FresnelZoneRadius(n,d,D) - radius of n-th Fresnel zone. MKS units.
/* On a direct line of sight path of length D.
/* Distance from source is d, not sanity checked.
/*--------------------------------------------------------------------------*/

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

/*--------------------------------------------------------------------------*/
/* FresnelZoneCutFromContext(Ls,Lz,Lt) -
/* On the base path Ls --> Lt, which zone intersected by Lz?
/* Calculation accounts for nonhorizontality of base path
/* Negative result is below base path
/*--------------------------------------------------------------------------*/

double FresnelZoneCutFromContext( L *Ls, L *Lz, L *Lt )
{
double ho, phz, ro,ro1,ro2, rr,rr1,rr2;
double Alpha, tanAlpha, delta, Ho, Hr, ozone, rzone;

printf("\n\tFZCFC( %s(%.2lf,%.2lf), %s(%.2lf,%.2lf), %s(%.2lf,%.2lf) ):\n"
	, Ls->name, Ls->d, Ls->h
	, Lz->name, Lz->d, Lz->h
	, Lt->name, Lt->d, Lt->h
	);

/* problem as yet unrotated onto base path */
ho = Ls->h - Lt->h; ro = Lt->d - Ls->d; tanAlpha = ho / ro;
ro1 = Lz->d - Ls->d; ro2 = Lt->d - Lz->d;

/* rotation parameters */
Alpha = atan2(ho,ro);
phz = Lt->h + ro2 * tanAlpha;
Ho = Lz->h - phz;
delta = Ho * sin(Alpha);
ozone = copysign(ro*Ho*Ho/(wavelength*ro1*ro2),Ho);

/* rotated problem */
Hr = Ho * cos(Alpha);
rr = ro / cos(Alpha);
rr1 = ro1 / cos(Alpha) + delta;
rr2 = ro2 / cos(Alpha) - delta;
rzone = copysign( rr * Hr * Hr / (wavelength * rr1 * rr2), Hr );

printf(
	"\t\tUnrotated: ho %.1lf, ro %.1lf, tanAlpha %.3lf,"
		" ro1 %.1lf, ro2 %.1lf\n"
	"\t\tZone %+.3lf\n"
	"\t\tAlpha %+.2le phz %.2lf, Ho %+.2lf = Hz - phz, delta-r %+.2lf\n"
	"\t\tRotated: Hr %.2lf rr %.2lf rr1 %.2lf rr2 %.2lf\n"
	"\t\tZone %+.3lf\n"
	, ho, ro, tanAlpha, ro1, ro2, ozone
	, Alpha, phz, Ho, delta
	, Hr, rr, rr1, rr2, rzone
	);

return(rzone);
}

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

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

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

#define BUFLEN	256
#define NZONES	8				/* number of fresnel zones to look at */

void main(void)
{
int i, j, k, nzones = NZONES, fx = -1, rx = -1;
L *Lp, *Lq, *Lr, *Lsource, *Ltarget, *Lpeak;
char buf[BUFLEN],*cp;
/* accept any elevation as maximum, not just one exceeding base LOS	*/
double fpeak = -PI2, rpeak = -PI2;
double IPSILON, tanIPSILON;
double x, h, r;
double pLoss, FreespaceLoss, dLoss, DiffractionLoss, TotalLoss;

/* Determine frequency for this problem ------------------------------------*/

wavelength = c / frequency;
printf("Welcome to analysis of path loss due to diffraction.\n");
printf("Supposed Picquenard method, implementation should be verified.\n");
printf("\nDefaults are Frequency = %.4leHz, Wavelength = %.4lfcm\n"
	,frequency,wavelength*100.0
	);

fprintf(stderr,"\nEnter alternative frequency or 0 to take default: ");
if(fgets(buf,BUFLEN,stdin))
	if( 1 == 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"))
				{ printf("Bad unit \"%s\", forget it.\n",cp); exit(0); }
		frequency = x;
		wavelength = c / frequency;
		printf("Now using Frequency = %.4leHz, Wavelength = %.4lfcm\n"
			,frequency,wavelength*100.0
			);
		}
	}
/*--------------------------------------------------------------------------*/

/*--------------------------------------------------------------------------*/
#if 0
printf("Enter a series of peaks...\n");
for( name = "A", nPeaks = 0, namespace = 0, Lq = 0; ; )
	{
	printf("\tEnter height(ft),distance(ft),name: ");
	bzero(buf,BUFLEN);
	if(feof(fp)) break;
	fscanf(fp,"%lf,%lf,%s",&h,&d,buf);
	if(i != 3) { printf("data error, read %d values. bye.\n",i); exit(1); }
	h *= F2M; d *= F2M;
	i = strlen(buf); *(buf + i - 1) = '\0';
	Lp = (L *)malloc(sizeof(L) + i); namespace += i;
	Lp->h = h; Lp->d = d;
	Lp->name = (char *)Lp + sizeof(L);
	strcpy(Lp->name,buf);
	if(Lq) Lq->succ = Lp; Lq = Lp;
	if( *name == "Z" ) *name = "a"; else (*name)++;
	nPeaks++;
	}
if(Lq) Lq->succ = 0;
Ldata = (L *)malloc(sizeof(L) + namespace);
if(!Ldata) { printf("Can't malloc array. bye.\n"); exit(1); }
for( cp = (char *)Ldata + nPeaks * sizeof(L), i = 0; i < nPeaks; i++ )
	{
	bcopy(Lp,Ldata+i,sizeof(L));
	Ldata[i].name = cp; strcpy(cp,Lp->name);
	cp += strlen(cp); *cp++ = '\0'; Lp = Lp->succ;
	}

/* Sort point-objects physically by distance */
qsort(Ldata,nPeaks,sizeof(L),compare_by_distance);
#endif
/*--------------------------------------------------------------------------*/

/* problem attributes ------------------------------------------------------*/

Lsource = Ldata;
Ltarget = Ldata + (Lcount - 1);
tanIPSILON = (Lsource->h - Ltarget->h) / (Ltarget->d - Lsource->d);
IPSILON = atan2( Lsource->h - Ltarget->h, Ltarget->d - Lsource->d );
k = Lcount * sizeof(double);
for( Lp = Lsource; Lp <= Ltarget; Lp++ )
	{
	/* Allocate space for Beta and Ipsilon arrays */
	Lp->Beta = (double *)malloc(k); bzero(Lp->Beta,k);
	Lp->Ipsilon = (double *)malloc(k); bzero(Lp->Ipsilon,k);
	}
/*--------------------------------------------------------------------------*/

#if 0
/* calculate attribute[i,j] for j > i --------------------------------------*/
for( Lp = Lsource; Lp < Ldata + Lcount - 1; Lp++ )
	for( Lq = Lp; ++Lq <= Ldata + Lcount; )
		{
		/* height of this point, with zero = height of source->target path */
		/* height-over-path of point */
		Lp->HeightAbovePath = Lp->h - (Ltarget->d - Lp->d) * tanIPSILON;

		/* angle from horizontal of line from point i to point j */
		Lp->Beta[j-i] = atan2( h[j] - h[i], d[j] - d[i] );

		/* angle from horizontal of line from point j to point i */
		Lp->Ipsilon[j-i] =
			PI2 - ( Lp->Beta[j-i] = atan2( h[i] - H, DT - d[i] ) );

		if(Lp != Ltarget)
			{
			*(Lp->Beta + (Lp-Lsource))
				= atan2( h[i+1] - h[i], d[i+1] - d[i] );
			fwd_elevation_angle[i][i+1] = atan2( d[i+1] - d[i], h[i+1] - h[i] );
			}
		}
/*--------------------------------------------------------------------------*/
#endif

#if 0
/* determine convex cover of path ------------------------------------------*/
for( i = 0; i < Lcount; i++ )
	{
	if( fwd_elevation_angle[i][0] > fpeak )
		{ fx = i; fpeak = fwd_elevation_angle[i][0]; }
	if( beta[i][0] > rpeak ) { rx = i; rpeak = beta[i][0]; }
	}
/*--------------------------------------------------------------------------*/
#endif

/* calculate attribute[i,j,k] for i < j < k --------------------------------*/
/*--------------------------------------------------------------------------*/

/*--------------------------------------------------------------------------*/
/* Picquenard method, deals with any number of obstructions.
/* Beginning at source, calculate loss due to first obstruction, then reset
/* path to be path from current obstruction to target and calculate next loss
/* according to height of next obstruction over new base path. Sum losses
/* until done, total loss is effective path loss from original source.
/*
/* For those points below base path, report path loss in any case. This can
/* actually be a gain! It does contribute, so report it now to consider
/* factoring it into the calculation.
/*
/*<> Knife-edge only. Shape of terrain not computed here.
/*<> Factor in path loss on each hop?
/*<> algorithm does not sum properly at last phase (check end condition)
/*<> does path = path over cover or sum(loss) peak to peak? ploss bugs...
/*--------------------------------------------------------------------------*/

TotalLoss = FreespaceLoss = DiffractionLoss = 0.0 dB;

for( Lpeak = Lp = Lsource; ++Lp < Ltarget; )
	{
	printf("\nPoint %s", Lp->name);
	pLoss = PathLoss(FreePathLength(Lpeak,Lp),wavelength);
	Lp->FresnelZoneCut = FresnelZoneCutFromContext(Lpeak,Lp,Ltarget);
	dLoss = exactFKloss(Lp->FresnelZoneCut);
	printf(" Zone %+5.3lf dLoss %+.3lf dB pLoss %+.3lf dB ...\n"
		, Lp->FresnelZoneCut, dLoss, pLoss
		);

	/* if zone < 0, report loss and proceed */
	if(Lp->FresnelZoneCut < 0.0)
		{
		/* printf("ignored.\n"); */
		DiffractionLoss += dLoss;
		TotalLoss += dLoss;
		printf("absorbed, total diffraction/freespace/path loss"
				" %+.3lf dB/%.3lf dB/%.3lf dB\n"
				"Retain source at point %s.\n"
			, DiffractionLoss, FreespaceLoss, TotalLoss, Lpeak->name
			);
		continue;
		}

	/* if zone >= 0, absorb loss and make this point the new source */
	DiffractionLoss += dLoss;
	FreespaceLoss += pLoss;
	TotalLoss += (dLoss + pLoss);
	Lpeak = Lp;
	printf("absorbed, total diffraction/freespace/path loss"
			" %+.3lf dB/%+.3lf dB/%+.3lf dB\n"
			"New source at point %s.\n"
		, DiffractionLoss, FreespaceLoss, TotalLoss, Lpeak->name
		);
	}

/* Clean up. Add in free path loss from last peak to target. */
pLoss = PathLoss(FreePathLength(Lpeak,Ltarget),wavelength);
FreespaceLoss += pLoss;
TotalLoss += pLoss;
printf("\nFinal peak point %s.\n"
		"Freespace loss %+5.3lf dB.\n"
		"Total diffraction loss %+5.3lf dB.\n"
		"Total freespace loss %+5.3lf dB.\n"
		"Total path loss %+5.3lf dB.\n"
	, Lpeak->name, pLoss, DiffractionLoss, FreespaceLoss, TotalLoss
	);

printf("\nCiao.\n");
}

/*--------------------------------------------------------------------------*/
/* end - fpath.c
/*--------------------------------------------------------------------------*/

