diff --git a/CHANGELOG b/CHANGELOG index 8e730db..338c08f 100644 --- a/CHANGELOG +++ b/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 diff --git a/Makefile b/Makefile index 828162f..10b5f16 100644 --- a/Makefile +++ b/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 diff --git a/data/terraindata.tgz b/data/terraindata.tgz new file mode 100644 index 0000000..afbda62 Binary files /dev/null and b/data/terraindata.tgz differ diff --git a/inputs.cc b/inputs.cc index cdf042b..0e76647 100644 --- a/inputs.cc +++ b/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 < min_west) - min_west = tiles[indx].min_west; + if (tiles[indx].min_west > 359.9) { + 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) - min_west = tiles[indx].min_west; + if (fabs(tiles[indx].min_west - min_west) < 180.0) { + 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) - min_west=0; - + if(eastF && fc>1 && min_west<=1.0001){ + 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; diff --git a/main.cc b/main.cc index ac2ed75..cfb5b15 100644 --- a/main.cc +++ b/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)); diff --git a/models/los.cc b/models/los.cc index 87b48f4..b6720a7 100644 --- a/models/los.cc +++ b/models/los.cc @@ -11,6 +11,7 @@ #include "sui.hh" #include "pel.hh" #include "egli.hh" +#include "soil.hh" #include #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); diff --git a/models/soil.cc b/models/soil.cc new file mode 100644 index 0000000..d824a2e --- /dev/null +++ b/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 + +// 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)); +} diff --git a/models/soil.hh b/models/soil.hh new file mode 100644 index 0000000..f1278f2 --- /dev/null +++ b/models/soil.hh @@ -0,0 +1,6 @@ +#ifndef _SOIL_HH_ +#define _SOIL_HH_ + +double SoilPathLoss(float f, float d, float t); + +#endif /* _SOIL_HH_ */ diff --git a/models/sui.cc b/models/sui.cc index a1b969b..c82e90d 100644 --- a/models/sui.cc +++ b/models/sui.cc @@ -2,44 +2,55 @@ #include #include + +// 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 - */ - d = d * 1000.0; // km to m + /* + 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 + */ + d *= 1e3; // km to m - // Urban (A1) is default - double a = 4.6; - double b = 0.0075; - double c = 12.6; - double s = 10.6; // Optional fading value. Max 10.6dB - double XhCF = -10.8; + // 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; - } - double d0 = 100; - double A = 20 * log10((4 * M_PI * d0) / (300.0 / f)); - double y = a - (b * TxH) + (c / TxH); - //Correction factors - double Xf = 6.0 * log10(f / 2000); - double Xh = XhCF * log10(RxH / 2000); - - return A + (10 * y * log10(d / d0)) + Xf + Xh + s; + if (mode == 2) { // Suburban + a = 4.0; + b = 0.0065; + c = 17.1; + XhCF = -10.8; + } + if (mode == 3) { // Rural + a = 3.6; + b = 0.005; + c = 20; + XhCF = -20; + } + float d0 = 100.0; + 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; } diff --git a/suitest.cc b/suitest.cc new file mode 100644 index 0000000..e075195 --- /dev/null +++ b/suitest.cc @@ -0,0 +1,64 @@ +#include +#include +#include + + +// 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)); +} diff --git a/test.sh b/test.sh index e0d52c7..b245ef5 100755 --- a/test.sh +++ b/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 -rm tests/lidar_50cm.ppm -rm tests/lidar_50cm.*cf +convert tests/1.ppm -transparent white -channel Alpha PNG32:tests/1.png + +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"