From 049e1daad6fb52b9beb6be40f9778c733b9755c1 Mon Sep 17 00:00:00 2001
From: Thibaut Paumard <thibaut.paumard@obspm.fr>
Date: Fri, 9 Aug 2019 09:37:16 +0200
Subject: [PATCH 01/17] python: install util.rayTrace as Scenery.rayTrace

---
 python/example.py        |  6 ++----
 python/gyoto.i           | 11 +++++++++++
 python/gyoto/__init__.py |  8 +++++++-
 python/gyoto/util.py     | 36 +++++++++++++++++++++++-------------
 4 files changed, 43 insertions(+), 18 deletions(-)

diff --git a/python/example.py b/python/example.py
index 99bb8fa9..da5507e6 100644
--- a/python/example.py
+++ b/python/example.py
@@ -30,8 +30,7 @@
 scr=gyoto.core.Screen()
 gg=gyoto.std.KerrBL()
 scr.metric(gg)
-pos=numpy.zeros(4, float)
-scr.getObserverPos(pos)
+pos=scr.getObserverPos()
 
 # Load Scenery
 
@@ -87,9 +86,8 @@
 # Ray-trace scenery
 
 # For that, we can use the short-hand:
-import gyoto.util
 sc.requestedQuantitiesString('Intensity EmissionTime MinDistance')
-results=gyoto.util.rayTrace(sc)
+results=sc.rayTrace()
 
 plt.imshow(results['Intensity'])
 plt.show()
diff --git a/python/gyoto.i b/python/gyoto.i
index 29f8ee8c..0a8cb3de 100644
--- a/python/gyoto.i
+++ b/python/gyoto.i
@@ -682,6 +682,17 @@ ExtendArrayNumPy(array_size_t, size_t);
 };
 %include "GyotoWorldline.h"
 
+%extend Gyoto::Screen {
+  // Support this syntax:
+  // pos = scr.getObserverPos()
+  void getObserverPos(double ARGOUT_ARRAY1[4]) {
+    ($self)->getObserverPos(ARGOUT_ARRAY1);
+  }
+  void getFourVel(double ARGOUT_ARRAY1[4]) {
+    ($self)->getFourVel(ARGOUT_ARRAY1);
+  }
+ };
+GyotoSmPtrClass(Screen)
 GyotoSmPtrClass(Screen)
 GyotoSmPtrClass(Scenery)
 %ignore Gyoto::Photon::Refined;
diff --git a/python/gyoto/__init__.py b/python/gyoto/__init__.py
index 06a65507..29ea0b2f 100644
--- a/python/gyoto/__init__.py
+++ b/python/gyoto/__init__.py
@@ -4,5 +4,11 @@
 future release. Please update your code to import gyoto.core instead.
 
 """
-
+# For backwards compatibility, expose gyoto.core as gyoto
 from gyoto.core import *
+
+# Provide a Pythonic wrapper around Scenery.rayTrace.
+# The underlying C++-like interface remains accessible.
+from gyoto import core, util
+core.Scenery.rayTrace = util.rayTrace
+core.Scenery.rayTrace.__doc__ += core._core.Scenery_rayTrace.__doc__
diff --git a/python/gyoto/util.py b/python/gyoto/util.py
index 6dff01dd..4345bf8a 100644
--- a/python/gyoto/util.py
+++ b/python/gyoto/util.py
@@ -34,24 +34,34 @@ def Coord1dSet(k, res, sz):
 def rayTrace(sc, coord2dset=None, width=None, height=None, fmt='\r j = '):
     '''Ray-trace scenery
 
-    Synopsis:
-    results=rayTrace(scenery)
+First form:
 
-    Input:
-    scenery -- a Gyoto Scenery object or XML file name.
+results=scenery.rayTrace([coord2dset [,width, height] [,fmt]])
 
-    Output:
-    results -- dict containing the various requested quantities as per
-               scenery.requestedQuantitiesString().
+optional parameters:
+coord2dset -- something to specify which part of the field to trace.
+           can be: None to trace all pixels (default)
+           a Coord2dSet instance
+width, height -- horizontal and vertical resolution (overrides what
+           is specified in scenery.screen().resolution()
+fmt     -- prefix to be written in front of the row number for
+           progress output
 
-    Keywords:
-    width, height -- width and height of the output image(s).
+Output:
+results -- dict containing the various requested quantities as per
+           scenery.requestedQuantitiesString().
 
-    TODO:
-    Support and debug various kinds of coord2dset to allow ray-tracing
-    only part of the Scenery.
+TODO:
+Support and debug various kinds of coord2dset to allow ray-tracing
+only part of the Scenery.
+
+
+Second form:
+'''
+    if isinstance(width, core.AstrobjProperties):
+        core._core.Scenery_rayTrace(sc, coord2dset, width, height)
+        return
 
-    '''
     # If needed, read sc
     if type(sc) is str:
         sc=core.Factory(sc).scenery()

From a0f2db907499a768e9d0c8322bd87386d5326b0a Mon Sep 17 00:00:00 2001
From: Thibaut Paumard <thibaut.paumard@obspm.fr>
Date: Sun, 11 Aug 2019 06:21:56 +0200
Subject: [PATCH 02/17] gyoto.util.rayTrace: accept paremeters i and j

---
 python/gyoto/util.py | 41 ++++++++++++++++++++++-------------------
 1 file changed, 22 insertions(+), 19 deletions(-)

diff --git a/python/gyoto/util.py b/python/gyoto/util.py
index 4345bf8a..5b9f7f1b 100644
--- a/python/gyoto/util.py
+++ b/python/gyoto/util.py
@@ -31,7 +31,7 @@ def Coord1dSet(k, res, sz):
             k=core.Angles(k)
     return k
 
-def rayTrace(sc, coord2dset=None, width=None, height=None, fmt='\r j = '):
+def rayTrace(sc, i=None, j=None, coord2dset=core.Grid, width=None, height=None, fmt='\r j = '):
     '''Ray-trace scenery
 
 First form:
@@ -58,8 +58,14 @@ def rayTrace(sc, coord2dset=None, width=None, height=None, fmt='\r j = '):
 
 Second form:
 '''
-    if isinstance(width, core.AstrobjProperties):
-        core._core.Scenery_rayTrace(sc, coord2dset, width, height)
+    if isinstance(j, core.AstrobjProperties):
+        ij=i
+        aop=j
+        if not isinstance(coord2dset, type):
+            ipct=coord2dset
+        else:
+            ipct=None
+        core._core.Scenery_rayTrace(sc, i, j, ipct)
         return
 
     # If needed, read sc
@@ -80,28 +86,25 @@ def rayTrace(sc, coord2dset=None, width=None, height=None, fmt='\r j = '):
 
     # Prepare coord2dset
     nx=None
-    if type(coord2dset) is tuple or coord2dset is None:
-        if type(coord2dset) is tuple:
-            if len(coord2dset) != 2:
-                raise ValueError('if coord2dset is a tuple, it must have 2 items (i, j)')
-            j=coord2dset[0]
-            i=coord2dset[1]
-        else:
-            j=None
-            i=None
+    if isinstance(coord2dset, type):
+        if not issubclass(coord2dset, core.Coord2dSet):
+            raise TypeError("when coord2dset is a type, it must be a subclass of gyoto.core.Coord2dSet")
         if not isinstance(i, core.Coord1dSet):
             i=Coord1dSet(i, res, width)
         if not isinstance(j, core.Coord1dSet):
             j=Coord1dSet(j, res, height)
-        coord2dset=core.Grid(i, j, fmt)
+        try:
+            coord2dset=coord2dset(i, j, fmt)
+        except TypeError:
+            coord2dset=coord2dset(i, j)
         nx=i.size()
-        ny=j.size()
-    elif not isinstance(coord2dset, core.Coord2dSet):
-        raise ValueError('this type of 2D set not yet implemented')
-
-    if nx is None:
-        nx=coord2dset.size()
+        ntot=coord2dset.size()
+        ny=ntot//nx
+    elif isinstance(coord2dset, core.Coord2dSet):
+        nx=ntot=coord2dset.size()
         ny=1
+    else:
+        raise TypeError('coord2dset must be a gyoto.core.Coord2dSet subclass or instance')
 
     # Prepare arrays to store results
     res = dict()

From 7d86c285b0b9fb095bbcb102fb39fdc0f96e076d Mon Sep 17 00:00:00 2001
From: Thibaut Paumard <thibaut.paumard@obspm.fr>
Date: Sun, 11 Aug 2019 07:21:59 +0200
Subject: [PATCH 03/17] gyoto.util.rayTrace: support indices, ranges and list
 of indices

---
 lib/Screen.C         |  1 +
 python/gyoto.i       |  7 +++++
 python/gyoto/util.py | 63 +++++++++++++++++++++++++++++++++++++-------
 3 files changed, 61 insertions(+), 10 deletions(-)

diff --git a/lib/Screen.C b/lib/Screen.C
index 201b608f..7ebea872 100644
--- a/lib/Screen.C
+++ b/lib/Screen.C
@@ -895,6 +895,7 @@ void Screen::fitsWriteMask(string const &fname) {
 #endif
 
 bool Screen::operator()(size_t i, size_t j) {
+  GYOTO_DEBUG << "i=" << i << ", j=" << j << endl;
   if ( i<=0 || i> npix_ || j>npix_ || j<=0) GYOTO_ERROR("wrong index");
   if (!mask_) return true;
   return mask_[i-1+npix_*(j-1)];
diff --git a/python/gyoto.i b/python/gyoto.i
index 0a8cb3de..c7006175 100644
--- a/python/gyoto.i
+++ b/python/gyoto.i
@@ -826,6 +826,12 @@ public:
   virtual double angle() const ;
   virtual Coord1dSet& operator++()=0;
 };
+%extend Coord1dSet {
+  // Get value of coord set
+  size_t value() const {
+    return **($self);
+  }
+};
 
 class Coord2dSet {
 public:
@@ -902,6 +908,7 @@ public:
   size_t size();
   Coord1dSet& operator++();
   size_t operator*() const ;
+  size_t index() const ;
 };
 %extend Indices {
   Indices (size_t DIM1, size_t *IN_ARRAY1) {
diff --git a/python/gyoto/util.py b/python/gyoto/util.py
index 5b9f7f1b..219f2bf2 100644
--- a/python/gyoto/util.py
+++ b/python/gyoto/util.py
@@ -12,24 +12,64 @@
 import numbers
 
 def Coord1dSet(k, res, sz):
-    '''Easily initialize a gyoto.core.Coord1dset'''
+    '''Easily initialize a gyoto.core.Coord1dset
+
+synopsis:
+set, data = Coord1dSet(k, res, sz)
+
+parameters:
+k   -- something to convert to Coord1dSet
+res -- Screen resolution
+sz  -- Screen size in that direction
+
+caveats:
+Gyoto indices start at 1. This function takes care of tranlating from
+the more standard Python 0-based indices to Gyoto 1-based indices.
+
+The instances of gyoto.core.Indices need to access a buffer that must
+be preallocated and survive as long as themselves. The second output
+(data) when ot None, is such a buffer. Make sure it remains in scope
+and is not destroyed as long as you need the first output (set).
+
+returns:
+if k is None:
+  data is None
+  set is a gyoto.core.Range() covering all pixels according to res and sz
+
+if k is a Python range:
+  data is None
+  set is the corresponding gyoto.core.Range
+
+if k is a scalar integer:
+  data is a NumPy array containing this scalar+1
+  set is a gyoto.core.Indices instances pointing to this array
+
+if k is an array-like continging only integers:
+  data is a NumPy array containing these integers+1
+  set is a gyoto.core.Indices instances pointing to this array
+
+    '''
+    data=None
     if (k is None):
         k=core.Range(res//2-sz//2+1, res//2-sz//2+sz, 1)
     elif type(k) is range:
-        k=core.Range(k.start, k.stop-1, k.step)
+        k=core.Range(k.start+1, k.stop, k.step)
     elif numpy.isscalar(k):
         if isinstance(k, numbers.Integral):
-            k=core.Indices([k])
+            data=numpy.array([k+1], numpy.uint64)
+            k=core.Indices(data)
         elif isinstance(k, numbers.Real):
             k=core.Angles([k])
         else:
             raise ValueError('unknown scalar type')
     else:
         if all([isinstance(n, numbers.Integral) for n in k]):
-            k=core.Indices(k)
+            data=numpy.asarray(k, dtype=numpy.uint64)
+            data += 1
+            k=core.Indices(data)
         else:
             k=core.Angles(k)
-    return k
+    return k, data
 
 def rayTrace(sc, i=None, j=None, coord2dset=core.Grid, width=None, height=None, fmt='\r j = '):
     '''Ray-trace scenery
@@ -39,9 +79,12 @@ def rayTrace(sc, i=None, j=None, coord2dset=core.Grid, width=None, height=None,
 results=scenery.rayTrace([coord2dset [,width, height] [,fmt]])
 
 optional parameters:
-coord2dset -- something to specify which part of the field to trace.
-           can be: None to trace all pixels (default)
-           a Coord2dSet instance
+i       -- horizontal specification of the part of the field to trace
+           (see gyoto.util.Coord1dSet)
+j       -- vertical specification of the part of the field to trace
+           (see gyoto.util.Coord1dSet)
+coord2dset -- a Coord2dSet subclass. Default: gyoto.core.Grid. The other
+           value that makes sense is gyoto.core.Bucket.
 width, height -- horizontal and vertical resolution (overrides what
            is specified in scenery.screen().resolution()
 fmt     -- prefix to be written in front of the row number for
@@ -90,9 +133,9 @@ def rayTrace(sc, i=None, j=None, coord2dset=core.Grid, width=None, height=None,
         if not issubclass(coord2dset, core.Coord2dSet):
             raise TypeError("when coord2dset is a type, it must be a subclass of gyoto.core.Coord2dSet")
         if not isinstance(i, core.Coord1dSet):
-            i=Coord1dSet(i, res, width)
+            i, idata=Coord1dSet(i, res, width)
         if not isinstance(j, core.Coord1dSet):
-            j=Coord1dSet(j, res, height)
+            j, jdata=Coord1dSet(j, res, height)
         try:
             coord2dset=coord2dset(i, j, fmt)
         except TypeError:

From 7852e47d204d4b6f87032abc36179a4cb23f479b Mon Sep 17 00:00:00 2001
From: Thibaut Paumard <thibaut.paumard@obspm.fr>
Date: Sun, 11 Aug 2019 07:45:42 +0200
Subject: [PATCH 04/17] gyoto.core.Scenery: add __getitem__ method

shortcut for rayTrace
---
 python/gyoto/__init__.py |  1 +
 python/gyoto/util.py     | 20 ++++++++++++++++++++
 2 files changed, 21 insertions(+)

diff --git a/python/gyoto/__init__.py b/python/gyoto/__init__.py
index 29ea0b2f..626d18ee 100644
--- a/python/gyoto/__init__.py
+++ b/python/gyoto/__init__.py
@@ -11,4 +11,5 @@
 # The underlying C++-like interface remains accessible.
 from gyoto import core, util
 core.Scenery.rayTrace = util.rayTrace
+core.Scenery.__getitem__ = util.Scenery_getitem
 core.Scenery.rayTrace.__doc__ += core._core.Scenery_rayTrace.__doc__
diff --git a/python/gyoto/util.py b/python/gyoto/util.py
index 219f2bf2..e8c87cc4 100644
--- a/python/gyoto/util.py
+++ b/python/gyoto/util.py
@@ -54,6 +54,22 @@ def Coord1dSet(k, res, sz):
         k=core.Range(res//2-sz//2+1, res//2-sz//2+sz, 1)
     elif type(k) is range:
         k=core.Range(k.start+1, k.stop, k.step)
+    elif type(k) is slice:
+        start=k.start
+        if start is None:
+            start=0
+        elif start < 0:
+            start=res+start
+        start += 1
+        stop=k.stop
+        if stop is None:
+            stop=res
+        elif stop < 0:
+            stop=res+stop
+        step=k.step
+        if step is None:
+            step=1
+        k=core.Range(start, stop, step)
     elif numpy.isscalar(k):
         if isinstance(k, numbers.Integral):
             data=numpy.array([k+1], numpy.uint64)
@@ -264,6 +280,10 @@ def rayTrace(sc, i=None, j=None, coord2dset=core.Grid, width=None, height=None,
 
     return res
 
+def Scenery_getitem(self, args):
+    '''Shortcut for Scenery.rayTrace(i, j)'''
+    self.rayTrace(args[0], args[1])
+
 def readScenery(filename):
     '''Read Scenery from XML file'''
     return core.Factory(filename).scenery()

From 18d2c3dfd3cd53e228b8ed0c91572a6232c07d04 Mon Sep 17 00:00:00 2001
From: Thibaut Paumard <thibaut.paumard@obspm.fr>
Date: Sun, 11 Aug 2019 07:53:12 +0200
Subject: [PATCH 05/17] gyoto.util: support angles in rayTrace

---
 python/gyoto/util.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/python/gyoto/util.py b/python/gyoto/util.py
index e8c87cc4..d4b15909 100644
--- a/python/gyoto/util.py
+++ b/python/gyoto/util.py
@@ -75,6 +75,7 @@ def Coord1dSet(k, res, sz):
             data=numpy.array([k+1], numpy.uint64)
             k=core.Indices(data)
         elif isinstance(k, numbers.Real):
+            data=numpy.array([k])
             k=core.Angles([k])
         else:
             raise ValueError('unknown scalar type')
@@ -84,6 +85,7 @@ def Coord1dSet(k, res, sz):
             data += 1
             k=core.Indices(data)
         else:
+            data=numpy.asarray(k)
             k=core.Angles(k)
     return k, data
 

From f92a846193e3bd3cfa2adc5a04371d82b5158156 Mon Sep 17 00:00:00 2001
From: Thibaut Paumard <thibaut.paumard@obspm.fr>
Date: Sun, 11 Aug 2019 14:00:19 +0200
Subject: [PATCH 06/17] gyoto.util.rayTrace/Scenery_getitem: improve
 dimensionality

- Scenery_getitem accepts all rayTrace parameters
- return 1D array when the Coord2dSet is 1D
---
 python/gyoto/util.py | 84 +++++++++++++++++++++++++-------------------
 1 file changed, 47 insertions(+), 37 deletions(-)

diff --git a/python/gyoto/util.py b/python/gyoto/util.py
index d4b15909..bf23cc5f 100644
--- a/python/gyoto/util.py
+++ b/python/gyoto/util.py
@@ -167,113 +167,122 @@ def rayTrace(sc, i=None, j=None, coord2dset=core.Grid, width=None, height=None,
     else:
         raise TypeError('coord2dset must be a gyoto.core.Coord2dSet subclass or instance')
 
+    if isinstance(coord2dset, core.Grid):
+        dims=(ny, nx)
+        array_double_fromnumpy1or2=core.array_double_fromnumpy2
+        array_double_fromnumpy2or3=core.array_double_fromnumpy3
+    else:
+        dims=(ntot,)
+        array_double_fromnumpy1or2=core.array_double_fromnumpy1
+        array_double_fromnumpy2or3=core.array_double_fromnumpy2
+
     # Prepare arrays to store results
     res = dict()
     aop=core.AstrobjProperties()
-    aop.offset=nx*ny
+    aop.offset=ntot
 
     if 'Spectrum' in sc.requestedQuantitiesString():
         nsamples=sc.screen().spectrometer().nSamples()
 
     if 'Intensity' in sc.requestedQuantitiesString():
-        intensity=numpy.zeros((ny, nx))
-        pintensity=core.array_double_fromnumpy2(intensity)
+        intensity=numpy.zeros(dims)
+        pintensity=array_double_fromnumpy1or2(intensity)
         aop.intensity=pintensity
         res['Intensity']=intensity
 
     if 'EmissionTime' in sc.requestedQuantitiesString():
-        time=numpy.zeros((ny, nx))
-        ptime=core.array_double_fromnumpy2(time)
+        time=numpy.zeros(dims)
+        ptime=array_double_fromnumpy1or2(time)
         aop.time=ptime
         res['EmissionTime'] = time
 
     if 'MinDistance' in sc.requestedQuantitiesString():
-        distance=numpy.zeros((ny, nx))
-        pdistance=core.array_double_fromnumpy2(distance)
+        distance=numpy.zeros(dims)
+        pdistance=array_double_fromnumpy1or2(distance)
         aop.distance=pdistance
         res['MinDistance'] = distance
 
     if 'FirstDistMin' in sc.requestedQuantitiesString():
-        first_dmin=numpy.zeros((ny, nx))
-        pfirst_dmin=core.array_double_fromnumpy2(first_dmin)
+        first_dmin=numpy.zeros(dims)
+        pfirst_dmin=array_double_fromnumpy1or2(first_dmin)
         aop.first_dmin=pfirst_dmin
         res['FirstDistMin'] = first_dmin
 
     if 'Redshift' in sc.requestedQuantitiesString():
-        redshift=numpy.zeros((ny, nx))
-        predshift=core.array_double_fromnumpy2(redshift)
+        redshift=numpy.zeros(dims)
+        predshift=array_double_fromnumpy1or2(redshift)
         aop.redshift=predshift
         res['Redshift'] = redshift
 
     if 'ImpactCoords' in sc.requestedQuantitiesString():
-        impactcoords=numpy.zeros((ny, nx, 16))
-        pimpactcoords=core.array_double_fromnumpy3(impactcoords)
+        impactcoords=numpy.zeros(dims+(16,))
+        pimpactcoords=array_double_fromnumpy2or3(impactcoords)
         aop.impactcoords=pimpactcoords
         res['ImpactCoords'] = impactcoords
 
     if 'User1' in sc.requestedQuantitiesString():
-        user1=numpy.zeros((ny, nx))
-        puser1=core.array_double_fromnumpy2(user1)
+        user1=numpy.zeros(dims)
+        puser1=array_double_fromnumpy1or2(user1)
         aop.user1=puser1
         res['User1'] = user1
 
     if 'User1' in sc.requestedQuantitiesString():
-        impactcoords=numpy.zeros((ny, nx))
-        pimpactcoords=core.array_double_fromnumpy2(impactcoords)
+        impactcoords=numpy.zeros(dims)
+        pimpactcoords=array_double_fromnumpy1or2(impactcoords)
         aop.impactcoords=pimpactcoords
         res['User1'] = impactcoords
 
     if 'User2' in sc.requestedQuantitiesString():
-        user2=numpy.zeros((ny, nx))
-        puser2=core.array_double_fromnumpy2(user2)
+        user2=numpy.zeros(dims)
+        puser2=array_double_fromnumpy1or2(user2)
         aop.user2=puser2
         res['User2'] = user2
 
     if 'User3' in sc.requestedQuantitiesString():
-        user3=numpy.zeros((ny, nx))
-        puser3=core.array_double_fromnumpy2(user3)
+        user3=numpy.zeros(dims)
+        puser3=array_double_fromnumpy1or2(user3)
         aop.user3=puser3
         res['User3'] = user3
 
     if 'User4' in sc.requestedQuantitiesString():
-        user4=numpy.zeros((ny, nx))
-        puser4=core.array_double_fromnumpy2(user4)
+        user4=numpy.zeros(dims)
+        puser4=array_double_fromnumpy1or2(user4)
         aop.user4=puser4
         res['User4'] = user4
 
     if 'User5' in sc.requestedQuantitiesString():
-        user5=numpy.zeros((ny, nx))
-        puser5=core.array_double_fromnumpy2(user5)
+        user5=numpy.zeros(dims)
+        puser5=array_double_fromnumpy1or2(user5)
         aop.user5=puser5
         res['User5'] = user5
 
     if 'Spectrum' in sc.requestedQuantitiesString():
-        spectrum=numpy.zeros((nsamples, ny, nx))
-        pspectrum=core.array_double_fromnumpy3(spectrum)
+        spectrum=numpy.zeros((nsamples,)+dims)
+        pspectrum=array_double_fromnumpy2or3(spectrum)
         aop.spectrum=pspectrum
         res['Spectrum'] = spectrum
 
     if 'SpectrumStokesQ' in sc.requestedQuantitiesString():
-        stokesQ=numpy.zeros((nsamples, ny, nx))
-        pstokesQ=core.array_double_fromnumpy3(stokesQ)
+        stokesQ=numpy.zeros((nsamples,)+dims)
+        pstokesQ=array_double_fromnumpy2or3(stokesQ)
         aop.stokesQ=pstokesQ
         res['SpectrumStokesQ'] = stokesQ
 
     if 'SpectrumStokesU' in sc.requestedQuantitiesString():
-        stokesU=numpy.zeros((nsamples, ny, nx))
-        pstokesU=core.array_double_fromnumpy3(stokesU)
+        stokesU=numpy.zeros((nsamples,)+dims)
+        pstokesU=array_double_fromnumpy2or3(stokesU)
         aop.stokesU=pstokesU
         res['SpectrumStokesU'] = stokesU
 
     if 'SpectrumStokesV' in sc.requestedQuantitiesString():
-        stokesV=numpy.zeros((nsamples, ny, nx))
-        pstokesV=core.array_double_fromnumpy3(stokesV)
+        stokesV=numpy.zeros((nsamples,)+dims)
+        pstokesV=array_double_fromnumpy2or3(stokesV)
         aop.stokesV=pstokesV
         res['SpectrumStokesV'] = stokesV
 
     if 'BinSpectrum' in sc.requestedQuantitiesString():
-        binspectrum=numpy.zeros((nsamples, ny, nx))
-        pbinspectrum=core.array_double_fromnumpy3(binspectrum)
+        binspectrum=numpy.zeros((nsamples,)+dims)
+        pbinspectrum=array_double_fromnumpy2or3(binspectrum)
         aop.binspectrum=pbinspectrum
         res['BinSpectrum'] = binspectrum
 
@@ -283,8 +292,9 @@ def rayTrace(sc, i=None, j=None, coord2dset=core.Grid, width=None, height=None,
     return res
 
 def Scenery_getitem(self, args):
-    '''Shortcut for Scenery.rayTrace(i, j)'''
-    self.rayTrace(args[0], args[1])
+    '''Shortcut for Scenery.rayTrace(i, j)
+'''+rayTrace.__doc__
+    return self.rayTrace(*args)
 
 def readScenery(filename):
     '''Read Scenery from XML file'''

From 6614509bccc4937123963c96011af8955fa04a51 Mon Sep 17 00:00:00 2001
From: Thibaut Paumard <thibaut.paumard@obspm.fr>
Date: Sun, 11 Aug 2019 14:07:43 +0200
Subject: [PATCH 07/17] gyoto.util: improve doc comments

---
 python/gyoto/__init__.py | 3 ++-
 python/gyoto/util.py     | 2 +-
 2 files changed, 3 insertions(+), 2 deletions(-)

diff --git a/python/gyoto/__init__.py b/python/gyoto/__init__.py
index 626d18ee..e6e4de37 100644
--- a/python/gyoto/__init__.py
+++ b/python/gyoto/__init__.py
@@ -11,5 +11,6 @@
 # The underlying C++-like interface remains accessible.
 from gyoto import core, util
 core.Scenery.rayTrace = util.rayTrace
-core.Scenery.__getitem__ = util.Scenery_getitem
 core.Scenery.rayTrace.__doc__ += core._core.Scenery_rayTrace.__doc__
+core.Scenery.__getitem__ = util.Scenery_getitem
+core.Scenery.__getitem__.__doc__ += core.Scenery.rayTrace.__doc__
diff --git a/python/gyoto/util.py b/python/gyoto/util.py
index bf23cc5f..f6cfd39e 100644
--- a/python/gyoto/util.py
+++ b/python/gyoto/util.py
@@ -293,7 +293,7 @@ def rayTrace(sc, i=None, j=None, coord2dset=core.Grid, width=None, height=None,
 
 def Scenery_getitem(self, args):
     '''Shortcut for Scenery.rayTrace(i, j)
-'''+rayTrace.__doc__
+'''
     return self.rayTrace(*args)
 
 def readScenery(filename):

From c4af46c430bb0cd73c6230571d3530146b0f126d Mon Sep 17 00:00:00 2001
From: Thibaut Paumard <thibaut.paumard@obspm.fr>
Date: Wed, 14 Aug 2019 14:31:42 +0200
Subject: [PATCH 08/17] Gyoto::Screen::Coord1dSet: copy beffur in Indices and
 Angles

Much safer in Python
---
 include/GyotoScreen.h |  6 ++++--
 lib/Screen.C          | 24 ++++++++++++++++++++----
 python/gyoto/util.py  | 27 +++++++++------------------
 3 files changed, 33 insertions(+), 24 deletions(-)

diff --git a/include/GyotoScreen.h b/include/GyotoScreen.h
index fd5ad5ad..a3604ee2 100644
--- a/include/GyotoScreen.h
+++ b/include/GyotoScreen.h
@@ -741,11 +741,12 @@ class Gyoto::Screen
     /// 1D specifier for an arbitrary pixel coordinate set.
     class Indices : public Coord1dSet {
     protected:
-      size_t const * const indices_;
+      size_t * indices_;
       size_t const sz_;
       size_t i_;
     public:
       Indices (size_t const*const buf, size_t sz);
+      ~Indices();
       void begin();
       bool valid();
       size_t size();
@@ -757,11 +758,12 @@ class Gyoto::Screen
     /// 1D specifier for an arbitrary angle coordinate set.
     class Angles : public Coord1dSet {
     protected:
-      double const * const buf_;
+      double * buf_;
       size_t const sz_;
       size_t i_;
     public:
       Angles (double const*const buf, size_t sz);
+      ~Angles();
       void begin();
       bool valid();
       size_t size();
diff --git a/lib/Screen.C b/lib/Screen.C
index 7ebea872..08a2fe00 100644
--- a/lib/Screen.C
+++ b/lib/Screen.C
@@ -1465,8 +1465,16 @@ size_t Screen::Range::index() const {return (cur_-mi_) / d_;}
 //////
 
 Screen::Indices::Indices (size_t const*const buf, size_t sz)
-  : Coord1dSet(pixel), indices_(buf), sz_(sz), i_(0)
-{}
+  : Coord1dSet(pixel), indices_(NULL), sz_(sz), i_(0)
+{
+  if (!buf) return;
+  indices_ = new size_t[sz];
+  memcpy(indices_, buf, sz*sizeof(size_t));
+}
+Screen::Indices::~Indices(){
+  if (!indices_) return;
+  delete[] indices_;
+}
 void Screen::Indices::begin() {i_=0;}
 bool Screen::Indices::valid() {return i_ < sz_;}
 size_t Screen::Indices::size(){return sz_;}
@@ -1478,8 +1486,16 @@ size_t Screen::Indices::index() const {return i_;}
 /////
 
 Screen::Angles::Angles (double const*const buf, size_t sz)
-  : Coord1dSet(Screen::angle), buf_(buf), sz_(sz), i_(0)
-{}
+  : Coord1dSet(Screen::angle), buf_(NULL), sz_(sz), i_(0)
+{
+  if (!buf) return;
+  buf_ = new double[sz];
+  memcpy(buf_, buf, sz*sizeof(double));
+}
+Screen::Angles::~Angles(){
+  if (!buf_) return;
+  delete[] buf_;
+}
 void Screen::Angles::begin() {i_=0;}
 bool Screen::Angles::valid() {return i_<sz_;}
 size_t Screen::Angles::size(){return sz_;}
diff --git a/python/gyoto/util.py b/python/gyoto/util.py
index f6cfd39e..9f032da6 100644
--- a/python/gyoto/util.py
+++ b/python/gyoto/util.py
@@ -15,7 +15,7 @@ def Coord1dSet(k, res, sz):
     '''Easily initialize a gyoto.core.Coord1dset
 
 synopsis:
-set, data = Coord1dSet(k, res, sz)
+set = Coord1dSet(k, res, sz)
 
 parameters:
 k   -- something to convert to Coord1dSet
@@ -26,27 +26,18 @@ def Coord1dSet(k, res, sz):
 Gyoto indices start at 1. This function takes care of tranlating from
 the more standard Python 0-based indices to Gyoto 1-based indices.
 
-The instances of gyoto.core.Indices need to access a buffer that must
-be preallocated and survive as long as themselves. The second output
-(data) when ot None, is such a buffer. Make sure it remains in scope
-and is not destroyed as long as you need the first output (set).
-
 returns:
 if k is None:
-  data is None
-  set is a gyoto.core.Range() covering all pixels according to res and sz
+  a gyoto.core.Range() covering all pixels according to res and sz
 
 if k is a Python range:
-  data is None
-  set is the corresponding gyoto.core.Range
+  the corresponding gyoto.core.Range
 
 if k is a scalar integer:
-  data is a NumPy array containing this scalar+1
-  set is a gyoto.core.Indices instances pointing to this array
+  a gyoto.core.Indices instance containing this value
 
-if k is an array-like continging only integers:
-  data is a NumPy array containing these integers+1
-  set is a gyoto.core.Indices instances pointing to this array
+if k is an array-like containing only integers:
+  the corresponding gyoto.core.Indices instance
 
     '''
     data=None
@@ -87,7 +78,7 @@ def Coord1dSet(k, res, sz):
         else:
             data=numpy.asarray(k)
             k=core.Angles(k)
-    return k, data
+    return k
 
 def rayTrace(sc, i=None, j=None, coord2dset=core.Grid, width=None, height=None, fmt='\r j = '):
     '''Ray-trace scenery
@@ -151,9 +142,9 @@ def rayTrace(sc, i=None, j=None, coord2dset=core.Grid, width=None, height=None,
         if not issubclass(coord2dset, core.Coord2dSet):
             raise TypeError("when coord2dset is a type, it must be a subclass of gyoto.core.Coord2dSet")
         if not isinstance(i, core.Coord1dSet):
-            i, idata=Coord1dSet(i, res, width)
+            i=Coord1dSet(i, res, width)
         if not isinstance(j, core.Coord1dSet):
-            j, jdata=Coord1dSet(j, res, height)
+            j=Coord1dSet(j, res, height)
         try:
             coord2dset=coord2dset(i, j, fmt)
         except TypeError:

From e1209889d5511fe9061f7c25155578a44f76e73c Mon Sep 17 00:00:00 2001
From: Thibaut Paumard <thibaut.paumard@obspm.fr>
Date: Wed, 14 Aug 2019 21:22:25 +0200
Subject: [PATCH 09/17] gyoto.util.rayTrace: reorder parameters

The order j, i is more Pythonic than i.j. This ways,
  sc[j, i]['Quantity'] == sc[:,:]['Quantity'][j,i].

Having height and width last allows writing:
  sc[j, i, core.Bucket, None]
skipping height and width to remove progress output.
---
 python/gyoto/util.py | 40 ++++++++++++++++++++++++----------------
 1 file changed, 24 insertions(+), 16 deletions(-)

diff --git a/python/gyoto/util.py b/python/gyoto/util.py
index 9f032da6..818cb044 100644
--- a/python/gyoto/util.py
+++ b/python/gyoto/util.py
@@ -80,44 +80,52 @@ def Coord1dSet(k, res, sz):
             k=core.Angles(k)
     return k
 
-def rayTrace(sc, i=None, j=None, coord2dset=core.Grid, width=None, height=None, fmt='\r j = '):
+def rayTrace(sc,
+             j=None, i=None,
+             coord2dset=core.Grid,
+             prefix='\r j = ',
+             height=None, width=None):
     '''Ray-trace scenery
 
 First form:
 
-results=scenery.rayTrace([coord2dset [,width, height] [,fmt]])
+results=scenery.rayTrace([j, i, [coord2dset [,prefix]]])
 
 optional parameters:
-i       -- horizontal specification of the part of the field to trace
-           (see gyoto.util.Coord1dSet)
 j       -- vertical specification of the part of the field to trace
            (see gyoto.util.Coord1dSet)
+i       -- horizontal specification of the part of the field to trace
+           (see gyoto.util.Coord1dSet)
 coord2dset -- a Coord2dSet subclass. Default: gyoto.core.Grid. The other
            value that makes sense is gyoto.core.Bucket.
-width, height -- horizontal and vertical resolution (overrides what
-           is specified in scenery.screen().resolution()
-fmt     -- prefix to be written in front of the row number for
+prefix  -- prefix to be written in front of the row number for
            progress output
+height, width -- vertical and horizontal resolution (overrides what
+           is specified in scenery.screen().resolution()
 
 Output:
 results -- dict containing the various requested quantities as per
            scenery.requestedQuantitiesString().
 
-TODO:
-Support and debug various kinds of coord2dset to allow ray-tracing
-only part of the Scenery.
+CAVEAT:
+This high level-wrapper is Pythonic and take the arguments as j, i,
+0-based indices whereas most Gyoto functions take them as i, j,
+1-based indices.
 
+TODO:
+Support impactcoords.
 
 Second form:
-'''
-    if isinstance(j, core.AstrobjProperties):
-        ij=i
-        aop=j
+
+    '''
+    if isinstance(i, core.AstrobjProperties):
+        ij=j
+        aop=i
         if not isinstance(coord2dset, type):
             ipct=coord2dset
         else:
             ipct=None
-        core._core.Scenery_rayTrace(sc, i, j, ipct)
+        core._core.Scenery_rayTrace(sc, ij, aop, ipct)
         return
 
     # If needed, read sc
@@ -146,7 +154,7 @@ def rayTrace(sc, i=None, j=None, coord2dset=core.Grid, width=None, height=None,
         if not isinstance(j, core.Coord1dSet):
             j=Coord1dSet(j, res, height)
         try:
-            coord2dset=coord2dset(i, j, fmt)
+            coord2dset=coord2dset(i, j, prefix)
         except TypeError:
             coord2dset=coord2dset(i, j)
         nx=i.size()

From 11b1a6016db823a39d42f993432a8271aee86ede Mon Sep 17 00:00:00 2001
From: Thibaut Paumard <thibaut.paumard@obspm.fr>
Date: Thu, 15 Aug 2019 09:01:36 +0200
Subject: [PATCH 10/17] gyoto.util.rayTrace: return arrays or the right
 dimension

if j and/or i are scalars, decrease dimension of output arrays
---
 python/gyoto/util.py | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/python/gyoto/util.py b/python/gyoto/util.py
index 818cb044..8f557e4f 100644
--- a/python/gyoto/util.py
+++ b/python/gyoto/util.py
@@ -145,6 +145,7 @@ def rayTrace(sc,
     sc.screen().resolution(res)
 
     # Prepare coord2dset
+    scalars=int(numpy.isscalar(i))+int(numpy.isscalar(j))
     nx=None
     if isinstance(coord2dset, type):
         if not issubclass(coord2dset, core.Coord2dSet):
@@ -166,7 +167,7 @@ def rayTrace(sc,
     else:
         raise TypeError('coord2dset must be a gyoto.core.Coord2dSet subclass or instance')
 
-    if isinstance(coord2dset, core.Grid):
+    if isinstance(coord2dset, core.Grid) and scalars is 0 :
         dims=(ny, nx)
         array_double_fromnumpy1or2=core.array_double_fromnumpy2
         array_double_fromnumpy2or3=core.array_double_fromnumpy3
@@ -288,6 +289,10 @@ def rayTrace(sc,
     # Perform the actual ray-tracing
     sc.rayTrace(coord2dset, aop)
 
+    if scalars is 2:
+        for key in res:
+            res[key]=res[key][0]
+
     return res
 
 def Scenery_getitem(self, args):

From 3277ebedd6d095a6f85c928ced59a680f8e90337 Mon Sep 17 00:00:00 2001
From: Thibaut Paumard <thibaut.paumard@obspm.fr>
Date: Thu, 15 Aug 2019 22:35:51 +0200
Subject: [PATCH 11/17] Python: make sure Worldline.getCartesian works both
 with NumPy and C arrays

---
 python/gyoto.i           | 22 ++++------------------
 python/gyoto/__init__.py |  3 +++
 python/gyoto/util.py     | 16 ++++++++++++++++
 python/gyoto_std.i       | 10 ++++++++++
 4 files changed, 33 insertions(+), 18 deletions(-)

diff --git a/python/gyoto.i b/python/gyoto.i
index c7006175..5db81725 100644
--- a/python/gyoto.i
+++ b/python/gyoto.i
@@ -588,10 +588,10 @@ ExtendArrayNumPy(array_size_t, size_t);
 %rename(Worldline__IntegState__Boost) Gyoto::Worldline::IntegState::Boost;
 %rename(Worldline__IntegState__Legacy) Gyoto::Worldline::IntegState::Legacy;
 %ignore Gyoto::Worldline::getCoord(double const * const dates, size_t const n_dates,
-		double * const x1dest,
-		double * const x2dest, double * const x3dest,
-		double * const x0dot=NULL,  double * const x1dot=NULL,
-		double * const x2dot=NULL,  double * const x3dot=NULL) ;
+                                   double * const x1dest,
+                                   double * const x2dest, double * const x3dest,
+                                   double * const x0dot=NULL,  double * const x1dot=NULL,
+                                   double * const x2dot=NULL,  double * const x3dot=NULL) ;
 %extend Gyoto::Worldline {
   void get_t(double * INPLACE_ARRAY1, size_t DIM1) {
     if (DIM1 != ($self)->get_nelements()) GYOTO_ERROR("wrong output array size");
@@ -613,20 +613,6 @@ ExtendArrayNumPy(array_size_t, size_t);
       GYOTO_ERROR("wrong size for output array");
     ($self)->getCoord(dates, n_dates, x1dest, x2dest, x3dest, x0dot, x1dot, x2dot, x3dot);
   }
-  void getCartesian(double * dates, size_t n_dates,
-		double * x1dest, size_t n1,
-		double * x2dest, size_t n2,
-                double * x3dest, size_t n3,
-                double * x1dot, size_t n1d,
-		double * x2dot, size_t n2d,
-                double * x3dot, size_t n3d) {
-    if (n1 != n_dates || n2 != n_dates || n3 != n_dates ||
-        (x1dot && n1d != n_dates) ||
-        (x2dot && n2d != n_dates) ||
-        (x3dot && n3d != n_dates))
-      GYOTO_ERROR("wrong size for output array");
-    ($self)->getCartesian(dates, n_dates, x1dest, x2dest, x3dest, x1dot, x2dot, x3dot);
-  }
   void getCoord(double * dates, size_t n_dates,
 		double * x1dest, size_t n1,
 		double * x2dest, size_t n2,
diff --git a/python/gyoto/__init__.py b/python/gyoto/__init__.py
index e6e4de37..b2c8e3a8 100644
--- a/python/gyoto/__init__.py
+++ b/python/gyoto/__init__.py
@@ -14,3 +14,6 @@
 core.Scenery.rayTrace.__doc__ += core._core.Scenery_rayTrace.__doc__
 core.Scenery.__getitem__ = util.Scenery_getitem
 core.Scenery.__getitem__.__doc__ += core.Scenery.rayTrace.__doc__
+
+core.Worldline.getCartesian =  util._Worldline_getCartesian
+core.Worldline.getCartesian.__doc__ = core._core.Worldline_getCartesian.__doc__
diff --git a/python/gyoto/util.py b/python/gyoto/util.py
index 8f557e4f..51e448f3 100644
--- a/python/gyoto/util.py
+++ b/python/gyoto/util.py
@@ -307,3 +307,19 @@ def readScenery(filename):
 def writeObject(obj, filename):
     '''Write Gyoto object (e.g. Scenery) to XML file'''
     core.Factory(obj).write(filename)
+
+### Pythonic extension of methods
+# Worldline
+# This version of  getCartesian accepts numpy arrays
+#   wl.getCartesian(t, x, y, z, [xprime, yprime, zprime])
+# in addition to to the raw level arrays and dimension
+def _Worldline_getCartesian(self, t, *arrays):
+    if isinstance(t, core.array_double):
+        core._core.Worldline_getCartesian(self, t, *arrays)
+    else:
+        core._core.Worldline_getCartesian(
+            self,
+            core.array_double_fromnumpy1(t),
+            t.size,
+            *[core.array_double_fromnumpy1(v) for v in arrays]
+            )
diff --git a/python/gyoto_std.i b/python/gyoto_std.i
index 679a054f..57a15263 100644
--- a/python/gyoto_std.i
+++ b/python/gyoto_std.i
@@ -40,6 +40,16 @@ using namespace Gyoto;
 #endif
  }
 
+// Worldline::getCartesian is enough
+%ignore Gyoto::Astrobj::Star::getCartesian(double const * const dates,
+                                           size_t const n_dates,
+                                           double * const x,
+                                           double * const y,
+                                           double * const z,
+                                           double * const xprime=NULL,
+                                           double * const yprime=NULL,
+                                           double * const zprime=NULL);
+
 // Typemaps to translate SmartPointer to specific classes
 GyotoSmPtrTypeMapClassDerived(Astrobj, UniformSphere)
 GyotoSmPtrTypeMapClassDerived(Astrobj, Star)

From 0dda0c96b3d35725f8221a40305c727128480e5c Mon Sep 17 00:00:00 2001
From: Thibaut Paumard <thibaut.paumard@obspm.fr>
Date: Fri, 16 Aug 2019 09:32:15 +0200
Subject: [PATCH 12/17] Python: make sure Worldline.getCoord works both with
 NumPy and C arrays

---
 python/gyoto.i           | 31 +------------------------------
 python/gyoto/__init__.py |  2 ++
 python/gyoto/util.py     | 18 ++++++++++++++++++
 python/gyoto_std.i       |  2 +-
 4 files changed, 22 insertions(+), 31 deletions(-)

diff --git a/python/gyoto.i b/python/gyoto.i
index 5db81725..ffdff5da 100644
--- a/python/gyoto.i
+++ b/python/gyoto.i
@@ -1,5 +1,5 @@
 /*
-    Copyright 2014-2018 Thibaut Paumard
+    Copyright 2014-2019 Thibaut Paumard
 
     This file is part of Gyoto.
 
@@ -587,40 +587,11 @@ ExtendArrayNumPy(array_size_t, size_t);
 %rename(Worldline__IntegState__Generic) Gyoto::Worldline::IntegState::Generic;
 %rename(Worldline__IntegState__Boost) Gyoto::Worldline::IntegState::Boost;
 %rename(Worldline__IntegState__Legacy) Gyoto::Worldline::IntegState::Legacy;
-%ignore Gyoto::Worldline::getCoord(double const * const dates, size_t const n_dates,
-                                   double * const x1dest,
-                                   double * const x2dest, double * const x3dest,
-                                   double * const x0dot=NULL,  double * const x1dot=NULL,
-                                   double * const x2dot=NULL,  double * const x3dot=NULL) ;
 %extend Gyoto::Worldline {
   void get_t(double * INPLACE_ARRAY1, size_t DIM1) {
     if (DIM1 != ($self)->get_nelements()) GYOTO_ERROR("wrong output array size");
     ($self)->get_t(INPLACE_ARRAY1);
   }
-  void getCoord(double * dates, size_t n_dates,
-		double * x1dest, size_t n1,
-		double * x2dest, size_t n2,
-                double * x3dest, size_t n3,
-		double * x0dot, size_t n0d,
-                double * x1dot, size_t n1d,
-		double * x2dot, size_t n2d,
-                double * x3dot, size_t n3d) {
-    if (n1 != n_dates || n2 != n_dates || n3 != n_dates ||
-        (x0dot && n0d != n_dates) ||
-        (x1dot && n1d != n_dates) ||
-        (x2dot && n2d != n_dates) ||
-        (x3dot && n3d != n_dates))
-      GYOTO_ERROR("wrong size for output array");
-    ($self)->getCoord(dates, n_dates, x1dest, x2dest, x3dest, x0dot, x1dot, x2dot, x3dot);
-  }
-  void getCoord(double * dates, size_t n_dates,
-		double * x1dest, size_t n1,
-		double * x2dest, size_t n2,
-                double * x3dest, size_t n3) {
-    if (n_dates != ($self)->get_nelements() || n1 != n_dates || n2 != n_dates || n3 != n_dates)
-      GYOTO_ERROR("wrong size for output array");
-    ($self)->getCoord(dates, x1dest, x2dest, x3dest);
-  }
   void get_xyz(
 		double * x1dest, size_t n1,
 		double * x2dest, size_t n2,
diff --git a/python/gyoto/__init__.py b/python/gyoto/__init__.py
index b2c8e3a8..7bc490bd 100644
--- a/python/gyoto/__init__.py
+++ b/python/gyoto/__init__.py
@@ -17,3 +17,5 @@
 
 core.Worldline.getCartesian =  util._Worldline_getCartesian
 core.Worldline.getCartesian.__doc__ = core._core.Worldline_getCartesian.__doc__
+core.Worldline.getCoord =  util._Worldline_getCoord
+core.Worldline.getCoord.__doc__ = core._core.Worldline_getCoord.__doc__
diff --git a/python/gyoto/util.py b/python/gyoto/util.py
index 51e448f3..fee4b5f1 100644
--- a/python/gyoto/util.py
+++ b/python/gyoto/util.py
@@ -317,9 +317,27 @@ def _Worldline_getCartesian(self, t, *arrays):
     if isinstance(t, core.array_double):
         core._core.Worldline_getCartesian(self, t, *arrays)
     else:
+        sizes=numpy.asarray([v.size for v in arrays])
+        if not (sizes == t.size).all():
+            raise ValueError('all arrays must be the same size')
         core._core.Worldline_getCartesian(
             self,
             core.array_double_fromnumpy1(t),
             t.size,
             *[core.array_double_fromnumpy1(v) for v in arrays]
             )
+
+# Same for getCoord
+def _Worldline_getCoord(self, t, *arrays):
+    if isinstance(t, core.array_double):
+        core._core.Worldline_getCoord(self, t, *arrays)
+    else:
+        sizes=numpy.asarray([v.size for v in arrays])
+        if not (sizes == t.size).all():
+            raise ValueError('all arrays must be the same size')
+        core._core.Worldline_getCoord(
+            self,
+            core.array_double_fromnumpy1(t),
+            t.size,
+            *[core.array_double_fromnumpy1(v) for v in arrays]
+            )
diff --git a/python/gyoto_std.i b/python/gyoto_std.i
index 57a15263..18bf5962 100644
--- a/python/gyoto_std.i
+++ b/python/gyoto_std.i
@@ -1,5 +1,5 @@
 /*
-    Copyright 2014-2018 Thibaut Paumard
+    Copyright 2014-2019 Thibaut Paumard
 
     This file is part of Gyoto.
 

From f277f6ebf3cc3810c7b56bd0b779540723272545 Mon Sep 17 00:00:00 2001
From: Thibaut Paumard <thibaut.paumard@obspm.fr>
Date: Sun, 25 Aug 2019 19:32:33 +0200
Subject: [PATCH 13/17] gyoto.util: make getCoord(double t, stat_t coord) work
 again

---
 python/gyoto/util.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/python/gyoto/util.py b/python/gyoto/util.py
index fee4b5f1..6dbbd24e 100644
--- a/python/gyoto/util.py
+++ b/python/gyoto/util.py
@@ -329,7 +329,7 @@ def _Worldline_getCartesian(self, t, *arrays):
 
 # Same for getCoord
 def _Worldline_getCoord(self, t, *arrays):
-    if isinstance(t, core.array_double):
+    if isinstance(t, core.array_double) or numpy.isscalar(t) :
         core._core.Worldline_getCoord(self, t, *arrays)
     else:
         sizes=numpy.asarray([v.size for v in arrays])

From a60e859cfc4ecbf7768b519c7b59e565ac5676c8 Mon Sep 17 00:00:00 2001
From: Thibaut Paumard <thibaut.paumard@obspm.fr>
Date: Fri, 27 Sep 2019 14:37:09 +0200
Subject: [PATCH 14/17] Fix example.py after changing Worldline.getCoord

---
 python/example.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/python/example.py b/python/example.py
index da5507e6..96a08d31 100644
--- a/python/example.py
+++ b/python/example.py
@@ -61,6 +61,7 @@
 phi=numpy.ndarray(n)
 
 # Call Gyoto method that takes these arrays as argument:
+ph.get_t(t)
 ph.getCoord(t, r, theta, phi)
 
 plt.plot(t, r)

From 3ee9b4e25f26c65ab162d64c61a5bcd66c86a78f Mon Sep 17 00:00:00 2001
From: Thibaut Paumard <thibaut.paumard@obspm.fr>
Date: Tue, 10 Sep 2019 10:25:36 +0200
Subject: [PATCH 15/17] GyotoDefs.h: give derived constants to machine
 precision

Some constants like GYOTO_G_OVER_C_SQUARED where only given to the same
precision as the constant they derive from (GYOTO_G in this case).

While it is correct that e.g. G is known only to 1e5, mixiin guses of GYOTO_G
and GYOTO_G_OVER_C_SQUARED both specified only to that precision leads to
additional errors. GYOTO_G_OVER_C_SQUARE, GYOTO_G_OVER_C_SQUARE_CGS,
GYOTO_PLANCK_OVER_C_SQUARE and GYOTO_PLANCK_OVER_BOLTZMANN are now given with
many more places so that GYOTO_G/GYOTO_G_OVER_C_SQUARE actually yields GYOTO_C^2.
---
 include/GyotoDefs.h      | 8 ++++----
 yorick/gyoto_constants.i | 2 +-
 2 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/include/GyotoDefs.h b/include/GyotoDefs.h
index c1718cd1..64c47e96 100644
--- a/include/GyotoDefs.h
+++ b/include/GyotoDefs.h
@@ -429,15 +429,15 @@ namespace Gyoto {
 /// \brief Gravitational constant (cgs: cm^3 * g^-1 * s-2)
 #define GYOTO_G_CGS 6.67428e-8
 /// \brief G/c^2=6.67428e-11/299792458.^2
-#define GYOTO_G_OVER_C_SQUARE 7.426138e-28
+#define GYOTO_G_OVER_C_SQUARE 7.4261380161175445989e-28
 /// \brief G/c^2 in cgs units
-#define GYOTO_G_OVER_C_SQUARE_CGS 7.426138e-29
+#define GYOTO_G_OVER_C_SQUARE_CGS 7.4261380161175445989e-29
 /// \brief Planck's constant (h) in SI (J.s=kg.m^2/s) 
 #define GYOTO_PLANCK 6.62606896e-34
 /// \brief Planck's constant (h) in c.g.s (g.cm^2/s) 
 #define GYOTO_PLANCK_CGS 6.62606896e-27
 /// \brief h/c^2 in SI (kg.s)
-#define GYOTO_PLANCK_OVER_C_SQUARE 7.372496e-51
+#define GYOTO_PLANCK_OVER_C_SQUARE 7.3724959997591407964e-51
 /// \brief Boltzmann's constant (k) in SI (J/K)
 #define GYOTO_BOLTZMANN 1.3806504e-23
 /// \brief Boltzmann's constant (k) in cgs (erg/K)
@@ -445,7 +445,7 @@ namespace Gyoto {
 /// \brief Stefan-Boltzmann's constant (sigma) in cgs (erg/cm2/s/K4)
 #define GYOTO_STEFANBOLTZMANN_CGS 5.670373e-5
 /// \brief h/k (K.s = K/Hz)
-#define GYOTO_PLANCK_OVER_BOLTZMANN 4.7992373e-11
+#define GYOTO_PLANCK_OVER_BOLTZMANN 4.7992373449498869688e-11
 /// \brief ideal gas constant R in SI
 #define GYOTO_GAS_CST 8.3144621
 /// \brief ideal gas constant R in erg/(K mol)
diff --git a/yorick/gyoto_constants.i b/yorick/gyoto_constants.i
index 54a09b9e..d44af163 100644
--- a/yorick/gyoto_constants.i
+++ b/yorick/gyoto_constants.i
@@ -20,7 +20,7 @@
 extern GYOTO_C, GYOTO_G, GYOTO_G_OVER_C_SQUARE, GYOTO_SUN_MASS, GYOTO_SUN_RADIUS, GYOTO_KPC;
 GYOTO_C               = 299792458.
 GYOTO_G               = 6.67428e-11
-GYOTO_G_OVER_C_SQUARE = 7.426138e-28
+GYOTO_G_OVER_C_SQUARE = 7.4261380161175445989e-28
 GYOTO_SUN_MASS        = 1.98843e30
 GYOTO_SUN_RADIUS      = 6.955e8
 GYOTO_KPC             = 3.08568025e19

From 703ead1f392b1ff2c92fbba04d566c4f5b6dd2da Mon Sep 17 00:00:00 2001
From: Thibaut Paumard <thibaut.paumard@obspm.fr>
Date: Fri, 27 Sep 2019 16:02:00 +0200
Subject: [PATCH 16/17] Fix copy constructor of RezzollaZhidenko and
 PolishDoughnut

---
 lib/PolishDoughnut.C   | 3 +++
 lib/RezzollaZhidenko.C | 2 +-
 2 files changed, 4 insertions(+), 1 deletion(-)

diff --git a/lib/PolishDoughnut.C b/lib/PolishDoughnut.C
index 51876ff0..e7c14fff 100644
--- a/lib/PolishDoughnut.C
+++ b/lib/PolishDoughnut.C
@@ -136,8 +136,11 @@ PolishDoughnut::PolishDoughnut(const PolishDoughnut& orig) :
   central_temperature_(orig.central_temperature_),
   beta_(orig.beta_),
   magnetizationParameter_(orig.magnetizationParameter_),
+  aa_(orig.aa_),
+  aa2_(orig.aa2_),
   spectral_oversampling_(orig.spectral_oversampling_),
   angle_averaged_(orig.angle_averaged_),
+  bremsstrahlung_(orig.bremsstrahlung_),
   deltaPL_(orig.deltaPL_),
   adaf_(orig.adaf_),
   ADAFtemperature_(orig.ADAFtemperature_),
diff --git a/lib/RezzollaZhidenko.C b/lib/RezzollaZhidenko.C
index d14ceb9c..8a56338f 100644
--- a/lib/RezzollaZhidenko.C
+++ b/lib/RezzollaZhidenko.C
@@ -99,7 +99,7 @@ Gyoto::Metric::RezzollaZhidenko::RezzollaZhidenko()
 }
 
 Gyoto::Metric::RezzollaZhidenko::RezzollaZhidenko(const RezzollaZhidenko & orig)
-  : Generic(GYOTO_COORDKIND_SPHERICAL, "RezzollaZhidenko"),
+  : Generic(orig),
     epsilon_(orig.epsilon_), rms_(orig.rms_), rmb_(orig.rms_), aparam_(NULL), bparam_(NULL)
 {
   GYOTO_DEBUG << endl;

From 904c849c86fd9552837d1bd55efad9fb5b88cd06 Mon Sep 17 00:00:00 2001
From: Thibaut Paumard <thibaut.paumard@obspm.fr>
Date: Fri, 27 Sep 2019 16:04:19 +0200
Subject: [PATCH 17/17] Releasing 1.3.6

---
 ChangeLog    | 129 ++++++++++++++++++++++++++++++++++++++++++++++++++-
 NEWS         |   4 ++
 configure    |  20 ++++----
 configure.ac |   2 +-
 4 files changed, 143 insertions(+), 12 deletions(-)

diff --git a/ChangeLog b/ChangeLog
index d0a3e095..995bffe8 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,6 +1,133 @@
 commit unknown
 Author: Thibaut Paumard <thibaut.paumard@obspm.fr>
-Date:   Fri Aug 3 09:10:24 2019 +0200
+Date:   Fri Sep 27 16:03:01 2019 +0200
+
+    Releasing 1.3.6
+
+commit 703ead1f392b1ff2c92fbba04d566c4f5b6dd2da
+Author: Thibaut Paumard <thibaut.paumard@obspm.fr>
+Date:   Fri Sep 27 16:02:00 2019 +0200
+
+    Fix copy constructor of RezzollaZhidenko and PolishDoughnut
+
+commit 3ee9b4e25f26c65ab162d64c61a5bcd66c86a78f
+Author: Thibaut Paumard <thibaut.paumard@obspm.fr>
+Date:   Tue Sep 10 10:25:36 2019 +0200
+
+    GyotoDefs.h: give derived constants to machine precision
+    
+    Some constants like GYOTO_G_OVER_C_SQUARED where only given to the same
+    precision as the constant they derive from (GYOTO_G in this case).
+    
+    While it is correct that e.g. G is known only to 1e5, mixiin guses of GYOTO_G
+    and GYOTO_G_OVER_C_SQUARED both specified only to that precision leads to
+    additional errors. GYOTO_G_OVER_C_SQUARE, GYOTO_G_OVER_C_SQUARE_CGS,
+    GYOTO_PLANCK_OVER_C_SQUARE and GYOTO_PLANCK_OVER_BOLTZMANN are now given with
+    many more places so that GYOTO_G/GYOTO_G_OVER_C_SQUARE actually yields GYOTO_C^2.
+
+commit a60e859cfc4ecbf7768b519c7b59e565ac5676c8
+Author: Thibaut Paumard <thibaut.paumard@obspm.fr>
+Date:   Fri Sep 27 14:37:09 2019 +0200
+
+    Fix example.py after changing Worldline.getCoord
+
+commit f277f6ebf3cc3810c7b56bd0b779540723272545
+Author: Thibaut Paumard <thibaut.paumard@obspm.fr>
+Date:   Sun Aug 25 19:32:33 2019 +0200
+
+    gyoto.util: make getCoord(double t, stat_t coord) work again
+
+commit 0dda0c96b3d35725f8221a40305c727128480e5c
+Author: Thibaut Paumard <thibaut.paumard@obspm.fr>
+Date:   Fri Aug 16 09:32:15 2019 +0200
+
+    Python: make sure Worldline.getCoord works both with NumPy and C arrays
+
+commit 3277ebedd6d095a6f85c928ced59a680f8e90337
+Author: Thibaut Paumard <thibaut.paumard@obspm.fr>
+Date:   Thu Aug 15 22:35:51 2019 +0200
+
+    Python: make sure Worldline.getCartesian works both with NumPy and C arrays
+
+commit 11b1a6016db823a39d42f993432a8271aee86ede
+Author: Thibaut Paumard <thibaut.paumard@obspm.fr>
+Date:   Thu Aug 15 09:01:36 2019 +0200
+
+    gyoto.util.rayTrace: return arrays or the right dimension
+    
+    if j and/or i are scalars, decrease dimension of output arrays
+
+commit e1209889d5511fe9061f7c25155578a44f76e73c
+Author: Thibaut Paumard <thibaut.paumard@obspm.fr>
+Date:   Wed Aug 14 21:22:25 2019 +0200
+
+    gyoto.util.rayTrace: reorder parameters
+    
+    The order j, i is more Pythonic than i.j. This ways,
+      sc[j, i]['Quantity'] == sc[:,:]['Quantity'][j,i].
+    
+    Having height and width last allows writing:
+      sc[j, i, core.Bucket, None]
+    skipping height and width to remove progress output.
+
+commit c4af46c430bb0cd73c6230571d3530146b0f126d
+Author: Thibaut Paumard <thibaut.paumard@obspm.fr>
+Date:   Wed Aug 14 14:31:42 2019 +0200
+
+    Gyoto::Screen::Coord1dSet: copy beffur in Indices and Angles
+    
+    Much safer in Python
+
+commit 6614509bccc4937123963c96011af8955fa04a51
+Author: Thibaut Paumard <thibaut.paumard@obspm.fr>
+Date:   Sun Aug 11 14:07:43 2019 +0200
+
+    gyoto.util: improve doc comments
+
+commit f92a846193e3bd3cfa2adc5a04371d82b5158156
+Author: Thibaut Paumard <thibaut.paumard@obspm.fr>
+Date:   Sun Aug 11 14:00:19 2019 +0200
+
+    gyoto.util.rayTrace/Scenery_getitem: improve dimensionality
+    
+    - Scenery_getitem accepts all rayTrace parameters
+    - return 1D array when the Coord2dSet is 1D
+
+commit 18d2c3dfd3cd53e228b8ed0c91572a6232c07d04
+Author: Thibaut Paumard <thibaut.paumard@obspm.fr>
+Date:   Sun Aug 11 07:53:12 2019 +0200
+
+    gyoto.util: support angles in rayTrace
+
+commit 7852e47d204d4b6f87032abc36179a4cb23f479b
+Author: Thibaut Paumard <thibaut.paumard@obspm.fr>
+Date:   Sun Aug 11 07:45:42 2019 +0200
+
+    gyoto.core.Scenery: add __getitem__ method
+    
+    shortcut for rayTrace
+
+commit 7d86c285b0b9fb095bbcb102fb39fdc0f96e076d
+Author: Thibaut Paumard <thibaut.paumard@obspm.fr>
+Date:   Sun Aug 11 07:21:59 2019 +0200
+
+    gyoto.util.rayTrace: support indices, ranges and list of indices
+
+commit a0f2db907499a768e9d0c8322bd87386d5326b0a
+Author: Thibaut Paumard <thibaut.paumard@obspm.fr>
+Date:   Sun Aug 11 06:21:56 2019 +0200
+
+    gyoto.util.rayTrace: accept paremeters i and j
+
+commit 049e1daad6fb52b9beb6be40f9778c733b9755c1
+Author: Thibaut Paumard <thibaut.paumard@obspm.fr>
+Date:   Fri Aug 9 09:37:16 2019 +0200
+
+    python: install util.rayTrace as Scenery.rayTrace
+
+commit bb7ce60adc3c9f476dafb0429b35df57ae99fe0a
+Author: Thibaut Paumard <thibaut.paumard@obspm.fr>
+Date:   Sat Aug 3 09:13:51 2019 +0200
 
     Releasing 1.3.5
 
diff --git a/NEWS b/NEWS
index 6ad66970..64367cac 100644
--- a/NEWS
+++ b/NEWS
@@ -1,3 +1,7 @@
+1.3.6 2019/09/27 FTR
+   * Python: improve interface in gyoto.util
+   * GyotoDefs.h: give derived constants to machine precision
+
 1.3.5 2019/08/03 BUG
    * Python: support Python2.7 in gyoto.util.rayTrace
 
diff --git a/configure b/configure
index fc77c4de..b5e24f0a 100755
--- a/configure
+++ b/configure
@@ -1,6 +1,6 @@
 #! /bin/sh
 # Guess values for system-dependent variables and create Makefiles.
-# Generated by GNU Autoconf 2.69 for Gyoto 1.3.5.
+# Generated by GNU Autoconf 2.69 for Gyoto 1.3.6.
 #
 # Report bugs to <gyoto@sympa.obspm.fr>.
 #
@@ -590,8 +590,8 @@ MAKEFLAGS=
 # Identity of this package.
 PACKAGE_NAME='Gyoto'
 PACKAGE_TARNAME='gyoto'
-PACKAGE_VERSION='1.3.5'
-PACKAGE_STRING='Gyoto 1.3.5'
+PACKAGE_VERSION='1.3.6'
+PACKAGE_STRING='Gyoto 1.3.6'
 PACKAGE_BUGREPORT='gyoto@sympa.obspm.fr'
 PACKAGE_URL=''
 
@@ -1488,7 +1488,7 @@ if test "$ac_init_help" = "long"; then
   # Omit some internal or obsolete options to make the list less imposing.
   # This message is too long to be a string in the A/UX 3.1 sh.
   cat <<_ACEOF
-\`configure' configures Gyoto 1.3.5 to adapt to many kinds of systems.
+\`configure' configures Gyoto 1.3.6 to adapt to many kinds of systems.
 
 Usage: $0 [OPTION]... [VAR=VALUE]...
 
@@ -1560,7 +1560,7 @@ fi
 
 if test -n "$ac_init_help"; then
   case $ac_init_help in
-     short | recursive ) echo "Configuration of Gyoto 1.3.5:";;
+     short | recursive ) echo "Configuration of Gyoto 1.3.6:";;
    esac
   cat <<\_ACEOF
 
@@ -1809,7 +1809,7 @@ fi
 test -n "$ac_init_help" && exit $ac_status
 if $ac_init_version; then
   cat <<\_ACEOF
-Gyoto configure 1.3.5
+Gyoto configure 1.3.6
 generated by GNU Autoconf 2.69
 
 Copyright (C) 2012 Free Software Foundation, Inc.
@@ -2508,7 +2508,7 @@ cat >config.log <<_ACEOF
 This file contains any messages produced by compilers while
 running configure, to aid debugging if configure makes a mistake.
 
-It was created by Gyoto $as_me 1.3.5, which was
+It was created by Gyoto $as_me 1.3.6, which was
 generated by GNU Autoconf 2.69.  Invocation command line was
 
   $ $0 $@
@@ -3484,7 +3484,7 @@ fi
 
 # Define the identity of the package.
  PACKAGE='gyoto'
- VERSION='1.3.5'
+ VERSION='1.3.6'
 
 
 cat >>confdefs.h <<_ACEOF
@@ -21955,7 +21955,7 @@ cat >>$CONFIG_STATUS <<\_ACEOF || ac_write_fail=1
 # report actual input values of CONFIG_FILES etc. instead of their
 # values after options handling.
 ac_log="
-This file was extended by Gyoto $as_me 1.3.5, which was
+This file was extended by Gyoto $as_me 1.3.6, which was
 generated by GNU Autoconf 2.69.  Invocation command line was
 
   CONFIG_FILES    = $CONFIG_FILES
@@ -22021,7 +22021,7 @@ _ACEOF
 cat >>$CONFIG_STATUS <<_ACEOF || ac_write_fail=1
 ac_cs_config="`$as_echo "$ac_configure_args" | sed 's/^ //; s/[\\""\`\$]/\\\\&/g'`"
 ac_cs_version="\\
-Gyoto config.status 1.3.5
+Gyoto config.status 1.3.6
 configured by $0, generated by GNU Autoconf 2.69,
   with options \\"\$ac_cs_config\\"
 
diff --git a/configure.ac b/configure.ac
index 5caf995e..63895f05 100644
--- a/configure.ac
+++ b/configure.ac
@@ -10,7 +10,7 @@ dnl without any warranty.
 
 
 AC_PREREQ([2.69])
-AC_INIT([Gyoto], [1.3.5], [gyoto@sympa.obspm.fr])
+AC_INIT([Gyoto], [1.3.6], [gyoto@sympa.obspm.fr])
 gyoto_test_CXXFLAGS=${CXXFLAGS+set}
 gyoto_test_CPPFLAGS=${CPPFLAGS+set}
 AC_CANONICAL_SYSTEM