Scenery.cpp

Go to the documentation of this file.

/*
 * Scenery.cpp
 *
 * Part of Fly! Legacy project
 *
 * Copyright 2003 Chris Wallace
 *
 * Fly! Legacy is free software; you can redistribute it and/or modify
 *   it under the terms of the GNU General Public License as published by
 *   the Free Software Foundation; either version 2 of the License, or
 *   (at your option) any later version.
 *
 * Fly! Legacy is distributed in the hope that it will be useful,
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *   GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 *   along with Fly! Legacy; if not, write to the Free Software
 *    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 */

#include <math.h>
#include "../Include/Utility.h"
#include "../Include/Globals.h"


//
// If this macro is defined, then a table of globe tile latitude boundaries is
//   dumped to a text file for debugging immediately after initialization
//
//#define DUMP_GLOBE_TILE_BOUNDARIES


//
// Allocate globe tile latitude lookup table
//
static double globe_tile_lat_table [129];

void init_globe_tile_table ()
{
  int i;

  double k = 360.0 / 256.0;
  globe_tile_lat_table[0] = 0.0;
  for (i=1; i<129; i++) {
    double prev = globe_tile_lat_table[i-1];
    double prevRad = DegToRad (prev);
    globe_tile_lat_table[i] = prev + k * cos(prevRad);
  }

#ifdef DUMP_GLOBE_TILE_BOUNDARIES
  // DEBUG : Dump table of globe tile latitude boundaries
  FILE* f = fopen ("Debug\\LatitudeTable.txt", "w");
  fprintf (f, "GT Table:\n\n");
  fprintf (f, "Z-Index     South Boundary      North Boundary\n");
  for (i=0; i<256; i++) {
    double n, s;
    globe_tile_lat_bounds (i, s, n);

    SPosition pos;
    pos.lon = pos.alt = 0;

    pos.lat = s * 3600;
    char temp[64];
    FormatPosition (pos, temp);
    char south[64];
    memset (south, 0, 64);
    strncpy (south, temp, 16);

    pos.lat = n * 3600;
    FormatPosition (pos, temp);
    char north[64];
    memset (north, 0, 64);
    strncpy (north, temp, 16);

    double refLat;
    if (i >= 128) {
      // Northern hemisphere, use southern boundary latitude as reference
      refLat = s * 3600;
    } else {
      // Southern hemisphere, use northern boundary latitude as reference
      refLat = n * 3600;
    }
    fprintf (f, " %3d       %s    %s\n", i, south, north);
  }

  // Dump QGT table
  fprintf (f, "\n\n");
  fprintf (f, "QGT Table:\n");
  fprintf (f, "----------\n");
  fprintf (f, "Z-Index     South Boundary      North Boundary\n");
  for (i=0; i<512; i++) {
    double n, s;
    qgt_lat_bounds (i, s, n);

    SPosition pos;
    pos.lon = pos.alt = 0;

    pos.lat = s * 3600;
    char temp[64];
    FormatPosition (pos, temp);
    char south[64];
    memset (south, 0, 64);
    strncpy (south, temp, 16);

    pos.lat = n * 3600;
    FormatPosition (pos, temp);
    char north[64];
    memset (north, 0, 64);
    strncpy (north, temp, 16);

    double refLat;
    if (i >= 128) {
      // Northern hemisphere, use southern boundary latitude as reference
      refLat = s * 3600;
    } else {
      // Southern hemisphere, use northern boundary latitude as reference
      refLat = n * 3600;
    }
    fprintf (f, " %3d       %s    %s\n", i, south, north);
  }

  fclose (f);
#endif
}


static double globe_tile_lon (int x)
{
  return (double)x * 360.0 / 256.0;
}

void qgt_position (SPosition pos, int &qx, int &qz, double &x, double &z)
{
  // Get boundaries of the position QGT and the reference QGT
  lat_lon_to_qgt (pos.lat, pos.lon, qx, qz);
  double w, e, s, n;
  qgt_lon_bounds (qx, w, e);
  qgt_lat_bounds (qz, s, n);
  w *= 3600.0;
  e *= 3600.0;
  n *= 3600.0;
  s *= 3600.0;

  double xSpan = (e - w) / 2.0;
  double zSpan = (n - s) / 2.0;
  double xMid = w + xSpan;
  double zMid = s + zSpan;
  x = -(pos.lon - xMid) / xSpan;
  z = -(pos.lat - zMid) / zSpan;
}

void delta_qgt (SPosition pos, int rx, int rz, double &x, double &z)
{
  // x,z indices of the position QGT
  int px, pz;
  qgt_position (pos, px, pz, x, z);
  // Adjust deltas by differences in indices
  x += (double)(rx - px);
  z += (double)(rz - pz);
}

void lat_lon_to_qgt (double latArg, double lonArg, int &x, int &z)
{
  // Convert from arcsec to degrees
  latArg /= 3600.0;
  lonArg /= 3600.0;

  int gtx = 0, gtz = 0;
  lat_lon_to_globe_tile (latArg, lonArg, gtx, gtz);

  // Determine longitude offset, 0 if the position is in the western half of the GT,
  //   1 if in the eastern half
  int xOffset = 0;
  double w, e;
  globe_tile_lon_bounds (gtx, w, e);

  if (lonArg > (w + ((e-w)/2.0))) {
    // Position is in the eastern half of the globe tile
    xOffset = 1;
  }

  // Determine latitude offset:
  //   0 if position is in the southern half of the GT
  //   1 if position is in the northern half of the GT
  int zOffset = 0;
  double n, s;
  globe_tile_lat_bounds (gtz, s, n);
  if (latArg > (s + ((n-s)/2))) {
    // Position is in the northern half of the globe tile
    zOffset = 1;
  }

  // Assign QGT indices
  x = (gtx * 2) + xOffset;
  z = (gtz * 2) + zOffset;
}

void lat_lon_to_globe_tile (double latArg, double lonArg, int &x, int &z)
{
  /*
   * Calculate longitude global and detail tile indices
   *
   * This is the easy calculation.  The first globe tile (index 0) has
   *   the Prime Meridian (longitude 0.0) as its western boundary. Each
   *   tile is equally spaced and numbered progressively eastward
   *   around the globe.
   *
   * The longitude argument is first normalized so that it represents
   *   the number of degrees East of the Prime Meridian, then the longitude
   *   globe tile and detail tile are simply calculated by dividing by 256
   */

  // Range-check latitude and longitude parameters
  double absLat = fabs(latArg);
  if ((lonArg < 0.0) || (lonArg > 360.0) || (absLat > globe_tile_lat_table[128])) {
    // Fatal error
    gtfo ("lat_lon_to_globe_tile : Illegal lat/lon lat=%f lon=%f", latArg, lonArg);
  }

  // Calculate longitude index
  x = (int)((lonArg / 360.0) * 256.0);

  // Binary search through latitude table for latitude index
  int min = 0;
  int max = 128;
  while ((max - min) > 1) {
    int mid = min + ((max - min) / 2);
    if (absLat > globe_tile_lat_table[mid]) {
      min = mid;
    } else {
      max = mid;
    }
  }

  if (latArg > 0.0) {
    // Northern hemisphere
    z = min + 128;
  } else {
    // Southern hemisphere
    z = 127 - min;
  }
}

void globe_tile_lon_bounds (int x, double &w, double &e)
{
  if ((x < 0) || (x > 255)) {
    // Fatal error
    gtfo ("globe_tile_lon_bounds : Illegal GT index x=%d", x);
  } else {
    w = globe_tile_lon (x);
    e = globe_tile_lon (x+1);
  }
}


void globe_tile_lat_bounds (int z, double &s, double &n)
{
  if (z <= 127) {
    // Southern hemisphere
    n = -globe_tile_lat_table[127 - z];
    s = -globe_tile_lat_table[128 - z];
  } else if (z <= 255) {
    // Northern hemisphere;
    n = globe_tile_lat_table[z - 127];
    s = globe_tile_lat_table[z - 128];
  } else {
    // Fatal error
    gtfo ("globe_tile_lat_bounds : Illegal GT index z=%d", z);
  }
}

void qgt_lon_bounds (int x, double &w, double &e)
{
  globe_tile_lon_bounds (x/2, w, e);
  // Globe tile bounds are valid
  double half = w + ((e - w) / 2.0);
  if ((x % 2) == 0) {
    // This is the western QGT within the GT
    e = half;
  } else {
    // This is the eastern QGT within the GT
    w = half;
  }
}

void qgt_lat_bounds (int z, double &s, double &n)
{
  globe_tile_lat_bounds (z/2, s, n);
  // Globe tile bounds are valid
  double half = s + ((n - s) / 2);
  if ((z % 2) == 0) {
    // This is the southern QGT within the GT
    n = half;
  } else {
    // This is the northern QGT within the GT
    s = half;
  }
}

#define ARCSEC2RAD    ((2.0 * PI) / 3600.0)

//
// Format a position struct into human-readable text
//
void FormatPosition (SPosition pos, char* s)
{
  double latSec = pos.lat;
  char latHemi;
  if (latSec >= 0) {
    latHemi = 'N';
  } else {
    latHemi = 'S';
    latSec = fabs(latSec);
  }

  double lonSec = pos.lon;
  char lonHemi;
  if (lonSec >= (180.0 * 3600)) {
    lonHemi = 'W';
    lonSec = (360 * 3600) - lonSec;
  } else {
    lonHemi = 'E';
  }

  double latDeg = floor (latSec / 3600);
  latSec -= (latDeg * 3600);
  double latMin = floor (latSec / 60);
  latSec -= (latMin * 60);

  double lonDeg = floor (lonSec / 3600);
  lonSec -= (lonDeg * 3600);
  double lonMin = floor (lonSec / 60);
  lonSec -= (lonMin * 60);

  sprintf (s, "%c %2.0f %2.0f'%7.4f\"  %c %2.0f %2.0f'%7.4f\"",
    latHemi, latDeg, latMin, latSec, lonHemi, lonDeg, lonMin, lonSec);
}


//
// Format Right Ascension/Declination (both in radians) to a human-readable string
//
void FormatRADec (double ra, double dec, char* s)
{
  ra = WrapTwoPi (ra) * 24.0 / (2 * SGD_PI);
  dec = RadToDeg (dec);

  // Convert RA to hrs/min/sec
  double ra_hr, ra_min, ra_sec;
  ra_hr = floor (ra);
  ra -= ra_hr;
  ra_min = floor (ra * 60);
  ra -= ra_min / 60;
  ra_sec = ra * 3600;

  // Convert Dec to deg/arcmin/arcsec
  double dec_deg, dec_min, dec_sec;
  dec_deg = floor (dec);
  dec -= dec_deg;
  dec_min = floor (dec * 60);
  dec -= dec_min / 60;
  dec_sec = dec * 3600;

  sprintf (s, "%2.0fh %2.0fm %7.3fs  %3.0fd %2.0fm %7.3fs",
    ra_hr, ra_min, ra_sec, dec_deg, dec_min, dec_sec);
}


//
// Format Sidereal Time (decimal hours) to a human-readable sting
//
void FormatSiderealTime (double st, char *s)
{
  double st_hr, st_min, st_sec;
  st_hr = floor (st);
  st -= st_hr;
  st_min = floor (st * 60);
  st -= st_min / 60;
  st_sec = st * 3600;

  sprintf (s, "%2.0f:%02.0f:%05.3fs", st_hr, st_min, st_sec);
}


//
// Convert CIE Yxy triplet to CIE XYZ
//
// Input parameters:
//  Y CIE Yxy luminance Y value   [0..100]
//  x CIE Yxy chrominance x value   [0..1]
//  y CIE Yxy chrominance y value   [0..1]
//
void CIE_Yxy_to_XYZ (const SCIE in, SCIE &out)
{
  out.XYZ.X = in.Yxy.x * (in.Yxy.Y / in.Yxy.y);
  out.XYZ.Y = in.Yxy.Y;
  out.XYZ.Z = (1 - in.Yxy.x - in.Yxy.y) * (in.Yxy.Y / in.Yxy.y);
}


//
// Convert CIE XYZ triplet to RGB colour using the D65 white point.
//
void CIE_XYZ_to_RGB_D65 (const SCIE XYZ, sgdVec3 &RGB)
{
  // Now multiply XYZ vector by D65 white point matrix
  float mWhite[3][3];
  mWhite[0][2] =  3.240479;
  mWhite[0][1] = -1.537150;
  mWhite[0][0] = -0.498535;

  mWhite[1][2] = -0.969256;
  mWhite[1][1] =  1.875992;
  mWhite[1][0] =  0.041556;

  mWhite[2][2] =  0.055648;
  mWhite[2][1] = -0.204043;
  mWhite[2][0] =  1.057311;

  RGB[0] = ((mWhite[0][2] * XYZ.XYZ.X) + (mWhite[0][1] * XYZ.XYZ.Y) + (mWhite[0][0] * XYZ.XYZ.Z)) / 100.0;
  RGB[1] = ((mWhite[1][2] * XYZ.XYZ.X) + (mWhite[1][1] * XYZ.XYZ.Y) + (mWhite[1][0] * XYZ.XYZ.Z)) / 100.0;
  RGB[2] = ((mWhite[2][2] * XYZ.XYZ.X) + (mWhite[2][1] * XYZ.XYZ.Y) + (mWhite[2][0] * XYZ.XYZ.Z)) / 100.0;
}


//
// Returns the cartesian coordinate scaling factor for terrain tiles at
//   a particular latitude.  See TerrainModel.doc for details.
//
float TerrainScale (SPosition pos)
{
  return cos(DegToRad(pos.lat / 3600.0));
}

//
// Returns the size, in feet, of a standardized globe tile
//
float GlobeTileDistance (void)
{
  return NmToFeet ((360.0 / 256.0) * 60.0);
}

//
// The following functions convert an SPosition struct into a set of cartesian
//   world coordinates for rendering.  See Terrain.doc for details.
//
// Find cartesian coordinates of a position relative to the center of the
//   given globe tile
//
SVector PosToFlatCartesian (SPosition pos, int x, int z)
{
  // Get lat and lon boundaries for this globe tile, convert to arcsec
  double s, n, w, e;
  globe_tile_lat_bounds (z, s, n);
  globe_tile_lon_bounds (x, w, e);
  double latScale = (360.0 / 256.0) / fabs(n - s);
  n *= 3600;
  s *= 3600;
  w *= 3600;
  e *= 3600;

  // Calculate midpoint lat and lon
  double midLon = w + ((e - w) / 2);
  double midLat;
  if (pos.lat >= 0) {
    // Northern hemisphere
    midLat = s + ((n - s) / 2);
  } else {
    // Southern hemisphere
    midLat = n + ((s - n) / 2);
  }

  // Calculate lon and lat arcsecs from the midpoint, adjusted to universal globe tile
  double dLon = pos.lon - midLon;
  double dLat = pos.lat - midLat;
  dLat *= latScale;

  // Convert to feet using 1 arcmin = 1 nm
  double easting = NmToFeet(dLon / 60.0);
  double northing = NmToFeet(dLat / 60.0);

  SVector v;
  v.x = easting;
  v.y = northing;
  v.z = pos.alt;

  return v;
}

SVector PosToFlatCartesian (SPosition pos)
{
  int gx, gz;
  lat_lon_to_globe_tile (pos.lat/3600.0, pos.lon/3600.0, gx, gz);
  return PosToFlatCartesian (pos, gx, gz);
}

SVector PosToScaledFlatCartesian (SPosition pos)
{
  SVector v;

  int gx, gz;
  lat_lon_to_globe_tile (pos.lat/3600.0, pos.lon/3600.0, gx, gz);
  v = PosToFlatCartesian (pos, gx, gz);

  float scale = TerrainScale (pos);
  v.x *= scale;
  v.y *= scale;

  return v;
}

//
// The normalized QGT size is 42.1875 nm square (equivalent to half of 1.40625 deg),
//   so if QGT center coordinate is 0,0 then the bounds are +/- 21.09375 nm.  This
//   macro defines the QGT bound in feet.
//
//#define QGT_FLAT_CARTESIAN_BOUND_FT   (1.28165E+5)
#define QGT_FLAT_CARTESIAN_BOUND_FT   (2.5633E+5)


SVector PosToFlatCartesianQgt (SPosition pos, int x, int z)
{
  // Get lat and lon boundaries for this QGT, convert to arcsec
  double s, n, w, e;
  qgt_lat_bounds (z, s, n);
  qgt_lon_bounds (x, w, e);
  
  // Latitude scale adjusts actual QGT size to normalized QGT size
  double latScale = (360.0 / 256.0) / (2.0 * fabs(n - s));
  
  // Convert boundary lat/lon to arcseconds
  n *= 3600;
  s *= 3600;
  w *= 3600;
  e *= 3600;

  // Calculate midpoint lat and lon
  double midLon = w + ((e - w) / 2.0);
  double midLat;
  if (pos.lat >= 0) {
    // Northern hemisphere
    midLat = s + ((n - s) / 2.0);
  } else {
    // Southern hemisphere
    midLat = n + ((s - n) / 2.0);
  }

  // Calculate lon and lat arcsecs from the midpoint, adjusted to universal QGT
  double dLon = pos.lon - midLon;
  double dLat = pos.lat - midLat;
  dLat *= latScale;

  // Convert to feet using 1 arcmin = 1 nm
  double easting = NmToFeet(dLon / 60.0);
  double northing = NmToFeet(dLat / 60.0);

  SVector v;
  v.x = easting;
  v.y = northing;
  v.z = pos.alt;

  return v;
}

SVector PosToFlatCartesianQgt (SPosition pos)
{
  int x, z;
  lat_lon_to_qgt (pos.lat, pos.lon, x, z);
  return PosToFlatCartesianQgt (pos, x, z);
}

SVector PosToScaledFlatCartesianQgt (SPosition pos)
{
  SVector v = PosToFlatCartesianQgt (pos);

  float scale = TerrainScale (pos);
  v.x *= scale;
  v.y *= scale;

  return v;
}

Documentation generated by