diff --git a/Cargo.lock b/Cargo.lock
index 969ffca..c585470 100644
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -121,9 +121,9 @@ checksum = "e78d4f1cc4ae33bbfc157ed5d5a5ef3bc29227303d595861deb238fcec4e9457"
 
 [[package]]
 name = "env_logger"
-version = "0.8.4"
+version = "0.9.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "a19187fea3ac7e84da7dacf48de0c45d63c6a76f9490dae389aead16c243fce3"
+checksum = "0b2cf0344971ee6c64c31be0d530793fba457d322dfec2810c453d0ef228f9c3"
 dependencies = [
  "atty",
  "humantime",
@@ -564,7 +564,7 @@ dependencies = [
 
 [[package]]
 name = "simquil"
-version = "0.1.1"
+version = "0.2.0"
 dependencies = [
  "anyhow",
  "env_logger",
diff --git a/Cargo.toml b/Cargo.toml
index ff0664b..b43bfd7 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -22,4 +22,4 @@ thiserror = "1.0.24"
 anyhow = "1.0.39"
 
 [dev-dependencies]
-pretty_assertions = "0.7.2"
+pretty_assertions = "0.7.2"
\ No newline at end of file
diff --git a/src/gates/mod.rs b/src/gates/mod.rs
index 461ec4b..d32c122 100644
--- a/src/gates/mod.rs
+++ b/src/gates/mod.rs
@@ -3,7 +3,7 @@ pub mod standard;
 use ndarray::{Array, Array2};
 use num::complex::Complex64;
 
-use crate::gates::standard::{ccnot, cnot, cz, h, i, rx, ry, rz, swap, x, z};
+use crate::gates::standard::swap;
 
 #[derive(Default, Clone, Debug)]
 pub struct QGate {
@@ -11,34 +11,12 @@ pub struct QGate {
     qubits: Vec<usize>,
 }
 
-pub fn gate_matrix(name: String, params: Vec<f64>, qubits: Vec<usize>) -> QGate {
-    match name.as_str() {
-        "I" => i(qubits[0]),
-        "X" => x(qubits[0]),
-        "Y" => x(qubits[0]),
-        "Z" => z(qubits[0]),
-        "RX" => rx(params[0], qubits[0]),
-        "RY" => ry(params[0], qubits[0]),
-        "RZ" => rz(params[0], qubits[0]),
-        "H" => h(qubits[0]),
-        "CZ" => cz(qubits[0], qubits[1]),
-        "CNOT" => cnot(qubits[0], qubits[1]),
-        "SWAP" => swap(qubits[0], qubits[1]),
-        "CCNOT" => ccnot(qubits[0], qubits[1], qubits[2]),
-        _ => todo!(),
-    }
-}
-
-// TODO Is it possible to define a SquareArray2 type (where dimensions are equal)
-// TODO Likewise, an array type which is both square and each dimension is 2^n
-
 impl QGate {
     fn lift_adjacent(&self, i: usize, n_qubits: usize) -> Array2<Complex64> {
         let gate_size = self.qubits.len();
         let n = num::pow(2, i);
         let bottom = Array::eye(n);
-        let top =
-            Array::eye(2u64.pow((n_qubits as u32) - (i as u32) - (gate_size as u32)) as usize);
+        let top = Array::eye(2u64.pow((n_qubits - i - gate_size) as u32) as usize);
 
         kron(&top, &kron(&self.matrix, &bottom))
     }
@@ -116,11 +94,12 @@ mod tests {
     use ndarray::arr2;
     use pretty_assertions::assert_eq;
 
+    use crate::gates::standard::{cnot, x};
     use crate::matrix::{C0, C1};
 
     #[test]
     fn lift_1q() {
-        let gate = super::x(0);
+        let gate = x(0);
         assert_eq!(
             gate.lift_adjacent(0, 2),
             arr2(&[
@@ -131,7 +110,7 @@ mod tests {
             ])
         );
 
-        let gate = super::x(1);
+        let gate = x(1);
         assert_eq!(
             gate.lift_adjacent(1, 2),
             arr2(&[
@@ -150,7 +129,7 @@ mod tests {
 
     #[test]
     fn lift_2q() {
-        let gate = super::cnot(0, 1);
+        let gate = cnot(0, 1);
         assert_eq!(
             gate.lift(2),
             arr2(&[
@@ -161,7 +140,7 @@ mod tests {
             ])
         );
 
-        let gate = super::cnot(1, 0);
+        let gate = cnot(1, 0);
         assert_eq!(
             gate.lift(2),
             arr2(&[
@@ -172,7 +151,7 @@ mod tests {
             ])
         );
 
-        let gate = super::cnot(0, 2);
+        let gate = cnot(0, 2);
         assert_eq!(
             gate.lift(3),
             arr2(&[
@@ -187,7 +166,7 @@ mod tests {
             ])
         );
 
-        let gate = super::cnot(2, 0);
+        let gate = cnot(2, 0);
         assert_eq!(
             gate.lift(3),
             arr2(&[
diff --git a/src/gates/standard.rs b/src/gates/standard.rs
index b289783..959478f 100644
--- a/src/gates/standard.rs
+++ b/src/gates/standard.rs
@@ -1,125 +1,250 @@
-use ndarray::arr2;
-
 use crate::matrix::{C0, C1, I1};
+use ndarray::arr2;
+use std::f64::consts::FRAC_PI_4;
+use thiserror::Error;
 
 use super::QGate;
 
-pub fn i(q: usize) -> QGate {
-    QGate {
-        matrix: arr2(&[[C1, C0], [C0, C1]]),
-        qubits: [q].to_vec(),
-    }
+#[derive(Error, Debug)]
+pub enum GateError {
+    #[error("Unknown gate {0}")]
+    UnknownGate(String),
+    #[error("Gate matrix missing qubit")]
+    GateMatrixMissingQubit,
+    #[error("Gate matrix missing parameter")]
+    GateMatrixMissingParameter,
 }
 
-pub fn x(q: usize) -> QGate {
-    QGate {
-        matrix: arr2(&[[C0, C1], [C1, C0]]),
-        qubits: [q].to_vec(),
-    }
+macro_rules! define_gate {
+    ($name:ident, [$($qubit:ident),+], $gate:expr) => {
+        pub fn $name($($qubit: usize),*) -> QGate {
+            QGate {
+                matrix: arr2(&$gate),
+                qubits: [$($qubit),*].to_vec(),
+            }
+        }
+    };
+    ($name:ident, [$($param:ident),+], [$($qubit:ident),+], $gate:expr) => {
+        pub fn $name($($param: f64),*,$($qubit: usize),*) -> QGate {
+            QGate {
+                matrix: arr2(&$gate),
+                qubits: [$($qubit),*].to_vec(),
+            }
+        }
+    };
 }
 
-pub fn y(q: usize) -> QGate {
-    QGate {
-        matrix: arr2(&[[C0, -I1], [I1, C0]]),
-        qubits: [q].to_vec(),
-    }
-}
+define_gate!(i, [q], [[C1, C0], [C0, C1]]);
+define_gate!(x, [q], [[C0, C1], [C1, C0]]);
+define_gate!(y, [q], [[C0, -I1], [I1, C0]]);
+define_gate!(z, [q], [[C1, C0], [C0, -C1]]);
+define_gate!(
+    h,
+    [q],
+    [
+        [1.0 / 2.0f64.sqrt() + C0, 1.0 / 2.0f64.sqrt() + C0],
+        [1.0 / 2.0f64.sqrt() + C0, -1.0 / 2.0f64.sqrt() + C0],
+    ]
+);
+define_gate!(s, [q], [[C1, C0], [C0, I1]]);
+define_gate!(t, [q], [[C1, C0], [C0, FRAC_PI_4.cos() + C0]]);
 
-pub fn z(q: usize) -> QGate {
-    QGate {
-        matrix: arr2(&[[C1, C0], [C0, -C1]]),
-        qubits: [q].to_vec(),
-    }
-}
-
-pub fn h(q: usize) -> QGate {
-    QGate {
-        matrix: arr2(&[
-            [1.0 / 2.0f64.sqrt() + C0, 1.0 / 2.0f64.sqrt() + C0],
-            [1.0 / 2.0f64.sqrt() + C0, -1.0 / 2.0f64.sqrt() + C0],
-        ]),
-        qubits: [q].to_vec(),
-    }
-}
-
-pub fn rx(param: f64, q: usize) -> QGate {
-    QGate {
-        matrix: arr2(&[
-            [(param / 2.0).cos() + C0, -I1 * (param / 2.0).sin()],
-            [-I1 * (param / 2.0).sin(), (param / 2.0).cos() + C0],
-        ]),
-        qubits: [q].to_vec(),
-    }
-}
+define_gate!(
+    rx,
+    [theta],
+    [q],
+    [
+        [(theta / 2.0).cos() + C0, -I1 * (theta / 2.0).sin()],
+        [-I1 * (theta / 2.0).sin(), (theta / 2.0).cos() + C0],
+    ]
+);
+define_gate!(
+    ry,
+    [theta],
+    [q],
+    [
+        [(theta / 2.0).cos() + C0, -(theta / 2.0).sin() + C0],
+        [(theta / 2.0).sin() + C0, (theta / 2.0).cos() + C0],
+    ]
+);
+define_gate!(
+    rz,
+    [theta],
+    [q],
+    [
+        [(-theta / 2.0).cos() + I1 * (-theta / 2.0).sin(), C0],
+        [C0, (-theta / 2.0).cos() + I1 * (theta / 2.0).sin()],
+    ]
+);
+define_gate!(
+    phase,
+    [theta],
+    [q],
+    [[C1, C0], [C0, theta.cos() + I1 * theta.sin()]]
+);
+define_gate!(
+    cphase,
+    [theta],
+    [q],
+    [
+        [C1, C0, C0, C0],
+        [C0, C1, C0, C0],
+        [C0, C0, C1, C0],
+        [C0, C0, C0, theta.cos() + I1 * theta.sin()],
+    ]
+);
 
-pub fn ry(param: f64, q: usize) -> QGate {
-    QGate {
-        matrix: arr2(&[
-            [(param / 2.0).cos() + C0, -(param / 2.0).sin() + C0],
-            [(param / 2.0).sin() + C0, (param / 2.0).cos() + C0],
-        ]),
-        qubits: [q].to_vec(),
-    }
-}
-
-pub fn rz(param: f64, q: usize) -> QGate {
-    QGate {
-        matrix: arr2(&[
-            [(-param / 2.0).cos() + I1 * (-param / 2.0).sin(), C0],
-            [C0, (-param / 2.0).cos() + I1 * (param / 2.0).sin()],
-        ]),
-        qubits: [q].to_vec(),
-    }
-}
-
-pub fn cz(control: usize, target: usize) -> QGate {
-    QGate {
-        matrix: arr2(&[
-            [C1, C0, C0, C0],
-            [C0, C1, C0, C0],
-            [C0, C0, C1, C0],
-            [C0, C0, C0, -C1],
-        ]),
-        qubits: [control, target].to_vec(),
-    }
-}
-
-pub fn cnot(control: usize, target: usize) -> QGate {
-    QGate {
-        matrix: arr2(&[
-            [C1, C0, C0, C0],
-            [C0, C1, C0, C0],
-            [C0, C0, C0, C1],
-            [C0, C0, C1, C0],
-        ]),
-        qubits: [control, target].to_vec(),
-    }
-}
-
-pub fn swap(q0: usize, q1: usize) -> QGate {
-    QGate {
-        matrix: arr2(&[
-            [C1, C0, C0, C0],
-            [C0, C0, C1, C0],
-            [C0, C1, C0, C0],
-            [C0, C0, C0, C1],
-        ]),
-        qubits: [q0, q1].to_vec(),
-    }
-}
+define_gate!(
+    cz,
+    [control, target],
+    [
+        [C1, C0, C0, C0],
+        [C0, C1, C0, C0],
+        [C0, C0, C1, C0],
+        [C0, C0, C0, -C1],
+    ]
+);
+define_gate!(
+    cnot,
+    [control, target],
+    [
+        [C1, C0, C0, C0],
+        [C0, C1, C0, C0],
+        [C0, C0, C0, C1],
+        [C0, C0, C1, C0],
+    ]
+);
+define_gate!(
+    swap,
+    [q0, q1],
+    [
+        [C1, C0, C0, C0],
+        [C0, C0, C1, C0],
+        [C0, C1, C0, C0],
+        [C0, C0, C0, C1],
+    ]
+);
+define_gate!(
+    iswap,
+    [q0, q1],
+    [
+        [C1, C0, C0, C0],
+        [C0, C0, I1, C0],
+        [C0, I1, C0, C0],
+        [C0, C0, C0, C1],
+    ]
+);
+define_gate!(
+    cswap,
+    [control, q0, q1],
+    [
+        [C1, C0, C0, C0, C0, C0, C0, C0],
+        [C0, C1, C0, C0, C0, C0, C0, C0],
+        [C0, C0, C1, C0, C0, C0, C0, C0],
+        [C0, C0, C0, C1, C0, C0, C0, C0],
+        [C0, C0, C0, C0, C1, C0, C0, C0],
+        [C0, C0, C0, C0, C0, C0, C1, C0],
+        [C0, C0, C0, C0, C0, C1, C0, C0],
+        [C0, C0, C0, C0, C0, C0, C0, C1],
+    ]
+);
+define_gate!(
+    ccnot,
+    [control0, control1, target],
+    [
+        [C1, C0, C0, C0, C0, C0, C0, C0],
+        [C0, C1, C0, C0, C0, C0, C0, C0],
+        [C0, C0, C1, C0, C0, C0, C0, C0],
+        [C0, C0, C0, C1, C0, C0, C0, C0],
+        [C0, C0, C0, C0, C1, C0, C0, C0],
+        [C0, C0, C0, C0, C0, C1, C0, C0],
+        [C0, C0, C0, C0, C0, C0, C0, C1],
+        [C0, C0, C0, C0, C0, C0, C1, C0],
+    ]
+);
 
-pub fn ccnot(control0: usize, control1: usize, target: usize) -> QGate {
-    QGate {
-        matrix: arr2(&[
-            [C1, C0, C0, C0, C0, C0, C0, C0],
-            [C0, C1, C0, C0, C0, C0, C0, C0],
-            [C0, C0, C1, C0, C0, C0, C0, C0],
-            [C0, C0, C0, C1, C0, C0, C0, C0],
-            [C0, C0, C0, C0, C1, C0, C0, C0],
-            [C0, C0, C0, C0, C0, C1, C0, C0],
-            [C0, C0, C0, C0, C0, C0, C0, C1],
-            [C0, C0, C0, C0, C0, C0, C1, C0],
-        ]),
-        qubits: [control0, control1, target].to_vec(),
+pub fn gate_matrix(name: String, params: Vec<f64>, qubits: Vec<usize>) -> Result<QGate, GateError> {
+    match name.as_str() {
+        "I" => {
+            let q = qubits.get(0).ok_or(GateError::GateMatrixMissingQubit)?;
+            Ok(i(*q))
+        }
+        "X" => {
+            let q = qubits.get(0).ok_or(GateError::GateMatrixMissingQubit)?;
+            Ok(x(*q))
+        }
+        "Y" => {
+            let q = qubits.get(0).ok_or(GateError::GateMatrixMissingQubit)?;
+            Ok(y(*q))
+        }
+        "Z" => {
+            let q = qubits.get(0).ok_or(GateError::GateMatrixMissingQubit)?;
+            Ok(z(*q))
+        }
+        "H" => {
+            let q = qubits.get(0).ok_or(GateError::GateMatrixMissingQubit)?;
+            Ok(h(*q))
+        }
+        "T" => {
+            let q = qubits.get(0).ok_or(GateError::GateMatrixMissingQubit)?;
+            Ok(t(*q))
+        }
+        "S" => {
+            let q = qubits.get(0).ok_or(GateError::GateMatrixMissingQubit)?;
+            Ok(s(*q))
+        }
+        "RX" => {
+            let q = qubits.get(0).ok_or(GateError::GateMatrixMissingQubit)?;
+            let p = params.get(0).ok_or(GateError::GateMatrixMissingParameter)?;
+            Ok(rx(*p, *q))
+        }
+        "RY" => {
+            let q = qubits.get(0).ok_or(GateError::GateMatrixMissingQubit)?;
+            let p = params.get(0).ok_or(GateError::GateMatrixMissingParameter)?;
+            Ok(ry(*p, *q))
+        }
+        "RZ" => {
+            let q = qubits.get(0).ok_or(GateError::GateMatrixMissingQubit)?;
+            let p = params.get(0).ok_or(GateError::GateMatrixMissingParameter)?;
+            Ok(rz(*p, *q))
+        }
+        "PHASE" => {
+            let q = qubits.get(0).ok_or(GateError::GateMatrixMissingQubit)?;
+            let p = params.get(0).ok_or(GateError::GateMatrixMissingParameter)?;
+            Ok(phase(*p, *q))
+        }
+        "CZ" => {
+            let control = qubits.get(0).ok_or(GateError::GateMatrixMissingQubit)?;
+            let target = qubits.get(1).ok_or(GateError::GateMatrixMissingQubit)?;
+            Ok(cz(*control, *target))
+        }
+        "CNOT" => {
+            let control = qubits.get(0).ok_or(GateError::GateMatrixMissingQubit)?;
+            let target = qubits.get(1).ok_or(GateError::GateMatrixMissingQubit)?;
+            Ok(cnot(*control, *target))
+        }
+        "SWAP" => {
+            let q0 = qubits.get(0).ok_or(GateError::GateMatrixMissingQubit)?;
+            let q1 = qubits.get(1).ok_or(GateError::GateMatrixMissingQubit)?;
+            Ok(cz(*q0, *q1))
+        }
+        "ISWAP" => {
+            let q0 = qubits.get(0).ok_or(GateError::GateMatrixMissingQubit)?;
+            let q1 = qubits.get(1).ok_or(GateError::GateMatrixMissingQubit)?;
+            Ok(iswap(*q0, *q1))
+        }
+        "CSWAP" => {
+            let control0 = qubits.get(0).ok_or(GateError::GateMatrixMissingQubit)?;
+            let control1 = qubits.get(1).ok_or(GateError::GateMatrixMissingQubit)?;
+            let target = qubits.get(2).ok_or(GateError::GateMatrixMissingQubit)?;
+            Ok(cswap(*control0, *control1, *target))
+        }
+        "CCNOT" => {
+            let control0 = qubits.get(0).ok_or(GateError::GateMatrixMissingQubit)?;
+            let control1 = qubits.get(1).ok_or(GateError::GateMatrixMissingQubit)?;
+            let target = qubits.get(1).ok_or(GateError::GateMatrixMissingQubit)?;
+            Ok(ccnot(*control0, *control1, *target))
+        }
+        _ => Err(GateError::UnknownGate(name)),
     }
 }
diff --git a/src/main.rs b/src/main.rs
index b831a76..cb1e4e8 100644
--- a/src/main.rs
+++ b/src/main.rs
@@ -24,8 +24,10 @@ pub mod wavefunction;
 
 #[derive(Error, Debug)]
 pub enum SimquilError {
-    #[error("failed to execute program")]
+    #[error("Failed to execute program")]
     ExecutionError(#[from] vm::VMError),
+    #[error("Failed to read program")]
+    ReadProgramError(#[from] io::Error),
 }
 
 #[derive(StructOpt)]
@@ -35,11 +37,13 @@ struct Cli {
 
 fn run(cli: Cli) -> Result<()> {
     let quil = match cli.quil_file {
-        Some(path) => fs::read_to_string(path).expect("bad read"),
+        Some(path) => fs::read_to_string(path).map_err(SimquilError::ReadProgramError)?,
         None => {
             let mut stdin = io::stdin();
             let mut buf = String::new();
-            stdin.read_to_string(&mut buf).expect("bad read from stdin");
+            stdin
+                .read_to_string(&mut buf)
+                .map_err(SimquilError::ReadProgramError)?;
             buf
         }
     };
@@ -57,12 +61,14 @@ fn run(cli: Cli) -> Result<()> {
     println!();
     println!("{:?}", vm);
 
-    println!("Classical memory:");
-    println!();
-    vm.memory
-        .iter()
-        .sorted_by_key(|x| x.0)
-        .for_each(|(mref, data)| println!("{}[0..{}]:\t {:?}", mref, data.len(), data));
+    if !vm.memory.is_empty() {
+        println!("Classical memory:");
+        println!();
+        vm.memory
+            .iter()
+            .sorted_by_key(|x| x.0)
+            .for_each(|(mref, data)| println!("{}[0..{}]:\t {:?}", mref, data.len(), data));
+    }
 
     Ok(())
 }
diff --git a/src/matrix.rs b/src/matrix.rs
index b458f85..7a9d2f1 100644
--- a/src/matrix.rs
+++ b/src/matrix.rs
@@ -2,7 +2,7 @@ use num::complex::Complex64;
 use std::collections::HashMap;
 use thiserror::Error;
 
-use crate::gates::{gate_matrix, QGate};
+use crate::gates::{standard::gate_matrix, standard::GateError, QGate};
 
 use crate::matrix::InstructionMatrixError::{InvalidInstruction, InvalidQubit};
 use quil::expression::EvaluationError;
@@ -14,12 +14,14 @@ pub const I1: Complex64 = Complex64::new(0.0, 1.0);
 
 #[derive(Error, Debug)]
 pub enum InstructionMatrixError {
-    #[error("invalid parameter")]
+    #[error("Invalid parameter")]
     InvalidParameter(#[from] EvaluationError),
-    #[error("cannot create matrix from gate with variable qubit {0}")]
+    #[error("Cannot create matrix from gate with variable qubit {0}")]
     InvalidQubit(String),
-    #[error("cannot create matrix from non-gate instruction {0}")]
+    #[error("Cannot create matrix from non-gate instruction {0}")]
     InvalidInstruction(String),
+    #[error("Gate matrix failed")]
+    InvalidGateError(#[from] GateError),
 }
 
 pub fn instruction_matrix(instruction: Instruction) -> Result<QGate, InstructionMatrixError> {
@@ -30,7 +32,7 @@ pub fn instruction_matrix(instruction: Instruction) -> Result<QGate, Instruction
             qubits,
             modifiers: _,
         } => {
-            let params: Result<Vec<f64>, EvaluationError> = parameters
+            let params: Vec<f64> = parameters
                 .iter()
                 .map(
                     |p| match p.to_owned().evaluate_to_complex(&HashMap::new()) {
@@ -38,15 +40,15 @@ pub fn instruction_matrix(instruction: Instruction) -> Result<QGate, Instruction
                         Err(e) => Err(e),
                     },
                 )
-                .collect();
-            let qubits: Result<Vec<_>, _> = qubits
+                .collect::<Result<Vec<_>, _>>()?;
+            let qubits: Vec<usize> = qubits
                 .iter()
                 .map(|q| match q {
                     Qubit::Fixed(i) => Ok(*i as usize),
                     Qubit::Variable(q) => Err(InvalidQubit(q.clone())),
                 })
-                .collect();
-            Ok(gate_matrix(name, params?, qubits?))
+                .collect::<Result<Vec<_>, _>>()?;
+            Ok(gate_matrix(name, params, qubits)?)
         }
         instruction => Err(InvalidInstruction(instruction.to_string())),
     }
@@ -54,10 +56,8 @@ pub fn instruction_matrix(instruction: Instruction) -> Result<QGate, Instruction
 
 #[cfg(test)]
 pub mod test {
-    use crate::gates::QGate;
-    use crate::matrix::{instruction_matrix, InstructionMatrixError};
+    use crate::matrix::instruction_matrix;
     use quil::expression::Expression;
-    use quil::expression::Expression::Address;
     use quil::instruction::{Instruction, Qubit};
 
     #[test]