static PyObject *
symbols(struct PyLibPlot *self, PyObject *args)
{
	PyObject *ox, *oy;
	PyObject *x, *y;
	double d0;
	int i0;
    npy_intp i, n;

	if ( !PyArg_ParseTuple( args, "OOid", &ox, &oy, &i0, &d0 ) )
		return NULL;

	x = PyArray_ContiguousFromAny( ox, NPY_DOUBLE, 1, 1 );
	y = PyArray_ContiguousFromAny( oy, NPY_DOUBLE, 1, 1 );

	if ( x == NULL || y == NULL )
		goto quit;

	n = BGL_MIN( PyArray_SIZE(x), PyArray_SIZE(y) );

	_symbol_begin( self->pl, i0, d0 );

	for ( i = 0; i < n; i++ )
		_symbol_draw( self->pl, BGL_DArray1(x,i), BGL_DArray1(y,i), i0, d0 );

	_symbol_end( self->pl, i0, d0 );

quit:
	Py_XDECREF(x);
	Py_XDECREF(y);
    Py_RETURN_NONE;
}

// create a new NDF (simple) structure
static PyObject*
pyndf_new(NDF *self, PyObject *args)
{
	// use ultracam defaults
	const char *ftype = "_REAL";
	int ndim;
	PyObject* lb;
	PyObject* ub;
	if(!PyArg_ParseTuple(args, "siOO:pyndf_new", &ftype, &ndim, &lb, &ub))
		return NULL;
	if(ndim < 0 || ndim > 7)
		return NULL;
	// TODO: check for ftype here
	int status = SAI__OK;
	PyArrayObject* lower = (PyArrayObject*) PyArray_FROM_OTF(lb, NPY_UINT, NPY_IN_ARRAY | NPY_FORCECAST);
	PyArrayObject* upper = (PyArrayObject*) PyArray_FROM_OTF(ub, NPY_UINT, NPY_IN_ARRAY | NPY_FORCECAST);
	if (!lower || !upper)
		return NULL;
	if(PyArray_SIZE(lower) != ndim || PyArray_SIZE(upper) != ndim)
		return NULL;
        errBegin(&status);
	ndfNew(ftype,ndim,(int*)PyArray_DATA(lower),(int*)PyArray_DATA(upper),&self->_place,&self->_ndfid,&status); // placeholder annulled by this routine
	Py_DECREF(lower);
	Py_DECREF(upper);
	if (raiseNDFException(&status))
		return NULL;
	Py_RETURN_NONE;
}

static PyObject *
curve(struct PyLibPlot *self, PyObject *args)
{
	PyObject *ox, *oy;
	PyObject *x, *y;
	npy_intp i, n;

	if ( !PyArg_ParseTuple( args, "OO", &ox, &oy ) )
		return NULL;

	x = PyArray_ContiguousFromAny( ox, NPY_DOUBLE, 1, 1 );
	y = PyArray_ContiguousFromAny( oy, NPY_DOUBLE, 1, 1 );

	if ( x == NULL || y == NULL )
		goto quit;

	n = BGL_MIN( PyArray_SIZE(x), PyArray_SIZE(y) );
	if ( n <= 0 )
		goto quit;

	pl_fmove_r( self->pl, BGL_DArray1(x,0), BGL_DArray1(y,0) );
	for ( i = 1; i < n; i++ )
		pl_fcont_r( self->pl, BGL_DArray1(x,i), BGL_DArray1(y,i) );
	pl_endpath_r( self->pl );

quit:
	Py_XDECREF(x);
	Py_XDECREF(y);
    Py_RETURN_NONE;
}

static PyObject* write_coords(PyObject* self, PyObject* args)
{
    oskar_MeasurementSet* h = 0;
    PyObject *capsule = 0;
    PyObject *obj[] = {0, 0, 0};
    PyArrayObject *uu = 0, *vv = 0, *ww = 0;
    int start_row = 0, num_baselines = 0;
    double exposure_sec = 0.0, interval_sec = 0.0, time_stamp = 0.0;
    if (!PyArg_ParseTuple(args, "OiiOOOddd", &capsule,
            &start_row, &num_baselines, &obj[0], &obj[1], &obj[2],
            &exposure_sec, &interval_sec, &time_stamp))
        return 0;
    if (!(h = (oskar_MeasurementSet*) get_handle(capsule, name))) return 0;

    /* Make sure input objects are arrays. Convert if required. */
    uu = (PyArrayObject*) PyArray_FROM_OF(obj[0], NPY_ARRAY_IN_ARRAY);
    vv = (PyArrayObject*) PyArray_FROM_OF(obj[1], NPY_ARRAY_IN_ARRAY);
    ww = (PyArrayObject*) PyArray_FROM_OF(obj[2], NPY_ARRAY_IN_ARRAY);
    if (!uu || !vv || !ww)
        goto fail;

    /* Check dimensions. */
    if (num_baselines != (int) PyArray_SIZE(uu) ||
            num_baselines != (int) PyArray_SIZE(vv) ||
            num_baselines != (int) PyArray_SIZE(ww))
    {
        PyErr_SetString(PyExc_RuntimeError, "Input data dimension mismatch.");
        goto fail;
    }

    /* Write the coordinates. */
    Py_BEGIN_ALLOW_THREADS
    if (PyArray_TYPE(uu) == NPY_DOUBLE)
        oskar_ms_write_coords_d(h, start_row, num_baselines,
                (const double*)PyArray_DATA(uu),
                (const double*)PyArray_DATA(vv),
                (const double*)PyArray_DATA(ww),
                exposure_sec, interval_sec, time_stamp);
    else
        oskar_ms_write_coords_f(h, start_row, num_baselines,
                (const float*)PyArray_DATA(uu),
                (const float*)PyArray_DATA(vv),
                (const float*)PyArray_DATA(ww),
                exposure_sec, interval_sec, time_stamp);
    Py_END_ALLOW_THREADS

    Py_XDECREF(uu);
    Py_XDECREF(vv);
    Py_XDECREF(ww);
    return Py_BuildValue("");

fail:
    Py_XDECREF(uu);
    Py_XDECREF(vv);
    Py_XDECREF(ww);
    return 0;
}

static PyArrayObject *mx2numeric(const mxArray *pArray)
{
  //current function returns PyArrayObject in c order currently
  mwSize nd;
  npy_intp  pydims[NPY_MAXDIMS];
  PyArrayObject *lRetval = NULL,*t=NULL;
  const double *lPR;
  const double *lPI;
  pyassert(PyArray_API,
           "Unable to perform this function without NumPy installed");

  nd = mxGetNumberOfDimensions(pArray);
  {
    const mwSize *dims;
    dims = mxGetDimensions(pArray);
    for (mwSize i=0; i != nd; i++){
        pydims[i] = static_cast<npy_intp>(dims[i]);
    }
  }
 //this function creates a fortran array
  t = (PyArrayObject *)
    PyArray_New(&PyArray_Type,static_cast<npy_intp>(nd), pydims,
                mxIsComplex(pArray) ? PyArray_CDOUBLE : PyArray_DOUBLE,
                NULL, // strides
                NULL, // data
                0,    //(ignored itemsize),
                NPY_F_CONTIGUOUS, 
                NULL); //  obj
  if (t == NULL) return NULL;
  
  lPR  = mxGetPr(pArray);
  if (mxIsComplex(pArray)) {
    double *lDst = (double *)PyArray_DATA(t);
    // AWMS unsigned int almost certainly can overflow on some platforms!
    npy_intp numberOfElements = PyArray_SIZE(t);
    lPI = mxGetPi(pArray);
    for (mwIndex i = 0; i != numberOfElements; i++) {
      *lDst++ = *lPR++;
      *lDst++ = *lPI++;
    }
  }
  else {
    double *lDst = (double *)PyArray_DATA(t);
    npy_intp numberOfElements = PyArray_SIZE(t);
    for (mwIndex i = 0; i != numberOfElements; i++) {
      *lDst++ = *lPR++;
    }
  }
  
  lRetval = (PyArrayObject *)PyArray_FromArray(t,NULL,NPY_C_CONTIGUOUS|NPY_ALIGNED|NPY_WRITEABLE);
  Py_DECREF(t);
  
  return lRetval;
  error_return:
  return NULL;
}

static PyObject *
clipped_colored_symbols(struct PyLibPlot *self, PyObject *args)
{
	PyObject *ox, *oy, *oc;
	PyObject *x, *y, *c;
	double xmin, xmax, ymin, ymax;
	double d0;
	int i0;
    npy_intp i, n;
	double px, py;
	int r, g, b;

	if ( !PyArg_ParseTuple( args, "OOOiddddd", &ox, &oy, &oc,
			&i0, &d0, &xmin, &xmax, &ymin, &ymax ) )
		return NULL;

	x = PyArray_ContiguousFromAny( ox, NPY_DOUBLE, 1, 1 );
	y = PyArray_ContiguousFromAny( oy, NPY_DOUBLE, 1, 1 );
	c = PyArray_ContiguousFromAny( oc, NPY_DOUBLE, 2, 2 );

	if ( x == NULL || y == NULL || c == NULL )
		goto quit;

	n = BGL_MIN( PyArray_SIZE(x), PyArray_SIZE(y) );
	n = BGL_MIN( n, PyArray_SIZE(c) );
	
	_symbol_begin( self->pl, i0, d0 );

	for ( i = 0; i < n; i++ )
	{
		px = BGL_DArray1(x,i);
		py = BGL_DArray1(y,i);

		if ( px >= xmin && px <= xmax &&
		     py >= ymin && py <= ymax ) {
			r = (int) floor( BGL_DArray2(c,i,0)*65535 );
			g = (int) floor( BGL_DArray2(c,i,1)*65535 );
			b = (int) floor( BGL_DArray2(c,i,2)*65535 );
			pl_fillcolor_r( self->pl, r, g, b );
			pl_pencolor_r(  self->pl, r, g, b );
		  
			_symbol_draw( self->pl, px, py, i0, d0 );
		}
	}

	_symbol_end( self->pl, i0, d0 );

quit:
	Py_XDECREF(x);
	Py_XDECREF(y);
	Py_XDECREF(c);
    Py_RETURN_NONE;
}

static PyObject *inverse_subtraction_inplace(PyObject *self, PyObject *args)
{
  PyObject* a = NULL;
  PyObject* b = NULL;
  npy_intp sz = 0, i;
  npy_complex64* tmp_sp = NULL;
  npy_complex64* a_data_sp = NULL;
  npy_float32* b_data_sp = NULL;
  npy_complex128* tmp_dp = NULL;
  npy_complex128* a_data_dp = NULL;
  npy_float64* b_data_dp = NULL;
  double c;
  if (!PyArg_ParseTuple(args, "OOd", &a, &b, &c))
    return NULL;
  if (!(PyArray_Check(a) && PyArray_Check(b)))
    {
      PyErr_SetString(PyExc_TypeError,"arguments must be array objects");
      return NULL;
    }
  sz = PyArray_SIZE(a);

  if (sz != PyArray_SIZE(b))
    {
      PyErr_SetString(PyExc_TypeError,"argument sizes must be equal");
      return NULL;
    }
  if ((PyArray_TYPE(a) == PyArray_COMPLEX64) && (PyArray_TYPE(b) == PyArray_FLOAT32))
    {
      a_data_sp = (npy_complex64*)PyArray_DATA(a);
      b_data_sp = (npy_float32*)PyArray_DATA(b);
      for (i=0; i<sz; ++i)
	{
	  tmp_sp = a_data_sp + i;
	  tmp_sp->real = b_data_sp[i] - tmp_sp->real * c; 
	}
    }
  else if ((PyArray_TYPE(a) == PyArray_COMPLEX128) && (PyArray_TYPE(b) == PyArray_FLOAT64))
    {
      a_data_dp = (npy_complex128*)PyArray_DATA(a);
      b_data_dp = (npy_float64*)PyArray_DATA(b);
      for (i=0; i<sz; ++i)
	{
	  tmp_dp = a_data_dp + i;
	  tmp_dp->real = b_data_dp[i] - tmp_dp->real * c; 
	}
    }
  else
    {
      PyErr_SetString(PyExc_TypeError,"argument types must be complex64 and float32");
      return NULL;
    }
  return Py_BuildValue("");
}

static inline int not_vect(PyObject *vv)
{
	PyArrayObject *v = (PyArrayObject*)vv;

	if (PyArray_NDIM(v) == 2 &&
	    (PyArray_SIZE(v) == PyArray_DIM(v, 0) ||
	     PyArray_SIZE(v) == PyArray_DIM(v, 1)))
		return 0;
	else if (PyArray_NDIM(v) == 1)
		return 0;
	else
		return 1;
}

static PyObject * merge(PyObject * self,
     PyObject * args, PyObject * kwds) {

    static char * kwlist[] = {
        "data", "A", "B", "out", NULL
    };

    PyArrayObject * data, * A, * B, * out;
    if(!PyArg_ParseTupleAndKeywords(args, kwds, 
        "O!O!O!O!", kwlist,
        &PyArray_Type, &data, 
        &PyArray_Type, &A, 
        &PyArray_Type, &B, 
        &PyArray_Type, &out)) return NULL;

    npy_intp i = 0;

    int (*compare)(void *, void *, void*) = PyArray_DESCR(data)->f->compare;

    size_t sizeA = PyArray_SIZE(A);
    size_t sizeB = PyArray_SIZE(B);
    npy_intp * Aptr = PyArray_DATA(A);
    npy_intp * Aend = Aptr + sizeA;

    npy_intp * Bptr = PyArray_DATA(B);
    npy_intp * Bend = Bptr + sizeB;
    npy_intp * Optr = PyArray_DATA(out);

    #define VA PyArray_GETPTR1(data, *Aptr)
    #define VB PyArray_GETPTR1(data, (*Bptr) + sizeA)

    Py_BEGIN_ALLOW_THREADS
    while(Aptr < Aend|| Bptr < Bend) {
        while(Aptr < Aend && (Bptr == Bend || compare(VA, VB, data) <= 0)) {
            *Optr = *Aptr;
            Aptr++;
            Optr++;
            //printf("adding from A, i = %ld, k = %ld v = %ld\n", i, k, v);
        }
        while(Bptr < Bend && (Aptr == Aend || compare(VA, VB, data) >= 0)) {
            *Optr = *Bptr + sizeA;
            Bptr++;
            Optr++;
            //printf("adding from B, j = %ld, k = %ld v = %ld\n", j, k, v);
        }
    }
    Py_END_ALLOW_THREADS
    Py_INCREF(Py_None);
    return Py_None;
}

static PyObject* gse_checksum(PyObject *dummy, PyObject *args) {
        
    int checksum, length, i;
    PyObject *array = NULL;
    PyArrayObject *carray = NULL;
    int *data;

    if (!PyArg_ParseTuple(args, "O", &array )) {
        PyErr_SetString(GSEError, "usage checksum(array)" );
        return NULL;
    }
    if (!PyArray_Check(array)) {
        PyErr_SetString(GSEError, "Data must be given as NumPy array." );
        return NULL;
    }
    if (PyArray_TYPE(array) != NPY_INT32) {
        PyErr_SetString(GSEError, "Data must be 32-bit integers.");
        return NULL;
    }
    
    carray = PyArray_GETCONTIGUOUS((PyArrayObject*)array);
    length = PyArray_SIZE(carray);
    data = (int*)PyArray_DATA(carray);
    
    checksum = 0;
    for (i=0; i<length; i++) {
        checksum += data[i] % MODULUS;
        checksum %= MODULUS;
    }
    
    return Py_BuildValue("i", abs(checksum));
}

/* Allocate line buffer data */
int NI_AllocateLineBuffer(PyArrayObject* array, int axis, npy_intp size1,
        npy_intp size2, npy_intp *lines, npy_intp max_size, double **buffer)
{
    npy_intp line_size, max_lines;

    /* the number of lines of the array is an upper limit for the
         number of lines in the buffer: */
    max_lines = PyArray_SIZE(array);
    if (PyArray_NDIM(array) > 0 && PyArray_DIM(array, axis) > 0) {
        max_lines /= PyArray_DIM(array, axis);
    }
    /* calculate the space needed for one line, including space to
         support the boundary conditions: */
    line_size = sizeof(double) * (PyArray_DIM(array, axis) + size1 + size2);
    /* if *lines < 1, no number of lines is proposed, so we calculate it
         from the maximum size allowed: */
    if (*lines < 1) {
        *lines = line_size > 0 ? max_size / line_size : 0;
        if (*lines < 1)
            *lines = 1;
    }
    /* no need to allocate too many lines: */
    if (*lines > max_lines)
        *lines = max_lines;
    /* allocate data for the buffer: */
    *buffer = malloc(*lines * line_size);
    if (!*buffer) {
        PyErr_NoMemory();
        return 0;
    }
    return 1;
}

NRT_adapt_ndarray_from_python(PyObject *obj, arystruct_t* arystruct) {
    PyArrayObject *ndary;
    int i, ndim;
    npy_intp *p;
    void *data;

    if (!PyArray_Check(obj)) {
        return -1;
    }

    ndary = (PyArrayObject*)obj;
    ndim = PyArray_NDIM(ndary);
    data = PyArray_DATA(ndary);

    arystruct->meminfo = meminfo_new_from_pyobject((void*)data, obj);
    arystruct->data = data;
    arystruct->nitems = PyArray_SIZE(ndary);
    arystruct->itemsize = PyArray_ITEMSIZE(ndary);
    arystruct->parent = obj;
    p = arystruct->shape_and_strides;
    for (i = 0; i < ndim; i++, p++) {
        *p = PyArray_DIM(ndary, i);
    }
    for (i = 0; i < ndim; i++, p++) {
        *p = PyArray_STRIDE(ndary, i);
    }

    NRT_Debug(nrt_debug_print("NRT_adapt_ndarray_from_python %p\n",
                              arystruct->meminfo));
    return 0;
}

// create a new NDF extension
static PyObject*
pyndf_xnew(NDF *self, PyObject *args)
{
	int ndim = 0;
	const char *xname, *type;
	PyObject *dim;
	if(!PyArg_ParseTuple(args, "ss|iO:pyndf_xnew", &xname, &type, &ndim, &dim))
		return NULL;
	int status = SAI__OK;
	HDSLoc *loc = NULL;
	// perform checks if we're not making an extension header
	if(ndim != 0) {
		// check for HDS types
		if (!checkHDStype(type))
			return NULL;
		// need dims if it's not an ext
		if(ndim < 1 || dim == NULL)
			return NULL;
		PyArrayObject *npydim = (PyArrayObject*) PyArray_FROM_OTF(dim,NPY_INT,NPY_IN_ARRAY|NPY_FORCECAST);
		if (PyArray_SIZE(npydim) != ndim)
			return NULL;
                errBegin(&status);
		ndfXnew(self->_ndfid,xname,type,ndim,(int*)PyArray_DATA(npydim),&loc,&status);
		Py_DECREF(npydim);
	} else {
		// making an ext/struct
                errBegin(&status);
		ndfXnew(self->_ndfid,xname,type,0,0,&loc,&status);
	}
        if (raiseNDFException(&status)) return NULL;
	PyObject* pobj = NpyCapsule_FromVoidPtr(loc, PyDelLoc);
	return Py_BuildValue("O",pobj);
}

static PyObject *
array_float(PyArrayObject *v)
{
    PyObject *pv, *pv2;
    if (PyArray_SIZE(v) != 1) {
        PyErr_SetString(PyExc_TypeError, "only length-1 arrays can "\
                        "be converted to Python scalars");
        return NULL;
    }
    pv = v->descr->f->getitem(v->data, v);
    if (pv == NULL) {
        return NULL;
    }
    if (pv->ob_type->tp_as_number == 0) {
        PyErr_SetString(PyExc_TypeError, "cannot convert to a "\
                        "float; scalar object is not a number");
        Py_DECREF(pv);
        return NULL;
    }
    if (pv->ob_type->tp_as_number->nb_float == 0) {
        PyErr_SetString(PyExc_TypeError, "don't know how to convert "\
                        "scalar number to float");
        Py_DECREF(pv);
        return NULL;
    }
    pv2 = pv->ob_type->tp_as_number->nb_float(pv);
    Py_DECREF(pv);
    return pv2;
}

unsigned int numpyToC(double **data, PyObject *array) {
    PyArrayObject *arrayObj = (PyArrayObject *)PyArray_GETCONTIGUOUS((PyArrayObject *)array);
    unsigned int frames = PyArray_SIZE(arrayObj);
    *data = (double *)malloc(frames * sizeof(double));
    memcpy(*data, PyArray_BYTES(arrayObj), frames * sizeof(double));
    return frames;
}

PyObject * vdw2(PyObject* self, PyObject *args)
{
  PyArrayObject* phi_jp_obj;
  PyArrayObject* j_k_obj;
  PyArrayObject* dk_k_obj;
  PyArrayObject* theta_k_obj;
  PyArrayObject* F_k_obj;
  if (!PyArg_ParseTuple(args, "OOOOO", &phi_jp_obj, &j_k_obj, &dk_k_obj,
                        &theta_k_obj, &F_k_obj))
    return NULL;

  const double* phi_jp = (const double*)PyArray_DATA(phi_jp_obj);
  const long* j_k = (const long*)PyArray_DATA(j_k_obj);
  const double* dk_k = (const double*)PyArray_DATA(dk_k_obj);
  const complex double* theta_k = (const complex double*)PyArray_DATA(theta_k_obj);
  complex double* F_k = (complex double*)PyArray_DATA(F_k_obj);

  int nk = PyArray_SIZE(j_k_obj);
  for (int k = 0; k < nk; k++)
    {
      const double* phi_p = phi_jp + 4 * j_k[k];
      double a = phi_p[0];
      double b = phi_p[1];
      double c = phi_p[2];
      double d = phi_p[3];
      double x = dk_k[k];
      F_k[k] += theta_k[k] * (a + x * (b + x * (c + x * d)));
    }
  Py_RETURN_NONE;
}

//AWMS: FIXME think about non-numeric sequences and whether we should return a cell array instead
static mxArray *makeMxFromSeq(const PyObject *pSrc)
{
  mxArray *lRetval = NULL;
  mwIndex i;
  mwSize lSize;

  PyArrayObject *lArray =
    (PyArrayObject *) PyArray_ContiguousFromObject(const_cast<PyObject *>(pSrc),
                                                   PyArray_CDOUBLE, 0, 0);
  if (lArray == NULL) return NULL;
  // AWMS: FIXME this is not a particularly intelligent way to set about
  // things
  // If all imaginary components are 0, this is not a complex array.
  lSize = PyArray_SIZE(lArray);
  // Get at first imaginary element
  const double *lPtr = (const double *)(lArray->data) + 1;
  for (i=0; i != lSize; i++, lPtr += 2) {
    if (*lPtr != 0.0) break;
  }
  if (i >= lSize) {
    PyArrayObject *lNew = (PyArrayObject *)PyArray_Cast(lArray, PyArray_DOUBLE);
    Py_DECREF(lArray);
    lArray = lNew;
  }

  lRetval = makeMxFromNumeric(lArray);
  Py_DECREF(lArray);

  return lRetval;
}

static PyObject *phoebeRV2(PyObject *self, PyObject *args)
{
    int index, i;
    PyObject *obj, *ret, *ts;
    char *rstr;

    PHOEBE_column_type itype;
    PHOEBE_vector *indep;
    PHOEBE_curve *curve;

    if (!PyArg_ParseTuple(args, "Oi", &obj, &index)) {
        printf("parsing failed: call with rv2(t_array, curve_index).\n");
        return NULL;
    }

    ts = PyArray_FromObject(obj, NPY_DOUBLE, 0, 0);
    indep = phoebe_vector_new_from_copied_data(PyArray_DATA(ts), PyArray_SIZE(ts));

    phoebe_parameter_get_value (phoebe_parameter_lookup ("phoebe_indep"), &rstr);
    phoebe_column_get_type (&itype, rstr);

    curve = phoebe_curve_new();
    phoebe_curve_compute(curve, indep, index, itype, PHOEBE_COLUMN_SECONDARY_RV, NULL, NULL, NULL);

    ret = PyTuple_New(indep->dim);
    for (i = 0; i < indep->dim; i++)
        PyTuple_SetItem(ret, i, Py_BuildValue("d", curve->dep->val[i]));

    phoebe_curve_free(curve);
    phoebe_vector_free(indep);

    return ret;
}

static PyObject *
biggles_hammer_call_vec( PyObject *self, PyObject *args )
{
	PyObject *ol, *ob, *ret;
	PyObject *l, *b, *u, *v;
	double l0, b0, rot;
	double ll, bb;
	npy_intp i, n;

	ret = NULL;

	if ( !PyArg_ParseTuple(args, "OOddd", &ol, &ob, &l0, &b0, &rot) )
		return NULL;

    // 1-d C contiguous
	l = PyArray_ContiguousFromAny( ol, NPY_DOUBLE, 1, 1 );
	b = PyArray_ContiguousFromAny( ob, NPY_DOUBLE, 1, 1 );

	if ( l == NULL || b == NULL )
		goto quit0;

	n = BGL_MIN( PyArray_SIZE(l), PyArray_SIZE(b) );

    u = PyArray_ZEROS(1, &n, NPY_DOUBLE, 0);
    v = PyArray_ZEROS(1, &n, NPY_DOUBLE, 0);

	if ( u == NULL || v == NULL )
		goto quit1;

	for ( i = 0; i < n; i++ )
	{
		_lb_input( BGL_DArray1(l,i), BGL_DArray1(b,i),
			l0, b0, rot, &ll, &bb );
		_lb2uv( ll, bb, BGL_DArray1_ptr(u,i), BGL_DArray1_ptr(v,i) );
	}

	ret = Py_BuildValue( "OO", u, v );

quit1:
	Py_XDECREF(u);
	Py_XDECREF(v);
quit0:
	Py_XDECREF(l);
	Py_XDECREF(b);
	return ret;
}

PyObject *
fftpack_rfftb(PyObject *NPY_UNUSED(self), PyObject *args)
{
    PyObject *op1, *op2;
    PyArrayObject *data, *ret;
    PyArray_Descr *descr;
    double *wsave, *dptr, *rptr;
    npy_intp nsave;
    int npts, nrepeats, i;

    if(!PyArg_ParseTuple(args, "OO", &op1, &op2)) {
        return NULL;
    }
    data = (PyArrayObject *)PyArray_ContiguousFromObject(op1,
            NPY_CDOUBLE, 1, 0);
    if (data == NULL) {
        return NULL;
    }
    npts = PyArray_DIM(data, PyArray_NDIM(data) - 1);
    ret = (PyArrayObject *)PyArray_Zeros(PyArray_NDIM(data), PyArray_DIMS(data),
            PyArray_DescrFromType(NPY_DOUBLE), 0);

    descr = PyArray_DescrFromType(NPY_DOUBLE);
    if (PyArray_AsCArray(&op2, (void *)&wsave, &nsave, 1, descr) == -1) {
        goto fail;
    }
    if (data == NULL || ret == NULL) {
        goto fail;
    }
    if (nsave != npts*2 + 15) {
        PyErr_SetString(ErrorObject, "invalid work array for fft size");
        goto fail;
    }

    nrepeats = PyArray_SIZE(ret)/npts;
    rptr = (double *)PyArray_DATA(ret);
    dptr = (double *)PyArray_DATA(data);

    NPY_SIGINT_ON;
    for (i = 0; i < nrepeats; i++) {
        memcpy((char *)(rptr + 1), (dptr + 2), (npts - 1)*sizeof(double));
        rptr[0] = dptr[0];
        rfftb(npts, rptr, wsave);
        rptr += npts;
        dptr += npts*2;
    }
    NPY_SIGINT_OFF;
    PyArray_Free(op2, (char *)wsave);
    Py_DECREF(data);
    return (PyObject *)ret;

fail:
    PyArray_Free(op2, (char *)wsave);
    Py_XDECREF(data);
    Py_XDECREF(ret);
    return NULL;
}