/* geostat.c */

/*
Copyright (C) 2006 Icare - http://www.icare.univ-lille1.fr
Fabrice Ducos, fabrice.ducos@icare.univ-lille1.fr

This program is a free software; you can redistribute it and/or
modify it under the terms of the CeCILL Public License
as published by www.cecill.info (License version 2 or later).

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
CeCILL Public License for more details.

You should have received a copy of the CeCILL Public License
along with this program; if not, please contact www.cecill.info

If you have any questions or concerns regarding this program,
please email to problem@icare.univ-lille1.fr

For general information about the ICARE Project, please email
to  icare-contact@univ-lille.fr
*/

#include <string.h> /* strlen */
#include <math.h> /* lround */
#include "geostat.h"
#include "fortran.h"
#include "msg_grid.h"

/* lround prototype is not included with math.h 
 * (maybe a bug in Linux implementation of math.h ?)
 */
extern long int lround(double x);

/* C interface for Fortran code */

#ifdef __cplusplus
# define C_LINK "C"
#else
# define C_LINK
#endif

extern C_LINK void FORTRAN_NAME(geolllp)(double *, double *, double *, double *, int *, int *);
extern C_LINK void FORTRAN_NAME(geolpll)(double *, double *, double *, double *, int *, int *);
extern C_LINK void FORTRAN_NAME(defgeo)(const char *, int *, int);
extern C_LINK void FORTRAN_NAME(geovis)(double *, double *, int *, double *, double *, int *);

geostat_err_t geostat_init(const char *satellite)
{
	int errcode;
	
	FORTRAN_NAME(defgeo)(satellite, &errcode, strlen(satellite));
	if (errcode) return GEOSTAT_ERR_BAD_INIT;
	return GEOSTAT_ERR_OK;
}

geostat_err_t geostat_latlon_to_rowcol_double(
	double latitude,
	double longitude,
	double *row,
	double *column)
{
	int errcode;
	const channel_t channel = 1; /* always 1 for MSG low resolution channels (3 km x 3km at subsatellite point) */

	FORTRAN_NAME(geolllp)((double *) &latitude, (double *) &longitude, row, column, (int *) &channel, &errcode);
	
	/* flips line and column returned by geolllp (output of geolllp gives line 1 at south and column 1 at east */
	*row = MSG_NROWS - *row + 1;
	*column = MSG_NCOLUMNS - *column + 1;
	
	return (geostat_err_t) errcode;

}

geostat_err_t geostat_rowcol_to_latlon_double(
	double row,
	double column,
	double *latitude,
	double *longitude)
{
	int errcode;
	const channel_t channel = 1; /* always 1 for MSG low resolution channels (3 km x 3km at subsatellite point) */

	/* flips line and column given geolllp (input of geolllp expects line 1 at south and column 1 at east */
	row = MSG_NROWS - row + 1;
	column = MSG_NCOLUMNS - column + 1;
	
	FORTRAN_NAME(geolpll)(&row, &column, latitude, longitude, (int *) &channel, &errcode);
	
	return (geostat_err_t) errcode;

}


geostat_err_t geostat_latlon_to_rowcol(
	const double latitude,
	const double longitude,
	unsigned long *row,
	unsigned long *column)
{
	int errcode;
	const channel_t channel = 1; /* always 1 for MSG low resolution channels (3 km x 3km at subsatellite point) */
	double dbl_row, dbl_column;
	
	FORTRAN_NAME(geolllp)((double *) &latitude, (double *) &longitude, &dbl_row, &dbl_column, (int *) &channel, &errcode);
	
	/* flips line and column returned by geolllp (output of geolllp gives line 1 at south and column 1 at east */
	*row = MSG_NROWS - lround(dbl_row) + 1;
	*column = MSG_NCOLUMNS - lround(dbl_column) + 1;
	
	return (geostat_err_t) errcode;

}


geostat_err_t geostat_rowcol_to_latlon(
	const unsigned long row,
	const unsigned long column,
	double *latitude,
	double *longitude)
{
	int errcode;
	const channel_t channel = 1; /* always 1 for MSG low resolution channels (3 km x 3km at subsatellite point) */

	/* flips line and column given geolllp (input of geolllp expects line 1 at south and column 1 at east */
	double dbl_row = (double) MSG_NROWS - row + 1;
	double dbl_column = (double) MSG_NCOLUMNS - column + 1;
	
	FORTRAN_NAME(geolpll)((double *) &dbl_row, (double *) &dbl_column, latitude, longitude, (int *) &channel, &errcode);
	
	return (geostat_err_t) errcode;

}

geostat_err_t geostat_get_azimuth_and_zenithal(
	const double latitude,
	const double longitude,
	double *azimuth,
	double *zenithal)
{
	const unsigned int npoints = 1;
	int errcode;
	
	FORTRAN_NAME(geovis)((double *) &latitude, (double *) &longitude, (int *) &npoints, azimuth, zenithal, &errcode);

	if (errcode) return GEOSTAT_ERR_BAD_INIT;
	
	return GEOSTAT_ERR_OK;

}


/*
 * These are old names for the previous functions. They are defined only for backward
 * compatibility purpose but should not be called in any new code.
 * They should be defined as inline functions in order to avoid call overload, but
 * this would limit their portability (inline is a C99 standard)
 */

geostat_err_t geostat_latlon_to_lincol_double(double latitude,
					      double longitude,
					      double *line,
					      double *column) {
  return geostat_latlon_to_rowcol_double(latitude, longitude, line, column);
}

geostat_err_t geostat_lincol_to_latlon_double(double line,
					      double column,
					      double *latitude,
					      double *longitude) {
  return geostat_rowcol_to_latlon_double(line, column, latitude, longitude);
}

geostat_err_t geostat_latlon_to_lincol(
	const double latitude,
	const double longitude,
	unsigned long *line,
	unsigned long *column)
{
  return geostat_latlon_to_rowcol(latitude, longitude, line, column);
}

geostat_err_t geostat_lincol_to_latlon(
	const unsigned long line,
	const unsigned long column,
	double *latitude,
	double *longitude)
{
  return geostat_rowcol_to_latlon(line, column, latitude, longitude);
}

geostat_err_t geostat_get_azimut_and_zenithal(
	const double latitude,
	const double longitude,
	double *azimuth,
	double *zenithal)
{
  return geostat_get_azimuth_and_zenithal(latitude, longitude, azimuth, zenithal);
}