Add example flow

Reviewed-on: #6

.gitignore

@@ -0,0 +1,136 @@
+# Logs
+logs
+*.log
+npm-debug.log*
+yarn-debug.log*
+yarn-error.log*
+lerna-debug.log*
+.pnpm-debug.log*
+
+# Diagnostic reports (https://nodejs.org/api/report.html)
+report.[0-9]*.[0-9]*.[0-9]*.[0-9]*.json
+
+# Runtime data
+pids
+*.pid
+*.seed
+*.pid.lock
+
+# Directory for instrumented libs generated by jscoverage/JSCover
+lib-cov
+
+# Coverage directory used by tools like istanbul
+coverage
+*.lcov
+
+# nyc test coverage
+.nyc_output
+
+# Grunt intermediate storage (https://gruntjs.com/creating-plugins#storing-task-files)
+.grunt
+
+# Bower dependency directory (https://bower.io/)
+bower_components
+
+# node-waf configuration
+.lock-wscript
+
+# Compiled binary addons (https://nodejs.org/api/addons.html)
+build/Release
+
+# Dependency directories
+node_modules/
+jspm_packages/
+
+# Snowpack dependency directory (https://snowpack.dev/)
+web_modules/
+
+# TypeScript cache
+*.tsbuildinfo
+
+# Optional npm cache directory
+.npm
+
+# Optional eslint cache
+.eslintcache
+
+# Optional stylelint cache
+.stylelintcache
+
+# Microbundle cache
+.rpt2_cache/
+.rts2_cache_cjs/
+.rts2_cache_es/
+.rts2_cache_umd/
+
+# Optional REPL history
+.node_repl_history
+
+# Output of 'npm pack'
+*.tgz
+
+# Yarn Integrity file
+.yarn-integrity
+
+# dotenv environment variable files
+.env
+.env.development.local
+.env.test.local
+.env.production.local
+.env.local
+
+# parcel-bundler cache (https://parceljs.org/)
+.cache
+.parcel-cache
+
+# Next.js build output
+.next
+out
+
+# Nuxt.js build / generate output
+.nuxt
+dist
+
+# Gatsby files
+.cache/
+# Comment in the public line in if your project uses Gatsby and not Next.js
+# https://nextjs.org/blog/next-9-1#public-directory-support
+# public
+
+# vuepress build output
+.vuepress/dist
+
+# vuepress v2.x temp and cache directory
+.temp
+.cache
+
+# vitepress build output
+**/.vitepress/dist
+
+# vitepress cache directory
+**/.vitepress/cache
+
+# Docusaurus cache and generated files
+.docusaurus
+
+# Serverless directories
+.serverless/
+
+# FuseBox cache
+.fusebox/
+
+# DynamoDB Local files
+.dynamodb/
+
+# TernJS port file
+.tern-port
+
+# Stores VSCode versions used for testing VSCode extensions
+.vscode-test
+
+# yarn v2
+.yarn/cache
+.yarn/unplugged
+.yarn/build-state.yml
+.yarn/install-state.gz
+.pnp.*

README.md

@@ -1,16 +1,20 @@
-# reactor
+# Reactor: Advanced Hydraulic Tank & Biological Process Simulator
 
-Reactor: Advanced Hydraulic Tank & Biological Process Simulator
+A comprehensive reactor class for wastewater treatment simulation featuring non-ideal plug flow hydraulics and ASM3 biological modeling.
-
+
-A comprehensive reactor class for wastewater treatment simulation featuring plug flow hydraulics, ASM1-ASM3 biological modeling, and multi-sectional concentration tracking. Implements hydraulic retention time calculations, dispersion modeling, and real-time biological reaction kinetics for accurate process simulation.
+## How to use this Node
-
+### Set Node Properties
-Key Features:
+- Set reactor type: Continuously Stirred Tank Reactor (CSTR) or Plug Flow Reactor (PFR)
-
+- Configure reactor sizing: Set reactor volume [ $m^3$ ] and (for PFRs) set the reactor length [ $m$ ]
-Plug Flow Hydraulics: Multi-section reactor with configurable sectioning factor and dispersion modeling
+- (For PFRs) set reactor spatial resolution: A value of 10 or 20 is good. A higher resolution means more accurate simulation, at higher computational expense. Note that connected reactors must have similar Δx values.
-ASM1 Integration: Complete biological nutrient removal modeling with 13 state variables (COD, nitrogen, phosphorus)
+- Set initial state of reactor: set the intial concentrations of all relevant reaction species.
-Dynamic Volume Control: Automatic section management with overflow handling and retention time calculations
+- (Optional) set $k_L a$ to calculate OTR internally, rather than providing it explicitly, using simple mass transfer model.
-Oxygen Transfer: Saturation-limited O2 transfer with Fick's law slowdown effects and solubility curves
+
-Real-time Kinetics: Continuous biological reaction rate calculations with configurable time acceleration
+### Accepted Node inputs
-Weighted Averaging: Volume-based concentration mixing for accurate mass balance calculations
+- \{ topic: clock, payload: \<timestamp [ $ms$ ]\> \} - **required** clock signal to make reactor update state.
-Child Registration: Integration with diffuser systems and upstream/downstream reactor networks
+- \{ topic: Fluent, payload: \{ F: \<flow rate [ $m^3 d^{-1}$ ]\>, C: \<array with concentrations\> \} \} - sets inflow composition and flow rate.
-Supports complex biological treatment train modeling with temperature compensation, sludge calculations, and comprehensive process monitoring for wastewater treatment plant optimization and regulatory compliance.
+- \{ topic: Dispersion, payload: \<dispersion coefficient in [ $m d^{-2}$ ]\> \} - sets PFR dispersion coefficient.
+- \{ topic: OTR, payload: \<oxygen transfer rate [ $ g d^{-1} m^{-3}$ ]\> \} - sets current oxygen transfer rate.
+
+## Troubleshooting
+Check for possible numerical warnings. These tell you which simulation parameters to change. If solutions appear to be oscillate, try reducing the time step. If solutions appear to be too dispersive, try increasing the reactor resolution.

example_flow/reactor_flows.json

@@ -0,0 +1,838 @@
+[
+    {
+        "id": "a6b85e226d144df1",
+        "type": "tab",
+        "label": "Flow 3",
+        "disabled": false,
+        "info": "",
+        "env": []
+    },
+    {
+        "id": "a94c85d65c71b66a",
+        "type": "inject",
+        "z": "a6b85e226d144df1",
+        "name": "Influx composition 1",
+        "props": [
+            {
+                "p": "payload"
+            },
+            {
+                "p": "topic",
+                "vt": "str"
+            },
+            {
+                "p": "timestamp",
+                "v": "",
+                "vt": "date"
+            }
+        ],
+        "repeat": "864",
+        "crontab": "",
+        "once": true,
+        "onceDelay": "5",
+        "topic": "Fluent",
+        "payload": "{\"inlet\":0,\"F\":6600,\"C\":[0,30,100,16,0,0,5,25,75,30,0,0,125]}",
+        "payloadType": "json",
+        "x": 260,
+        "y": 240,
+        "wires": [
+            [
+                "b726c5be41c24dcb"
+            ]
+        ]
+    },
+    {
+        "id": "0ebfbbf57bba79f1",
+        "type": "inject",
+        "z": "a6b85e226d144df1",
+        "name": "",
+        "props": [
+            {
+                "p": "timestamp",
+                "v": "",
+                "vt": "date"
+            },
+            {
+                "p": "topic",
+                "vt": "str"
+            }
+        ],
+        "repeat": "1",
+        "crontab": "",
+        "once": true,
+        "onceDelay": 0.1,
+        "topic": "clock",
+        "x": 300,
+        "y": 280,
+        "wires": [
+            [
+                "b726c5be41c24dcb"
+            ]
+        ]
+    },
+    {
+        "id": "cd0e4a78d1a59a6e",
+        "type": "inject",
+        "z": "a6b85e226d144df1",
+        "name": "",
+        "props": [
+            {
+                "p": "payload"
+            },
+            {
+                "p": "topic",
+                "vt": "str"
+            }
+        ],
+        "repeat": "",
+        "crontab": "",
+        "once": true,
+        "onceDelay": 0.1,
+        "topic": "Dispersion",
+        "payload": "100",
+        "payloadType": "num",
+        "x": 260,
+        "y": 320,
+        "wires": [
+            [
+                "b726c5be41c24dcb",
+                "b0819ff9a4010227"
+            ]
+        ]
+    },
+    {
+        "id": "b726c5be41c24dcb",
+        "type": "reactor",
+        "z": "a6b85e226d144df1",
+        "name": "Anoxic 1",
+        "reactor_type": "PFR",
+        "volume": "730",
+        "length": "10",
+        "resolution_L": "20",
+        "kla": "",
+        "S_O_init": 0,
+        "S_I_init": 30,
+        "S_S_init": 100,
+        "S_NH_init": 16,
+        "S_N2_init": 0,
+        "S_NO_init": 0,
+        "S_HCO_init": 5,
+        "X_I_init": 25,
+        "X_S_init": 75,
+        "X_H_init": 30,
+        "X_STO_init": 0,
+        "X_A_init": "30",
+        "X_TS_init": "132",
+        "timeStep": "2",
+        "enableLog": true,
+        "logLevel": "info",
+        "positionVsParent": "upstream",
+        "x": 540,
+        "y": 240,
+        "wires": [
+            [
+                "057ab2dcd4739aef"
+            ],
+            [],
+            [
+                "b0819ff9a4010227",
+                "214aff7330e81d3f",
+                "2ec38c4ae9aa6e7e"
+            ]
+        ]
+    },
+    {
+        "id": "b0819ff9a4010227",
+        "type": "reactor",
+        "z": "a6b85e226d144df1",
+        "name": "Aerobic 1",
+        "reactor_type": "PFR",
+        "volume": "1460",
+        "length": "20",
+        "resolution_L": "40",
+        "kla": "2400",
+        "S_O_init": 0,
+        "S_I_init": 30,
+        "S_S_init": 100,
+        "S_NH_init": 16,
+        "S_N2_init": 0,
+        "S_NO_init": 0,
+        "S_HCO_init": 5,
+        "X_I_init": 25,
+        "X_S_init": 75,
+        "X_H_init": "500",
+        "X_STO_init": 0,
+        "X_A_init": "30",
+        "X_TS_init": "132",
+        "timeStep": "2",
+        "enableLog": true,
+        "logLevel": "info",
+        "positionVsParent": "upstream",
+        "x": 940,
+        "y": 240,
+        "wires": [
+            [
+                "e0049d66fefeb3e6"
+            ],
+            [],
+            [
+                "47120bc82aa7bf49"
+            ]
+        ]
+    },
+    {
+        "id": "214aff7330e81d3f",
+        "type": "rotatingMachine",
+        "z": "a6b85e226d144df1",
+        "speed": 1,
+        "startup": 0,
+        "warmup": 0,
+        "shutdown": 0,
+        "cooldown": 0,
+        "machineCurve": {},
+        "flowNumber": "1",
+        "uuid": "",
+        "supplier": "Hidrostal",
+        "category": "Pumps",
+        "assetType": "Centrifugal",
+        "model": "hidrostal-H05K-S03R",
+        "unit": "l/s",
+        "enableLog": true,
+        "logLevel": "info",
+        "positionVsParent": "downstream",
+        "positionIcon": "→",
+        "hasDistance": false,
+        "distance": "",
+        "x": 740,
+        "y": 340,
+        "wires": [
+            [],
+            [],
+            [
+                "b0819ff9a4010227"
+            ]
+        ]
+    },
+    {
+        "id": "0be74d5c77febec1",
+        "type": "measurement",
+        "z": "a6b85e226d144df1",
+        "name": "sensor",
+        "scaling": false,
+        "i_min": 0,
+        "i_max": 0,
+        "i_offset": 0,
+        "o_min": 900,
+        "o_max": 1200,
+        "simulator": true,
+        "smooth_method": "",
+        "count": 10,
+        "uuid": "",
+        "supplier": "Vega",
+        "category": "Sensor",
+        "assetType": "Flow",
+        "model": "VegaFlow 10",
+        "unit": "m³/h",
+        "enableLog": true,
+        "logLevel": "error",
+        "positionVsParent": "atEquipment",
+        "positionIcon": "⊥",
+        "hasDistance": false,
+        "distance": "",
+        "x": 600,
+        "y": 460,
+        "wires": [
+            [],
+            [],
+            [
+                "214aff7330e81d3f"
+            ]
+        ]
+    },
+    {
+        "id": "668f41a05698f21b",
+        "type": "inject",
+        "z": "a6b85e226d144df1",
+        "name": "",
+        "props": [
+            {
+                "p": "payload"
+            },
+            {
+                "p": "topic",
+                "vt": "str"
+            }
+        ],
+        "repeat": "",
+        "crontab": "",
+        "once": true,
+        "onceDelay": 0.1,
+        "topic": "execSequence",
+        "payload": "{\"source\":\"parent\",\"action\":\"execSequence\",\"parameter\":\"startup\"}",
+        "payloadType": "json",
+        "x": 330,
+        "y": 400,
+        "wires": [
+            [
+                "214aff7330e81d3f",
+                "2ec38c4ae9aa6e7e"
+            ]
+        ]
+    },
+    {
+        "id": "47120bc82aa7bf49",
+        "type": "settler",
+        "z": "a6b85e226d144df1",
+        "name": "",
+        "model": "mb-model",
+        "enableLog": true,
+        "logLevel": "info",
+        "positionVsParent": "atEquipment",
+        "x": 1100,
+        "y": 400,
+        "wires": [
+            [
+                "b1af49d1d0eb6783",
+                "6baffb7954cf0cce"
+            ],
+            [],
+            []
+        ]
+    },
+    {
+        "id": "2ec38c4ae9aa6e7e",
+        "type": "rotatingMachine",
+        "z": "a6b85e226d144df1",
+        "speed": 1,
+        "startup": 0,
+        "warmup": 0,
+        "shutdown": 0,
+        "cooldown": 0,
+        "machineCurve": {},
+        "flowNumber": "2",
+        "uuid": "",
+        "supplier": "Hidrostal",
+        "category": "Pumps",
+        "assetType": "Centrifugal",
+        "model": "hidrostal-H05K-S03R",
+        "unit": "l/s",
+        "enableLog": true,
+        "logLevel": "info",
+        "positionVsParent": "downstream",
+        "positionIcon": "→",
+        "hasDistance": false,
+        "distance": "",
+        "x": 860,
+        "y": 460,
+        "wires": [
+            [],
+            [],
+            [
+                "47120bc82aa7bf49"
+            ]
+        ]
+    },
+    {
+        "id": "21eff612371519c5",
+        "type": "measurement",
+        "z": "a6b85e226d144df1",
+        "name": "sensor",
+        "scaling": false,
+        "i_min": 0,
+        "i_max": 0,
+        "i_offset": 0,
+        "o_min": 500,
+        "o_max": 660,
+        "simulator": true,
+        "smooth_method": "",
+        "count": 10,
+        "uuid": "",
+        "supplier": "Vega",
+        "category": "Sensor",
+        "assetType": "Flow",
+        "model": "VegaFlow 10",
+        "unit": "m³/h",
+        "enableLog": true,
+        "logLevel": "error",
+        "positionVsParent": "atEquipment",
+        "positionIcon": "⊥",
+        "hasDistance": false,
+        "distance": "",
+        "x": 720,
+        "y": 580,
+        "wires": [
+            [],
+            [],
+            [
+                "2ec38c4ae9aa6e7e"
+            ]
+        ]
+    },
+    {
+        "id": "057ab2dcd4739aef",
+        "type": "function",
+        "z": "a6b85e226d144df1",
+        "name": "Data_converter",
+        "func": "if (msg.topic != \"Fluent\") {\n    return;\n}\n\nconst [S_O, S_I, S_S, S_NH, S_N2, S_NO, S_HCO, X_I, X_S, X_H, X_STO, X_A, X_TS] = msg.payload.C;\n\nmsg = {payload: [\n    { \"Series\": \"S_O\",      \"Y\": S_O},\n    { \"Series\": \"S_I\",      \"Y\": S_I},\n    { \"Series\": \"S_S\",      \"Y\": S_S},\n    { \"Series\": \"S_NH\",     \"Y\": S_NH},\n    { \"Series\": \"S_N2\",     \"Y\": S_N2},\n    { \"Series\": \"S_NO\",     \"Y\": S_NO},\n    { \"Series\": \"S_HCO\",    \"Y\": S_HCO},\n    { \"Series\": \"X_I\",      \"Y\": X_I},\n    { \"Series\": \"X_S\",      \"Y\": X_S},\n    { \"Series\": \"X_H\",      \"Y\": X_H},\n    { \"Series\": \"X_STO\",    \"Y\": X_STO},\n    { \"Series\": \"X_A\",      \"Y\": X_A},\n    { \"Series\": \"X_TS\",     \"Y\": X_TS}\n    ]};\n\nreturn msg;",
+        "outputs": 1,
+        "timeout": 0,
+        "noerr": 0,
+        "initialize": "",
+        "finalize": "",
+        "libs": [],
+        "x": 720,
+        "y": 180,
+        "wires": [
+            [
+                "02a932317c76721e"
+            ]
+        ]
+    },
+    {
+        "id": "02a932317c76721e",
+        "type": "ui-chart",
+        "z": "a6b85e226d144df1",
+        "group": "ae38454098a37db0",
+        "name": "Anoxic reactor",
+        "label": "Anoxic reactor",
+        "order": 9007199254740991,
+        "chartType": "line",
+        "category": "Series",
+        "categoryType": "property",
+        "xAxisLabel": "",
+        "xAxisProperty": "",
+        "xAxisPropertyType": "timestamp",
+        "xAxisType": "time",
+        "xAxisFormat": "",
+        "xAxisFormatType": "auto",
+        "xmin": "",
+        "xmax": "",
+        "yAxisLabel": "",
+        "yAxisProperty": "Y",
+        "yAxisPropertyType": "property",
+        "ymin": "",
+        "ymax": "",
+        "bins": 10,
+        "action": "append",
+        "stackSeries": false,
+        "pointShape": "circle",
+        "pointRadius": 4,
+        "showLegend": true,
+        "removeOlder": "8",
+        "removeOlderUnit": "3600",
+        "removeOlderPoints": "2000",
+        "colors": [
+            "#0095ff",
+            "#ff0000",
+            "#ff7f0e",
+            "#2ca02c",
+            "#a347e1",
+            "#d62728",
+            "#ff9896",
+            "#9467bd",
+            "#c5b0d5"
+        ],
+        "textColor": [
+            "#666666"
+        ],
+        "textColorDefault": true,
+        "gridColor": [
+            "#e5e5e5"
+        ],
+        "gridColorDefault": true,
+        "width": 6,
+        "height": 8,
+        "className": "",
+        "interpolation": "linear",
+        "x": 920,
+        "y": 180,
+        "wires": [
+            []
+        ]
+    },
+    {
+        "id": "e0049d66fefeb3e6",
+        "type": "function",
+        "z": "a6b85e226d144df1",
+        "name": "Data_converter",
+        "func": "if (msg.topic != \"Fluent\") {\n    return;\n}\n\nconst [S_O, S_I, S_S, S_NH, S_N2, S_NO, S_HCO, X_I, X_S, X_H, X_STO, X_A, X_TS] = msg.payload.C;\n\nmsg = {payload: [\n    { \"Series\": \"S_O\",      \"Y\": S_O},\n    { \"Series\": \"S_I\",      \"Y\": S_I},\n    { \"Series\": \"S_S\",      \"Y\": S_S},\n    { \"Series\": \"S_NH\",     \"Y\": S_NH},\n    { \"Series\": \"S_N2\",     \"Y\": S_N2},\n    { \"Series\": \"S_NO\",     \"Y\": S_NO},\n    { \"Series\": \"S_HCO\",    \"Y\": S_HCO},\n    { \"Series\": \"X_I\",      \"Y\": X_I},\n    { \"Series\": \"X_S\",      \"Y\": X_S},\n    { \"Series\": \"X_H\",      \"Y\": X_H},\n    { \"Series\": \"X_STO\",    \"Y\": X_STO},\n    { \"Series\": \"X_A\",      \"Y\": X_A},\n    { \"Series\": \"X_TS\",     \"Y\": X_TS}\n    ]};\n\nreturn msg;",
+        "outputs": 1,
+        "timeout": 0,
+        "noerr": 0,
+        "initialize": "",
+        "finalize": "",
+        "libs": [],
+        "x": 1140,
+        "y": 220,
+        "wires": [
+            [
+                "70eeb4c8caa2db77"
+            ]
+        ]
+    },
+    {
+        "id": "70eeb4c8caa2db77",
+        "type": "ui-chart",
+        "z": "a6b85e226d144df1",
+        "group": "ae38454098a37db0",
+        "name": "Aerobic reactor",
+        "label": "Aerobic reactor / recirculation",
+        "order": 9007199254740991,
+        "chartType": "line",
+        "category": "Series",
+        "categoryType": "property",
+        "xAxisLabel": "",
+        "xAxisProperty": "",
+        "xAxisPropertyType": "timestamp",
+        "xAxisType": "time",
+        "xAxisFormat": "",
+        "xAxisFormatType": "auto",
+        "xmin": "",
+        "xmax": "",
+        "yAxisLabel": "",
+        "yAxisProperty": "Y",
+        "yAxisPropertyType": "property",
+        "ymin": "",
+        "ymax": "",
+        "bins": 10,
+        "action": "append",
+        "stackSeries": false,
+        "pointShape": "circle",
+        "pointRadius": 4,
+        "showLegend": true,
+        "removeOlder": "8",
+        "removeOlderUnit": "3600",
+        "removeOlderPoints": "2000",
+        "colors": [
+            "#0095ff",
+            "#ff0000",
+            "#ff7f0e",
+            "#2ca02c",
+            "#a347e1",
+            "#d62728",
+            "#ff9896",
+            "#9467bd",
+            "#c5b0d5"
+        ],
+        "textColor": [
+            "#666666"
+        ],
+        "textColorDefault": true,
+        "gridColor": [
+            "#e5e5e5"
+        ],
+        "gridColorDefault": true,
+        "width": 6,
+        "height": 8,
+        "className": "",
+        "interpolation": "linear",
+        "x": 1340,
+        "y": 220,
+        "wires": [
+            []
+        ]
+    },
+    {
+        "id": "b1af49d1d0eb6783",
+        "type": "function",
+        "z": "a6b85e226d144df1",
+        "name": "Data_converter",
+        "func": "if (msg.topic != \"Fluent\" || msg.payload.inlet == 0) {\n    return;\n}\n\nconst [S_O, S_I, S_S, S_NH, S_N2, S_NO, S_HCO, X_I, X_S, X_H, X_STO, X_A, X_TS] = msg.payload.C;\n\nmsg = {payload: [\n    { \"Series\": \"X_I\",      \"Y\": X_I},\n    { \"Series\": \"X_S\",      \"Y\": X_S},\n    { \"Series\": \"X_H\",      \"Y\": X_H},\n    { \"Series\": \"X_STO\",    \"Y\": X_STO},\n    { \"Series\": \"X_A\",      \"Y\": X_A},\n    { \"Series\": \"X_TS\",     \"Y\": X_TS}\n    ]};\n\nreturn msg;",
+        "outputs": 1,
+        "timeout": 0,
+        "noerr": 0,
+        "initialize": "",
+        "finalize": "",
+        "libs": [],
+        "x": 1320,
+        "y": 380,
+        "wires": [
+            [
+                "0a3ffa69046a4844"
+            ]
+        ]
+    },
+    {
+        "id": "0a3ffa69046a4844",
+        "type": "ui-chart",
+        "z": "a6b85e226d144df1",
+        "group": "de8b029d69f26c0e",
+        "name": "Sludge composition",
+        "label": "Sludge composition",
+        "order": 9007199254740991,
+        "chartType": "line",
+        "category": "Series",
+        "categoryType": "property",
+        "xAxisLabel": "",
+        "xAxisProperty": "",
+        "xAxisPropertyType": "timestamp",
+        "xAxisType": "time",
+        "xAxisFormat": "",
+        "xAxisFormatType": "auto",
+        "xmin": "",
+        "xmax": "",
+        "yAxisLabel": "",
+        "yAxisProperty": "Y",
+        "yAxisPropertyType": "property",
+        "ymin": "",
+        "ymax": "",
+        "bins": 10,
+        "action": "append",
+        "stackSeries": false,
+        "pointShape": "circle",
+        "pointRadius": 4,
+        "showLegend": true,
+        "removeOlder": "8",
+        "removeOlderUnit": "3600",
+        "removeOlderPoints": "2000",
+        "colors": [
+            "#0095ff",
+            "#ff0000",
+            "#ff7f0e",
+            "#2ca02c",
+            "#a347e1",
+            "#d62728",
+            "#ff9896",
+            "#9467bd",
+            "#c5b0d5"
+        ],
+        "textColor": [
+            "#666666"
+        ],
+        "textColorDefault": true,
+        "gridColor": [
+            "#e5e5e5"
+        ],
+        "gridColorDefault": true,
+        "width": 6,
+        "height": 8,
+        "className": "",
+        "interpolation": "linear",
+        "x": 1530,
+        "y": 380,
+        "wires": [
+            []
+        ]
+    },
+    {
+        "id": "272b4b0050479d13",
+        "type": "measurement",
+        "z": "a6b85e226d144df1",
+        "name": "sensor",
+        "scaling": false,
+        "i_min": 0,
+        "i_max": 0,
+        "i_offset": 0,
+        "o_min": 3200,
+        "o_max": 4000,
+        "simulator": true,
+        "smooth_method": "",
+        "count": 10,
+        "uuid": "",
+        "supplier": "Vega",
+        "category": "Sensor",
+        "assetType": "Quantity (TSS)",
+        "model": "VegaSolidsProbe",
+        "unit": "g/m³",
+        "enableLog": true,
+        "logLevel": "error",
+        "positionVsParent": "atEquipment",
+        "positionIcon": "⊥",
+        "hasDistance": false,
+        "distance": "",
+        "x": 930,
+        "y": 580,
+        "wires": [
+            [],
+            [],
+            [
+                "47120bc82aa7bf49"
+            ]
+        ]
+    },
+    {
+        "id": "6baffb7954cf0cce",
+        "type": "function",
+        "z": "a6b85e226d144df1",
+        "name": "Data_converter",
+        "func": "if (msg.topic != \"Fluent\" || msg.payload.inlet == 1 || msg.payload.inlet == 2) {\n    return;\n}\n\nconst [S_O, S_I, S_S, S_NH, S_N2, S_NO, S_HCO, X_I, X_S, X_H, X_STO, X_A, X_TS] = msg.payload.C;\n\nmsg = {payload: [\n    { \"Series\": \"S_O\",      \"Y\": S_O},\n    { \"Series\": \"S_I\",      \"Y\": S_I},\n    { \"Series\": \"S_S\",      \"Y\": S_S},\n    { \"Series\": \"S_NH\",     \"Y\": S_NH},\n    { \"Series\": \"S_N2\",     \"Y\": S_N2},\n    { \"Series\": \"S_NO\",     \"Y\": S_NO},\n    { \"Series\": \"S_HCO\",    \"Y\": S_HCO}\n    ]};\n\nreturn msg;",
+        "outputs": 1,
+        "timeout": 0,
+        "noerr": 0,
+        "initialize": "",
+        "finalize": "",
+        "libs": [],
+        "x": 1320,
+        "y": 420,
+        "wires": [
+            [
+                "df288814a2a9a2b1"
+            ]
+        ]
+    },
+    {
+        "id": "df288814a2a9a2b1",
+        "type": "ui-chart",
+        "z": "a6b85e226d144df1",
+        "group": "de8b029d69f26c0e",
+        "name": "Effluent",
+        "label": "Effluent",
+        "order": 9007199254740991,
+        "chartType": "line",
+        "category": "Series",
+        "categoryType": "property",
+        "xAxisLabel": "",
+        "xAxisProperty": "",
+        "xAxisPropertyType": "timestamp",
+        "xAxisType": "time",
+        "xAxisFormat": "",
+        "xAxisFormatType": "auto",
+        "xmin": "",
+        "xmax": "",
+        "yAxisLabel": "",
+        "yAxisProperty": "Y",
+        "yAxisPropertyType": "property",
+        "ymin": "",
+        "ymax": "",
+        "bins": 10,
+        "action": "append",
+        "stackSeries": false,
+        "pointShape": "circle",
+        "pointRadius": 4,
+        "showLegend": true,
+        "removeOlder": "8",
+        "removeOlderUnit": "3600",
+        "removeOlderPoints": "2000",
+        "colors": [
+            "#0095ff",
+            "#ff0000",
+            "#ff7f0e",
+            "#2ca02c",
+            "#a347e1",
+            "#d62728",
+            "#ff9896",
+            "#9467bd",
+            "#c5b0d5"
+        ],
+        "textColor": [
+            "#666666"
+        ],
+        "textColorDefault": true,
+        "gridColor": [
+            "#e5e5e5"
+        ],
+        "gridColorDefault": true,
+        "width": 6,
+        "height": 8,
+        "className": "",
+        "interpolation": "linear",
+        "x": 1520,
+        "y": 420,
+        "wires": [
+            []
+        ]
+    },
+    {
+        "id": "ae38454098a37db0",
+        "type": "ui-group",
+        "name": "Group 3",
+        "page": "ca564642bfc5606c",
+        "width": 6,
+        "height": 1,
+        "order": -1,
+        "showTitle": true,
+        "className": "",
+        "visible": "true",
+        "disabled": "false",
+        "groupType": "default"
+    },
+    {
+        "id": "de8b029d69f26c0e",
+        "type": "ui-group",
+        "name": "Group 4",
+        "page": "ca564642bfc5606c",
+        "width": 6,
+        "height": 1,
+        "order": -1,
+        "showTitle": true,
+        "className": "",
+        "visible": "true",
+        "disabled": "false",
+        "groupType": "default"
+    },
+    {
+        "id": "ca564642bfc5606c",
+        "type": "ui-page",
+        "name": "PFR",
+        "ui": "90eb5f47d95b4087",
+        "path": "/page2",
+        "icon": "home",
+        "layout": "grid",
+        "theme": "2c8bcaa0046b4323",
+        "breakpoints": [
+            {
+                "name": "Default",
+                "px": "0",
+                "cols": "3"
+            },
+            {
+                "name": "Tablet",
+                "px": "576",
+                "cols": "6"
+            },
+            {
+                "name": "Small Desktop",
+                "px": "768",
+                "cols": "9"
+            },
+            {
+                "name": "Desktop",
+                "px": "1024",
+                "cols": "12"
+            }
+        ],
+        "order": -1,
+        "className": "",
+        "visible": "true",
+        "disabled": "false"
+    },
+    {
+        "id": "90eb5f47d95b4087",
+        "type": "ui-base",
+        "name": "Dashboard",
+        "path": "/dashboard",
+        "appIcon": "",
+        "includeClientData": true,
+        "acceptsClientConfig": [
+            "ui-notification",
+            "ui-control"
+        ],
+        "showPathInSidebar": false,
+        "headerContent": "page",
+        "navigationStyle": "default",
+        "titleBarStyle": "default",
+        "showReconnectNotification": true,
+        "notificationDisplayTime": 1,
+        "showDisconnectNotification": true,
+        "allowInstall": true
+    },
+    {
+        "id": "2c8bcaa0046b4323",
+        "type": "ui-theme",
+        "name": "Default",
+        "colors": {
+            "surface": "#ffffff",
+            "primary": "#0094ce",
+            "bgPage": "#eeeeee",
+            "groupBg": "#ffffff",
+            "groupOutline": "#cccccc"
+        },
+        "sizes": {
+            "density": "default",
+            "pagePadding": "12px",
+            "groupGap": "12px",
+            "groupBorderRadius": "4px",
+            "widgetGap": "12px"
+        }
+    }
+]

package.json

@@ -0,0 +1,32 @@
+{
+  "name": "reactor",
+  "version": "0.0.1",
+  "description": "Implementation of the asm3 model for Node-Red",
+  "repository": {
+    "type": "git",
+    "url": "https://gitea.centraal.wbd-rd.nl/RnD/reactor.git"
+  },
+  "keywords": [
+    "asm3",
+    "activated sludge",
+    "wastewater",
+    "biological model",
+    "EVOLV",
+    "node-red"
+  ],
+  "license": "SEE LICENSE",
+  "author": "P.R. van der Wilt",
+  "main": "reactor.js",
+  "scripts": {
+    "test": "node reactor.js"
+  },
+  "node-red": {
+    "nodes": {
+      "reactor": "reactor.js"
+    }
+  },
+  "dependencies": {
+    "generalFunctions": "git+https://gitea.centraal.wbd-rd.nl/RnD/generalFunctions.git",
+    "mathjs": "^14.5.2"
+  }
+}

reactor.html

@@ -0,0 +1,239 @@
+<script src="/reactor/menu.js"></script>
+
+<script type="text/javascript">
+    RED.nodes.registerType("reactor", {
+        category: "EVOLV",
+        color: "#c4cce0",
+        defaults: {
+            name: { value: "" },
+            reactor_type: { value: "CSTR", required: true },
+            volume: { value: 0., required: true },
+            length: { value: 0.},
+            resolution_L: { value: 0.},
+            kla: { value: null },
+            
+            S_O_init: { value: 0., required: true },
+            S_I_init: { value: 30., required: true },
+            S_S_init: { value: 100., required: true },
+            S_NH_init: { value: 16., required: true },
+            S_N2_init: { value: 0., required: true },
+            S_NO_init: { value: 0., required: true },
+            S_HCO_init: { value: 5., required: true },
+            X_I_init: { value: 25., required: true },
+            X_S_init: { value: 75., required: true },
+            X_H_init: { value: 30., required: true },
+            X_STO_init: { value: 0., required: true },
+            X_A_init: { value: 0.001, required: true },
+            X_TS_init: { value: 125.0009, required: true },
+
+            timeStep: { value: 1, required: true },
+
+            enableLog: { value: false },
+            logLevel: { value: "error" },
+
+            positionVsParent: { value: "" }
+        },
+        inputs: 1,
+        outputs: 3,
+        inputLabels: ["input"],
+        outputLabels: ["process", "dbase", "parent"],
+        icon: "font-awesome/fa-recycle",
+        label: function() {
+            return this.name || "Reactor";
+        },
+        oneditprepare: function() {
+            // wait for the menu scripts to load
+            const waitForMenuData = () => {
+                if (window.EVOLV?.nodes?.reactor?.initEditor) {
+                    window.EVOLV.nodes.reactor.initEditor(this);
+                } else {
+                    setTimeout(waitForMenuData, 50);
+                }
+            };
+            waitForMenuData();
+
+            $("#node-input-volume").typedInput({
+                type:"num",
+                types:["num"]
+            });
+            $("#node-input-length").typedInput({
+                type:"num",
+                types:["num"]
+            });
+            $("#node-input-resolution_L").typedInput({
+                type:"num",
+                types:["num"]
+            });
+            $("#node-input-kla").typedInput({
+                type:"num",
+                types:["num"]
+            });
+            $(".concentrations").typedInput({
+                type:"num",
+                types:["num"]
+            });
+            $("#node-input-reactor_type").typedInput({
+                types: [
+                    {
+                        value: "CSTR",
+                        options: [
+                            { value: "CSTR", label: "CSTR"},
+                            { value: "PFR", label: "PFR"}
+                        ]
+                    }
+                ]
+            })
+            $("#node-input-reactor_type").on("change", function() {
+                const type = $("#node-input-reactor_type").typedInput("value");
+                if (type === "CSTR") {
+                    $(".PFR").hide();
+                } else {
+                    $(".PFR").show();
+                }
+            });
+            $("#node-input-timeStep").typedInput({
+                type:"num",
+                types:["num"]
+            })
+            // Set initial visibility on dialog open
+            const initialType = $("#node-input-reactor_type").typedInput("value");
+            if (initialType === "CSTR") {
+                $(".PFR").hide();
+            } else {
+                $(".PFR").show();
+            }
+
+            const updateDx = () => {
+                const length = parseFloat($("#node-input-length").val()) || 0;
+                const resolution = parseFloat($("#node-input-resolution_L").val()) || 1;
+                const dx = resolution > 0 ? (length / resolution).toFixed(6) : "N/A";
+                $("#dx-output").text(dx + " m");
+            };
+
+            // Set up event listeners for real-time updates
+            $("#node-input-length, #node-input-resolution_L").on("change keyup", updateDx);
+            
+            // Initial calculation
+            updateDx();
+        },
+        oneditsave: function() {
+            // save logger fields
+            if (window.EVOLV?.nodes?.reactor?.loggerMenu?.saveEditor) {
+                window.EVOLV.nodes.reactor.loggerMenu.saveEditor(this);
+            }
+
+            // save position field
+            if (window.EVOLV?.nodes?.measurement?.positionMenu?.saveEditor) {
+                window.EVOLV.nodes.rotatingMachine.positionMenu.saveEditor(this);
+            }
+
+            let volume = parseFloat($("#node-input-volume").typedInput("value"));
+            if (isNaN(volume) || volume <= 0) {
+                RED.notify("Fluid volume not set correctly", {type: "error"});
+            }
+        }
+    });
+</script>
+
+<script type="text/html" data-template-name="reactor">
+    <div class="form-row">
+        <label for="node-input-name"><i class="fa fa-tag"></i> Name</label>
+        <input type="text" id="node-input-name" placeholder="Name">
+    </div>
+    <h2> Reactor properties </h2>
+    <div class="form-row">
+        <label for="node-input-reactor_type"><i class="fa fa-tag"></i> Reactor type</label>
+        <input type="text" id="node-input-reactor_type">
+    </div>
+    <div class="form-row">
+        <label for="node-input-volume"><i class="fa fa-tag"></i> Fluid volume [m3]</label>
+        <input type="text" id="node-input-volume" placeholder="m3">
+    </div>
+    <div class="form-row PFR">
+        <label for="node-input-length"><i class="fa fa-tag"></i> Reactor length [m]</label>
+        <input type="text" id="node-input-length" placeholder="m">
+    </div>
+    <h2> Simulation parameters </h2>
+    <div class="form-row">
+        <label for="node-input-timeStep"><i class="fa fa-tag"></i> Time step [s]</label>
+        <input type="text" id="node-input-timeStep" placeholder="s">
+    </div>
+    <div class="form-row PFR">
+        <label for="node-input-resolution_L"><i class="fa fa-tag"></i> Spatial resolution</label>
+        <input type="text" id="node-input-resolution_L" placeholder="#">
+    </div>
+    <div class="form-row PFR">
+        <label for="node-input-dx"><i class="fa fa-tag"></i> Δx (length / resolution) [m]</label>
+        <span id="dx-output" style="display: inline-block; padding: 8px; font-weight: bold;">--</span>
+    </div>
+    <h3> Internal mass transfer calculation (optional) </h3>
+    <div class="form-row">
+        <label for="node-input-kla"><i class="fa fa-tag"></i> kLa [d-1]</label>
+        <input type="text" id="node-input-kla" placeholder="d-1">
+    </div>
+    <h2> Dissolved components </h2>
+    <div class="form-row">
+        <label for="node-input-S_O_init"><i class="fa fa-tag"></i> Initial dissolved oxygen [g O2 m-3]</label>
+        <input type="text" id="node-input-S_O_init" class="concentrations">
+    </div>
+    <div class="form-row">
+        <label for="node-input-S_I_init"><i class="fa fa-tag"></i> Initial soluble inert  organics [g COD m-3]</label>
+        <input type="text" id="node-input-S_I_init" class="concentrations">
+    </div>
+    <div class="form-row">
+        <label for="node-input-S_S_init"><i class="fa fa-tag"></i> Initial readily biodegrable substrates [g COD m-3]</label>
+        <input type="text" id="node-input-S_S_init" class="concentrations">
+    </div>
+    <div class="form-row">
+        <label for="node-input-S_NH_init"><i class="fa fa-tag"></i> Initial ammonium / ammonia [g N m-3]</label>
+        <input type="text" id="node-input-S_NH_init" class="concentrations">
+    </div>
+    <div class="form-row">
+        <label for="node-input-S_N2_init"><i class="fa fa-tag"></i> Initial dinitrogen, released by denitrification [g N m-3]</label>
+        <input type="text" id="node-input-S_N2_init" class="concentrations">
+    </div>
+    <div class="form-row">
+        <label for="node-input-S_NO_init"><i class="fa fa-tag"></i> Initial nitrite + nitrate [g N m-3]</label>
+        <input type="text" id="node-input-S_NO_init" class="concentrations">
+    </div>
+    <div class="form-row">
+        <label for="node-input-S_HCO_init"><i class="fa fa-tag"></i> Initial alkalinity, bicarbonate [mole HCO3- m-3]</label>
+        <input type="text" id="node-input-S_HCO_init" class="concentrations">
+    </div>
+    <h2> Particulate components </h2>
+    <div class="form-row">
+        <label for="node-input-X_I_init"><i class="fa fa-tag"></i> Initial inert particulate organics [g COD m-3]</label>
+        <input type="text" id="node-input-X_I_init" class="concentrations">
+    </div>
+    <div class="form-row">
+        <label for="node-input-X_S_init"><i class="fa fa-tag"></i> Initial slowly biodegrable substrates [g COD m-3]</label>
+        <input type="text" id="node-input-X_S_init" class="concentrations">
+    </div>
+    <div class="form-row">
+        <label for="node-input-X_H_init"><i class="fa fa-tag"></i> Initial heterotrophic biomass [g COD m-3]</label>
+        <input type="text" id="node-input-X_H_init" class="concentrations">
+    </div>
+    <div class="form-row">
+        <label for="node-input-X_STO_init"><i class="fa fa-tag"></i> Initial Organics stored by heterotrophs [g COD m-3]</label>
+        <input type="text" id="node-input-X_STO_init" class="concentrations">
+    </div>
+    <div class="form-row">
+        <label for="node-input-X_A_init"><i class="fa fa-tag"></i> Initial autotrophic, nitrifying biomass [g COD m-3]</label>
+        <input type="text" id="node-input-X_A_init" class="concentrations">
+    </div>
+    <div class="form-row">
+        <label for="node-input-X_TS_init"><i class="fa fa-tag"></i> Initial total suspended solids [g TSS m-3]</label>
+        <input type="text" id="node-input-X_TS_init" class="concentrations">
+    </div>
+
+    <!-- Logger fields injected here -->
+	<div id="logger-fields-placeholder"></div>
+
+    <!-- Position fields will be injected here -->
+    <div id="position-fields-placeholder"></div>
+
+</script>
+
+<script type="text/html" data-help-name="reactor">
+    <p>New reactor node</p>
+</script>

reactor.js

@@ -0,0 +1,26 @@
+const nameOfNode = "reactor"; // name of the node, should match file name and node type in Node-RED
+const nodeClass = require('./src/nodeClass.js'); // node class
+const { MenuManager } = require('generalFunctions');
+
+
+module.exports = function (RED) {
+    // Register the node type
+    RED.nodes.registerType(nameOfNode, function (config) {
+        // Initialize the Node-RED node first
+        RED.nodes.createNode(this, config);
+        // Then create your custom class and attach it
+        this.nodeClass = new nodeClass(config, RED, this, nameOfNode);
+    });
+
+    const menuMgr = new MenuManager();
+
+    // Serve /advancedReactor/menu.js
+    RED.httpAdmin.get(`/${nameOfNode}/menu.js`, (req, res) => {
+        try {
+            const script = menuMgr.createEndpoint(nameOfNode, ['logger', 'position']);
+            res.type('application/javascript').send(script);
+        } catch (err) {
+            res.status(500).send(`// Error generating menu: ${err.message}`);
+        }
+    });
+};

src/nodeClass.js

@@ -0,0 +1,159 @@
+const { Reactor_CSTR, Reactor_PFR } = require('./specificClass.js');
+
+
+class nodeClass {
+    /**
+    * Node-RED node class for advanced-reactor.
+    * @param {object} uiConfig - Node-RED node configuration
+    * @param {object} RED - Node-RED runtime API
+    * @param {object} nodeInstance - Node-RED node instance
+    * @param {string} nameOfNode - Name of the node
+    */
+    constructor(uiConfig, RED, nodeInstance, nameOfNode) {
+        // Preserve RED reference for HTTP endpoints if needed
+        this.node = nodeInstance;
+        this.RED = RED;
+        this.name = nameOfNode;
+        this.source = null;
+
+        this._loadConfig(uiConfig)
+        this._setupClass();
+
+        this._attachInputHandler();
+        this._registerChild();
+        this._startTickLoop();
+        this._attachCloseHandler();
+    }
+
+    /**
+     * Handle node-red input messages
+     */
+    _attachInputHandler() {
+        this.node.on('input', (msg, send, done) => {
+
+            switch (msg.topic) {
+                case "clock":
+                    this.source.updateState(msg.timestamp);
+                    break;
+                case "Fluent":
+                    this.source.setInfluent = msg;
+                    break;
+                case "OTR":
+                    this.source.setOTR = msg;
+                    break;
+                case "Dispersion":
+                    this.source.setDispersion = msg;
+                    break;
+                case 'registerChild':
+                    // Register this node as a parent of the child node
+                    const childId = msg.payload;
+                    const childObj = this.RED.nodes.getNode(childId); 
+                    this.source.childRegistrationUtils.registerChild(childObj.source, msg.positionVsParent);
+                    break;
+                default:
+                    console.log("Unknown topic: " + msg.topic);
+            }
+
+            if (done) {
+                done();
+            }
+        });
+    }
+
+    /**
+     * Parse node configuration
+     * @param {object} uiConfig Config set in UI in node-red
+     */
+    _loadConfig(uiConfig) {
+        this.config = {
+            general: {
+                name: uiConfig.name || this.name,
+                id: this.node.id,
+                unit: null,
+                logging: {
+                    enabled: uiConfig.enableLog,
+                    logLevel: uiConfig.logLevel
+                }
+            },
+            functionality: {
+                positionVsParent: uiConfig.positionVsParent || 'atEquipment', // Default to 'atEquipment' if not specified
+                softwareType: "reactor" // should be set in config manager
+            },
+            reactor_type: uiConfig.reactor_type,
+            volume: parseFloat(uiConfig.volume),
+            length: parseFloat(uiConfig.length),
+            resolution_L: parseInt(uiConfig.resolution_L),
+            kla: parseFloat(uiConfig.kla),
+            initialState: [
+                parseFloat(uiConfig.S_O_init),
+                parseFloat(uiConfig.S_I_init),
+                parseFloat(uiConfig.S_S_init),
+                parseFloat(uiConfig.S_NH_init),
+                parseFloat(uiConfig.S_N2_init),
+                parseFloat(uiConfig.S_NO_init),
+                parseFloat(uiConfig.S_HCO_init),
+                parseFloat(uiConfig.X_I_init),
+                parseFloat(uiConfig.X_S_init),
+                parseFloat(uiConfig.X_H_init),
+                parseFloat(uiConfig.X_STO_init),
+                parseFloat(uiConfig.X_A_init),
+                parseFloat(uiConfig.X_TS_init)
+            ],
+            timeStep: parseFloat(uiConfig.timeStep)
+        }
+    }
+
+    /**
+     * Register this node as a child upstream and downstream.
+     * Delayed to avoid Node-RED startup race conditions.
+     */
+    _registerChild() {
+        setTimeout(() => {
+        this.node.send([
+            null,
+            null,
+            { topic: 'registerChild', payload: this.node.id, positionVsParent: this.config?.functionality?.positionVsParent || 'atEquipment' }
+        ]);
+        }, 100);
+    }
+
+    /**
+     * Setup reactor class based on config
+     */
+    _setupClass() {
+        let new_reactor;
+
+        switch (this.config.reactor_type) {
+            case "CSTR":
+                new_reactor = new Reactor_CSTR(this.config);
+                break;
+            case "PFR":
+                new_reactor = new Reactor_PFR(this.config);
+                break;
+            default:
+                console.warn("Unknown reactor type: " + uiConfig.reactor_type);
+        }
+
+        this.source = new_reactor; // protect from reassignment
+        this.node.source = this.source;
+    }
+
+    _startTickLoop() {
+        setTimeout(() => {
+            this._tickInterval = setInterval(() => this._tick(), 1000);
+        }, 1000);
+    }
+
+    _tick(){
+        this.node.send([this.source.getEffluent, null, null]);
+    }
+
+    _attachCloseHandler() {
+        this.node.on('close', (done) => {
+            clearInterval(this._tickInterval);
+            done();
+        });
+    }
+}
+
+module.exports = nodeClass;

src/reaction_modules/asm3_class Koch.js

@@ -0,0 +1,222 @@
+const math = require('mathjs');
+
+const ASM_CONSTANTS = {
+  S_O_INDEX: 0,
+  S_NH_INDEX: 3,
+  S_NO_INDEX: 5,
+  NUM_SPECIES: 13
+};
+
+const KINETIC_CONSTANTS = {
+  // Hydrolysis
+  k_H:      9.,   // hydrolysis rate constant       [g X_S g-1 X_H d-1]
+  K_X:      1.,   // hydrolysis saturation constant [g X_S g-1 X_H] 
+  // Heterotrophs
+  k_STO:    12.,  // storage rate constant          [g S_S g-1 X_H d-1]
+  nu_NO:    0.5,  // anoxic reduction factor        [-]
+  K_O:      0.2,  // saturation constant S_0        [g O2 m-3]
+  K_NO:     0.5,  // saturation constant S_NO       [g NO3-N m-3]
+  K_S:      10.,  // saturation constant S_s        [g COD m-3]
+  K_STO:    0.1,  // saturation constant X_STO      [g X_STO g-1 X_H]
+  mu_H_max: 3.,   // maximum specific growth rate   [d-1]
+  K_NH:     0.01, // saturation constant S_NH3      [g NH3-N m-3]
+  K_HCO:    0.1,  // saturation constant S_HCO      [mole HCO3 m-3]
+  b_H_O:    0.3,  // aerobic respiration rate       [d-1]
+  b_H_NO:   0.15, // anoxic respiration rate        [d-1]
+  b_STO_O:  0.3,  // aerobic respitation rate X_STO [d-1]
+  b_STO_NO: 0.15, // anoxic respitation rate X_STO  [d-1]
+  // Autotrophs
+  mu_A_max: 1.3,  // maximum specific growth rate   [d-1]
+  K_A_NH:   1.4,  // saturation constant S_NH3      [g NH3-N m-3] 
+  K_A_O:    0.5,  // saturation constant S_0        [g O2 m-3]
+  K_A_HCO:  0.5,  // saturation constant S_HCO      [mole HCO3 m-3]
+  b_A_O:    0.20, // aerobic respiration rate       [d-1]
+  b_A_NO:   0.10  // anoxic respiration rate        [d-1]
+};
+
+const STOICHIOMETRIC_CONSTANTS = {
+  // Fractions
+  f_SI:     0.,   // fraction S_I from hydrolysis   [g S_I g-1 X_S]
+  f_XI:     0.2,  // fraction X_I from decomp X_H   [g X_I g-1 X_H]
+  // Yields
+  Y_STO_O:  0.80, // aerobic yield X_STO per S_S    [g X_STO g-1 S_S]
+  Y_STO_NO: 0.70, // anoxic yield X_STO per S_S     [g X_STO g-1 S_S]
+  Y_H_O:    0.80, // aerobic yield X_H per X_STO    [g X_H g-1 X_STO]
+  Y_H_NO:   0.65, // anoxic yield X_H per X_STO     [g X_H g-1 X_STO]
+  Y_A:      0.24, // anoxic yield X_A per S_NO      [g X_A g-1 NO3-N]
+  // Composition (COD via DoR)
+  i_CODN:  -1.71, // COD content (DoR)              [g COD g-1 N2-N]
+  i_CODNO: -4.57, // COD content (DoR)              [g COD g-1 NO3-N]
+  // Composition (nitrogen)
+  i_NSI:    0.01, // nitrogen content S_I           [g N g-1 S_I]
+  i_NSS:    0.03, // nitrogen content S_S           [g N g-1 S_S]
+  i_NXI:    0.04, // nitrogen content X_I           [g N g-1 X_I]
+  i_NXS:    0.03, // nitrogen content X_S           [g N g-1 X_S]
+  i_NBM:    0.07, // nitrogen content X_H / X_A     [g N g-1 X_H / X_A]
+  // Composition (TSS)
+  i_TSXI:   0.75, // TSS content X_I                [g TS g-1 X_I]
+  i_TSXS:   0.75, // TSS content X_S                [g TS g-1 X_S]
+  i_TSBM:   0.90, // TSS content X_H / X_A          [g TS g-1 X_H / X_A]
+  i_TSSTO:  0.60, // TSS content X_STO (PHB based)  [g TS g-1 X_STO]
+  // Composition (charge)
+  i_cNH:    1/14, // charge per S_NH                [mole H+ g-1 NH3-N]
+  i_cNO:   -1/14  // charge per S_NO                [mole H+ g-1 NO3-N]
+};
+
+/**
+ * ASM3 class for the Activated Sludge Model No. 3 (ASM3). Using Koch et al. 2000 parameters.
+ */
+class ASM3 {
+
+  constructor() {
+    /**
+     * Kinetic parameters for ASM3 at 20 C. Using Koch et al. 2000 parameters.
+     * @property {Object} kin_params - Kinetic parameters
+     */
+    this.kin_params = KINETIC_CONSTANTS;
+
+    /**
+     * Stoichiometric and composition parameters for ASM3. Using Koch et al. 2000 parameters.
+     * @property {Object} stoi_params - Stoichiometric parameters
+     */
+    this.stoi_params = STOICHIOMETRIC_CONSTANTS;
+
+    /**
+     * Temperature theta parameters for ASM3. Using Koch et al. 2000 parameters.
+     * These parameters are used to adjust reaction rates based on temperature.
+     * @property {Object} temp_params - Temperature theta parameters
+     */
+    this.temp_params = {
+      // Hydrolysis
+      theta_H:        0.04,
+      // Heterotrophs
+      theta_STO:      0.07,
+      theta_mu_H:     0.07,
+      theta_b_H_O:    0.07,
+      theta_b_H_NO:   0.07,
+      theta_b_STO_O:  this._compute_theta(0.1, 0.3, 10, 20),
+      theta_b_STO_NO: this._compute_theta(0.05, 0.15, 10, 20),
+      // Autotrophs
+      theta_mu_A:     0.105,
+      theta_b_A_O:    0.105,
+      theta_b_A_NO:   0.105
+    };
+
+    this.stoi_matrix = this._initialise_stoi_matrix();
+  }
+
+  /**
+   * Initialises the stoichiometric matrix for ASM3.
+   * @returns {Array} - The stoichiometric matrix for ASM3. (2D array)
+   */
+  _initialise_stoi_matrix() { // initialise stoichiometric matrix
+    const { f_SI, f_XI, Y_STO_O, Y_STO_NO, Y_H_O, Y_H_NO, Y_A, i_CODN, i_CODNO, i_NSI, i_NSS, i_NXI, i_NXS, i_NBM, i_TSXI, i_TSXS, i_TSBM, i_TSSTO, i_cNH, i_cNO } = this.stoi_params;
+
+    const stoi_matrix = Array(12);
+    // S_O, S_I, S_S, S_NH, S_N2, S_NO, S_HCO, X_I, X_S, X_H, X_STO, X_A, X_TS
+    stoi_matrix[0]  = [0., f_SI, 1.-f_SI, i_NXS-(1.-f_SI)*i_NSS-f_SI*i_NSI, 0., 0., (i_NXS-(1.-f_SI)*i_NSS-f_SI*i_NSI)*i_cNH, 0., -1., 0., 0., 0., -i_TSXS];
+    stoi_matrix[1]  = [-(1.-Y_STO_O), 0, -1., i_NSS, 0., 0., i_NSS*i_cNH, 0., 0., 0., Y_STO_O, 0., Y_STO_O*i_TSSTO];
+    stoi_matrix[2]  = [0., 0., -1., i_NSS, -(1.-Y_STO_NO)/(i_CODNO-i_CODN), (1.-Y_STO_NO)/(i_CODNO-i_CODN), i_NSS*i_cNH + (1.-Y_STO_NO)/(i_CODNO-i_CODN)*i_cNO, 0., 0., 0., Y_STO_NO, 0., Y_STO_NO*i_TSSTO];
+    stoi_matrix[3]  = [-(1.-Y_H_O)/Y_H_O, 0., 0., -i_NBM, 0., 0., -i_NBM*i_cNH, 0., 0., 1., -1./Y_H_O, 0., i_TSBM-i_TSSTO/Y_H_O];
+    stoi_matrix[4]  = [0., 0., 0., -i_NBM, -(1.-Y_H_NO)/(Y_H_NO*(i_CODNO-i_CODN)), (1.-Y_H_NO)/(Y_H_NO*(i_CODNO-i_CODN)), -i_NBM*i_cNH+(1.-Y_H_NO)/(Y_H_NO*(i_CODNO-i_CODN))*i_cNO, 0., 0., 1., -1./Y_H_NO, 0., i_TSBM-i_TSSTO/Y_H_NO];
+    stoi_matrix[5]  = [f_XI-1., 0., 0., i_NBM-f_XI*i_NXI, 0., 0., (i_NBM-f_XI*i_NXI)*i_cNH, f_XI, 0., -1., 0., 0., f_XI*i_TSXI-i_TSBM];
+    stoi_matrix[6]  = [0., 0., 0., i_NBM-f_XI*i_NXI, -(1.-f_XI)/(i_CODNO-i_CODN), (1.-f_XI)/(i_CODNO-i_CODN), (i_NBM-f_XI*i_NXI)*i_cNH+(1-f_XI)/(i_CODNO-i_CODN)*i_cNO, f_XI, 0., -1., 0., 0., f_XI*i_TSXI-i_TSBM];
+    stoi_matrix[7]  = [-1., 0., 0., 0., 0., 0., 0., 0., 0., 0., -1., 0., -i_TSSTO];
+    stoi_matrix[8]  = [0., 0., 0., 0., -1./(i_CODNO-i_CODN), 1./(i_CODNO-i_CODN), i_cNO/(i_CODNO-i_CODN), 0., 0., 0., -1., 0., -i_TSSTO];
+    stoi_matrix[9]  = [1.+i_CODNO/Y_A, 0., 0., -1./Y_A-i_NBM, 0., 1./Y_A, (-1./Y_A-i_NBM)*i_cNH+i_cNO/Y_A, 0., 0., 0., 0., 1., i_TSBM];
+    stoi_matrix[10] = [f_XI-1., 0., 0., i_NBM-f_XI*i_NXI, 0., 0., (i_NBM-f_XI*i_NXI)*i_cNH, f_XI, 0., 0., 0., -1., f_XI*i_TSXI-i_TSBM];
+    stoi_matrix[11] = [0., 0., 0., i_NBM-f_XI*i_NXI, -(1.-f_XI)/(i_CODNO-i_CODN), (1.-f_XI)/(i_CODNO-i_CODN), (i_NBM-f_XI*i_NXI)*i_cNH+(1-f_XI)/(i_CODNO-i_CODN)*i_cNO, 0., 0., 0., 0., -1., f_XI*i_TSXI-i_TSBM];
+
+    return stoi_matrix[0].map((col, i) => stoi_matrix.map(row => row[i])); // transpose matrix
+  }
+
+  /**
+   * Computes the Monod equation rate value for a given concentration and half-saturation constant.
+   * @param {number} c - Concentration of reaction species.
+   * @param {number} K - Half-saturation constant for the reaction species.
+   * @returns {number} - Monod equation rate value for the given concentration and half-saturation constant.
+   */
+  _monod(c, K) {
+    return c / (K + c);
+  }
+
+  /**
+   * Computes the inverse Monod equation rate value for a given concentration and half-saturation constant. Used for inhibition.
+   * @param {number} c - Concentration of reaction species.
+   * @param {number} K - Half-saturation constant for the reaction species. 
+   * @returns {number} - Inverse Monod equation rate value for the given concentration and half-saturation constant.
+   */
+  _inv_monod(c, K) {
+    return K / (K + c);
+  }
+
+  /**
+   * Adjust the rate parameter for temperature T using simplied Arrhenius equation based on rate constant at 20 degrees Celsius and theta parameter.
+   * @param {number} k - Rate constant at 20 degrees Celcius.
+   * @param {number} theta - Theta parameter.
+   * @param {number} T - Temperature in Celcius.
+   * @returns {number} - Adjusted rate parameter at temperature T based on the Arrhenius equation.
+   */
+  _arrhenius(k, theta, T) {
+    return k * Math.exp(theta*(T-20));
+  }
+
+  /**
+   * Computes the temperature theta parameter based on two rate constants and their corresponding temperatures.
+   * @param {number} k1 - Rate constant at temperature T1.
+   * @param {number} k2 - Rate constant at temperature T2.
+   * @param {number} T1 - Temperature T1 in Celcius.
+   * @param {number} T2 - Temperature T2 in Celcius.
+   * @returns {number} - Theta parameter.
+   */
+  _compute_theta(k1, k2, T1, T2) {
+    return Math.log(k1/k2)/(T1-T2);
+  }
+
+  /**
+   * Computes the reaction rates for each process reaction based on the current state and temperature.
+   * @param {Array} state - State vector containing concentrations of reaction species.
+   * @param {number} [T=20] - Temperature in degrees Celsius (default is 20).
+   * @returns {Array} - Reaction rates for each process reaction.
+   */
+  compute_rates(state, T = 20) {
+    // state: S_O, S_I, S_S, S_NH, S_N2, S_NO, S_HCO, X_I, X_S, X_H, X_STO, X_A, X_TS
+    const rates = Array(12);
+    const [S_O, S_I, S_S, S_NH, S_N2, S_NO, S_HCO, X_I, X_S, X_H, X_STO, X_A, X_TS] = state;
+    const { k_H, K_X, k_STO, nu_NO, K_O, K_NO, K_S, K_STO, mu_H_max, K_NH, K_HCO, b_H_O, b_H_NO, b_STO_O, b_STO_NO, mu_A_max, K_A_NH, K_A_O, K_A_HCO, b_A_O, b_A_NO } = this.kin_params;
+    const { theta_H, theta_STO, theta_mu_H, theta_b_H_O, theta_b_H_NO, theta_b_STO_O, theta_b_STO_NO, theta_mu_A, theta_b_A_O, theta_b_A_NO } = this.temp_params;
+
+    // Hydrolysis
+    rates[0] = X_H == 0 ? 0 : this._arrhenius(k_H, theta_H, T) * this._monod(X_S / X_H, K_X) * X_H;
+
+    // Heterotrophs
+    rates[1] = this._arrhenius(k_STO, theta_STO, T) * this._monod(S_O, K_O) * this._monod(S_S, K_S) * X_H;
+    rates[2] = this._arrhenius(k_STO, theta_STO, T) * nu_NO * this._inv_monod(S_O, K_O) * this._monod(S_NO, K_NO) * this._monod(S_S, K_S) * X_H;
+    rates[3] = X_H == 0 ? 0 : this._arrhenius(mu_H_max, theta_mu_H, T) * this._monod(S_O, K_O) * this._monod(S_NH, K_NH) * this._monod(S_HCO, K_HCO) * this._monod(X_STO/X_H, K_STO) * X_H;
+    rates[4] = X_H == 0 ? 0 : this._arrhenius(mu_H_max, theta_mu_H, T) * nu_NO * this._inv_monod(S_O, K_O) * this._monod(S_NO, K_NO) * this._monod(S_NH, K_NH) * this._monod(S_HCO, K_HCO) * this._monod(X_STO/X_H, K_STO) * X_H;
+    rates[5] = this._arrhenius(b_H_O, theta_b_H_O, T) * this._monod(S_O, K_O) * X_H;
+    rates[6] = this._arrhenius(b_H_NO, theta_b_H_NO, T) * this._inv_monod(S_O, K_O) * this._monod(S_NO, K_NO) * X_H;
+    rates[7] = this._arrhenius(b_STO_O, theta_b_STO_O, T) * this._monod(S_O, K_O) * X_H;
+    rates[8] = this._arrhenius(b_STO_NO, theta_b_STO_NO, T) * this._inv_monod(S_O, K_O) * this._monod(S_NO, K_NO) * X_STO;
+    
+    // Autotrophs
+    rates[9] = this._arrhenius(mu_A_max, theta_mu_A, T) * this._monod(S_O, K_A_O) * this._monod(S_NH, K_A_NH) * this._monod(S_HCO, K_A_HCO) * X_A;
+    rates[10] = this._arrhenius(b_A_O, theta_b_A_O, T) * this._monod(S_O, K_O) * X_A;
+    rates[11] = this._arrhenius(b_A_NO, theta_b_A_NO, T) * this._inv_monod(S_O, K_A_O) * this._monod(S_NO, K_NO) * X_A;
+    
+    return rates;
+  }
+
+  /**
+   * Computes the change in concentrations of reaction species based on the current state and temperature.
+   * @param {Array} state - State vector containing concentrations of reaction species.
+   * @param {number} [T=20] - Temperature in degrees Celsius (default is 20).
+   * @returns {Array} - Change in reaction species concentrations.
+   */
+  compute_dC(state, T = 20) { // compute changes in concentrations
+    // state: S_O, S_I, S_S, S_NH, S_N2, S_NO, S_HCO, X_I, X_S, X_H, X_STO, X_A, X_TS
+    return math.multiply(this.stoi_matrix, this.compute_rates(state, T));
+  }
+}
+
+module.exports = { ASM3, ASM_CONSTANTS, KINETIC_CONSTANTS, STOICHIOMETRIC_CONSTANTS };

src/reaction_modules/asm3_class.js

@@ -0,0 +1,222 @@
+const math = require('mathjs');
+
+const ASM_CONSTANTS = {
+  S_O_INDEX: 0,
+  S_NH_INDEX: 3,
+  S_NO_INDEX: 5,
+  NUM_SPECIES: 13
+};
+
+const KINETIC_CONSTANTS = {
+  // Hydrolysis
+  k_H:      3.,   // hydrolysis rate constant       [g X_S g-1 X_H d-1]
+  K_X:      1.,   // hydrolysis saturation constant [g X_S g-1 X_H] 
+  // Heterotrophs
+  k_STO:    5.,   // storage rate constant          [g S_S g-1 X_H d-1]
+  nu_NO:    0.6,  // anoxic reduction factor        [-]
+  K_O:      0.2,  // saturation constant S_0        [g O2 m-3]
+  K_NO:     0.5,  // saturation constant S_NO       [g NO3-N m-3]
+  K_S:      2.,   // saturation constant S_s        [g COD m-3]
+  K_STO:    1.,   // saturation constant X_STO      [g X_STO g-1 X_H]
+  mu_H_max: 2.,   // maximum specific growth rate   [d-1]
+  K_NH:     0.01, // saturation constant S_NH3      [g NH3-N m-3]
+  K_HCO:    0.1,  // saturation constant S_HCO      [mole HCO3 m-3]
+  b_H_O:    0.2,  // aerobic respiration rate       [d-1]
+  b_H_NO:   0.1,  // anoxic respiration rate        [d-1]
+  b_STO_O:  0.2,  // aerobic respitation rate X_STO [d-1]
+  b_STO_NO: 0.1,  // anoxic respitation rate X_STO  [d-1]
+  // Autotrophs
+  mu_A_max: 1.0,  // maximum specific growth rate   [d-1]
+  K_A_NH:   1.,   // saturation constant S_NH3      [g NH3-N m-3] 
+  K_A_O:    0.5,  // saturation constant S_0        [g O2 m-3]
+  K_A_HCO:  0.5,  // saturation constant S_HCO      [mole HCO3 m-3]
+  b_A_O:    0.15, // aerobic respiration rate       [d-1]
+  b_A_NO:   0.05  // anoxic respiration rate        [d-1]
+};
+
+const STOICHIOMETRIC_CONSTANTS = {
+  // Fractions
+  f_SI:     0.,   // fraction S_I from hydrolysis   [g S_I g-1 X_S]
+  f_XI:     0.2,  // fraction X_I from decomp X_H   [g X_I g-1 X_H]
+  // Yields
+  Y_STO_O:  0.85, // aerobic yield X_STO per S_S    [g X_STO g-1 S_S]
+  Y_STO_NO: 0.80, // anoxic yield X_STO per S_S     [g X_STO g-1 S_S]
+  Y_H_O:    0.63, // aerobic yield X_H per X_STO    [g X_H g-1 X_STO]
+  Y_H_NO:   0.54, // anoxic yield X_H per X_STO     [g X_H g-1 X_STO]
+  Y_A:      0.24, // anoxic yield X_A per S_NO      [g X_A g-1 NO3-N]
+  // Composition (COD via DoR)
+  i_CODN:  -1.71, // COD content (DoR)              [g COD g-1 N2-N]
+  i_CODNO: -4.57, // COD content (DoR)              [g COD g-1 NO3-N]
+  // Composition (nitrogen)
+  i_NSI:    0.01, // nitrogen content S_I           [g N g-1 S_I]
+  i_NSS:    0.03, // nitrogen content S_S           [g N g-1 S_S]
+  i_NXI:    0.02, // nitrogen content X_I           [g N g-1 X_I]
+  i_NXS:    0.04, // nitrogen content X_S           [g N g-1 X_S]
+  i_NBM:    0.07, // nitrogen content X_H / X_A     [g N g-1 X_H / X_A]
+  // Composition (TSS)
+  i_TSXI:   0.75, // TSS content X_I                [g TS g-1 X_I]
+  i_TSXS:   0.75, // TSS content X_S                [g TS g-1 X_S]
+  i_TSBM:   0.90, // TSS content X_H / X_A          [g TS g-1 X_H / X_A]
+  i_TSSTO:  0.60, // TSS content X_STO (PHB based)  [g TS g-1 X_STO]
+  // Composition (charge)
+  i_cNH:    1/14, // charge per S_NH                [mole H+ g-1 NH3-N]
+  i_cNO:   -1/14  // charge per S_NO                [mole H+ g-1 NO3-N]
+};
+
+/**
+ * ASM3 class for the Activated Sludge Model No. 3 (ASM3).
+ */
+class ASM3 {
+
+  constructor() {
+    /**
+     * Kinetic parameters for ASM3 at 20 C.
+     * @property {Object} kin_params - Kinetic parameters
+     */
+    this.kin_params = KINETIC_CONSTANTS;
+
+    /**
+     * Stoichiometric and composition parameters for ASM3.
+     * @property {Object} stoi_params - Stoichiometric parameters
+     */
+    this.stoi_params = STOICHIOMETRIC_CONSTANTS;
+
+    /**
+     * Temperature theta parameters for ASM3.
+     * These parameters are used to adjust reaction rates based on temperature.
+     * @property {Object} temp_params - Temperature theta parameters
+     */
+    this.temp_params = {
+      // Hydrolysis
+      theta_H:        this._compute_theta(2, 3, 10, 20),
+      // Heterotrophs
+      theta_STO:      this._compute_theta(2.5, 5, 10, 20),
+      theta_mu_H:     this._compute_theta(1, 2, 10, 20),
+      theta_b_H_O:    this._compute_theta(0.1, 0.2, 10, 20),
+      theta_b_H_NO:   this._compute_theta(0.05, 0.1, 10, 20),
+      theta_b_STO_O:  this._compute_theta(0.1, 0.2, 10, 20),
+      theta_b_STO_NO: this._compute_theta(0.05, 0.1, 10, 20),
+      // Autotrophs
+      theta_mu_A:     this._compute_theta(0.35, 1, 10, 20),
+      theta_b_A_O:    this._compute_theta(0.05, 0.15, 10, 20),
+      theta_b_A_NO:   this._compute_theta(0.02, 0.05, 10, 20)
+    };
+
+    this.stoi_matrix = this._initialise_stoi_matrix();
+  }
+
+  /**
+   * Initialises the stoichiometric matrix for ASM3.
+   * @returns {Array} - The stoichiometric matrix for ASM3. (2D array)
+   */
+  _initialise_stoi_matrix() { // initialise stoichiometric matrix
+    const { f_SI, f_XI, Y_STO_O, Y_STO_NO, Y_H_O, Y_H_NO, Y_A, i_CODN, i_CODNO, i_NSI, i_NSS, i_NXI, i_NXS, i_NBM, i_TSXI, i_TSXS, i_TSBM, i_TSSTO, i_cNH, i_cNO } = this.stoi_params;
+
+    const stoi_matrix = Array(12);
+    // S_O, S_I, S_S, S_NH, S_N2, S_NO, S_HCO, X_I, X_S, X_H, X_STO, X_A, X_TS
+    stoi_matrix[0]  = [0., f_SI, 1.-f_SI, i_NXS-(1.-f_SI)*i_NSS-f_SI*i_NSI, 0., 0., (i_NXS-(1.-f_SI)*i_NSS-f_SI*i_NSI)*i_cNH, 0., -1., 0., 0., 0., -i_TSXS];
+    stoi_matrix[1]  = [-(1.-Y_STO_O), 0, -1., i_NSS, 0., 0., i_NSS*i_cNH, 0., 0., 0., Y_STO_O, 0., Y_STO_O*i_TSSTO];
+    stoi_matrix[2]  = [0., 0., -1., i_NSS, -(1.-Y_STO_NO)/(i_CODNO-i_CODN), (1.-Y_STO_NO)/(i_CODNO-i_CODN), i_NSS*i_cNH + (1.-Y_STO_NO)/(i_CODNO-i_CODN)*i_cNO, 0., 0., 0., Y_STO_NO, 0., Y_STO_NO*i_TSSTO];
+    stoi_matrix[3]  = [-(1.-Y_H_O)/Y_H_O, 0., 0., -i_NBM, 0., 0., -i_NBM*i_cNH, 0., 0., 1., -1./Y_H_O, 0., i_TSBM-i_TSSTO/Y_H_O];
+    stoi_matrix[4]  = [0., 0., 0., -i_NBM, -(1.-Y_H_NO)/(Y_H_NO*(i_CODNO-i_CODN)), (1.-Y_H_NO)/(Y_H_NO*(i_CODNO-i_CODN)), -i_NBM*i_cNH+(1.-Y_H_NO)/(Y_H_NO*(i_CODNO-i_CODN))*i_cNO, 0., 0., 1., -1./Y_H_NO, 0., i_TSBM-i_TSSTO/Y_H_NO];
+    stoi_matrix[5]  = [f_XI-1., 0., 0., i_NBM-f_XI*i_NXI, 0., 0., (i_NBM-f_XI*i_NXI)*i_cNH, f_XI, 0., -1., 0., 0., f_XI*i_TSXI-i_TSBM];
+    stoi_matrix[6]  = [0., 0., 0., i_NBM-f_XI*i_NXI, -(1.-f_XI)/(i_CODNO-i_CODN), (1.-f_XI)/(i_CODNO-i_CODN), (i_NBM-f_XI*i_NXI)*i_cNH+(1-f_XI)/(i_CODNO-i_CODN)*i_cNO, f_XI, 0., -1., 0., 0., f_XI*i_TSXI-i_TSBM];
+    stoi_matrix[7]  = [-1., 0., 0., 0., 0., 0., 0., 0., 0., 0., -1., 0., -i_TSSTO];
+    stoi_matrix[8]  = [0., 0., 0., 0., -1./(i_CODNO-i_CODN), 1./(i_CODNO-i_CODN), i_cNO/(i_CODNO-i_CODN), 0., 0., 0., -1., 0., -i_TSSTO];
+    stoi_matrix[9]  = [1.+i_CODNO/Y_A, 0., 0., -1./Y_A-i_NBM, 0., 1./Y_A, (-1./Y_A-i_NBM)*i_cNH+i_cNO/Y_A, 0., 0., 0., 0., 1., i_TSBM];
+    stoi_matrix[10] = [f_XI-1., 0., 0., i_NBM-f_XI*i_NXI, 0., 0., (i_NBM-f_XI*i_NXI)*i_cNH, f_XI, 0., 0., 0., -1., f_XI*i_TSXI-i_TSBM];
+    stoi_matrix[11] = [0., 0., 0., i_NBM-f_XI*i_NXI, -(1.-f_XI)/(i_CODNO-i_CODN), (1.-f_XI)/(i_CODNO-i_CODN), (i_NBM-f_XI*i_NXI)*i_cNH+(1-f_XI)/(i_CODNO-i_CODN)*i_cNO, 0., 0., 0., 0., -1., f_XI*i_TSXI-i_TSBM];
+
+    return stoi_matrix[0].map((col, i) => stoi_matrix.map(row => row[i])); // transpose matrix
+  }
+
+  /**
+   * Computes the Monod equation rate value for a given concentration and half-saturation constant.
+   * @param {number} c - Concentration of reaction species.
+   * @param {number} K - Half-saturation constant for the reaction species.
+   * @returns {number} - Monod equation rate value for the given concentration and half-saturation constant.
+   */
+  _monod(c, K) {
+    return c / (K + c);
+  }
+
+  /**
+   * Computes the inverse Monod equation rate value for a given concentration and half-saturation constant. Used for inhibition.
+   * @param {number} c - Concentration of reaction species.
+   * @param {number} K - Half-saturation constant for the reaction species. 
+   * @returns {number} - Inverse Monod equation rate value for the given concentration and half-saturation constant.
+   */
+  _inv_monod(c, K) {
+    return K / (K + c);
+  }
+
+  /**
+   * Adjust the rate parameter for temperature T using simplied Arrhenius equation based on rate constant at 20 degrees Celsius and theta parameter.
+   * @param {number} k - Rate constant at 20 degrees Celcius.
+   * @param {number} theta - Theta parameter.
+   * @param {number} T - Temperature in Celcius.
+   * @returns {number} - Adjusted rate parameter at temperature T based on the Arrhenius equation.
+   */
+  _arrhenius(k, theta, T) {
+    return k * Math.exp(theta*(T-20));
+  }
+
+  /**
+   * Computes the temperature theta parameter based on two rate constants and their corresponding temperatures.
+   * @param {number} k1 - Rate constant at temperature T1.
+   * @param {number} k2 - Rate constant at temperature T2.
+   * @param {number} T1 - Temperature T1 in Celcius.
+   * @param {number} T2 - Temperature T2 in Celcius.
+   * @returns {number} - Theta parameter.
+   */
+  _compute_theta(k1, k2, T1, T2) {
+    return Math.log(k1/k2)/(T1-T2);
+  }
+
+  /**
+   * Computes the reaction rates for each process reaction based on the current state and temperature.
+   * @param {Array} state - State vector containing concentrations of reaction species.
+   * @param {number} [T=20] - Temperature in degrees Celsius (default is 20).
+   * @returns {Array} - Reaction rates for each process reaction.
+   */
+  compute_rates(state, T = 20) {
+    // state: S_O, S_I, S_S, S_NH, S_N2, S_NO, S_HCO, X_I, X_S, X_H, X_STO, X_A, X_TS
+    const rates = Array(12);
+    const [S_O, S_I, S_S, S_NH, S_N2, S_NO, S_HCO, X_I, X_S, X_H, X_STO, X_A, X_TS] = state;
+    const { k_H, K_X, k_STO, nu_NO, K_O, K_NO, K_S, K_STO, mu_H_max, K_NH, K_HCO, b_H_O, b_H_NO, b_STO_O, b_STO_NO, mu_A_max, K_A_NH, K_A_O, K_A_HCO, b_A_O, b_A_NO } = this.kin_params;
+    const { theta_H, theta_STO, theta_mu_H, theta_b_H_O, theta_b_H_NO, theta_b_STO_O, theta_b_STO_NO, theta_mu_A, theta_b_A_O, theta_b_A_NO } = this.temp_params;
+
+    // Hydrolysis
+    rates[0] = X_H == 0 ? 0 : this._arrhenius(k_H, theta_H, T) * this._monod(X_S / X_H, K_X) * X_H;
+
+    // Heterotrophs
+    rates[1] = this._arrhenius(k_STO, theta_STO, T) * this._monod(S_O, K_O) * this._monod(S_S, K_S) * X_H;
+    rates[2] = this._arrhenius(k_STO, theta_STO, T) * nu_NO * this._inv_monod(S_O, K_O) * this._monod(S_NO, K_NO) * this._monod(S_S, K_S) * X_H;
+    rates[3] = X_H == 0 ? 0 : this._arrhenius(mu_H_max, theta_mu_H, T) * this._monod(S_O, K_O) * this._monod(S_NH, K_NH) * this._monod(S_HCO, K_HCO) * this._monod(X_STO/X_H, K_STO) * X_H;
+    rates[4] = X_H == 0 ? 0 : this._arrhenius(mu_H_max, theta_mu_H, T) * nu_NO * this._inv_monod(S_O, K_O) * this._monod(S_NO, K_NO) * this._monod(S_NH, K_NH) * this._monod(S_HCO, K_HCO) * this._monod(X_STO/X_H, K_STO) * X_H;
+    rates[5] = this._arrhenius(b_H_O, theta_b_H_O, T) * this._monod(S_O, K_O) * X_H;
+    rates[6] = this._arrhenius(b_H_NO, theta_b_H_NO, T) * this._inv_monod(S_O, K_O) * this._monod(S_NO, K_NO) * X_H;
+    rates[7] = this._arrhenius(b_STO_O, theta_b_STO_O, T) * this._monod(S_O, K_O) * X_H;
+    rates[8] = this._arrhenius(b_STO_NO, theta_b_STO_NO, T) * this._inv_monod(S_O, K_O) * this._monod(S_NO, K_NO) * X_STO;
+    
+    // Autotrophs
+    rates[9] = this._arrhenius(mu_A_max, theta_mu_A, T) * this._monod(S_O, K_A_O) * this._monod(S_NH, K_A_NH) * this._monod(S_HCO, K_A_HCO) * X_A;
+    rates[10] = this._arrhenius(b_A_O, theta_b_A_O, T) * this._monod(S_O, K_O) * X_A;
+    rates[11] = this._arrhenius(b_A_NO, theta_b_A_NO, T) * this._inv_monod(S_O, K_A_O) * this._monod(S_NO, K_NO) * X_A;
+    
+    return rates;
+  }
+
+  /**
+   * Computes the change in concentrations of reaction species based on the current state and temperature.
+   * @param {Array} state - State vector containing concentrations of reaction species.
+   * @param {number} [T=20] - Temperature in degrees Celsius (default is 20).
+   * @returns {Array} - Change in reaction species concentrations.
+   */
+  compute_dC(state, T = 20) { // compute changes in concentrations
+    // state: S_O, S_I, S_S, S_NH, S_N2, S_NO, S_HCO, X_I, X_S, X_H, X_STO, X_A, X_TS
+    return math.multiply(this.stoi_matrix, this.compute_rates(state, T));
+  }
+}
+
+module.exports = { ASM3, ASM_CONSTANTS, KINETIC_CONSTANTS, STOICHIOMETRIC_CONSTANTS };

src/specificClass.js

@@ -0,0 +1,488 @@
+const { ASM3, ASM_CONSTANTS } = require('./reaction_modules/asm3_class.js');
+const { create, all, isArray } = require('mathjs');
+const { assertNoNaN } = require('./utils.js');
+const { childRegistrationUtils, logger, MeasurementContainer } = require('generalFunctions');
+const EventEmitter = require('events');
+
+const mathConfig = {
+  matrix: 'Array' // use Array as the matrix type
+};
+
+const math = create(all, mathConfig);
+
+const BC_PADDING = 2; // Boundary Condition padding for open boundaries in extendedState variable
+const DEBUG = false;
+const DAY2MS = 1000 * 60 * 60 * 24; // convert between days and milliseconds
+
+class Reactor {
+  /**
+   * Reactor base class.
+   * @param {object} config - Configuration object containing reactor parameters.
+   */
+  constructor(config) {
+    this.config = config;
+    // EVOLV stuff
+    this.logger = new logger(this.config.general.logging.enabled, this.config.general.logging.logLevel, config.general.name);
+    this.emitter = new EventEmitter();
+    this.measurements = new MeasurementContainer();
+    this.childRegistrationUtils = new childRegistrationUtils(this); // child registration utility
+
+    // placeholder variables for children and parents
+    this.upstreamReactor = null;
+    this.downstreamReactor = null;
+    this.returnPump = null;
+
+    this.asm = new ASM3(); // Reaction model
+
+    this.volume = config.volume; // fluid volume reactor [m3]
+
+    this.Fs = [0]; // fluid debits per inlet [m3 d-1]
+    this.Cs_in = [Array(ASM_CONSTANTS.NUM_SPECIES).fill(0)]; // composition influents
+    this.OTR = 0.0; // oxygen transfer rate [g O2 d-1 m-3]
+    this.temperature = 20; // temperature [C]
+
+    this.kla = config.kla; // if NaN, use externaly provided OTR [d-1]
+
+    this.currentTime = null; // milliseconds since epoch [ms]
+    this.timeStep = 1 / (24*60*60) * this.config.timeStep; // time step in seconds, converted to days.
+    this.speedUpFactor = 1; // speed up factor for simulation, 60 means 1 minute per simulated second
+  }
+
+  /**
+   * Setter for influent data.
+   * @param {object} input - Input object (msg) containing payload with inlet index, flow rate, and concentrations.
+  */
+  set setInfluent(input) {
+    const i_in = input.payload.inlet;
+    if (this.Fs.length <= i_in) {
+      this.logger.debug(`Adding new inlet index ${i_in}.`);
+      this.Fs.push(0);
+      this.Cs_in.push(Array(ASM_CONSTANTS.NUM_SPECIES).fill(0));
+      this.setInfluent = input;
+    }
+    this.Fs[i_in] = input.payload.F;
+    this.Cs_in[i_in] = input.payload.C;
+  }
+
+  /**
+   * Setter for OTR (Oxygen Transfer Rate).
+   * @param {object} input - Input object (msg) containing payload with OTR value [g O2 d-1 m-3].
+   */
+  set setOTR(input) {
+    this.OTR = input.payload;
+  }
+
+  /**
+   * Getter for effluent data.
+   * @returns {object} Effluent data object (msg).
+   */
+  get getEffluent() {
+    const Cs = isArray(this.state.at(-1)) ? this.state.at(-1) : this.state;
+    const effluent = [{ topic: "Fluent", payload: { inlet: 0, F: math.sum(this.Fs), C: Cs }, timestamp: this.currentTime }];
+    if (this.returnPump) {
+      const recirculationFlow = this.returnPump.measurements.type("flow").variant("measured").position("atEquipment").getCurrentValue();
+      // constrain flow to prevent negatives
+      const F_main = Math.max(effluent[0].payload.F - recirculationFlow, 0);
+      const F_sidestream = effluent[0].payload.F < recirculationFlow ? effluent[0].payload.F : recirculationFlow;
+      effluent[0].payload.F = F_main;
+      effluent.push({ topic: "Fluent", payload: { inlet: 1, F: F_sidestream, C: Cs }, timestamp: this.currentTime });
+    }
+    return effluent;
+  }
+
+  /**
+   * Calculate the oxygen transfer rate (OTR) based on the dissolved oxygen concentration and temperature.
+   * @param {number} S_O - Dissolved oxygen concentration [g O2 m-3].
+   * @param {number} T - Temperature in Celsius, default to 20 C.
+   * @returns {number} - Calculated OTR [g O2 d-1 m-3].
+   */
+  _calcOTR(S_O, T = 20.0) { // caculate the OTR using basic correlation, default to temperature: 20 C
+    const S_O_sat = 14.652 - 4.1022e-1 * T + 7.9910e-3 * T*T + 7.7774e-5 * T*T*T;
+    return this.kla * (S_O_sat - S_O);
+  }
+
+  /**
+   * Clip values in an array to zero.
+   * @param {Array} arr - Array of values to clip.
+   * @returns {Array} - New array with values clipped to zero.
+   */
+  _arrayClip2Zero(arr) {
+    if (Array.isArray(arr)) {
+      return arr.map(x => this._arrayClip2Zero(x));
+    } else {
+      return arr < 0 ? 0 : arr;
+    }
+  }
+
+  /**
+   * Register child function required for child registration.
+   * @param {object} child 
+   * @param {string} softwareType 
+   */
+  registerChild(child, softwareType) {
+    if(!child) {
+      this.logger.error(`Invalid ${softwareType} child provided.`);
+      return;
+    }
+
+    switch (softwareType) {
+      case "measurement":
+        this.logger.debug(`Registering measurement child...`);
+        this._connectMeasurement(child);
+        break;
+      case "reactor":
+        this.logger.debug(`Registering reactor child...`);
+        this._connectReactor(child);
+        break;
+      case "machine":
+        this.logger.debug(`Registering machine child...`);
+        this._connectMachine(child);
+        break;
+
+      default:
+        this.logger.error(`Unrecognized softwareType: ${softwareType}`);
+    }
+  }
+
+  _connectMeasurement(measurementChild) {
+    const position = measurementChild.config.functionality.positionVsParent;
+    const measurementType = measurementChild.config.asset.type;
+    const eventName = `${measurementType}.measured.${position}`;
+
+    // Register event listener for measurement updates
+    measurementChild.measurements.emitter.on(eventName, (eventData) => {
+      this.logger.debug(`${position} ${measurementType} from ${eventData.childName}: ${eventData.value} ${eventData.unit}`);
+
+      // Store directly in parent's measurement container
+      this.measurements
+        .type(measurementType)
+        .variant("measured")
+        .position(position)
+        .value(eventData.value, eventData.timestamp, eventData.unit);
+      
+      this._updateMeasurement(measurementType, eventData.value, position, eventData);
+    });
+  }
+
+
+  _connectReactor(reactorChild) {
+    if (reactorChild.config.functionality.positionVsParent != "upstream") {
+      this.logger.warn("Reactor children of other reactors should always be upstream!");
+    }
+
+    // set upstream and downstream reactor variable in current and child nodes respectively for easy access
+    this.upstreamReactor = reactorChild;
+    reactorChild.downstreamReactor = this;
+
+    reactorChild.emitter.on("stateChange", (eventData) => { // Triggers state update in downstream reactor.
+      this.logger.debug(`State change of upstream reactor detected.`);
+      this.updateState(eventData);
+    });
+  }
+
+  _connectMachine(machineChild) {
+    if (machineChild.config.functionality.positionVsParent == "downstream") {
+      machineChild.upstreamSource = this;
+      this.returnPump = machineChild;
+    }
+  }
+  
+  _updateMeasurement(measurementType, value, position, context) {
+    this.logger.debug(`---------------------- updating ${measurementType} ------------------ `);
+    switch (measurementType) {
+      case "temperature":
+        if (position == "atEquipment") {
+          this.temperature = value;
+        }
+        break;
+      default:
+        this.logger.error(`Type '${measurementType}' not recognized for measured update.`);
+        return;
+    }
+  }
+
+  /**
+   * Update the reactor state based on the new time.
+   * @param {number} newTime - New time to update reactor state to, in milliseconds since epoch.
+   */
+  updateState(newTime) {
+    if (!this.currentTime) { // initialise currentTime variable
+      this.currentTime = newTime;
+      return;
+    }
+
+    if (this.upstreamReactor) { // grab main effluent upstream reactor
+      this.setInfluent = this.upstreamReactor.getEffluent[0];
+    }
+
+    const n_iter = Math.floor(this.speedUpFactor * (newTime-this.currentTime) / (this.timeStep*DAY2MS));
+
+    if (n_iter == 0) { // no update required, change in currentTime smaller than time step
+      return;
+    }
+
+    let n = 0;
+    while (n < n_iter) {
+      this.tick(this.timeStep);
+      n += 1;
+    }
+    this.currentTime += n_iter * this.timeStep * DAY2MS / this.speedUpFactor;
+    this.emitter.emit("stateChange", this.currentTime); // update downstream reactors
+
+    if (this.returnPump) { // update recirculation pump state
+      this.returnPump.updateSourceSink();
+    }
+  }
+}
+
+class Reactor_CSTR extends Reactor {
+  /**
+   * Reactor_CSTR class for Continuous Stirred Tank Reactor.
+   * @param {object} config - Configuration object containing reactor parameters.
+   */
+  constructor(config) {
+      super(config);
+      this.state = config.initialState;
+  }
+
+  _updateMeasurement(measurementType, value, position, context) {
+
+    switch(measurementType) {
+      case "quantity (oxygen)":
+        this.state[ASM_CONSTANTS.S_O_INDEX] = value;
+        break;
+      case "quantity (ammonium)":
+        this.state[ASM_CONSTANTS.S_NH_INDEX] = value;
+        break;
+      case "quantity (nox)":
+        this.state[ASM_CONSTANTS.S_NO_INDEX] = value;
+        break;
+      default:
+        super._updateMeasurement(measurementType, value, position, context);
+    }
+  }
+
+  /**
+   * Tick the reactor state using the forward Euler method.
+   * @param {number} time_step - Time step for the simulation [d].
+   * @returns {Array} - New reactor state.
+   */
+  tick(time_step) { // tick reactor state using forward Euler method
+    const inflow = math.multiply(math.divide([this.Fs], this.volume), this.Cs_in)[0];
+    const outflow = math.multiply(-1 * math.sum(this.Fs) / this.volume, this.state);
+    const reaction = this.asm.compute_dC(this.state, this.temperature);
+    const transfer = Array(ASM_CONSTANTS.NUM_SPECIES).fill(0.0);
+    transfer[ASM_CONSTANTS.S_O_INDEX] = isNaN(this.kla) ? this.OTR : this._calcOTR(this.state[ASM_CONSTANTS.S_O_INDEX], this.temperature); // calculate OTR if kla is not NaN, otherwise use externaly calculated OTR
+
+    const dC_total = math.multiply(math.add(inflow, outflow, reaction, transfer), time_step)
+    this.state = this._arrayClip2Zero(math.add(this.state, dC_total)); // clip value element-wise to avoid negative concentrations
+    if(DEBUG){
+      assertNoNaN(dC_total, "change in state");
+      assertNoNaN(this.state, "new state");
+    }
+    return this.state;
+  }
+}
+
+class Reactor_PFR extends Reactor {
+  /**
+   * Reactor_PFR class for Plug Flow Reactor.
+   * @param {object} config - Configuration object containing reactor parameters.
+   */
+  constructor(config) {
+    super(config);
+
+    this.length = config.length; // reactor length [m]
+    this.n_x = config.resolution_L; // number of slices
+
+    this.d_x = this.length / this.n_x;
+    this.A = this.volume / this.length; // crosssectional area [m2]
+
+    this.state = Array.from(Array(this.n_x), () => config.initialState.slice());
+    this.extendedState = Array.from(Array(this.n_x + 2*BC_PADDING), () => new Array(ASM_CONSTANTS.NUM_SPECIES).fill(0));
+
+    // initialise extended state
+    this.state.forEach((row, i) => this.extendedState[i+BC_PADDING] = row);
+
+    this.D = 0.0; // axial dispersion [m2 d-1]
+
+    this.D_op = this._makeDoperator();
+    assertNoNaN(this.D_op, "Derivative operator");
+
+    this.D2_op = this._makeD2operator();
+    assertNoNaN(this.D2_op, "Second derivative operator");
+  }
+
+  /**
+   * Setter for axial dispersion.
+   * @param {object} input - Input object (msg) containing payload with dispersion value [m2 d-1].
+   */
+  set setDispersion(input) {
+    this.D = this._constrainDispersion(input.payload);
+  }
+
+  _connectReactor(reactorChild) {
+    if (math.abs(reactorChild.d_x - this.d_x) / this.d_x < 0.025) {
+      this.logger.warn("Significant grid sizing discrepancies between adjacent reactors! Change resolutions to match reactors grid step, or implement boundary value interpolation.");
+    }
+    super._connectReactor(reactorChild);
+  }
+
+  /**
+   * Update the reactor state based on the new time. Performs checks specific to PFR.
+   * @param {number} newTime - New time to update reactor state to, in milliseconds since epoch.
+   */
+  updateState(newTime) {
+    super.updateState(newTime);
+    
+    this.D = this._constrainDispersion(this.D); // constrains D to minimum dispersion, so that local Péclet number is always above 2
+    const Co_D = this.D*this.timeStep/(this.d_x*this.d_x);
+
+    (Co_D >= 0.5) && this.logger.warn(`Courant number (${Co_D}) is too high! Reduce time step size.`);
+
+    if(DEBUG) {
+      console.log("Inlet state max " + math.max(this.state[0]))
+      console.log("Pe total " + this.length*math.sum(this.Fs)/(this.D*this.A));
+      console.log("Pe local " + Pe_local);
+      console.log("Co ad " + math.sum(this.Fs)*this.timeStep/(this.A*this.d_x));
+      console.log("Co D " + Co_D);
+    }
+  }
+
+  /**
+   * Tick the reactor state using explicit finite difference method.
+   * @param {number} time_step - Time step for the simulation [d].
+   * @returns {Array} - New reactor state.
+   */
+  tick(time_step) {
+    this._applyBoundaryConditions();
+
+    const dispersion = math.multiply(this.D / (this.d_x*this.d_x), this.D2_op, this.extendedState);
+    const advection = math.multiply(-1 * math.sum(this.Fs) / (this.A*this.d_x), this.D_op, this.extendedState);
+    const reaction = this.extendedState.map((state_slice) => this.asm.compute_dC(state_slice, this.temperature));
+    const transfer = Array.from(Array(this.n_x+2*BC_PADDING), () => new Array(ASM_CONSTANTS.NUM_SPECIES).fill(0));
+
+    if (isNaN(this.kla)) { // calculate OTR if kla is not NaN, otherwise use externally calculated OTR
+      for (let i = BC_PADDING+1; i < BC_PADDING+this.n_x - 1; i++) {
+        transfer[i][ASM_CONSTANTS.S_O_INDEX] = this.OTR * this.n_x/(this.n_x-2);
+      }
+    } else {
+      for (let i = BC_PADDING+1; i < BC_PADDING+this.n_x - 1; i++) {
+        transfer[i][ASM_CONSTANTS.S_O_INDEX] = this._calcOTR(this.extendedState[i][ASM_CONSTANTS.S_O_INDEX], this.temperature);
+      }
+    }
+
+    const dC_total = math.multiply(math.add(dispersion, advection, reaction, transfer).slice(BC_PADDING, this.n_x+BC_PADDING), time_step);
+
+    const stateNew = math.add(this.state, dC_total);
+
+    if (DEBUG) {
+      assertNoNaN(dispersion, "dispersion");
+      assertNoNaN(advection, "advection");
+      assertNoNaN(reaction, "reaction");
+      assertNoNaN(dC_total, "change in state");
+      assertNoNaN(stateNew, "new state post BC");
+    }
+
+    this.state = this._arrayClip2Zero(stateNew);
+    this.state.forEach((row, i) => this.extendedState[i+BC_PADDING] = row);
+    return stateNew;
+  }
+
+  _updateMeasurement(measurementType, value, position, context) {
+    const grid_pos = Math.round(context.distance / this.config.length * this.n_x);
+
+    // naive approach for reconciling measurements and simulation
+    // could benefit from Kalman filter?
+    switch(measurementType) {
+      case "quantity (oxygen)":
+        this.state[grid_pos][ASM_CONSTANTS.S_O_INDEX] = value;
+        break;
+      case "quantity (ammonium)":
+        this.state[grid_pos][ASM_CONSTANTS.S_NH_INDEX] = value;
+        break;
+      case "quantity (nox)":
+        this.state[grid_pos][ASM_CONSTANTS.S_NO_INDEX] = value;
+        break;
+      default:
+        super._updateMeasurement(measurementType, value, position, context);
+    }
+  }
+
+  /**
+   * Apply boundary conditions to the reactor state.
+   *  for inlet, apply generalised Danckwerts BC, if there is not flow, apply Neumann BC with no flux
+   *  for outlet, apply regular Danckwerts BC (Neumann BC with no flux)
+   */
+  _applyBoundaryConditions() {
+    // Upstream BC
+    if (this.upstreamReactor && this.upstreamReactor.config.reactor_type == "PFR") {
+      // Open boundary, if upstream reactor is PFR
+      this.extendedState.splice(0, BC_PADDING, ...this.upstreamReactor.state.slice(-BC_PADDING)); 
+    } else {
+      if (math.sum(this.Fs) > 0) {
+        // Danckwerts BC
+        const BC_C_in = math.multiply(1 / math.sum(this.Fs), [this.Fs], this.Cs_in)[0];
+        const BC_dispersion_term = this.D*this.A/(math.sum(this.Fs)*this.d_x);
+        this.extendedState[BC_PADDING] = math.multiply(1/(1+BC_dispersion_term), math.add(BC_C_in, math.multiply(BC_dispersion_term, this.extendedState[BC_PADDING+1])));
+        // Numerical boundary condition (first-order accurate)
+        this.extendedState[BC_PADDING-1] = math.add(math.multiply(2, this.extendedState[BC_PADDING]), math.multiply(-2, this.extendedState[BC_PADDING+2]), this.extendedState[BC_PADDING+3]);
+      } else {
+        // Neumann BC (no flux)
+        this.extendedState.fill(this.extendedState[BC_PADDING], 0, BC_PADDING);
+      }
+    }
+
+    // Downstream BC
+    if (this.downstreamReactor && this.downstreamReactor.config.reactor_type == "PFR") {
+      // Open boundary, if downstream reactor is PFR
+      this.extendedState.splice(this.n_x+BC_PADDING, BC_PADDING, ...this.downstreamReactor.state.slice(0, BC_PADDING));
+    } else {
+      // Neumann BC (no flux)
+      this.extendedState.fill(this.extendedState.at(-1-BC_PADDING), BC_PADDING+this.n_x);
+    }
+  }
+
+  /**
+   * Create finite difference first derivative operator.
+   */
+  _makeDoperator() { // create gradient operator
+    const D_size = this.n_x+2*BC_PADDING;
+    const I = math.resize(math.diag(Array(D_size).fill(1/12), -2), [D_size, D_size]);
+    const A = math.resize(math.diag(Array(D_size).fill(-2/3), -1), [D_size, D_size]);
+    const B = math.resize(math.diag(Array(D_size).fill(2/3), 1), [D_size, D_size]);
+    const C = math.resize(math.diag(Array(D_size).fill(-1/12), 2), [D_size, D_size]);
+    const D = math.add(I, A, B, C);
+    // set by BCs elsewhere
+    D.forEach((row, i) => i < BC_PADDING || i >= this.n_x+BC_PADDING ? row.fill(0) : row);
+    return D;
+  }
+
+  /**
+   * Create central finite difference second derivative operator.
+   */
+  _makeD2operator() { // create the central second derivative operator
+    const D_size = this.n_x+2*BC_PADDING;
+    const I = math.diag(Array(D_size).fill(-2), 0);
+    const A = math.resize(math.diag(Array(D_size).fill(1), 1), [D_size, D_size]);
+    const B = math.resize(math.diag(Array(D_size).fill(1), -1), [D_size, D_size]);
+    const D2 = math.add(I, A, B);
+    // set by BCs elsewhere
+    D2.forEach((row, i) => i < BC_PADDING || i >= this.n_x+BC_PADDING ? row.fill(0) : row);
+    return D2;
+  }
+
+  /**
+   * Constrains dispersion so that local Péclet number stays below 2. Needed for stable central differencing method.
+   */
+  _constrainDispersion(D) {
+    const Dmin = math.sum(this.Fs) * this.d_x / (1.999 * this.A);
+    if (D < Dmin) {
+      this.logger.warn(`Local Péclet number too high! Constraining given dispersion (${D}) to minimal value (${Dmin}).`);
+      return Dmin;
+    }
+    return D;
+  }
+}
+
+module.exports = { Reactor_CSTR, Reactor_PFR };

src/utils.js

@@ -0,0 +1,18 @@
+/**
+ * Assert that no NaN values are present in an array.
+ * @param {Array} arr 
+ * @param {string} label 
+ */
+function assertNoNaN(arr, label = "array") {
+    if (Array.isArray(arr)) {
+        for (const el of arr) {
+            assertNoNaN(el, label);
+        }
+    } else {
+        if (Number.isNaN(arr)) {
+            throw new Error(`NaN detected in ${label}!`);
+        }
+    }
+}
+
+module.exports = { assertNoNaN };