#ifndef INTERPOLATE_H
#define INTERPOLATE_H

#include <stdbool.h>

#include "grid_points.h"
#include "modtran_utils.h"


/* Defines index locations for the parameters in the at_height array */
typedef enum
{
    AHP_TRANSMISSION,
    AHP_UPWELLED_RADIANCE,
    AHP_DOWNWELLED_RADIANCE,
    AHP_NUM_PARAMETERS
} AT_HEIGHT_PARAMETERS;

/* Defines index locations in the vertices array for the current cell to be
   used for interpolation of the pixel */
typedef enum
{
    LL_POINT,
    UL_POINT,
    UR_POINT,
    LR_POINT,
    NUM_CELL_POINTS
} CELL_POINTS;

void interpolate_parameters
(
    MODEL_POINT *model_points,     /* I: results from model runs */
    const GRID_POINTS *points,           /* I: coordinate points */
    int *cell_vertices,            /* I: current cell vertices */
    double interpolate_height,     /* I: current landsat pixel height */
    double interpolate_easting,    /* I: interpolate to easting */
    double interpolate_northing,   /* I: interpolate to northing */
    double *parameters             /* O: interpolated pixel atmospheric 
                                         parameters */
);

int dec_binsearch
(
    double *d,    /* I: Array to search in */
    double v,     /* I: Element to search for */
    int imin,     /* I: Minimum of range to search */
    int imax      /* I: Maximum of range to search */
);

int inc_binsearch
(
    double *d,    /* I: Array to search in */
    double v,     /* I: Element to search for */
    int imin,     /* I: Minimum of range to search */
    int imax      /* I: Maximum of range to search */
);

int findIx
(
    double qtarg, /* I: Element to search for */
    double *samp, /* I: Array to search in */
    int nsamps    /* I: Size of search array */
);

double dval2
(
    const double *d, /* I: Look up value in this flattened 2-D array */
    int x,           /* I: X axis location */
    int y,           /* I: Y axis location */
    int ysize        /* I: Size of Y axis */
);

void setup_smooth_interpolate
(
    int num_points,                    /* I: Number of reanalysis points */
    int nsampsx,                       /* I: Number of reanalysis points in X
                                             direction */
    int nsampsy,                       /* I: Number of reanalysis points in Y
                                             direction */
    const MODEL_POINTS *model_results, /* I: results from model runs */
    bool antimeridian_crossing,        /* I: 0=doesn't cross, 1=does cross */
    double *upwelled_radiance,         /* O: Upwelled radiance */
    double *downwelled_radiance,       /* O: Downwelled radiance */
    double *transmission,              /* O: Transmission */
    double *xsamp,                     /* O: X linear coordinate vector */
    double *ysamp                      /* O: Y linear coordinate vector */
);

double smooth_interpolate_parameters
(
    double pixel_lat,        /* I: Latitude of the current pixel */
    double pixel_lon,        /* I: Longitude of the current pixel */
    const double *parameter, /* I: Upwelled radiance, downwelled radiance, or
                                   atmospheric transmittance */
    const double *xsamp,     /* I: Reanalysis latitudes along x axis */
    const double *ysamp,     /* I: Reanalysis longitudes along y axis*/
    int ix,                  /* I: X index into reanalysis grid */
    int iy,                  /* I: Y index into reanalysis grid */
    int nsampsy              /* I: Number of reanalysis points along y axis */
);

#endif