Plots for https://lisakov.com/blog/air-resistance/

air-length.py

import numpy as np
from sympy import *
from scipy import interpolate
import matplotlib.pyplot as plt
from matplotlib import rc, rcParams

# Plot and save trajectories y(x) for fixed v0, k and varying
# angle α.
# Last line prints an angle corresponding to the maximum length


# Triginometry {{{
def cos(x): return round( np.cos( np.radians(x) ), 10 )
def sin(x): return round( np.sin( np.radians(x) ), 10 )
def sins(x): return round( np.sin( np.radians(x) )**2, 10 )
def asin(x): return np.degrees( np.arcsin(x) )
def tg(x): return round( np.tan( np.radians(x) ), 10 )
def atg(x): return np.degrees( np.arctan(x) )
#}}}
# Plot setup {{{
rcParams['font.size'] = 16.
rcParams['text.usetex'] = True
#rcParams['font.sans-serif'] = ['Computer Modern']
rc('text.latex', preamble=r'\usepackage[T2A]{fontenc} \usepackage[utf8]{inputenc} \usepackage{bm}')

rc('axes', linewidth=1.2)
rc('xtick.major', size=4, width=1.1)
rc('ytick.major', size=4, width=1.1)
fig, ax = plt.subplots()
fig.subplots_adjust(left=.15, bottom=.15, right=.95, top=.85)
#}}}

v0 = 20
m  = 1
k  = 0.3
g  = 10
e  = np.exp(1)

times = np.linspace(0,5,1000)
angles = np.arange(0.0001, 90, 0.1) # Set larger step instead of 0.1 to reduce images number

ls=[]
for a in angles:
    # Sans l'air
    x=[]
    for t in times:
        x = np.append( x, v0*cos(a)*t)
    y=[]
    for t in times:
        y = np.append( y, v0*sin(a)*t - g*t**2/2)
    ax.plot(x, y, ls=':', label='$F_\\textrm{сопр}=0$')

    # Avec l'air
    x=[]
    for t in times:
        x = np.append( x, -v0*cos(a)*m/k * (np.exp(-k*t/m)-1) )
    y=[]
    for t in times:
        y = np.append( y, -(v0*sin(a)+m/k*g) * m/k * (np.exp(-k*t/m)-1) - g*m/k*t )
    ax.plot(x, y, label='$F_\\textrm{сопр}=kv$')

    ax.set_xlabel('$x$ [\\textrm{м}]')
    ax.set_ylabel('$y$ [\\textrm{м}]')

    ax.set_xlim([0,42])
    ax.set_ylim([0,21])

    ax.plot([], [], ' ', label='$k={}~$'.format( "{:.1f}".format(k) )+'\\textrm{кг/с}')
    ax.plot([], [], ' ', label='$m={}~$'.format( "{:.0f}".format(m) )+'\\textrm{кг}')
    ax.plot([], [], ' ', label='$v_0={}~$'.format( "{:.0f}".format(v0) )+'\\textrm{м/с}')
    ax.plot([], [], ' ', label='$\\alpha={}$'.format( "{:.1f}".format(a) )+'$^\circ$')
    ax.legend(markerscale=0, handletextpad=0.4, loc='center',
                bbox_to_anchor=(0.5, 1.15), ncol=3, columnspacing=0.0)
    ax.grid()

    #plt.savefig('length/' + str("{:02.3f}".format(a))+'.png', bbox_inches='tight')
    plt.cla()

    # Find numerically the time of flight (tau) when y=0
    t, y = symbols('t y') 
    y    =  -(v0*sin(a)+m/k*g) * m/k * (e**(-k*t/m)-1) - g*m/k*t
    tau  = nsolve(Eq(y, 0), 2) # time of flight
    def x(t):
        return -v0*cos(a)*m/k * (e**(-k*t/m)-1) #np.exp won't work here
    l = x(tau)
    ls = np.append (ls, l)

# Interpolate angles and find the angle corresponding to the maximum distance
# Not sure that interpolation is needed
f = interpolate.interp1d(ls, angles, kind='linear', bounds_error=False, fill_value=0.0)

for l,a in zip(ls,angles):
    print ("{:>8.2f}".format(a), "{:>7.6}".format(str(l)))

print ( 'Coefficient k is', k )
print ( 'Max distance is', np.amax(ls) )
print ( 'Angle, corresponding to the maximum distance, is', f(float(np.amax(ls))) )

### angle(k) corresponding to maximum length
#  k    angle
#  0.1  42.5
#  0.2  40.5
#  0.3  38.6
#  0.4  36.9
#  0.5  35.7

air-x.py

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rc, rcParams

# Plot and save x(t) for fixed v0, α and varying
# resistance coefficient k.


# Triginometry {{{
def cos(x): return round( np.cos( np.radians(x) ), 10 )
def sin(x): return round( np.sin( np.radians(x) ), 10 )
def sins(x): return round( np.sin( np.radians(x) )**2, 10 )
def asin(x): return np.degrees( np.arcsin(x) )
def tg(x): return round( np.tan( np.radians(x) ), 10 )
def atg(x): return np.degrees( np.arctan(x) )
#}}}
# Plot setup {{{
rcParams['font.size'] = 16.
rcParams['text.usetex'] = True
rcParams['font.sans-serif'] = ['Computer Modern']
rc('text.latex', preamble=r'\usepackage[T2A]{fontenc} \usepackage[utf8]{inputenc} \usepackage{bm}')

rc('axes', linewidth=1.2)
rc('xtick.major', size=4, width=1.1)
rc('ytick.major', size=4, width=1.1)
fig, ax = plt.subplots()
fig.subplots_adjust(left=.15, bottom=.15, right=.95, top=.75)
#}}}

v0 = 20
a  = 30
m  = 1
#k  = 0.4
g  = 10

times = np.linspace(0,10,1000)
ks = np.arange(0.0001,1.1,0.005) # Set larger step instead of 0.005 to reduce images number

for k in ks:
    x=[]
    for t in times:
        x = np.append( x, v0*cos(a)*t )
    ax.plot(times, x, ls=':', label='$F_\\textrm{сопр}=0$')

    x=[]
    for t in times:
        x = np.append( x, -v0*cos(a)*m/k * (np.exp(-k*t/m)-1) )
    ax.plot(times, x, label='$F_\\textrm{сопр}=kv$')

    ax.set_xlabel('$t$ [\\textrm{с}]')
    ax.set_ylabel('$x$ [\\textrm{м}]')

    ax.plot([], [], ' ', label='$k={}~$'.format( "{:.2f}".format(k) )+'\\textrm{кг/с}')
    ax.plot([], [], ' ', label='$m={}~$'.format( "{:.0f}".format(m) )+'\\textrm{кг}')
    ax.plot([], [], ' ', label='$v_0={}~$'.format( "{:.0f}".format(v0) )+'\\textrm{м/с}')
    ax.plot([], [], ' ', label='$\\alpha={}$'.format( "{:.0f}".format(a) )+'$^\circ$')
    ax.legend(markerscale=0, handletextpad=0.4, loc='center', bbox_to_anchor=(0.5, 1.15), ncol=3, columnspacing=0.0)
    ax.grid()

    plt.savefig('x/' + str(k)+'.png', bbox_inches='tight')
    plt.cla()

air-y.py

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rc, rcParams

# Plot and save y(t) for fixed v0, α and varying
# resistance coefficient k.

# Triginometry {{{
def cos(x): return round( np.cos( np.radians(x) ), 10 )
def sin(x): return round( np.sin( np.radians(x) ), 10 )
def sins(x): return round( np.sin( np.radians(x) )**2, 10 )
def asin(x): return np.degrees( np.arcsin(x) )
def tg(x): return round( np.tan( np.radians(x) ), 10 )
def atg(x): return np.degrees( np.arctan(x) )
#}}}
# Plot setup {{{
rcParams['font.size'] = 16.
rcParams['text.usetex'] = True
rc('text.latex', preamble=r'\usepackage[T2A]{fontenc} \usepackage[utf8]{inputenc} \usepackage{bm}')

rc('axes', linewidth=1.2)
rc('xtick.major', size=4, width=1.1)
rc('ytick.major', size=4, width=1.1)
fig, ax = plt.subplots()
fig.subplots_adjust(left=.15, bottom=.15, right=.95, top=.85)
#}}}

v0 = 20
a  = 60
m  = 1
#k  = 0.4
g  = 10

times = np.linspace(0,5,1000)
ks = np.arange(0.0001,1.1,0.005) # Set larger step instead of 0.005 to reduce images number

for k in ks:
    y=[]
    for t in times:
        y = np.append( y, v0*sin(a)*t - g*t**2/2)
    ax.plot(times, y, ls=':', label='$F_\\textrm{сопр}=0$')

    y=[]
    for t in times:
        y = np.append( y, -(v0*sin(a)+m/k*g) * m/k * (np.exp(-k*t/m)-1) - g*m/k*t )
    ax.plot(times, y, label='$F_\\textrm{сопр}=kv$')

    ax.set_xlabel('$t$ [\\textrm{с}]')
    ax.set_ylabel('$y$ [\\textrm{м}]')

    ax.plot([], [], ' ', label='$k={}~$'.format( "{:.2f}".format(k) )+'\\textrm{кг/с}')
    ax.plot([], [], ' ', label='$m={}~$'.format( "{:.0f}".format(m) )+'\\textrm{кг}')
    ax.plot([], [], ' ', label='$v_0={}~$'.format( "{:.0f}".format(v0) )+'\\textrm{м/с}')
    ax.plot([], [], ' ', label='$\\alpha={}$'.format( "{:.0f}".format(a) )+'$^\circ$')
    ax.legend(markerscale=0, handletextpad=0.4, loc='center',
                bbox_to_anchor=(0.5, 1.15), ncol=3, columnspacing=0.0)
    ax.grid()

    plt.savefig('y/' + str(k)+'.png', bbox_inches='tight')
    plt.cla()

air-yx.py

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rc, rcParams

# Plot and save trajectory y(x) for fixed v0, α and varying
# resistance coefficient k.

# Triginometry {{{
def cos(x): return round( np.cos( np.radians(x) ), 10 )
def sin(x): return round( np.sin( np.radians(x) ), 10 )
def sins(x): return round( np.sin( np.radians(x) )**2, 10 )
def asin(x): return np.degrees( np.arcsin(x) )
def tg(x): return round( np.tan( np.radians(x) ), 10 )
def atg(x): return np.degrees( np.arctan(x) )
#}}}
# Plot setup {{{
rcParams['font.size'] = 16.
rcParams['text.usetex'] = True
#rcParams['font.sans-serif'] = ['Computer Modern']
rc('text.latex', preamble=r'\usepackage[T2A]{fontenc} \usepackage[utf8]{inputenc} \usepackage{bm}')

rc('axes', linewidth=1.2)
rc('xtick.major', size=4, width=1.1)
rc('ytick.major', size=4, width=1.1)
fig, ax = plt.subplots()
fig.subplots_adjust(left=.15, bottom=.15, right=.95, top=.85)
#}}}

v0 = 20
a  = 45
#k = 0.3
m  = 1
g  = 10

times = np.linspace(0,5,1000)
ks = np.arange(0.0001,1.1,0.005) # Set larger step instead of 0.005 to reduce images number

for k in ks:
    # Air resistance off
    x=[]
    for t in times:
        x = np.append( x, v0*cos(a)*t)
    y=[]
    for t in times:
        y = np.append( y, v0*sin(a)*t - g*t**2/2)
    ax.plot(x, y, ls=':', label='$F_\\textrm{сопр}=0$')

    # Air resistance on
    x=[]
    for t in times:
        x = np.append( x, -v0*cos(a)*m/k * (np.exp(-k*t/m)-1) )
    y=[]
    for t in times:
        y = np.append( y, -(v0*sin(a)+m/k*g) * m/k * (np.exp(-k*t/m)-1) - g*m/k*t )
    ax.plot(x, y, label='$F_\\textrm{сопр}=kv$')

    ax.set_xlabel('$x$ [\\textrm{м}]')
    ax.set_ylabel('$y$ [\\textrm{м}]')

    ax.plot([], [], ' ', label='$k={}~$'.format( "{:.2f}".format(k) )+'\\textrm{кг/с}')
    ax.plot([], [], ' ', label='$m={}~$'.format( "{:.0f}".format(m) )+'\\textrm{кг}')
    ax.plot([], [], ' ', label='$v_0={}~$'.format( "{:.0f}".format(v0) )+'\\textrm{м/с}')
    ax.plot([], [], ' ', label='$\\alpha={}$'.format( "{:.0f}".format(a) )+'$^\circ$')
    ax.legend(markerscale=0, handletextpad=0.4, loc='center',
                bbox_to_anchor=(0.5, 1.15), ncol=3, columnspacing=0.0)
    ax.grid()

    plt.savefig('yx/' + str(k)+'.png', bbox_inches='tight')
    plt.cla()