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

/***********************/
int gfits_set_bintable_column (Header *header, FTable *table, char *label, void *data, off_t Nrow) {

  off_t Nx, Ny, nbytes;
  int i, N, Nfields;
  int 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; strcasecmp (label, tlabel) && (i < Nfields + 1); i++) {
    sprintf (field, "TTYPE%d", i);
    gfits_scan (header, field, "%s", 1, tlabel);
  }
  if (strcasecmp (label, tlabel)) return (FALSE);
  N = i - 1;

  /* interpret format */
  sprintf (field, "TSCAL%d", N);
  if (!gfits_scan (header, field, "%lf", 1, &Bscale)) {
    Bscale = 1.0;
  }
  sprintf (field, "TZERO%d", N);
  if (!gfits_scan (header, field, "%lf", 1, &Bzero)) {
    Bzero = 0.0;
  }
  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", OFF_T_FMT, 1,  &Nx);
  gfits_scan (header, "NAXIS2", OFF_T_FMT, 1,  &Ny);
  if (Ny == 0) { 
    Ny = Nrow;
    header[0].Naxis[1] = Ny;
    gfits_modify (header, "NAXIS2", OFF_T_FMT, 1,  Ny);

    nbytes = gfits_data_size (header);
    REALLOCATE (table[0].buffer, char, MAX (nbytes, 1));
    bzero (table[0].buffer, nbytes);
    table[0].datasize = nbytes;
  }
  if (Ny != Nrow) return (FALSE);

  /* 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;
  }

  /* make duplicate of data with correct type
     byte swap and Bzero/Bscale */
  ALLOCATE (array, char, Nbytes*Nval*Nrow);
  Pin = data;
  Pout = array;
  // does it makes sense to scale 'char' data?
  if (!strcmp (type, "char")) {
    for (i = 0; i < Nval*Nrow; i++, Pin+=Nbytes, Pout+=Nbytes) {
      *(char *)Pout = (*(char *)Pin - Bzero) / Bscale;
    }
  }
  if (!strcmp (type, "byte")) {
    for (i = 0; i < Nval*Nrow; i++, Pin+=Nbytes, Pout+=Nbytes) {
      *(char *)Pout = (*(char *)Pin - Bzero) / Bscale;
    }
  }
  if (!strcmp (type, "short")) {
    for (i = 0; i < Nval*Nrow; i++, Pin+=Nbytes, Pout+=Nbytes) {
      *(short *)Pout = (*(short *)Pin - Bzero) / Bscale;
# ifdef BYTE_SWAP
      SWAP_BYTE;
# endif
    }  
  }
  if (!strcmp (type, "int")) {
    for (i = 0; i < Nval*Nrow; i++, Pin+=Nbytes, Pout+=Nbytes) {
      *(int *)Pout = (*(int *)Pin - Bzero) / Bscale;
# ifdef BYTE_SWAP
      SWAP_WORD;
# endif
    }
  }
  if (!strcmp (type, "int64_t")) {
    // XXX 64 bit int operations with Bzero & Bscale are inaccurate even with doubles
    if ((Bzero == 0.0) && (Bscale == 1.0)) {
      for (i = 0; i < Nval*Nrow; i++, Pin+=Nbytes, Pout+=Nbytes) {
	*(int64_t *)Pout = *(int64_t *)Pin;
# ifdef BYTE_SWAP
	SWAP_DBLE;
# endif
      }
    } else {
      for (i = 0; i < Nval*Nrow; i++, Pin+=Nbytes, Pout+=Nbytes) {
	*(int64_t *)Pout = (*(int64_t *)Pin - Bzero) / Bscale;
# ifdef BYTE_SWAP
	SWAP_DBLE;
# endif
      }
    }
  }
  if (!strcmp (type, "float")) {
    for (i = 0; i < Nval*Nrow; i++, Pin+=Nbytes, Pout+=Nbytes) {
      *(float *)Pout = (*(float *)Pin - Bzero) / Bscale;
# ifdef BYTE_SWAP
      SWAP_WORD;
# endif
    }
  }
  if (!strcmp (type, "double")) {
    for (i = 0; i < Nval*Nrow; i++, Pin+=Nbytes, Pout+=Nbytes) {
      *(double *)Pout = (*(double *)Pin - Bzero) / Bscale;
# ifdef BYTE_SWAP
      SWAP_DBLE;
# endif
    }
  }

  /* check array space */
  if (Nx*Ny < Nx*(Nrow - 1) + Nstart + Nval*Nbytes) {
    fprintf (stderr, "mismatch in array sizes\n");
    return (FALSE);
  }

  /* insert bytes from array into appropriate section of buffer */
  Pout = table[0].buffer + Nstart;
  Pin  = array;
  for (i = 0; i < Nrow; i++, Pout += Nx, Pin += Nval*Nbytes) {
    memcpy (Pout, Pin, Nval*Nbytes);
  }

  free (array);
  return (TRUE);
}

/***********************/
int gfits_set_bintable_column_reformat (Header *header, FTable *table, char *label, char *intype, void *data, off_t Nrow) {

  off_t Nx, Ny, nbytes;
  int i, N, Nfields;
  int Nval, NbytesOut, Nstart, Nv, Nb;
  char tlabel[80], field[80], format[80], outtype[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; strcasecmp (label, tlabel) && (i < Nfields + 1); i++) {
    sprintf (field, "TTYPE%d", i);
    gfits_scan (header, field, "%s", 1, tlabel);
  }
  if (strcasecmp (label, tlabel)) return (FALSE);
  N = i - 1;

  /* interpret format */
  sprintf (field, "TSCAL%d", N);
  if (!gfits_scan (header, field, "%lf", 1, &Bscale)) {
    Bscale = 1.0;
  }
  sprintf (field, "TZERO%d", N);
  if (!gfits_scan (header, field, "%lf", 1, &Bzero)) {
    Bzero = 0.0;
  }
  sprintf (field, "TFORM%d", N);
  gfits_scan (header, field, "%s", 1, format);

  if (!gfits_bintable_format (format, outtype, &Nval, &NbytesOut)) return (FALSE);
  
  /* check existing table dimensions */
  gfits_scan (header, "NAXIS1", OFF_T_FMT, 1,  &Nx);
  gfits_scan (header, "NAXIS2", OFF_T_FMT, 1,  &Ny);
  if (Ny == 0) { 
    Ny = Nrow;
    header[0].Naxis[1] = Ny;
    gfits_modify (header, "NAXIS2", OFF_T_FMT, 1,  Ny);

    nbytes = gfits_data_size (header);
    REALLOCATE (table[0].buffer, char, MAX (nbytes, 1));
    bzero (table[0].buffer, nbytes);
    table[0].datasize = nbytes;
  }
  if (Ny != Nrow) return (FALSE);

  /* 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;
  }

  /* make duplicate of data with correct type
     byte swap and Bzero/Bscale */
  ALLOCATE (array, char, NbytesOut*Nval*Nrow);
  Pin = data;
  Pout = array;

  // # define ASSIGN_DATA(ITYPE,INAME,ISIZE,OTYPE,ONAME)
  // if (!strcmp (intype, #ITYPE)) {

  /** input == char **/
  if (!strcmp (outtype, "char") && !strcmp (intype, "char")) {
    int NbytesIn = 1;
    for (i = 0; i < Nval*Nrow; i++, Pin+=NbytesIn, Pout+=NbytesOut) {
      *(char *)Pout = (*(char *)Pin - Bzero) / Bscale;
    }
  }
  if (!strcmp (outtype, "byte") && !strcmp (intype, "char")) {
    int NbytesIn = 1;
    for (i = 0; i < Nval*Nrow; i++, Pin+=NbytesIn, Pout+=NbytesOut) {
      *(char *)Pout = (*(char *)Pin - Bzero) / Bscale;
    }
  }
  if (!strcmp (outtype, "short") && !strcmp (intype, "char")) {
    int NbytesIn = 1;
    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
      *(short *)Pout = (*(char *)Pin - Bzero) / Bscale;
# ifdef BYTE_SWAP
      SWAP_BYTE;
# endif
    }  
  }
  if (!strcmp (outtype, "int") && !strcmp (intype, "char")) {
    int NbytesIn = 1;
    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
      *(int *)Pout = (*(char *)Pin - Bzero) / Bscale;
# ifdef BYTE_SWAP
      SWAP_WORD;
# endif
    }
  }
  if (!strcmp (outtype, "int64_t") && !strcmp (intype, "char")) {
    int NbytesIn = 1;
    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
      *(int64_t *)Pout = (*(char *)Pin - Bzero) / Bscale;
# ifdef BYTE_SWAP
      SWAP_DBLE;
# endif
    }
  }
  if (!strcmp (outtype, "float") && !strcmp (intype, "char")) {
    int NbytesIn = 1;
    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
      *(float *)Pout = (*(char *)Pin - Bzero) / Bscale;
# ifdef BYTE_SWAP
      SWAP_WORD;
# endif
    }
  }
  if (!strcmp (outtype, "double") && !strcmp (intype, "char")) {
    int NbytesIn = 1;
    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
      *(double *)Pout = (*(char *)Pin - Bzero) / Bscale;
# ifdef BYTE_SWAP
      SWAP_DBLE;
# endif
    }
  }

  /** input == byte **/
  if (!strcmp (outtype, "char") && !strcmp (intype, "byte")) {
    int NbytesIn = 1;
    for (i = 0; i < Nval*Nrow; i++, Pin+=NbytesIn, Pout+=NbytesOut) {
      *(char *)Pout = (*(char *)Pin - Bzero) / Bscale;
    }
  }
  if (!strcmp (outtype, "byte") && !strcmp (intype, "byte")) {
    int NbytesIn = 1;
    for (i = 0; i < Nval*Nrow; i++, Pin+=NbytesIn, Pout+=NbytesOut) {
      *(char *)Pout = (*(char *)Pin - Bzero) / Bscale;
    }
  }
  if (!strcmp (outtype, "short") && !strcmp (intype, "byte")) {
    int NbytesIn = 1;
    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
      *(short *)Pout = (*(char *)Pin - Bzero) / Bscale;
# ifdef BYTE_SWAP
      SWAP_BYTE;
# endif
    }  
  }
  if (!strcmp (outtype, "int") && !strcmp (intype, "byte")) {
    int NbytesIn = 1;
    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
      *(int *)Pout = (*(char *)Pin - Bzero) / Bscale;
# ifdef BYTE_SWAP
      SWAP_WORD;
# endif
    }
  }
  if (!strcmp (outtype, "int64_t") && !strcmp (intype, "byte")) {
    int NbytesIn = 1;
    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
      *(int64_t *)Pout = (*(char *)Pin - Bzero) / Bscale;
# ifdef BYTE_SWAP
      SWAP_DBLE;
# endif
    }
  }
  if (!strcmp (outtype, "float") && !strcmp (intype, "byte")) {
    int NbytesIn = 1;
    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
      *(float *)Pout = (*(char *)Pin - Bzero) / Bscale;
# ifdef BYTE_SWAP
      SWAP_WORD;
# endif
    }
  }
  if (!strcmp (outtype, "double") && !strcmp (intype, "byte")) {
    int NbytesIn = 1;
    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
      *(double *)Pout = (*(char *)Pin - Bzero) / Bscale;
# ifdef BYTE_SWAP
      SWAP_DBLE;
# endif
    }
  }

  /** input == short **/
  if (!strcmp (outtype, "char") && !strcmp (intype, "short")) {
    int NbytesIn = 2;
    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
      *(char *)Pout = (*(short *)Pin - Bzero) / Bscale;
    }
  }
  if (!strcmp (outtype, "byte") && !strcmp (intype, "short")) {
    int NbytesIn = 2;
    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
      *(char *)Pout = (*(short *)Pin - Bzero) / Bscale;
    }
  }
  if (!strcmp (outtype, "short") && !strcmp (intype, "short")) {
    int NbytesIn = 2;
    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
      *(short *)Pout = (*(short *)Pin - Bzero) / Bscale;
# ifdef BYTE_SWAP
      SWAP_BYTE;
# endif
    }  
  }
  if (!strcmp (outtype, "int") && !strcmp (intype, "short")) {
    int NbytesIn = 2;
    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
      *(int *)Pout = (*(short *)Pin - Bzero) / Bscale;
# ifdef BYTE_SWAP
      SWAP_WORD;
# endif
    }
  }
  if (!strcmp (outtype, "int64_t") && !strcmp (intype, "short")) {
    int NbytesIn = 2;
    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
      *(int64_t *)Pout = (*(short *)Pin - Bzero) / Bscale;
# ifdef BYTE_SWAP
      SWAP_DBLE;
# endif
    }
  }
  if (!strcmp (outtype, "float") && !strcmp (intype, "short")) {
    int NbytesIn = 2;
    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
      *(float *)Pout = (*(short *)Pin - Bzero) / Bscale;
# ifdef BYTE_SWAP
      SWAP_WORD;
# endif
    }
  }
  if (!strcmp (outtype, "double") && !strcmp (intype, "short")) {
    int NbytesIn = 2;
    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
      *(double *)Pout = (*(short *)Pin - Bzero) / Bscale;
# ifdef BYTE_SWAP
      SWAP_DBLE;
# endif
    }
  }

  /** input == int **/
  if (!strcmp (outtype, "char") && !strcmp (intype, "int")) {
    int NbytesIn = 4;
    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
      *(char *)Pout = (*(int *)Pin - Bzero) / Bscale;
    }
  }
  if (!strcmp (outtype, "byte") && !strcmp (intype, "int")) {
    int NbytesIn = 4;
    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
      *(char *)Pout = (*(int *)Pin - Bzero) / Bscale;
    }
  }
  if (!strcmp (outtype, "short") && !strcmp (intype, "int")) {
    int NbytesIn = 4;
    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
      *(short *)Pout = (*(int *)Pin - Bzero) / Bscale;
# ifdef BYTE_SWAP
      SWAP_BYTE;
# endif
    }  
  }
  if (!strcmp (outtype, "int") && !strcmp (intype, "int")) {
    int NbytesIn = 4;
    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
      *(int *)Pout = (*(int *)Pin - Bzero) / Bscale;
# ifdef BYTE_SWAP
      SWAP_WORD;
# endif
    }
  }
  if (!strcmp (outtype, "int64_t") && !strcmp (intype, "int")) {
    int NbytesIn = 4;
    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
      *(int64_t *)Pout = (*(int *)Pin - Bzero) / Bscale;
# ifdef BYTE_SWAP
      SWAP_DBLE;
# endif
    }
  }
  if (!strcmp (outtype, "float") && !strcmp (intype, "int")) {
    int NbytesIn = 4;
    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
      *(float *)Pout = (*(int *)Pin - Bzero) / Bscale;
# ifdef BYTE_SWAP
      SWAP_WORD;
# endif
    }
  }
  if (!strcmp (outtype, "double") && !strcmp (intype, "int")) {
    int NbytesIn = 4;
    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
      *(double *)Pout = (*(int *)Pin - Bzero) / Bscale;
# ifdef BYTE_SWAP
      SWAP_DBLE;
# endif
    }
  }

  /** input == int64_t **/
  if (!strcmp (outtype, "char") && !strcmp (intype, "int64_t")) {
    int NbytesIn = 8;
    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
      *(char *)Pout = (*(int64_t *)Pin - Bzero) / Bscale;
    }
  }
  if (!strcmp (outtype, "byte") && !strcmp (intype, "int64_t")) {
    int NbytesIn = 8;
    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
      *(char *)Pout = (*(int64_t *)Pin - Bzero) / Bscale;
    }
  }
  if (!strcmp (outtype, "short") && !strcmp (intype, "int64_t")) {
    int NbytesIn = 8;
    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
      *(short *)Pout = (*(int64_t *)Pin - Bzero) / Bscale;
# ifdef BYTE_SWAP
      SWAP_BYTE;
# endif
    }  
  }
  if (!strcmp (outtype, "int") && !strcmp (intype, "int64_t")) {
    int NbytesIn = 8;
    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
      *(int *)Pout = (*(int64_t *)Pin - Bzero) / Bscale;
# ifdef BYTE_SWAP
      SWAP_WORD;
# endif
    }
  }
  if (!strcmp (outtype, "int64_t") && !strcmp (intype, "int64_t")) {
    int NbytesIn = 8;
    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
      *(int64_t *)Pout = (*(int64_t *)Pin - Bzero) / Bscale;
# ifdef BYTE_SWAP
      SWAP_DBLE;
# endif
    }
  }
  if (!strcmp (outtype, "float") && !strcmp (intype, "int64_t")) {
    int NbytesIn = 8;
    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
      *(float *)Pout = (*(int64_t *)Pin - Bzero) / Bscale;
# ifdef BYTE_SWAP
      SWAP_WORD;
# endif
    }
  }
  if (!strcmp (outtype, "double") && !strcmp (intype, "int64_t")) {
    int NbytesIn = 8;
    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
      *(double *)Pout = (*(int64_t *)Pin - Bzero) / Bscale;
# ifdef BYTE_SWAP
      SWAP_DBLE;
# endif
    }
  }

  /** input == float **/
  if (!strcmp (outtype, "char") && !strcmp (intype, "float")) {
    int NbytesIn = 4;
    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
      *(char *)Pout = (*(float *)Pin - Bzero) / Bscale;
    }
  }
  if (!strcmp (outtype, "byte") && !strcmp (intype, "float")) {
    int NbytesIn = 4;
    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
      *(char *)Pout = (*(float *)Pin - Bzero) / Bscale;
    }
  }
  if (!strcmp (outtype, "short") && !strcmp (intype, "float")) {
    int NbytesIn = 4;
    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
      *(short *)Pout = (*(float *)Pin - Bzero) / Bscale;
# ifdef BYTE_SWAP
      SWAP_BYTE;
# endif
    }  
  }
  if (!strcmp (outtype, "int") && !strcmp (intype, "float")) {
    int NbytesIn = 4;
    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
      *(int *)Pout = (*(float *)Pin - Bzero) / Bscale;
# ifdef BYTE_SWAP
      SWAP_WORD;
# endif
    }
  }
  if (!strcmp (outtype, "int64_t") && !strcmp (intype, "float")) {
    int NbytesIn = 4;
    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
      *(int64_t *)Pout = (*(float *)Pin - Bzero) / Bscale;
# ifdef BYTE_SWAP
      SWAP_DBLE;
# endif
    }
  }
  if (!strcmp (outtype, "float") && !strcmp (intype, "float")) {
    int NbytesIn = 4;
    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
      *(float *)Pout = (*(float *)Pin - Bzero) / Bscale;
# ifdef BYTE_SWAP
      SWAP_WORD;
# endif
    }
  }
  if (!strcmp (outtype, "double") && !strcmp (intype, "float")) {
    int NbytesIn = 4;
    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
      *(double *)Pout = (*(float *)Pin - Bzero) / Bscale;
# ifdef BYTE_SWAP
      SWAP_DBLE;
# endif
    }
  }

  /** input == double **/
  if (!strcmp (outtype, "char") && !strcmp (intype, "double")) {
    int NbytesIn = 8;
    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
      *(char *)Pout = (*(double *)Pin - Bzero) / Bscale;
    }
  }
  if (!strcmp (outtype, "byte") && !strcmp (intype, "double")) {
    int NbytesIn = 8;
    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
      *(char *)Pout = (*(double *)Pin - Bzero) / Bscale;
    }
  }
  if (!strcmp (outtype, "short") && !strcmp (intype, "double")) {
    int NbytesIn = 8;
    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
      *(short *)Pout = (*(double *)Pin - Bzero) / Bscale;
# ifdef BYTE_SWAP
      SWAP_BYTE;
# endif
    }  
  }
  if (!strcmp (outtype, "int") && !strcmp (intype, "double")) {
    int NbytesIn = 8;
    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
      *(int *)Pout = (*(double *)Pin - Bzero) / Bscale;
# ifdef BYTE_SWAP
      SWAP_WORD;
# endif
    }
  }
  if (!strcmp (outtype, "int64_t") && !strcmp (intype, "double")) {
    int NbytesIn = 8;
    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
      *(int64_t *)Pout = (*(double *)Pin - Bzero) / Bscale;
# ifdef BYTE_SWAP
      SWAP_DBLE;
# endif
    }
  }
  if (!strcmp (outtype, "float") && !strcmp (intype, "double")) {
    int NbytesIn = 8;
    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
      *(float *)Pout = (*(double *)Pin - Bzero) / Bscale;
# ifdef BYTE_SWAP
      SWAP_WORD;
# endif
    }
  }
  if (!strcmp (outtype, "double") && !strcmp (intype, "double")) {
    int NbytesIn = 8;
    for (i = 0; i < Nval*Nrow; i++, Pin += NbytesIn, Pout += NbytesOut) {
      *(double *)Pout = (*(double *)Pin - Bzero) / Bscale;
# ifdef BYTE_SWAP
      SWAP_DBLE;
# endif
    }
  }

  /* check array space */
  if (Nx*Ny < Nx*(Nrow - 1) + Nstart + Nval*NbytesOut) {
    fprintf (stderr, "mismatch in array sizes\n");
    return (FALSE);
  }

  /* insert bytes from array into appropriate section of buffer */
  Pout = table[0].buffer + Nstart;
  Pin  = array;
  for (i = 0; i < Nrow; i++, Pout += Nx, Pin += Nval*NbytesOut) {
    memcpy (Pout, Pin, Nval*NbytesOut);
  }

  free (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_set_table_column (Header *header, FTable *table, char *label, void *data, off_t Nrow) {

  off_t Nx, Ny, nbytes;
  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);

  /* 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", OFF_T_FMT, 1,  &Nx);
  gfits_scan (header, "NAXIS2", OFF_T_FMT, 1,  &Ny);
  if (Ny == 0) { 
    Ny = Nrow;
    header[0].Naxis[1] = Ny;
    gfits_modify (header, "NAXIS2", OFF_T_FMT, 1,  Ny);

    nbytes = gfits_data_size (header);
    REALLOCATE (table[0].buffer, char, MAX (nbytes, 1));
    bzero (table[0].buffer, nbytes);
    table[0].datasize = nbytes;
  }
  if (Ny != Nrow) return (FALSE);

  /* 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;
  }

  /* print line with appropriate formatting */
  ALLOCATE (array, char, Nbytes*Nval*Nrow);
  ALLOCATE (line, char, Nval+1);
  Pin = data;
  Pout = array;
  if (!strcmp (type, "char")) {
    for (i = 0; i < Nval*Nrow; i++, Pin++, Pout++) {
      *Pout = *Pin;
    }
  }
  if (!strcmp (type, "int")) {
    for (i = 0; i < Nrow; i++, Pin+=4, Pout+=Nval) {
      snprintf (line, Nval + 1, cformat, *(int *)Pin);
      memcpy (Pout, line, Nval);
    }
  }
  if (!strcmp (type, "float")) {
    for (i = 0; i < Nrow; i++, Pin+=4, Pout+=Nval) {
      snprintf (line, Nval + 1, cformat, *(float *)Pin);
      memcpy (Pout, line, Nval);
    }
  }

  /* check array space */
  if (Nx*Ny < Nx*(Nrow - 1) + Nstart + Nval*Nbytes) {
    fprintf (stderr, "mismatch in array sizes\n");
    return (FALSE);
  }

  /* insert bytes from array into appropriate section of buffer */
  Pout = table[0].buffer + Nstart;
  Pin  = array;
  for (i = 0; i < Nrow; i++, Pout += Nx, Pin += Nval*Nbytes) {
    memcpy (Pout, Pin, Nval*Nbytes);
  }

  free (line);
  free (array);
  return (TRUE);
}

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

