3.10 More tests and 360 min_west fix

CHANGELOG

@@ -1,5 +1,8 @@
 SIGNAL SERVER CHANGELOG 
 
+3.10 - 15 Apr 2018
+360 min_west fix 
+
 3.09 - 7 Feb 2018
 Meridian LIDAR cropping fix
 

Makefile

@@ -8,7 +8,7 @@ LIBS		= -lm -lpthread -ldl
 
 VPATH		= models
 objects 	= main.o cost.o ecc33.o ericsson.o fspl.o hata.o itwom3.0.o \
-		  los.o sui.o pel.o egli.o inputs.o outputs.o image.o image-ppm.o tiles.o
+		  los.o sui.o pel.o egli.o soil.o inputs.o outputs.o image.o image-ppm.o tiles.o
 
 GCC_MAJOR	:= $(shell $(CXX) -dumpversion 2>&1 | cut -d . -f 1)
 GCC_MINOR	:= $(shell $(CXX) -dumpversion 2>&1 | cut -d . -f 2)
@@ -43,14 +43,14 @@ main.o: main.cc common.h inputs.hh outputs.hh itwom3.0.hh los.hh
 inputs.o: inputs.cc common.h main.hh
 
 outputs.o: outputs.cc common.h inputs.hh main.hh cost.hh ecc33.hh ericsson.hh \
-	   fspl.hh hata.hh itwom3.0.hh sui.hh pel.hh egli.hh
+	   fspl.hh hata.hh itwom3.0.hh sui.hh pel.hh egli.hh soil.hh
 
 image.o: image.cc image-ppm.o
 
 image-ppm.o: image-ppm.cc
 
 los.o: los.cc common.h main.hh cost.hh ecc33.hh ericsson.hh fspl.hh hata.hh \
-       itwom3.0.hh sui.hh pel.hh egli.hh
+       itwom3.0.hh sui.hh pel.hh egli.hh soil.hh
 
 tiles.o: tiles.cc tiles.hh common.h
 

inputs.cc

@@ -217,7 +217,7 @@ int loadLIDAR(char *filenames, int resample)
 
 		if (tiles[indx].max_west > max_west)
 			max_west = tiles[indx].max_west;
-	
+
 		if (tiles[indx].min_west < min_west)
 			min_west = tiles[indx].min_west;
 
@@ -225,10 +225,6 @@ int loadLIDAR(char *filenames, int resample)
 			eastF=1;
 		if (tiles[indx].max_west < 2.0)
 			westF=1;
-		if (tiles[indx].min_west > 359){
-			westF=0;
-			eastF=1;
-		}
 
 		if (max_west == -1) {
 			max_west = tiles[indx].max_west;
@@ -242,33 +238,40 @@ int loadLIDAR(char *filenames, int resample)
 			}
 		}
 
-		
 
-		if (fabs(tiles[indx].min_west - min_west) < 180.0) {
+		if (tiles[indx].min_west > 359.9) {
-			if (tiles[indx].min_west < min_west)
+			min_west = tiles[indx].min_west;
-				min_west = tiles[indx].min_west;
+			eastF=1;
+			if(debug)
+				fprintf(stderr,"min_west set to 360\n");
 		} else {
-			if (tiles[indx].min_west > min_west)
+			if (fabs(tiles[indx].min_west - min_west) < 180.0) {
-				min_west = tiles[indx].min_west;
+				if (tiles[indx].min_west < min_west)
+					min_west = tiles[indx].min_west;
+			} else {
+				if (tiles[indx].min_west > min_west)
+					min_west = tiles[indx].min_west;
+			}
 		}
 
-
 	}
 
 	// Meridian fix
 	if(eastF && westF){
 		max_west=0; //reset
 		for (indx = 0; indx < fc; indx++) {
-			if(tiles[indx].max_west<2.0 && tiles[indx].max_west > max_west)
+			if(tiles[indx].max_west<=2.0 && tiles[indx].max_west > max_west)
 				max_west=tiles[indx].max_west;
 		}
 		max_west*=1.0;  // WGS84 to westing. -1.5 = 1.5
-		
+
 	}
 	// -1 fix
-	if(eastF && fc>1 && min_west<=1.0)
+	if(eastF && fc>1 && min_west<=1.0001){
-		min_west=0;
+		min_west=0.0;
-	
+		if(debug)
+			fprintf(stderr,"min_west=0\n");
+	}
 	/* Iterate through all of the tiles to find the smallest resolution. We will
 	 * need to rescale every tile from here on out to this value */
 	float smallest_res = 0;
@@ -288,7 +291,7 @@ int loadLIDAR(char *filenames, int resample)
 
 	// Don't resize large 1 deg tiles in large multi-degree plots as it gets messy
 	if(tiles[0].width != 3600){
-	
+
 		for (size_t i = 0; i< fc; i++) {
 			float rescale = tiles[i].resolution / (float)desired_resolution;
 			if(debug)
@@ -306,7 +309,9 @@ int loadLIDAR(char *filenames, int resample)
 
 
 	/* Now we work out the size of the giant lidar tile. */
-	fprintf(stderr,"mw:%lf Mnw:%lf\n", max_west, min_west);
+	if(debug){
+		fprintf(stderr,"mw:%lf Mnw:%lf\n", max_west, min_west);
+	}
 	double total_width = max_west - min_west >= 0 ? max_west - min_west : max_west + (360 - min_west);
 	double total_height = max_north - min_north;
 	if (debug) {
@@ -317,7 +322,7 @@ int loadLIDAR(char *filenames, int resample)
 
 	//detect problematic layouts eg. vertical rectangles
 	// 1x2
-	if(fc >= 2 && desired_resolution < 15 && total_height > total_width*1.2){
+	if(fc >= 2 && desired_resolution < 28 && total_height > total_width*1.2){
 		tiles[fc].max_north=max_north;
 		tiles[fc].min_north=min_north;
 		westoffset=westoffset-(total_height-total_width); // WGS84 for stdout only
@@ -338,7 +343,7 @@ int loadLIDAR(char *filenames, int resample)
 		}	
 	}
 	// 2x1
-	if(fc >= 2 && desired_resolution < 15 && total_width > total_height*1.2){
+	if(fc >= 2 && desired_resolution < 28 && total_width > total_height*1.2){
 		tiles[fc].max_north=max_north+(total_width-total_height);
 		tiles[fc].min_north=max_north;
 		tiles[fc].max_west=max_west; // Positive westing
@@ -385,6 +390,11 @@ int loadLIDAR(char *filenames, int resample)
 	if (debug)
 		fprintf(stderr,"Lidar tile dimensions w:%lf(%zu) h:%lf(%zu)\n", total_width, new_width, total_height, new_height);
 
+	//sanity check!
+	if(new_width > 50e3 || new_height > 50e3){
+		fprintf(stdout,"Not processing a tile with these dimensions: %zu x %zu\n",new_width,new_height);
+		exit(1);
+	}
 	/* ...If we wanted a value other than sea level here, we would need to initialize the array... */
 	size_t prevPixelOffsetW=0;
 	size_t prevPixelOffsetN=0;

main.cc

@@ -1,4 +1,4 @@
-double version = 3.09;
+double version = 3.10;
 /****************************************************************************\
 *  Signal Server: Radio propagation simulator by Alex Farrant QCVS, 2E0TDW   *
 ******************************************************************************
@@ -1105,7 +1105,8 @@ int main(int argc, char *argv[])
 		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, 10: Plane earth, 11: Egli VHF/UHF\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");
@@ -1746,7 +1747,7 @@ int main(int argc, char *argv[])
 			fprintf(stderr,"%.4f,%.4f,%.4f,%.4f,%.4f,%.4f\n",max_north,min_west,min_north,max_west,max_lon,min_lon);
 		}
 
-		max_lon-=3;
+		//max_lon-=3;
 
 		if( (result = LoadTopoData(max_lon, min_lon, max_lat, min_lat)) != 0 ){
 			// This only fails on errors loading SDF tiles
@@ -1844,7 +1845,7 @@ int main(int argc, char *argv[])
 			}
 		}
 		ppd=(double)ippd;
-		//yppd=ppd; 
+		yppd=ppd; 
 
 		width = (unsigned)(ippd * ReduceAngle(max_west - min_west));
 		height = (unsigned)(ippd * ReduceAngle(max_north - min_north));

models/los.cc

@@ -11,6 +11,7 @@
 #include "sui.hh"
 #include "pel.hh"
 #include "egli.hh"
+#include "soil.hh"
 #include <pthread.h>
 
 #define NUM_SECTIONS 4
@@ -186,7 +187,7 @@ static double incidenceAngle(double opp, double adj)
  */
 static double ked(double freq, double rxh, double dkm)
 {
-	double obh, obd, rxobaoi = 0, d, dipheight = 25;
+	double obh, obd, rxobaoi = 0, d;
 
 	obh = 0;		// Obstacle height
 	obd = 0;		// Obstacle distance
@@ -199,7 +200,7 @@ static double ked(double freq, double rxh, double dkm)
 		d = (n - 2) * elev[1];	// no of points * delta = km
 
 		//Find dip(s)
-		if (elev[n] < (obh + dipheight)) {
+		if (elev[n] < obh) {
 
 			// Angle from Rx point to obstacle
 			rxobaoi =
@@ -218,9 +219,9 @@ static double ked(double freq, double rxh, double dkm)
 	}
 
 	if (rxobaoi >= 0) {
-		return rxobaoi / (300 / freq);	// Diffraction angle divided by wavelength (m)
+		return (rxobaoi / (300 / freq))+3;	// Diffraction angle divided by wavelength (m)
 	} else {
-		return 0;
+		return 1;
 	}
 }
 
@@ -530,6 +531,12 @@ void PlotPropPath(struct site source, struct site destination,
 				// Egli VHF/UHF
 				loss = EgliPathLoss(LR.frq_mhz, source.alt * METERS_PER_FOOT, (path.elevation[y] * METERS_PER_FOOT) + (destination.alt * METERS_PER_FOOT),dkm);
 				break;
+                        case 12:
+                                // Soil
+                                loss = SoilPathLoss(LR.frq_mhz, dkm, LR.eps_dielect);
+                                break;
+
+
 			default:
 				point_to_point_ITM(source.alt * METERS_PER_FOOT,
 						   destination.alt *
@@ -543,8 +550,8 @@ void PlotPropPath(struct site source, struct site destination,
 
 			}
 
-		
+
-			if (knifeedge == 1) {
+			if (knifeedge == 1 && propmodel > 1) {
 				diffloss =
 				    ked(LR.frq_mhz,
 					destination.alt * METERS_PER_FOOT, dkm);

models/soil.cc

@@ -0,0 +1,33 @@
+/*****************************************************************************
+*  Soil Path Loss model for Signal Server by Alex Farrant                    *
+*  21 February 2018                                                          *
+*                                                                            *
+*  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.                                                         *
+*
+* Frequency: Any MHz
+* Distance: Any Km
+* Terrain permittivity: 1 - 15 (Bad to Good)
+*/
+
+#include <math.h>
+
+// use call with log/ln as this may be faster
+// use constant of value 20.0/log(10.0)
+static __inline float _20log10f(float x)
+{
+  return(8.685889f*logf(x));
+}
+
+double SoilPathLoss(float f, float d, float terdic)
+{
+  float soil = (120/terdic);
+  return(6.4 + _20log10f(d) + _20log10f(f)+(8.69*soil));
+}

models/soil.hh

@@ -0,0 +1,6 @@
+#ifndef _SOIL_HH_
+#define _SOIL_HH_
+
+double SoilPathLoss(float f, float d, float t);
+
+#endif /* _SOIL_HH_ */

models/sui.cc

@@ -2,44 +2,55 @@
 #include <stdlib.h>
 #include <math.h>
 
+
+// use call with log/ln as this may be faster
+// use constant of value 20.0/log(10.0)
+static __inline float _20log10f(float x)
+{
+  return(8.685889f*logf(x));
+}
+
+
 double SUIpathLoss(double f, double TxH, double RxH, double d, int mode)
 {
-	/*
+        /*
-	   f = Frequency (MHz) 1900 to 11000
+           f = Frequency (MHz) 1900 to 11000
-	   TxH =  Transmitter height (m)
+           TxH =  Transmitter height (m)
-	   RxH = Receiver height (m)
+           RxH = Receiver height (m)
-	   d = distance (km)
+           d = distance (km)
-	   mode A1 = URBAN / OBSTRUCTED
+           mode A1 = URBAN / OBSTRUCTED
-	   mode B2 = SUBURBAN / PARTIALLY OBSTRUCTED
+           mode B2 = SUBURBAN / PARTIALLY OBSTRUCTED
-	   mode C3 = RURAL / OPEN
+           mode C3 = RURAL / OPEN
-	   http://www.cl.cam.ac.uk/research/dtg/lce-pub/public/vsa23/VTC05_Empirical.pdf
+           http://www.cl.cam.ac.uk/research/dtg/lce-pub/public/vsa23/VTC05_Empirical.pdf
-	 */
+         */
-	d = d * 1000.0;		// km to m
+        d *= 1e3;               // km to m
 
-	// Urban (A1) is default
+        // Urban (A1) is default
-	double a = 4.6;
+        float a = 4.6;
-	double b = 0.0075;
+        float b = 0.0075;
-	double c = 12.6;
+        float c = 12.6;
-	double s = 10.6; // Optional fading value. Max 10.6dB
+        float s = 10.6; // Optional fading value. Max 10.6dB
-	double XhCF = -10.8;
+        float XhCF = -10.8;
 
-	if (mode == 2) { // Suburban
+        if (mode == 2) { // Suburban
-		a = 4.0;
+                a = 4.0;
-		b = 0.0065;
+                b = 0.0065;
-		c = 17.1;
+                c = 17.1;
-	}
+		XhCF = -10.8;
-	if (mode == 3) { // Rural
+        }
-		a = 3.6;
+        if (mode == 3) { // Rural
-		b = 0.005;
+                a = 3.6;
-		c = 20;
+                b = 0.005;
-		XhCF = -20;
+                c = 20;
-	}
+                XhCF = -20;
-	double d0 = 100;
+        }
-	double A = 20 * log10((4 * M_PI * d0) / (300.0 / f));
+        float d0 = 100.0;
-	double y = a - (b * TxH) + (c / TxH);
+        float A = _20log10f((4 * M_PI * d0) / (300.0 / f));
-	//Correction factors
+        float y = a - (b * TxH) + (c / TxH);
-	double Xf = 6.0 * log10(f / 2000);
+        //Correction factors
-	double Xh = XhCF * log10(RxH / 2000);
+        float Xf = 6.0 * log10(f / 2000.0);
-
+        float Xh = XhCF * log10(RxH / 2000.0);
-	return A + (10 * y * log10(d / d0)) + Xf + Xh + s;
+        Xh *=-1;
+        //fprintf(stdout,"A %lf y %lf Xf %lf Xh %lf\n",A,y,Xf,Xh);
+        return A + (10 * y * log10(d / d0)) + Xf + Xh + s;
 }

suitest.cc

@@ -0,0 +1,64 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+
+// use call with log/ln as this may be faster
+// use constant of value 20.0/log(10.0)
+static __inline float _20log10f(float x)
+{
+  return(8.685889f*logf(x));
+}
+
+double SUIpathLoss(double f, double TxH, double RxH, double d, int mode)
+{
+        /*
+           f = Frequency (MHz) 1900 to 11000
+           TxH =  Transmitter height (m)
+           RxH = Receiver height (m)
+           d = distance (km)
+           mode A1 = URBAN / OBSTRUCTED
+           mode B2 = SUBURBAN / PARTIALLY OBSTRUCTED
+           mode C3 = RURAL / OPEN
+           http://www.cl.cam.ac.uk/research/dtg/lce-pub/public/vsa23/VTC05_Empirical.pdf
+         */
+	//f /=1e3;
+        d *= 1e3;               // km to m
+
+        // Urban (A1) is default
+        float a = 4.6;
+        float b = 0.0075;
+        float c = 12.6;
+        float s = 10.6; // Optional fading value. Max 10.6dB
+        float XhCF = -10.8;
+
+        if (mode == 2) { // Suburban
+                a = 4.0;
+                b = 0.0065;
+                c = 17.1;
+        }
+        if (mode == 3) { // Rural
+                a = 3.6;
+                b = 0.005;
+                c = 20;
+                XhCF = -20;
+        }
+        float d0 = 100;
+        float A = _20log10f((4 * M_PI * d0) / (300.0 / f));
+        float y = a - (b * TxH) + (c / TxH);
+        //Correction factors
+        float Xf = 6.0 * log10(f / 2000.0);
+        float Xh = XhCF * log10(RxH / 2000.0);
+	Xh *=-1;
+	fprintf(stdout,"A %lf y %lf Xf %lf Xh %lf\n",A,y,Xf,Xh);
+        return A + (10 * y * log10(d / d0)) + Xf + Xh + s;
+
+
+}
+
+int main(){
+
+	fprintf(stdout,"%lf\n",SUIpathLoss(3500.0,15.0,10.0,1,1));
+	fprintf(stdout,"%lf\n",SUIpathLoss(3500.0,15.0,10.0,1,2));
+	fprintf(stdout,"%lf\n",SUIpathLoss(3500.0,15.0,10.0,1,3));
+}

test.sh

@@ -1,16 +1,20 @@
 #!/bin/bash
 mkdir tests
 RAD=5
-MAXRAD=20
+MAXRAD=50
 FRQ=446
 ERP=25
 
 echo "Running 50cm LIDAR test..."
-./signalserverLIDAR -lid data/sk3587_50cm.asc -lat 53.383 -lon -1.468 -txh 8 -f 446 -erp 1 -rxh 2 -m -dbm -rt -95 -o tests/lidar_50cm -R 0.5 -t
+./signalserverLIDAR -lid data/sk3587_50cm.asc -lat 53.383 -lon -1.468 -txh 8 -f $FRQ -erp 1 -rxh 2 -m -dbm -rt -90 -o tests/1 -R 0.5 -t
 echo "Converting to PNG..."
-convert tests/lidar_50cm.ppm -transparent white -channel Alpha PNG32:tests/lidar_50cm.png
+convert tests/1.ppm -transparent white -channel Alpha PNG32:tests/1.png
-rm tests/lidar_50cm.ppm
+
-rm tests/lidar_50cm.*cf
+echo "Running 30m Meridian test..."
+./signalserverLIDAR -lid data/N051E000_AVE_DSM.tif.asc,data/N051W001_AVE_DSM.tif.asc -lat 51.472 -lon 0.0096 -txh 12 -f $FRQ -erp 5 -rxh 2 -m -dbm -rt -100 -o tests/2 -R 10 -t
+echo "Converting to PNG..."
+convert tests/2.ppm -transparent white -channel Alpha PNG32:tests/2.png
+
 
 echo "Running soak test out to $MAXRAD"