grid.c

#if defined (HAVE_CONFIG_H)
#  include "config.h"
#endif

#include <string.h>

#include "dmemory.h"
#include "cdi.h"
#include "stream_int.h"
#include "grid.h"
#include "gaussgrid.h"


#ifndef  RAD2DEG
#define  RAD2DEG  (180./M_PI)   /* conversion for rad to deg */
#endif

#ifndef  DEG2RAD
#define  DEG2RAD  (M_PI/180.)   /* conversion for deg to rad */
#endif


char *Grids[] = {
  /*  0 */  "undefined",
  /*  1 */  "generic",
  /*  2 */  "gaussian",
  /*  3 */  "gaussian reduced",
  /*  4 */  "lonlat",
  /*  5 */  "spectral",
  /*  6 */  "fourier",
  /*  7 */  "gme",
  /*  8 */  "trajectory",
  /*  9 */  "unstructured",
  /* 10 */  "curvilinear",
  /* 11 */  "lcc",
  /* 12 */  "lcc2",
  /* 13 */  "laea",
  /* 14 */  "sinusoidal",
  /* 15 */  "reference",
  /* 16 */  "projection",
};


static int  GRID_Debug = 0;   /* If set to 1, debugging */

static int _grid_max = MAX_GRIDS;

static void grid_initialize(void);

static int _grid_init = FALSE;

#if  defined  (HAVE_LIBPTHREAD)
#  include <pthread.h>

static pthread_once_t _grid_init_thread = PTHREAD_ONCE_INIT;
static pthread_mutex_t _grid_mutex;

#  define GRID_LOCK()         pthread_mutex_lock(&_grid_mutex)
#  define GRID_UNLOCK()       pthread_mutex_unlock(&_grid_mutex)
#  define GRID_INIT()        \
   if ( _grid_init == FALSE ) pthread_once(&_grid_init_thread, grid_initialize)

#else

#  define GRID_LOCK()
#  define GRID_UNLOCK()
#  define GRID_INIT()        \
   if ( _grid_init == FALSE ) grid_initialize()

#endif


typedef struct _gridPtrToIdx {
  int idx;
  grid_t *ptr;
  struct _gridPtrToIdx *next;
} gridPtrToIdx;


static gridPtrToIdx *_gridList  = NULL;
static gridPtrToIdx *_gridAvail = NULL;


static
void grid_list_new(void)
{
  assert(_gridList == NULL);

  _gridList = (gridPtrToIdx *) malloc(_grid_max*sizeof(gridPtrToIdx));
}

static
void grid_list_delete(void)
{
  if ( _gridList ) free(_gridList);
}

static
void grid_init_pointer(void)
{
  int  i;
  
  for ( i = 0; i < _grid_max; i++ )
    {
      _gridList[i].next = _gridList + i + 1;
      _gridList[i].idx  = i;
      _gridList[i].ptr  = 0;
    }

  _gridList[_grid_max-1].next = 0;

  _gridAvail = _gridList;
}

static
grid_t *grid_to_pointer(int idx)
{
  grid_t *gridptr = NULL;

  GRID_INIT();

  if ( idx >= 0 && idx < _grid_max )
    {
      GRID_LOCK();

      gridptr = _gridList[idx].ptr;

      GRID_UNLOCK();
    }
  else
    Error("grid index %d undefined!", idx);

  return (gridptr);
}

/* Create an index from a pointer */
static
int grid_from_pointer(grid_t *ptr)
{
  int      idx = -1;
  gridPtrToIdx *newptr;

  if ( ptr )
    {
      GRID_LOCK();

      if ( _gridAvail )
	{
	  newptr       = _gridAvail;
	  _gridAvail   = _gridAvail->next;
	  newptr->next = 0;
	  idx	       = newptr->idx;
	  newptr->ptr  = ptr;
      
	  if ( GRID_Debug )
	    Message("Pointer %p has idx %d from grid list", ptr, idx);
	}
      else
	Warning("Too many open grids (limit is %d)!", _grid_max);

      GRID_UNLOCK();
    }
  else
    Error("Internal problem (pointer %p undefined)", ptr);

  return (idx);
}


void grid_init(grid_t *gridptr)
{
  gridptr->self         = CDI_UNDEFID;
  gridptr->type         = CDI_UNDEFID;
  gridptr->proj         = CDI_UNDEFID;
  gridptr->mask         = NULL;
  gridptr->mask_gme     = NULL;
  gridptr->xvals        = NULL;
  gridptr->yvals        = NULL;
  gridptr->area         = NULL;
  gridptr->xbounds      = NULL;
  gridptr->ybounds      = NULL;
  gridptr->rowlon       = NULL;
  gridptr->nrowlon      = 0;
  gridptr->xinc         = 0.0;
  gridptr->yinc         = 0.0;
  gridptr->lcc_originLon = 0.0;
  gridptr->lcc_originLat = 0.0;
  gridptr->lcc_lonParY  = 0.0;
  gridptr->lcc_lat1     = 0.0;
  gridptr->lcc_lat2     = 0.0;
  gridptr->lcc_xinc     = 0.0;
  gridptr->lcc_yinc     = 0.0;
  gridptr->lcc_projflag = 0;
  gridptr->lcc_scanflag = 0;
  gridptr->lcc_defined  = FALSE;
  gridptr->lcc2_lon_0   = 0.0;
  gridptr->lcc2_lat_0   = 0.0;
  gridptr->lcc2_lat_1   = 0.0;
  gridptr->lcc2_lat_2   = 0.0;
  gridptr->lcc2_a       = 0.0;
  gridptr->lcc2_defined = FALSE;
  gridptr->laea_lon_0   = 0.0;
  gridptr->laea_lat_0   = 0.0;
  gridptr->laea_a       = 0.0;
  gridptr->laea_defined = FALSE;
  gridptr->trunc        = 0;
  gridptr->nvertex      = 0;
  gridptr->nd           = 0;
  gridptr->ni           = 0;
  gridptr->ni2          = 0;
  gridptr->ni3          = 0;
  gridptr->number       = 0;
  gridptr->position     = 0;
  gridptr->reference    = NULL;
  gridptr->prec         = 0;
  gridptr->size         = 0;
  gridptr->xsize        = 0;
  gridptr->ysize        = 0;
  gridptr->np           = 0;
  gridptr->xdef         = 0;
  gridptr->ydef         = 0;
  gridptr->isCyclic     = CDI_UNDEFID;
  gridptr->isRotated    = FALSE;
  gridptr->xpole        = 0.0;
  gridptr->ypole        = 0.0;
  gridptr->angle        = 0.0;
  gridptr->locked       = FALSE;
  gridptr->lcomplex     = 0;
  gridptr->xname[0]     = 0;
  gridptr->yname[0]     = 0;
  gridptr->xlongname[0] = 0;
  gridptr->ylongname[0] = 0;
  gridptr->xunits[0]    = 0;
  gridptr->yunits[0]    = 0;
  gridptr->xstdname[0]  = 0;
  gridptr->ystdname[0]  = 0;
  gridptr->uuid[0]      = 0;
  gridptr->name         = NULL;
}


void grid_free(grid_t *gridptr)
{
  if ( gridptr->mask      ) free(gridptr->mask);
  if ( gridptr->mask_gme  ) free(gridptr->mask_gme);
  if ( gridptr->xvals     ) free(gridptr->xvals);
  if ( gridptr->yvals     ) free(gridptr->yvals);
  if ( gridptr->area      ) free(gridptr->area);
  if ( gridptr->xbounds   ) free(gridptr->xbounds);
  if ( gridptr->ybounds   ) free(gridptr->ybounds);
  if ( gridptr->rowlon    ) free(gridptr->rowlon);
  if ( gridptr->reference ) free(gridptr->reference);
  if ( gridptr->name      ) free(gridptr->name);

  grid_init(gridptr);
}

static
void grid_init_entry(grid_t *gridptr)
{
  grid_init(gridptr);

  gridptr->self = grid_from_pointer(gridptr);
}

static
grid_t *grid_new_entry(void)
{
  grid_t *gridptr;

  gridptr = (grid_t *) malloc(sizeof(grid_t));

  if ( gridptr ) grid_init_entry(gridptr);

  return (gridptr);
}

static
void grid_delete_entry(grid_t *gridptr)
{
  int idx;

  idx = gridptr->self;

  GRID_LOCK();

  free(gridptr);

  _gridList[idx].next = _gridAvail;
  _gridList[idx].ptr  = 0;
  _gridAvail          = &_gridList[idx];

  GRID_UNLOCK();

  if ( GRID_Debug )
    Message("Removed idx %d from grid list", idx);
}

static
void grid_initialize(void)
{
  char *env;

#if  defined  (HAVE_LIBPTHREAD)
  /* initialize global API mutex lock */
  pthread_mutex_init(&_grid_mutex, NULL);
#endif

  env = getenv("GRID_DEBUG");
  if ( env ) GRID_Debug = atoi(env);

  grid_list_new();
  atexit(grid_list_delete);

  GRID_LOCK();

  grid_init_pointer();

  GRID_UNLOCK();

  _grid_init = TRUE;
}

static
void grid_copy(grid_t *gridptr2, grid_t *gridptr1)
{
  int gridID2;

  gridID2 = gridptr2->self;
  memcpy(gridptr2, gridptr1, sizeof(grid_t));
  gridptr2->self = gridID2;
}

static
void gridCheckPtr(const char *caller, int gridID, grid_t *gridptr)
{
  if ( gridptr == NULL )
    Errorc("grid %d undefined!", gridID);
}

#define  grid_check_ptr(gridID, gridptr)  gridCheckPtr(__func__, gridID, gridptr)

int gridSize(void)
{
  int gridsize = 0;
  long i;
  
  GRID_INIT();

  GRID_LOCK();

  for ( i = 0; i < _grid_max; i++ )
    if ( _gridList[i].ptr ) gridsize++;

  GRID_UNLOCK();

  return (gridsize);
}


void gridGenXvals(int xsize, double xfirst, double xlast, double xinc, double *xvals)
{
  long i;

  if ( (! (fabs(xinc) > 0)) && xsize > 1 )
    {
      if ( xfirst >= xlast )
	{
	  while ( xfirst >= xlast ) xlast += 360;
	  xinc = (xlast-xfirst)/(xsize);
	}
      else
	{
	  xinc = (xlast-xfirst)/(xsize-1);
	}
    }

  for ( i = 0; i < xsize; i++ )
    xvals[i] = xfirst + i*xinc;
}

static
void calc_gaussgrid(double *yvals, int ysize, double yfirst, double ylast)
{
  double *yw;
  long yhsize;
  long i;

  yw = (double *) malloc(ysize*sizeof(double));
  gaussaw(yvals, yw, ysize);
  free(yw);
  for ( i = 0; i < ysize; i++ )
    yvals[i] = asin(yvals[i])/M_PI*180.0;

  if ( yfirst < ylast && yfirst > -90.0 && ylast < 90.0 )
    {
      double ytmp;
      yhsize = ysize/2;
      for ( i = 0; i < yhsize; i++ )
	{
	  ytmp = yvals[i];
	  yvals[i] = yvals[ysize-i-1];
	  yvals[ysize-i-1] = ytmp;
	}
    }
}


void gridGenYvals(int gridtype, int ysize, double yfirst, double ylast, double yinc, double *yvals)
{
  long i;
  double deleps = 0.002;

  if ( gridtype == GRID_GAUSSIAN || gridtype == GRID_GAUSSIAN_REDUCED )
    {
      if ( ysize > 2 )
	{
	  calc_gaussgrid(yvals, ysize, yfirst, ylast);

	  if ( ! (IS_EQUAL(yfirst, 0) && IS_EQUAL(ylast, 0)) )
	    if ( fabs(yvals[0] - yfirst) > deleps || fabs(yvals[ysize-1] - ylast) > deleps )
	      {
		double yinc = fabs(ylast-yfirst)/(ysize-1);
		double *ytmp = NULL;
		int nstart, lfound = 0;
		int ny = (int) (180./yinc + 0.5);
		ny -= ny%2;
		/* printf("%g %g %g %g %g %d\n", ylast, yfirst, ylast-yfirst,yinc, 180/yinc, ny); */
		if ( ny > ysize && ny < 4096 )
		  {
		    ytmp = (double *) malloc(ny*sizeof(double));
		    calc_gaussgrid(ytmp, ny, yfirst, ylast);
		    for ( i = 0; i < (ny-ysize); i++ )
		      if ( fabs(ytmp[i] - yfirst) < deleps ) break;

		    nstart = i;

		    if ( (nstart+ysize-1) < ny )
		      if ( fabs(ytmp[nstart+ysize-1] - ylast) < deleps ) lfound = 1;
		  }

		if ( lfound )
		  {
		    for ( i = 0; i < ysize; i++ ) yvals[i] = ytmp[i+nstart];
		  }
		else
		  {
		    Warning("Cannot calculate gaussian latitudes for lat1 = %g latn = %g!", yfirst, ylast);
		    for ( i = 0; i < ysize; i++ ) yvals[i] = 0;
		    yvals[0] = yfirst;
		    yvals[ysize-1] = ylast;
		  }

		if ( ytmp ) free(ytmp);
	      }
	}
      else
	{
	  yvals[0] = yfirst;
	  yvals[ysize-1] = ylast;
	}
    }
  /*     else if ( gridtype == GRID_LONLAT || gridtype == GRID_GENERIC ) */
  else
    {
      if ( (! (fabs(yinc) > 0)) && ysize > 1 )
	{
	  if ( IS_EQUAL(yfirst, ylast) && IS_NOT_EQUAL(yfirst, 0) ) ylast *= -1;

	  if ( yfirst > ylast )
	    yinc = (yfirst-ylast)/(ysize-1);
	  else if ( yfirst < ylast )
	    yinc = (ylast-yfirst)/(ysize-1);
	  else
	    {
	      if ( ysize%2 != 0 )
		{
		  yinc = 180.0/(ysize-1);
		  yfirst = -90;
		}
	      else
		{
		  yinc = 180.0/ysize;
		  yfirst = -90 + yinc/2;
		}
	    }
	}

      if ( yfirst > ylast && yinc > 0 ) yinc = -yinc;

      for ( i = 0; i < ysize; i++ )
	yvals[i] = yfirst + i*yinc;
    }
  /*
    else
    Error("unable to calculate values for %s grid!", gridNamePtr(gridtype));
  */
}


/*
@Function  gridCreate
@Title     Create a horizontal Grid

@Prototype int gridCreate(int gridtype, int size)
@Parameter
    @Item  gridtype  The type of the grid, one of the set of predefined CDI grid types.
                     The valid CDI grid types are @func{GRID_GENERIC}, @func{GRID_GAUSSIAN},
                     @func{GRID_LONLAT}, @func{GRID_LCC}, @func{GRID_SPECTRAL},
                     @func{GRID_GME}, @func{GRID_CURVILINEAR}, @func{GRID_UNSTRUCTURED} and
                     @func{GRID_REFERENCE}.
    @Item  size      Number of gridpoints.

@Description
The function @func{gridCreate} creates a horizontal Grid.

@Result
@func{gridCreate} returns an identifier to the Grid.

@Example
Here is an example using @func{gridCreate} to create a regular lon/lat Grid:

@Source
#include "cdi.h"
   ...
#define  nlon  12
#define  nlat   6
   ...
double lons[nlon] = {0, 30, 60, 90, 120, 150, 180, 210, 240, 270, 300, 330};
double lats[nlat] = {-75, -45, -15, 15, 45, 75};
int gridID;
   ...
gridID = gridCreate(GRID_LONLAT, nlon*nlat);
gridDefXsize(gridID, nlon);
gridDefYsize(gridID, nlat);
gridDefXvals(gridID, lons);
gridDefYvals(gridID, lats);
   ...
@EndSource
@EndFunction
*/
int gridCreate(int gridtype, int size)
{
  int gridID;
  grid_t *gridptr;

  if ( CDI_Debug )
    Message("gridtype: %d size: %d", gridtype, size);

  GRID_INIT();

  gridptr = grid_new_entry();
  if ( ! gridptr ) Error("No memory");

  gridID = gridptr->self;

  if ( CDI_Debug ) Message("gridID: %d", gridID);

  gridptr->type = gridtype;
  gridptr->size = size;

  /*  if ( gridtype == GRID_GENERIC )     gridptr->xsize = size; */
  if ( gridtype == GRID_UNSTRUCTURED )  gridptr->xsize = size;
  if ( gridtype == GRID_CURVILINEAR  )  gridptr->nvertex = 4;

  switch (gridtype)
    {
    case GRID_LONLAT:
    case GRID_GAUSSIAN:
    case GRID_GAUSSIAN_REDUCED:
    case GRID_CURVILINEAR:
    case GRID_TRAJECTORY:
      {
	if ( gridtype == GRID_TRAJECTORY )
	  {
	    gridDefXname(gridID, "tlon");
	    gridDefYname(gridID, "tlat");
	  }
	else
	  {
	    gridDefXname(gridID, "lon");
	    gridDefYname(gridID, "lat");
	  }
	gridDefXlongname(gridID, "longitude");
	gridDefYlongname(gridID, "latitude");
	/*
	if ( gridtype == GRID_CURVILINEAR )
	  {
	    strcpy(gridptr->xstdname, "grid_longitude");
	    strcpy(gridptr->ystdname, "grid_latitude");
	    gridDefXunits(gridID, "degrees");
	    gridDefYunits(gridID, "degrees");
	  }
	else
	*/
	  {
	    strcpy(gridptr->xstdname, "longitude");
	    strcpy(gridptr->ystdname, "latitude");
	    gridDefXunits(gridID, "degrees_east");
	    gridDefYunits(gridID, "degrees_north");
	  }

	break;
      }
    case GRID_GME:
    case GRID_UNSTRUCTURED:
      {
	gridDefXname(gridID, "lon");
	gridDefYname(gridID, "lat");
	strcpy(gridptr->xstdname, "longitude");
	strcpy(gridptr->ystdname, "latitude");
	gridDefXunits(gridID, "degrees_east");
	gridDefYunits(gridID, "degrees_north");
	break;
      }
    case GRID_GENERIC:
      {
	gridDefXname(gridID, "x");
	gridDefYname(gridID, "y");
	strcpy(gridptr->xstdname, "grid_longitude");
	strcpy(gridptr->ystdname, "grid_latitude");
	gridDefXunits(gridID, "degrees");
	gridDefYunits(gridID, "degrees");
	break;
      }
    case GRID_LCC2:
    case GRID_SINUSOIDAL:
    case GRID_LAEA:
      {
	gridDefXname(gridID, "x");
	gridDefYname(gridID, "y");
	strcpy(gridptr->xstdname, "projection_x_coordinate");
	strcpy(gridptr->ystdname, "projection_y_coordinate");
	gridDefXunits(gridID, "m");
	gridDefYunits(gridID, "m");
	break;
      }
    }

  return (gridID);
}


/*
@Function  gridDestroy
@Title     Destroy a horizontal Grid

@Prototype void gridDestroy(int gridID)
@Parameter
    @Item  gridID   Grid ID, from a previous call to @fref{gridCreate}.

@EndFunction
*/
void gridDestroy(int gridID)
{
  grid_t *gridptr;

  gridptr = grid_to_pointer(gridID);
  
  grid_check_ptr(gridID, gridptr);

  if ( gridptr->mask      ) free(gridptr->mask);
  if ( gridptr->mask_gme  ) free(gridptr->mask_gme);
  if ( gridptr->xvals     ) free(gridptr->xvals);
  if ( gridptr->yvals     ) free(gridptr->yvals);
  if ( gridptr->area      ) free(gridptr->area);
  if ( gridptr->xbounds   ) free(gridptr->xbounds);
  if ( gridptr->ybounds   ) free(gridptr->ybounds);
  if ( gridptr->rowlon    ) free(gridptr->rowlon);
  if ( gridptr->reference ) free(gridptr->reference);

  grid_delete_entry(gridptr);
}


char *gridNamePtr(int gridtype)
{
  char *name;
  int size = (int) (sizeof(Grids)/sizeof(char *));

  if ( gridtype >= 0 && gridtype < size )
    name = Grids[gridtype];
  else
    name = Grids[GRID_GENERIC];

  return (name);
}


void gridName(int gridtype, char *gridname)
{
  strcpy(gridname, gridNamePtr(gridtype));
}


/*
@Function  gridDefXname
@Title     Define the name of a X-axis

@Prototype void gridDefXname(int gridID, const char *name)
@Parameter
    @Item  gridID   Grid ID, from a previous call to @fref{gridCreate}.
    @Item  name     Name of the X-axis.

@Description
The function @func{gridDefXname} defines the name of a X-axis.

@EndFunction
*/
void gridDefXname(int gridID, const char *xname)
{
  grid_t *gridptr;

  gridptr = grid_to_pointer(gridID);

  grid_check_ptr(gridID, gridptr);

  if ( xname )
    strcpy(gridptr->xname, xname);
}


/*
@Function  gridDefXlongname
@Title     Define the longname of a X-axis

@Prototype void gridDefXlongname(int gridID, const char *longname)
@Parameter
    @Item  gridID   Grid ID, from a previous call to @fref{gridCreate}.
    @Item  longname Longname of the X-axis.

@Description
The function @func{gridDefXlongname} defines the longname of a X-axis.

@EndFunction
*/
void gridDefXlongname(int gridID, const char *xlongname)
{
  grid_t *gridptr;

  gridptr = grid_to_pointer(gridID);

  if ( xlongname )
    strcpy(gridptr->xlongname, xlongname);
}


/*
@Function  gridDefXunits
@Title     Define the units of a X-axis

@Prototype void gridDefXunits(int gridID, const char *units)
@Parameter
    @Item  gridID   Grid ID, from a previous call to @fref{gridCreate}.
    @Item  units    Units of the X-axis.

@Description
The function @func{gridDefXunits} defines the units of a X-axis.

@EndFunction
*/
void gridDefXunits(int gridID, const char *xunits)
{
  grid_t *gridptr;

  gridptr = grid_to_pointer(gridID);

  grid_check_ptr(gridID, gridptr);

  if ( xunits )
    strcpy(gridptr->xunits, xunits);
}


/*
@Function  gridDefYname
@Title     Define the name of a Y-axis

@Prototype void gridDefYname(int gridID, const char *name)
@Parameter
    @Item  gridID   Grid ID, from a previous call to @fref{gridCreate}.
    @Item  name     Name of the Y-axis.

@Description
The function @func{gridDefYname} defines the name of a Y-axis.

@EndFunction
*/
void gridDefYname(int gridID, const char *yname)
{
  grid_t *gridptr;

  gridptr = grid_to_pointer(gridID);

  grid_check_ptr(gridID, gridptr);

  if ( yname )
    strcpy(gridptr->yname, yname);
}


/*
@Function  gridDefYlongname
@Title     Define the longname of a Y-axis

@Prototype void gridDefYlongname(int gridID, const char *longname)
@Parameter
    @Item  gridID   Grid ID, from a previous call to @fref{gridCreate}.
    @Item  longname Longname of the Y-axis.

@Description
The function @func{gridDefYlongname} defines the longname of a Y-axis.

@EndFunction
*/
void gridDefYlongname(int gridID, const char *ylongname)
{
  grid_t *gridptr;

  gridptr = grid_to_pointer(gridID);

  grid_check_ptr(gridID, gridptr);

  if ( ylongname )
    strcpy(gridptr->ylongname, ylongname);
}


/*
@Function  gridDefYunits
@Title     Define the units of a Y-axis

@Prototype void gridDefYunits(int gridID, const char *units)
@Parameter
    @Item  gridID   Grid ID, from a previous call to @fref{gridCreate}.
    @Item  units    Units of the Y-axis.

@Description
The function @func{gridDefYunits} defines the units of a Y-axis.

@EndFunction
*/
void gridDefYunits(int gridID, const char *yunits)
{
  grid_t *gridptr;

  gridptr = grid_to_pointer(gridID);

  grid_check_ptr(gridID, gridptr);

  if ( yunits )
    strcpy(gridptr->yunits, yunits);
}


/*
@Function  gridInqXname
@Title     Get the name of a X-axis

@Prototype void gridInqXname(int gridID, char *name)
@Parameter
    @Item  gridID   Grid ID, from a previous call to @fref{gridCreate}.
    @Item  name     Name of the X-axis. The caller must allocate space for the 
                    returned string. The maximum possible length, in characters, of
                    the string is given by the predefined constant @func{CDI_MAX_NAME}.

@Description
The function @func{gridInqXname} returns the name of a X-axis.

@Result
@func{gridInqXname} returns the name of the X-axis to the parameter name.

@EndFunction
*/
void gridInqXname(int gridID, char *xname)
{
  grid_t *gridptr;

  gridptr = grid_to_pointer(gridID);

  grid_check_ptr(gridID, gridptr);

  strcpy(xname, gridptr->xname);
}


/*
@Function  gridInqXlongname
@Title     Get the longname of a X-axis

@Prototype void gridInqXlongname(int gridID, char *longname)
@Parameter
    @Item  gridID   Grid ID, from a previous call to @fref{gridCreate}.
    @Item  longname Longname of the X-axis. The caller must allocate space for the 
                    returned string. The maximum possible length, in characters, of
                    the string is given by the predefined constant @func{CDI_MAX_NAME}.

@Description
The function @func{gridInqXlongname} returns the longname of a X-axis.

@Result
@func{gridInqXlongname} returns the longname of the X-axis to the parameter longname.

@EndFunction
*/
void gridInqXlongname(int gridID, char *xlongname)
{
  grid_t *gridptr;

  gridptr = grid_to_pointer(gridID);

  grid_check_ptr(gridID, gridptr);

  strcpy(xlongname, gridptr->xlongname);
}


/*
@Function  gridInqXunits
@Title     Get the units of a X-axis

@Prototype void gridInqXunits(int gridID, char *units)
@Parameter
    @Item  gridID   Grid ID, from a previous call to @fref{gridCreate}.
    @Item  units    Units of the X-axis. The caller must allocate space for the 
                    returned string. The maximum possible length, in characters, of
                    the string is given by the predefined constant @func{CDI_MAX_NAME}.

@Description
The function @func{gridInqXunits} returns the units of a X-axis.

@Result
@func{gridInqXunits} returns the units of the X-axis to the parameter units.

@EndFunction
*/
void gridInqXunits(int gridID, char *xunits)
{
  grid_t *gridptr;

  gridptr = grid_to_pointer(gridID);

  grid_check_ptr(gridID, gridptr);

  strcpy(xunits, gridptr->xunits);
}


void gridInqXstdname(int gridID, char *xstdname)
{
  grid_t *gridptr;

  gridptr = grid_to_pointer(gridID);

  grid_check_ptr(gridID, gridptr);

  strcpy(xstdname, gridptr->xstdname);
}


/*
@Function  gridInqYname
@Title     Get the name of a Y-axis

@Prototype void gridInqYname(int gridID, char *name)
@Parameter
    @Item  gridID   Grid ID, from a previous call to @fref{gridCreate}.
    @Item  name     Name of the Y-axis. The caller must allocate space for the 
                    returned string. The maximum possible length, in characters, of
                    the string is given by the predefined constant @func{CDI_MAX_NAME}.

@Description
The function @func{gridInqYname} returns the name of a Y-axis.

@Result
@func{gridInqYname} returns the name of the Y-axis to the parameter name.

@EndFunction
*/
void gridInqYname(int gridID, char *yname)
{
  grid_t *gridptr;

  gridptr = grid_to_pointer(gridID);

  grid_check_ptr(gridID, gridptr);

  strcpy(yname, gridptr->yname);
}


/*
@Function  gridInqYlongname
@Title     Get the longname of a Y-axis

@Prototype void gridInqXlongname(int gridID, char *longname)
@Parameter
    @Item  gridID   Grid ID, from a previous call to @fref{gridCreate}.
    @Item  longname Longname of the Y-axis. The caller must allocate space for the 
                    returned string. The maximum possible length, in characters, of
                    the string is given by the predefined constant @func{CDI_MAX_NAME}.

@Description
The function @func{gridInqYlongname} returns the longname of a Y-axis.

@Result
@func{gridInqYlongname} returns the longname of the Y-axis to the parameter longname.

@EndFunction
*/
void gridInqYlongname(int gridID, char *ylongname)
{
  grid_t *gridptr;

  gridptr = grid_to_pointer(gridID);

  grid_check_ptr(gridID, gridptr);

  strcpy(ylongname, gridptr->ylongname);
}


/*
@Function  gridInqYunits
@Title     Get the units of a Y-axis

@Prototype void gridInqYunits(int gridID, char *units)
@Parameter
    @Item  gridID   Grid ID, from a previous call to @fref{gridCreate}.
    @Item  units    Units of the Y-axis. The caller must allocate space for the 
                    returned string. The maximum possible length, in characters, of
                    the string is given by the predefined constant @func{CDI_MAX_NAME}.

@Description
The function @func{gridInqYunits} returns the units of a Y-axis.

@Result
@func{gridInqYunits} returns the units of the Y-axis to the parameter units.

@EndFunction
*/
void gridInqYunits(int gridID, char *yunits)
{
  grid_t *gridptr;

  gridptr = grid_to_pointer(gridID);

  grid_check_ptr(gridID, gridptr);

  strcpy(yunits, gridptr->yunits);
}

void gridInqYstdname(int gridID, char *ystdname)
{
  grid_t *gridptr;

  gridptr = grid_to_pointer(gridID);

  grid_check_ptr(gridID, gridptr);

  strcpy(ystdname, gridptr->ystdname);
}


/*
@Function  gridInqType
@Title     Get the type of a Grid

@Prototype int gridInqType(int gridID)
@Parameter
    @Item  gridID   Grid ID, from a previous call to @fref{gridCreate}.

@Description
The function @func{gridInqType} returns the type of a Grid.

@Result
@func{gridInqType} returns the type of the grid,
one of the set of predefined CDI grid types.
The valid CDI grid types are @func{GRID_GENERIC}, @func{GRID_GAUSSIAN},
@func{GRID_LONLAT}, @func{GRID_LCC}, @func{GRID_SPECTRAL}, @func{GRID_GME},
@func{GRID_CURVILINEAR}, @func{GRID_UNSTRUCTURED} and @func{GRID_REFERENCE}.

@EndFunction
*/
int gridInqType(int gridID)
{
  grid_t *gridptr;

  gridptr = grid_to_pointer(gridID);

  grid_check_ptr(gridID, gridptr);

  return (gridptr->type);
}


/*
@Function  gridInqSize
@Title     Get the size of a Grid

@Prototype int gridInqSize(int gridID)
@Parameter
    @Item  gridID   Grid ID, from a previous call to @fref{gridCreate}.

@Description
The function @func{gridInqSize} returns the size of a Grid.

@Result
@func{gridInqSize} returns the number of grid points of a Grid.

@EndFunction
*/
int gridInqSize(int gridID)
{
  int size = 0;
  grid_t *gridptr;

  gridptr = grid_to_pointer(gridID);

  grid_check_ptr(gridID, gridptr);

  size = gridptr->size;

  if ( ! size )
    {
      int xsize, ysize;

      xsize = gridptr->xsize;
      ysize = gridptr->ysize;

      if ( ysize )
	size = xsize *ysize;
      else
	size = xsize;

      gridptr->size = size;  
    }

  return (size);
}

static
int nsp2trunc(int nsp)
{
  /*  nsp = (trunc+1)*(trunc+1)              */
  /*      => trunc^2 + 3*trunc - (x-2) = 0   */
  /*                                         */
  /*  with:  y^2 + p*y + q = 0               */
  /*         y = -p/2 +- sqrt((p/2)^2 - q)   */
  /*         p = 3 and q = - (x-2)           */
  int trunc;

  trunc = (int) (sqrt(nsp*4 + 1.) - 3) / 2;

  return (trunc);
}


int gridInqTrunc(int gridID)
{
  grid_t *gridptr;

  gridptr = grid_to_pointer(gridID);

  grid_check_ptr(gridID, gridptr);

  if ( gridptr->trunc == 0 )
    {
      if ( gridptr->type == GRID_SPECTRAL )
	gridptr->trunc = nsp2trunc(gridptr->size);
      /*
      else if      ( gridptr->type == GRID_GAUSSIAN )
	gridptr->trunc = nlat2trunc(gridptr->ysize);
      */
    }

  return (gridptr->trunc);
}


void gridDefTrunc(int gridID, int trunc)
{
  grid_t *gridptr;

  gridptr = grid_to_pointer(gridID);

  grid_check_ptr(gridID, gridptr);

  gridptr->trunc = trunc;
}


/*
@Function  gridDefXsize
@Title     Define the number of values of a X-axis

@Prototype void gridDefXsize(int gridID, int xsize)
@Parameter
    @Item  gridID   Grid ID, from a previous call to @fref{gridCreate}.
    @Item  xsize    Number of values of a X-axis.

@Description
The function @func{gridDefXsize} defines the number of values of a X-axis.

@EndFunction
*/
void gridDefXsize(int gridID, int xsize)
{
  grid_t *gridptr;

  gridptr = grid_to_pointer(gridID);

  grid_check_ptr(gridID, gridptr);

  if ( xsize > gridInqSize(gridID) )
    Error("xsize %d is greater then gridsize %d", xsize, gridInqSize(gridID));

  if ( gridInqType(gridID) == GRID_UNSTRUCTURED && xsize != gridInqSize(gridID) )
    Error("xsize %d must be equal gridsize %d for gridtype CELL", xsize, gridInqSize(gridID));

  gridptr->xsize = xsize;

  if ( gridInqType(gridID) != GRID_UNSTRUCTURED &&
       gridptr->xsize*gridptr->ysize > gridInqSize(gridID) )
    Error("inconsistent grid declaration! (xsize %d ysize %d gridsize %d)",
	  gridptr->xsize, gridptr->ysize, gridInqSize(gridID));
}


/*
@Function
@Title

@Prototype
@Parameter
    @Item  Grid identifier

@EndFunction
*/
void gridDefPrec(int gridID, int prec)
{
  grid_t *gridptr;

  gridptr = grid_to_pointer(gridID);

  grid_check_ptr(gridID, gridptr);

  gridptr->prec = prec;
}


/*
@Function  
@Title     

@Prototype 
@Parameter
    @Item  Grid identifier

@EndFunction
*/
int gridInqPrec(int gridID)
{
  grid_t *gridptr;

  gridptr = grid_to_pointer(gridID);

  grid_check_ptr(gridID, gridptr);

  return (gridptr->prec);
}


/*
@Function  gridInqXsize
@Title     Get the number of values of a X-axis

@Prototype int gridInqXsize(int gridID)
@Parameter