reworked coord transformations to reduce roundoff errors

readhorizon.c

@@ -48,10 +48,7 @@ static int doNoMercator = 0;
 static GLint uniform_view_z;
 static GLint uniform_demfileN, uniform_demfileW;
 static GLint uniform_WDEM;
-static GLint uniform_sin_view_lon;
-static GLint uniform_cos_view_lon;
-static GLint uniform_sin_view_lat;
-static GLint uniform_cos_view_lat;
+static GLint uniform_view_lon, uniform_view_lat;
 static GLint uniform_aspect;
 
 #define WDEM        1201
@@ -314,20 +311,17 @@ static void loadGeometry(void)
     uniform_demfileN     = glGetUniformLocation(program, "demfileN"    ); assert( glGetError() == GL_NO_ERROR );
     uniform_demfileW     = glGetUniformLocation(program, "demfileW"    ); assert( glGetError() == GL_NO_ERROR );
     uniform_WDEM         = glGetUniformLocation(program, "WDEM"        ); assert( glGetError() == GL_NO_ERROR );
-    uniform_sin_view_lon = glGetUniformLocation(program, "sin_view_lon"); assert( glGetError() == GL_NO_ERROR );
-    uniform_cos_view_lon = glGetUniformLocation(program, "cos_view_lon"); assert( glGetError() == GL_NO_ERROR );
-    uniform_sin_view_lat = glGetUniformLocation(program, "sin_view_lat"); assert( glGetError() == GL_NO_ERROR );
-    uniform_cos_view_lat = glGetUniformLocation(program, "cos_view_lat"); assert( glGetError() == GL_NO_ERROR );
+    uniform_view_lat     = glGetUniformLocation(program, "view_lat"    ); assert( glGetError() == GL_NO_ERROR );
+    uniform_view_lon     = glGetUniformLocation(program, "view_lon"    ); assert( glGetError() == GL_NO_ERROR );
     uniform_aspect       = glGetUniformLocation(program, "aspect"      ); assert( glGetError() == GL_NO_ERROR );
 
     glUniform1f( uniform_view_z,       getHeight(view_i, view_j));
     glUniform1i( uniform_demfileN,     demfileN);
     glUniform1i( uniform_demfileW,     demfileW);
     glUniform1i( uniform_WDEM,         WDEM);
-    glUniform1f( uniform_sin_view_lon, sinf( M_PI / 180.0f * view_lon ));
-    glUniform1f( uniform_cos_view_lon, cosf( M_PI / 180.0f * view_lon ));
-    glUniform1f( uniform_sin_view_lat, sinf( M_PI / 180.0f * view_lat ));
-    glUniform1f( uniform_cos_view_lat, cosf( M_PI / 180.0f * view_lat ));
+
+    glUniform1f( uniform_view_lon,     view_lon * M_PI / 180.0f );
+    glUniform1f( uniform_view_lat,     view_lat * M_PI / 180.0f );
   }
 }
 

vertex.glsl

@@ -22,14 +22,19 @@ void main(void)
 
   float lon = radians( float(-demfileW)     + vin.x/float(WDEM-1) );
   float lat = radians( float( demfileN + 1) - vin.y/float(WDEM-1) );
-  float sin_lon  = sin( lon );
-  float cos_lon  = cos( lon );
+
+  // Here I compute 4 sin/cos. Previously I was sending sincos( view_lat/lon) as
+  // a uniform, so no trig was needed here. I think this may have been causing
+  // roundoff issues, so I'm not doing that anymore. Specifically, sin(+eps) was
+  // being slightly negative
+  float sin_dlat = sin( lat - view_lat );
+  float cos_dlat = cos( lat - view_lat );
+  float sin_dlon = sin( lon - view_lon );
+  float cos_dlon = cos( lon - view_lon );
+
   float sin_lat  = sin( lat );
   float cos_lat  = cos( lat );
-  float sin_dlat = sin_lat * cos_view_lat - cos_lat * sin_view_lat;
-  float cos_dlat = cos_lat * cos_view_lat + sin_lat * sin_view_lat;
-  float sin_dlon = sin_lon * cos_view_lon - cos_lon * sin_view_lon;
-  float cos_dlon = cos_lon * cos_view_lon + sin_lon * sin_view_lon;
+
 
   vec3 v = vec3( (Rearth + vin.z) * ( cos_lat * sin_dlon ),
                  (Rearth + vin.z) * ( cos_dlat*cos_dlon + sin_lat*sin_view_lat*(1.0 - cos_dlon) ),

vertex_header.glsl

@@ -4,10 +4,8 @@
 uniform float view_z;
 uniform int   demfileN, demfileW;
 uniform int   WDEM;
-uniform float sin_view_lon;
-uniform float cos_view_lon;
-uniform float sin_view_lat;
-uniform float cos_view_lat;
+uniform float view_lon;
+uniform float view_lat;
 
 uniform float aspect;
 varying float red;