/*
 * cactus_total.c
 * Copyright (c) 2004-2008 Royal Observatory of Belgium
 * Author: Bogdan Nicula
 *
 * v1.1 Feb 2008 - DLM, bounds checking
 * v1.0 Aug 2004 - initial implemenatation
 */

#include <stdio.h>
#include "idl_export.h"

#define _MIN(a,b) (((a) < (b)) ? (a) : (b))
#define _MAX(a,b) (((a) > (b)) ? (a) : (b))

/*
	0 velmap
	1 myxt
	2 vectors
	3 min
	4 max
*/

static void the_total(int argc, IDL_VPTR argv[], char *argk)
{
	int i, j, k, d0, d1, d2, st, min, max;
	unsigned int inlen, idx;
	const unsigned int *vec, *ovec;
	const float *in;
	float *out;
	double s;

	IDL_VPTR velmap = argv[0];
	IDL_VPTR myxt = argv[1];
	IDL_VPTR vectors = argv[2];

	IDL_EXCLUDE_CONST(velmap);
	IDL_ENSURE_ARRAY(velmap);
	IDL_ENSURE_ARRAY(myxt);
	IDL_ENSURE_ARRAY(vectors);

	if (velmap->value.arr->n_dim != 2)
		IDL_Message(IDL_M_NAMED_GENERIC, IDL_MSG_LONGJMP,
					"first argument (velmap) must be a 2-D array");
	if (myxt->value.arr->n_dim != 2)
		IDL_Message(IDL_M_NAMED_GENERIC, IDL_MSG_LONGJMP,
					"second argument (myxt) must be a 2-D array");
	if (vectors->value.arr->n_dim != 3)
		IDL_Message(IDL_M_NAMED_GENERIC, IDL_MSG_LONGJMP,
					"third argument (vectors) must be a 3-D array");
/* */
	if (velmap->type != IDL_TYP_FLOAT)
		IDL_Message(IDL_M_NAMED_GENERIC, IDL_MSG_LONGJMP,
					"first argument (velmap) must be an array of floats");
	if (myxt->type != IDL_TYP_FLOAT)
		IDL_Message(IDL_M_NAMED_GENERIC, IDL_MSG_LONGJMP,
					"second argument (myxt) must be an array of floats");
/* */

	if (vectors->type != IDL_TYP_ULONG)
		IDL_Message(IDL_M_NAMED_GENERIC, IDL_MSG_LONGJMP,
					"third argument (vectors) must be an array of unsigned longs");

	d0 = vectors->value.arr->dim[0];
	d1 = vectors->value.arr->dim[1];
	d2 = vectors->value.arr->dim[2];
	st = d0 * d1;

	min = IDL_LongScalar(argv[3]);
	max = IDL_LongScalar(argv[4]);
	min = _MAX(min, 0);
	max = _MIN(max, d2 - 1);
	if (min > max)
		min = max;

	if (velmap->value.arr->dim[0] != d2)
		IDL_Message(IDL_M_NAMED_GENERIC, IDL_MSG_LONGJMP,
					"velmap dim[0] != vectors dim[2]");
	if (velmap->value.arr->dim[1] != d1)
		IDL_Message(IDL_M_NAMED_GENERIC, IDL_MSG_LONGJMP,
					"velmap dim[1] != vectors dim[1]");

	inlen = myxt->value.arr->dim[0] * myxt->value.arr->dim[1] - 1;

	out = (float *) velmap->value.arr->data;
	in = (const float *) myxt->value.arr->data;
	vec = (const unsigned int *) vectors->value.arr->data;

	out += min;
	vec += min * st;

	for (k = min; k <= max; ++k) {
		ovec = vec;
		for (j = 0; j < d1; ++j) {
			s = 0.;
			for (i = 0; i < d0; ++i) {
				idx = *ovec++;
				if (idx > inlen)
					IDL_Message(IDL_M_NAMED_GENERIC, IDL_MSG_LONGJMP,
								"index out of bounds");
				s += in[idx];
			}
			out[j * d2] = s;
		}
		vec += st;
		++out;
	}
}

int IDL_Load(void)
{
	static IDL_SYSFUN_DEF2 proc[] =
		{ {(IDL_SYSRTN_UNION) the_total, "CACTUS_TOTAL", 5, 5, 0, NULL}, };
	return IDL_SysRtnAdd(proc, FALSE, IDL_CARRAY_ELTS(proc));
}