/**
* Copyright 2005-2007 ECMWF
*
* Licensed under the GNU Lesser General Public License which
* incorporates the terms and conditions of version 3 of the GNU
* General Public License.
* See LICENSE and gpl-3.0.txt for details.
*/

/**
* Author: Enrico Fucile
* date: 31/10/2007
*
*/

#include "grib_api_internal.h"


int grib_nearest_find(grib_nearest *nearest,grib_handle* h,double inlat,double inlon,
    unsigned long flags,double* outlats,double* outlons,
    double* values,double* distances, int* indexes, size_t *len)
{
  grib_nearest_class *c = nearest->cclass;
  while(c)
  {
    grib_nearest_class *s = c->super ? *(c->super) : NULL;
    if(c->find) {
		int ret=c->find(nearest,h,inlat,inlon,flags,outlats,outlons,values,distances,indexes,len);
		if (ret!=GRIB_SUCCESS) {
			if (inlon>0) inlon-=360;
			else inlon+=360;
			ret=c->find(nearest,h,inlat,inlon,flags,outlats,outlons,values,distances,indexes,len);
		}
		return ret;
	}
    c = s;
  }
  Assert(0);
  return 0;
}


/* For this one, ALL init are called */

static int init_nearest(grib_nearest_class* c,grib_nearest* i, grib_handle *h, grib_arguments* args)
{

  if(c) {
    int ret = GRIB_SUCCESS;
    grib_nearest_class *s = c->super ? *(c->super) : NULL;
    if(!c->inited)
    {
      if(c->init_class) c->init_class(c);
      c->inited = 1;
    }
    if(s) ret = init_nearest(s,i,h,args);

    if(ret != GRIB_SUCCESS) return ret;

    if(c->init) return c->init(i,h, args);
  }
  return GRIB_INTERNAL_ERROR;
}

int grib_nearest_init(grib_nearest* i, grib_handle *h, grib_arguments* args)
{
  return init_nearest(i->cclass,i,h,args);
}

/* For this one, ALL destroy are called */

int grib_nearest_delete(grib_nearest *i)
{
  grib_nearest_class *c = i->cclass;
  while(c)
  {
    grib_nearest_class *s = c->super ? *(c->super) : NULL;
    if(c->destroy) c->destroy(i);
    c = s;
  }
  return 0;
}


void grib_binary_search(double xx[], const unsigned long n, double x,
    int *ju,int *jl)
{
  size_t jm=0;
  int ascending=0;
  *jl=0;
  *ju=n;
  ascending=(xx[n] >= xx[0]);
  while (*ju-*jl > 1) {
    jm=(*ju+*jl) >> 1;
    if ((x >= xx[jm]) == ascending) *jl=jm;
    else *ju=jm;
  }
}
		
#define RADIAN(x) ((x) * acos(0.0) / 90.0)

double grib_nearest_distance(double radius,double lon1, double lat1, double lon2, double lat2){
    double rlat1=RADIAN(lat1);
    double rlat2=RADIAN(lat2);
    double rlon1=lon1;
    double rlon2=lon2;
	double a;

    if (rlon1 >=360) rlon1-=360.0;
    rlon1=RADIAN(rlon1);
    if (rlon2 >=360) rlon2-=360.0;
    rlon2=RADIAN(rlon2);

	a = sin(rlat1)*sin(rlat2)+ cos(rlat1)*cos(rlat2)*cos(rlon2-rlon1);

	if (a > 1 || a < -1 ) a=(int)a;

	return radius*acos(a);
}

int grib_nearest_find_multiple(grib_handle* h,int is_lsm,
    double* inlats,double* inlons,long npoints,
    double* outlats,double* outlons,
    double* values,double* distances, int* indexes)
{
  grib_nearest* nearest=0;
  double* pdistances=distances;
  double* poutlats=outlats;
  double* poutlons=outlons;
  double* pvalues=values;
  int* pindexes=indexes;
  int idx=0,ii=0;
  double max,min;
  double qdistances[4]={0,};
  double qoutlats[4]={0,};
  double qoutlons[4]={0,};
  double qvalues[4]={0,};
  int qindexes[4]={0,};
  int ret=0;
  long i=0;
  size_t len=4;
  int flags=GRIB_NEAREST_SAME_GRID | GRIB_NEAREST_SAME_DATA;

  nearest=grib_nearest_new(h,&ret);
  if (ret!=GRIB_SUCCESS) return ret;

  if (is_lsm) {
    int noland=1;
    for (i=0;i<npoints;i++) {
      ret=grib_nearest_find(nearest,h,inlats[i],inlons[i],flags,qoutlats,qoutlons,
                         qvalues,qdistances,qindexes,&len);
      max=qdistances[0];
      for (ii=0;ii<4;ii++) {
        if (max<qdistances[ii]) {max=qdistances[ii];idx=ii;}
        if (qvalues[ii] >= 0.5) noland=0;
      }
      min=max;
      for (ii=0;ii<4;ii++) {
        if ((min >= qdistances[ii]) && ( noland || (qvalues[ii] >= 0.5)) ) {
          min = qdistances[ii];
          idx=ii;
        }
      }
      *poutlats=qoutlats[idx];poutlats++;
      *poutlons=qoutlons[idx];poutlons++;
      *pvalues=qvalues[idx];pvalues++;
      *pdistances=qdistances[idx];pdistances++;
      *pindexes=qindexes[idx];pindexes++;

    }
  } else {
    for (i=0;i<npoints;i++) {
       ret=grib_nearest_find(nearest,h,inlats[i],inlons[i],flags,qoutlats,qoutlons,
                         qvalues,qdistances,qindexes,&len);
      min=qdistances[0];
      for (ii=0;ii<4;ii++) {
        if ((min >= qdistances[ii])) {
          min = qdistances[ii];
          idx=ii;
        }
      }
      *poutlats=qoutlats[idx];poutlats++;
      *poutlons=qoutlons[idx];poutlons++;
      *pvalues=qvalues[idx];pvalues++;
      *pdistances=qdistances[idx];pdistances++;
      *pindexes=qindexes[idx];pindexes++;
    }
  }

  grib_nearest_delete(nearest);

  return ret;
}