/*****************************************************************************
FILE: modtran_utils.c

PURPOSE:  This file contains routines for working with MODTRAN data.
*****************************************************************************/
#include <stdio.h>
#include "espa_geoloc.h"
#include "const.h"
#include "utilities.h"

#include "modtran_utils.h"

/*****************************************************************************
Method:  load_elevations

Description:  Loads the grid elevations into a data structure.

Notes:
    1. The grid elevation file must be present in the current working directory.
    2. The grid elevation entries should be in sync with the grid file.

RETURN: SUCCESS
        FAILURE
*****************************************************************************/
static int load_elevations
(
    MODEL_POINTS *model_points
)
{
    const char *FUNC_NAME = "load_elevations";

    FILE *elevation_fd = NULL;

    int status;
    int index;   /* Index into point structure */

    char elevation_filename[] = "grid_elevations.txt";
    char errmsg[PATH_MAX];

    snprintf(errmsg, sizeof(errmsg), "Failed reading %s", elevation_filename);

    elevation_fd = fopen(elevation_filename, "r");
    if (elevation_fd == NULL)
    {
        RETURN_ERROR(errmsg, FUNC_NAME, FAILURE);
    }

    /* Read the elevations into the 0 elevation positions in the model 
       point structure.  The file and structure should have the same order. */
    for (index = 0; index < model_points->count; index++)
    {
        MODEL_POINT *model_ptr = &model_points->points[index];

        /* Keep looking for a model point that was actually run. */
        if (model_ptr->ran_model == 0)
        {
            continue; 
        }

        status = fscanf(elevation_fd, "%lf %lf\n", 
            &(model_ptr->elevations->elevation),
            &(model_ptr->elevations->elevation_directory));
        if (status != 2)
        {
            RETURN_ERROR(errmsg, FUNC_NAME, FAILURE);
        }
    }

    fclose(elevation_fd);

    return SUCCESS;
}


/*****************************************************************************
Method:  initialize_model_points

Description:  Allocate the memory need to hold the model results and
              initialize known values.

RETURN: SUCCESS
FAILURE
*****************************************************************************/
int initialize_model_points
(
    GRID_POINTS *grid_points,    /* I: The coordinate points */
    MODEL_POINTS *model_points,  /* O: Memory Allocated */
    bool modtran_flag            /* I: 0=RTTOV, 1=MODTRAN */
)
{
    const char FUNC_NAME[] = "initialize_model_points";

    double gndalt[MAX_NUM_ELEVATIONS];

    int index;
    int num_elevations;            /* Number of elevations actually used */
    int elevation_index;           /* Index into elevations */
    int status;                    /* Function return status */
    bool antimeridian_crossing = 0; /* 0=doesn't cross, 1=does cross */

    if (modtran_flag)
    {
        FILE *model_elevation_fd = NULL;

        char model_elevation_filename[] = "modtran_elevations.txt";
        char errmsg[PATH_MAX];

        snprintf(errmsg, sizeof(errmsg), "Failed reading %s", 
            model_elevation_filename);

        model_elevation_fd = fopen(model_elevation_filename, "r");
        if (model_elevation_fd == NULL)
        {
            RETURN_ERROR(errmsg, FUNC_NAME, FAILURE);
        }

        status = fscanf(model_elevation_fd, "%d\n", &num_elevations);
        if (status != 1)
        {
            RETURN_ERROR(errmsg, FUNC_NAME, FAILURE);
        }

        /* Read the elevations into the gndalt structure. */ 
        for (index = 0; index < num_elevations; index++)
        {
            status = fscanf(model_elevation_fd, "%lf\n", &gndalt[index]);
            if (status <= 0)
            {
                RETURN_ERROR(errmsg, FUNC_NAME, FAILURE);
            }
        }

        fclose(model_elevation_fd);
    }
    else
    {
        /* RTTOV is only run at 1 elevation, the elevation of the reanalysis
           pixel */
        num_elevations = 1;
    }

    model_points->count = grid_points->count;

    model_points->points = malloc(model_points->count *
                                    sizeof(MODEL_POINT));
    if (model_points->points == NULL)
    {
        RETURN_ERROR("Failed allocating memory for model points",
                     FUNC_NAME, FAILURE);
    }

    /* Check if the grid points cross the antimeridian */
    if (!modtran_flag)
    {
        antimeridian_crossing = check_grid_antimeridian(grid_points);
    }

    for (index = 0; index < model_points->count; index++)
    {
        /* convenience pointers */
        MODEL_POINT *model_ptr = &model_points->points[index];
        GRID_POINT *grid_ptr = &grid_points->points[index];
 
        model_ptr->count = num_elevations;
        model_ptr->ran_model = grid_ptr->run_model;
        model_ptr->row = grid_ptr->row;
        model_ptr->col = grid_ptr->col;
        model_ptr->reanalysis_row = grid_ptr->reanalysis_row;
        model_ptr->reanalysis_col = grid_ptr->reanalysis_col;
        model_ptr->lon = grid_ptr->lon;
        /* Handle the antimeridian case.  This shift is only needed for RTTOV,
           since the interpolation in RTTOV cases has strict requirements */
        if (!modtran_flag && antimeridian_crossing)
        {
            if (grid_ptr->lon < 0)
            {
                model_ptr->lon += 360.0;
            }
        }
        model_ptr->lat = grid_ptr->lat;
        model_ptr->map_x = grid_ptr->map_x;
        model_ptr->map_y = grid_ptr->map_y;

        model_ptr->elevations = malloc(num_elevations*
                                       sizeof(MODEL_ELEVATION));
        if (model_ptr->elevations == NULL)
        {
            RETURN_ERROR("Failed allocating memory for model point"
                         " elevations", FUNC_NAME, FAILURE);
        }

        /* If MODTRAN is used, multiple elevations are processed.  There is
           no else clause for RTTOV because we don't make a set of elevation
           subdirectories for RTTOV, and the elevation comes from the DEM */
        if (modtran_flag)
        {
            /* Iterate over the elevations and assign the elevation values. */
            for (elevation_index = 0; elevation_index < num_elevations; 
                elevation_index++)
            {
                MODEL_ELEVATION *model_elev =
                                        &model_ptr->elevations[elevation_index];

                model_elev->elevation = gndalt[elevation_index];
                model_elev->elevation_directory = gndalt[elevation_index];
            }
        }
    }

    if (modtran_flag)
    {
        /* Load the first elevation values if needed. */
        if (load_elevations(model_points) != SUCCESS)
        {
            RETURN_ERROR("calling load_elevations", FUNC_NAME, EXIT_FAILURE);
        }
    }

    return SUCCESS;
}


/*****************************************************************************
 * Method:  free_model_points
 *
 * Description:  Free allocated memory for the model points.
 * *****************************************************************************/
void free_model_points
(
    MODEL_POINTS *model_points
)
{
    int index;     /* Index into model points structure */

    for (index = 0; index < model_points->count; index++)
    {
        free(model_points->points[index].elevations);
        model_points->points[index].elevations = NULL;
    }

    free(model_points->points);
    model_points->points = NULL;
}