# include <ohana.h>
# include <gfitsio.h>
#include <inttypes.h>
# define SWAP_BYTE { \
  char tmp; \
  tmp = Pin[0]; Pin[0] = Pin[1]; Pin[1] = tmp; }
# define SWAP_WORD { \
  char tmp; \
  tmp = Pin[0]; Pin[0] = Pin[3]; Pin[3] = tmp; \
  tmp = Pin[1]; Pin[1] = Pin[2]; Pin[2] = tmp; }
# define SWAP_DBLE { \
  char tmp; \
  tmp = Pin[0]; Pin[0] = Pin[7]; Pin[7] = tmp; \
  tmp = Pin[1]; Pin[1] = Pin[6]; Pin[6] = tmp; \
  tmp = Pin[2]; Pin[2] = Pin[5]; Pin[5] = tmp; \
  tmp = Pin[3]; Pin[3] = Pin[4]; Pin[4] = tmp; }

void *gfits_get_bintable_column_data (Header *header, FTable *table, char *label, char *type, off_t *Nrow, int *Ncol) {

  off_t Nx, Ny;
  int i, N, Nfields, Nval, Nbytes, Nstart, Nv, Nb;
  char tlabel[80], field[80], format[80], tmpline[16];
  char *Pin, *Pout, *array;
  double Bscale, Bzero;

  if (label == (char *) NULL) return (NULL);
  if (label[0] == 0) return (NULL);

  /* find label in header */
  tlabel[0] = 0;
  if (!gfits_scan (header, "TFIELDS", "%d", 1, &Nfields)) return (NULL);
  for (i = 1; strcmp (label, tlabel) && (i < Nfields + 1); i++) {
    sprintf (field, "TTYPE%d", i);
    gfits_scan (header, field, "%s", 1, tlabel);
  }
  if (strcmp (label, tlabel)) return (NULL);
  N = i - 1;

  Bscale = 1; 
  Bzero  = 0;

  /* interpret format */
  sprintf (field, "TSCAL%d", N);
  gfits_scan (header, field, "%lf", 1, &Bscale);
  sprintf (field, "TZERO%d", N);
  gfits_scan (header, field, "%lf", 1, &Bzero);
  sprintf (field, "TFORM%d", N);
  gfits_scan (header, field, "%s", 1, format);

  if (!gfits_bintable_format (format, type, &Nval, &Nbytes)) return (NULL);
  
  /* check existing table dimensions */
  gfits_scan (header, "NAXIS1",  "%lld", 1, (long long *) &Nx);
  gfits_scan (header, "NAXIS2",  "%lld", 1, (long long *) &Ny);

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

  /* extract bytes from table into array */
  ALLOCATE (array, char, Nbytes*Nval*Ny);
  Pin  = table[0].buffer + Nstart;
  Pout = array;
  for (i = 0; i < Ny; i++, Pin += Nx, Pout += Nval*Nbytes) {
    memcpy (Pout, Pin, Nval*Nbytes);
  }

  /* convert data in-situ with correct type, byte swap and Bzero/Bscale */
  Pin  = array;
  Pout = array;
  if (!strcmp (type, "char")) {
    for (i = 0; i < Nval*Ny; i++, Pin+=Nbytes, Pout+=Nbytes) {
      *(char *)Pout = *(char *)Pin*Bscale + Bzero;
    }
  }
  if (!strcmp (type, "short")) {
    for (i = 0; i < Nval*Ny; i++, Pin+=Nbytes, Pout+=Nbytes) {
# ifdef BYTE_SWAP
      SWAP_BYTE;
# endif
      *(short *)Pout = *(short *)Pin*Bscale + Bzero;
    }  
  }
  if (!strcmp (type, "int")) {
    for (i = 0; i < Nval*Ny; i++, Pin+=Nbytes, Pout+=Nbytes) {
# ifdef BYTE_SWAP
      SWAP_WORD;
# endif
      *(int *)Pout = *(int *)Pin*Bscale + Bzero;
    }
  }
  if (!strcmp (type, "int64_t")) {
    for (i = 0; i < Nval*Ny; i++, Pin+=Nbytes, Pout+=Nbytes) {
# ifdef BYTE_SWAP
      SWAP_DBLE;
# endif
      *(int64_t *)Pout = *(int64_t *)Pin*Bscale + Bzero;
    }
  }
  if (!strcmp (type, "float")) {
    for (i = 0; i < Nval*Ny; i++, Pin+=Nbytes, Pout+=Nbytes) {
# ifdef BYTE_SWAP
      SWAP_WORD;
# endif
      *(float *)Pout = *(float *)Pin*Bscale + Bzero;
    }
  }
  if (!strcmp (type, "double")) {
    for (i = 0; i < Nval*Ny; i++, Pin+=Nbytes, Pout+=Nbytes) {
# ifdef BYTE_SWAP
      SWAP_DBLE;
# endif
      *(double *)Pout = *(double *)Pin*Bscale + Bzero;
    }
  }

  *Ncol = Nval;
  *Nrow = Ny;
  return (array);
}

/***********************/
int gfits_get_bintable_column_type (Header *header, char *label, char *type, int *Nval) {

  int i, N, Nfields, Nbytes;
  char tlabel[80], field[80], format[80];

  if (label == (char *) NULL) return (FALSE);
  if (label[0] == 0) return (FALSE);

  /* find label in header */
  tlabel[0] = 0;
  gfits_scan (header, "TFIELDS", "%d", 1, &Nfields);
  for (i = 1; strcmp (label, tlabel) && (i < Nfields + 1); i++) {
    sprintf (field, "TTYPE%d", i);
    gfits_scan (header, field, "%s", 1, tlabel);
  }
  if (strcmp (label, tlabel)) return (FALSE);
  N = i - 1;

  sprintf (field, "TFORM%d", N);
  gfits_scan (header, field, "%s", 1, format);

  if (!gfits_bintable_format (format, type, Nval, &Nbytes)) return (FALSE);
  return (TRUE);
}

/***********************/
int gfits_get_bintable_column (Header *header, FTable *table, char *label, void **data) {

  off_t Nx, Ny;
  int i, N, Nfields, Nval, Nbytes, Nstart, Nv, Nb;
  char tlabel[80], field[80], format[80], type[16], tmpline[16];
  char *Pin, *Pout, *array;
  double Bscale, Bzero;

  if (label == (char *) NULL) return (FALSE);
  if (label[0] == 0) return (FALSE);

  /* find label in header */
  tlabel[0] = 0;
  gfits_scan (header, "TFIELDS", "%d", 1, &Nfields);
  for (i = 1; strcmp (label, tlabel) && (i < Nfields + 1); i++) {
    sprintf (field, "TTYPE%d", i);
    gfits_scan (header, field, "%s", 1, tlabel);
  }
  if (strcmp (label, tlabel)) return (FALSE);
  N = i - 1;

  Bscale = 1; 
  Bzero  = 0;

  /* interpret format */
  sprintf (field, "TSCAL%d", N);
  gfits_scan (header, field, "%lf", 1, &Bscale);
  sprintf (field, "TZERO%d", N);
  gfits_scan (header, field, "%lf", 1, &Bzero);
  sprintf (field, "TFORM%d", N);
  gfits_scan (header, field, "%s", 1, format);

  if (!gfits_bintable_format (format, type, &Nval, &Nbytes)) return (FALSE);
  
  /* check existing table dimensions */
  gfits_scan (header, "NAXIS1",  "%lld", 1, (long long *) &Nx);
  gfits_scan (header, "NAXIS2",  "%lld", 1, (long long *) &Ny);

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

  /* extract bytes from table into temporary array */
  ALLOCATE (array, char, Nbytes*Nval*Ny);
  Pin  = table[0].buffer + Nstart;
  Pout = array;
  for (i = 0; i < Ny; i++, Pin += Nx, Pout += Nval*Nbytes) {
    memcpy (Pout, Pin, Nval*Nbytes);
  }

  /* convert data in-situ with correct type, byte swap and Bzero/Bscale */
  Pin  = array;
  Pout = array;
  if (!strcmp (type, "char")) {
    for (i = 0; i < Nval*Ny; i++, Pin+=Nbytes, Pout+=Nbytes) {
      *(char *)Pout = *(char *)Pin*Bscale + Bzero;
    }
  }
  if (!strcmp (type, "short")) {
    for (i = 0; i < Nval*Ny; i++, Pin+=Nbytes, Pout+=Nbytes) {
# ifdef BYTE_SWAP
      SWAP_BYTE;
# endif
      *(short *)Pout = *(short *)Pin*Bscale + Bzero;
    }  
  }
  if (!strcmp (type, "int")) {
    for (i = 0; i < Nval*Ny; i++, Pin+=Nbytes, Pout+=Nbytes) {
# ifdef BYTE_SWAP
      SWAP_WORD;
# endif
      *(int *)Pout = *(int *)Pin*Bscale + Bzero;
    }
  }
  if (!strcmp (type, "int64_t")) {
    for (i = 0; i < Nval*Ny; i++, Pin+=Nbytes, Pout+=Nbytes) {
# ifdef BYTE_SWAP
      SWAP_DBLE;
# endif
      *(int64_t *)Pout = *(int64_t *)Pin*Bscale + Bzero;
    }
  }
  if (!strcmp (type, "float")) {
    for (i = 0; i < Nval*Ny; i++, Pin+=Nbytes, Pout+=Nbytes) {
# ifdef BYTE_SWAP
      SWAP_WORD;
# endif
      *(float *)Pout = *(float *)Pin*Bscale + Bzero;
    }
  }
  if (!strcmp (type, "double")) {
    for (i = 0; i < Nval*Ny; i++, Pin+=Nbytes, Pout+=Nbytes) {
# ifdef BYTE_SWAP
      SWAP_DBLE;
# endif
      *(double *)Pout = *(double *)Pin*Bscale + Bzero;
    }
  }

  *data = array;
  return (TRUE);
}

/* 
   valid BINTABLE column formats:
   L - logical
   X - bit
   I - 16 bit int
   J - 32 bit int
   A - char
   E - float 
   D - double
   B - unsigned bytes
   C - complex float
   M - complex double
   P - var lenght array descrpt (64 bit)
   
   all can be preceeded by integer N which specified number
   of entries in array
   
*/

/***********************/
int gfits_get_table_column_type (Header *header, char *label, char *type) {

  int i, N, Nfields, Nval, Nbytes;
  char tlabel[80], field[80], format[80];

  if (label == (char *) NULL) return (FALSE);
  if (label[0] == 0) return (FALSE);

  /* find label in header */
  tlabel[0] = 0;
  gfits_scan (header, "TFIELDS", "%d", 1, &Nfields);
  for (i = 1; strcmp (label, tlabel) && (i < Nfields + 1); i++) {
    sprintf (field, "TTYPE%d", i);
    gfits_scan (header, field, "%s", 1, tlabel);
  }
  if (strcmp (label, tlabel)) return (FALSE);
  N = i - 1;

  /* interpret format */
  sprintf (field, "TFORM%d", N);
  gfits_scan (header, field, "%s", 1, format);

  if (!gfits_table_format (format, type, &Nval, &Nbytes)) return (FALSE);
  return (TRUE);
}

/***********************/
int gfits_get_table_column (Header *header, FTable *table, char *label, void **data) {

  off_t Nx, Ny;
  int i, N, Nfields, Nval, Nbytes, Nstart, Nv, Nb;
  char tlabel[80], field[80], format[80], cformat[80], type[16], tmp[16];
  char *array, *Pin, *Pout, *line;

  if (label == (char *) NULL) return (FALSE);
  if (label[0] == 0) return (FALSE);
  array = NULL;

  /* find label in header */
  tlabel[0] = 0;
  gfits_scan (header, "TFIELDS", "%d", 1, &Nfields);
  for (i = 1; strcmp (label, tlabel) && (i < Nfields + 1); i++) {
    sprintf (field, "TTYPE%d", i);
    gfits_scan (header, field, "%s", 1, tlabel);
  }
  if (strcmp (label, tlabel)) return (FALSE);
  N = i - 1;

  /* interpret format */
  sprintf (field, "TFORM%d", N);
  gfits_scan (header, field, "%s", 1, format);

  if (!gfits_table_format (format, type, &Nval, &Nbytes)) return (FALSE);
  strcpy (cformat, format);
  
  /* check existing table dimensions */
  gfits_scan (header, "NAXIS1",  "%lld", 1, (long long *) &Nx);
  gfits_scan (header, "NAXIS2",  "%lld", 1, (long long *) &Ny);

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

  /* allocate temporary line, init pointers */
  ALLOCATE (line, char, Nval+1);
  bzero (line, Nval+1);
  Pin  = table[0].buffer + Nstart;

  /* allocate output array, copy/scan line as needed */
  if (!strcmp (type, "char")) {
    ALLOCATE (array, char, Ny*Nval);
    Pout = array;
    for (i = 0; i < Ny; i++, Pin+=Nx, Pout+=Nval*Nbytes) {
      memcpy (Pout, Pin, Nval*Nbytes);
    }
  }
  if (!strcmp (type, "int")) {
    ALLOCATE (array, char, Ny*4);
    Pout = array;
    for (i = 0; i < Ny; i++, Pin+=Nx, Pout+=4) {
      memcpy (line, Pin, Nval*Nbytes);
      sscanf (line, "%d", (int *)Pout);
    }
  }
  if (!strcmp (type, "int64_t")) {
    ALLOCATE (array, char, Ny*8);
    Pout = array;
    for (i = 0; i < Ny; i++, Pin+=Nx, Pout+=8) {
      memcpy (line, Pin, Nval*Nbytes);
      sscanf (line, "%" PRId64, (int64_t *)Pout);
    }
  }
  if (!strcmp (type, "float")) {
    ALLOCATE (array, char, Ny*4);
    Pout = array;
    for (i = 0; i < Ny; i++, Pin+=Nx, Pout+=4) {
      memcpy (line, Pin, Nval*Nbytes);
      sscanf (line, "%f", (float *)Pout);
    }
  }
  *data = array;
  free (line);

  return (TRUE);
}

/*
  valid TABLE column formats:
  FN.N  - floating point
  INN   - integer
  ANN   - NN char string
  
*/

