+// Calculate shortest distance from point to line
+static double point_line_dist_squared( const Point3D& tc, const Point3D& vp,
+ MAT3vec d )
+{
+ MAT3vec p, p0;
+ double dist;
+
+ p[0] = tc.x(); p[1] = tc.y(); p[2] = tc.z();
+ p0[0] = vp.x(); p0[1] = vp.y(); p0[2] = vp.z();
+
+ dist = fgPointLineSquared(p, p0, d);
+
+ // cout << "dist = " << dist << endl;
+
+ return(dist);
+}
+
+
+// Determine scenery altitude. Normally this just happens when we
+// render the scene, but we'd also like to be able to do this
+// explicitely. lat & lon are in radians. abs_view_pos in meters.
+// Returns result in meters.
+double
+fgTileMgrCurElev( const fgBUCKET& p ) {
+ fgTILE *t;
+ fgFRAGMENT *frag_ptr;
+ Point3D abs_view_pos = current_view.abs_view_pos;
+ Point3D earth_center(0.0);
+ Point3D result;
+ MAT3vec local_up;
+ double dist, lat_geod, alt, sea_level_r;
+ int index;
+
+ local_up[0] = abs_view_pos.x();
+ local_up[1] = abs_view_pos.y();
+ local_up[2] = abs_view_pos.z();
+
+ // Find current translation offset
+ // fgBucketFind(lon * RAD_TO_DEG, lat * RAD_TO_DEG, &p);
+ index = global_tile_cache.exists(p);
+ if ( index < 0 ) {
+ FG_LOG( FG_TERRAIN, FG_WARN, "Tile not found" );
+ return 0.0;
+ }
+
+ t = global_tile_cache.get_tile(index);
+
+ scenery.next_center = t->center;
+
+ // earth_center = Point3D(0.0, 0.0, 0.0);
+
+ FG_LOG( FG_TERRAIN, FG_DEBUG,
+ "Current bucket = " << p << " Index = " << fgBucketGenIndex(&p) );
+
+ // calculate tile offset
+ // x = (t->offset.x = t->center.x - scenery.center.x);
+ // y = (t->offset.y = t->center.y - scenery.center.y);
+ // z = (t->offset.z = t->center.z - scenery.center.z);
+
+ // calc current terrain elevation calculate distance from
+ // vertical tangent line at current position to center of
+ // tile.
+
+ /* printf("distance squared = %.2f, bounding radius = %.2f\n",
+ point_line_dist_squared(&(t->offset), &(v->view_pos),
+ v->local_up), t->bounding_radius); */
+
+ dist = point_line_dist_squared( t->center, abs_view_pos, local_up );
+ if ( dist < FG_SQUARE(t->bounding_radius) ) {
+
+ // traverse fragment list for tile
+ fgTILE::FragmentIterator current = t->begin();
+ fgTILE::FragmentIterator last = t->end();
+
+ for ( ; current != last; ++current ) {
+ frag_ptr = &(*current);
+ /* printf("distance squared = %.2f, bounding radius = %.2f\n",
+ point_line_dist_squared( &(frag_ptr->center),
+ &abs_view_pos), local_up),
+ frag_ptr->bounding_radius); */
+
+ dist = point_line_dist_squared( frag_ptr->center,
+ abs_view_pos,
+ local_up);
+ if ( dist <= FG_SQUARE(frag_ptr->bounding_radius) ) {
+ if ( frag_ptr->intersect( abs_view_pos,
+ earth_center, 0, result ) ) {
+ FG_LOG( FG_TERRAIN, FG_DEBUG, "intersection point " <<
+ result );
+ // compute geocentric coordinates of tile center
+ Point3D pp = fgCartToPolar3d(result);
+ FG_LOG( FG_TERRAIN, FG_DEBUG, " polar form = " << pp );
+ // convert to geodetic coordinates
+ fgGeocToGeod(pp.lat(), pp.radius(), &lat_geod,
+ &alt, &sea_level_r);
+
+ // printf("alt = %.2f\n", alt);
+ // exit since we found an intersection
+ if ( alt > -9999.0 ) {
+ // printf("returning alt\n");
+ return alt;
+ } else {
+ // printf("returning 0\n");
+ return 0.0;
+ }
+ }
+ }
+ }
+ }
+
+ cout << "no terrain intersection found\n";
+ return 0.0;
+}
+
+
+// Determine scenery altitude. Normally this just happens when we
+// render the scene, but we'd also like to be able to do this
+// explicitely. lat & lon are in radians. abs_view_pos in meters.
+// Returns result in meters.
+double
+fgTileMgrCurElevOLD( double lon, double lat, const Point3D& abs_view_pos ) {
+ fgTILECACHE *c;
+ fgTILE *t;
+ // fgVIEW *v;
+ fgFRAGMENT *frag_ptr;
+ fgBUCKET p;
+ Point3D earth_center(0.0);
+ Point3D result;
+ MAT3vec local_up;
+ double dist, lat_geod, alt, sea_level_r;
+ // double x, y, z;
+ int index;
+
+ c = &global_tile_cache;
+ // v = ¤t_view;
+
+ local_up[0] = abs_view_pos.x();
+ local_up[1] = abs_view_pos.y();
+ local_up[2] = abs_view_pos.z();
+
+ // Find current translation offset
+ fgBucketFind(lon * RAD_TO_DEG, lat * RAD_TO_DEG, &p);
+ index = c->exists(p);
+ if ( index < 0 ) {
+ FG_LOG( FG_TERRAIN, FG_WARN, "Tile not found" );
+ return 0.0;
+ }
+
+ t = c->get_tile(index);
+
+ scenery.next_center = t->center;
+
+ // earth_center = Point3D(0.0, 0.0, 0.0);
+
+ FG_LOG( FG_TERRAIN, FG_DEBUG,
+ "Pos = (" << lon * RAD_TO_DEG << ", " << lat * RAD_TO_DEG
+ << ") Current bucket = " << p
+ << " Index = " << fgBucketGenIndex(&p) );
+
+ // calculate tile offset
+ // x = (t->offset.x = t->center.x - scenery.center.x);
+ // y = (t->offset.y = t->center.y - scenery.center.y);
+ // z = (t->offset.z = t->center.z - scenery.center.z);
+
+ // calc current terrain elevation calculate distance from
+ // vertical tangent line at current position to center of
+ // tile.
+
+ /* printf("distance squared = %.2f, bounding radius = %.2f\n",
+ point_line_dist_squared(&(t->offset), &(v->view_pos),
+ v->local_up), t->bounding_radius); */
+
+ dist = point_line_dist_squared( t->center, abs_view_pos, local_up );
+ if ( dist < FG_SQUARE(t->bounding_radius) ) {
+
+ // traverse fragment list for tile
+ fgTILE::FragmentIterator current = t->begin();
+ fgTILE::FragmentIterator last = t->end();
+
+ for ( ; current != last; ++current ) {
+ frag_ptr = &(*current);
+ /* printf("distance squared = %.2f, bounding radius = %.2f\n",
+ point_line_dist_squared( &(frag_ptr->center),
+ &abs_view_pos), local_up),
+ frag_ptr->bounding_radius); */
+
+ dist = point_line_dist_squared( frag_ptr->center,
+ abs_view_pos,
+ local_up);
+ if ( dist <= FG_SQUARE(frag_ptr->bounding_radius) ) {
+ if ( frag_ptr->intersect( abs_view_pos,
+ earth_center, 0, result ) ) {
+ FG_LOG( FG_TERRAIN, FG_DEBUG, "intersection point " <<
+ result );
+ // compute geocentric coordinates of tile center
+ Point3D pp = fgCartToPolar3d(result);
+ FG_LOG( FG_TERRAIN, FG_DEBUG, " polar form = " << pp );
+ // convert to geodetic coordinates
+ fgGeocToGeod(pp.lat(), pp.radius(), &lat_geod,
+ &alt, &sea_level_r);
+
+ // printf("alt = %.2f\n", alt);
+ // exit since we found an intersection
+ if ( alt > -9999.0 ) {
+ // printf("returning alt\n");
+ return alt;
+ } else {
+ // printf("returning 0\n");
+ return 0.0;
+ }
+ }
+ }
+ }
+ }
+
+ cout << "no terrain intersection found\n";
+ return 0.0;
+}
+
+