/*
 * renderer.cpp
 *
 * $Id: renderer.cpp,v 0.9 1999/12/06 16:13:26 thorsten Exp thorsten $
 *
 * This file is part of XGlobe. See README for details.
 *
 * Copyright (C) 1998 Thorsten Scheuermann
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public Licenses as published by
 * the Free Software Foundation.
 *
 * $Log: renderer.cpp,v $
 * Revision 0.9  1999/12/06 16:13:26  thorsten
 * -added -rot switch
 *
 * Revision 0.8  1999/12/05 19:18:30  thorsten
 * by Andrew:
 * - added cloud layer rendering
 *
 * Revision 0.7  1999/09/19 19:45:53  thorsten
 * by Andrew Sumner:
 * - added -shade_area switch
 *
 * Revision 0.6  1999/08/06 16:14:46  thorsten
 * Forgot to initialize a variable, which could cause a segfault.
 *
 * Revision 0.5  1999/07/19 12:44:39  thorsten
 * - works with QT 2
 * - added -term switch
 * - nicer output when using night maps
 *
 * Revision 0.4  1999/07/13 17:40:02  thorsten
 * 8 bit images can now be used without increasing colordepth to 32 bit
 * can now draw stars in the background
 * by Andrew Sumner:
 * can now display a background image
 * globe does no longer need to be centered
 *
 * Revision 0.2  1999/01/08 17:25:28  thorsten
 * added grid drawing code
 *
 * Revision 0.2  1999/01/06 13:21:33  thorsten
 * added optional display of a grid
 *
 * Revision 0.1  1998/12/10 20:08:57  thorsten
 * initial revision
 *
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <qapplication.h>
#include <qpainter.h>
#include <qpixmap.h>
#include <qdatetime.h>
#include <qfile.h>
#include <string.h>
#include "config.h"
#include "renderer.h"
#include "sunpos.h"

#include "marker.xpm"

#define MAX(x, y) ((x) < (y) ? (y) : (x))
#define MIN(x, y) ((x) > (y) ? (y) : (x))

/* ------------------------------------------------------------------------*/

Renderer::Renderer(const QSize &size, const char *mapfile)
{
  markerlist = NULL;
  show_markers = FALSE;
  use_nightmap = FALSE;
  use_cloudmap = FALSE;
  use_backimage = FALSE;
  map = NULL;
  mapnight = NULL;
  
  if(!mapfile)
  {
    mapFileName = XGLOBE_LIB_DIR;
    mapFileName += "/map.bmp";
  }
  else
  {
    if(QFile::exists(mapfile))
      mapFileName = mapfile;
    else
    {
      mapFileName = XGLOBE_LIB_DIR;
      mapFileName += "/";
      mapFileName += mapfile;
      if(!QFile::exists(mapFileName))
      {
        fprintf(stderr, "File \"%s\" not found!\n", (const char *)mapfile);
        ::exit(1);
      }
    }
  }

  renderedImage = new QImage(size, 32);
  if(!renderedImage)
  {
    fprintf(stderr, "Not enough memory for offscreen image buffer!\n");
    ::exit(1);
  }

  map = loadImage((const char *)mapFileName);
  
  // fprintf(stderr, "Map size: %dx%d\n", map->width(), map->height());
 
  this->radius = 1000.;
  this->view_long = 0.;
  this->view_lat = 0.;
  this->sun_long = 0.;
  this->sun_lat = 0.;
  this->fov = 0.5*PI/180.;
  this->zoom = 0.9;
  this->ambientRed = 0.15;
  this->ambientGreen = 0.15;
  this->ambientBlue = 0.15;
  this->markerpixmap = new QPixmap((const char **)marker_xpm);
  ASSERT(markerpixmap != NULL);
  this->show_label = TRUE;
  this->gridtype = NO_GRID;
  this->d_gridline = 15.0*PI/180.;
  this->d_griddot = PI/180.;
  this->cloud_filter = 110;
  this->shift_x = 0;
  this->shift_y = 0;
  this->show_stars = true;
  this->star_frequency = 0.002;
  this->trans = 0.0;
  this->rot = 0.0;
  
  calcDistance();
}

/* ------------------------------------------------------------------------*/

QImage *Renderer::loadImage(const QString &name)
{
  QImage *m;

  m = new QImage();
  if(!m)
  {
    fprintf(stderr, "Not enough memory for map!\n");
    ::exit(1);
  }

  if(!m->load(name))
  {
    fprintf(stderr, "Error while opening map \"%s\"!\n",
            (const char *)name);
    ::exit(1);
  }

  if(m->depth() < 8)
    *m = m->convertDepth(8);
  /*
  if(m->depth() != 32)
    *m = m->convertDepth(32);
  */
  
  return m;
}

/* ------------------------------------------------------------------------*/

int Renderer::loadNightMap(const char *nmapfile)
{
  if(use_nightmap)                // we already have a night map!
    return 1;
  
  if(nmapfile)
  {
    if(QFile::exists(nmapfile))
      mapNightFileName = nmapfile;
    else
    {
      mapNightFileName = XGLOBE_LIB_DIR;
      mapNightFileName += "/";
      mapNightFileName += nmapfile;
      if(!QFile::exists(mapNightFileName))
      {
        fprintf(stderr, "File \"%s\" not found!\n", (const char *)nmapfile);
        ::exit(1);
      }
    }
  }
  else
  {
    mapNightFileName = XGLOBE_LIB_DIR;
    mapNightFileName += "/mapnight.bmp";
  }
  
  mapnight = loadImage(mapNightFileName);

  use_nightmap = TRUE;
  return 1;
}

/* ------------------------------------------------------------------------*/

int Renderer::loadCloudMap(const char *cmapfile, int cf)
{
  if(use_cloudmap)                // we already have a cloud map!
    return 1;

  cloud_filter = cf;
  
  if(cmapfile)
  {
    if(QFile::exists(cmapfile))
      mapCloudFileName = cmapfile;
    else
    {
      mapCloudFileName = XGLOBE_LIB_DIR;
      mapCloudFileName += "/";
      mapCloudFileName += cmapfile;
      if(!QFile::exists(mapCloudFileName))
      {
        fprintf(stderr, "File \"%s\" not found!\n", (const char *)cmapfile);
        ::exit(1);
      }
    }
  }
  else
    return 1;
  
  mapcloud = loadImage(mapCloudFileName);
  if (!mapcloud)
  {
    fprintf(stderr, "Error loading cloud mapfile \"%s\"\n", (const char *)cmapfile);
    ::exit(1);
  }
  int endy = mapcloud->height();
  int endx = mapcloud->width();

  if (mapcloud->depth() == 8)
    *mapcloud = mapcloud->convertDepth(32);

  int sb, sg, sr;
  int lb=0, lg=0, lr=0;
  unsigned char *p1;
  unsigned int *q1;

  /* filter out the pink continent outlines */
  for(int py=0; py<endy; py++)
  {
    q1 = (unsigned int *)mapcloud->scanLine(py);
    for (int px=0; px<endx; px++)
    {
      p1 = (unsigned char *)mapcloud->scanLine(py) + 4*px;
      //q1 = (unsigned int *)mapcloud->scanLine(py) + 4*px;
  
      sb = qBlue(*(unsigned int *)p1);
      sg = qGreen(*(unsigned int *)p1);
      sr = qRed(*(unsigned int *)p1);

      if (sb == 255 && sr == 255 && sg == 0)
      {
        /* continent outline */
        sb = lb;
        sg = lg;
        sr = lr;
        *q1 = qRgb((int)sr,(int)sg,(int)sb);
      }
      else if (sb > sr && sb > sg)
      {
        /* polar area, no cloud data */
        sb = sr = sg = 0;
        *q1 = qRgb((int)sr,(int)sg,(int)sb);
      }
      /*
      else if (sb < cloud_filter || sg < cloud_filter || sr < cloud_filter)
      {
        // cloud is too dark to show 
        sb = sr = sg = 0;
        *q1 = qRgb((int)sr,(int)sg,(int)sb);
      }
      */

      lb = sb;
      lr = sr;
      lg = sg;
      q1++;
    }
  }

  use_cloudmap = TRUE;
  return 1;
}

/* ------------------------------------------------------------------------*/

int Renderer::loadBackImage(const char *imagefile, bool tld)
{
  QImage *bi;

  if (use_backimage == 1)
    return 1;

  tiled = tld;
  if(imagefile)
  {
    if(QFile::exists(imagefile))
      backImageFileName = imagefile;
    else
    {
      backImageFileName = XGLOBE_LIB_DIR;
      backImageFileName += "/";
      backImageFileName += imagefile;
      if(!QFile::exists(backImageFileName))
      {
        fprintf(stderr, "File \"%s\" not found!\n", (const char *)imagefile);
        ::exit(1);
      }
    }
  }
  else
  {
    backImageFileName = XGLOBE_LIB_DIR;
    backImageFileName += "/back.bmp";
  }
  
  bi = loadImage(backImageFileName);

  if (tiled == FALSE)
    backImage = bi->smoothScale(renderedImage->width(), renderedImage->height());
  else
    backImage = bi->copy();
  bi->~QImage();
  
  use_backimage = TRUE;
  return 1;
}

/* ------------------------------------------------------------------------*/

Renderer::~Renderer()
{
  if(renderedImage)
    delete renderedImage;
  if(map)
    delete map;

  if(mapnight)
    delete mapnight;

  delete markerpixmap;
}

/* ------------------------------------------------------------------------*/

void Renderer::setViewPos(double lat, double lon)
{
  while(lat >= 360.)
    lat -= 360.;
  while(lat <= -360.)
    lat += 360.;
  if(lat > 90.)
  {
    lat = 90. - (lat - 90.);
    lon += 180.;
  }
  if(lat < -90.)
  {
    lat = -90. + (-lat - 90.);
    lon += 180.;
  }
  
  while(lon >= 360.)
    lon -= 360.;
  while(lon <= -360.)
    lon += 360.;
  if(lon > 180.)
    lon = -180. + (lon - 180.);
  if(lon < -180.)
    lon = 180. + (lon + 180.);
  
  view_lat = lat*PI/180.;
  view_long = lon*PI/180.;
}

/* ------------------------------------------------------------------------*/

double Renderer::getViewLat()
{
  return view_lat*180./PI;
}

/* ------------------------------------------------------------------------*/

double Renderer::getViewLong()
{
  return view_long*180./PI;
}

/* ------------------------------------------------------------------------*/

void Renderer::setRotation(double r)
{
  rot = r*PI/180.;
}

/* ------------------------------------------------------------------------*/

double Renderer::getRotation()
{
  return rot*180./PI;
}

/* ------------------------------------------------------------------------*/

double Renderer::getSunLat()
{
  return sun_lat*180./PI;
}

/* ------------------------------------------------------------------------*/

double Renderer::getSunLong()
{
  return sun_long*180./PI;
}

/* ------------------------------------------------------------------------*/

void Renderer::setZoom(double z)
{
  zoom = z;
  calcDistance();
}

/* ------------------------------------------------------------------------*/

double Renderer::getZoom()
{
  return zoom;
}

/* ------------------------------------------------------------------------*/

void Renderer::setMarkerList(MarkerList *l)
{
  ASSERT(l != NULL);
  markerlist = l;
}

/* ------------------------------------------------------------------------*/

void Renderer::showMarkers(bool show)
{
  show_markers = show;
}

/* ------------------------------------------------------------------------*/

void Renderer::showLabel(bool show)
{
  show_label = show;
}

/* ------------------------------------------------------------------------*/

void Renderer::setLabelPos(int x, int y)
{
  label_x = x;
  label_y = y;
}

/* ------------------------------------------------------------------------*/

void Renderer::setShadeArea(double area)
{
  shade_area = area;
}

/* ------------------------------------------------------------------------*/

void Renderer::setAmbientRGB(double ambient_red, double ambient_green,
                             double ambient_blue)
{
  ambientRed = ambient_red;
  ambientGreen = ambient_green;
  ambientBlue = ambient_blue;
}

/* ------------------------------------------------------------------------*/

void Renderer::setFov(double fov)
{
  this->fov = fov*PI/180.;
  calcDistance();
}

/* ------------------------------------------------------------------------*/

void Renderer::setTime(time_t t)
{
  time_to_render = t;
  calcLightVector();  // calc. current sun position
}

/* ------------------------------------------------------------------------*/

time_t Renderer::getTime()
{
  return time_to_render;
}

/* ------------------------------------------------------------------------*/

void Renderer::setNumGridLines(int num)
{
  d_gridline = PI/(2.0*num);
}

/* ------------------------------------------------------------------------*/

int Renderer::getNumGridLines()
{
  return (int)(PI/2*d_gridline);
}
  
/* ------------------------------------------------------------------------*/

void Renderer::setNumGridDots(int num)
{
  d_griddot = 2.0*PI/num;
}

/* ------------------------------------------------------------------------*/

int Renderer::getNumGridDots()
{
  return (int)(2.0*PI/d_griddot);
}

/* ------------------------------------------------------------------------*/

void Renderer::setGridType(int type)
{
  gridtype = type;
  if((type < 0) || (type > NICE_GRID))
    gridtype = NO_GRID;
}

/* ------------------------------------------------------------------------*/

int Renderer::getGridType()
{
  return gridtype;
}

/* ------------------------------------------------------------------------*/

void Renderer::setShift(int x, int y)
{
  shift_x = x;
  shift_y = y;
}

/* ------------------------------------------------------------------------*/

int Renderer::getShiftX()
{
  return shift_x;
}

/* ------------------------------------------------------------------------*/

int Renderer::getShiftY()
{
  return shift_y;
}

/* ------------------------------------------------------------------------*/

void Renderer::showStars(bool show)
{
  show_stars = show;
}

/* ------------------------------------------------------------------------*/

void Renderer::setStarFrequency(double f)
{
  star_frequency = f;
}

/* ------------------------------------------------------------------------*/

double Renderer::getStarFrequency()
{
  return star_frequency;
}

/* ------------------------------------------------------------------------*/

void Renderer::setTransition(double t)
{
  trans = t;
  if(trans >= 1.0)
    trans = 0.9999;
  else if(trans < 0.0)
    trans = 0.0;
}

/* ------------------------------------------------------------------------*/

double Renderer::getTransition()
{
  return trans;
}

/* ------------------------------------------------------------------------*/

void Renderer::calcDistance()
{
  double x;
  double tan_a;
  
  // distance of camera to projection plane
  proj_dist = MIN(renderedImage->width(), renderedImage->height())/tan(fov);

  x = zoom*MIN(renderedImage->width(), renderedImage->height())/2.;
  tan_a = x/proj_dist;
  // distance of camera camera to center of earth ( = coordinate origin)
  center_dist = radius/sin(atan(tan_a));
}

/* ------------------------------------------------------------------------*/

void Renderer::renderFrame()
{
  double dir_x, dir_y, dir_z;      // direction of cast ray
  double hit_x, hit_y, hit_z;      // hit position on earth surface
  double hit2_x, hit2_y, hit2_z;   // mirrored hit position on earth surface
  double sp_x, sp_y, sp_z;         // intersection point of globe and ray
  double a, b, c;                  // coeff. of quadratic equation
  double radikand;
  double wurzel;
  double r;                        // r'
  double radiusq = radius*radius;
  double s1, s2, s;                // distance between intersections and
                                   // camera position
  double longitude, latitude;      // coordinates of hit position
  double light_angle;              // cosine of angle between sunlight and
                                   // surface normal
  int startx, endx;                // the region to be painted
  int starty, endy;
  int real_startx, real_endx;
  int temp;

  unsigned int *p;                 // pointer to current pixel
  unsigned int *q;

  int half_width = renderedImage->width()/2 + renderedImage->width()%2 -1;
  
  // clear image
  for(int i=0; i < renderedImage->height(); i++)
  {
    p = (unsigned int *)renderedImage->scanLine(i);
    memset(p, 0, renderedImage->bytesPerLine());
  }

  if (use_backimage)
    copyBackImage();

  if(show_stars)
    drawStars();
  
  // rotation matrix
  double m11 = cos(rot)*cos(view_long)-sin(view_lat)*sin(view_long)*sin(rot);
  double m12 = -sin(rot)*cos(view_long)-sin(view_lat)*sin(view_long)*cos(rot);
  double m13 = cos(view_lat)*sin(view_long);
  double m21 = sin(rot)*cos(view_lat);
  double m22 = cos(view_lat)*cos(rot);
  double m23 = sin(view_lat);
  double m31 = -sin(view_long)*cos(rot)-sin(view_lat)*cos(view_long)*sin(rot);
  double m32 = sin(view_long)*sin(rot)-sin(view_lat)*cos(view_long)*cos(rot);
  double m33 = cos(view_lat)*cos(view_long);

  
  dir_z = -proj_dist;

  // indifferent coeff.
  c = center_dist*center_dist - radiusq;

  // calc. radius of projected sphere
  b = 2*center_dist*dir_z;
  radius_proj = (int)sqrt(b*b/(4*c) - dir_z*dir_z);
  
  startx = (renderedImage->width()/2 - radius_proj -1);
  startx = (startx < 0) ? 0 : startx;
  endx = renderedImage->width() - startx-1;
  starty = (renderedImage->height()/2 - radius_proj -1);
  starty = (starty < 0) ? 0 : starty;
  endy = renderedImage->height() - starty-1;

  for(int py=starty; py<=endy; py++)
  {
    // handle any paint events waiting in the queue
    qApp->processEvents();

    temp = radius_proj*radius_proj - (py-renderedImage->height()/2) *
           (py-renderedImage->height()/2);

    if(temp >= 0)
      startx = (renderedImage->width()/2 - (int)sqrt(temp));
    else
      startx = (renderedImage->width()/2);
    
    startx = (startx < 0) ? 0 : startx;
    endx = renderedImage->width() - startx-1;

    // calculate offset into image data
    if (py + shift_y < 0 || py + shift_y >= renderedImage->height())
      continue;
    real_startx = startx + shift_x;
    real_startx = (real_startx < 0 ? 0 : real_startx);
    real_startx = (real_startx >= renderedImage->width() ?
                   renderedImage->width()-1 : real_startx);

    real_endx = endx + shift_x;
    real_endx = (real_endx < 0 ? 0 : real_endx);
    real_endx = (real_endx >= renderedImage->width() ?
                 renderedImage->width()-1 : real_endx);

    p = (unsigned int *)renderedImage->scanLine(py + shift_y)+real_startx;
    q = (unsigned int *)renderedImage->scanLine(py + shift_y)+real_endx;

    if(rot == 0.)           // optimization when using no rotation
      endx = half_width;

    for(int px=startx; px<=endx; px++)
    {
      dir_x =  (px - renderedImage->width()/2);
      dir_y = (-py + renderedImage->height()/2);

      a = dir_x*dir_x + dir_y*dir_y + dir_z*dir_z;
      b = 2*center_dist*dir_z;
      // c constant, see above

      radikand = b*b-4*a*c;			// what's under the sq.root when solving the
                                // quadratic equation
      if(radikand >= 0.)				// solution exists <=> intersection
      {
        wurzel = sqrt(radikand);
        s1 = (-b + wurzel)/(2.*a);
        s2 = (-b - wurzel)/(2.*a);
        s = (s1<s2) ? s1 : s2;    // smaller solution belongs to nearer 
                                  // intersection
        sp_x = s * dir_x;         // sp = camera pos + s*dir
        sp_y = s * dir_y;
        sp_z = center_dist + s * dir_z;

        hit_x =  m11*sp_x + m12*sp_y + m13*sp_z;
        hit_y =  m21*sp_x + m22*sp_y + m23*sp_z;
        hit_z =  m31*sp_x + m32*sp_y + m33*sp_z;

        if(rot == 0.)             // optimization when using no rotation
        {
          hit2_x = -m11*sp_x + m12*sp_y + m13*sp_z;
          hit2_y = -m21*sp_x + m22*sp_y + m23*sp_z;
          hit2_z = -m31*sp_x + m32*sp_y + m33*sp_z;
        }
        
        longitude = atan(hit_x/hit_z);
        if(hit_z < 0.)
          longitude = PI + longitude;
        
        r = (double)sqrt(hit_x*hit_x + hit_z*hit_z);
        latitude = atan(-hit_y/r);
      
        light_angle = (light_x*hit_x + light_y*hit_y + light_z*hit_z)/radius;
        light_angle = pow(light_angle, 1.0-trans);
        
        // Set pixel in image
        *p++ = getPixelColor(longitude, latitude, light_angle);

        // only when using no rotation:
        // mirror the left half-circle of the globe: we need a new position
        // and have to recalculate the light intensity
        if(rot == 0.)
        {
          light_angle = (light_x*hit2_x + light_y*hit2_y + light_z*hit2_z)/
                        radius;
          light_angle = pow(light_angle, 1.0-trans);
          *q-- = getPixelColor(2*view_long - longitude, latitude, light_angle);
        }
      }
      else
      {
        p++;
        q--;
      }   
		}
	}

  if(gridtype != NO_GRID)
    drawGrid();
  
  if(show_markers)
    drawMarkers();

  if(show_label)
    drawLabel();
}

/* ------------------------------------------------------------------------*/

void Renderer::copyBackImage()
{
  unsigned int *p, *bp;
  unsigned char *c_bp;
  unsigned int y, x, by, bx;
  unsigned int mywidth = renderedImage->width(), myheight = renderedImage->height();
  unsigned int bwidth = backImage.width(), bheight = backImage.height();
  
  for (y=0, by=0; y<myheight; y++, by++)
  {
    if (by >= bheight)
      by = 0;
    
    p = (unsigned int *)renderedImage->scanLine(y);
    bp = (unsigned int *)backImage.scanLine(by);
    c_bp = (unsigned char *)backImage.scanLine(by);

    if(backImage.depth() == 32)
    {
      for (x=0, bx=0; x<mywidth; x++, bx++)
      {
        if (bx >= bwidth)
        {
          bx = 0;
          bp = (unsigned int *) backImage.scanLine(by);
        }
        *p++ = *bp++;
      }
    }
    else
    {
      for (x=0, bx=0; x<mywidth; x++, bx++)
      {
        if (bx >= bwidth)
        {
          bx = 0;
          c_bp = (unsigned char *) backImage.scanLine(by);
        }
        *p++ = backImage.color(*c_bp++);
      }
    }  
  }
}


/* ------------------------------------------------------------------------*/

void Renderer::calcLightVector()
{
  SunPos::GetSunPos(time_to_render, &sun_lat, &sun_long);

  light_x = cos(sun_lat)*sin(sun_long);
  light_y = sin(sun_lat);
  light_z = cos(sun_lat)*cos(sun_long);
}

/* ------------------------------------------------------------------------*/

unsigned int Renderer::getPixelColor(double longitude, double latitude,
                                     double angle)
{
  float r, g, b, shade_angle;
  int dr, dg, db;                  // rgb values of day and night pixel
  int nr, ng, nb;  
  int cr, cg, cb;
  bool cloud_used = FALSE;

  if (shade_area)
    shade_angle = angle / shade_area;
  else
    shade_angle = 1.0;

  getMapColorLinear(map, longitude, latitude,
                    &dr, &dg, &db);
  if(use_nightmap)
    getMapColorLinear(mapnight, longitude, latitude,
                      &nr, &ng, &nb);
  if(use_cloudmap && mapcloud)
  {
    getMapColorLinear(mapcloud, longitude, latitude,
                      &cr, &cg, &cb);
    if (cr>cloud_filter && cb>cloud_filter && cg>cloud_filter)
    {
      if (angle > 0.0 || !use_nightmap)
      {
        if (dr < cr)
          dr += (int) ((double) (cr - dr) / 1.5);
        else
          dr = cr;
        if (dg < cg)
          dg += (int) ((double) (cg - dg) / 1.2);
        else
          dg = cr;
        if (db < cb)
          db += (int) ((double) (cb - db) / 1.5);
        else
          dg = cr;
      }
      if (use_nightmap)
      {
        cr /= 5;
        cb /= 5;
        cg /= 5;
        
        if (nr < cr)
          nr = (int) ((cr - ((double) (cr - nr) / 2)) / (1.-ambientRed));
        else
          nr = (int) ((cr + ((double) (nr - cr) / 4)) / (1.-ambientRed));
        if (ng < cg)
          ng = (int) ((cg - ((double) (cg - ng) / 2)) / (1.-ambientGreen));
        else
          ng = (int) ((cg + ((double) (ng - cg) / 4)) / (1.-ambientGreen));
        if (nb < cb)
          nb = (int) ((cb - ((double) (cb - nb) / 2)) / (1.-ambientBlue));
        else
          nb = (int) ((cb + ((double) (nb - cb) / 4)) / (1.-ambientBlue));
      }
      cloud_used = TRUE;
    }
  }

  if(use_nightmap)
  {
    if (angle > shade_area)
    {
      r = dr*ambientRed + dr*(1.-ambientRed);
      g = dg*ambientGreen + dg*(1.-ambientGreen);
      b = db*ambientBlue + db*(1.-ambientBlue);
    }
    else if(angle > 0.1) 
    {
      r = dr*ambientRed + dr*shade_angle*(1.-ambientRed);
      g = dg*ambientGreen + dg*shade_angle*(1.-ambientGreen);
      b = db*ambientBlue + db*shade_angle*(1.-ambientBlue);
    }
    else if(fabs(angle) <= 0.1) 
    {
      double x;

      if(angle > 0.)
      {
        r = dr*ambientRed + dr*shade_angle*(1.-ambientRed);
        g = dg*ambientGreen + dg*shade_angle*(1.-ambientGreen);
        b = db*ambientBlue + db*shade_angle*(1.-ambientBlue);
      }
      else
      {
        r = dr*ambientRed;
        g = dg*ambientGreen;
        b = db*ambientBlue;
      }
      
      x = -5.0*angle+0.5;
      r = x*nr + (1.0-x)*r;
      g = x*ng + (1.0-x)*g;
      b = x*nb + (1.0-x)*b;
    }
    else
    {
      r = nr;
      g = ng;
      b = nb;
    }
  }
  else
  {
    if (angle > shade_area)
    {
      r = dr*ambientRed + dr*(1.-ambientRed);
      g = dg*ambientGreen + dg*(1.-ambientGreen);
      b = db*ambientBlue + db*(1.-ambientBlue);
    }
    else if(angle > 0.)
    {
      r = dr*ambientRed + dr*shade_angle*(1.-ambientRed);
      g = dg*ambientGreen + dg*shade_angle*(1.-ambientGreen);
      b = db*ambientBlue + db*shade_angle*(1.-ambientBlue);
    }
    else
    {
      r = dr*ambientRed;
      g = dg*ambientGreen;
      b = db*ambientBlue;
    }
  }

  return qRgb((int)r,(int)g,(int)b);
}

/* ------------------------------------------------------------------------*/

void Renderer::getMapColorLinear(QImage *m, double longitude, double latitude,
                                 int *r, int *g, int *b)
{
  unsigned char *p1;
  unsigned char *p2;

  if(longitude < 0.)
  	longitude += 2*PI;

  latitude += PI/2;

  longitude += PI;
  if(longitude >= 2*PI)
  	longitude -= 2*PI;

  int x11 = (int)(longitude*m->width()/(2*PI));
  int y11 = (int)(latitude*m->height()/PI);
  int x12 = (x11+1) % m->width();
  int y12 = y11;
  int x21 = x11;
  int y21 = (y11+1) % m->height();
  int x22 = x12;
  int y22 = y21;
  double dx = (longitude*m->width()/(2*PI))-x11;
  double dy = (latitude*m->height()/PI)-y11;
  x11 %= m->width();
  y11 %= m->height();             // this isn't really correct, but
                                    // nobody will see the difference...
  double r1,g1,b1, r2,g2,b2;

  int sb, sg, sr;
  int tb, tg, tr;

  // offset into map pixel data
  if(m->depth() == 32)
  {
    p1 = (unsigned char *)m->scanLine(y11) + 4*x11;
    p2 = (unsigned char *)m->scanLine(y12) + 4*x12;
  }
  else //m->depth() == 8
  {
    p1 = (unsigned char *)m->scanLine(y11) + x11;
    p2 = (unsigned char *)m->scanLine(y12) + x12;
  }
  
  if(m->depth() == 32)
  {
    sb = qBlue(*(unsigned int *)p1);
    sg = qGreen(*(unsigned int *)p1);
    sr = qRed(*(unsigned int *)p1);
    tb = qBlue(*(unsigned int *)p2);
    tg = qGreen(*(unsigned int *)p2);
    tr = qRed(*(unsigned int *)p2);
  }
  else
  {
    QRgb c = m->color(*p1);
    sb = qBlue(c);
    sg = qGreen(c);
    sr = qRed(c);
    c = m->color(*p2);
    tb = qBlue(c);
    tg = qGreen(c);
    tr = qRed(c);
  }
  
  r1 = (tr-sr)*dx+sr;
  g1 = (tg-sg)*dx+sg;
  b1 = (tb-sb)*dx+sb;
  
  if(m->depth() == 32)
  {
    p1 = (unsigned char *)m->scanLine(y21) + 4*x21;
    p2 = (unsigned char *)m->scanLine(y22) + 4*x22;
  }
  else
  {
    p1 = (unsigned char *)m->scanLine(y21) + x21;
    p2 = (unsigned char *)m->scanLine(y22) + x22;
  }

    
  if(m->depth() == 32)
  {
    sb = qBlue(*(unsigned int *)p1);
    sg = qGreen(*(unsigned int *)p1);
    sr = qRed(*(unsigned int *)p1);
    tb = qBlue(*(unsigned int *)p2);
    tg = qGreen(*(unsigned int *)p2);
    tr = qRed(*(unsigned int *)p2);
  }
  else
  {
    QRgb c = m->color(*p1);
    sb = qBlue(c);
    sg = qGreen(c);
    sr = qRed(c);
    c = m->color(*p2);
    tb = qBlue(c);
    tg = qGreen(c);
    tr = qRed(c);
  }
  r2 = (tr-sr)*dx+sr;
  g2 = (tg-sg)*dx+sg;
  b2 = (tb-sb)*dx+sb;

  *r = (int)((r2-r1)*dy+r1);
  *g = (int)((g2-g1)*dy+g1);
  *b = (int)((b2-b1)*dy+b1);
}

/* ------------------------------------------------------------------------*/

void Renderer::drawMarkers()
{
  Location *l;
  double lon, lat;
  double s_x, s_y, s_z;
  double loc_x, loc_y, loc_z;
  int screen_x, screen_y;
  double m11, m12, m13,
         m21, m22, m23,
         m31, m32, m33;
  double visible_angle;

  int i, num;
  Location **visible_locations;

  visible_locations = new Location*[markerlist->count()];
  ASSERT(visible_locations != NULL);

  // Matrix M of renderFrame, but transposed
  m11 = cos(rot)*cos(view_long)-sin(view_lat)*sin(view_long)*sin(rot);
  m12 = sin(rot)*cos(view_lat);
  m13 = -sin(view_long)*cos(rot)-sin(view_lat)*cos(view_long)*sin(rot);
  m21 = -sin(rot)*cos(view_long)-sin(view_lat)*sin(view_long)*cos(rot);
  m22 = cos(view_lat)*cos(rot);
  m23 = sin(view_long)*sin(rot)-sin(view_lat)*cos(view_long)*cos(rot);
  m31 = cos(view_lat)*sin(view_long);
  m32 = sin(view_lat);
  m33 = cos(view_lat)*cos(view_long);

  visible_angle = radius/center_dist;
  
  for(i=0, l = markerlist->first(); l != NULL; l = markerlist->next())
  {
    lon = l->getLongitude()*PI/180.;
    lat = l->getLatitude()*PI/180.;

    s_x = radius*cos(lat)*sin(lon);
    s_y = radius*sin(lat);
    s_z = radius*cos(lat)*cos(lon);
    loc_x =  m11*s_x + m12*s_y + m13*s_z;
    loc_y =  m21*s_x + m22*s_y + m23*s_z;
    loc_z =  m31*s_x + m32*s_y + m33*s_z;
    
    l->cos_angle = loc_z/radius;
    
    if(l->cos_angle < visible_angle)
      // location lies on the other side
      continue;

    loc_y = -loc_y;
    loc_z = center_dist - loc_z;
    screen_x = (int)(loc_x*proj_dist/loc_z);
    screen_y = (int)(loc_y*proj_dist/loc_z);
    screen_x += renderedImage->width()/2;
    screen_y += renderedImage->height()/2;

    if((screen_x < 0) || (screen_x >= renderedImage->width()))
      // location out of bounds
      continue;
    if((screen_y <0) || (screen_y >= renderedImage->height()))
      continue;

    l->x = screen_x + shift_x;
    l->y = screen_y + shift_y;

    visible_locations[i] = l;
    i++;
  }

  num = i;

  // sort the markers according to depth
  qsort(visible_locations, num, sizeof(Location *),
        Renderer::compareLocations);

  for(i=0; i<num; i++)
    paintMarker(visible_locations[i]->x,
                visible_locations[i]->y, visible_locations[i]);

  delete[] visible_locations;
}

/* ------------------------------------------------------------------------*/

int Renderer::compareLocations(const void *l1, const void *l2)
{
  double c1, c2;
  
  c1 = (*((Location **)l1))->cos_angle;
  c2 = (*((Location **)l2))->cos_angle;
  
  if(c1 > c2)
    return 1;
  if(c1 < c2)
    return -1;
  return 0;
}

/* ------------------------------------------------------------------------*/

void Renderer::drawGrid()
{
  double lon, lat;
  double s_x, s_y, s_z;
  double loc_x, loc_y, loc_z;
  int screen_x, screen_y;
  double m11, m12, m13,
         m21, m22, m23,
         m31, m32, m33;
  double visible_angle;
  double cos_angle;
  double light_angle;
  double temp;
  
  unsigned int *p;
  unsigned int pixel;
  int r, g, b;
  
  // Matrix M of renderFrame, but transposed
  m11 = cos(rot)*cos(view_long)-sin(view_lat)*sin(view_long)*sin(rot);
  m12 = sin(rot)*cos(view_lat);
  m13 = -sin(view_long)*cos(rot)-sin(view_lat)*cos(view_long)*sin(rot);
  m21 = -sin(rot)*cos(view_long)-sin(view_lat)*sin(view_long)*cos(rot);
  m22 = cos(view_lat)*cos(rot);
  m23 = sin(view_long)*sin(rot)-sin(view_lat)*cos(view_long)*cos(rot);
  m31 = cos(view_lat)*sin(view_long);
  m32 = sin(view_lat);
  m33 = cos(view_lat)*cos(view_long);

  visible_angle = radius/center_dist;

  temp = PI/2.0 - d_gridline;
  
  for(lat=-temp; lat<=temp+0.01; lat+=d_gridline)
  {
    s_y = radius*sin(lat);

    for(lon=-PI; lon<PI; lon+=d_griddot)
    {
      s_x = radius*cos(lat)*sin(lon);
      s_z = radius*cos(lat)*cos(lon);
      loc_x =  m11*s_x + m12*s_y + m13*s_z;
      loc_y =  m21*s_x + m22*s_y + m23*s_z;
      loc_z =  m31*s_x + m32*s_y + m33*s_z;
    
      cos_angle = loc_z/radius;
      light_angle = (light_x*s_x + light_y*s_y + light_z*s_z)/radius;
    
      if(cos_angle < visible_angle)
      // location lies on the other side
      continue;
      
      loc_z = center_dist - loc_z;
      screen_x = (int)(loc_x*proj_dist/loc_z);
      screen_y = (int)(-loc_y*proj_dist/loc_z);
      screen_x += renderedImage->width()/2 + shift_x;
      screen_y += renderedImage->height()/2 + shift_y;

      if((screen_x < 0) || (screen_x >= renderedImage->width()))
        // location out of bounds
        continue;
      if((screen_y <0) || (screen_y >= renderedImage->height()))
        continue;

      p = (unsigned int *)renderedImage->scanLine(screen_y) + screen_x;

      // Set pixel in image
      if(gridtype == NICE_GRID)
      {
        pixel = getPixelColor(lon, -lat, light_angle);
        r = qRed(pixel)*3;
        g = qGreen(pixel)*3;
        b = qBlue(pixel)*3;

        r = (r>255) ? 255 : r;
        g = (g>255) ? 255 : g;
        b = (b>255) ? 255 : b;
      
        *p = qRgb(r, g, b);
      }
      else
        *p = qRgb(255, 255, 255);
    }
  }

  for(lon=-PI; lon<PI; lon+=d_gridline)
  {
    for(lat=-temp; lat<=temp; lat+=d_griddot)
    {
      s_x = radius*cos(lat)*sin(lon);
      s_y = radius*sin(lat);
      s_z = radius*cos(lat)*cos(lon);
      loc_x =  m11*s_x + m12*s_y + m13*s_z;
      loc_y =  m21*s_x + m22*s_y + m23*s_z;
      loc_z =  m31*s_x + m32*s_y + m33*s_z;
    
      cos_angle = loc_z/radius;
      light_angle = (light_x*s_x + light_y*s_y + light_z*s_z)/radius;
    
      if(cos_angle < visible_angle)
      // location lies on the other side
      continue;
      
      loc_z = center_dist - loc_z;
      screen_x = (int)(loc_x*proj_dist/loc_z);
      screen_y = (int)(-loc_y*proj_dist/loc_z);
      screen_x += renderedImage->width()/2 + shift_x;
      screen_y += renderedImage->height()/2 + shift_y;

      if((screen_x < 0) || (screen_x >= renderedImage->width()))
        // location out of bounds
        continue;
      if((screen_y <0) || (screen_y >= renderedImage->height()))
        continue;

      p = (unsigned int *)renderedImage->scanLine(screen_y) + screen_x;

      // Set pixel in image
      if(gridtype == NICE_GRID)
      {
        pixel = getPixelColor(lon, -lat, light_angle);
        r = qRed(pixel)*3;
        g = qGreen(pixel)*3;
        b = qBlue(pixel)*3;

        r = (r>255) ? 255 : r;
        g = (g>255) ? 255 : g;
        b = (b>255) ? 255 : b;
      
        *p = qRgb(r, g, b);
      }
      else
        *p = qRgb(255, 255, 255);
    }
  }
}

/* ------------------------------------------------------------------------*/

QImage *Renderer::getImage()
{
  QImage *clonedImage = NULL;

  clonedImage = new QImage(*renderedImage);
  ASSERT(clonedImage != NULL);
  return clonedImage;
}

/* ------------------------------------------------------------------------*/

void Renderer::paintMarker(int x, int y, Location *l)
{
  QPainter p;
  int wx, wy;
  unsigned int *src, *dest;
  unsigned int markercolor;
  
  QFontMetrics fm(QFont("helvetica", 12, QFont::Bold));
  
  QRect br = fm.boundingRect(l->getName());
  QPixmap pm(markerpixmap->width()+2 + br.width()+2,
             MAX(markerpixmap->height(), br.height()+2));
  
  p.begin(&pm);
  p.setFont(QFont("helvetica", 12, QFont::Bold));
  p.fillRect(0,0,pm.width(), pm.height(), QColor(0, 0, 255));

  // draw a black seam around the text
#if QT_VERSION >= 200  
  p.setPen(Qt::black);
#else
  p.setPen(black);
#endif
  wx = -br.x()+markerpixmap->width()+2;
  wy = -br.y();
  p.drawText(wx, wy+1, l->getName());
  p.drawText(wx+1, wy, l->getName());
  p.drawText(wx+1, wy+2, l->getName());
  p.drawText(wx+2, wy+1, l->getName());

  markercolor = l->getColor().rgb();
  
  p.setPen(QColor(255, 0, 0));
  p.drawText(wx+1, wy+1, l->getName());
  
  p.drawPixmap(0, (pm.height() - markerpixmap->height())/2,
               *markerpixmap);
  
  p.end();
  
  QImage labelimage = pm.convertToImage();
  if(labelimage.depth() != 32)
    labelimage = labelimage.convertDepth(32);

  x -= markerpixmap->width()/2;
  y -= labelimage.height()/2;

  QRect screenrect(0, 0, renderedImage->width(), renderedImage->height());
  QRect labelrect(x, y, pm.width(), pm.height());
  QRect visiblerect = screenrect.intersect(labelrect);

  if(visiblerect.isEmpty())
    // the label is not visible
    return;
  
  for(wy=0 ; wy<visiblerect.height(); wy++)
  {
    dest = (unsigned int *)renderedImage->scanLine(visiblerect.y()+wy) +
           visiblerect.x();
    src = (unsigned int *)labelimage.scanLine(visiblerect.y()-y+wy) +
          (visiblerect.x()-x);
    
    for(wx=0 ; wx<visiblerect.width(); wx++)
    {
      switch(*src++)
      {
        case 0x00000000:
          *dest++ = 0;
          break;

#if QT_VERSION >= 200          
        case 0x00FF0000:
#else
        case 0x000000FF:
#endif
          *dest++ = markercolor;
          break;

        default:
          dest++;
          break;
      }
    }
  }  
}

/* ------------------------------------------------------------------------*/

void Renderer::drawLabel()
{  
  QDateTime    dt;
  QString      labelstring;
  QPainter     p;
  unsigned int transparentcolor = qRgb(255, 255, 0);
  unsigned int pixel;
  int          x, y;
  int          wx, wy;
  unsigned int *src, *dest;
  double vlon, vlat;
  double slon, slat;
  struct tm    *tm;
  
  dt.setTime_t(time_to_render);
  tm = localtime(&time_to_render);
  
  vlon = view_long*180./PI;
  vlat = view_lat*180./PI;
  slon = sun_long*180./PI;
  slat = sun_lat*180./PI;
  
  labelstring.sprintf("%s, %s %d. %d, %d:%02d %s\n"
                      "View pos %2.2f° %c %2.2f° %c\n"
                      "Sun pos %2.2f° %c %2.2f° %c",
                      dt.date().dayName(dt.date().dayOfWeek()),
                      dt.date().monthName(dt.date().month()),
                      dt.date().day(), dt.date().year(),
                      dt.time().hour(), dt.time().minute(),
                      tzname[tm->tm_isdst],
                      fabs(vlat), (vlat < 0.) ? 'S' : 'N',
                      fabs(vlon), (vlon < 0.) ? 'W' : 'E',
                      fabs(slat), (slat < 0.) ? 'S' : 'N',
                      fabs(slon), (slon < 0.) ? 'W' : 'E');

  QFontMetrics fm(QFont("helvetica", 12, QFont::Bold));
  
#if QT_VERSION >= 200
  QRect br = fm.boundingRect(0, 0, 0, 0, Qt::AlignLeft|Qt::AlignTop,
                             labelstring);
#else
  QRect br = fm.boundingRect(0, 0, 0, 0, AlignLeft|AlignTop, labelstring);
#endif
  QPixmap pm(br.width()+10, br.height()+10);
  
  p.begin(&pm);
  p.setFont(QFont("helvetica", 12, QFont::Bold));
  p.fillRect(0,0,pm.width(), pm.height(), QColor(255, 255, 0));

#if QT_VERSION >= 200
  p.setPen(Qt::white);
  p.drawText(5, 5, br.width(), br.height(), Qt::AlignLeft|Qt::AlignTop,
             labelstring);
#else
  p.setPen(white);
  p.drawText(5, 5, br.width(), br.height(), AlignLeft|AlignTop, labelstring);
#endif
  p.end();
  
  QImage labelimage = pm.convertToImage();
  if(labelimage.depth() != 32)
    labelimage = labelimage.convertDepth(32);

  if(label_x > 0)
    x = label_x;
  else
    x = renderedImage->width()-labelimage.width()+label_x;
  if(label_y > 0)
    y = label_y;
  else
    y = renderedImage->height()-labelimage.height()+label_y;
  
  for(wy=0 ; wy<labelimage.height(); wy++)
  {
    dest = (unsigned int *)renderedImage->scanLine(y+wy) + x;
    src = (unsigned int *)labelimage.scanLine(wy);
    
    for(wx=0 ; wx<labelimage.width(); wx++)
    {      
      if(*src != transparentcolor)
        *dest++ = *src++;
      else
      {
        pixel = *dest;
        src++;
        *dest++ = qRgb(qRed(pixel)/2, qGreen(pixel)/2, qBlue(pixel)/2);
      }      
    }
  }
} 

/* ------------------------------------------------------------------------*/

void Renderer::drawStars()
{
  int x, y, brightness;
  unsigned int *p;
  int numstars;

  numstars = (int)(renderedImage->width() * renderedImage->height() *
                   star_frequency);
  
  for(int i=0; i<numstars; i++)
  {
    x = rand() % renderedImage->width();
    y = rand() % renderedImage->height();
    brightness = 150+ (rand() % 106);
    p = (unsigned int *)renderedImage->scanLine(y);
    p += x;
    *p = qRgb(brightness, brightness, brightness);
  }
}

/* ------------------------------------------------------------------------*/