Metrologist: MQL recipes

Recipes

GUM uncertainty budget of a pendulum measurement of g
Gauge block calibration with exact decimals
Calibration curve of a Pt100 sensor
Interlaboratory comparison: reference value, En numbers and Birge ratio
Type A and Type B uncertainty of a resistance measurement
From older CODATA values to the exact SI 2019 constants
Josephson and quantum Hall standards
Calibrating a psi gauge against a kPa reference
Lab bias on a certified reference material
Floating-point resolution of time stamps and frequency offsets

executed on MentDB

GUM uncertainty budget of a pendulum measurement of g

Measure g = 4 pi^2 L / T^2 with a length known to 0.3 mm and a period known to 0.12 ms, propagate the uncertainties, compute the effective degrees of freedom, the 95 % coverage factor from the Student law and the normalised error against standard gravity.

Pipeline

01 · unitsLength in munit convert→02 · propagateCombined uncertainty and variance sharesnum uncertainty→03 · dofWelch-Satterthwaitemath calc→04 · coverageStudent quantile with beta_incnum solve→05 · referenceStandard gravity g0const get

Script

admin

#1. Pendulum: L = 1000.4 mm, u 0.3 mm type B 50 dof, T = 2.00658 s, u 0.00012 s from 10 readings;
-> "[L]" (unit convert 1000.4 "mm" "m");
-> "[uL]" (unit convert 0.3 "mm" "m");
-> "[r]" (num uncertainty "4*pi^2*L/T^2" (concat "{\"L\": [" [L] ", " [uL] "], \"T\": [2.00658, 0.00012]}"));
json load "rx_g" [r];
#2. Effective degrees of freedom, Welch-Satterthwaite, from the variance shares;
-> "[nu]" (math calc (concat "1/(" (json select "rx_g" /variance_share/L) "^2/50+" (json select "rx_g" /variance_share/T) "^2/9)"));
#3. Coverage factor for 95 percent: Student law written with the incomplete beta function;
-> "[k]" (num solve (concat "0.5*beta_inc(" [nu] "/2, 0.5, " [nu] "/(" [nu] "+t^2))-0.025") "t" 1 10);
-> "[U]" (math calc (concat [k] "*" (json select "rx_g" /uncertainty)));
#4. Normalised error En against the standard gravity g0, an exact constant;
-> "[g0]" (const get "g0");
-> "[exact]" (json load select "rx_g0" (const info "g0") /exact);
concat "g = " (json select "rx_g" /text) " m/s2, share L " (math round (json select "rx_g" /variance_share/L) 3) ", nu_eff = " (math round [nu] 1) ", k95 = " (math round [k] 3) ", U = " (math round [U] 4) " m/s2, En = " (math round (math calc (concat "abs(" (json select "rx_g" /value) "-" [g0] ")/" [U])) 2) " vs g0 = " [g0] " (exact: " [exact] ")"

Result

mentdb

g = 9.8089 ± 0.0032 m/s2, share L 0.863, nu_eff = 58.9, k95 = 2.001, U = 0.0063 m/s2, En = 0.36 vs g0 = 9.80665 (exact: true)

Go further

Add the amplitude correction T0 = T / (1 + theta^2/16) to the model
Read the uncertainty of G or h with const info

executed on MentDB

Gauge block calibration with exact decimals

Six comparator readings of a 25 mm gauge block are averaged exactly, corrected for the reference block and for the thermal expansion (an exact fraction), rounded half-even to the 0.001 um resolution and compared with the grade 0 tolerance.

Pipeline

01 · meanExact meanbigdec mean→02 · thermalL alpha dT as a fractionarith frac_mul→03 · correctionsExact sum of correctionsbigdec sum→04 · roundHALF_EVEN to 0.001 umbigdec round→05 · toleranceCompare with 0.14 umbigdec compare

Script

admin

#1. Six comparator readings of a 25 mm gauge block against the reference, in um;
-> "[dev]" "[0.120, 0.150, 0.110, 0.140, 0.130, 0.125]";
-> "[m]" (bigdec mean [dev] 20);
#2. Thermal correction 25000 um x 1.15e-5 /K x 0.15 K as an exact fraction;
-> "[f]" (arith frac_mul (arith frac_mul 25000 0.0000115) 0.15);
-> "[th]" (bigdec div (atom get [f] 1 "/") (atom get [f] 2 "/") 20);
#3. Deviation = mean + reference block 0.045 um - thermal;
-> "[d]" (bigdec sum (concat "[" [m] ", 0.045, " (bigdec negate [th]) "]"));
#4. Round to the resolution 0.001 um, HALF_EVEN, and check grade 0 tolerance 0.14 um;
-> "[dr]" (bigdec round [d] 3 HALF_EVEN);
-> "[ok]" (bigdec compare (bigdec abs [dr]) 0.14);
concat "mean = " (bigdec precision [m] 6) " um, thermal = " [f] " = " (bigdec strip [th]) " um, deviation = " [dr] " um = " (unit prefix [dr] u n) " nm, |dev| vs 0.14 um: " [ok] ", calibrated length = " (bigdec sum (concat "[25, " (unit prefix [dr] u m) "]")) " mm"

Result

mentdb

mean = 0.129167 um, thermal = 69/1600 = 0.043125 um, deviation = 0.131 um = 131 nm, |dev| vs 0.14 um: -1, calibrated length = 25.000131 mm

Go further

Convert the deviation to inches with unit convert
Keep 50 digits in a long chain with bigdec div and bigdec precision

executed on MentDB

Calibration curve of a Pt100 sensor

Fit the Callendar-Van Dusen equation R = R0 (1 + A T + B T^2) to six calibration points by least squares on a scaled design matrix, then draw the points with the fitted curve and the pchip interpolation.

Pipeline

01 · least squaresOverdetermined system (QR)matrix solve→02 · coefficientsR0, A, Bjson select→03 · fitPoints and fitted modelplot fit→04 · interpolationOverlay the pchip curveplot interpolate

Script

admin

#1. Pt100 calibration points: T in degC and R in ohm;
-> "[T]" "[0, 50, 100, 150, 200, 250]";
-> "[R]" "[100.004, 119.394, 138.509, 157.322, 175.859, 194.096]";
#2. Least squares on 1, t, t^2 with t = T/100 to keep the matrix well conditioned;
-> "[X]" "[[1, 0, 0], [1, 0.5, 0.25], [1, 1, 1], [1, 1.5, 2.25], [1, 2, 4], [1, 2.5, 6.25]]";
-> "[c]" (matrix solve [X] [R]);
json load "rx_cal" [c];
#3. Callendar-Van Dusen R = R0 1 + A T + B T^2;
-> "[R0]" (json select "rx_cal" /[0]);
-> "[A]" (math calc (concat (json select "rx_cal" /[1]) "/(100*" [R0] ")"));
-> "[B]" (math calc (concat (json select "rx_cal" /[2]) "/(10000*" [R0] ")"));
#4. Draw the points with the fitted model, then the pchip interpolation over them;
plot fit [T] [R] (concat [R0] "*(1+" [A] "*x+" [B] "*x^2)") "x";
plot add (plot interpolate [T] [R] "pchip")

Result

mentdb

The chart is shown in the Plot tab of the editor.

Go further

Invert the curve for a reading of 150 ohm with num poly_roots
Check the conditioning of the design matrix with matrix cond and compare with matrix pinv

executed on MentDB

Interlaboratory comparison: reference value, En numbers and Birge ratio

Five laboratories calibrate the same 10 kOhm travelling standard. The reference value is their mean weighted by 1/u^2; the chi2 and the Birge ratio test the consistency, and the En number of each lab (k = 2, corrected for its correlation with the reference) flags the outliers.

Pipeline

01 · weightsw = 1/u^2 built in a JSON arrayjson iarray→02 · referenceWeighted mean and 1/sqrt(sum w)stat weighted_mean→03 · consistencychi2, Birge ratio, p-valuestat dist_cdf→04 · En(xi - xr)/2.sqrt(ui^2 - ur^2)math calc

Script

admin

#1. Interlaboratory comparison of a 10 kOhm standard: deviations from nominal in ppm and standard uncertainties;
-> "[x]" "[1.8, 2.6, 1.1, 3.9, 2.2]"; -> "[u]" "[0.4, 0.6, 0.3, 0.9, 0.5]";
json load "rx_tu1_x" [x]; json load "rx_tu1_u" [u];
#2. Weights 1/u^2, reference value = weighted mean, its uncertainty 1/sqrt of the sum of the weights;
-> "[w]" "[]"; json load "rx_tu1_w" [w];
for (-> "[i]" 0) (< [i] 5) (++ "[i]") {
	json iarray "rx_tu1_w" / (math calc (concat "1/" (json select "rx_tu1_u" (concat "/[" [i] "]")) "^2")) NUM;
};
-> "[w]" (json doc "rx_tu1_w");
-> "[xr]" (stat weighted_mean [x] [w]);
-> "[ur]" (math calc (concat "1/sqrt(" (stat sum [w]) ")"));
#3. Consistency chi2 and Birge ratio, then En number of each lab with k = 2;
-> "[chi2]" 0; -> "[en]" "";
for (-> "[i]" 0) (< [i] 5) (++ "[i]") {
	-> "[xi]" (json select "rx_tu1_x" (concat "/[" [i] "]"));
	-> "[ui]" (json select "rx_tu1_u" (concat "/[" [i] "]"));
	-> "[chi2]" (math calc "[chi2]+(([xi]-[xr])/[ui])^2");
	-> "[en]" (concat [en] " " (math round (math calc "([xi]-[xr])/(2*sqrt([ui]^2-[ur]^2))") 2));
};
concat "reference = " (math round [xr] 3) " ppm, u = " (math round [ur] 3) " ppm; chi2 = " (math round [chi2] 2) ", Birge ratio = " (math round (math sqrt (/ [chi2] 4)) 3) ", p = " (math round (- 1 (stat dist_cdf "chi2" "[4]" [chi2])) 3) "; En =" [en]

Result

mentdb

reference = 1.746 ppm, u = 0.199 ppm; chi2 = 13.23, Birge ratio = 1.819, p = 0.01; En = 0.08 0.75 -1.44 1.23 0.49

Go further

Remove the inconsistent labs (|En| > 1) and compute the reference again
Use the median and a robust scale (stat median, stat quartiles) when the comparison has many labs

executed on MentDB

Type A and Type B uncertainty of a resistance measurement

Ten readings of a 100 Ohm resistor with a 6.5-digit multimeter: normality check, Type A uncertainty of the mean, Type B from the calibration certificate and the resolution, combined uncertainty, Welch-Satterthwaite degrees of freedom and expanded uncertainty.

Pipeline

01 · normalityKurtosis and Kolmogorov-Smirnovstat ks_test→02 · type As/sqrt(n) from the sample variancestat variance→03 · type BCertificate k = 2, resolution rectangularmath calc→04 · coveragenu_eff and Student k at 95 %stat dist_inv

Script

admin

#1. Ten repeated readings of a 100 Ohm resistor with a 6.5-digit DMM, in Ohm;
-> "[r]" "[100.0213, 100.0208, 100.0219, 100.0211, 100.0215, 100.0209, 100.0217, 100.0212, 100.0214, 100.0210]";
#2. Normality check before the Type A evaluation: kurtosis, Kolmogorov-Smirnov against the fitted normal law;
-> "[m]" (stat mean [r]); -> "[s]" (math sqrt (stat variance [r]));
-> "[ks]" (json load select "rx_tu1_ks" (stat ks_test [r] "normal" (concat "[" [m] ", " [s] "]")) "/p_value");
#3. Type A: standard deviation of the mean. Type B: certificate 8 ppm k = 2 and resolution 0.1 mOhm rectangular;
-> "[uA]" (math calc "[s]/sqrt(10)");
-> "[uB1]" (math calc "8e-6*100/2");
-> "[uB2]" (math calc "0.0001/2/sqrt(3)");
-> "[uc]" (math calc "sqrt([uA]^2+[uB1]^2+[uB2]^2)");
#4. Welch-Satterthwaite effective degrees of freedom, Type B taken as infinite, and coverage factor at 95 %;
-> "[nu]" (math calc "[uc]^4/([uA]^4/9)");
-> "[k]" (stat dist_inv "student" (concat "[" (math floor [nu]) "]") 0.975);
concat "R = " (math round [m] 5) " Ohm, median " (stat median [r]) ", excess kurtosis " (math round (stat kurtosis [r]) 3) ", KS p = " (math round [ks] 3) "; uA = " (math round_sig [uA] 2) ", uB = " (math round_sig [uB1] 2) " and " (math round_sig [uB2] 2) ", uc = " (math round_sig [uc] 3) " Ohm, nu_eff = " (math round [nu] 0) ", k = " (math round [k] 3) ", U = " (math round_sig (* [k] [uc]) 2) " Ohm"

Result

mentdb

R = 100.02128 Ohm, median 100.02125, excess kurtosis -0.612, KS p = 1; uA = 0.00011, uB = 0.00040 and 0.000029, uc = 0.000416 Ohm, nu_eff = 1757, k = 1.961, U = 0.00082 Ohm

Go further

Propagate the result through a formula with num uncertainty
Show the readings with plot histogram and look for a drift with stat correlation against the reading number

executed on MentDB

From older CODATA values to the exact SI 2019 constants

Since 2019, h, e, k and N_A are exact. The script reads the older values (with their uncertainties) from the constant family, the exact values from const get, and gives the shift and the old relative uncertainty in ppb.

Pipeline

01 · oldOlder CODATA value with uncertaintyconstant physics_planck→02 · evalSame for e, k, N_Aeval→03 · exactExact SI valueconst get→04 · compareShift and old u_r in ppbmath round_sig

Script

admin

#1. The constant family keeps the older CODATA values, const get gives the exact SI 2019 values;
-> "[out]" "";
-> "[old]" "physics_planck,physics_e,physics_k,physics_n";
-> "[new]" "h,e,k,NA";
for (-> "[i]" 1) (<= [i] 4) (++ "[i]") {
	-> "[txt]" (eval (concat "constant " (atom get [old] [i] ",")));
	-> "[v]" (string replace (atom get [txt] 1 " ") "(" "");
	-> "[u]" (string replace (atom get [txt] 3 " ") ")" "");
	-> "[ex]" (const get (atom get [new] [i] ","));
	#2. Shift of the value and old relative standard uncertainty, in parts per billion;
	-> "[out]" (concat [out] (atom get [new] [i] ",") ": " [v] " -> " [ex] ", shift " (math round (math calc (concat "(" [ex] "-" [v] ")/" [ex] "*1e9")) 1) " ppb, old u_r " (math round_sig (math calc (concat [u] "/" [v] "*1e9")) 2) " ppb; ");
};
[out]

Result

mentdb

h: 6.626069E-34 -> 6.62607015E-34, shift 173.6 ppb, old u_r 170 ppb; e: 1.6021765E-19 -> 1.602176634E-19, shift 83.6 ppb, old u_r 87 ppb; k: 1.380651E-23 -> 1.380649E-23, shift -1448.6 ppb, old u_r 1700 ppb; NA: 6.022142E23 -> 6.02214076E23, shift -205.9 ppb, old u_r 170 ppb;

Go further

List all the exact constants with const list and a filter on "exact"
Propagate the old uncertainty of k into a gas thermometer with num uncertainty

executed on MentDB

Josephson and quantum Hall standards

A Josephson array gives 10 V at 70 GHz and a quantum Hall bar gives RK/2. The script counts the junctions, computes RK from the conductance quantum and from Z0/(2 alpha), and the ratio of a 10 kOhm transfer standard.

Pipeline

01 · fluxMagnetic flux quantumconstant physics_magn_flux_quantum→02 · hallRK = 2/G0constant physics_g0→03 · impedanceRK = Z0 / 2 alphaconstant physics_z0→04 · alphaFine structure constantconstant physics_fine_structure→05 · exactKJ and RK of the SIconst get

Script

admin

#1. Josephson voltage standard: 10 V from an array driven at 70 GHz, V = n f / KJ, KJ = 1 / Phi0;
-> "[phi0]" (string replace (atom get (constant physics_magn_flux_quantum) 1 " ") "(" "");
-> "[KJ]" (const get "KJ");
-> "[n]" (math calc "10*[KJ]/70e9");
#2. Quantum Hall resistance RK = h/e^2 = 2/G0 and the link with the vacuum impedance RK = Z0 / 2 alpha;
-> "[G0]" (string replace (atom get (constant physics_g0) 1 " ") "(" "");
-> "[Z0]" (string replace (atom get (constant physics_z0) 1 " ") "(" "");
-> "[alpha]" (string replace (atom get (constant physics_fine_structure) 1 " ") "(" "");
-> "[RK1]" (math calc "2/[G0]");
-> "[RK2]" (math calc "[Z0]/(2*[alpha])");
#3. A 10 kOhm transfer standard against the plateau i = 2: R/RK2 ratio;
-> "[ratio]" (math calc (concat "10000/(" (const get "RK") "/2)"));
concat "n = " (math ceil [n]) " junctions for 10 V (KJ = " [KJ] " Hz/V, 1/Phi0 from the constant = " (math round_sig (math calc "1/[phi0]") 8) "); RK = " (math round [RK1] 3) " ohm from G0, " (math round [RK2] 3) " ohm from Z0/2alpha, exact " (math round (const get "RK") 3) " ohm; 10 kOhm / RK/2 = " (math round [ratio] 6)

Result

mentdb

n = 69086 junctions for 10 V (KJ = 4.835978484169836E14 Hz/V, 1/Phi0 from the constant = 483597850000000); RK = 25812.807 ohm from G0, 25812.807 ohm from Z0/2alpha, exact 25812.807 ohm; 10 kOhm / RK/2 = 0.774809

Go further

Compute the Josephson step voltage n f / KJ for n = 1 at 70 GHz in microvolts
Compare the old conventional values KJ-90 = 483597.9 GHz/V and RK-90 = 25812.807 ohm with the SI ones

executed on MentDB

Calibrating a psi gauge against a kPa reference

A pressure gauge graduated in psi is checked against a reference in kPa. The script reads the conversion factors from the unit catalogue, converts the readings, computes the errors and the maximal error in % of full scale.

Pipeline

01 · catalogueUnits of pressure with their SI factorunit list→02 · countNumber of unitsjson count→03 · convertReadings in kPavector scale→04 · errorsGauge - referencevector sub→05 · full scale100 psi in kPaunit convert

Script

admin

#1. Pressure gauge graduated in psi, reference in kPa: factors from the unit catalogue;
json load "rx_p4b_units" (unit list "pressure");
-> "[k_psi]" ""; -> "[k_kpa]" "";
for (-> "[i]" 0) (< [i] (json count "rx_p4b_units" /units)) (++ "[i]") {
	-> "[u]" (json select "rx_p4b_units" (concat "/units[" [i] "]/unit"));
	if (equal [u] "psi") { -> "[k_psi]" (json select "rx_p4b_units" (concat "/units[" [i] "]/to_si")); };
	if (equal [u] "kPa") { -> "[k_kpa]" (json select "rx_p4b_units" (concat "/units[" [i] "]/to_si")); };
};
-> "[nu]" (json count "rx_p4b_units" /units);
json unload "rx_p4b_units";
#2. Calibration points: gauge reading in psi, reference in kPa;
-> "[g]" "[0, 25, 50, 75, 100]";
-> "[r]" "[0, 172.9, 345.1, 516.4, 688.2]";
#3. Error of the gauge in kPa at each point and the maximal error in % of full scale;
-> "[gk]" (vector scale [g] (math calc "[k_psi]/[k_kpa]"));
-> "[err]" (vector sub [gk] [r]);
-> "[fs]" (unit convert 100 "psi" "kPa");
concat [nu] " pressure units known, 1 psi = " (math round (math calc "[k_psi]/[k_kpa]") 6) " kPa; errors (kPa) " [err] "; max error " (math round (math calc (concat "max(abs(" (stat max [err]) "), abs(" (stat min [err]) "))/" [fs] "*100")) 3) " % of full scale " (math round [fs] 2) " kPa"

Result

mentdb

11 pressure units known, 1 psi = 6.894757 kPa; errors (kPa) [0,-0.5310676707910034,-0.3621353415820181,0.7067969876270581,1.2757293168359638]; max error 0.185 % of full scale 689.48 kPa

Go further

Fit the error curve with stat fit and correct the gauge
Show all the quantities known by unit list without argument

executed on MentDB

Lab bias on a certified reference material

A lab measures lead in a soil CRM certified at 50.0 mg/kg with a known reproducibility of 0.8 mg/kg. The script runs a z-test, rebuilds the p-value with erfc and gives the coverage of k = 2 and the k for 99 %.

Pipeline

01 · testz-test with known sigmastat ztest→02 · p-valuep = erfc(|z|/sqrt 2)math erfc→03 · coverageerf(k/sqrt 2)math erf→04 · k99k = sqrt 2 erf_inv(0.99)math erf_inv→05 · spreadPopulation variancestat pvariance

Script

admin

#1. Lead in a soil CRM certified at 50.0 mg/kg, reproducibility sigma 0.8 mg/kg, 7 results of the lab;
-> "[x]" "[50.9, 51.3, 50.6, 51.8, 50.2, 51.1, 50.8]";
#2. z-test of the mean against the certified value with the known sigma;
json load "rx_p4b_z" (stat ztest [x] 50.0 0.8);
-> "[z]" (json select "rx_p4b_z" /z);
-> "[p]" (json select "rx_p4b_z" /p_value);
json unload "rx_p4b_z";
#3. Same p-value from the complementary error function: p = erfc |z| / sqrt 2;
-> "[p2]" (math erfc (math calc "abs([z])/sqrt(2)"));
#4. Coverage probability of k = 2, and the k for 99 %;
-> "[cov]" (math erf (math calc "2/sqrt(2)"));
-> "[k99]" (math calc (concat "sqrt(2)*" (math erf_inv 0.99)));
concat "bias " (math round (math calc (concat (stat mean [x]) "-50")) 3) " mg/kg, z = " (math round [z] 3) ", p = " (math round_sig [p] 3) " = " (math round_sig [p2] 3) " from erfc: significant bias; spread of the lab " (math round (math sqrt (stat pvariance [x])) 3) " mg/kg; k = 2 covers " (math round (math calc "[cov]*100") 2) " %, 99 % needs k = " (math round [k99] 3)

Result

mentdb

bias 0.957 mg/kg, z = 3.165, p = 0.00155 = 0.00155 from erfc: significant bias; spread of the lab 0.475 mg/kg; k = 2 covers 95.45 %, 99 % needs k = 2.576

Go further

Compute the En number with the uncertainty of the certificate
Use stat ttest_1 when sigma is not known

executed on MentDB

Floating-point resolution of time stamps and frequency offsets

A data logger stores time stamps as doubles in seconds since 1970, and a quartz has a fractional frequency offset of 1e-13. The script shows the ulp at that magnitude, what is left of the fraction of second, and why log1p and expm1 are needed for tiny offsets.

Pipeline

01 · ulpSpacing of doubles at tmath ulp→02 · fractionDecimal fraction of secondmath frac_part→03 · naiveln(1 + y) computed directlymath calc→04 · accuratelog1p and expm1math log1p

Script

admin

#1. Time stamps of a data logger stored as doubles in seconds since 1970;
-> "[t]" 1761955200.123456789;
#2. Smallest step a double can hold at this magnitude, and the fraction of second kept;
-> "[ulp]" (math ulp [t]);
-> "[frac]" (math frac_part [t]);
-> "[kept]" (math calc "[t]-1761955200");
#3. Fractional frequency offset of a quartz y = 1e-13: ln 1+y naive against log1p, and e^y - 1;
-> "[naive]" (math calc "log(1+1e-13)");
-> "[good]" (math log1p 1e-13);
-> "[em1]" (math expm1 1e-13);
concat "ulp at t = " [ulp] " s: the fraction " [frac] " s is stored as " [kept] " s, nanoseconds are lost; ln(1+y): naive " [naive] " (error " (math round (math calc (concat "abs(" [naive] "-1e-13)/1e-13*100")) 3) " %), log1p " [good] ", expm1 " [em1] "; ulp of 1.0 = " (math ulp 1)

Result

mentdb

ulp at t = 2.384185791015625E-7 s: the fraction 0.123456789 s is stored as 0.12345671653747559 s, nanoseconds are lost; ln(1+y): naive 4.3394735994895746E-14 (error 56.605 %), log1p 9.9999999999995E-14, expm1 1.00000000000005E-13; ulp of 1.0 = 2.220446049250313E-16

Go further