diff --git a/CHANGELOG.md b/CHANGELOG.md
index 7c3919e..c9a8926 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,6 +2,7 @@ QuantumAnnealing.jl Change Log
 ==============================
 
 ### Staged
+- Add tools for working with classical Ising models (#8)
 - Add data processing tools (#6)
 
 ### v0.0.1
diff --git a/docs/src/api.md b/docs/src/api.md
index c3d119a..9295455 100644
--- a/docs/src/api.md
+++ b/docs/src/api.md
@@ -32,6 +32,14 @@ AS_CIRCULAR
 AS_DW_QUADRATIC
 ```
 
+## Ising Functions
+```@autodocs
+Modules = [QuantumAnnealing]
+Pages   = ["ising.jl"]
+Order   = [:function]
+Private  = true
+```
+
 # D-Wave Functions
 ```@autodocs
 Modules = [QuantumAnnealing]
diff --git a/src/QuantumAnnealing.jl b/src/QuantumAnnealing.jl
index ea47414..92a8d81 100644
--- a/src/QuantumAnnealing.jl
+++ b/src/QuantumAnnealing.jl
@@ -7,6 +7,7 @@ module QuantumAnnealing
     import DifferentialEquations
 
     include("base.jl")
+    include("ising.jl")
 
     include("simulate.jl")
     include("simulate_de.jl")
diff --git a/src/dwave.jl b/src/dwave.jl
index d0e61b9..0a8aea9 100644
--- a/src/dwave.jl
+++ b/src/dwave.jl
@@ -274,19 +274,6 @@ function read_bqpjson(bqpjson_file::String)
 end
 
 
-function _eval_state_energy(state, ising_model)
-    energy = 0.0
-    for (ids,v) in ising_model
-        val = v
-        for qid in ids
-            val *= state[qid]
-        end
-        energy += val
-    end
-    return energy
-end
-
-
 function write_dwisc(outfile::String, ρ, ising_model, qubit_ids; simulated_num_reads=1e17, annealing_time=1000)
     dwisc_data = Dict()
     dwisc_data["variable_ids"] = qubit_ids
@@ -303,7 +290,7 @@ function write_dwisc(outfile::String, ρ, ising_model, qubit_ids; simulated_num_
     for state_int in 0:(2^n-1)
         prob = probs[state_int+1]
         spin_vector = int2spin(state_int, pad=n)
-        energy = _eval_state_energy(spin_vector, ising_model)
+        energy = eval_ising_state_energy(spin_vector, ising_model)
 
         sol_data = Dict(
             "energy" => energy,
diff --git a/src/ising.jl b/src/ising.jl
new file mode 100644
index 0000000..9abc0b4
--- /dev/null
+++ b/src/ising.jl
@@ -0,0 +1,79 @@
+### Helper Functions for Classical Ising Models ###
+
+"""
+given a state vector of spin values and an Ising model
+computes the energy of that spin configuration
+"""
+function eval_ising_state_energy(spin_state::Vector, ising_model::Dict)
+    energy = 0.0
+    for (ids,v) in ising_model
+        val = v
+        for qid in ids
+            val *= spin_state[qid]
+        end
+        energy += val
+    end
+    return energy
+end
+
+"""
+given an Ising model computes a mapping from state integers to energy values
+"""
+function compute_ising_state_energies(ising_model::Dict)
+    n = _check_ising_model_ids(ising_model)
+
+    state_energy = Dict{Int,Float64}()
+    for state_int in 0:(2^n-1)
+        spin_vector = binary2spin(int2binary(state_int, pad=n))
+        energy = eval_ising_state_energy(spin_vector, ising_model)
+        state_energy[state_int] = energy
+    end
+
+    return state_energy
+end
+
+"""
+given an Ising model computes a mapping from energy values to collections of state integers
+"""
+function compute_ising_energy_levels(ising_model::Dict)
+    state_energies = compute_ising_state_energies(ising_model)
+
+    energy_levels = Dict{Float64,Set{Int}}()
+    for (state_id, energy) in state_energies
+        if !haskey(energy_levels, energy)
+            energy_levels[energy] = Set{Int}()
+        end
+        push!(energy_levels[energy], state_id)
+    end
+
+    return energy_levels
+end
+
+"""
+given an Ising model prints state strings by ascending energy levels.
+`limit` is used the stop the printing after a number of energy levels have been
+presented (`limit` <= 0 will print all states).
+"""
+function print_ising_energy_levels(ising_model::Dict; limit=3)
+    n = maximum(max(k...) for (k,v) in ising_model)
+
+    energy_levels = compute_ising_energy_levels(ising_model)
+
+    energies = sort(collect(keys(energy_levels)))
+
+    for (i,energy) in enumerate(energies)
+        state_ids = energy_levels[energy]
+
+        println("\033[1menergy: $(energy)\033[0m")
+        for state_id in sort(collect(state_ids))
+            state = binary2spin(int2binary(state_id, pad=n))
+            state_string = spin2braket(state)
+            println("  $(state_string)")
+        end
+
+        if limit > 0 && i >= limit
+            println("first $(i) of $(length(energy_levels)) energy levels shown")
+            break
+        end
+    end
+end
diff --git a/test/base.jl b/test/base.jl
index a50464a..3fb51cb 100644
--- a/test/base.jl
+++ b/test/base.jl
@@ -69,6 +69,42 @@
 
 end
 
+@testset "ising energy computations" begin
+
+    @testset "single qubit, single state" begin
+        @test isapprox(eval_ising_state_energy([1],Dict((1,) => 1)), 1)
+        @test isapprox(eval_ising_state_energy([-1],Dict((1,) => 1)), -1)
+    end
+
+    @testset "single qubit, all states" begin
+        energies = compute_ising_state_energies(Dict((1,) => 1))
+        @test isapprox(energies[0], 1)
+        @test isapprox(energies[1], -1)
+    end
+
+    @testset "two qubit, single state" begin
+        ising_model = Dict((1,) => 1, (2,) => 1, (1,2) => -1)
+        @test isapprox(eval_ising_state_energy([1, 1], ising_model), 1)
+        @test isapprox(eval_ising_state_energy([1, -1], ising_model), 1)
+        @test isapprox(eval_ising_state_energy([-1, 1], ising_model), 1)
+        @test isapprox(eval_ising_state_energy([-1, -1], ising_model), -3)
+    end
+
+    @testset "two qubit, all state energies" begin
+        energies = compute_ising_state_energies(Dict((1,) => 1, (1,2) => -1))
+        @test isapprox(energies[0], 0)
+        @test isapprox(energies[1], 0)
+        @test isapprox(energies[2], 2)
+        @test isapprox(energies[3], -2)
+    end
+
+    @testset "two qubit, energy levels" begin
+        energy_levels = compute_ising_energy_levels(Dict((1,2) => -1))
+        @test energy_levels[-1.0] == Set([0,3])
+        @test energy_levels[1.0] == Set([1,2])
+    end
+
+end
 
 @testset "csv annealing schedules" begin
     s_100 = range(0, 1, length=100)
diff --git a/test/simulate.jl b/test/simulate.jl
index 8ffe2e1..df62b1b 100644
--- a/test/simulate.jl
+++ b/test/simulate.jl
@@ -446,8 +446,8 @@ end
         steps = 100
         numshots = 10
 
-        #random numbers between -1.0 and 1.0
-        z_bias = (Random.rand(numshots) .- 0.5)*2
+        # random numbers between -0.75 and 0.75
+        z_bias = (Random.rand(numshots) .- 0.5)*1.5
 
         bqpjson_file = "data/bqpjson_1q.json"
         dwisc_file = "tmp.json"
@@ -467,9 +467,9 @@ end
         steps = 100
         numshots = 10
 
-        #random numbers between -1.0 and 1.0
-        x_bias = (Random.rand(numshots) .- 0.5)*2
-        z_bias = (Random.rand(numshots) .- 0.5)*2
+        # random numbers between -0.75 and 0.75
+        x_bias = (Random.rand(numshots) .- 0.5)*1.5
+        z_bias = (Random.rand(numshots) .- 0.5)*1.5
 
         bqpjson_file = "data/bqpjson_1q.json"
         dwisc_file = "tmp.json"