diff --git a/cython_numpy/erfa.pyx b/cython_numpy/erfa.pyx index 6504f38..f3943e8 100644 --- a/cython_numpy/erfa.pyx +++ b/cython_numpy/erfa.pyx @@ -3,7 +3,8 @@ cimport numpy as np np.import_array() -__all__ = ['atco13', 'd2dtf'] +__all__ = ['atco13', 'd2dtf', 'aper'] + cdef extern from "erfa.h": int eraAtco13(double rc, double dc, @@ -15,28 +16,86 @@ cdef extern from "erfa.h": double *dob, double *rob, double *eo) int eraD2dtf(const char *scale, int ndp, double d1, double d2, int *iy, int *im, int *id, int ihmsf[4]) + void eraAper(double theta, eraASTROM *astrom) + +cdef struct eraASTROM: + double pmt + double eb[3] + double eh[3] + double em + double v[3] + double bm1 + double bpn[3][3] + double along + double phi + double xpl + double ypl + double sphi + double cphi + double diurab + double eral + double refa + double refb + +dt_eraASTROM = np.dtype([('pmt','d'), + ('eb','d',(3,)), + ('eh','d',(3,)), + ('em','d'), + ('v','d',(3,)), + ('bm1 ','d'), + ('bpn','d',(3,3)), + ('along','d'), + ('phi','d'), + ('xpl','d'), + ('ypl','d'), + ('sphi','d'), + ('cphi','d'), + ('diurab','d'), + ('eral','d'), + ('refa','d'), + ('refb','d')], align=True) + #Note: the pattern used here follows https://github.com/cython/cython/wiki/tutorials-numpy#dimensionally-simple-functions + def atco13(rc, dc, pr, pd, px, rv, utc1, utc2, dut1, elong, phi, hm, xp, yp, phpa, tk, rh, wl): shape = np.broadcast(rc, dc, pr, pd, px, rv, utc1, utc2, dut1, elong, phi, hm, xp, yp, phpa, tk, rh, wl).shape - aob = np.empty(shape, dtype=np.double) - zob = np.empty(shape, dtype=np.double) - hob = np.empty(shape, dtype=np.double) - dob = np.empty(shape, dtype=np.double) - rob = np.empty(shape, dtype=np.double) - eo = np.empty(shape, dtype=np.double) + aob_out = np.empty(shape, dtype=np.double) + zob_out = np.empty(shape, dtype=np.double) + hob_out = np.empty(shape, dtype=np.double) + dob_out = np.empty(shape, dtype=np.double) + rob_out = np.empty(shape, dtype=np.double) + eo_out = np.empty(shape, dtype=np.double) - cdef np.broadcast it = np.broadcast(rc, dc, pr, pd, px, rv, utc1, utc2, dut1, elong, phi, hm, xp, yp, phpa, tk, rh, wl, aob, zob, hob, dob, rob, eo) + cdef np.broadcast it = np.broadcast(rc, dc, pr, pd, px, rv, utc1, utc2, dut1, elong, phi, hm, xp, yp, phpa, tk, rh, wl, aob_out, zob_out, hob_out, dob_out, rob_out, eo_out) - cdef double _aob - cdef double _zob - cdef double _hob - cdef double _dob - cdef double _rob - cdef double _eo + cdef double _rc + cdef double _dc + cdef double _pr + cdef double _pd + cdef double _px + cdef double _rv + cdef double _utc1 + cdef double _utc2 + cdef double _dut1 + cdef double _elong + cdef double _phi + cdef double _hm + cdef double _xp + cdef double _yp + cdef double _phpa + cdef double _tk + cdef double _rh + cdef double _wl + cdef double *_aob + cdef double *_zob + cdef double *_hob + cdef double *_dob + cdef double *_rob + cdef double *_eo while np.PyArray_MultiIter_NOTDONE(it): @@ -58,47 +117,76 @@ def atco13(rc, dc, pr, pd, px, rv, utc1, utc2, dut1, elong, phi, hm, xp, yp, php _tk = (np.PyArray_MultiIter_DATA(it, 15))[0] _rh = (np.PyArray_MultiIter_DATA(it, 16))[0] _wl = (np.PyArray_MultiIter_DATA(it, 17))[0] + _aob = (np.PyArray_MultiIter_DATA(it, 18)) + _zob = (np.PyArray_MultiIter_DATA(it, 19)) + _hob = (np.PyArray_MultiIter_DATA(it, 20)) + _dob = (np.PyArray_MultiIter_DATA(it, 21)) + _rob = (np.PyArray_MultiIter_DATA(it, 22)) + _eo = (np.PyArray_MultiIter_DATA(it, 23)) - ret = eraAtco13(_rc, _dc, _pr, _pd, _px, _rv, _utc1, _utc2, _dut1, _elong, _phi, _hm, _xp, _yp, _phpa, _tk, _rh, _wl, &_aob, &_zob, &_hob, &_dob, &_rob, &_eo) - - (np.PyArray_MultiIter_DATA(it, 18))[0] = _aob - (np.PyArray_MultiIter_DATA(it, 19))[0] = _zob - (np.PyArray_MultiIter_DATA(it, 20))[0] = _hob - (np.PyArray_MultiIter_DATA(it, 21))[0] = _dob - (np.PyArray_MultiIter_DATA(it, 22))[0] = _rob - (np.PyArray_MultiIter_DATA(it, 23))[0] = _eo + ret = eraAtco13(_rc, _dc, _pr, _pd, _px, _rv, _utc1, _utc2, _dut1, _elong, _phi, _hm, _xp, _yp, _phpa, _tk, _rh, _wl, _aob, _zob, _hob, _dob, _rob, _eo) np.PyArray_MultiIter_NEXT(it) - return aob, zob, hob, dob, rob, eo + return aob_out, zob_out, hob_out, dob_out, rob_out, eo_out + + def d2dtf(scale, ndp, d1, d2): - shape = np.broadcast(d1, d2).shape - iy = np.empty(shape, dtype=np.int) - im = np.empty(shape, dtype=np.int) - id = np.empty(shape, dtype=np.int) - ihmsf = np.empty(shape, dtype=[('h','i'),('m','i'),('s','i'),('f','i')]) + shape = np.broadcast(scale, ndp, d1, d2).shape + iy_out = np.empty(shape, dtype=np.int) + im_out = np.empty(shape, dtype=np.int) + id_out = np.empty(shape, dtype=np.int) + ihmsf_out = np.empty(shape, dtype=[('h','i'),('m','i'),('s','i'),('f','i')]) - cdef np.broadcast it = np.broadcast(d1, d2, iy, im, id, ihmsf) + cdef np.broadcast it = np.broadcast(scale, ndp, d1, d2, iy_out, im_out, id_out, ihmsf_out) - cdef int _iy - cdef int _im - cdef int _id - cdef int _ihmsf[4] + cdef char *_scale + cdef int _ndp + cdef double _d1 + cdef double _d2 + cdef int *_iy + cdef int *_im + cdef int *_id + cdef int *_ihmsf while np.PyArray_MultiIter_NOTDONE(it): - _d1 = (np.PyArray_MultiIter_DATA(it, 0))[0] - _d2 = (np.PyArray_MultiIter_DATA(it, 1))[0] + _scale = ( np.PyArray_MultiIter_DATA(it, 0)) + _ndp = ( np.PyArray_MultiIter_DATA(it, 1))[0] + _d1 = (np.PyArray_MultiIter_DATA(it, 2))[0] + _d2 = (np.PyArray_MultiIter_DATA(it, 3))[0] + _iy = ( np.PyArray_MultiIter_DATA(it, 4)) + _im = ( np.PyArray_MultiIter_DATA(it, 5)) + _id = ( np.PyArray_MultiIter_DATA(it, 6)) + _ihmsf = ( np.PyArray_MultiIter_DATA(it, 7)) + + ret = eraD2dtf(_scale, _ndp, _d1, _d2, _iy, _im, _id, _ihmsf) + + np.PyArray_MultiIter_NEXT(it) + + return iy_out, im_out, id_out, ihmsf_out + + + +def aper(theta, astrom): + + shape = np.broadcast(theta, astrom).shape + astrom_out = np.array(np.broadcast_arrays(theta, astrom)[1], dtype=dt_eraASTROM) + + cdef np.broadcast it = np.broadcast(theta, astrom_out) + + cdef double _theta + cdef eraASTROM *_astrom + + while np.PyArray_MultiIter_NOTDONE(it): - ret = eraD2dtf(scale, ndp, _d1, _d2, &_iy, &_im, &_id, _ihmsf) + _theta = ( np.PyArray_MultiIter_DATA(it, 0))[0] + _astrom = (np.PyArray_MultiIter_DATA(it, 1)) - (np.PyArray_MultiIter_DATA(it, 2))[0] = _iy - (np.PyArray_MultiIter_DATA(it, 3))[0] = _im - (np.PyArray_MultiIter_DATA(it, 4))[0] = _id - (np.PyArray_MultiIter_DATA(it, 5))[0:4] = _ihmsf + eraAper(_theta, _astrom) np.PyArray_MultiIter_NEXT(it) - return iy, im, id, ihmsf + return astrom_out \ No newline at end of file diff --git a/tests/test_erfa.py b/tests/test_erfa.py index 80653bc..9bdd66d 100644 --- a/tests/test_erfa.py +++ b/tests/test_erfa.py @@ -22,3 +22,9 @@ print(iy.shape, ihmsf.shape, ihmsf.dtype) print(ihmsf) +astrom = np.zeros([2],dtype=erfa.dt_eraASTROM) +theta = np.arange(0,10.0) +print(theta.shape) +print(astrom.shape) +astrom = erfa.aper(theta[:,None], astrom[None,:]) +print(astrom)