Init commit

.gitignore

@@ -0,0 +1,133 @@
+# Byte-compiled / optimized / DLL files
+__pycache__/
+*.py[cod]
+*$py.class
+
+# C extensions
+*.so
+
+# Distribution / packaging
+.Python
+build/
+develop-eggs/
+dist/
+downloads/
+eggs/
+.eggs/
+lib/
+lib64/
+parts/
+sdist/
+var/
+wheels/
+pip-wheel-metadata/
+share/python-wheels/
+*.egg-info/
+.installed.cfg
+*.egg
+MANIFEST
+
+# PyInstaller
+#  Usually these files are written by a python script from a template
+#  before PyInstaller builds the exe, so as to inject date/other infos into it.
+*.manifest
+*.spec
+
+# Installer logs
+pip-log.txt
+pip-delete-this-directory.txt
+
+# Unit test / coverage reports
+htmlcov/
+.tox/
+.nox/
+.coverage
+.coverage.*
+.cache
+nosetests.xml
+coverage.xml
+*.cover
+*.py,cover
+.hypothesis/
+.pytest_cache/
+
+# Translations
+*.mo
+*.pot
+
+# Django stuff:
+*.log
+local_settings.py
+db.sqlite3
+db.sqlite3-journal
+
+# Flask stuff:
+instance/
+.webassets-cache
+
+# Scrapy stuff:
+.scrapy
+
+# Sphinx documentation
+docs/_build/
+
+# PyBuilder
+target/
+
+# Jupyter Notebook
+.ipynb_checkpoints
+
+# IPython
+profile_default/
+ipython_config.py
+
+# pyenv
+.python-version
+
+# pipenv
+#   According to pypa/pipenv#598, it is recommended to include Pipfile.lock in version control.
+#   However, in case of collaboration, if having platform-specific dependencies or dependencies
+#   having no cross-platform support, pipenv may install dependencies that don't work, or not
+#   install all needed dependencies.
+#Pipfile.lock
+
+# PEP 582; used by e.g. github.com/David-OConnor/pyflow
+__pypackages__/
+
+# Celery stuff
+celerybeat-schedule
+celerybeat.pid
+
+# SageMath parsed files
+*.sage.py
+
+# Environments
+.env
+.venv
+env/
+venv/
+ENV/
+env.bak/
+venv.bak/
+
+# Spyder project settings
+.spyderproject
+.spyproject
+
+# Rope project settings
+.ropeproject
+
+# mkdocs documentation
+/site
+
+# mypy
+.mypy_cache/
+.dmypy.json
+dmypy.json
+
+# Pyre type checker
+.pyre/
+
+.idea/
+
+*.ipynb_checkpoints/

app/Disc Golf Simulator.py

@@ -0,0 +1,117 @@
+from logging import getLogger
+from pathlib import Path
+proj_dir = Path(__file__).parents[1]
+import sys
+sys.path.append(str(proj_dir))
+
+import numpy as np
+import streamlit as st
+
+from shotshaper.projectile import DiscGolfDisc
+from visualize import get_plot, get_stl, get_subplots, visualize_disc
+
+
+# Define the default values
+default_U = 24.2
+default_omega = 116.8
+default_z0 = 1.3
+default_pitch = 15.5
+default_nose = 0.0
+default_roll = 14.7
+
+
+def main():
+    tab1, tab2 = st.tabs(['Simulator', 'FAQ'])
+    with tab1:
+        disc_names = {
+            'Innova Wraith': 'dd2',
+            'Innova Firebird': 'cd1',
+            'Innova Roadrunner': 'cd5',
+            'Innova Fairway Driver': 'fd2',
+            }
+        disc_selected = st.sidebar.selectbox("Disc Selection", disc_names.keys())
+        disc_name = disc_names[disc_selected]
+
+        # Create the sliders with the default values
+        U = st.sidebar.slider("Throwing Velocity (m/s)", min_value=0.0, max_value=40.0, value=default_U, step=0.1,
+                              help='Fastest Throw on record is ~40m/s by Simon Lizotte')
+        omega = st.sidebar.slider("Omega", min_value=0.0, max_value=200.0, value=default_omega, step=0.1)
+        z0 = st.sidebar.slider("Release Height (m)", min_value=0.0, max_value=2.0, value=default_z0, step=0.1)
+        pitch = st.sidebar.slider("Pitch Angle (deg) | Release angle", min_value=0.0, max_value=90.0, value=default_pitch,
+                                  step=0.1)
+        nose = st.sidebar.slider("Nose Angle (deg) | Up/Down", min_value=0.0, max_value=90.0, value=default_nose,
+                                 step=0.1)
+        roll = st.sidebar.slider("Roll Angle (deg) | Tilt Left/Right", min_value=-90.0, max_value=90.0, value=default_roll,
+                                 step=0.1)
+
+        pos = np.array((0, 0, default_z0))
+        disc_dict = DiscGolfDisc(disc_name)
+
+        stl_mesh = get_stl(proj_dir / 'shotshaper' / 'discs' / (disc_name + '.stl'))
+        fig = visualize_disc(stl_mesh, nose=nose, roll=roll)
+
+        st.markdown("""## Disc orientation""")
+        st.plotly_chart(fig)
+        st.markdown("""## Flight Path""")
+        shot = disc_dict.shoot(speed=U, omega=omega, pitch=pitch,
+                               position=pos, nose_angle=nose, roll_angle=roll)
+
+        # Plot trajectory
+        x, y, z = shot.position
+        x_new, y_new = -1 * y, x
+
+        # Reversed x and y to mimic a throw
+        fig = get_plot(x_new, y_new, z)
+        st.plotly_chart(fig, True)
+
+        st.markdown(
+                f"""
+        **Arrows in Blue** show you where your *s-turn* is.
+
+        **Arrows in Red** show you your *max height* and *lateral deviance*.
+
+        Hit Play to watch your animated throw.
+
+        | Metric       | Value  |
+        |--------------|--------|
+        | Drift Left   | {round(min(x_new), 2)} |
+        | Drift Right  | {round(max(x_new), 2)} |
+        | Max Height   | {round(max(z), 2)}    |
+        | Distance     | {round(max(y_new), 2)} |
+
+        """
+                )
+
+        arc, alphas, betas, lifts, drags, moms, rolls = disc_dict.post_process(shot, omega)
+        fig = get_subplots(arc, alphas, lifts, drags, moms, rolls, shot.velocity)
+        st.plotly_chart(fig, True)
+
+    with tab2:
+        st.markdown("""
+        # Motivation
+        I saw some great work by [kegiljarhus](https://github.com/kegiljarhus) [repo](https://github.com/kegiljarhus/shotshaper)
+        and wanted to make this available as an app so people could learn more about disc golf. I really want to commend
+        the amazing idea of writing a [scientific article](https://link.springer.com/article/10.1007/s12283-022-00390-5)
+         AND releasing code, and actually executing it well. This is what gets people excited about STEM.
+        
+        I originally saw this 
+        [reddit post](https://www.reddit.com/r/discgolf/comments/yyhbcj/wrote_a_scientific_article_on_disc_golf_flight/)
+        which really piqued my interest.
+        
+        # Questions
+        - I imagine some of you will want to add your disc here, if you can convert your disc into an `.stl` then I will 
+        add it to the database. If this gets common enough I will add an option to upload your own.
+            - I imagine there will be a barrier to entry to do this.
+            - If you have any ideas, just let me know in a discussion or in a pull request
+        """)
+
+
+if __name__ == "__main__":
+    # Setting up Logger and proj_dir
+    logger = getLogger(__name__)
+    proj_dir = Path(__file__).parents[1]
+
+    st.title("Disc Golf Simulator")
+
+    # initialize_state()
+    main()

app/extrema.py

@@ -0,0 +1,32 @@
+import numpy as np
+from scipy.signal import argrelextrema
+
+def find_extrema(x, y):
+    # Find the indices of the local maxima and minima
+    maxima_idx = argrelextrema(x, np.greater)[0]
+    minima_idx = argrelextrema(x, np.less)[0]
+
+    # Get the x and y values at the extrema
+    maxima_x = x[maxima_idx]
+    maxima_y = y[maxima_idx]
+    minima_x = x[minima_idx]
+    minima_y = y[minima_idx]
+
+    # Combine the maxima and minima into a single array
+    extrema_x = np.concatenate((maxima_x, minima_x))
+    extrema_y = np.concatenate((maxima_y, minima_y))
+
+    # Determine whether each extrema is a maximum or minimum
+    is_maxima = np.zeros(len(extrema_x), dtype=bool)
+    is_maxima[:len(maxima_x)] = True
+
+    # Sort the extrema by x-value
+    idx = np.argsort(extrema_x)
+    extrema_x = extrema_x[idx]
+    extrema_y = extrema_y[idx]
+    is_maxima = is_maxima[idx]
+
+    # Convert the boolean array to a list of strings
+    extrema_type = ["arrow-left" if is_max else "arrow-right" for is_max in is_maxima]
+
+    return extrema_x, extrema_y, extrema_type

app/get_disc.py

@@ -0,0 +1,101 @@
+import requests
+
+headers = {
+    'authority': 'alldiscs.com',
+    'accept': 'application/json, text/javascript, */*; q=0.01',
+    'accept-language': 'en-US,en;q=0.6',
+    'content-type': 'application/x-www-form-urlencoded; charset=UTF-8',
+    'origin': 'https://alldiscs.com',
+    'referer': 'https://alldiscs.com/',
+    'sec-fetch-dest': 'empty',
+    'sec-fetch-mode': 'cors',
+    'sec-fetch-site': 'same-origin',
+    'sec-gpc': '1',
+    'user-agent': 'Mozilla/5.0 (Windows NT 10.0; Win64; x64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/107.0.0.0 Safari/537.36',
+    'x-requested-with': 'XMLHttpRequest',
+}
+
+params = {
+    'action': 'get_wdtable',
+    'table_id': '5',
+}
+
+data = {
+    'draw': '4',
+    'columns[0][data]': '0',
+    'columns[0][name]': 'wdt_ID',
+    'columns[0][searchable]': 'true',
+    'columns[0][orderable]': 'true',
+    'columns[0][search][value]': '',
+    'columns[0][search][regex]': 'false',
+    'columns[1][data]': '1',
+    'columns[1][name]': 'brand',
+    'columns[1][searchable]': 'true',
+    'columns[1][orderable]': 'true',
+    'columns[1][search][value]': '',
+    'columns[1][search][regex]': 'false',
+    'columns[2][data]': '2',
+    'columns[2][name]': 'mold',
+    'columns[2][searchable]': 'true',
+    'columns[2][orderable]': 'true',
+    'columns[2][search][value]': '',
+    'columns[2][search][regex]': 'false',
+    'columns[3][data]': '3',
+    'columns[3][name]': 'type',
+    'columns[3][searchable]': 'true',
+    'columns[3][orderable]': 'true',
+    'columns[3][search][value]': 'Distance|Fairway|Midrange|Putter',
+    'columns[3][search][regex]': 'false',
+    'columns[4][data]': '4',
+    'columns[4][name]': 'speed',
+    'columns[4][searchable]': 'true',
+    'columns[4][orderable]': 'true',
+    'columns[4][search][value]': '1|15',
+    'columns[4][search][regex]': 'false',
+    'columns[5][data]': '5',
+    'columns[5][name]': 'glide',
+    'columns[5][searchable]': 'true',
+    'columns[5][orderable]': 'true',
+    'columns[5][search][value]': '1|7',
+    'columns[5][search][regex]': 'false',
+    'columns[6][data]': '6',
+    'columns[6][name]': 'turn',
+    'columns[6][searchable]': 'true',
+    'columns[6][orderable]': 'true',
+    'columns[6][search][value]': '-5|1',
+    'columns[6][search][regex]': 'false',
+    'columns[7][data]': '7',
+    'columns[7][name]': 'fade',
+    'columns[7][searchable]': 'true',
+    'columns[7][orderable]': 'true',
+    'columns[7][search][value]': '0|5',
+    'columns[7][search][regex]': 'false',
+    'columns[8][data]': '8',
+    'columns[8][name]': 'inproduction',
+    'columns[8][searchable]': 'true',
+    'columns[8][orderable]': 'true',
+    'columns[8][search][value]': 'Coming Soon|Yes',
+    'columns[8][search][regex]': 'false',
+    'columns[9][data]': '9',
+    'columns[9][name]': 'dateapproved',
+    'columns[9][searchable]': 'true',
+    'columns[9][orderable]': 'true',
+    'columns[9][search][value]': '|',
+    'columns[9][search][regex]': 'false',
+    'columns[10][data]': '10',
+    'columns[10][name]': 'link',
+    'columns[10][searchable]': 'true',
+    'columns[10][orderable]': 'true',
+    'columns[10][search][value]': '',
+    'columns[10][search][regex]': 'false',
+    'order[0][column]': '0',
+    'order[0][dir]': 'asc',
+    'start': '0',
+    'length': '10',
+    'search[value]': 'wraith',
+    'search[regex]': 'false',
+    'wdtNonce': '511bd3400c',
+    'sRangeSeparator': '|',
+}
+
+response = requests.post('https://alldiscs.com/wp-admin/admin-ajax.php', params=params, headers=headers, data=data)

app/visualize.py

@@ -0,0 +1,272 @@
+import math
+
+import numpy as np
+import plotly.graph_objects as go
+from plotly.colors import sequential
+from stl.mesh import Mesh
+
+from extrema import find_extrema
+
+
+def get_stl(stl_file):
+    """
+    Taken from https://community.plotly.com/t/view-3d-cad-data/16920/9
+    """
+    stl_mesh = Mesh.from_file(stl_file)
+    return stl_mesh
+
+
+def visualize_disc(stl_mesh, nose, roll):
+    """
+    Taken from https://community.plotly.com/t/view-3d-cad-data/16920/9
+    """
+    stl_mesh.rotate([1, 0, 0], math.radians(-1*nose))
+    stl_mesh.rotate([0, 1, 0], math.radians(roll))
+    # stl_mesh.rotate([0, 0, 1], math.radians(z_angle))
+
+    p, q, r = stl_mesh.vectors.shape  # (p, 3, 3)
+    # the array stl_mesh.vectors.reshape(p*q, r) can contain multiple copies of the same vertex;
+    # extract unique vertices from all mesh triangles
+    vertices, ixr = np.unique(stl_mesh.vectors.reshape(p * q, r), return_inverse=True, axis=0)
+    I = np.take(ixr, [3 * k for k in range(p)])
+    J = np.take(ixr, [3 * k + 1 for k in range(p)])
+    K = np.take(ixr, [3 * k + 2 for k in range(p)])
+
+    x, y, z = vertices.T
+    trace = go.Mesh3d(x=x, y=y, z=z, i=I, j=J, k=K)
+    # optional parameters to make it look nicer
+    trace.update(flatshading=True, lighting_facenormalsepsilon=0, lighting_ambient=0.7)
+
+    fig = go.Figure(trace)
+
+    # Add camera controls to the plot
+    camera = dict(
+            eye=dict(x=0, y=-3, z=0)
+            )
+    fig.update_layout(scene_camera=camera)
+
+    fig.update_layout(
+            scene=dict(
+                    xaxis=dict(nticks=4, range=[-0.11, 0.11], ),
+                    yaxis=dict(nticks=4, range=[-0.11, 0.11], ),
+                    zaxis=dict(nticks=4, range=[-0.11, 0.11], ), )
+            )
+    return fig
+
+
+import plotly.subplots as sp
+
+
+def get_plot(x, y, z):
+    xm = np.min(x) - 1.5
+    xM = np.max(x) + 1.5
+    ym = -5
+    yM = np.max(y) + 1.5
+    zm = np.min(z)
+    zM = np.max(z)
+    N = len(x)
+    category = 'Height'
+    x_extrema, y_extrema, extrema_type = find_extrema(x, y)
+
+    xM_abs = max(abs(xm), abs(xM))
+    xm_abs = -1 * xM_abs
+
+    carats_v = go.Scatter(
+            x=[category, category],
+            y=[min(z), max(z)],
+            mode='markers',
+            showlegend=False,
+            marker=dict(symbol=['arrow-up', 'arrow-down'], size=20, color=['red', 'red']),
+            name='Carets',
+            )
+    carats_h = go.Scatter(
+            x=[min(x), max(x)],
+            y=['', ''],
+            mode='markers',
+            showlegend=False,
+            marker=dict(symbol=['arrow-right', 'arrow-left'], size=20, color=['red', 'red']),
+            name='Carets',
+            )
+
+    extrema = go.Scatter(
+            x=x_extrema,
+            y=y_extrema,
+            mode='markers',
+            showlegend=False,
+            marker=dict(symbol=extrema_type, size=20, color='blue'),
+            name='Extrema',
+            )
+
+    # Create figure with subplots
+    fig = sp.make_subplots(rows=2, cols=2, subplot_titles=("Flight Path", "Height", "Lateral Deviance"),
+                           specs=[[{"rowspan": 2}, {}], [None, {}]], row_heights=[0.5, 0.5])
+
+    # Add traces to the main plot
+    fig.add_trace(
+            go.Scatter(x=x, y=y,
+                       mode="lines",
+                       showlegend=False,
+                       line=dict(width=1, color='black')),
+            row=1, col=1
+            )
+
+    fig.add_trace(
+            go.Scatter(x=x, y=y,
+                       showlegend=False,
+                       mode="markers", marker_colorscale=sequential.Peach,
+                       marker=dict(color=z, size=3, showscale=True)),
+            row=1, col=1
+            )
+
+    fig.add_trace(
+            extrema,
+            row=1, col=1
+            )
+
+    # Add trace for the subplot
+    fig.add_trace(carats_v,
+                  row=1, col=2
+                  )
+    fig.add_trace(
+            go.Bar(
+                    x=[],
+                    y=[],
+                    showlegend=False,
+                    ),
+            row=1, col=2
+            )
+    # Add trace for the subplot
+    fig.add_trace(carats_h,
+                  row=2, col=2
+                  )
+    fig.add_trace(
+            go.Bar(
+                    x=[],
+                    y=[],
+                    showlegend=False,
+                    orientation='h',
+                    ),
+            row=2, col=2
+            )
+
+    # Update layout
+    fig.update_layout(
+            xaxis=dict(range=[xm, xM], autorange=False, zeroline=False),
+            yaxis=dict(range=[ym, yM], autorange=False, zeroline=False),
+            title_text="Flight Path",
+            hovermode="closest",
+            updatemenus=[
+                dict(
+                        type="buttons",
+                        buttons=[
+                            dict(
+                                    label="Play",
+                                    method="animate",
+                                    args=[None, {"frame": {"duration": 30, "redraw": True}, "fromcurrent": True}]
+                                    ),
+                            dict(
+                                    label="Pause",
+                                    method="animate",
+                                    args=[[None], {"frame": {"duration": 0, "redraw": False}, "mode": "immediate",
+                                                   "transition": {"duration": 0}}]
+                                    )
+                            ],
+                        showactive=False,
+                        x=0.05,
+                        y=0.05
+                        )
+                ],
+            )
+
+    # Create frames for the main plot and subplot
+    all_frames = [
+        go.Frame(data=[
+            go.Scatter(
+                    x=[x[k]],
+                    y=[y[k]],
+                    mode="markers",
+                    showlegend=False,
+                    marker=dict(color="red", size=10)),
+            go.Scatter(x=x, y=y,
+                       mode="markers", marker_colorscale=sequential.Peach,
+                       marker=dict(color=z, size=3, showscale=True)),
+            extrema,
+            carats_v,
+            go.Bar(
+                    x=['Height'],
+                    y=[z[k]],
+                    showlegend=False,
+                    name='Value',
+                    marker=dict(color='orange', line=dict(width=1))
+                    # Set color of value bar to orange and line width to 1
+                    ),
+            carats_h,
+            go.Bar(
+                    x=[x[k]],
+                    y=[''],
+                    showlegend=False,
+                    name='Value',
+                    orientation='h',
+                    marker=dict(color='orange', line=dict(width=1))
+                    # Set color of value bar to orange and line width to 1
+                    )
+            ])
+        for k in range(N)
+        ]
+
+    # Combine frames for the main plot and subplot
+    fig.frames = all_frames
+    fig.update_yaxes(scaleanchor="x", scaleratio=1, row=1, col=1)
+    fig.update_yaxes(range=[zm - 2, zM + 2], fixedrange=True, row=1, col=2)
+    fig.update_xaxes(range=[xm_abs, xM_abs], fixedrange=True, row=2, col=2)
+
+    # # Add green rectangle at the bottom of the plot
+    fig.update_layout(
+            shapes=[dict(type="rect", xref="x", yref="y",
+                         x0=-1, y0=0, x1=1, y1=-4, fillcolor="gray",
+                         opacity=1, layer="below")],
+            plot_bgcolor="green",
+            )
+
+    return fig
+
+
+import plotly.graph_objs as go
+from plotly.subplots import make_subplots
+
+
+def get_subplots(arc, alphas, lifts, drags, moms, rolls, velocity):
+    fig = make_subplots(rows=2, cols=3, specs=[[{}, {}, {}], [{}, {}, {}]],
+                        subplot_titles=("Lift force (N)", "Drag force (N)", "Moment (Nm)",
+                                        "Angle of attack (deg)", "Velocities (m/s)", "Roll rate (rad/s)"),
+                        shared_xaxes=True)
+
+    fig.add_trace(go.Scatter(x=arc, y=lifts, name="Lift force (N)"), row=1, col=1)
+    fig.update_xaxes(title_text="Distance (m)", row=1, col=1)
+    fig.update_yaxes(title_text="Lift force (N)", row=1, col=1)
+
+    fig.add_trace(go.Scatter(x=arc, y=drags, name="Drag force (N)"), row=1, col=2)
+    fig.update_xaxes(title_text="Distance (m)", row=1, col=2)
+    fig.update_yaxes(title_text="Drag force (N)", row=1, col=2)
+
+    fig.add_trace(go.Scatter(x=arc, y=moms, name="Moment (Nm)"), row=1, col=3)
+    fig.update_xaxes(title_text="Distance (m)", row=1, col=3)
+    fig.update_yaxes(title_text="Moment (Nm)", row=1, col=3)
+
+    fig.add_trace(go.Scatter(x=arc, y=alphas, name="Angle of attack (deg)"), row=2, col=1)
+    fig.update_xaxes(title_text="Distance (m)", row=2, col=1)
+    fig.update_yaxes(title_text="Angle of attack (deg)", row=2, col=1)
+
+    fig.add_trace(go.Scatter(x=arc, y=velocity[0, :], name="u"), row=2, col=2)
+    fig.add_trace(go.Scatter(x=arc, y=velocity[1, :], name="v"), row=2, col=2)
+    fig.add_trace(go.Scatter(x=arc, y=velocity[2, :], name="w"), row=2, col=2)
+    fig.update_xaxes(title_text="Distance (m)", row=2, col=2)
+    fig.update_yaxes(title_text="Velocities (m/s)", row=2, col=2)
+    fig.update_traces(mode='lines', row=2, col=2)
+
+    fig.add_trace(go.Scatter(x=arc, y=rolls, name="Roll rate (rad/s)"), row=2, col=3)
+    fig.update_xaxes(title_text="Distance (m)", row=2, col=3)
+    fig.update_yaxes(title_text="Roll rate (rad/s)", row=2, col=3)
+
+    fig.update_layout(height=600, width=1000, title_text="Plotly Subplots", hovermode='x')
+    return fig

bootstrap.py

@@ -0,0 +1,12 @@
+from pathlib import Path
+import streamlit.web.bootstrap
+from streamlit import config as _config
+
+proj_dir = Path(__file__).parent
+filename = proj_dir / "app" / "app.py"
+
+_config.set_option("server.headless", True)
+args = []
+
+# streamlit.cli.main_run(filename, args)
+streamlit.web.bootstrap.run(str(filename), "", args, "")

examples/disc_experimental_trajectory_comparison.py

@@ -0,0 +1,106 @@
+# -*- coding: utf-8 -*-
+
+import sys
+from shotshaper.projectile import DiscGolfDisc
+import matplotlib.pyplot as pl
+import numpy as np
+import shotshaper.environment as env
+from shotshaper.transforms import T_12
+from random import uniform
+
+
+def rotz(pos,betad):
+    beta = np.radians(betad)
+    TZ =  np.array([[np.cos(beta), -np.sin(beta), 0],
+                    [np.sin(beta),  np.cos(beta), 0],
+                    [0,             0,            1]])
+
+    return np.matmul(TZ,pos)
+
+throws=[1,6,15]
+
+nthrow = len(throws)
+        # pitch, roll, nose,speed,spin,  yaw,  wind, length
+params = [[15.5, 21.8, 0.0, 24.7, 138, -31.6, 4.8, 89.7],
+          [12.3, 14.7, 0.8, 24.2, 128.5, -9.60, 4.8, 87.0],
+          [5.20,-0.70, 0.0, 24.5, 147.7,  8.00, 4.8, 106.6]]
+         
+
+d = DiscGolfDisc('dd2')
+fig1, ax1 = pl.subplots( )
+    
+fig1.set_figheight(4)
+fig1.set_figwidth(6)
+
+for i in range(nthrow):
+    p = params[i]
+    t = throws[i]
+    U = p[3]
+    omega = p[4]
+    z0 = 1.5
+    pos = np.array((0,0,z0))
+    pitch = p[0]
+    yaw = p[5]
+    nose = p[2]
+    roll = p[1]
+    
+    env.Uref = p[6]
+    env.winddir = np.array((1,0,0))
+
+    # Currently handle yaw by rotating the position after the throw, hence
+    # also need to rotate the wind vector accordingly
+    env.winddir = rotz(env.winddir, -yaw)
+
+    s = d.shoot(speed=U, omega=omega, pitch=pitch, position=pos, nose_angle=nose, roll_angle=roll)
+
+    pos = s.position
+    for j in range(len(pos[0,:])):
+        pos[:,j] = rotz(pos[:,j], yaw)
+    x,y,z = pos
+    arc,alphas,betas,lifts,drags,moms,rolls = d.post_process(s, omega)
+    
+    # Plot trajectory
+    ax1.plot(x,y,f'C{i}-')
+    
+    # Experiment
+    te,xe,ye,ve = np.loadtxt(f'data/throw{t}',skiprows=2,unpack=True)
+    ax1.plot(xe,ye,f'C{i}--')
+
+    ax1.set_xlabel('Distance (m)')
+    ax1.set_ylabel('Drift (m)')
+    ax1.axis('equal')
+    
+    # Plot other parameters
+    
+    # axes[0,0].plot(arc, lifts)
+    # axes[0,0].set_xlabel('Distance (m)')
+    # axes[0,0].set_ylabel('Lift force (N)')
+    
+    # axes[0,1].plot(arc, drags)
+    # axes[0,1].set_xlabel('Distance (m)')
+    # axes[0,1].set_ylabel('Drag force (N)')
+    
+    # axes[0,2].plot(arc, moms)
+    # axes[0,2].set_xlabel('Distance (m)')
+    # axes[0,2].set_ylabel('Moment (Nm)')
+    
+    # axes[1,0].plot(arc, alphas)
+    # axes[1,0].set_xlabel('Distance (m)')
+    # axes[1,0].set_ylabel('Angle of attack (deg)')
+    
+    # axes[1,1].plot(arc, s.velocity[0,:])
+    # axes[1,1].plot(arc, s.velocity[1,:])
+    # axes[1,1].plot(arc, s.velocity[2,:])
+    # axes[1,1].set_xlabel('Distance (m)')
+    # axes[1,1].set_ylabel('Velocities (m/s)')
+    
+    # axes[1,2].plot(arc, rolls)
+    # axes[1,2].set_xlabel('Distance (m)')
+    # axes[1,2].set_ylabel('Roll rate (rad/s)')
+
+pl.tight_layout()
+pl.show()
+
+
+
+

examples/disc_golf_throw.py

@@ -0,0 +1,66 @@
+# -*- coding: utf-8 -*-
+"""
+Example showing a single disc throw.
+"""
+
+from shotshaper.projectile import DiscGolfDisc
+import matplotlib.pyplot as pl
+import numpy as np
+import plotly.express as px
+
+d = DiscGolfDisc('dd2')
+U = 24.2
+omega = 116.8
+z0 = 1.3
+pos = np.array((0,0,z0))
+pitch = 15.5
+nose = 0.0
+roll = 14.7
+
+shot = d.shoot(speed=U, omega=omega, pitch=pitch, 
+               position=pos, nose_angle=nose, roll_angle=roll)
+
+# Plot trajectory
+pl.figure(1)
+x,y,z = shot.position
+pl.plot(x,y)
+fig = px.scatter(x,y)
+fig.show()
+
+pl.xlabel('Distance (m)')
+pl.ylabel('Drift (m)')
+pl.axis('equal')
+
+# Plot other parameters
+arc,alphas,betas,lifts,drags,moms,rolls = d.post_process(shot, omega)
+fig, axes = pl.subplots(nrows=2, ncols=3, dpi=80,figsize=(13,5))
+
+axes[0,0].plot(arc, lifts)
+axes[0,0].set_xlabel('Distance (m)')
+axes[0,0].set_ylabel('Lift force (N)')
+
+axes[0,1].plot(arc, drags)
+axes[0,1].set_xlabel('Distance (m)')
+axes[0,1].set_ylabel('Drag force (N)')
+
+axes[0,2].plot(arc, moms)
+axes[0,2].set_xlabel('Distance (m)')
+axes[0,2].set_ylabel('Moment (Nm)')
+
+axes[1,0].plot(arc, alphas)
+axes[1,0].set_xlabel('Distance (m)')
+axes[1,0].set_ylabel('Angle of attack (deg)')
+
+axes[1,1].plot(arc, shot.velocity[0,:])
+axes[1,1].plot(arc, shot.velocity[1,:])
+axes[1,1].plot(arc, shot.velocity[2,:])
+axes[1,1].set_xlabel('Distance (m)')
+axes[1,1].set_ylabel('Velocities (m/s)')
+axes[1,1].legend(('u','v','w'))
+
+axes[1,2].plot(arc, rolls)
+axes[1,2].set_xlabel('Distance (m)')
+axes[1,2].set_ylabel('Roll rate (rad/s)')
+pl.tight_layout()
+
+pl.show()

examples/disc_gui2d.py

@@ -0,0 +1,105 @@
+# -*- coding: utf-8 -*-
+"""
+Graphical user interface to explore the influence
+of disc throw parameters on the trajectory
+"""
+
+import numpy as np
+import matplotlib.pyplot as pl
+from matplotlib.widgets import Slider, TextBox
+from shotshaper.projectile import DiscGolfDisc
+
+name = 'dd2'
+mass = 0.175
+
+d = DiscGolfDisc(name, mass=mass)
+speed = 24
+omega = d.empirical_spin(speed)
+z0 = 1.3
+pos = np.array((0,0,z0))
+pitch = 10
+nose = 0.0
+roll = 15.0  
+yaw = 0
+adjust_axes = False
+     
+s = d.shoot(speed=speed, omega=omega, pitch=pitch, position=pos, nose_angle=nose, roll_angle=roll,yaw=yaw)
+
+x,y,z = s.position
+
+# Creating figure
+fig = pl.figure(1,figsize=(13, 6), dpi=80)
+ax1 = pl.subplot(2,3,2)
+ax2 = pl.subplot(2,3,5)
+ax3 = pl.subplot(1,3,3)
+fac = 1.6
+ylim = fac*max(abs(min(y)),abs(max(y)))
+ax1.axis((min(x),fac*max(x),-ylim,ylim))
+ax2.axis((min(x),fac*max(x),min(z),fac*max(z)))
+ax3.axis((-ylim,ylim,min(z),fac*max(z)))
+        
+ax3.invert_xaxis()
+
+l1, = ax1.plot(x,y,lw=2)
+ax1.set_xlabel('Distance (m)')
+ax1.set_ylabel('Drift (m)')
+l2, = ax2.plot(x,z,lw=2)
+ax2.set_xlabel('Distance (m)')
+ax2.set_ylabel('Height (m)')
+
+l3, = ax3.plot(y,z,lw=2)
+ax3.set_xlabel('Drift (m)')
+ax3.set_ylabel('Height (m)')
+
+xax = 0.07
+ax4  = pl.axes([xax, 0.80, 0.25, 0.03], facecolor='lightgrey')
+ax5  = pl.axes([xax, 0.75, 0.25, 0.03], facecolor='lightgrey')
+ax6  = pl.axes([xax, 0.70, 0.25, 0.03], facecolor='lightgrey')
+#ax7  = pl.axes([xax, 0.65, 0.25, 0.03], facecolor='lightgrey')
+ax8  = pl.axes([xax, 0.65, 0.25, 0.03], facecolor='lightgrey')
+ax9  = pl.axes([xax, 0.60, 0.25, 0.03], facecolor='lightgrey')
+ax11 = pl.axes([xax, 0.55, 0.25, 0.03], facecolor='lightgrey')
+
+s1 = Slider(ax=ax4, label='Speed (m/s)', valmin=15,  valmax=35, valinit=speed)
+s2 = Slider(ax=ax5, label='Roll (deg)',   valmin=-110, valmax=110, valinit=roll)
+s3 = Slider(ax=ax6, label='Pitch (deg)', valmin=-10,   valmax=50, valinit=pitch)
+s5 = Slider(ax=ax8, label='Nose (deg)',   valmin=-5, valmax=5, valinit=nose)
+s7 = Slider(ax=ax9, label='Mass (kg)',   valmin=0.140, valmax=0.200, valinit=mass)
+s6 = Slider(ax=ax11, label='Spin (-)',   valmin=0, valmax=2, valinit=1.0)
+
+def update(x):
+    speed = s1.val
+    roll = s2.val
+    pitch = s3.val
+    nose = s5.val
+    spin = s6.val
+    mass = s7.val 
+    d = DiscGolfDisc(name,mass=mass)
+    
+    omega = spin*d.empirical_spin(speed)
+    s = d.shoot(speed=speed, omega=omega, pitch=pitch, position=pos, nose_angle=nose, roll_angle=roll)
+    x,y,z = s.position
+    
+    l1.set_xdata(x)
+    l1.set_ydata(y)
+    l2.set_xdata(x)
+    l2.set_ydata(z)
+    l3.set_xdata(y)
+    l3.set_ydata(z)
+    
+    if adjust_axes:
+        ax1.axis((min(x),max(x),min(y),max(y)))
+        ax2.axis((min(x),max(x),min(z),max(z)))
+        ax3.axis((min(y),max(y),min(z),max(z)))
+        
+    fig.canvas.draw_idle()
+
+s1.on_changed(update)
+s2.on_changed(update)
+s3.on_changed(update)
+s5.on_changed(update)
+s6.on_changed(update)
+s7.on_changed(update)
+
+pl.show()
+    

examples/driver.py

@@ -0,0 +1,35 @@
+# -*- coding: utf-8 -*-
+"""
+Simple example comparing different projectiles.
+"""
+
+from shotshaper.projectile import _Particle, ShotPutBall, SoccerBall
+import matplotlib.pyplot as pl
+import numpy as np
+
+U = 10.0
+angle = 20.0
+
+p = _Particle()
+shot = p.shoot(speed=U, pitch=angle)
+pl.plot(shot.position[0,:],shot.position[2,:])
+
+p = ShotPutBall('M')
+shot = p.shoot(speed=U, pitch=angle)
+pl.plot(shot.position[0,:],shot.position[2,:])
+
+p = SoccerBall()
+
+spin = np.array((0,0,0))
+shot = p.shoot(speed=U, pitch=angle, spin=spin)
+pl.plot(shot.position[0,:],shot.position[2,:])
+
+spin = np.array((0,-10,0))
+shot = p.shoot(speed=U, pitch=angle, spin=spin)
+pl.plot(shot.position[0,:],shot.position[2,:])
+
+pl.legend(('vacuum','air', 'no spin','spin'))
+
+pl.show()
+
+

examples/validation_balls.py

@@ -0,0 +1,50 @@
+# -*- coding: utf-8 -*-
+"""
+This program calculates trajectories that were experimentally conducted in
+
+Mencke, J. E., Salewski, M., Trinhammer, O. L., & Adler, A. T. (2020). 
+Flight and bounce of spinning sports balls. American Journal of Physics, 
+88(11), 934-947.
+
+Specifically, this includes data for:
+    - a shot put throw, representing a heavy projectile where the 
+      gravitational forces dominate over the aerodynamic forces
+    - a soccer ball, with slightly more influence of air drag
+    - a table tennis ball, significantly impacted by air drag due
+      to low weight
+"""
+
+from shotshaper.projectile import SoccerBall, ShotPutBall, TableTennisBall
+import matplotlib.pyplot as pl
+import numpy as np
+
+cases = ('Shot','Soccer ball','Table tennis ball')
+projectiles = (ShotPutBall('M'), SoccerBall(), TableTennisBall())
+ux = (9.0, 11.0, 11.8)
+uy = (7.1, 10.2, 12.0)
+z0 = (2.4, 0.0, 1.0)
+spin = (0,0,0)
+
+f, ax = pl.subplots(1, 1, figsize=(6,4))
+for i,c in enumerate(cases):
+    speed = np.sqrt(ux[i]**2 + uy[i]**2)
+    pitch = np.degrees(np.arctan2(uy[i],ux[i]))
+    position = (0, 0, z0[i])
+    
+    s = projectiles[i].shoot(speed=speed,pitch=pitch,position=position,spin=spin)
+
+    x,y,z = s.position
+
+    ax.plot(x,z,'C0-')
+    ax.text(x[0]-0.5, z[0], c, fontsize=9, ha='right')
+    
+    x,z = np.loadtxt(f'data/{c}_trajectory.dat', unpack=True, delimiter=';')
+    ax.plot(x,z,'C1--')
+    
+ax.legend(('Simulation','Experiment'))
+ax.axis((-10,25,-1,6))
+    
+pl.xlabel('Length (m)')
+pl.ylabel('Height (m)')
+
+pl.show()

examples/validation_spinning_ball.py

@@ -0,0 +1,48 @@
+# -*- coding: utf-8 -*-
+"""
+This program calculates the trajectory of a spinning soccer ball
+that was experimentally investigated in
+
+Mencke, J. E., Salewski, M., Trinhammer, O. L., & Adler, A. T. (2020). 
+Flight and bounce of spinning sports balls. American Journal of Physics, 
+88(11), 934-947.
+
+"""
+
+from shotshaper.projectile import SoccerBall
+import matplotlib.pyplot as pl
+import numpy as np
+from scipy.linalg import norm
+
+ball = SoccerBall()
+u = np.array([20.0, 2.5, 4.9])
+# Note that spin is here negative along the z-axis, as opposed
+# to positive around the y-axis, since we define z as upwards,
+# while Mencke et al. define y as upwards 
+spin = (0,0,-46) 
+speed = norm(u)
+pitch = np.degrees(np.arctan2(u[2],u[0]))
+yaw = -np.degrees(np.arctan2(u[1],u[0]))
+    
+s = ball.shoot(speed=speed,pitch=pitch,yaw=yaw,spin=spin)
+
+x,y,z = s.position
+
+f, (ax1, ax2) = pl.subplots(2, 1, figsize=(5,7))
+
+ax1.plot(x,z,'C0-')
+ax2.plot(x,y,'C0-')    
+x,z = np.loadtxt('data/Spin_z_trajectory.dat', unpack=True, delimiter=';')
+ax1.plot(x,z,'C1--')
+x,y = np.loadtxt('data/Spin_y_trajectory.dat', unpack=True, delimiter=';')
+ax2.plot(x,y,'C1--')    
+
+ax1.legend(('Simulation','Experiment'))
+ax2.legend(('Simulation','Experiment'))
+    
+ax1.set_xlabel('Length (m)')
+ax1.set_ylabel('Height (m)')
+ax2.set_xlabel('Length (m)')
+ax2.set_ylabel('Drift (m)')
+
+pl.show()

notebooks/disc_plot.ipynb

@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:1006c45576f250234568583529a2f594a1d6fd5400aeb877b2dbee7dd770cde9
+size 11383386

requirements.txt

@@ -0,0 +1,5 @@
+numpy>=1.18.1
+scipy>=1.4.1
+PyYAML==6.0
+numpy_stl==3.0.1
+plotly==5.13.1

setup.py

@@ -0,0 +1,7 @@
+from setuptools import setup, find_packages
+setup(
+    name="shotshaper",
+    version="0.1.0",
+    packages=find_packages(),
+    include_package_data=True
+)

shotshaper/__init__.py

@@ -0,0 +1,3 @@
+# -*- coding: utf-8 -*-
+
+

shotshaper/discs/.gitattributes

@@ -0,0 +1,4 @@
+cd1.stl filter=lfs diff=lfs merge=lfs -text
+cd5.stl filter=lfs diff=lfs merge=lfs -text
+dd2.stl filter=lfs diff=lfs merge=lfs -text
+fd2.stl filter=lfs diff=lfs merge=lfs -text

shotshaper/discs/cd1.stl

@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:750a9ec26e520f1d662453952e44a20c2b02eeb4d7f52a1d702c7a5dff335ad0
+size 3869884

shotshaper/discs/cd1.yaml

@@ -0,0 +1,24 @@
+# Overstable control driver, modelled after an Innova Firebird
+# Data taken from CFD simulations
+
+diameter: 0.211
+J_xy: 3.759e-03
+J_z: 7.485e-03
+alpha: [-90,-80,-70,-60,-50,-40,-30,-20,-15,
+        -10,-8,-6,-4,-2,0,2,4,6,8,10,12,14,16,18,20,
+        25,30,40,50,60,70,80,90]
+Cd: [ 0.9641927, 1.040903, 1.014881, 0.8892227, 0.79079, 0.674083, 0.5343097, 0.272792, 0.1632757, 
+
+ 0.0813, 0.0663, 0.0554, 0.0474, 0.0474, 0.0458, 0.057, 0.0668, 0.0846, 0.109, 0.137, 0.1729, 0.2182, 0.2672, 0.3268, 0.3924,
+
+.5411498, 0.7411681, 1.111394, 1.115411, 1.080624, 1.186235, 1.023045, 1.029478 ]
+
+Cl: [ 0.0001273178, -0.1869159, -0.3665411, -0.5029187, -0.6435312, -0.7736677, -0.8636631, -0.6659618, -0.4805856,
+
+-0.29, -0.219, -0.1522, -0.0802, -0.0026, 0.0674, 0.1646, 0.2486, 0.3467, 0.4523, 0.5633, 0.6728, 0.7904, 0.8948, 1.0084, 1.0936,
+1.199855, 1.309407, 1.322162, 0.9143433, 0.6057466, 0.4181427, 0.1679042, -0.001581172 ]
+
+Cm: [ 0.0003381886, -0.03136248, -0.04697123, -0.06378013, -0.07869955, -0.08915614, -0.100739, -0.07532892, -0.08049738, 
+ -0.075, -0.0617, -0.0484, -0.0346, -0.0238, -0.0099, 0.002, 0.0112, 0.0211, 0.0323, 0.0461, 0.0601, 0.0749, 0.0929, 0.1087, 0.1159,
+
+0.1514894, 0.1528616, 0.135503, 0.09038871, 0.06759896, 0.07489785, 0.04788332, -0.002841088 ]

shotshaper/discs/cd5.stl

@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:ffdbc36b511ef0a98ab99ab247c32353f0e15475d2ca2ef8fbf87f5d6e3c188d
+size 6317884

shotshaper/discs/cd5.yaml

@@ -0,0 +1,23 @@
+# Understable control driver, modelled after an Innova Roadrunner
+# Data taken from CFD simulations
+
+diameter: 0.211
+J_xy: 3.829e-03
+J_z: 7.616e-03
+alpha: [-90,-80,-70,-60,-50,-40,-30,-20,-15,
+        -10,-8,-6,-4,-2,0,2,4,6,8,10,12,14,16,18,20,
+        25,30,40,50,60,70,80,90]
+Cd: [ 0.9179461, 0.9709634, 0.9750627, 0.8941224, 0.7740944, 0.611287, 0.4877146, 0.2566395, 0.1423688, 
+
+ 0.0773, 0.0606, 0.0514, 0.0468, 0.0439, 0.0528, 0.0631, 0.0763, 0.0968, 0.1171, 0.1417, 0.1776, 0.2136, 0.2587, 0.2986, 0.3612,
+0.5332736, 0.7454067, 1.100109, 1.231819, 1.080511, 1.06314, 1.025655, 1.000069 ]
+
+Cl: [ -0.003535795, -0.1650353, -0.3380664, -0.4879364, -0.6047266, -0.6677451, -0.7633718, -0.6095486, -0.4310734, 
+                                             #
+ -0.2642, -0.1891, -0.1239, -0.0469, 0.0208, 0.1081, 0.1917, 0.2859, 0.4015, 0.4846, 0.5882, 0.6957, 0.7988, 0.9093, 1.0012, 1.1181,
+1.246332, 1.367988, 1.333375, 1.038364, 0.6176756, 0.4386948, 0.1819935, 0.0001378738 ]
+
+Cm: [ -0.004521051, -0.03329007, -0.05428122, -0.072994, -0.09090575, -0.09778813, -0.1002716, -0.08264331, -0.09134467, 
+ -0.0783, -0.0634, -0.0539, -0.0436, -0.0304, -0.0229, -0.012, -0.0013, 0.005, 0.019, 0.0302, 0.042, 0.0559, 0.0692, 0.0901, 0.1032,
+
+0.1538899, 0.1695394, 0.1369221, 0.1079746, 0.06659547, 0.07765108, 0.04283614, -2.493789e-05 ]

shotshaper/discs/dd2.stl

@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:0d9226b43dba147376cb204532fa9547c6729d61d086a51803a5769e56016d67
+size 3673484

shotshaper/discs/dd2.yaml

@@ -0,0 +1,21 @@
+diameter: 0.211
+J_xy: 3.755e-03
+J_z: 7.465e-03
+
+alpha: [-90,-80,-70,-60,-50,-40,-30,-20,-15,-10,-8,-6,-4,-2,0,2,4,6,8,10,12,14,
+         16,18,20,25,30,40,50,60,70,80,90]
+         
+Cd: [1.0156, 1.0413, 0.995, 0.8904, 0.7415, 0.6219, 0.4648, 0.2499, 0.1556, 
+     0.0769, 0.0617, 0.0509, 0.0432, 0.0429, 0.0417, 0.0578, 0.0664, 0.0825, 
+     0.1058, 0.1281, 0.1624, 0.2015, 0.2462, 0.2962, 0.3591, 0.561, 0.8141, 
+     1.2255, 1.2382, 1.0793, 1.1797, 1.199, 1.1867]
+
+Cl: [0.0001, -0.2141, -0.4086, -0.5726, -0.6814, -0.7987, -0.8534, -0.694, 
+    -0.4917, -0.2137, -0.1607, -0.0952, -0.0401, 0.0347, 0.1091, 0.2013, 
+     0.2761, 0.3666, 0.477, 0.5604, 0.6684, 0.7773, 0.888, 1.0013, 1.1157, 
+     1.67, 1.8447, 1.8161, 1.2563, 0.7495, 0.5236, 0.2646, 0.0033]
+
+Cm: [-0.0007, -0.0326, -0.0544, -0.0712, -0.085, -0.0985, -0.0953, -0.0728, 
+     -0.0742, -0.0966, -0.076, -0.0621, -0.0455, -0.0342, -0.0175, -0.0122, 
+     0.0013, 0.0128, 0.0217, 0.0364, 0.0486, 0.0618, 0.0758, 0.0925, 0.1081, 
+     0.1568, 0.1712, 0.1419, 0.0995, 0.0731, 0.0692, 0.0421, 0.0016]

shotshaper/discs/fd2.stl

@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:3f4a12f3fc29c8ddf1661c92ddd750e49ce2fb5f056fe3508b919fb39f9172cc
+size 4492684

shotshaper/discs/fd2.yaml

@@ -0,0 +1,27 @@
+# Stable fairway driver, modelled after an Innova Teebird
+# Data taken from CFD simulations (-90 to -15, and 25 to 90 copied from cd1)
+
+diameter: 0.211
+J_xy: 3.809e-03
+J_z: 7.577e-03
+alpha: [-90,-80,-70,-60,-50,-40,-30,-20,-15,
+        -10,-8,-6,-4,-2,0,2,4,6,8,10,12,14,16,18,20,
+        25,30,40,50,60,70,80,90]
+Cd: [ 0.9641927, 1.040903, 1.014881, 0.8892227, 0.79079, 0.674083, 0.5343097, 0.272792, 0.1632757, 
+
+#0.0913, 0.0644, 0.0543, 0.0502, 0.0458, 0.0436, 0.0594, 0.0671, 0.0875, 0.1111, 0.1415, 0.1749, 0.2152, 0.2613, 0.3132, 0.3829,
+  0.08, 0.0674, 0.0555, 0.0496, 0.046, 0.0454, 0.0488, 0.0701, 0.0872, 0.1096, 0.1336, 0.1688, 0.212, 0.2604, 0.3172, 0.387,
+.5411498, 0.7411681, 1.111394, 1.115411, 1.080624, 1.186235, 1.023045, 1.029478 ]
+
+Cl: [ 0.0001273178, -0.1869159, -0.3665411, -0.5029187, -0.6435312, -0.7736677, -0.8636631, -0.6659618, -0.4805856,
+
+#-0.2833, -0.1839, -0.1311, -0.0441, 0.0169, 0.0734, 0.1942, 0.2687, 0.3722, 0.478, 0.5899, 0.6983, 0.8132, 0.9264, 1.0354, 1.1698,
+-0.2624, -0.1941, -0.1464, -0.0673, -0.0005, 0.0672, 0.1408, 0.2622, 0.3575, 0.4587, 0.5593, 0.673, 0.7957, 0.9149, 1.0373, 1.148,
+
+1.199855, 1.309407, 1.322162, 0.9143433, 0.6057466, 0.4181427, 0.1679042, -0.001581172 ]
+
+Cm: [ 0.0003381886, -0.03136248, -0.04697123, -0.06378013, -0.07869955, -0.08915614, -0.100739, -0.07532892, -0.08049738, 
+
+# -0.0763, -0.0653, -0.0475, -0.0432, -0.027, -0.0128, -0.0112, 0.0016, 0.0095, 0.019, 0.0292, 0.042, 0.055, 0.0701, 0.0888, 0.1095,
+ -0.0694, -0.0567, -0.0429, -0.0361, -0.0233, -0.0112, 0.0006, 0.0026, 0.0125, 0.0224, 0.0344, 0.0452, 0.0566, 0.0703, 0.0886, 0.1093,
+0.1514894, 0.1528616, 0.135503, 0.09038871, 0.06759896, 0.07489785, 0.04788332, -0.002841088 ]

shotshaper/environment.py

@@ -0,0 +1,27 @@
+# -*- coding: utf-8 -*-
+"""
+"""
+import numpy as np
+
+g = -9.81
+# Air
+rho = 1.225
+mu = 1.81e-5
+
+winddir = np.array((1,0,0))
+z0 = 0.1
+Uref = 0.0
+zref = 1.5
+kappa = 0.41
+
+def wind_abl(z):
+    # For a constant wind:
+    # return Uref*winddir
+
+    if z < 0.0:
+        z = 0.0
+    
+    ustar = Uref*kappa/(np.log((zref+z0)/z0))
+    u = ustar/kappa*np.log((z+z0)/z0)
+    
+    return u*winddir

shotshaper/projectile.py

@@ -0,0 +1,607 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Tue Oct 19 18:29:10 2021
+
+@author: 2913452
+
+TODO:
+    - height of athlete, optimal angle shot put
+        - note that biomech can influence
+          how much force an athlete can emit for each angle
+    - 
+"""
+
+from abc import ABC, abstractmethod
+from scipy.integrate import solve_ivp
+from scipy.interpolate import interp1d
+from .transforms import T_12, T_23, T_34, T_14, T_41, T_31
+import matplotlib.pyplot as pl
+from numpy import exp,matmul,pi,sqrt,arctan2,radians,degrees,sin,cos,array,concatenate,linspace,zeros_like,cross,zeros,argmin
+from numpy.linalg import norm
+from . import environment
+import os
+import yaml
+
+T_END = 60
+N_STEP = 200
+
+def hit_ground(t, y, *args): 
+    return y[2]
+
+def stopped(t, y, *args):
+    U = norm(y[3:6])
+    return U - 1e-4
+
+class Shot:
+    def __init__(self,t,x,v,att=None):
+        self.time = t
+        self.position = x
+        self.velocity = v
+        if att is not None:
+            self.attitude = att
+        
+
+class _Projectile(ABC):
+    def __init__(self):
+        pass
+   
+    def _shoot(self, advance_function, y0, *args):
+        hit_ground.terminal = True
+        hit_ground.direction = -1
+        stopped.terminal = True
+        stopped.direction = -1
+        
+        sol = solve_ivp(advance_function,[0,T_END],y0,
+                        dense_output=True,args=args,
+                        method='RK45',
+                        events=(hit_ground,stopped))
+        
+        t = linspace(0,sol.t[-1],N_STEP)
+        
+        f = sol.sol(t)
+        pos = array([f[0],f[1],f[2]])
+        vel = array([f[3],f[4],f[5]])
+        
+        if len(f) <= 6:
+            shot = Shot(t, pos, vel)
+        else:    
+            att = array([f[6],f[7],f[8]])
+            shot = Shot(t, pos, vel, att)
+        
+        return shot
+         
+    @abstractmethod
+    def advance(self,t,vec,*args):
+        """
+        :param float T: Thrust
+        :param float Q: Torque
+        :param float P: Power
+        :return: Right hand side of kinematic equations for a projectile
+        :rtype: array
+        """
+        
+class _Particle(_Projectile):
+    def __init__(self):
+        super().__init__()
+        
+        self.g = environment.g
+        
+    def initialize_shot(self, **kwargs):
+        kwargs.setdefault('yaw', 0.0) 
+        
+        pitch = radians(kwargs["pitch"])
+        yaw = radians(kwargs["yaw"])
+        U = kwargs["speed"]
+        xy = cos(pitch)
+        u = U*xy*cos(yaw)
+        w = U*sin(pitch)
+        v = U*xy*sin(-yaw)
+        if "position" in kwargs:
+            x,y,z = kwargs["position"]
+        else:
+            x = 0.
+            y = 0.
+            z = 0.
+        
+        y0 = array((x,y,z,u,v,w))
+        return y0
+            
+    def shoot(self, **kwargs):
+
+        y0 = self.initialize_shot(**kwargs)
+        shot = self._shoot(self.advance, y0)
+        
+        return shot
+        
+    def gravity_force(self, x=None):
+        if x is None:
+            return array((0,0,environment.g))
+        else:
+            # Messy way to return g array...
+            f = zeros_like(x)
+            f[0,:] = 0
+            f[1,:] = 0
+            f[2,:] = environment.g
+            return f
+        
+    def advance(self, t, vec, *args):
+        # x, y, z, u, v, w = vec
+        x = vec[0:3]
+        u = vec[3:6]
+        
+        f = self.gravity_force()
+        
+        return concatenate((u,f))
+        
+    
+class _SphericalParticleAirResistance(_Particle):
+    def __init__(self, mass, diameter):
+        super().__init__()
+        
+        self.mass = mass
+        self.diameter = diameter
+        self.radius = 0.5*diameter
+        self.area = 0.25*pi*diameter**2
+        self.volume = 4./3.*pi*self.radius**3
+        
+    
+    def air_resistance_force(self, U, Cd):
+        
+        f = -0.5*environment.rho*self.area*Cd*norm(U)*U/self.mass
+        #f = -0.5*environment.rho*self.area*Cd*Umag*U/self.mass
+        
+        return f
+    
+    def advance(self, t, vec, *args):
+        x = vec[0:3]
+        u = vec[3:6]
+        
+        Cd = self.drag_coefficient(norm(u))
+        
+        f = self.air_resistance_force(u, Cd) \
+          + self.gravity_force()
+        
+        return concatenate((u,f))
+       
+        
+    def reynolds_number(self, velocity):
+        """
+        Reynolds number, non-dimensional number giving the 
+        ratio of inertial forces to viscous forces. Used
+        for calculating the drag coefficient.
+        
+        :param float velocity: Velocity seen by particle
+        :return: Reynolds number
+        :rtype: float
+        
+        """
+        return environment.rho*velocity*self.diameter/environment.mu
+    
+    def drag_coefficient(self, velocity):
+        """
+        Drag coefficient for sphere, empirical curve fit
+        taken from:
+        
+        F. A. Morrison, An Introduction to Fluid Mechanics, (Cambridge
+        University Press, New York, 2013). This correlation appears in
+        Figure 8.13 on page 625. 
+
+        The full formula is:
+
+        .. math::
+            F = \\frac{2}{\\pi}\\cos^{-1}e^{-f} \\\\
+            f = \\frac{B}{2}\\frac{R-r}{r\\sin\\phi}
+
+
+        :param float velocity: Velocity seen by particle
+        :return: Drag coefficient
+        :rtype: float
+        """
+    
+        Re = self.reynolds_number(velocity)
+        
+        if Re <= 0:
+            return 1e30
+        
+        tmp1 = Re/5.0
+        tmp2 = Re/2.63e5
+        tmp3 = Re/1e6
+        
+        Cd = 24.0/Re \
+           + 2.6*tmp1/(1 + tmp1**1.52) \
+           + 0.411*tmp2**-7.94/(1 + tmp2**-8) \
+           + 0.25*tmp3/(1 + tmp3) 
+           
+        return Cd
+    
+
+class _SphericalParticleAirResistanceSpin(_SphericalParticleAirResistance):
+    def __init__(self, mass, diameter):
+        super().__init__(mass, diameter)
+        
+        
+    def lift_coefficient(self, Umag, omega):
+        # TODO - complex dependency on Re. For now,
+        #        assume constant
+        return 0.9
+        
+    def shoot(self, **kwargs):
+        y0 = self.initialize_shot(**kwargs)
+        spin = array((kwargs["spin"]))
+        
+        shot = self._shoot(self.advance, y0, spin)
+        
+        return shot        
+    
+    def spin_force(self,U,spin):
+        
+        Umag = norm(U)
+        omega = norm(spin)
+        
+        Cl = self.lift_coefficient(Umag, omega)
+        
+        if U.ndim == 1:
+            f = Cl*pi*self.radius**3*environment.rho*cross(spin, U)/self.mass
+        else:
+            # Messy way to return spin array for post-processing
+            f = zeros_like(U)
+            for i in range(U.shape[1]):
+                f[:,i] = Cl*pi*self.radius**3*environment.rho*cross(spin, U[:,i])/self.mass
+        
+        return f
+    
+    def advance(self, t, vec, spin):
+        x = vec[0:3]
+        u = vec[3:6]
+        
+        Cd = self.drag_coefficient(norm(u), norm(spin))
+        
+        f = self.air_resistance_force(u, Cd) \
+          + self.gravity_force() \
+          + self.spin_force(u,spin)
+        
+        return concatenate((u,f))
+
+
+class ShotPutBall(_SphericalParticleAirResistance):
+    """
+    Note that diameter can vary 110 mm to 130mm
+    and 95 mm to 110 mm
+    """
+    def __init__(self, weight_class):
+        
+        if weight_class == 'M':
+            mass = 7.26
+            diameter = 0.11
+        elif weight_class == 'F':
+            mass = 4.0
+            diameter = 0.095
+            
+        super().__init__(mass, diameter)
+        
+        
+class SoccerBall(_SphericalParticleAirResistanceSpin):
+    """
+    Note that diameter can vary 110 mm to 130mm
+    and 95 mm to 110 mm
+    """
+    def __init__(self, mass=0.430, diameter=0.22):
+                    
+        super(SoccerBall, self).__init__(mass, diameter)
+    
+    def drag_coefficient(self, velocity, omega):
+        # Texture, sewing pattern and spin will alter
+        # the drag coefficient.
+        # Here, use correlation from
+        
+        # Goff, J. E., & Carré, M. J. (2010). Soccer ball lift 
+        # coefficients via trajectory analysis. 
+        # European Journal of Physics, 31(4), 775.
+        
+        
+        vc = 12.19
+        vs = 1.309
+        
+        S = omega*self.radius/velocity;
+        if S > 0.05 and velocity > vc:
+            Cd = 0.4127*S**0.3056;
+        else:
+            Cd = 0.155 + 0.346 / (1 + exp((velocity - vc)/vs))
+        
+        return Cd
+    
+    def lift_coefficient(self, Umag, omega):
+        # TODO - complex dependency on Re and spin, skin texture etc
+        return 0.9
+ 
+class TableTennisBall(SoccerBall):
+    """
+    
+    """
+    def __init__(self):
+        
+        mass = 2.7e-3   
+        diameter = 40e-3 
+            
+        super(TableTennisBall, self).__init__(mass, diameter)
+        
+        
+class DiscGolfDisc(_Projectile):
+    def __init__(self, name, mass=0.175):
+        this_dir = os.path.dirname(os.path.abspath(__file__))
+        path = os.path.join(this_dir, 'discs', name + '.yaml')
+    
+        self.name = name
+        
+        with open(path, 'r') as f:
+            data = yaml.load(f, Loader=yaml.FullLoader)
+            
+            self.diameter = data['diameter']
+            self.mass = mass
+            self.weight = environment.g*mass
+            self.area = pi*self.diameter**2/4.0
+            self.I_xy = mass*data['J_xy']
+            self.I_z = mass*data['J_z']
+            
+            a = array(data['alpha'])
+            cl = array(data['Cl'])
+            cd = array(data['Cd'])
+            cm = array(data['Cm'])
+            
+        
+        self._alpha,self._Cl,self._Cd,self._Cm = self._flip(a,cl,cd,cm)
+        kind = 'linear'
+        self.Cl_func = interp1d(self._alpha, self._Cl, kind=kind)
+        self.Cd_func = interp1d(self._alpha, self._Cd, kind=kind)
+        self.Cm_func = interp1d(self._alpha, self._Cm, kind=kind)
+        
+    def _flip(self,a,cl,cd,cm):
+        """
+        Data given from -90 deg to 90 deg.
+        Expand to -180 to 180 using symmetry considerations.
+        """
+        n = len(a)
+
+        idx = argmin(abs(a))
+        a2 = zeros(2*n)
+        cl2 = zeros(2*n)
+        cd2 = zeros(2*n)
+        cm2 = zeros(2*n)
+        
+        a2[idx:idx+n] = a[:]
+        cl2[idx:idx+n] = cl[:]
+        cd2[idx:idx+n] = cd[:]
+        cm2[idx:idx+n] = cm[:]
+        for i in range(idx):
+            a2[i] = -(180 + a[idx-i])
+            cl2[i] = -cl[idx-i]
+            cd2[i] =  cd[idx-i]
+            cm2[i] = -cm[idx-i]
+        
+        for i in range(idx+n,2*n):
+            a2[i] = 180 - a[idx+n-i-2]
+            cl2[i] = -cl[idx+n-i-2]
+            cd2[i] =  cd[idx+n-i-2]
+            cm2[i] = -cm[idx+n-i-2]
+        
+        return a2,cl2,cd2,cm2
+        
+    def _normalize_angle(self, alpha):
+        """
+        Ensure that the angle fulfils :math:`-\\pi < \\alpha < \\pi`
+
+        :param float alpha: Angle in radians
+        :return: Normalized angle
+        :rtype: float
+        """
+
+        return arctan2(sin(alpha), cos(alpha))
+    
+    def Cd(self, alpha): 
+        """
+        Provide drag coefficent for a given angle of attack.
+
+        :param float alpha: Angle in radians
+        :return: Drag coefficient
+        :rtype: float
+        """
+        
+        # NB! The stored data uses degrees for the angle
+        return self.Cd_func(degrees(self._normalize_angle(alpha)))
+
+    def Cl(self, alpha): 
+        """
+        Provide drag coefficent for a given angle of attack.
+
+        :param float alpha: Angle in radians
+        :return: Drag coefficient
+        :rtype: float
+        """
+        
+        # NB! The stored data uses degrees for the angle
+        return self.Cl_func(degrees(self._normalize_angle(alpha)))
+
+    def Cm(self, alpha): 
+        """
+        Provide coefficent of moment for a given angle of attack.
+
+        :param float alpha: Angle in radians
+        :return: Coefficient of moment
+        :rtype: float
+        """
+    
+        # NB! The stored data uses degrees for the angle
+        return self.Cm_func(degrees(self._normalize_angle(alpha)))
+
+
+    def plot_coeffs(self, color='k'):
+        """
+        Utility function to quickly explore disc coefficients.
+
+        :param string color: Matplotlib color key. Default value is k, i.e. black.
+        """
+        pl.plot(self._alpha, self._Cl, 'C0-o',label='$C_L$')
+        pl.plot(self._alpha, self._Cd, 'C1-o',label='$C_D$')
+        pl.plot(self._alpha, 3*self._Cm, 'C2-o',label='$C_M$')
+        
+        a = linspace(-pi,pi,200)
+        #pl.plot(degrees(a), self.Cl(a), 'C0-',label='$C_L$')
+        #pl.plot(degrees(a), self.Cd(a), 'C1-',label='$C_D$')
+        #pl.plot(degrees(a), 3*self.Cm(a), 'C2-',label='$C_M$')
+        
+        pl.xlabel('Angle of attack ($^\circ$)')
+        pl.ylabel('Aerodynamic coefficients (-)')
+        pl.legend(loc='upper left')
+        ax = pl.gca()
+        ax2 = pl.gca().twinx()
+        ax2.set_ylabel("Aerodynamic efficiency, $C_L/C_D$")
+        pl.plot(self._alpha, self._Cl/self._Cd, 'C3-.',label='$C_L/C_D$')
+        ax2.legend(loc='upper right')
+        
+        return ax,ax2
+    
+    
+    def empirical_spin(self, speed):
+        # Simple empirical formula for spin rate, based on curve-fitting
+        # data from:
+        # https://www.dgcoursereview.com/dgr/forums/viewtopic.php?f=2&t=7097
+        #omega = -0.257*speed**2 + 15.338*speed
+
+        # Alternatively, experiments indicate a linear relationship,
+        omega = 5.2*speed
+
+        return omega
+
+            
+    
+    def initialize_shot(self, **kwargs):
+        U = kwargs["speed"]
+        
+        kwargs.setdefault('yaw', 0.0) 
+        #kwargs.setdefault('omega', self.empirical_spin(U)) 
+        
+        pitch = radians(kwargs["pitch"])
+        yaw = radians(kwargs["yaw"])
+        omega = kwargs["omega"]
+        
+        # phi, theta
+        roll_angle = radians(kwargs["roll_angle"]) # phi
+        nose_angle = radians(kwargs["nose_angle"]) # theta
+        # psi, rotation around z irrelevant for starting position
+        #      since the disc is symmetric
+        
+        # Initialize position
+        if "position" in kwargs:
+            x,y,z = kwargs["position"]
+        else:
+            x = 0.
+            y = 0.
+            z = 0.
+        
+        # Initialize velocity
+        xy = cos(pitch)
+        u = U*xy*cos(yaw)
+        v = U*xy*sin(-yaw)
+        w = U*sin(pitch)
+        
+        # Initialize angles
+        attitude = array([roll_angle, nose_angle, 0])
+        # The initial orientation of the disc must also account for the
+        # angle of the throw itself, i.e. the launch angle. 
+        attitude += matmul(T_12(attitude), array((0, pitch, 0)))
+        
+        #attitude = matmul(T_23(yaw),attitude)
+        #attitude += matmul(T_12(attitude), array((0, pitch, 0)))
+        phi, theta, psi = attitude
+        y0 = array((x,y,z,u,v,w,phi,theta,psi))
+        return y0, omega
+            
+    def shoot(self, **kwargs):
+
+        y0, omega = self.initialize_shot(**kwargs)
+               
+        shot = self._shoot(self.advance, y0, omega)
+        
+        return shot
+    
+    def post_process(self, s, omega):
+        n = len(s.time)
+        alphas = zeros(n)
+        betas = zeros(n)
+        lifts = zeros(n)
+        drags = zeros(n)
+        moms = zeros(n)
+        rolls = zeros(n)
+        for i in range(n):
+            x = s.position[:,i]
+            u = s.velocity[:,i]
+            a = s.attitude[:,i]
+            
+            alpha, beta, Fd, Fl, M, g4 = self.forces(x, u, a, omega)
+            
+            alphas[i] = alpha
+            betas[i] = beta
+            lifts[i] = Fl
+            drags[i] = Fd
+            moms[i] = M
+            rolls[i] = -M/(omega*(self.I_xy - self.I_z))
+        
+        arc_length = norm(s.position, axis=0)
+        return arc_length,degrees(alphas),degrees(betas),lifts,drags,moms,degrees(rolls)
+            
+    def forces(self, x, u, a, omega):
+        # Velocity in body axes
+        urel = u - environment.wind_abl(x[2])
+        u2 = matmul(T_12(a), urel)
+        # Side slip angle is the angle between the x and y velocity
+        beta = -arctan2(u2[1], u2[0])
+        # Velocity in zero side slip axes
+        u3 = matmul(T_23(beta), u2)
+        # Angle of attack is the angle between 
+        # vertical and horizontal velocity
+        alpha = -arctan2(u3[2], u3[0])
+        # Velocity in wind system, where forces are to be calculated
+        u4 = matmul(T_34(alpha), u3)
+        
+        # Convert gravitational force from Earth to Wind axes
+        g = array((0, 0, self.mass*environment.g))
+        g4 = T_14(g, a, beta, alpha)
+        
+        # Aerodynamic forces
+        q = 0.5*environment.rho*u4[0]**2
+        S = self.area
+        D = self.diameter
+        
+        Fd = q*S*self.Cd(alpha)
+        Fl = q*S*self.Cl(alpha)
+        M  = q*S*D*self.Cm(alpha)
+        
+        return alpha, beta, Fd, Fl, M, g4
+        
+    def advance(self, t, vec, omega):
+        x = vec[0:3]
+        u = vec[3:6]
+        a = vec[6:9]
+        
+        alpha, beta, Fd, Fl, M, g4 = self.forces(x, u, a, omega)
+        
+        m = self.mass
+        # Calculate accelerations
+        dudt = (-Fd + g4[0])/m
+        dvdt =        g4[1]/m
+        dwdt = ( Fl + g4[2])/m
+        acc4 = array((dudt,dvdt,dwdt))
+        # Roll rate acts around x-axis (in axes 3: zero side slip axes)
+        dphidt = -M/(omega*(self.I_xy - self.I_z))
+        # Other angular rotations are ignored, assume zero wobble
+        angvel3 = array((dphidt, 0, 0))
+        
+        acc1 = T_41(acc4, a, beta, alpha)
+        angvel1 = T_31(angvel3, a, beta)
+        
+        return concatenate((u,acc1,angvel1)) 
+
+    
+    
+
+    

shotshaper/transforms.py

@@ -0,0 +1,68 @@
+# -*- coding: utf-8 -*-
+"""
+
+The axes are denoted:
+1: Earth
+2: Body
+3: Zero side slip - rotation around z, to apply roll
+4: Wind axes - to replicate the CFD condition with wind coming straight onto the body
+
+In Crowther & Potts, z-axis is defined pointing downwards. Here, it points upwards
+leading to some different signs.
+"""
+import numpy as np
+from numpy import cos,sin,matmul
+
+def T_12(attitude):
+    """
+    Transform from Earth axes to Body axes
+    """
+    phi, theta, psi = attitude
+                                                                    
+    return np.array([[cos(theta)*cos(psi), sin(phi)*sin(theta)*cos(psi) - cos(phi)*sin(psi), cos(phi)*sin(theta)*cos(psi) + sin(phi)*sin(psi)],
+                     [cos(theta)*sin(psi), sin(phi)*sin(theta)*sin(psi) + cos(phi)*cos(psi), cos(phi)*sin(theta)*sin(psi) - sin(phi)*cos(psi)],
+                     [-sin(theta),         sin(phi)*cos(theta),                              cos(phi)*cos(theta)                             ]])
+
+def T_23(beta):
+    """
+    Transform from Body axes to Zero side slip axes,
+    Rotation around z-axis by the side-slip angle
+    """
+    return np.array([[cos(beta), -sin(beta), 0],
+                     [sin(beta),  cos(beta), 0],
+                     [0,          0,         1]])
+
+
+def T_34(alpha):
+    """
+    Transform from Zero side slip axes to Wind axes
+    Rotation around y-axis by the angle of attack.
+    """
+    return np.array([[cos(alpha), 0, -sin(alpha)],
+                     [0,          1,  0         ],
+                     [sin(alpha), 0,  cos(alpha)]])
+
+def T_14(vec, attitude, beta, alpha):
+    return matmul(T_34(alpha), matmul(T_23(beta), matmul(T_12(attitude), vec)))
+
+def T_21(attitude):
+    """
+    Transform from Body axes to Earth axes.
+    Done by transposing the opposite 
+    """
+    return np.transpose(T_12(attitude))
+
+def T_32(beta):
+    return np.transpose(T_23(beta))
+
+def T_43(alpha):
+    return np.transpose(T_34(alpha))
+
+def T_41(vec, attitude, beta, alpha):
+     return matmul(T_21(attitude), matmul(T_32(beta), matmul(T_43(alpha), vec)))
+
+def T_31(vec, attitude, beta):
+    return matmul(T_21(attitude), matmul(T_32(beta), vec))
+
+
+                   

utils/disc_geometric_properties.py

@@ -0,0 +1,28 @@
+# -*- coding: utf-8 -*-
+"""
+Calculate moments of inertia for a disc from STL file.
+"""
+
+import numpy as np
+import trimesh
+import sys
+import os 
+
+path = os.path.dirname(os.path.realpath(__file__))
+# attach to logger so trimesh messages will be printed to console
+#trimesh.util.attach_to_log()
+
+name = sys.argv[-1]
+
+m = trimesh.load(os.path.join(path, 'discs', name + '.stl'))
+trimesh.repair.fix_inversion(m)
+trimesh.repair.fix_normals(m)
+trimesh.repair.fix_winding(m)
+
+if m.is_watertight and m.is_winding_consistent and m.is_volume:
+    V = m.volume
+    J = m.principal_inertia_components/V
+    print('Volume: ', V)
+    print('J_xy: %4.3e' % J[0])
+    print('J_z: %4.3e' % J[2])
+