From b546ba6eaa57c7308ed08ca62683800d17b6db9e Mon Sep 17 00:00:00 2001
From: David A Minton <daminton@purdue.edu>
Date: Mon, 9 Jan 2023 08:35:42 -0500
Subject: [PATCH] Added Merger origin_type

---
 examples/Chambers2013/aescattermovie.py | 50 +++++++++++--------------
 src/collision/collision_resolve.f90     |  6 ++-
 2 files changed, 27 insertions(+), 29 deletions(-)

diff --git a/examples/Chambers2013/aescattermovie.py b/examples/Chambers2013/aescattermovie.py
index 1695edd90..40a8e7202 100755
--- a/examples/Chambers2013/aescattermovie.py
+++ b/examples/Chambers2013/aescattermovie.py
@@ -7,13 +7,15 @@
 
 titletext = "Chambers (2013)"
 radscale = 2000
+AU = 1.0
 xmin = 0.0
-xmax = 2.20
-ymin = 0.0 
+xmax = 2.25
+ymin = 0.0
 ymax = 1.0
-framejump = 1
+framejump = 10
+ncutoff = 1e20
 
-class AnimatedScatter():
+class AnimatedScatter(object):
     """An animated scatter plot using matplotlib.animations.FuncAnimation."""
     def __init__(self, ds, param):
 
@@ -23,7 +25,8 @@ def __init__(self, ds, param):
         self.param = param
         self.Rcb = self.ds['radius'].sel(name="Sun").isel(time=0).values[()]
 
-        self.clist = {'Initial conditions' : 'xkcd:faded blue',
+        self.clist = {'Initial conditions' : 'k',
+                      'Merger' : 'xkcd:faded blue',
                       'Disruption' : 'xkcd:marigold',
                       'Supercatastrophic' : 'xkcd:shocking pink',
                       'Hit and run fragmentation' : 'xkcd:baby poop green'}
@@ -37,7 +40,7 @@ def __init__(self, ds, param):
         self.ax.set_ylim(ymin, ymax)
         fig.add_axes(self.ax)
         self.ani = animation.FuncAnimation(fig, self.update, interval=1, frames=nframes, init_func=self.setup_plot, blit=True)
-        self.ani.save('aescatter.mp4', fps=60, dpi=300, extra_args=['-vcodec', 'libx264'])
+        self.ani.save('aescatter.mp4', fps=30, dpi=300, extra_args=['-vcodec', 'libx264'])
         print('Finished writing aescattter.mp4')
 
     def scatters(self, pl, radmarker, origin):
@@ -58,26 +61,22 @@ def setup_plot(self):
         self.ax.set_xlabel("Semimajor Axis (AU)", fontsize='16', labelpad=1)
         self.ax.set_ylabel("Eccentricity", fontsize='16', labelpad=1)
 
-        self.title = self.ax.text(0.50, 1.05, "", bbox={'facecolor': 'w', 'alpha': 0.5, 'pad': 5}, transform=self.ax.transAxes,
+        title = self.ax.text(0.50, 1.05, "", bbox={'facecolor': 'w', 'alpha': 0.5, 'pad': 5}, transform=self.ax.transAxes,
                         ha="center")
 
-        self.title.set_text(f"{titletext} -  Time = ${t*1e-6:6.2f}$ My with ${npl:4.0f}$ particles")
-        slist = self.scatters(pl, radmarker, origin)
-        self.s0 = slist[0]
-        self.s1 = slist[1]
-        self.s2 = slist[2]
-        self.s3 = slist[3]
+        title.set_text(f"{titletext} -  Time = ${t*1e-6:6.2f}$ My with ${npl:4.0f}$ particles")
+        self.slist = self.scatters(pl, radmarker, origin)
+        self.slist.append(title)
         leg = plt.legend(loc="upper right", scatterpoints=1, fontsize=10)
         for i,l in enumerate(leg.legendHandles):
            leg.legendHandles[i]._sizes = [20]
-        return self.s0, self.s1, self.s2, self.s3, self.title
+        return self.slist
 
     def data_stream(self, frame=0):
         while True:
             d = self.ds.isel(time = frame)
-            name_good = d.name.where(d['status'] != 1, drop=True)
-            name_good = name_good.where(name_good != "Sun", drop=True)
-            d = d.sel(name=name_good) 
+            name_good = d['name'].where(d['name'] != "Sun", drop=True)
+            d = d.sel(name=name_good)
             d['radmarker'] = (d['radius'] / self.Rcb) * radscale
             radius = d['radmarker'].values
 
@@ -98,22 +97,17 @@ def update(self,frame):
         """Update the scatter plot."""
         t, name, Gmass, radius, npl, pl, radmarker, origin = next(self.data_stream(framejump * frame))
 
-        self.title.set_text(f"{titletext} - Time = ${t*1e-6:6.3f}$ My with ${npl:4.0f}$ particles")
+        self.slist[-1].set_text(f"{titletext} - Time = ${t*1e-6:6.3f}$ My with ${npl:4.0f}$ particles")
 
         # We need to return the updated artist for FuncAnimation to draw..
         # Note that it expects a sequence of artists, thus the trailing comma.
-        s = [self.s0, self.s1, self.s2, self.s3]
         for i, (key, value) in enumerate(self.clist.items()):
             idx = origin == key
-            s[i].set_sizes(radmarker[idx])
-            s[i].set_offsets(pl[idx,:])
-            s[i].set_facecolor(value)
-
-        self.s0 = s[0]
-        self.s1 = s[1]
-        self.s2 = s[2]
-        self.s3 = s[3]
-        return self.s0, self.s1, self.s2, self.s3, self.title,
+            self.slist[i].set_sizes(radmarker[idx])
+            self.slist[i].set_offsets(pl[idx,:])
+            self.slist[i].set_facecolor(value)
+
+        return self.slist
 
 sim = swiftest.Simulation(read_old_output=True)
 print('Making animation')
diff --git a/src/collision/collision_resolve.f90 b/src/collision/collision_resolve.f90
index 4b88971cb..b2becbc90 100644
--- a/src/collision/collision_resolve.f90
+++ b/src/collision/collision_resolve.f90
@@ -409,7 +409,11 @@ module subroutine collision_resolve_mergeaddsub(nbody_system, param, t, status)
                end do 
             case(MERGED)
                call plnew%info(1)%copy(pl%info(ibiggest))
-               plnew%status(1) = OLD_PARTICLE
+               write(origin_type,*) "Merger"
+               plnew%status(1) = NEW_PARTICLE
+               call plnew%info(1)%set_value(origin_type=origin_type, origin_time=t,&
+                                            origin_rh=plnew%rh(:,1), origin_vh=plnew%vh(:,1), &
+                                            collision_id=param%maxid_collision)
                do i = 1, nimpactors
                   if (impactors%id(i) == ibiggest) cycle