From cc4108e7763bee277337f8b315a92d34f21b3d69 Mon Sep 17 00:00:00 2001
From: Ben Mather <brmather1@gmail.com>
Date: Wed, 16 Aug 2023 14:57:53 +1000
Subject: [PATCH] Fix bug in motion paths

---
 gplately/reconstruction.py | 33 ++++++++++++++++++++++-----------
 1 file changed, 22 insertions(+), 11 deletions(-)

diff --git a/gplately/reconstruction.py b/gplately/reconstruction.py
index 70ecda5f..6ddcabfa 100644
--- a/gplately/reconstruction.py
+++ b/gplately/reconstruction.py
@@ -743,10 +743,14 @@ def create_motion_path(self, lons, lats, time_array, plate_id=None, anchor_plate
         # rates of motion (if requested)
         rlons = np.empty((len(time_array), len(lons)))
         rlats = np.empty((len(time_array), len(lons)))
-        StepTimes = np.empty(((len(time_array)-1)*2, len(lons)))
-        StepRates = np.empty(((len(time_array)-1)*2, len(lons)))
+        
+        if plate_id is None:
+            query_plate_id = True
+        else:
+            query_plate_id = False
+            plate_ids = np.ones(len(lons), dtype=int) * plate_id
 
-        seed_points = list(zip(lats,lons))
+        seed_points = zip(lats,lons)
         for i, lat_lon in enumerate(seed_points):
 
             seed_points_at_digitisation_time = pygplates.PointOnSphere(
@@ -754,37 +758,43 @@ def create_motion_path(self, lons, lats, time_array, plate_id=None, anchor_plate
             )
             # Allocate the present-day plate ID to the PointOnSphere if 
             # it was not given.
-            if plate_id is None:
+            if query_plate_id:
                 plate_id = _tools.plate_partitioner_for_point(
                     lat_lon,
                     self.topology_features,
                     self.rotation_model
                 )
+            else:
+                plate_id = plate_ids[i]
+                
             # Create the motion path feature. enforce float and int for C++ signature.
             motion_path_feature = pygplates.Feature.create_motion_path(
                 seed_points_at_digitisation_time, 
                 time_array, 
                 valid_time=(time_array.max(), time_array.min()),
-                relative_plate=int(anchor_plate_id),
-                reconstruction_plate_id=int(plate_id))
+                relative_plate=int(plate_id),
+                reconstruction_plate_id=int(anchor_plate_id))
 
             reconstructed_motion_paths = self.reconstruct(
                 motion_path_feature, 
                 to_time=0, 
                 reconstruct_type=pygplates.ReconstructType.motion_path,
-                anchor_plate_id=int(anchor_plate_id),
+                anchor_plate_id=int(plate_id),
             )
             # Turn motion paths in to lat-lon coordinates
             for reconstructed_motion_path in reconstructed_motion_paths:
                 trail = reconstructed_motion_path.get_motion_path().to_lat_lon_array()
 
             lon, lat = np.flipud(trail[:, 1]), np.flipud(trail[:, 0])
-
+            
             rlons[:,i] = lon
             rlats[:,i] = lat
 
             # Obtain step-plot coordinates for rate of motion
             if return_rate_of_motion is True:
+
+                StepTimes = np.empty(((len(time_array)-1)*2, len(lons)))
+                StepRates = np.empty(((len(time_array)-1)*2, len(lons)))
                 
                 # Get timestep
                 TimeStep = []
@@ -1490,14 +1500,15 @@ def motion_path(self, time_array, anchor_plate_id=0, return_rate_of_motion=False
                 valid_time=(time_array.max(), time_array.min()),
                 #relative_plate=int(self.plate_id[i]),
                 #reconstruction_plate_id=int(anchor_plate_id))
-                relative_plate=int(anchor_plate_id),
-                reconstruction_plate_id=int(self.plate_id[i]))
+                relative_plate=int(self.plate_id[i]),
+                reconstruction_plate_id=int(anchor_plate_id))
 
             reconstructed_motion_paths = self.plate_reconstruction.reconstruct(
                 motion_path_feature, 
                 to_time=0, 
                 #from_time=0, 
-                reconstruct_type=pygplates.ReconstructType.motion_path)
+                reconstruct_type=pygplates.ReconstructType.motion_path,
+                anchor_plate_id=int(anchor_plate_id))
 
             # Turn motion paths in to lat-lon coordinates
             for reconstructed_motion_path in reconstructed_motion_paths: