diff --git a/CHANGELOG b/CHANGELOG index 1a0aad0..3ffa688 100644 --- a/CHANGELOG +++ b/CHANGELOG @@ -1,4 +1,9 @@ SIGNAL SERVER CHANGELOG +2.7 - 03 Jan 2016 +Added support for LIDAR data in ASCII grid format with WGS84 bounds +Increased resolution to +54000 pixels per degree (2m) +Removed support for Windows SDF filenames without:colons. + 2.63 - 10 Nov 2015 Added sanity check and handicap to ECC33 model when used with low Tx heights in hilly areas = sea of red @@ -11,10 +16,10 @@ Credit to Nils Lofstad for helping nail this down. FSPL model floor reduced from 150MHz to 20MHz. 2.6 - 9 June 2015 -Multithreading support added by Michael Ramnarine +Multi-threading support added by Michael Ramnarine PlotPropagation() and PlotLOSMap() use four threads by default Feature can be disabled with -nothreads flag -Static and global variables have been made threadsafe +Static and global variables have been made thread safe 2.5 - 27 May 2015 Code refactored by Andrew Clayton / ac000 with header files @@ -109,7 +114,7 @@ New Earth Dielectric range (Permittivity): 80 to 0.1 New Conductivity range (Siemens/m): 0.01 to 0.000001 1.0 - 19 November 2011 -SS released. +Signal Server forked from SPLAT! diff --git a/Makefile b/Makefile index 8ae3401..ea51ead 100644 --- a/Makefile +++ b/Makefile @@ -35,7 +35,9 @@ signalserver: $(objects) @$(CXX) $(objects) -o $@ ${LIBS} @echo -e " SYMLNK\tsignalserverHD -> $@" @ln -sf $@ signalserverHD - + @echo -e " SYMLNK\tsignalserverLIDAR -> $@" + @ln -sf $@ signalserverLIDAR + main.o: main.cc common.h inputs.hh outputs.hh itwom3.0.hh los.hh inputs.o: inputs.cc common.h main.hh @@ -48,4 +50,4 @@ los.o: los.cc common.h main.hh cost.hh ecc33.hh ericsson.hh fspl.hh hata.hh \ .PHONY: clean clean: - rm -f $(objects) signalserver signalserverHD + rm -f $(objects) signalserver signalserverHD signalserverLIDAR diff --git a/README.md b/README.md index c348498..d796d0e 100644 --- a/README.md +++ b/README.md @@ -1,11 +1,11 @@ /****************************************************************************\ -* Signal Server: Server optimised SPLAT! by Alex Farrant * +* Signal Server: Server optimised SPLAT! by Alex Farrant, M6ZUJ * ****************************************************************************** * 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. * +* 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 * @@ -20,7 +20,7 @@ * * \****************************************************************************/ - -- Signal Server 2.6 -- + -- Signal Server -- Compiled for 64 tiles at 1200 pixels/degree -d Directory containing .sdf tiles @@ -53,5 +53,17 @@ -ng Normalise Path Profile graph -haf Halve 1 or 2 (optional) -nothreads Turn off threaded processing (optional) + + Example usage: + + INPUTS: 900MHz tower at 25m AGL with 5W ERP + OUTPUTS: 1200 resolution, 30km radius, -90dBm receiver threshold, Longley Rice model + ./signalserver -d /data/SRTM3 -lat 51.849 -lon -2.2299 -txh 25 -f 900 -erp 5 -rxh 2 -rt -90 -dbm -m -o test1 -R 30 -res 1200 -pm 1 + INPUTS: 450MHz tower at 25f AGL with 20W ERP + OUTPUTS: 3600 resolution, 30km radius, 10dBuV receiver threshold, Hata model + ./signalserverHD -d /data/SRTM1 -lat 51.849 -lon -2.2299 -txh 25 -f 450 -erp 20 -rxh 2 -rt 10 -o test2 -R 30 -res 3600 -pm 3 + INPUTS: 1800MHz tower at 15m AGL with 1W ERP + OUTPUTS: 2m LIDAR resolution, 5km radius, -90dBm receiver threshold, Longley Rice model + ./signalserverLIDAR -lid /data/LIDAR2m/Gloucester_2m.asc -lat 51.849 -lon -2.2299 -txh 15 -f 1800 -erp 1 -rxh 2 -rt -90 -dbm -m -o test3 -R 30 -res 5000 -pm 1 \ No newline at end of file diff --git a/common.h b/common.h index e98e500..61feb60 100644 --- a/common.h +++ b/common.h @@ -23,10 +23,10 @@ #define FOUR_THIRDS 1.3333333333333 struct dem { - int min_north; - int max_north; - int min_west; - int max_west; + float min_north; + float max_north; + float min_west; + float max_west; int max_el; int min_el; short **data; @@ -73,10 +73,10 @@ extern int MAXPAGES; extern int ARRAYSIZE; extern int IPPD; -extern int min_north; -extern int max_north; -extern int min_west; -extern int max_west; +extern double min_north; +extern double max_north; +extern double min_west; +extern double max_west; extern int ippd; extern int mpi; extern int max_elevation; @@ -84,6 +84,8 @@ extern int min_elevation; extern int contour_threshold; extern int loops; extern int jgets; +extern int width; +extern int height; extern double earthradius; extern double north; @@ -93,6 +95,7 @@ extern double west; extern double max_range; extern double dpp; extern double ppd; +extern double yppd; extern double fzone_clearance; extern double clutter; extern double dBm; diff --git a/inputs.cc b/inputs.cc index f9cb009..7149f1c 100644 --- a/inputs.cc +++ b/inputs.cc @@ -3,11 +3,144 @@ #include #include #include - #include "common.h" #include "main.hh" -int LoadSDF_SDF(char *name, int winfiles) +int loadLIDAR(char *filename) +{ + /* This function reads a single LIDAR tile of n rows and n columns in ASCII grid format. + The tile must have WGS84 bounds in the header in the order: WEST,SOUTH,EAST,NORTH + ncols 5000 + nrows 5000 + xllcorner -2.291359 + yllcorner 51.788295 + xurcorner -2.146674 + yurcorner 51.878474 + cellsize 2 + NODATA_value -9999 + */ + + int x, y, width, height, cellsize; + double xll, yll, xur, yur; + char found, free_page = 0, line[50000], jline[20], lid_file[255], + path_plus_name[255], *junk = NULL; + FILE *fd; + char * pch; + fd = fopen(filename, "rb"); + if (fd != NULL) { + if (debug) { + fprintf(stdout,"Loading \"%s\"...\n", filename); + fflush(stdout); + } + if (fgets(line, 20, fd) != NULL) { + pch = strtok (line," "); + pch = strtok (NULL, " "); + width=atoi(pch); + } + if (fgets(line, 20, fd) != NULL) { + height=atoi(pch); + } + fgets(line, 21, fd); // + + if (fgets(line, 22, fd) != NULL) { + xll=atof(pch); + } + + fgets(line, 21, fd); // + + if (fgets(line, 22, fd) != NULL) { + yll=atof(pch); + } + + fgets(line, 21, fd); // + + if (fgets(line, 22, fd) != NULL) { + xur=atof(pch); + } + + fgets(line, 21, fd); // + + if (fgets(line, 22, fd) != NULL) { + yur=atof(pch); + } + + fgets(line, 21, fd); // + + if (fgets(line, 21, fd) != NULL) { + cellsize=atoi(pch); + } + + // Transform WGS84 longitudes into 'west' values as society finishes east of Greenwich ;) + if(xll > 0){ + xll=360-xll; + } + if(xur > 0){ + xur=360-xur; + } + if(xll < 0){ + xll=xll*-1; + } + if(xur < 0){ + xur=xur*-1; + } + dem[0].min_north=yll; + min_north=yll; + dem[0].max_north=yur; + max_north=yur; + dem[0].min_west=xur; + min_west=xur; + dem[0].max_west=xll; + max_west=xll; + if(debug){ + fprintf(stdout,"yll %.2f yur %.2f xur %.2f xll %.2f\n",yll,yur,xur,xll); + } + if(width!=height){ + fprintf(stdout,"LIDAR tile is not a square. Rows != Columns\n"); + return 0; + } + + fgets(line, 30, fd); // NODATA + for (y = height-1; y > -1; y--) { + x=width-1; + if (fgets(line, 50000,fd) != NULL) { + pch = strtok (line, " "); // 500 + while(pch != NULL){ + if(atoi(pch)<-999){ + pch="0"; + } + + dem[0].data[y][x]=atoi(pch); + dem[0].signal[x][y] = 0; + dem[0].mask[x][y] = 0; + + if (atoi(pch) > dem[0].max_el){ + dem[0].max_el = atoi(pch); + max_elevation = atoi(pch); + } + + if (atoi(pch) < dem[0].min_el){ + dem[0].min_el = atoi(pch); + min_elevation = dem[0].min_el; + } + + + x--; + pch = strtok (NULL, " "); // 500 + } + }else{ + fprintf(stdout,"LIDAR error @ line %d\n",x); + return 0; + } + } + + fclose(fd); + return 0; + } + else + return -1; +} + +int LoadSDF_SDF(char *name) { /* This function reads uncompressed ss Data Files (.sdf) containing digital elevation model data into memory. @@ -27,13 +160,10 @@ int LoadSDF_SDF(char *name, int winfiles) sdf_file[x] = 0; /* Parse filename for minimum latitude and longitude values */ - if (winfiles == 1) { - sscanf(sdf_file, "%d=%d=%d=%d", &minlat, &maxlat, &minlon, - &maxlon); - } else { + sscanf(sdf_file, "%d:%d:%d:%d", &minlat, &maxlat, &minlon, &maxlon); - } + sdf_file[x] = '.'; sdf_file[x + 1] = 's'; @@ -87,19 +217,19 @@ int LoadSDF_SDF(char *name, int winfiles) } if (fgets(line, 19, fd) != NULL) { - sscanf(line, "%d", &dem[indx].max_west); + sscanf(line, "%f", &dem[indx].max_west); } if (fgets(line, 19, fd) != NULL) { - sscanf(line, "%d", &dem[indx].min_north); + sscanf(line, "%f", &dem[indx].min_north); } if (fgets(line, 19, fd) != NULL) { - sscanf(line, "%d", &dem[indx].min_west); + sscanf(line, "%f", &dem[indx].min_west); } if (fgets(line, 19, fd) != NULL) { - sscanf(line, "%d", &dem[indx].max_north); + sscanf(line, "%f", &dem[indx].max_north); } /* Here X lines of DEM will be read until IPPD is reached. @@ -205,7 +335,7 @@ int LoadSDF_SDF(char *name, int winfiles) return 0; } -char LoadSDF(char *name, int winfiles) +char LoadSDF(char *name) { /* This function loads the requested SDF file from the filesystem. It first tries to invoke the LoadSDF_SDF() function to load an @@ -220,19 +350,16 @@ char LoadSDF(char *name, int winfiles) char found, free_page = 0; int return_value = -1; - return_value = LoadSDF_SDF(name, winfiles); + return_value = LoadSDF_SDF(name); /* If neither format can be found, then assume the area is water. */ if (return_value == 0 || return_value == -1) { - if (winfiles == 1) { - sscanf(name, "%d=%d=%d=%d", &minlat, &maxlat, &minlon, - &maxlon); - } else { + sscanf(name, "%d:%d:%d:%d", &minlat, &maxlat, &minlon, &maxlon); - } + /* Is it already in memory? */ for (indx = 0, found = 0; indx < MAXPAGES && found == 0; indx++) { @@ -1138,8 +1265,7 @@ void LoadDBMColors(struct site xmtr) } } -void LoadTopoData(int max_lon, int min_lon, int max_lat, int min_lat, - int winfiles) +void LoadTopoData(int max_lon, int min_lon, int max_lat, int min_lat) { /* This function loads the SDF files required to cover the limits of the region specified. */ @@ -1167,17 +1293,7 @@ void LoadTopoData(int max_lon, int min_lon, int max_lat, int min_lat, while (ymax >= 360) ymax -= 360; - if (winfiles == 1) { - if (ippd == 3600) - snprintf(string, 19, - "%d=%d=%d=%d=hd", x, - x + 1, ymin, ymax); - else - snprintf(string, 16, - "%d=%d=%d=%d", x, - x + 1, ymin, ymax); - } else { if (ippd == 3600) snprintf(string, 19, "%d:%d:%d:%d-hd", x, @@ -1186,9 +1302,7 @@ void LoadTopoData(int max_lon, int min_lon, int max_lat, int min_lat, snprintf(string, 16, "%d:%d:%d:%d", x, x + 1, ymin, ymax); - } - - LoadSDF(string, winfiles); + LoadSDF(string); } } @@ -1211,17 +1325,6 @@ void LoadTopoData(int max_lon, int min_lon, int max_lat, int min_lat, while (ymax >= 360) ymax -= 360; - if (winfiles == 1) { - if (ippd == 3600) - snprintf(string, 19, - "%d=%d=%d=%d=hd", x, - x + 1, ymin, ymax); - else - snprintf(string, 16, - "%d=%d=%d=%d", x, - x + 1, ymin, ymax); - - } else { if (ippd == 3600) snprintf(string, 19, "%d:%d:%d:%d-hd", x, @@ -1230,9 +1333,7 @@ void LoadTopoData(int max_lon, int min_lon, int max_lat, int min_lat, snprintf(string, 16, "%d:%d:%d:%d", x, x + 1, ymin, ymax); - } - - LoadSDF(string, winfiles); + LoadSDF(string); } } } diff --git a/inputs.hh b/inputs.hh index 64ddb96..8a19945 100644 --- a/inputs.hh +++ b/inputs.hh @@ -9,8 +9,8 @@ void LoadPAT(char *filename); void LoadSignalColors(struct site xmtr); void LoadLossColors(struct site xmtr); void LoadDBMColors(struct site xmtr); -void LoadTopoData(int max_lon, int min_lon, int max_lat, int min_lat, - int winfiles); +void LoadTopoData(int max_lon, int min_lon, int max_lat, int min_lat); void LoadUDT(char *filename); +int loadLIDAR(char *filename); #endif /* _INPUTS_HH_ */ diff --git a/main.cc b/main.cc index 26069e5..c074d87 100644 --- a/main.cc +++ b/main.cc @@ -1,6 +1,6 @@ -double version = 2.63; +double version = 2.7; /****************************************************************************\ -* Signal Server: Server optimised SPLAT! by Alex Farrant * +* Signal Server: Server optimised SPLAT! by Alex Farrant, M6ZUJ * ****************************************************************************** * SPLAT! Project started in 1997 by John A. Magliacane, KD2BD * * * @@ -34,7 +34,6 @@ double version = 2.63; #include "models/itwom3.0.hh" #include "models/los.hh" -// 90m mode (default) int MAXPAGES = 64; int ARRAYSIZE = 76810; int IPPD = 1200; @@ -42,14 +41,15 @@ int IPPD = 1200; 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, +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; + north, east, south, west, dBm, loss, field_strength, + min_north = 90, max_north = -90, min_west = 360, max_west = -1; -int min_north = 90, max_north = -90, min_west = 360, max_west = -1, ippd, mpi, +int ippd, mpi, max_elevation = -32768, min_elevation = 32768, bzerror, contour_threshold, pred, pblue, pgreen, ter, multiplier = 256, debug = 0, loops = 64, jgets = - 0, MAXRAD, hottest = 10; + 0, MAXRAD, hottest = 10, height, width; unsigned char got_elevation_pattern, got_azimuth_pattern, metric = 0, dbm = 0; @@ -161,7 +161,7 @@ int PutMask(double lat, double lon, int value) for (indx = 0, found = 0; indx < MAXPAGES && found == 0;) { x = (int)rint(ppd * (lat - dem[indx].min_north)); - y = mpi - (int)rint(ppd * (LonDiff(dem[indx].max_west, lon))); + y = mpi - (int)rint(yppd * (LonDiff(dem[indx].max_west, lon))); if (x >= 0 && x <= mpi && y >= 0 && y <= mpi) found = 1; @@ -191,8 +191,9 @@ int OrMask(double lat, double lon, int value) for (indx = 0, found = 0; indx < MAXPAGES && found == 0;) { x = (int)rint(ppd * (lat - dem[indx].min_north)); - y = mpi - (int)rint(ppd * (LonDiff(dem[indx].max_west, lon))); + y = mpi - (int)rint(yppd * (LonDiff(dem[indx].max_west, lon))); + if (x >= 0 && x <= mpi && y >= 0 && y <= mpi) found = 1; else @@ -226,14 +227,14 @@ int PutSignal(double lat, double lon, unsigned char signal) 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(ppd * (LonDiff(dem[indx].max_west, lon))); + y = mpi - (int)rint(yppd * (LonDiff(dem[indx].max_west, lon))); if (x >= 0 && x <= mpi && y >= 0 && y <= mpi) found = 1; @@ -263,7 +264,7 @@ unsigned char GetSignal(double lat, double lon) for (indx = 0, found = 0; indx < MAXPAGES && found == 0;) { x = (int)rint(ppd * (lat - dem[indx].min_north)); - y = mpi - (int)rint(ppd * (LonDiff(dem[indx].max_west, lon))); + y = mpi - (int)rint(yppd * (LonDiff(dem[indx].max_west, lon))); if (x >= 0 && x <= mpi && y >= 0 && y <= mpi) found = 1; @@ -290,7 +291,7 @@ double GetElevation(struct site location) for (indx = 0, found = 0; indx < MAXPAGES && found == 0;) { x = (int)rint(ppd * (location.lat - dem[indx].min_north)); y = mpi - - (int)rint(ppd * + (int)rint(yppd * (LonDiff(dem[indx].max_west, location.lon))); if (x >= 0 && x <= mpi && y >= 0 && y <= mpi) @@ -319,7 +320,7 @@ int AddElevation(double lat, double lon, double height) for (indx = 0, found = 0; indx < MAXPAGES && found == 0;) { x = (int)rint(ppd * (lat - dem[indx].min_north)); - y = mpi - (int)rint(ppd * (LonDiff(dem[indx].max_west, lon))); + y = mpi - (int)rint(yppd * (LonDiff(dem[indx].max_west, lon))); if (x >= 0 && x <= mpi && y >= 0 && y <= mpi) found = 1; @@ -455,9 +456,7 @@ void ReadPath(struct site source, struct site 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; } @@ -986,25 +985,32 @@ 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, winfiles, knifeedge = 0, ppa = - 0, normalise = 0, haf = 0, pmenv = 1; + nortRxHin, nortRxHax, propmodel, knifeedge = 0, ppa = + 0, normalise = 0, haf = 0, pmenv = 1, lidar=0; 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], udt_file[255], ano_filename[255]; + char mapfile[255], udt_file[255], ano_filename[255], lidar_tile[255]; 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 = 14844; + ARRAYSIZE = 32410; IPPD = 3600; } + if (strstr(argv[0], "signalserverLIDAR")) { + MAXPAGES = 1; // 10km2 for now :) + IPPD = 5000; // // 2m resolution + ARRAYSIZE = 5010; + ppd=IPPD; + } + strncpy(ss_name, "Signal Server\0", 14); if (argc == 1) { @@ -1013,7 +1019,8 @@ int main(int argc, char *argv[]) fprintf(stdout, "\tSet for %d tiles at %d pixels/degree\n\n", MAXPAGES, IPPD); - fprintf(stdout, " -d Directory containing .sdf tiles\n"); + fprintf(stdout, " -sdf Directory containing .sdf tiles\n"); + fprintf(stdout, " -lid LIDAR ASCII tile with WGS84 bounds\n"); fprintf(stdout, " -lat Tx Latitude (decimal degrees) -70/+70\n"); fprintf(stdout, @@ -1046,7 +1053,7 @@ int main(int argc, char *argv[]) fprintf(stdout, " -o Filename. Required. \n"); fprintf(stdout, " -R Radius (miles/kilometers)\n"); fprintf(stdout, - " -res Pixels per degree. 300/600/1200(default)/3600 (optional)\n"); + " -res Pixels per tile. 300/600/1200/3600/5000/10000 (optional)\n"); fprintf(stdout, " -t Terrain background\n"); fprintf(stdout, " -pm Prop model. 1: ITM, 2: LOS, 3: Hata, 4: ECC33,\n"); @@ -1074,7 +1081,6 @@ int main(int argc, char *argv[]) alloc_path(); y = argc - 1; - kml = 0; geo = 0; dbm = 0; @@ -1096,8 +1102,6 @@ int main(int argc, char *argv[]) earthradius = EARTHRADIUS; max_range = 1.0; propmodel = 1; //ITM - winfiles = 0; - lat = 0; lon = 0; txh = 0; @@ -1146,25 +1150,32 @@ int main(int argc, char *argv[]) switch (ippd) { case 300: - MAXRAD = 300; - jgets = 3; + MAXRAD = 500; + jgets = 3; // 3 dummy reads break; - case 600: - MAXRAD = 150; + MAXRAD = 500; jgets = 1; break; - + case 1200: + MAXRAD = 200; + ippd = 1200; + break; case 3600: MAXRAD = 100; ippd = 3600; - jgets = 0; break; - + case 5000: // LIDAR 2m + MAXRAD = 10; + ippd = 5000; + break; + case 10000: // LIDAR 1m + MAXRAD = 5; + ippd = 10000; + break; default: - MAXRAD = 100; + MAXRAD = 200; ippd = 1200; - jgets = 0; break; } } @@ -1220,13 +1231,20 @@ int main(int argc, char *argv[]) if (strcmp(argv[x], "-dbm") == 0) dbm = 1; - if (strcmp(argv[x], "-d") == 0) { + 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_tile, argv[z], 253); + } + if (strcmp(argv[x], "-lat") == 0) { z = x + 1; @@ -1382,10 +1400,7 @@ int main(int argc, char *argv[]) debug = 1; } - ppd = (double)ippd; /* pixels per degree (double) */ - dpp = 1.0 / ppd; /* degrees per pixel */ - mpi = ippd - 1; /* maximum pixel index per degree */ - + /*UDT*/ if (strcmp(argv[x], "-udt") == 0) { z = x + 1; @@ -1416,11 +1431,7 @@ int main(int argc, char *argv[]) z = x + 1; knifeedge = 1; } - //Windows friendly SDF filenames - if (strcmp(argv[x], "-wf") == 0) { - z = x + 1; - winfiles = 1; - } + //Normalise Path Profile chart if (strcmp(argv[x], "-ng") == 0) { z = x + 1; @@ -1486,7 +1497,7 @@ int main(int argc, char *argv[]) "ERROR: Rx altitude above ground was too high!"); exit(0); } - if (ippd < 300 || ippd > 3600) { + if (ippd < 300 || ippd > 10000) { fprintf(stdout, "ERROR: resolution out of range!"); exit(0); } @@ -1562,96 +1573,110 @@ int main(int argc, char *argv[]) max_lon = rxlon; } - /* Load the required SDF files */ + /* Load the required tiles */ + if(lidar){ + loadLIDAR(lidar_tile); + ppd=rint(ippd / (max_north-min_north)); + yppd=rint(ippd / (max_west-min_west)); + height=ippd; + width=ippd; + }else{ + LoadTopoData(max_lon, min_lon, max_lat, min_lat); + 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. */ - LoadTopoData(max_lon, min_lon, max_lat, min_lat, winfiles); + tx_range = + sqrt(1.5 * + (tx_site[z].alt + GetElevation(tx_site[z]))); - 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. */ + if (LRmap) + rx_range = sqrt(1.5 * altitudeLR); + else + rx_range = sqrt(1.5 * altitude); - tx_range = - sqrt(1.5 * - (tx_site[z].alt + GetElevation(tx_site[z]))); + /* deg_range determines the maximum + amount of topo data we read */ - if (LRmap) - rx_range = sqrt(1.5 * altitudeLR); - else - rx_range = sqrt(1.5 * altitude); + deg_range = (tx_range + rx_range) / 57.0; - /* deg_range determines the maximum - amount of topo data we read */ + /* 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. */ - deg_range = (tx_range + rx_range) / 57.0; + if (max_range == 0.0) + max_range = tx_range + rx_range; - /* 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. */ + deg_range = max_range / 57.0; - if (max_range == 0.0) - max_range = tx_range + rx_range; + // No more than 8 degs + deg_limit = 3.5; - deg_range = max_range / 57.0; + 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); - // No more than 8 degs - deg_limit = 3.5; + /* Correct for squares in degrees not being square in miles */ - 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); + if (deg_range > deg_limit) + deg_range = deg_limit; - /* Correct for squares in degrees not being square in miles */ + if (deg_range_lon > deg_limit) + deg_range_lon = deg_limit; - if (deg_range > deg_limit) - deg_range = deg_limit; + nortRxHin = (int)floor(tx_site[z].lat - deg_range); + nortRxHax = (int)floor(tx_site[z].lat + deg_range); - if (deg_range_lon > deg_limit) - deg_range_lon = deg_limit; + west_min = (int)floor(tx_site[z].lon - deg_range_lon); - nortRxHin = (int)floor(tx_site[z].lat - deg_range); - nortRxHax = (int)floor(tx_site[z].lat + deg_range); + while (west_min < 0) + west_min += 360; - west_min = (int)floor(tx_site[z].lon - deg_range_lon); + while (west_min >= 360) + west_min -= 360; - while (west_min < 0) - west_min += 360; + west_max = (int)floor(tx_site[z].lon + deg_range_lon); - while (west_min >= 360) - west_min -= 360; + while (west_max < 0) + west_max += 360; - west_max = (int)floor(tx_site[z].lon + deg_range_lon); + while (west_max >= 360) + west_max -= 360; - while (west_max < 0) - west_max += 360; + if (nortRxHin < min_lat) + min_lat = nortRxHin; - while (west_max >= 360) - west_max -= 360; + if (nortRxHax > max_lat) + max_lat = nortRxHax; - if (nortRxHin < min_lat) - min_lat = nortRxHin; + if (LonDiff(west_min, min_lon) < 0.0) + min_lon = west_min; - if (nortRxHax > max_lat) - max_lat = nortRxHax; + if (LonDiff(west_max, max_lon) >= 0.0) + max_lon = west_max; + } - if (LonDiff(west_min, min_lon) < 0.0) - min_lon = west_min; + /* Load any additional SDF files, if required */ - if (LonDiff(west_max, max_lon) >= 0.0) - max_lon = west_max; + LoadTopoData(max_lon, min_lon, max_lat, min_lat); } - - /* Load any additional SDF files, if required */ - - LoadTopoData(max_lon, min_lon, max_lat, min_lat, winfiles); + ppd=(double)ippd; + yppd=ppd; + width = (unsigned)(ippd * ReduceAngle(max_west - min_west)); + height = (unsigned)(ippd * ReduceAngle(max_north - min_north)); } - // UDT clutter + + dpp = 1 / ppd; + mpi = ippd-1; + + // User defined clutter file LoadUDT(udt_file); if (ppa == 0) { @@ -1663,21 +1688,21 @@ int main(int argc, char *argv[]) PlotPropagation(tx_site[0], altitudeLR, ano_filename, propmodel, knifeedge, haf, pmenv, use_threads); - // Near field bugfix - //PutSignal(tx_site[0].lat, tx_site[0].lon, hottest); - - - - for (lat = tx_site[0].lat - 0.001; - lat <= tx_site[0].lat + 0.001; - lat = lat + 0.0001) { - for (lon = tx_site[0].lon - 0.001; - lon <= tx_site[0].lon + 0.001; - lon = lon + 0.0001) { - PutSignal(lat, lon, hottest); + if(!lidar){ + // nearfield bugfix + for (lat = tx_site[0].lat - 0.001; + lat <= tx_site[0].lat + 0.001; + lat = lat + 0.0001) { + for (lon = tx_site[0].lon - 0.001; + lon <= tx_site[0].lon + 0.001; + lon = lon + 0.0001) { + PutSignal(lat, lon, hottest); + } } - } + } + + // Write bitmap if (LR.erp == 0.0) DoPathLoss(mapfile, geo, kml, ngs, tx_site, txsites); @@ -1687,11 +1712,14 @@ int main(int argc, char *argv[]) else DoSigStr(mapfile, geo, kml, ngs, tx_site, txsites); + } - fprintf(stdout, "|%.5f", north); - fprintf(stdout, "|%.5f", east); - fprintf(stdout, "|%.5f", south); - fprintf(stdout, "|%.5f|", west); + + // Print WGS84 bounds + fprintf(stdout, "|%.6f", north); + fprintf(stdout, "|%.6f", east); + fprintf(stdout, "|%.6f", south); + fprintf(stdout, "|%.6f|", west); } else { strncpy(tx_site[0].name, "Tx", 3); @@ -1704,9 +1732,10 @@ int main(int argc, char *argv[]) } fflush(stdout); - free_elev(); + // deprecated garbage collection + /*free_elev(); free_dem(); free_path(); - + */ return 0; } diff --git a/models/los.cc b/models/los.cc index 0264654..870f202 100644 --- a/models/los.cc +++ b/models/los.cc @@ -1,6 +1,5 @@ #include #include - #include "../main.hh" #include "los.hh" #include "cost.hh" @@ -13,7 +12,7 @@ #include #define MAXPAGES 64 -#define IPPD 3600 +#define IPPD 5000 #define NUM_SECTIONS 4 namespace { @@ -25,7 +24,7 @@ namespace { bool has_init_processed = false; struct propagationRange { - int min_west, max_west, min_north, max_north; + double min_west, max_west, min_north, max_north; double altitude; bool eastwest, los, use_threads; site source; @@ -105,7 +104,7 @@ namespace { for (indx = 0, found = 0; indx < MAXPAGES && found == 0;) { x = (int)rint(ppd * (lat - dem[indx].min_north)); - y = mpi - (int)rint(ppd * (LonDiff(dem[indx].max_west, lon))); + y = mpi - (int)rint(yppd * (LonDiff(dem[indx].max_west, lon))); if (x >= 0 && x <= mpi && y >= 0 && y <= mpi) found = 1; @@ -310,6 +309,7 @@ void PlotPropPath(struct site source, struct site destination, elevation[x]) * METERS_PER_FOOT); + /* Copy ending points without clutter */ elev[2] = path.elevation[0] * METERS_PER_FOOT; @@ -332,7 +332,7 @@ void PlotPropPath(struct site source, struct site destination, y++) { /* Process this point only if it has not already been processed. */ - + if ( (GetMask(path.lat[y], path.lon[y]) & 248) != (mask_value << 3) && can_process(path.lat[y], path.lon[y])) { @@ -435,6 +435,7 @@ void PlotPropPath(struct site source, struct site destination, dkm = (elev[1] * elev[0]) / 1000; // km + switch (propmodel) { case 1: // Longley Rice ITM @@ -682,17 +683,17 @@ void PlotLOSMap(struct site source, double altitude, char *plo_filename, if (fd != NULL) { fprintf(fd, - "%d, %d\t; max_west, min_west\n%d, %d\t; max_north, min_north\n", + "%.3f, %.3f\t; max_west, min_west\n%.3f, %.3f\t; max_north, min_north\n", max_west, min_west, max_north, min_north); } // Four sections start here // Process north edge east/west, east edge north/south, // south edge east/west, west edge north/south - int range_min_west[] = {min_west, min_west, min_west, max_west}; - int range_min_north[] = {max_north, min_north, min_north, min_north}; - int range_max_west[] = {max_west, min_west, max_west, max_west}; - int range_max_north[] = {max_north, max_north, min_north, max_north}; + double range_min_west[] = {min_west, min_west, min_west, max_west}; + double range_min_north[] = {max_north, min_north, min_north, min_north}; + double range_max_west[] = {max_west, min_west, max_west, max_west}; + double range_max_north[] = {max_north, max_north, min_north, max_north}; propagationRange* r[NUM_SECTIONS]; for(int i = 0; i < NUM_SECTIONS; ++i) { @@ -773,27 +774,23 @@ void PlotPropagation(struct site source, double altitude, char *plo_filename, metric ? clutter * METERS_PER_FOOT : clutter, metric ? "meters" : "feet"); - if (debug) { - fprintf(stdout, "...\n\n 0%c to 25%c ", 37, 37); - fflush(stdout); - } - if (plo_filename[0] != 0) fd = fopen(plo_filename, "wb"); if (fd != NULL) { fprintf(fd, - "%d, %d\t; max_west, min_west\n%d, %d\t; max_north, min_north\n", + "%.3f, %.3f\t; max_west, min_west\n%.3f, %.3f\t; max_north, min_north\n", max_west, min_west, max_north, min_north); } + // Four sections start here // Process north edge east/west, east edge north/south, // south edge east/west, west edge north/south - int range_min_west[] = {min_west, min_west, min_west, max_west}; - int range_min_north[] = {max_north, min_north, min_north, min_north}; - int range_max_west[] = {max_west, min_west, max_west, max_west}; - int range_max_north[] = {max_north, max_north, min_north, max_north}; + double range_min_west[] = {min_west, min_west, min_west, max_west}; + double range_min_north[] = {max_north, min_north, min_north, min_north}; + double range_max_west[] = {max_west, min_west, max_west, max_west}; + double range_max_north[] = {max_north, max_north, min_north, max_north}; propagationRange* r[NUM_SECTIONS]; for(int i = 0; i < NUM_SECTIONS; ++i) { diff --git a/outputs.cc b/outputs.cc index 44dba1d..0f90687 100644 --- a/outputs.cc +++ b/outputs.cc @@ -25,7 +25,7 @@ void DoPathLoss(char *filename, unsigned char geo, unsigned char kml, points up and east points right in the image generated. */ char mapfile[255]; - unsigned width, height, red, green, blue, terrain = 0; + unsigned red, green, blue, terrain = 0; unsigned char found, mask, cityorcounty; int indx, x, y, z, x0 = 0, y0 = 0, loss, match; double lat, lon, conversion, one_over_gamma, minwest; @@ -36,9 +36,6 @@ void DoPathLoss(char *filename, unsigned char geo, unsigned char kml, 255.0 / pow((double)(max_elevation - min_elevation), one_over_gamma); - width = (unsigned)(ippd * ReduceAngle(max_west - min_west)); - height = (unsigned)(ippd * ReduceAngle(max_north - min_north)); - LoadLossColors(xmtr[0]); if (filename[0] == 0) { @@ -117,11 +114,6 @@ void DoPathLoss(char *filename, unsigned char geo, unsigned char kml, loss = (dem[indx].signal[x0][y0]); cityorcounty = 0; - if (debug) { - fprintf(stdout, "\n%d\t%d\t%d\t%d", - loss, indx, x0, y0); - fflush(stdout); - } match = 255; red = 0; @@ -257,7 +249,7 @@ void DoSigStr(char *filename, unsigned char geo, unsigned char kml, points up and east points right in the image generated. */ char mapfile[255]; - unsigned width, height, terrain, red, green, blue; + unsigned terrain, red, green, blue; unsigned char found, mask, cityorcounty; int indx, x, y, z = 1, x0 = 0, y0 = 0, signal, match; double conversion, one_over_gamma, lat, lon, minwest; @@ -268,9 +260,6 @@ void DoSigStr(char *filename, unsigned char geo, unsigned char kml, 255.0 / pow((double)(max_elevation - min_elevation), one_over_gamma); - width = (unsigned)(ippd * ReduceAngle(max_west - min_west)); - height = (unsigned)(ippd * ReduceAngle(max_north - min_north)); - LoadSignalColors(xmtr[0]); if (filename[0] == 0) { @@ -344,13 +333,6 @@ void DoSigStr(char *filename, unsigned char geo, unsigned char kml, mask = dem[indx].mask[x0][y0]; signal = (dem[indx].signal[x0][y0]) - 100; cityorcounty = 0; - - if (debug) { - fprintf(stdout, "\n%d\t%d\t%d\t%d", - signal, indx, x0, y0); - fflush(stdout); - } - match = 255; red = 0; @@ -500,20 +482,18 @@ void DoRxdPwr(char *filename, unsigned char geo, unsigned char kml, points up and east points right in the image generated. */ char mapfile[255]; - unsigned width, height, terrain, red, green, blue; + unsigned terrain, red, green, blue; unsigned char found, mask, cityorcounty; int indx, x, y, z = 1, x0 = 0, y0 = 0, dBm, match; double conversion, one_over_gamma, lat, lon, minwest; FILE *fd; + one_over_gamma = 1.0 / GAMMA; conversion = 255.0 / pow((double)(max_elevation - min_elevation), one_over_gamma); - width = (unsigned)(ippd * ReduceAngle(max_west - min_west)); - height = (unsigned)(ippd * ReduceAngle(max_north - min_north)); - LoadDBMColors(xmtr[0]); if (filename[0] == 0) { @@ -553,15 +533,14 @@ void DoRxdPwr(char *filename, unsigned char geo, unsigned char kml, fd = fopen(mapfile, "wb"); - fprintf(fd, "P6\n%u %u\n255\n", width, (kml ? height : height + 30)); + fprintf(fd, "P6\n%u %u\n255\n", width, (kml ? height : height)); if (debug) { fprintf(stdout, "\nWriting \"%s\" (%ux%u pixmap image)... ", - mapfile, width, (kml ? height : height + 30)); + mapfile, width, (kml ? height : height)); fflush(stdout); } - // WriteKML() - //writeKML(xmtr,filename); + // Draw image of x by y pixels for (y = 0, lat = north; y < (int)height; y++, lat = north - (dpp * (double)y)) { for (x = 0, lon = max_west; x < (int)width; @@ -571,14 +550,13 @@ void DoRxdPwr(char *filename, unsigned char geo, unsigned char kml, for (indx = 0, found = 0; indx < MAXPAGES && found == 0;) { - x0 = (int)rint(ppd * - (lat - - (double)dem[indx].min_north)); + x0 = (int)rint((ppd * + (lat - + (double)dem[indx].min_north))); // +4549 fix y0 = mpi - - (int)rint(ppd * + (int)rint(ppd * // ONLY exception. All others are yppd (LonDiff - ((double)dem[indx].max_west, - lon))); + ((double)dem[indx].max_west,lon))); if (x0 >= 0 && x0 <= mpi && y0 >= 0 && y0 <= mpi) @@ -591,13 +569,6 @@ void DoRxdPwr(char *filename, unsigned char geo, unsigned char kml, mask = dem[indx].mask[x0][y0]; dBm = (dem[indx].signal[x0][y0]) - 200; cityorcounty = 0; - - if (debug) { - fprintf(stdout, "\n%d\t%d\t%d\t%d", dBm, - indx, x0, y0); - fflush(stdout); - } - match = 255; red = 0; @@ -639,9 +610,7 @@ void DoRxdPwr(char *filename, unsigned char geo, unsigned char kml, else if (mask & 4) { /* County Boundaries: Black */ - fprintf(fd, "%c%c%c", 0, 0, 0); - cityorcounty = 1; } @@ -690,7 +659,7 @@ void DoRxdPwr(char *filename, unsigned char geo, unsigned char kml, "%c%c%c", 255, 255, - 255); + 255); // WHITE else { if (dem[indx]. data[x0][y0] @@ -700,7 +669,7 @@ void DoRxdPwr(char *filename, unsigned char geo, unsigned char kml, "%c%c%c", 0, 0, - 170); + 170); // BLUE else { /* Elevation: Greyscale */ terrain @@ -732,6 +701,7 @@ void DoRxdPwr(char *filename, unsigned char geo, unsigned char kml, } } + fclose(fd); } @@ -746,7 +716,7 @@ void DoLOS(char *filename, unsigned char geo, unsigned char kml, points up and east points right in the image generated. */ char mapfile[255]; - unsigned width, height, terrain; + unsigned terrain; unsigned char found, mask; int indx, x, y, x0 = 0, y0 = 0; double conversion, one_over_gamma, lat, lon, minwest; @@ -757,9 +727,6 @@ void DoLOS(char *filename, unsigned char geo, unsigned char kml, 255.0 / pow((double)(max_elevation - min_elevation), one_over_gamma); - width = (unsigned)(ippd * ReduceAngle(max_west - min_west)); - height = (unsigned)(ippd * ReduceAngle(max_north - min_north)); - if (filename[0] == 0) { strncpy(filename, xmtr[0].filename, 254); filename[strlen(filename) - 4] = 0; /* Remove .qth */ @@ -803,8 +770,6 @@ void DoLOS(char *filename, unsigned char geo, unsigned char kml, mapfile, width, (kml ? height : height + 30)); fflush(stdout); } - // WriteKML() - //writeKML(xmtr,filename); for (y = 0, lat = north; y < (int)height; y++, lat = north - (dpp * (double)y)) {