Scipy - solve_ivp - 解决角度的弹道问题

如何解决Scipy - solve_ivp - 解决角度的弹道问题

我需要描述飞碟的投掷。投掷从坐标 x0 = 60 和 y0 = 0 开始，速度 = 40 ms^-1。也有相反方向的风，风速为 3.5 ms^-1，并且飞碟在 x 和 y 方向具有不同的阻力。

picture - skeet throw

现在我需要计算角度，使飞碟降落在坐标 x=0 和 y=0 上。我的问题是我不知道如何在没有定义角度的情况下进行整合。我希望有人可以帮助我。

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

v_wind = 3.5 #m / s in positive x direction
rho_l = 1.2 # kg/m3

# skeet params
cw_skeet_top = 1.11
cw_skeet_side = 0.5
A_skeet_top = 0.0086
A_skeet_side = 0.00432
m_skeet = 0.105 #kg

v0_skeet = 40
x0_skeet = 60
y0_skeet = 0

# rocket params
cw_rocket_side = 1
cw_rocket_front = 0.5
A_rocket_side = 0.08
A_rocket_front = 0.0312
m_rocket = 2

v0_rocket = 0
x0_rocket = 20
y0_rocket = 0
F_accel_rocket = 120 #in Newton,divide by m_rocket to find accel in m/s^2
t_accel_end = 0.5

# grenade parameters
cw_grenade = 0.5
A_grenade = np.pi*0.1**2
m_grenade = 1 #kg
t_boom = 5

# Function that calculates drag force as a function of relative veLocity and c_w
def drag(v_rel,rho_l,cw,A):
    Fw = rho_l * cw * A * v_rel**2/2.0
    if (v_rel < 0):
        Fw = -Fw
    return Fw
def skeet(t,y,A,m):
    # Zur Übersichtlichkeit:
    # y[0] - x-position
    # y[1] - x-veLocity
    # y[2] - y-position
    # y[3] - y-veLocity

    d2xdt2 = -drag(y[1]-v_wind(t),cw_skeet_side,A_skeet_side)/m_skeet   
    dxdt = y[1]    
    d2ydt2 = -9.81 - drag(y[3],cw_skeet_top,A_skeet_top)/m_skeet
    dydt = y[3]
    
    return dxdt,d2xdt2,dydt,d2ydt2

def hit_ground(t,*args):
    return y[2]
hit_ground.terminal = True
hit_ground.direction = -1

x0 = 60 #m
y0 = 0
v0 = 40 #meters per second
angle = 
tspan = [0,20]
maxStep=tspan[1]/1000 

start_vec = [x0,v0*np.cos(angle/180*np.pi),y0,v0*np.sin(angle/180*np.pi)]
sol = solve_ivp(skeet,tspan,start_vec,args=(rho_l,m),max_step=maxStep,events=(hit_ground))