diff --git a/fuse/abb_test.go b/fuse/abb_test.go
new file mode 100644
index 0000000..9097e23
--- /dev/null
+++ b/fuse/abb_test.go
@@ -0,0 +1,51 @@
+package fuse
+
+import (
+ "math"
+ "testing"
+)
+
+// Every gG band against every vertex of the ABB curve it was fitted to,
+// not only the eleven fitted points: a three-body ladder cannot follow the
+// curve's whole shape (see gg), so the bound is loose, 0.5 decades, and
+// only catches a regression.
+func TestGGBandsAgainstTheWholeABBCurve(t *testing.T) {
+ p := profile(t, IecGg)
+ for _, b := range p.Bands {
+ vs := abbVertices[b.RatingA]
+ if len(vs) != 48 {
+ t.Fatalf("%g A: %d ABB vertices", b.RatingA, len(vs))
+ }
+ var sq, worst, worstAt float64
+ for _, v := range vs {
+ got := b.MeltSeconds(v[0], RefAmbientC)
+ d := math.Log10(got / v[1])
+ if math.IsInf(got, 1) {
+ d = math.Inf(1)
+ }
+ sq += d * d
+ if math.Abs(d) > math.Abs(worst) {
+ worst, worstAt = d, v[0]
+ }
+ }
+ rms := math.Sqrt(sq / float64(len(vs)))
+ t.Logf("%g A: rms %.3f decades, worst %+.2f at %.2f ×", b.RatingA, rms, worst, worstAt)
+ if math.Abs(worst) > 0.5 {
+ t.Errorf("%g A: %+.2f decades from ABB at %.2f ×, want within 0.5", b.RatingA, worst, worstAt)
+ }
+ }
+}
+
+// abbVertices are the 48 vertices of each ABB OFAF 000 gG pre-arcing curve
+// (1SCC317001C0201 page 14), as multiples of In and seconds, read from the
+// PDF's vector paths.
+var abbVertices = map[float64][][2]float64{
+ 16: {{1.4324, 7990.5}, {1.4564, 5987.4}, {1.4718, 5000}, {1.4898, 3994.1}, {1.5147, 2993.7}, {1.5517, 1997}, {1.5777, 1496.8}, {1.6139, 998.51}, {1.6334, 801.12}, {1.6608, 600.3}, {1.681, 499.25}, {1.7039, 400.56}, {1.7351, 300.15}, {1.7857, 200.28}, {1.8239, 150.07}, {1.8855, 100.14}, {1.9228, 80.014}, {1.9758, 59.956}, {2.0118, 50.068}, {2.0612, 39.996}, {2.1308, 29.978}, {2.243, 20.003}, {2.3329, 14.989}, {2.478, 9.9988}, {2.5657, 7.9895}, {2.6886, 6.0112}, {2.7752, 4.9993}, {2.8909, 3.9947}, {3.0522, 3.0056}, {3.3065, 1.9973}, {3.5125, 1.5028}, {3.8396, 0.99866}, {4.0422, 0.80124}, {4.3263, 0.60039}, {4.5201, 0.49932}, {4.7722, 0.40051}, {5.1313, 0.30019}, {5.6945, 0.20031}, {6.1324, 0.15009}, {6.8258, 0.10012}, {7.251, 0.080004}, {7.8318, 0.059965}, {8.2196, 0.050076}, {8.7308, 0.040002}, {9.4583, 0.029982}, {10.593, 0.020001}, {11.511, 0.014987}, {13.008, 0.01}},
+ 20: {{1.4136, 7990.5}, {1.435, 5987.4}, {1.4481, 5000}, {1.4657, 3994.1}, {1.488, 2993.7}, {1.5221, 1997}, {1.5452, 1496.8}, {1.5807, 998.51}, {1.6022, 801.12}, {1.6316, 600.3}, {1.6513, 499.25}, {1.6766, 400.56}, {1.7124, 300.15}, {1.7675, 200.28}, {1.8107, 150.07}, {1.8833, 100.14}, {1.9264, 80.014}, {1.9914, 59.956}, {2.0341, 50.068}, {2.0932, 39.996}, {2.1771, 29.978}, {2.3092, 20.003}, {2.4198, 14.989}, {2.5899, 9.9988}, {2.6978, 7.9895}, {2.8442, 6.0112}, {2.9491, 4.9993}, {3.0812, 3.9947}, {3.2681, 3.0056}, {3.5564, 1.9973}, {3.7893, 1.5028}, {4.1486, 0.99866}, {4.367, 0.80124}, {4.6672, 0.60039}, {4.8763, 0.49932}, {5.1409, 0.40051}, {5.5192, 0.30019}, {6.1062, 0.20031}, {6.5656, 0.15009}, {7.2975, 0.10012}, {7.7521, 0.080004}, {8.3722, 0.059965}, {8.8003, 0.050076}, {9.362, 0.040002}, {10.173, 0.029982}, {11.463, 0.020001}, {12.511, 0.014987}, {14.225, 0.01}},
+ 25: {{1.397, 7990.5}, {1.4204, 5987.4}, {1.4355, 5000}, {1.4551, 3994.1}, {1.4795, 2993.7}, {1.5134, 1997}, {1.5412, 1496.8}, {1.5788, 998.51}, {1.6028, 801.12}, {1.6321, 600.3}, {1.6544, 499.25}, {1.6797, 400.56}, {1.7156, 300.15}, {1.7734, 200.28}, {1.8195, 150.07}, {1.8895, 100.14}, {1.9358, 80.014}, {1.9951, 59.956}, {2.0377, 50.068}, {2.0971, 39.996}, {2.1745, 29.978}, {2.3028, 20.003}, {2.406, 14.989}, {2.5674, 9.9988}, {2.6703, 7.9895}, {2.8109, 6.0112}, {2.9059, 4.9993}, {3.0314, 3.9947}, {3.2104, 3.0056}, {3.499, 1.9973}, {3.7223, 1.5028}, {4.0879, 0.99866}, {4.3164, 0.80124}, {4.6335, 0.60039}, {4.8485, 0.49932}, {5.1426, 0.40051}, {5.554, 0.30019}, {6.2112, 0.20031}, {6.7488, 0.15009}, {7.6035, 0.10012}, {8.1387, 0.080004}, {8.9105, 0.059965}, {9.4511, 0.050076}, {10.147, 0.040002}, {11.159, 0.029982}, {12.803, 0.020001}, {14.124, 0.014987}, {16.277, 0.01}},
+ 32: {{1.3939, 7990.5}, {1.4216, 5987.4}, {1.4367, 5000}, {1.4586, 3994.1}, {1.4851, 2993.7}, {1.5261, 1997}, {1.554, 1496.8}, {1.5969, 998.51}, {1.6211, 801.12}, {1.6558, 600.3}, {1.6784, 499.25}, {1.7065, 400.56}, {1.7456, 300.15}, {1.8074, 200.28}, {1.857, 150.07}, {1.9315, 100.14}, {1.9786, 80.014}, {2.0455, 59.956}, {2.0892, 50.068}, {2.1501, 39.996}, {2.2328, 29.978}, {2.3683, 20.003}, {2.4744, 14.989}, {2.6483, 9.9988}, {2.7544, 7.9895}, {2.9039, 6.0112}, {3.0064, 4.9993}, {3.141, 3.9947}, {3.3316, 3.0056}, {3.6363, 1.9973}, {3.8748, 1.5028}, {4.2614, 0.99866}, {4.4996, 0.80124}, {4.8376, 0.60039}, {5.0694, 0.49932}, {5.3769, 0.40051}, {5.8166, 0.30019}, {6.5236, 0.20031}, {7.1101, 0.15009}, {8.0599, 0.10012}, {8.6521, 0.080004}, {9.5155, 0.059965}, {10.124, 0.050076}, {10.934, 0.040002}, {12.116, 0.029982}, {14.092, 0.020001}, {15.734, 0.014987}, {18.493, 0.01}},
+ 35: {{1.3931, 7990.5}, {1.4208, 5987.4}, {1.4381, 5000}, {1.46, 3994.1}, {1.4889, 2993.7}, {1.5299, 1997}, {1.5602, 1496.8}, {1.6057, 998.51}, {1.6301, 801.12}, {1.6649, 600.3}, {1.6876, 499.25}, {1.7186, 400.56}, {1.758, 300.15}, {1.82, 200.28}, {1.8675, 150.07}, {1.9451, 100.14}, {1.9897, 80.014}, {2.057, 59.956}, {2.1007, 50.068}, {2.162, 39.996}, {2.2452, 29.978}, {2.3814, 20.003}, {2.488, 14.989}, {2.6629, 9.9988}, {2.7696, 7.9895}, {2.9199, 6.0112}, {3.023, 4.9993}, {3.1584, 3.9947}, {3.35, 3.0056}, {3.6568, 1.9973}, {3.8961, 1.5028}, {4.2849, 0.99866}, {4.5244, 0.80124}, {4.8643, 0.60039}, {5.0973, 0.49932}, {5.4066, 0.40051}, {5.8487, 0.30019}, {6.5697, 0.20031}, {7.1493, 0.15009}, {8.1044, 0.10012}, {8.6999, 0.080004}, {9.568, 0.059965}, {10.18, 0.050076}, {10.995, 0.040002}, {12.183, 0.029982}, {14.17, 0.020001}, {15.82, 0.014987}, {18.595, 0.01}},
+ 40: {{1.3987, 7990.5}, {1.4263, 5987.4}, {1.4437, 5000}, {1.4634, 3994.1}, {1.4925, 2993.7}, {1.5335, 1997}, {1.564, 1496.8}, {1.6071, 998.51}, {1.6315, 801.12}, {1.6664, 600.3}, {1.6891, 499.25}, {1.7201, 400.56}, {1.7595, 300.15}, {1.8244, 200.28}, {1.8747, 150.07}, {1.9527, 100.14}, {2.0005, 80.014}, {2.0681, 59.956}, {2.1154, 50.068}, {2.177, 39.996}, {2.264, 29.978}, {2.4051, 20.003}, {2.5203, 14.989}, {2.7016, 9.9988}, {2.8142, 7.9895}, {2.9715, 6.0112}, {3.0764, 4.9993}, {3.2191, 3.9947}, {3.4193, 3.0056}, {3.7382, 1.9973}, {3.9948, 1.5028}, {4.3934, 0.99866}, {4.6461, 0.80124}, {4.9951, 0.60039}, {5.2425, 0.49932}, {5.5605, 0.40051}, {6.0238, 0.30019}, {6.7665, 0.20031}, {7.3642, 0.15009}, {8.3471, 0.10012}, {8.9752, 0.080004}, {9.8708, 0.059965}, {10.485, 0.050076}, {11.342, 0.040002}, {12.549, 0.029982}, {14.573, 0.020001}, {16.246, 0.014987}, {19.039, 0.01}},
+ 50: {{1.4162, 7990.5}, {1.4464, 5987.4}, {1.464, 5000}, {1.4884, 3994.1}, {1.518, 2993.7}, {1.5645, 1997}, {1.598, 1496.8}, {1.6469, 998.51}, {1.6745, 801.12}, {1.7155, 600.3}, {1.7416, 499.25}, {1.7761, 400.56}, {1.8224, 300.15}, {1.8953, 200.28}, {1.9533, 150.07}, {2.0439, 100.14}, {2.0971, 80.014}, {2.1744, 59.956}, {2.2277, 50.068}, {2.2995, 39.996}, {2.3988, 29.978}, {2.5597, 20.003}, {2.6865, 14.989}, {2.8928, 9.9988}, {3.0177, 7.9895}, {3.1958, 6.0112}, {3.3137, 4.9993}, {3.4728, 3.9947}, {3.7001, 3.0056}, {4.051, 1.9973}, {4.3291, 1.5028}, {4.7758, 0.99866}, {5.0427, 0.80124}, {5.4215, 0.60039}, {5.6812, 0.49932}, {6.0259, 0.40051}, {6.5086, 0.30019}, {7.2892, 0.20031}, {7.908, 0.15009}, {8.9377, 0.10012}, {9.5659, 0.080004}, {10.489, 0.059965}, {11.142, 0.050076}, {12.016, 0.040002}, {13.274, 0.029982}, {15.369, 0.020001}, {17.109, 0.014987}, {20.017, 0.01}},
+ 63: {{1.317, 7990.5}, {1.3512, 5987.4}, {1.3718, 5000}, {1.3991, 3994.1}, {1.4354, 2993.7}, {1.4861, 1997}, {1.5247, 1496.8}, {1.5858, 998.51}, {1.6197, 801.12}, {1.6695, 600.3}, {1.7025, 499.25}, {1.7467, 400.56}, {1.8058, 300.15}, {1.8979, 200.28}, {1.971, 150.07}, {2.0811, 100.14}, {2.1482, 80.014}, {2.2409, 59.956}, {2.3027, 50.068}, {2.3805, 39.996}, {2.491, 29.978}, {2.662, 20.003}, {2.7938, 14.989}, {2.9993, 9.9988}, {3.124, 7.9895}, {3.2936, 6.0112}, {3.4049, 4.9993}, {3.5519, 3.9947}, {3.7562, 3.0056}, {4.075, 1.9973}, {4.3223, 1.5028}, {4.7108, 0.99866}, {4.9517, 0.80124}, {5.2759, 0.60039}, {5.5037, 0.49932}, {5.8024, 0.40051}, {6.2293, 0.30019}, {6.9131, 0.20031}, {7.4784, 0.15009}, {8.3874, 0.10012}, {8.964, 0.080004}, {9.799, 0.059965}, {10.378, 0.050076}, {11.158, 0.040002}, {12.291, 0.029982}, {14.165, 0.020001}, {15.719, 0.014987}, {18.282, 0.01}},
+}
diff --git a/fuse/fuse.go b/fuse/fuse.go
new file mode 100644
index 0000000..4fd40ab
--- /dev/null
+++ b/fuse/fuse.go
@@ -0,0 +1,604 @@
+// Package fuse models the fuse on each phase of a site's main supply as a
+// thermal system: the current heats the fuse element (I², so export heats
+// it as much as import), the heat flows through the fuse body and sand to
+// the holder and the cabinet, and the element melts when it reaches its
+// melting temperature. A fuse that has carried its rating for an hour is
+// hot and blows at an overload sooner than a cold one; short overloads add
+// up; a fuse cools along the sum of its bodies' exponentials.
+//
+// Each phase has three thermal bodies in a ladder:
+//
+// heat I² → element → fuse body and sand → holder and cabinet → ambient
+//
+// In per-unit terms the heat into the element is u² (u = I/In) and the
+// three thermal resistances add up to 1, so at a steady current the element
+// sits u² above ambient. It melts at UMin², UMin being the steady current
+// (as a multiple of In) that just melts it at the reference ambient.
+//
+// The gG bands keep IEC 60269-1's conventional currents, gates and
+// pre-arcing I²t window at 0.01 s, and within them follow a manufacturer's
+// published curves (see gg). aM keeps its gates. The other profiles are gG
+// with scaled times and are not calibrated. The model is a central
+// estimate of one maker's fuse-links, not a worst case, and it is
+// independent of any EMS's own fuse protection: it is the judge.
+//
+// It has no dependencies, so an EMS's own fast simulation can import it.
+package fuse
+
+import (
+ "fmt"
+ "math"
+ "sort"
+)
+
+// Profile ids.
+const (
+ Manual = "manual" // never blows by itself
+ IecGg = "iec_gg"
+ IecGm = "iec_gm"
+ IecAm = "iec_am"
+ IecGr = "iec_gr"
+ IecGs = "iec_gs"
+)
+
+const (
+ // MeltC is the element's melting temperature: tin, 231.93 °C, the
+ // solder spot (M-effect) at which a gG element melts under overload.
+ MeltC = 232.0
+ // RefAmbientC is the ambient the calibration holds at: IEC 60269 gives
+ // time-current characteristics at a reference ambient of 20 °C,
+ // approximately valid at 30 °C (Enel GSCL004 rev. 0, 2016, §6.1).
+ RefAmbientC = 20.0
+ // DefaultAmbientC is a meter cabinet's air when a phase gives none.
+ DefaultAmbientC = 25.0
+)
+
+// Point is one point of a time-current curve from cold: at Multiple × the
+// rating the fuse melts after Seconds.
+type Point struct {
+ Multiple float64
+ Seconds float64
+}
+
+// Gate bounds the melting time from cold at a multiple of the rating: with
+// Min the fuse must not melt sooner than Seconds, otherwise it must melt
+// within Seconds.
+type Gate struct {
+ Multiple float64
+ Seconds float64
+ Min bool
+}
+
+// Thermal is a band's three bodies. Re, Rb and Rh are the shares of the
+// steady temperature rise across element-to-body, body-to-holder and
+// holder-to-ambient (they add up to 1). Each Tau is that body's heat
+// capacity times its own outgoing resistance; they are parameters, not the
+// system's modes, which come out of the coupled ladder (at 25 A 0.12 s,
+// 9.4 s and 1,308 s).
+type Thermal struct {
+ UMin float64 // steady melting multiple at the reference ambient
+ Re, TauE float64 // element
+ Rb, TauB float64 // fuse body and sand
+ Rh, TauH float64 // holder and cabinet
+ TimeFactor float64 // scales every time constant; 0 or 1 as calibrated
+}
+
+// Band is a fuse type at one rating: its thermal bodies and the points it
+// is checked against.
+type Band struct {
+ RatingA float64
+ Thermal Thermal
+ // Gates are the standard's bounds; I2tMin and I2tMax bound the
+ // pre-arcing I²t at 0.01 s (A²s). Reference is the manufacturer's curve
+ // the band is fitted to.
+ Gates []Gate
+ I2tMin, I2tMax float64
+ Reference []Point
+
+ sys system
+}
+
+// Profile is a fuse type over its rating bands.
+type Profile struct {
+ ID string
+ Label string
+ Source string // where its numbers come from
+ Bands []Band
+ // Curve is the reference band's cold curve, sampled, for display.
+ Curve []Point
+ // RefRatingA is the band TripSeconds and Curve use.
+ RefRatingA float64
+}
+
+// gg are the gG bands, 16 to 63 A.
+//
+// Gates, I²t windows and conventional currents are IEC 60269-1's, as tabled
+// in Schneider Electric, "Complementary technical information" CA908036E
+// v3.2 (2018); the 25, 40, 50 and 63 A rows are also in Enel GSCL004 rev. 0
+// (2016) Tables 6, 7 and 14, and 25 A in Mersen EduPack TM-104 (2012) slide
+// 6. 35 A is not in those tables: its gates and I²t are interpolated on log
+// axes between 32 A and 40 A. Conventional currents: Inf = 1.25 In must not
+// melt and If = 1.6 In must melt within the conventional time, 1 h for
+// 16 ≤ In ≤ 63 A.
+//
+// The gates leave a wide tube: no gate constrains 1.4-2 ×, where home
+// overloads land. Within it each band follows ABB's pre-arcing curves for
+// OFAF gG fuse-links, size 000, 500 V (ABB "Fusegear DIN-type HRC-fuse
+// links" 1SCC317001C0201, page 14), digitized from the PDF's vector paths
+// (48 vertices a curve; read-off error under 2 % in time). All six thermal
+// parameters are fitted, least squares on log time from cold at 20 °C, to
+// ABB's melting times at 1.5, 1.6, 2, 2.5, 3, 4, 5, 8 and 10 × In, the
+// current that melts in 0.01 s and the curve's start at 7,990 s; the fit
+// keeps every gate. At those eleven points the residuals are rms
+// 0.08-0.095 decades, at most 0.17. Against all 48 vertices of each ABB
+// curve they are rms 0.13-0.20, at most 0.22-0.46: the model is too slow
+// at 1.6-1.95 × (16 A at 1.82 ×: ABB 150 s, model 447 s) and too fast at
+// 2.0-2.5 × (down to -0.27). Where the element's quasi-steady rise over
+// the body reaches the melting rise, u = UMin/√(Re + Rb) = 2.02-2.09 × In,
+// the body stage takes over from the holder and the curve drops steeply:
+// its log-log slope between 1.9 and 2.3 × reaches -18 to -52, ABB's -7 to
+// -11. Fitting all 48 vertices reaches only rms 0.12-0.15, so it is the
+// three-body ladder's shape. A simulation judging overloads at 1.9-2.1 ×
+// In sits on that step: a small change of current there moves the melting
+// time a lot.
+// Letting the element's resistance rise with temperature (silver, 0.0038
+// per K) made the fit worse (rms 0.11-0.14), so it is left out.
+var gg = []struct {
+ ratingA float64
+ imin10, imax5, imin01, imax1 float64 // A
+ i2tMin, i2tMax float64 // A²s at 0.01 s
+ uMin, re, tauE, rh float64 // fitted
+ tauB, tauH float64 // fitted, s
+ abb [9]float64
+ abb001A, abbStartA float64 // A melting in 0.01 s; A at the curve's 7,990 s start
+}{
+ {16, 33, 65, 85, 150, 300, 1000, 1.429959, 0.073378, 0.074249, 0.500174, 7.03102, 1604.18,
+ [9]float64{3546.8, 1165.2, 53.1, 9.445, 3.29, 0.8381, 0.3328, 0.0554, 0.0246}, 208.1, 22.918},
+ {20, 42, 85, 110, 200, 500, 1800, 1.412945, 0.070246, 0.091392, 0.514825, 9.47660, 1251.63,
+ [9]float64{2592.9, 819.5, 57.8, 12.3, 4.58, 1.177, 0.4498, 0.0711, 0.0318}, 284.5, 28.272},
+ {25, 52, 110, 150, 260, 1000, 3000, 1.396178, 0.073599, 0.128139, 0.527645, 8.81978, 1295.49,
+ [9]float64{2341.3, 821.6, 58.7, 11.8, 4.222, 1.098, 0.445, 0.0847, 0.0419}, 406.9, 34.925},
+ {32, 75, 150, 200, 350, 1800, 5000, 1.392001, 0.075676, 0.167797, 0.532203, 10.60002, 1500.19,
+ [9]float64{2581.2, 970.5, 72.9, 14.1, 5.056, 1.311, 0.5272, 0.1026, 0.0519}, 591.8, 44.605},
+ {35, 82.47, 164.94, 222.22, 387.17, 2209.86, 6331.18, 1.390453, 0.075385, 0.169990, 0.533283, 10.92408, 1601.99, // gates interpolated
+ [9]float64{2680.3, 1049.8, 76.5, 14.6, 5.206, 1.342, 0.5387, 0.1044, 0.0527}, 650.8, 48.759},
+ {40, 95, 190, 260, 450, 3000, 9000, 1.396104, 0.074093, 0.177481, 0.529719, 11.81891, 1595.82,
+ [9]float64{2778.0, 1066.7, 80.2, 15.8, 5.714, 1.494, 0.5982, 0.1149, 0.0577}, 761.5, 55.948},
+ {50, 125, 250, 350, 610, 5000, 16000, 1.409864, 0.071730, 0.197218, 0.525250, 16.77090, 1948.86,
+ [9]float64{3566.2, 1471.4, 121.6, 23.2, 8.242, 2.115, 0.8293, 0.1445, 0.0696}, 1000.8, 70.810},
+ {63, 160, 320, 450, 820, 9000, 27000, 1.314326, 0.063439, 0.165090, 0.604208, 19.24866, 1554.48,
+ [9]float64{1798.4, 910.0, 134.6, 29.3, 9.987, 2.193, 0.7666, 0.1183, 0.0563}, 1151.8, 82.971},
+}
+
+var abbMultiples = [9]float64{1.5, 1.6, 2, 2.5, 3, 4, 5, 8, 10}
+
+func ggBands(timeFactor float64, calibrated bool) []Band {
+ var out []Band
+ for _, r := range gg {
+ th := Thermal{UMin: r.uMin, Re: r.re, TauE: r.tauE, Rb: 1 - r.re - r.rh, TauB: r.tauB,
+ Rh: r.rh, TauH: r.tauH, TimeFactor: timeFactor}
+ b := Band{RatingA: r.ratingA, Thermal: th}
+ if calibrated {
+ in := r.ratingA
+ b.Gates = []Gate{
+ {1.25, 3600, true}, {1.6, 3600, false},
+ {r.imin10 / in, 10, true}, {r.imax5 / in, 5, false},
+ {r.imin01 / in, 0.1, true}, {r.imax1 / in, 0.1, false},
+ }
+ b.I2tMin, b.I2tMax = r.i2tMin, r.i2tMax
+ b.Reference = append(b.Reference, Point{r.abbStartA / in, 7990})
+ for i, m := range abbMultiples {
+ b.Reference = append(b.Reference, Point{m, r.abb[i]})
+ }
+ b.Reference = append(b.Reference, Point{r.abb001A / in, 0.01})
+ }
+ out = append(out, b)
+ }
+ return out
+}
+
+// amHolderTauS is aM's holder and cabinet, as gG's were before they were
+// fitted: a modelling choice.
+const amHolderTauS = 1200.0
+
+// amBand is aM at any rating. Gates from Mersen EduPack TM-104 (2012) slide
+// 7, "as per IEC 60269", for any current rating: not melting within 60 s at
+// 4 In, 0.5 s at 8 In and 0.2 s at 10 In; melting within 60 s at 6.3 In,
+// 0.5 s at 12.5 In and 0.1 s at 19 In. Fitted to √(4 × 6.3) In at 60 s,
+// 10 In at √(0.2 × 0.5) s and 19 In at 0.05 s, with UMin = 4 (an
+// assumption: aM protects against short circuits, not overloads) and the
+// holder at a 0.2 share, which the fit needs.
+func amBand() Band {
+ return Band{RatingA: 0,
+ Thermal: Thermal{UMin: 4, Re: 0.2005, TauE: 0.2001, Rb: 0.5995, TauB: 45.9302, Rh: 0.2, TauH: amHolderTauS},
+ Gates: []Gate{{4, 60, true}, {6.3, 60, false}, {8, 0.5, true}, {10, 0.2, true},
+ {12.5, 0.5, false}, {19, 0.1, false}},
+ Reference: []Point{{math.Sqrt(4 * 6.3), 60}, {10, math.Sqrt(0.2 * 0.5)}, {19, 0.05}},
+ }
+}
+
+var profiles = func() map[string]Profile {
+ m := map[string]Profile{
+ Manual: {ID: Manual, Label: "Manual", Source: "never blows by itself"},
+ IecGg: {ID: IecGg, Label: "IEC 60269 gG/gL (general purpose)", RefRatingA: 25, Bands: ggBands(1, true),
+ Source: "IEC 60269-1 gates, I²t and conventional currents; ABB OFAF 000 gG curves within them (see gg)"},
+ IecAm: {ID: IecAm, Label: "IEC 60269 aM (motor, short-circuit)", RefRatingA: 25, Bands: []Band{amBand()},
+ Source: "IEC 60269 aM gates (see amBand)"},
+ // Not calibrated: gG with every time scaled, as the simulator's
+ // earlier curves had them against gG.
+ IecGm: {ID: IecGm, Label: "IEC 60269 gM (motor, dual-rated) - UK common", RefRatingA: 25, Bands: ggBands(2, false),
+ Source: "not calibrated: gG with times × 2 (a gM is a gG of its characteristic current, which the site does not give)"},
+ IecGr: {ID: IecGr, Label: "IEC 60269 gR (semiconductor, fast)", RefRatingA: 25, Bands: ggBands(0.3, false),
+ Source: "not calibrated: gG with times × 0.3"},
+ IecGs: {ID: IecGs, Label: "IEC 60269 gS (semiconductor, full-range)", RefRatingA: 25, Bands: ggBands(0.5, false),
+ Source: "not calibrated: gG with times × 0.5"},
+ }
+ for id, p := range m {
+ for i := range p.Bands {
+ p.Bands[i].sys = newSystem(p.Bands[i].Thermal)
+ }
+ if b, err := p.band(p.RefRatingA); err == nil {
+ for _, x := range []float64{1.5, 1.6, 2, 3, 5, 10, 20} {
+ if t := b.MeltSeconds(x, RefAmbientC); !math.IsInf(t, 1) {
+ p.Curve = append(p.Curve, Point{x, t})
+ }
+ }
+ }
+ m[id] = p
+ }
+ return m
+}()
+
+// Lookup returns the profile with this id.
+func Lookup(id string) (Profile, bool) {
+ p, ok := profiles[id]
+ return p, ok
+}
+
+// All returns every profile, ordered by id.
+func All() []Profile {
+ out := make([]Profile, 0, len(profiles))
+ for _, p := range profiles {
+ out = append(out, p)
+ }
+ sort.Slice(out, func(i, j int) bool { return out[i].ID < out[j].ID })
+ return out
+}
+
+// Ratings are the ratings a gG-based profile is calibrated for; between
+// them a rating takes the nearest band on log axes.
+func (p Profile) Ratings() []float64 {
+ var out []float64
+ for _, b := range p.Bands {
+ if b.RatingA > 0 {
+ out = append(out, b.RatingA)
+ }
+ }
+ return out
+}
+
+// band is the band for ratingA: the nearest on log axes between the first
+// and the last rating, any rating for a rating-independent band (aM), an
+// error outside the bands or for a profile without any.
+func (p Profile) band(ratingA float64) (*Band, error) {
+ if len(p.Bands) == 0 {
+ return nil, fmt.Errorf("%s has no fuse model", p.ID)
+ }
+ if !(ratingA > 0) || math.IsInf(ratingA, 0) {
+ return nil, fmt.Errorf("fuse rating %v A is not a rating", ratingA)
+ }
+ if len(p.Bands) == 1 && p.Bands[0].RatingA == 0 {
+ return &p.Bands[0], nil
+ }
+ lo, hi := p.Bands[0].RatingA, p.Bands[len(p.Bands)-1].RatingA
+ if ratingA < lo || ratingA > hi {
+ return nil, fmt.Errorf("%s is modelled for %g-%g A, not %g A", p.ID, lo, hi, ratingA)
+ }
+ best, dist := 0, math.Inf(1)
+ for i, b := range p.Bands {
+ if d := math.Abs(math.Log(ratingA / b.RatingA)); d < dist {
+ best, dist = i, d
+ }
+ }
+ return &p.Bands[best], nil
+}
+
+// Band returns the band a fuse rated ratingA uses, and false when the
+// profile has none for it.
+func (p Profile) Band(ratingA float64) (Band, bool) {
+ b, err := p.band(ratingA)
+ if err != nil {
+ return Band{}, false
+ }
+ return *b, true
+}
+
+// Supports says whether the profile models a fuse rated ratingA.
+func (p Profile) Supports(ratingA float64) bool {
+ _, err := p.band(ratingA)
+ return err == nil
+}
+
+// TripSeconds is the time to melt from cold at multiple × the rating, at the
+// reference band and ambient; +Inf if it never melts there.
+func (p Profile) TripSeconds(multiple float64) float64 {
+ return p.MeltSeconds(multiple, p.RefRatingA, RefAmbientC)
+}
+
+// MeltSeconds is the time to melt from cold at multiple × a rating of
+// ratingA in ambientC; +Inf if it never melts or the rating has no band.
+func (p Profile) MeltSeconds(multiple, ratingA, ambientC float64) float64 {
+ b, err := p.band(ratingA)
+ if err != nil {
+ return math.Inf(1)
+ }
+ return b.MeltSeconds(multiple, ambientC)
+}
+
+// MeltSeconds is the band's time to melt from cold.
+func (b *Band) MeltSeconds(multiple, ambientC float64) float64 {
+ return b.sys.timeToMelt([3]float64{}, multiple, b.sys.meltRise(ambientC))
+}
+
+// Phase is the state of one phase's fuse.
+type Phase struct {
+ // Theta is each body's rise over ambient, in per-unit (element, body,
+ // holder); zero is cold.
+ Theta [3]float64
+ // AmbientC is the cabinet's air; nil means DefaultAmbientC. 0 is 0 °C.
+ AmbientC *float64
+ // Damage is the element's way from ambient to melting, 0 to 1; Step
+ // sets it.
+ Damage float64
+ Blown bool
+}
+
+// Ambient returns a pointer to c, for Phase.AmbientC.
+func Ambient(c float64) *float64 { return &c }
+
+func (ph *Phase) ambient() float64 {
+ if ph.AmbientC == nil {
+ return DefaultAmbientC
+ }
+ return *ph.AmbientC
+}
+
+// Status describes a phase after a step.
+type Status struct {
+ Ratio float64 `json:"ratio"` // current over the rating
+ TripSeconds float64 `json:"trip_s"` // time to melt at this ratio from cold; 0 if never
+ RemainingSeconds float64 `json:"remaining_s"` // time to melt at this ratio from now; 0 if never
+ Damage float64 `json:"damage"` // the element's way from ambient to melting
+ Active bool `json:"active"` // it will melt if the current stays
+ Invalid string `json:"invalid,omitempty"` // why the step was refused; the state is unchanged
+}
+
+// Step advances one phase by dt seconds at amps through a fuse rated
+// ratingA. The sign of amps does not matter, and an infinite current blows
+// the fuse at once. The step is exact for a current held over dt, so the
+// result does not depend on the step size; dt = 0 changes nothing and asks
+// how long the fuse has left at that current.
+//
+// Step refuses, leaving the phase as it was and saying why in
+// Status.Invalid: a NaN current, a negative, NaN or infinite dt, an ambient
+// that is not finite or not below the melting point, and a rating the
+// profile does not model (the profile without a model, manual, refuses
+// nothing and never blows).
+func (p Profile) Step(ph *Phase, amps, ratingA, dt float64) Status {
+ if len(p.Bands) == 0 {
+ return Status{Damage: ph.Damage}
+ }
+ b, err := p.band(ratingA)
+ switch {
+ case err != nil:
+ return Status{Damage: ph.Damage, Invalid: err.Error()}
+ case math.IsNaN(amps):
+ return Status{Damage: ph.Damage, Invalid: "current is NaN"}
+ case !(dt >= 0) || math.IsInf(dt, 1):
+ return Status{Damage: ph.Damage, Invalid: fmt.Sprintf("step of %v s", dt)}
+ }
+ ambient := ph.ambient()
+ if math.IsNaN(ambient) || math.IsInf(ambient, 0) || ambient >= MeltC {
+ return Status{Damage: ph.Damage, Invalid: fmt.Sprintf("ambient of %v °C", ambient)}
+ }
+ if ph.Blown {
+ ph.Damage = 1
+ return Status{Damage: 1}
+ }
+ u := math.Abs(amps) / ratingA
+ st := Status{Ratio: u}
+ if math.IsInf(u*u, 1) {
+ // A current whose square overflows (or an infinite one) melts the
+ // fuse at once; the ratio stays finite for JSON.
+ ph.Blown, ph.Damage = true, 1
+ st.Ratio, st.Damage = math.MaxFloat64, 1
+ return st
+ }
+ melt := b.sys.meltRise(ambient)
+ if cold := b.sys.timeToMelt([3]float64{}, u, melt); !math.IsInf(cold, 1) {
+ st.TripSeconds = cold
+ }
+ left := b.sys.timeToMelt(ph.Theta, u, melt)
+ if left <= dt {
+ ph.Theta = b.sys.advance(ph.Theta, u, left)
+ ph.Blown, ph.Damage = true, 1
+ st.Damage = 1
+ return st
+ }
+ ph.Theta = b.sys.advance(ph.Theta, u, dt)
+ ph.Damage = math.Max(0, math.Min(1, ph.Theta[0]/melt))
+ st.Damage = ph.Damage
+ if !math.IsInf(left, 1) {
+ st.RemainingSeconds, st.Active = left-dt, true
+ }
+ return st
+}
+
+// Preheat puts the phase in the steady state of a current held for ever in
+// its ambient, as a fuse that has carried the evening's load: a scenario
+// that starts mid-evening starts warm, not cold. It refuses a blown fuse
+// (clear Blown and Theta to replace one), a current that would melt the
+// fuse before reaching that state, and anything Step refuses, leaving the
+// phase as it was.
+func (p Profile) Preheat(ph *Phase, amps, ratingA float64) error {
+ if len(p.Bands) == 0 {
+ return fmt.Errorf("%s has no fuse model", p.ID)
+ }
+ if ph.Blown {
+ return fmt.Errorf("the fuse is blown")
+ }
+ if st := p.Step(&Phase{AmbientC: ph.AmbientC}, amps, ratingA, 0); st.Invalid != "" {
+ return fmt.Errorf("%s", st.Invalid)
+ }
+ b, _ := p.band(ratingA)
+ u := math.Abs(amps) / ratingA
+ inf, _ := b.sys.modes([3]float64{}, u)
+ if inf[0] >= b.sys.meltRise(ph.ambient()) {
+ return fmt.Errorf("%.3g × %g A melts the fuse before it settles", u, ratingA)
+ }
+ ph.Theta = inf
+ ph.Damage = inf[0] / b.sys.meltRise(ph.ambient())
+ return nil
+}
+
+// system is a band's ladder solved once: C dθ/dt = -G θ + (u², 0, 0).
+// With y = θ/s (s = 1/√C), dy/dt = -M y + ..., M = s G s symmetric, so
+// θ(t) = θ∞ + P e^{-Λt} Q (θ0 - θ∞), P = diag(s) V, Q = Vᵀ diag(1/s).
+type system struct {
+ uMin2 float64
+ inf [3]float64 // θ∞ at u = 1
+ lambda [3]float64
+ p, q [3][3]float64
+}
+
+func newSystem(t Thermal) system {
+ f := t.TimeFactor
+ if f <= 0 {
+ f = 1
+ }
+ c := [3]float64{f * t.TauE / t.Re, f * t.TauB / t.Rb, f * t.TauH / t.Rh}
+ g1, g2, g3 := 1/t.Re, 1/t.Rb, 1/t.Rh
+ g := [3][3]float64{{g1, -g1, 0}, {-g1, g1 + g2, -g2}, {0, -g2, g2 + g3}}
+ var s [3]float64
+ for i := range s {
+ s[i] = 1 / math.Sqrt(c[i])
+ }
+ var m [3][3]float64
+ for i := 0; i < 3; i++ {
+ for j := 0; j < 3; j++ {
+ m[i][j] = s[i] * g[i][j] * s[j]
+ }
+ }
+ lam, v := eigen3(m)
+ sys := system{uMin2: t.UMin * t.UMin, lambda: lam,
+ inf: [3]float64{t.Re + t.Rb + t.Rh, t.Rb + t.Rh, t.Rh}}
+ for i := 0; i < 3; i++ {
+ for j := 0; j < 3; j++ {
+ sys.p[i][j] = s[i] * v[i][j]
+ sys.q[i][j] = v[j][i] / s[j]
+ }
+ }
+ return sys
+}
+
+// meltRise is the element's melting rise over ambient, in per-unit: UMin²
+// at the reference ambient, less in a warmer cabinet.
+func (s *system) meltRise(ambientC float64) float64 {
+ return s.uMin2 * (MeltC - ambientC) / (MeltC - RefAmbientC)
+}
+
+// modes returns θ∞ at u and the modal amplitudes of θ0 - θ∞.
+func (s *system) modes(theta [3]float64, u float64) (inf, amp [3]float64) {
+ u2 := u * u
+ var d [3]float64
+ for i := range inf {
+ inf[i] = s.inf[i] * u2
+ d[i] = theta[i] - inf[i]
+ }
+ for k := 0; k < 3; k++ {
+ for j := 0; j < 3; j++ {
+ amp[k] += s.q[k][j] * d[j]
+ }
+ }
+ return inf, amp
+}
+
+func (s *system) advance(theta [3]float64, u, dt float64) [3]float64 {
+ inf, amp := s.modes(theta, u)
+ var out [3]float64
+ for i := 0; i < 3; i++ {
+ out[i] = inf[i]
+ for k := 0; k < 3; k++ {
+ out[i] += s.p[i][k] * math.Exp(-s.lambda[k]*dt) * amp[k]
+ }
+ }
+ return out
+}
+
+// timeToMelt is when the element first reaches melt at u from theta; +Inf
+// if it never does (within ten million seconds).
+func (s *system) timeToMelt(theta [3]float64, u, melt float64) float64 {
+ inf, amp := s.modes(theta, u)
+ f := func(t float64) float64 {
+ x := inf[0] - melt
+ for k := 0; k < 3; k++ {
+ x += s.p[0][k] * math.Exp(-s.lambda[k]*t) * amp[k]
+ }
+ return x
+ }
+ if f(0) >= 0 {
+ return 0
+ }
+ prev := 0.0
+ for t := 1e-6; t <= 1e7; t *= 1.2 {
+ if f(t) >= 0 {
+ lo, hi := prev, t
+ for i := 0; i < 200 && hi-lo > 1e-12*hi; i++ {
+ mid := (lo + hi) / 2
+ if f(mid) >= 0 {
+ hi = mid
+ } else {
+ lo = mid
+ }
+ }
+ return hi
+ }
+ prev = t
+ }
+ return math.Inf(1)
+}
+
+// eigen3 diagonalises a symmetric 3×3 matrix by Jacobi rotations: the
+// eigenvalues and the eigenvectors as columns.
+func eigen3(a [3][3]float64) ([3]float64, [3][3]float64) {
+ v := [3][3]float64{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}
+ for sweep := 0; sweep < 50; sweep++ {
+ off := a[0][1]*a[0][1] + a[0][2]*a[0][2] + a[1][2]*a[1][2]
+ if off < 1e-30 {
+ break
+ }
+ for _, pq := range [][2]int{{0, 1}, {0, 2}, {1, 2}} {
+ p, q := pq[0], pq[1]
+ if math.Abs(a[p][q]) < 1e-300 {
+ continue
+ }
+ theta := (a[q][q] - a[p][p]) / (2 * a[p][q])
+ t := math.Copysign(1, theta) / (math.Abs(theta) + math.Sqrt(theta*theta+1))
+ c := 1 / math.Sqrt(t*t+1)
+ sn := t * c
+ for k := 0; k < 3; k++ {
+ akp, akq := a[k][p], a[k][q]
+ a[k][p], a[k][q] = c*akp-sn*akq, sn*akp+c*akq
+ }
+ for k := 0; k < 3; k++ {
+ apk, aqk := a[p][k], a[q][k]
+ a[p][k], a[q][k] = c*apk-sn*aqk, sn*apk+c*aqk
+ }
+ for k := 0; k < 3; k++ {
+ vkp, vkq := v[k][p], v[k][q]
+ v[k][p], v[k][q] = c*vkp-sn*vkq, sn*vkp+c*vkq
+ }
+ }
+ }
+ return [3]float64{a[0][0], a[1][1], a[2][2]}, v
+}
diff --git a/fuse/fuse_test.go b/fuse/fuse_test.go
new file mode 100644
index 0000000..8ed3439
--- /dev/null
+++ b/fuse/fuse_test.go
@@ -0,0 +1,381 @@
+package fuse
+
+import (
+ "encoding/json"
+ "math"
+ "testing"
+)
+
+func profile(t *testing.T, id string) Profile {
+ t.Helper()
+ p, ok := Lookup(id)
+ if !ok {
+ t.Fatalf("no %s profile", id)
+ }
+ return p
+}
+
+func ref() *Phase { return &Phase{AmbientC: Ambient(RefAmbientC)} }
+
+// hold runs a phase at multiple × rating for seconds in steps of dt and
+// says how long it took to blow (+Inf if it held).
+func hold(p Profile, ph *Phase, multiple, rating, seconds, dt float64) float64 {
+ for t := 0.0; t < seconds; t += dt {
+ p.Step(ph, multiple*rating, rating, dt)
+ if ph.Blown {
+ return t + dt
+ }
+ }
+ return math.Inf(1)
+}
+
+// Every gG band, from cold at the reference ambient, keeps IEC 60269-1's
+// gates and its pre-arcing I²t window at 0.01 s, and sits within 0.2
+// decades of the ABB curve it was fitted to.
+func TestGGBandsMeetTheGatesAndFollowABB(t *testing.T) {
+ p := profile(t, IecGg)
+ for _, b := range p.Bands {
+ for _, g := range b.Gates {
+ got := b.MeltSeconds(g.Multiple, RefAmbientC)
+ if g.Min && got < g.Seconds || !g.Min && got > g.Seconds {
+ kind := "within"
+ if g.Min {
+ kind = "not sooner than"
+ }
+ t.Errorf("%g A at %.2f×: melts after %.3g s, want %s %g s", b.RatingA, g.Multiple, got, kind, g.Seconds)
+ }
+ }
+ u := multipleAt(&b, 0.01)
+ if i2t := math.Pow(u*b.RatingA, 2) * 0.01; i2t < b.I2tMin || i2t > b.I2tMax {
+ t.Errorf("%g A: pre-arcing I²t at 0.01 s %.0f A²s, want %g-%g", b.RatingA, i2t, b.I2tMin, b.I2tMax)
+ }
+ if len(b.Reference) != 11 {
+ t.Fatalf("%g A: %d reference points", b.RatingA, len(b.Reference))
+ }
+ for _, pt := range b.Reference {
+ got := b.MeltSeconds(pt.Multiple, RefAmbientC)
+ if d := math.Abs(math.Log10(got / pt.Seconds)); d > 0.2 {
+ t.Errorf("%g A at %.2f×: %.3g s, %.2f decades from ABB's %.3g s", b.RatingA, pt.Multiple, got, d, pt.Seconds)
+ }
+ }
+ }
+}
+
+// multipleAt is the multiple at which the band melts after seconds.
+func multipleAt(b *Band, seconds float64) float64 {
+ lo, hi := 1.0, 1000.0
+ for i := 0; i < 100; i++ {
+ mid := math.Sqrt(lo * hi)
+ if b.MeltSeconds(mid, RefAmbientC) > seconds {
+ lo = mid
+ } else {
+ hi = mid
+ }
+ }
+ return hi
+}
+
+// Through Step, as a site runs it: 1.25 × In holds for the conventional
+// hour, 1.6 × In blows within it.
+func TestConventionalCurrents(t *testing.T) {
+ p := profile(t, IecGg)
+ for _, b := range p.Bands {
+ if got := hold(p, ref(), 1.25, b.RatingA, 3600, 1); !math.IsInf(got, 1) {
+ t.Errorf("%g A: 1.25× blew after %.0f s, want it to hold 1 h", b.RatingA, got)
+ }
+ if got := hold(p, ref(), 1.6, b.RatingA, 3600, 1); got > 3600 {
+ t.Errorf("%g A: 1.6× held 1 h, want it to blow", b.RatingA)
+ }
+ }
+}
+
+// A fuse that has carried its rating for an hour blows at 2 × in under
+// half its cold time; at 1.5 × in under three quarters.
+func TestPreloadBlowsSooner(t *testing.T) {
+ p := profile(t, IecGg)
+ for _, rating := range []float64{16, 25, 63} {
+ for _, c := range []struct{ multiple, factor float64 }{{2, 0.5}, {1.5, 0.75}} {
+ cold := hold(p, ref(), c.multiple, rating, 7200, 0.5)
+ warm := ref()
+ hold(p, warm, 1, rating, 3600, 1)
+ hot := hold(p, warm, c.multiple, rating, 7200, 0.5)
+ t.Logf("%g A, %g×: cold %.1f s, after 1 h at In %.1f s", rating, c.multiple, cold, hot)
+ if hot >= c.factor*cold {
+ t.Errorf("%g A, %g× after 1 h at In: %.1f s, want under %.0f %% of the cold %.1f s", rating, c.multiple, hot, 100*c.factor, cold)
+ }
+ }
+ }
+}
+
+// Preheat sets the steady state of a held current: the same as a long soak,
+// and a preheated fuse blows sooner than a cold one. It refuses a current
+// that melts the fuse before settling.
+func TestPreheat(t *testing.T) {
+ p := profile(t, IecGg)
+ warm := ref()
+ if err := p.Preheat(warm, 25, 25); err != nil {
+ t.Fatal(err)
+ }
+ soaked := ref()
+ hold(p, soaked, 1, 25, 6*3600, 10)
+ if math.Abs(warm.Damage-soaked.Damage) > 0.01 {
+ t.Fatalf("preheated at In damage %.3f, after a 6 h soak %.3f", warm.Damage, soaked.Damage)
+ }
+ if hot, cold := hold(p, warm, 2, 25, 7200, 0.5), p.MeltSeconds(2, 25, RefAmbientC); hot >= cold/2 {
+ t.Errorf("2× from preheated: %.1f s, want under half the cold %.1f s", hot, cold)
+ }
+ for _, amps := range []float64{1.6 * 25, math.NaN(), math.Inf(1)} {
+ ph := ref()
+ if err := p.Preheat(ph, amps, 25); err == nil || ph.Theta != ([3]float64{}) {
+ t.Errorf("preheat at %v A: err %v, state %v; want refused, cold", amps, err, ph.Theta)
+ }
+ }
+ // A blown fuse stays blown: preheating does not replace it.
+ blown := ref()
+ hold(p, blown, 3, 25, 600, 1)
+ before := *blown
+ if err := p.Preheat(blown, 25, 25); err == nil || !blown.Blown || blown.Theta != before.Theta {
+ t.Errorf("preheat on a blown fuse: err %v, blown %v; want refused, still blown", err, blown.Blown)
+ }
+}
+
+// 30 s at 1.5 × every 2 min (0.5 × between) heats a fuse further each
+// cycle: the element ends each burst hotter than the last, and after ten
+// cycles a held 1.5 × blows sooner than from cold.
+func TestShortOverloadsAccumulate(t *testing.T) {
+ p := profile(t, IecGg)
+ ph := ref()
+ var peaks []float64
+ for cycle := 0; cycle < 10; cycle++ {
+ hold(p, ph, 1.5, 25, 30, 1)
+ peaks = append(peaks, ph.Damage)
+ hold(p, ph, 0.5, 25, 90, 1)
+ }
+ for i := 1; i < len(peaks); i++ {
+ if peaks[i] <= peaks[i-1] {
+ t.Fatalf("burst %d left damage %.3f, not above the last %.3f", i+1, peaks[i], peaks[i-1])
+ }
+ }
+ t.Logf("damage after each burst: %.3f", peaks)
+ after := hold(p, ph, 1.5, 25, 7200, 1)
+ cold := p.MeltSeconds(1.5, 25, RefAmbientC)
+ t.Logf("1.5× after ten bursts blows in %.0f s, from cold %.0f s", after, cold)
+ if after >= 0.9*cold {
+ t.Errorf("1.5× after ten bursts: %.0f s, want under 90 %% of the cold %.0f s", after, cold)
+ }
+}
+
+// With no current a hot fuse cools along its slowest mode: after an hour,
+// each further 10 minutes take off the same share.
+func TestCoolingIsExponential(t *testing.T) {
+ p := profile(t, IecGg)
+ ph := ref()
+ hold(p, ph, 1.3, 25, 7200, 1)
+ hold(p, ph, 0, 25, 3600, 1)
+ var ratios []float64
+ for i := 0; i < 4; i++ {
+ before := ph.Theta[0]
+ hold(p, ph, 0, 25, 600, 1)
+ ratios = append(ratios, ph.Theta[0]/before)
+ }
+ for _, r := range ratios[1:] {
+ if math.Abs(r/ratios[0]-1) > 0.01 {
+ t.Fatalf("10-minute cooling ratios %.4f, want one constant share", ratios)
+ }
+ }
+ b, _ := p.Band(25)
+ slow := math.Min(b.sys.lambda[0], math.Min(b.sys.lambda[1], b.sys.lambda[2]))
+ if want := math.Exp(-600 * slow); math.Abs(ratios[0]-want) > 0.01 {
+ t.Errorf("cooling ratio %.4f per 10 min, want e^(-600/τ) = %.4f for the slowest mode", ratios[0], want)
+ }
+}
+
+func TestExportHeatsAsImport(t *testing.T) {
+ p := profile(t, IecGg)
+ var in, out Phase
+ for s := 0; s < 120; s++ {
+ a := p.Step(&in, 40, 25, 1)
+ b := p.Step(&out, -40, 25, 1)
+ if a != b || in.Theta != out.Theta || in.Damage != out.Damage {
+ t.Fatalf("second %d: import %+v, export %+v", s, a, b)
+ }
+ }
+}
+
+// The step is exact: steps of 0.5, 1 and 2 s reach the same state, and the
+// time left (a dt = 0 probe) drops by exactly the time stepped.
+func TestStepIsExact(t *testing.T) {
+ p := profile(t, IecGg)
+ for _, m := range []float64{1.7, 2.5, 6} {
+ cold := p.MeltSeconds(m, 25, RefAmbientC)
+ var states [][3]float64
+ for _, dt := range []float64{0.5, 1, 2} {
+ ph := ref()
+ elapsed := 0.0
+ for elapsed+dt < math.Min(cold, 60) {
+ p.Step(ph, m*25, 25, dt)
+ elapsed += dt
+ }
+ left := p.Step(ph, m*25, 25, 0).RemainingSeconds
+ if math.Abs(left-(cold-elapsed)) > 1e-9*cold {
+ t.Errorf("%g× in %g s steps: %.12g s left after %g s, want %.12g", m, dt, left, elapsed, cold-elapsed)
+ }
+ ph2 := ref()
+ for i := 0; i < int(math.Round(4/dt)); i++ {
+ p.Step(ph2, m*25, 25, dt)
+ }
+ states = append(states, ph2.Theta)
+ }
+ for _, s := range states[1:] {
+ for i := range s {
+ if math.Abs(s[i]-states[0][i]) > 1e-9*math.Abs(states[0][i]) {
+ t.Errorf("%g× after 4 s: state %v, with 0.5 s steps %v", m, s, states[0])
+ }
+ }
+ }
+ }
+}
+
+// Status: trip_s is the cold time, remaining_s the time from now, damage
+// the element's way from ambient to melting.
+func TestStatusDescribesThePhase(t *testing.T) {
+ p := profile(t, IecGg)
+ ph := ref()
+ first := p.Step(ph, 50, 25, 1)
+ second := p.Step(ph, 50, 25, 1)
+ cold := p.MeltSeconds(2, 25, RefAmbientC)
+ if !first.Active || math.Abs(first.TripSeconds-cold) > 1e-9*cold || math.Abs(first.RemainingSeconds-(cold-1)) > 1e-9*cold {
+ t.Fatalf("first second at 2×: %+v, want active, trip %.2f s, remaining %.2f s", first, cold, cold-1)
+ }
+ if math.Abs(second.RemainingSeconds-(first.RemainingSeconds-1)) > 1e-9*cold || second.Damage <= first.Damage {
+ t.Fatalf("second second: %+v after %+v", second, first)
+ }
+ idle := p.Step(ref(), 25, 25, 1)
+ if idle.Active || idle.TripSeconds != 0 || idle.RemainingSeconds != 0 || idle.Damage <= 0 {
+ t.Fatalf("at the rating: %+v, want warming but never melting", idle)
+ }
+ if hold(p, ph, 2, 25, 7200, 1); !ph.Blown || ph.Damage != 1 {
+ t.Fatalf("held at 2×: %+v, want blown", ph)
+ }
+ if st := p.Step(ph, 0, 25, 1); st.Damage != 1 || !ph.Blown {
+ t.Fatalf("a blown fuse: %+v", st)
+ }
+}
+
+// 0 °C is a temperature: a cold cabinet, not the default.
+func TestZeroAmbientIsZero(t *testing.T) {
+ p := profile(t, IecGg)
+ at0 := p.Step(&Phase{AmbientC: Ambient(0)}, 1.5*25, 25, 0).TripSeconds
+ unset := p.Step(&Phase{}, 1.5*25, 25, 0).TripSeconds
+ if want := p.MeltSeconds(1.5, 25, 0); at0 != want {
+ t.Errorf("Step at 0 °C: %.0f s, MeltSeconds at 0 °C %.0f s", at0, want)
+ }
+ if want := p.MeltSeconds(1.5, 25, DefaultAmbientC); unset != want || at0 <= unset {
+ t.Errorf("unset ambient %.0f s, want the default's %.0f s and slower at 0 °C (%.0f s)", unset, want, at0)
+ }
+}
+
+// Input that is not a current, a step or a cabinet leaves the phase as it
+// was and says so; an infinite current blows the fuse.
+func TestInvalidInputLeavesTheState(t *testing.T) {
+ p := profile(t, IecGg)
+ warm := ref()
+ hold(p, warm, 1.2, 25, 600, 1)
+ for _, c := range []struct {
+ name string
+ amps, rating, dt float64
+ ambient *float64
+ }{
+ {"NaN current", math.NaN(), 25, 1, nil},
+ {"NaN step", 30, 25, math.NaN(), nil},
+ {"infinite step", 30, 25, math.Inf(1), nil},
+ {"negative step", 30, 25, -1, nil},
+ {"NaN ambient", 30, 25, 1, Ambient(math.NaN())},
+ {"infinite ambient", 30, 25, 1, Ambient(math.Inf(1))},
+ {"ambient at the melting point", 30, 25, 1, Ambient(MeltC)},
+ {"no rating", 30, 0, 1, nil},
+ {"NaN rating", 30, math.NaN(), 1, nil},
+ {"10 A, below the bands", 15, 10, 1, nil},
+ {"80 A, above them", 100, 80, 1, nil},
+ } {
+ ph := *warm
+ if c.ambient != nil {
+ ph.AmbientC = c.ambient
+ }
+ before := ph
+ st := p.Step(&ph, c.amps, c.rating, c.dt)
+ if st.Invalid == "" || ph.Theta != before.Theta || ph.Damage != before.Damage || ph.Blown {
+ t.Errorf("%s: %+v, state %v → %v; want refused, unchanged", c.name, st, before.Theta, ph.Theta)
+ }
+ }
+ for _, amps := range []float64{math.Inf(-1), 1e160 * 25} {
+ ph := *warm
+ st := p.Step(&ph, amps, 25, 1)
+ _, err := json.Marshal(st)
+ if !ph.Blown || st.Damage != 1 || st.Invalid != "" || math.IsInf(st.Ratio, 0) || err != nil {
+ t.Errorf("%g A: %+v (json: %v), want blown with a finite ratio", amps, st, err)
+ }
+ for _, x := range ph.Theta {
+ if math.IsNaN(x) || math.IsInf(x, 0) {
+ t.Errorf("%g A left the state %v", amps, ph.Theta)
+ }
+ }
+ }
+ if p.Supports(10) || p.Supports(80) || !p.Supports(16) || !p.Supports(30) || !p.Supports(63) {
+ t.Error("gG supports 16-63 A only")
+ }
+}
+
+// A warmer cabinet leaves less margin to the melting point.
+func TestWarmCabinetBlowsSooner(t *testing.T) {
+ p := profile(t, IecGg)
+ at20 := p.MeltSeconds(1.5, 25, 20)
+ at40 := p.MeltSeconds(1.5, 25, 40)
+ if at40 >= at20 {
+ t.Errorf("1.5× at 40 °C %.0f s, at 20 °C %.0f s: want sooner when warmer", at40, at20)
+ }
+ // IEC 60269's curves hold approximately at 30 °C as at 20 °C.
+ if at30 := p.MeltSeconds(5, 25, 30); math.Abs(math.Log10(at30/p.MeltSeconds(5, 25, 20))) > 0.1 {
+ t.Errorf("5× at 30 °C %.2f s, more than 0.1 decades from 20 °C", at30)
+ }
+}
+
+// aM keeps its gates at any rating.
+func TestAMGates(t *testing.T) {
+ p := profile(t, IecAm)
+ for _, rating := range []float64{10, 16, 63, 125} {
+ b, ok := p.Band(rating)
+ if !ok {
+ t.Fatalf("aM has no band for %g A", rating)
+ }
+ for _, g := range b.Gates {
+ got := b.MeltSeconds(g.Multiple, RefAmbientC)
+ if g.Min && got < g.Seconds || !g.Min && got > g.Seconds {
+ t.Errorf("aM at %g×: %.3g s, gate %g s (min %v)", g.Multiple, got, g.Seconds, g.Min)
+ }
+ }
+ }
+}
+
+func TestManualNeverBlows(t *testing.T) {
+ p := profile(t, Manual)
+ ph := &Phase{}
+ if hold(p, ph, 20, 25, 600, 1); ph.Blown {
+ t.Fatal("a manual fuse blew")
+ }
+}
+
+// Every profile but manual has a curve and a source and holds at its rating.
+func TestEveryProfileHasACurve(t *testing.T) {
+ for _, p := range All() {
+ if p.ID == Manual {
+ continue
+ }
+ if len(p.Curve) == 0 || len(p.Bands) == 0 || p.Source == "" {
+ t.Errorf("%s: %d curve points, %d bands, source %q", p.ID, len(p.Curve), len(p.Bands), p.Source)
+ }
+ if !math.IsInf(p.MeltSeconds(1, 25, RefAmbientC), 1) {
+ t.Errorf("%s melts at its rating", p.ID)
+ }
+ }
+}
diff --git a/internal/mcp/server.go b/internal/mcp/server.go
index 5bb7491..2e246af 100644
--- a/internal/mcp/server.go
+++ b/internal/mcp/server.go
@@ -689,7 +689,7 @@ func buildToolSpecs() ([]Tool, map[string]ToolSpec) {
"id": schemaInteger("Site ID"),
"name": schemaString("Site name"),
"base_load_w": schemaNumber("Base load in watts"),
- "fuse_size": schemaInteger("Fuse size"),
+ "fuse_size": schemaInteger("Fuse size in A: 16, 20, 25, 32, 35, 40, 50 or 63"),
"fuse_profile": schemaString("Fuse profile"),
"phase_factors": schemaObject(nil, nil),
"fuse_blown": schemaObject(nil, nil),