/* planet2.c */
/* planet generating program */
/* Copyright 1993 Torben AE. Mogensen */
/* this version uses integers */

/* See Usage section further below, for how to use.

/* The method used is a spatial subdivision of a box containing the globe */

/* I have tried to avoid using machine specific features, so it should */
/* be easy to port the program to any machine. Beware, though that due */
/* to different precision on different machines, the same seed numbers */
/* can yield very different planets. */
/* The very primitive user interface is a result of portability concerns */


/* slightly modified by Jouni Aro <jaro@hut.fi> in December 1993 */

/* I added commandline options and changed the output format to */
/* raw altitude data, which can be easily processed further. */
/* Look for the utility pla2gif that converts the output to gif format */

/* This one compiles for sure with Borland Compilers, probably in */
/* most msdos and also unix environments */

/* All comments and complaints are welcome */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#ifdef __MSDOS__
#include <conio.h>
#include <io.h>
#include <fcntl.h>

extern int getopt(int ac, char **av, char *optstr);

#else
#include <unistd.h>
#endif

#define MAXALT   32767
#define MINALT  -32768

#define WIDTH 320
#define HEIGHT 200

/* these three values can be changed to change world characteristica */
int     M    = -246;   /* initial altitude (slightly below sea level) */
int     Dda0 =  3277;  /* initial altitude variation */
int     Ddda =  13926; /* rate of decrease of altitude variation */

#define Da0     32768   /* initial scale (size of box) */
#define DEG2RAD 0.017453292 /* pi/180 */
#define PI      3.1415926535897

#define Maxlevel 64

/* prototypes */
void fatal(char *msg);
void usage(void);
int planet0(double x,double y,double longi,double lat);
long planet1(long x,long y,long z, int level);
long planet(long x,long y,long z, int level);
long rand2(long p,long q);
double log2(double x);
extern void alt_encode(int alt, FILE *fp); /* in alt_code.c */
extern char *VersionStr;                   /* this also     */

int Depth; /* depth of subdivisions */

long t1,t2,t3;
long r1,r2,r3,r4,r5,r6,r7,r8; /* seeds */
double longi,lat,scale;

typedef enum { Mercator, Global, Scaled } ViewType;

struct plan {long a,b,c,d,e,f,g,h;         /* altitudes of the eight corners */
	     long as,bs,cs,ds,es,fs,gs,hs; /* seeds of the eight corners */
	     long ax,ay,az;                /* bottom left closest corner */
	     long ae,bf,cg,dh,aes,bfs,cgs,dhs;
	     int direction;}
       levels[Maxlevel];

long scales[Maxlevel], dalt[Maxlevel];

char *prog_name;

void fatal(char *msg)
{
  fprintf(stderr, "\n%s: %s\n", prog_name, msg);
  exit(1);
}

/*** USAGE ***/
void usage()
{
  fprintf(stderr,
"Usage: %s [-s width height] [-v view [long lat scale]]\n\
        [-q quietmode] [-h] seed [filename]\n\n", prog_name);
  fprintf(stderr,
"     seed is a long integer value\n\
     default map size is 320x200 \n\
     view is m (or 0), g (or 1)  or s (or 2). \n\
     - m gives a Mercator projection world map \n\
     - g gives an orthogonal projection from given angle \n\
       extra parameters long,lat are used as \n\
       longitude and latitude (in degrees) for centre of picture \n\
     - s means orthogonal projection from a given point to a given scale \n\
       extra parameters long,lat,scale are used as \n\
       longitude and latitude (in degrees) for centre of picture \n\
       and a scaling factor (1.0 ~ full view)\n");
  exit(1);
}

/*** MAIN ***/
int main(ac,av)
int ac;
char **av;
{
  char line[80];
  int i,j, width, height, c, longitude, latitude;
  ViewType view;
  double x,y;
  long seed;
  long tm;
  FILE *out;

  extern int optind;
  extern char *optarg;
  int opt;

#ifdef __MSDOS__
  int quiet = 0;
  if ((prog_name = strrchr(av[0], '\\')) == NULL)
    prog_name = av[0];
  else
    prog_name++;
#else
  int quiet = 1;
  prog_name = av[0];
#endif
  
  width = WIDTH;
  height = HEIGHT;
  view = Mercator;
  longitude = latitude = 0;
  scale = 1.0;

  do {
    opt = getopt(ac, av, "s:v:q:h");

    switch(opt) {
    case -1 : break;
    case 's': 
      width = atoi(optarg); 
      height = atoi(av[optind++]); 
      break;
    case 'v':
      switch (*optarg) {
      case '0':
      case 'm':
        view = Mercator;
        scale = 1.0;
        break;
      case '1':
      case 'g':
        view = Global;
        longitude = atoi(av[optind++]);
	latitude = atoi(av[optind++]);
	scale = 1.5; /* I think this is best <jaro> */
        break;
      case '2':
      case 's':
        view = Scaled;
        longitude = atoi(av[optind++]);
	latitude = atoi(av[optind++]);
	scale = atof(av[optind++]);
        break;
      default:
        usage();
      }
      break;
    case 'q': 
      quiet = atoi(optarg); 
      break;
    case 'h':
    default: 
      usage();
    }
  } while (opt != -1);


  if (optind < ac)
    seed = atof(av[optind++]);   /* seed for random number generator */
  else 
    usage();

  if (optind < ac) {
    if((out = fopen(av[optind++], "wb")) == NULL)
      fatal("Couldn't open outfile.");
  }
  else {
    out = stdout;
#ifdef __MSDOS__
/* Make sure it's in binary form */
    setmode(1, O_BINARY); 
#endif
  }

  if (view == Mercator)
    sprintf(line,"%s %ld %dx%d Mercator", 
            VersionStr, seed, width, height);
  else if (view == Global)
    sprintf(line,"%s %ld %dx%d Global from (%d,%d)", 
            VersionStr, seed, width, height,
            longitude, latitude);
  else /* view == Scaled */
    sprintf(line,"%s %ld %dx%d Scaled %.2fX at (%d,%d)", 
            VersionStr, seed, width, height, 
            scale, longitude, latitude);
  fputs(line, out);
  fputc('\n', out);

  if (!quiet) fprintf(stderr, "\n%s:     ", line);

  longi = longitude*DEG2RAD;
  lat   = latitude*DEG2RAD;

  Depth = 3*((int)(log2(scale*HEIGHT)))+3;

  r1 = rand2(seed,seed);
  r2 = rand2(r1,r1);
  r3 = rand2(r1,r2);
  r4 = rand2(r2,r3);
  r5 = rand2(r4,r3);
  r6 = rand2(r4,r5);
  r7 = rand2(r6,r5);
  r8 = rand2(r6,r7);
  t1 = 16384 * exp(log(2.)/3.);    /* dimension of box = (scale,t1*scale,t2*scale) */
  t2 = 16384 * exp(2.*log(2.)/3.); /* so halving will keep the same relative scale */
  t3 = t2/2;                      /* like A4 paper (but in 3 dimensions) */

  tm = Da0;
  for (i=0; i<Maxlevel; i++) {
    tm = tm*t3>>14;
    scales[i] = tm;
  }

  tm = Dda0;
  for (i=0; i<Maxlevel; i++) {
    dalt[i] = tm;
    tm = tm*Ddda>>14;
  }

  levels->a = levels->b = levels->c = levels->d =
  levels->e = levels->f = levels->g = levels->h = M;

  levels->as = r1; levels->bs = r2; levels->cs = r3; levels->ds = r4;
  levels->es = r5; levels->fs = r6; levels->gs = r7; levels->hs = r8;

  levels->ax = -16384; levels->ay = -t1; levels->az = -t2;

  levels->direction = 0;

  for (j = 0; j < height; j++) {

    for (i = 0; i < width ; i++) {
      switch (view) {

        case Mercator: /* Mercator projection */
	  y = PI*1.025*(2.0*j-height)/width;
	  y = exp(2.*y);
	  y = (y-1.)/(y+1.);
	  /* old version: y = y/sqrt(1.0+y*y); */
	  scale = (double)width/height/sqrt(1.0-y*y)/PI/1.05;
	  Depth = 3*((int)(log2(scale*height)))+3;
	  c = planet0(0.0,y,PI*1.025*(2.0*i-width)/width,0.0);
          break;

        case Global: /* global view */
          x = 2*(i-(double)height/2)/height; 
          y = 2*(j-(double)height/2)/height;
          c = planet0(x,y,longi,lat);
          break;
      
        case Scaled: /* use longitude and latitude */
	  c = planet0((2.0*i-width)/height/scale,
		      (2.0*j-height)/height/scale,longi,lat);
          break;

      }
      alt_encode(c, out);

#ifdef __MSDOS__
      if(kbhit() && getch() == 'q') {
/*	fprintf(out,"*break");*/
	exit(1);	
      }
#endif
    }
    if(!quiet) fprintf(stderr, "\b\b\b\b%-4d", j);
  }
  /* end picture */
  fclose(out);
  if(!quiet) fputc('\n', stderr);
  
  return 0;
}
	
int planet0(double x,double y,double longi,double lat)
{
  long x1,y1,z1,alt;
  double theta1,theta2,cos2;
  if (x*x+y*y >= 1.0) 
    alt = MAXALT;
  else {
    theta2 = asin(y)-lat;
    if (fabs(y) > 0.999) theta1 = 0;
    else theta1 = acos(x/sqrt(1.0-y*y))-longi;
    cos2 = cos(theta2);
    x1 = 16384 * cos(theta1)*cos2;
    y1 = 16384 * sin(theta1)*cos2;
    z1 = 16384 * sin(theta2);
  
    alt = planet1(x1,y1,z1,
	             /* coordinates of point we want colour of */
		 0);
    if (alt >= MAXALT) /* not in int limits */
      alt = MAXALT-1;
    else if (alt <= MINALT) 
      alt = MINALT;
  }

  return (int) alt;
}

long planet1(long x,long y,long z, int level)
             /* goal point */
{
  if (level == Depth)
    return(planet(x,y,z,level));
  else
    if (levels[level].az+scales[level] >= z)
      if (levels[level].direction == 1)
	return(planet1(z,x,y,level+1));
      else
	return(planet(x,y,z,level));
    else
      if (levels[level].direction == 2)
	return(planet1(z,x,y,level+1));
      else
	return(planet(x,y,z,level));
}

long planet(long x,long y,long z, int level)
           /* goal point */
{
  struct plan *level1, *level2;

  level1 = levels+level;
  if (level == Depth) {  
    level1->direction = 0;
    return((level1->a+level1->b+level1->c+level1->d+
	    level1->e+level1->f+level1->g+level1->h)/8);
  }
  else { /* subdivide over longest axis */
        /* find seeds for new points */
    level1->aes = rand2(level1->as,level1->es);
    level1->bfs = rand2(level1->bs,level1->fs);
    level1->cgs = rand2(level1->cs,level1->gs);
    level1->dhs = rand2(level1->ds,level1->hs);
        /* find altitudes as mean+variation */
    level1->ae = (level1->a+level1->e+(dalt[level]*level1->aes>>14))/2;
    level1->bf = (level1->b+level1->f+(dalt[level]*level1->bfs>>14))/2;
    level1->cg = (level1->c+level1->g+(dalt[level]*level1->cgs>>14))/2;
    level1->dh = (level1->d+level1->h+(dalt[level]*level1->dhs>>14))/2;
    level2 = level1+1;
	/* find which halfbox goal point is in and continue in that */
    if (level1->az+scales[level] >= z) {
      level1->direction = 1;
      level2->a = level1->a;
      level2->b = level1->ae;
      level2->c = level1->bf;
      level2->d = level1->b;
      level2->e = level1->d;
      level2->f = level1->dh;
      level2->g = level1->cg;
      level2->h = level1->c;
      level2->as = level1->as;
      level2->bs = level1->aes;
      level2->cs = level1->bfs;
      level2->ds = level1->bs;
      level2->es = level1->ds;
      level2->fs = level1->dhs;
      level2->gs = level1->cgs;
      level2->hs = level1->cs;
      level2->ax = level1->az;
      level2->ay = level1->ax;
      level2->az = level1->ay;
      return(planet(z,x,y, level+1));
    }
    else {
      level1->direction = 2;
      level2->a = level1->ae;
      level2->b = level1->e;
      level2->c = level1->f;
      level2->d = level1->bf;
      level2->e = level1->dh;
      level2->f = level1->h;
      level2->g = level1->g;
      level2->h = level1->cg;
      level2->as = level1->aes;
      level2->bs = level1->es;
      level2->cs = level1->fs;
      level2->ds = level1->bfs;
      level2->es = level1->dhs;
      level2->fs = level1->hs;
      level2->gs = level1->gs;
      level2->hs = level1->cgs;
      level2->ax = level1->az+scales[level];
      level2->ay = level1->ax;
      level2->az = level1->ay;
      return(planet(z,x,y, level+1));
    }
  }
}

long rand2(long p,long q)
/* random number generator taking two seeds */
/* rand2(p,q) = rand2(q,p) is important     */
{
  long r;
  r = (p+q+37415)<<1;
  if (!r || r&0x2000) r^=51473;
  return((r&32767)-16384);
}

double log2(double x)
{ return(log(x)/log(2)); }