Pythons snippets computing the time needed for a bowling ball to reach the bottom of the mariana trench (XKCD)

BowlingBall.py

#!/usr/bin/python
# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt

#Fonction donnat la friction pour un reynolds donné
re_friction = lambda x:(np.sqrt(24/x)+0.5406)**2

#Fonction trouvant la vélocité d'une bulle de diamètre X dans
#l'eau

def find_velocity(weight): #weight in kg

    #Définition des variables nécessaires
    grav = 9.81 # m/sec²
    diam = 0.026 # m
    density_water = 1.000*(10**3) # kg/m³
    density_air   = weight/((4*3.1416*(diam/2)**3)/3) # kg/m³
    print(density_air)
    viscosity_water = 1*(10**-3) # Pascal / sec


    #Définition des fonctions utilisé lors
    #de la méthode ittérative
    #Fonction donnant la vitesse a partir d'un reynolds
    reynolds= lambda v:((diam*v*density_water)/ viscosity_water)
  
    #Fonction donnant le reynolds a partir d'une vitesse
    velocity_rey = lambda rey:np.sqrt((4*grav*diam*np.abs((density_air-density_water)))/(3*density_water*re_friction(rey)))

    #Vélocité initiale (n'as pas vraiment d'importance)
    velocity = velocity_rey(reynolds(0.001))
    
    #Méthode ittérative
    precision=50
    itteration=0
    while itteration < precision:
        velocity = velocity_rey(reynolds(velocity))
        itteration+=1

    return [velocity,reynolds(velocity)]
    
#Création d'un tableau de valeur entre 0.1 mm et 1 cm
d = np.arange(6,120,0.1) # cm, cm, cm/sép

#Trouver la vélocité pour chaque éléments du tableau
y_axe=find_velocity(d)[0] # m/sec²
#Trouver le reynolds pour chaque éléments du tableau
z_axe=find_velocity(d)[1] # nSu


#Paramètrer l'affichage
fig = plt.figure()
ax = fig.add_subplot(2,1,1)
bx = fig.add_subplot(2,1,2)

ax.set_title("Speed of sinking a bowling ball and its time \n to hit the bottom of the mariana trench")
ax.set_xlabel("Mass (kg)")
ax.set_ylabel("Speed (m/sec)")
ax.grid(True)
bx.set_xlabel("Mass (kg)")
bx.set_ylabel("Time (min)")
bx.grid(True)

find_time=lambda x:(11000/x)/60
name=["regular","gold","osmium"]
it=0
for pos in [7,62,105]:
	weight=pos
	speed=find_velocity(pos)[0]
	time=find_time(find_velocity(pos)[0])
	ax.annotate("{}\n{}kg , {:.1f} m/sec".format(name[it],pos,speed),xy=(pos,speed))
	bx.annotate("{}\n{}kg , {:.1f} min".format(name[it],pos,time),xy=(pos,time))
	it+=1

p = ax.plot(d,y_axe)
p = bx.plot(d,find_time(y_axe))

#Montrer le graphique
plt.show()