From 0eba95c483f2db216116ed8ebc99848415d04b3f Mon Sep 17 00:00:00 2001
From: sesteves <sroesteves@gmail.com>
Date: Sun, 10 Jul 2016 16:40:03 +0100
Subject: [PATCH] multiple improvements
---
.../AtmosphericRefraction.java | 3 +-
.../MultiLayerModel.java | 90 +++++++++++++------
2 files changed, 67 insertions(+), 26 deletions(-)
diff --git a/src/main/java/org/orekit/rugged/atmosphericrefraction/AtmosphericRefraction.java b/src/main/java/org/orekit/rugged/atmosphericrefraction/AtmosphericRefraction.java
index 4e5e0ba0..837f259a 100644
--- a/src/main/java/org/orekit/rugged/atmosphericrefraction/AtmosphericRefraction.java
+++ b/src/main/java/org/orekit/rugged/atmosphericrefraction/AtmosphericRefraction.java
@@ -19,6 +19,7 @@ package org.orekit.rugged.atmosphericrefraction;
import org.apache.commons.math3.geometry.euclidean.threed.Vector3D;
import org.orekit.rugged.raster.Tile;
+import org.orekit.rugged.utils.NormalizedGeodeticPoint;
/**
* Interface for atmospheric refraction.
@@ -26,6 +27,6 @@ import org.orekit.rugged.raster.Tile;
*/
public interface AtmosphericRefraction {
- double getDeviation(Vector3D pos, Vector3D los, Vector3D zenith, double altitude, Tile tile);
+ NormalizedGeodeticPoint getPointOnGround(Vector3D pos, Vector3D los, Vector3D zenith, double altitude, Tile tile);
}
diff --git a/src/main/java/org/orekit/rugged/atmosphericrefraction/MultiLayerModel.java b/src/main/java/org/orekit/rugged/atmosphericrefraction/MultiLayerModel.java
index 1ebaa007..a7d7a1bc 100644
--- a/src/main/java/org/orekit/rugged/atmosphericrefraction/MultiLayerModel.java
+++ b/src/main/java/org/orekit/rugged/atmosphericrefraction/MultiLayerModel.java
@@ -18,13 +18,17 @@ package org.orekit.rugged.atmosphericrefraction;
import org.apache.commons.math3.geometry.euclidean.threed.Vector3D;
import org.apache.commons.math3.util.FastMath;
+import org.orekit.bodies.OneAxisEllipsoid;
+import org.orekit.errors.OrekitException;
+import org.orekit.frames.FramesFactory;
+import org.orekit.rugged.errors.RuggedException;
import org.orekit.rugged.raster.Tile;
import org.orekit.rugged.utils.ExtendedEllipsoid;
import org.orekit.rugged.utils.NormalizedGeodeticPoint;
+import org.orekit.utils.Constants;
+import org.orekit.utils.IERSConventions;
-import java.util.Collections;
-import java.util.Map;
-import java.util.TreeMap;
+import java.util.*;
/**
* Multi layer model for atmospheric refraction.
@@ -33,7 +37,9 @@ import java.util.TreeMap;
public class MultiLayerModel implements AtmosphericRefraction {
// maps altitude (lower bound) to refraction index
- private static final Map<Double, Double> meanAtmoshpericRefractions;
+ private static final Map<Double, Double> meanAtmosphericRefractions;
+
+ private static final Map<Double, ExtendedEllipsoid> atmosphericEllipsoids;
static {
Map<Double, Double> refractions = new TreeMap(Collections.reverseOrder());
@@ -52,35 +58,69 @@ public class MultiLayerModel implements AtmosphericRefraction {
refractions.put(40000.00000000000, 1.00000100000);
refractions.put(50000.00000000000, 1.00000000000);
refractions.put(100000.00000000000, 1.00000000000);
- meanAtmoshpericRefractions = Collections.unmodifiableMap(refractions);
+ meanAtmosphericRefractions = Collections.unmodifiableMap(refractions);
+
+ Map<Double, ExtendedEllipsoid> ellipsoids = new HashMap<Double, ExtendedEllipsoid>();
+ try {
+ for (Double altitude : refractions.keySet()) {
+ OneAxisEllipsoid ellipsoid = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS + altitude,
+ Constants.WGS84_EARTH_FLATTENING, FramesFactory.getITRF(IERSConventions.IERS_2010, true));
+ ellipsoid = new ExtendedEllipsoid(ellipsoid.getEquatorialRadius(), ellipsoid.getFlattening(),
+ ellipsoid.getBodyFrame());
+ ellipsoids.put(altitude, (ExtendedEllipsoid) ellipsoid);
+ }
+ } catch (OrekitException e) {
+ e.printStackTrace();
+ }
+ atmosphericEllipsoids = Collections.unmodifiableMap(ellipsoids);
}
@Override
- public double getDeviation(Vector3D pos, Vector3D los, Vector3D zenith, double altitude, Tile tile) {
+ public NormalizedGeodeticPoint getPointOnGround(Vector3D initialPos, Vector3D initialLos, Vector3D initialZenith,
+ double altitude, Tile tile) {
+ NormalizedGeodeticPoint newGeodeticPoint = null;
+ try {
- new ExtendedEllipsoid(ellipsoid.getEquatorialRadius(), ellipsoid.getFlattening(),
- ellipsoid.getBodyFrame());
+ Vector3D pos = initialPos;
+ Vector3D los = initialLos;
+ Vector3D zenith = initialZenith;
+ double theta1 = Vector3D.angle(los, zenith), theta2;
+ double previousRefractionIndex = -1;
+ NormalizedGeodeticPoint gp = null;
+ for (Map.Entry<Double, Double> entry : meanAtmosphericRefractions.entrySet()) {
+ if (pos.getZ() < entry.getKey()) {
+ continue;
+ }
- double incidenceAngleSin = FastMath.sin(Vector3D.angle(los, zenith));
- double previousRefractionIndex = -1;
- double xDistance = 0;
- for(Map.Entry<Double, Double> entry : meanAtmoshpericRefractions.entrySet()) {
- if(pos.getZ() < entry.getKey()) {
- continue;
- }
- if(previousRefractionIndex > 0) {
- incidenceAngleSin = previousRefractionIndex * incidenceAngleSin / entry.getValue();
- }
- xDistance += (pos.getZ() - entry.getKey()) * FastMath.tan(incidenceAngleSin);
+ if (previousRefractionIndex > 0) {
+ theta2 = previousRefractionIndex * FastMath.sin(theta1) / entry.getValue();
- if(altitude > entry.getKey()) {
- break;
+ // get current los
+ double a = FastMath.sqrt((1 - FastMath.pow(FastMath.cos(theta2), 2)) /
+ (1 - FastMath.pow(FastMath.cos(theta1), 2)));
+ double b = a * FastMath.cos(theta1) - FastMath.cos(theta2);
+ los = new Vector3D(a, los, b, zenith);
+
+ theta1 = theta2;
+ }
+
+ // get intersection point
+ ExtendedEllipsoid ellipsoid = atmosphericEllipsoids.get(entry.getKey());
+ gp = ellipsoid.pointOnGround(pos, los, 0.0);
+ pos = ellipsoid.transform(gp);
+ zenith = gp.getZenith();
+
+ if (altitude > entry.getKey()) {
+ break;
+ }
+ previousRefractionIndex = entry.getValue();
}
- previousRefractionIndex = entry.getValue();
- }
- NormalizedGeodeticPoint geodeticPoint = tile.cellIntersection(pos, los, 0, 0);
+ newGeodeticPoint = tile.cellIntersection(gp, los, 0, 0);
+ } catch (RuggedException e) {
+ e.printStackTrace();
+ }
- return pos.getX() + xDistance;
+ return newGeodeticPoint;
}
}
--
GitLab