/*****************************************************************************
 *
 * MODULE:       Grass PDE Numerical Library
 * AUTHOR(S):    Soeren Gebbert, Berlin (GER) Dec 2006
 *                 soerengebbert <at> gmx <dot> de
 *
 * PURPOSE:      gradient management functions
 *               part of the gpde library
 *
 * COPYRIGHT:    (C) 2000 by the GRASS Development Team
 *
 *               This program is free software under the GNU General Public
 *               License (>=v2). Read the file COPYING that comes with GRASS
 *               for details.
 *
 *****************************************************************************/

#include <grass/N_pde.h>

/*!
 * \brief Allocate a N_gradient_2d structure
 *
 * \return N_gradient_2d *
 *
 * */
N_gradient_2d *N_alloc_gradient_2d(void)
{
    N_gradient_2d *grad;

    grad = (N_gradient_2d *)G_calloc(1, sizeof(N_gradient_2d));

    return grad;
}

/*!
 * \brief Free's a N_gradient_2d structure
 *
 * \return void
 *
 * */
void N_free_gradient_2d(N_gradient_2d *grad)
{
    G_free(grad);
    grad = NULL;

    return;
}

/*!
 * \brief allocate and initialize a N_gradient_2d structure
 *
 * \param NC double - the gradient between northern and center cell
 * \param SC double - the gradient between southern and center cell
 * \param WC double - the gradient between western and center cell
 * \param EC double - the gradient between eastern and center cell
 * \return N_gradient_2d *
 *
 * */
N_gradient_2d *N_create_gradient_2d(double NC, double SC, double WC, double EC)
{
    N_gradient_2d *grad;

    G_debug(5, "N_create_gradient_2d: create N_gradient_2d");

    grad = N_alloc_gradient_2d();

    grad->NC = NC;
    grad->SC = SC;
    grad->WC = WC;
    grad->EC = EC;

    return grad;
}

/*!
 * \brief copy a N_gradient_2d structure
 *
 * \param source - the source N_gradient_2d struct
 * \param target - the target N_gradient_2d struct
 * \return int - 1 success, 0 failure while copying
 *
 * */
int N_copy_gradient_2d(N_gradient_2d *source, N_gradient_2d *target)
{
    G_debug(5, "N_copy_gradient_2d: copy N_gradient_2d");

    if (!source || !target)
        return 0;

    target->NC = source->NC;
    target->SC = source->SC;
    target->WC = source->WC;
    target->EC = source->EC;

    return 1;
}

/*!
 * \brief Return a N_gradient_2d structure calculated from the input gradient
 * field at position [row][col]
 *
 *  This function returns the gradient of a cell at position [row][col] from the
 * input gradient field. Returned is a new structure of type N_gradient_2d.
 *
 *  \param field N_gradient_field_2d * - A two dimensional gradient field
 *  \param gradient N_gradient_2d * - the gradient structure which should be
 * filled with data, if a NULL pointer is given, a new structure will be created
 *  \param col int
 *  \param row int
 *  \return N_gradient_2d * - the new or filled gradient structure
 *
 *
 * */
N_gradient_2d *N_get_gradient_2d(N_gradient_field_2d *field,
                                 N_gradient_2d *gradient, int col, int row)
{
    double NC = 0, SC = 0, WC = 0, EC = 0;
    N_gradient_2d *grad = gradient;

    NC = N_get_array_2d_d_value(field->y_array, col, row);
    SC = N_get_array_2d_d_value(field->y_array, col, row + 1);
    WC = N_get_array_2d_d_value(field->x_array, col, row);
    EC = N_get_array_2d_d_value(field->x_array, col + 1, row);

    G_debug(
        5, "N_get_gradient_2d: calculate N_gradient_2d NC %g SC %g WC %g EC %g",
        NC, SC, WC, EC);

    /*if gradient is a NULL pointer, create a new one */
    if (!grad) {
        grad = N_create_gradient_2d(NC, SC, WC, EC);
    }
    else {
        grad->NC = NC;
        grad->SC = SC;
        grad->WC = WC;
        grad->EC = EC;
    }

    return grad;
}

/*!
 * \brief Allocate a N_gradient_3d structure
 *
 * \return N_gradient_3d *
 *
 * */
N_gradient_3d *N_alloc_gradient_3d(void)
{
    N_gradient_3d *grad;

    grad = (N_gradient_3d *)G_calloc(1, sizeof(N_gradient_3d));

    return grad;
}

/*!
 * \brief Free's a N_gradient_3d structure
 *
 * \return void
 *
 * */
void N_free_gradient_3d(N_gradient_3d *grad)
{
    G_free(grad);
    grad = NULL;

    return;
}

/*!
 * \brief allocate and initialize a N_gradient_3d structure
 *
 * \param NC double - the gradient between northern and center cell
 * \param SC double - the gradient between southern and center cell
 * \param WC double - the gradient between western and center cell
 * \param EC double - the gradient between eastern and center cell
 * \param TC double - the gradient between top and center cell
 * \param BC double - the gradient between bottom and center cell
 * \return N_gradient_3d *
 *
 * */
N_gradient_3d *N_create_gradient_3d(double NC, double SC, double WC, double EC,
                                    double TC, double BC)
{
    N_gradient_3d *grad;

    G_debug(5, "N_create_gradient_3d: create N_gradient_3d");

    grad = N_alloc_gradient_3d();

    grad->NC = NC;
    grad->SC = SC;
    grad->WC = WC;
    grad->EC = EC;
    grad->TC = TC;
    grad->BC = BC;

    return grad;
}

/*!
 * \brief copy a N_gradient_3d structure
 *
 * \param source - the source N_gradient_3d struct
 * \param target - the target N_gradient_3d struct
 * \return int - 1 success, 0 failure while copying
 *
 * */
int N_copy_gradient_3d(N_gradient_3d *source, N_gradient_3d *target)
{
    G_debug(5, "N_copy_gradient_3d: copy N_gradient_3d");

    if (!source || !target)
        return 0;

    target->NC = source->NC;
    target->SC = source->SC;
    target->WC = source->WC;
    target->EC = source->EC;
    target->TC = source->TC;
    target->BC = source->BC;

    return 1;
}

/*!
 * \brief Return a N_gradient_3d structure calculated from the input gradient
 * field at position [depth][row][col]
 *
 *  This function returns the gradient of a 3d cell at position
 * [depth][row][col] from the input gradient field. Returned is a new structure
 * of type N_gradient_3d.
 *
 *  \param field N_gradient_field_3d * - A three dimensional gradient field
 *  \param gradient N_gradient_3d * - an existing gradient structure or a NULL
 * pointer, if a NULL pointer is providet a new structure will be returned
 *  \param col int
 *  \param row int
 *  \param depth int
 *  \return N_gradient_3d *
 *
 *
 * */
N_gradient_3d *N_get_gradient_3d(N_gradient_field_3d *field,
                                 N_gradient_3d *gradient, int col, int row,
                                 int depth)
{
    double NC, SC, WC, EC, TC, BC;
    N_gradient_3d *grad = gradient;

    NC = N_get_array_3d_d_value(field->y_array, col, row, depth);
    SC = N_get_array_3d_d_value(field->y_array, col, row + 1, depth);
    WC = N_get_array_3d_d_value(field->x_array, col, row, depth);
    EC = N_get_array_3d_d_value(field->x_array, col + 1, row, depth);
    BC = N_get_array_3d_d_value(field->z_array, col, row, depth);
    TC = N_get_array_3d_d_value(field->z_array, col, row, depth + 1);

    G_debug(6,
            "N_get_gradient_3d: calculate N_gradient_3d NC %g SC %g WC %g EC "
            "%g TC %g BC %g",
            NC, SC, WC, EC, TC, BC);

    /*if gradient is a NULL pointer, create a new one */
    if (!grad) {
        grad = N_create_gradient_3d(NC, SC, WC, EC, TC, BC);
    }
    else {
        grad->NC = NC;
        grad->SC = SC;
        grad->WC = WC;
        grad->EC = EC;
        grad->BC = BC;
        grad->TC = TC;
    }

    return grad;
}

/*!
 * \brief Allocate a N_gradient_neighbours_x structure
 *
 * This structure contains all neighbour gradients in x direction of one cell
 *
 * \return N_gradient_neighbours_x  *
 *
 * */
N_gradient_neighbours_x *N_alloc_gradient_neighbours_x(void)
{
    N_gradient_neighbours_x *grad;

    grad =
        (N_gradient_neighbours_x *)G_calloc(1, sizeof(N_gradient_neighbours_x));

    return grad;
}

/*!
 * \brief Free's a N_gradient_neighbours_x structure
 *
 * \return void
 *
 * */
void N_free_gradient_neighbours_x(N_gradient_neighbours_x *grad)
{
    G_free(grad);
    grad = NULL;

    return;
}

/*!
 * \brief Allocate and initialize a N_gradient_neighbours_x structure
 *
 * \param NWN double - the gradient between north-west and northern cell
 * \param NEN double - the gradient between north-east and northern cell
 * \param WC double - the gradient between western and center cell
 * \param EC double - the gradient between eastern and center cell
 * \param SWS double - the gradient between south-west and southern cell
 * \param SES double - the gradient between south-east and southern cell
 * \return N_gradient_neighbours_x *

 *
 * */
N_gradient_neighbours_x *N_create_gradient_neighbours_x(double NWN, double NEN,
                                                        double WC, double EC,
                                                        double SWS, double SES)
{
    N_gradient_neighbours_x *grad;

    G_debug(6,
            "N_create_gradient_neighbours_x: create N_gradient_neighbours_x");

    grad = N_alloc_gradient_neighbours_x();

    grad->NWN = NWN;
    grad->NEN = NEN;
    grad->WC = WC;
    grad->EC = EC;
    grad->SWS = SWS;
    grad->SES = SES;

    return grad;
}

/*!
 * \brief copy a N_gradient_neighbours_x structure
 *
 * \param source - the source N_gradient_neighbours_x struct
 * \param target - the target N_gradient_neighbours_x struct
 * \return int - 1 success, 0 failure while copying
 *
 * */
int N_copy_gradient_neighbours_x(N_gradient_neighbours_x *source,
                                 N_gradient_neighbours_x *target)
{
    G_debug(6, "N_copy_gradient_neighbours_x: copy N_gradient_neighbours_x");

    if (!source || !target)
        return 0;

    target->NWN = source->NWN;
    target->NEN = source->NEN;
    target->WC = source->WC;
    target->EC = source->EC;
    target->SWS = source->SWS;
    target->SES = source->SES;

    return 1;
}

/*!
 * \brief Allocate a N_gradient_neighbours_y structure
 *
 * This structure contains all neighbour gradients in y direction of one cell
 *
 * \return N_gradient_neighbours_y  *
 *
 * */
N_gradient_neighbours_y *N_alloc_gradient_neighbours_y(void)
{
    N_gradient_neighbours_y *grad;

    grad =
        (N_gradient_neighbours_y *)G_calloc(1, sizeof(N_gradient_neighbours_y));

    return grad;
}

/*!
 * \brief Free's a N_gradient_neighbours_y structure
 *
 * \return void
 *
 * */
void N_free_gradient_neighbours_y(N_gradient_neighbours_y *grad)
{
    G_free(grad);
    grad = NULL;

    return;
}

/*!
 * \brief Allocate and initialize a N_gradient_neighbours_y structure
 *
 * \param NWW double - the gradient between north-west and western cell
 * \param NEE double - the gradient between north-east and eastern cell
 * \param NC double - the gradient between northern and center cell
 * \param SC double - the gradient between southern and center cell
 * \param SWW double - the gradient between south-west and western cell
 * \param SEE double - the gradient between south-east and eastern cell
 * \return N_gradient_neighbours_y *

 *
 * */
N_gradient_neighbours_y *N_create_gradient_neighbours_y(double NWW, double NEE,
                                                        double NC, double SC,
                                                        double SWW, double SEE)
{
    N_gradient_neighbours_y *grad;

    G_debug(6,
            "N_create_gradient_neighbours_y: create N_gradient_neighbours_y");

    grad = N_alloc_gradient_neighbours_y();

    grad->NWW = NWW;
    grad->NEE = NEE;
    grad->NC = NC;
    grad->SC = SC;
    grad->SWW = SWW;
    grad->SEE = SEE;

    return grad;
}

/*!
 * \brief copy a N_gradient_neighbours_y structure
 *
 * \param source - the source N_gradient_neighbours_y struct
 * \param target - the target N_gradient_neighbours_y struct
 * \return int - 1 success, 0 failure while copying
 *
 * */
int N_copy_gradient_neighbours_y(N_gradient_neighbours_y *source,
                                 N_gradient_neighbours_y *target)
{
    G_debug(6, "N_copy_gradient_neighbours_y: copy N_gradient_neighbours_y");

    if (!source || !target)
        return 0;

    target->NWW = source->NWW;
    target->NEE = source->NEE;
    target->NC = source->NC;
    target->SC = source->SC;
    target->SWW = source->SWW;
    target->SEE = source->SEE;

    return 1;
}

/*!
 * \brief Allocate a N_gradient_neighbours_z structure
 *
 * This structure contains all neighbour gradients in z direction of one cell
 *
 * \return N_gradient_neighbours_z  *
 *
 * */
N_gradient_neighbours_z *N_alloc_gradient_neighbours_z(void)
{
    N_gradient_neighbours_z *grad;

    grad =
        (N_gradient_neighbours_z *)G_calloc(1, sizeof(N_gradient_neighbours_z));

    return grad;
}

/*!
 * \brief Free's a N_gradient_neighbours_z structure
 *
 * \return void
 *
 * */
void N_free_gradient_neighbours_z(N_gradient_neighbours_z *grad)
{
    G_free(grad);
    grad = NULL;

    return;
}

/*!
 * \brief Allocate and initialize a N_gradient_neighbours_z structure
 *
 * \param NWZ double - the gradient between upper and lower north-western cells
 * \param NZ double - the gradient between upper and lower northern cells
 * \param NEZ double - the gradient between upper and lower north-eastern cells
 * \param WZ double - the gradient between upper and lower western cells
 * \param CZ double - the gradient between upper and lower center cells
 * \param EZ double - the gradient between upper and lower eastern cells
 * \param SWZ double - the gradient between upper and lower south-western cells
 * \param SZ double - the gradient between upper and lower southern cells
 * \param SEZ double - the gradient between upper and lower south-eastern cells
 * \return N_gradient_neighbours_z *

 *
 * */
N_gradient_neighbours_z *N_create_gradient_neighbours_z(double NWZ, double NZ,
                                                        double NEZ, double WZ,
                                                        double CZ, double EZ,
                                                        double SWZ, double SZ,
                                                        double SEZ)
{
    N_gradient_neighbours_z *grad;

    G_debug(6,
            "N_create_gradient_neighbours_z: create N_gradient_neighbours_z");

    grad = N_alloc_gradient_neighbours_z();

    grad->NWZ = NWZ;
    grad->NZ = NZ;
    grad->NEZ = NEZ;
    grad->WZ = WZ;
    grad->CZ = CZ;
    grad->EZ = EZ;
    grad->SWZ = SWZ;
    grad->SZ = SZ;
    grad->SEZ = SEZ;

    return grad;
}

/*!
 * \brief copy a N_gradient_neighbours_z structure
 *
 * \param source - the source N_gradient_neighbours_z struct
 * \param target - the target N_gradient_neighbours_z struct
 * \return int - 1 success, 0 failure while copying
 *
 * */
int N_copy_gradient_neighbours_z(N_gradient_neighbours_z *source,
                                 N_gradient_neighbours_z *target)
{
    G_debug(6, "N_copy_gradient_neighbours_z: copy N_gradient_neighbours_z");

    if (!source || !target)
        return 0;

    target->NWZ = source->NWZ;
    target->NZ = source->NZ;
    target->NEZ = source->NEZ;
    target->WZ = source->WZ;
    target->CZ = source->CZ;
    target->EZ = source->EZ;
    target->SWZ = source->SWZ;
    target->SZ = source->SZ;
    target->SEZ = source->SEZ;

    return 1;
}

/*!
 * \brief Allocate a N_gradient_neighbours_2d structure
 *
 * This structure contains all neighbour gradients in all directions of one cell
 * in a 2d raster layer
 *
 * \return N_gradient_neighbours_2d *
 *
 * */
N_gradient_neighbours_2d *N_alloc_gradient_neighbours_2d(void)
{
    N_gradient_neighbours_2d *grad;

    grad = (N_gradient_neighbours_2d *)G_calloc(
        1, sizeof(N_gradient_neighbours_2d));

    grad->x = N_alloc_gradient_neighbours_x();
    grad->y = N_alloc_gradient_neighbours_y();

    return grad;
}

/*!
 * \brief Free's a N_gradient_neighbours_2d structure
 *
 * \return void
 *
 * */
void N_free_gradient_neighbours_2d(N_gradient_neighbours_2d *grad)
{

    N_free_gradient_neighbours_x(grad->x);
    N_free_gradient_neighbours_y(grad->y);

    G_free(grad);
    grad = NULL;

    return;
}

/*!
 * \brief Allocate and initialize a N_gradient_neighbours_2d structure
 *
 * The parameter N_gradient_neighbours x and y are copied into the new allocated
 * structure and can be deleted after the initializing
 *
 * \return N_gradient_neighbours_2d * -- if failure NULL is returned
 *
 * */
N_gradient_neighbours_2d *
N_create_gradient_neighbours_2d(N_gradient_neighbours_x *x,
                                N_gradient_neighbours_y *y)
{
    N_gradient_neighbours_2d *grad;
    int fail = 0;

    G_debug(5,
            "N_create_gradient_neighbours_2d: create N_gradient_neighbours_2d");

    grad = N_alloc_gradient_neighbours_2d();

    if (!N_copy_gradient_neighbours_x(x, grad->x))
        fail++;
    if (!N_copy_gradient_neighbours_y(y, grad->y))
        fail++;

    if (fail > 0) {
        N_free_gradient_neighbours_2d(grad);
        grad = NULL;
    }

    return grad;
}

/*!
 * \brief copy a N_gradient_neighbours_2d structure
 *
 * \param source - the source N_gradient_neighbours_2d struct
 * \param target - the target N_gradient_neighbours_2d struct
 * \return int - 1 success, 0 failure while copying
 *
 * */
int N_copy_gradient_neighbours_2d(N_gradient_neighbours_2d *source,
                                  N_gradient_neighbours_2d *target)
{
    int fail = 0;

    G_debug(5, "N_copy_gradient_neighbours_2d: copy N_gradient_neighbours_2d");

    if (!source || !target)
        return 0;

    if (!(N_copy_gradient_neighbours_x(source->x, target->x)))
        fail++;
    if (!(N_copy_gradient_neighbours_y(source->y, target->y)))
        fail++;

    if (fail > 0) {
        return 0;
    }

    return 1;
}

/*!
 * \brief Return a N_gradient_neighbours_2d structure calculated from the input
 * gradient field at position [row][col]
 *
 *  This function returns the gradient neighbours in x and y dierection
 *  of a cell at position [row][col] from the input gradient field.
 *  Returned is a pointer to a structure of type N_gradient_neighbours_2d.
 *
 *  \param field N_gradient_field_2d * - A two dimensional gradient field
 *  \param gradient N_gradient_neighbours_2d * - the gradient structure which
 * should be filled with data, if a NULL pointer is given, a new structure will
 * be created \param col int \param row int \return N_gradient_neighbours_2d * -
 * the new or filled gradient structure
 *
 *
 * */
N_gradient_neighbours_2d *
N_get_gradient_neighbours_2d(N_gradient_field_2d *field,
                             N_gradient_neighbours_2d *gradient, int col,
                             int row)
{
    double NWN, NEN, WC, EC, SWS, SES;
    double NWW, NEE, NC, SC, SWW, SEE;
    N_gradient_neighbours_2d *grad = NULL;
    N_gradient_neighbours_x *grad_x = NULL;
    N_gradient_neighbours_y *grad_y = NULL;

    NWN = N_get_array_2d_d_value(field->x_array, col, row - 1);
    NEN = N_get_array_2d_d_value(field->x_array, col + 1, row - 1);
    WC = N_get_array_2d_d_value(field->x_array, col, row);
    EC = N_get_array_2d_d_value(field->x_array, col + 1, row);
    SWS = N_get_array_2d_d_value(field->x_array, col, row + 1);
    SES = N_get_array_2d_d_value(field->x_array, col + 1, row + 1);

    NWW = N_get_array_2d_d_value(field->y_array, col - 1, row);
    NEE = N_get_array_2d_d_value(field->y_array, col + 1, row);
    NC = N_get_array_2d_d_value(field->y_array, col, row);
    SC = N_get_array_2d_d_value(field->y_array, col, row + 1);
    SWW = N_get_array_2d_d_value(field->y_array, col - 1, row + 1);
    SEE = N_get_array_2d_d_value(field->y_array, col + 1, row + 1);

    grad_x = N_create_gradient_neighbours_x(NWN, NEN, WC, EC, SWS, SES);
    grad_y = N_create_gradient_neighbours_y(NWW, NEE, NC, SC, SWW, SEE);

    G_debug(5,
            "N_get_gradient_neighbours_2d: calculate N_gradient_neighbours_x "
            "NWN %g NEN %g WC %g EC %g SWS %g SES %g",
            NWN, NEN, WC, EC, SWS, SES);

    G_debug(5,
            "N_get_gradient_neighbours_2d: calculate N_gradient_neighbours_y "
            "NWW %g NEE %g NC %g SC %g SWW %g SEE %g",
            NWW, NEE, NC, SC, SWW, SEE);

    /*if gradient is a NULL pointer, create a new one */
    if (!gradient) {
        grad = N_create_gradient_neighbours_2d(grad_x, grad_y);
        gradient = grad;
    }
    else {
        grad = N_create_gradient_neighbours_2d(grad_x, grad_y);
        N_copy_gradient_neighbours_2d(grad, gradient);
        N_free_gradient_neighbours_2d(grad);
    }

    N_free_gradient_neighbours_x(grad_x);
    N_free_gradient_neighbours_y(grad_y);

    return gradient;
}

/*!
 * \brief Allocate a N_gradient_neighbours_3d structure
 *
 * This structure contains all neighbour gradients in all directions of one cell
 * in a 3d raster layer
 *
 * \return N_gradient_neighbours_3d *
 *
 * */
N_gradient_neighbours_3d *N_alloc_gradient_neighbours_3d(void)
{
    N_gradient_neighbours_3d *grad;

    grad = (N_gradient_neighbours_3d *)G_calloc(
        1, sizeof(N_gradient_neighbours_3d));

    grad->xt = N_alloc_gradient_neighbours_x();
    grad->xc = N_alloc_gradient_neighbours_x();
    grad->xb = N_alloc_gradient_neighbours_x();
    grad->yt = N_alloc_gradient_neighbours_y();
    grad->yc = N_alloc_gradient_neighbours_y();
    grad->yb = N_alloc_gradient_neighbours_y();
    grad->zt = N_alloc_gradient_neighbours_z();
    grad->zb = N_alloc_gradient_neighbours_z();

    return grad;
}

/*!
 * \brief Free's a N_gradient_neighbours_3d structure
 *
 * \return void
 *
 * */
void N_free_gradient_neighbours_3d(N_gradient_neighbours_3d *grad)
{

    N_free_gradient_neighbours_x(grad->xt);
    N_free_gradient_neighbours_x(grad->xc);
    N_free_gradient_neighbours_x(grad->xb);
    N_free_gradient_neighbours_y(grad->yt);
    N_free_gradient_neighbours_y(grad->yc);
    N_free_gradient_neighbours_y(grad->yb);
    N_free_gradient_neighbours_z(grad->zt);
    N_free_gradient_neighbours_z(grad->zb);

    G_free(grad);
    grad = NULL;

    return;
}

/*!
 * \brief Allocate and initialize a N_gradient_neighbours_3d structure
 *
 * The parameter N_gradient_neighbours x(tcb) and y(tcb) and z(tb) are copied
 into the new allocated structure
 * and can be deleted after the initializing
 *
 * \return N_gradient_neighbours_3d * -- if failure NULL is returned

 *
 * */
N_gradient_neighbours_3d *N_create_gradient_neighbours_3d(
    N_gradient_neighbours_x *xt, N_gradient_neighbours_x *xc,
    N_gradient_neighbours_x *xb, N_gradient_neighbours_y *yt,
    N_gradient_neighbours_y *yc, N_gradient_neighbours_y *yb,
    N_gradient_neighbours_z *zt, N_gradient_neighbours_z *zb)
{
    N_gradient_neighbours_3d *grad;
    int fail = 0;

    G_debug(5,
            "N_create_gradient_neighbours_3d: create N_gradient_neighbours_3d");

    grad = N_alloc_gradient_neighbours_3d();

    if (!(N_copy_gradient_neighbours_x(xt, grad->xt)))
        fail++;
    if (!(N_copy_gradient_neighbours_x(xc, grad->xc)))
        fail++;
    if (!(N_copy_gradient_neighbours_x(xb, grad->xb)))
        fail++;
    if (!(N_copy_gradient_neighbours_y(yt, grad->yt)))
        fail++;
    if (!(N_copy_gradient_neighbours_y(yc, grad->yc)))
        fail++;
    if (!(N_copy_gradient_neighbours_y(yb, grad->yb)))
        fail++;
    if (!(N_copy_gradient_neighbours_z(zt, grad->zt)))
        fail++;
    if (!(N_copy_gradient_neighbours_z(zb, grad->zb)))
        fail++;

    if (fail > 0) {
        return NULL;
    }

    return grad;
}

/*!
 * \brief copy a N_gradient_neighbours_3d structure
 *
 * \param source - the source N_gradient_neighbours_3d struct
 * \param target - the target N_gradient_neighbours_3d struct
 * \return int - 1 success, 0 failure while copying
 *
 * */
int N_copy_gradient_neighbours_3d(N_gradient_neighbours_3d *source,
                                  N_gradient_neighbours_3d *target)
{
    int fail = 0;

    G_debug(5, "N_copy_gradient_neighbours_3d: copy N_gradient_neighbours_3d");

    if (!source || !target)
        return 0;

    if (!(N_copy_gradient_neighbours_x(source->xt, target->xt)))
        fail++;
    if (!(N_copy_gradient_neighbours_x(source->xc, target->xc)))
        fail++;
    if (!(N_copy_gradient_neighbours_x(source->xb, target->xb)))
        fail++;
    if (!(N_copy_gradient_neighbours_y(source->yt, target->yt)))
        fail++;
    if (!(N_copy_gradient_neighbours_y(source->yc, target->yc)))
        fail++;
    if (!(N_copy_gradient_neighbours_y(source->yb, target->yb)))
        fail++;
    if (!(N_copy_gradient_neighbours_z(source->zt, target->zt)))
        fail++;
    if (!(N_copy_gradient_neighbours_z(source->zb, target->zb)))
        fail++;

    if (fail > 0) {
        return 0;
    }

    return 1;
}

/*!
 * \brief Allocate a N_gradient_field_2d
 *
 * The field arrays are of type DCELL.
 *
 * \param rows - number of rows of the 2d array from which the gradient should
 * be calculated \param cols - number of cols of the 2d array from which the
 * gradient should be calculated \return N_gradient_field_2d *
 *
 * */
N_gradient_field_2d *N_alloc_gradient_field_2d(int cols, int rows)
{
    N_gradient_field_2d *field;

    G_debug(5,
            "N_alloc_gradient_field_2d: allocate a N_gradient_field_2d struct");

    field = (N_gradient_field_2d *)G_calloc(1, sizeof(N_gradient_field_2d));

    field->x_array = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);
    field->y_array = N_alloc_array_2d(cols, rows, 1, DCELL_TYPE);

    field->cols = cols;
    field->rows = rows;

    return field;
}

/*!
 * \brief Free's a N_gradient_neighbours_2d structure
 *
 * \return void
 *
 * */
void N_free_gradient_field_2d(N_gradient_field_2d *field)
{

    N_free_array_2d(field->x_array);
    N_free_array_2d(field->y_array);

    G_free(field);

    field = NULL;

    return;
}

/*!
 * \brief Copy N_gradient_field_2d structure from source to target
 *
 * \param source - the source N_gradient_field_2d struct
 * \param target - the target N_gradient_field_2d struct
 * \return int - 1 success, 0 failure while copying
 *
 * */
int N_copy_gradient_field_2d(N_gradient_field_2d *source,
                             N_gradient_field_2d *target)
{
    G_debug(3, "N_copy_gradient_field_2d: copy N_gradient_field_2d");

    if (!source || !target)
        return 0;

    N_copy_array_2d(source->x_array, target->x_array);
    N_copy_array_2d(source->y_array, target->y_array);

    return 1;
}

/*! \brief Print gradient field information to stdout
 *
 * \param field N_gradient_2d_field *
 * \return void
 *
 * */
void N_print_gradient_field_2d_info(N_gradient_field_2d *field)
{
    fprintf(stdout, "N_gradient_field_2d \n");
    fprintf(stdout, "Cols %i\n", field->cols);
    fprintf(stdout, "Rows: %i\n", field->rows);
    fprintf(stdout, "X array pointer: %p\n", (void *)field->x_array);
    fprintf(stdout, "Y array pointer: %p\n", (void *)field->y_array);
    fprintf(stdout, "Min %g\n", field->min);
    fprintf(stdout, "Max %g\n", field->max);
    fprintf(stdout, "Sum %g\n", field->sum);
    fprintf(stdout, "Mean %g\n", field->mean);
    fprintf(stdout, "Nonull %i\n", field->nonull);
    fprintf(stdout, "X array info \n");
    N_print_array_2d_info(field->x_array);
    fprintf(stdout, "Y array info \n");
    N_print_array_2d_info(field->y_array);

    return;
}

/*!
 * \brief Allocate a N_gradient_field_3d
 *
 * The field arrays are always of type DCELL_TYPE.
 *
 * \param cols - number of cols of the 3d array from which the gradient should
 * be calculated \param rows - number of rows of the 3d array from which the
 * gradient should be calculated \param depths - number of depths of the 3d
 * array from which the gradient should be calculated \return
 * N_gradient_field_3d *
 *
 * */
N_gradient_field_3d *N_alloc_gradient_field_3d(int cols, int rows, int depths)
{
    N_gradient_field_3d *field;

    G_debug(5,
            "N_alloc_gradient_field_3d: allocate a N_gradient_field_3d struct");

    field = (N_gradient_field_3d *)G_calloc(1, sizeof(N_gradient_field_3d));

    field->x_array = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
    field->y_array = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);
    field->z_array = N_alloc_array_3d(cols, rows, depths, 1, DCELL_TYPE);

    field->cols = cols;
    field->rows = rows;
    field->depths = depths;

    return field;
}

/*!
 * \brief Free's a N_gradient_neighbours_3d structure
 *
 * \return void
 *
 * */
void N_free_gradient_field_3d(N_gradient_field_3d *field)
{

    N_free_array_3d(field->x_array);
    N_free_array_3d(field->y_array);
    N_free_array_3d(field->z_array);

    G_free(field);

    field = NULL;

    return;
}

/*!
 * \brief Copy N_gradient_field_3d structure from source to target
 *
 * \param source - the source N_gradient_field_3d struct
 * \param target - the target N_gradient_field_3d struct
 * \return int - 1 success, 0 failure while copying
 *
 * */
int N_copy_gradient_field_3d(N_gradient_field_3d *source,
                             N_gradient_field_3d *target)
{
    G_debug(3, "N_copy_gradient_field_3d: copy N_gradient_field_3d");

    if (!source || !target)
        return 0;

    N_copy_array_3d(source->x_array, target->x_array);
    N_copy_array_3d(source->y_array, target->y_array);
    N_copy_array_3d(source->z_array, target->z_array);

    return 1;
}

/*! \brief Print gradient field information to stdout
 *
 * \param field N_gradient_3d_field *
 * \return void
 *
 * */
void N_print_gradient_field_3d_info(N_gradient_field_3d *field)
{

    fprintf(stdout, "N_gradient_field_3d \n");
    fprintf(stdout, "Cols %i\n", field->cols);
    fprintf(stdout, "Rows: %i\n", field->rows);
    fprintf(stdout, "Depths %i\n", field->depths);
    fprintf(stdout, "X array pointer: %p\n", (void *)field->x_array);
    fprintf(stdout, "Y array pointer: %p\n", (void *)field->y_array);
    fprintf(stdout, "Z array pointer: %p\n", (void *)field->z_array);
    fprintf(stdout, "Min %g\n", field->min);
    fprintf(stdout, "Max %g\n", field->max);
    fprintf(stdout, "Sum %g\n", field->sum);
    fprintf(stdout, "Mean %g\n", field->mean);
    fprintf(stdout, "Nonull %i\n", field->nonull);
    fprintf(stdout, "X array info \n");
    N_print_array_3d_info(field->x_array);
    fprintf(stdout, "Y array info \n");
    N_print_array_3d_info(field->y_array);
    fprintf(stdout, "Z array info \n");
    N_print_array_3d_info(field->z_array);

    return;
}