diff --git a/.gitignore b/.gitignore
index 37905002381e4d47e31fcc450d81bf593ee71a2d..29cf53ef64723c8b6822c87dc308132d623475bf 100644
--- a/.gitignore
+++ b/.gitignore
@@ -10,3 +10,4 @@ dist
 
 tests/resources/test_output
 
+pydoubletcalc_install
\ No newline at end of file
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 89e92217fd584f5b4897e288639f6b52b0cce93a..e61c1f5a94a752043177e063d0c84559f862de0f 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -17,6 +17,7 @@ Run tests:
     - export PATH=$JAVA_HOME/bin:$PATH
     - export THERMOGIS_JAR=$CI_PROJECT_DIR/resources/thermogis_jar/thermogis-1.7.0-shaded.jar
     - pixi install
+    - pixi run pip install pydoubletcalc --index-url https://gitlab-ci-token:$PYDOUBLETCALC_READ@ci.tno.nl/gitlab/api/v4/projects/19719/packages/pypi/simple
     - pixi run pytest -rx tests/
     - rm -rf /var/lib/apt/lists/*
   only:
diff --git a/README.md b/README.md
index 9153aea91330a90c9bb7ca77092d0a1a2903874e..32e887b2335775a2e616fd331f718cff1af76da5 100644
--- a/README.md
+++ b/README.md
@@ -18,6 +18,27 @@ PyThermoGIS has been designed to be used as a python package you import into you
 It works by creating a python API access to the ThermoGIS techno-economic application, which is written in Java. 
 Because of this dependency you need to Install a Java 17 VM and store the ThermoGIS Jar (a Java executable file) on your computer:
 
+## Installing pydoubletcalc
+
+[pydoubletcalc](https://ci.tno.nl/gitlab/AGS/Geothermal/thermogis/pydoubletcalc) is the python backend we are writing to replace the java backend for pythermogis.
+
+It is not yet public, until then you can install the latest version of the code using the following command:
+
+```bash
+pixi run pip install pydoubletcalc --index-url https://token_name:token_password@ci.tno.nl/gitlab/api/v4/projects/19719/packages/pypi/simple
+```
+replace token_name and token_password with the credentials from a Personal Access Token, with read access.
+(If you don't have a PAT then: Login to gitlab > preferences > Access Tokens)
+
+To update pydoubletcalc you will first need to uninstall it and then reinstall using the above method.
+
+
+## Installing pydoubletcalc locally
+
+```bash
+ pixi add --editable pydoubletcalc@file:C:path\to\pydoubletcalc --pypi
+```
+
 ### 1. Install Java 17 and Download the ThermoGIS JAR
 
 This package requires a Java 17 VM (we recommend using [Amazon Corretto 17](https://docs.aws.amazon.com/corretto/latest/corretto-17-ug/downloads-list.html)) and a [ThermoGIS Jar file](https://ci.tno.nl/gitlab/ags_public/pythermogis/-/blob/main/resources/thermogis_jar/thermogis-1.7.0-shaded.jar?ref_type=heads) (Version >=1.7).
diff --git a/pixi.lock b/pixi.lock
index f59970bbb164a4fed72bc559d1f84602778805f4..41e6e6b2dbad86d663e87745fb000a8b1fccdc3e 100644
--- a/pixi.lock
+++ b/pixi.lock
@@ -252,7 +252,6 @@ environments:
       - pypi: https://files.pythonhosted.org/packages/8f/e9/6a7d025d8da8c4931522922cd706105aa32b3291d1add8c5427cdcd66e63/kiwisolver-1.4.8-cp313-cp313-manylinux_2_17_x86_64.manylinux2014_x86_64.whl
       - pypi: https://files.pythonhosted.org/packages/f5/64/41c4367bcaecbc03ef0d2a3ecee58a7065d0a36ae1aa817fe573a2da66d4/matplotlib-3.10.3-cp313-cp313-manylinux_2_17_x86_64.manylinux2014_x86_64.whl
       - pypi: https://files.pythonhosted.org/packages/d1/80/b9c19f1bb4ac6c5fa6f94a4f278bc68a778473d1814a86a375d7cffa193a/netCDF4-1.7.2-cp313-cp313-manylinux_2_17_x86_64.manylinux2014_x86_64.whl
-      - pypi: https://files.pythonhosted.org/packages/78/f9/690a8600b93c332de3ab4a344a4ac34f00c8f104917061f779db6a918ed6/pathlib-1.0.1-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/02/65/ad2bc85f7377f5cfba5d4466d5474423a3fb7f6a97fd807c06f92dd3e721/plotly-6.0.1-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/c9/58/c3bc54c0fad9a82899e9a2703e04ee6e8eaa76caa90c0689fd1b468a4427/pygridsio-1.0-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/ea/00/d815833441d8c52bf4a6930952e77d3de77d0bf67b3202ccc12dabdae279/pykrige-1.7.2.tar.gz
@@ -264,9 +263,7 @@ environments:
       - pypi: https://files.pythonhosted.org/packages/2a/2f/63d2cacc0e525f8e3398bcf32bd3620385f22cd1600834ec49d7f3597a7b/rioxarray-0.19.0-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/21/f6/4bfb5695d8941e5c570a04d9fcd0d36bce7511b7d78e6e75c8f9791f82d0/scipy-1.16.3-cp313-cp313-manylinux2014_x86_64.manylinux_2_17_x86_64.whl
       - pypi: https://files.pythonhosted.org/packages/d1/a7/5c9cb413e4e2ce52c16be717e94abd40ce91b1f8974624d5d56154c5d40b/shapely-2.1.0-cp313-cp313-manylinux_2_17_x86_64.manylinux2014_x86_64.whl
-      - pypi: https://files.pythonhosted.org/packages/d0/30/dc54f88dd4a2b5dc8a0279bdd7270e735851848b762aeb1c1184ed1f6b14/tqdm-4.67.1-py3-none-any.whl
       - pypi: ./
-      - pypi: C:/Users/knappersfy/work/thermogis/pydoubletcalc
       win-64:
       - conda: https://conda.anaconda.org/conda-forge/win-64/_openmp_mutex-4.5-2_gnu.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/_python_abi3_support-1.0-hd8ed1ab_2.conda
@@ -489,7 +486,6 @@ environments:
       - pypi: https://files.pythonhosted.org/packages/d0/dc/c1abe38c37c071d0fc71c9a474fd0b9ede05d42f5a458d584619cfd2371a/kiwisolver-1.4.8-cp313-cp313-win_amd64.whl
       - pypi: https://files.pythonhosted.org/packages/b1/0f/eed564407bd4d935ffabf561ed31099ed609e19287409a27b6d336848653/matplotlib-3.10.3-cp313-cp313-win_amd64.whl
       - pypi: https://files.pythonhosted.org/packages/66/b5/e04550fd53de57001dbd5a87242da7ff784c80790adc48897977b6ccf891/netCDF4-1.7.2-cp313-cp313-win_amd64.whl
-      - pypi: https://files.pythonhosted.org/packages/78/f9/690a8600b93c332de3ab4a344a4ac34f00c8f104917061f779db6a918ed6/pathlib-1.0.1-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/02/65/ad2bc85f7377f5cfba5d4466d5474423a3fb7f6a97fd807c06f92dd3e721/plotly-6.0.1-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/c9/58/c3bc54c0fad9a82899e9a2703e04ee6e8eaa76caa90c0689fd1b468a4427/pygridsio-1.0-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/ea/00/d815833441d8c52bf4a6930952e77d3de77d0bf67b3202ccc12dabdae279/pykrige-1.7.2.tar.gz
@@ -501,9 +497,7 @@ environments:
       - pypi: https://files.pythonhosted.org/packages/2a/2f/63d2cacc0e525f8e3398bcf32bd3620385f22cd1600834ec49d7f3597a7b/rioxarray-0.19.0-py3-none-any.whl
       - pypi: https://files.pythonhosted.org/packages/cd/01/1204382461fcbfeb05b6161b594f4007e78b6eba9b375382f79153172b4d/scipy-1.16.3-cp313-cp313-win_amd64.whl
       - pypi: https://files.pythonhosted.org/packages/e3/f0/9f8cdf2258d7aed742459cea51c70d184de92f5d2d6f5f7f1ded90a18c31/shapely-2.1.0-cp313-cp313-win_amd64.whl
-      - pypi: https://files.pythonhosted.org/packages/d0/30/dc54f88dd4a2b5dc8a0279bdd7270e735851848b762aeb1c1184ed1f6b14/tqdm-4.67.1-py3-none-any.whl
       - pypi: ./
-      - pypi: C:/Users/knappersfy/work/thermogis/pydoubletcalc
 packages:
 - conda: https://conda.anaconda.org/conda-forge/linux-64/_libgcc_mutex-0.1-conda_forge.tar.bz2
   sha256: fe51de6107f9edc7aa4f786a70f4a883943bc9d39b3bb7307c04c41410990726
@@ -4477,10 +4471,6 @@ packages:
   - pkg:pypi/partd?source=hash-mapping
   size: 20884
   timestamp: 1715026639309
-- pypi: https://files.pythonhosted.org/packages/78/f9/690a8600b93c332de3ab4a344a4ac34f00c8f104917061f779db6a918ed6/pathlib-1.0.1-py3-none-any.whl
-  name: pathlib
-  version: 1.0.1
-  sha256: f35f95ab8b0f59e6d354090350b44a80a80635d22efdedfa84c7ad1cf0a74147
 - conda: https://conda.anaconda.org/conda-forge/noarch/pathspec-0.12.1-pyhd8ed1ab_1.conda
   sha256: 9f64009cdf5b8e529995f18e03665b03f5d07c0b17445b8badef45bde76249ee
   md5: 617f15191456cc6a13db418a275435e5
@@ -4762,20 +4752,6 @@ packages:
   - pkg:pypi/pycparser?source=hash-mapping
   size: 110100
   timestamp: 1733195786147
-- pypi: C:/Users/knappersfy/work/thermogis/pydoubletcalc
-  name: pydoubletcalc
-  version: 0.0.1
-  sha256: af6852e94daff598f382da9a04e74eb8649592424d163653d66a833bc140c065
-  requires_dist:
-  - pandas
-  - matplotlib
-  - pathlib
-  - xarray
-  - tqdm>=4.67.1
-  - scipy>=1.16.2
-  - numba>=0.62.1
-  requires_python: '>=3.13,<3.14'
-  editable: true
 - conda: https://conda.anaconda.org/conda-forge/noarch/pygments-2.19.1-pyhd8ed1ab_0.conda
   sha256: 28a3e3161390a9d23bc02b4419448f8d27679d9e2c250e29849e37749c8de86b
   md5: 232fb4577b6687b2d503ef8e254270c9
@@ -4943,7 +4919,7 @@ packages:
 - pypi: ./
   name: pythermogis
   version: 1.2.0
-  sha256: f25432777c09250f39cf8019497063cc9860feeb7207723ab5ec7783781125d6
+  sha256: f484b4cf939f0f4b0cccabf53e301e4ab4b85f1f4d96b533dbee0e63dfeb1415
   requires_dist:
   - jpype1>=1.5.2,<2
   - xarray==2024.9.0.*
@@ -5640,22 +5616,6 @@ packages:
   - pkg:pypi/tornado?source=hash-mapping
   size: 878044
   timestamp: 1748003914685
-- pypi: https://files.pythonhosted.org/packages/d0/30/dc54f88dd4a2b5dc8a0279bdd7270e735851848b762aeb1c1184ed1f6b14/tqdm-4.67.1-py3-none-any.whl
-  name: tqdm
-  version: 4.67.1
-  sha256: 26445eca388f82e72884e0d580d5464cd801a3ea01e63e5601bdff9ba6a48de2
-  requires_dist:
-  - colorama ; sys_platform == 'win32'
-  - pytest>=6 ; extra == 'dev'
-  - pytest-cov ; extra == 'dev'
-  - pytest-timeout ; extra == 'dev'
-  - pytest-asyncio>=0.24 ; extra == 'dev'
-  - nbval ; extra == 'dev'
-  - requests ; extra == 'discord'
-  - slack-sdk ; extra == 'slack'
-  - requests ; extra == 'telegram'
-  - ipywidgets>=6 ; extra == 'notebook'
-  requires_python: '>=3.7'
 - conda: https://conda.anaconda.org/conda-forge/noarch/twine-6.1.0-pyh29332c3_0.conda
   sha256: c5b373f6512b96324c9607d7d91a76bb53c1056cb1012b4f9c86900c6b7f8898
   md5: d319066fad04e07a0223bf9936090161
diff --git a/pyproject.toml b/pyproject.toml
index be7171d83404499a18b36fc0b8448a88dc52bf70..281dff81bfc1a613ce9e68882bb7bfac7ce7397b 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -36,13 +36,12 @@ markers = [
 [tool.setuptools]
 package-dir = {"" = "src"}
 
-[tool.pixi.project]
+[tool.pixi.workspace]
 channels = ["conda-forge"]
 platforms = ["win-64", "linux-64"]
 
 [tool.pixi.pypi-dependencies]
 pythermogis = { path = ".", editable = true }
-pydoubletcalc = { path = 'C:\Users\knappersfy\work\thermogis\pydoubletcalc', editable = true }
 
 [tool.pixi.tasks]
 test = "PYTHONPATH=src/pythermogis pytest -s tests/"
diff --git a/src/pythermogis/thermogis_classes/doublet.py b/src/pythermogis/thermogis_classes/doublet.py
index 3360a824e7ba0769e48cfff1785798bdf95501d1..50d31633c4abef871eecaac9fed89b9e0e2ed200 100644
--- a/src/pythermogis/thermogis_classes/doublet.py
+++ b/src/pythermogis/thermogis_classes/doublet.py
@@ -1,6 +1,9 @@
 import numpy as np
 import xarray as xr
 from jpype import JClass
+from pythermogis.workflow.utc.doublet import calculate_doublet_performance, DoubletInput, DoubletOutput
+from pythermogis.workflow.utc.utc_properties import UTCConfiguration
+
 
 def simulate_doublet(
     output_data: xr.Dataset,
@@ -102,27 +105,64 @@ def calculate_performance_of_single_location(mask: float, depth: float, thicknes
     if not np.isnan(mask) or temperature < utc_properties.minProdTemp():
         return (mask_value,) * 14
 
-    # Instantiate ThermoGIS doublet
-    doublet = instantiate_thermogis_doublet(utc_properties, rng_seed)
-
-    DoubletInput = JClass("thermogis.calc.utc.doublet.records.DoubletInput")
-    input = DoubletInput(
-            -9999.0, # unknowninput
-            thickness,
-            transmissivity,
-            transmissivity_with_ntg,
-            ntg,
-            depth,
-            porosity,
-            temperature,
-            None, # ates input
+    use_java_backend = True
+    if use_java_backend:
+        doublet = instantiate_thermogis_doublet(utc_properties, rng_seed)
+        JavaDoubletInput = JClass("thermogis.calc.utc.doublet.records.DoubletInput")
+        input = JavaDoubletInput(
+                -9999.0, # unknowninput
+                thickness,
+                transmissivity,
+                transmissivity_with_ntg,
+                ntg,
+                depth,
+                porosity,
+                temperature,
+                None, # ates input
+            )
+        # The Java routine which calculates DoubletPerformance, for more detail on the
+        # simulation inspect the Java source code
+        results_java = doublet.calculateDoubletPerformance(input)
+        results = DoubletOutput(
+            results_java.power(),
+            results_java.hppower(),
+            results_java.capex(),
+            0.0,
+            0.0,
+            results_java.opex(),
+            results_java.utc(),
+            results_java.npv(),
+            results_java.hprod(),
+            results_java.cop(),
+            results_java.cophp(),
+            results_java.pres(),
+            results_java.flow(),
+            results_java.welld(),
+            results_java.breakthrough(),
+            results_java.cooling(),
+            results_java.productionTemp(),
+            results_java.injectionTemp(),
         )
-
-    # The Java routine which calculates DoubletPerformance, for more detail on the simulation inspect the Java source code
-    results = doublet.calculateDoubletPerformance(input)
+    else:
+        # Arrange: instantiate default UTCConfiguration
+        props = UTCConfiguration(segment_length=1)
+        # Create a minimal valid DoubletInput
+        input_data = DoubletInput(
+            unknown_input_value=-9999.0,
+            thickness=thickness,
+            transmissivity=transmissivity,
+            transmissivity_with_ntg=transmissivity_with_ntg,
+            ntg=ntg,
+            depth=depth,
+            porosity=porosity,
+            temperature=temperature,
+        )
+        # Act
+        #perform one simulation to compile all the numba functions before timing
+        results = calculate_doublet_performance(props, input_data)
 
     # If calculation was not successful, return mask value
-    if results.utc() == -9999.0:
+    if results.utc == -9999.0:
         return (mask_value,) * 14
 
     # calculate net-present-value using the utc-cutoffs
@@ -131,24 +171,23 @@ def calculate_performance_of_single_location(mask: float, depth: float, thicknes
     else:
         utc_cut = utc_properties.utcCutoff()
 
-    hprod = results.hprod()
-    npv = 1e-6 * (utc_cut - results.utc()) * 3.6 * hprod * (1 - utc_properties.taxRate())
+    npv = 1e-6 * (utc_cut - results.utc) * 3.6 * results.hprod * (1 - utc_properties.taxRate())
 
     # get values from doublet
-    output_values = {"power": results.power(),
-                     "heat_pump_power": results.hppower(),
-                     "capex": results.capex(),
-                     "opex": results.opex(),
-                     "utc": results.utc(),
+    output_values = {"power": results.power,
+                     "heat_pump_power": results.hppower,
+                     "capex": results.capex,
+                     "opex": results.opex,
+                     "utc": results.utc,
                      "npv": npv,
-                     "hprod": hprod,
-                     "cop": results.cop(),
-                     "cophp": results.cophp(),
-                     "pres": results.pres(),
-                     "flow_rate": results.flow(),
-                     "welld": results.welld(),
-                     "inj_temp": results.injectionTemp(),
-                     "prd_temp": results.productionTemp(),
+                     "hprod": results.hprod,
+                     "cop": results.cop,
+                     "cophp": results.cophp,
+                     "pres": results.pres,
+                     "flow_rate": results.flow,
+                     "welld": results.welld,
+                     "inj_temp": results.injection_temp,
+                     "prd_temp": results.production_temp,
                      }
 
     # Reset doublet variables for next calculation
diff --git a/src/pythermogis/utils/__init__.py b/src/pythermogis/utils/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/src/pythermogis/utils/timer.py b/src/pythermogis/utils/timer.py
new file mode 100644
index 0000000000000000000000000000000000000000..818a092f468448ca36da6a0b37deb72d220ab74c
--- /dev/null
+++ b/src/pythermogis/utils/timer.py
@@ -0,0 +1,28 @@
+import sys
+import timeit
+
+
+def print_time(start, message: str = "", verbose: bool = True):
+    if not verbose:
+        return
+
+    time_taken_seconds = timeit.default_timer() - start
+    time_taken_message = format_time(time_taken_seconds)
+    print(f"{message}{time_taken_message}", flush=True, file=sys.stdout)
+
+    return timeit.default_timer()
+
+
+def format_time(seconds: float) -> str:
+    hours = seconds // 3600
+    minutes = (seconds % 3600) // 60
+    seconds = seconds % 60
+
+    if hours > 0:
+        return f"{hours:.0f}h {minutes:.0f}m {seconds:.0f}s"
+    elif minutes > 0:
+        return f"{minutes:.0f}m {seconds:.0f}s"
+    elif seconds > 10:
+        return f"{seconds:.1f}s"
+    else:
+        return f"{seconds:.3f}s"
diff --git a/src/pythermogis/workflow/utc/cooling_temp.py b/src/pythermogis/workflow/utc/cooling_temp.py
index 1a5be51032d74782c6c2088bbcd791b55978e540..bd62df841e0f1064524435d0ac2239597b7c69bb 100644
--- a/src/pythermogis/workflow/utc/cooling_temp.py
+++ b/src/pythermogis/workflow/utc/cooling_temp.py
@@ -1,4 +1,5 @@
 from pythermogis.workflow.utc.flow import calculate_volumetric_flow
+from pythermogis.workflow.utc.doublet_utils import calc_lifetime
 
 
 def calculate_cooling_temperature(
diff --git a/src/pythermogis/workflow/utc/doublet.py b/src/pythermogis/workflow/utc/doublet.py
index 006a2fe6dd36fee5210e68512bb43533b2d37d87..71debd9c7b4969ff5d2e3fa9da4f07be49c7494d 100644
--- a/src/pythermogis/workflow/utc/doublet.py
+++ b/src/pythermogis/workflow/utc/doublet.py
@@ -1,10 +1,16 @@
+import timeit
 from dataclasses import dataclass
 
-from pythermogis.workflow.utc.utc_properties import UTCConfiguration
-from pythermogis.workflow.utc.doublet_utils import calculate_injection_temp_with_heat_pump, \
-    calc_lifetime
-from pythermogis.workflow.utc.pressure import calculate_max_pressure, optimize_pressure
+import numpy as np
+
+from pythermogis.utils.timer import print_time
 from pythermogis.workflow.utc.cooling_temp import calculate_cooling_temperature
+from pythermogis.workflow.utc.doublet_utils import (
+    calc_lifetime,
+    calculate_injection_temp_with_heat_pump,
+)
+from pythermogis.workflow.utc.pressure import calculate_max_pressure, optimize_pressure
+from pythermogis.workflow.utc.utc_properties import UTCConfiguration
 from pythermogis.workflow.utc.well_distance import optimize_well_distance
 
 EUR_PER_CT_PER_KWH = 0.36
@@ -47,31 +53,31 @@ class DoubletOutput:
     production_temp: float
     injection_temp: float
 
-def calculate_doublet_performance(props: UTCConfiguration, input: DoubletInput) -> DoubletOutput | None:
+def calculate_doublet_performance(props: UTCConfiguration, input: DoubletInput, verbose: bool = False) -> DoubletOutput | None:
+    timer = timeit.default_timer()
+    # determine initial well distance
     well_distance = (
-        (props.optim_dist_well_dist_min + props.optim_dist_well_dist_max) / 2
+        np.mean([props.optim_dist_well_dist_min + props.optim_dist_well_dist_max])
         if props.optim_well_dist
         else props.default_well_distance
     )
 
-    injection_temperature = (
-        max(input.temperature - props.max_cooling_temp_range,
-            props.hp_minimum_injection_temperature)
-        if props.use_heat_pump
-        else max(input.temperature - props.max_cooling_temp_range,
-                 props.dh_return_temp)
-    )
-
-    if props.use_heat_pump and props.calculate_cop and not props.hp_application_mode:
+    # determine injection temperature
+    if props.use_heat_pump:
         injection_temperature = calculate_injection_temp_with_heat_pump(
             input.temperature,
             props.hp_direct_heat_input_temp,
             props.dh_return_temp,
             input.temperature,
             props.max_cooling_temp_range,
-            props.hp_minimum_injection_temperature
+            props.hp_minimum_injection_temperature,
         )
+    else:
+        injection_temperature = max(input.temperature - props.max_cooling_temp_range,
+            props.dh_return_temp)
+    timer = print_time(timer, "injection temperature: ", verbose=verbose)
 
+    # calculate maximum pressure
     drawdown_pressure = calculate_max_pressure(
         props,
         input,
@@ -79,10 +85,11 @@ def calculate_doublet_performance(props: UTCConfiguration, input: DoubletInput)
         well_distance,
         injection_temperature,
     )
-
     if drawdown_pressure == 0:
         return None
+    timer = print_time(timer, "maximum pressure: ", verbose=verbose)
 
+    # cooling temperature and well distance optimization
     if props.optim_well_dist:
         end_temperature_p = (
                 props.optim_dist_cooling_fraction *
@@ -96,6 +103,7 @@ def calculate_doublet_performance(props: UTCConfiguration, input: DoubletInput)
             well_distance,
             injection_temperature,
         )
+    timer = print_time(timer, "cooling temperature: ", verbose=verbose)
 
     if props.optim_well_dist:
         well_distance = optimize_well_distance(
@@ -104,13 +112,16 @@ def calculate_doublet_performance(props: UTCConfiguration, input: DoubletInput)
             drawdown_pressure,
             injection_temperature,
         )
+    timer = print_time(timer, "well distance optimizer: ", verbose=verbose)
 
+    # stimulation capex
     stimulation_capex = (
         0.0
         if (not props.use_stimulation or input.transmissivity_with_ntg > props.stim_kh_max)
         else props.stimulation_capex
     )
 
+    # pressure optimizer
     pressure_results = optimize_pressure(
         props=props,
         input=input,
@@ -119,10 +130,11 @@ def calculate_doublet_performance(props: UTCConfiguration, input: DoubletInput)
         injection_temp=injection_temperature,
         stimulation_capex=stimulation_capex,
     )
-
+    timer = print_time(timer, "pressure optimizer: ", verbose=verbose)
     if pressure_results is None:
         return None
 
+    # everything underneath here is fast, no point in optimizing
     heat_power_per_doublet = pressure_results.heat_power_per_doublet
     flowrate = pressure_results.flowrate
     discounted_heat_produced_p = pressure_results.discounted_heat_produced
@@ -179,7 +191,7 @@ def calculate_doublet_performance(props: UTCConfiguration, input: DoubletInput)
             * (1 - props.tax_rate)
     )
     opex_first_prod_year = total_opex_ts[props.drilling_time]
-    hp_cop = 3.0
+    print_time(timer, "economics: ", verbose=verbose)
 
     return DoubletOutput(
         power=heat_power_per_doublet,
@@ -192,7 +204,7 @@ def calculate_doublet_performance(props: UTCConfiguration, input: DoubletInput)
         npv=npv,
         hprod=discounted_heat_produced_p,
         cop=cop,
-        cophp=hp_cop,
+        cophp=pressure_results.heat_pump_cop,
         pres=drawdown_pressure / 1e5,
         flow=flowrate,
         welld=well_distance,
diff --git a/src/pythermogis/workflow/utc/doublet_utils.py b/src/pythermogis/workflow/utc/doublet_utils.py
index 7eba85a8ae08db55971ea9e3c9e912bb84c30d9a..12ac94c9449a8f04cb81e65d1745e06d6cd8634b 100644
--- a/src/pythermogis/workflow/utc/doublet_utils.py
+++ b/src/pythermogis/workflow/utc/doublet_utils.py
@@ -59,7 +59,7 @@ def get_cop_carnot(eta: float, Tout: float, Tin: float) -> float:
 
     return eta * (Tcond + TKELVIN) / (Tcond - Tevap)
 
-
+@njit
 def calc_lifetime(
     well_distance: float,
     thickness: float,
@@ -113,7 +113,7 @@ def calc_lifetime(
     return Eseg / Eflowseg
 
 
-
+@njit
 def get_along_hole_length(
     true_vertical_depth: float,
     well_distance: float,
diff --git a/src/pythermogis/workflow/utc/doubletcalc.py b/src/pythermogis/workflow/utc/doubletcalc.py
index d83ed5529ebe621013ac7d1095201d77e024baed..548a168bccf39ca60972f766e72cd4aa745f1325 100644
--- a/src/pythermogis/workflow/utc/doubletcalc.py
+++ b/src/pythermogis/workflow/utc/doubletcalc.py
@@ -1,12 +1,10 @@
 import math
-
 from dataclasses import dataclass
 from pythermogis.workflow.utc.utc_properties import UTCConfiguration
 from pydoubletcalc import Aquifer, Doublet, Well, WellPipeSegment
 from pythermogis.workflow.utc.water import get_salinity
 from pythermogis.workflow.utc.rock import get_geothermal_gradient
-
-
+from numba import njit
 INCH_SI = 0.0254
 
 @dataclass
@@ -24,7 +22,7 @@ def doubletcalc(
         drawdown_pressure: float,
         well_distance: float,
         injection_temp: float,
-) -> Doublet1DResults:
+) -> Doublet1DResults | None:
     aquifer = Aquifer(
         permeability=input.permeability,
         porosity=input.porosity,
@@ -42,12 +40,19 @@ def doubletcalc(
     injector = Well(
         aquifer=aquifer,
         well_type="injector",
-        pipe_segments=[WellPipeSegment(
-            ah_depth=get_along_hole_length(input.depth, well_distance, props.well_curv_scaling, props.max_tvd_stepout_factor),
-            tv_depth=input.depth,
-            inner_diameter=props.inner_diameter * INCH_SI,
-            roughness=props.roughness * 1e-3 * INCH_SI
-        )],
+        pipe_segments=[
+            WellPipeSegment(
+                ah_depth=get_along_hole_length(
+                    input.depth,
+                    well_distance,
+                    props.well_curv_scaling,
+                    props.max_tvd_stepout_factor,
+                ),
+                tv_depth=input.depth,
+                inner_diameter=props.inner_diameter * INCH_SI,
+                roughness=props.roughness * 1e-3 * INCH_SI,
+            )
+        ],
         aquifer_top_depth=input.depth,
         pipe_scaling=0,
         target_segment_length=props.segment_length,
@@ -61,6 +66,7 @@ def doubletcalc(
         viscosity_mode=props.viscosity_mode,
         heat_exchanger_exit_temperature=injection_temp,
     )
+
     producer = Well(
         aquifer=aquifer,
         well_type="producer",
@@ -94,7 +100,13 @@ def doubletcalc(
         viscosity_mode=props.viscosity_mode,
         heat_exchanger_exit_temperature=injection_temp,
     )
-    doublet.simulate()
+
+    # don't calculate lifetime as it is unneded for Doublet1DResults
+    doublet.calculate_mass_rate()
+    if not doublet.simulation_success:
+        return None
+    doublet.calculate_geothermal_power()
+    doublet.calculate_power()
 
     return Doublet1DResults(
         geothermal_powers=doublet.geothermal_power,
@@ -124,13 +136,13 @@ def get_total_skin_production(props, input):
 def get_pump_production_depth(props, depth: float) -> float:
     return min(props.pump_depth, depth / 2)
 
+@njit
 def get_along_hole_length(
     true_vertical_depth: float,
     well_distance: float,
     curve_scaling_factor: float,
     max_true_vertical_depth_stepout_factor: float,
 ) -> float:
-
     if curve_scaling_factor == 0:
         # Vertical well: along-hole length equals TVD
         return true_vertical_depth
@@ -143,6 +155,7 @@ def get_along_hole_length(
 
     stepout -= horizontal_distance
 
-    oblique_distance = math.sqrt(stepout**2 + true_vertical_depth**2) * curve_scaling_factor
-
+    oblique_distance = (
+        math.sqrt(stepout**2 + true_vertical_depth**2) * curve_scaling_factor
+    )
     return oblique_distance + horizontal_distance
\ No newline at end of file
diff --git a/src/pythermogis/workflow/utc/flow.py b/src/pythermogis/workflow/utc/flow.py
index 2c5d3a4f3acd2fee54bceac986cde738f2c79af1..97411a17ea82fc6a67b2188f8ab1233cc6e382f6 100644
--- a/src/pythermogis/workflow/utc/flow.py
+++ b/src/pythermogis/workflow/utc/flow.py
@@ -1,13 +1,16 @@
 from dataclasses import dataclass
+import numpy as np
 
 from pythermogis.workflow.utc.heatpump import calculate_heat_pump_performance
 from pythermogis.workflow.utc.doublet_utils import get_orc_efficiency
 from pythermogis.workflow.utc.doubletcalc import doubletcalc
 
+
 @dataclass
 class VolumetricFlowResults:
     hp_cop: float
     hp_added_power: float
+    hp_elec_consumption: float
     heat_power_per_doublet: float
     cop: float
     flowrate: float
@@ -21,85 +24,74 @@ def calculate_volumetric_flow(
     input_data,
     original_pressure: float,
     well_distance: float,
-    injection_temp: float
+    injection_temp: float,
 ):
-    STEPS = [0, 1, -1, 2, -2, 3, -3, 4, -4, 5, -5, 6, -6]
-    BAR_SI = 1e5  # Units.BAR_SI
-    hp_cop = 3.0
-    hp_added_power = 0.0
-
-    d1d_results = None
-    iter_index = 0
-
-    while d1d_results is None and iter_index < len(STEPS):
-        drawdown_pressure = original_pressure + STEPS[iter_index] * BAR_SI
-
+    for step in [0, 1e5, -1e5, 2e5, -2e5, 3e5, -3e5]:
         d1d_results = doubletcalc(
-            props, input_data, drawdown_pressure, well_distance, injection_temp
+            props,
+            input_data,
+            original_pressure + step,
+            well_distance,
+            injection_temp
         )
-        iter_index += 1
+        if d1d_results is not None:
+            break
 
+    # if no valid output results return None
     if d1d_results is None:
         return None
 
-    geothermal_powers = d1d_results.geothermal_powers
-    cop = d1d_results.cop
-    flowrate = d1d_results.flowrate
-    pump_power_required = d1d_results.pump_power_required
-    production_temp = d1d_results.production_temp
-    heat_power_produced = d1d_results.heat_power_produced
-
-    he2 = props.heat_exchanger_efficiency
-
     if props.use_orc:
-        Ts = props.heat_exchanger_basetemp
-        he2 = get_orc_efficiency(
-            production_temp, Ts, props.heat_exchanger_efficiency
+        heat_exchanger_efficiency = get_orc_efficiency(
+            d1d_results.production_temp,
+            props.heat_exchanger_basetemp,
+            props.heat_exchanger_efficiency,
         )
-
-    heat_power_per_doublet = max(1e-6, geothermal_powers * he2)
-
-    for i in range(len(heat_power_produced)):
-        heat_power_produced[i] = max(1e-6, heat_power_produced[i] * he2)
-
-    ignore_subsurface = (
-        props.well_cost_scaling
-        + props.well_cost_const
-        + props.well_cost_z
-        + props.well_cost_z2
-    ) < 1e-3
-
-    if ignore_subsurface:
-        cop = 1e4
-        pump_power_required = 0.0
-
-    power_consumption = ((heat_power_per_doublet / he2) / cop) + \
-                        (heat_power_per_doublet * props.heat_exchanger_parasitic)
-
+    else:
+        heat_exchanger_efficiency = props.heat_exchanger_efficiency
+
+    heat_power_produced = np.clip(
+        np.array(d1d_results.heat_power_produced) * heat_exchanger_efficiency,
+        1e-6,
+        None,
+    )
+
+    # Calculate cop
+    heat_power_per_doublet = max(
+        1e-6, d1d_results.geothermal_powers * heat_exchanger_efficiency
+    )
+    power_consumption = (
+        (heat_power_per_doublet / heat_exchanger_efficiency) / d1d_results.cop
+    ) + (heat_power_per_doublet * props.heat_exchanger_parasitic)
     if props.use_orc:
         heat_power_per_doublet -= power_consumption
-
     cop = heat_power_per_doublet / power_consumption
 
+    # Calculate heat pump performance
     if props.use_heat_pump:
         hp_results = calculate_heat_pump_performance(
             props,
             input_data,
-            production_temp,
+            d1d_results.production_temp,
             injection_temp,
-            flowrate
+            d1d_results.flowrate,
         )
-
         hp_cop = hp_results.hp_cop
         hp_added_power = hp_results.hp_added_power
+        hp_elec_consumption = hp_results.hp_elec_consumption
+    else:
+        hp_cop = 0.0
+        hp_added_power = 0.0
+        hp_elec_consumption = 0.0
 
     return VolumetricFlowResults(
         hp_cop,
         hp_added_power,
+        hp_elec_consumption,
         heat_power_per_doublet,
         cop,
-        flowrate,
-        pump_power_required,
-        production_temp,
-        heat_power_produced
+        d1d_results.flowrate,
+        d1d_results.pump_power_required,
+        d1d_results.production_temp,
+        heat_power_produced,
     )
\ No newline at end of file
diff --git a/src/pythermogis/workflow/utc/heatpump.py b/src/pythermogis/workflow/utc/heatpump.py
index 7e1b6b11f6560a629774d5b861a4966c832dfedd..1043160a8663e271d05d4d559e36e5028f65e304 100644
--- a/src/pythermogis/workflow/utc/heatpump.py
+++ b/src/pythermogis/workflow/utc/heatpump.py
@@ -7,6 +7,7 @@ from pythermogis.workflow.utc.water import get_hydrostatic_pressure, get_salinit
 class HeatPumpPerformanceResults:
     hp_cop: float
     hp_added_power: float
+    hp_elec_consumption: float
 
 
 def calculate_heat_pump_performance(
@@ -35,28 +36,31 @@ def calculate_heat_pump_performance(
         hydrostatic_pressure, props.surface_temperature, salinity
     )
 
-    cp_water = heat_capacity(
-        hydrostatic_pressure, props.surface_temperature, salinity
-    )
-
-    if props.hp_calc_cop:
-        Tout = props.hp_direct_heat_input_temp
-        Tin = min(injection_temp, props.dh_return_temp)
-        hp_cop = get_cop_carnot(ETA_CARNOT, Tout, Tin)
-    else:
-        hp_cop = props.hp_capex
+    cp_water = heat_capacity(hydrostatic_pressure, props.surface_temperature, salinity)
 
+    # calculate the power used by the heat pump
     waste_temp = min(heat_pump_start_temp, props.dh_return_temp)
     delta_temp = max(0.0, waste_temp - injection_temp)
-
     hp_added_power = (flowrate / 3600.0) * delta_temp * rho_water * cp_water * 1e-6
-
     if hp_added_power == 0:
-        hp_cop = 0.0
+        return HeatPumpPerformanceResults(
+            hp_cop=0,
+            hp_added_power=0,
+            hp_elec_consumption=0,
+        )
+
+    # calculate heat pump coefficient of performance
+    Tout = props.hp_direct_heat_input_temp
+    Tin = min(injection_temp, props.dh_return_temp)
+    hp_cop = get_cop_carnot(ETA_CARNOT, Tout, Tin)
+
+    # calculate hp_elec_consumption
+    hp_elec_consumption = hp_added_power / (hp_cop -1)
 
     return HeatPumpPerformanceResults(
         hp_cop=hp_cop,
         hp_added_power=hp_added_power,
+        hp_elec_consumption=hp_elec_consumption
     )
 
 
@@ -64,11 +68,10 @@ def calculate_heat_pump_start_temp(props, production_temp: float, injection_temp
     """
     Python version of calculateHeatPumpStartTemp(...)
     """
-
     delta_temp_geothermal = production_temp - injection_temp
     delta_temp_for_hex = 0.0
 
-    if props.hp_application_mode and production_temp > props.dh_return_temp:
+    if production_temp > props.dh_return_temp:
         delta_temp_for_hex = production_temp - props.dh_return_temp
         delta_temp_for_hex = max(0.0, min(delta_temp_for_hex, delta_temp_geothermal))
 
diff --git a/src/pythermogis/workflow/utc/pressure.py b/src/pythermogis/workflow/utc/pressure.py
index 2303f67a5efa18f9e141549c6deee4b1546dbe4a..40179fc41694b14df50264fcc128b084941eb74c 100644
--- a/src/pythermogis/workflow/utc/pressure.py
+++ b/src/pythermogis/workflow/utc/pressure.py
@@ -1,7 +1,9 @@
 
 from dataclasses import dataclass
-from pythermogis.workflow.utc.flow import calculate_volumetric_flow
+
 from pythermogis.workflow.utc.economics import calculate_economics
+from pythermogis.workflow.utc.flow import calculate_volumetric_flow
+
 
 def calculate_max_pressure(
     props,
@@ -10,16 +12,15 @@ def calculate_max_pressure(
     well_distance: float,
     injection_temp: float
 ) -> float:
-
     if use_olsthoorn_max_pressure:
         max_pres = 0.2 * 0.1 * input_data.depth * 100000
     else:
         max_pres = input_data.depth * (0.135 - props.hy_gradient) * 100000
 
-    pres = min(props.max_pump * 1e5, max_pres)
+    pressure = min(props.max_pump * 1e5, max_pres)
 
     results = calculate_volumetric_flow(
-        props, input_data, pres, well_distance, injection_temp
+        props, input_data, pressure, well_distance, injection_temp
     )
 
     if results is None:
@@ -27,51 +28,52 @@ def calculate_max_pressure(
         iter_count = 0
 
         pres_min = 1e5
-        pres_max = pres
+        pres_max = pressure
 
         while iter_count < 1000 and abs(pres_max - pres_min) > pres_tol:
             iter_count += 1
-            pres = 0.5 * (pres_min + pres_max)
+            pressure = 0.5 * (pres_min + pres_max)
 
             results = calculate_volumetric_flow(
-                props, input_data, pres, well_distance, injection_temp
+                props, input_data, pressure, well_distance, injection_temp
             )
 
             if results is not None:
-                pres_min = pres
+                pres_min = pressure
             else:
-                pres_max = pres
+                pres_max = pressure
 
         if results is None:
-            pres -= pres_tol
+            pressure -= pres_tol
             results = calculate_volumetric_flow(
-                props, input_data, pres, well_distance, injection_temp
+                props, input_data, pressure, well_distance, injection_temp
             )
 
         if results is None or iter_count >= 1000:
             return 0.0
 
-        return pres
+        return pressure
 
     if results.heat_power_per_doublet < 0 and not props.is_ates:
-        pres /= 2.0
+        pressure /= 2.0
 
         results = calculate_volumetric_flow(
-            props, input_data, pres, well_distance, injection_temp
+            props, input_data, pressure, well_distance, injection_temp
         )
 
         if results.heatPowerPerDoublet() < 0:
             return 0.0
         else:
-            pres *= 2.0
+            pressure *= 2.0
 
-    return pres
+    return pressure
 
 @dataclass
 class PressureOptimizerResults:
     drawdown_pressure: float
     flowrate: float
     heat_power_per_doublet: float
+    heat_pump_cop: float
     cop: float
     utc: float
     sum_capex: float
@@ -104,9 +106,7 @@ def optimize_pressure(
     if flow_results is None:
         return None
 
-
     if flow_results.production_temp > injection_temp:
-
         iter_count = 0
         pres_min = 1e5 * props.min_pump
         pres_max = drawdown_pressure
@@ -151,11 +151,9 @@ def optimize_pressure(
                     pres_max = pres_max * 1.01
 
     if flow_results is not None and flow_results.flowrate > props.max_flow:
-
         iter_count = 0
         pres_min = 0.0
         pres_max = drawdown_pressure
-
         while iter_count < 1000 and abs(pres_max - pres_min) > pres_tol:
             iter_count += 1
             pres = 0.5 * (pres_min + pres_max)
@@ -177,29 +175,20 @@ def optimize_pressure(
         # Final adjustment if flow still too high
         if flow_results.flowrate > props.max_flow:
             drawdown_pressure -= pres_tol
+            flow_results = calculate_volumetric_flow(
+                props=props,
+                input_data=input,
+                original_pressure=drawdown_pressure,
+                well_distance=well_distance,
+                injection_temp=injection_temp,
+            )
 
-
-    flow_results = calculate_volumetric_flow(
-        props=props,
-        input_data=input,
-        original_pressure=drawdown_pressure,
-        well_distance=well_distance,
-        injection_temp=injection_temp,
-    )
     if flow_results is None or flow_results.flowrate > props.max_flow:
         return None
 
-
-    hp_added_power = flow_results.hp_added_power
-    if hp_added_power > 0:
-        hp_elec_consumption = hp_added_power / (flow_results.hp_cop - 1.0)
-    else:
-        hp_elec_consumption = 0.0
-
-
     cop = (
-        (flow_results.heat_power_per_doublet + hp_elec_consumption) /
-        (flow_results.heat_power_per_doublet / flow_results.cop + hp_elec_consumption)
+        (flow_results.heat_power_per_doublet + flow_results.hp_elec_consumption) /
+        (flow_results.heat_power_per_doublet / flow_results.cop + flow_results.hp_elec_consumption)
     )
 
     econ = calculate_economics(
@@ -220,6 +209,7 @@ def optimize_pressure(
         drawdown_pressure=drawdown_pressure,
         flowrate=flow_results.flowrate,
         heat_power_per_doublet=flow_results.heat_power_per_doublet,
+        heat_pump_cop=flow_results.hp_cop,
         cop=cop,
         utc=econ.utc.utc,
         sum_capex=econ.capex.sum_capex,
diff --git a/src/pythermogis/workflow/utc/well_distance.py b/src/pythermogis/workflow/utc/well_distance.py
index 01249f6deafa58f4bc181ef03b827c31a1cf816d..0b9be4e3de48ce5609b9c1586d4c2377809fb9d4 100644
--- a/src/pythermogis/workflow/utc/well_distance.py
+++ b/src/pythermogis/workflow/utc/well_distance.py
@@ -1,42 +1,58 @@
+import numpy as np
+from scipy.optimize import brentq
+
 from pythermogis.workflow.utc.doublet_utils import calc_lifetime
 from pythermogis.workflow.utc.flow import calculate_volumetric_flow
 
 
-def optimize_well_distance(
+def optimize_well_distance_original(
     props,
     input,
     drawdown_pressure: float,
     injection_temp: float,
 ) -> float:
+    """
+    The original well distance optimization algorithm. This is kept for testing purposes
+    but has been replaced with the scipy implementation of the brenth algorithm which
+    does the same thing more accurately.
+
+    Parameters
+    ----------
+    props
+    input
+    drawdown_pressure
+    injection_temp
+
+    Returns
+    -------
+
+    """
 
     dist_min = props.optim_dist_well_dist_min
     dist_max = props.optim_dist_well_dist_max
-    dist = 0.5 * (dist_min + dist_max)
+    well_distance = np.mean([dist_min, dist_max])
 
     for iter_count in range(1000):
         if abs(dist_max - dist_min) <= 10.0:
-            break
-
-        dist = 0.5 * (dist_min + dist_max)
+            return well_distance
+        well_distance = np.mean([dist_min, dist_max])
 
         # --- Compute flow for this distance ---
         results = calculate_volumetric_flow(
             props=props,
             input_data=input,
             original_pressure=drawdown_pressure,
-            well_distance=dist,
+            well_distance=well_distance,
             injection_temp=injection_temp,
         )
 
-        flowrate = min(results.flowrate, props.max_flow)
-
         # --- Compute lifetime for this distance ---
         lifetime = calc_lifetime(
-            well_distance=dist,
+            well_distance=well_distance,
             thickness=input.thickness * input.ntg,
             delta_temp_fraction=props.optim_dist_cooling_fraction,
             porosity=input.porosity,
-            flowrate=flowrate,
+            flowrate=min(results.flowrate, props.max_flow),
             depth=input.depth,
             reservoir_temp=input.temperature,
             salinity_surface=props.salinity_surface,
@@ -47,12 +63,60 @@ def optimize_well_distance(
 
         # --- Bisection rule ---
         if lifetime < props.optim_dist_lifetime:
-            dist_min = dist
+            dist_min = well_distance
         else:
-            dist_max = dist
-
+            dist_max = well_distance
     # If no convergence in 1000 iterations
     else:
-        print(f"WARNING: Well distance optimization failed to converge. Final dist={dist}")
+        print(
+            f"WARNING: Well distance optimization failed to converge. Final dist={well_distance}"
+        )
+
+    return well_distance
+
 
-    return dist
+def f1(
+    well_distance, props, input, drawdown_pressure: float, injection_temp: float
+) -> float:
+    # --- Compute flow for this distance ---
+    results = calculate_volumetric_flow(
+        props=props,
+        input_data=input,
+        original_pressure=drawdown_pressure,
+        well_distance=well_distance,
+        injection_temp=injection_temp,
+    )
+
+    # --- Compute lifetime for this distance ---
+    return calc_lifetime(
+        well_distance=well_distance,
+        thickness=input.thickness * input.ntg,
+        delta_temp_fraction=props.optim_dist_cooling_fraction,
+        porosity=input.porosity,
+        flowrate=min(results.flowrate, props.max_flow),
+        depth=input.depth,
+        reservoir_temp=input.temperature,
+        salinity_surface=props.salinity_surface,
+        salinity_gradient=props.salinity_gradient,
+        cp_rock=props.optim_dist_cp_rock,
+        rho_rock=props.optim_dist_rho_rock,
+    ) - props.optim_dist_lifetime
+
+
+def optimize_well_distance(
+    props,
+    input,
+    drawdown_pressure: float,
+    injection_temp: float,
+) -> float:
+    # find the well distance between the min and max which comes closest to the optimal
+    # doublet lifetime, as defined in props.
+    well_distance = brentq(
+        f1,
+        props.optim_dist_well_dist_min,
+        props.optim_dist_well_dist_max,
+        xtol=10,
+        maxiter=100,
+        args=(props, input, drawdown_pressure, injection_temp),
+    )
+    return well_distance
\ No newline at end of file
diff --git a/tests/test_dask_parralelization.py b/tests/test_dask_parralelization.py
index 9573cf868cd31705b5b1368b2fe92351353e44a1..138d424cc81b2ea255204e1f79b7fd3294b0c6ab 100644
--- a/tests/test_dask_parralelization.py
+++ b/tests/test_dask_parralelization.py
@@ -2,15 +2,17 @@ from os import path
 from pathlib import Path
 
 import numpy as np
+import pytest
 import xarray as xr
 
 from pythermogis import auto_chunk_dataset, assess_optimal_chunk_size
 
+# @pytest.mark.skip("Useful for making a plot, but not needed for pipeline testing")
 def test_dask_parralelization():
     output_data_path = Path(path.dirname(__file__), "resources") / "test_output" / "parallelization"
     output_data_path.mkdir(parents=True, exist_ok=True)
 
-    assess_optimal_chunk_size(n_simulations = 100, chunk_step_size=50, plot_outfile = output_data_path / "parallelization.png")
+    assess_optimal_chunk_size(n_simulations = 100, chunk_step_size=10, plot_outfile = output_data_path / "parallelization.png")
 
     assert output_data_path.exists()
 
diff --git a/tests/test_doc_examples.py b/tests/test_doc_examples.py
index faa28fd92f430322920f52bcfc6dd84e5a46defa..23a2b626576df0e0760ff1fd207c90d27d771a0d 100644
--- a/tests/test_doc_examples.py
+++ b/tests/test_doc_examples.py
@@ -280,7 +280,7 @@ def test_example8():
     output_data_path.mkdir(parents=True, exist_ok=True)
 
     # generate simulation samples across desired reservoir properties
-    Nsamples = 100
+    Nsamples = 1000
     thickness_samples = np.random.uniform(low=150, high=300, size=Nsamples)
     porosity_samples = np.random.uniform(low=0.5, high=0.8, size=Nsamples)
     ntg_samples = np.random.uniform(low=0.25, high=0.5, size=Nsamples)
@@ -298,7 +298,7 @@ def test_example8():
     )
 
     # For every sample, run a doublet simulation store the output values
-    simulations = calculate_doublet_performance(reservoir_properties, print_execution_duration=True)
+    simulations = calculate_doublet_performance(reservoir_properties, chunk_size=100, print_execution_duration=True)
 
     # post process the samples to calculate the percentiles of each variable
     percentiles = np.arange(1,99)
diff --git a/tests/test_doublet_speed.py b/tests/test_doublet_speed.py
index 36be0e513366e375dc6c38fae7c426740c70226f..c93073de07bd61d9c2028e0a9ba9916b4b07235c 100644
--- a/tests/test_doublet_speed.py
+++ b/tests/test_doublet_speed.py
@@ -27,15 +27,14 @@ class PyThermoGIS(TestCase):
         non_nan_data = data[~np.isnan(data)]
         print(f"Running Performance calculation for ROSL_ROSLU, Pvalues: {p_values}, The number of Non Nan Cells: {len(non_nan_data)}")
 
-        # Run calculation across all dimensions of input_grids, and all provided P_values: The Equivalent calculation with 10 cores takes 43 seconds in the Java code
+        # Run calculation across all dimensions of input_grids, and all provided P_values
         start = timeit.default_timer()
-        calculate_doublet_performance_stochastic(input_grids, chunk_size=100, rng_seed=123, p_values=p_values)
+        calculate_doublet_performance_stochastic(input_grids, chunk_size=10, rng_seed=123, p_values=p_values)
         time_elapsed = timeit.default_timer() - start
         print(f"Python calculation took {time_elapsed:.1f} seconds.")
 
-
     def read_input_grids(self):
-        new_cellsize=5000 # in m
+        new_cellsize=10e3 # in m
         input_grids = resample_grid(read_grid(self.test_files_out_path / "ROSL_ROSLU__thick.zmap"), new_cellsize=new_cellsize).to_dataset(name="thickness_mean")
         input_grids["thickness_sd"] = resample_grid(read_grid(self.test_files_out_path / "ROSL_ROSLU__thick_sd.zmap"), new_cellsize=new_cellsize)
         input_grids["ntg"] = resample_grid(read_grid(self.test_files_out_path / "ROSL_ROSLU__ntg.zmap"), new_cellsize=new_cellsize)
diff --git a/tests/utc/test_calculate_volumetric_flow.py b/tests/utc/test_calculate_volumetric_flow.py
index afed8c710d7615c16619fc410187a4915c124f53..a6b93e39f0f856a9fca1ae723f758409c485cf3b 100644
--- a/tests/utc/test_calculate_volumetric_flow.py
+++ b/tests/utc/test_calculate_volumetric_flow.py
@@ -1,8 +1,9 @@
-from pythermogis.workflow.utc.utc_properties import UTCConfiguration
+import numpy as np
+
 from pythermogis.workflow.utc.doublet import DoubletInput
 from pythermogis.workflow.utc.flow import calculate_volumetric_flow
+from pythermogis.workflow.utc.utc_properties import UTCConfiguration
 
-import numpy as np
 
 def test_calculate_volumetric_flow():
     props = UTCConfiguration()
@@ -20,7 +21,7 @@ def test_calculate_volumetric_flow():
 
     flow_results = calculate_volumetric_flow(props, input_data, 600000, 1550, 40)
 
-    assert np.isclose(flow_results.hp_cop, 3.0, rtol=0.01)
+    assert np.isclose(flow_results.hp_cop, 0.0, rtol=0.01)
     assert np.isclose(flow_results.hp_added_power, 0.0, rtol=0.01)
     assert np.isclose(flow_results.heat_power_per_doublet, 0.027926914290870505, rtol=0.01)
     assert np.isclose(flow_results.cop, 5.396068601571214, rtol=0.01)
diff --git a/tests/utc/test_doublet.py b/tests/utc/test_doublet.py
index 61eaed357c7fa8b47a275323abe7557c62255d30..785a5f405fdc1311d6877697ee06f9009a1d017a 100644
--- a/tests/utc/test_doublet.py
+++ b/tests/utc/test_doublet.py
@@ -1,15 +1,59 @@
+import timeit
+
 import numpy as np
 
-from pythermogis.workflow.utc.doublet import calculate_doublet_performance, DoubletInput, DoubletOutput
+from pythermogis.workflow.utc.doublet import DoubletInput, calculate_doublet_performance
 from pythermogis.workflow.utc.utc_properties import UTCConfiguration
 
 
+def test_calculate_doublet_performance_precise():
+    # Arrange: instantiate default UTCConfiguration
+    props = UTCConfiguration(
+        viscosity_mode="kestin",
+        dh_return_temp=40,
+    )
 
-def test_calculate_doublet_performance():
+    # Create a minimal valid DoubletInput
+    input_data = DoubletInput(
+        unknown_input_value=-999.0,
+        thickness=100.0,
+        transmissivity=17500.0,
+        transmissivity_with_ntg=0.0,
+        ntg=1.0,
+        depth=2000,
+        porosity=0.0,
+        temperature=76.0,
+    )
+
+    # Act
+    # perform one simulation to compile all the numba functions before timing
+    calculate_doublet_performance(props, input_data)
+
+    start = timeit.default_timer()
+    result = calculate_doublet_performance(props, input_data, verbose=True)
+    time_elapsed = timeit.default_timer() - start
+
+    print(f"1 simulation took: {time_elapsed:.1f} seconds\n"
+          f"{1/time_elapsed:.1f} simulations per second")
+
+
+    # Assert
+    rtol = 0.005 # accurate to 0.5%
+    assert np.isclose(result.flow, 227.2757568359375, rtol=rtol)
+    assert np.isclose(result.pres, 60, rtol=rtol)
+    assert np.isclose(result.utc, 6.616096470753937, rtol=rtol)
+    assert np.isclose(result.welld, 1159.17968, rtol=rtol)
+    assert np.isclose(result.power, 8.636903762817383, rtol=rtol)
+    assert np.isclose(result.cop, 13.627557754516602, rtol=rtol)
+    assert np.isclose(result.var_opex, -7.510325908660889, rtol=rtol)
+    assert np.isclose(result.fixed_opex, -11.227973937988281, rtol=rtol)
+
+def test_calculate_doublet_performance_approximate():
     # Arrange: instantiate default UTCConfiguration
     props = UTCConfiguration(
         viscosity_mode="kestin",
         dh_return_temp=40,
+        segment_length=1.0,
     )
 
     # Create a minimal valid DoubletInput
@@ -25,17 +69,28 @@ def test_calculate_doublet_performance():
     )
 
     # Act
+    #perform one simulation to compile all the numba functions before timing
     result = calculate_doublet_performance(props, input_data)
 
-    # Assert
-    assert np.isclose(result.flow, 227.2757568359375, atol=1.5)
-    assert np.isclose(result.pres, 60, atol=0.001)
-    assert np.isclose(result.utc, 6.616096470753937, atol=0.001)
-    assert np.isclose(result.welld, 1159.17968, atol=0.001)
-    assert np.isclose(result.power, 8.636903762817383, atol=0.01)
-    assert np.isclose(result.cop, 13.627557754516602, atol=0.01)
-    assert np.isclose(result.var_opex, -7.510325908660889, atol=0.01)
-    assert np.isclose(result.fixed_opex, -11.227973937988281, atol=0.01)
+    start = timeit.default_timer()
+    n_sims = 500
+    for _ in range(n_sims):
+        result = calculate_doublet_performance(props, input_data)
+    time_elapsed = timeit.default_timer() - start
 
+    print(f"{n_sims} simulations took: {time_elapsed:.1f} seconds\n"
+          f"{n_sims/time_elapsed:.1f} simulations per second")
+
+
+    # Assert
+    rtol = 0.01 # accurate to 1%
+    assert np.isclose(result.flow, 227.2757568359375, rtol=rtol)
+    assert np.isclose(result.pres, 60, rtol=rtol)
+    assert np.isclose(result.utc, 6.616096470753937, rtol=rtol)
+    assert np.isclose(result.welld, 1159.17968, rtol=rtol)
+    assert np.isclose(result.power, 8.636903762817383, rtol=rtol)
+    assert np.isclose(result.cop, 13.627557754516602, rtol=rtol)
+    assert np.isclose(result.var_opex, -7.510325908660889, rtol=rtol)
+    assert np.isclose(result.fixed_opex, -11.227973937988281, rtol=rtol)
 
 
diff --git a/tests/utc/test_doubletcalc.py b/tests/utc/test_doubletcalc.py
index 770b7d8e79f4b2c6b9a2f7170e572dd6ab0dfc37..1f0f7a475f7bb76f1695a908cf4b035e55b7b60e 100644
--- a/tests/utc/test_doubletcalc.py
+++ b/tests/utc/test_doubletcalc.py
@@ -1,7 +1,6 @@
 from pythermogis.workflow.utc.utc_properties import UTCConfiguration
 from pythermogis.workflow.utc.doublet import DoubletInput
 from pythermogis.workflow.utc.doubletcalc import doubletcalc
-
 import numpy as np
 
 def test_doubletcalc():
@@ -11,7 +10,7 @@ def test_doubletcalc():
         unknown_input_value=-999.0,
         thickness=100.0,
         transmissivity=17500.0,
-        transmissivity_with_ntg=0.0,
+        transmissivity_with_ntg=17500.0,
         ntg=1.0,
         depth=2000,
         porosity=0.0,
diff --git a/tests/utc/test_properties.py b/tests/utc/test_properties.py
index 074468d4733cafe18b80ed92783a1313a0b40101..ab98eb17c08247970d6bd8f91131fd2721919c52 100644
--- a/tests/utc/test_properties.py
+++ b/tests/utc/test_properties.py
@@ -11,7 +11,7 @@ def test_utc_configuration_instantiation():
     assert cfg.n_threads == 1
     assert cfg.is_ates is False
     assert cfg.scenario == "basecase"
-    assert cfg.viscosity_mode == "BATZLEWANG"
+    assert cfg.viscosity_mode == "batzlewang"
     assert cfg.surface_temperature == 10.0
     assert cfg.temp_gradient == 31.0
     assert 30.0 in cfg.p_values
diff --git a/tests/utc/test_well_distance_optimizer.py b/tests/utc/test_well_distance_optimizer.py
index ee8eda2cc687ed52df40b322be1ee8b7a6e3a8ed..00401228f3d88d7f9051a2bcfe6dba0af8f2f36c 100644
--- a/tests/utc/test_well_distance_optimizer.py
+++ b/tests/utc/test_well_distance_optimizer.py
@@ -1,8 +1,12 @@
-from pythermogis.workflow.utc.utc_properties import UTCConfiguration
+import numpy as np
+
 from pythermogis.workflow.utc.doublet import DoubletInput
-from pythermogis.workflow.utc.well_distance import optimize_well_distance
+from pythermogis.workflow.utc.utc_properties import UTCConfiguration
+from pythermogis.workflow.utc.well_distance import (
+    optimize_well_distance,
+    optimize_well_distance_original,
+)
 
-import numpy as np
 
 def test_well_distance_optimizer():
     props = UTCConfiguration()
@@ -18,6 +22,25 @@ def test_well_distance_optimizer():
         temperature=50.0,
     )
 
+    well_distance = optimize_well_distance_original(props, input_data, 6000000, 40)
+
+    assert np.isclose(well_distance, 1011.9140625, rtol=0.015)
+
+
+def test_well_distance_optimizer2():
+    props = UTCConfiguration()
+
+    input_data = DoubletInput(
+        unknown_input_value=-999.0,
+        thickness=100.0,
+        transmissivity=17500.0,
+        transmissivity_with_ntg=0.0,
+        ntg=1.0,
+        depth=2000,
+        porosity=0.0,
+        temperature=50.0,
+    )
+
     well_distance = optimize_well_distance(props, input_data, 6000000, 40)
 
-    assert np.isclose(well_distance, 1011.9140625, rtol=0.000001)
\ No newline at end of file
+    assert np.isclose(well_distance, 1011.9140625, rtol=0.015)
\ No newline at end of file