text:
#include "uavsar.h"
#include "asf_meta.h"
#include "asf_endian.h"
#include <ctype.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_multiroots.h>

#define SQR(x) (x*x)
#define FLOAT_COMPARE_TOLERANCE(a, b, t) (fabs (a - b) <= t ? 1: 0)
#define ASF_EXPORT_FLOAT_MICRON 0.000000001
#define FLOAT_EQUIVALENT(a, b) (FLOAT_COMPARE_TOLERANCE \
                                (a, b, ASF_EXPORT_FLOAT_MICRON))
#define MAX_LINE 512

char **get_uavsar_products(const char *data_type, char *type, int *num_product)
{
  char *rest, *token;
  int ii, product_count;
  char *tmp = (char *) MALLOC(sizeof(char)*60);
  strcpy(tmp, data_type);

  if (strcmp_case(type, "POLSAR") == 0) {
    product_count = 5;
    if (strcmp_case(tmp, "ALL") == 0)
      sprintf(tmp, "SLC,MLC,DAT,GRD,HGT");
  }
  else if (strcmp_case(type, "INSAR") == 0) {
    product_count = 9;
    if (strcmp_case(tmp, "ALL") == 0)
      sprintf(tmp, "AMP,INT,UNW,COR,AMP_GRD,INT_GRD,UNW_GRD,COR_GRD,HGT_GRD");
  }

  char **product = (char **) MALLOC(sizeof(char*) * product_count);
  for (ii = 0; ii < product_count; ii++) {
    product[ii] = (char *) MALLOC(sizeof(char)*10);
    strcpy(product[ii], "");
  }

  ii = 0;
  while ((token = strtok_r(tmp, ",", &rest))) {
    strcpy(product[ii], token);
    tmp = rest;
    ii++;
  }

  *num_product = ii;

  return product;
}

int parse_annotation_line(char *line, char *key, char *value)
{
  char *ks = line, *ke = line, *vs = line, *ve = line, *b = line;

  while(isspace(*b)) ks = ++b; // Move forward past the whitespace to the beginning of the key

  if(*b == ';' || *b == '\0') {
    // This line has nothing interesting, return nothing
    *key = *value = '\0';
    return 1;
  }

  while(*b != '(' && *b != '\0' && *b != '=') ++b; // Advance to a delimiter; either an equals sign, a '(' (to denote the beginning of the unit specification) or the end of the line

  if(*b == '\0') return 0; // Unexpected end-of-line

  ke = b - 1;
  while(isspace(*ke)) --ke; // Back up until we hit the end of the key

  while(*b != '\0' && *b != '=') b++; // Advance to the equals sign

  if(*b == '\0') return 0; // Unexpected end-of-line

  b++;
  while(isspace(*b) && *b != '\0') b++; // Move forward past the whitespace to the beginning of the value

  if(*b == '\0') return 0; // Unexpected end-of-line

  vs = b++;
  while(*b != ';' && *b != '\0') ++b; // Move forward to the end of the line or the beginning of a comment

  ve = b - 1;
  while(isspace(*ve)) --ve; // Back up until we hit the end of the key

  strncpy(key, ks, ke - ks + 1);
  key[ke-ks+1] = '\0';
  strncpy(value, vs, ve - vs + 1);
  value[ve-vs+1] = '\0';

  return 1;
}

int check_file(const char *path, char *line, char **fileName)
{
  char *p, *r;
  p = strchr(line, '=');
  char *file = (char *) MALLOC(sizeof(char)*255);
  strcpy(file, p+1);
  r = strchr(file, ';');
  if (r)
    r[0] = '\0';
  *fileName = MALLOC(sizeof(char)*(strlen(path) + strlen(file) + 32));
  strcpy(*fileName, path);
  char *trimmed_file = trim_spaces(file);
  strcat(*fileName, trimmed_file);
  FREE(file);
  FREE(trimmed_file);
  p = strchr(line, ';');
  long long size;
  sscanf(p+12, "%lld", &size);
  if (fileExists(*fileName) && fileSize(*fileName) == size)
    return TRUE;
  else
    return FALSE;
}

void
get_uavsar_file_names(const char *dataFile, uavsar_type_t type,
                      char ***pDataName, char ***pElement,
                      int **pDataType, int *nBands)
{
  int ii, slc = 0, mlc = 0, dat = 0, grd = 0, hgt = 0;
  int igram = 0, unw = 0, cor = 0, amp = 0;
  int igram_grd = 0, unw_grd = 0, cor_grd = 0, amp_grd = 0, hgt_grd = 0;

  char *file = (char *) MALLOC(sizeof(char) * 1024);
  char **dataName = (char **) MALLOC(6 * sizeof(char *));
  char **element = (char **) MALLOC(6 * sizeof(char *));
  int *dataType = (int *) MALLOC(6 * sizeof(int));
  for (ii = 0; ii < 6; ii++) {
    dataName[ii] = (char *) MALLOC(sizeof(char) * 512);
    element[ii] = (char *) MALLOC(sizeof(char) * 10);
  }
  *pDataName = dataName;
  *pElement = element;
  *pDataType = dataType;
  *nBands = 0;

  char *path = get_dirname(dataFile);

  char line[MAX_LINE], key[MAX_LINE], value[MAX_LINE];
  FILE *fp = FOPEN(dataFile, "r");
  while (fgets(line, MAX_LINE, fp)) {
    if(!parse_annotation_line(line, key, value)) {
      asfPrintWarning("Unable to parse line in annotation file: %s", line);
      continue;
    }
    if (!strcmp(key, ""))
      continue;
    if (type == POLSAR_SLC) {
      if (!strcmp(key, "slcHH")) {
        if (check_file(path, line, &file)) {
          strcpy(dataName[slc], file);
          strcpy(element[slc], "HH");
          dataType[slc] = 1;
          slc++;
        }
      }
      else if (!strcmp(key, "slcHV")) {
        if (check_file(path, line, &file)) {
          strcpy(dataName[slc], file);
          strcpy(element[slc], "HV");
          dataType[slc] = 1;
          slc++;
        }
      }
      else if (!strcmp(key, "slcVH")) {
        if (check_file(path, line, &file)) {
          strcpy(dataName[slc], file);
          strcpy(element[slc], "VH");
          dataType[slc] = 1;
          slc++;
        }
      }
      else if (!strcmp(key, "slcVV")) {
        if (check_file(path, line, &file)) {
          strcpy(dataName[slc], file);
          strcpy(element[slc], "VV");
          dataType[slc] = 1;
          slc++;
        }
      }
      *nBands = slc;
    }
    else if (type == POLSAR_MLC) {
      if (!strcmp(key, "mlcHHHH")) {
        if (check_file(path, line, &file)) {
          strcpy(dataName[mlc], file);
          strcpy(element[mlc], "C11");
          dataType[mlc] = 0;
          mlc++;
        }
      }
      else if (!strcmp(key, "mlcHVHV")) {
        if (check_file(path, line, &file)) {
          strcpy(dataName[mlc], file);
          strcpy(element[mlc], "C22");
          dataType[mlc] = 0;
          mlc++;
        }
      }
      else if (!strcmp(key, "mlcVVVV")) {
        if (check_file(path, line, &file)) {
          strcpy(dataName[mlc], file);
          strcpy(element[mlc], "C33");
          dataType[mlc] = 0;
          mlc++;
        }
      }
      else if (!strcmp(key, "mlcHHHV")) {
        if (check_file(path, line, &file)) {
          strcpy(dataName[mlc], file);
          strcpy(element[mlc], "C12");
          dataType[mlc] = 1;
          mlc++;
        }
      }
      else if (!strcmp(key, "mlcHHVV")) {
        if (check_file(path, line, &file)) {
          strcpy(dataName[mlc], file);
          strcpy(element[mlc], "C13");
          dataType[mlc] = 1;
          mlc++;
        }
      }
      else if (!strcmp(key, "mlcHVVV")) {
        if (check_file(path, line, &file)) {
          strcpy(dataName[mlc], file);
          strcpy(element[mlc], "C23");
          dataType[mlc] = 1;
          mlc++;
        }
      }
      *nBands = mlc;
    }
    else if (type == POLSAR_DAT) {
      if (!strcmp(key, "dat")) {
        if (check_file(path, line, &file)) {
          strcpy(dataName[dat], file);
          dat++;
        }
      }
      *nBands = dat;
    }
    else if (type == POLSAR_GRD) {
      if (!strcmp(key, "grdHHHH")) {
        if (check_file(path, line, &file)) {
          strcpy(dataName[grd], file);
          strcpy(element[grd], "C11");
          dataType[grd] = 0;
          grd++;
        }
      }
      else if (!strcmp(key, "grdHVHV")) {
        if (check_file(path, line, &file)) {
          strcpy(dataName[grd], file);
          strcpy(element[grd], "C22");
          dataType[grd] = 0;
          grd++;
        }
      }
      else if (!strcmp(key, "grdVVVV")) {
        if (check_file(path, line, &file)) {
          strcpy(dataName[grd], file);
          strcpy(element[grd], "C33");
          dataType[grd] = 0;
          grd++;
        }
      }
      else if (!strcmp(key, "grdHHHV")) {
        if (check_file(path, line, &file)) {
          strcpy(dataName[grd], file);
          strcpy(element[grd], "C12");
          dataType[grd] = 1;
          grd++;
        }
      }
      else if (!strcmp(key, "grdHHVV")) {
        if (check_file(path, line, &file)) {
          strcpy(dataName[grd], file);
          strcpy(element[grd], "C13");
          dataType[grd] = 1;
          grd++;
        }
      }
      else if (!strcmp(key, "grdHVVV")) {
        if (check_file(path, line, &file)) {
          strcpy(dataName[grd], file);
          strcpy(element[grd], "C23");
          dataType[grd] = 1;
          grd++;
        }
      }
      *nBands = grd;
    }
    else if (type == POLSAR_HGT) {
      if (!strcmp(key, "hgt")) {
        if (check_file(path, line, &file)) {
          strcpy(dataName[hgt], file);
          strcpy(element[hgt], "HH");
          dataType[hgt] = 0;
          hgt++;
        }
      }
      *nBands = hgt;
    }
    else if (type == INSAR_AMP) {
      if (!strcmp(key, "Slant Range Amplitude of Pass 1")) {
        if (check_file(path, line, &file)) {
          strcpy(dataName[amp], file);
          strcpy(element[amp], "HH");
          dataType[amp] = 0;
          amp++;
        }
      }
      else if (!strcmp(key, "Slant Range Amplitude of Pass 2")) {
        if (check_file(path, line, &file)) {
          strcpy(dataName[amp], file);
          strcpy(element[amp], "HH");
          dataType[amp] = 0;
          amp++;
        }
      }
      *nBands = amp;
    }
    else if (type == INSAR_AMP_GRD) {
      if (!strcmp(key, "Ground Range Amplitude of Pass 1")) {
        if (check_file(path, line, &file)) {
          strcpy(dataName[amp_grd], file);
          strcpy(element[amp_grd], "HH");
          dataType[amp_grd] = 0;
          amp_grd++;
        }
      }
      else if (!strcmp(key, "Ground Range Amplitude of Pass 2")) {
        if (check_file(path, line, &file)) {
          strcpy(dataName[amp_grd], file);
          strcpy(element[amp_grd], "HH");
          dataType[amp_grd] = 0;
          amp_grd++;
        }
      }
      *nBands = amp_grd;
    }
    else if (type == INSAR_INT) {
      if (!strcmp(key, "Slant Range Interferogram")) {
        if (check_file(path, line, &file)) {
          strcpy(dataName[igram], file);
          strcpy(element[igram], "HH");
          dataType[igram] = 1;
          igram++;
        }
      }
      *nBands = igram;
    }
    else if (type == INSAR_INT_GRD) {
      if (!strcmp(key, "Ground Range Interferogram")) {
        if (check_file(path, line, &file)) {
          strcpy(dataName[igram_grd], file);
          strcpy(element[igram_grd], "HH");
          dataType[igram_grd] = 1;
          igram_grd++;
        }
      }
      *nBands = igram_grd;
    }
    else if (type == INSAR_UNW) {
      if (!strcmp(key, "Slant Range Unwrapped Phase")) {
        if (check_file(path, line, &file)) {
          strcpy(dataName[unw], file);
          strcpy(element[unw], "HH");
          dataType[unw] = 0;
          unw++;
        }
      }
      *nBands = unw;
    }
    else if (type == INSAR_UNW_GRD) {
      if (!strcmp(key, "Ground Range Unwrapped Phase")) {
        if (check_file(path, line, &file)) {
          strcpy(dataName[unw_grd], file);
          strcpy(element[unw_grd], "HH");
          dataType[unw_grd] = 0;
          unw_grd++;
        }
      }
      *nBands = unw_grd;
    }
    else if (type == INSAR_COR) {
      if (!strcmp(key, "Slant Range Correlation")) {
        if (check_file(path, line, &file)) {
          strcpy(dataName[cor], file);
          strcpy(element[cor], "HH");
          dataType[cor] = 0;
          cor++;
        }
      }
      *nBands = cor;
    }
    else if (type == INSAR_COR_GRD) {
      if (!strcmp(key, "Ground Range Correlation")) {
        if (check_file(path, line, &file)) {
          strcpy(dataName[cor_grd], file);
          strcpy(element[cor_grd], "HH");
          dataType[cor_grd] = 0;
          cor_grd++;
        }
      }
      *nBands = cor_grd;
    }
    else if (type == INSAR_HGT_GRD) {
      if (!strcmp(key, "DEM Used in Ground Projection")) {
        if (check_file(path, line, &file)) {
          strcpy(dataName[hgt_grd], file);
          strcpy(element[hgt_grd], "HH");
          dataType[hgt_grd] = 0;
          hgt_grd++;
        }
      }
      *nBands = hgt_grd;
    }
  }

  FCLOSE(fp);
  FREE(file);
}

uavsar_polsar *
read_uavsar_polsar_params(const char *dataFile, uavsar_type_t type)
{
  uavsar_polsar *params = (uavsar_polsar *) MALLOC(sizeof(uavsar_polsar));

  // Determine ID
  char *dirName = (char *) MALLOC(sizeof(char) * 1024);
  char *fileName = (char *) MALLOC(sizeof(char) * 1024);
  split_dir_and_file(dataFile, dirName, fileName);
  sprintf(params->id, "%s", stripExt(fileName));

  // Read annotation file
  char line[MAX_LINE], key[MAX_LINE], value[MAX_LINE];
  FILE *fp = FOPEN(dataFile, "r");
  while (fgets(line, MAX_LINE, fp)) {
    if(!parse_annotation_line(line, key, value)) {
      asfPrintWarning("Unable to parse line in annotation file: %s", line);
      continue;
    }
    if (!strcmp(key, ""))
      continue;
    if (!strcmp(key, "Site Description"))
      strcpy(params->site, value);
    if (!strcmp(key, "Acquisition Mode"))
      strcpy(params->acquisition_mode, value);
    if (type == POLSAR_SLC) {
      params->type = POLSAR_SLC;
      if (!strcmp(key, "slc_mag.set_rows"))
        params->row_count = atoi(value);
      else if (!strcmp(key, "slc_mag.set_cols"))
        params->column_count = atoi(value);
      else if (!strcmp(key, "slc_mag.set_proj"))
        strcpy(params->projection, value);
      else if (!strcmp(key, "slc_mag.row_addr"))
        params->along_track_offset = atof(value);
      else if (!strcmp(key, "slc_mag.col_addr"))
        params->cross_track_offset = atof(value);
      else if (!strcmp(key, "slc_mag.row_mult"))
        params->azimuth_pixel_spacing = atof(value);
      else if (!strcmp(key, "slc_mag.col_mult"))
        params->range_pixel_spacing = atof(value);
      else if (!strcmp(key, "slc_mag.val_size"))
        params->bytes_per_pixel = atoi(value);
      else if (!strcmp(key, "slc_mag.val_frmt"))
        strcpy(params->value_format, value);
      params->range_look_count = 1;
      params->azimuth_look_count = 1;
      if (!strcmp(key, "SLC Data Units"))
        strcpy(params->data_units, value);
    }
    else if (type == POLSAR_MLC) {
      params->type = POLSAR_MLC;
      if (!strcmp(key, "mlc_pwr.set_rows"))
        params->row_count = atoi(value);
      else if (!strcmp(key, "mlc_pwr.set_cols"))
        params->column_count = atoi(value);
      else if (!strcmp(key, "mlc_pwr.set_proj"))
        strcpy(params->projection, value);
      else if (!strcmp(key, "mlc_pwr.row_addr"))
        params->along_track_offset = atof(value);
      else if (!strcmp(key, "mlc_pwr.col_addr"))
        params->cross_track_offset = atof(value);
      else if (!strcmp(key, "mlc_pwr.row_mult"))
        params->azimuth_pixel_spacing = atof(value);
      else if (!strcmp(key, "mlc_pwr.col_mult"))
        params->range_pixel_spacing = atof(value);
      else if (!strcmp(key, "mlc_pwr.val_size"))
        params->bytes_per_pixel = atoi(value);
      else if (!strcmp(key, "mlc_pwr.val_frmt"))
        strcpy(params->value_format, value);
      else if (!strcmp(key, "Number of Range Looks in MLC"))
        params->range_look_count = atoi(value);
      else if (!strcmp(key, "Number of Azimuth Looks in MLC"))
        params->azimuth_look_count = atoi(value);
      else if (!strcmp(key, "MLC Data Units"))
        strcpy(params->data_units, value);
      params->range_look_count = 1;
      params->azimuth_look_count = 1;
    }
    else if (type == POLSAR_DAT) {
      params->type = POLSAR_DAT;
      if (!strcmp(key, "dat.set_rows"))
        params->row_count = atoi(value);
      else if (!strcmp(key, "dat.set_cols"))
        params->column_count = atoi(value);
      else if (!strcmp(key, "dat.set_proj"))
        strcpy(params->projection, value);
      else if (!strcmp(key, "dat.row_addr"))
        params->along_track_offset = atof(value);
      else if (!strcmp(key, "dat.col_addr"))
        params->cross_track_offset = atof(value);
      else if (!strcmp(key, "dat.row_mult"))
        params->azimuth_pixel_spacing = atof(value);
      else if (!strcmp(key, "dat.col_mult"))
        params->range_pixel_spacing = atof(value);
      else if (!strcmp(key, "dat.val_size"))
        params->bytes_per_pixel = atoi(value);
      else if (!strcmp(key, "dat.val_frmt"))
        strcpy(params->value_format, value);
      params->range_look_count = 1;
      params->azimuth_look_count = 1;
      strcpy(params->data_units, MAGIC_UNSET_STRING);
    }
    else if (type == POLSAR_GRD) {
      params->type = POLSAR_GRD;
      if (!strcmp(key, "grd_pwr.set_rows"))
        params->row_count = atoi(value);
      else if (!strcmp(key, "grd_pwr.set_cols"))
        params->column_count = atoi(value);
      else if (!strcmp(key, "grd_pwr.set_proj"))
        strcpy(params->projection, value);
      else if (!strcmp(key, "grd_pwr.row_addr"))
        params->along_track_offset = atof(value);
      else if (!strcmp(key, "grd_pwr.col_addr"))
        params->cross_track_offset = atof(value);
      else if (!strcmp(key, "grd_pwr.row_mult"))
        params->azimuth_pixel_spacing = atof(value);
      else if (!strcmp(key, "grd_pwr.col_mult"))
        params->range_pixel_spacing = atof(value);
      else if (!strcmp(key, "grd_pwr.val_size"))
        params->bytes_per_pixel = atoi(value);
      else if (!strcmp(key, "grd_pwr.val_frmt"))
        strcpy(params->value_format, value);
      else if (!strcmp(key, "Number of Range Looks in MLC"))
        params->range_look_count = atoi(value);
      else if (!strcmp(key, "Number of Azimuth Looks in MLC"))
        params->azimuth_look_count = atoi(value);
      else if (!strcmp(key, "GRD Data Units"))
        strcpy(params->data_units, value);
    }
    else if (type == POLSAR_HGT) {
      params->type = POLSAR_HGT;
      if (!strcmp(key, "hgt.set_rows"))
        params->row_count = atoi(value);
      else if (!strcmp(key, "hgt.set_cols"))
        params->column_count = atoi(value);
      else if (!strcmp(key, "hgt.set_proj"))
        strcpy(params->projection, value);
      else if (!strcmp(key, "hgt.row_addr"))
        params->along_track_offset = atof(value);
      else if (!strcmp(key, "hgt.col_addr"))
        params->cross_track_offset = atof(value);
      else if (!strcmp(key, "hgt.row_mult"))
        params->azimuth_pixel_spacing = atof(value);
      else if (!strcmp(key, "hgt.col_mult"))
        params->range_pixel_spacing = atof(value);
      else if (!strcmp(key, "hgt.val_size"))
        params->bytes_per_pixel = atoi(value);
      else if (!strcmp(key, "hgt.val_frmt"))
        strcpy(params->value_format, value);
      params->range_look_count = 1;
      params->azimuth_look_count = 1;
      if (!strcmp(key, "HGT Data Units"))
        strcpy(params->data_units, value);
    }
    if (!strcmp(key, "set_hddr"))
      params->header_bytes = atoi(value);
    else if (!strcmp(key, "set_tail"))
      params->tail_bytes = atoi(value);
    else if (!strcmp(key, "set_plat"))
      params->lat_peg_point = atof(value);
    else if (!strcmp(key, "set_plon"))
      params->lon_peg_point = atof(value);
    else if (!strcmp(key, "set_phdg"))
      params->head_peg_point = atof(value);
    else if (!strcmp(key, "val_endi"))
      strcpy(params->endianess, value);
    else if (!strcmp(key, "val_mult"))
      params->data_scale = atof(value);
    else if (!strcmp(key, "val_addr"))
      params->data_shift = atof(value);
    else if (!strcmp(key, "val_minv"))
      params->min_value = atof(value);
    else if (!strcmp(key, "val_maxv"))
      params->max_value = atof(value);
    else if (!strcmp(key, "Center Wavelength"))
      params->wavelength = atof(value);
    else if (!strcmp(key, "Ellipsoid Semi-major Axis"))
      params->semi_major = atof(value);
    else if (!strcmp(key, "Ellipsoid Eccentricity Squared"))
      params->eccentricity = atof(value);
    else if (!strcmp(key, "Look Direction"))
      strcpy(params->look_direction, value);
    else if (!strcmp(key, "Range Spacing per Bin"))
      params->range_spacing = atof(value);
    else if (!strcmp(key, "Azimuth Spacing"))
      params->azimuth_spacing = atof(value);
    else if (!strcmp(key, "Image Starting Range"))
      params->slant_range_first_pixel = atof(value);
    else if (!strcmp(key, "Global Average Yaw"))
      params->yaw = atof(value);
    else if (!strcmp(key, "Global Average Pitch"))
      params->pitch = atof(value);
    else if (!strcmp(key, "Global Average Roll"))
      params->roll = atof(value);
    else if (!strcmp(key, "Global Average Altitude"))
      params->altitude = atof(value);
    else if (!strcmp(key, "Average GPS Altitude"))
      params->altitude = atof(value);
    else if (!strcmp(key, "Global Average Terrain Height"))
      params->terrain_height = atof(value);
    else if (!strcmp(key, "Global Average Squint Angle"))
      params->squint_angle = atof(value);
    else if (!strcmp(key, "Pulse Length"))
      params->pulse_length = atof(value);
    else if (!strcmp(key, "Steering Angle"))
      params->steering_angle = atof(value);
    else if (!strcmp(key, "Bandwidth"))
      params->bandwidth = atof(value);
    else if (!strcmp(key, "Approximate Upper Left Latitude"))
      params->lat_upper_left = atof(value);
    else if (!strcmp(key, "Approximate Upper Left Longitude"))
      params->lon_upper_left = atof(value);
    else if (!strcmp(key, "Approximate Upper Right Latitude"))
      params->lat_upper_right = atof(value);
    else if (!strcmp(key, "Approximate Upper Right Longitude"))
      params->lon_upper_right = atof(value);
    else if (!strcmp(key, "Approximate Lower Left Latitude"))
      params->lat_lower_left = atof(value);
    else if (!strcmp(key, "Approximate Lower Left Longitude"))
      params->lon_lower_left = atof(value);
    else if (!strcmp(key, "Approximate Lower Right Latitude"))
      params->lat_lower_right = atof(value);
    else if (!strcmp(key, "Approximate Lower Right Longitude"))
      params->lon_lower_right = atof(value);
    else if (!strcmp(key, "Date of Acquisition"))
      strcpy(params->acquisition_date, value);
    else if (!strcmp(key, "Processor Version Number"))
      strcpy(params->processor, value);
  }
  FCLOSE(fp);

  return params;
}

uavsar_insar *
read_uavsar_insar_params(const char *dataFile, uavsar_type_t type)
{
  uavsar_insar *params = (uavsar_insar *) MALLOC(sizeof(uavsar_insar));
  char time1[50], time2[50];

  // Determine ID
  char *dirName = (char *) MALLOC(sizeof(char) * 1024);
  char *fileName = (char *) MALLOC(sizeof(char) * 1024);
  split_dir_and_file(dataFile, dirName, fileName);
  sprintf(params->id, "%s", stripExt(fileName));

  // Read annotation file
  char line[MAX_LINE], key[MAX_LINE], value[MAX_LINE];
  FILE *fp = FOPEN(dataFile, "r");
  while (fgets(line, MAX_LINE, fp)) {
    if(!parse_annotation_line(line, key, value)) {
      asfPrintWarning("Unable to parse line in annotation file: %s", line);
      continue;
    }
    if(!strcmp(key, ""))
      continue;
    if (!strcmp(key, "Site Description"))
      strcpy(params->site, value);
    if (!strcmp(key, "Processing Mode"))
      strcpy(params->processing_mode, value);
    if (!strcmp(key, "Polarization"))
      strcpy(params->polarization, value);
    if (!strcmp(key, "Number of Looks in Range"))
      params->range_look_count = atoi(value);
    if (!strcmp(key, "Number of Looks in Azimuth"))
      params->azimuth_look_count = atoi(value);
    if (type == INSAR_INT) {
      params->type = INSAR_INT;
      if (!strcmp(key, "Interferogram Bytes Per Pixel"))
        params->bytes_per_pixel = atoi(value);
      else if (!strcmp(key, "Interferogram Pixel Format"))
        strcpy(params->value_format, value);
      else if (!strcmp(key, "Interferogram Units"))
        strcpy(params->data_units, value);
    }
    else if (type == INSAR_UNW) {
      params->type = INSAR_UNW;
      if (!strcmp(key, "Unwrapped Phase Bytes Per Pixel"))
        params->bytes_per_pixel = atoi(value);
      else if (!strcmp(key, "Unwrapped Phase Pixel Format"))
        strcpy(params->value_format, value);
      else if (!strcmp(key, "Unwrapped Phase Units"))
        strcpy(params->data_units, value);
    }
    else if (type == INSAR_COR) {
      params->type = INSAR_COR;
      if (!strcmp(key, "Correlation Bytes Per Pixel"))
        params->bytes_per_pixel = atoi(value);
      else if (!strcmp(key, "Correlation Pixel Format"))
        strcpy(params->value_format, value);
      else if (!strcmp(key, "Correlation Units"))
        strcpy(params->data_units, value);
    }
    else if (type == INSAR_AMP) {
      params->type = INSAR_AMP;
      if (!strcmp(key, "Amplitude Bytes Per Pixel"))
        params->bytes_per_pixel = atoi(value);
      else if (!strcmp(key, "Amplitude Pixel Format"))
        strcpy(params->value_format, value);
      else if (!strcmp(key, "Amplitude Units"))
        strcpy(params->data_units, value);
    }
    if (type >= INSAR_AMP && type <= INSAR_COR) {
      if (!strcmp(key, "Slant Range Data Azimuth Lines"))
        params->row_count = atoi(value);
      else if (!strcmp(key, "Slant Range Data Range Samples"))
        params->column_count = atoi(value);
      else if (!strcmp(key, "slt.set_proj"))
        strcpy(params->projection, value);
      else if (!strcmp(key, "slt.row_addr"))
        params->along_track_offset = atof(value);
      else if (!strcmp(key, "slt.col_addr"))
        params->cross_track_offset = atof(value);
      else if (!strcmp(key, "Slant Range Data at Near Range")) {
        params->slant_range_first_pixel = atof(value);
        params->slant_range_first_pixel /= 1000.0;
      }
      else if (!strcmp(key, "Slant Range Data Azimuth Spacing"))
        params->azimuth_pixel_spacing = atof(value);
      else if (!strcmp(key, "Slant Range Data Range Spacing"))
        params->range_pixel_spacing = atof(value);
    }
    else if (type >= INSAR_AMP_GRD && type <= INSAR_HGT_GRD) {
      if (!strcmp(key, "Ground Range Data Latitude Lines"))
        params->row_count = atoi(value);
      else if (!strcmp(key, "Ground Range Data Longitude Samples"))
        params->column_count = atoi(value);
      else if (!strcmp(key, "grd.set_proj"))
        strcpy(params->projection, value);
      else if (!strcmp(key, "grd.row_addr"))
        params->along_track_offset = atof(value);
      else if (!strcmp(key, "grd.col_addr"))
        params->cross_track_offset = atof(value);
      else if (!strcmp(key, "Ground Range Data Latitude Spacing"))
        params->azimuth_pixel_spacing = atof(value);
      else if (!strcmp(key, "Ground Range Data Longitude Spacing"))
        params->range_pixel_spacing = atof(value);
    }
    if (type == INSAR_INT_GRD) {
      params->type = INSAR_INT_GRD;
      if (!strcmp(key, "Interferogram Bytes Per Pixel"))
        params->bytes_per_pixel = atoi(value);
      else if (!strcmp(key, "Interferogram Pixel Format"))
        strcpy(params->value_format, value);
      else if (!strcmp(key, "Interferogram Units"))
        strcpy(params->data_units, value);
    }
    else if (type == INSAR_UNW_GRD) {
      params->type = INSAR_UNW_GRD;
      if (!strcmp(key, "Unwrapped Phase Bytes Per Pixel"))
        params->bytes_per_pixel = atoi(value);
      else if (!strcmp(key, "Unwrapped Phase Pixel Format"))
        strcpy(params->value_format, value);
      else if (!strcmp(key, "Unwrapped Phase Units"))
        strcpy(params->data_units, value);
    }
    else if (type == INSAR_COR_GRD) {
      params->type = INSAR_COR_GRD;
      if (!strcmp(key, "Correlation Bytes Per Pixel"))
        params->bytes_per_pixel = atoi(value);
      else if (!strcmp(key, "Correlation Pixel Format"))
        strcpy(params->value_format, value);
      else if (!strcmp(key, "Correlation Units"))
        strcpy(params->data_units, value);
    }
    else if (type == INSAR_AMP_GRD) {
      params->type = INSAR_AMP_GRD;
      if (!strcmp(key, "Amplitude Bytes Per Pixel"))
        params->bytes_per_pixel = atoi(value);
      else if (!strcmp(key, "Amplitude Pixel Format"))
        strcpy(params->value_format, value);
      else if (!strcmp(key, "Amplitude Units"))
        strcpy(params->data_units, value);
    }
    else if (type == INSAR_HGT_GRD) {
      params->type = INSAR_HGT_GRD;
      if (!strcmp(key, "DEM Bytes Per Pixel"))
        params->bytes_per_pixel = atoi(value);
      else if (!strcmp(key, "DEM Pixel Format"))
        strcpy(params->value_format, value);
      else if (!strcmp(key, "DEM Units"))
        strcpy(params->data_units, value);
    }
    if (!strcmp(key, "set_hddr"))
      params->header_bytes = atoi(value);
    else if (!strcmp(key, "set_tail"))
      params->tail_bytes = atoi(value);
    else if (!strcmp(key, "set_plat"))
      params->lat_peg_point = atof(value);
    else if (!strcmp(key, "set_plon"))
      params->lon_peg_point = atof(value);
    else if (!strcmp(key, "set_phdg"))
      params->head_peg_point = atof(value);
    else if (!strcmp(key, "val_endi"))
      strcpy(params->endianess, value);
    else if (!strcmp(key, "val_mult"))
      params->data_scale = atof(value);
    else if (!strcmp(key, "val_addr"))
      params->data_shift = atof(value);
    else if (!strcmp(key, "val_minv"))
      params->min_value = atof(value);
    else if (!strcmp(key, "val_maxv"))
      params->max_value = atof(value);
    else if (!strcmp(key, "Center Wavelength"))
      params->wavelength = atof(value);
    else if (!strcmp(key, "Ellipsoid Semi-major Axis"))
      params->semi_major = atof(value);
    else if (!strcmp(key, "Ellipsoid Eccentricity Squared"))
      params->eccentricity = atof(value);
    else if (!strcmp(key, "Look Direction"))
      strcpy(params->look_direction, value);
    else if (!strcmp(key, "Range Spacing per Bin"))
      params->range_spacing = atof(value);
    else if (!strcmp(key, "Azimuth Spacing"))
      params->azimuth_spacing = atof(value);
    else if (!strcmp(key, "Image Starting Range"))
      params->slant_range_first_pixel = atof(value);
    else if (!strcmp(key, "Global Average Yaw"))
      params->yaw = atof(value);
    else if (!strcmp(key, "Global Average Pitch"))
      params->pitch = atof(value);
    else if (!strcmp(key, "Global Average Roll"))
      params->roll = atof(value);
    else if (!strcmp(key, "Global Average Altitude"))
      params->altitude = atof(value);
    else if (!strcmp(key, "Average GPS Altitude"))
      params->altitude = atof(value);
    else if (!strcmp(key, "Global Average Terrain Height"))
      params->terrain_height = atof(value);
    else if (!strcmp(key, "Global Average Squint Angle"))
      params->squint_angle = atof(value);
    else if (!strcmp(key, "Pulse Length"))
      params->pulse_length = atof(value);
    else if (!strcmp(key, "Steering Angle"))
      params->steering_angle = atof(value);
    else if (!strcmp(key, "Bandwidth"))
      params->bandwidth = atof(value);
    else if (!strcmp(key, "Approximate Upper Left Latitude")) {
      if (strcmp_case(value, "N/A") == 0)
        params->lat_upper_left = MAGIC_UNSET_DOUBLE;
      else
        params->lat_upper_left = atof(value);
    }
    else if (!strcmp(key, "Approximate Upper Left Longitude")) {
      if (strcmp_case(value, "N/A") == 0)
        params->lon_upper_left = MAGIC_UNSET_DOUBLE;
      else
        params->lon_upper_left = atof(value);
    }
    else if (!strcmp(key, "Approximate Upper Right Latitude")) {
      if (strcmp_case(value, "N/A") == 0)
        params->lat_upper_right = MAGIC_UNSET_DOUBLE;
      else
        params->lat_upper_right = atof(value);
    }
    else if (!strcmp(key, "Approximate Upper Right Longitude")) {
      if (strcmp_case(value, "N/A") == 0)
        params->lon_upper_right = MAGIC_UNSET_DOUBLE;
      else
        params->lon_upper_right = atof(value);
    }
    else if (!strcmp(key, "Approximate Lower Left Latitude")) {
      if (strcmp_case(value, "N/A") == 0)
        params->lat_lower_left = MAGIC_UNSET_DOUBLE;
      else
        params->lat_lower_left = atof(value);
    }
    else if (!strcmp(key, "Approximate Lower Left Longitude")) {
      if (strcmp_case(value, "N/A") == 0)
        params->lon_lower_left = MAGIC_UNSET_DOUBLE;
      else
        params->lon_lower_left = atof(value);
    }
    else if (!strcmp(key, "Approximate Lower Right Latitude")) {
      if (strcmp_case(value, "N/A") == 0)
        params->lat_lower_right = MAGIC_UNSET_DOUBLE;
      else
        params->lat_lower_right = atof(value);
    }
    else if (!strcmp(key, "Approximate Lower Right Longitude")) {
      if (strcmp_case(value, "N/A") == 0)
        params->lon_lower_right = MAGIC_UNSET_DOUBLE;
      else
        params->lon_lower_right = atof(value);
    }
    else if (!strcmp(key, "Time of Acquisition for Pass 1"))
      strcpy(time1, value);
    else if (!strcmp(key, "Time of Acquisition for Pass 2"))
      strcpy(time2, value);
    else if (!strcmp(key, "Processor Version Number"))
      strcpy(params->processor, value);
  }
  FCLOSE(fp);

  sprintf(params->acquisition_date, "%s, %s", time1, time2);

  return params;
}

int sign(char byteBuf)
{
  if (byteBuf < 0)
    return -1;
  else
    return 1;
}

char *check_data_type(const char *inFileName)
{
  char line[MAX_LINE], key[MAX_LINE], value[MAX_LINE];
  char *type = (char *) MALLOC(sizeof(char)*25);
  FILE *fp = FOPEN(inFileName, "r");
  while (fgets(line, MAX_LINE, fp)) {
    if(!parse_annotation_line(line, key, value)) {
      asfPrintWarning("Unable to parse line in annotation file: %s", line);
      continue;
    }
    if (!strcmp(key, ""))
      continue;
    if (!strcmp(key, "Acquisition Mode"))
      strcpy(type, value);
    else if (!strcmp(key, "Processing Mode"))
      strcpy(type, value);
    else if (strcmp_case(type, "RPI") == 0)
      sprintf(type, "InSAR");
  }
  FCLOSE(fp);

  return type;
}

void import_uavsar(const char *inFileName, int line, int sample, int width,
		   int height, radiometry_t radiometry,
		   const char *data_type, const char *outBaseName) {

  // UAVSAR comes in two flavors: InSAR and PolSAR
  // Things look similar to AirSAR data, just organized a little different.
  // There does not seem to be a consistent identifier in the annotation file,
  // that would allow us to easily identify the data set as UAVSAR. No mention
  // of UAVSAR whatsoever.

  // The data can come in a large variety of flavors (single look versus multi-
  // look, derived magnitude and phase, etc.). I assume we can take anything or
  // nothing from this menu. The different files have different dimensions,
  // so we will need to generate several output images to accommodate that.

  // PolSAR data
  // slc - Single look complex slant range data
  // mlc - Multilook cross product slant range data
  // dat - Compressed Stokes matrix of multilooked data
  // grd - Ground range projected (equi-rectangular) and multilooked data
  // hgt - Digital elevation model projected in projection
  // slc_mag and slc_phase - 8 bytes per pixel derived from slc
  // mlc_mag and mlc_phase - 8 bytes per pixel derived from mlc

  // InSAR data
  // int - Slant range interferogram
  // unw - Slant range unwrapped phase
  // cor - Slant range correlation
  // amp - Slant range amplitudes
  // int_grd - Ground range interferogram
  // unw_grd - Ground range unwrapped phase
  // cor_grd - Ground range correlation
  // amp_grd - Ground range amplitudes
  // hgt_grd - Digital elevation model in ground projection

  FILE *fpIn, *fpOut;
  int ii, kk, ll, nn, pp, nBands, ns, *dataType, product_count;
  int multi = FALSE;
  float *floatAmp, *floatPhase, *floatAmpBuf, *amp, re, im;
  float *floatComplexReal, *floatComplexImag;
  float *floatComplexBuf;
  char **dataName, **element, **product, tmp[50];
  char *type;
  char *outName = (char *) MALLOC(sizeof(char)*(strlen(outBaseName)+15));
  uavsar_polsar *polsar_params;
  uavsar_insar *insar_params;
  meta_parameters *metaIn, *metaOut;

  type = check_data_type(inFileName);
  asfPrintStatus("   Data type: %s\n", type);
  product = get_uavsar_products(data_type, type, &product_count);
  if (product_count > 1)
    multi = TRUE;

  for (pp = 0; pp < product_count; pp++) {

    // InSAR data
    // Ground range interferogram
    if (strcmp_case(type, "InSAR") == 0 &&
	strcmp_case(product[pp], "INT_GRD") == 0) {
      get_uavsar_file_names(inFileName, INSAR_INT_GRD, &dataName, &element,
			    &dataType, &nBands);
      if (!nBands) {
	asfPrintWarning("Ground range interferogram does not exist. "
			"Will skip the ingest.\n");
	continue;
      }
      insar_params =
	read_uavsar_insar_params(inFileName, INSAR_INT_GRD);
      metaIn = uavsar_insar2meta(insar_params);
      metaOut = uavsar_insar2meta(insar_params);
      ns = metaIn->general->sample_count;
      nn = 0;
      floatAmp = (float *) CALLOC(ns, sizeof(float));
      floatPhase = (float *) CALLOC(ns, sizeof(float));
      floatComplexBuf = (float *) CALLOC(2*ns, sizeof(float));
      outName = (char *) MALLOC(sizeof(char)*(strlen(outBaseName)+15));
      metaOut->general->band_count = 2;
      if (multi)
	outName = appendToBasename(outBaseName, "_int_grd.img");
      else
	outName = appendExt(outBaseName, ".img");
      asfPrintStatus("\nGround range interferogram:\n");
      fpOut = FOPEN(outName, "wb");
      char *filename = get_filename(dataName[nn]);
      asfPrintStatus("Ingesting %s ...\n", filename);
      FREE(filename);
      fpIn = FOPEN(dataName[nn], "rb");
      sprintf(metaOut->general->bands, "INTERFEROGRAM_AMP,INTERFEROGRAM_PHASE");
      for (ii=0; ii<metaIn->general->line_count; ii++) {
	metaIn->general->sample_count = 2*ns;
	get_float_line(fpIn, metaIn, ii, floatComplexBuf);
	for (kk=0; kk<ns; kk++) {
	  re = floatComplexBuf[kk*2];
	  im = floatComplexBuf[kk*2+1];
	  ieee_big32(re);
	  ieee_big32(im);
	  floatAmp[kk] = hypot(re,im);
	  floatPhase[kk] = atan2_check(im,re);
	}
	put_band_float_line(fpOut, metaOut, 0, ii, floatAmp);
	put_band_float_line(fpOut, metaOut, 1, ii, floatPhase);
	asfLineMeter(ii, metaIn->general->line_count);
      }
      FCLOSE(fpIn);
      FCLOSE(fpOut);
      meta_write(metaOut, outName);
      FREE(floatAmp);
      FREE(floatPhase);
      FREE(floatComplexBuf);
      FREE(outName);
      meta_free(metaIn);
      meta_free(metaOut);
      FREE(insar_params);
    }

    // Ground range unwrapped phase
    if (strcmp_case(type, "InSAR") == 0 &&
	strcmp_case(product[pp], "UNW_GRD") == 0) {
      get_uavsar_file_names(inFileName, INSAR_UNW_GRD, &dataName, &element,
			    &dataType, &nBands);
      if (!nBands) {
	asfPrintWarning("Ground range unwrapped phase does not exist. "
			"Will skip the ingest.\n");
	continue;
      }
      insar_params =
	read_uavsar_insar_params(inFileName, INSAR_UNW_GRD);
      metaIn = uavsar_insar2meta(insar_params);
      metaOut = uavsar_insar2meta(insar_params);
      ns = metaOut->general->sample_count;
      nn = 0;
      floatAmpBuf = (float *) MALLOC(sizeof(float)*ns);
      outName = (char *) MALLOC(sizeof(char)*(strlen(outBaseName)+15));
      if (multi)
	outName = appendToBasename(outBaseName, "_unw_grd.img");
      else
	outName = appendExt(outBaseName, ".img");
      asfPrintStatus("\nGround range unwrapped phase:\n");
      char *filename = get_filename(dataName[nn]);
      asfPrintStatus("Ingesting %s ...\n", filename);
      FREE(filename);
      fpIn = FOPEN(dataName[nn], "rb");
      fpOut = FOPEN(outName, "wb");
      strcpy(metaOut->general->bands, "UNWRAPPED_PHASE");
      for (ii=0; ii<metaIn->general->line_count; ii++) {
	get_float_line(fpIn, metaIn, ii, floatAmpBuf);
	for (kk=0; kk<metaIn->general->sample_count; kk++)
	  ieee_big32(floatAmpBuf[kk]);
	put_float_line(fpOut, metaOut, ii, floatAmpBuf);
      asfLineMeter(ii, metaIn->general->line_count);
      }
      FCLOSE(fpIn);
      FCLOSE(fpOut);
      meta_write(metaOut, outName);
      FREE(floatAmpBuf);
      FREE(outName);
      meta_free(metaIn);
      meta_free(metaOut);
      FREE(insar_params);
    }

    // Ground range correlation image
    if (strcmp_case(type, "InSAR") == 0 &&
	strcmp_case(product[pp], "COR_GRD") == 0) {
      get_uavsar_file_names(inFileName, INSAR_COR_GRD, &dataName, &element,
			    &dataType, &nBands);
      if (!nBands) {
	asfPrintWarning("Ground range correlation image does not exist. "
			"Will skip the ingest.\n");
	continue;
      }
      insar_params =
	read_uavsar_insar_params(inFileName, INSAR_COR_GRD);
      metaIn = uavsar_insar2meta(insar_params);
      metaOut = uavsar_insar2meta(insar_params);
      ns = metaOut->general->sample_count;
      nn = 0;
      floatAmpBuf = (float *) MALLOC(sizeof(float)*ns);
      if (multi)
	outName = appendToBasename(outBaseName, "_cor_grd.img");
      else
	outName = appendExt(outBaseName, ".img");
      asfPrintStatus("\nGround range correlation image:\n");
      char *filename = get_filename(dataName[nn]);
      asfPrintStatus("Ingesting %s ...\n", filename);
      FREE(filename);
      fpIn = FOPEN(dataName[nn], "rb");
      fpOut = FOPEN(outName, "wb");
      strcpy(metaOut->general->bands, "COHERENCE");
      for (ii=0; ii<metaIn->general->line_count; ii++) {
	get_float_line(fpIn, metaIn, ii, floatAmpBuf);
	for (kk=0; kk<metaIn->general->sample_count; kk++)
	  ieee_big32(floatAmpBuf[kk]);
	put_float_line(fpOut, metaOut, ii, floatAmpBuf);
	asfLineMeter(ii, metaIn->general->line_count);
      }
      FCLOSE(fpIn);
      FCLOSE(fpOut);
      meta_write(metaOut, outName);
      FREE(floatAmpBuf);
      FREE(outName);
      meta_free(metaIn);
      meta_free(metaOut);
      FREE(insar_params);
    }

    // Ground range amplitude images
    if (strcmp_case(type, "InSAR") == 0 &&
	strcmp_case(product[pp], "AMP_GRD") == 0) {
      get_uavsar_file_names(inFileName, INSAR_AMP_GRD, &dataName, &element,
			    &dataType, &nBands);
      if (!nBands) {
	asfPrintWarning("Ground range amplitude images does not exist. "
			"Will skip the ingest.\n");
	continue;
      }
      insar_params =
	read_uavsar_insar_params(inFileName, INSAR_AMP_GRD);
      metaIn = uavsar_insar2meta(insar_params);
      metaOut = uavsar_insar2meta(insar_params);
      metaOut->general->band_count = 2;
      ns = metaOut->general->sample_count;
      nn = 0;
      floatAmpBuf = (float *) MALLOC(sizeof(float)*ns);
      outName = (char *) MALLOC(sizeof(char)*(strlen(outBaseName)+15));
      if (multi)
	outName = appendToBasename(outBaseName, "_amp_grd.img");
      else
	outName = appendExt(outBaseName, ".img");
      fpOut = FOPEN(outName, "wb");
      strcpy(metaOut->general->bands, "AMP1,AMP2");
      asfPrintStatus("\nGround range amplitude images:\n");
      for (nn=0; nn<nBands; nn++) {
        char *filename = get_filename(dataName[nn]);
        asfPrintStatus("Ingesting %s ...\n", filename);
        FREE(filename);
	fpIn = FOPEN(dataName[nn], "rb");
	for (ii=0; ii<metaIn->general->line_count; ii++) {
	  get_float_line(fpIn, metaIn, ii, floatAmpBuf);
	  for (kk=0; kk<metaIn->general->sample_count; kk++)
	    ieee_big32(floatAmpBuf[kk]);
	  put_band_float_line(fpOut, metaOut, nn, ii, floatAmpBuf);
	  asfLineMeter(ii, metaIn->general->line_count);
	}
	FCLOSE(fpIn);
      }
      FCLOSE(fpOut);
      meta_write(metaOut, outName);
      FREE(floatAmpBuf);
      FREE(outName);
      meta_free(metaIn);
      meta_free(metaOut);
      FREE(insar_params);
    }

    // Ground range digital elevation model
    if (strcmp_case(type, "InSAR") == 0 &&
	strcmp_case(product[pp], "HGT_GRD") == 0) {
      get_uavsar_file_names(inFileName, INSAR_HGT_GRD, &dataName, &element,
			    &dataType, &nBands);
      if (!nBands) {
	asfPrintWarning("Ground range digital elevation model does not exist. "
			"Will skip the ingest.\n");
	continue;
      }
      insar_params =
	read_uavsar_insar_params(inFileName, INSAR_HGT_GRD);
      metaIn = uavsar_insar2meta(insar_params);
      metaOut = uavsar_insar2meta(insar_params);
      ns = metaOut->general->sample_count;
      nn = 0;
      floatAmpBuf = (float *) MALLOC(sizeof(float)*ns);
      outName = (char *) MALLOC(sizeof(char)*(strlen(outBaseName)+15));
      if (multi)
	outName = appendToBasename(outBaseName, "_hgt_grd.img");
      else
	outName = appendExt(outBaseName, ".img");
      fpOut = FOPEN(outName, "wb");
      strcpy(metaOut->general->bands, "HEIGHT");
      asfPrintStatus("\nGround range digital elevation model:\n");
      for (nn=0; nn<nBands; nn++) {
        char *filename = get_filename(dataName[nn]);
        asfPrintStatus("Ingesting %s ...\n", filename);
        FREE(filename);
	fpIn = FOPEN(dataName[nn], "rb");
	for (ii=0; ii<metaIn->general->line_count; ii++) {
	  get_float_line(fpIn, metaIn, ii, floatAmpBuf);
	  for (kk=0; kk<metaIn->general->sample_count; kk++)
	    ieee_big32(floatAmpBuf[kk]);
	  put_band_float_line(fpOut, metaOut, nn, ii, floatAmpBuf);
	  asfLineMeter(ii, metaIn->general->line_count);
	}
	FCLOSE(fpIn);
      }
      FCLOSE(fpOut);
      meta_write(metaOut, outName);
      FREE(floatAmpBuf);
      FREE(outName);
      meta_free(metaIn);
      meta_free(metaOut);
      FREE(insar_params);
    }

    // Slant range interferogram
    if (strcmp_case(type, "InSAR") == 0 &&
	strcmp_case(product[pp], "INT") == 0) {
      get_uavsar_file_names(inFileName, INSAR_INT, &dataName, &element,
			    &dataType, &nBands);
      if (!nBands) {
	asfPrintWarning("Slant range interferogram does not exist. "
			"Will skip the ingest.\n");
	continue;
      }
      insar_params =
	read_uavsar_insar_params(inFileName, INSAR_INT);
      metaIn = uavsar_insar2meta(insar_params);
      metaOut = uavsar_insar2meta(insar_params);
      ns = metaIn->general->sample_count;
      nn = 0;
      floatAmp = (float *) CALLOC(ns, sizeof(float));
      floatPhase = (float *) CALLOC(ns, sizeof(float));
      floatComplexBuf = (float *) CALLOC(2*ns, sizeof(float));
      outName = (char *) MALLOC(sizeof(char)*(strlen(outBaseName)+15));
      if (multi)
	outName = appendToBasename(outBaseName, "_int.img");
      else
	outName = appendExt(outBaseName, ".img");
      metaOut->general->band_count = 2;
      if (multi)
	outName = appendToBasename(outBaseName, "_int.img");
      else
	outName = appendExt(outBaseName, ".img");
      asfPrintStatus("\nSlant range interferogram:\n");
      fpOut = FOPEN(outName, "wb");
      char *filename = get_filename(dataName[nn]);
      asfPrintStatus("Ingesting %s ...\n", filename);
      FREE(filename);
      fpIn = FOPEN(dataName[nn], "rb");
      sprintf(metaOut->general->bands, "INTERFEROGRAM_AMP,INTERFEROGRAM_PHASE");
      for (ii=0; ii<metaIn->general->line_count; ii++) {
	metaIn->general->sample_count = 2*ns;
	get_float_line(fpIn, metaIn, ii, floatComplexBuf);
	for (kk=0; kk<ns; kk++) {
	  re = floatComplexBuf[kk*2];
	  im = floatComplexBuf[kk*2+1];
	  ieee_big32(re);
	  ieee_big32(im);
	  floatAmp[kk] = hypot(re,im);
	  floatPhase[kk] = atan2_check(im,re);
	}
	put_band_float_line(fpOut, metaOut, 0, ii, floatAmp);
	put_band_float_line(fpOut, metaOut, 1, ii, floatPhase);
	asfLineMeter(ii, metaIn->general->line_count);
      }
      FCLOSE(fpIn);
      FCLOSE(fpOut);
      meta_write(metaOut, outName);
      FREE(floatAmp);
      FREE(floatPhase);
      FREE(floatComplexBuf);
      FREE(outName);
      meta_free(metaIn);
      meta_free(metaOut);
      FREE(insar_params);
    }

    // Slant range unwrapped phase
    if (strcmp_case(type, "InSAR") == 0 &&
	strcmp_case(product[pp], "UNW") == 0) {
      get_uavsar_file_names(inFileName, INSAR_UNW, &dataName, &element,
			    &dataType, &nBands);
      if (!nBands) {
	asfPrintWarning("Slant range unwrapped phase does not exist. "
			"Will skip the ingest.\n");
	continue;
      }
      insar_params =
	read_uavsar_insar_params(inFileName, INSAR_UNW);
      metaIn = uavsar_insar2meta(insar_params);
      metaOut = uavsar_insar2meta(insar_params);
      ns = metaOut->general->sample_count;
      nn = 0;
      floatAmpBuf = (float *) MALLOC(sizeof(float)*ns);
      outName = (char *) MALLOC(sizeof(char)*(strlen(outBaseName)+15));
      if (multi)
	outName = appendToBasename(outBaseName, "_unw.img");
      else
	outName = appendExt(outBaseName, ".img");
      asfPrintStatus("\nSlant range unwrapped phase:\n");
      char *filename = get_filename(dataName[nn]);
      asfPrintStatus("Ingesting %s ...\n", filename);
      FREE(filename);
      fpIn = FOPEN(dataName[nn], "rb");
      fpOut = FOPEN(outName, "wb");
      strcpy(metaOut->general->bands, "UNWRAPPED_PHASE");
      for (ii=0; ii<metaIn->general->line_count; ii++) {
	get_float_line(fpIn, metaIn, ii, floatAmpBuf);
	for (kk=0; kk<metaIn->general->sample_count; kk++)
	  ieee_big32(floatAmpBuf[kk]);
	put_float_line(fpOut, metaOut, ii, floatAmpBuf);
	asfLineMeter(ii, metaIn->general->line_count);
      }
      FCLOSE(fpIn);
      FCLOSE(fpOut);
      meta_write(metaOut, outName);
      FREE(floatAmpBuf);
      FREE(outName);
      meta_free(metaIn);
      meta_free(metaOut);
      FREE(insar_params);
    }

    // Slant range correlation image
    if (strcmp_case(type, "InSAR") == 0 &&
	strcmp_case(product[pp], "COR") == 0) {
      get_uavsar_file_names(inFileName, INSAR_COR, &dataName, &element,
			    &dataType, &nBands);
      if (!nBands) {
	asfPrintWarning("Slant range correlation image does not exist. "
			"Will skip the ingest.\n");
	continue;
      }
      insar_params =
	read_uavsar_insar_params(inFileName, INSAR_COR);
      metaIn = uavsar_insar2meta(insar_params);
      metaOut = uavsar_insar2meta(insar_params);
      ns = metaOut->general->sample_count;
      nn = 0;
      floatAmpBuf = (float *) MALLOC(sizeof(float)*ns);
      outName = (char *) MALLOC(sizeof(char)*(strlen(outBaseName)+15));
      if (multi)
	outName = appendToBasename(outBaseName, "_cor.img");
      else
	outName = appendExt(outBaseName, ".img");
      asfPrintStatus("\nSlant range correlation image:\n");
      char *filename = get_filename(dataName[nn]);
      asfPrintStatus("Ingesting %s ...\n", filename);
      FREE(filename);
      fpIn = FOPEN(dataName[nn], "rb");
      fpOut = FOPEN(outName, "wb");
      strcpy(metaOut->general->bands, "COHERENCE");
      for (ii=0; ii<metaIn->general->line_count; ii++) {
	get_float_line(fpIn, metaIn, ii, floatAmpBuf);
	for (kk=0; kk<metaIn->general->sample_count; kk++)
	  ieee_big32(floatAmpBuf[kk]);
	put_float_line(fpOut, metaOut, ii, floatAmpBuf);
	asfLineMeter(ii, metaIn->general->line_count);
      }
      FCLOSE(fpIn);
      FCLOSE(fpOut);
      meta_write(metaOut, outName);
      FREE(floatAmpBuf);
      FREE(outName);
      meta_free(metaIn);
      meta_free(metaOut);
      FREE(insar_params);
    }

    // Slant range amplitude images
    if (strcmp_case(type, "InSAR") == 0 &&
	strcmp_case(product[pp], "AMP") == 0) {
      get_uavsar_file_names(inFileName, INSAR_AMP, &dataName, &element,
			    &dataType, &nBands);
      if (!nBands) {
	asfPrintWarning("Slant range amplitude images does not exist. "
			"Will skip the ingest.\n");
	continue;
      }
      insar_params =
	read_uavsar_insar_params(inFileName, INSAR_AMP);
      metaIn = uavsar_insar2meta(insar_params);
      metaOut = uavsar_insar2meta(insar_params);
      metaOut->general->band_count = 2;
      ns = metaOut->general->sample_count;
      nn = 0;
      floatAmpBuf = (float *) MALLOC(sizeof(float)*ns);
      outName = (char *) MALLOC(sizeof(char)*(strlen(outBaseName)+15));
      if (multi)
	outName = appendToBasename(outBaseName, "_amp.img");
      else
	outName = appendExt(outBaseName, ".img");
      fpOut = FOPEN(outName, "wb");
      strcpy(metaOut->general->bands, "AMP1,AMP2");
      asfPrintStatus("\nSlant range amplitude images:\n");
      for (nn=0; nn<nBands; nn++) {
        char *filename = get_filename(dataName[nn]);
        asfPrintStatus("Ingesting %s ...\n", filename);
        FREE(filename);
	fpIn = FOPEN(dataName[nn], "rb");
	for (ii=0; ii<metaIn->general->line_count; ii++) {
	  get_float_line(fpIn, metaIn, ii, floatAmpBuf);
	  for (kk=0; kk<metaIn->general->sample_count; kk++)
	    ieee_big32(floatAmpBuf[kk]);
	  put_band_float_line(fpOut, metaOut, nn, ii, floatAmpBuf);
	  asfLineMeter(ii, metaIn->general->line_count);
	}
	FCLOSE(fpIn);
      }
      FCLOSE(fpOut);
      meta_write(metaOut, outName);
      FREE(floatAmpBuf);
      FREE(outName);
      meta_free(metaIn);
      meta_free(metaOut);
      FREE(insar_params);
    }

    // PolSAR data
    // Single look complex data
    if (strcmp_case(type, "PolSAR") == 0 &&
	strcmp_case(product[pp], "SLC") == 0) {
      asfPrintWarning("Ingest of SLC data is currently not supported!\n");
      /*
      get_uavsar_file_names(inFileName, POLSAR_SLC, &dataName, &element,
			    &dataType, &nBands);
      polsar_params =
	read_uavsar_polsar_params(inFileName, POLSAR_SLC);
      metaIn = uavsar_polsar2meta(polsar_params);
      metaOut = uavsar_polsar2meta(polsar_params);
      if (multi)
	outName = appendToBasename(outBaseName, "_slc.img");
      else
	outName = appendExt(outBaseName, ".img");
      for (ii=0; ii<nBands; ii++)
	printf("file: %s\n", dataName[ii]);
      //meta_write(metaOut, outName);
      meta_free(metaIn);
      meta_free(metaOut);
      FREE(polsar_params);
      */
    }

    // Multilooked data
    if (strcmp_case(type, "PolSAR") == 0  &&
	strcmp_case(product[pp], "MLC") == 0) {
      get_uavsar_file_names(inFileName, POLSAR_MLC, &dataName, &element,
			    &dataType, &nBands);
      if (!nBands) {
	asfPrintWarning("Slant range multilooked data do not exist. "
			"Will skip the ingest.\n");
	continue;
      }
      polsar_params =
	read_uavsar_polsar_params(inFileName, POLSAR_MLC);
      metaIn = uavsar_polsar2meta(polsar_params);
      metaOut = uavsar_polsar2meta(polsar_params);
      ns = metaIn->general->sample_count;
      amp = (float *) MALLOC(sizeof(float)*ns);
      floatAmp = (float *) MALLOC(sizeof(float)*ns);
      floatAmpBuf = (float *) MALLOC(sizeof(float)*ns);
      floatComplexReal = (float *) MALLOC(sizeof(float)*ns);
      floatComplexImag = (float *) MALLOC(sizeof(float)*ns);
      floatComplexBuf = (float *) MALLOC(sizeof(float)*2*ns);
      outName = (char *) MALLOC(sizeof(char)*(strlen(outBaseName)+15));
      metaOut->general->band_count = ll = 1;
      if (multi)
	outName = appendToBasename(outBaseName, "_mlc.img");
      else
	outName = appendExt(outBaseName, ".img");
      asfPrintStatus("\nMultilooked data:\n");
      fpOut = FOPEN(outName, "wb");
      for (nn=0; nn<nBands; nn++) {
        char *filename = get_filename(dataName[nn]);
        asfPrintStatus("Ingesting %s ...\n", filename);
        FREE(filename);
	if (dataType[nn] == 0)
	  metaOut->general->band_count += 1;
	else
	  metaOut->general->band_count += 2;
	fpIn = FOPEN(dataName[nn], "rb");
	if (nn == 0)
	  sprintf(metaOut->general->bands, "AMP,%s", element[0]);
	else {
	  if (dataType[nn])
	    sprintf(tmp, ",%s_real,%s_imag", element[nn], element[nn]);
	  else
	    sprintf(tmp, ",%s", element[nn]);
	  strcat(metaOut->general->bands, tmp);
	}
	for (ii=0; ii<metaIn->general->line_count; ii++) {
	  if (dataType[nn]) {
	    metaIn->general->sample_count = 2*ns;
	    get_float_line(fpIn, metaIn, ii, floatComplexBuf);
	    for (kk=0; kk<ns; kk++) {
	      floatComplexReal[kk] = floatComplexBuf[kk*2];
	      floatComplexImag[kk] = floatComplexBuf[kk*2+1];
	      ieee_big32(floatComplexReal[kk]);
	      ieee_big32(floatComplexImag[kk]);
	    }
	    put_band_float_line(fpOut, metaOut, ll, ii, floatComplexReal);
	    put_band_float_line(fpOut, metaOut, ll+1, ii, floatComplexImag);
	  }
	  else {
	    metaIn->general->sample_count = ns;
	    get_float_line(fpIn, metaIn, ii, floatAmpBuf);
	    if (nn == 0) {
	      for (kk=0; kk<ns; kk++) {
		ieee_big32(floatAmpBuf[kk]);
		floatAmp[kk] = sqrt(floatAmpBuf[kk]);
	      }
	    }
	    else {
	      for (kk=0; kk<ns; kk++)
		ieee_big32(floatAmpBuf[kk]);
	    }
	    put_band_float_line(fpOut, metaOut, ll, ii, floatAmpBuf);
	  }
	  asfLineMeter(ii, metaIn->general->line_count);
	}
	if (dataType[nn])
	  ll += 2;
	else
	  ll++;
	FCLOSE(fpIn);
      }
      FCLOSE(fpOut);
      meta_write(metaOut, outName);
      FREE(amp);
      FREE(floatAmp);
      FREE(floatAmpBuf);
      FREE(floatComplexReal);
      FREE(floatComplexImag);
      FREE(floatComplexBuf);
      FREE(outName);
      meta_free(metaIn);
      meta_free(metaOut);
      FREE(polsar_params);
    }

    // Compressed Stokes matrix
    if (strcmp_case(type, "PolSAR") == 0 &&
	strcmp_case(product[pp], "DAT") == 0) {
      get_uavsar_file_names(inFileName, POLSAR_DAT, &dataName, &element,
			    &dataType, &nBands);
      if (!nBands) {
	asfPrintWarning("Slant range Stokes matrix does not exist. "
			"Will skip the ingest.\n");
	continue;
      }
      polsar_params =
	read_uavsar_polsar_params(inFileName, POLSAR_DAT);
      metaIn = uavsar_polsar2meta(polsar_params);
      metaOut = uavsar_polsar2meta(polsar_params);
      outName = (char *) MALLOC(sizeof(char)*(strlen(outBaseName)+15));
      if (multi)
	outName = appendToBasename(outBaseName, "_dat.img");
      else
	outName = appendExt(outBaseName, ".img");
      metaOut->general->band_count = 9;
      asfPrintStatus("\nCompressed Stokes matrix:\n");
      char *filename = get_filename(dataName[0]);
      asfPrintStatus("Ingesting %s ...\n", filename);
      FREE(filename);
      if (radiometry == r_AMP)
	strcpy(metaOut->general->bands,
	       "AMP,AMP_HH,PHASE_HH,AMP_HV,PHASE_HV,AMP_VH,PHASE_VH,"	\
	       "AMP_VV,PHASE_VV");
      else if (radiometry == r_SIGMA)
	strcpy(metaOut->general->bands,
	       "AMP,SIGMA-AMP-HH,SIGMA-PHASE-HH,SIGMA-AMP-HV,SIGMA-PHASE-HV," \
	       "SIGMA-AMP-VH,SIGMA-PHASE-VH,SIGMA-AMP-VV,SIGMA-PHASE-VV");
      else if (radiometry == r_SIGMA_DB)
	strcpy(metaOut->general->bands,
	       "AMP,SIGMA_DB-AMP-HH,SIGMA_DB-PHASE-HH,SIGMA_DB-AMP-HV,"	\
	       "SIGMA_DB-PHASE-HV,SIGMA_DB-AMP-VH,SIGMA_DB-PHASE-VH,"	\
	       "SIGMA_DB-AMP-VV,SIGMA_DB-PHASE-VV");
      int ns = metaOut->general->sample_count;
      float total_power, ysca, amp, phase;
      complexFloat cpx;
      float *power = (float *) MALLOC(sizeof(float)*ns);
      float *shh_amp = (float *) MALLOC(sizeof(float)*ns);
      float *shh_phase = (float *) MALLOC(sizeof(float)*ns);
      float *shv_amp = (float *) MALLOC(sizeof(float)*ns);
      float *shv_phase = (float *) MALLOC(sizeof(float)*ns);
      float *svh_amp = (float *) MALLOC(sizeof(float)*ns);
      float *svh_phase = (float *) MALLOC(sizeof(float)*ns);
      float *svv_amp = (float *) MALLOC(sizeof(float)*ns);
      float *svv_phase = (float *) MALLOC(sizeof(float)*ns);
      char *byteBuf = (char *) MALLOC(sizeof(char)*10);
      fpIn = FOPEN(dataName[0], "rb");
      fpOut = FOPEN(outName, "wb");
      for (ii=0; ii<metaOut->general->line_count; ii++) {
	for (kk=0; kk<metaOut->general->sample_count; kk++) {
	  FREAD(byteBuf, sizeof(char), 10, fpIn);
          //float m11, m12, m13, m14, m22, m23, m24, m33, m34, m44;
	  // Scale is always 1.0 according to Bruce Chapman
	  //m11 = ((float)byteBuf[1]/254.0 + 1.5) * pow(2, byteBuf[0]);
	  //m12 = (float)byteBuf[2] * m11 / 127.0;
	  //m13 = sign(byteBuf[3]) * SQR((float)byteBuf[3] / 127.0) * m11;
	  //m14 = sign(byteBuf[4]) * SQR((float)byteBuf[4] / 127.0) * m11;
	  //m23 = sign(byteBuf[5]) * SQR((float)byteBuf[5] / 127.0) * m11;
	  //m24 = sign(byteBuf[6]) * SQR((float)byteBuf[6] / 127.0) * m11;
	  //m33 = (float)byteBuf[7] * m11 / 127.0;
	  //m34 = (float)byteBuf[8] * m11 / 127.0;
	  //m44 = (float)byteBuf[9] * m11 / 127.0;
	  //m22 = 1 - m33 -m44;
	  total_power =
	    ((float)byteBuf[1]/254.0 + 1.5) * pow(2, byteBuf[0]);
	  ysca = 2.0 * sqrt(total_power);
	  power[kk] = sqrt(total_power);
	  cpx.real = (float)byteBuf[2] * ysca / 127.0;
	  cpx.imag = (float)byteBuf[3] * ysca / 127.0;
	  amp = sqrt(cpx.real*cpx.real + cpx.imag*cpx.imag);
	  phase = atan2(cpx.imag, cpx.real);
	  if (radiometry == r_AMP) {
	    shh_amp[kk] = amp;
	    shh_phase[kk] = phase;
	  }
	  else if (radiometry == r_SIGMA) {
	    shh_amp[kk] = amp*amp;
	    shh_phase[kk] = phase;
	  }
	  else if (radiometry == r_SIGMA_DB) {
	    shh_amp[kk] = amp;
	    shh_phase[kk] = phase;
	  }
	  cpx.real = (float)byteBuf[4] * ysca / 127.0;
	  cpx.imag = (float)byteBuf[5] * ysca / 127.0;
	  amp = sqrt(cpx.real*cpx.real + cpx.imag*cpx.imag);
	  phase = atan2(cpx.imag, cpx.real);
	  if (radiometry == r_AMP) {
	    shv_amp[kk] = amp;
	    shv_phase[kk] = phase;
	  }
	  else if (radiometry == r_SIGMA) {
	    shv_amp[kk] = amp*amp;
	    shv_phase[kk] = phase;
	  }
	  else if (radiometry == r_SIGMA_DB) {
	    shv_amp[kk] = amp;
	    shv_phase[kk] = phase;
	  }
	  cpx.real = (float)byteBuf[6] * ysca / 127.0;
	  cpx.imag = (float)byteBuf[7] * ysca / 127.0;
	  amp = sqrt(cpx.real*cpx.real + cpx.imag*cpx.imag);
	  phase = atan2(cpx.imag, cpx.real);
	  if (radiometry == r_AMP) {
	    svh_amp[kk] = amp;
	    svh_phase[kk] = phase;
	  }
	  else if (radiometry == r_SIGMA) {
	    svh_amp[kk] = amp*amp;
	    svh_phase[kk] = phase;
	  }
	  else if (radiometry == r_SIGMA_DB) {
	    svh_amp[kk] = amp;
	    svh_phase[kk] = phase;
	  }
	  cpx.real = (float)byteBuf[8] * ysca / 127.0;
	  cpx.imag = (float)byteBuf[9] * ysca / 127.0;
	  amp = sqrt(cpx.real*cpx.real + cpx.imag*cpx.imag);
	  phase = atan2(cpx.imag, cpx.real);
	  if (radiometry == r_AMP) {
	    svv_amp[kk] = amp;
	    svv_phase[kk] = phase;
	  }
	  else if (radiometry == r_SIGMA) {
	    svv_amp[kk] = amp*amp;
	    svv_phase[kk] = phase;
	  }
	  else if (radiometry == r_SIGMA_DB) {
	    svv_amp[kk] = amp;
	    svv_phase[kk] = phase;
	  }
	}
	put_band_float_line(fpOut, metaOut, 0, ii, power);
	put_band_float_line(fpOut, metaOut, 1, ii, shh_amp);
	put_band_float_line(fpOut, metaOut, 2, ii, shh_phase);
	put_band_float_line(fpOut, metaOut, 3, ii, shv_amp);
	put_band_float_line(fpOut, metaOut, 4, ii, shv_phase);
	put_band_float_line(fpOut, metaOut, 5, ii, svh_amp);
	put_band_float_line(fpOut, metaOut, 6, ii, svh_phase);
	put_band_float_line(fpOut, metaOut, 7, ii, svv_amp);
	put_band_float_line(fpOut, metaOut, 8, ii, svv_phase);
	asfLineMeter(ii, metaOut->general->line_count);
      }
      FCLOSE(fpIn);
      FCLOSE(fpOut);
      FREE(power);
      FREE(shh_amp);
      FREE(shh_phase);
      FREE(shv_amp);
      FREE(shv_phase);
      FREE(svh_amp);
      FREE(svh_phase);
      FREE(svv_amp);
      FREE(svv_phase);
      meta_write(metaOut, outName);
      meta_free(metaIn);
      meta_free(metaOut);
      FREE(outName);
      FREE(polsar_params);
    }

    // Ground range projected data
    if (strcmp_case(type, "PolSAR") == 0 &&
	strcmp_case(product[pp], "GRD") == 0) {
      get_uavsar_file_names(inFileName, POLSAR_GRD, &dataName, &element,
			    &dataType, &nBands);
      if (!nBands) {
	asfPrintWarning("Ground range projected data do not exist. "
			"Will skip the ingest.\n");
	continue;
      }
      polsar_params =
	read_uavsar_polsar_params(inFileName, POLSAR_GRD);
      metaIn = uavsar_polsar2meta(polsar_params);
      metaOut = uavsar_polsar2meta(polsar_params);
      ns = metaIn->general->sample_count;
      floatAmpBuf = (float *) CALLOC(ns, sizeof(float));
      floatComplexReal = (float *) CALLOC(ns, sizeof(float));
      floatComplexImag = (float *) CALLOC(ns, sizeof(float));
      floatComplexBuf = (float *) CALLOC(2*ns, sizeof(float));
      outName = (char *) MALLOC(sizeof(char)*(strlen(outBaseName)+15));
      metaOut->general->band_count = ll = 0;
      if (multi)
	outName = appendToBasename(outBaseName, "_grd.img");
      else
	outName = appendExt(outBaseName, ".img");
      asfPrintStatus("\nGround range projected data:\n");
      fpOut = FOPEN(outName, "wb");
      for (nn=0; nn<nBands; nn++) {
        char *filename = get_filename(dataName[nn]);
        asfPrintStatus("Ingesting %s ...\n", filename);
        FREE(filename);
	if (dataType[nn])
	  metaOut->general->band_count += 2;
	else
	  metaOut->general->band_count += 1;
	fpIn = FOPEN(dataName[nn], "rb");
	if (nn == 0)
	  sprintf(metaOut->general->bands, "%s", element[0]);
	else {
	  if (dataType[nn])
	    sprintf(tmp, ",%s_real,%s_imag", element[nn], element[nn]);
	  else
	    sprintf(tmp, ",%s", element[nn]);
	  strcat(metaOut->general->bands, tmp);
	}
	if (dataType[nn]) {
	  for (ii=0; ii<metaIn->general->line_count; ii++) {
	    metaIn->general->sample_count = 2*ns;
	    get_float_line(fpIn, metaIn, ii, floatComplexBuf);
	    for (kk=0; kk<ns; kk++) {
	      floatComplexReal[kk] = floatComplexBuf[kk*2];
	      floatComplexImag[kk] = floatComplexBuf[kk*2+1];
	      ieee_big32(floatComplexReal[kk]);
	      ieee_big32(floatComplexImag[kk]);
	    }
	    put_band_float_line(fpOut, metaOut, ll, ii, floatComplexReal);
	    put_band_float_line(fpOut, metaOut, ll+1, ii, floatComplexImag);
	    asfLineMeter(ii, metaIn->general->line_count);
	  }
	}
	else {
	  for (ii=0; ii<metaIn->general->line_count; ii++) {
	    metaIn->general->sample_count = ns;
	    get_float_line(fpIn, metaIn, ii, floatAmpBuf);
	    for (kk=0; kk<ns; kk++)
	      ieee_big32(floatAmpBuf[kk]);
	    put_band_float_line(fpOut, metaOut, ll, ii, floatAmpBuf);
	    asfLineMeter(ii, metaIn->general->line_count);
	  }
	}
	FCLOSE(fpIn);
	if (dataType[nn])
	  ll += 2;
	else
	  ll++;
      }
      FCLOSE(fpOut);
      meta_write(metaOut, outName);
      FREE(floatAmpBuf);
      FREE(floatComplexReal);
      FREE(floatComplexImag);
      FREE(floatComplexBuf);
      FREE(outName);
      meta_free(metaIn);
      meta_free(metaOut);
      FREE(polsar_params);
    }

    // Digital elevation model
    if (strcmp_case(type, "PolSAR") == 0 &&
	strcmp_case(product[pp], "HGT") == 0) {
      get_uavsar_file_names(inFileName, POLSAR_HGT, &dataName, &element,
			    &dataType, &nBands);
      if (!nBands) {
	asfPrintWarning("Digital elevation model does not exist. "
			"Will skip the ingest.\n");
	continue;
      }
      polsar_params =
	read_uavsar_polsar_params(inFileName, POLSAR_HGT);
      metaIn = uavsar_polsar2meta(polsar_params);
      metaOut = uavsar_polsar2meta(polsar_params);
      ns = metaOut->general->sample_count;
      floatAmpBuf = (float *) MALLOC(sizeof(float)*ns);
      outName = (char *) MALLOC(sizeof(char)*(strlen(outBaseName)+15));
      if (multi)
	outName = appendToBasename(outBaseName, "_hgt.img");
      else
	outName = appendExt(outBaseName, ".img");
      asfPrintStatus("\nDigital elevation model:\n");
      nn=0;
      char *filename = get_filename(dataName[nn]);
      asfPrintStatus("Ingesting %s ...\n", filename);
      FREE(filename);
      fpIn = FOPEN(dataName[nn], "rb");
      fpOut = FOPEN(outName, "wb");
      strcpy(metaOut->general->bands, "HEIGHT");
      for (ii=0; ii<metaIn->general->line_count; ii++) {
	get_float_line(fpIn, metaIn, ii, floatAmpBuf);
	for (kk=0; kk<metaIn->general->sample_count; kk++)
	  ieee_big32(floatAmpBuf[kk]);
	put_float_line(fpOut, metaOut, ii, floatAmpBuf);
	asfLineMeter(ii, metaIn->general->line_count);
      }
      FCLOSE(fpIn);
      FCLOSE(fpOut);
      meta_write(metaOut, outName);
      FREE(floatAmpBuf);
      FREE(outName);
      meta_free(metaIn);
      meta_free(metaOut);
      FREE(polsar_params);
    }
  }
  if (strcmp_case(type, "PolSAR") == 0)
    product_count = 5;
  else if (strcmp_case(type, "InSAR") == 0)
    product_count = 9;
  for (pp=0; pp<product_count; pp++)
    FREE(product[pp]);
  FREE(product);
  FREE(type);
}

void read_meta_uavsar(const char *inFileName, const char *outBaseName)
{
  int pp, product_count, *dataType, nBands;
  char *type = check_data_type(inFileName);
  char **product = get_uavsar_products("ALL", type, &product_count);
  char **dataName, **element;
  char *outName = (char *) MALLOC(sizeof(char)*255);
  uavsar_polsar *polsar_params;
  uavsar_insar *insar_params;
  meta_parameters *meta;

  for (pp = 0; pp < product_count; pp++) {

    // InSAR
    if (strcmp_case(type, "InSAR") == 0 &&
	strcmp_case(product[pp], "INT_GRD") == 0) {
      get_uavsar_file_names(inFileName, INSAR_INT_GRD, &dataName, &element,
			    &dataType, &nBands);
      if (!nBands)
	continue;
      insar_params =
	read_uavsar_insar_params(inFileName, INSAR_INT_GRD);
      meta = uavsar_insar2meta(insar_params);
      sprintf(outName, "%s_int_grd.xml", outBaseName);
      meta_write_xml(meta, outName);
      FREE(outName);
      meta_free(meta);
      FREE(insar_params);
      FREE(dataType);
    }
    if (strcmp_case(type, "InSAR") == 0 &&
	strcmp_case(product[pp], "UNW_GRD") == 0) {
      get_uavsar_file_names(inFileName, INSAR_UNW_GRD, &dataName, &element,
			    &dataType, &nBands);
      if (!nBands)
	continue;
      insar_params =
	read_uavsar_insar_params(inFileName, INSAR_UNW_GRD);
      meta = uavsar_insar2meta(insar_params);
      sprintf(outName, "%s_unw_grd.xml", outBaseName);
      meta_write_xml(meta, outName);
      FREE(outName);
      meta_free(meta);
      FREE(insar_params);
      FREE(dataType);
    }
    if (strcmp_case(type, "InSAR") == 0 &&
	strcmp_case(product[pp], "COR_GRD") == 0) {
      get_uavsar_file_names(inFileName, INSAR_COR_GRD, &dataName, &element,
			    &dataType, &nBands);
      if (!nBands)
	continue;
      insar_params =
	read_uavsar_insar_params(inFileName, INSAR_COR_GRD);
      meta = uavsar_insar2meta(insar_params);
      sprintf(outName, "%s_cor_grd.xml", outBaseName);
      meta_write_xml(meta, outName);
      FREE(outName);
      meta_free(meta);
      FREE(insar_params);
      FREE(dataType);
    }
    if (strcmp_case(type, "InSAR") == 0 &&
	strcmp_case(product[pp], "AMP_GRD") == 0) {
      get_uavsar_file_names(inFileName, INSAR_AMP_GRD, &dataName, &element,
			    &dataType, &nBands);
      if (!nBands)
	continue;
      insar_params =
	read_uavsar_insar_params(inFileName, INSAR_AMP_GRD);
      meta = uavsar_insar2meta(insar_params);
      sprintf(outName, "%s_amp_grd.xml", outBaseName);
      meta_write_xml(meta, outName);
      FREE(outName);
      meta_free(meta);
      FREE(insar_params);
      FREE(dataType);
    }
    if (strcmp_case(type, "InSAR") == 0 &&
	strcmp_case(product[pp], "HGT_GRD") == 0) {
      get_uavsar_file_names(inFileName, INSAR_HGT_GRD, &dataName, &element,
			    &dataType, &nBands);
      if (!nBands)
	continue;
      insar_params =
	read_uavsar_insar_params(inFileName, INSAR_HGT_GRD);
      meta = uavsar_insar2meta(insar_params);
      sprintf(outName, "%s_hgt_grd.xml", outBaseName);
      meta_write_xml(meta, outName);
      FREE(outName);
      meta_free(meta);
      FREE(insar_params);
      FREE(dataType);
    }
    if (strcmp_case(type, "InSAR") == 0 &&
	strcmp_case(product[pp], "INT") == 0) {
      get_uavsar_file_names(inFileName, INSAR_INT, &dataName, &element,
			    &dataType, &nBands);
      if (!nBands)
	continue;
      insar_params =
	read_uavsar_insar_params(inFileName, INSAR_INT);
      meta = uavsar_insar2meta(insar_params);
      sprintf(outName, "%s_int.xml", outBaseName);
      meta_write_xml(meta, outName);
      FREE(outName);
      meta_free(meta);
      FREE(insar_params);
      FREE(dataType);
    }
    if (strcmp_case(type, "InSAR") == 0 &&
	strcmp_case(product[pp], "UNW") == 0) {
      get_uavsar_file_names(inFileName, INSAR_UNW, &dataName, &element,
			    &dataType, &nBands);
      if (!nBands)
	continue;
      insar_params =
	read_uavsar_insar_params(inFileName, INSAR_UNW);
      meta = uavsar_insar2meta(insar_params);
      sprintf(outName, "%s_unw.xml", outBaseName);
      meta_write_xml(meta, outName);
      FREE(outName);
      meta_free(meta);
      FREE(insar_params);
      FREE(dataType);
    }
    if (strcmp_case(type, "InSAR") == 0 &&
	strcmp_case(product[pp], "COR") == 0) {
      get_uavsar_file_names(inFileName, INSAR_COR, &dataName, &element,
			    &dataType, &nBands);
      if (!nBands)
	continue;
      insar_params =
	read_uavsar_insar_params(inFileName, INSAR_COR);
      meta = uavsar_insar2meta(insar_params);
      sprintf(outName, "%s_cor.xml", outBaseName);
      meta_write_xml(meta, outName);
      FREE(outName);
      meta_free(meta);
      FREE(insar_params);
      FREE(dataType);
    }
    if (strcmp_case(type, "InSAR") == 0 &&
	strcmp_case(product[pp], "AMP") == 0) {
      get_uavsar_file_names(inFileName, INSAR_AMP, &dataName, &element,
			    &dataType, &nBands);
      if (!nBands)
	continue;
      insar_params =
	read_uavsar_insar_params(inFileName, INSAR_AMP);
      meta = uavsar_insar2meta(insar_params);
      sprintf(outName, "%s_amp.xml", outBaseName);
      meta_write_xml(meta, outName);
      FREE(outName);
      meta_free(meta);
      FREE(insar_params);
      FREE(dataType);
    }
    // PolSAR
    if (strcmp_case(type, "PolSAR") == 0  &&
	strcmp_case(product[pp], "MLC") == 0) {
      get_uavsar_file_names(inFileName, POLSAR_MLC, &dataName, &element,
			    &dataType, &nBands);
      if (!nBands)
	continue;
      polsar_params =
	read_uavsar_polsar_params(inFileName, POLSAR_MLC);
      meta = uavsar_polsar2meta(polsar_params);
      sprintf(outName, "%s_mlc.xml", outBaseName);
      meta_write_xml(meta, outName);
      FREE(outName);
      meta_free(meta);
      FREE(insar_params);
      FREE(dataType);
    }
    if (strcmp_case(type, "PolSAR") == 0 &&
	strcmp_case(product[pp], "DAT") == 0) {
      get_uavsar_file_names(inFileName, POLSAR_DAT, &dataName, &element,
			    &dataType, &nBands);
      if (!nBands)
	continue;
      polsar_params =
	read_uavsar_polsar_params(inFileName, POLSAR_DAT);
      meta = uavsar_polsar2meta(polsar_params);
      sprintf(outName, "%s_dat.xml", outBaseName);
      meta_write_xml(meta, outName);
      FREE(outName);
      meta_free(meta);
      FREE(insar_params);
      FREE(dataType);
    }
    if (strcmp_case(type, "PolSAR") == 0 &&
	strcmp_case(product[pp], "GRD") == 0) {
      get_uavsar_file_names(inFileName, POLSAR_GRD, &dataName, &element,
			    &dataType, &nBands);
      if (!nBands)
	continue;
      polsar_params =
	read_uavsar_polsar_params(inFileName, POLSAR_GRD);
      meta = uavsar_polsar2meta(polsar_params);
      sprintf(outName, "%s_grd.xml", outBaseName);
      meta_write_xml(meta, outName);
      FREE(outName);
      meta_free(meta);
      FREE(insar_params);
      FREE(dataType);
    }
    if (strcmp_case(type, "PolSAR") == 0 &&
	strcmp_case(product[pp], "HGT") == 0) {
      get_uavsar_file_names(inFileName, POLSAR_HGT, &dataName, &element,
			    &dataType, &nBands);
      if (!nBands)
	continue;
      polsar_params =
	read_uavsar_polsar_params(inFileName, POLSAR_HGT);
      meta = uavsar_polsar2meta(polsar_params);
      sprintf(outName, "%s_hgt.xml", outBaseName);
      meta_write_xml(meta, outName);
      FREE(outName);
      meta_free(meta);
      FREE(insar_params);
      FREE(dataType);
    }
  }
  for (pp=0; pp<product_count; pp++)
    FREE(product[pp]);
  FREE(product);
  FREE(type);
}

max_stars_repo_name: glshort/MapReady

max_issues_repo_name: glshort/MapReady

max_forks_repo_name: glshort/MapReady


text:
/*****************************************************************************/
/* Solvers for the Least squares problem min_x || B - A*x ||                 */
/*****************************************************************************/
#include <lah.h>

/* lah_solveLS
 * ------------
 * Solves the least squares problem by applying the Lapack routine gels
 * which uses QR or LQ factorization.
 * Params:  - A, B, X: Matrices of the LS problem min_x || B - A*x ||
*/
#ifdef HAVE_LAPACK /* Use a LAPACK package to do the heavy work */
#include <lapacke.h>
lah_Return lah_solveLS(koop_mat *A, *B, *X)
{
    lah result = 0;
    gels('N', A->nR, A->nC, NRHS, A->data, LDA, B, LDB, WORK,
                        LWORK, result);

    /* Alternative: used for ill conditioned problems
     * ca. 5 times slower. Could check residuum output to determine which precision to use*/
    /* gelsd () */
    return result;
}
#endif

ext: c


text:
//
// ePDF
//

#pragma once

#include <gsl/gsl_integration.h>

namespace ePDF
{
  struct int_params
  {
    double z;
    int nl;
  };

  /**
   * @brief The "NumericIntegrals" class
   */
  class NumericIntegrals
  {
  public:
    /**
     * @brief The "NumericIntegrals" constructor.
     * @param flv: 0 = singlet, 1 = photon, 2 = nonsinglet
     * @param nl: number of leptons (needed for singlet)
     */
    NumericIntegrals(int const& iflv, int const& nl);

    /**
     * @brief The "NumericIntegrals" destructor.
     */
    ~NumericIntegrals();

    /**
     * @brief Perform the integration
     * @param z: momentum fraction
     */
    double integrate(double const& z) const;

  private:
    int    const _iflv;
    int    const _nl;
    double const _epsabs;
    double const _epsrel;
    gsl_integration_workspace *_w;
  };
}

ext: h

max_stars_repo_name: vbertone/ePDF

max_issues_repo_name: gstagnit/ePDF

max_forks_repo_name: gstagnit/ePDF


text:
#pragma once

#include <memory>
#include <gsl/span>

#include "halley/time/halleytime.h"
#include "../session/network_session.h"
#include "entity_network_remote_peer.h"
#include "halley/bytes/serialization_dictionary.h"

namespace Halley {
	class EntityFactory;
	class Resources;
	class World;
	class NetworkSession;

	class EntitySessionSharedData : public SharedData {
	public:

	};

	class EntityClientSharedData : public SharedData {
	public:
		std::optional<Rect4i> viewRect;

		void serialize(Serializer& s) const override;
		void deserialize(Deserializer& s) override;
	};

	class EntityNetworkSession : NetworkSession::IListener, NetworkSession::ISharedDataHandler {
    public:
		class IEntityNetworkSessionListener {
		public:
			virtual ~IEntityNetworkSessionListener() = default;
			virtual void onRemoteEntityCreated(EntityRef entity, NetworkSession::PeerId peerId) {}
			virtual void onPreSendDelta(EntityDataDelta& delta) {}
			virtual bool isEntityInView(EntityRef entity, const EntityClientSharedData& clientData) = 0;
		};

		EntityNetworkSession(std::shared_ptr<NetworkSession> session, Resources& resources, std::set<String> ignoreComponents, IEntityNetworkSessionListener* listener);
		~EntityNetworkSession() override;

		void setWorld(World& world);

		void sendUpdates(Time t, Rect4i viewRect, gsl::span<const std::pair<EntityId, uint8_t>> entityIds); // Takes pairs of entity id and owner peer id
		void receiveUpdates();

		World& getWorld() const;
		EntityFactory& getFactory() const;
		NetworkSession& getSession() const;
		bool hasWorld() const;

		const EntityFactory::SerializationOptions& getEntitySerializationOptions() const;
		const EntityDataDelta::Options& getEntityDeltaOptions() const;
		const SerializerOptions& getByteSerializationOptions() const;

		Time getMinSendInterval() const;

		void onRemoteEntityCreated(EntityRef entity, NetworkSession::PeerId peerId);
		void onPreSendDelta(EntityDataDelta& delta);

		bool isReadyToStart() const;
		bool isEntityInView(EntityRef entity, const EntityClientSharedData& clientData) const;

		std::vector<Rect4i> getRemoteViewPorts() const;

	protected:
		void onStartSession(NetworkSession::PeerId myPeerId) override;
		void onPeerConnected(NetworkSession::PeerId peerId) override;
		void onPeerDisconnected(NetworkSession::PeerId peerId) override;
		std::unique_ptr<SharedData> makeSessionSharedData() override;
		std::unique_ptr<SharedData> makePeerSharedData() override;

	private:
		struct QueuedPacket {
			NetworkSession::PeerId fromPeerId;
			EntityNetworkHeaderType type;
			InboundNetworkPacket packet;
		};

		Resources& resources;
		std::shared_ptr<EntityFactory> factory;
		IEntityNetworkSessionListener* listener = nullptr;

		EntityFactory::SerializationOptions entitySerializationOptions;
		EntityDataDelta::Options deltaOptions;
		SerializerOptions byteSerializationOptions;
		SerializationDictionary serializationDictionary;

		std::shared_ptr<NetworkSession> session;
		std::vector<EntityNetworkRemotePeer> peers;

		std::vector<QueuedPacket> queuedPackets;

		bool readyToStart = false;

		void onReceiveEntityUpdate(NetworkSession::PeerId fromPeerId, EntityNetworkHeaderType type, InboundNetworkPacket packet);
		void onReceiveReady(NetworkSession::PeerId fromPeerId);

		void setupDictionary();
	};
}

max_stars_repo_name: JoelOtter/halley

max_issues_repo_name: JoelOtter/halley

max_forks_repo_name: JoelOtter/halley


text:
/*
* -----------------------------------------------------------------
*  Binary Search Tree Library --- bst_lib.h
*  Version: 1.6180
*  Date: Feb 19, 2010
* -----------------------------------------------------------------
*  Programmer: Americo Barbosa da Cunha Junior
*              americo.cunhajr@gmail.com
* -----------------------------------------------------------------
*  Copyright (c) 2010 by Americo Barbosa da Cunha Junior
*
*  This program is free software: you can redistribute it and/or
*  modify it under the terms of the GNU General Public License as
*  published by the Free Software Foundation, either version 3 of
*  the License, or (at your option) any later version.
*
*  This program is distributed in the hope that it will be useful,
*  but WITHOUT ANY WARRANTY; without even the implied warranty of
*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
*  GNU General Public License for more details.
*
*  A copy of the GNU General Public License is available in
*  LICENSE.txt or http://www.gnu.org/licenses/.
* -----------------------------------------------------------------
*  This is the header file for a library with
*  to operates a binary search tree.
* -----------------------------------------------------------------
*/

#ifndef __BST_LIB_H__
#define __BST_LIB_H__

#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>

#undef  RIGHT
#define RIGHT 1

#undef  LEFT
#define LEFT 0

/*
*------------------------------------------------------------
*   structure of binary search tree leaf
*
*   last update: Feb 2, 2009
*------------------------------------------------------------
*/

typedef struct leaf
{
    gsl_vector *phi;   /* composition */
    gsl_vector *Rphi;  /* reaction mapping */
    gsl_matrix *A;     /* mapping gradient matrix */
    gsl_matrix *L;     /* ellipsoid Cholesky matrix */
} bst_leaf;
/*------------------------------------------------------------*/

/*
*------------------------------------------------------------
*   structure of binary search tree node
*
*   last update: Feb 2, 2009
*------------------------------------------------------------
*/

typedef struct node
{
    gsl_vector *v;          /* cutting plane normal vector */
    double a;               /* scalar a */
    struct node *r_node;    /* right node */
    struct node *l_node;    /* left node */
    struct leaf *r_leaf;    /* right leaf */
    struct leaf *l_leaf;    /* left leaf */
} bst_node;
/*------------------------------------------------------------*/

/*
*------------------------------------------------------------
* function prototypes
*------------------------------------------------------------
*/

bst_leaf *bst_leaf_alloc();

bst_node* bst_node_alloc();

void bst_leaf_free(void **leaf_bl);

void bst_node_free(void **node_bl);

int bst_leaf_set(gsl_vector *phi,
                 gsl_vector *Rphi,
                 gsl_matrix *A,
                 gsl_matrix *L,
                 bst_leaf *leaf);

double bst_cutplane(gsl_vector *phi_l,
                    gsl_vector *phi_r,
                    gsl_vector *v);

int bst_node_set(bst_leaf *l_leaf,
                 bst_leaf *r_leaf,
                 bst_node *node);

int bst_node_add(int side,
                 bst_node *old_node,
                 bst_node *new_node);

int bst_search(bst_node *root,
               gsl_vector *phi,
               bst_node **end_node,
               bst_leaf **end_leaf);

int bst_height(bst_node *root);

#endif /* __BST_LIB_H__ */

max_stars_repo_name: americocunhajr/CRFlowLib

max_issues_repo_name: americocunhajr/CRFlowLib

max_forks_repo_name: americocunhajr/CRFlowLib


text:
#ifndef _RMSD_H_
#define _RMSD_H_

#include <algorithm>
#include <map>
#include <cmath>
#include <iostream>
#include <vector>
#include <gsl/gsl_matrix.h>
#include "pdb.h"
#include "atom.h"
#include "align.h"
#include "contact.h"

// Indicates that, when calculating the ligand rmsd between a receptor/ligand pair
// only atoms (in the ligand) closer than 10 angstroms to any receptor atom should
// be considered
#define LIGAND_RMSD_THRESHOLD 10

/*
 * This version of the RMSD calculation method assumes that both vectors
 * are the same size and that each position in each array corresponds to
 * the actual alignment of atoms
 */
double calculate_allatom_rmsd(vector<atom>& S, vector<atom>& T);
double calculate_lrms(vector<atom>&, vector<atom>&, map<int, int>&);
double calculate_irmsd(vector<atom>&, vector<atom>&, vector<atom>&,
    vector<atom>&, map<int, int>&, map<int, int>&, int, int);
void get_fnat(vector<string>& target, vector<string>& predicted, double& fnat, double& fnon_nat);
// overloaded version that not only calculates fnat fot a target/prediction, but for each of the chains too
// according to the per-chain information given to it
void get_fnat(vector<string>& target, vector<string>& predicted, vector< vector<string> >& target_per_chain, vector< vector<string> >& predicted_per_chain,
                         double& fnat, double& fnon_nat, double* fnat_per_chain, double* fnon_nat_per_chain);
// This function takes a multiple chain target and its prediction and calculates the global c-alpha rmsd
// between them. It directly returns the total rmsd but it also returns the specific rmsd for each chain
// in the output parameter named rmsd_by_chain. The memory for this parameter should have been allocated previously
double calculate_multiple_rmsd(vector< vector<atom> >& target, vector< vector<atom> >& predicted, double* rmsd_by_chain, bool interface_only);
#endif

max_stars_repo_name: tecdatalab/legacy

max_issues_repo_name: tecdatalab/legacy

max_forks_repo_name: tecdatalab/legacy


text:
/* specfunc/transport.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

/* Author:  G. Jungman */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_transport.h>

#include "error.h"
#include "check.h"

#include "chebyshev.h"
#include "cheb_eval.c"

static double transport2_data[18] = {
   1.671760446434538503,
  -0.147735359946794490,
   0.148213819946936338e-01,
  -0.14195330326305613e-02,
   0.1306541324415708e-03,
  -0.117155795867579e-04,
   0.10333498445756e-05,
  -0.901911304223e-07,
   0.78177169833e-08,
  -0.6744565684e-09,
   0.579946394e-10,
  -0.49747619e-11,
   0.425961e-12,
  -0.36422e-13,
   0.3111e-14,
  -0.265e-15,
   0.23e-16,
  -0.19e-17
};
static cheb_series transport2_cs = {
  transport2_data,
  17,
  -1, 1,
  9
};

static double transport3_data[18] = {
   0.762012543243872007,
  -0.105674387705058533,
   0.119778084819657810e-01,
  -0.12144015203698307e-02,
   0.1155099769392855e-03,
  -0.105815992124423e-04,
   0.9474663385302e-06,
  -0.836221212858e-07,
   0.73109099278e-08,
  -0.6350594779e-09,
   0.549118282e-10,
  -0.47321395e-11,
   0.4067695e-12,
  -0.348971e-13,
   0.29892e-14,
  -0.256e-15,
   0.219e-16,
  -0.19e-17
};
static cheb_series transport3_cs = {
  transport3_data,
  17,
  -1, 1,
  9
};

static double transport4_data[18] = {
  0.4807570994615110579,
 -0.8175378810321083956e-01,
  0.1002700665975162973e-01,
 -0.10599339359820151e-02,
  0.1034506245030405e-03,
 -0.96442705485899e-05,
  0.8745544408515e-06,
 -0.779321207981e-07,
  0.68649886141e-08,
 -0.5999571076e-09,
  0.521366241e-10,
 -0.45118382e-11,
  0.3892159e-12,
 -0.334936e-13,
  0.28767e-14,
 -0.2467e-15,
  0.211e-16,
 -0.18e-17
};
static cheb_series transport4_cs = {
  transport4_data,
  17,
  -1, 1,
  9
};

static double transport5_data[18] = {
   0.347777777133910789,
  -0.66456988976050428e-01,
   0.8611072656883309e-02,
  -0.9396682223755538e-03,
   0.936324806081513e-04,
  -0.88571319340833e-05,
   0.811914989145e-06,
  -0.72957654233e-07,
   0.646971455e-08,
  -0.568490283e-09,
   0.49625598e-10,
  -0.4310940e-11,
   0.373100e-12,
  -0.32198e-13,
   0.2772e-14,
  -0.238e-15,
   0.21e-16,
  -0.18e-17
};
static cheb_series transport5_cs = {
  transport5_data,
  17,
  -1, 1,
  9
};

static
double
transport_sumexp(const int numexp, const int order, const double t, double x)
{
  double rk = (double)numexp;
  double sumexp = 0.0;
  int k;
  for(k=1; k<=numexp; k++) {
    double sum2 = 1.0;
    double xk  = 1.0/(rk*x);
    double xk1 = 1.0;
    int j;
    for(j=1; j<=order; j++) {
      sum2 = sum2*xk1*xk + 1.0;
      xk1 += 1.0;
    }
    sumexp *= t;
    sumexp += sum2;
    rk -= 1.0;
  }
  return sumexp;
}

/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/

int
gsl_sf_transport_2_e(const double x, gsl_sf_result * result)
{
  const double val_infinity = 3.289868133696452873;

  /* CHECK_POINTER(result) */

  if(x < 0.0) {
    DOMAIN_ERROR(result);
  }
  else if(x < 3.0*GSL_SQRT_DBL_EPSILON) {
    result->val = x;
    result->err = GSL_DBL_EPSILON*fabs(x) + x*x/2.0;
    return GSL_SUCCESS;
  }
  else if(x <= 4.0) {
    double t = (x*x/8.0 - 0.5) - 0.5;
    gsl_sf_result result_c;
    cheb_eval_e(&transport2_cs, t, &result_c);
    result->val  = x * result_c.val;
    result->err  = x * result_c.err;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(x < -GSL_LOG_DBL_EPSILON) {
    const int    numexp = (int)((-GSL_LOG_DBL_EPSILON)/x) + 1;
    const double sumexp = transport_sumexp(numexp, 2, exp(-x), x);
    const double t = 2.0 * log(x) - x + log(sumexp);
    if(t < GSL_LOG_DBL_EPSILON) {
      result->val  = val_infinity;
      result->err  = 2.0 * GSL_DBL_EPSILON * val_infinity;
    }
    else {
      const double et = exp(t);
      result->val = val_infinity - et;
      result->err = 2.0 * GSL_DBL_EPSILON * (val_infinity + fabs(t) * et);
    }
    return GSL_SUCCESS;
  }
  else if(x < 2.0/GSL_DBL_EPSILON) {
    const int    numexp = 1;
    const double sumexp = transport_sumexp(numexp, 2, 1.0, x);
    const double t = 2.0 * log(x) - x + log(sumexp);
    if(t < GSL_LOG_DBL_EPSILON) {
      result->val = val_infinity;
      result->err = 2.0 * GSL_DBL_EPSILON * val_infinity;
    }
    else {
      const double et = exp(t);
      result->val = val_infinity - et;
      result->err = 2.0 * GSL_DBL_EPSILON * (val_infinity + (fabs(t)+1.0) * et);
    }
    return GSL_SUCCESS;
  }
  else {
    const double t = 2.0 * log(x) - x;
    if(t < GSL_LOG_DBL_EPSILON) {
      result->val = val_infinity;
      result->err = 2.0 * GSL_DBL_EPSILON * val_infinity;
    }
    else {
      const double et = exp(t);
      result->val = val_infinity - et;
      result->err = 2.0 * GSL_DBL_EPSILON * (val_infinity + (fabs(t)+1.0) * et);
    }
    return GSL_SUCCESS;
  }
}

int
gsl_sf_transport_3_e(const double x, gsl_sf_result * result)
{
  const double val_infinity = 7.212341418957565712;

  /* CHECK_POINTER(result) */

  if(x < 0.0) {
    DOMAIN_ERROR(result);
  }
  else if(x == 0.0) {
    result->val = 0.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else if(x < 3.0*GSL_SQRT_DBL_EPSILON) {
    result->val = 0.5*x*x;
    result->err = 2.0 * GSL_DBL_EPSILON * result->val;
    CHECK_UNDERFLOW(result);
    return GSL_SUCCESS;
  }
  else if(x <= 4.0) {
    const double x2 = x*x;
    const double t = (x2/8.0 - 0.5) - 0.5;
    gsl_sf_result result_c;
    cheb_eval_e(&transport3_cs, t, &result_c);
    result->val  = x2 * result_c.val;
    result->err  = x2 * result_c.err;
    result->err += GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(x < -GSL_LOG_DBL_EPSILON) {
    const int    numexp = (int)((-GSL_LOG_DBL_EPSILON)/x) + 1;
    const double sumexp = transport_sumexp(numexp, 3, exp(-x), x);
    const double t = 3.0 * log(x) - x + log(sumexp);
    if(t < GSL_LOG_DBL_EPSILON) {
      result->val = val_infinity;
      result->err = 2.0 * GSL_DBL_EPSILON * val_infinity;
    }
    else {
      const double et = exp(t);
      result->val = val_infinity - et;
      result->err = 2.0 * GSL_DBL_EPSILON * (val_infinity + fabs(t) * et);
    }
    return GSL_SUCCESS;
  }
  else if(x < 3.0/GSL_DBL_EPSILON) {
    const int    numexp = 1;
    const double sumexp = transport_sumexp(numexp, 3, 1.0, x);
    const double t = 3.0 * log(x) - x + log(sumexp);
    if(t < GSL_LOG_DBL_EPSILON) {
      result->val = val_infinity;
      result->err = 2.0 * GSL_DBL_EPSILON * val_infinity;
    }
    else {
      const double et = exp(t);
      result->val = val_infinity - et;
      result->err = 2.0 * GSL_DBL_EPSILON * (val_infinity + (fabs(t)+1.0) * et);
    }
    return GSL_SUCCESS;
  }
  else {
    const double t = 3.0 * log(x) - x;
    if(t < GSL_LOG_DBL_EPSILON) {
      result->val = val_infinity;
      result->err = 2.0 * GSL_DBL_EPSILON * val_infinity;
    }
    else {
      const double et = exp(t);
      result->val = val_infinity - et;
      result->err = 2.0 * GSL_DBL_EPSILON * (val_infinity + (fabs(t)+1.0) * et);
    }
    return GSL_SUCCESS;
  }
}

int
gsl_sf_transport_4_e(const double x, gsl_sf_result * result)
{
  const double val_infinity = 25.97575760906731660;

  /* CHECK_POINTER(result) */

  if(x < 0.0) {
    DOMAIN_ERROR(result);
  }
  else if(x == 0.0) {
    result->val = 0.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else if(x < 3.0*GSL_SQRT_DBL_EPSILON) {
    result->val = x*x*x/3.0;
    result->err = 3.0 * GSL_DBL_EPSILON * result->val;
    CHECK_UNDERFLOW(result);
    return GSL_SUCCESS;
  }
  else if(x <= 4.0) {
    const double x2 = x*x;
    const double t = (x2/8.0 - 0.5) - 0.5;
    gsl_sf_result result_c;
    cheb_eval_e(&transport4_cs, t, &result_c);
    result->val  = x2*x * result_c.val;
    result->err  = x2*x * result_c.err;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(x < -GSL_LOG_DBL_EPSILON) {
    const int    numexp = (int)((-GSL_LOG_DBL_EPSILON)/x) + 1;
    const double sumexp = transport_sumexp(numexp, 4, exp(-x), x);
    const double t = 4.0 * log(x) - x + log(sumexp);
    if(t < GSL_LOG_DBL_EPSILON) {
      result->val = val_infinity;
      result->err = 2.0 * GSL_DBL_EPSILON * val_infinity;
    }
    else {
      const double et = exp(t);
      result->val = val_infinity - et;
      result->err = 2.0 * GSL_DBL_EPSILON * (val_infinity + (fabs(t)+1.0) * et);
    }
    return GSL_SUCCESS;
  }
  else if(x < 3.0/GSL_DBL_EPSILON) {
    const int    numexp = 1;
    const double sumexp = transport_sumexp(numexp, 4, 1.0, x);
    const double t = 4.0 * log(x) - x + log(sumexp);
    if(t < GSL_LOG_DBL_EPSILON) {
      result->val = val_infinity;
      result->err = 2.0 * GSL_DBL_EPSILON * val_infinity;
    }
    else {
      const double et = exp(t);
      result->val = val_infinity - et;
      result->err = 2.0 * GSL_DBL_EPSILON * (val_infinity + (fabs(t)+1.0) * et);
    }
    return GSL_SUCCESS;
  }
  else {
    const double t = 4.0 * log(x) - x;
    if(t < GSL_LOG_DBL_EPSILON) {
      result->val = val_infinity;
      result->err = 2.0 * GSL_DBL_EPSILON * val_infinity;
    }
    else {
      const double et = exp(t);
      result->val = val_infinity - et;
      result->err = 2.0 * GSL_DBL_EPSILON * (val_infinity + (fabs(t)+1.0) * et);
    }
    return GSL_SUCCESS;
  }
}

int
gsl_sf_transport_5_e(const double x, gsl_sf_result * result)
{
  const double val_infinity = 124.4313306172043912;

  /* CHECK_POINTER(result) */

  if(x < 0.0) {
    DOMAIN_ERROR(result);
  }
  else if(x == 0.0) {
    result->val = 0.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else if(x < 3.0*GSL_SQRT_DBL_EPSILON) {
    result->val = x*x*x*x/4.0;
    result->err = 4.0 * GSL_DBL_EPSILON * result->val;
    CHECK_UNDERFLOW(result);
    return GSL_SUCCESS;
  }
  else if(x <= 4.0) {
    const double x2 = x*x;
    const double t = (x2/8.0 - 0.5) - 0.5;
    gsl_sf_result result_c;
    cheb_eval_e(&transport5_cs, t, &result_c);
    result->val  = x2*x2 * result_c.val;
    result->err  = x2*x2 * result_c.err;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(x < -GSL_LOG_DBL_EPSILON) {
    const int    numexp = (int)((-GSL_LOG_DBL_EPSILON)/x) + 1;
    const double sumexp = transport_sumexp(numexp, 5, exp(-x), x);
    const double t = 5.0 * log(x) - x + log(sumexp);
    if(t < GSL_LOG_DBL_EPSILON) {
      result->val = val_infinity;
      result->err = 2.0 * GSL_DBL_EPSILON * val_infinity;
    }
    else {
      const double et = exp(t);
      result->val = val_infinity - et;
      result->err = 2.0 * GSL_DBL_EPSILON * (val_infinity + (fabs(t)+1.0) * et);
    }
    return GSL_SUCCESS;
  }
  else if(x < 3.0/GSL_DBL_EPSILON) {
    const int    numexp = 1;
    const double sumexp = transport_sumexp(numexp, 5, 1.0, x);
    const double t = 5.0 * log(x) - x + log(sumexp);
    if(t < GSL_LOG_DBL_EPSILON) {
      result->val = val_infinity;
      result->err = 2.0 * GSL_DBL_EPSILON * val_infinity;
    }
    else {
      const double et = exp(t);
      result->val = val_infinity - et;
      result->err = 2.0 * GSL_DBL_EPSILON * (val_infinity + (fabs(t)+1.0) * et);
    }
    return GSL_SUCCESS;
  }
  else {
    const double t = 5.0 * log(x) - x;
    if(t < GSL_LOG_DBL_EPSILON) {
      result->val = val_infinity;
      result->err = 2.0 * GSL_DBL_EPSILON * val_infinity;
    }
    else {
      const double et = exp(t);
      result->val = val_infinity - et;
      result->err = 2.0 * GSL_DBL_EPSILON * (val_infinity + (fabs(t)+1.0) * et);
    }
    return GSL_SUCCESS;
  }
}

/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/

#include "eval.h"

double gsl_sf_transport_2(const double x)
{
  EVAL_RESULT(gsl_sf_transport_2_e(x, &result));
}

double gsl_sf_transport_3(const double x)
{
  EVAL_RESULT(gsl_sf_transport_3_e(x, &result));
}

double gsl_sf_transport_4(const double x)
{
  EVAL_RESULT(gsl_sf_transport_4_e(x, &result));
}

double gsl_sf_transport_5(const double x)
{
  EVAL_RESULT(gsl_sf_transport_5_e(x, &result));
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE