Untitled

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# kastparabel, ger refernslinje i grafen,


def kastparabel(x, v0, grader):
    g = 9.81
    alpha = np.radians(grader)
    y = x * np.tan(alpha) - (g * x**2) / (2 * v0**2 * np.cos(alpha)**2)
    return y

# rörelseekvationer, första ordningens diff, ode använder följnade funktion flera gånger


def derivator(t, s, m, Cd, A, rho, g, k, c, mu, vvind_x, vvind_y):
    x = s[0]
    y = s[1]
    vx = s[2]
    vy = s[3]

    # boll i luften tyngdkraft och luftmotstånd
    if y > 0:
        Ftyngd_y = -m * g  # tyngdkraften

        # relativ hastighet enligt docs
        vrel_x = vx - vvind_x
        vrel_y = vy - vvind_y
        vrel_storlek = np.sqrt(vrel_x**2 + vrel_y**2)

        Fd = 0.5 * Cd * A * rho * vrel_storlek**2  # storhet på luftmotstånd

        # riktning för luftmostånd
        if vrel_storlek > 0:
            Fluft_x = Fd * (-vrel_x / vrel_storlek)
            Fluft_y = Fd * (-vrel_y / vrel_storlek)
        else:  # undviker div på 0
            Fluft_x = 0.0
            Fluft_y = 0.0

        # acc a = F/m
        ax = Fluft_x / m
        ay = (Ftyngd_y + Fluft_y) / m

    # bollen studsar, fjäder och dämpare verkar uppåt och rullmotstånd horisontellt
    else:
        # fjäder och dämpare, fs=k*(-y) + c*(-vy), -vy bromsar studsen medans -y positiv bollen trycks ned kraft uppåt
        Fs = k * (-y) + c * (-vy)
        # rullmoståndet i x-led, vx styr riktning, (bromsande)
        Fr = mu * m * g * vx

        ax = -Fr / m
        ay = (Fs / m) - g

    return [vx, vy, ax, ay]

# Numerisk lösning


def sim_golfboll(v0, grader, vvind_x=0.0, vvind_y=0.0):

    # parametarar
    m = 0.046  # massan
    Cd = 0.25  # luftmotståndkoefficient
    A = 1430e-6  # area
    rho = 1.225  # densitet på luft
    g = 9.81  # tyngdacceleration
    k = 1000.0  # styvhet på fjäder
    c = 2.0  # dämpningskonstant
    mu = 0.15  # rullfriktionskoefficient

    # begynnelsevilkoren
    alpha = np.radians(grader)

    s0 = [0.0,
          0.001,              # y0 = 0.001 så studsar ej vid start, sedan vx0, vy0
          v0 * np.cos(alpha),
          v0 * np.sin(alpha)]

    # simulering 15 sekunder
    t_eval = np.linspace(0, 15, 1000)

    resultat = solve_ivp(
        fun=derivator,
        t_span=(0, 15),
        y0=s0,
        t_eval=t_eval,
        args=(m, Cd, A, rho, g, k, c, mu, vvind_x, vvind_y),
        method='RK45',
        max_step=0.01
    )

    return resultat.y[0], resultat.y[1], resultat.y[2], resultat.y[3]


# simulering och analys
def main():

    print("Golfboll-simulering  –  MMS135\n")
    # första analys, v0=50, a=15, ingen vind
    v0 = 50.0
    grader = 15.0

    x, y, vx, vy = sim_golfboll(v0, grader)

    # här används kastparabeln som referns
    x_kast = np.linspace(0, 220, 500)
    y_kast = kastparabel(x_kast, v0, grader)
    y_kast = np.where(y_kast < 0, np.nan, y_kast)

    # stoppdelen, där bollen är helt stilla på marken
    x_slut = x[-1]
    print(f"Första analys (a=15 grader, utan vind ger resultat: )")
    print(f" x = {x_slut:.1f} m\n (bollen stannat)")

    # Graf 1, banan i xy, och hastigheter
    fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8))
    fig.suptitle(f"Första analys v0={v0} m/s, a={grader} grader, utan vind")

    ax1.plot(x, y, 'b-', linewidth=1.5,
             label='Simulerad bana med luftmotståndet')
    ax1.plot(x_kast, y_kast, 'k--', linewidth=1.2,
             label='kastparabel som referns')
    ax1.set_xlabel("x (m)")
    ax1.set_ylabel("y (m)")
    ax1.set_title("Bollens bana i xy-planet")
    ax1.legend()
    ax1.set_ylim(bottom=-0.5)
    ax1.set_xlim(left=0)
    ax1.grid(True, alpha=0.3)

    ax2.plot(x, vx, 'b-', linewidth=1.5, label='vₓ (horisontell)')
    ax2.plot(x, vy, 'r-', linewidth=1.5, label='vᵧ (vertikal)')
    ax2.axhline(0, color='black', linewidth=0.8, linestyle='--')
    ax2.set_xlabel("x (m)")
    ax2.set_ylabel("Hastighet (m/s)")
    ax2.set_title("Bollens hastigheter som funktion av x")
    ax2.legend()
    ax2.set_xlim(left=0)
    ax2.grid(True, alpha=0.3)

    plt.tight_layout()
    plt.savefig("grundfall.png", dpi=120)
    plt.show()

    # analys del 2 för den optimala utslagsvinkeln (resultat printas)
    print("Andra analys, optimal utslagsvinkel (v0=50 m/s, utan vind)")

    vinklar = np.arange(5, 65, 1)
    slutpunkter = []

    for alfa in vinklar:
        x_i, y_i, vx_i, vy_i = sim_golfboll(50.0, alfa)
        slutpunkter.append(x_i[-1])

    slutpunkter = np.array(slutpunkter)
    basta_index = np.argmax(slutpunkter)
    print(
        f"Längsta skottet: {slutpunkter[basta_index]:.1f} m vid a = {vinklar[basta_index]} grader\n")

    # Analys del 3, 5m/s som luftmoståendet, andra värden maneullt valda efter testning
    print("Analys 3, hole-in-one (5 m/s luftmotstånd, hål x=190 m)")

    v0_val = 50.0
    alpha_val = 39.0

    x_b, y_b, vx_b, vy_b = sim_golfboll(v0_val, alpha_val, vvind_x=-5.0)

    # kontroll av hastighten och höjden vid hålet
    idx = np.argmin(np.abs(x_b - 190.0))
    vx_hal = abs(vx_b[idx])
    y_hal = abs(y_b[idx])

    print(f"v0={v0_val} m/s, a={alpha_val} grader")
    print(f"|vx| x=190 m: {vx_hal:.2f} m/s ,måste vara lägre än 0.2m/s")
    print(f"vid x=190 m:  {y_hal:.3f} m  , måste vara nära 0")

    if vx_hal < 0.2 and y_hal < 0.1:
        print("hole in one lyckat")
    else:
        print("justera parameterar")

    # graf för hole in one
    fig2, ax = plt.subplots(figsize=(12, 5))
    ax.plot(x_b, y_b, 'b-', linewidth=1.5,
            label=f'v0={v0_val} m/s, a={alpha_val} grader')
    ax.set_xlabel("x (m)")
    ax.set_ylabel("y (m)")
    ax.set_title(
        f"Hole-in-one med 5 m/s motvind  (v0={v0_val} m/s, a={alpha_val})")
    ax.legend()
    ax.set_ylim(bottom=-0.5)
    ax.set_xlim(left=0)
    ax.grid(True, alpha=0.3)
    plt.tight_layout()
    plt.savefig("hole_in_one.png", dpi=120)
    plt.show()


if __name__ == "__main__":
    main()