n_arrays.c

Go to the documentation of this file.

/*****************************************************************************
*
* MODULE:       Grass PDE Numerical Library
* AUTHOR(S):    Soeren Gebbert, Berlin (GER) Dec 2006
*               soerengebbert <at> gmx <dot> de
*
* PURPOSE:      Array 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 <math.h>

#include <grass/N_pde.h>
#include <grass/raster.h>
#include <grass/glocale.h>


/* ******************** 2D ARRAY FUNCTIONS *********************** */

/*!
 * \brief Allocate memory for a N_array_2d data structure.
 *
 * This function allocates memory for an array of type N_array_2d
 * and returns the pointer to the new allocated memory.
 * <br><br>
 * The data type of this array is set by "type" and must be
 * CELL_TYPE, FCELL_TYPE or DCELL_TYPE accordingly to the raster map data types.
 * The offset sets the number of boundary cols and rows.
 * This option is useful to generate homogeneous Neumann boundary conditions around
 * an array or to establish overlapping boundaries. The array is initialized with 0 by default.
 * <br><br>
 * If the offset is greater then 0, negative indices are possible.
 * <br><br>
 *
 * The data structure of a array with 3 rows and cols and an offset of 1
 * will looks like this:
 * <br><br>
 *
 \verbatim
 0 0 0 0 0
 0 0 1 2 0
 0 3 4 5 0
 0 6 7 8 0
 0 0 0 0 0
 \endverbatim
 *
 * 0 is the boundary.
 * <br><br>
 * Internal a one dimensional array is allocated to save memory and to speed up the memory access.
 * To access the one dimensional array with a two dimensional index use the provided
 * get and put functions. The internal representation of the above data will look like this:
 *
 \verbatim
 0 0 0 0 0 0 0 1 2 0 0 3 4 5 0 0 6 7 8 0 0 0 0 0 0
 \endverbatim
 *
 * \param cols int
 * \param rows int
 * \param offset int
 * \param type int
 * \return N_array_2d *
 *
 * */
N_array_2d *N_alloc_array_2d(int cols, int rows, int offset, int type)
{
    N_array_2d *data = NULL;

    if (rows < 1 || cols < 1)
        G_fatal_error("N_alloc_array_2d: cols and rows should be > 0");

    if (type != CELL_TYPE && type != FCELL_TYPE && type != DCELL_TYPE)
        G_fatal_error
            ("N_alloc_array_2d: Wrong data type, should be CELL_TYPE, FCELL_TYPE or DCELL_TYPE");

    data = (N_array_2d *) G_calloc(1, sizeof(N_array_2d));

    data->cols = cols;
    data->rows = rows;
    data->type = type;
    data->offset = offset;
    data->rows_intern = rows + 2 * offset;      /*offset position at booth sides */
    data->cols_intern = cols + 2 * offset;      /*offset position at booth sides */
    data->cell_array = NULL;
    data->fcell_array = NULL;
    data->dcell_array = NULL;

    if (data->type == CELL_TYPE) {
        data->cell_array =
            (CELL *) G_calloc((size_t) data->rows_intern * data->cols_intern,
                              sizeof(CELL));
        G_debug(3,
                "N_alloc_array_2d: CELL array allocated rows_intern %i cols_intern %i offset %i",
                data->rows_intern, data->cols_intern, data->offset = offset);
    }
    else if (data->type == FCELL_TYPE) {
        data->fcell_array =
            (FCELL *) G_calloc((size_t) data->rows_intern * data->cols_intern,
                               sizeof(FCELL));
        G_debug(3,
                "N_alloc_array_2d: FCELL array allocated rows_intern %i cols_intern %i offset %i",
                data->rows_intern, data->cols_intern, data->offset = offset);

    }
    else if (data->type == DCELL_TYPE) {
        data->dcell_array =
            (DCELL *) G_calloc((size_t) data->rows_intern * data->cols_intern,
                               sizeof(DCELL));
        G_debug(3,
                "N_alloc_array_2d: DCELL array allocated rows_intern %i cols_intern %i offset %i",
                data->rows_intern, data->cols_intern, data->offset = offset);
    }

    return data;
}

/*!
 * \brief Release the memory of a N_array_2d structure
 *
 * \param data N_array_2d *
 * \return void
 * */
void N_free_array_2d(N_array_2d * data)
{

    if (data != NULL) {
        G_debug(3, "N_free_array_2d: free N_array_2d");

        if (data->type == CELL_TYPE && data->cell_array != NULL) {
            G_free(data->cell_array);
        }
        else if (data->type == FCELL_TYPE && data->fcell_array != NULL) {
            G_free(data->fcell_array);

        }
        else if (data->type == DCELL_TYPE && data->dcell_array != NULL) {
            G_free(data->dcell_array);
        }

        G_free(data);
        data = NULL;

    }

    return;
}


/*!
 * \brief Return the data type of the N_array_2d struct
 *
 * The data type can be CELL_TYPE, FCELL_TYPE or DCELL_TYPE accordingly to the raster map data types.
 *
 * \param array N_array_2d *
 * \return type int
 * */
int N_get_array_2d_type(N_array_2d * array)
{
    return array->type;
}

/*!
 * \brief Write the value of the N_array_2d struct at position col, row to value
 *
 * The value must be of the same type as the array. Otherwise you will risk data losses.
 *
 * \param data N_array_2d *
 * \param col int
 * \param row int
 * \param value void * - this variable contains the array value at col, row position
 * \return void
 * */

void N_get_array_2d_value(N_array_2d * data, int col, int row, void *value)
{

    if (data->offset == 0) {
        if (data->type == CELL_TYPE && data->cell_array != NULL) {
            *((CELL *) value) =
                data->cell_array[row * data->cols_intern + col];
        }
        else if (data->type == FCELL_TYPE && data->fcell_array != NULL) {
            *((FCELL *) value) =
                data->fcell_array[row * data->cols_intern + col];
        }
        else if (data->type == DCELL_TYPE && data->dcell_array != NULL) {
            *((DCELL *) value) =
                data->dcell_array[row * data->cols_intern + col];
        }
    }
    else {
        if (data->type == CELL_TYPE && data->cell_array != NULL) {
            *((CELL *) value) =
                data->cell_array[(row + data->offset) * data->cols_intern +
                                 col + data->offset];
        }
        else if (data->type == FCELL_TYPE && data->fcell_array != NULL) {
            *((FCELL *) value) =
                data->fcell_array[(row + data->offset) * data->cols_intern +
                                  col + data->offset];
        }
        else if (data->type == DCELL_TYPE && data->dcell_array != NULL) {
            *((DCELL *) value) =
                data->dcell_array[(row + data->offset) * data->cols_intern +
                                  col + data->offset];
        }
    }

    return;
}

/*!
 * \brief Returns 1 if the value of N_array_2d struct at position col, row
 * is of type null, otherwise 0
 *
 * This function checks automatically the type of the array and checks for the
 * data type null value.
 *
 * \param data N_array_2d *
 * \param col int
 * \param row int
 * \return int - 1 = is null, 0 otherwise
 * */
int N_is_array_2d_value_null(N_array_2d * data, int col, int row)
{

    if (data->offset == 0) {
        if (data->type == CELL_TYPE && data->cell_array != NULL) {
            G_debug(6,
                    "N_is_array_2d_value_null: null value is of type CELL at pos [%i][%i]",
                    col, row);
            return Rast_is_null_value((void *)
                                   &(data->
                                     cell_array[row * data->cols_intern +
                                                col]), CELL_TYPE);
        }
        else if (data->type == FCELL_TYPE && data->fcell_array != NULL) {
            G_debug(6,
                    "N_is_array_2d_value_null: null value is of type FCELL at pos [%i][%i]",
                    col, row);
            return Rast_is_null_value((void *)
                                   &(data->
                                     fcell_array[row * data->cols_intern +
                                                 col]), FCELL_TYPE);
        }
        else if (data->type == DCELL_TYPE && data->dcell_array != NULL) {
            G_debug(6,
                    "N_is_array_2d_value_null: null value is of type DCELL at pos [%i][%i]",
                    col, row);
            return Rast_is_null_value((void *)
                                   &(data->
                                     dcell_array[row * data->cols_intern +
                                                 col]), DCELL_TYPE);
        }
    }
    else {
        if (data->type == CELL_TYPE && data->cell_array != NULL) {
            G_debug(6,
                    "N_is_array_2d_value_null: null value is of type CELL at pos [%i][%i]",
                    col, row);
            return Rast_is_null_value((void *)
                                   &(data->
                                     cell_array[(row +
                                                 data->offset) *
                                                data->cols_intern + col +
                                                data->offset]), CELL_TYPE);
        }
        else if (data->type == FCELL_TYPE && data->fcell_array != NULL) {
            G_debug(6,
                    "N_is_array_2d_value_null: null value is of type FCELL at pos [%i][%i]",
                    col, row);
            return Rast_is_null_value((void *)
                                   &(data->
                                     fcell_array[(row +
                                                  data->offset) *
                                                 data->cols_intern + col +
                                                 data->offset]), FCELL_TYPE);
        }
        else if (data->type == DCELL_TYPE && data->dcell_array != NULL) {
            G_debug(6,
                    "N_is_array_2d_value_null: null value is of type DCELL at pos [%i][%i]",
                    col, row);
            return Rast_is_null_value((void *)
                                   &(data->
                                     dcell_array[(row +
                                                  data->offset) *
                                                 data->cols_intern + col +
                                                 data->offset]), DCELL_TYPE);
        }
    }

    return 0;
}


/*!
 * \brief Returns the value of type CELL at position col, row
 *
 * The data array can be of type CELL, FCELL or DCELL, the value will be casted to the CELL type.
 *
 * \param data N_array_2d *
 * \param col int
 * \param row int
 * \return CELL
 *
 * */
CELL N_get_array_2d_c_value(N_array_2d * data, int col, int row)
{
    CELL value = 0;
    FCELL fvalue = 0.0;
    DCELL dvalue = 0.0;

    switch (data->type) {
    case CELL_TYPE:
        N_get_array_2d_value(data, col, row, (void *)&value);
        return (CELL) value;
    case FCELL_TYPE:
        N_get_array_2d_value(data, col, row, (void *)&fvalue);
        return (CELL) fvalue;
    case DCELL_TYPE:
        N_get_array_2d_value(data, col, row, (void *)&dvalue);
        return (CELL) dvalue;
    }

    return value;
}

/*!
 * \brief Returns the value of type FCELL at position col, row
 *
 * The data array can be of type CELL, FCELL or DCELL, the value will be casted to the FCELL type.
 *
 * \param data N_array_2d *
 * \param col int
 * \param row int
 * \return FCELL

 * */
FCELL N_get_array_2d_f_value(N_array_2d * data, int col, int row)
{
    CELL value = 0;
    FCELL fvalue = 0.0;
    DCELL dvalue = 0.0;

    switch (data->type) {
    case CELL_TYPE:
        N_get_array_2d_value(data, col, row, (void *)&value);
        return (FCELL) value;
    case FCELL_TYPE:
        N_get_array_2d_value(data, col, row, (void *)&fvalue);
        return (FCELL) fvalue;
    case DCELL_TYPE:
        N_get_array_2d_value(data, col, row, (void *)&dvalue);
        return (FCELL) dvalue;
    }

    return fvalue;
}

/*!
 * \brief Returns the value of type DCELL at position col, row
 *
 * The data array can be of type CELL, FCELL or DCELL, the value will be casted to the DCELL type.
 *
 * \param data N_array_2d *
 * \param col int
 * \param row int
 * \return DCELL
 *
 * */
DCELL N_get_array_2d_d_value(N_array_2d * data, int col, int row)
{
    CELL value = 0;
    FCELL fvalue = 0.0;
    DCELL dvalue = 0.0;

    switch (data->type) {
    case CELL_TYPE:
        N_get_array_2d_value(data, col, row, (void *)&value);
        return (DCELL) value;
    case FCELL_TYPE:
        N_get_array_2d_value(data, col, row, (void *)&fvalue);
        return (DCELL) fvalue;
    case DCELL_TYPE:
        N_get_array_2d_value(data, col, row, (void *)&dvalue);
        return (DCELL) dvalue;
    }

    return dvalue;

}

/*!
 * \brief Writes a value to the N_array_2d struct at position col, row
 *
 * The value will be automatically cast to the array type.
 *
 * \param data N_array_2d *
 * \param col int
 * \param row int
 * \param value char *
 * \return void
 * */
void N_put_array_2d_value(N_array_2d * data, int col, int row, char *value)
{

    G_debug(6, "N_put_array_2d_value: put value to array");

    if (data->offset == 0) {
        if (data->type == CELL_TYPE && data->cell_array != NULL) {
            data->cell_array[row * data->cols_intern + col] =
                *((CELL *) value);
        }
        else if (data->type == FCELL_TYPE && data->fcell_array != NULL) {
            data->fcell_array[row * data->cols_intern + col] =
                *((FCELL *) value);
        }
        else if (data->type == DCELL_TYPE && data->dcell_array != NULL) {
            data->dcell_array[row * data->cols_intern + col] =
                *((DCELL *) value);
        }
    }
    else {
        if (data->type == CELL_TYPE && data->cell_array != NULL) {
            data->cell_array[(row + data->offset) * data->cols_intern + col +
                             data->offset] = *((CELL *) value);
        }
        else if (data->type == FCELL_TYPE && data->fcell_array != NULL) {
            data->fcell_array[(row + data->offset) * data->cols_intern + col +
                              data->offset] = *((FCELL *) value);
        }
        else if (data->type == DCELL_TYPE && data->dcell_array != NULL) {
            data->dcell_array[(row + data->offset) * data->cols_intern + col +
                              data->offset] = *((DCELL *) value);
        }
    }

    return;
}

/*!
 * \brief Writes the null value to the N_array_2d struct at position col, row
 *
 * The null value will be automatically set to the array data type (CELL, FCELL or DCELL).
 *
 * \param data N_array_2d *
 * \param col int
 * \param row int
 * \return void
 * */
void N_put_array_2d_value_null(N_array_2d * data, int col, int row)
{

    G_debug(6,
            "N_put_array_2d_value_null: put null value to array pos [%i][%i]",
            col, row);

    if (data->offset == 0) {
        if (data->type == CELL_TYPE && data->cell_array != NULL) {
            Rast_set_c_null_value((void *)
                               &(data->
                                 cell_array[row * data->cols_intern + col]),
                               1);
        }
        else if (data->type == FCELL_TYPE && data->fcell_array != NULL) {
            Rast_set_f_null_value((void *)
                               &(data->
                                 fcell_array[row * data->cols_intern + col]),
                               1);
        }
        else if (data->type == DCELL_TYPE && data->dcell_array != NULL) {
            Rast_set_d_null_value((void *)
                               &(data->
                                 dcell_array[row * data->cols_intern + col]),
                               1);
        }
    }
    else {
        if (data->type == CELL_TYPE && data->cell_array != NULL) {
            Rast_set_c_null_value((void *)
                               &(data->
                                 cell_array[(row +
                                             data->offset) *
                                            data->cols_intern + col +
                                            data->offset]), 1);
        }
        else if (data->type == FCELL_TYPE && data->fcell_array != NULL) {
            Rast_set_f_null_value((void *)
                               &(data->
                                 fcell_array[(row +
                                              data->offset) *
                                             data->cols_intern + col +
                                             data->offset]), 1);
        }
        else if (data->type == DCELL_TYPE && data->dcell_array != NULL) {
            Rast_set_d_null_value((void *)
                               &(data->
                                 dcell_array[(row +
                                              data->offset) *
                                             data->cols_intern + col +
                                             data->offset]), 1);
        }
    }

    return;
}

/*!
 * \brief Writes a CELL value to the N_array_2d struct at position col, row
 *
 * \param data N_array_2d *
 * \param col int
 * \param row int
 * \param value CELL
 * \return void
 * */
void N_put_array_2d_c_value(N_array_2d * data, int col, int row, CELL value)
{
    FCELL fvalue;
    DCELL dvalue;

    switch (data->type) {
    case FCELL_TYPE:
        fvalue = (FCELL) value;
        N_put_array_2d_value(data, col, row, (char *)&fvalue);
        return;
    case DCELL_TYPE:
        dvalue = (DCELL) value;
        N_put_array_2d_value(data, col, row, (char *)&dvalue);
        return;
    }

    N_put_array_2d_value(data, col, row, (char *)&value);

    return;
}

/*!
 * \brief Writes a FCELL value to the N_array_2d struct at position col, row
 *
 * \param data N_array_2d *
 * \param col int
 * \param row int
 * \param value FCELL
 * \return void
 * */
void N_put_array_2d_f_value(N_array_2d * data, int col, int row, FCELL value)
{
    CELL cvalue;
    DCELL dvalue;

    switch (data->type) {
    case CELL_TYPE:
        cvalue = (CELL) value;
        N_put_array_2d_value(data, col, row, (char *)&cvalue);
        return;
    case DCELL_TYPE:
        dvalue = (DCELL) value;
        N_put_array_2d_value(data, col, row, (char *)&dvalue);
        return;
    }

    N_put_array_2d_value(data, col, row, (char *)&value);

    return;
}

/*!
 * \brief Writes a DCELL value to the N_array_2d struct at position col, row
 *
 * \param data N_array_2d *
 * \param col int
 * \param row int
 * \param value DCELL
 * \return void
 * */
void N_put_array_2d_d_value(N_array_2d * data, int col, int row, DCELL value)
{
    CELL cvalue;
    FCELL fvalue;

    switch (data->type) {
    case CELL_TYPE:
        cvalue = (CELL) value;
        N_put_array_2d_value(data, col, row, (char *)&cvalue);
        return;
    case FCELL_TYPE:
        fvalue = (FCELL) value;
        N_put_array_2d_value(data, col, row, (char *)&fvalue);
        return;
    }

    N_put_array_2d_value(data, col, row, (char *)&value);

    return;
}

/*!
 * \brief This function writes the data info of the array data to stdout
 *
 * \param data N_array_2d *
 * \return void
 * */
void N_print_array_2d_info(N_array_2d * data)
{

    fprintf(stdout, "N_array_2d \n");
    fprintf(stdout, "Cols %i\n", data->cols);
    fprintf(stdout, "Rows: %i\n", data->rows);
    fprintf(stdout, "Array type: %i\n", data->type);
    fprintf(stdout, "Offset: %i\n", data->offset);
    fprintf(stdout, "Internal cols: %i\n", data->cols_intern);
    fprintf(stdout, "Internal rows: %i\n", data->rows_intern);
    fprintf(stdout, "CELL array pointer: %p\n", data->cell_array);
    fprintf(stdout, "FCELL array pointer: %p\n", data->fcell_array);
    fprintf(stdout, "DCELL array pointer: %p\n", data->dcell_array);


    return;
}

/*!
 * \brief Write info and content of the N_array_2d struct to stdout
 *
 * Offsets are ignored
 *
 * \param data N_array_2d *
 * \return void
 * */
void N_print_array_2d(N_array_2d * data)
{
    int i, j;

    N_print_array_2d_info(data);

    for (j = 0 - data->offset; j < data->rows + data->offset; j++) {
        for (i = 0 - data->offset; i < data->cols + data->offset; i++) {
            if (data->type == CELL_TYPE)
                fprintf(stdout, "%6d ", N_get_array_2d_c_value(data, i, j));
            else if (data->type == FCELL_TYPE)
                fprintf(stdout, "%6.6f ", N_get_array_2d_f_value(data, i, j));
            else if (data->type == DCELL_TYPE)
                printf("%6.6f ", N_get_array_2d_d_value(data, i, j));
        }
        fprintf(stdout, "\n");
    }
    fprintf(stdout, "\n");

    return;
}


/* ******************** 3D ARRAY FUNCTIONS *********************** */

/*!
 * \brief Allocate memory for a N_array_3d data structure.
 *
 * This functions allocates an array of type N_array_3d and returns a pointer
 * to the new allocated memory.
 * <br><br>
 * The data type of this array set by "type" must be
 * FCELL_TYPE or DCELL_TYPE accordingly to the raster3d map data types.
 * The offsets sets the number of boundary cols, rows and depths.
 * This option is useful to generate homogeneous Neumann boundary conditions around
 * an array or to establish overlapping boundaries. The arrays are initialized with 0 by default.
 * <br><br>
 * If the offset is greater then 0, negative indices are possible.
 * The data structure of a array with 3 depths, rows and cols and an offset of 1
 * will looks like this:
 *
 \verbatim
 0  0  0  0  0
 0  0  0  0  0
 0  0  0  0  0
 0  0  0  0  0
 0  0  0  0  0

 0  0  0  0  0
 0  0  1  2  0
 0  3  4  5  0
 0  6  7  8  0
 0  0  0  0  0

 0  0  0  0  0
 0  9 10 11  0
 0 12 13 14  0
 0 15 16 17  0
 0  0  0  0  0

 0  0  0  0  0
 0 18 19 20  0
 0 21 22 23  0
 0 24 25 26  0
 0  0  0  0  0

 0  0  0  0  0
 0  0  0  0  0
 0  0  0  0  0
 0  0  0  0  0
 0  0  0  0  0

 \endverbatim

 The depth counts from the bottom to the top.

 * <br><br>
 * Internal a one dimensional array is allocated to speed up the memory access.
 * To access the dimensional array with a three dimensional indexing use the provided
 * get and put functions.
 *
 * \param cols int
 * \param rows int
 * \param depths int
 * \param offset int
 * \param type int
 * \return N_array_3d *
 *
 * */
N_array_3d *N_alloc_array_3d(int cols, int rows, int depths, int offset,
                             int type)
{
    N_array_3d *data = NULL;

    if (rows < 1 || cols < 1 || depths < 1)
        G_fatal_error
            ("N_alloc_array_3d: depths, cols and rows should be > 0");

    if (type != DCELL_TYPE && type != FCELL_TYPE)
        G_fatal_error
            ("N_alloc_array_3d: Wrong data type, should be FCELL_TYPE or DCELL_TYPE");

    data = (N_array_3d *) G_calloc(1, sizeof(N_array_3d));

    data->cols = cols;
    data->rows = rows;
    data->depths = depths;
    data->type = type;
    data->offset = offset;
    data->rows_intern = rows + 2 * offset;
    data->cols_intern = cols + 2 * offset;
    data->depths_intern = depths + 2 * offset;
    data->fcell_array = NULL;
    data->dcell_array = NULL;

    if (data->type == FCELL_TYPE) {
        data->fcell_array =
            (float *)G_calloc((size_t) data->depths_intern * data->rows_intern *
                              data->cols_intern, sizeof(float));
        G_debug(3,
                "N_alloc_array_3d: float array allocated rows_intern %i cols_intern %i depths_intern %i offset %i",
                data->rows_intern, data->cols_intern, data->depths_intern,
                data->offset = offset);
    }
    else if (data->type == DCELL_TYPE) {
        data->dcell_array =
            (double *)G_calloc((size_t) data->depths_intern * data->rows_intern *
                               data->cols_intern, sizeof(double));
        G_debug(3,
                "N_alloc_array_3d: double array allocated rows_intern %i cols_intern %i depths_intern %i offset %i",
                data->rows_intern, data->cols_intern, data->depths_intern,
                data->offset = offset);
    }

    return data;
}

/*!
 * \brief Release the memory of a N_array_3d
 *
 * \param data N_array_3d *
 * \return void
 * */
void N_free_array_3d(N_array_3d * data)
{

    if (data != NULL) {
        G_debug(3, "N_free_array_3d: free N_array_3d");

        if (data->type == FCELL_TYPE && data->fcell_array != NULL) {
            G_free(data->fcell_array);
        }
        else if (data->type == DCELL_TYPE && data->dcell_array != NULL) {
            G_free(data->dcell_array);
        }

        G_free(data);
        data = NULL;

    }

    return;
}

/*!
 * \brief Return the data type of the N_array_3d
 *
 * The data type can be FCELL_TYPE and DCELL_TYPE accordingly to the raster map data types.
 *
 * \param array N_array_3d *
 * \return type int -- FCELL_TYPE or DCELL_TYPE
 * */
int N_get_array_3d_type(N_array_3d * array)
{
    return array->type;
}


/*!
 * \brief This function writes the value of N_array_3d data at position col, row, depth
 *        to the variable value
 *
 * The value must be from the same type as the array. Otherwise you will risk data losses.
 *
 * \param data N_array_3d *
 * \param col int
 * \param row int
 * \param depth int
 * \param value void *
 * \return void
 * */
void
N_get_array_3d_value(N_array_3d * data, int col, int row, int depth,
                     void *value)
{

    if (data->offset == 0) {
        if (data->type == FCELL_TYPE && data->fcell_array != NULL) {
            *((float *)value) =
                data->fcell_array[depth *
                                  (data->rows_intern * data->cols_intern) +
                                  row * data->cols_intern + col];
        }
        else if (data->type == DCELL_TYPE && data->dcell_array != NULL) {
            *((double *)value) =
                data->dcell_array[depth *
                                  (data->rows_intern * data->cols_intern) +
                                  row * data->cols_intern + col];
        }
    }
    else {
        if (data->type == FCELL_TYPE && data->fcell_array != NULL) {
            *((float *)value) =
                data->fcell_array[(depth + data->offset) *
                                  (data->rows_intern * data->cols_intern) +
                                  (row + data->offset) * data->cols_intern +
                                  (col + data->offset)];

        }
        else if (data->type == DCELL_TYPE && data->dcell_array != NULL) {
            *((double *)value) =
                data->dcell_array[(depth + data->offset) *
                                  (data->rows_intern * data->cols_intern) +
                                  (row + data->offset) * data->cols_intern +
                                  (col + data->offset)];
        }
    }

    return;
}

/*!
 * \brief This function returns 1 if value of N_array_3d data at position col, row, depth
 * is of type null, otherwise 0
 *
 * This function checks automatically the type of the array and checks for the
 * data type null value.
 *
 * \param data N_array_3d *
 * \param col int
 * \param row int
 * \param depth int
 * \return void
 * */
int N_is_array_3d_value_null(N_array_3d * data, int col, int row, int depth)
{

    if (data->offset == 0) {
        if (data->type == FCELL_TYPE && data->fcell_array != NULL) {
            G_debug(6,
                    "N_is_array_3d_value_null: null value is of type DCELL_TYPE at pos [%i][%i][%i]",
                    depth, row, col);
            return Rast3d_is_null_value_num((void *)
                                      &(data->
                                        fcell_array[depth *
                                                    (data->rows_intern *
                                                     data->cols_intern) +
                                                    row * data->cols_intern +
                                                    col]), FCELL_TYPE);
        }
        else if (data->type == DCELL_TYPE && data->dcell_array != NULL) {
            G_debug(6,
                    "N_is_array_3d_value_null: null value is of type DCELL_TYPE at pos [%i][%i][%i]",
                    depth, row, col);
            return Rast3d_is_null_value_num((void *)
                                      &(data->
                                        dcell_array[depth *
                                                    (data->rows_intern *
                                                     data->cols_intern) +
                                                    row * data->cols_intern +
                                                    col]), DCELL_TYPE);
        }
    }
    else {
        if (data->type == FCELL_TYPE && data->fcell_array != NULL) {
            G_debug(6,
                    "N_is_array_3d_value_null: null value is of type DCELL_TYPE at pos [%i][%i][%i]",
                    depth, row, col);
            return Rast3d_is_null_value_num((void *)
                                      &(data->
                                        fcell_array[(depth +
                                                     data->offset) *
                                                    (data->rows_intern *
                                                     data->cols_intern) +
                                                    (row + data->offset)
                                                    * data->cols_intern +
                                                    (col + data->offset)]),
                                      FCELL_TYPE);

        }
        else if (data->type == DCELL_TYPE && data->dcell_array != NULL) {
            G_debug(6,
                    "N_is_array_3d_value_null: null value is of type DCELL_TYPE at pos [%i][%i][%i]",
                    depth, row, col);
            return Rast3d_is_null_value_num((void *)
                                      &(data->
                                        dcell_array[(depth +
                                                     data->offset) *
                                                    (data->rows_intern *
                                                     data->cols_intern) +
                                                    (row +
                                                     data->offset) *
                                                    data->cols_intern + (col +
                                                                         data->
                                                                         offset)]),
                                      DCELL_TYPE);
        }
    }

    return 0;
}

/*!
 * \brief This function returns the value of type float at position col, row, depth
 *
 * The data type can be FCELL_TYPE or DCELL_TYPE accordingly to the raster map data types.
 *
 * \param data N_array_3d *
 * \param col int
 * \param row int
 * \param depth int
 * \return float
 *
 * */
float N_get_array_3d_f_value(N_array_3d * data, int col, int row, int depth)
{
    float fvalue = 0.0;
    double dvalue = 0.0;

    switch (data->type) {
    case FCELL_TYPE:
        N_get_array_3d_value(data, col, row, depth, (void *)&fvalue);
        return (float)fvalue;
    case DCELL_TYPE:
        N_get_array_3d_value(data, col, row, depth, (void *)&dvalue);
        return (float)dvalue;
    }

    return fvalue;
}

/*!
 * \brief This function returns the value of type float at position col, row, depth
 *
 * The data type can be FCELL_TYPE or DCELL_TYPE accordingly to the raster map data types.
 *
 * \param data N_array_3d *
 * \param col int
 * \param row int
 * \param depth int
 * \return double
 *
 * */
double N_get_array_3d_d_value(N_array_3d * data, int col, int row, int depth)
{
    float fvalue = 0.0;
    double dvalue = 0.0;

    switch (data->type) {

    case FCELL_TYPE:
        N_get_array_3d_value(data, col, row, depth, (void *)&fvalue);
        return (double)fvalue;
    case DCELL_TYPE:
        N_get_array_3d_value(data, col, row, depth, (void *)&dvalue);
        return (double)dvalue;
    }

    return dvalue;
}

/*!
 * \brief This function writes a value to the N_array_3d data at position col, row, depth
 *
 * The value will be automatically cast to the array type.
 *
 * \param data N_array_3d *
 * \param col int
 * \param row int
 * \param depth int
 * \param value cahr *
 * \return void
 * */
void
N_put_array_3d_value(N_array_3d * data, int col, int row, int depth,
                     char *value)
{

    G_debug(6, "N_put_array_3d_value: put value to array at pos [%i][%i][%i]",
            depth, row, col);

    if (data->offset == 0) {
        if (data->type == FCELL_TYPE && data->fcell_array != NULL) {
            data->fcell_array[depth *
                              (data->rows_intern * data->cols_intern) +
                              row * data->cols_intern + col]
                = *((float *)value);
        }
        else if (data->type == DCELL_TYPE && data->dcell_array != NULL) {

            data->dcell_array[depth *
                              (data->rows_intern * data->cols_intern) +
                              row * data->cols_intern + col]
                = *((double *)value);
        }
    }
    else {
        if (data->type == FCELL_TYPE && data->fcell_array != NULL) {
            data->fcell_array[(depth + data->offset) *
                              (data->rows_intern * data->cols_intern) + (row +
                                                                         data->
                                                                         offset)
                              * data->cols_intern + (col + data->offset)] =
                *((float *)value);
        }
        else if (data->type == DCELL_TYPE && data->dcell_array != NULL) {
            data->dcell_array[(depth + data->offset) *
                              (data->rows_intern * data->cols_intern) + (row +
                                                                         data->
                                                                         offset)
                              * data->cols_intern + (col + data->offset)] =
                *((double *)value);
        }
    }

    return;
}

/*!
 * \brief This function writes a null value to the N_array_3d data at position col, row, depth
 *
 * The null value will be automatically set to the array type.
 *
 * \param data N_array_3d *
 * \param col int
 * \param row int
 * \param depth int
 * \return void
 * */
void N_put_array_3d_value_null(N_array_3d * data, int col, int row, int depth)
{

    G_debug(6,
            "N_put_array_3d_value_null: put null value to array at pos [%i][%i][%i]",
            depth, row, col);

    if (data->offset == 0) {
        if (data->type == FCELL_TYPE && data->fcell_array != NULL) {
            Rast3d_set_null_value((void *)
                             &(data->
                               fcell_array[depth *
                                           (data->rows_intern *
                                            data->cols_intern) +
                                           row * data->cols_intern + col]), 1,
                             FCELL_TYPE);
        }
        else if (data->type == DCELL_TYPE && data->dcell_array != NULL) {
            Rast3d_set_null_value((void *)
                             &(data->
                               dcell_array[depth *
                                           (data->rows_intern *
                                            data->cols_intern) +
                                           row * data->cols_intern + col]), 1,
                             DCELL_TYPE);
        }
    }
    else {
        if (data->type == FCELL_TYPE && data->fcell_array != NULL) {
            Rast3d_set_null_value((void *)
                             &(data->
                               fcell_array[(depth +
                                            data->offset) *
                                           (data->rows_intern *
                                            data->cols_intern) + (row +
                                                                  data->
                                                                  offset) *
                                           data->cols_intern + (col +
                                                                data->
                                                                offset)]), 1,
                             FCELL_TYPE);
        }
        else if (data->type == DCELL_TYPE && data->dcell_array != NULL) {
            Rast3d_set_null_value((void *)
                             &(data->
                               dcell_array[(depth +
                                            data->offset) *
                                           (data->rows_intern *
                                            data->cols_intern) + (row +
                                                                  data->
                                                                  offset) *
                                           data->cols_intern + (col +
                                                                data->
                                                                offset)]), 1,
                             DCELL_TYPE);
        }
    }

    return;
}

/*!
 * \brief This function writes a float value to the N_array_3d data at position col, row, depth
 *
 * \param data N_array_3d *
 * \param col int
 * \param row int
 * \param depth int
 * \param value float
 * \return void
 * */
void
N_put_array_3d_f_value(N_array_3d * data, int col, int row, int depth,
                       float value)
{
    double dval;

    if (data->type == DCELL_TYPE) {
        dval = (double)value;
        N_put_array_3d_value(data, col, row, depth, (void *)&dval);
    }
    else {
        N_put_array_3d_value(data, col, row, depth, (void *)&value);
    }

    return;
}

/*!
 * \brief Writes a double value to the N_array_3d struct at position col, row, depth
 *
 * \param data N_array_3d *
 * \param col int
 * \param row int
 * \param depth int
 * \param value double
 * \return void
 * */
void
N_put_array_3d_d_value(N_array_3d * data, int col, int row, int depth,
                       double value)
{
    float fval;

    if (data->type == FCELL_TYPE) {
        fval = (double)value;
        N_put_array_3d_value(data, col, row, depth, (void *)&fval);
    }
    else {
        N_put_array_3d_value(data, col, row, depth, (void *)&value);
    }

    return;
}

/*!
 * \brief Write the info of the array to stdout
 *
 * \param data N_array_3d *
 * \return void
 * */
void N_print_array_3d_info(N_array_3d * data)
{

    fprintf(stdout, "N_array_3d \n");
    fprintf(stdout, "Cols %i\n", data->cols);
    fprintf(stdout, "Rows: %i\n", data->rows);
    fprintf(stdout, "Depths: %i\n", data->depths);
    fprintf(stdout, "Array type: %i\n", data->type);
    fprintf(stdout, "Offset: %i\n", data->offset);
    fprintf(stdout, "Internal cols: %i\n", data->cols_intern);
    fprintf(stdout, "Internal rows: %i\n", data->rows_intern);
    fprintf(stdout, "Internal depths: %i\n", data->depths_intern);
    fprintf(stdout, "FCELL array pointer: %p\n", data->fcell_array);
    fprintf(stdout, "DCELL array pointer: %p\n", data->dcell_array);

    return;
}

/*!
 * \brief Write info and content of the array data to stdout
 *
 * Offsets are ignored
 *
 * \param data N_array_2d *
 * \return void
 * */
void N_print_array_3d(N_array_3d * data)
{
    int i, j, k;

    N_print_array_3d_info(data);

    for (k = 0; k < data->depths; k++) {
        for (j = 0; j < data->rows; j++) {
            for (i = 0; i < data->cols; i++) {
                if (data->type == FCELL_TYPE)
                    printf("%6.6f ", N_get_array_3d_f_value(data, i, j, k));
                else if (data->type == DCELL_TYPE)
                    printf("%6.6f ", N_get_array_3d_d_value(data, i, j, k));
            }
            printf("\n");
        }
        printf("\n");
    }
    printf("\n");

    return;
}