/*
 * xyz2rgb.c - convert an XYZ file (or a concatenation of several XYZ files)
 *             to a ImageMagick-readable bitmap file (or several bitmaps)
 *
 * Orientation: when longitude and latitude are both zero (the default case)
 * the camera lies on the x axis. The camera always faces the origin. The z
 * axis always points up, and for zero longitude, the y axis points to the
 * right. Increasing longitude will swing the camera around toward the y axis.
 * There are two other parameters for controlling the camera, which are the
 * width of the viewing area (an angle) and the distance from the camera to
 * the origin.
 *
 * The viewing area is always assumed to have an aspect ratio of 4:3.
 *
 * $Id: xyz2rgb.c,v 1.5 1997/05/22 01:44:01 wware Exp $
 */

#include 
#include 
#include 
#include 

char usage_string[] =
"Usage: xyz2rgb.c [options] foo.xyz\n"
"Options:\n"
"    -longitude a    specify horizontal angle of camera\n"
"    -latitude a     specify vertical angle of camera\n"
"    -distance x     specify distance to origin of camera\n"
"    -width a        specify wide-angleness of camera\n"
"    -hsize n        specify number of horizontal pixels\n"
"    -vsize n        specify number of vertical pixels\n"
"    -radii [c,v]    select covalent/VanDerWaals atomic radii\n"
"    -o filename     output filename, files will have numbers\n"
"                     and \".rgb\" appended to their names\n"
;

typedef struct
  {
    double longitude, latitude, distance, width;
  }
camera;

typedef struct
  {
    unsigned char r, g, b;
    double radius;
    char name;
  }
species;

typedef struct
  {
    double x, y, z, radius;
    species *spc;
  }
atom;

#define MAXHSIZE 1000
unsigned char rbuf[MAXHSIZE], gbuf[MAXHSIZE], bbuf[MAXHSIZE];
double xbuf[MAXHSIZE];

species known_species[] =
{
  {
    255, 255, 255, 1.5, 'H'
  }
  ,
  {
    128, 128, 128, 1.94, 'C'
  }
  ,
  {
    255, 0, 0, 1.82, 'O'
  }
  ,
  {
    0, 0, 255, 1.74, 'N'
  }
};

int num_known_species = sizeof (known_species) / sizeof (species);

#ifdef DEBUG
#define H    &known_species[0]
#define C    &known_species[1]
#define O    &known_species[2]
#define N    &known_species[3]
atom atoms[] =
{
  {
    17.047, 14.099, 3.625, 0.0, N
  }
  ,
  {
    16.967, 12.784, 4.338, 0.0, C
  }
  ,
  {
    15.685, 12.755, 5.133, 0.0, C
  }
  ,
  {
    15.268, 13.825, 5.594, 0.0, O
  }
  ,
  {
    18.170, 12.703, 5.337, 0.0, H
  }
  ,
  {
    19.334, 12.829, 4.463, 0.0, O
  }
  ,
  {
    18.150, 11.546, 6.304, 0.0, H
  }
  ,
  {
    15.115, 11.555, 5.265, 0.0, N
  }
  ,
  {
    13.856, 11.469, 6.066, 0.0, H
  }
  ,
  {
    14.164, 10.785, 7.379, 0.0, C
  }
  ,
  {
    14.993, 9.862, 7.443, 0.0, O
  }
  ,
  {
    12.732, 10.711, 5.261, 0.0, C
  }
  ,
  {
    13.308, 9.439, 4.926, 0.0, O
  }
  ,
  {
    12.484, 11.442, 3.895, 0.0, C
  }
  ,
  {
    13.488, 11.241, 8.417, 0.0, N
  }
  ,
  {
    13.660, 10.707, 9.787, 0.0, H
  }
  ,
  {
    10.646, 8.991, 11.408, 0.0, C
  }
  ,
  {
    10.654, 8.793, 12.919, 0.0, C
  }
  ,
  {
    11.659, 8.296, 13.491, 0.0, O
  }
  ,
  {
    10.057, 7.752, 10.682, 0.0, C
  }
  ,
  {
    9.837, 8.018, 8.904, 0.0, O
  }
  ,
  {
    9.561, 9.108, 13.563, 0.0, N
  }
  ,
  {
    9.448, 9.034, 15.012, 0.0, C
  }
};
#else
atom *atoms;
#endif

static void
transform_coords (atom * atm, camera * cam)
{
  double angle, tmp, pfactor;

  /* rotate by longitude */
  angle = cam->longitude;
  tmp = atm->x * cos (angle) + atm->y * sin (angle);
  atm->y = atm->y * cos (angle) - atm->x * sin (angle);
  atm->x = tmp;

  /* rotate by longitude */
  angle = cam->latitude;
  tmp = atm->x * cos (angle) + atm->z * sin (angle);
  atm->z = atm->z * cos (angle) - atm->x * sin (angle);
  atm->x = tmp;

  /* apply perspective */
  pfactor = (cam->distance - atm->x) / cam->distance;
  pfactor = (pfactor < 0.01) ? 100.0 : (1.0 / pfactor);
  atm->y *= pfactor;
  atm->z *= pfactor;
  atm->radius *= pfactor;
}

static void
render_atom (atom * a, double z, double ymax, int half_hsize)
{
  double local_radius, rsq, zsq, y1, y2;
  int i, i1, i2;

  rsq = a->radius * a->radius;
  zsq = (z - a->z) * (z - a->z);
  local_radius = rsq - zsq;
  if (local_radius < 0.0)
    return;
  local_radius = sqrt (local_radius);
  y1 = a->y - local_radius;
  if (y1 > ymax)
    return;
  y2 = a->y + local_radius;
  if (y2 < -ymax)
    return;
  i1 = (int) (half_hsize * (y1 / ymax));
  i2 = (int) (half_hsize * (y2 / ymax));
  for (i = i1 + half_hsize; i <= i2 + half_hsize; i++)
    {
      double y, ysq, xblip;
      y = (i - half_hsize) * (ymax / half_hsize);
      ysq = (y - a->y) * (y - a->y);
      xblip = rsq - ysq - zsq;
      if (xblip < 0.0)
	continue;
      xblip = sqrt (xblip);
      if (xbuf[i] < a->x + xblip)
	{
	  double color_intensity = xblip / a->radius;
	  xbuf[i] = a->x + xblip;
	  rbuf[i] = (unsigned char) (color_intensity * a->spc->r);
	  gbuf[i] = (unsigned char) (color_intensity * a->spc->g);
	  bbuf[i] = (unsigned char) (color_intensity * a->spc->b);
	}
    }
}

int hsize, vsize, num_atoms;
double ymax, zmax;
camera cam;

void
render_structure (FILE * outf)
{
  int i, h, v;

  /* center the structure */
  {
    double xc, yc, zc;
    xc = yc = zc = 0.0;
    for (i = 0; i < num_atoms; i++)
      {
	xc += atoms[i].x;
	yc += atoms[i].y;
	zc += atoms[i].z;
      }
    xc /= num_atoms;
    yc /= num_atoms;
    zc /= num_atoms;
    for (i = 0; i < num_atoms; i++)
      {
	atoms[i].x -= xc;
	atoms[i].y -= yc;
	atoms[i].z -= zc;
      }
  }

  for (i = 0; i < num_atoms; i++)
    {
      atoms[i].radius = atoms[i].spc->radius;
      transform_coords (&atoms[i], &cam);
    }

  for (v = 0; v < vsize; v++)
    {
      for (h = 0; h < hsize; h++)
	{
	  /* sky blue background */
	  rbuf[h] = 192;
	  gbuf[h] = 192;
	  bbuf[h] = 255;
	  xbuf[h] = -10000000000000.0;
	}

      /* draw atoms */
      for (i = 0; i < num_atoms; i++)
	{
	  double z = (-2.0 * zmax / vsize) * v + zmax;
	  render_atom (&atoms[i], z, ymax, hsize / 2);
	}

      for (i = 0; i < hsize; i++)
	fputc (rbuf[i], outf);

      for (i = 0; i < hsize; i++)
	fputc (gbuf[i], outf);

      for (i = 0; i < hsize; i++)
	fputc (bbuf[i], outf);
    }
}

int
ouch (char *s)
{
  fprintf (stderr, "%s", s);
  return 1;
}

char line[200], infname[40], fname_prefix[40], fname[40];

int
main (int argc, char *argv[])
{
  int i, j, fn;
  FILE *inf, *outf;

  /* initialize all options to default values */
  cam.longitude = 0.0;
  cam.latitude = 0.0;
  cam.distance = 30.0;
  cam.width = 30.0 * (3.14159 / 180.0);
  hsize = 640;
  vsize = 480;
  infname[0] = '\0';
  fname_prefix[0] = '\0';

  if (argc < 2)
    return ouch (usage_string);

  /* collect command line options */
  for (i = 1; i < argc; i++)
    {
      if (strcmp (argv[i], "-longitude") == 0)
	{
	  double tmp;
	  if (sscanf (argv[++i], "%lf", &tmp) == 1)
	    cam.longitude = tmp * (3.14159 / 180.0);
	}
      else if (strcmp (argv[i], "-latitude") == 0)
	{
	  double tmp;
	  if (sscanf (argv[++i], "%lf", &tmp) == 1)
	    cam.latitude = tmp * (3.14159 / 180.0);
	}
      else if (strcmp (argv[i], "-distance") == 0)
	{
	  double tmp;
	  if (sscanf (argv[++i], "%lf", &tmp) == 1)
	    cam.distance = tmp;
	}
      else if (strcmp (argv[i], "-width") == 0)
	{
	  double tmp;
	  if (sscanf (argv[++i], "%lf", &tmp) == 1)
	    cam.width = tmp * (3.14159 / 180.0);
	}
      else if (strcmp (argv[i], "-hsize") == 0)
	{
	  int tmp;
	  if (sscanf (argv[++i], "%d", &tmp) == 1)
	    hsize = tmp;
	}
      else if (strcmp (argv[i], "-vsize") == 0)
	{
	  int tmp;
	  if (sscanf (argv[++i], "%d", &tmp) == 1)
	    vsize = tmp;
	}
      else if (strcmp (argv[i], "-o") == 0)
	{
	  strcpy (fname_prefix, argv[++i]);
	}
      else
	{
	  strcpy (infname, argv[i]);
	}
    }

  ymax = cam.distance * tan (cam.width);
  zmax = (3.0 / 4.0) * ymax;

#ifdef DEBUG
  num_atoms = sizeof (atoms) / sizeof (atom);
  outf = stdout;
  render_structure (outf);
  if (outf != stdout)
    fclose (outf);
  return 0;
#else
  if (strlen(infname) == 0)
    inf = stdin;
  else
    inf = fopen(infname, "r");
  if (inf == NULL)
    return ouch ("Can't open input file\n");

  /* find out how many atoms to grab */
  fn = 0;
  while (fgets (line, 200, inf) != NULL)
    {
      char *s;
      int found;

      if (sscanf (line, "%d", &num_atoms) < 1)
	return 0;
      atoms = malloc (num_atoms * sizeof (atom));
      if (atoms == NULL)
	return ouch ("Not enough memory for so many atoms!\n");
      fgets (line, 200, inf);
      for (i = 0; i < num_atoms; i++)
	{
	  fgets (line, 200, inf);
	  s = line;
	  while (*s == ' ' || *s == '\t')
	    s++;
	  for (j = found = 0; !found && j < num_known_species; j++)
	    if (*s == known_species[j].name)
	      {
		found = 1;
		atoms[i].spc = &known_species[j];
	      }
	  if (!found)
	    atoms[i].spc = &known_species[0];
	  s++;
	  while (*s == ' ' || *s == '\t')
	    s++;
	  sscanf(s, "%lf %lf %lf", &atoms[i].x, &atoms[i].y, &atoms[i].z);
	}

      if (strlen(fname_prefix) > 0)
	{
	  sprintf(fname, "%s%03d.rgb", fname_prefix, fn++);
	  outf = fopen(fname, "w");
	}
      else
	outf = stdout;
      if (outf == NULL)
	return ouch ("Can't open output file\n");
      render_structure (outf);
      if (outf != stdout)
	fclose (outf);
    }

  if (inf != stdin)
    fclose (inf);
  return 0;
#endif
}