1 
2 /*****************************************************************************
3 *
4 * MODULE:       Grass PDE Numerical Library
5 * AUTHOR(S):    Soeren Gebbert, Berlin (GER) Dec 2006
6 * 		soerengebbert <at> gmx <dot> de
7 *
8 * PURPOSE:      This file contains definitions of variables and data types
9 *
10 * COPYRIGHT:    (C) 2000 by the GRASS Development Team
11 *
12 *               This program is free software under the GNU General Public
13 *               License (>=v2). Read the file COPYING that comes with GRASS
14 *               for details.
15 *
16 *****************************************************************************/
17 
18 #include <grass/gis.h>
19 #include <grass/raster3d.h>
20 #include <grass/glocale.h>
21 #include <grass/gmath.h>
22 
23 #ifndef _N_PDE_H_
24 #define _N_PDE_H_
25 
26 #define N_NORMAL_LES 0
27 #define N_SPARSE_LES 1
28 /*!
29  * Boundary conditions for cells
30  */
31 #define N_CELL_INACTIVE 0
32 #define N_CELL_ACTIVE 1
33 #define N_CELL_DIRICHLET 2
34 #define N_CELL_TRANSMISSION 3
35 /*!
36  * \brief the maximum number of available cell states (eg: boundary condition, inactiven active)
37  * */
38 #define N_MAX_CELL_STATE 20
39 
40 #define N_5_POINT_STAR 0
41 #define N_7_POINT_STAR 1
42 #define N_9_POINT_STAR 2
43 #define N_27_POINT_STAR 3
44 
45 #define N_MAXIMUM_NORM 0
46 #define N_EUKLID_NORM 1
47 
48 #define N_ARRAY_SUM 0		/* summ two arrays */
49 #define N_ARRAY_DIF 1		/* calc the difference between two arrays */
50 #define N_ARRAY_MUL 2		/* multiply two arrays */
51 #define N_ARRAY_DIV 3		/* array division, if div with 0 the NULL value is set */
52 
53 #define N_UPWIND_FULL 0		/*full upwinding stabilization */
54 #define N_UPWIND_EXP  1		/*exponential upwinding stabilization */
55 #define N_UPWIND_WEIGHT 2	/*weighted upwinding stabilization */
56 
57 
58 
59 /* *************************************************************** */
60 /* *************** LINEARE EQUATION SYSTEM PART ****************** */
61 /* *************************************************************** */
62 
63 /*!
64  * \brief The linear equation system (les) structure
65  *
66  * This structure manages the Ax = b system.
67  * It manages regular quadratic matrices or
68  * sparse matrices. The vector b and x are normal one dimensional
69  * memory structures of type double. Also the number of rows
70  * and the matrix type are stored in this structure.
71  * */
72 typedef struct
73 {
74     double *x;			/*the value vector */
75     double *b;			/*the right side of Ax = b */
76     double **A;			/*the normal quadratic matrix */
77     G_math_spvector **Asp;	/*the sparse matrix */
78     int rows;			/*number of rows */
79     int cols;			/*number of cols */
80     int quad;			/*is the matrix quadratic (1-quadratic, 0 not) */
81     int type;			/*the type of the les, normal == 0, sparse == 1 */
82 } N_les;
83 
84 extern N_les *N_alloc_les_param(int cols, int rows, int type, int param);
85 extern N_les *N_alloc_les(int rows, int type);
86 extern N_les *N_alloc_les_A(int rows, int type);
87 extern N_les *N_alloc_les_Ax(int rows, int type);
88 extern N_les *N_alloc_les_Ax_b(int rows, int type);
89 extern N_les *N_alloc_nquad_les(int cols, int rows, int type);
90 extern N_les *N_alloc_nquad_les_A(int cols, int rows, int type);
91 extern N_les *N_alloc_nquad_les_Ax(int cols, int rows, int type);
92 extern N_les *N_alloc_nquad_les_Ax_b(int cols, int rows, int type);
93 extern void N_print_les(N_les * les);
94 extern void N_free_les(N_les * les);
95 
96 /* *************************************************************** */
97 /* *************** GEOMETRY INFORMATION ************************** */
98 /* *************************************************************** */
99 
100 /*!
101  * \brief Geometric information about the structured grid
102  * */
103 typedef struct
104 {
105     int planimetric;		/*If the projection is not planimetric (0), the array calculation is different for each row */
106     double *area;		/* the vector of area values for non-planimetric projection for each row */
107     int dim;			/* 2 or 3 */
108 
109     double dx;
110     double dy;
111     double dz;
112 
113     double Az;
114 
115     int depths;
116     int rows;
117     int cols;
118 
119 } N_geom_data;
120 
121 extern N_geom_data *N_alloc_geom_data(void);
122 extern void N_free_geom_data(N_geom_data * geodata);
123 extern N_geom_data *N_init_geom_data_3d(RASTER3D_Region * region3d, N_geom_data * geodata);
124 extern N_geom_data *N_init_geom_data_2d(struct Cell_head *region, N_geom_data * geodata);
125 extern double N_get_geom_data_area_of_cell(N_geom_data * geom, int row);
126 
127 /* *************************************************************** */
128 /* *************** READING RASTER AND VOLUME DATA **************** */
129 /* *************************************************************** */
130 
131 typedef struct
132 {
133     int type;			/* which raster type CELL_TYPE, FCELL_TYPE, DCELL_TYPE */
134     int rows, cols;
135     int rows_intern, cols_intern;
136     int offset;			/*number of cols/rows offset at each boundary */
137     CELL *cell_array;		/*The data is stored in an one dimensional array internally */
138     FCELL *fcell_array;		/*The data is stored in an one dimensional array internally */
139     DCELL *dcell_array;		/*The data is stored in an one dimensional array internally */
140 } N_array_2d;
141 
142 extern N_array_2d *N_alloc_array_2d(int cols, int rows, int offset, int type);
143 extern void N_free_array_2d(N_array_2d * data_array);
144 extern int N_get_array_2d_type(N_array_2d * array2d);
145 extern void N_get_array_2d_value(N_array_2d * array2d, int col, int row, void *value);
146 extern CELL N_get_array_2d_c_value(N_array_2d * array2d, int col, int row);
147 extern FCELL N_get_array_2d_f_value(N_array_2d * array2d, int col, int row);
148 extern DCELL N_get_array_2d_d_value(N_array_2d * array2d, int col, int row);
149 extern void N_put_array_2d_value(N_array_2d * array2d, int col, int row, char *value);
150 extern void N_put_array_2d_c_value(N_array_2d * array2d, int col, int row, CELL value);
151 extern void N_put_array_2d_f_value(N_array_2d * array2d, int col, int row, FCELL value);
152 extern void N_put_array_2d_d_value(N_array_2d * array2d, int col, int row, DCELL value);
153 extern int N_is_array_2d_value_null(N_array_2d * array2d, int col, int row);
154 extern void N_put_array_2d_value_null(N_array_2d * array2d, int col, int row);
155 extern void N_print_array_2d(N_array_2d * data);
156 extern void N_print_array_2d_info(N_array_2d * data);
157 extern void N_copy_array_2d(N_array_2d * source, N_array_2d * target);
158 extern double N_norm_array_2d(N_array_2d * array1, N_array_2d * array2, int type);
159 extern N_array_2d *N_math_array_2d(N_array_2d * array1, N_array_2d * array2, N_array_2d * result, int type);
160 extern int N_convert_array_2d_null_to_zero(N_array_2d * a);
161 extern N_array_2d *N_read_rast_to_array_2d(char *name, N_array_2d * array);
162 extern void N_write_array_2d_to_rast(N_array_2d * array, char *name);
163 extern void N_calc_array_2d_stats(N_array_2d * a, double *min, double *max, double *sum, int *nonzero, int withoffset);
164 
165 typedef struct
166 {
167     int type;			/* which raster type FCELL_TYPE, DCELL_TYPE */
168     int rows, cols, depths;
169     int rows_intern, cols_intern, depths_intern;
170     int offset;			/*number of cols/rows/depths offset at each boundary */
171     float *fcell_array;		/*The data is stored in an one dimensional array internally */
172     double *dcell_array;	/*The data is stored in an one dimensional array internally */
173 } N_array_3d;
174 
175 extern N_array_3d *N_alloc_array_3d(int cols, int rows, int depths, int offset, int type);
176 extern void N_free_array_3d(N_array_3d * data_array);
177 extern int N_get_array_3d_type(N_array_3d * array3d);
178 extern void N_get_array_3d_value(N_array_3d * array3d, int col, int row, int depth, void *value);
179 extern float N_get_array_3d_f_value(N_array_3d * array3d, int col, int row, int depth);
180 extern double N_get_array_3d_d_value(N_array_3d * array3d, int col, int row, int depth);
181 extern void N_put_array_3d_value(N_array_3d * array3d, int col, int row, int depth, char *value);
182 extern void N_put_array_3d_f_value(N_array_3d * array3d, int col, int row, int depth, float value);
183 extern void N_put_array_3d_d_value(N_array_3d * array3d, int col, int row, int depth, double value);
184 extern int N_is_array_3d_value_null(N_array_3d * array3d, int col, int row, int depth);
185 extern void N_put_array_3d_value_null(N_array_3d * array3d, int col, int row, int depth);
186 extern void N_print_array_3d(N_array_3d * data);
187 extern void N_print_array_3d_info(N_array_3d * data);
188 extern void N_copy_array_3d(N_array_3d * source, N_array_3d * target);
189 extern double N_norm_array_3d(N_array_3d * array1, N_array_3d * array2, int type);
190 extern N_array_3d *N_math_array_3d(N_array_3d * array1, N_array_3d * array2, N_array_3d * result, int type);
191 extern int N_convert_array_3d_null_to_zero(N_array_3d * a);
192 extern N_array_3d *N_read_rast3d_to_array_3d(char *name, N_array_3d * array, int mask);
193 extern void N_write_array_3d_to_rast3d(N_array_3d * array, char *name, int mask);
194 extern void N_calc_array_3d_stats(N_array_3d * a, double *min, double *max, double *sum, int *nonzero, int withoffset);
195 
196 /* *************************************************************** */
197 /* *************** MATRIX ASSEMBLING METHODS ********************* */
198 /* *************************************************************** */
199 /*!
200  * \brief Matrix entries for a mass balance 5/7/9 star system
201  *
202  * Matrix entries for the mass balance of a 5 star system
203  *
204  * The entries are center, east, west, north, south and the
205  * right side vector b of Ax = b. This system is typically used in 2d.
206 
207  \verbatim
208  N
209  |
210  W-- C --E
211  |
212  S
213  \endverbatim
214 
215  * Matrix entries for the mass balance of a 7 star system
216  *
217  * The entries are center, east, west, north, south, top, bottom and the
218  * right side vector b of Ax = b. This system is typically used in 3d.
219 
220  \verbatim
221  T N
222  |/
223  W-- C --E
224  /|
225  S B
226  \endverbatim
227 
228  * Matrix entries for the mass balance of a 9 star system
229  *
230  * The entries are center, east, west, north, south, north-east, south-east,
231  * north-wast, south-west and the
232  * right side vector b of Ax = b. This system is typically used in 2d.
233 
234  \verbatim
235  NW  N  NE
236  \ | /
237  W-- C --E
238  / | \
239  SW  S  SE
240  \endverbatim
241 
242  * Matrix entries for the mass balance of a 27 star system
243  *
244  * The entries are center, east, west, north, south, north-east, south-east,
245  * north-wast, south-west, same for top and bottom and the
246  * right side vector b of Ax = b. This system is typically used in 2d.
247 
248  \verbatim
249  top:
250  NW_T N_Z NE_T
251  \ | /
252  W_T-- T --E_T
253  / | \
254  SW_T S_T SE_T
255 
256  center:
257  NW  N  NE
258  \ | /
259  W-- C --E
260  / | \
261  SW  S  SE
262 
263  bottom:
264  NW_B N_B NE_B
265  \ | /
266  W_B-- B --E_B
267  / | \
268  SW_B S_B SE_B
269  \endverbatim
270 
271  */
272 typedef struct
273 {
274     int type;
275     int count;
276     double C, W, E, N, S, NE, NW, SE, SW, V;
277     /*top part */
278     double T, W_T, E_T, N_T, S_T, NE_T, NW_T, SE_T, SW_T;
279     /*bottom part */
280     double B, W_B, E_B, N_B, S_B, NE_B, NW_B, SE_B, SW_B;
281 } N_data_star;
282 
283 /*!
284  * \brief callback structure for 3d matrix assembling
285  * */
286 typedef struct
287 {
288     N_data_star *(*callback) ();
289 } N_les_callback_3d;
290 
291 /*!
292  * \brief callback structure for 2d matrix assembling
293  * */
294 typedef struct
295 {
296     N_data_star *(*callback) ();
297 } N_les_callback_2d;
298 
299 
300 extern void N_set_les_callback_3d_func(N_les_callback_3d * data, N_data_star * (*callback_func_3d) ());
301 extern void N_set_les_callback_2d_func(N_les_callback_2d * data, N_data_star * (*callback_func_2d) ());
302 extern N_les_callback_3d *N_alloc_les_callback_3d(void);
303 extern N_les_callback_2d *N_alloc_les_callback_2d(void);
304 extern N_data_star *N_alloc_5star(void);
305 extern N_data_star *N_alloc_7star(void);
306 extern N_data_star *N_alloc_9star(void);
307 extern N_data_star *N_alloc_27star(void);
308 extern N_data_star *N_create_5star(double C, double W, double E, double N,
309 				   double S, double V);
310 extern N_data_star *N_create_7star(double C, double W, double E, double N,
311 				   double S, double T, double B, double V);
312 extern N_data_star *N_create_9star(double C, double W, double E, double N,
313 				   double S, double NW, double SW, double NE,
314 				   double SE, double V);
315 extern N_data_star *N_create_27star(double C, double W, double E, double N,
316 				    double S, double NW, double SW, double NE,
317 				    double SE, double T, double W_T,
318 				    double E_T, double N_T, double S_T,
319 				    double NW_T, double SW_T, double NE_T,
320 				    double SE_T, double B, double W_B,
321 				    double E_B, double N_B, double S_B,
322 				    double NW_B, double SW_B, double NE_B,
323 				    double SE_B, double V);
324 extern N_data_star *N_callback_template_3d(void *data, N_geom_data * geom, int col, int row, int depth);
325 extern N_data_star *N_callback_template_2d(void *data, N_geom_data * geom, int col, int row);
326 extern N_les *N_assemble_les_3d(int les_type, N_geom_data * geom, N_array_3d * status, N_array_3d * start_val, void *data, N_les_callback_3d * callback);
327 extern N_les *N_assemble_les_3d_active(int les_type, N_geom_data * geom, N_array_3d * status, N_array_3d * start_val, void *data, N_les_callback_3d * callback);
328 extern N_les *N_assemble_les_3d_dirichlet(int les_type, N_geom_data * geom, N_array_3d * status, N_array_3d * start_val, void *data, N_les_callback_3d * callback);
329 extern N_les *N_assemble_les_3d_param(int les_type, N_geom_data * geom, N_array_3d * status, N_array_3d * start_val, void *data, N_les_callback_3d * callback, int cell_type);
330 extern N_les *N_assemble_les_2d(int les_type, N_geom_data * geom, N_array_2d * status, N_array_2d * start_val, void *data, N_les_callback_2d * callback);
331 extern N_les *N_assemble_les_2d_active(int les_type, N_geom_data * geom, N_array_2d * status, N_array_2d * start_val, void *data, N_les_callback_2d * callback);
332 extern N_les *N_assemble_les_2d_dirichlet(int les_type, N_geom_data * geom, N_array_2d * status, N_array_2d * start_val, void *data, N_les_callback_2d * callback);
333 extern N_les *N_assemble_les_2d_param(int les_type, N_geom_data * geom, N_array_2d * status, N_array_2d * start_val, void *data, N_les_callback_2d * callback, int cell_Type);
334 extern int N_les_pivot_create(N_les * les);
335 int N_les_integrate_dirichlet_2d(N_les * les, N_geom_data * geom, N_array_2d * status, N_array_2d * start_val);
336 int N_les_integrate_dirichlet_3d(N_les * les, N_geom_data * geom, N_array_3d * status, N_array_3d * start_val);
337 
338 /* *************************************************************** */
339 /* *************** GPDE STANDARD OPTIONS ************************* */
340 /* *************************************************************** */
341 
342 /*! \brief Standard options of the gpde library
343  * */
344 typedef enum
345 {
346     N_OPT_SOLVER_SYMM,		/*! solver for symmetric, positive definite linear equation systems */
347     N_OPT_SOLVER_UNSYMM,	/*! solver for unsymmetric linear equation systems */
348     N_OPT_MAX_ITERATIONS,	/*! Maximum number of iteration used to solver the linear equation system */
349     N_OPT_ITERATION_ERROR,	/*! Error break criteria for the iterative solver (jacobi, sor, cg or bicgstab) */
350     N_OPT_SOR_VALUE,		/*! The relaxation parameter used by the jacobi and sor solver for speedup or stabilizing */
351     N_OPT_CALC_TIME		/*! The calculation time in seconds */
352 } N_STD_OPT;
353 
354 extern struct Option *N_define_standard_option(int opt);
355 
356 /* *************************************************************** */
357 /* *************** GPDE MATHEMATICAL TOOLS *********************** */
358 /* *************************************************************** */
359 
360 extern double N_calc_arith_mean(double a, double b);
361 extern double N_calc_arith_mean_n(double *a, int size);
362 extern double N_calc_geom_mean(double a, double b);
363 extern double N_calc_geom_mean_n(double *a, int size);
364 extern double N_calc_harmonic_mean(double a, double b);
365 extern double N_calc_harmonic_mean_n(double *a, int size);
366 extern double N_calc_quad_mean(double a, double b);
367 extern double N_calc_quad_mean_n(double *a, int size);
368 
369 /* *************************************************************** */
370 /* *************** UPWIND STABILIZATION ALGORITHMS *************** */
371 /* *************************************************************** */
372 
373 extern double N_full_upwinding(double sprod, double distance, double D);
374 extern double N_exp_upwinding(double sprod, double distance, double D);
375 
376 
377 /* *************************************************************** */
378 /* *************** METHODS FOR GRADIENT CALCULATION ************** */
379 /* *************************************************************** */
380 /*!
381    \verbatim
382 
383    ______________
384    |    |    |    |
385    |    |    |    |
386    |----|-NC-|----|
387    |    |    |    |
388    |   WC    EC   |
389    |    |    |    |
390    |----|-SC-|----|
391    |    |    |    |
392    |____|____|____|
393 
394 
395    |  /
396    TC NC
397    |/
398    --WC-----EC--
399    /|
400    SC BC
401    /  |
402 
403    \endverbatim
404 
405  */
406 
407 /*! \brief Gradient between the cells in X and Y direction */
408 typedef struct
409 {
410 
411     double NC, SC, WC, EC;
412 
413 } N_gradient_2d;
414 
415 /*! \brief Gradient between the cells in X, Y and Z direction */
416 typedef struct
417 {
418 
419     double NC, SC, WC, EC, TC, BC;
420 
421 } N_gradient_3d;
422 
423 
424 /*!
425    \verbatim
426 
427    Gradient in X direction between the cell neighbours
428    ____ ____ ____
429    |    |    |    |
430    |   NWN  NEN   |
431    |____|____|____|
432    |    |    |    |
433    |   WN    EN   |
434    |____|____|____|
435    |    |    |    |
436    |   SWS  SES   |
437    |____|____|____|
438 
439    Gradient in Y direction between the cell neighbours
440    ______________
441    |    |    |    |
442    |    |    |    |
443    |NWW-|-NC-|-NEE|
444    |    |    |    |
445    |    |    |    |
446    |SWW-|-SC-|-SEE|
447    |    |    |    |
448    |____|____|____|
449 
450    Gradient in Z direction between the cell neighbours
451    /______________/
452    /|    |    |    |
453    | NWZ| NZ | NEZ|
454    |____|____|____|
455    /|    |    |    |
456    | WZ | CZ | EZ |
457    |____|____|____|
458    /|    |    |    |
459    | SWZ| SZ | SEZ|
460    |____|____|____|
461    /____/____/____/
462 
463 
464    \endverbatim
465  */
466 
467 /*! \brief Gradient between the cell neighbours in X direction */
468 typedef struct
469 {
470 
471     double NWN, NEN, WC, EC, SWS, SES;
472 
473 } N_gradient_neighbours_x;
474 
475 /*! \brief Gradient between the cell neighbours in Y direction */
476 typedef struct
477 {
478 
479     double NWW, NEE, NC, SC, SWW, SEE;
480 
481 } N_gradient_neighbours_y;
482 
483 /*! \brief Gradient between the cell neighbours in Z direction */
484 typedef struct
485 {
486 
487     double NWZ, NZ, NEZ, WZ, CZ, EZ, SWZ, SZ, SEZ;
488 
489 } N_gradient_neighbours_z;
490 
491 /*! \brief Gradient between the cell neighbours in X and Y direction */
492 typedef struct
493 {
494 
495     N_gradient_neighbours_x *x;
496     N_gradient_neighbours_y *y;
497 
498 } N_gradient_neighbours_2d;
499 
500 
501 /*! \brief Gradient between the cell neighbours in X, Y and Z direction */
502 typedef struct
503 {
504 
505     N_gradient_neighbours_x *xt;	/*top values */
506     N_gradient_neighbours_x *xc;	/*center values */
507     N_gradient_neighbours_x *xb;	/*bottom values */
508 
509     N_gradient_neighbours_y *yt;	/*top values */
510     N_gradient_neighbours_y *yc;	/*center values */
511     N_gradient_neighbours_y *yb;	/*bottom values */
512 
513     N_gradient_neighbours_z *zt;	/*top-center values */
514     N_gradient_neighbours_z *zb;	/*bottom-center values */
515 
516 } N_gradient_neighbours_3d;
517 
518 
519 /*! Two dimensional gradient field */
520 typedef struct
521 {
522 
523     N_array_2d *x_array;
524     N_array_2d *y_array;
525     int cols, rows;
526     double min, max, mean, sum;
527     int nonull;
528 
529 } N_gradient_field_2d;
530 
531 /*! Three dimensional gradient field */
532 typedef struct
533 {
534 
535     N_array_3d *x_array;
536     N_array_3d *y_array;
537     N_array_3d *z_array;
538     int cols, rows, depths;
539     double min, max, mean, sum;
540     int nonull;
541 
542 } N_gradient_field_3d;
543 
544 
545 extern N_gradient_2d *N_alloc_gradient_2d(void);
546 extern void N_free_gradient_2d(N_gradient_2d * grad);
547 extern N_gradient_2d *N_create_gradient_2d(double NC, double SC, double WC, double EC);
548 extern int N_copy_gradient_2d(N_gradient_2d * source, N_gradient_2d * target);
549 extern N_gradient_2d *N_get_gradient_2d(N_gradient_field_2d * field, N_gradient_2d * gradient, int col, int row);
550 extern N_gradient_3d *N_alloc_gradient_3d(void);
551 extern void N_free_gradient_3d(N_gradient_3d * grad);
552 extern N_gradient_3d *N_create_gradient_3d(double NC, double SC, double WC, double EC, double TC, double BC);
553 extern int N_copy_gradient_3d(N_gradient_3d * source, N_gradient_3d * target);
554 extern N_gradient_3d *N_get_gradient_3d(N_gradient_field_3d * field, N_gradient_3d * gradient, int col, int row, int depth);
555 extern N_gradient_neighbours_x *N_alloc_gradient_neighbours_x(void);
556 extern void N_free_gradient_neighbours_x(N_gradient_neighbours_x * grad);
557 extern N_gradient_neighbours_x *N_create_gradient_neighbours_x(double NWN,
558 							       double NEN,
559 							       double WC,
560 							       double EC,
561 							       double SWS,
562 							       double SES);
563 extern int N_copy_gradient_neighbours_x(N_gradient_neighbours_x * source, N_gradient_neighbours_x * target);
564 extern N_gradient_neighbours_y *N_alloc_gradient_neighbours_y(void);
565 extern void N_free_gradient_neighbours_y(N_gradient_neighbours_y * grad);
566 extern N_gradient_neighbours_y *N_create_gradient_neighbours_y(double NWW,
567 							       double NEE,
568 							       double NC,
569 							       double SC,
570 							       double SWW,
571 							       double SEE);
572 extern int N_copy_gradient_neighbours_y(N_gradient_neighbours_y * source, N_gradient_neighbours_y * target);
573 extern N_gradient_neighbours_z *N_alloc_gradient_neighbours_z(void);
574 extern void N_free_gradient_neighbours_z(N_gradient_neighbours_z * grad);
575 extern N_gradient_neighbours_z *N_create_gradient_neighbours_z(double NWZ,
576 							       double NZ,
577 							       double NEZ,
578 							       double WZ,
579 							       double CZ,
580 							       double EZ,
581 							       double SWZ,
582 							       double SZ,
583 							       double SEZ);
584 extern int N_copy_gradient_neighbours_z(N_gradient_neighbours_z * source, N_gradient_neighbours_z * target);
585 extern N_gradient_neighbours_2d *N_alloc_gradient_neighbours_2d(void);
586 extern void N_free_gradient_neighbours_2d(N_gradient_neighbours_2d * grad);
587 extern N_gradient_neighbours_2d * N_create_gradient_neighbours_2d(N_gradient_neighbours_x * x, N_gradient_neighbours_y * y);
588 extern int N_copy_gradient_neighbours_2d(N_gradient_neighbours_2d * source, N_gradient_neighbours_2d * target);
589 extern N_gradient_neighbours_2d * N_get_gradient_neighbours_2d(N_gradient_field_2d * field, N_gradient_neighbours_2d * gradient, int col, int row);
590 extern N_gradient_neighbours_3d *N_alloc_gradient_neighbours_3d(void);
591 extern void N_free_gradient_neighbours_3d(N_gradient_neighbours_3d * grad);
592 extern N_gradient_neighbours_3d
593     * N_create_gradient_neighbours_3d(N_gradient_neighbours_x * xt,
594 				      N_gradient_neighbours_x * xc,
595 				      N_gradient_neighbours_x * xb,
596 				      N_gradient_neighbours_y * yt,
597 				      N_gradient_neighbours_y * yc,
598 				      N_gradient_neighbours_y * yb,
599 				      N_gradient_neighbours_z * zt,
600 				      N_gradient_neighbours_z * zb);
601 extern int N_copy_gradient_neighbours_3d(N_gradient_neighbours_3d * source, N_gradient_neighbours_3d * target);
602 extern void N_print_gradient_field_2d_info(N_gradient_field_2d * field);
603 extern void N_calc_gradient_field_2d_stats(N_gradient_field_2d * field);
604 extern N_gradient_field_2d *N_alloc_gradient_field_2d(int cols, int rows);
605 extern void N_free_gradient_field_2d(N_gradient_field_2d * field);
606 extern int N_copy_gradient_field_2d(N_gradient_field_2d * source, N_gradient_field_2d * target);
607 extern N_gradient_field_2d *N_compute_gradient_field_2d(N_array_2d * pot,
608 							N_array_2d * weight_x,
609 							N_array_2d * weight_y,
610 							N_geom_data * geom,
611 							N_gradient_field_2d *
612 							gradfield);
613 extern void N_compute_gradient_field_components_2d(N_gradient_field_2d * field, N_array_2d * x_comp, N_array_2d * y_comp);
614 extern void N_print_gradient_field_3d_info(N_gradient_field_3d * field);
615 extern void N_calc_gradient_field_3d_stats(N_gradient_field_3d * field);
616 extern N_gradient_field_3d *N_alloc_gradient_field_3d(int cols, int rows, int depths);
617 extern void N_free_gradient_field_3d(N_gradient_field_3d * field);
618 extern int N_copy_gradient_field_3d(N_gradient_field_3d * source, N_gradient_field_3d * target);
619 extern N_gradient_field_3d *N_compute_gradient_field_3d(N_array_3d * pot,
620 							N_array_3d * weight_x,
621 							N_array_3d * weight_y,
622 							N_array_3d * weight_z,
623 							N_geom_data * geom,
624 							N_gradient_field_3d *
625 							gradfield);
626 extern void N_compute_gradient_field_components_3d(N_gradient_field_3d * field, N_array_3d * x_comp, N_array_3d * y_comp, N_array_3d * z_comp);
627 
628 #endif
629