Physicist: MQL recipes

Recipes

Range of a baseball with air drag
Period of a pendulum at large amplitude
Logarithmic decrement and Q of a damped oscillator
Quantum scattering by a hard sphere
Total internal reflection and the evanescent wave
Relativistic velocity addition with rapidities
Skydiver with quadratic drag: tanh and acosh
Do the stored constants agree with each other?
Fraction of fast molecules in a Maxwell gas
Is our pendulum g compatible with standard gravity?

executed on MentDB

Range of a baseball with air drag

A baseball leaves the bat at 40 m/s and 35 degrees. Integrate the equations of motion with quadratic drag (k = rho*Cd*A/2m = 0.0062 1/m) until it lands, and compare with the range in vacuum.

Pipeline

01 · launchVelocity componentsunit calc→02 · integrateEquations of motion, stop at y = 0num ode→03 · readLanding pointjson select→04 · compareRange in vacuumunit calc

Script

admin

#1. Baseball hit at 40 m/s and 35 degrees, drag factor k = rho Cd A / 2m in 1/m;
-> "[vx0]" (unit calc "40 m/s * cos(35 deg)" "m/s");
-> "[vy0]" (unit calc "40 m/s * sin(35 deg)" "m/s");
-> "[k]" 0.0062;
#2. Trajectory with quadratic drag, stop when y crosses 0 going down;
-> "[traj]" (num ode "{\"x\": \"vx\", \"y\": \"vy\", \"vx\": \"-[k]*sqrt(vx^2+vy^2)*vx\", \"vy\": \"-_g0-[k]*sqrt(vx^2+vy^2)*vy\"}" (concat "{\"x\": 0, \"y\": 0, \"vx\": " [vx0] ", \"vy\": " [vy0] "}") 0 10 "{\"points\": 11, \"event\": \"y\", \"direction\": -1}");
-> "[range]" (json load select "rx_proj" [traj] /event/x);
-> "[tf]" (json select "rx_proj" /event/t);
#3. Compare with the range in vacuum v0^2 sin 2a / g;
concat "range with drag = " (math round [range] 1) " m, flight time " (math round [tf] 2) " s; in vacuum = " (math round (unit calc "(40 m/s)^2*sin(70 deg)/_g0" "m") 1) " m"

Result

mentdb

range with drag = 94.8 m, flight time 4.06 s; in vacuum = 153.3 m

Go further

Draw the trajectory with plot ode
Find the launch angle that maximises the range with drag using num maximize around the ODE

executed on MentDB

Period of a pendulum at large amplitude

A 2 m pendulum is released at 60 degrees. Compare the small-angle period 2*pi*sqrt(L/g) with the exact period, which uses the complete elliptic integral K(k), k = sin(theta0/2).

Pipeline

01 · kModulus k = sin(theta0/2)math eval→02 · T0Small-angle periodunit calc→03 · KElliptic integralnum integrate→04 · compareRelative errormath calc

Script

admin

#1. Pendulum L = 2 m released at 60 degrees;
-> "[L]" 2; -> "[theta0]" 60;
-> "[k]" (math eval "sin(t/2)" (concat "{\"t\": " [theta0] ", \"_angle\": \"deg\"}"));
#2. Small-angle period 2 pi sqrt L/g;
-> "[T0]" (unit calc (concat "2*pi*sqrt(" [L] " m/_g0)") "s");
#3. Exact period T = T0 * 2K/pi, the integral from 0 to pi gives 2K;
-> "[K2]" (num integrate (concat "1/sqrt(1-" [k] "^2*sin(p)^2)") "p" 0 (math pi));
-> "[T]" (math calc "[T0]*[K2]/pi");
concat "k = " [k] ", T0 = " (math round [T0] 4) " s, exact T = " (math round [T] 4) " s: the small-angle formula is " (math round (math calc "([T]-[T0])/[T]*100") 2) " % too short"

Result

mentdb

k = 0.5, T0 = 2.8375 s, exact T = 3.0451 s: the small-angle formula is 6.82 % too short

Go further

Check by integrating theta'' = -g/L sin(theta) with num ode and an event on the angular velocity
Tabulate T/T0 for 10 to 90 degrees with num table

executed on MentDB

Logarithmic decrement and Q of a damped oscillator

Simulate a mass-spring-damper (m = 0.5 kg, k = 200 N/m, b = 0.4 kg/s) released from 5 cm, detect the successive maxima, and get the damped period, the logarithmic decrement and the quality factor, compared with theory.

Pipeline

01 · simulateFree response x(t)num ode→02 · peaksSuccessive maximadsp peaks→03 · decrementdelta = ln(A1/A2), Q = pi/deltamath calc

Script

admin

#1. Mass-spring-damper: m = 0.5 kg, k = 200 N/m, b = 0.4 kg/s;
-> "[m]" 0.5; -> "[k]" 200; -> "[b]" 0.4;
#2. Free response from x = 5 cm, 2001 points over 2 s;
-> "[sol]" (json load select "rx_osc" (num ode "{\"x\": \"v\", \"v\": \"-([b]*v+[k]*x)/[m]\"}" "{\"x\": 0.05, \"v\": 0}" 0 2 "{\"points\": 2001}") /x);
#3. Successive maxima;
json load "rx_pk" (dsp peaks [sol]);
-> "[A1]" (json select "rx_pk" /values[0]); -> "[A2]" (json select "rx_pk" /values[1]);
-> "[Td]" (math calc (concat "(" (json select "rx_pk" /indices[1]) "-" (json select "rx_pk" /indices[0]) ")*0.001"));
#4. Logarithmic decrement and quality factor;
-> "[delta]" (math calc "ln([A1]/[A2])");
concat "Td = " [Td] " s, delta = " (math round [delta] 4) " (theory b*Td/2m = " (math round (math calc "[b]*[Td]/(2*[m])") 4) "), Q = pi/delta = " (math round (math calc "pi/[delta]") 1) " (theory sqrt(k*m)/b = " (math round (math calc "sqrt([k]*[m])/[b]") 1) ")"

Result

mentdb

Td = 0.314 s, delta = 0.1257 (theory b*Td/2m = 0.1256), Q = pi/delta = 25 (theory sqrt(k*m)/b = 25)

Go further

Fit the envelope A*exp(-t/tau) on the peaks with stat fit_model
Add a sinusoidal force and find the resonance with num maximize

executed on MentDB

Quantum scattering by a hard sphere

For a hard sphere of radius a, the phase shift of the l-th partial wave is tan(delta_l) = j_l(ka)/y_l(ka). The script checks delta_0 = -ka, sums the partial waves into the total cross section at ka = 1.5, and recovers the low-energy limit 4 pi a^2.

Pipeline

01 · phasedelta_0 from j_0 and y_0math calc→02 · sumPartial wave sum to l = 30num sum→03 · limitLow energy checksph_bessel_y(n, x)

Script

admin

#1. Hard sphere scattering: radius a, wave number k, ka = 1.5;
-> "[ka]" 1.5;
#2. Phase shift of the l wave: tan delta_l = j_l of ka / y_l of ka, here l = 0;
-> "[d0]" (math calc "atan(sph_bessel_j(0, [ka])/sph_bessel_y(0, [ka]))");
#3. Partial wave sum: sigma / pi a^2 = 4/ka^2 * sum of 2l+1 sin^2 delta_l;
-> "[sig]" (math calc (concat "4/[ka]^2*" (num sum "(2*l+1)*sin(atan(sph_bessel_j(l, 1.5)/sph_bessel_y(l, 1.5)))^2" "l" 0 30)));
#4. Low energy limit ka = 0.01: four times the geometric cross section;
-> "[sig0]" (math calc (concat "4/0.01^2*" (num sum "(2*l+1)*sin(atan(sph_bessel_j(l, 0.01)/sph_bessel_y(l, 0.01)))^2" "l" 0 10)));
concat "{\"delta0_rad\": " [d0] ", \"minus_ka\": " (math calc "-[ka]") ", \"sigma_over_pi_a2_ka1.5\": " [sig] ", \"sigma_over_pi_a2_ka0.01\": " [sig0] "}"

Result

mentdb

{"delta0_rad": -1.5, "minus_ka": -1.5, "sigma_over_pi_a2_ka1.5": 3.152393933314819, "sigma_over_pi_a2_ka0.01": 3.9998666817761683}

Go further

Tabulate sigma against ka with num table and draw it with plot xy: it tends to 2 pi a^2 at high energy.
The spherical Bessel functions sph_bessel_j and sph_bessel_y also give the modes of a spherical cavity.

executed on MentDB

Total internal reflection and the evanescent wave

A HeNe beam hits a glass-air interface at 60 degrees, beyond the critical angle. Snell law gives a complex refraction angle; the Fresnel coefficient has modulus 1 and a phase shift, and the evanescent intensity decays over a few tens of nanometres.

Pipeline

01 · criticalCritical angle asin n2/n1math asin→02 · snellComplex angle tt = asin of sin > 1complex asin→03 · fresnelrs = (n1 cos ti - n2 cos tt)/(sum)complex div→04 · evanescentcos tt = i.kappa, decay lengthcomplex sqrt

Script

admin

#1. Glass n1 = 1.5 to air n2 = 1, incidence 60 degrees, HeNe laser 633 nm;
-> "[n1]" 1.5; -> "[n2]" 1; -> "[ti]" (math deg_to_rad 60); -> "[lambda]" 633;
-> "[tc]" (math rad_to_deg (math asin (/ [n2] [n1])));
#2. Snell law beyond the critical angle: sin tt > 1, the refraction angle is complex;
-> "[s]" (math calc "[n1]/[n2]*sin([ti])");
-> "[tt]" (complex asin [s]);
-> "[ct]" (complex cos [tt]);
#3. Fresnel coefficient rs = n1 cos ti - n2 cos tt over n1 cos ti + n2 cos tt;
-> "[a]" (math calc "[n1]*cos([ti])");
-> "[rs]" (complex div (complex sub [a] (complex mul [n2] [ct])) (complex add [a] (complex mul [n2] [ct])));
#4. Evanescent field: cos tt = i.kappa, penetration depth of the intensity;
-> "[kappa]" (complex im (complex sqrt (math calc "1-[s]^2")));
-> "[dp]" (math calc "[lambda]/(4*pi*[n2]*[kappa])");
concat "critical angle = " (math round [tc] 2) " deg, sin tt = " (math round [s] 4) ", tt = " [tt] ", |rs| = " (math round (complex abs [rs]) 6) ", phase shift of the s wave = " (math round (math rad_to_deg (complex arg [rs])) 2) " deg, intensity decay length = " (math round [dp] 1) " nm"

Result

mentdb

critical angle = 41.81 deg, sin tt = 1.299, tt = 1.5707963267948966-0.7552738752366037i, |rs| = 1, phase shift of the s wave = -95.74 deg, intensity decay length = 60.8 nm

Go further

Do the p polarization: rp with n2 cos ti and n1 cos tt, and the phase difference used in a Fresnel rhomb
Scan the incidence angle with a for loop and draw the phase with plot xy

executed on MentDB

Relativistic velocity addition with rapidities

A ship moves at 0.8c and fires a probe at 0.9c in its own frame. Rapidities add (phi = atanh beta), so the probe speed is tanh(phi1 + phi2), its Lorentz factor cosh(phi) and its beta-gamma sinh(phi).

Pipeline

01 · cSpeed of lightconstant physics_c→02 · rapidityphi = atanh betamath atanh→03 · addbeta = tanh of the summath tanh→04 · gammagamma = cosh, beta gamma = sinhmath cosh

Script

admin

-> "[c]" (string regexp_extract (constant physics_c) "[-0-9.]+E?[-0-9]*" 0);
#1. Rapidities of the ship and of the probe in the ship frame;
-> "[p1]" (math atanh 0.8);
-> "[p2]" (math atanh 0.9);
#2. Rapidities add;
-> "[p]" (math calc "[p1]+[p2]");
-> "[b]" (math tanh [p]);
#3. Check with the textbook formula u+v over 1+uv;
concat "phi = " (math round [p1] 4) " + " (math round [p2] 4) " = " (math round [p] 4) "; beta = " [b] " = " (math calc "(0.8+0.9)/(1+0.8*0.9)") ", speed = " (math round_sig (math calc "[b]*[c]") 6) " m/s; gamma = cosh phi = " (math round (math cosh [p]) 4) ", beta gamma = sinh phi = " (math round (math sinh [p]) 4)

Result

mentdb

phi = 1.0986 + 1.4722 = 2.5708; beta = 0.9883720930232558 = 0.9883720930232558, speed = 296306000 m/s; gamma = cosh phi = 6.5766, beta gamma = sinh phi = 6.5001

Go further

A constant proper acceleration a gives phi = a tau / c: the rapidity grows linearly with the traveller's time.
Compare with unit calc for the kinetic energy (gamma - 1) m c^2.

executed on MentDB

Skydiver with quadratic drag: tanh and acosh

An 80 kg skydiver (Cd.A = 0.7 m^2, air 1.2 kg/m^3) falls from rest. With quadratic drag v(t) = vt tanh(g t/vt) and y(t) = vt^2/g ln cosh(g t/vt): find the terminal speed, the time to fall 1000 m and the speed then.

Pipeline

01 · gStandard gravityconstant physics_g→02 · vtTerminal speedmath sqrt→03 · timeTime for 1000 m by acoshmath acosh→04 · speedSpeed at that timemath tanh

Script

admin

-> "[g]" (string regexp_extract (constant physics_g) "[-0-9.]+E?[-0-9]*" 0);
-> "[m]" 80;
-> "[k]" (math calc "1.2*0.7/2");
#1. Terminal speed where drag = weight;
-> "[vt]" (math sqrt (math calc "[m]*[g]/[k]"));
#2. Invert y = vt^2/g ln cosh: t = vt/g acosh of exp of g h / vt^2;
-> "[t]" (math calc (concat "[vt]/[g]*" (math acosh (math exp (math calc "[g]*1000/[vt]^2")))));
-> "[v]" (math calc (concat "[vt]*" (math tanh (math calc "[g]*[t]/[vt]"))));
concat "vt = " (math round [vt] 2) " m/s = " (math round (unit convert [vt] "m/s" "km/h") 1) " km/h; 1000 m in " (math round [t] 1) " s, speed " (math round [v] 2) " m/s; without drag: " (math round (math sqrt (math calc "2*1000/[g]")) 1) " s"

Result

mentdb

vt = 43.22 m/s = 155.6 km/h; 1000 m in 26.2 s, speed 43.22 m/s; without drag: 14.3 s

Go further

Integrate the same fall with num ode and compare with the closed form.
With the head-down position Cd.A drops to about 0.2 m^2: recompute vt.

executed on MentDB

Do the stored constants agree with each other?

Check the defining relations R = N_A k and F = N_A e with the CODATA 2022 values of the 'constant' family, then compare them with 'const get' (k, e and N_A are exact since the 2019 SI redefinition, so R and F are exact too).

Pipeline

01 · exactExact N_A, k, econstant physics_n→02 · RR and F storedconstant physics_r→03 · faradayFaraday constantconstant physics_f→04 · codataSame values in const getconst get

Script

admin

#1. Numbers out of the constants, exact ones have no uncertainty;
-> "[n]" (string regexp_extract (constant physics_n) "[-0-9.]+E?[-0-9]*" 0);
-> "[k]" (string regexp_extract (constant physics_k) "[-0-9.]+E?[-0-9]*" 0);
-> "[e]" (string regexp_extract (constant physics_e) "[-0-9.]+E?[-0-9]*" 0);
-> "[R]" (string regexp_extract (constant physics_r) "[-0-9.]+E?[-0-9]*" 0);
-> "[F]" (string regexp_extract (constant physics_f) "[-0-9.]+E?[-0-9]*" 0);
#2. Defining relations;
-> "[dR]" (math round_sig (math calc "[n]*[k]/[R]-1") 2);
-> "[dF]" (math round_sig (math calc "[n]*[e]/[F]-1") 2);
#3. Gap with the values of const get;
concat "R = " (constant physics_r) ", N.k/R - 1 = " [dR] ", N.e/F - 1 = " [dF] "; const get: R = " (const get "R") ", relative gap of R " (math round_sig (math calc (concat "[R]/" (const get "R") "-1")) 2) ", of e " (math round_sig (math calc (concat "[e]/" (const get "e") "-1")) 2)

Result

mentdb

R = 8.31446261815324 J/(mol·K), N.k/R - 1 = 0, N.e/F - 1 = 0; const get: R = 8.31446261815324, relative gap of R 0, of e 0

Go further

const info gives the uncertainty and whether a constant is exact.
Since 2019, h, e, k and N_A are fixed: R = N_A k and F = N_A e are exact, so both gaps are exactly 0.

executed on MentDB

Fraction of fast molecules in a Maxwell gas

In nitrogen at 300 K, what fraction of the molecules move faster than 1000 m/s? With the most probable speed vp = sqrt(2kT/m) and x = v/vp, the Maxwell tail is erfc(x) + 2x exp(-x^2)/sqrt(pi).

Pipeline

01 · constantsk and the atomic mass unitconstant physics_amu→02 · vpMost probable speedmath sqrt→03 · tailComplementary error functionmath erfc→04 · checkBelow vp: erf partmath erf

Script

admin

-> "[k]" (string regexp_extract (constant physics_k) "[-0-9.]+E?[-0-9]*" 0);
-> "[m]" (math calc (concat "28.014*" (string regexp_extract (constant physics_amu) "[-0-9.]+E?[-0-9]*" 0)));
-> "[T]" 300;
#1. Most probable speed;
-> "[vp]" (math sqrt (math calc "2*[k]*[T]/[m]"));
#2. Tail above 1000 m/s;
-> "[x]" (math calc "1000/[vp]");
-> "[F]" (math calc (concat (math erfc [x]) "+2*[x]*exp(-[x]^2)/sqrt(pi)"));
#3. Fraction below vp, x = 1;
-> "[low]" (math calc (concat (math erf 1) "-2*exp(-1)/sqrt(pi)"));
concat "vp = " (math round [vp] 1) " m/s, mean speed = " (math round (math calc "2*[vp]/sqrt(pi)") 1) " m/s; faster than 1000 m/s: " (math round (math calc "100*[F]") 3) " %; slower than vp: " (math round (math calc "100*[low]") 2) " %"

Result

mentdb

vp = 422 m/s, mean speed = 476.2 m/s; faster than 1000 m/s: 1.054 %; slower than vp: 42.76 %

Go further

Escape of hydrogen from a planet: redo it with m = 2 u and v = 11.2 km/s.
Integrate the Maxwell density with num integrate to check the erfc formula.

executed on MentDB

Is our pendulum g compatible with standard gravity?

Eight groups timed a 1.000 m pendulum and got g values. The known timing scatter is 0.03 m/s^2. A z-test against standard gravity tells whether the setup has a systematic error; the mean period is deduced with T = 2 pi sqrt(L/g).

Pipeline

01 · referenceStandard gravityconstant physics_g→02 · meanMean of the groupsmath avg→03 · testz-test, sigma knownstat ztest→04 · periodPeriod from pi^2constant math_pi_square

Script

admin

-> "[g0]" (string regexp_extract (constant physics_g) "[-0-9.]+E?[-0-9]*" 0);
-> "[g]" "[9.79, 9.76, 9.82, 9.77, 9.78, 9.75, 9.80, 9.78]";
#1. Mean and z-test with known sigma;
-> "[gm]" (math avg [g]);
-> "[z]" (stat ztest [g] [g0] 0.03);
#2. Period of the 1 m pendulum from g = 4 pi^2 L / T^2;
-> "[pi2]" (string regexp_extract (constant math_pi_square) "[-0-9.]+E?[-0-9]*" 0);
concat "mean g = " [gm] " vs " [g0] "; z-test " [z] "; implied period " (math round (math sqrt (math calc "4*[pi2]*1/[gm]")) 4) " s vs " (math round (math sqrt (math calc "4*[pi2]/[g0]")) 4) " s"

Result

mentdb

mean g = 9.78125 vs 9.80665; z-test {"z": -2.3947349656183867, "p_value": 0.01663238288346709}; implied period 2.009 s vs 2.0064 s

Go further

A p-value below 0.05 points to a bias: the large-angle correction or the length to the centre of mass.
Local g also depends on latitude: about 9.78 at the equator, 9.83 at the poles.