diff --git a/myps/test-scripts/factory/furnace/furnace2.myps b/myps/test-scripts/factory/furnace/furnace2.myps
new file mode 100644
index 0000000..7b62b30
--- /dev/null
+++ b/myps/test-scripts/factory/furnace/furnace2.myps
@@ -0,0 +1,162 @@
+# NOTE: To start operation, write a "value" greater than zero to the housing. This value is taken to
+# be a code for the target temperature and pressure with the format `PPPTTT`, where tT=TTT*10$ and
+# pT=PPP*100.
+#
+# e.g. 500050 -> tT = 050*10 = 500 K, pT = 500*100 = 50,000 kPa = 50 MPa
+#
+# IMPORTANT: Validation for this input is not done here, and the program will attempt to use
+# whatever it is given. If the target temperature is less-than the cold source temperature or
+# greater than the hot source temperature then wackiness will certainly occurr!
+# (I may change this point in the future)
+furnace    = d0
+tankC      = d1
+tankH      = d2
+nCPump     = d3
+nHPump     = d4
+nIAnalyzer = d5
+
+# Constants
+def vF = 1000             # furnace volume
+def R  = 8.31446261815324 # universal gas consant
+def dt = 0.5              # time step
+
+# (NOTE: ADJUST PER SETUP)
+# Source volume-mole factor
+# e.g. 6100  = 1 furnace (6000) + 2 pipes (200) L
+# (TODO: Need to implement defines calculated from other defines)
+def nCfactor = (6100) / 6000
+def nHfactor = (6100) / 6000
+
+# (NOTE: MAY NEED TO ADJUST)
+# PD (proportional, derivative) term coefficients
+def kPF = 0.01
+def kDF = 0
+def kPC = 0.1
+def kDC = 0
+def kPH = 0.1
+def kDH = 0
+
+# Acceptable error amounts for
+# (a) the removal of gas from the furnace,
+# (b) the transfer of cold CO2 (i.e. removal from source)
+# (c) the transfer of hot CO2
+def ERRORF = 0.25
+def ERRORC = 0.25
+def ERRORH = 0.25
+
+db.Setting = -1
+furnace.SettingInput = 0
+furnace.SettingOutput = 0
+loop:
+    yield()
+    input = db.Setting
+    if input <= 0:
+        db.Setting = -1
+    else:
+        db.Setting = 0
+        fix rF, rC, rH, eFPrev, eCPrev, eHPrev, eTotal0
+        tag calcMolesTargets:
+            # (a) unpack input target temperature and pressure (tT/pT)
+            input = trunc(input)
+            tT = (input % 1000) * 10
+            pT = trunc(input / 1000) * 100
+
+            # (b) read furnace moles and temperature (tF/pF)
+            nF = furnace.TotalMoles
+            tF = furnace.Temperature
+
+            # (c) read tank temperatures (tC/tH)
+            tC = tankC.Temperature
+            tH = tankH.Temperature
+
+            # (d) calculate moles to remove and add from C and H (nR/nCI/nHI)
+            nT = pT * vF / (R * tT)
+            nRC = nT * (tT - tH) / (tH - tF) + nF
+            nRH = nT * (tT - tC) / (tC - tF) + nF
+            nR = max(0, max(nRC, nRH))
+            nI = nT - nF + nR
+            tI = (tT * nT - tF * (nF - nR)) / nI
+            nHI = nI * (tC - tI) / (tC - tH)
+            nCI = nI - nHI
+
+            # (e) calculate target moles for each source (rF/rC/rH),
+            # initial errors (eFPrev/eCPrev/eHPrev), and total initial error (eTotal0)
+            rF = nF - nR
+            eFPrev = nF - rF
+
+            nC = nCfactor * tankC.TotalMoles
+            rC = nC - nCI
+            eCPrev = nC - rC
+
+            nH = nHfactor * tankH.TotalMoles
+            rH = nH - nHI
+            eHPrev = nH - rH
+
+            eTotal0 = eFPrev + eCPrev + eHPrev
+
+        # (f) reach target moles for each source using PD (proportional, derivative) controls
+        # NOTE: currently some pumps are bugged in that their Settings are unbounded outside
+        # their normal range [0,100] when set with logic (i.e. it's possible to set negative values
+        # and larger than 100 values). Those control values are for now temporarily bounded.
+        tag MolesLoop:
+            nF = furnace.TotalMoles
+            nF - nR
+
+            
+
+        tag PDControl:
+            tag PDLoop:
+                yield()
+                # (f.1) Calculate the furnace output control value (uF) and update error (eFPrev)
+                nF = furnace.TotalMoles
+                eF = nF - rF
+                #uF = max(0, min(100, (kPF * eF) + (kDF * (eF - eFPrev) / dt))) # (kPF, kDF)
+                uF = (kPF * eF) + (kDF * (eF - eFPrev) / dt) # (kPF, kDF)
+                eFPrev = eF
+
+                # (f.2) Calculate the cold source control value (uC) and update error (eCPrev)
+                nC = nCfactor * tankC.TotalMoles
+                eC = nC - rC
+                #uC = max(0, min(100, (kPC * eC) + (kDC * (eC - eCPrev) / dt))) # (kPC, kDC)
+                uC = (kPC * eC) + (kDC * (eC - eCPrev) / dt) # (kPC, kDC)
+                eCPrev = eC
+
+                # (f.3) Calculate the hot source control value (uH) and update error (eHPrev)
+                nH = nHfactor * tankH.TotalMoles
+                eH = nH - rH
+                #uH = max(0, min(100, (kPH * eH) + (kDH * (eH - eHPrev) / dt)))  # (kPH, kDH)
+                uH = (kPH * eH) + (kDH * (eH - eHPrev) / dt) # (kPH, kDH)
+                eHPrev = eH
+
+                # (f.4) Set housing to current progress (total error / initial error)
+                eTotal = (eF + eC + eH) / eTotal0
+                db.Setting = 1.05 - eTotal
+
+                # (f.5) Set rates from control values, and loop if any error is above limit
+                FNotDone = (eF > ERRORF)
+                CNotDone = (eC > ERRORC)
+                HNotDone = (eH > ERRORH)
+
+                furnace.SettingOutput = FNotDone ? uF : 0
+                #furnace.SettingInput  = FNotDone ? 0  : 1000
+                nCPump.Setting        = CNotDone ? uC : 0
+                nHPump.Setting        = HNotDone ? uH : 0
+
+                nCPump.On = CNotDone
+                nHPump.On = HNotDone
+
+                bnez(FNotDone, PDLoop)
+                bnez(CNotDone, PDLoop)
+                bnez(HNotDone, PDLoop)
+
+        # (g) input the nC and nH gas
+        # NOTE: As mentioned previously, furnace.SettingInput is currently unbounded when set by
+        # logic. Using this to cheat here and instantly fill the furnace volume from the mixer
+        # system by setting the input rate to be equal to the mixer system volume.
+        furnace.SettingInput  = 1000
+        while nIAnalyzer.TotalMoles > 0:
+            yield()
+        furnace.SettingInput = 0
+
+        # (?) enjoy
+        db.Setting = -1