/***************************************************************************
 *   Copyright (C) 2005 by Nicolas PASCAL   *
 *   nicolas.pascal@icare.univ-lille1.fr   *
 *                                                                         *
 *   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 2 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.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/
#ifndef PARASOLFILEDATA_H
#define PARASOLFILEDATA_H

#include "filedata.h"
#include "parasolleader.h"
#include "parasoldata.h"
#include "satellitefiledata.h"

#include <cmath>
#include <cfloat>
#include "hdfi.h"

/**
@brief manages the reading of PARASOL files.
@author Nicolas PASCAL
*/
class PARASOLFileData : public SatelliteFileData
{
public:
    static const int nb_parasol_direction;
    // altitude of the satellite in m. It is a virtual mean altitude. The real movement is an ellipse, not a perfect circle
    static const double satellite_altitude;
    // size in m of a NADIR pixel border (6km)
    static const double pix_sz_nadir;
    // maximum size in m of an observed pixel border, ie a pixel observed with a zenithal angle of 60 degrees
    static const double pix_sz_max;
    // maximum angle between viewibng direction and satellite nadir : 60 degrees
    static const double max_sat_nadir_angle;
    // angle in radians at satellite position on an observed nadir pixel at full resolution (= 1/16 degree earth angle)
    static const double pix_obs_nadir_angle;
    // maximum number of acquisition sequences
    static const int nb_acq_seq_max;
    // number of images during an acquisition sequence
    static const int nb_acq_img_max;

protected:
    /**
    * The leader reader
    */
    PARASOLLeader* leader;
    /**
    * The Data reader
    */
    PARASOLData* data;
private:
    template <typename NumType>
    static bool greater_comparator (NumType i,NumType j) { return (i > j); };
    /**
    defines the possible processing lines of PARASOL L2 products
    */
    enum ProcessingLine{
        UNDEFINED,
        BASIC,
        LAND_SURFACE,
        OCEAN_COLOR,
        RADIATION_BUDGET
    };

    /**
    * lat2time_coef is a constant coefficient that makes the links the latitude of a record and the acquisition time, using the equation :
    * time=time(ascending_node)+latitude*lat2time_coef, where lat2time_coef, time(ascending_node) are constant
    * lat2time_coef can be decomposed as : 16/233 * 1/360 * (24*3600)
    */
    static const double lat2time_coef;

    /**
    * the POLDER instrument used : 1 if POLDER1, 2 if POLDER2, -1 if unknown
    */
    int instrument;
    /**
    * The file's level : 1,2 or 3. -1 if unknown
    * @warning level 3 isn't treated in a good way, because the file name convention is different than level 1 and 2.
    */
    int level;
    /**
    * The processing line
    */
    ProcessingLine processing_line;
    /**
    * The product : A, B or C. Space character if unknown
    */
    char product;
    /**
    * The orbit cycle number. -1 if undefined.
    * @warning This attribute is unuseful for level 3 files
    */
    int orbit_cycle_nb;
    /**
    * The orbit number in one cycle. -1 if undefined.
    * @warning This attribute is unuseful for level 3 files
    */
    int orbit_nb;
    /**
    * The reprocessing number : a character from 'A' to 'Z'.  Space character if undefined.
    */
    char reprocessing_nb;
    /**
    * The grid reduction factor. It is defined as :
     * <i>grid_factor = full_resolution_grid_x_max/product_grid_x_max = full_resolution_grid_y_max/product_grid_y_max</i>
    */
    int grid_factor;
    /** size of the as [nrow, ncol]. ncol is the maximum number of columns (at the equator) */
    int sz_grid [2];
    /** size of a cell side, in degrees */
    double sz_cell;
    /** version of the software used to process the file */
    string software_version_number;
    /**
    * the line indexes
    */
    unsigned short *line_data;
    /**
    * the column indexes
    */
    unsigned short *col_data;
    /**
    * geolocation frame
    */
    float lat_min,lat_max,lon_min,lon_max;
    // inner data buffers used to speed up viewing directions computations
    /** viewed pixel altitude data buffer */
    short * v_alt_view;
    /** solar azimuth angles data buffer. 1 per record */
    double * v_saa;
    /** viewing zenith angles data buffer. N_directions  per record */
    double * v_vza;
    /** relative azimuth angles data buffer. N_directions  per record */
    double * v_raa;
    /** number of available directions data buffer */
    unsigned char * v_ni;
    /** acquisition sequence number */
    unsigned char * v_sn;
    /** X, Y, Z coordinates of the satellite in ECR system */
    double * v_x_sat;
    double * v_y_sat;
    double * v_z_sat;

    /**
     * @brief preload the requested data buffers for sppeding up the computations of pixels viewing directions
     */
    void init_viewing_direction();

    void init();
    void parse_filename(const string& short_filename);
    bool check_filename(const string& short_filename) const;
    /**
     * @brief check if the given grid indices are valid
     * @throw g_exception if indices are out of range. All indices start from 1
     * @param irow grid row indice
     * @param icol grid column indice
     */
    void check_grid_coord (const int irow, const int icol) const;
    /**
     * @brief set the PARASOL grid reduction factor depending of the product.
     * The definition of the reduction factor is :
     * <i>factor = full_resolution_grid_x_max/product_grid_x_max = full_resolution_grid_y_max/product_grid_y_max</i>
     */
    void set_grid_factor();
    /**
     * @brief initialize if needed the reading range of the read_data/read_scaled_data/read_count_data methods
     * @warning If start and edges are NULL, the method does the allocation itself. In this case, the caller must free the memory using the delete[] operator. Also, if rank is set to -1, the method sets it to the good value.
     * @param entry_idx index of the entry to be read
     * @param start the start indexes of the data to read.
     * @param edges the number of data to read along each direction.
     * @param rank the size of the start and edges arrays
     * @return true if the method does the allocations
     */
    bool init_read_data_range(const int entry_idx, int * &start, int * &edges, int &rank);
    /**
     * @brief fill the buffer @a data with the values read in the PARASOL file for the given parameter (using its index).
     * The filled buffer has only one dimension. So, for accessing to the value of a given direction (if you have load all directions values) you must use the formula :
     * <b>data[ record_index*nb_directions + direction_index ]</b>, record_index and direction_index starting at 0.
     * @param data the buffer to fill. If NULL, the method does the needed allocations
     * @param var_idx the variable index
     * @param start the start indexes of the data to read.
     * @param edges the number of data to read along each direction.
     * @param rank the size of the start and edges arrays
     * @return the pointer to the data buffer
     */
    void* read_count_data(void* data, const int var_idx = -1, int * start = NULL, int * edges = NULL, int rank = -1);
    /**
     * @brief read the calibration of a parameter, using its index
     * @warning DON'T USE THIS METHOD for accessing the record number, the record length, line nb, col nb, Pixel Altitude or Land/Water/Mixed indicator. They are not supposed to have a number. Read them using their name instead. Actually, this method is not supposed to be used as it. Prefer the get_scaling one that uses the parameter name.
     * @param var_idx the index of the parameter. This index is called the "parametre number" in the Parasol DPC
     * @param slope the scaling slope
     * @param offset the scaling offset
     * @param nb_bytes the number of bytes used by this parameter count
     */
    void get_scaling(const int var_idx, double &slope, double &offset, unsigned short& nb_bytes);

    /**
     * @brief retrieve the entry format for the given variable, using its index
     * @param var_idx the index of the variable
     * @return the entry format, or NULL if no entry with this name can be found
     */
    const EntryFormat *get_entry_format(const int &var_idx);
protected:
    /**
     * @brief access the size of the given parametre data
     * @param var_idx the index of the variable (also known as parametre)
     * @return a vector that contains the dimensions of the data in each direction (convention [Z,Y,X]). Use return_vector.size() to know the number of dimensions
     */
    vector<int> get_data_dimension(const int var_idx);
public:
    /**
     * Constructor
     * @param name the name of the file to load. Can be either the leader filename, or the data file one.
     * @warning it assumes the data and leader files are in the same directory.
     * @param mode opening mode. At this time, only "r" is allowed
     */
    PARASOLFileData(const string& name, const string& mode="r");
    /**
     * Destructor
     */
    ~PARASOLFileData();
    /**
     * @brief access to the radix of the file.
     * The 'Radix' is the filename with leader or data indicator removed (ie the ending 'L' or 'D' character).
     * @return the radix of the file
     */
    const string get_radix();
    /**
     * @brief access to the instrument code. -1 if unknown
     * @return the instrument
     */
    const int get_instrument() const {
        return instrument;
    }
    /**
     * @brief access to the product level. -1 if unknown
     * @return the product level
     */
    const int get_level() const {
        return level;
    }
    /**
     * @brief access to the processing line.
     * @return the processing line
     */
    const ProcessingLine get_processing_line() const {
        return processing_line;
    }
    /**
     * @brief access to the processing line, using a string representation, like "ocean color", "radiation budget"...
     * @return a string representing the processing line
     */
    const string get_processing_line_string() const;
    /**
     * @brief access to the processing line, using the PARASOL character representation, like 'O' for "ocean color", 'R' for "radiation budget"...
     * return a space character if undefined.
     * @return a character representing the processing line
     */
    const char get_processing_line_char() const;
    /**
     * @brief access to the product code (A,B,C or space character if unknown)
     * @return the product code
     */
    const char get_product() const {
        return product;
    }
    /**
     * @brief access to the orbit cycle number. -1 if unknown
     * @return the orbit cycle number
     */
    const int get_orbit_cycle_nb() const {
        return orbit_cycle_nb;
    }
    /**
     * @brief access to the orbit number of the cycle. -1 if unknown
     * @return the orbit number
     */
    const int get_orbit_nb() const {
        return orbit_nb;
    }
    /**
     * @brief access to the reprocessing number. -1 if unknown
     * @return the reprocessing number
     */
    const char get_reprocessing_nb() const {
        return reprocessing_nb;
    }
    /**
     * @brief get the grid indices corresponding to a point given by its (lat,lon) coordinates
     * @param lat_lon [IN] the geolocation coordinates
     * @param grid_idx [OUT] the indices in the POLDER grid. Format [irow,icol]. Indices start from 1
     */
    void geolocation_to_grid(const float * lat_lon, int * grid_idx);
    /**
     * @brief get the grid indices corresponding to a point given by its (lat,lon) coordinates
     * @throw g_exception if (lat, lon) not valid
     * Indices start from 1
     * @param lat [IN] latitude
     * @param lon [IN] longitude
     * @param irow [OUT] grid row indice. Starts from 1
     * @param icol [OUT] grid column indice. Starts from 1
     */
    void geolocation_to_grid(const float lat, const float lon, int & irow, int & icol);
    /**
     * @brief get the lat,lon extrema for a grid point
     * if the indexes are out of the grid or not valid, lat_minmax and lon_minmax are set to -1.
     * @param igrid the indexes in the POLDER grid (line,col). Starts at 1
     * @param lat_minmax 2 elements array containing the latitude extrema. Mustn't be NULL
     * @param lon_minmax 2 elements array containing the longitude extrema. Can be NULL. In this case, only lat_minmax will be computed
     */
    void grid_to_geolocation (const int igrid[2], float lat_minmax[2], float lon_minmax[2]);
    /**
     * @brief get the geolocation for a grid point center
     * @param igrid indices in the POLDER grid (line,col). Starts at 1
     * @param pos (lat, lon) at center of grid cell
    */
    void grid_to_geolocation (const int igrid[2], float pos[2]);
    /**
     * @brief get the geolocation for a grid point center
     * @param igrid indices in the POLDER grid (line,col). Starts at 1
     * @param lat latitude at center of grid cell
     * @param lon longitude at center of grid cell
     */
    void grid_to_geolocation (const int igrid[2], float & lat, float & lon);
    /**
     * @brief get the geolocation for a grid point center
     * All indices start from 1
     * @param irow grid row indice
     * @param icol grid column indice
     * @param lat latitude at center of grid cell
     * @param lon longitude at center of grid cell
     */
    void grid_to_geolocation (const int irow, const int icol, float & lat, float & lon);
    /**
     * @brief access to the number of lines of the PARASOL grid for this product
     * @return the number of lines
     */
    const int get_nrow() const;
    /**
     * @brief access to the number of valid columns of the PARASOL grid for this product
     * @param irow the indice of the row to consider. If not given, return the maximum number of colums. Start at 1
     * @return the number of columns
     */
    const int get_ncol(const int irow = -1) const;
    /**
     * @brief access to the number of valid columns of the PARASOL grid for this product at a given latitude
     * @param lat a latitude
     * @return the number of columns
     */
    const int get_ncol(const double lat) const;
    /**
     * @brief test if the given grid indices represent an earth pixel (is it inside the sinusoidal grid)
     * All indices start from 1
     * @param irow grid row indice
     * @param icol grid column indice
     */
    bool is_earth_grid (const int irow, const int icol) const;
    /**
     * @brief return the latitude of the middle of the given row in the PARASOL grid
     * @param irow the indice of the row to consider. If not given, return the maximum number of colums. Start at 1
     * @return the latitude of the middle of the row
     */
    const double get_grid_lat (const int irow) const;
    /**
     * @brief retrieve the coordinates of the record at index @a irec
     * The output coordinates are given for the center of the pixel
     * @param irec [IN] indice of the record. Be careful starts at 0 here but at 2 in the PARASOL data
     * @param lat [OUT] the latitude of the record
     * @param lon [OUT] the longitude of the record
     */
    void get_record_coord(const int irec, float &lat, float &lon);
    /**
     * @brief maps a record indice to its grid coordinates
     * @param irec [IN] indice of the record. Be careful starts at 0 here but at 2 in the PARASOL data
     * @param igrid [OUT] 2D grid coordinates [irow, icol]
     */
    void idata2igrid (const int irec, int igrid[2]);
    /**
     * @brief access the size of the given parametre data
     * @param var_name the name of the variable (also known as parametre)
     * @return a vector that contains the dimensions of the data in each direction (convention [Z,Y,X]). Use return_vector.size() to know the number of dimensions
     */
    vector<int> get_data_dimension(const string &var_name);
    /**
     * @brief fill the buffer @a data with the values read in the PARASOL file for the given parameter (using its index).
     * The filled buffer has only one dimension. So, for accessing to the value of a given direction (if you have load all directions values) you must use the formula :
     * <b>data[ record_index*nb_directions + direction_index ]</b>, record_index and direction_index starting at 0.


     *  The returned buffer is a 1D linearized array, even if the read data have 2D. To know the real dimensions, you can use the "get_data_dimension" method that will return a vector containing the real dimension of the data, and access to one value can be done as y_index*X_DIM+x_index. @n Let's see a little example : reading the dqx of a L1_B file @n
     * @code
     * string param="dqx";
     * unsigned short* bit_data=NULL; // The dqx bitfield is coded on 2 bytes. So on, use an unsigned short array to read it. Set it to NULL to delegate the buffer allocation to the read_data method
     * // read the dimensions of the variable data : for this parametre, has 2 ones, [nb_data_record,NB_PARASOL_DIRECTIONS]
     * vector<int> dimensions=pfd.get_data_dimension(param);
     * int nb_data_record=dimensions[0];
     * int nb_viewing_directions=dimensions[1];
     * // read the data. The dqx is a bitfield, so use the read_count_data method : scaling is not needed (and will produce invalid data)
     * bit_data=static_cast<unsigned short*>(pfd.read_count_data(bit_data,param.c_str()));
     * // print out the values
     * if (bit_data!=NULL) {
     *   if (show_val_display) {
     *       cout<<"Record\tbits"<<endl;
     *       for (int i = 0 ; i<nb_data;++i) {
     *         cout<<i;
     *         for (int j=0;j < nb_viewing_directions;++j)
     *           cout<<"\t"<<MyTools::to_bit(&bit_data[i*nb_viewing_directions+j],2); // Use the MyTools::to_bit method to print out bitfield data
     *         cout<<endl;
     *       }
     *   }
     * }
     * // free the buffer once used
     * delete[] bit_data;
     * bit_data=NULL;
     * @endcode
     * @param data the buffer to fill. If NULL, the method does the needed allocations
     * @param var_name the variable name
     * @param start the start indexes of the data to read.
     * @param edges the number of data to read along each direction.
     * @param rank the size of the start and edges arrays
     * @return the pointer to the data buffer
     */
    void* read_count_data(void * data, const char* var_name = "alt", int * start = NULL, int * edges = NULL, int rank = -1);
    /**
     * @brief read the data and apply the scaling to it.
     * @warning all scaled data are returned in a double format.
     * @warning at this time, @a start, @a stride and @a edges are ignored. You can only read the whole data of a variable
     * The filled buffer has only one dimension. So, for accessing to the value of a given direction, you must use the formula :
     * <b>data[ record_index*nb_directions + direction_index ]</b>, record_index and direction_index starting at 0. In this case, for one record, there always will be 16 directionnal values (16 is the maximal number of available PARASOL directions), and the unavailable directions values will be set to -1. You can also use the "get_data_dimension" method that will return a vector containing the real dimension of the data, and access to one value can be done as y_index*X_DIM+x_index.@n
     * Let's see a little example : read the sequence number of a L1_B file, and print them out @n
     * @code
     * // Initialize the PARASOL file instance
     * PARASOLFileData pfd(filename);
     * string param="seq_nb";
     * cout<<"---------- Extract Parameter "<<param<<" ----------"<<endl;
     * // read the dimensions of the variable data : for this parametre, has 2 ones [nb_data_record],NB_PARASOL_DIRECTIONS]
     * vector<int> dimensions=pfd.get_data_dimension(param);
     * int nb_data_record=dimensions[0];
     * int nb_viewing_directions=dimensions[1];
     * cout<<"Dimensions : ["<<dimensions[0]<<","<<dimensions[1]<<"]"<<endl;
     * cout<<"---------------------------------------------------"<<endl;
     * // read the data
     * signed short* data=NULL; // Set the data buffer to NULL to delegate its allocation to the read_data method
     * data=static_cast<signed short*>(pfd.read_scaled_data(data,param.c_str()));
     * // print out the values
     * if (data!=NULL) {
     *     if (show_val_display) {
     *         for (int j = 0 ; j<nb_data_record;++j) {
     *             cout<<"Record #"<<j<<endl;
     *             for (int i = 0 ; i<nb_viewing_directions;++i) {
     *                 cout<<"\t"<<data[j*nb_viewing_directions+i];
     *             }
     *             cout<<endl;
     *         }
     *     }
     * }
     * // free the buffer once used
     * delete[] data;
     * data=NULL;
     * @endcode
     * @param data a buffer of double values that will contain the data. If NULL, the method does the needed allocations.
     * @param var_name the variable name
     * @param start the start index. If NULL will be 0
     * @param edges the number of values to read. If NULL will be the whole data
     * @return the pointer to data
     */
    double* read_scaled_data(double * data, const char * var_name = "alt", int * start = NULL, int * edges = NULL, int rank = -1);
    /**
     * @brief read the data and apply a scaling to them. This method implements the interface of filedata
     * If the read data is a physical value, it will return it as double. If the variable is a bitfield, it will return it as count (unscaled). You can used unsigned char to decode it.
     * @warning at this time, @a stride parameter is ignored.
      *  The returned buffer is a 1D linearized array, even if the read data have 2D. To know the real dimensions, you can use the "get_data_dimension" method that will return a vector containing the real dimension of the data, and access to one value can be done as y_index*X_DIM+x_index.@n
     * Let's see a little example : reading the sequence number of a L1_B file, and print them out @n
     * @code
     * // Initialize the PARASOL file instance
     * PARASOLFileData pfd(filename);
     * string param="seq_nb";
     * cout<<"---------- Extract Parameter "<<param<<" ----------"<<endl;
     * // read the dimensions of the variable data : for this parametre, has 2 ones [nb_data_record],NB_PARASOL_DIRECTIONS]
     * vector<int> dimensions=pfd.get_data_dimension(param);
     * int nb_data_record=dimensions[0];
     * int nb_viewing_directions=dimensions[1];
     * cout<<"Dimensions : ["<<dimensions[0]<<","<<dimensions[1]<<"]"<<endl;
     * cout<<"---------------------------------------------------"<<endl;
     * // read the data
     * signed short* data=NULL; // Set the data buffer to NULL to delegate its allocation to the read_data method
     * data=static_cast<signed short*>(pfd.read_data(data,param.c_str()));
     * // print out the values
     * if (data!=NULL) {
     *     if (show_val_display) {
     *         for (int j = 0 ; j<nb_data_record;++j) {
     *             cout<<"Record #"<<j<<endl;
     *             for (int i = 0 ; i<nb_viewing_directions;++i) {
     *                 cout<<"\t"<<data[j*nb_viewing_directions+i];
     *             }
     *             cout<<endl;
     *         }
     *     }
     * }
     * // free the buffer once used
     * delete[] data;
     * data=NULL;
     * @endcode
     *
     * Another example : reading the dqx of the first viewing direction of a L1_B file, and print them out @n
     * @code
     * // Initialize the PARASOL file instance
     * PARASOLFileData pfd(filename);
     * string param="dqx";
     * cout<<"---------- Extract Parameter "<<param<<" ----------"<<endl;
     * // read the dimensions of the variable data : for this parametre, has 2 ones [nb_data_record],NB_PARASOL_DIRECTIONS]
     * vector<int> dimensions=pfd.get_data_dimension(param);
     * int nb_data_record=dimensions[0];
     * int nb_viewing_directions=dimensions[1];
     * cout<<"Dimensions : ["<<dimensions[0]<<","<<dimensions[1]<<"]"<<endl;
     * cout<<"---------------------------------------------------"<<endl;
     * // read the data -> set the range of the data to read
     * int start[]={0,4}; // means : start to read the data from the first record (0), and from the 5th viewing direction (4)
     * int edges[]={-1,1}; // means : read all the records (-1 is a kind of metaselector), and read only one viewing direction value (1)
     * signed short* data=NULL; // Set the data buffer to NULL to delegate its allocation to the read_data method
     * data=static_cast<signed short*>(pfd.read_data(data,param.c_str()));
     * // print out the values
     * if (data!=NULL) {
     *     if (show_val_display) {
     *         cout<<"Record\tbits"<<endl;
     *         for (int i = 0 ; i<nb_data_record];++i)
                 cout<<i<<"\t"<<MyTools::to_bit(&ushort_count_data[i],2)<<endl;
     *     }
     * }
     * // free the buffer once used
     * delete[] data;
     * data=NULL;
     * @endcode
     * @param data a buffer that will contain the data.
     * @param var_name the variable name
     * @param start the start indexes of the data to read. Indexes of -1 are considered as metacharacters and will be replaced by 0. Actually, if you read a directionnal parametre, start is [start_record_index, start_direction]
     * @param stride step between 2 read values (ignored at this time : 1 used by default)0.
     * @param edges the number of data to read along each direction. Indexes of -1 are considered as metacharacters and will be replaced by the maximum number of values. Actually, if you read a directionnal parametre, edges is [nb_record_to_read, nb_direction_to_read]
     * @param rank the size of the start and edges arrays : 1 for undirectionnal parametres, and 2 for directionnal ones. This parametre is not mandatory : if not specified, the method computes it.
     * @return the pointer to the data buffer
     * @return the pointer to data
     */
    void * read_data(void * data, const char * var_name, int * start = NULL, int * stride = NULL, int * edges = NULL, int rank = -1);
    /**
     * @brief access to the number of the data record contained in the file
     * @return the number of data record
     */
    const int get_nb_data();
    /**
     * @brief load the sequences of lines and columns indexes
     */
    virtual void load_geolocation_data();
    /**
     * @brief free the sequences of lines and columns indexes
     */
    virtual void free_geolocation_data();
    /**
     * @brief load the list data pixels
     */
    void load_v_pixel() ;
    /**
     * @brief check if the lines and columns indexes have already been loaded
     * @return true if lines and columns indexes are loaded
     */
    virtual const bool is_geolocation_data_loaded() const;
    /**
     * check if this file has data coincident with (lat,lon,time)
     * @param lat the latitude of the event
     * @param lon the longitude of the event
     * @param time the time of the event. If -1., time test will be ignored
     * @return true if this is the case
     */
    const bool contain_data(const float &lat,const float &lon, const double &time);
     /**
     * @brief check if this file has data coincident with (lat,lon,time)
     * @param lat the latitude of the event
     * @param lon the longitude of the event
     * @param time the time of the event. If -1., time test will be ignored
     * @param tolerance colocation tolerance.
     * @return true if this is the case
     */
    const bool contain_data(const float &lat,const float &lon, const double &time, const double &tolerance);
    /**
     * @brief check if this file has EVENTUALLY data coincident with (lat,lon)
     * Actually, it only tests if (lat,lon) is contained in the data's bounding rectangle.
     * @param lat the latitude of the event
     * @param lon the longitude of the event
     * @return true if this is the case
     */
    const bool contain_location(const float &lat,const float &lon);
    /**
     * @brief check if this file has EVENTUALLY data coincident with (lat,lon)
     * Actually, it only tests if (lat,lon) is contained in the data's bounding rectangle.
     * @param lat the latitude of the event
     * @param lon the longitude of the event
     * @param tolerance colocation tolerance.
     * @return true if this is the case
     */
    const bool contain_location(const float &lat,const float &lon, const double &tolerance);
    /**
     * @brief check if the parasol data contain the given pixel
     * @param pix_idx the pixel indexes {line_index,column_index}
     * @return true if the pixel is contained in the data
     */
    const bool contain_pixel(const int* pix_idx);
    /**
     * @brief read the calibration of a parameter, using its name
     * @param var_name the name of the parameter
     * @param slope the scaling slope
     * @param offset the scaling offset
     * @param nb_bytes the number of bytes used by this parameter count
     */
    void get_scaling(const string &var_name, double &slope, double &offset, unsigned short& nb_bytes);
    /**
     * @brief retrieve the entry format for the given variable
     * @param var_name the name of the variable
     * @return the entry format, or NULL if no entry with this name can be found
     */
    const EntryFormat *get_entry_format(const string &var_name);
    /**
     * @brief opens the data for reading
     */
    void open_data_file(){
        assert(data!=NULL);
        data->open_file();
    };
    void close_data_file(){
        assert(data!=NULL);
        data->close_file();
    };
    const bool is_data_file_opened(){
        return (data!=NULL && data->is_file_loaded());
    };
    unsigned short* get_line_data() const {
        assert(line_data!=NULL); // TODO do load geolocation_data before if NULL
        return line_data;
    }
    unsigned short* get_col_data() const {
        assert(col_data!=NULL); // TODO do load geolocation_data before if NULL
        return col_data;
    }
    /**
     * @brief retrieve the record index (from 0) that contains the (lat,lon) point.
     * If the given point isn't contained in the file, the method will return false, and idx will be set to -1
     * @param lat [IN] the latitude
     * @param lon [IN] the longitude
     * @param irec [OUT] indice of the record. Be careful starts at 0 here but at 2 in the PARASOL data
     * @return true if the point has been found. false else
     */
    const bool get_index(const float& lat, const float& lon, int &irec);
    /**
     * @brief retrieve the record index (from 0) that contains the (lat,lon) point.
     * If the given point isn't contained in the file, the method will return false, and idx will be set to -1
     * @param lat [IN] the latitude
     * @param lon [IN] the longitude
     * @param irec [OUT] indice of the record. Be careful starts at 0 here but at 2 in the PARASOL data
     * @param tolerance colocation tolerance.
     * @return true if the point has been found. false else
     */
    const bool get_index(const float &lat, const float& lon, int &irec, const float colocation_tolerance);
    /**
     * @brief retrieve the record index (from 0) that contains the pixel pix_idx.
     * If the given pixel isn't contained in the file, the method will return false, and irec will be set to -1
     * @param ipix [IN] the indexes (line,column) of the searched parasol pixel
     * @param irec [OUT] indice of the record. Be careful starts at 0 here but at 2 in the PARASOL data
     * @return true if the pixel has been found. false else
     */
    const bool get_index(const int ipix[2], int & irec);
    /**
     * @brief maps grid coordinates to the matching record indice
     * @param igrid [IN] 2D grid coordinates [irow, icol]
     * @param irec [OUT] indice of the record. Be careful starts at 0 here but at 2 in the PARASOL data
     */
    const bool get_index(const vector <int> & ipix, int & irec);
    /**
     * @brief maps grid coordinates to the matching record indice
     * @param igrid [IN] 2D grid coordinates [irow, icol]
     * @param irec [OUT] indice of the record. Be careful starts at 0 here but at 2 in the PARASOL data
     */
    void igrid2idata (const int igrid[2], int & irec);
    /**
     * @brief maps grid coordinates to the matching record indice
     * @param igrid [IN] 2D grid coordinates [irow, icol]
     * @param irec [OUT] indice of the record. Be careful starts at 0 here but at 2 in the PARASOL data
     */
    void igrid2idata (const vector<int> & igrid, int & irec);
    /**
     * @brief build the list of indices of pixels that are in colocation tolerance, sorted by increasing distance to (lat,lon)
     * If (lat,lon) is not found or out of the colocalisation_frame, returns an empty vector
     * @param lat the latitude
     * @param lon the longitude
     * @param colocation_tolerance the acceptable bias (in km or degrees. Supposed to be in, a plane approximation) between [lat,lon] and the nearest data point.
     * @return the list of indices of pixels
     */
     virtual void get_vindex(vector < vector < int > > &v_index, const float &lat, const float& lon,
                                            const float colocation_tolerance) {
        UnimplementedMethod e(__FILE__,__LINE__,__PRETTY_FUNCTION__);
        throw e;
     };

    const float get_nearest_point_distance(const float & lat, const float & lon, const float coloc_tolerance);

    /**
     * @brief check if a file (wether leader or data) is currently opened
     * @return true if the leader and/or data file is opened
     */
    const bool is_file_loaded();
    /**
     * @brief close an opened file (wether leader or data)
     */
    void close_file();
    /**
     * @brief access to the maximal number of directions that can be found in a PARASOL file
     * @return the max number of PARASOL directions
     */
    static const int get_max_direction();
    /**
     * @brief read the pixel (lin,col) of the given record
     * @param irec [IN] indice of the record. Be careful starts at 0 here but at 2 in the PARASOL data
     * @param ipix [OUT] 2D grid coordinates [irow, icol]
     */
    void get_pixel(const int &irec, int ipix[2]);
    /**
     * @brief print out the scaling factors of all parametres
     */
    void print_scaling_factors(){
        leader->print_scaling_factors();
    };
    /**
     * @brief print out the spatio temporal characteristics as set in the leader file
     */
    void print_spatio_temp_char(){
        leader->print_spatio_temp_char();
    };
    /**
     * @brief print out the technological paramaters as set in the leader file
     */
    void print_techno_param(){
        leader->print_techno_param();
    };
    /**
     * @brief access the measure time of the record @a irec
     * @param irec [IN] indice of the record. Be careful starts at 0 here but at 2 in the PARASOL data
     * @return the TAI 93 time
     */
    inline const double get_record_time(const int &irec) {
      assert(leader->is_record_loaded(SPATIO_TEMPORAL_CHARACTERISTICS));
      int pixel[]={-1,-1};
      // record to pixel mapping
      get_pixel(irec,pixel);
      // pixel to geolocation mapping
      float latminmax[]={-1.,-1.};
      float lonminmax[]={-1.,-1.};
      grid_to_geolocation(pixel,latminmax,lonminmax);
      float lat=(latminmax[1]-latminmax[0])/2.;
      // Ascending node time
      string s_asc_node_time(leader->spatio_temp_char->desc_node_date);
      s_asc_node_time=s_asc_node_time.substr(0,s_asc_node_time.size()-2);
      Date d; d.set_date_str(s_asc_node_time,"%Y%m%d%H%M%S");
      return d.get_TAI93_time()+lat*lat2time_coef;
    };

    /**
     * @brief compute the view zenith angle for any filter, using the parametres stored in the PARASOL L1B files
     * This method is based on the formula given in the Appendix C of the PARASOL L1B Product Guide
     * For PARASOL, the @a xj values are : -6:490P, -4:443NP, -3:1020NP, -2:565NP, 0:670P, 2:763NP, 3:765NP, 4:910NP, 6:865P
     * @param view_zenith_angle_filter8 the view zenith angle for filter #8
     * @param relative_azimuth_angle_filter8 the relative azimuth angle for filter #8
     * @param delta_thetav_cosphi the Theta_v*Cos_Phi relative variation of viewing geometry between the filters in degrees
     * @param delta_thetav_sinphi the Theta_v*Sin_Phi relative variation of viewing geometry between the filters in degrees
     * @param xj index of desired filter. See Appendix C of the PARASOL L1B Product Guide for details
     * @return the view zenith angle for the given filter in degrees
     */
    static inline const double get_filter_view_zenith_angle(const double &view_zenith_angle_filter8,
                                                             const double &relative_azimuth_angle_filter8,
                                                             const double &delta_thetav_cosphi,
                                                             const double &delta_thetav_sinphi,
                                                             const int    xj) {
        double a1=view_zenith_angle_filter8*cos(relative_azimuth_angle_filter8) +
                      (double)(xj)*delta_thetav_cosphi;
        double a2=view_zenith_angle_filter8*sin(relative_azimuth_angle_filter8) +
                      (double)(xj)*delta_thetav_sinphi;
        return sqrt(pow (a1,2) + pow(a2,2));
    }
    /**
     * @brief compute the relative azimuth angle for any filter, using the parametres stored in the PARASOL L1B files
     * This method is based on the formula given in the Appendix C of the PARASOL L1B Product Guide
     * For PARASOL, the @a xj values are : -6:490P, -4:443NP, -3:1020NP, -2:565NP, 0:670P, 2:763NP, 3:765NP, 4:910NP, 6:865P
     * @param view_zenith_angle_filter8 the view zenith angle for filter #8
     * @param relative_azimuth_angle_filter8 the relative azimuth angle for filter #8
     * @param delta_thetav_cosphi the Theta_v*Cos_Phi relative variation of viewing geometry between the filters in degrees
     * @param delta_thetav_sinphi the Theta_v*Sin_Phi relative variation of viewing geometry between the filters in degrees
     * @param xj index of desired filter. See Appendix C of the PARASOL L1B Product Guide for details
     * @return the relative azimuth angle for the given filter in degrees
     */
    static inline const double get_filter_relative_azimuth_angle(const double &view_zenith_angle_filter8,
                                                                  const double &relative_azimuth_angle_filter8,
                                                                  const double &delta_thetav_cosphi,
                                                                  const double &delta_thetav_sinphi,
                                                                  const int    xj) {
        double a1=view_zenith_angle_filter8*sin(relative_azimuth_angle_filter8) +
                      (double)(xj)*delta_thetav_sinphi;
        double a2=view_zenith_angle_filter8*cos(relative_azimuth_angle_filter8) +
                      (double)(xj)*delta_thetav_cosphi;
        double offset=0.;
        if (a2<0)
            offset+=180.;
        return offset + atan (a1/a2) ;
    }
    /**
     * @brief retrieve the coordinates of a pixel using its index
     * @param ipix [IN] index of the pixel
     * @param lat [OUT] latitude of the pixel
     * @param lon [OUT] longitude of the pixel
     * @param time [OUT] timestamp of the pixel
     */
    virtual void get_pixel_coord (const vector < int > & ipix, float &lat, float &lon, double &time) {
        UnimplementedMethod e(__FILE__,__LINE__,__PRETTY_FUNCTION__);
        throw e;
    };
    /**
     * @brief accessor to the version of the software used to process the file
     * @return the software version number as a 6 characters string
     */
    const string get_software_version_number() const {
        return software_version_number;
    }
    /**
     * @brief build the earth viewing directions, as a segment from satellite position to viewed pixel center. All positions are stored in Earth Center Rotating carthesian coordinates.
     * The method used is different for L1 and L2 products (L3 nbot supported at this time), because the available parametres are different. In L1, the satellite position can directly be read in the technological parametres. In L2, it must deduced from the viewing directions
     * For directionnal products (L1), a pixel is identified by the 3 indices : [irow, icolum, idirection]. For non directionnal ones, only [irow, icolum] is used
     */
//     virtual void load_viewing_directions ();
    /**
     * @brief constucts the viewing directions observations for the given grid pixel
     * @param ipix 2D grid indices.
     * @param v_obs vector of observation(s). Can contains more than one for directional products
     */
    virtual void get_viewing_directions (const vector <int> & ipix, vector <Observation> & v_obs);
    /**
     * @brief test if the data requested for computing the viewing directions has been loaded
     */
    virtual bool is_viewing_directions_data_loaded();
    /**
     * @brief load the data requested for computing the viewing directions
     */
    virtual void load_viewing_directions_data ();
    /**
     * @brief free the data requested for computing the viewing directions
     */
    virtual void free_viewing_directions_data ();
    /**
     * @brief build the table that maps the gridded pixels indices to their data representation
     *  Most of time, it will be exactly the same, but for instance, in the PARASOL data, it will maps the gridded pixels indices to their matching record number
     */
    virtual void load_pix2data_map ();
    /**
     * @brief computes the list of grid coordinates of the pixels that are neighbours of [irow, icol], in the area +/- [d_irow, d_icol]
     * All indices in this method start at 1
     * @param irow [IN] grid row indice
     * @param icol [IN] grid column indice
     * @param d_irow [IN] number of neighbours pixels to consider along rows
     * @param d_icol [IN] number of neighbours pixels to consider along columns
     * @param v_irow [OUT] indices of the neighbour pixels along rows (sorted increasing)
     * @param v_icol [OUT] indices of the neighbour pixels along columns (sorted increasing)
     */
//     void get_neighbours_coord (const int irow, const int icol, const int d_irow, const int d_icol, vector<int> & v_irow, vector<int> & v_icol);
    /**
     * @brief computes the list of grid coordinates of the pixels that are neighbours of [irow, icol], in the area [d_irow, d_icol]
     * All indices in this method start at 1.
     * @warning as the lines in the [irow + dirow_min -> irow + d_irow_max] can have different number of valid columns in an sinusoidal grid, @a v_icol is returned with the maximum number of columns. As a consequence, some returned [irow, icol] pixels can be invalid.
     * For example, in the PARASOL full resolution grid, if the 2 first rows are requested for the central column 3240, with dicol_min = -3, dicol_max = 3, it will return v_irow = [1, 2] and v_icol = [3237, 3238, 3239, 3240, 3241, 3242, 3243]. First line has 4 valid pixels and second has 10. The biggest one is onsidered
     * @param irow [IN] grid row indice
     * @param icol [IN] grid column indice
     * @param dirow_min [IN] number of north rows to consider (<=0)
     * @param dirow_max [IN] number of south rows to consider
     * @param dicol_min [IN] number of east columns to consider (<=0)
     * @param dicol_max [IN] number of west columns to consider
     * @param v_irow [OUT] indices of the neighbour pixels along rows (sorted increasing)
     * @param v_icol [OUT] indices of the neighbour pixels along columns (sorted increasing)
     */
    void get_neighbours_coord (const int irow, const int icol, int dirow_min, int dirow_max, int dicol_min, int dicol_max, vector<int> & v_irow, vector<int> & v_icol);
    /**
     * @brief computes the list of grid coordinates of the pixels that are neighbours of [irow, icol], in the area +/- [d_irow, d_icol]
     * All indices in this method start at 1
     * @param irow [IN] grid row indice
     * @param icol [IN] grid column indice
     * @param dicol_min [IN] number of east columns to consider (<=0)
     * @param dicol_max [IN] number of west columns to consider
     * @param v_icol [OUT] a sorted list of unique columns indices, mostly [icol + dicol_min -> icol + dicol_max], but also treats the case of crossing the change date meridian
     */
    void get_v_icol_neighbours (const int irow, const int icol, int dicol_min, int dicol_max, vector <int> & v_icol);

    /**
     * @brief accessor to the number of available directions data buffer
     * @return the number of available directions data buffer
     */
    unsigned char * get_ndir_data() const;
private:
    /**
     * @brief constucts the viewing directions observations for the given grid pixel for L1 products
     * @param ipix 2D grid indices
     * @param irec record indice matching ipix
     * @param v_obs vector of observation(s). Can contains more than one for directional products
     */
    void _get_viewing_directions_l1(const vector< int > & ipix, const int irec, vector< Observation > & v_obs);
    /**
     * @brief constucts the viewing directions observations for the given grid pixel for L23 products
     * @param ipix 2D grid indices
     * @param irec record indice matching ipix
     * @param v_obs vector of observation(s). Can contains more than one for directional products
     */
    void _get_viewing_directions_l2(const vector< int > & ipix, const int irec, vector< Observation > & v_obs);

    /**
     * @brief build the earth viewing directions of L1 products, as a segment from satellite position to viewed pixel center. All positions are stored in Earth Center Rotating carthesian coordinates.
     */
//     virtual void _load_viewing_directions_l1 ();
    /**
     * @brief build the pixels viewing directions of L2 products, as a segment from satellite position to viewed pixel center. All positions are stored in Earth Center Rotating carthesian coordinates.
     * For directionnal products (L1), a pixel is identified by the 3 indices : [irow, icolum, idirection]. For non directionnal ones, only [irow, icolum] is used
     */
//     virtual void _load_viewing_directions_l2 ();
};

#endif