HFNOC/Signal-Server
2
0
Fork 0
forked from ExternalVendorCode/Signal-Server
Code Pull Requests Activity
Files
147617bb5dba49e91b6605f42018a77d61477551
Signal-Server/main.cc
root 635b328871 Meridian cleanup and tested for 30m.
2018-09-25 21:37:12 +01:00

1969 lines
50 KiB
C++
Raw Blame History

double version = 3.10;
/****************************************************************************\
* Signal Server: Radio propagation simulator by Alex Farrant QCVS, 2E0TDW *
******************************************************************************
* SPLAT! Project started in 1997 by John A. Magliacane, KD2BD *
******************************************************************************
* Please consult the SPLAT! documentation for a complete list of *
* individuals who have contributed to this project. *
******************************************************************************
* *
* This program is free software; you can redistribute it and/or modify it *
* under the terms of the GNU General Public License as published by the *
* Free Software Foundation; either version 2 of the License or any later *
* version. *
* *
* This program is distributed in the hope that it will useful, but WITHOUT *
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or *
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License *
* for more details. *
* *
\****************************************************************************/
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <unistd.h>
#include <errno.h>
#include <limits.h>
#include "common.h"
#include "inputs.hh"
#include "outputs.hh"
#include "models/itwom3.0.hh"
#include "models/los.hh"
#include "models/pel.hh"
#include "image.hh"
int MAXPAGES = 10*10;
int IPPD = 1200;
int ARRAYSIZE = (MAXPAGES * IPPD) + 10;
char string[255], sdf_path[255], udt_file[255], opened = 0, gpsav =
0, ss_name[16], dashes[80];
double earthradius, max_range = 0.0, forced_erp, dpp, ppd, yppd,
fzone_clearance = 0.6, forced_freq, clutter, lat, lon, txh, tercon, terdic,
north, east, south, west, dBm, loss, field_strength,
min_north = 90, max_north = -90, min_west = 360, max_west = -1, westoffset=180, eastoffset=-180, delta=0, rxGain=0,
cropLat=-70, cropLon=0,cropLonNeg=0;
int ippd, mpi,
max_elevation = -32768, min_elevation = 32768, bzerror, contour_threshold,
pred, pblue, pgreen, ter, multiplier = 256, debug = 0, loops = 100, jgets =
0, MAXRAD, hottest = 0, height, width, resample = 0;
unsigned char got_elevation_pattern, got_azimuth_pattern, metric = 0, dbm = 0;
bool to_stdout = false, cropping = true;
__thread double *elev;
__thread struct path path;
struct site tx_site[2];
struct dem *dem;
struct LR LR;
struct region region;
double arccos(double x, double y)
{
/* This function implements the arc cosine function,
returning a value between 0 and TWOPI. */
double result = 0.0;
if (y > 0.0)
result = acos(x / y);
if (y < 0.0)
result = PI + acos(x / y);
return result;
}
int ReduceAngle(double angle)
{
/* This function normalizes the argument to
an integer angle between 0 and 180 degrees */
double temp;
temp = acos(cos(angle * DEG2RAD));
return (int)rint(temp / DEG2RAD);
}
double LonDiff(double lon1, double lon2)
{
/* This function returns the short path longitudinal
difference between longitude1 and longitude2
as an angle between -180.0 and +180.0 degrees.
If lon1 is west of lon2, the result is positive.
If lon1 is east of lon2, the result is negative. */
double diff;
diff = lon1 - lon2;
if (diff <= -180.0)
diff += 360.0;
if (diff >= 180.0)
diff -= 360.0;
return diff;
}
char *dec2dms(double decimal)
{
/* Converts decimal degrees to degrees, minutes, seconds,
(DMS) and returns the result as a character string. */
char sign;
int degrees, minutes, seconds;
double a, b, c, d;
if (decimal < 0.0) {
decimal = -decimal;
sign = -1;
}
else
sign = 1;
a = floor(decimal);
b = 60.0 * (decimal - a);
c = floor(b);
d = 60.0 * (b - c);
degrees = (int)a;
minutes = (int)c;
seconds = (int)d;
if (seconds < 0)
seconds = 0;
if (seconds > 59)
seconds = 59;
string[0] = 0;
snprintf(string, 250, "%d%c %d\' %d\"", degrees * sign, 176, minutes,
seconds);
return (string);
}
int PutMask(double lat, double lon, int value)
{
/* Lines, text, markings, and coverage areas are stored in a
mask that is combined with topology data when topographic
maps are generated by ss. This function sets and resets
bits in the mask based on the latitude and longitude of the
area pointed to. */
int x = 0, y = 0, indx;
char found;
for (indx = 0, found = 0; indx < MAXPAGES && found == 0;) {
x = (int)rint(ppd * (lat - dem[indx].min_north));
y = mpi - (int)rint(yppd * (LonDiff(dem[indx].max_west, lon)));
if (x >= 0 && x <= mpi && y >= 0 && y <= mpi)
found = 1;
else
indx++;
}
if (found) {
dem[indx].mask[x][y] = value;
return ((int)dem[indx].mask[x][y]);
}
else
return -1;
}
int OrMask(double lat, double lon, int value)
{
/* Lines, text, markings, and coverage areas are stored in a
mask that is combined with topology data when topographic
maps are generated by ss. This function sets bits in
the mask based on the latitude and longitude of the area
pointed to. */
int x = 0, y = 0, indx;
char found;
for (indx = 0, found = 0; indx < MAXPAGES && found == 0;) {
x = (int)rint(ppd * (lat - dem[indx].min_north));
y = mpi - (int)rint(yppd * (LonDiff(dem[indx].max_west, lon)));
if (x >= 0 && x <= mpi && y >= 0 && y <= mpi)
found = 1;
else
indx++;
}
if (found) {
dem[indx].mask[x][y] |= value;
return ((int)dem[indx].mask[x][y]);
}
else
return -1;
}
int GetMask(double lat, double lon)
{
/* This function returns the mask bits based on the latitude
and longitude given. */
return (OrMask(lat, lon, 0));
}
int PutSignal(double lat, double lon, unsigned char signal)
{
/* This function writes a signal level (0-255)
at the specified location for later recall. */
char dotfile[255];
snprintf(dotfile, 80, "%s.dot%c", tx_site[0].filename, 0);
snprintf(dotfile, 80, "%s.dot%c", tx_site[0].filename, 0);
int x = 0, y = 0, indx;
char found;
if (signal > hottest) // dBm, dBuV
hottest = signal;
//lookup x/y for this co-ord
for (indx = 0, found = 0; indx < MAXPAGES && found == 0;) {
x = (int)rint(ppd * (lat - dem[indx].min_north));
y = mpi - (int)rint(yppd * (LonDiff(dem[indx].max_west, lon)));
if (x >= 0 && x <= mpi && y >= 0 && y <= mpi)
found = 1;
else
indx++;
}
if (found) {
// Write values to file
dem[indx].signal[x][y] = signal;
return (dem[indx].signal[x][y]);
}
else
return 0;
}
unsigned char GetSignal(double lat, double lon)
{
/* This function reads the signal level (0-255) at the
specified location that was previously written by the
complimentary PutSignal() function. */
int x = 0, y = 0, indx;
char found;
for (indx = 0, found = 0; indx < MAXPAGES && found == 0;) {
x = (int)rint(ppd * (lat - dem[indx].min_north));
y = mpi - (int)rint(yppd * (LonDiff(dem[indx].max_west, lon)));
if (x >= 0 && x <= mpi && y >= 0 && y <= mpi)
found = 1;
else
indx++;
}
if (found)
return (dem[indx].signal[x][y]);
else
return 0;
}
double GetElevation(struct site location)
{
/* This function returns the elevation (in feet) of any location
represented by the digital elevation model data in memory.
Function returns -5000.0 for locations not found in memory. */
char found;
int x = 0, y = 0, indx;
double elevation;
for (indx = 0, found = 0; indx < MAXPAGES && found == 0;) {
x = (int)rint(ppd * (location.lat - dem[indx].min_north));
y = mpi -
(int)rint(yppd *
(LonDiff(dem[indx].max_west, location.lon)));
if (x >= 0 && x <= mpi && y >= 0 && y <= mpi)
found = 1;
else
indx++;
}
if (found)
elevation = 3.28084 * dem[indx].data[x][y];
else
elevation = -5000.0;
return elevation;
}
int AddElevation(double lat, double lon, double height, int size)
{
/* This function adds a user-defined terrain feature
(in meters AGL) to the digital elevation model data
in memory. Does nothing and returns 0 for locations
not found in memory. */
char found;
int i,j,x = 0, y = 0, indx;
for (indx = 0, found = 0; indx < MAXPAGES && found == 0;) {
x = (int)rint(ppd * (lat - dem[indx].min_north));
y = mpi - (int)rint(yppd * (LonDiff(dem[indx].max_west, lon)));
if (x >= 0 && x <= mpi && y >= 0 && y <= mpi)
found = 1;
else
indx++;
}
if (found && size<2)
dem[indx].data[x][y] += (short)rint(height);
// Make surrounding area bigger for wide area landcover. Should enhance 3x3 pixels including c.p
if (found && size>1){
for(i=size*-1; i <= size; i=i+1){
for(j=size*-1; j <= size; j=j+1){
if(x+j >= 0 && x+j <=IPPD && y+i >= 0 && y+i <=IPPD)
dem[indx].data[x+j][y+i] += (short)rint(height);
}
}
}
return found;
}
double dist(double lat1, double lon1, double lat2, double lon2)
{
//ENHANCED HAVERSINE FORMULA WITH RADIUS SLIDER
double dx, dy, dz;
int polarRadius=6357;
int equatorRadius=6378;
int delta = equatorRadius-polarRadius; // 21km
float earthRadius = equatorRadius - ((lat1/100) * delta);
lon1 -= lon2;
lon1 *= DEG2RAD, lat1 *= DEG2RAD, lat2 *= DEG2RAD;
dz = sin(lat1) - sin(lat2);
dx = cos(lon1) * cos(lat1) - cos(lat2);
dy = sin(lon1) * cos(lat1);
return asin(sqrt(dx * dx + dy * dy + dz * dz) / 2) * 2 * earthRadius;
}
double Distance(struct site site1, struct site site2)
{
/* This function returns the great circle distance
in miles between any two site locations. */
double lat1, lon1, lat2, lon2, distance;
lat1 = site1.lat * DEG2RAD;
lon1 = site1.lon * DEG2RAD;
lat2 = site2.lat * DEG2RAD;
lon2 = site2.lon * DEG2RAD;
distance =
3959.0 * acos(sin(lat1) * sin(lat2) +
cos(lat1) * cos(lat2) * cos((lon1) - (lon2)));
return distance;
}
double Azimuth(struct site source, struct site destination)
{
/* This function returns the azimuth (in degrees) to the
destination as seen from the location of the source. */
double dest_lat, dest_lon, src_lat, src_lon,
beta, azimuth, diff, num, den, fraction;
dest_lat = destination.lat * DEG2RAD;
dest_lon = destination.lon * DEG2RAD;
src_lat = source.lat * DEG2RAD;
src_lon = source.lon * DEG2RAD;
/* Calculate Surface Distance */
beta =
acos(sin(src_lat) * sin(dest_lat) +
cos(src_lat) * cos(dest_lat) * cos(src_lon - dest_lon));
/* Calculate Azimuth */
num = sin(dest_lat) - (sin(src_lat) * cos(beta));
den = cos(src_lat) * sin(beta);
fraction = num / den;
/* Trap potential problems in acos() due to rounding */
if (fraction >= 1.0)
fraction = 1.0;
if (fraction <= -1.0)
fraction = -1.0;
/* Calculate azimuth */
azimuth = acos(fraction);
/* Reference it to True North */
diff = dest_lon - src_lon;
if (diff <= -PI)
diff += TWOPI;
if (diff >= PI)
diff -= TWOPI;
if (diff > 0.0)
azimuth = TWOPI - azimuth;
return (azimuth / DEG2RAD);
}
double ElevationAngle(struct site source, struct site destination)
{
/* This function returns the angle of elevation (in degrees)
of the destination as seen from the source location.
A positive result represents an angle of elevation (uptilt),
while a negative result represents an angle of depression
(downtilt), as referenced to a normal to the center of
the earth. */
register double a, b, dx;
a = GetElevation(destination) + destination.alt + earthradius;
b = GetElevation(source) + source.alt + earthradius;
dx = FEET_PER_MILE * Distance(source, destination);
/* Apply the Law of Cosines */
return ((180.0 *
(acos(((b * b) + (dx * dx) - (a * a)) / (2.0 * b * dx))) /
PI) - 90.0);
}
void ReadPath(struct site source, struct site destination)
{
/* This function generates a sequence of latitude and
longitude positions between source and destination
locations along a great circle path, and stores
elevation and distance information for points
along that path in the "path" structure. */
int c;
double azimuth, distance, lat1, lon1, beta, den, num,
lat2, lon2, total_distance, dx, dy, path_length,
miles_per_sample, samples_per_radian = 68755.0;
struct site tempsite;
lat1 = source.lat * DEG2RAD;
lon1 = source.lon * DEG2RAD;
lat2 = destination.lat * DEG2RAD;
lon2 = destination.lon * DEG2RAD;
samples_per_radian = ppd * 57.295833;
azimuth = Azimuth(source, destination) * DEG2RAD;
total_distance = Distance(source, destination);
if (total_distance > (30.0 / ppd)) {
dx = samples_per_radian * acos(cos(lon1 - lon2));
dy = samples_per_radian * acos(cos(lat1 - lat2));
path_length = sqrt((dx * dx) + (dy * dy));
miles_per_sample = total_distance / path_length;
}
else {
c = 0;
dx = 0.0;
dy = 0.0;
path_length = 0.0;
miles_per_sample = 0.0;
total_distance = 0.0;
lat1 = lat1 / DEG2RAD;
lon1 = lon1 / DEG2RAD;
path.lat[c] = lat1;
path.lon[c] = lon1;
path.elevation[c] = GetElevation(source);
path.distance[c] = 0.0;
}
for (distance = 0.0, c = 0;
(total_distance != 0.0 && distance <= total_distance
&& c < ARRAYSIZE); c++, distance = miles_per_sample * (double)c) {
beta = distance / 3959.0;
lat2 =
asin(sin(lat1) * cos(beta) +
cos(azimuth) * sin(beta) * cos(lat1));
num = cos(beta) - (sin(lat1) * sin(lat2));
den = cos(lat1) * cos(lat2);
if (azimuth == 0.0 && (beta > HALFPI - lat1))
lon2 = lon1 + PI;
else if (azimuth == HALFPI && (beta > HALFPI + lat1))
lon2 = lon1 + PI;
else if (fabs(num / den) > 1.0)
lon2 = lon1;
else {
if ((PI - azimuth) >= 0.0)
lon2 = lon1 - arccos(num, den);
else
lon2 = lon1 + arccos(num, den);
}
while (lon2 < 0.0)
lon2 += TWOPI;
while (lon2 > TWOPI)
lon2 -= TWOPI;
lat2 = lat2 / DEG2RAD;
lon2 = lon2 / DEG2RAD;
path.lat[c] = lat2;
path.lon[c] = lon2;
tempsite.lat = lat2;
tempsite.lon = lon2;
path.elevation[c] = GetElevation(tempsite);
// fix for tile gaps in multi-tile LIDAR plots
if(path.elevation[c]==0 && path.elevation[c-1] > 10)
path.elevation[c]=path.elevation[c-1];
path.distance[c] = distance;
}
/* Make sure exact destination point is recorded at path.length-1 */
if (c < ARRAYSIZE) {
path.lat[c] = destination.lat;
path.lon[c] = destination.lon;
path.elevation[c] = GetElevation(destination);
path.distance[c] = total_distance;
c++;
}
if (c < ARRAYSIZE)
path.length = c;
else
path.length = ARRAYSIZE - 1;
}
double ElevationAngle2(struct site source, struct site destination, double er)
{
/* This function returns the angle of elevation (in degrees)
of the destination as seen from the source location, UNLESS
the path between the sites is obstructed, in which case, the
elevation angle to the first obstruction is returned instead.
"er" represents the earth radius. */
int x;
char block = 0;
double source_alt, destination_alt, cos_xmtr_angle,
cos_test_angle, test_alt, elevation, distance,
source_alt2, first_obstruction_angle = 0.0;
struct path temp;
temp = path;
ReadPath(source, destination);
distance = FEET_PER_MILE * Distance(source, destination);
source_alt = er + source.alt + GetElevation(source);
destination_alt = er + destination.alt + GetElevation(destination);
source_alt2 = source_alt * source_alt;
/* Calculate the cosine of the elevation angle of the
destination (receiver) as seen by the source (transmitter). */
cos_xmtr_angle =
((source_alt2) + (distance * distance) -
(destination_alt * destination_alt)) / (2.0 * source_alt *
distance);
/* Test all points in between source and destination locations to
see if the angle to a topographic feature generates a higher
elevation angle than that produced by the destination. Begin
at the source since we're interested in identifying the FIRST
obstruction along the path between source and destination. */
for (x = 2, block = 0; x < path.length && block == 0; x++) {
distance = FEET_PER_MILE * path.distance[x];
test_alt =
earthradius + (path.elevation[x] ==
0.0 ? path.elevation[x] : path.elevation[x] +
clutter);
cos_test_angle =
((source_alt2) + (distance * distance) -
(test_alt * test_alt)) / (2.0 * source_alt * distance);
/* Compare these two angles to determine if
an obstruction exists. Since we're comparing
the cosines of these angles rather than
the angles themselves, the sense of the
following "if" statement is reversed from
what it would be if the angles themselves
were compared. */
if (cos_xmtr_angle >= cos_test_angle) {
block = 1;
first_obstruction_angle =
((acos(cos_test_angle)) / DEG2RAD) - 90.0;
}
}
if (block)
elevation = first_obstruction_angle;
else
elevation = ((acos(cos_xmtr_angle)) / DEG2RAD) - 90.0;
path = temp;
return elevation;
}
double ReadBearing(char *input)
{
/* This function takes numeric input in the form of a character
string, and returns an equivalent bearing in degrees as a
decimal number (double). The input may either be expressed
in decimal format (40.139722) or degree, minute, second
format (40 08 23). This function also safely handles
extra spaces found either leading, trailing, or
embedded within the numbers expressed in the
input string. Decimal seconds are permitted. */
double seconds, bearing = 0.0;
char string[20];
int a, b, length, degrees, minutes;
/* Copy "input" to "string", and ignore any extra
spaces that might be present in the process. */
string[0] = 0;
length = strlen(input);
for (a = 0, b = 0; a < length && a < 18; a++) {
if ((input[a] != 32 && input[a] != '\n')
|| (input[a] == 32 && input[a + 1] != 32
&& input[a + 1] != '\n' && b != 0)) {
string[b] = input[a];
b++;
}
}
string[b] = 0;
/* Count number of spaces in the clean string. */
length = strlen(string);
for (a = 0, b = 0; a < length; a++)
if (string[a] == 32)
b++;
if (b == 0) /* Decimal Format (40.139722) */
sscanf(string, "%lf", &bearing);
if (b == 2) { /* Degree, Minute, Second Format (40 08 23.xx) */
sscanf(string, "%d %d %lf", &degrees, &minutes, &seconds);
bearing = fabs((double)degrees);
bearing += fabs(((double)minutes) / 60.0);
bearing += fabs(seconds / 3600.0);
if ((degrees < 0) || (minutes < 0) || (seconds < 0.0))
bearing = -bearing;
}
/* Anything else returns a 0.0 */
if (bearing > 360.0 || bearing < -360.0)
bearing = 0.0;
return bearing;
}
void ObstructionAnalysis(struct site xmtr, struct site rcvr, double f,
FILE *outfile)
{
/* Perform an obstruction analysis along the
path between receiver and transmitter. */
int x;
struct site site_x;
double h_r, h_t, h_x, h_r_orig, cos_tx_angle, cos_test_angle,
cos_tx_angle_f1, cos_tx_angle_fpt6, d_tx, d_x,
h_r_f1, h_r_fpt6, h_f, h_los, lambda = 0.0;
char string[255], string_fpt6[255], string_f1[255];
ReadPath(xmtr, rcvr);
h_r = GetElevation(rcvr) + rcvr.alt + earthradius;
h_r_f1 = h_r;
h_r_fpt6 = h_r;
h_r_orig = h_r;
h_t = GetElevation(xmtr) + xmtr.alt + earthradius;
d_tx = FEET_PER_MILE * Distance(rcvr, xmtr);
cos_tx_angle =
((h_r * h_r) + (d_tx * d_tx) - (h_t * h_t)) / (2.0 * h_r * d_tx);
cos_tx_angle_f1 = cos_tx_angle;
cos_tx_angle_fpt6 = cos_tx_angle;
if (f)
lambda = 9.8425e8 / (f * 1e6);
if (clutter > 0.0) {
fprintf(outfile, "Terrain has been raised by");
if (metric)
fprintf(outfile, " %.2f meters",
METERS_PER_FOOT * clutter);
else
fprintf(outfile, " %.2f feet", clutter);
fprintf(outfile, " to account for ground clutter.\n\n");
}
/* At each point along the path calculate the cosine
of a sort of "inverse elevation angle" at the receiver.
From the antenna, 0 deg. looks at the ground, and 90 deg.
is parallel to the ground.
Start at the receiver. If this is the lowest antenna,
then terrain obstructions will be nearest to it. (Plus,
that's the way ppa!'s original los() did it.)
Calculate cosines only. That's sufficient to compare
angles and it saves the extra computational burden of
acos(). However, note the inverted comparison: if
acos(A) > acos(B), then B > A. */
for (x = path.length - 1; x > 0; x--) {
site_x.lat = path.lat[x];
site_x.lon = path.lon[x];
site_x.alt = 0.0;
h_x = GetElevation(site_x) + earthradius + clutter;
d_x = FEET_PER_MILE * Distance(rcvr, site_x);
/* Deal with the LOS path first. */
cos_test_angle =
((h_r * h_r) + (d_x * d_x) -
(h_x * h_x)) / (2.0 * h_r * d_x);
if (cos_tx_angle > cos_test_angle) {
if (h_r == h_r_orig)
fprintf(outfile,
"Between %s and %s, obstructions were detected at:\n\n",
rcvr.name, xmtr.name);
if (site_x.lat >= 0.0) {
if (metric)
fprintf(outfile,
" %8.4f N,%9.4f W, %5.2f kilometers, %6.2f meters AMSL\n",
site_x.lat, site_x.lon,
KM_PER_MILE * (d_x / FEET_PER_MILE),
METERS_PER_FOOT * (h_x -
earthradius));
else
fprintf(outfile,
" %8.4f N,%9.4f W, %5.2f miles, %6.2f feet AMSL\n",
site_x.lat, site_x.lon,
d_x / FEET_PER_MILE,
h_x - earthradius);
}
else {
if (metric)
fprintf(outfile,
" %8.4f S,%9.4f W, %5.2f kilometers, %6.2f meters AMSL\n",
-site_x.lat, site_x.lon,
KM_PER_MILE * (d_x / FEET_PER_MILE),
METERS_PER_FOOT * (h_x -
earthradius));
else
fprintf(outfile,
" %8.4f S,%9.4f W, %5.2f miles, %6.2f feet AMSL\n",
-site_x.lat, site_x.lon,
d_x / FEET_PER_MILE,
h_x - earthradius);
}
}
while (cos_tx_angle > cos_test_angle) {
h_r += 1;
cos_test_angle =
((h_r * h_r) + (d_x * d_x) -
(h_x * h_x)) / (2.0 * h_r * d_x);
cos_tx_angle =
((h_r * h_r) + (d_tx * d_tx) -
(h_t * h_t)) / (2.0 * h_r * d_tx);
}
if (f) {
/* Now clear the first Fresnel zone... */
cos_tx_angle_f1 =
((h_r_f1 * h_r_f1) + (d_tx * d_tx) -
(h_t * h_t)) / (2.0 * h_r_f1 * d_tx);
h_los =
sqrt(h_r_f1 * h_r_f1 + d_x * d_x -
2 * h_r_f1 * d_x * cos_tx_angle_f1);
h_f = h_los - sqrt(lambda * d_x * (d_tx - d_x) / d_tx);
while (h_f < h_x) {
h_r_f1 += 1;
cos_tx_angle_f1 =
((h_r_f1 * h_r_f1) + (d_tx * d_tx) -
(h_t * h_t)) / (2.0 * h_r_f1 * d_tx);
h_los =
sqrt(h_r_f1 * h_r_f1 + d_x * d_x -
2 * h_r_f1 * d_x * cos_tx_angle_f1);
h_f =
h_los -
sqrt(lambda * d_x * (d_tx - d_x) / d_tx);
}
/* and clear the 60% F1 zone. */
cos_tx_angle_fpt6 =
((h_r_fpt6 * h_r_fpt6) + (d_tx * d_tx) -
(h_t * h_t)) / (2.0 * h_r_fpt6 * d_tx);
h_los =
sqrt(h_r_fpt6 * h_r_fpt6 + d_x * d_x -
2 * h_r_fpt6 * d_x * cos_tx_angle_fpt6);
h_f =
h_los -
fzone_clearance * sqrt(lambda * d_x * (d_tx - d_x) /
d_tx);
while (h_f < h_x) {
h_r_fpt6 += 1;
cos_tx_angle_fpt6 =
((h_r_fpt6 * h_r_fpt6) + (d_tx * d_tx) -
(h_t * h_t)) / (2.0 * h_r_fpt6 * d_tx);
h_los =
sqrt(h_r_fpt6 * h_r_fpt6 + d_x * d_x -
2 * h_r_fpt6 * d_x *
cos_tx_angle_fpt6);
h_f =
h_los -
fzone_clearance * sqrt(lambda * d_x *
(d_tx - d_x) / d_tx);
}
}
}
if (h_r > h_r_orig) {
if (metric)
snprintf(string, 150,
"\nAntenna at %s must be raised to at least %.2f meters AGL\nto clear all obstructions detected.\n",
rcvr.name,
METERS_PER_FOOT * (h_r - GetElevation(rcvr) -
earthradius));
else
snprintf(string, 150,
"\nAntenna at %s must be raised to at least %.2f feet AGL\nto clear all obstructions detected.\n",
rcvr.name,
h_r - GetElevation(rcvr) - earthradius);
}
else
snprintf(string, 150,
"\nNo obstructions to LOS path due to terrain were detected\n");
if (f) {
if (h_r_fpt6 > h_r_orig) {
if (metric)
snprintf(string_fpt6, 150,
"\nAntenna at %s must be raised to at least %.2f meters AGL\nto clear %.0f%c of the first Fresnel zone.\n",
rcvr.name,
METERS_PER_FOOT * (h_r_fpt6 -
GetElevation(rcvr) -
earthradius),
fzone_clearance * 100.0, 37);
else
snprintf(string_fpt6, 150,
"\nAntenna at %s must be raised to at least %.2f feet AGL\nto clear %.0f%c of the first Fresnel zone.\n",
rcvr.name,
h_r_fpt6 - GetElevation(rcvr) -
earthradius, fzone_clearance * 100.0,
37);
}
else
snprintf(string_fpt6, 150,
"\n%.0f%c of the first Fresnel zone is clear.\n",
fzone_clearance * 100.0, 37);
if (h_r_f1 > h_r_orig) {
if (metric)
snprintf(string_f1, 150,
"\nAntenna at %s must be raised to at least %.2f meters AGL\nto clear the first Fresnel zone.\n",
rcvr.name,
METERS_PER_FOOT * (h_r_f1 -
GetElevation(rcvr) -
earthradius));
else
snprintf(string_f1, 150,
"\nAntenna at %s must be raised to at least %.2f feet AGL\nto clear the first Fresnel zone.\n",
rcvr.name,
h_r_f1 - GetElevation(rcvr) -
earthradius);
}
else
snprintf(string_f1, 150,
"\nThe first Fresnel zone is clear.\n");
}
fprintf(outfile, "%s", string);
if (f) {
fprintf(outfile, "%s", string_f1);
fprintf(outfile, "%s", string_fpt6);
}
}
static void free_dem(void)
{
int i;
int j;
for (i = 0; i < MAXPAGES; i++) {
for (j = 0; j < IPPD; j++) {
delete [] dem[i].data[j];
delete [] dem[i].mask[j];
delete [] dem[i].signal[j];
}
delete [] dem[i].data;
delete [] dem[i].mask;
delete [] dem[i].signal;
}
delete [] dem;
}
void free_elev(void) {
delete [] elev;
}
void free_path(void)
{
delete [] path.lat;
delete [] path.lon;
delete [] path.elevation;
delete [] path.distance;
}
void alloc_elev(void)
{
elev = new double[ARRAYSIZE + 10];
}
static void alloc_dem(void)
{
int i;
int j;
dem = new struct dem[MAXPAGES];
for (i = 0; i < MAXPAGES; i++) {
dem[i].data = new short *[IPPD];
dem[i].mask = new unsigned char *[IPPD];
dem[i].signal = new unsigned char *[IPPD];
for (j = 0; j < IPPD; j++) {
dem[i].data[j] = new short[IPPD];
dem[i].mask[j] = new unsigned char[IPPD];
dem[i].signal[j] = new unsigned char[IPPD];
}
}
}
void alloc_path(void)
{
path.lat = new double[ARRAYSIZE];
path.lon = new double[ARRAYSIZE];
path.elevation = new double[ARRAYSIZE];
path.distance = new double[ARRAYSIZE];
}
void do_allocs(void)
{
int i;
alloc_elev();
alloc_dem();
alloc_path();
for (i = 0; i < MAXPAGES; i++) {
dem[i].min_el = 32768;
dem[i].max_el = -32768;
dem[i].min_north = 90;
dem[i].max_north = -90;
dem[i].min_west = 360;
dem[i].max_west = -1;
}
}
int main(int argc, char *argv[])
{
int x, y, z = 0, min_lat, min_lon, max_lat, max_lon,
rxlat, rxlon, txlat, txlon, west_min, west_max,
nortRxHin, nortRxHax, propmodel, knifeedge = 0, ppa =
0, normalise = 0, haf = 0, pmenv = 1, lidar=0, cropped, result;
bool use_threads = true;
unsigned char LRmap = 0, txsites = 0, topomap = 0, geo = 0, kml =
0, area_mode = 0, max_txsites, ngs = 0;
char mapfile[255], ano_filename[255], lidar_tiles[27000], clutter_file[255];
char *az_filename, *el_filename, *udt_file = NULL;
double altitude = 0.0, altitudeLR = 0.0, tx_range = 0.0,
rx_range = 0.0, deg_range = 0.0, deg_limit = 0.0, deg_range_lon;
if (strstr(argv[0], "signalserverHD")) {
MAXPAGES = 9;
ARRAYSIZE = 32410;
IPPD = 3600;
}
if (strstr(argv[0], "signalserverLIDAR")) {
MAXPAGES = 100; // 10x10
lidar = 1;
IPPD = 6000; // will be overridden based upon file header...
}
strncpy(ss_name, "Signal Server\0", 14);
if (argc == 1) {
fprintf(stdout, "Version: %s %.2f (Built for %d DEM tiles at %d pixels)\n", ss_name, version,MAXPAGES, IPPD);
fprintf(stdout, "License: GNU General Public License (GPL) version 2\n\n");
fprintf(stdout, "Radio propagation simulator by Alex Farrant QCVS, 2E0TDW\n");
fprintf(stdout, "Based upon SPLAT! by John Magliacane, KD2BD\n\n");
fprintf(stdout, "Usage: signalserver [data options] [input options] [output options] -o outputfile\n\n");
fprintf(stdout, "Data:\n");
fprintf(stdout, " -sdf Directory containing SRTM derived .sdf DEM tiles\n");
fprintf(stdout, " -lid ASCII grid tile (LIDAR) with dimensions and resolution defined in header\n");
fprintf(stdout, " -udt User defined point clutter as decimal co-ordinates: 'latitude,longitude,height'\n");
fprintf(stdout, " -clt MODIS 17-class wide area clutter in ASCII grid format\n");
fprintf(stdout, "Input:\n");
fprintf(stdout, " -lat Tx Latitude (decimal degrees) -70/+70\n");
fprintf(stdout, " -lon Tx Longitude (decimal degrees) -180/+180\n");
fprintf(stdout, " -txh Tx Height (above ground)\n");
fprintf(stdout, " -rla (Optional) Rx Latitude for PPA (decimal degrees) -70/+70\n");
fprintf(stdout, " -rlo (Optional) Rx Longitude for PPA (decimal degrees) -180/+180\n");
fprintf(stdout, " -f Tx Frequency (MHz) 20MHz to 100GHz (LOS after 20GHz)\n");
fprintf(stdout, " -erp Tx Total Effective Radiated Power in Watts (dBd) inc Tx+Rx gain. 2.14dBi = 0dBd\n");
fprintf(stdout, " -rxh Rx Height(s) (optional. Default=0.1)\n");
fprintf(stdout, " -rxg Rx gain dBd (optional for PPA text report)\n");
fprintf(stdout, " -hp Horizontal Polarisation (default=vertical)\n");
fprintf(stdout, " -gc Random ground clutter (feet/meters)\n");
fprintf(stdout, " -m Metric units of measurement\n");
fprintf(stdout, " -te Terrain code 1-6 (optional)\n");
fprintf(stdout, " -terdic Terrain dielectric value 2-80 (optional)\n");
fprintf(stdout, " -tercon Terrain conductivity 0.01-0.0001 (optional)\n");
fprintf(stdout, " -cl Climate code 1-6 (optional)\n");
fprintf(stdout, " -rel Reliability for ITM model 50 to 99 (optional)\n");
fprintf(stdout, " -resample Reduce Lidar resolution by specified factor (2 = 50%)\n");
fprintf(stdout, "Output:\n");
fprintf(stdout, " -dbm Plot Rxd signal power instead of field strength\n");
fprintf(stdout, " -rt Rx Threshold (dB / dBm / dBuV/m)\n");
fprintf(stdout, " -o Filename. Required. \n");
fprintf(stdout, " -R Radius (miles/kilometers)\n");
fprintf(stdout, " -res Pixels per tile. 300/600/1200/3600 (Optional. LIDAR res is within the tile)\n");
fprintf(stdout, " -pm Propagation model. 1: ITM, 2: LOS, 3: Hata, 4: ECC33,\n");
fprintf(stdout, " 5: SUI, 6: COST-Hata, 7: FSPL, 8: ITWOM, 9: Ericsson,\n");
fprintf(stdout, " 10: Plane earth, 11: Egli VHF/UHF, 12: Soil\n");
fprintf(stdout, " -pe Propagation model mode: 1=Urban,2=Suburban,3=Rural\n");
fprintf(stdout, " -ked Knife edge diffraction (Already on for ITM)\n");
fprintf(stdout, "Debugging:\n");
fprintf(stdout, " -t Terrain greyscale background\n");
fprintf(stdout, " -dbg Verbose debug messages\n");
fprintf(stdout, " -ng Normalise Path Profile graph\n");
fprintf(stdout, " -haf Halve 1 or 2 (optional)\n");
fprintf(stdout, " -nothreads Turn off threaded processing\n");
fflush(stdout);
return 1;
}
/*
* If we're not called as signalserverLIDAR we can allocate various
* memory now. For LIDAR we need to wait until we've parsed
* the headers in the .asc file to know how much memory to allocate...
*/
if (!lidar)
do_allocs();
y = argc - 1;
kml = 0;
geo = 0;
dbm = 0;
gpsav = 0;
metric = 0;
string[0] = 0;
mapfile[0] = 0;
clutter_file[0] = 0;
clutter = 0.0;
forced_erp = -1.0;
forced_freq = 0.0;
sdf_path[0] = 0;
udt_file = NULL;
path.length = 0;
max_txsites = 30;
fzone_clearance = 0.6;
contour_threshold = 0;
resample = 0;
ano_filename[0] = 0;
earthradius = EARTHRADIUS;
max_range = 1.0;
propmodel = 1; //ITM
lat = 0;
lon = 0;
txh = 0;
ngs = 1; // no terrain background
kml = 1;
LRmap = 1;
area_mode = 1;
ippd = IPPD; // default resolution
sscanf("0.1", "%lf", &altitudeLR);
// Defaults
LR.eps_dielect = 15.0; // Farmland
LR.sgm_conductivity = 0.005; // Farmland
LR.eno_ns_surfref = 301.0;
LR.frq_mhz = 19.0; // Deliberately too low
LR.radio_climate = 5; // continental
LR.pol = 1; // vert
LR.conf = 0.50;
LR.rel = 0.50;
LR.erp = 0.0; // will default to Path Loss
tx_site[0].lat = 91.0;
tx_site[0].lon = 361.0;
tx_site[1].lat = 91.0;
tx_site[1].lon = 361.0;
/* Scan for command line arguments */
for (x = 1; x <= y; x++) {
if (strcmp(argv[x], "-R") == 0) {
z = x + 1;
if (z <= y && argv[z][0] && argv[z][0] != '-') {
sscanf(argv[z], "%lf", &max_range);
}
}
if (strcmp(argv[x], "-gc") == 0) {
z = x + 1;
if (z <= y && argv[z][0] && argv[z][0] != '-') {
sscanf(argv[z], "%lf", &clutter);
if (clutter < 0.0)
clutter = 0.0;
}
}
if (strcmp(argv[x], "-clt") == 0) {
z = x + 1;
if (z <= y && argv[z][0] && argv[z][0] != '-') {
strncpy(clutter_file, argv[z], 253);
}
}
if (strcmp(argv[x], "-o") == 0) {
z = x + 1;
if (z <= y && argv[z][0] && argv[z][0] != '-') {
strncpy(mapfile, argv[z], 253);
strncpy(tx_site[0].name, "Tx", 2);
strncpy(tx_site[0].filename, argv[z], 253);
/* Antenna pattern files have the same basic name as the output file
* but with a different extension. If they exist, load them now */
if( (az_filename = (char*) calloc(strlen(argv[z]) + strlen(AZ_FILE_SUFFIX) + 1, sizeof(char))) == NULL )
return ENOMEM;
strcpy(az_filename, argv[z]);
strcat(az_filename, AZ_FILE_SUFFIX);
if( (el_filename = (char*) calloc(strlen(argv[z]) + strlen(EL_FILE_SUFFIX) + 1, sizeof(char))) == NULL ){
free(az_filename);
return ENOMEM;
}
strcpy(el_filename, argv[z]);
strcat(el_filename, EL_FILE_SUFFIX);
if( (result = LoadPAT(az_filename,el_filename)) != 0 ){
fprintf(stderr,"Permissions error reading antenna pattern file\n");
free(az_filename);
free(el_filename);
exit(result);
}
free(az_filename);
free(el_filename);
} else if (z <= y && argv[z][0] && argv[z][0] == '-' && argv[z][1] == '\0' ) {
/* Handle writing image data to stdout */
to_stdout = true;
mapfile[0] = '\0';
strncpy(tx_site[0].name, "Tx", 2);
tx_site[0].filename[0] = '\0';
fprintf(stderr,"Writing to stdout\n");
}
}
if (strcmp(argv[x], "-so") == 0) {
z = x + 1;
if(image_set_library(argv[z]) != 0){
fprintf(stderr,"Error configuring image processor\n");
exit(EINVAL);
}
}
if (strcmp(argv[x], "-rt") == 0) {
z = x + 1;
if (z <= y && argv[z][0]) /* A minus argument is legal here */
sscanf(argv[z], "%d", &contour_threshold);
}
if (strcmp(argv[x], "-m") == 0) {
metric = 1;
}
if (strcmp(argv[x], "-t") == 0) {
ngs = 0; // greyscale background
}
if (strcmp(argv[x], "-dbm") == 0)
dbm = 1;
if (strcmp(argv[x], "-sdf") == 0) {
z = x + 1;
if (z <= y && argv[z][0] && argv[z][0] != '-')
strncpy(sdf_path, argv[z], 253);
}
if (strcmp(argv[x], "-lid") == 0) {
z = x + 1;
lidar=1;
if (z <= y && argv[z][0] && argv[z][0] != '-')
strncpy(lidar_tiles, argv[z], 27000); // 900 tiles!
}
if (strcmp(argv[x], "-res") == 0) {
z = x + 1;
if (!lidar &&
z <= y &&
argv[z][0] &&
argv[z][0] != '-') {
sscanf(argv[z], "%d", &ippd);
switch (ippd) {
case 300:
MAXRAD = 500;
jgets = 3; // 3 dummy reads
break;
case 600:
MAXRAD = 500;
jgets = 1;
break;
case 1200:
MAXRAD = 200;
ippd = 1200;
break;
case 3600:
MAXRAD = 100;
ippd = 3600;
break;
default:
MAXRAD = 200;
ippd = 1200;
break;
}
}
}
if (strcmp(argv[x], "-resample") == 0) {
z = x + 1;
if(!lidar){
fprintf(stderr, "Error, this should only be used with LIDAR tiles.\n");
return -1;
}
sscanf(argv[z], "%d", &resample);
}
if (strcmp(argv[x], "-lat") == 0) {
z = x + 1;
if (z <= y && argv[z][0]) {
tx_site[0].lat = ReadBearing(argv[z]);
}
}
if (strcmp(argv[x], "-lon") == 0) {
z = x + 1;
if (z <= y && argv[z][0]) {
tx_site[0].lon = ReadBearing(argv[z]);
tx_site[0].lon *= -1;
if (tx_site[0].lon < 0.0)
tx_site[0].lon += 360.0;
}
}
//Switch to Path Profile Mode if Rx co-ords specified
if (strcmp(argv[x], "-rla") == 0) {
z = x + 1;
if (z <= y && argv[z][0]) {
ppa = 1;
tx_site[1].lat = ReadBearing(argv[z]);
}
}
if (strcmp(argv[x], "-rlo") == 0) {
z = x + 1;
if (z <= y && argv[z][0]) {
tx_site[1].lon = ReadBearing(argv[z]);
tx_site[1].lon *= -1;
if (tx_site[1].lon < 0.0)
tx_site[1].lon += 360.0;
}
}
if (strcmp(argv[x], "-txh") == 0) {
z = x + 1;
if (z <= y && argv[z][0] && argv[z][0] != '-') {
sscanf(argv[z], "%f", &tx_site[0].alt);
}
txsites = 1;
}
if (strcmp(argv[x], "-rxh") == 0) {
z = x + 1;
if (z <= y && argv[z][0] && argv[z][0] != '-') {
sscanf(argv[z], "%lf", &altitudeLR);
sscanf(argv[z], "%f", &tx_site[1].alt);
}
}
if (strcmp(argv[x], "-rxg") == 0) {
z = x + 1;
if (z <= y && argv[z][0] && argv[z][0] != '-') {
sscanf(argv[z], "%lf", &rxGain);
}
}
if (strcmp(argv[x], "-f") == 0) {
z = x + 1;
if (z <= y && argv[z][0] && argv[z][0] != '-') {
sscanf(argv[z], "%lf", &LR.frq_mhz);
}
}
if (strcmp(argv[x], "-erp") == 0) {
z = x + 1;
if (z <= y && argv[z][0] && argv[z][0] != '-') {
sscanf(argv[z], "%lf", &LR.erp);
}
}
if (strcmp(argv[x], "-cl") == 0) {
z = x + 1;
if (z <= y && argv[z][0] && argv[z][0] != '-') {
sscanf(argv[z], "%d", &LR.radio_climate);
}
}
if (strcmp(argv[x], "-te") == 0) {
z = x + 1;
if (z <= y && argv[z][0] && argv[z][0] != '-') {
sscanf(argv[z], "%d", &ter);
switch (ter) {
case 1: // Water
terdic = 80;
tercon = 0.010;
break;
case 2: // Marsh
terdic = 12;
tercon = 0.007;
break;
case 3: // Farmland
terdic = 15;
tercon = 0.005;
break;
case 4: //Mountain
terdic = 13;
tercon = 0.002;
break;
case 5: //Desert
terdic = 13;
tercon = 0.002;
break;
case 6: //Urban
terdic = 5;
tercon = 0.001;
break;
}
LR.eps_dielect = terdic;
LR.sgm_conductivity = tercon;
}
}
if (strcmp(argv[x], "-terdic") == 0) {
z = x + 1;
if (z <= y && argv[z][0] && argv[z][0] != '-') {
sscanf(argv[z], "%lf", &terdic);
LR.eps_dielect = terdic;
}
}
if (strcmp(argv[x], "-tercon") == 0) {
z = x + 1;
if (z <= y && argv[z][0] && argv[z][0] != '-') {
sscanf(argv[z], "%lf", &tercon);
LR.sgm_conductivity = tercon;
}
}
if (strcmp(argv[x], "-hp") == 0) {
// Horizontal polarisation (0)
LR.pol = 0;
}
if (strcmp(argv[x], "-dbg") == 0) {
debug = 1;
}
/*UDT*/
if (strcmp(argv[x], "-udt") == 0) {
z = x + 1;
if (z <= y && argv[z][0]) {
udt_file = (char*) calloc(PATH_MAX+1, sizeof(char));
if( udt_file == NULL )
return ENOMEM;
strncpy(udt_file, argv[z], 253);
}
}
/*Prop model */
if (strcmp(argv[x], "-pm") == 0) {
z = x + 1;
if (z <= y && argv[z][0]) {
sscanf(argv[z], "%d", &propmodel);
}
}
// Prop model variant eg. urban/suburban
if (strcmp(argv[x], "-pe") == 0) {
z = x + 1;
if (z <= y && argv[z][0]) {
sscanf(argv[z], "%d", &pmenv);
}
}
//Knife edge diffraction
if (strcmp(argv[x], "-ked") == 0) {
z = x + 1;
knifeedge = 1;
}
//Normalise Path Profile chart
if (strcmp(argv[x], "-ng") == 0) {
z = x + 1;
normalise = 1;
}
//Halve the problem
if (strcmp(argv[x], "-haf") == 0) {
z = x + 1;
if (z <= y && argv[z][0]) {
sscanf(argv[z], "%d", &haf);
}
}
//Disable threads
if (strcmp(argv[x], "-nothreads") == 0) {
z = x + 1;
use_threads = false;
}
// Reliability % for ITM model
if (strcmp(argv[x], "-rel") == 0) {
z = x + 1;
if (z <= y && argv[z][0]) {
sscanf(argv[z], "%lf", &LR.rel);
sscanf(argv[z], "%lf", &LR.conf);
LR.rel=LR.rel/100;
LR.conf=LR.conf/100;
}
}
}
/* ERROR DETECTION */
if (tx_site[0].lat > 90 || tx_site[0].lat < -90) {
fprintf(stderr,
"ERROR: Either the lat was missing or out of range!");
exit(EINVAL);
}
if (tx_site[0].lon > 360 || tx_site[0].lon < 0) {
fprintf(stderr,
"ERROR: Either the lon was missing or out of range!");
exit(EINVAL);
}
if (LR.frq_mhz < 20 || LR.frq_mhz > 100000) {
fprintf(stderr,
"ERROR: Either the Frequency was missing or out of range!");
exit(EINVAL);
}
if (LR.erp > 500000000) {
fprintf(stderr, "ERROR: Power was out of range!");
exit(EINVAL);
}
if (LR.eps_dielect > 80 || LR.eps_dielect < 0.1) {
fprintf(stderr, "ERROR: Ground Dielectric value out of range!");
exit(EINVAL);
}
if (LR.sgm_conductivity > 0.01 || LR.sgm_conductivity < 0.000001) {
fprintf(stderr, "ERROR: Ground conductivity out of range!");
exit(EINVAL);
}
if (tx_site[0].alt < 0 || tx_site[0].alt > 60000) {
fprintf(stderr,
"ERROR: Tx altitude above ground was too high: %f",
tx_site[0].alt);
exit(EINVAL);
}
if (altitudeLR < 0 || altitudeLR > 60000) {
fprintf(stderr,
"ERROR: Rx altitude above ground was too high!");
exit(EINVAL);
}
if(!lidar){
if (ippd < 300 || ippd > 10000) {
fprintf(stderr, "ERROR: resolution out of range!");
exit(EINVAL);
}
}
if (contour_threshold < -200 || contour_threshold > 240) {
fprintf(stderr,
"ERROR: Receiver threshold out of range (-200 / +240)");
exit(EINVAL);
}
if (propmodel > 2 && propmodel < 7 && LR.frq_mhz < 150) {
fprintf(stderr,
"ERROR: Frequency too low for Propagation model");
exit(EINVAL);
}
if (to_stdout == true && ppa != 0) {
fprintf(stderr,
"ERROR: Cannot write to stdout in ppa mode");
exit(EINVAL);
}
if(resample > 10){
fprintf(stderr,
"ERROR: Cannot resample higher than a factor of 10");
exit(EINVAL);
}
if (metric) {
altitudeLR /= METERS_PER_FOOT; /* 10ft * 0.3 = 3.3m */
max_range /= KM_PER_MILE; /* 10 / 1.6 = 7.5 */
altitude /= METERS_PER_FOOT;
tx_site[0].alt /= METERS_PER_FOOT; /* Feet to metres */
tx_site[1].alt /= METERS_PER_FOOT; /* Feet to metres */
clutter /= METERS_PER_FOOT; /* Feet to metres */
}
/* Ensure a trailing '/' is present in sdf_path */
if (sdf_path[0]) {
x = strlen(sdf_path);
if (sdf_path[x - 1] != '/' && x != 0) {
sdf_path[x] = '/';
sdf_path[x + 1] = 0;
}
}
x = 0;
y = 0;
min_lat = 70;
max_lat = -70;
min_lon = (int)floor(tx_site[0].lon);
max_lon = (int)floor(tx_site[0].lon);
txlat = (int)floor(tx_site[0].lat);
txlon = (int)floor(tx_site[0].lon);
if (txlat < min_lat)
min_lat = txlat;
if (txlat > max_lat)
max_lat = txlat;
if (LonDiff(txlon, min_lon) < 0.0)
min_lon = txlon;
if (LonDiff(txlon, max_lon) >= 0.0)
max_lon = txlon;
if (ppa == 1) {
rxlat = (int)floor(tx_site[1].lat);
rxlon = (int)floor(tx_site[1].lon);
if (rxlat < min_lat)
min_lat = rxlat;
if (rxlat > max_lat)
max_lat = rxlat;
if (LonDiff(rxlon, min_lon) < 0.0)
min_lon = rxlon;
if (LonDiff(rxlon, max_lon) >= 0.0)
max_lon = rxlon;
}
/* Load the required tiles */
if(lidar){
if( (result = loadLIDAR(lidar_tiles, resample)) != 0 ){
fprintf(stderr, "Couldn't find one or more of the "
"lidar files. Please ensure their paths are "
"correct and try again.\n");
fprintf(stderr, "Error %d: %s\n", result, strerror(result));
exit(result);
}
ppd=(double) (height / (max_north-min_north));
yppd=ppd;
if(debug){
fprintf(stderr,"ppd %lf, yppd %lf, %.4f,%.4f,%.4f,%.4f,%d x %d\n",ppd,yppd,max_north,min_west,min_north,max_west,width,height);
}
if(yppd<ppd/4){
fprintf(stderr,"yppd is bad! Check longitudes\n");
exit(1);
}
if(delta>0){
tx_site[0].lon+=delta;
}
}else{
// DEM first
if(debug){
fprintf(stderr,"%.4f,%.4f,%.4f,%.4f,%.4f,%.4f\n",max_north,min_west,min_north,max_west,max_lon,min_lon);
}
//max_lon-=3;
if( (result = LoadTopoData(max_lon, min_lon, max_lat, min_lat)) != 0 ){
// This only fails on errors loading SDF tiles
fprintf(stderr, "Error loading topo data\n");
return result;
}
if (area_mode || topomap) {
for (z = 0; z < txsites && z < max_txsites; z++) {
/* "Ball park" estimates used to load any additional
SDF files required to conduct this analysis. */
tx_range =
sqrt(1.5 *
(tx_site[z].alt + GetElevation(tx_site[z])));
if (LRmap)
rx_range = sqrt(1.5 * altitudeLR);
else
rx_range = sqrt(1.5 * altitude);
/* deg_range determines the maximum
amount of topo data we read */
deg_range = (tx_range + rx_range) / 57.0;
/* max_range regulates the size of the
analysis. A small, non-zero amount can
be used to shrink the size of the analysis
and limit the amount of topo data read by
ss A large number will increase the
width of the analysis and the size of
the map. */
if (max_range == 0.0)
max_range = tx_range + rx_range;
deg_range = max_range / 57.0;
// No more than 8 degs
deg_limit = 3.5;
if (fabs(tx_site[z].lat) < 70.0)
deg_range_lon =
deg_range / cos(DEG2RAD * tx_site[z].lat);
else
deg_range_lon = deg_range / cos(DEG2RAD * 70.0);
/* Correct for squares in degrees not being square in miles */
if (deg_range > deg_limit)
deg_range = deg_limit;
if (deg_range_lon > deg_limit)
deg_range_lon = deg_limit;
nortRxHin = (int)floor(tx_site[z].lat - deg_range);
nortRxHax = (int)floor(tx_site[z].lat + deg_range);
west_min = (int)floor(tx_site[z].lon - deg_range_lon);
while (west_min < 0)
west_min += 360;
while (west_min >= 360)
west_min -= 360;
west_max = (int)floor(tx_site[z].lon + deg_range_lon);
while (west_max < 0)
west_max += 360;
while (west_max >= 360)
west_max -= 360;
if (nortRxHin < min_lat)
min_lat = nortRxHin;
if (nortRxHax > max_lat)
max_lat = nortRxHax;
if (LonDiff(west_min, min_lon) < 0.0)
min_lon = west_min;
if (LonDiff(west_max, max_lon) >= 0.0)
max_lon = west_max;
}
/* Load any additional SDF files, if required */
if( (result = LoadTopoData(max_lon, min_lon, max_lat, min_lat)) != 0 ){
// This only fails on errors loading SDF tiles
fprintf(stderr, "Error loading topo data\n");
return result;
}
}
ppd=(double)ippd;
yppd=ppd;
width = (unsigned)(ippd * ReduceAngle(max_west - min_west));
height = (unsigned)(ippd * ReduceAngle(max_north - min_north));
}
dpp = 1 / ppd;
mpi = ippd-1;
// User defined clutter file
if( udt_file != NULL && (result = LoadUDT(udt_file)) != 0 ){
fprintf(stderr, "Error loading clutter file\n");
return result;
}
// Enrich with Clutter
if(strlen(clutter_file) > 1){
/*
Clutter tiles cover 16 x 12 degs but we only need a fraction of that area.
Limit by max_range / miles per degree (at equator)
*/
if( (result = loadClutter(clutter_file,max_range/45,tx_site[0])) != 0 ){
fprintf(stderr, "Error, invalid or clutter file not found\n");
return result;
}
}
if(max_range>100 || LR.frq_mhz==446.446){
cropping=false;
}
if (ppa == 0) {
if (propmodel == 2) {
cropping = false;
PlotLOSMap(tx_site[0], altitudeLR, ano_filename, use_threads);
DoLOS(mapfile, geo, kml, ngs, tx_site, txsites);
} else {
// 90% of effort here
PlotPropagation(tx_site[0], altitudeLR, ano_filename,
propmodel, knifeedge, haf, pmenv, use_threads);
if(debug)
fprintf(stderr,"Finished PlotPropagation()\n");
// nearfield void
for (float x=-0.001; x<0.001;x=x+0.0001){
for (float y=-0.001; y<0.001;y=y+0.0001){
if(GetSignal(tx_site[0].lat+y, tx_site[0].lon+x)<=0){
PutSignal(tx_site[0].lat+y, tx_site[0].lon+x, hottest);
}
}
}
if(cropping){
// CROPPING. Factor is determined in propPathLoss().
// cropLon is the circle radius in pixels at it's widest (east/west)
cropLon*=dpp; // pixels to degrees
max_north=cropLat; // degrees
max_west=cropLon+tx_site[0].lon; // degrees west (positive)
cropLat-=tx_site[0].lat; // angle from tx to edge
if(debug)
fprintf(stderr,"Cropping 1: max_west: %.4f cropLat: %.4f cropLon: %.4f longitude: %.5f dpp %.7f\n",max_west,cropLat,cropLon,tx_site[0].lon,dpp);
width=(int)((cropLon*ppd)*2);
height=(int)((cropLat*ppd)*2);
if(debug)
fprintf(stderr,"Cropping 2: max_west: %.4f cropLat: %.4f cropLon: %.7f longitude: %.5f width %d\n",max_west,cropLat,cropLon,tx_site[0].lon,width);
if(width>3600*10 || cropLon < 0){
fprintf(stderr,"FATAL BOUNDS! max_west: %.4f cropLat: %.4f cropLon: %.7f longitude: %.5f\n",max_west,cropLat,cropLon,tx_site[0].lon);
return 0;
}
}
// Write bitmap
if (LR.erp == 0.0)
DoPathLoss(mapfile, geo, kml, ngs, tx_site,
txsites);
else if (dbm)
DoRxdPwr((to_stdout == true ? NULL : mapfile), geo, kml, ngs, tx_site,
txsites);
else
if( (result = DoSigStr(mapfile, geo, kml, ngs, tx_site,txsites)) != 0 )
return result;
}
/*if(lidar){
east=eastoffset;
west=westoffset;
}*/
if (tx_site[0].lon > 0.0){
tx_site[0].lon *= -1;
}
if (tx_site[0].lon < -180.0){
tx_site[0].lon += 360;
}
if (cropping) {
fprintf(stderr, "|%.6f", tx_site[0].lat+cropLat);
fprintf(stderr, "|%.6f", tx_site[0].lon+cropLon);
fprintf(stderr, "|%.6f", tx_site[0].lat-cropLat);
fprintf(stderr, "|%.6f|",tx_site[0].lon-cropLon);
}else{
fprintf(stderr, "|%.6f", max_north);
fprintf(stderr, "|%.6f", east);
fprintf(stderr, "|%.6f", min_north);
fprintf(stderr, "|%.6f|",west);
}
fprintf(stderr, "\n");
} else {
strncpy(tx_site[0].name, "Tx", 3);
strncpy(tx_site[1].name, "Rx", 3);
PlotPath(tx_site[0], tx_site[1], 1);
PathReport(tx_site[0], tx_site[1], tx_site[0].filename, 0,
propmodel, pmenv, rxGain);
// Order flipped for benefit of graph. Makes no difference to data.
SeriesData(tx_site[1], tx_site[0], tx_site[0].filename, 1,
normalise);
}
fflush(stderr);
return 0;
}
View Git Blame Copy Permalink