# include <ohana.h>
# include <gfitsio.h>

char *print_table_row (char *row, Header *header);
FILE *load_extension (char *file, int Nextend, char *Extname, Header *header);
void print_column (FTable *table, int Column, char *Colname);
void usage();
void list_extnames (char *file);
void print_layout (Header *header);
int   ByteSwap (char *ptr, int size, int nitems, char *type);
int Binary;

int main (int argc, char **argv) {

  int i, N, Nx, Ny, Nbytes, Nread;
  int Nextend, Column, Row, ListExtname, Layout;
  char *Extname, *Colname, *line, ttype[80];
  FTable table;
  Header header;
  FILE *f;

  if (get_argument (argc, argv, "-h")) usage ();
  if (get_argument (argc, argv, "--help")) usage ();

  Nextend = 0;
  if ((N = get_argument (argc, argv, "-x"))) {
    remove_argument (N, &argc, argv);
    Nextend = atoi (argv[N]);
    remove_argument (N, &argc, argv);
  }

  Extname = (char *) NULL;
  if ((N = get_argument (argc, argv, "-n"))) {
    if (Nextend) usage ();
    remove_argument (N, &argc, argv);
    Extname = strcreate (argv[N]);
    remove_argument (N, &argc, argv);
  }

  Row = 0;
  if ((N = get_argument (argc, argv, "-row"))) {
    remove_argument (N, &argc, argv);
    Row = atof (argv[N]);
    remove_argument (N, &argc, argv);
  }

  Column = 0;
  if ((N = get_argument (argc, argv, "-ncolumn"))) {
    remove_argument (N, &argc, argv);
    Column = atof (argv[N]);
    remove_argument (N, &argc, argv);
  }

  Colname = 0;
  if ((N = get_argument (argc, argv, "-column"))) {
    if (Column) usage ();
    remove_argument (N, &argc, argv);
    Colname = strcreate (argv[N]);
    remove_argument (N, &argc, argv);
  }

  ListExtname = FALSE;
  if ((N = get_argument (argc, argv, "-list"))) {
    remove_argument (N, &argc, argv);
    ListExtname = TRUE;
  }

  Layout = FALSE;
  if ((N = get_argument (argc, argv, "-layout"))) {
    remove_argument (N, &argc, argv);
    Layout = TRUE;
  }

  if (argc != 2) usage ();

  if (ListExtname) list_extnames (argv[1]);

  /* load header */
  table.header = &header;
  f = load_extension (argv[1], Nextend, Extname, table.header);

  if (Layout) print_layout (table.header);

  Binary = FALSE;
  gfits_scan (table.header, "XTENSION", "%s", 1, ttype);
  if (!strcmp (ttype, "BINTABLE")) Binary = TRUE;

  /* load table data array */
  Nbytes = gfits_matrix_size (table.header);
  ALLOCATE (table.buffer, char, Nbytes);
  Nread = fread (table.buffer, sizeof (char), Nbytes, f);
  if (Nread != Nbytes) {
    fprintf (stderr, "failed to read all table data\n");
    exit (1);
  }
  table.size = Nbytes;

  gfits_scan (table.header, "NAXIS1",  "%d", 1, &Nx);
  gfits_scan (table.header, "NAXIS2",  "%d", 1, &Ny);

  /* print a column */
  if (Column || (Colname != (char *) NULL)) print_column (&table, Column, Colname);

  /* print a row */
  if (Row) {
    line = print_table_row (&table.buffer[Nx*Row], table.header);
    fprintf (stdout, "%s\n", line);
    free (line);
    exit (0);
  }

  /* print complete table */
  for (i = 0; i < Ny; i++) {
    line = print_table_row (&table.buffer[Nx*i], table.header);
    fprintf (stdout, "%s\n", line);
    free (line);
  }    
  exit (0);
}

/* print an ASCII table to a row with single spaces separating value */
char *print_table_row (char *row, Header *header) {
  
  int i, j, Nx, Nfields, Nbytes, Nvals, Oout, Oin;
  char field[16], type[16], format[16], *line;

  gfits_scan (header, "NAXIS1",  "%d", 1, &Nx);
  gfits_scan (header, "TFIELDS", "%d", 1, &Nfields);

  /* assume we have one space per byte column */
  ALLOCATE (line, char, 2*Nx + 1);

  Oin = Oout = 0;
  for (i = 1; i <= Nfields; i++) {
    sprintf (field, "TFORM%d", i);
    gfits_scan (header, field, "%s", 1, format); /* get field format */
    gfits_table_format (format, type, &Nvals, &Nbytes);    /* convert to c-style */
    memcpy (&line[Oout], &row[Oin], Nvals*Nbytes);
    for (j = 0; j < Nvals*Nbytes; j++) if (line[Oout+j] == 0) line[Oout+j] = ' ';
    line[Oout+Nvals*Nbytes] = ' ';
    Oout += Nvals*Nbytes + 1;
    Oin += Nvals*Nbytes;
  }

  return (line);
}

void usage () {
    fprintf (stderr, "USAGE: [-h] [-N N] (table.fits)\n");
    fprintf (stderr, " -x N:       operate on extension N (0 is default)\n");
    fprintf (stderr, " -n EXTNAME: operate on named extension (incompatible with -x)\n");
    fprintf (stderr, " -row N:     print row number N\n");
    fprintf (stderr, " -column n:  print column named n\n");
    fprintf (stderr, " -ncolumn N: print column number N\n");
    fprintf (stderr, " -list:      print extension names\n");
    exit (2);
}

void list_extnames (char *file) {

  FILE *f;
  int i, Naxis, Nelem, Ncomp, extend, Nbytes, status;
  char extname[82], exttype[82], axisname[32];
  Header header;

  f = fopen (file, "r");
  if (f == (FILE *) NULL) {
    fprintf (stderr, "can't open file %s\n", file);
    exit (1);
  }

  /* 
     if (!gfits_fread_header (f, &header)) {
     fprintf (stderr, "can't read header from %s\n", file);
     exit (1);
     }
     Nbytes = gfits_matrix_size (&header);
     fseek (f, Nbytes, SEEK_CUR);
  */

  fprintf (stdout, "%-30s %-15s NAXIS NAXIS(i)...\n", "extname", "datatype"); 

  Ncomp = 0;
  while (gfits_fread_header (f, &header)) {
    /* extract the EXTNAME for this component (set to PHU for 0th component) */
    status = gfits_scan (&header, "EXTNAME", "%s", 1, extname);
    if (!status) {
      if (Ncomp == 0) {
	strcpy (extname, "PHU");
      } else {
	strcpy (extname, "UNKNOWN");
      }
    }
    fprintf (stdout, "%-30s ", extname);

    /* extract the datatype for this component (IMAGE for 0th component) */
    if (Ncomp == 0) {
      strcpy (exttype, "IMAGE");
    } else {
      status = gfits_scan (&header, "XTENSION", "%s", 1, exttype);
      if (!status) {
	strcpy (exttype, "UNKNOWN");
      }
    }
    fprintf (stdout, "%-15s ", exttype);

    /* extract the rank of the component */
    status = gfits_scan (&header, "NAXIS",  "%d", 1, &Naxis);
    if (!status) {
      fprintf (stderr, "component %d is missing Naxis!\n", Ncomp);
      Ncomp ++;
      continue;
    }
    fprintf (stdout, " %4d", Naxis);

    /* extract the individual axes */
    for (i = 0; i < Naxis; i++) {
      sprintf (axisname, "NAXIS%d", i+1);
      status = gfits_scan (&header, axisname,  "%d", 1, &Nelem);
      if (!status) {
	fprintf (stderr, "missing %s\n", axisname);
      }
      fprintf (stdout, " %7d", Nelem);
    }
    fprintf (stdout, "\n");

    /* are extensions identified? (we will scan for them anyway) */
    if (Ncomp == 0) {
      extend = FALSE;
      gfits_scan (&header, "EXTEND", "%t", 1, &extend);
      if (!extend) {
	fprintf (stderr, "no extensions listed in file\n");
      }
    }

    Nbytes = gfits_matrix_size (&header);
    fseek (f, Nbytes, SEEK_CUR);

    Ncomp ++;
  }
  fclose (f);
  exit (0);
}

FILE *load_extension (char *file, int Nextend, char *Extname, Header *header) {

  int i, extend, Nbytes;
  char extname[82];
  FILE *f;

  /* open file */
  f = fopen (file, "r");
  if (f == (FILE *) NULL) {
    fprintf (stderr, "can't open file %s\n", file);
    exit (1);
  }

  /* read PHU */
  if (!gfits_fread_header (f, header)) {
    fprintf (stderr, "can't read header from %s\n", file);
    exit (1);
  }

  /* check for existence of extensions */
  extend = FALSE;
  gfits_scan (header, "EXTEND", "%t", 1, &extend);
  if (!extend) {
    fprintf (stderr, "no extensions listed in file\n");
  }

  /* skip first data array */
  Nbytes = gfits_matrix_size (header);
  fseek (f, Nbytes, SEEK_CUR);

  /* search for extension of interest */
  for (i = 0; gfits_fread_Theader (f, header); i++) {
    if ((Extname == NULL) && (Nextend == i)) return (f);

    gfits_scan (header, "EXTNAME", "%s", 1, extname);
    if ((Extname != NULL) && (!strcmp (Extname, extname))) return (f);

    Nbytes = gfits_matrix_size (header);
    fseek (f, Nbytes, SEEK_CUR);
    gfits_free_header (header);
  }
  fclose (f);

  fprintf (stderr, "failed to load extension of interest\n");
  exit (1);
}

void print_column (FTable *table, int Column, char *Colname) {
  
  int i, j, Nfields, Nstart, Nv, Nb, Nx, Ny;
  Header *header;
  char format[16], field[16], type[16], *line, *data;
  int *Ti;
  float *Tf;

  header =  table[0].header;
  data   =  table[0].buffer;

  gfits_scan (header, "TFIELDS", "%d", 1, &Nfields);
  gfits_scan (header, "NAXIS1",  "%d", 1, &Nx);
  gfits_scan (header, "NAXIS2",  "%d", 1, &Ny);

  if (Colname != (char *) NULL) {
    /* find matching column entry */
    for (i = 1; i <= Nfields; i++) {
      sprintf (field, "TTYPE%d", i);
      gfits_scan (header, field, "%s", 1, type);
      if (!strcmp (type, Colname)) {
	Column = i;
	break;
      }
    }
    if (!Column) {
      fprintf (stderr, "column %s not found\n", Colname);
      exit (1);
    }
  } else {
    /* check if Column is in range */
    if (Column > Nfields) {
      fprintf (stderr, "-ncolumn %d too large, only %d columns\n", Column, Nfields);
      exit (1);
    }
  }

  /* scan columns to find insert point */
  Nstart = 0;
  for (i = 1; i < Column; i++) {
    sprintf (field, "TFORM%d", i);
    gfits_scan (header, field, "%s", 1, format);
    if (Binary) 
      gfits_table_format (format, type, &Nv, &Nb);
    else 
      gfits_bintable_format (format, type, &Nv, &Nb);
    Nstart += Nv*Nb;
  }

  sprintf (field, "TFORM%d", Column);
  gfits_scan (header, field, "%s", 1, format);
  if (Binary) 
    gfits_bintable_format (format, type, &Nv, &Nb);    /* convert to c-style */
  else 
    gfits_table_format (format, type, &Nv, &Nb);    /* convert to c-style */

  ALLOCATE (line, char, Nv*Nb + 1);

  if (Binary) {
    if (!gfits_get_bintable_column_type (header, Colname, type, &Nv)) return;
    if (!strcmp (type, "char")) return;
    if (!gfits_get_bintable_column (header, table, Colname, (void **)&data)) return;

    for (i = 0; i < Ny; i++) {
      if (!strcmp (type, "char")) {
	memcpy (line, &data[i*Nx + Nstart], Nv*Nb);
	fprintf (stdout, "%s\n", line);
      } else {
	for (j = 0; j < Nv; j++) {
	  if (!strcmp (type, "int")) {
	    memcpy (line, &data[i*Nx + Nstart + Nb*j], Nb);
	    ByteSwap (line, Nb, 1, "int");
	    Ti = (int *)line;
	    fprintf (stdout, "%d ", Ti[0]);
	  }
	  if (!strcmp (type, "float")) {
	    memcpy (line, &data[i*Nx + Nstart + Nb*j], Nb);
	    ByteSwap (line, Nb, 1, "float");
	    Tf = (float *)line;
	    fprintf (stdout, "%f ", Tf[0]);
	  }
	}
	fprintf (stdout, "\n");
      }
    }
  } else {
    for (i = 0; i < Ny; i++) {
      memcpy (line, &data[i*Nx + Nstart], Nv*Nb);
      sprintf (format, "%%%ds\n", Nv*Nb);
      fprintf (stdout, format, line);
    }
  }
  exit (0);
}

void print_layout (Header *header) {

  int i, Nfields;
  char field[16], type[80], comment[80], format[80], unit[80], null[80], nval[80];

  gfits_scan (header, "TFIELDS", "%d", 1, &Nfields);

  for (i = 1; i <= Nfields; i++) {
    sprintf (field, "TTYPE%d", i);
    gfits_scan (header, field, "%s", 1, type);
    gfits_scan (header, field, "%C", 1, comment);
    sprintf (field, "TFORM%d", i);
    gfits_scan (header, field, "%s", 1, format);
    sprintf (field, "TUNIT%d", i);
    gfits_scan (header, field, "%s", 1, unit);
    sprintf (field, "TNULL%d", i);
    if (!gfits_scan (header, field, "%s", 1, null)) strcpy (null, "none");
    sprintf (field, "TNVAL%d", i);
    if (!gfits_scan (header, field, "%s", 1, nval)) strcpy (nval, "none");

    fprintf (stdout, "%-18s %-32s %-18s %-6s %-8s %-8s\n", 
	  type, comment, unit, format, null, nval); 
  }
  exit (0);
}
