#include <stdio.h>		/* use printf() */
#include <stdlib.h>		/* use malloc(), realloc() */

#ifdef _WIN32
#include <fcntl.h>		/* use setmode()'s parameter macros */
#include <io.h>			/* use setmode() */
#include <excpt.h>
#endif

#ifdef MAC
#include <signal.h>		/* use sigaction() */
#endif

#include "mslib.h"
#include "msmsAllProto.h"	/* Use MS_getTri() */

#undef APP_DOES_FREE		/* undefined means quicker OS deallocation */

static int read_binary_xyzr(FILE *f, float ***xyzr, int *n);
static int compute_surface_areas(MOLSRF *surface);
static int write_binary_surface(FILE *f, MOLSRF *surface);
#ifndef APP_DOES_FREE
static void free_coords(float **coords, int nc);
#endif

int main(int argc, char **argv)
{
	int	num_atoms;
	float	**coords;
	char	**atom_names = NULL;
	int	status;
	MOLSRF	*surface;
	MOLSRF_ROOT *msr;
	const char *filename = "<stdin>";
	float probe_radius, vertex_density;
	int all_components;

#ifdef MAC
	/* Suppress crash reporter dialog on Mac 10.5. */
	struct sigaction sa;
	sa.sa_flags = 0;
	sigemptyset(&sa.sa_mask);
	sa.sa_handler = exit;
	sigaction(SIGBUS, &sa, NULL);
	sigaction(SIGSEGV, &sa, NULL);
#endif

#ifdef _WIN32
	/*
	 * Need to change stdin from text to binary mode, otherwise ctrl-z
	 * in input terminates the input stream.
	 */
	setmode(fileno(stdin), O_BINARY);
	setmode(fileno(stdout), O_BINARY);
#endif

	if (argc != 4)
	  {
	    fprintf(stderr, "Syntax: mscalc <probe-radius> <vertex-density> <all-components>\n");
	    return 1;
	  }
	probe_radius = atof(argv[1]);
	vertex_density = atof(argv[2]);
	all_components = atoi(argv[3]);

	if (!read_binary_xyzr(stdin, &coords, &num_atoms))
	  {
	    fprintf(stderr, "Failed reading xyzr coordinates\n");
	    return 1;
	  }

#ifdef _WIN32
	__try {
#endif
	MS_init_msmslib();
	msr = MS_init_molsrf_root();
	surface = MS_add_molsrf(msr, (char *)filename, coords, num_atoms,
				0, atom_names);
	surface->all_components = all_components;
	status = MS_compute_surface(surface, probe_radius, vertex_density);
#ifdef _WIN32
	} __except (EXCEPTION_EXECUTE_HANDLER) {
	  /* Catch structured exceptions so the user doesn't get a dialog
	   * asking to send a bug report to Microsoft */
	  _exit(5);
	}
#endif
	if (status != MS_OK)
	  {
	    fprintf(stderr, "Failed calculating surface.\n");
	    return 2;
	  }
	if (!compute_surface_areas(surface))
	  {
	    fprintf(stderr, "Failed computing per-atom areas.\n");
	    return 3;
	  }
	if (!write_binary_surface(stdout, surface))
	  {
	    fprintf(stderr, "Failed writing triangulated surfaces.\n");
	    return 4;
	  }

#ifdef APP_DOES_FREE
	free_coords(coords, num_atoms);
#endif
	return 0;
}

static int read_binary_xyzr(FILE *f, float ***xyzr, int *n)
{
  int nxyzr, c;
  float **coords;

  if (fread(&nxyzr, sizeof(int), 1, f) != 1)
    return 0;

  coords = (float **) malloc(nxyzr * sizeof(float *));
  for (c = 0 ; c < nxyzr ; ++c)
    {
      coords[c] = (float *) malloc(4*sizeof(float));
      if (fread(coords[c], sizeof(float), 4, f) != 4)
	return 0;	/* TODO: free coordinates. */
    }
  *xyzr = coords;
  *n = nxyzr;
  return 1;
}

static int compute_surface_areas(MOLSRF *surface)
{
  int nc, c;
  RS *rs;

  nc = surface->rsr->nb;
  for (c = 0 ; c < nc ; ++c)
    {
      rs = MS_find_rs_component_by_num(surface, c);
      if (MS_compute_numerical_area_vol(surface, rs->ses, MS_SEMI_ANALYTICAL)
	  != MS_OK)
	return 0;
    }
  if (MS_compute_SES_area(surface) != MS_OK)
    return 0;
  return 1;
}

static int write_binary_surface(FILE *f, MOLSRF *surface)
{
  float *v, *area;
  int nc, c, i, nv, nf, *tri, *ind, selnum=3, base=0, natm;
  unsigned char *flag;
  RS *rs;
  RSV *atm;

  natm = surface->nbat;
  area = (float *)malloc(2*natm * sizeof(float));
  flag = (unsigned char *)malloc(natm * sizeof(unsigned char));
  for (i = 0 ; i < natm ; i++) flag[i] = 1;

  /* fprintf(stderr, "%d components\n", nc); */
  nc = surface->rsr->nb;
  fwrite(&nc, sizeof(int), 1, f);
  atm = surface->rsr->atm;
  for (c = 0 ; c < nc ; ++c)
    {
      rs = MS_find_rs_component_by_num(surface, c);
      /* fprintf(stderr, "component %d has %d triangles\n", c, rs->ses->nbtri); */
      /* allocate array for atom flags */
      if (MS_getTri( flag, rs->ses, base, selnum, &nv, &v, &ind, &nf, &tri)
	  != MS_OK)
	return 0;
      fwrite((void *)&nv, sizeof(int), 1, f);
      fwrite((void *)v, sizeof(float), nv*8, f);
      fwrite((void *)&nv, sizeof(int), 1, f);
      fwrite((void *)ind, sizeof(int), nv*3, f);
      fwrite((void *)&nf, sizeof(int), 1, f);
      fwrite((void *)tri, sizeof(int), nf*5, f);
      for (i = 0 ; i < natm ; ++i)
	{
	  area[2*i] = atm[i].ses_area[c];
	  area[2*i+1] = atm[i].sas_area[c];
	}
      fwrite((void *)&natm, sizeof(int), 1, f);
      fwrite((void *)area, sizeof(float), natm*2, f);
    }

#ifdef APP_DOES_FREE
  free(flag);
  free(area);
#endif

  return 1;
}

#ifdef APP_DOES_FREE
static void free_coords(float **coords, int nc)
{
  int c;
  for (c = 0 ; c < nc ; ++c)
    free(coords[c]);
  free(coords);
}
#endif