# include "dvo.h"

AstromOffsetTable *AstromOffsetMapLoad (char *filename, int VERBOSE) {

  off_t Nmap;

  Header header;
  Matrix matrix;
  Header theader;
  FTable ftable;

  FILE *f = fopen (filename, "r");
  if (f == NULL) {
    if (VERBOSE) fprintf (stderr, "can't find Astrom Offset Map file %s\n", filename);
    return (NULL);
  }

  /* load in table data */
  ftable.header = &theader;
  if (!gfits_fread_header (f, &header)) {
    if (VERBOSE) fprintf (stderr, "can't read Astrom Offset Map header\n");
    fclose (f);
    return (NULL);
  }
  if (!gfits_fread_matrix (f, &matrix, &header)) {
    if (VERBOSE) fprintf (stderr, "can't read Astrom Offset Map matrix\n");
    gfits_free_header (&header);
    fclose (f);
    return (NULL);
  }

  // for now, we only have one flavor (6x6)

  if (!gfits_fread_ftable (f, &ftable, "ASTROM_OFFSET_MAP_DISK_6x6")) {
    if (VERBOSE) fprintf (stderr, "can't read Astrom Offset Map table\n");
    gfits_free_header (&header);
    gfits_free_matrix (&matrix);
    fclose (f);
    return (NULL);
  }
  AstromOffsetMap_Disk_6x6 *map_disk = gfits_table_get_AstromOffsetMap_Disk_6x6 (&ftable, &Nmap, NULL);
  if (!map_disk) {
    fprintf (stderr, "ERROR: failed to read Astrom Offset Map\n");
    exit (2);
  }

  // AstromOffsetMap_Disk_6x6 *map_disk is an external format stored as an array of maps.
  // Convert the disk array of maps to then internal format in a rich structure:
  AstromOffsetTable *table = AstromOffsetMapToTable (map_disk, Nmap);

  gfits_free_header (&theader);
  gfits_free_matrix (&matrix);
  gfits_free_table  (&ftable);

  return (table);
}

int AstromOffsetMapSave (AstromOffsetTable *table, char *filename) {

  off_t Nmap;

  Header header;
  Matrix matrix;
  Header theader;
  FTable ftable;

  FILE *f;

  /* make phu header (no matrix needed) */
  gfits_init_header (&header);
  header.extend = TRUE;
  gfits_create_header (&header);
  gfits_create_matrix (&header, &matrix);

  // AstromOffsetMap_Disk_6x6 *map_disk is an external format stored as an array of maps.
  // Convert the internal format in a rich structure into a disk array of maps.
  AstromOffsetMap_Disk_6x6 *map_disk = AstromOffsetTableToMap (table, &Nmap); 

  ftable.header = &theader;
  gfits_table_set_AstromOffsetMap_Disk_6x6 (&ftable, map_disk, Nmap);

  f = fopen (filename, "w");
  if (f == (FILE *) NULL) { 
    fprintf (stderr, "cannot open %s for output\n", filename);
    return (FALSE);
  }
  
  gfits_fwrite_header  (f, &header);
  gfits_fwrite_matrix  (f, &matrix);
  gfits_fwrite_Theader (f, &theader);
  gfits_fwrite_table  (f, &ftable);
  fclose (f);

  return (TRUE);
}

AstromOffsetTable *AstromOffsetMapToTable(AstromOffsetMap_Disk_6x6 *map_disk, off_t Nmap) {

  int i, j, k;

  AstromOffsetTable *table = NULL;
  ALLOCATE (table, AstromOffsetTable, 1);

  table->Nmap = Nmap;
  ALLOCATE (table->map, AstromOffsetMap *, Nmap);

  // find the max value of imageID
  int MaxTableID = 0;
  int MaxImageID = 0;
  for (i = 0; i < Nmap; i++) {
    MaxTableID = MAX(map_disk[i].tableID, MaxTableID);
    MaxImageID = MAX(map_disk[i].imageID, MaxImageID);
  }
  table->MaxTableID = MaxTableID;
  table->MaxImageID = MaxImageID;

  // generate the index and init values to -1
  ALLOCATE (table->imageIDtoTableSeq, int, MaxImageID + 1);
  for (i = 0; i <= MaxImageID; i++) {
    table->imageIDtoTableSeq[i] = -1;
  }

  // assign the ID values
  for (i = 0; i < Nmap; i++) {
    int ImageID = map_disk[i].imageID;
    myAssert (table->imageIDtoTableSeq[ImageID] == -1, "oops, duplicate image IDs");
    table->imageIDtoTableSeq[ImageID] = i;
  }

  // assign the map values (this allocates just the area needed for each image, not the
  // full 6x6, saving some memory while doing the analysis)
  for (i = 0; i < Nmap; i++) {
    ALLOCATE (table->map[i], AstromOffsetMap, 1);

    table->map[i][0].Nx 	  = map_disk[i].Nx;
    table->map[i][0].Ny 	  = map_disk[i].Ny;
    table->map[i][0].tableID = map_disk[i].tableID;
    table->map[i][0].imageID = map_disk[i].imageID;
    
    // STORE THESE VALUES?
    table->map[i][0].dY 	  = map_disk[i].dX;
    table->map[i][0].dX 	  = map_disk[i].dY;

    // since this was on disk, we obviously keep it
    table->map[i][0].keep 	  = TRUE;

    ALLOCATE (table->map[i][0].dXv, float *, map_disk[i].Nx);
    ALLOCATE (table->map[i][0].dYv, float *, map_disk[i].Nx);

    for (j = 0; j < map_disk[i].Nx; j++) {
      ALLOCATE (table->map[i][0].dXv[j], float, map_disk[i].Ny);
      ALLOCATE (table->map[i][0].dYv[j], float, map_disk[i].Ny);

      for (k = 0; k < map_disk[i].Ny; k++) {
	table->map[i][0].dXv[j][k] = map_disk[i].dXv[j][k];
	table->map[i][0].dYv[j][k] = map_disk[i].dYv[j][k];
      }
    }
  }
  return table;
}

AstromOffsetMap_Disk_6x6 *AstromOffsetTableToMap(AstromOffsetTable *table, off_t *Nmap) {

  int i, j, k;

  AstromOffsetMap_Disk_6x6 *map_disk = NULL;
  ALLOCATE (map_disk, AstromOffsetMap_Disk_6x6, table->Nmap);

  // assign the map values (this allocates just the area needed for each image, not the
  // full 6x6, saving some memory while doing the analysis)
  // some maps in the table should be skipped because their image no longer uses them
  int Ndisk = 0;
  for (i = 0; i < table->Nmap; i++) {
    if (!table->map[i][0].keep) continue;
    map_disk[Ndisk].Nx 	    = table->map[i][0].Nx;
    map_disk[Ndisk].Ny 	    = table->map[i][0].Ny;
    map_disk[Ndisk].tableID = table->map[i][0].tableID;
    map_disk[Ndisk].imageID = table->map[i][0].imageID;
    
    map_disk[Ndisk].dX 	    = table->map[i][0].dX;
    map_disk[Ndisk].dY 	    = table->map[i][0].dY;

    

    for (j = 0; j < map_disk[Ndisk].Nx; j++) {
      for (k = 0; k < map_disk[Ndisk].Ny; k++) {
	map_disk[Ndisk].dXv[j][k] = table->map[i][0].dXv[j][k];
	map_disk[Ndisk].dYv[j][k] = table->map[i][0].dYv[j][k];
      }
      for (k = map_disk[Ndisk].Ny; k < 6; k++) {
	map_disk[Ndisk].dXv[j][k] = 0.0;
	map_disk[Ndisk].dYv[j][k] = 0.0;
      }
    }
    for (j = map_disk[Ndisk].Nx; j < 6; j++) {
      for (k = 0; k < 6; k++) {
	map_disk[Ndisk].dXv[j][k] = 0.0;
	map_disk[Ndisk].dYv[j][k] = 0.0;
      }
    }
    Ndisk ++;
  }
  *Nmap = Ndisk;
  return map_disk;
}

