diff --git a/qualtran/surface_code/azure_cost_model.ipynb b/qualtran/surface_code/azure_cost_model.ipynb
new file mode 100644
index 000000000..7b9423d12
--- /dev/null
+++ b/qualtran/surface_code/azure_cost_model.ipynb
@@ -0,0 +1,711 @@
+{
+ "cells": [
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Azure Cost Model\n",
+    "In this notebook, we go through the cost model in \"Assessing requirements to scale to practical quantum advantage\" by [Beverland et al](https://arxiv.org/abs/2211.07629) to reproduce the estimated costs in Appendix F.\n",
+    "\n",
+    "The paper describes the formulas used for estimating cost in the various appendices. The final estimation procedure is put together in Appendix E. Throughout, we reproduce the costs from the paper; in particular Appendix F."
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The Azure cost model has 5 inputs:\n",
+    "1. The physical assumptions about the hardware (e.g. error rate, the latency of Clifford and measurement operations, etc).\n",
+    "1. A summary of the circuit/algorithm to execute (e.g. number of T gates, etc).\n",
+    "1. The magic state factory (e.g. number of qubits consumed by it, its error rate, etc).\n",
+    "1. Cost model of Approximating rotations using T operations up to error $\\epsilon$.\n",
+    "1. The quantum error correction scheme.\n",
+    "\n",
+    "We take a look at each of these and then reproduce the results of Appendix F."
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## The Inputs\n",
+    "### Physical Parameters\n",
+    "These are assumptions about the quantum hardware and are:\n",
+    "- `t_gate_ns`: Latency time of Clifford gates. \n",
+    "- `t_meas_ns`: Latency time of Measurement operations.\n",
+    "- `physical_error_rate`: Physical error rate ($p$).\n",
+    "\n",
+    "In Qualtran these are represented by the `PhysicalParameters` dataclass."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from qualtran.surface_code.physical_parameters import BEVERLAND_PARAMS\n",
+    "\n",
+    "BEVERLAND_PARAMS"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Algorithm Summary\n",
+    "This is a summary of the circuit or algorithm to execute. This summary is several simple counts:\n",
+    "- `algorithm_qubits` is the number of algorithm qubits.\n",
+    "- `measurements` is the number of Clifford measurements ($M_R$).\n",
+    "- `t_gates` is the number of T gates ($M_T$).\n",
+    "- `toffoli_gates` is the number of Toffoli gates ($M_{Tof}$).\n",
+    "- `rotation_gates` is the number of rotations ($M_R$).\n",
+    "- `rotation_circuit_depth` is the depth of rotation circuit  ($D_R$).\n",
+    "\n",
+    "Note: The symbol in parentheses corresponds to the notation in the paper\n",
+    "\n",
+    "In Qualtran the algorithm specs are represented by the data class `qualtran.surface_code.AlgorithmSummary`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from qualtran.surface_code.algorithm_summary import AlgorithmSummary\n",
+    "\n",
+    "# A circuit that applies a T then measures the qubit.\n",
+    "circuit_summary = AlgorithmSummary(algorithm_qubits=1, t_gates=1, measurements=1)\n",
+    "circuit_summary"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Magic State Factory\n",
+    "\n",
+    "The magic state factory in our case is a T-state factory. The paper describes 15-to-1 factories in Appendix C, but only the overall summary in Table VII in terms of physical qubit count and generation time is used in estimation in Appendix F."
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Rotation Approximation Model\n",
+    "\n",
+    "This is a model that approximates the number of T gates needed to implement a rotation up to error $\\epsilon$. In the paper they use the approximation model from [Kliuchnikov et al](https://arxiv.org/abs/2203.10064). The formula for the number of T gates used to approximate an angle up to error $\\epsilon$ is given by:\n",
+    "$$\n",
+    "a \\log_2{\\frac{1}{\\epsilon}} + b\n",
+    "$$\n",
+    "\n",
+    "Where $a$ and $b$ are constants that depend on the gate set and the approximation protocol. Table 1 of [Kliuchnikov et al](https://arxiv.org/abs/2203.10064) gives estimates for those constants for different combinations of gate sets and protocols.\n",
+    "\n",
+    "In the paper, they use $a=0.53, b=5.3$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from qualtran.surface_code.rotation_cost_model import BeverlandEtAlRotationCost\n",
+    "\n",
+    "BeverlandEtAlRotationCost"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Quantum Error Correction Scheme\n",
+    "The quantum error correction scheme determines three things:\n",
+    "1. The logical error rate given a code distance and physical error rate: $P(d)$.\n",
+    "1. The number of physical qubits needed per logical qubit: $n(d)$.\n",
+    "1. The number of logical time steps: $\\tau(d)$.\n",
+    "\n",
+    "Table V of the paper lists how these are related to the QEC scheme.\n",
+    "\n",
+    "In the paper, they use gate-based QEC which has $P(d), n(d), \\tau(d)$ as:\n",
+    "$$\n",
+    "P(d) = 0.03 \\left ( \\frac{p}{0.01} \\right) ^ \\frac{d+1}{2}\\\\\n",
+    "n(d) = 2 d^2\\\\\n",
+    "\\tau(d) = \\textrm{\\{single stabilizer time\\}} \\cdot d\\\\\n",
+    "$$\n",
+    "\n",
+    "The error detection circuit time depends on several factors physical factors including the time to apply a Clifford, measurement and reset operations as well as classical processing.\n",
+    "\n",
+    "In Table V they don't take into account the classical processing part and assume that a reset takes the same time as a measurement leading to the formula:\n",
+    "$$\n",
+    "\\textrm{\\{single stabilizer time\\}} = 4t_\\textrm{gate} + 2t_\\textrm{meas}\n",
+    "$$\n",
+    "\n",
+    "Other authors (e.g. [Fowler, Gidney](https://arxiv.org/abs/1808.06709)) assume that the entire process takes a specific time (e.g. $\\approx 1\\mu s$).\n",
+    "\n",
+    "In Qualtran, we use \"single stabilizer time\" rather than fall to the low-level hardware parameters. We also default to the Fowler, Gidney parameters (i.e. `FowlerSuperconductingQubits`) when none are given."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from qualtran.surface_code.quantum_error_correction_scheme_summary import (\n",
+    "    BeverlandSuperconductingQubits,\n",
+    ")\n",
+    "\n",
+    "qec = BeverlandSuperconductingQubits\n",
+    "qec"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Resource Estimation\n",
+    "Now we move to reproduce the results in Appendix F."
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Quantum Dynamics\n",
+    "The algorithm specs of this circuit are given as:\n",
+    "- number of algorithm qubits: $100$\n",
+    "- number of rotation gates: $30,100$\n",
+    "- number of measurements: $1.4 \\times 10^6$\n",
+    "- number of T gates: $0$\n",
+    "- number of Toffoli gates: $0$\n",
+    "- depth of rotation circuit: $501$\n",
+    "\n",
+    "with an error budget $\\epsilon$ of $0.001$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "quantum_dynamics_specs = AlgorithmSummary(\n",
+    "    algorithm_qubits=100,\n",
+    "    rotation_gates=30100,\n",
+    "    # Note in the paper the number of measurements\n",
+    "    # has an extra zero which we assume to be a typo.\n",
+    "    measurements=1.4e5,\n",
+    "    rotation_circuit_depth=501,\n",
+    ")\n",
+    "quantum_dynamics_specs"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from qualtran.surface_code import azure_cost_model\n",
+    "\n",
+    "# First calculate the number of logical qubits (Eq. D1)\n",
+    "logical_qubits = azure_cost_model.logical_qubits(quantum_dynamics_specs)\n",
+    "print('Q =', logical_qubits)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Calculate the minimum number of logical time steps (Eq D3)\n",
+    "error_budget = 0.001\n",
+    "c_min = azure_cost_model.minimum_time_steps(\n",
+    "    error_budget=error_budget, alg=quantum_dynamics_specs, rotation_model=BeverlandEtAlRotationCost\n",
+    ")\n",
+    "print('C_min = %e' % c_min)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# And the number of needed T operations (Eq D4)\n",
+    "t_operations = azure_cost_model.t_states(\n",
+    "    error_budget=error_budget, alg=quantum_dynamics_specs, rotation_model=BeverlandEtAlRotationCost\n",
+    ")\n",
+    "print('M = %e' % t_operations)"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Comparing our esimates of $Q = 230, C_{min} = 1.8 \\times 10^5, M = 6 \\times 10^5$ to the paper estimates of \n",
+    "$Q = 230, C_{min} = 1.5 \\times 10^5, M = 2.4 \\times 10^6$. We find a match for $Q$ and $C_{min}$ however we are off by 4x for the number of T gates $M$.\n",
+    "\n",
+    "D4 gives the formula for $M$ as $$M = M_T + 4 M_{Tof} + M_R \\lceil A \\log_2{M_R/\\epsilon_{syn}} + B\\rceil$$\n",
+    "\n",
+    "Since $M_T = M_{Tof} = 0$, the only contribution to $M$ comes from rotations where $M_R = 30100, \\epsilon_{syn} = \\frac{10^{-3}}{3} , A = 0.53, B = 5.3$ which give our value of $6.02e5$."
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Next, we estimate the code distance $d$. $d$ satisfies \n",
+    "$$\n",
+    "P(d) = a \\left ( \\frac{p}{p^*} \\right )^\\frac{d+1}{2} \n",
+    "$$\n",
+    "subject ot the constraint on the logical error rate of $Q \\cdot C \\cdot P(d) = \\frac{\\epsilon}{3}$. Where $p$ is the physical error rate. $a$ and $p^*$ are constants determined by the QEC scheme."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "d = azure_cost_model.code_distance(\n",
+    "    error_budget=error_budget,\n",
+    "    time_steps=c_min,\n",
+    "    alg=quantum_dynamics_specs,\n",
+    "    qec=BeverlandSuperconductingQubits,\n",
+    "    physical_error_rate=BEVERLAND_PARAMS.physical_error_rate,\n",
+    ")\n",
+    "print(f'{d=}')"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Matching the paper's code distance of $d = 9$. This leads to a total run time (Eq. E3) of 0.65s which is close to the time in the paper of 0.55s"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "t_s = qec.error_detection_circuit_time_us(d) * 1e-6 * c_min\n",
+    "f'algorithm run time of {t_s:g} seconds'.format()"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Next, we examine the magic state factories. In the paper, for the quantum dynamics example, we are given $199$ factories each producing one T state every $46.8 \\mu s$ at an error rate of $5.6e-11$ while consuming $3,240$ qubits. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "num_factories = 199\n",
+    "factory_qubits = 3240"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Leading to a total number of physical qubits (from E6)\n",
+    "distillation_qubits = num_factories * factory_qubits\n",
+    "q = distillation_qubits + logical_qubits * qec.physical_qubits(d)\n",
+    "print('total number of physical qubits:', q)\n",
+    "print('percentage of distillation qubits: {}%'.format(round(distillation_qubits / q * 100, 1)))"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Our estimate of 0.68M physical qubits with 94.5% of them being consumed by the T states factories match the numbers in the paper."
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Quantum Chemistry\n",
+    "The algorithm specs of this circuit are given as:\n",
+    "- number of algorithm qubits: 1318\n",
+    "- number of rotation gates: $2.06 \\times 10^8$\n",
+    "- number of measurements: $1.37 \\times 10^9$\n",
+    "- number of T gates: $5.53 \\times 10^7$\n",
+    "- number of Toffoli gates: $1.35 \\times 10^{11}$\n",
+    "- depth of rotation circuit: $2.05 \\times 10^8$\n",
+    "\n",
+    "with an error budget of 0.01"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "quantum_chemistry_specs = AlgorithmSummary(\n",
+    "    algorithm_qubits=1318,\n",
+    "    rotation_gates=2.06e8,\n",
+    "    measurements=1.37e9,\n",
+    "    toffoli_gates=1.35e11,\n",
+    "    t_gates=5.53e7,\n",
+    "    rotation_circuit_depth=2.05e8,\n",
+    ")\n",
+    "quantum_chemistry_specs"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "logical_qubits = azure_cost_model.logical_qubits(quantum_chemistry_specs)\n",
+    "print('Q =', logical_qubits)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Calculate the minimum number of logical time steps (Eq D3)\n",
+    "error_budget = 0.01\n",
+    "c_min = azure_cost_model.minimum_time_steps(\n",
+    "    error_budget=error_budget, alg=quantum_chemistry_specs, rotation_model=BeverlandEtAlRotationCost\n",
+    ")\n",
+    "print('C_{min} = %g' % c_min)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# And the number of needed T operations (Eq D4)\n",
+    "t_operations = azure_cost_model.t_states(\n",
+    "    error_budget=error_budget, alg=quantum_chemistry_specs, rotation_model=BeverlandEtAlRotationCost\n",
+    ")\n",
+    "print('M = %g' % t_operations)"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Note that our estimates match nicely to the paper estimates of \n",
+    "$$Q = 2740\\\\\n",
+    "C_{min} = 4.10 \\times 10^{11}\\\\\n",
+    "M = 5.44 \\times 10^{11}$$"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Next, we estimate the code distance."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "d = azure_cost_model.code_distance(\n",
+    "    error_budget=error_budget,\n",
+    "    time_steps=c_min,\n",
+    "    alg=quantum_chemistry_specs,\n",
+    "    qec=qec,\n",
+    "    physical_error_rate=BEVERLAND_PARAMS.physical_error_rate,\n",
+    ")\n",
+    "print(f'{d=}')"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Again we match the code distance from the paper of $d = 17$. This leads to a total run time (Eq. E3) of"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "total_seconds = qec.error_detection_circuit_time_us(d) * 1e-6 * c_min\n",
+    "total_days = total_seconds / 3600 / 24\n",
+    "'algorithm run time of %g days' % (total_days)"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We get a run time estimate of 32.4 days. In the paper, it says the run time is 1 month and 1 day."
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Next we examine the magic state factories. In the paper, for the quantum chemistry example, we are given $17$ factories each producing one T state every $83.2\\mu s$ at an error rate of $2.13e-15$ while consuming $16,000$ qubits. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "num_factories = 17\n",
+    "factory_qubits = 16000"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Leading to a total number of physical qubits (from E6)\n",
+    "distillation_qubits = num_factories * factory_qubits\n",
+    "q = distillation_qubits + logical_qubits * qec.physical_qubits(d)\n",
+    "print('total number of physical qubits: %g M' % round(q * 1e-6, 2))\n",
+    "print('percentage of distillation qubits: {}%'.format(round(distillation_qubits / q * 100, 1)))"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Our estimate of 1.86M physical qubits matches the paper's estimate."
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Shor Factoring\n",
+    "The algorithm specs of this circuit are given as:\n",
+    "- number of algorithm qubits: $12,581$\n",
+    "- number of rotation gates: $12$\n",
+    "- number of measurements: $1.08 \\times 10^9$\n",
+    "- number of T gates: 12\n",
+    "- number of Toffoli gates: $3.73 \\times 10^{10}$\n",
+    "- depth of rotation circuit: $12$\n",
+    "\n",
+    "with an error budget of $1/3$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "shor_specs = AlgorithmSummary(\n",
+    "    algorithm_qubits=12581,\n",
+    "    rotation_gates=12,\n",
+    "    measurements=1.08e9,\n",
+    "    rotation_circuit_depth=12,\n",
+    "    # Note in the paper the number of Toffoli operations is 3.73e10.\n",
+    "    # However we assume that the exponent has a typo and that the number is 3.73e9.\n",
+    "    toffoli_gates=3.73e9,\n",
+    "    t_gates=12,\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "logical_qubits = azure_cost_model.logical_qubits(shor_specs)\n",
+    "print('Q =', logical_qubits)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Calculate the minimum number of logical time steps (Eq D3)\n",
+    "error_budget = 1 / 3\n",
+    "c_min = azure_cost_model.minimum_time_steps(\n",
+    "    error_budget=error_budget, alg=shor_specs, rotation_model=BeverlandEtAlRotationCost\n",
+    ")\n",
+    "print('C_min = %e' % c_min)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# And the number of needed T operations (Eq D4)\n",
+    "t_operations = azure_cost_model.t_states(\n",
+    "    error_budget=error_budget, alg=shor_specs, rotation_model=BeverlandEtAlRotationCost\n",
+    ")\n",
+    "print('M = %e' % t_operations)"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Our estimates of $Q = 25481, C_{min} = 1.23 \\cdot 10^{10}, M = 1.49 \\cdot 10^{10}$ match the estimates of the paper."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "d = azure_cost_model.code_distance(\n",
+    "    error_budget=error_budget,\n",
+    "    time_steps=c_min,\n",
+    "    alg=shor_specs,\n",
+    "    qec=qec,\n",
+    "    physical_error_rate=BEVERLAND_PARAMS.physical_error_rate,\n",
+    ")\n",
+    "print(f'{d=}')"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Matching the code distance of $d = 13$ in the paper."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "total_seconds = qec.error_detection_circuit_time_us(d) * 1e-6 * c_min\n",
+    "total_hours = total_seconds / 3600\n",
+    "h = int(total_hours)\n",
+    "m = (total_hours - h) * 60\n",
+    "'algorithm run time of %d hours %d' % (h, m)"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Our estimate runtime of 17h43m matches with the estimate of the paper."
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Next we examine the magic state factories. In the paper, for the quantum chemistry example, we are given $18$ factories each producing one T state every $72.8 \\mu s$ at an error rate of $5.51e-13$ while consuming $5,760$ qubits. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "num_factories = 18\n",
+    "factory_qubits = 5760"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Leading to a total number of physical qubits (from E6)\n",
+    "distillation_qubits = num_factories * factory_qubits\n",
+    "q = distillation_qubits + logical_qubits * qec.physical_qubits(d)\n",
+    "print('total number of physical qubits: %g M' % round(q * 1e-6, 2))\n",
+    "print('percentage of distillation qubits: {}%'.format(round(distillation_qubits / q * 100, 1)))"
+   ]
+  },
+  {
+   "attachments": {},
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Our estimate of 8.72M physical qubits matches the paper's estimate."
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "qualtran",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.10.13"
+  },
+  "orig_nbformat": 4
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/qualtran/surface_code/azure_cost_model.py b/qualtran/surface_code/azure_cost_model.py
new file mode 100644
index 000000000..afa18b99d
--- /dev/null
+++ b/qualtran/surface_code/azure_cost_model.py
@@ -0,0 +1,127 @@
+#  Copyright 2023 Google LLC
+#
+#  Licensed under the Apache License, Version 2.0 (the "License");
+#  you may not use this file except in compliance with the License.
+#  You may obtain a copy of the License at
+#
+#      https://www.apache.org/licenses/LICENSE-2.0
+#
+#  Unless required by applicable law or agreed to in writing, software
+#  distributed under the License is distributed on an "AS IS" BASIS,
+#  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#  See the License for the specific language governing permissions and
+#  limitations under the License.
+
+import math
+
+from qualtran.surface_code.algorithm_summary import AlgorithmSummary
+from qualtran.surface_code.quantum_error_correction_scheme_summary import (
+    QuantumErrorCorrectionSchemeSummary,
+)
+from qualtran.surface_code.rotation_cost_model import RotationCostModel
+
+
+def logical_qubits(algorithm_specs: AlgorithmSummary) -> int:
+    r"""Number of logical qubits needed for the algorithm.
+
+    Equals:
+    $$
+        2 Q_\mathrm{alg} + \lceil \sqrt{8 Q_\mathrm{alg}} \rceil + 1
+    $$
+
+    Source: Equation D1 in https://arxiv.org/abs/2211.07629.
+
+    Args:
+        algorithm_specs: A summary of an algorithm/circuit.
+    """
+    q_alg = algorithm_specs.algorithm_qubits
+    return math.ceil(2 * q_alg + math.sqrt(8 * q_alg) + 1)
+
+
+def minimum_time_steps(
+    error_budget: float, alg: AlgorithmSummary, rotation_model: RotationCostModel
+) -> int:
+    r"""Minimum number of time steps needed for the algorithm.
+
+    Equals
+    $$
+        M_\mathrm{meas} + M_R + M_T + 3 M_\mathrm{Tof} + D_R \textrm{rotation cost}
+    $$
+    Where:
+        $M_\mathrm{meas}$ is the number of measurements.
+        $M_R$ is the number of rotations.
+        $M_T$ is the number of T operations.
+        $M_mathrm{Tof}$ is the number of toffoli operations.
+        $D_R$ is the depth of the rotation circuit.
+        $\textrm{rotation cost}$ is the number of T operations needed to approximate a rotation to $\epsilon/(3*M_R)$.
+    Source: Equation D3 in https://arxiv.org/abs/2211.07629.
+
+    Args:
+        error_budget: Error Budget.
+        alg: A summary of an algorithm/circuit.
+        rotation_model: Cost model used to compute the number of T gates
+            needed to approximate rotations.
+    """
+    c_min = math.ceil(alg.measurements + alg.rotation_gates + alg.t_gates + 3 * alg.toffoli_gates)
+    eps_syn = error_budget / 3
+    c_min += math.ceil(
+        alg.rotation_circuit_depth
+        * rotation_model.rotation_cost(eps_syn / alg.rotation_gates).t_gates
+    )
+    return c_min
+
+
+def code_distance(
+    error_budget: float,
+    time_steps: float,
+    alg: AlgorithmSummary,
+    qec: QuantumErrorCorrectionSchemeSummary,
+    physical_error_rate: float,
+) -> int:
+    r"""Minimum code distance needed to run the algorithm within the error budget.
+
+    This is the code distance $d$ that satisfies $QCP = \epsilon/3$. Where:
+        $\epsilon$ is the error budget.
+        Q is the number of logical qubits.
+        C is the number of time steps.
+        P(d) is the logical error rate.
+
+    Args:
+        error_budget: Error Budget.
+        time_steps: Number of time steps used to run the algorithm.
+        alg: A summary of an algorithm/circuit.
+        qec: Quantum Error Correction Scheme.
+        physical_error_rate: The physical error rate of the device.
+    """
+    q = logical_qubits(alg)
+    return qec.code_distance_from_budget(physical_error_rate, error_budget / (3 * q * time_steps))
+
+
+def t_states(
+    error_budget: float, alg: AlgorithmSummary, rotation_model: RotationCostModel
+) -> float:
+    r"""Total number of T states consumed by the algorithm.
+
+    Equals
+    $$
+       M_T + 4 M_\mathrm{Tof} + M_R \textrm{rotation cost}
+    $$
+    Where:
+        $M_R$ is the number of rotations.
+        $M_T$ is the number of T operations.
+        $M_mathrm{Tof}$ is the number of toffoli operations.
+        $\textrm{rotation cost}$ is the number of T operations needed to approximate a rotation to $\epsilon/(3*M_R)$.
+    Source: D4 in https://arxiv.org/abs/2211.07629.
+
+    Args:
+        error_budget: Error Budget.
+        alg: A summary of an algorithm/circuit.
+        rotation_model: Cost model used to compute the number of T gates
+            needed to approximate rotations.
+    """
+    eps_syn = error_budget / 3
+    return (
+        alg.t_gates
+        + 4 * alg.toffoli_gates
+        + alg.rotation_gates * rotation_model.rotation_cost(eps_syn / alg.rotation_gates).t_gates
+    )
diff --git a/qualtran/surface_code/azure_cost_model_test.py b/qualtran/surface_code/azure_cost_model_test.py
new file mode 100644
index 000000000..3800cdc69
--- /dev/null
+++ b/qualtran/surface_code/azure_cost_model_test.py
@@ -0,0 +1,120 @@
+#  Copyright 2023 Google LLC
+#
+#  Licensed under the Apache License, Version 2.0 (the "License");
+#  you may not use this file except in compliance with the License.
+#  You may obtain a copy of the License at
+#
+#      https://www.apache.org/licenses/LICENSE-2.0
+#
+#  Unless required by applicable law or agreed to in writing, software
+#  distributed under the License is distributed on an "AS IS" BASIS,
+#  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#  See the License for the specific language governing permissions and
+#  limitations under the License.
+
+import pytest
+from attrs import frozen
+
+from qualtran.surface_code import azure_cost_model
+from qualtran.surface_code.algorithm_summary import AlgorithmSummary
+from qualtran.surface_code.quantum_error_correction_scheme_summary import (
+    BeverlandSuperconductingQubits,
+)
+from qualtran.surface_code.rotation_cost_model import BeverlandEtAlRotationCost
+
+
+@frozen
+class Test:
+    alg: AlgorithmSummary
+    q_alg: float
+
+    error_budget: float
+    c_min: float
+
+    time_steps: float
+    code_distance: float
+
+    t_states: float
+
+
+_TESTS = [
+    Test(
+        alg=AlgorithmSummary(
+            algorithm_qubits=100,
+            rotation_gates=30100,
+            measurements=1.4 * 10**6,
+            rotation_circuit_depth=501,
+        ),
+        q_alg=230,
+        error_budget=1e-3,
+        c_min=1.5e6,
+        time_steps=1.5e5,
+        code_distance=9,
+        t_states=602000,
+    ),
+    Test(
+        alg=AlgorithmSummary(
+            algorithm_qubits=1318,
+            t_gates=5.53e7,
+            rotation_circuit_depth=2.05e8,
+            rotation_gates=2.06e8,
+            toffoli_gates=1.35e11,
+            measurements=1.37e9,
+        ),
+        q_alg=2740,
+        error_budget=1e-2,
+        c_min=4.1e11,
+        time_steps=4.1e11,
+        code_distance=17,
+        t_states=5.44e11,
+    ),
+    Test(
+        alg=AlgorithmSummary(
+            algorithm_qubits=12581,
+            t_gates=12,
+            rotation_circuit_depth=12,
+            rotation_gates=12,
+            toffoli_gates=3.73e9,
+            measurements=1.08e9,
+        ),
+        q_alg=25481,
+        error_budget=1 / 3,
+        c_min=1.23e10,
+        time_steps=1.23e10,
+        code_distance=13,
+        t_states=1.49e10,
+    ),
+]
+
+
+@pytest.mark.parametrize('test', _TESTS)
+def test_logical_qubits(test: Test):
+    assert azure_cost_model.logical_qubits(test.alg) == test.q_alg
+
+
+@pytest.mark.parametrize('test', _TESTS)
+def test_minimum_time_step(test: Test):
+    got = azure_cost_model.minimum_time_steps(
+        test.error_budget, test.alg, rotation_model=BeverlandEtAlRotationCost
+    )
+    assert got == pytest.approx(test.c_min, rel=0.1)
+
+
+@pytest.mark.parametrize('test', _TESTS)
+def test_code_distance(test: Test):
+    got = azure_cost_model.code_distance(
+        test.error_budget,
+        test.time_steps,
+        test.alg,
+        qec=BeverlandSuperconductingQubits,
+        physical_error_rate=1e-4,
+    )
+    assert got == test.code_distance
+
+
+@pytest.mark.parametrize('test', _TESTS)
+def test_t_states(test: Test):
+    got = azure_cost_model.t_states(
+        test.error_budget, test.alg, rotation_model=BeverlandEtAlRotationCost
+    )
+    assert got == pytest.approx(test.t_states, rel=0.1)
diff --git a/qualtran/surface_code/physical_parameters.py b/qualtran/surface_code/physical_parameters.py
index 1d4f8a745..9944de8af 100644
--- a/qualtran/surface_code/physical_parameters.py
+++ b/qualtran/surface_code/physical_parameters.py
@@ -35,7 +35,7 @@ class PhysicalParameters:
     reference: Optional[str] = None
 
 
-BeverlandEtAl = PhysicalParameters(
+BEVERLAND_PARAMS = PhysicalParameters(
     t_gate_ns=50,  # 50ns
     t_meas_ns=100,  # 100ns
     physical_error_rate=1e-4,