/* 
   EarthView
   Copyright Jean-Pierre Demailly <demailly@ujf-grenoble.fr>
   Released under the LGPL

   GenEarth displays a color map of the Earth, using info from
   topographic/bathymetric or population data. It can also convert or
   compress these data in various forms.

   This program is published without any warranty
*/

#include <unistd.h>
#include <stdio.h>
#include <stdlib.h>
#include <ctype.h>
#include <string.h>
#include <math.h>
#include <sys/stat.h>
#include <zlib.h>
#include <bzlib.h>
#include <X11/X.h>
#include <X11/Xlib.h>
#include <X11/Xatom.h>
#include <X11/xpm.h>
#include <X11/keysym.h>
#include <X11/keysymdef.h>

typedef struct {
int data, w_data, h_data;
char *name;
double lon1, lon2, lat1, lat2;
} record;

Display * dpy = NULL;
Visual * visual;
Colormap cmap;
GC gc;
Atom wm_delete_window, wm_protocols;
int scr;
int depth;
Window root, win;
Pixmap pixmap;
Pixel black, white, *pix=NULL;

/* input modes */
enum {UNDEF=-1, POP=0, ALT, AZ2, PZ2, CPS, GTOPO30, MIXED };

/* output modes */
enum {RAW=0, PACK, STAT, PPM, XWIN };

int verbose=0;
int input=UNDEF, output=XWIN;
int header_data=0, header_out=1, aspect=1, listheader=0, num_scheme=1;
int w_data, h_data, im_data, jm_data, xm_data, ym_data;
int width, height;
int x_orig, y_orig, rectangle=0;
int new_colors=1;
int num_cps;

double londata1, londata2, latdata1, latdata2;
double lat1, lat2, lon1, lon2;
double lonbak1, lonbak2, latbak1=120.0, latbak2;
double gamma_f, speed = 2.0;
char **buf_out=NULL;
record *cps=NULL;

/* Default Digital Elevation Model file */
char *default_alt_file;
char *default_pop_file;

char *data_dir = NULL;
char *data_file = NULL;
char *cps_file, *cps_dir;
FILE *fd = NULL;
gzFile *gzd = NULL;

/* Predefined color schemes */
/* You can of course add as many as you like... */

char *color_scheme_alt[] = {
/* Topography color schemes */
 "-11000,0x000080,0,0x80C0FF,1,0x008000,200,0x00FF00,400,0xBBBB00,1000,0xFFFF00,2000,0xFF9000,3000,0xAA6000,3750,0x994000,4200,0x882000,4500,0x771000,5000,0x660000,8000,0x000000", 
 "-11000,0x000080,0,0x80C0FF,1,0x00AA00,200,0x00FF00,400,0xFFFF00,800,0xFF9000,1200,0xAE00AE,2000,0xA0A0A0,8000,0xFFFFFF",
 "-11000,0x000080,0,0x80C0FF,1,0x008000,200,0x00FF00,400,0xFFFF00,1000,0xFF9000,2000,0xAE20AE,3000,0xAE00AE,5000,0x61009E,9000,0x230046",
 "-11000,0x000080,0,0x80C0FF,1,0x00AA00,200,0x00FF00,400,0xFFFF00,800,0xFF9000,1200,0xD000D0,2000,0xB0A0B0,2500,0xA080A0,4000,0x808080,8000,0x000000"
};

char *color_scheme_pop[] = {
/* Population color schemes */
"-12000,0x80C0FF,-1,0x80C0FF,0,0xFFD7A6,1280,0xFF96AE,2560,0xCFCFCF,3840,0x9E9E9E,5120,0xAE69AE,6400,0xAE20AE,7680,0xAE00AE,8960,0x61009E,10240,0x230046",
 "-12000,0x80C0FF,-1,0x80C0FF,0,0xFFFF7F,2560,0x7FFF7F,5120,0x7F7F7F,7680,0xFF007F,10240,0xFF0000"
};
char *range = NULL;
int max_scheme_alt, max_scheme_pop;

int endian_data, endian_out, min_alt, max_alt;
double invcosphi = 1.0, min_den, max_den;
short int ocean_bit=0, in_ocean=0;

/* Red, Green, Blue values */
unsigned char *red=NULL, *green=NULL, *blue=NULL;

/* Deta files for GTOPO30 */
char * gtopo30_blocks[] = {
  "W180N90.DEM",
  "W140N90.DEM",
  "W100N90.DEM",
  "W060N90.DEM",
  "W020N90.DEM",
  "E020N90.DEM",
  "E060N90.DEM",
  "E100N90.DEM",
  "E140N90.DEM",
  "W180N40.DEM",
  "W140N40.DEM",
  "W100N40.DEM",
  "W060N40.DEM",
  "W020N40.DEM",
  "E020N40.DEM",
  "E060N40.DEM",
  "E100N40.DEM",
  "E140N40.DEM",
  "W180S10.DEM",
  "W140S10.DEM",
  "W100S10.DEM",
  "W060S10.DEM",
  "W020S10.DEM",
  "E020S10.DEM",
  "E060S10.DEM",
  "E100S10.DEM",
  "E140S10.DEM",
  "W180S60.DEM",
  "W120S60.DEM",
  "W060S60.DEM",
  "W000S60.DEM",
  "E060S60.DEM",
  "E120S60.DEM"
};
FILE *gfd[33];

void usage()
{ 
    printf("\n"
      "GenEarth displays a color map of the Earth, using info from\n"
      "topographic/bathymetric or population data. It can also convert or\n"
      "compress these data in various forms.\n\n"
      "Usage: earthview [-help] [-verbose] [-population]\n"
      "         [-nointerp] [-noaspect] [-listheader] [-oceanbit]\n"
      "         [-w width] [-h height] [-wdata width] [-hdata height]\n"
      "         [-lon val1 val2] [-lat val1 val2]\n"
      "         [-londata val1 val2] [-latdata val1 val2]\n"
      "         [-scheme num] [-range i1,color1,i2,color2...]\n"
      "         [-data {pop|alt|az2|pz2|gtopo30|gtopo30+}]\n"
      "         [-endian {data+|data-|out+|out-}]\n"
      "         [-header {data+|data-|out+|out-}]\n"
      "         [-output {xwin|raw|ppm|stat}]\n"
      "         [-pack num1 num2] [-gamma value(0.0...2.0)]\n"
      "         [-datadir directory] [-file filename] [> outfile.xxx]\n"
      "There are %d predefined color schemes (-scheme 1,..,-scheme %d)\n"
      "for topographic data, and %d schemes for the density of population.\n"
      "Color specifications should be in the RGB hexadecimal form 0xUVWXYZ\n"
      "The color range is meant to specify density colors for\n"
      "i=-12000,...i=+12000, indices i1,i2,... should therefore be in\n"
      "that interval.\n"
      "GenEarth automatically interpolates the other color values.\n\n",
        max_scheme_alt, max_scheme_alt, max_scheme_pop);
    exit(0);
}

int nfread(FILE * fd)
{
unsigned char c;
int n;
    c = fgetc(fd);
    if (feof(fd)) return -1;
    n = c;
    c = fgetc(fd);
    if (feof(fd)) return -1;
    n = (n<<8) | c;
    c = fgetc(fd);
    if (feof(fd)) return -1;
    n = (n<<8) | c;
    c = fgetc(fd);
    if (feof(fd)) return -1;
    n = (n<<8) | c;
    return n;
}

int ngzread(gzFile *gzd)
{
unsigned char c;
int n;
    c = gzgetc(gzd);
    if (gzeof(gzd)) return -2;
    n = c;
    c = gzgetc(gzd);
    if (gzeof(gzd)) return -2;
    n = (n<<8) | c;
    c = gzgetc(gzd);
    if (gzeof(gzd)) return -2;
    n = (n<<8) | c;
    c = gzgetc(gzd);
    if (gzeof(gzd)) return -2;
    n = (n<<8) | c;
    return n;
}

void nprintf(FILE * fo, int n)
{
    fprintf(fo, "%c", (n>>24) & 0xFF);
    fprintf(fo, "%c", (n>>16) & 0xFF);
    fprintf(fo, "%c", (n>>8)  & 0xFF);
    fprintf(fo, "%c",    n    & 0xFF);
}

#define NOVALUE -12000
int read_input(int i)
{
    double den;
    int n;
    int ip, j, k, l, u, v, x, y, xp, yp, s1, s2;
    unsigned char c1, c2, c3;
    char *buf_in;
    unsigned int destlen;

    if ((input|1)==PZ2) {
      x = i%w_data;
      y = i/w_data;
      u = x/xm_data;
      v = y/ym_data;
      xp = x%xm_data;
      yp = y%ym_data;
      k = u+v*im_data;
      /* allocate new output buffer if this is not yet done */
      if (!buf_out[k]) {
        /* before that, free all buffers that are no longer needed */
        if (v>0 && buf_out[k-im_data]) {
          for (l=0; l<im_data; l++) {
            ip = l+(v-1)*im_data;
            if (buf_out[ip]) {
              free(buf_out[ip]);
              buf_out[ip] = NULL;
            }
          }
        }
        if (verbose)
          fprintf(stderr, "\nReading sector [%d,%d]\n", u, v);
        fseek(fd, 4*k+header_data, SEEK_SET);
        s1 = nfread(fd);
        s2 = nfread(fd);
        buf_in = malloc((s2-s1)*sizeof(char));
	l = 2+(input==PZ2);
        buf_out[k] = malloc(l*xm_data*ym_data*sizeof(char));
        fseek(fd, header_data+4*(im_data*jm_data+1)+s1, SEEK_SET);
        for (l=0; l<s2-s1; l++)
          buf_in[l] = fgetc(fd);
        destlen = (unsigned int)(l*xm_data*ym_data);
        if (buf_out[k])
          l = BZ2_bzBuffToBuffDecompress(buf_out[k], &destlen,
                                         buf_in, (unsigned int) s2-s1, 0, 0);
        else
          l = BZ_OUTBUFF_FULL;
        free(buf_in);
        if (l==BZ_OUTBUFF_FULL) {
          fprintf(stderr, "Memory buffer[%d,%d] full !!\n", u, v);
          exit(-1);
        }
      }
      if (input==AZ2) {
        l = 2*(xp + yp*xm_data);
	/* I've decided that AZ2 is always smallendian !! */
        c1 = buf_out[k][l]; 
        c2 = buf_out[k][l+1]; 
        n = (c1<<8)|c2;
        /* high bit might be "ocean bit"; next bit is, in any case, sign bit */
        n = n&0x7FFF;
        if (n&0x4000) n |= 0xFFFF8000;
        if (n<min_alt) min_alt = n;
        if (n>max_alt) max_alt = n;
        in_ocean = c1>>7;
        return n;
      } else {
        l = 3*(xp + yp*xm_data);
	/* I've decided that PZ2 is always smallendian !! */
        c1 = buf_out[k][l]; 
        c2 = buf_out[k][l+1]; 
        c3 = buf_out[k][l+2];
        k = (int)((c1<<16) | (c2<<8) | c3);
        if (output<=PACK) return k;
	if (c1==255) return NOVALUE;
        den = ((double)k)*invcosphi;
        if (den<min_den) min_den = den;
        if (den>max_den) max_den = den;
        n = (int) (913.036*(log(den+40)-3.6888));
        if (n>10240) n=10240;
	return n;
      }
    }

    if (input>=GTOPO30) {
      x = i%43200;
      y = i/43200;
      if (y<18000) {
        u = x/4800;
        v = y/6000;
        j = u+9*v;
        xp = x-u*4800;
        yp = y-v*6000; 
        ip = xp+yp*4800;
      } else {
        u = x/7200;
        j = 27+u;
        xp = x-u*7200;
        yp = y-18000;
        ip = xp+yp*7200;
      }
      fseek(gfd[j], 2*ip, SEEK_SET);
      if (feof(gfd[j])) return NOVALUE;
      c1 = fgetc(gfd[j]);
      if (feof(gfd[j])) return NOVALUE;
      c2 = fgetc(gfd[j]);
      n = (c1<<8) | c2;
      if (n&0x8000) n |= 0xFFFF0000;
      if (n<min_alt) min_alt = n;
      if (n>max_alt) max_alt = n;
      /* This corrects a bug of GTOPO30 */
      if (y==21599) {
        n = read_input(x+(x==0)+(y-1)*43200);
	return n;
      } else
      if (x==0 && y>=18000) {
        n = read_input(1+y*43200);
        return n;
      }
      /* end of GTOPO30 bug fix */
      in_ocean = (n==-9999);
      if (in_ocean && input==MIXED) {
        /* We start reading in ETOPO2 for bathymetry */
        int n1, n2, n3, n4, xpp, ypp;
        double fx, fy;
        input = ALT;
        xp = x/4;
        fx = 0.25*(x%4);
        xpp = xp+1;
        if (xpp==10800) xpp=0;
        yp = y/4;
        fy = 0.25*(y%4);
        ypp = yp+1;
        if (ypp==5400) ypp=5399;
        n1 = read_input(xp+yp*10800);
        n2 = read_input(xpp+yp*10800);
        n3 = read_input(xp+ypp*10800);
        n4 = read_input(xpp+ypp*10800);
        n = (int)((1.0-fy)*((1.0-fx)*n1+fx*n2)+fy*((1.0-fx)*n3+fx*n4)-0.5);
        input = MIXED;
      }
      return n;
    }

    if (input==ALT) {
      gzseek(gzd, 2*i+header_data, SEEK_SET);
      if (gzeof(gzd)) return NOVALUE;
      if (endian_data) {
         c2 = gzgetc(gzd);
         if (gzeof(gzd)) return NOVALUE;
         c1 = gzgetc(gzd);
      } else {
         c1 = gzgetc(gzd);
         if (gzeof(gzd)) return NOVALUE;
         c2 = gzgetc(gzd);
      }
      n = (c1<<8) | c2;
      /* high bit might be "ocean bit"; next bit is, in any case, sign bit */
      n = n&0x7FFF;
      if (n&0x4000) n |= 0xFFFF8000;
      if (n<min_alt) min_alt = n;
      if (n>max_alt) max_alt = n;
      in_ocean = c1>>7;
      return n;
    }

    if (input==POP) {
      gzseek(gzd, 3*i+header_data, SEEK_SET);
      if (gzeof(gzd)) return NOVALUE;
      if (endian_data) {
         c3 = gzgetc(gzd);
         if (gzeof(gzd)) return NOVALUE;
         c2 = gzgetc(gzd);
         if (gzeof(gzd)) return NOVALUE;
         c1 = gzgetc(gzd);
      } else {
         c1 = gzgetc(gzd);
         if (gzeof(gzd)) return NOVALUE;
         c2 = gzgetc(gzd);
         if (gzeof(gzd)) return NOVALUE;
         c3 = gzgetc(gzd);
      }
      k = (int)((c1<<16) | (c2<<8) | c3);
      if (output<=PACK) return k;
      if (c1==255) return NOVALUE;
      den = ((double)k)*invcosphi;
      if (den<min_den) min_den = den;
      if (den>max_den) max_den = den;
      n = (int) (913.036*(log(den+40)-3.6888));
      if (n>10240) n=10240;
    }

    return n;
}

void do_output(int n, int x, int y)
{
    char c1, c2, c3;
    int l; 
    if (output==XWIN) {
      XSetForeground(dpy, gc, pix[n+12000]);
      XDrawPoint(dpy, win, gc, x, y); 
      XDrawPoint(dpy, pixmap, gc, x, y); 
    } else
    if (output==PPM) {
      if (n<-12000 || n>12000) {
        fprintf(stderr, "Strange value %d\n!!!", n);
        printf("%c%c%c", 0, 0, 0);
      } else {
        l = n+12000;
        printf("%c%c%c", red[l], green[l], blue[l]);
      }
    } else
    if (output==RAW && endian_out==0) {
      if (input==POP || input==PZ2) {
	if (n<0) n = 0xFF0000;
        c1 = (n & 0xFF0000)>>16;
        c2 = (n & 0xFF00)>>8;
        c3 = n & 0xFF;
        printf("%c%c%c", c1, c2, c3);
      } else {
        c1 = (n & 0xFF00)>>8;
        if (ocean_bit) {
	  if (in_ocean) 
	    c1 |= 0x80;
	  else
	    c1 &= 0x7F;
	}
        c2 = n & 0xFF;
        printf("%c%c", c1, c2);
      }
    } else
    if (output==RAW && endian_out==1) {
      if (input==POP || input==PZ2) {
	if (n<0) n = 0xFF0000;
        c1 = (n & 0xFF0000)>>16;
        c2 = (n & 0xFF00)>>8;
        c3 = n & 0xFF;
        printf("%c%c%c", c3, c2, c1);
      } else {
        c1 = (n & 0xFF00)>>8;
        if (ocean_bit) {
	  if (in_ocean) 
	    c1 |= 0x80;
	  else
	    c1 &= 0x7F;
	}
        c2 = n & 0xFF;
        printf("%c%c", c2, c1);
      }
    }
}

void init_colors(char *range)
{
    int i, j, k, boo;
    char *ptr1, *ptr2;
    unsigned int col;
    unsigned char c1, c2, c3, u, v;
    double f;

    if (!range) {
      if (input==POP || input==PZ2) {
	if (num_scheme>max_scheme_pop) num_scheme = max_scheme_pop;
        range = strdup(color_scheme_pop[num_scheme-1]);
      } else {
	if (num_scheme>max_scheme_alt) num_scheme = max_scheme_alt;
        range = strdup(color_scheme_alt[num_scheme-1]);
      }
    }
    if (verbose)
      fprintf(stderr, "Computing color range (scheme %s%d)...\n", 
        ((input==POP||input==PZ2)?"POP":"ALT"), num_scheme);

    ptr1 = ptr2 = range;
    if (!red)
      red = (unsigned char *)malloc(24002*sizeof(unsigned char));
    if (!green)
      green = (unsigned char *)malloc(24002*sizeof(unsigned char));
    if (!blue)
      blue = (unsigned char *)malloc(24002*sizeof(unsigned char));

    i = -1;
    boo = 1;
 begrange:
    while (*ptr2 && *ptr2!=',') ++ptr2;
    if (*ptr2==',') {
       *ptr2 = '\0';
       j = atoi(ptr1)+12000;
    } else
       j = -1;
    if (j<0 || j>24000) {
       fprintf(stderr, "Erroneous color range !!\n\n");
       exit(-1);
    }
    ++ptr2;
    ptr1 = ptr2;
    while (*ptr2 && *ptr2!=',') ++ptr2;
    if (!*ptr2) boo=0;
    *ptr2 = '\0';
    sscanf(ptr1, "%x", &col);
    red[j] = col>>16;
    green[j] = (col>>8)&255;
    blue[j] = col&255;
    if (i==-1) {
       c1 = red[j];
       c2 = green[j];
       c3 = blue[j];
    } else {
       c1 = red[i];
       c2 = green[i];
       c3 = blue[i];
    }
    for (k=i+1; k<j; k++) {
        f = ((double)(k-i))/((double)(j-i));
        red[k] = (unsigned char) ((1.0-f)*((double)c1)+f*((double)red[j])+0.5);
        green[k] = (unsigned char) ((1.0-f)*((double)c2)+f*((double)green[j])+0.5);
        blue[k] = (unsigned char) ((1.0-f)*((double)c3)+f*((double)blue[j])+0.5);
    } 
    if (boo) {
       ++ptr2;
       ptr1 = ptr2;
       i = j;
       goto begrange;
    }
    for (k=j+1; k<=24000; k++) {
        red[k] = red[j];
        green[k] = green[j];
        blue[k] = blue[j];
    }
    if (gamma_f<1.0) {
      for (k=0; k<=24000; k++) {
        red[k] = (unsigned char) (gamma_f*red[k]+0.5);
        green[k] = (unsigned char) (gamma_f*green[k]+0.5);
        blue[k] = (unsigned char) (gamma_f*blue[k]+0.5);

      }
    } else 
    if (gamma_f>1.0) {
      gamma_f = gamma_f - 1.0;
      for (k=0; k<=24000; k++) {
        red[k] = (unsigned char) ((1.0-gamma_f)*red[k]+gamma_f*255.0+0.5);
        green[k] = (unsigned char) ((1.0-gamma_f)*green[k]+gamma_f*255.0+0.5);
        blue[k] = (unsigned char) ((1.0-gamma_f)*blue[k]+gamma_f*255.0+0.5);
      }
      gamma_f = gamma_f + 1.0;
    }
    if (verbose)
      fprintf(stderr, "colors[-12000...12000] written to palette.\n\n");

    if (output==XWIN) {
       if (!pix)
          pix = (Pixel *)malloc(24002*sizeof(Pixel));
       if (verbose)
          fprintf(stderr, "Indexing colors for the X Window system\n");
       for (k=0; k<=24000; k++) {
          if (verbose && k%500==0) {
             fprintf(stderr, "[%d]", k);
             if ((k%500)%10==9) 
                fprintf(stderr, "\n");
             else
                fprintf(stderr, " ");
          }
          if (depth>=24)
            pix[k] = (red[k]<<16) | (green[k]<<8) | blue[k];
          else
          if (depth==16) {
             u = (blue[k]&248)>>3 | (green[k]&28)<<3;
             v = (green[k]&224)>>5 | (red[k]&248);
             pix[k] = (v<<8)|u;
          } else
          if (depth==15) {
             u = (blue[k]&248)>>3 | (green[k]&56)<<2;
             v = (green[k]&192)>>6 | (red[k]&248)>>1;
             pix[k] = (v<<8)|u;
          }
       }
       if (verbose)
          fprintf(stderr, "\n\n");
    }
    new_colors = 0;
}

void correct_aspect()
{
  double ratio, center_w, center_h, fw, fh;
  if (!aspect) return;
  ratio = sqrt(((double)width)*(lat2-lat1)/(((double)height)*(lon2-lon1))/
               cos(M_PI*(lat1+lat2)/360.0));
  center_w = 0.5*(lon1+lon2);
  center_h = 0.5*(lat1+lat2);
  fw = 0.5*(lon2-lon1)*ratio;
  fh = 0.5*(lat2-lat1)/ratio;
  lon1 = center_w-fw;
  lon2 = center_w+fw;
  lat1 = center_h-fh;
  lat2 = center_h+fh;

  while (lon1>180.0) {
    lon1 -= 360.00;
    lon2 -= 360.00;
  }
  while (lon1<-180.0) {
    lon1 += 360.00;
    lon2 += 360.00;
  }
}

Bool evpred(d, e, a)
Display *              d;
XEvent *               e;
XPointer               a;
{
        return (True);
}

void save_latlon()
{
    lonbak1 = lon1;
    lonbak2 = lon2;
    latbak1 = lat1;
    latbak2 = lat2;
}

int manage_event(XEvent event)
{
    double f, cor;
    int i;

    if (event.type == VisibilityNotify || 
        event.type == Expose) {
      XSetWindowBackgroundPixmap(dpy, win, pixmap);
      XClearWindow(dpy, win);
      XSync(dpy, True);   
      return 0;
    }
    if ((event.xclient.message_type == wm_protocols &&
         event.xclient.format == 32 &&
         event.xclient.data.l[0] == wm_delete_window))
      return 2;

    if (event.type == ConfigureNotify) {
      Window root;
      int i, j, x, y, w, h;
      unsigned int b, d;
      XGetGeometry(dpy, win, &root, &x, &y, &w, &h, &b, &d);
      if (w!=width && h!=height) {
	width = w;
        height = h;
        XFreePixmap(dpy, pixmap);
        pixmap = XCreatePixmap(dpy, win, width, height, depth);
	XSetForeground(dpy, gc, black);
        for (j=0; j<height; j++)
        for (i=0; i<width; i++)
	  XDrawPoint(dpy, pixmap, gc, i, j);
        correct_aspect();
        return 1;
      } else
	return 0;
    }
    if (event.type == ButtonPress) {
      x_orig = event.xbutton.x;
      y_orig = event.xbutton.y;
      XSetWindowBackgroundPixmap(dpy, win, pixmap);
      rectangle = 1;
      return 0;
    }
    if (event.type == ButtonRelease) {
      int x, y, w, h;
      double center_x, center_y, fw, fh;
      if (event.xbutton.x<0 || event.xbutton.x>width) return 0;
      if (event.xbutton.y<0 || event.xbutton.y>height) return 0;
      if (event.xbutton.x>x_orig) x=x_orig; else x=event.xbutton.x;
      if (event.xbutton.y>y_orig) y=y_orig; else y=event.xbutton.y;
      w = abs(event.xbutton.x-x_orig);
      h = abs(event.xbutton.y-y_orig);
      rectangle = 0;
      XClearWindow(dpy, win);
      center_x = lon1+(lon2-lon1)*((double)(x+w/2))/((double)width);
      center_y = lat2-(lat2-lat1)*((double)(y+h/2))/((double)height);
      if (w<=2 && h<=2) {
	if (event.xbutton.button!=3) return 0;
	/* click almost same point --> translate view rectangle */
        fw = 0.5*(lon2-lon1);
	fh = 0.5*(lat2-lat1);
      } else {
	/* click other point --> zoom scene */
        fw = 0.5*(lon2-lon1)*((double)w)/((double)width);
        fh = 0.5*(lat2-lat1)*((double)h)/((double)height);
      }
      save_latlon();
      lon1 = center_x-fw;
      lon2 = center_x+fw;
      lat1 = center_y-fh;
      lat2 = center_y+fh;
      correct_aspect();
      return 1;
    }
    if (event.type == MotionNotify && rectangle) {
      int x, y, w, h;
      XClearWindow(dpy, win);
      if (event.xbutton.x>x_orig) x=x_orig; else x=event.xbutton.x;
      if (event.xbutton.y>y_orig) y=y_orig; else y=event.xbutton.y;
      w = abs(event.xbutton.x-x_orig);
      h = abs(event.xbutton.y-y_orig);
      XSetForeground(dpy, gc, white);
      XDrawRectangle(dpy, win, gc, x, y, w, h);
      return 0;
    }
    if (event.type == KeyPress) {
      KeySym keysym;
      char buffer[1];
      XLookupString((XKeyEvent *) &event, buffer, 1, &keysym, NULL);
      if (keysym==XK_BackSpace || keysym==XK_Delete) {
	 if (latbak1>100.0) return 0;
         lon1 = lonbak1;
         lon2 = lonbak2;
         lat1 = latbak1;
         lat2 = latbak2;
	 return 1;
      }
      if (((i=(int)(keysym-XK_0))>=1 && i<=9) ||
          ((i=(int)(keysym-XK_KP_0))>=1 && i<=9)) {
	 if (i==num_scheme) return 0;
	 num_scheme = i;
	 if (verbose) {
           fprintf(stderr, "Changing to color scheme %d...\n", num_scheme);
	 }
	 free(range);
	 range = NULL;
	 new_colors = 1;
	 return 1;
      }
      if (keysym==XK_equal) {
         save_latlon();
         lon1 = londata1;
         lon2 = londata2;
         lat1 = latdata1;
         lat2 = latdata2;
	 return 1;
      }
      if (keysym==XK_Page_Up) {
	speed = speed*speed;
	return 0;
      }
      if (keysym==XK_Page_Down) {
	speed = sqrt(speed);
	return 0;
      }
      if (keysym==XK_c || keysym==XK_C) {
        if (cps && cps_file==default_alt_file) {
	  data_file = default_pop_file;
	  for (i=0; i<num_cps; i++) free(cps[i].name);
	  free(cps);
          if (cps_dir) free(cps_dir);
	  cps_dir = NULL;
          cps = NULL;
	  num_cps = 0;
	  input = UNDEF;
	  free(range);
	  range = NULL;
	  new_colors = 1;
	  return 1;
	} else
        if (cps && cps_file==default_pop_file) {
	  data_file = default_alt_file;
	  for (i=0; i<num_cps; i++) free(cps[i].name);
	  free(cps);
          if (cps_dir) free(cps_dir);
	  cps_dir = NULL;
	  num_cps = 0;
          cps = NULL;
	  input = UNDEF;
	  free(range);
	  range = NULL;
	  new_colors = 1;
	  return 1;
	} else
	  return 0;
      }
      if (keysym==XK_d || keysym==XK_D) {
	if (gamma_f<=1.2/1.1)
	  gamma_f = gamma_f/1.2;
	else
	  gamma_f = 2.0 - (2.0-gamma_f)*1.2;
	if (verbose)
	  fprintf(stderr, "New gamma factor %.3f\n", gamma_f);
	new_colors = 1;
	return 1;
      }
      if (keysym==XK_l || keysym==XK_L) {
	if (gamma_f<=1.0/1.1)
	  gamma_f = gamma_f*1.2;
	else
	  gamma_f = 2.0 - (2.0-gamma_f)/1.2;
	if (verbose)
	  fprintf(stderr, "New gamma factor %.3f\n", gamma_f);
	new_colors = 1;
	return 1;
      }
      if (keysym==XK_plus || keysym==XK_KP_Add) {
        save_latlon();
	cor = 0.5*(1.0-sqrt(1.0/speed));
        f = cor*(lon2-lon1);
	lon1 += f;
	lon2 -= f;
        f = cor*(lat2-lat1);
	lat1 += f;
	lat2 -= f;
	correct_aspect();
	return 1;
      }
      if (keysym==XK_minus || keysym==XK_KP_Subtract) {
        save_latlon();
	cor = 0.5*(sqrt(speed)-1);
        f = cor*(lon2-lon1);
	lon1 -= f;
	lon2 += f;
        f = cor*(lat2-lat1);
	lat1 -= f;
        lat2 += f;
	correct_aspect();
	return 1;
      }
      if (keysym==XK_Right) {
        save_latlon();
        f = (lon2-lon1)*speed/(1.0+speed);
	lon1 += f;
	lon2 += f;
	return 1;
      }
      if (keysym==XK_Left) {
        save_latlon();
        f = (lon2-lon1)*speed/(1.0+speed);
	lon1 -= f;
	lon2 -= f;
	correct_aspect();
	return 1;
      }
      if (keysym==XK_Up) {
        save_latlon();
        f = (lat2-lat1)*speed/(1.0+speed);
	lat1 += f;
	lat2 += f;
	correct_aspect();
	return 1;
      }
      if (keysym==XK_Down) {
        save_latlon();
        f = (lat2-lat1)*speed/(1.0+speed);
	lat1 -= f;
	lat2 -= f;
	correct_aspect();
	return 1;
      }
      if (keysym=='q' || keysym=='Q') return 2;
    }
    return 0;
}

int main(int argc, char **argv)
{
    double alon, blon, alat, blat, fcopy, fxp, fyp;
    int i, j, k, l, dil_x, dil_y, alt, n_alt, o_alt, 
      im_out, jm_out, xm_out, ym_out, 
      x, y, xp, yp, n_xp, n_yp,
      interp=1, reiterate, do_pop = 0;
    int *n_mem_alt = NULL, *o_mem_alt = NULL;
    XGCValues gcv;
    XEvent event;
    char title[256];
    char buf_read[256];

    endian_data = 0; /* smallendian is default for data */
    endian_out = 0;  /* smallendian is default for output */
    width = 0;
    height = 0;
    londata1 = -180.0;
    londata2 = +180.0;
    latdata1 = -90.0;
    latdata2 = +90.0;
    lon1 = -200.0;
    lon2 = -200.0;
    lat1 = 100.0;
    lat2 = 100.0;
    gamma_f = 1.0;
    max_scheme_alt = sizeof(color_scheme_alt)/sizeof(char *);
    max_scheme_pop = sizeof(color_scheme_pop)/sizeof(char *);
    w_data = -1;
    h_data = -1;

    for (i=1; i<argc; i++) {
       if (!strcasecmp(argv[i], "-help")) usage();
       else
       if (!strcasecmp(argv[i], "-verbose")) verbose = 1;
       else
       if (!strcasecmp(argv[i], "-nointerp")) interp = 0;
       else
       if (!strcasecmp(argv[i], "-noaspect")) aspect = 0;
       else
       if (!strcasecmp(argv[i], "-oceanbit")) ocean_bit = 1;
       else
       if (!strncasecmp(argv[i], "-population", 4)) do_pop = 1;
       else
       if (!strcasecmp(argv[i], "-listheader")) {
	 listheader = 1;
	 verbose = 1;
       }
       else
       if (i<argc-1 && !strcasecmp(argv[i], "-scheme")) {
         num_scheme = atoi(argv[++i]);
         if (num_scheme<=0) num_scheme = 1;
       }
       else
       if (i<argc-1 && !strcasecmp(argv[i], "-header")) {
         char *arg = argv[++i];
         if (!strcasecmp(arg, "data+")) header_data = 1;
         else
         if (!strcasecmp(arg, "data-")) header_data = 0;
         else
         if (!strcasecmp(arg, "out+")) header_out = 1;
         else
         if (!strcasecmp(arg, "out-")) header_out = 0;
         else
           usage();
       }
       else
       if (i<argc-1 && !strcasecmp(argv[i], "-endian")) {
         char *arg = argv[++i];
         if (!strcasecmp(arg, "data+")) endian_data = 1;
         else
         if (!strcasecmp(arg, "data-")) endian_data = 0;
         else
         if (!strcasecmp(arg, "out+")) endian_out = 1;
         else
         if (!strcasecmp(arg, "out-")) endian_out = 0;
         else
           usage();
       }
       else
       if (i<argc-1 && !strcasecmp(argv[i], "-data")) {
         char *arg = argv[++i];
         if (!strcasecmp(arg, "pop")) input = POP;
         else
         if (!strcasecmp(arg, "alt")) input = ALT;
         else
         if (!strcasecmp(arg, "az2")) input = AZ2;
         else
         if (!strcasecmp(arg, "pz2")) input = PZ2;
         else
         if (!strcasecmp(arg, "gtopo30")) input = GTOPO30;
         else
         if (!strcasecmp(arg, "gtopo30+")) input = MIXED;
         else
           usage();
       }
       else
       if (i<argc-1 && !strcasecmp(argv[i], "-output")) {
         char *arg = argv[++i];
         if (!strcasecmp(arg, "stat")) {
           output = STAT;
           verbose = 1;
         }
         else
         if (!strcasecmp(arg, "raw")) output = RAW;
         else
         if (!strcasecmp(arg, "ppm")) output = PPM;
         else
         if (!strcasecmp(arg, "xwin")) output = XWIN;
         else
           usage();
       }
       else
       if (i<argc-1 && !strcasecmp(argv[i], "-wdata")) {
         w_data = atoi(argv[++i]);
         if (w_data<100) exit(1);
       }
       else
       if (i<argc-1 && !strcasecmp(argv[i], "-hdata")) {
         h_data = atoi(argv[++i]);
         if (h_data<50) exit(1);
       }
       else
       if (i<argc-1 && !strncasecmp(argv[i], "-w", 2)) {
         width = atoi(argv[++i]);
         if (width<=10) width = 10;
       }
       else
       if (i<argc-1 && !strncasecmp(argv[i], "-h", 2)) {
         height = atoi(argv[++i]);
         if (height<=5) height = 5;
       }
       else
       if (i<argc-1 && !strcasecmp(argv[i], "-gamma")) {
         gamma_f = atof(argv[++i]);
         if (gamma_f<0) gamma_f = 0.0;
         if (gamma_f>2.0) gamma_f = 2.0;
       }
       else
       if (i<argc-1 && !strcasecmp(argv[i], "-file"))
         data_file = argv[++i];
       else
       if (i<argc-1 && !strcasecmp(argv[i], "-datadir"))
         data_dir = argv[++i];
       else
       if (i<argc-1 && !strcasecmp(argv[i], "-range")) 
         range = strdup(argv[++i]);
       else
       if (i<argc-2 && !strcasecmp(argv[i], "-lon")) {
         lon1 = atof(argv[++i]);
         lon2 = atof(argv[++i]);
       }
       else
       if (i<argc-2 && !strcasecmp(argv[i], "-lat")) {
         lat1 = atof(argv[++i]);
         lat2 = atof(argv[++i]);
       }
       else
       if (i<argc-2 && !strcasecmp(argv[i], "-londata")) {
         londata1 = atof(argv[++i]);
         londata2 = atof(argv[++i]);
       }
       else
       if (i<argc-2 && !strcasecmp(argv[i], "-lat")) {
         latdata1 = atof(argv[++i]);
         latdata2 = atof(argv[++i]);
       }
       else
       if (i<argc-2 && !strcasecmp(argv[i], "-pack")) {
         output = PACK;
         im_out = atoi(argv[++i]);
         jm_out = atoi(argv[++i]);
       }
       else
       if (argv[i][0] != '-')
	 data_file = argv[i];
       else
         usage();
    }

    if (data_dir==NULL)
      data_dir = DATADIR;

    default_alt_file = (char*)malloc((strlen(data_dir)+50)*sizeof(char));
    default_pop_file = (char*)malloc((strlen(data_dir)+50)*sizeof(char));
    sprintf(default_alt_file, "%s/%s", data_dir, "topography/topography.cps");
    sprintf(default_pop_file, "%s/%s", data_dir, "population/population.cps");

    if (data_file==NULL)
      data_file = do_pop? default_pop_file : default_alt_file;

 initialize:
    reiterate = 0;

    if (lon1==-200.0 && lon2==-200.0 && lat1==100.0 && lat2==100.0) {
      lon1 = londata1;
      lon2 = londata2;
      lat1 = latdata1;
      lat2 = latdata2;
    }

    if (input >= GTOPO30) {
       char *ptr1, *ptr2;
       w_data = 43200;
       h_data = 21600;
       ptr1 = strdup(data_file);
       ptr2 = index(ptr1, ':');
       if (ptr2) {
	  *ptr2 = '\0';
	  ++ptr2;
       }
       endian_data = 0;
       for (i=0; i<33; i++) {
          char fname[512];
          sprintf(fname, "%s/%s", ptr1, gtopo30_blocks[i]);
          gfd[i] = fopen(fname, "r");
          if (!gfd[i]) {
            fprintf(stderr, "Directory '%s' does not exist or\n"
                            "File '%s' cannot be read !!\n\n", ptr1, fname);
            exit(-1);
          }
       }
       if (input==MIXED) {
         if (!ptr2) {
           fprintf(stderr, "etopo2 data missing!!\n");
           exit(-1);
         }
         gzd = gzopen(ptr2, "r");
         if (!gzd) {
           fprintf(stderr, "File '%s' cannot be read !!\n\n", ptr2);
         exit(-1);
         }
	 header_data = 0;
         gzgets(gzd, buf_read, 255);
         if (!strncmp(buf_read, "@#!ALT", 6)) {
	   i = 0;
	   j = 1;
           do {
	     gzgets(gzd, buf_read, 255);
	     ++i;
	   }
	   while (i<=20 && !gzeof(gzd) && (j=strncmp(buf_read, "size:", 5)));
	   if (j) {
             fprintf(stderr, 
               "\nCouldn't find start of data section in '%s' !!\n", ptr2);
	     exit(-1);
	   } else {
	     header_data = gztell(gzd);
 	     if (verbose)
               fprintf(stderr, 
                 "\nFound start of data section in '%s'\n"
                 "after header of length %d ...\n", 
                 ptr2, header_data);
	   }
	 }
       }
       free(ptr1);
       goto processdata;
    }

    if (verbose)
      fprintf(stderr, "Opening file '%s' ...\n", data_file);
    gzd = gzopen(data_file, "r");
    if (!gzd) {
      fprintf(stderr, "File '%s' cannot be read !!\n\n", data_file);
      exit(-1);
    }
    gzgets(gzd, buf_read, 255);
    i = UNDEF;
    if (!strncmp(buf_read, "@#!ALT", 6)) i = ALT;
    else
    if (!strncmp(buf_read, "@#!POP", 6)) i = POP;
    else
    if (!strncmp(buf_read, "@#!AZ2", 6)) i = AZ2;
    else
    if (!strncmp(buf_read, "@#!PZ2", 6)) i = PZ2;
    else
    if (!strncmp(buf_read, "@#!CPS", 6)) i = CPS;
    if (i!=UNDEF) {
      header_data = 1;
      if (input==UNDEF) 
        input = i;
      else
      if (input!=i) {
        fprintf(stderr, 
          "The specified data type doesn't match the actual data.\n"
	  "Aborting!!\n");
        exit(-1);
      }
    }
    if (input==UNDEF) {
      fprintf(stderr, 
        "Cannot determine the data type: file %s has no legible header.\n"
	"You may try to set it manually with the option -data XXX\n"
	"Aborting!!\n", data_file);
      exit(-1);
    }
    if (input==AZ2 || input==PZ2) header_data = 1;
    if (input==CPS && !cps) {
      char data[10], name[220];
      char *ptr;
      struct stat buf;
      /* composite data */
      header_data = 0;
      cps_file = data_file;
      ptr = rindex(cps_file, '/');
      if (ptr) {
	i = ptr-cps_file+1;
	cps_dir = (char *)malloc((i+1)*sizeof(char));
	strncpy(cps_dir, cps_file, i);
	cps_dir[i] = '\0';
      } else
	cps_dir = NULL;
      num_cps = 0;
      while (!gzeof(gzd)) {
        gzgets(gzd, buf_read, 255);
	if (!strncmp(buf_read, "%end", 4)) break;
	if (*buf_read=='%' || isspace(*buf_read)) continue;
	cps = realloc(cps, (num_cps+1)*sizeof(record));
	i = sscanf(buf_read, "%s %d %d %s %lg %lg %lg %lg",
	  data, &cps[num_cps].w_data, &cps[num_cps].h_data, name, 
          &cps[num_cps].lon1, &cps[num_cps].lon2,
          &cps[num_cps].lat1, &cps[num_cps].lat2);
	if (i<8) {
	  fprintf(stderr, "Incorrect syntax in %s :\n%s\n",
	    data_file, buf_read);
	  continue;
	}
	i = UNDEF;
        if (!strcasecmp(data, "ALT")) i = ALT;
        else
        if (!strcasecmp(data, "POP")) i = POP;
        else
        if (!strcasecmp(data, "AZ2")) i = AZ2;
        else
        if (!strcasecmp(data, "PZ2")) i = PZ2;
        else
        if (!strcasecmp(data, "GTOPO30")) i = GTOPO30;
        else
        if (!strcasecmp(data, "GTOPO30+")) i = MIXED;
	if (i==UNDEF) {
	  fprintf(stderr, "Wrong data type definition in %s :\n%s\n",
	    data_file, buf_read);
	  continue;
	}
	if (*name=='/' || !cps_dir)
	  /* path is absolute path or current directory is that of .cps file */
	  cps[num_cps].name = strdup(name);
	else {
	  /* relative path */
	  cps[num_cps].name = 
            malloc((strlen(name)+strlen(cps_dir)+1)*sizeof(char));
	  sprintf(cps[num_cps].name, "%s%s", cps_dir, name);
	}
	j = stat(cps[num_cps].name, &buf);
	if (j!=0 || !S_ISREG(buf.st_mode)) {
	  if (verbose)
	    fprintf(stderr, "!!cancelled: %s\n", cps[num_cps].name);
	  free(cps[num_cps].name);
	  continue;
	}
	if (verbose)
	  fprintf(stderr, "%s", buf_read);
	cps[num_cps].data = i;
	++num_cps;
      }
      if (num_cps==0) {
	  fprintf(stderr, "Composite file %s contains no valid record.\n"
		          "Aborting!!\n", data_file);
	  exit(-1);
      }
      londata1 = cps[0].lon1;
      londata2 = cps[0].lon2;
      latdata1 = cps[0].lat1;
      latdata2 = cps[0].lat2;
      w_data = cps[0].w_data;
      h_data = cps[0].h_data;
      if (gzd) {
	gzclose(gzd);
	gzd = NULL;
      }
    }

    if (header_data) {
      if (verbose)
	fprintf(stderr, "\nReading header data from '%s'\n\n", data_file);
      if (verbose) fprintf(stderr, buf_read); /* print magic string */
      gzgets(gzd, buf_read, 255);  /* source (just for info) */
      if (verbose) fprintf(stderr, buf_read);
      gzgets(gzd, buf_read, 255);  /* reading lon/lat of data */
      i = sscanf(buf_read, "lon/lat: %lg %lg %lg %lg",
            &londata1, &londata2, &latdata1, &latdata2);
      if (i<4) {
        fprintf(stderr, 
          "Header of input data %s has incorrect lon/lat values:\n  %s\n",
          data_file, buf_read);
	exit(-1);
      }
      if (verbose)
        fprintf(stderr, "lon/lat: %.4f %.4f %.4f %.4f\n",
	  londata1, londata2, latdata1, latdata2);
      gzgets(gzd, buf_read, 255);  /* reading w_data h_data */
      i = sscanf(buf_read, "size: %d %d", &w_data, &h_data);
      if (i<2) {
        fprintf(stderr, 
          "Header of input data %s has incorrect size:\n  %s\n",
          data_file, buf_read);
	exit(-1);
      }
      if (verbose)
        fprintf(stderr, "size: %d %d\n", w_data, h_data);
      if (input == AZ2 || input == PZ2) {
        gzgets(gzd, buf_read, 255);  /* reading im_data jm_data */
        i = sscanf(buf_read, "pack: %d %d", &im_data, &jm_data);
        if (i<2) {
          fprintf(stderr, 
            "Header of input data %s has incorrect pack values:\n  %s\n",
            data_file, buf_read);
	  exit(-1);
	}
        if (verbose)
          fprintf(stderr, "pack: %d %d\n", im_data, jm_data);
	header_data = gztell(gzd);
        gzclose(gzd);
	gzd = NULL;
	fd = fopen(data_file, "r"); 
	xm_data = w_data/im_data;
	ym_data = h_data/jm_data;
        buf_out = (char **)realloc(buf_out, im_data*jm_data*sizeof(char *));
        memset(buf_out, 0, im_data*jm_data*sizeof(char *));
      } else
	header_data = gztell(gzd);
    }
    if (verbose && header_data)
      fprintf(stderr, "header length: %d\n\n", header_data);
    if (listheader) {
      if (fd) fclose(fd);
      if (gzd) gzclose(gzd);
      exit(0);
    }

 processdata:
    if (w_data==-1 || h_data==-1) {
      fprintf(stderr, 
         "Dimensions (w=%d,h=%d) of input incorrectly specified !!\n",
            w_data, h_data);
      exit(-1);
    }
    if (w_data<=10 || h_data<=5) {
      fprintf(stderr, "data_file file '%s' seems to be corrupted !!\n", 
         data_file);
      exit(-1);
    }

    if (londata1==londata2 || latdata1==latdata2) {
       fprintf(stderr, "Specified data region empty !!\n");
       exit(-1);
    }

    if (lon1==lon2 || lat1==lat2) {
       fprintf(stderr, "Specified output region empty !!\n");
       exit(-1);
    }

    if (londata1>londata2) {
       fcopy = londata1;
       londata1 = londata2;
       londata2 = fcopy;
    }
    if (latdata1>latdata2) {
       fcopy = latdata1;
       latdata1 = latdata2;
       latdata2 = fcopy;
    }

    if (lon1>lon2) {
       fcopy = lon1;
       lon1 = lon2;
       lon2 = fcopy;
    }
    if (lat1>lat2) {
       fcopy = lat1;
       lat1 = lat2;
       lat2 = fcopy;
    }

    if (lat1<latdata1) lat1 = latdata1;
    if (lat2>latdata2) lat2 = latdata2;

    while (lon1>180.0) {
      lon1 -= 360.00;
      lon2 -= 360.00;
    }

    while (lon1<-180.0) {
      lon1 += 360.00;
      lon2 += 360.00;
    }

    if (fabs(londata2-londata1-360.0)>1E-5) {
      if (lon1<londata1) lon1 = londata1;
      if (lon2>londata2) lon2 = londata2;
    }

    if (width==0 || height==0) {
      fcopy = (lon2-lon1)/(lat2-lat1);
      if (aspect) fcopy *= cos(M_PI*(lat1+lat2)/360.0);
      if (width==0 && height==0) {
        width = (int)(720 * sqrt(fcopy)+0.5);
        height = (int)(720 / sqrt(fcopy)+0.5);
      } else
      if (width!=0)
        height = (int) (width/fcopy+0.5);
      else
      if (height!=0)
        width = (int) (height*fcopy+0.5);
    }

 cpsinit:
    if (input==CPS && cps) {
      j = 0;
      for (i=0; i<num_cps; i++) {
        if (((lon2-lon1)/(londata2-londata1))*(double)cps[i].w_data
               <(double)width-0.5 ||
            ((lat2-lat1)/(latdata2-latdata1))*(double)cps[i].h_data
               <(double)height-0.5) {
	  j = 1;
	  break;
	}
      }
      if (j==0) i=num_cps-1;
      if (j>0 && i>0) --i;
      data_file = cps[i].name;
      input = cps[i].data;
      londata1 = cps[i].lon1;
      londata2 = cps[i].lon2;
      latdata1 = cps[i].lat1;
      latdata2 = cps[i].lat2;
      w_data = cps[i].w_data;
      h_data = cps[i].h_data;
      goto initialize;
    }

    dil_x = 0;
    dil_y = 0;    
    alon = (lon2-lon1)*((double)w_data)/((double)width*(londata2-londata1));
    blon = (lon1-londata1)*((double)w_data)/(londata2-londata1) + 0.0001;
    alat = (lat2-lat1)*((double)h_data)/((double)height*(latdata2-latdata1));
    blat = (lat1-latdata1)*((double)h_data)/(latdata2-latdata1) - 0.0001;

    if (interp) {
      if (alon<0.99999) dil_x = 1;
      if (alat<0.99999) dil_y = 1;
    }

    if (output==XWIN && dpy==NULL) {
       dpy = XOpenDisplay(NULL);
       scr = DefaultScreen(dpy);
       visual = DefaultVisual(dpy, scr);
       root = DefaultRootWindow(dpy);
       black = BlackPixel(dpy, scr);
       white = WhitePixel(dpy, scr);
       gcv.foreground = black;
       gc = XCreateGC(dpy, root, GCForeground, &gcv);
       depth = DefaultDepth(dpy, scr);
       if (depth<=8) {
          fprintf(stderr, 
                 "Can't use X Window output in depth <= 8. Aborting!!\n");
          exit(1);
       }
       else
          cmap = DefaultColormap(dpy, scr);
       win = XCreateSimpleWindow(dpy, root,
                 10, 10, width, height, 0, white, black);
       pixmap = XCreatePixmap(dpy, win, width, height, depth);
       for (j=0; j<height; j++)
       for (i=0; i<width; i++)
	 XDrawPoint(dpy, pixmap, gc, i, j);
       if (depth == 16 && visual->green_mask==992) depth = 15;
       wm_protocols = XInternAtom(dpy, "WM_PROTOCOLS", False);
       wm_delete_window = XInternAtom(dpy, "WM_DELETE_WINDOW", False);
       XSetWMProtocols(dpy, win, &wm_delete_window, 1);
       XSelectInput(dpy, win, 
         KeyPressMask|VisibilityChangeMask|ExposureMask|PointerMotionMask|
	 StructureNotifyMask|ButtonPressMask|ButtonReleaseMask);
       XMapWindow(dpy, win);
       XFlush(dpy);
       usleep(10000);
    }

    if (output==XWIN && dpy!=NULL) {
      XResizeWindow(dpy, win, width, height);
      XSync(dpy, True);
    }
    if (output==XWIN) {
      sprintf(title, "%s (lon %.2f %.2f / lat %.2f %.2f) scheme %d [%dx%d%s]", 
        ((input==POP||input==PZ2)?"Population":"Topography"), 
	lon1, lon2, lat1, lat2, 
        num_scheme, w_data, h_data, ((dil_x || dil_y)?",i":""));
      XStoreName(dpy, win, title);
    }

    if (new_colors)
      init_colors(range);

    if (verbose) {
      if (input==POP||input==PZ2)
        fprintf(stderr, "Reading population densities from '%s':\n", 
          data_file);
      else
        fprintf(stderr, "Reading elevation values from '%s':\n", data_file);
      fprintf(stderr, "  scanning %dx%d grid at resolution %dx%d (%s)\n",
        w_data, h_data, width, height,
        ((output==STAT)? "statistics only":"image output"));
      if (output!=STAT)
          fprintf(stderr, "Output image %dx%d : %.3f <= lon <= %.3f, "
                       "%.3f <= lat <= %.3f\n", width, height, 
                       lon1, lon2, lat1, lat2);
    }

    max_alt = -12000;
    min_alt = 10000;
    min_den = 1E9;
    max_den = -12000.0;

    if (output==PACK) {
      int total, mode;
      FILE *fo;
      char *source, *dest;
      xm_out = width/im_out;
      ym_out = height/jm_out;
      if (xm_out*im_out!=width || ym_out*jm_out!=height) {
        fprintf(stderr, 
          "Output size %dx%d not compatible with (%d,%d) packing\n"
	  "Aborting!!\n", width, height, im_out, jm_out);
	exit(-1);
      }
      /* create all necessary buffers for the .bz2 chunks */
      mode = (input==POP||input==PZ2);
      source = (unsigned char *)
        malloc((2+mode)*xm_out*ym_out*sizeof(unsigned char));
      dest = (unsigned char *)
        malloc(((2+mode)*xm_out*ym_out+100)*sizeof(unsigned char));
      if (!source || !dest) {
	fprintf(stderr, "Memory allocation failed!!\n");
	exit(-1);
      }
      /* create output file */
      sprintf(buf_read, "%dx%d-%dx%d.%s", 
        width, height, im_out, jm_out, 
        ((mode)? "pz2":"az2"));
      fo = fopen(buf_read, "w");
      if (!fo) {
        fprintf(stderr, "File '%s' cannot be created\n", buf_read);
        exit(-1);
      }
      fprintf(stderr, "Writing file '%s' ...\n", buf_read);
      /* 8 char identification string for new formats "packed BZ2" */
      if (mode)
        sprintf(buf_read, "%s", "@#!PZ2\n");
      else
        sprintf(buf_read, "%s", "@#!AZ2\n");
      fprintf(fo,"%s", buf_read);
      header_out = strlen(buf_read);
      sprintf(buf_read, "source: %dx%d, %s\n", w_data, h_data, data_file);
      fprintf(fo,"%s", buf_read);
      header_out += strlen(buf_read);
      sprintf(buf_read, "lon/lat: %.6f %.6f %.6f %.6f\n", 
        lon1, lon2, lat1, lat2);
      fprintf(fo,"%s", buf_read);
      header_out += strlen(buf_read);
      sprintf(buf_read, "size: %d %d\n", width, height);
      fprintf(fo,"%s", buf_read);
      header_out += strlen(buf_read);
      sprintf(buf_read, "pack: %d %d\n", im_out, jm_out);
      fprintf(fo,"%s", buf_read);
      header_out += strlen(buf_read);
      /* this finishes the header of AZ2/PZ2 format */
      /* just fill the addresses with zeroes as they are still unknown */
      for (i=0; i<=im_out*jm_out; i++) nprintf(fo, total);
      if (verbose)
        fprintf(stderr, 
          "Starting to pack %dx%d chunks of size %dx%d\n",
          im_out, jm_out, xm_out, ym_out);
      total = 0;
      for (j=0; j<jm_out; j++) {
        for (i=0; i<im_out; i++) {
	  l = 0;
          for (y=0; y<ym_out;y++) {
            yp = (int)
              ((double)(h_data-1)-(double)(height-1-j*ym_out-y)*alat-blat);
            if (mode)
            for (x=0; x<xm_out;x++) {
	      xp = (int)((double)(i*xm_out+x)*alon+blon);
              alt = read_input(xp+yp*w_data);
	      source[l] = (alt & 0xFF0000)>>16;
              ++l;
	      source[l] = (alt & 0xFF00)>>8;
	      ++l;
	      source[l] = alt & 0xFF;
	      ++l;
            } else 
            for (x=0; x<xm_out;x++) {
	      xp = (int)((double)(i*xm_out+x)*alon+blon);
              alt = read_input(xp+yp*w_data);
	      source[l] = (alt & 0xFF00)>>8;
	      if (ocean_bit) {
		if (in_ocean)
		  source[l] |= 0x80;
		else
		  source[l] &= 0x7F;
	      }
              ++l;
	      source[l] = alt & 0xFF; 
	      ++l;
            } 
          }
	  l = (2+mode)*xm_out*ym_out+100;
          k = BZ2_bzBuffToBuffCompress(dest, &l,
                source, (2+mode)*xm_out*ym_out, 9, 0, 30);
	  if (k==BZ_OUTBUFF_FULL) {
            fprintf(stderr, "Bzip2 compression failed!! Bailing out...\n");
	    exit(-1);
	  }
          /* write addresses of bz2 chunks immediately after header */
          fseek(fo, header_out+4*(i+im_out*j), SEEK_SET);
          nprintf(fo, total);
          nprintf(fo, total+l);
	  /* write .bz2 compressed buffer */
	  fseek(fo, header_out+4*(im_out*jm_out+1)+total, SEEK_SET);
	  for (k=0; k<l; k++) putc(dest[k], fo);
          total += l;
	  if (verbose)
            fprintf(stderr, "%02x_%02x --> "
                            "%02x_%02x.bz2, size = %d, total = %d\n", 
              i, j, i, j, l, total);
        }
      }
      fclose(fo);
      free(source);
      free(dest);
      exit(0);
    }

    /* write header output if required */
    if (output==PPM) {
      printf("P6\n%d %d\n255\n", width, height);
    } else
    if (output==RAW && header_out) {
      if (input==POP || input==PZ2)
        printf("%s", "@#!POP\n");
      else
        printf("%s", "@#!ALT\n");
      printf("source: %dx%d, %s\n", w_data, h_data, data_file);
      printf("lon/lat: %.6f %.6f %.6f %.6f\n", lon1, lon2, lat1, lat2);
      printf("size: %d %d\n", width, height);
    }

    n_mem_alt = (int *)realloc(n_mem_alt, width*sizeof(int));
    if (dil_y)
      o_mem_alt = (int *)realloc(o_mem_alt, width*sizeof(int));

    i = 0;
    yp = -1;
    for (y=-dil_y; y<height; y++) {
      /* report progress if verbose mode */
      if (verbose && (y%50==0 || y==height-1)) {
        fprintf(stderr, "[%d]", y);
        ++i;
        if (i%10==0)
          fprintf(stderr, "\n");
        else
          fprintf(stderr, " ");
      }
      if (y<0) {
        yp = (int)((double)(h_data-1)-(double)(height-1)*alat-blat);
        if (yp<0) yp = 0;
           if (yp>=h_data) yp = h_data-1;
      } else {
        fyp = (double)(h_data-1)-(double)(height-1-y)*alat-blat;
        n_yp = (int)fyp;
        fyp = fyp-(double)n_yp;
        n_yp = n_yp+dil_y;
        if (n_yp<0) n_yp=0;
        if (n_yp>=h_data) n_yp = h_data-1;
        if (dil_y) {
          if (n_yp==yp) {
            for (x=0; x<width; x++) {
              alt = (int)((1.0-fyp)*o_mem_alt[x]+fyp*n_mem_alt[x]+0.5);
              do_output(alt, x, y);
            }
            continue;
          } else 
            memcpy(o_mem_alt, n_mem_alt, width*sizeof(int));
        } else {
          if (n_yp==yp) {
            for (x=0; x<width; x++) {
              alt = n_mem_alt[x];
              do_output(alt, x, y);
            }
            continue;
          }
        }
        yp = n_yp;
      }
      if (input==POP || input==PZ2)
        invcosphi = 1.0/sin(M_PI*((yp+0.5)/(double)h_data));
      if (dil_x) {
        xp = ((int)blon)%w_data;
        n_alt = read_input(xp+yp*w_data);
      }
      for (x=0; x<width; x++) {
        if (output==XWIN && XCheckIfEvent(dpy, &event, evpred, (XPointer)0)) {
          k = manage_event(event);
          if (k==1) {
	    reiterate = 1;
	    goto freebuffers;
	  }
          if (k==2) goto finish;
        }
        fxp = ((double)x)*alon+blon;
        n_xp = (int)fxp;
        fxp = fxp - (double)n_xp;
        n_xp = (n_xp+dil_x)%w_data;
        if (n_xp==xp) goto samevalue;
        xp = n_xp;
        o_alt = n_alt;
        n_alt = read_input(xp+yp*w_data);
      samevalue:
        if (dil_x)
          alt = (int)((1.0-fxp)*o_alt+fxp*n_alt+0.5);
        else
          alt = n_alt;
        n_mem_alt[x] = alt;
        if (y<0) continue;
        if (dil_y)
          alt = (int)((1.0-fyp)*o_mem_alt[x]+fyp*alt+0.5);
        do_output(alt, x, y);
      }
    }

  freebuffers:
    /* free temporary buffers */
    if (dil_y) {
      free(o_mem_alt);
      o_mem_alt = NULL;
    }
    free(n_mem_alt);
    n_mem_alt = NULL;
    if (input==AZ2 || input==PZ2) {
      fclose(fd);
      for (i=0; i<im_data; i++)
      for (j=0; j<jm_data; j++)
        if (buf_out[i+j*im_data]) {
	  free(buf_out[i+j*im_data]);
	  buf_out[i+j*im_data] = 0;
	}
    } else
    if (input==GTOPO30)
      for (i=0; i<33; i++) fclose(gfd[i]);
    else {
      if (gzd) gzclose(gzd);
      gzd = NULL;
    }
 restart:
    if (cps && reiterate) {
      input = CPS;
      header_data = 0;
      goto cpsinit;
    }
    if (reiterate) goto initialize;

    if (output==XWIN) {
      XSetWindowBackgroundPixmap(dpy, win, pixmap);
      XClearWindow(dpy, win);
      while(1) {
        XNextEvent(dpy, &event);
	i = manage_event(event);
	if (i==1) {
	  reiterate = 1;
          goto restart;
	}
	if (i==2) break;
      }
    }

 finish:
    /* free X11 objects */
    if (output==XWIN) {
      XDestroyWindow(dpy, win);
      XFreePixmap(dpy, pixmap);
      XFreeGC(dpy, gc);
    }

    /* free pixel and color buffers */
    if (pix) free(pix);
    free(red);
    free(green);
    free(blue);

    if (input==AZ2 || input==PZ2) {
      /* free all buffers that are still open */
      for (i=0; i<im_data; i++)
      for (j=0; j<jm_data; j++)
        if (buf_out[i+j*im_data])
	  free(buf_out[i+j*im_data]);
      free(buf_out);
    }

    if (cps) {
      if (cps_dir) free(cps_dir);
      for (i=0; i<num_cps; i++) free(cps[i].name);
      free(cps);
    }

    if (verbose) {
      if (input==POP || input==PZ2)
        fprintf(stderr, 
          "\n\nMinimum/maximum population density %.3f %.3f\n\n", 
          min_den, max_den);
      else
        fprintf(stderr, 
          "\n\nMinimum/maximum altitude %d %d\n\n", min_alt, max_alt);
    }
    exit(0);
}