From fad65d8b65673c66554bbf4f85fe46bef7e39542 Mon Sep 17 00:00:00 2001
From: Johannes Hofmann <Johannes.Hofmann@gmx.de>
Date: Tue, 8 Apr 2008 16:16:12 +0200
Subject: add basic refraction correction
---
src/Panorama.H | 3 ++-
src/Panorama.cxx | 19 ++++++++++++++++---
2 files changed, 18 insertions(+), 4 deletions(-)
(limited to 'src')
diff --git a/src/Panorama.H b/src/Panorama.H
index aea48b8..36003d9 100644
--- a/src/Panorama.H
+++ b/src/Panorama.H
@@ -35,6 +35,7 @@ class Panorama {
double distance(double phi, double lam);
double alpha(const Hill *m);
double nick(const Hill *m);
+ double refraction_change(const Hill *m);
double comp_center_angle(double alph_a, double alph_b, double d1, double d2);
double comp_scale(double alph_a, double alph_b, double d1, double d2);
int get_matrix(double m[]);
@@ -73,7 +74,7 @@ class Panorama {
double get_view_long();
double get_view_height();
double get_earth_radius(double latitude);
- double get_real_distance(Hill *m);
+ double get_real_distance(const Hill *m);
int comp_params(Hills *h);
ProjectionLSQ::Projection_t get_projection();
void set_projection(ProjectionLSQ::Projection_t p);
diff --git a/src/Panorama.cxx b/src/Panorama.cxx
index 32e14f4..49010f5 100644
--- a/src/Panorama.cxx
+++ b/src/Panorama.cxx
@@ -428,15 +428,28 @@ Panorama::alpha(const Hill *m) {
return fmod(atan2(sin_alph, cos_alph) + 2.0 * pi_d, 2.0 * pi_d);
}
+double
+Panorama::refraction_change(const Hill *m) {
+ double a, b, c, alpha = 6.5 / 1000.0, T0 = 10.0;
+
+ a = 2.9e-4 * exp (-view_height / 10000.0) / (1.0 + 2.9 * T0 / 760.0);
+ b = 2.9 * alpha / (760.0 * (1.0 + 2.9 * T0 / 760.0));
+ c = a * (b - 1.0 / 10000.0);
+
+ return (c * get_real_distance(m) / (2.0 * (1.0 + a)));
+}
double
Panorama::nick(const Hill *m) {
- double b, c;
+ double b, c, theta = refraction_change(m);
b = m->height + get_earth_radius(m->phi);
c = view_height + get_earth_radius(view_phi);
- return atan((cos(m->dist) * b - c) / (sin(m->dist) * b));
+if (get_real_distance(m) < 200000)
+fprintf(stderr, "=== %s, %g\n", m->name, theta);
+
+ return atan((cos(m->dist) * b - c) / (sin(m->dist) * b)) - theta;
}
// return local distance to center of WGS84 ellipsoid
@@ -452,7 +465,7 @@ Panorama::get_earth_radius(double phi) {
}
double
-Panorama::get_real_distance(Hill *m) {
+Panorama::get_real_distance(const Hill *m) {
double a, b, c, gam;
a = view_height + get_earth_radius(view_phi);
--
cgit v1.2.3