비 올 때 뛰면 덜 젖을까? 같은 100 m를 걷고 뛰는 빗물 모형

"""Fixed-distance rain interception, not clothing absorption or human measurements.
Python 3.12; numpy, matplotlib, Pillow. MP4 additionally needs ffmpeg on PATH.
Run: python rain_experiment.py. Outputs are written beside this script.
"""
from pathlib import Path
import json, sys, shutil
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation, FFMpegWriter
from matplotlib import font_manager

OUT = Path(__file__).resolve().parent
# English labels keep the script portable without requiring a Korean font.
D = 100.0  # m, same trip for every speed
WIDTH, DEPTH, HEIGHT = 0.45, 0.25, 1.70  # m, illustrative cuboid
AF, AS, AT = WIDTH*HEIGHT, DEPTH*HEIGHT, WIDTH*DEPTH  # m^2
R = 10.0  # mm/hour of vertical rainfall, illustrative, not measured
W = 5.0  # m/s downward rain speed, fixed for every scenario
RHO = 1000.0*(R/1000.0/3600.0)/W  # kg of liquid water per m^3 of air

def parts(v, ux=0.0, uy=0.0):
    v = np.asarray(v, dtype=float)
    if np.any(v <= 0):
        raise ValueError('Travel speed must be positive.')
    top = RHO*AT*W*D/v
    front_back = RHO*AF*np.abs(v-ux)*D/v
    side = RHO*AS*abs(uy)*D/v
    return top, front_back, side

def mass(v, ux=0.0, uy=0.0):
    return sum(parts(v, ux, uy))

def face_integral(v, ux, uy, n):
    # Independently sum inward flux on all six faces, then midpoint time steps.
    normals = np.array([[1,0,0],[-1,0,0],[0,1,0],[0,-1,0],[0,0,1],[0,0,-1]])
    areas = np.array([AF,AF,AS,AS,AT,AT])
    relative = np.array([ux-v,uy,-W])
    rate = RHO*np.sum(areas*np.maximum(0.0,-normals@relative))
    return np.sum(np.full(n,rate*D/v/n))

plt.rcParams.update({'font.size':12,'axes.spines.top':False,'axes.spines.right':False,
                     'figure.facecolor':'#faf8f3','axes.facecolor':'#faf8f3','savefig.facecolor':'#faf8f3'})
COLORS=['#147d92','#d26643','#7161a8','#498663']
v = np.linspace(0.5,6,551)
fig, ax = plt.subplots(figsize=(9,5.5),layout='constrained')
for (label,ux,uy),color in zip([('No wind',0,0),('Headwind: -3 m/s',-3,0),('Tailwind: +3 m/s',3,0)],COLORS):
    ax.plot(v,mass(v,ux,uy)*1000,label=label,color=color,lw=3)
ax.set(xlabel='Travel speed (m/s)',ylabel='Intercepted water (g)',title='Same 100 m trip, different rain directions',ylim=(0,400))
ax.legend(frameon=False);ax.grid(alpha=.18)
fig.savefig(OUT/'rain-speed-curves.png',dpi=170);plt.close(fig)
fig,ax=plt.subplots(figsize=(9,5.5),layout='constrained')
speeds=np.array([1.5,3.0,5.0]);top,front,side=parts(speeds)
ax.bar(['1.5 m/s','3 m/s','5 m/s'],front*1000,label='Front: same swept volume',color=COLORS[0])
ax.bar(['1.5 m/s','3 m/s','5 m/s'],top*1000,bottom=front*1000,label='Top: less time in rain',color=COLORS[1])
for i,m in enumerate(mass(speeds)*1000):ax.text(i,m+2,f'{m:.1f} g',ha='center')
ax.set(ylabel='Intercepted water (g)',title='Vertical rain: faster reduces the top contribution',ylim=(0,80));ax.legend(frameon=False);ax.grid(axis='y',alpha=.18)
fig.savefig(OUT/'rain-top-front.png',dpi=170);plt.close(fig)
fig,ax=plt.subplots(figsize=(9,5.5),layout='constrained')
for uy,color in zip([0,5],COLORS):ax.plot(v,mass(v,3,uy)*1000,lw=3,label=f'Tailwind +3 m/s; crosswind {uy} m/s',color=color)
ax.axvline(3,color='#555555',ls=':',label='Speed matches tailwind')
ax.set(xlabel='Travel speed (m/s)',ylabel='Intercepted water (g)',title='A crosswind can change the best speed',ylim=(0,550));ax.legend(frameon=False);ax.grid(alpha=.18)
fig.savefig(OUT/'rain-crosswind.png',dpi=170);plt.close(fig)

errors=[]
for speed in [0.5,1.5,3,5,6]:
    for ux,uy in [(0,0),(-3,0),(3,0),(3,5)]:
        for n in [10,100,1000]:errors.append(abs(face_integral(speed,ux,uy,n)-float(mass(speed,ux,uy))))
assert max(errors)<1e-12
assert np.all(np.diff(mass(v,0,0))<0)
assert np.all(np.diff(mass(v,-3,0))<0)
assert abs(v[np.argmin(mass(v,3,0))]-3)<1e-10
assert np.all(np.diff(mass(v,3,5))<0)
assert np.allclose(parts(speeds)[1],RHO*D*AF)
try: mass(0)
except ValueError: pass
else: raise AssertionError('Zero speed was not rejected')
summary={'python':sys.version.split()[0],'numpy':np.__version__,'matplotlib':matplotlib.__version__,
         'distance_m':D,'rain_mm_hr':R,'rain_down_m_s':W,'water_in_air_kg_m3':RHO,
         'areas_m2':{'front':AF,'side':AS,'top':AT},'max_abs_error_kg':max(errors),
         'cases':[{'speed_m_s':float(s),'time_s':D/s,'no_wind_g':float(mass(s)*1000),
                   'headwind_g':float(mass(s,-3)*1000),'tailwind_g':float(mass(s,3)*1000)} for s in speeds]}
(OUT/'results.json').write_text(json.dumps(summary,indent=2),encoding='utf-8')
print(json.dumps(summary,indent=2))
# Animation is a time-compressed accumulation chart, not a measured rain movie.
if shutil.which('ffmpeg'):
    fig,axes=plt.subplots(1,2,figsize=(10,4.5),layout='constrained')
    tmax=D/1.5; ts=np.linspace(0,tmax,121)
    lines=[]
    for ax,title,ylab,limit in zip(axes,['Distance travelled','Water intercepted, vertical rain'],['Distance (m)','Water (g)'],[105,75]):
        ax.set(xlim=(0,tmax),ylim=(0,limit),xlabel='Elapsed time (s)',ylabel=ylab,title=title);ax.grid(alpha=.2)
        for s,c in zip(speeds,COLORS):lines.append(ax.plot([],[],color=c,lw=3,label=f'{s:g} m/s')[0])
        ax.legend(frameon=False,loc='upper left')
    def update(k):
        t=ts[:k+1]
        for j,s in enumerate(speeds):
            active=np.minimum(t,D/s)
            lines[j].set_data(t,s*active)
            lines[3+j].set_data(t,float(mass(s))*1000*active/(D/s))
        return lines
    ani=FuncAnimation(fig,update,frames=len(ts),interval=100,blit=False)
    ani.save(OUT/'rain-fixed-distance.mp4',writer=FFMpegWriter(fps=10,codec='libx264',extra_args=['-pix_fmt','yuv420p']))
    update(120);fig.savefig(OUT/'video-final-frame.png',dpi=130);plt.close(fig)
else:
    print('MP4 skipped: ffmpeg unavailable; figures and tests still produced.')