Minor updates for devel compatibility

Binary/README.md

@@ -1,7 +1,25 @@
+# Binary star example
+
 A simple simulation of a binary star system, to exercise the code and give some experience with the file formats. The example takes two stars and evolves them under direct interaction. The initial conditions for the stars are specified in `binary.bods`, which follows standard `EXP` input format: the first line is the number of bodies (followed by two integers that specify the number of extra integer and double fields -- here set to be zero). Lines after the first line specify, for each body, the mass, xyz position, and xyz velocity.
 
+## How to run
+
 Run as follows.  The YAML configuration file is set for the Docker container. If you are using a native build and installation, change the ldlibdir parameter to point at your library install directory.  Then, execute the command: `mpirun -np 1 exp binary.yml`.
 
+## A fixed potential example
+
 There is also an example where the stars no longer feel mutual self-gravity, but instead are evolved in an external logarithmic potential. This may be run using `mpirun -np 1 exp logpot.yml`. This example demonstrates `noforce`.
 
+## Viewing results
+
+Each example creates a trajectory file named `ORBTRACE.runX` where `X`
+is either `1` or `2`, for `binary.yml` or `logpot.yml`
+respectively. These may be plotted using the command:
+
+``` bash
+python3 plotOrbit.py --input=ORBTRACE.runX
+```
+
+## Additional exploration
+
 A first exercise to extend these examples would be to instead use the `usermw` module to simulate the orbit of the sun in the MW.

Binary/logpot.yml

@@ -15,7 +15,7 @@ Global:
   multistep  : 12
   dynfracA   : 0.01
   dynfracV   : 0.01
-  VERBOSE    : 4
+  VERBOSE    : 1
 # ------------------------------------------------------------------------
 # This is a sequence of components.  The parameters for the force are
 # now included as a parameter map, rather than a separate file.

Binary/plotOrbit.py

@@ -0,0 +1,57 @@
+import matplotlib.pyplot as plt
+import numpy as np
+import argparse
+
+def plot_orbit(positions, title="Orbital Trajectory", save_path=None):
+    """
+    Plots two orbital trajectories given a list of 3D positions.
+
+    Parameters:
+    positions: The x, y, z positions for N times both trajectories as an ndarray of shape (N, 6).
+    title (str): The title of the plot.
+    save_path (str): If provided, the plot will be saved to this path.
+    """
+
+    # Extract x, y, z coordinates
+    x1 = positions[:, 0]
+    y1 = positions[:, 1]
+    z1 = positions[:, 2]
+    x2 = positions[:, 3]
+    y2 = positions[:, 4]
+    z2 = positions[:, 5]
+
+    # Create a 3D plot
+    fig = plt.figure()
+    ax = fig.add_subplot(111, projection='3d')
+
+    # Plot the trajectory
+    ax.plot(x1, y1, z1, '-')
+    ax.plot(x2, y2, z2, '-')
+
+    # Set labels and title
+    ax.set_xlabel('X Position')
+    ax.set_ylabel('Y Position')
+    ax.set_zlabel('Z Position')
+    ax.set_title(title)
+
+    # Show grid
+    ax.grid(True)
+
+    # Save or show the plot
+    if save_path:
+        plt.savefig(save_path)
+    else:
+        plt.show()
+
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(description="Plot orbital trajectory from 3D positions.")
+    parser.add_argument('--input', type=str,
+                        default='ORBTRACE.run0',
+                        help="Path to the input file containing 3D positions (one per line as x,y,z).")
+
+    args = parser.parse_args()
+
+    data = np.loadtxt(args.input)
+
+    plot_orbit(data[:, [1,2,3,11,12,13]], title="Binary orbit")

DiskHalo/config.yml

@@ -39,7 +39,7 @@ Components:
     bodyfile   : disk.bod
     force :
       id : cylinder
-      parameters : {acyl: 1.0, nmaxfid: 32, mmax: 2, hcyl: 0.05, nmax: 10, ncylodd: 2, ncylnx: 128, ncylny: 64, cachename: disk.cache}
+      parameters : {acyl: 1.0, lmaxfid: 64, mmax: 2, hcyl: 0.05, sech2: false, nmax: 10, ncylodd: 2, ncylnx: 200, ncylny: 100, cachename: disk.cache}
 #
 # The parameters could be expressed like this, equivalently:
 #

Nbody/README.md

@@ -48,3 +48,13 @@ testing.
 
 The directory `data` contains some snapshots and coefficients for use
 with the pyEXP tutorials.
+
+## Viewing results
+
+This example creates a trajectory file with selected disk orbits named
+`ORBTRACE.run0`.  These orbits may be plotted using the command:
+
+``` bash
+python3 plotOrbit.py
+```
+

Nbody/plotOrbit.py

@@ -0,0 +1,56 @@
+import matplotlib.pyplot as plt
+import numpy as np
+import argparse
+
+def plot_orbit(positions, title="Orbital Trajectory", save_path=None):
+    """
+    Plots two orbital trajectories given a list of 3D positions.
+
+    Parameters:
+    positions: A list of 3d particle trajectories (x, y, z) positions for N times, as a list of M ndarrays of shape (N, 3).
+    title (str): The title of the plot.
+    save_path (str): If provided, the plot will be saved to this path.
+    """
+
+    # Create a 3D plot
+    fig = plt.figure()
+    ax = fig.add_subplot(111, projection='3d')
+
+    # Plot the trajectories
+    for n in range(len(positions)):
+        ax.plot(positions[n][:, 0], positions[n][:, 1], positions[n][:, 2], '-')
+
+    # Set labels and title
+    ax.set_xlabel('X Position')
+    ax.set_ylabel('Y Position')
+    ax.set_zlabel('Z Position')
+    ax.set_title(title)
+
+    # Show grid
+    ax.grid(True)
+
+    # Save or show the plot
+    if save_path:
+        plt.savefig(save_path)
+    else:
+        plt.show()
+
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(description="Plot orbital trajectory from 3D positions.")
+    parser.add_argument('--input', type=str,
+                        default='ORBTRACE.run0',
+                        help="Path to the input file containing 3D positions (one per line as x,y,z).")
+
+    args = parser.parse_args()
+
+    data = np.loadtxt(args.input)
+
+    pos = []
+    pos.append(data[:, 1:4])    # First object positions
+    pos.append(data[:, 7:10])   # Second object positions
+    pos.append(data[:, 13:16])  # Third object positions
+    pos.append(data[:, 19:22])  # Fourth object positions
+    pos.append(data[:, 25:28])  # Fifth object positions
+
+    plot_orbit(pos, title="Selected disk orbits")