Skip to content
GitLab
Commit 16daa957 authored by Stephan Korsholm's avatar Stephan Korsholm
Browse files

Merge branch '146-modifications-to-radiometeorological-data-maps-in-pmp-p-452-16-17' into 'master' 

Resolve "Modifications to radiometeorological Data Maps in PMP P.452-16/17"

Closes #146

See merge request SEAM/seamcat!47
parents 5b2023c5 6127981e
Hide whitespace changes
Inline Side-by-side
Showing
with 769 additions and 380 deletions
+769 -380
core/src/main/java/org/seamcat/Seamcat.java
+ 1
- 2
View file @ 16daa957
......@@ -126,11 +126,10 @@ public class Seamcat {
SandboxInitializer.initializeSandbox();
// initialize Factory
Factory.initialize(new DistributionFactoryImpl(), new PropagationModelFactoryImpl(new ClimateData()),
Factory.initialize(new DistributionFactoryImpl(), new PropagationModelFactoryImpl(new ClimateDataImpl()),
new BuildersImpl(), new AntennaGainFactoryImpl(), new FunctionFactoryImpl(),
new CoverageRadiusFactoryImpl(), new EventProcessingFactoryImpl(), new CorrelationModeFactoryImpl(),
new PluginFactoryImpl(), new ResultFactoryImpl(), new SystemHandlerImpl());
ClimateData.loadITUdata();
// default ui factory with limited capabilities
Factory.initializeUI(new LoggingUIFactory());
......
core/src/main/java/org/seamcat/model/propagation/ClimateDataImpl.java 0 → 100644
+ 121
- 0
View file @ 16daa957
package org.seamcat.model.propagation;
import java.io.BufferedReader;
import java.io.IOException;
import java.io.InputStream;
import java.io.InputStreamReader;
import java.util.ArrayList;
import java.util.List;
import org.seamcat.model.plugin.propagation.ClimateData;
import org.seamcat.model.plugin.propagation.TerrainCoordinate;
public class ClimateDataImpl implements ClimateData {
private List<String> deltaNData;
private List<String> n0Data;
/**
* loads the data from ITU-R F452-16 additional file (deltaN and N0)
*/
public ClimateDataImpl() {
// TODO get deltaN and N0 from ITU data and set to PMP
deltaNData = new ArrayList<>();
n0Data = new ArrayList<>();
String sDN50 = "/p452/DN50.TXT";
String sN0 = "/p452/N050.TXT";
readSource(sDN50, deltaNData);
readSource(sN0, n0Data);
}
private static void readSource(String name, List<String> target) {
InputStream is = ClimateDataImpl.class.getResourceAsStream(name);
try (BufferedReader br = new BufferedReader(new InputStreamReader(is))) {
String line = br.readLine();
while (line != null) {
target.add(line);
line = br.readLine();
}
} catch (IOException e) {
}
}
public double[] getClimateData(TerrainCoordinate point) {
ArrayList<Double> posLonInMap = new ArrayList<>();
ArrayList<Double> posLatInMap = new ArrayList<>();
int lon = (int) point.getLon();
int lat = (int) point.getLat();
for (double d = -90; d <= 90; d += 1.5) {
posLatInMap.add(d);
}
for (double l = 0; l <= 360; l += 1.5) {
posLonInMap.add(l);
}
int index = -1; // takes account of the first ++
while (lat <= posLatInMap.get(++index)) {
}
int indexLat0 = index;
int indexLat1 = index + 1;
if (lat < 0)
indexLat1 = index - 1;
index = -1;
if (lon < 0)
lon = 360 + lon;
while (lon > posLonInMap.get(++index)) {
}
int indexLon0 = index;
int indexLon1 = index + 1;
if (lon < 0)
indexLon1 = index - 1;
double[] bilinear = new double[4];
int offset = 9 * indexLon0; // each entry has together with its leading spaces 9 digits
if (offset + 9 > 2169) {
offset = 0; // TODO was at the end of the line at 359 degree -> take the value at zero degree
}
bilinear[0] = Double.parseDouble(deltaNData.get(indexLat0).substring(offset, offset + 9)); // C / R
bilinear[2] = Double.parseDouble(deltaNData.get(indexLat1).substring(offset, offset + 9)); // C / R + 1
offset = 9 * indexLon1;
if (offset + 9 > 2169) {
offset = 0; // TODO was at the end of the line
}
bilinear[1] = Double.parseDouble(deltaNData.get(indexLat0).substring(offset, offset + 9)); // C + 1 / R
bilinear[3] = Double.parseDouble(deltaNData.get(indexLat1).substring(offset, offset + 9)); // C + 1 / R + 1
double rDeltaN = doBilinear(bilinear, point, lat, lon);
offset = 9 * indexLon0; // each entry has together with its leading spaces 9 digits
if (offset + 9 > 2169) {
offset = 0; // TODO was at the end of the line at 359 degree -> take the value at zero degree
}
bilinear[0] = Double.parseDouble(n0Data.get(indexLat0).substring(offset, offset + 9)); // C / R
bilinear[2] = Double.parseDouble(n0Data.get(indexLat1).substring(offset, offset + 9)); // C / R + 1
offset = 9 * indexLon1;
if (offset + 9 > 2169) {
offset = 0; // TODO was at the end of the line at 359 degree -> take the value at zero degree
}
bilinear[1] = Double.parseDouble(n0Data.get(indexLat0).substring(offset, offset + 9)); // C + 1 / R
bilinear[3] = Double.parseDouble(n0Data.get(indexLat1).substring(offset, offset + 9)); // C + 1 / R + 1
double rN0 = doBilinear(bilinear, point, lat, lon);
return new double[] {rDeltaN, rN0};
}
private double doBilinear(double[] bilinear, TerrainCoordinate point, int lat, int lon) {
int C = lon;
int C1 = C + 1;
int R = lat;
int R1 = R + 1;
double p0 = point.getLat();
double p1 = point.getLon();
return bilinear[0] * (R1 - p0) * (C1 - p1) + bilinear[2] * (p0 - R) * (C1 - p1)
+ bilinear[1] * (R1 - p0) * (p1 - C) + bilinear[3] * (p0 - R) * (p1 - C);
}
}
core/src/main/java/org/seamcat/model/propagation/DigitalMap.java 0 → 100644
+ 103
- 0
View file @ 16daa957
package org.seamcat.model.propagation;
import java.io.BufferedReader;
import java.io.InputStream;
import java.io.InputStreamReader;
import java.util.ArrayList;
import java.util.List;
public class DigitalMap {
double[][] _map;
int _sizeX, _sizeY;
private double _spacing;
public DigitalMap(String path, boolean lastAndFirstColumnEqual) {
String name = path;
try {
InputStream inputStream = getClass().getResourceAsStream(name);
List<String> lines = new ArrayList<>();
BufferedReader br = new BufferedReader(new InputStreamReader(inputStream));
String line;
while (null != (line = br.readLine())) {
lines.add(line);
}
_sizeY = lines.size();
_sizeX = ParseLine(lines.get(0)).length;
_map = new double[_sizeY][_sizeX];
for (int i = 0; i < _sizeY; i++) { /* DO */
double[] data = ParseLine(lines.get(i));
for (int j = 0; j < _sizeX; j++) {
_map[i][j] = data[j];
}
}
} catch (Exception ex) {
throw new IllegalArgumentException("Could not load map: '" + name + "'");
}
if (lastAndFirstColumnEqual) {
_spacing = 360.0 / (_map[1].length - 1);
} else {
_spacing = 360.0 / (_map[1].length);
}
}
private double[] ParseLine(String line) {
String[] parts = line.trim().split("\\s+");
;
double[] data = new double[parts.length];
for (int i = 0; i < parts.length; i++) {
data[i] = Double.parseDouble(parts[i]);
// System.out.printf( "parts.length = %f\n" ,data[i]);
}
return data;
}
public double GetClosestGridPointValue(double longitude, double latitude) {
double longitudeOffset = longitude + 180.0;
double latitudeOffset = 90.0 - latitude;
int latitudeIndex = (int) (latitudeOffset / _spacing);
int longitudeIndex = (int) (longitudeOffset / _spacing);
latitudeIndex %= _sizeY;
longitudeIndex %= _sizeX;
double val = _map[latitudeIndex][longitudeIndex];
return val;
}
public double GetInterpolatedValue(double longitude, double latitude) {
double longitudeOffset = longitude;
if (longitude < 0.0) {
longitudeOffset = longitude + 360.0;
}
double latitudeOffset = 90.0 - latitude;
int latitudeIndex = (int) (latitudeOffset / _spacing);
int longitudeIndex = (int) (longitudeOffset / _spacing);
// System.out.printf( "latitudeIndex = %d\n" ,latitudeIndex);
// System.out.printf( "longitudeIndex = %d\n" ,longitudeIndex);
double latitudeFraction = (latitudeOffset / _spacing) - latitudeIndex;
double longitudeFraction = (longitudeOffset / _spacing) - longitudeIndex;
double value_ul = _map[latitudeIndex][longitudeIndex];
double value_ur = _map[latitudeIndex][(longitudeIndex + 1) % _sizeX];
double value_ll = _map[(latitudeIndex + 1) % _sizeY][longitudeIndex];
double value_lr = _map[(latitudeIndex + 1) % _sizeY][(longitudeIndex + 1) % _sizeX];
double interpolatedHeight1 = (longitudeFraction * (value_ur - value_ul)) + value_ul;
double interpolatedHeight2 = (longitudeFraction * (value_lr - value_ll)) + value_ll;
double interpolatedHeight3 =
latitudeFraction * (interpolatedHeight2 - interpolatedHeight1) + interpolatedHeight1;
return interpolatedHeight3;
}
}
\ No newline at end of file
core/src/main/java/org/seamcat/model/propagation/P2001ver2DigitalMaps.java
+ 17
- 116
View file @ 16daa957
package org.seamcat.model.propagation;
import java.io.BufferedReader;
import java.io.InputStream;
import java.io.InputStreamReader;
import java.util.ArrayList;
import java.util.List;
import org.seamcat.model.propagation.DigitalMap;
public class P2001ver2DigitalMaps {
private DigitalMap _mapSdN;
......@@ -23,20 +19,22 @@ public class P2001ver2DigitalMaps {
private DigitalMap _mapTropoClim;
public P2001ver2DigitalMaps() {
_mapSdN = new DigitalMap("DN_Median.txt", true);
_mapSdNsup = new DigitalMap("DN_SupSlope.txt", true);
_mapSdNsub = new DigitalMap("DN_SubSlope.txt", true);
_mapDndz_01 = new DigitalMap("dndz_01.txt", true);
_mapRho_sur = new DigitalMap("surfwv_50_fixed.txt", true);
_mapEsarain_Pr6 = new DigitalMap("Esarain_Pr6_v5.txt", true);
_mapEsarain_Mt = new DigitalMap("Esarain_Mt_v5.txt", true);
_mapEsarain_Beta = new DigitalMap("Esarain_Beta_v5.txt", true);
_mapData_h0 = new DigitalMap("h0.txt", true);
_mapFoEs50 = new DigitalMap("FoEs50.txt", true);
_mapFoEs10 = new DigitalMap("FoEs10.txt", true);
_mapFoEs01 = new DigitalMap("FoEs01.txt", true);
_mapFoEs0_1 = new DigitalMap("FoEs0.1.txt", true);
_mapTropoClim = new DigitalMap("TropoClim.txt", false);
String rec = "/p2001/";
_mapSdN = new DigitalMap(rec + "DN_Median.txt", true);
_mapSdNsup = new DigitalMap(rec + "DN_SupSlope.txt", true);
_mapSdNsub = new DigitalMap(rec + "DN_SubSlope.txt", true);
_mapDndz_01 = new DigitalMap(rec + "dndz_01.txt", true);
_mapRho_sur = new DigitalMap(rec + "surfwv_50_fixed.txt", true);
_mapEsarain_Pr6 = new DigitalMap(rec + "Esarain_Pr6_v5.txt", true);
_mapEsarain_Mt = new DigitalMap(rec + "Esarain_Mt_v5.txt", true);
_mapEsarain_Beta = new DigitalMap(rec + "Esarain_Beta_v5.txt", true);
_mapData_h0 = new DigitalMap(rec + "h0.txt", true);
_mapFoEs50 = new DigitalMap(rec + "FoEs50.txt", true);
_mapFoEs10 = new DigitalMap(rec + "FoEs10.txt", true);
_mapFoEs01 = new DigitalMap(rec + "FoEs01.txt", true);
_mapFoEs0_1 = new DigitalMap(rec + "FoEs0.1.txt", true);
_mapTropoClim = new DigitalMap(rec + "TropoClim.txt", false);
}
public double GetSdN(double lon, double lat) {
......@@ -94,101 +92,4 @@ public class P2001ver2DigitalMaps {
public int GetClimZone(double lon, double lat) {
return (int) _mapTropoClim.GetClosestGridPointValue(lon, lat);
}
class DigitalMap {
double[][] _map;
int _sizeX, _sizeY;
private double _spacing;
public DigitalMap(String path, boolean lastAndFirstColumnEqual) {
String name = "/p2001/" + path;
try {
InputStream inputStream = getClass().getResourceAsStream(name);
List<String> lines = new ArrayList<>();
BufferedReader br = new BufferedReader(new InputStreamReader(inputStream));
String line;
while (null != (line = br.readLine())) {
lines.add(line);
}
_sizeY = lines.size();
_sizeX = ParseLine(lines.get(0)).length;
_map = new double[_sizeY][_sizeX];
for (int i = 0; i < _sizeY; i++) { /* DO */
double[] data = ParseLine(lines.get(i));
for (int j = 0; j < _sizeX; j++) {
_map[i][j] = data[j];
}
}
} catch (Exception ex) {
throw new IllegalArgumentException("Could not load map: '" + name + "'");
}
if (lastAndFirstColumnEqual) {
_spacing = 360.0 / (_map[1].length - 1);
} else {
_spacing = 360.0 / (_map[1].length);
}
}
private double[] ParseLine(String line) {
String[] parts = line.trim().split("\\s+");
;
double[] data = new double[parts.length];
for (int i = 0; i < parts.length; i++) {
data[i] = Double.parseDouble(parts[i]);
// System.out.printf( "parts.length = %f\n" ,data[i]);
}
return data;
}
public double GetClosestGridPointValue(double longitude, double latitude) {
double longitudeOffset = longitude + 180.0;
double latitudeOffset = 90.0 - latitude;
int latitudeIndex = (int) (latitudeOffset / _spacing);
int longitudeIndex = (int) (longitudeOffset / _spacing);
latitudeIndex %= _sizeY;
longitudeIndex %= _sizeX;
double val = _map[latitudeIndex][longitudeIndex];
return val;
}
public double GetInterpolatedValue(double longitude, double latitude) {
double longitudeOffset = longitude;
if (longitude < 0.0) {
longitudeOffset = longitude + 360.0;
}
double latitudeOffset = 90.0 - latitude;
int latitudeIndex = (int) (latitudeOffset / _spacing);
int longitudeIndex = (int) (longitudeOffset / _spacing);
// System.out.printf( "latitudeIndex = %d\n" ,latitudeIndex);
// System.out.printf( "longitudeIndex = %d\n" ,longitudeIndex);
double latitudeFraction = (latitudeOffset / _spacing) - latitudeIndex;
double longitudeFraction = (longitudeOffset / _spacing) - longitudeIndex;
double value_ul = _map[latitudeIndex][longitudeIndex];
double value_ur = _map[latitudeIndex][(longitudeIndex + 1) % _sizeX];
double value_ll = _map[(latitudeIndex + 1) % _sizeY][longitudeIndex];
double value_lr = _map[(latitudeIndex + 1) % _sizeY][(longitudeIndex + 1) % _sizeX];
double interpolatedHeight1 = (longitudeFraction * (value_ur - value_ul)) + value_ul;
double interpolatedHeight2 = (longitudeFraction * (value_lr - value_ll)) + value_ll;
double interpolatedHeight3 =
latitudeFraction * (interpolatedHeight2 - interpolatedHeight1) + interpolatedHeight1;
return interpolatedHeight3;
}
}
}
\ No newline at end of file
core/src/main/java/org/seamcat/model/propagation/P452ver16PropagationModel.java
+ 186
- 21
View file @ 16daa957
......@@ -25,8 +25,10 @@ public class P452ver16PropagationModel
// Removed consistency checks and exceptions related to this
// distance
// v3 19JUL22 Ivica Stevanovic, OFCOM Introduced polarization as input argument in PMP UI
// v4 09MAR23 Ivica Stevanovic, OFCOM Introduced DN and N0 maps and midpoint calculation for terrain
// profile
//
// Copyright (c) 2016-2017, Ivica Stevanovic
// Copyright (c) 2016-, Ivica Stevanovic
// All rights reserved.
//
// Redistribution and use in source and binary forms, with or without
......@@ -56,6 +58,8 @@ public class P452ver16PropagationModel
// THE AUTHORS AND OFCOM (CH) DO NOT PROVIDE ANY SUPPORT FOR THIS SOFTWARE
////
private P452ver17DigitalMaps maps;
@Override
public void consistencyCheck(ConsistencyCheckContext context, P452ver16Input input) {
Distribution frequency = context.getFrequency();
......@@ -108,6 +112,9 @@ public class P452ver16PropagationModel
@Override
public double evaluate(LinkResult linkResult, boolean variations, P452ver16Input input) {
if (maps == null)
maps = new P452ver17DigitalMaps();
double rFreq = linkResult.getFrequency();
double rHTx = linkResult.txAntenna().getHeight();
double rHRx = linkResult.rxAntenna().getHeight();
......@@ -148,10 +155,16 @@ public class P452ver16PropagationModel
double[] rDisti;
double[] rHighti;
int[] iZonei;
double refractionIndex;
double seaLevelSurfaceRefractivity;
// double refractionIndex;
// double seaLevelSurfaceRefractivity;
double N0;
double DN;
double txLatitude;
double rxLatitude;
double txLongitude;
double rxLongitude;
double lon_mid;
double lat_mid;
TerrainData terrainData = linkResult.getTerrainData();
if (terrainData != null) {
......@@ -168,10 +181,31 @@ public class P452ver16PropagationModel
iZonei[i] = dataPoint.getZone();
}
refractionIndex = terrainData.getRefractionIndex();
seaLevelSurfaceRefractivity = terrainData.getSeaLevelSurfaceRefractivity();
// refractionIndex = terrainData.getRefractionIndex();
// seaLevelSurfaceRefractivity = terrainData.getSeaLevelSurfaceRefractivity();
rxLatitude = terrainData.getRxLatitude();
txLatitude = terrainData.getTxLatitude();
rxLongitude = terrainData.getTo().getLon();
txLongitude = terrainData.getFrom().getLon();
// Calculate the longitude and latitude of the mid-point o the path, lon_mid,
// and lat_mid for dpnt = 0.5
// Use the method from P.2001 Annex H for computation of path center along the great circle path
double Re = 6371;
double dpnt = 0.5;
double[] gcp = great_circle_path(rxLongitude, txLongitude, rxLatitude, txLatitude, Re, dpnt);
lon_mid = gcp[0];
lat_mid = gcp[1];
// The following input arguments related to meteorological data:
// DN: The average radio-refractivity lapse-rate through the lowest 1 km of the atmosphere at the path
// center N0: The sea-level surface refractivity at the path center
DN = maps.GetDN50(lon_mid, lat_mid); // interpolate the value of DN50 at the point given by (lon, lat)
N0 = maps.GetN050(lon_mid, lat_mid); // interpolate the value of N050 at the point given by (lon,
} else {
double rDist = linkResult.getTxRxDistance();
......@@ -206,10 +240,14 @@ public class P452ver16PropagationModel
iZonei[i] = 2;
}
refractionIndex = input.refractionIndex();
seaLevelSurfaceRefractivity = input.seaLevelSurfaceRefractivity();
txLatitude = input.latitude();
rxLatitude = input.latitude();
// refractionIndex = input.refractionIndex();
// seaLevelSurfaceRefractivity = input.seaLevelSurfaceRefractivity();
DN = input.refractionIndex();
N0 = input.seaLevelSurfaceRefractivity();
// txLatitude = input.latitude();
// rxLatitude = input.latitude();
lat_mid = input.latitude();
}
double rpol = 1; //(1 = horizontal polarization, 2 = vertical);
......@@ -222,10 +260,9 @@ public class P452ver16PropagationModel
}
double dct = 500; // for land only paths
double dcr = 500;
rL = tl_p452(rFreq / 1000, rP, rDisti, rHighti, iZonei, rHTx, rHRx, txLatitude, rxLatitude,
linkResult.txAntenna().getGain(), linkResult.rxAntenna().getGain(), rpol, dct, dcr, refractionIndex,
seaLevelSurfaceRefractivity, input.surfacePressure(), input.surfaceTemperature(), ha_t, ha_r, dk_t, dk_r,
input.diffractionOnly());
rL = tl_p452(rFreq / 1000, rP, rDisti, rHighti, iZonei, rHTx, rHRx, lat_mid, linkResult.txAntenna().getGain(),
linkResult.rxAntenna().getGain(), rpol, dct, dcr, DN, N0, input.surfacePressure(),
input.surfaceTemperature(), ha_t, ha_r, dk_t, dk_r, input.diffractionOnly());
// calculation of the standard deviation
if (variations) {
......@@ -286,12 +323,12 @@ public class P452ver16PropagationModel
return rValue;
}
public double tl_p452(double f, double p, double[] d, double[] h, int[] zone, double htg, double hrg, double phi_t,
double phi_r, double Gt, double Gr, double pol, double dct, double dcr, double DN, double N0, double press,
public double tl_p452(double f, double p, double[] d, double[] h, int[] zone, double htg, double hrg,
double phi_path, double Gt, double Gr, double pol, double dct, double dcr, double DN, double N0, double press,
double temp, double ha_t, double ha_r, double dk_t, double dk_r, boolean diffractionOnly) {
// tl_p452 basic transmission loss according to P.452-16
// Lb = tl_p452(f, p, d, h, zone, htg, hrg, phi_t, phi_r, Gt, Gr, pol, dct, dcr, DN, N0, press, temp, ha_t,
// ha_r, dk_t, dk_r, diffraction_only )
// Lb = tl_p452(f, p, d, h, zone, htg, hrg, phi_path, Gt, Gr, pol, dct, dcr, DN, N0, press, temp, ha_t, ha_r,
// dk_t, dk_r, diffraction_only )
//
// This is the MAIN function that computes the basic transmission loss not exceeded for p% of time
// as defined in ITU-R P.452-16 (Section 4.6).
......@@ -306,8 +343,7 @@ public class P452ver16PropagationModel
// zone - Zone type: Coastal land (1), Inland (2) or Sea (3)
// htg - Tx Antenna center heigth above ground level (m)
// hrg - Rx Antenna center heigth above ground level (m)
// phi_t - Latitude of Tx station (degrees)
// phi_r - Latitude of Rx station (degrees)
// phi_path- Latitude of the mid point between Tx and Rx stations along the great circle (degrees)
// Gt, Gr - Antenna gain in the direction of the horizon along the
// great-circle interference path (dBi)
// pol - polarization of the signal (1) horizontal, (2) vertical
......@@ -363,7 +399,7 @@ public class P452ver16PropagationModel
// Compute the path profile parameters
// Path center latitude
double phi_path = (phi_t + phi_r) / 2;
// double phi_path = (phi_t + phi_r) / 2;
// Compute dtm - the longest continuous land (inland + coastal) section of the great-circle path (km)
zone_r = 12;
......@@ -2281,6 +2317,136 @@ public class P452ver16PropagationModel
return out;
}
public double[] great_circle_path(double Phire, double Phite, double Phirn, double Phitn, double Re, double dpnt) {
// great_circle_path Great-circle path calculations according to Attachment H
// This function computes the great-circle intermediate points on the
// radio path as defined in ITU-R P.2001-3 Attachment H
//
// Input parameters:
// Phire - Receiver longitude, positive to east (deg)
// Phite - Transmitter longitude, positive to east (deg)
// Phirn - Receiver latitude, positive to north (deg)
// Phitn - Transmitter latitude, positive to north (deg)
// Re - Average Earth radius (km)
// dpnt - Distance from the transmitter to the intermediate point expressed as a fraction of great
// circle distance (dgc)
//
// Output parameters:
// Phipnte - Longitude of the intermediate point (deg)
// Phipntn - Latitude of the intermediate point (deg)
// Bt2r - Bearing of the great-circle path from Tx towards the Rx (deg)
// dgc - Great-circle path length (km)
//
//
// Rev Date Author Description
// -------------------------------------------------------------------------------
// v0 12JUL16 Ivica Stevanovic, OFCOM Initial version
// v1 13JUL17 Ivica Stevanovic, OFCOM Initial version in Java
// v2 27JAN23 Ivica Stevanovic, OFCOM modified dpnt to be a fraction od dgc (and not
// absolute distance)
//// H.2 Path length and bearing
// Difference (deg) in longitude between the terminals (H.2.1)
double Dlon = Phire - Phite;
// Calculate quantity r (H.2.2)
double r = sind(Phitn) * sind(Phirn) + cosd(Phitn) * cosd(Phirn) * cosd(Dlon);
// Calculate the path length as the angle subtended at the center of
// average-radius Earth (H.2.3)
double Phid = Math.acos(r); // radians
// Calculate the great-circle path length (H.2.4)
double dgc = Phid * Re; // km
// Calculate the quantity x1 (H.2.5a)
double x1 = sind(Phirn) - r * sind(Phitn);
// Calculate the quantity y1 (H.2.5b)
double y1 = cosd(Phitn) * cosd(Phirn) * sind(Dlon);
// Calculate the bearing of the great-circle path for Tx to Rx (H.2.6)
double Bt2r = Phire;
if (Math.abs(x1) < 1e-9 && Math.abs(y1) < 1e-9) {
Bt2r = Phire;
} else {
Bt2r = atan2d(y1, x1);
}
//// H.3 Calculation of intermediate path point
// Calculate the distance to the point as the angle subtended at the center
// of average-radius Earth (H.3.1)
double Phipnt = dpnt * dgc / Re; // radians
// Calculate quantity s (H.3.2)
double s = sind(Phitn) * Math.cos(Phipnt) + cosd(Phitn) * Math.sin(Phipnt) * cosd(Bt2r);
// The latitude of the intermediate point is now given by (H.3.3)
double Phipntn = asind(s); // degs
// Calculate the quantity x2 (H.3.4a)
double x2 = Math.cos(Phipnt) - s * sind(Phitn);
// Calculate the quantity y2 (H.3.4b)
double y2 = cosd(Phitn) * Math.sin(Phipnt) * sind(Bt2r);
// Calculate the longitude of the intermediate point Phipnte (H.3.5)
double Phipnte = Bt2r;
if (x2 < 1e-9 && y2 < 1e-9) {
Phipnte = Bt2r;
} else {
Phipnte = Phite + atan2d(y2, x2);
}
double[] out = new double[4];
out[0] = Phipnte;
out[1] = Phipntn;
out[2] = Bt2r;
out[3] = dgc;
return out;
}
public double sind(double theta_deg) {
double theta_rad = theta_deg * Math.PI / 180.0;
double y = Math.sin(theta_rad);
return y;
}
public double cosd(double theta_deg) {
double theta_rad = theta_deg * Math.PI / 180.0;
double y = Math.cos(theta_rad);
return y;
}
public double atan2d(double y, double x) {
double res = Math.atan2(y, x) * 180.0 / Math.PI;
return res;
}
public double asind(double y) {
double x = Math.asin(y) * 180.0 / Math.PI;
return x;
}
@Override
public Description description() {
String text = "<html>" + DescriptionTags.startTag.div();
......@@ -2300,7 +2466,6 @@ public class P452ver16PropagationModel
"<br><b><u>Note 3:</u></b> SRTM DEMs do not contain only the bare terrain profile, but may include any objects on the ground (clutter).";
text +=
"<br><b><u>Note 4:</u></b> Some tiles of SRTM DEM may contain voids in regions outside the US. When constructing terrain profile in SEAMCAT, heights in the profile having void are replaced with the last valid height from the profile.";
text += "<br><b><u>Note 5:</u></b> Implementation uses Inland zone type for Radio climatic zones settings.";
text += DescriptionTags.endTagNotes.div();
text += DescriptionTags.endTag.div() + "</html>";
return new DescriptionImpl("ITU-R P.452-16", text);
......
core/src/main/java/org/seamcat/model/propagation/P452ver17DigitalMaps.java 0 → 100644
+ 21
- 0
View file @ 16daa957
package org.seamcat.model.propagation;
public class P452ver17DigitalMaps {
private DigitalMap _mapDN50;
private DigitalMap _mapN050;
public P452ver17DigitalMaps() {
String rec = "/p452/";
_mapN050 = new DigitalMap(rec + "N050.TXT", true);
_mapDN50 = new DigitalMap(rec + "DN50.TXT", true);
}
public double GetDN50(double lon, double lat) {
return _mapDN50.GetInterpolatedValue(lon, lat);
}
public double GetN050(double lon, double lat) {
return _mapN050.GetInterpolatedValue(lon, lat);
}
}
\ No newline at end of file
core/src/main/java/org/seamcat/model/propagation/P452ver17PropagationModel.java
+ 187
- 26
View file @ 16daa957
......@@ -25,6 +25,8 @@ public class P452ver17PropagationModel
// Removed consistency checks and exceptions related to this
// distance
// v3 25OCT21 Ivica Stevanovic, OFCOM Aligned to P.452-17 (free space alignment with P.525)
// v4 27JAN23 Ivica Stevanovic, OFCOM Included path center calculation according to ITU-R P.2001
// Annex H
//
// Copyright (c) 2017- , Ivica Stevanovic
// All rights reserved.
......@@ -56,6 +58,8 @@ public class P452ver17PropagationModel
// THE AUTHORS AND OFCOM (CH) DO NOT PROVIDE ANY SUPPORT FOR THIS SOFTWARE
////
private P452ver17DigitalMaps maps;
@Override
public void consistencyCheck(ConsistencyCheckContext context, P452ver17Input input) {
Distribution frequency = context.getFrequency();
......@@ -108,6 +112,9 @@ public class P452ver17PropagationModel
@Override
public double evaluate(LinkResult linkResult, boolean variations, P452ver17Input input) {
if (maps == null)
maps = new P452ver17DigitalMaps();
double rFreq = linkResult.getFrequency();
double rDist = linkResult.getTxRxDistance();
double rHTx = linkResult.txAntenna().getHeight();
......@@ -149,10 +156,16 @@ public class P452ver17PropagationModel
double[] rDisti;
double[] rHighti;
int[] iZonei;
double refractionIndex;
double seaLevelSurfaceRefractivity;
// double refractionIndex;
// double seaLevelSurfaceRefractivity;
double N0;
double DN;
double txLatitude;
double rxLatitude;
double txLongitude;
double rxLongitude;
double lon_mid;
double lat_mid;
TerrainData terrainData = linkResult.getTerrainData();
if (terrainData != null) {
......@@ -169,10 +182,29 @@ public class P452ver17PropagationModel
iZonei[i] = dataPoint.getZone();
}
refractionIndex = terrainData.getRefractionIndex();
seaLevelSurfaceRefractivity = terrainData.getSeaLevelSurfaceRefractivity();
// refractionIndex = terrainData.getRefractionIndex();
// seaLevelSurfaceRefractivity = terrainData.getSeaLevelSurfaceRefractivity();
rxLatitude = terrainData.getRxLatitude();
txLatitude = terrainData.getTxLatitude();
rxLongitude = terrainData.getTo().getLon();
txLongitude = terrainData.getFrom().getLon();
// Use the method from P.2001 Annex H for computation of path center along the great circle path
double Re = 6371;
double dpnt = 0.5;
double[] gcp = great_circle_path(rxLongitude, txLongitude, rxLatitude, txLatitude, Re, dpnt);
lon_mid = gcp[0];
lat_mid = gcp[1];
// The following input arguments related to meteorological data:
// DN: The average radio-refractivity lapse-rate through the lowest 1 km of the atmosphere at the path
// center N0: The sea-level surface refractivity at the path center
DN = maps.GetDN50(lon_mid, lat_mid); // interpolate the value of DN50 at the point given by (lon, lat)
N0 = maps.GetN050(lon_mid, lat_mid); // interpolate the value of N050 at the point given by (lon, lat)
} else {
int rTotalLength = Math.max(5, (int) Math.round(rDist / rStep));
rStep = rDist / (rTotalLength - 1); // adjust rStep so that d(0)=0 and d(end)=rDist !!!
......@@ -180,10 +212,15 @@ public class P452ver17PropagationModel
rHighti = new double[rTotalLength];
iZonei = new int[rTotalLength];
refractionIndex = input.refractionIndex();
seaLevelSurfaceRefractivity = input.seaLevelSurfaceRefractivity();
txLatitude = input.latitude();
rxLatitude = input.latitude();
// refractionIndex = input.refractionIndex();
// seaLevelSurfaceRefractivity = input.seaLevelSurfaceRefractivity();
DN = input.refractionIndex();
N0 = input.seaLevelSurfaceRefractivity();
// txLatitude = input.latitude();
// rxLatitude = input.latitude();
lat_mid = input.latitude();
rDisti[0] = 0;
rHighti[0] = 0; // for flat only terrain
......@@ -223,10 +260,9 @@ public class P452ver17PropagationModel
}
double dct = 500; // for land only paths
double dcr = 500;
rL = tl_p452(rFreq / 1000, rP, rDisti, rHighti, iZonei, rHTx, rHRx, txLatitude, rxLatitude,
linkResult.txAntenna().getGain(), linkResult.rxAntenna().getGain(), rpol, dct, dcr, refractionIndex,
seaLevelSurfaceRefractivity, input.surfacePressure(), input.surfaceTemperature(), ha_t, ha_r, dk_t, dk_r,
input.diffractionOnly());
rL = tl_p452(rFreq / 1000, rP, rDisti, rHighti, iZonei, rHTx, rHRx, lat_mid, linkResult.txAntenna().getGain(),
linkResult.rxAntenna().getGain(), rpol, dct, dcr, DN, N0, input.surfacePressure(),
input.surfaceTemperature(), ha_t, ha_r, dk_t, dk_r, input.diffractionOnly());
// calculation of the standard deviation
if (variations) {
......@@ -287,12 +323,12 @@ public class P452ver17PropagationModel
return rValue;
}
public double tl_p452(double f, double p, double[] d, double[] h, int[] zone, double htg, double hrg, double phi_t,
double phi_r, double Gt, double Gr, double pol, double dct, double dcr, double DN, double N0, double press,
public double tl_p452(double f, double p, double[] d, double[] h, int[] zone, double htg, double hrg,
double phi_path, double Gt, double Gr, double pol, double dct, double dcr, double DN, double N0, double press,
double temp, double ha_t, double ha_r, double dk_t, double dk_r, boolean diffractionOnly) {
// tl_p452 basic transmission loss according to P.452-17
// Lb = tl_p452(f, p, d, h, zone, htg, hrg, phi_t, phi_r, Gt, Gr, pol, dct, dcr, DN, N0, press, temp, ha_t,
// ha_r, dk_t, dk_r, diffraction_only )
// Lb = tl_p452(f, p, d, h, zone, htg, hrg, phi_path, Gt, Gr, pol, dct, dcr, DN, N0, press, temp, ha_t, ha_r,
// dk_t, dk_r, diffraction_only )
//
// This is the MAIN function that computes the basic transmission loss not exceeded for p% of time
// as defined in ITU-R P.452-17 (Section 4.6).
......@@ -307,8 +343,7 @@ public class P452ver17PropagationModel
// zone - Zone type: Coastal land (1), Inland (2) or Sea (3)
// htg - Tx Antenna center heigth above ground level (m)
// hrg - Rx Antenna center heigth above ground level (m)
// phi_t - Latitude of Tx station (degrees)
// phi_r - Latitude of Rx station (degrees)
// phi_path- Latitude of the mid point between Tx and Rx stations along great circle path (degrees)
// Gt, Gr - Antenna gain in the direction of the horizon along the
// great-circle interference path (dBi)
// pol - polarization of the signal (1) horizontal, (2) vertical
......@@ -335,8 +370,8 @@ public class P452ver17PropagationModel
// Lb - basic transmission loss according to P.452-16
//
// Example:
// Lb = tl_p452(f, p, d, h, zone, htg, hrg, phi_t, phi_r, Gt, Gr, pol, dct, dcr, DN, N0, press, temp, ha_t,
// ha_r, dk_t, dk_r, diffraction_only)
// Lb = tl_p452(f, p, d, h, zone, htg, hrg, phi_t, phi_r, lon_t, lon_r, Gt, Gr, pol, dct, dcr, DN, N0,
// press, temp, ha_t, ha_r, dk_t, dk_r, diffraction_only)
// Rev Date Author Description
// -------------------------------------------------------------------------------
......@@ -347,7 +382,8 @@ public class P452ver17PropagationModel
// v4 17NOV16 Ivica Stevanovic, OFCOM Translated to Java
// v5 01MAY17 Ivica Stevanovic, OFCOM corrected bug in dl_bull (Srim computation)
// v6 20JUL18 Ivica Stevanovic, OFCOM introduced 3-D path length for shorter distances in
// free space model v7 19JUL22 Ivica Stevanovic, OFCOM Introduced polarization as input
// free space model v7 19JUL22 Ivica Stevanovic, OFCOM Introduced polarization as input
// argument v8 08MAR23 Ivica Stevanovic, OFCOM Introduced mid point latitude as input
// argument
//
......@@ -363,10 +399,6 @@ public class P452ver17PropagationModel
//
int zone_r;
// Compute the path profile parameters
// Path center latitude
double phi_path = (phi_t + phi_r) / 2;
// Compute dtm - the longest continuous land (inland + coastal) section of the great-circle path (km)
zone_r = 12;
double dtm = longest_cont_dist(d, zone, zone_r);
......@@ -2285,6 +2317,136 @@ public class P452ver17PropagationModel
return out;
}
public double[] great_circle_path(double Phire, double Phite, double Phirn, double Phitn, double Re, double dpnt) {
// great_circle_path Great-circle path calculations according to Attachment H
// This function computes the great-circle intermediate points on the
// radio path as defined in ITU-R P.2001-3 Attachment H
//
// Input parameters:
// Phire - Receiver longitude, positive to east (deg)
// Phite - Transmitter longitude, positive to east (deg)
// Phirn - Receiver latitude, positive to north (deg)
// Phitn - Transmitter latitude, positive to north (deg)
// Re - Average Earth radius (km)
// dpnt - Distance from the transmitter to the intermediate point expressed as a fraction of great
// circle distance (dgc)
//
// Output parameters:
// Phipnte - Longitude of the intermediate point (deg)
// Phipntn - Latitude of the intermediate point (deg)
// Bt2r - Bearing of the great-circle path from Tx towards the Rx (deg)
// dgc - Great-circle path length (km)
//
//
// Rev Date Author Description
// -------------------------------------------------------------------------------
// v0 12JUL16 Ivica Stevanovic, OFCOM Initial version
// v1 13JUL17 Ivica Stevanovic, OFCOM Initial version in Java
// v2 27JAN23 Ivica Stevanovic, OFCOM modified dpnt to be a fraction od dgc (and not
// absolute distance)
//// H.2 Path length and bearing
// Difference (deg) in longitude between the terminals (H.2.1)
double Dlon = Phire - Phite;
// Calculate quantity r (H.2.2)
double r = sind(Phitn) * sind(Phirn) + cosd(Phitn) * cosd(Phirn) * cosd(Dlon);
// Calculate the path length as the angle subtended at the center of
// average-radius Earth (H.2.3)
double Phid = Math.acos(r); // radians
// Calculate the great-circle path length (H.2.4)
double dgc = Phid * Re; // km
// Calculate the quantity x1 (H.2.5a)
double x1 = sind(Phirn) - r * sind(Phitn);
// Calculate the quantity y1 (H.2.5b)
double y1 = cosd(Phitn) * cosd(Phirn) * sind(Dlon);
// Calculate the bearing of the great-circle path for Tx to Rx (H.2.6)
double Bt2r = Phire;
if (Math.abs(x1) < 1e-9 && Math.abs(y1) < 1e-9) {
Bt2r = Phire;
} else {
Bt2r = atan2d(y1, x1);
}
//// H.3 Calculation of intermediate path point
// Calculate the distance to the point as the angle subtended at the center
// of average-radius Earth (H.3.1)
double Phipnt = dpnt * dgc / Re; // radians
// Calculate quantity s (H.3.2)
double s = sind(Phitn) * Math.cos(Phipnt) + cosd(Phitn) * Math.sin(Phipnt) * cosd(Bt2r);
// The latitude of the intermediate point is now given by (H.3.3)
double Phipntn = asind(s); // degs
// Calculate the quantity x2 (H.3.4a)
double x2 = Math.cos(Phipnt) - s * sind(Phitn);
// Calculate the quantity y2 (H.3.4b)
double y2 = cosd(Phitn) * Math.sin(Phipnt) * sind(Bt2r);
// Calculate the longitude of the intermediate point Phipnte (H.3.5)
double Phipnte = Bt2r;
if (x2 < 1e-9 && y2 < 1e-9) {
Phipnte = Bt2r;
} else {
Phipnte = Phite + atan2d(y2, x2);
}
double[] out = new double[4];
out[0] = Phipnte;
out[1] = Phipntn;
out[2] = Bt2r;
out[3] = dgc;
return out;
}
public double sind(double theta_deg) {
double theta_rad = theta_deg * Math.PI / 180.0;
double y = Math.sin(theta_rad);
return y;
}
public double cosd(double theta_deg) {
double theta_rad = theta_deg * Math.PI / 180.0;
double y = Math.cos(theta_rad);
return y;
}
public double atan2d(double y, double x) {
double res = Math.atan2(y, x) * 180.0 / Math.PI;
return res;
}
public double asind(double y) {
double x = Math.asin(y) * 180.0 / Math.PI;
return x;
}
@Override
public Description description() {
String text = "<html>" + DescriptionTags.startTag.div();
......@@ -2304,7 +2466,6 @@ public class P452ver17PropagationModel
"<br><b><u>Note 3:</u></b> SRTM DEMs do not contain only the bare terrain profile, but may include any objects on the ground (clutter).";
text +=
"<br><b><u>Note 4:</u></b> Some tiles of SRTM DEM may contain voids in regions outside the US. When constructing terrain profile in SEAMCAT heights in the profile having void are replaced with the last valid height from the profile.";
text += "<br><b><u>Note 5:</u></b> Implementation uses Inland zone type for Radio climatic zones settings.";
text += DescriptionTags.endTagNotes.div();
text += DescriptionTags.endTag.div() + "</html>";
return new DescriptionImpl("ITU-R P.452-17", text);
......
core/src/test/java/org/seamcat/model/propagation/P452ver16Test.java
+ 38
- 33
View file @ 16daa957
......@@ -17,7 +17,8 @@ public class P452ver16Test {
// v1 01MAY17 Ivica Stevanovic, OFCOM introduced Test8 with land-only terrain profile
// v2 16MAY17 Ivica Stevanovic, OFCOM introduced Test9 and Test10 with ASTER and SRTM terrain
// profile for CRAF example v3 20JUL18 Ivica Stevanovic, OFCOM introduced modifications related
// with 3-D short distances
// with 3-D short distances v4 08MAR23 Ivica Stevanovic, OFCOM introduced changes in function
// calls related to latitude of the mid point
TestUtil util;
......@@ -34,8 +35,8 @@ public class P452ver16Test {
double press = 1013.25;
double DN = 40;
double temp = 15;
double phi_t = 45;
double phi_r = 45;
double phi_path = 45;
double Gt = 0;
double Gr = 0;
double N0 = 325;
......@@ -105,8 +106,8 @@ public class P452ver16Test {
iZonei[i] = 2;
}
double result = calculator.tl_p452(f, p, rDisti, rHighti, iZonei, htg, hrg, phi_t, phi_r, Gt, Gr, pol, dct,
dcr, DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
double result = calculator.tl_p452(f, p, rDisti, rHighti, iZonei, htg, hrg, phi_path, Gt, Gr, pol, dct, dcr,
DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
util.assertDoubleEquals(expectedResult[ii], result);
}
}
......@@ -118,8 +119,8 @@ public class P452ver16Test {
double press = 1013.25;
double DN = 40;
double temp = 15;
double phi_t = 45;
double phi_r = 45;
double phi_path = 45;
double Gt = 0;
double Gr = 0;
double N0 = 325;
......@@ -189,8 +190,8 @@ public class P452ver16Test {
iZonei[i] = 2;
}
double result = calculator.tl_p452(f, p, rDisti, rHighti, iZonei, htg, hrg, phi_t, phi_r, Gt, Gr, pol, dct,
dcr, DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
double result = calculator.tl_p452(f, p, rDisti, rHighti, iZonei, htg, hrg, phi_path, Gt, Gr, pol, dct, dcr,
DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
util.assertDoubleEquals(expectedResult[ii], result);
}
}
......@@ -202,8 +203,8 @@ public class P452ver16Test {
double press = 1013.25;
double DN = 40;
double temp = 15;
double phi_t = 45;
double phi_r = 45;
double phi_path = 45;
double Gt = 0;
double Gr = 0;
double N0 = 325;
......@@ -273,8 +274,8 @@ public class P452ver16Test {
iZonei[i] = 2;
}
double result = calculator.tl_p452(f, p, rDisti, rHighti, iZonei, htg, hrg, phi_t, phi_r, Gt, Gr, pol, dct,
dcr, DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
double result = calculator.tl_p452(f, p, rDisti, rHighti, iZonei, htg, hrg, phi_path, Gt, Gr, pol, dct, dcr,
DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
util.assertDoubleEquals(expectedResult[ii], result);
}
}
......@@ -286,8 +287,8 @@ public class P452ver16Test {
double press = 1013.25;
double DN = 40;
double temp = 15;
double phi_t = 45;
double phi_r = 45;
double phi_path = 45;
double Gt = 0;
double Gr = 0;
double N0 = 325;
......@@ -363,8 +364,8 @@ public class P452ver16Test {
iZonei[i] = 2;
}
double result = calculator.tl_p452(f, p, rDisti, rHighti, iZonei, htg, hrg, phi_t, phi_r, Gt, Gr, pol, dct,
dcr, DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
double result = calculator.tl_p452(f, p, rDisti, rHighti, iZonei, htg, hrg, phi_path, Gt, Gr, pol, dct, dcr,
DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
util.assertDoubleEquals(expectedResult[ii], result);
}
}
......@@ -376,8 +377,8 @@ public class P452ver16Test {
double press = 1013.25;
double DN = 40;
double temp = 15;
double phi_t = 45;
double phi_r = 45;
double phi_path = 45;
double Gt = 0;
double Gr = 0;
double N0 = 325;
......@@ -447,8 +448,8 @@ public class P452ver16Test {
iZonei[i] = 2;
}
double result = calculator.tl_p452(f, p, rDisti, rHighti, iZonei, htg, hrg, phi_t, phi_r, Gt, Gr, pol, dct,
dcr, DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
double result = calculator.tl_p452(f, p, rDisti, rHighti, iZonei, htg, hrg, phi_path, Gt, Gr, pol, dct, dcr,
DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
util.assertDoubleEquals(expectedResult[ii], result);
}
}
......@@ -460,8 +461,8 @@ public class P452ver16Test {
double press = 1013.25;
double DN = 40;
double temp = 15;
double phi_t = 45;
double phi_r = 45;
double phi_path = 45;
double Gt = 0;
double Gr = 0;
double N0 = 325;
......@@ -531,8 +532,8 @@ public class P452ver16Test {
iZonei[i] = 2;
}
double result = calculator.tl_p452(f, p, rDisti, rHighti, iZonei, htg, hrg, phi_t, phi_r, Gt, Gr, pol, dct,
dcr, DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
double result = calculator.tl_p452(f, p, rDisti, rHighti, iZonei, htg, hrg, phi_path, Gt, Gr, pol, dct, dcr,
DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
util.assertDoubleEquals(expectedResult[ii], result);
}
}
......@@ -606,8 +607,8 @@ public class P452ver16Test {
double press = 1013.000000;
double DN = 50.000000;
double temp = 15.000000;
double phi_t = 50.525000;
double phi_r = 50.525000;
double phi_path = 50.525000;
double Gt = 12.000000;
double Gr = 22.000000;
double N0 = 301.000000;
......@@ -1927,8 +1928,8 @@ public class P452ver16Test {
rHighti[i] = rHeight[i];
iZonei[i] = iZone[i];
}
double result = calculator.tl_p452(f, p, rDisti, rHighti, iZonei, htg, hrg, phi_t, phi_r, Gt, Gr, pol, dct,
dcr, DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
double result = calculator.tl_p452(f, p, rDisti, rHighti, iZonei, htg, hrg, phi_path, Gt, Gr, pol, dct, dcr,
DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
util.assertDoubleEquals(expectedResult[ii - 1], result);
......@@ -1943,6 +1944,7 @@ public class P452ver16Test {
double DN = 50.000000;
double temp = 15.000000;
double phi_t = 50.525000;
double phi_path = 0.0;
double Gt = 12.000000;
double Gr = 22.000000;
......@@ -3764,8 +3766,9 @@ public class P452ver16Test {
iZonei[i] = iZone[i];
}
rphi_r = phi_r[5 + ii - 2];
double result = calculator.tl_p452(f, p, rDisti, rHighti, iZonei, htg, hrg, phi_t, rphi_r, Gt, Gr, pol, dct,
dcr, DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
phi_path = (phi_t + rphi_r) / 2.0;
double result = calculator.tl_p452(f, p, rDisti, rHighti, iZonei, htg, hrg, phi_path, Gt, Gr, pol, dct, dcr,
DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
util.assertDoubleEquals(expectedResult[ii - 1], result);
}
}
......@@ -3777,6 +3780,7 @@ public class P452ver16Test {
double DN = 50.000000;
double temp = 15.000000;
double phi_t = 50.525000;
double phi_path = 0.0;
double Gt = 12.000000;
double Gr = 22.000000;
......@@ -9126,8 +9130,9 @@ public class P452ver16Test {
iZonei[i] = iZone[i];
}
rphi_r = phi_r[5 + ii - 2];
double result = calculator.tl_p452(f, p, rDisti, rHighti, iZonei, htg, hrg, phi_t, rphi_r, Gt, Gr, pol, dct,
dcr, DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
phi_path = (phi_t + rphi_r) / 2.0;
double result = calculator.tl_p452(f, p, rDisti, rHighti, iZonei, htg, hrg, phi_path, Gt, Gr, pol, dct, dcr,
DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
util.assertDoubleEquals(expectedResult[ii - 1], result);
}
}
......
core/src/test/java/org/seamcat/model/propagation/P452ver17Test.java
+ 80
- 58
View file @ 16daa957
......@@ -18,6 +18,8 @@ public class P452ver17Test {
// v2 16MAY17 Ivica Stevanovic, OFCOM introduced Test9 and Test10 with ASTER and SRTM terrain
// profile for CRAF example v3 20JUL18 Ivica Stevanovic, OFCOM introduced modifications related
// with 3-D short distances v4 25OCT21 Ivica Stevanovic, OFCOM Implementation for ITU-R P.452-17
// v5 08MAR23 Ivica Stevanovic, OFCOM Updated to take into account additional arguments lon_t and
// lon_r in function call
TestUtil util;
......@@ -34,8 +36,7 @@ public class P452ver17Test {
double press = 1013.25;
double DN = 40;
double temp = 15;
double phi_t = 45;
double phi_r = 45;
double phi_path = 45;
double Gt = 0;
double Gr = 0;
double N0 = 325;
......@@ -105,8 +106,8 @@ public class P452ver17Test {
iZonei[i] = 2;
}
double result = calculator.tl_p452(f, p, rDisti, rHighti, iZonei, htg, hrg, phi_t, phi_r, Gt, Gr, pol, dct,
dcr, DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
double result = calculator.tl_p452(f, p, rDisti, rHighti, iZonei, htg, hrg, phi_path, Gt, Gr, pol, dct, dcr,
DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
util.assertDoubleEquals(expectedResult[ii], result);
}
}
......@@ -118,8 +119,7 @@ public class P452ver17Test {
double press = 1013.25;
double DN = 40;
double temp = 15;
double phi_t = 45;
double phi_r = 45;
double phi_path = 45;
double Gt = 0;
double Gr = 0;
double N0 = 325;
......@@ -189,8 +189,8 @@ public class P452ver17Test {
iZonei[i] = 2;
}
double result = calculator.tl_p452(f, p, rDisti, rHighti, iZonei, htg, hrg, phi_t, phi_r, Gt, Gr, pol, dct,
dcr, DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
double result = calculator.tl_p452(f, p, rDisti, rHighti, iZonei, htg, hrg, phi_path, Gt, Gr, pol, dct, dcr,
DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
util.assertDoubleEquals(expectedResult[ii], result);
}
}
......@@ -202,8 +202,7 @@ public class P452ver17Test {
double press = 1013.25;
double DN = 40;
double temp = 15;
double phi_t = 45;
double phi_r = 45;
double phi_path = 45;
double Gt = 0;
double Gr = 0;
double N0 = 325;
......@@ -273,8 +272,8 @@ public class P452ver17Test {
iZonei[i] = 2;
}
double result = calculator.tl_p452(f, p, rDisti, rHighti, iZonei, htg, hrg, phi_t, phi_r, Gt, Gr, pol, dct,
dcr, DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
double result = calculator.tl_p452(f, p, rDisti, rHighti, iZonei, htg, hrg, phi_path, Gt, Gr, pol, dct, dcr,
DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
util.assertDoubleEquals(expectedResult[ii], result);
}
}
......@@ -286,8 +285,7 @@ public class P452ver17Test {
double press = 1013.25;
double DN = 40;
double temp = 15;
double phi_t = 45;
double phi_r = 45;
double phi_path = 45;
double Gt = 0;
double Gr = 0;
double N0 = 325;
......@@ -363,8 +361,8 @@ public class P452ver17Test {
iZonei[i] = 2;
}
double result = calculator.tl_p452(f, p, rDisti, rHighti, iZonei, htg, hrg, phi_t, phi_r, Gt, Gr, pol, dct,
dcr, DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
double result = calculator.tl_p452(f, p, rDisti, rHighti, iZonei, htg, hrg, phi_path, Gt, Gr, pol, dct, dcr,
DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
util.assertDoubleEquals(expectedResult[ii], result);
}
}
......@@ -376,8 +374,7 @@ public class P452ver17Test {
double press = 1013.25;
double DN = 40;
double temp = 15;
double phi_t = 45;
double phi_r = 45;
double phi_path = 45;
double Gt = 0;
double Gr = 0;
double N0 = 325;
......@@ -447,8 +444,8 @@ public class P452ver17Test {
iZonei[i] = 2;
}
double result = calculator.tl_p452(f, p, rDisti, rHighti, iZonei, htg, hrg, phi_t, phi_r, Gt, Gr, pol, dct,
dcr, DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
double result = calculator.tl_p452(f, p, rDisti, rHighti, iZonei, htg, hrg, phi_path, Gt, Gr, pol, dct, dcr,
DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
util.assertDoubleEquals(expectedResult[ii], result);
}
}
......@@ -460,8 +457,7 @@ public class P452ver17Test {
double press = 1013.25;
double DN = 40;
double temp = 15;
double phi_t = 45;
double phi_r = 45;
double phi_path = 45;
double Gt = 0;
double Gr = 0;
double N0 = 325;
......@@ -531,8 +527,8 @@ public class P452ver17Test {
iZonei[i] = 2;
}
double result = calculator.tl_p452(f, p, rDisti, rHighti, iZonei, htg, hrg, phi_t, phi_r, Gt, Gr, pol, dct,
dcr, DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
double result = calculator.tl_p452(f, p, rDisti, rHighti, iZonei, htg, hrg, phi_path, Gt, Gr, pol, dct, dcr,
DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
util.assertDoubleEquals(expectedResult[ii], result);
}
}
......@@ -607,8 +603,7 @@ public class P452ver17Test {
double press = 1013.000000;
double DN = 53.000000;
double temp = 15.000000;
double phi_t = 40.250000;
double phi_r = 40.250000;
double phi_path = 40.250000;
double Gt = 20.000000;
double Gr = 5.000000;
double N0 = 328.000000;
......@@ -3707,14 +3702,14 @@ public class P452ver17Test {
expectedResult[34] = 272.839206;
for (int ii = 1; ii <= 18; ii++) {
double pfix = 49;
double result = calculator.tl_p452(ff[ii - 1], pfix, rDist, rHeight, iZone, htg, hrg, phi_t, phi_r, Gt, Gr,
pol, dct, dcr, DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
double result = calculator.tl_p452(ff[ii - 1], pfix, rDist, rHeight, iZone, htg, hrg, phi_path, Gt, Gr, pol,
dct, dcr, DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
util.assertDoubleEquals(expectedResult[ii - 1], result);
}
for (int ii = 1; ii <= 17; ii++) {
double ffix = 2;
double result = calculator.tl_p452(ffix, pp[ii - 1], rDist, rHeight, iZone, htg, hrg, phi_t, phi_r, Gt, Gr,
pol, dct, dcr, DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
double result = calculator.tl_p452(ffix, pp[ii - 1], rDist, rHeight, iZone, htg, hrg, phi_path, Gt, Gr, pol,
dct, dcr, DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
util.assertDoubleEquals(expectedResult[ii - 1 + 18], result);
}
}
......@@ -3742,9 +3737,13 @@ public class P452ver17Test {
double dk_t = 0.000000;
double ha_r = 0.000000;
double dk_r = 0.000000;
double[] rDist = new double[101];
double[] rHeight = new double[101];
int[] iZone = new int[101];
double phi_path = (phi_t + phi_r) / 2.0;
rDist[0] = 0.000000;
rDist[1] = 1.000000;
rDist[2] = 2.000000;
......@@ -4126,14 +4125,14 @@ public class P452ver17Test {
expectedResult[34] = 193.141029;
for (int ii = 1; ii <= 18; ii++) {
double pfix = 49;
double result = calculator.tl_p452(ff[ii - 1], pfix, rDist, rHeight, iZone, htg, hrg, phi_t, phi_r, Gt, Gr,
pol, dct, dcr, DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
double result = calculator.tl_p452(ff[ii - 1], pfix, rDist, rHeight, iZone, htg, hrg, phi_path, Gt, Gr, pol,
dct, dcr, DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
util.assertDoubleEquals(expectedResult[ii - 1], result);
}
for (int ii = 1; ii <= 17; ii++) {
double ffix = 2;
double result = calculator.tl_p452(ffix, pp[ii - 1], rDist, rHeight, iZone, htg, hrg, phi_t, phi_r, Gt, Gr,
pol, dct, dcr, DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
double result = calculator.tl_p452(ffix, pp[ii - 1], rDist, rHeight, iZone, htg, hrg, phi_path, Gt, Gr, pol,
dct, dcr, DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
util.assertDoubleEquals(expectedResult[ii - 1 + 18], result);
}
}
......@@ -4161,9 +4160,13 @@ public class P452ver17Test {
double dk_t = 0.000000;
double ha_r = 0.000000;
double dk_r = 0.000000;
double[] rDist = new double[2002];
double[] rHeight = new double[2002];
int[] iZone = new int[2002];
double phi_path = (phi_t + phi_r) / 2.0;
rDist[0] = 0.000000;
rDist[1] = 0.034953;
rDist[2] = 0.069905;
......@@ -10248,14 +10251,14 @@ public class P452ver17Test {
expectedResult[34] = 194.045255;
for (int ii = 1; ii <= 18; ii++) {
double pfix = 10;
double result = calculator.tl_p452(ff[ii - 1], pfix, rDist, rHeight, iZone, htg, hrg, phi_t, phi_r, Gt, Gr,
pol, dct, dcr, DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
double result = calculator.tl_p452(ff[ii - 1], pfix, rDist, rHeight, iZone, htg, hrg, phi_path, Gt, Gr, pol,
dct, dcr, DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
util.assertDoubleEquals(expectedResult[ii - 1], result);
}
for (int ii = 1; ii <= 17; ii++) {
double ffix = 2;
double result = calculator.tl_p452(ffix, pp[ii - 1], rDist, rHeight, iZone, htg, hrg, phi_t, phi_r, Gt, Gr,
pol, dct, dcr, DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
double result = calculator.tl_p452(ffix, pp[ii - 1], rDist, rHeight, iZone, htg, hrg, phi_path, Gt, Gr, pol,
dct, dcr, DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
util.assertDoubleEquals(expectedResult[ii - 1 + 18], result);
}
}
......@@ -10286,6 +10289,9 @@ public class P452ver17Test {
double[] rDist = new double[110];
double[] rHeight = new double[110];
int[] iZone = new int[110];
double phi_path = (phi_t + phi_r) / 2.0;
rDist[0] = 0.000000;
rDist[1] = 1.000000;
rDist[2] = 2.000000;
......@@ -10694,14 +10700,14 @@ public class P452ver17Test {
expectedResult[34] = 159.795515;
for (int ii = 1; ii <= 18; ii++) {
double pfix = 0.1;
double result = calculator.tl_p452(ff[ii - 1], pfix, rDist, rHeight, iZone, htg, hrg, phi_t, phi_r, Gt, Gr,
pol, dct, dcr, DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
double result = calculator.tl_p452(ff[ii - 1], pfix, rDist, rHeight, iZone, htg, hrg, phi_path, Gt, Gr, pol,
dct, dcr, DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
util.assertDoubleEquals(expectedResult[ii - 1], result);
}
for (int ii = 1; ii <= 17; ii++) {
double ffix = 0.2;
double result = calculator.tl_p452(ffix, pp[ii - 1], rDist, rHeight, iZone, htg, hrg, phi_t, phi_r, Gt, Gr,
pol, dct, dcr, DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
double result = calculator.tl_p452(ffix, pp[ii - 1], rDist, rHeight, iZone, htg, hrg, phi_path, Gt, Gr, pol,
dct, dcr, DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
util.assertDoubleEquals(expectedResult[ii - 1 + 18], result);
}
}
......@@ -10729,9 +10735,13 @@ public class P452ver17Test {
double dk_t = 0.000000;
double ha_r = 0.000000;
double dk_r = 0.000000;
double[] rDist = new double[501];
double[] rHeight = new double[501];
int[] iZone = new int[501];
double phi_path = (phi_t + phi_r) / 2.0;
rDist[0] = 0.000000;
rDist[1] = 0.010000;
rDist[2] = 0.020000;
......@@ -12313,14 +12323,14 @@ public class P452ver17Test {
expectedResult[34] = 112.419795;
for (int ii = 1; ii <= 18; ii++) {
double pfix = 49;
double result = calculator.tl_p452(ff[ii - 1], pfix, rDist, rHeight, iZone, htg, hrg, phi_t, phi_r, Gt, Gr,
pol, dct, dcr, DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
double result = calculator.tl_p452(ff[ii - 1], pfix, rDist, rHeight, iZone, htg, hrg, phi_path, Gt, Gr, pol,
dct, dcr, DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
util.assertDoubleEquals(expectedResult[ii - 1], result);
}
for (int ii = 1; ii <= 17; ii++) {
double ffix = 2;
double result = calculator.tl_p452(ffix, pp[ii - 1], rDist, rHeight, iZone, htg, hrg, phi_t, phi_r, Gt, Gr,
pol, dct, dcr, DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
double result = calculator.tl_p452(ffix, pp[ii - 1], rDist, rHeight, iZone, htg, hrg, phi_path, Gt, Gr, pol,
dct, dcr, DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
util.assertDoubleEquals(expectedResult[ii - 1 + 18], result);
}
}
......@@ -12348,9 +12358,13 @@ public class P452ver17Test {
double dk_t = 0.050000;
double ha_r = 20.000000;
double dk_r = 0.050000;
double[] rDist = new double[501];
double[] rHeight = new double[501];
int[] iZone = new int[501];
double phi_path = (phi_t + phi_r) / 2.0;
rDist[0] = 0.000000;
rDist[1] = 0.010000;
rDist[2] = 0.020000;
......@@ -13932,14 +13946,14 @@ public class P452ver17Test {
expectedResult[34] = 143.469590;
for (int ii = 1; ii <= 18; ii++) {
double pfix = 49;
double result = calculator.tl_p452(ff[ii - 1], pfix, rDist, rHeight, iZone, htg, hrg, phi_t, phi_r, Gt, Gr,
pol, dct, dcr, DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
double result = calculator.tl_p452(ff[ii - 1], pfix, rDist, rHeight, iZone, htg, hrg, phi_path, Gt, Gr, pol,
dct, dcr, DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
util.assertDoubleEquals(expectedResult[ii - 1], result);
}
for (int ii = 1; ii <= 17; ii++) {
double ffix = 2;
double result = calculator.tl_p452(ffix, pp[ii - 1], rDist, rHeight, iZone, htg, hrg, phi_t, phi_r, Gt, Gr,
pol, dct, dcr, DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
double result = calculator.tl_p452(ffix, pp[ii - 1], rDist, rHeight, iZone, htg, hrg, phi_path, Gt, Gr, pol,
dct, dcr, DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
util.assertDoubleEquals(expectedResult[ii - 1 + 18], result);
}
}
......@@ -13967,9 +13981,13 @@ public class P452ver17Test {
double dk_t = 0.020000;
double ha_r = 12.000000;
double dk_r = 0.020000;
double[] rDist = new double[501];
double[] rHeight = new double[501];
int[] iZone = new int[501];
double phi_path = (phi_t + phi_r) / 2.0;
rDist[0] = 0.000000;
rDist[1] = 0.010000;
rDist[2] = 0.020000;
......@@ -15551,14 +15569,14 @@ public class P452ver17Test {
expectedResult[34] = 114.738442;
for (int ii = 1; ii <= 18; ii++) {
double pfix = 49;
double result = calculator.tl_p452(ff[ii - 1], pfix, rDist, rHeight, iZone, htg, hrg, phi_t, phi_r, Gt, Gr,
pol, dct, dcr, DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
double result = calculator.tl_p452(ff[ii - 1], pfix, rDist, rHeight, iZone, htg, hrg, phi_path, Gt, Gr, pol,
dct, dcr, DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
util.assertDoubleEquals(expectedResult[ii - 1], result);
}
for (int ii = 1; ii <= 17; ii++) {
double ffix = 2;
double result = calculator.tl_p452(ffix, pp[ii - 1], rDist, rHeight, iZone, htg, hrg, phi_t, phi_r, Gt, Gr,
pol, dct, dcr, DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
double result = calculator.tl_p452(ffix, pp[ii - 1], rDist, rHeight, iZone, htg, hrg, phi_path, Gt, Gr, pol,
dct, dcr, DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
util.assertDoubleEquals(expectedResult[ii - 1 + 18], result);
}
}
......@@ -15586,9 +15604,13 @@ public class P452ver17Test {
double dk_t = 0.020000;
double ha_r = 25.000000;
double dk_r = 0.020000;
double[] rDist = new double[501];
double[] rHeight = new double[501];
int[] iZone = new int[501];
double phi_path = (phi_t + phi_r) / 2.0;
rDist[0] = 0.000000;
rDist[1] = 0.010000;
rDist[2] = 0.020000;
......@@ -17170,14 +17192,14 @@ public class P452ver17Test {
expectedResult[34] = 149.347206;
for (int ii = 1; ii <= 18; ii++) {
double pfix = 49;
double result = calculator.tl_p452(ff[ii - 1], pfix, rDist, rHeight, iZone, htg, hrg, phi_t, phi_r, Gt, Gr,
pol, dct, dcr, DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
double result = calculator.tl_p452(ff[ii - 1], pfix, rDist, rHeight, iZone, htg, hrg, phi_path, Gt, Gr, pol,
dct, dcr, DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
util.assertDoubleEquals(expectedResult[ii - 1], result);
}
for (int ii = 1; ii <= 17; ii++) {
double ffix = 2;
double result = calculator.tl_p452(ffix, pp[ii - 1], rDist, rHeight, iZone, htg, hrg, phi_t, phi_r, Gt, Gr,
pol, dct, dcr, DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
double result = calculator.tl_p452(ffix, pp[ii - 1], rDist, rHeight, iZone, htg, hrg, phi_path, Gt, Gr, pol,
dct, dcr, DN, N0, press, temp, ha_t, ha_r, dk_t, dk_r, false);
util.assertDoubleEquals(expectedResult[ii - 1 + 18], result);
}
}
......
core/src/test/java/org/seamcat/model/propagation/TerrainHelperTest.java
+ 5
- 0
View file @ 16daa957
......@@ -210,6 +210,11 @@ public class TerrainHelperTest extends DeterministicWorkspaceLoader {
testUtil.assertDoubleEquals(IRSSB_EXPECTED, iRssB);
}
@Test
public void testPMP452() throws Exception {
loadWorkspace("/Test_terrain_RadioMet__146.sws");
}
// Test of results and PL with Terrain ISL Cellular & Generic -> VSL Generic systems
// Uses PMPs ITU-r P.452-16, P.452-17, P.2001-4
@Test
......
core/src/test/resources/integrationtests/Test_terrain_RadioMet__146.sws 0 → 100644
+ 0
- 0
View file @ 16daa957
File added
model/src/main/java/org/seamcat/model/plugin/propagation/ClimateData.java
+ 2
- 116
View file @ 16daa957
package org.seamcat.model.plugin.propagation;
import java.io.BufferedReader;
import java.io.IOException;
import java.io.InputStream;
import java.io.InputStreamReader;
import java.util.ArrayList;
import java.util.List;
public class ClimateData {
public static List<String> deltaNData;
public static List<String> n0Data;
/**
* loads the data from ITU-R F452-16 additional file (deltaN and N0)
*/
public static void loadITUdata() {
// TODO get deltaN and N0 from ITU data and set to PMP
deltaNData = new ArrayList<>();
n0Data = new ArrayList<>();
String sDN50 = "/p2001/DN50.TXT";
String sN0 = "/p2001/N050.TXT";
readSource(sDN50, deltaNData);
readSource(sN0, n0Data);
}
private static void readSource(String name, List<String> target) {
InputStream is = ClimateData.class.getResourceAsStream(name);
try (BufferedReader br = new BufferedReader(new InputStreamReader(is))) {
String line = br.readLine();
while (line != null) {
target.add(line);
line = br.readLine();
}
} catch (IOException e) {
}
}
public double[] getClimateData(TerrainCoordinate point) {
ArrayList<Double> posLonInMap = new ArrayList<>();
ArrayList<Double> posLatInMap = new ArrayList<>();
int lon = (int) point.getLon();
int lat = (int) point.getLat();
for (double d = -90; d <= 90; d += 1.5) {
posLatInMap.add(d);
}
for (double l = 0; l <= 360; l += 1.5) {
posLonInMap.add(l);
}
int index = -1; // takes account of the first ++
while (lat <= posLatInMap.get(++index)) {
}
int indexLat0 = index;
int indexLat1 = index + 1;
if (lat < 0)
indexLat1 = index - 1;
index = -1;
if (lon < 0)
lon = 360 + lon;
while (lon > posLonInMap.get(++index)) {
}
int indexLon0 = index;
int indexLon1 = index + 1;
if (lon < 0)
indexLon1 = index - 1;
double[] bilinear = new double[4];
int offset = 9 * indexLon0; // each entry has together with its leading spaces 9 digits
if (offset + 9 > 2169) {
offset = 0; // TODO was at the end of the line at 359 degree -> take the value at zero degree
}
bilinear[0] = Double.parseDouble(deltaNData.get(indexLat0).substring(offset, offset + 9)); // C / R
bilinear[2] = Double.parseDouble(deltaNData.get(indexLat1).substring(offset, offset + 9)); // C / R + 1
offset = 9 * indexLon1;
if (offset + 9 > 2169) {
offset = 0; // TODO was at the end of the line
}
bilinear[1] = Double.parseDouble(deltaNData.get(indexLat0).substring(offset, offset + 9)); // C + 1 / R
bilinear[3] = Double.parseDouble(deltaNData.get(indexLat1).substring(offset, offset + 9)); // C + 1 / R + 1
double rDeltaN = doBilinear(bilinear, point, lat, lon);
offset = 9 * indexLon0; // each entry has together with its leading spaces 9 digits
if (offset + 9 > 2169) {
offset = 0; // TODO was at the end of the line at 359 degree -> take the value at zero degree
}
bilinear[0] = Double.parseDouble(n0Data.get(indexLat0).substring(offset, offset + 9)); // C / R
bilinear[2] = Double.parseDouble(n0Data.get(indexLat1).substring(offset, offset + 9)); // C / R + 1
offset = 9 * indexLon1;
if (offset + 9 > 2169) {
offset = 0; // TODO was at the end of the line at 359 degree -> take the value at zero degree
}
bilinear[1] = Double.parseDouble(n0Data.get(indexLat0).substring(offset, offset + 9)); // C + 1 / R
bilinear[3] = Double.parseDouble(n0Data.get(indexLat1).substring(offset, offset + 9)); // C + 1 / R + 1
double rN0 = doBilinear(bilinear, point, lat, lon);
return new double[] {rDeltaN, rN0};
}
private double doBilinear(double[] bilinear, TerrainCoordinate point, int lat, int lon) {
int C = lon;
int C1 = C + 1;
int R = lat;
int R1 = R + 1;
double p0 = point.getLat();
double p1 = point.getLon();
return bilinear[0] * (R1 - p0) * (C1 - p1) + bilinear[2] * (p0 - R) * (C1 - p1)
+ bilinear[1] * (R1 - p0) * (p1 - C) + bilinear[3] * (p0 - R) * (p1 - C);
}
public interface ClimateData {
double[] getClimateData(TerrainCoordinate point);
}
model/src/main/java/org/seamcat/model/plugin/propagation/P452ver16Input.java
+ 4
- 4
View file @ 16daa957
......@@ -18,7 +18,7 @@ public interface P452ver16Input {
boolean variations();
boolean variations = false; // REMARK: added by KK
@Config(order = 2, name = "Use Terrain Profile Data")
@Config(order = 2, name = "Use Terrain Profile Data", defineGroup = "useTerrain")
boolean terrain();
boolean terrain = false;
......@@ -31,16 +31,16 @@ public interface P452ver16Input {
@Config(order = 5, name = "Surface pressure", unit = Unit.hPa)
double surfacePressure();
@Config(order = 6, name = "Refraction index gradient", unit = Unit.perKm)
@Config(order = 6, name = "Refractivity gradient", unit = Unit.perKm, invertedGroup = "useTerrain")
double refractionIndex();
@Config(order = 7, name = "Surface temperature", unit = Unit.degC)
double surfaceTemperature();
@Config(order = 8, name = "Latitude", unit = Unit.deg)
@Config(order = 8, name = "Latitude", unit = Unit.deg, invertedGroup = "useTerrain")
double latitude();
@Config(order = 9, name = "Sea level surface refractivity")
@Config(order = 9, name = "Sea level surface refractivity", invertedGroup = "useTerrain")
double seaLevelSurfaceRefractivity();
@Config(order = 10, name = "Time percentage", unit = Unit.percent,
......
model/src/main/java/org/seamcat/model/plugin/propagation/P452ver17Input.java
+ 4
- 4
View file @ 16daa957
......@@ -18,7 +18,7 @@ public interface P452ver17Input {
boolean variations();
boolean variations = false; // REMARK: added by KK
@Config(order = 2, name = "Use Terrain Profile Data")
@Config(order = 2, name = "Use Terrain Profile Data", defineGroup = "useTerrain")
boolean terrain();
boolean terrain = false;
......@@ -31,16 +31,16 @@ public interface P452ver17Input {
@Config(order = 5, name = "Surface pressure", unit = Unit.hPa)
double surfacePressure();
@Config(order = 6, name = "Refraction index gradient", unit = Unit.perKm)
@Config(order = 6, name = "Refractivity gradient", unit = Unit.perKm, invertedGroup = "useTerrain")
double refractionIndex();
@Config(order = 7, name = "Surface temperature", unit = Unit.degC)
double surfaceTemperature();
@Config(order = 8, name = "Latitude", unit = Unit.deg)
@Config(order = 8, name = "Latitude", unit = Unit.deg, invertedGroup = "useTerrain")
double latitude();
@Config(order = 9, name = "Sea level surface refractivity")
@Config(order = 9, name = "Sea level surface refractivity", invertedGroup = "useTerrain")
double seaLevelSurfaceRefractivity();
@Config(order = 10, name = "Time percentage", unit = Unit.percent,
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment