diff --git a/.gitignore b/.gitignore
index 94f1119..c491b10 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,2 +1,3 @@
 .DS_Store
 .vscode
+output/
diff --git a/Manifest.toml b/Manifest.toml
new file mode 100644
index 0000000..c4a3262
--- /dev/null
+++ b/Manifest.toml
@@ -0,0 +1,788 @@
+# This file is machine-generated - editing it directly is not advised
+
+[[AbstractFFTs]]
+deps = ["LinearAlgebra"]
+git-tree-sha1 = "485ee0867925449198280d4af84bdb46a2a404d0"
+uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c"
+version = "1.0.1"
+
+[[Adapt]]
+deps = ["LinearAlgebra"]
+git-tree-sha1 = "84918055d15b3114ede17ac6a7182f68870c16f7"
+uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
+version = "3.3.1"
+
+[[ArgTools]]
+uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
+
+[[ArnoldiMethod]]
+deps = ["LinearAlgebra", "Random", "StaticArrays"]
+git-tree-sha1 = "f87e559f87a45bece9c9ed97458d3afe98b1ebb9"
+uuid = "ec485272-7323-5ecc-a04f-4719b315124d"
+version = "0.1.0"
+
+[[Artifacts]]
+uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
+
+[[AxisAlgorithms]]
+deps = ["LinearAlgebra", "Random", "SparseArrays", "WoodburyMatrices"]
+git-tree-sha1 = "a4d07a1c313392a77042855df46c5f534076fab9"
+uuid = "13072b0f-2c55-5437-9ae7-d433b7a33950"
+version = "1.0.0"
+
+[[Base64]]
+uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
+
+[[CSVFiles]]
+deps = ["CodecZlib", "DataValues", "FileIO", "HTTP", "IterableTables", "IteratorInterfaceExtensions", "TableShowUtils", "TableTraits", "TableTraitsUtils", "TextParse"]
+git-tree-sha1 = "d4dd66b73d3c811daa67587980bf45a179d16983"
+uuid = "5d742f6a-9f54-50ce-8119-2520741973ca"
+version = "1.0.1"
+
+[[ChainRulesCore]]
+deps = ["Compat", "LinearAlgebra", "SparseArrays"]
+git-tree-sha1 = "dcc25ff085cf548bc8befad5ce048391a7c07d40"
+uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
+version = "0.10.11"
+
+[[Classes]]
+deps = ["DataStructures", "InteractiveUtils", "MacroTools"]
+git-tree-sha1 = "96a8a88e6a59e870a47f415ad1519b6b4cb0fc02"
+uuid = "1a9c1350-211b-5766-99cd-4544d885a0d1"
+version = "1.4.1"
+
+[[CodecZlib]]
+deps = ["TranscodingStreams", "Zlib_jll"]
+git-tree-sha1 = "ded953804d019afa9a3f98981d99b33e3db7b6da"
+uuid = "944b1d66-785c-5afd-91f1-9de20f533193"
+version = "0.7.0"
+
+[[ColorTypes]]
+deps = ["FixedPointNumbers", "Random"]
+git-tree-sha1 = "32a2b8af383f11cbb65803883837a149d10dfe8a"
+uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
+version = "0.10.12"
+
+[[Colors]]
+deps = ["ColorTypes", "FixedPointNumbers", "Reexport"]
+git-tree-sha1 = "417b0ed7b8b838aa6ca0a87aadf1bb9eb111ce40"
+uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
+version = "0.12.8"
+
+[[Combinatorics]]
+git-tree-sha1 = "08c8b6831dc00bfea825826be0bc8336fc369860"
+uuid = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
+version = "1.0.2"
+
+[[Compat]]
+deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
+git-tree-sha1 = "dc7dedc2c2aa9faf59a55c622760a25cbefbe941"
+uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
+version = "3.31.0"
+
+[[CompilerSupportLibraries_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
+
+[[Compose]]
+deps = ["Base64", "Colors", "DataStructures", "Dates", "IterTools", "JSON", "LinearAlgebra", "Measures", "Printf", "Random", "Requires", "Statistics", "UUIDs"]
+git-tree-sha1 = "c6461fc7c35a4bb8d00905df7adafcff1fe3a6bc"
+uuid = "a81c6b42-2e10-5240-aca2-a61377ecd94b"
+version = "0.9.2"
+
+[[ConstructionBase]]
+deps = ["LinearAlgebra"]
+git-tree-sha1 = "f74e9d5388b8620b4cee35d4c5a618dd4dc547f4"
+uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
+version = "1.3.0"
+
+[[Crayons]]
+git-tree-sha1 = "3f71217b538d7aaee0b69ab47d9b7724ca8afa0d"
+uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
+version = "4.0.4"
+
+[[DataAPI]]
+git-tree-sha1 = "ee400abb2298bd13bfc3df1c412ed228061a2385"
+uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
+version = "1.7.0"
+
+[[DataFrames]]
+deps = ["Compat", "DataAPI", "Future", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrettyTables", "Printf", "REPL", "Reexport", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"]
+git-tree-sha1 = "1dadfca11c0e08e03ab15b63aaeda55266754bad"
+uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+version = "1.2.0"
+
+[[DataStructures]]
+deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
+git-tree-sha1 = "4437b64df1e0adccc3e5d1adbc3ac741095e4677"
+uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
+version = "0.18.9"
+
+[[DataValueInterfaces]]
+git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
+uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
+version = "1.0.0"
+
+[[DataValues]]
+deps = ["DataValueInterfaces", "Dates"]
+git-tree-sha1 = "d88a19299eba280a6d062e135a43f00323ae70bf"
+uuid = "e7dc6d0d-1eca-5fa6-8ad6-5aecde8b7ea5"
+version = "0.4.13"
+
+[[Dates]]
+deps = ["Printf"]
+uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
+
+[[DelimitedFiles]]
+deps = ["Mmap"]
+uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"
+
+[[Distributed]]
+deps = ["Random", "Serialization", "Sockets"]
+uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
+
+[[Distributions]]
+deps = ["FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SparseArrays", "SpecialFunctions", "Statistics", "StatsBase", "StatsFuns"]
+git-tree-sha1 = "a837fdf80f333415b69684ba8e8ae6ba76de6aaa"
+uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
+version = "0.24.18"
+
+[[DocStringExtensions]]
+deps = ["LibGit2"]
+git-tree-sha1 = "a32185f5428d3986f47c2ab78b1f216d5e6cc96f"
+uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+version = "0.8.5"
+
+[[DoubleFloats]]
+deps = ["GenericLinearAlgebra", "LinearAlgebra", "Polynomials", "Printf", "Quadmath", "Random", "Requires", "SpecialFunctions"]
+git-tree-sha1 = "1c962cf7e75c09a5f1fbf504df7d6a06447a1129"
+uuid = "497a8b3b-efae-58df-a0af-a86822472b78"
+version = "1.1.23"
+
+[[Downloads]]
+deps = ["ArgTools", "LibCURL", "NetworkOptions"]
+uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
+
+[[Electron]]
+deps = ["Base64", "FilePaths", "JSON", "Pkg", "Sockets", "URIParser", "UUIDs"]
+git-tree-sha1 = "a53025d3eabe23659065b3c5bba7b4ffb1327aa0"
+uuid = "a1bb12fb-d4d1-54b4-b10a-ee7951ef7ad3"
+version = "3.1.2"
+
+[[ExprTools]]
+git-tree-sha1 = "b7e3d17636b348f005f11040025ae8c6f645fe92"
+uuid = "e2ba6199-217a-4e67-a87a-7c52f15ade04"
+version = "0.1.6"
+
+[[FFTW]]
+deps = ["AbstractFFTs", "FFTW_jll", "LinearAlgebra", "MKL_jll", "Preferences", "Reexport"]
+git-tree-sha1 = "f985af3b9f4e278b1d24434cbb546d6092fca661"
+uuid = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
+version = "1.4.3"
+
+[[FFTW_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "3676abafff7e4ff07bbd2c42b3d8201f31653dcc"
+uuid = "f5851436-0d7a-5f13-b9de-f02708fd171a"
+version = "3.3.9+8"
+
+[[FileIO]]
+deps = ["Pkg", "Requires", "UUIDs"]
+git-tree-sha1 = "256d8e6188f3f1ebfa1a5d17e072a0efafa8c5bf"
+uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
+version = "1.10.1"
+
+[[FilePaths]]
+deps = ["FilePathsBase", "MacroTools", "Reexport", "Requires"]
+git-tree-sha1 = "919d9412dbf53a2e6fe74af62a73ceed0bce0629"
+uuid = "8fc22ac5-c921-52a6-82fd-178b2807b824"
+version = "0.8.3"
+
+[[FilePathsBase]]
+deps = ["Dates", "Mmap", "Printf", "Test", "UUIDs"]
+git-tree-sha1 = "0f5e8d0cb91a6386ba47bd1527b240bd5725fbae"
+uuid = "48062228-2e41-5def-b9a4-89aafe57970f"
+version = "0.9.10"
+
+[[FillArrays]]
+deps = ["LinearAlgebra", "Random", "SparseArrays"]
+git-tree-sha1 = "693210145367e7685d8604aee33d9bfb85db8b31"
+uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
+version = "0.11.9"
+
+[[FixedPointNumbers]]
+deps = ["Statistics"]
+git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
+uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
+version = "0.8.4"
+
+[[Formatting]]
+deps = ["Printf"]
+git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
+uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
+version = "0.4.2"
+
+[[Future]]
+deps = ["Random"]
+uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"
+
+[[GenericLinearAlgebra]]
+deps = ["LinearAlgebra", "Printf", "Random"]
+git-tree-sha1 = "ff291c1827030ffaacaf53e3c83ed92d4d5e6fb6"
+uuid = "14197337-ba66-59df-a3e3-ca00e7dcff7a"
+version = "0.2.5"
+
+[[GlobalSensitivityAnalysis]]
+deps = ["DataFrames", "DataStructures", "Distributed", "Distributions", "KernelDensity", "NumericalIntegration", "ProgressMeter", "Random", "Sobol", "Statistics", "StatsBase", "VegaLite"]
+git-tree-sha1 = "3fb26b1bf1ee9f897341ee1a6ff186da514e8666"
+uuid = "1b10255b-6da3-57ce-9089-d24e8517b87e"
+version = "1.1.2"
+
+[[GraphPlot]]
+deps = ["ArnoldiMethod", "ColorTypes", "Colors", "Compose", "DelimitedFiles", "LightGraphs", "LinearAlgebra", "Random", "SparseArrays"]
+git-tree-sha1 = "dd8f15128a91b0079dfe3f4a4a1e190e54ac7164"
+uuid = "a2cc645c-3eea-5389-862e-a155d0052231"
+version = "0.4.4"
+
+[[HTTP]]
+deps = ["Base64", "Dates", "IniFile", "Logging", "MbedTLS", "NetworkOptions", "Sockets", "URIs"]
+git-tree-sha1 = "c6a1fff2fd4b1da29d3dccaffb1e1001244d844e"
+uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3"
+version = "0.9.12"
+
+[[Inflate]]
+git-tree-sha1 = "f5fc07d4e706b84f72d54eedcc1c13d92fb0871c"
+uuid = "d25df0c9-e2be-5dd7-82c8-3ad0b3e990b9"
+version = "0.1.2"
+
+[[IniFile]]
+deps = ["Test"]
+git-tree-sha1 = "098e4d2c533924c921f9f9847274f2ad89e018b8"
+uuid = "83e8ac13-25f8-5344-8a64-a9f2b223428f"
+version = "0.5.0"
+
+[[IntelOpenMP_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "d979e54b71da82f3a65b62553da4fc3d18c9004c"
+uuid = "1d5cc7b8-4909-519e-a0f8-d0f5ad9712d0"
+version = "2018.0.3+2"
+
+[[InteractiveUtils]]
+deps = ["Markdown"]
+uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
+
+[[Interpolations]]
+deps = ["AxisAlgorithms", "ChainRulesCore", "LinearAlgebra", "OffsetArrays", "Random", "Ratios", "Requires", "SharedArrays", "SparseArrays", "StaticArrays", "WoodburyMatrices"]
+git-tree-sha1 = "1470c80592cf1f0a35566ee5e93c5f8221ebc33a"
+uuid = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
+version = "0.13.3"
+
+[[Intervals]]
+deps = ["Dates", "Printf", "RecipesBase", "Serialization", "TimeZones"]
+git-tree-sha1 = "323a38ed1952d30586d0fe03412cde9399d3618b"
+uuid = "d8418881-c3e1-53bb-8760-2df7ec849ed5"
+version = "1.5.0"
+
+[[InvertedIndices]]
+deps = ["Test"]
+git-tree-sha1 = "15732c475062348b0165684ffe28e85ea8396afc"
+uuid = "41ab1584-1d38-5bbf-9106-f11c6c58b48f"
+version = "1.0.0"
+
+[[IterTools]]
+git-tree-sha1 = "05110a2ab1fc5f932622ffea2a003221f4782c18"
+uuid = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
+version = "1.3.0"
+
+[[IterableTables]]
+deps = ["DataValues", "IteratorInterfaceExtensions", "Requires", "TableTraits", "TableTraitsUtils"]
+git-tree-sha1 = "70300b876b2cebde43ebc0df42bc8c94a144e1b4"
+uuid = "1c8ee90f-4401-5389-894e-7a04a3dc0f4d"
+version = "1.0.0"
+
+[[IteratorInterfaceExtensions]]
+git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
+uuid = "82899510-4779-5014-852e-03e436cf321d"
+version = "1.0.0"
+
+[[JLD2]]
+deps = ["DataStructures", "FileIO", "MacroTools", "Mmap", "Pkg", "Printf", "Reexport", "TranscodingStreams", "UUIDs"]
+git-tree-sha1 = "4813826871754cf52607e76ad37acb36ccf52719"
+uuid = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
+version = "0.4.11"
+
+[[JLLWrappers]]
+deps = ["Preferences"]
+git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e"
+uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
+version = "1.3.0"
+
+[[JSON]]
+deps = ["Dates", "Mmap", "Parsers", "Unicode"]
+git-tree-sha1 = "81690084b6198a2e1da36fcfda16eeca9f9f24e4"
+uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
+version = "0.21.1"
+
+[[JSONSchema]]
+deps = ["HTTP", "JSON", "ZipFile"]
+git-tree-sha1 = "b84ab8139afde82c7c65ba2b792fe12e01dd7307"
+uuid = "7d188eb4-7ad8-530c-ae41-71a32a6d4692"
+version = "0.3.3"
+
+[[KernelDensity]]
+deps = ["Distributions", "DocStringExtensions", "FFTW", "Interpolations", "StatsBase"]
+git-tree-sha1 = "591e8dc09ad18386189610acafb970032c519707"
+uuid = "5ab0869b-81aa-558d-bb23-cbf5423bbe9b"
+version = "0.6.3"
+
+[[LazyArtifacts]]
+deps = ["Artifacts", "Pkg"]
+uuid = "4af54fe1-eca0-43a8-85a7-787d91b784e3"
+
+[[LibCURL]]
+deps = ["LibCURL_jll", "MozillaCACerts_jll"]
+uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"
+
+[[LibCURL_jll]]
+deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
+uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"
+
+[[LibGit2]]
+deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
+uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"
+
+[[LibSSH2_jll]]
+deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
+uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"
+
+[[Libdl]]
+uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
+
+[[LightGraphs]]
+deps = ["ArnoldiMethod", "DataStructures", "Distributed", "Inflate", "LinearAlgebra", "Random", "SharedArrays", "SimpleTraits", "SparseArrays", "Statistics"]
+git-tree-sha1 = "432428df5f360964040ed60418dd5601ecd240b6"
+uuid = "093fc24a-ae57-5d10-9952-331d41423f4d"
+version = "1.3.5"
+
+[[LinearAlgebra]]
+deps = ["Libdl"]
+uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+
+[[LogExpFunctions]]
+deps = ["DocStringExtensions", "LinearAlgebra"]
+git-tree-sha1 = "7bd5f6565d80b6bf753738d2bc40a5dfea072070"
+uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
+version = "0.2.5"
+
+[[Logging]]
+uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
+
+[[MKL_jll]]
+deps = ["Artifacts", "IntelOpenMP_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "Pkg"]
+git-tree-sha1 = "c253236b0ed414624b083e6b72bfe891fbd2c7af"
+uuid = "856f044c-d86e-5d09-b602-aeab76dc8ba7"
+version = "2021.1.1+1"
+
+[[MacroTools]]
+deps = ["Markdown", "Random"]
+git-tree-sha1 = "6a8a2a625ab0dea913aba95c11370589e0239ff0"
+uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
+version = "0.5.6"
+
+[[Markdown]]
+deps = ["Base64"]
+uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
+
+[[MbedTLS]]
+deps = ["Dates", "MbedTLS_jll", "Random", "Sockets"]
+git-tree-sha1 = "1c38e51c3d08ef2278062ebceade0e46cefc96fe"
+uuid = "739be429-bea8-5141-9913-cc70e7f3736d"
+version = "1.0.3"
+
+[[MbedTLS_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"
+
+[[Measures]]
+git-tree-sha1 = "e498ddeee6f9fdb4551ce855a46f54dbd900245f"
+uuid = "442fdcdd-2543-5da2-b0f3-8c86c306513e"
+version = "0.3.1"
+
+[[MetaGraphs]]
+deps = ["JLD2", "LightGraphs", "Random"]
+git-tree-sha1 = "9cc756ced4fc5ecfff62605462c1001eec1cf324"
+uuid = "626554b9-1ddb-594c-aa3c-2596fe9399a5"
+version = "0.6.7"
+
+[[Mimi]]
+deps = ["CSVFiles", "Classes", "DataFrames", "DataStructures", "Dates", "DelimitedFiles", "Distributions", "Electron", "FileIO", "FilePaths", "GlobalSensitivityAnalysis", "GraphPlot", "IterTools", "IteratorInterfaceExtensions", "JSON", "LightGraphs", "LinearAlgebra", "Logging", "MacroTools", "MetaGraphs", "NamedArrays", "Pkg", "ProgressMeter", "Random", "Serialization", "Statistics", "StatsBase", "StringBuilders", "TableTraits", "VegaLite"]
+git-tree-sha1 = "265766488a3c7c4e45a2dbeaea0c2de9e02208f8"
+uuid = "e4e893b0-ee5e-52ea-8111-44b3bdec128c"
+version = "1.3.0"
+
+[[Missings]]
+deps = ["DataAPI"]
+git-tree-sha1 = "f8c673ccc215eb50fcadb285f522420e29e69e1c"
+uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
+version = "0.4.5"
+
+[[Mmap]]
+uuid = "a63ad114-7e13-5084-954f-fe012c677804"
+
+[[Mocking]]
+deps = ["ExprTools"]
+git-tree-sha1 = "916b850daad0d46b8c71f65f719c49957e9513ed"
+uuid = "78c3b35d-d492-501b-9361-3d52fe80e533"
+version = "0.7.1"
+
+[[MozillaCACerts_jll]]
+uuid = "14a3606d-f60d-562e-9121-12d972cd8159"
+
+[[MutableArithmetics]]
+deps = ["LinearAlgebra", "SparseArrays", "Test"]
+git-tree-sha1 = "3927848ccebcc165952dc0d9ac9aa274a87bfe01"
+uuid = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
+version = "0.2.20"
+
+[[NamedArrays]]
+deps = ["Combinatorics", "DataStructures", "DelimitedFiles", "InvertedIndices", "LinearAlgebra", "Random", "Requires", "SparseArrays", "Statistics"]
+git-tree-sha1 = "9ba8ddb0c06a08b1bad81b7120d13288e5d766fa"
+uuid = "86f7a689-2022-50b4-a561-43c23ac3c673"
+version = "0.9.5"
+
+[[NetworkOptions]]
+uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
+
+[[NodeJS]]
+deps = ["Pkg"]
+git-tree-sha1 = "905224bbdd4b555c69bb964514cfa387616f0d3a"
+uuid = "2bd173c7-0d6d-553b-b6af-13a54713934c"
+version = "1.3.0"
+
+[[Nullables]]
+git-tree-sha1 = "8f87854cc8f3685a60689d8edecaa29d2251979b"
+uuid = "4d1e1d77-625e-5b40-9113-a560ec7a8ecd"
+version = "1.0.0"
+
+[[NumericalIntegration]]
+deps = ["Interpolations", "LinearAlgebra", "Logging"]
+git-tree-sha1 = "2a4ef5fc235053f9747d59cfdee19bcb8ba1e833"
+uuid = "e7bfaba1-d571-5449-8927-abc22e82249b"
+version = "0.3.3"
+
+[[OffsetArrays]]
+deps = ["Adapt"]
+git-tree-sha1 = "2bf78c5fd7fa56d2bbf1efbadd45c1b8789e6f57"
+uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
+version = "1.10.2"
+
+[[OpenSpecFun_jll]]
+deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
+uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
+version = "0.5.5+0"
+
+[[OrderedCollections]]
+git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
+uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
+version = "1.4.1"
+
+[[PDMats]]
+deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse"]
+git-tree-sha1 = "4dd403333bcf0909341cfe57ec115152f937d7d8"
+uuid = "90014a1f-27ba-587c-ab20-58faa44d9150"
+version = "0.11.1"
+
+[[Parsers]]
+deps = ["Dates"]
+git-tree-sha1 = "c8abc88faa3f7a3950832ac5d6e690881590d6dc"
+uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
+version = "1.1.0"
+
+[[Pkg]]
+deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
+uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+
+[[Polynomials]]
+deps = ["Intervals", "LinearAlgebra", "MutableArithmetics", "RecipesBase"]
+git-tree-sha1 = "0bbfdcd8cda81b8144de4be8a67f5717e959a005"
+uuid = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
+version = "2.0.14"
+
+[[PooledArrays]]
+deps = ["DataAPI", "Future"]
+git-tree-sha1 = "cde4ce9d6f33219465b55162811d8de8139c0414"
+uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720"
+version = "1.2.1"
+
+[[Preferences]]
+deps = ["TOML"]
+git-tree-sha1 = "00cfd92944ca9c760982747e9a1d0d5d86ab1e5a"
+uuid = "21216c6a-2e73-6563-6e65-726566657250"
+version = "1.2.2"
+
+[[PrettyTables]]
+deps = ["Crayons", "Formatting", "Markdown", "Reexport", "Tables"]
+git-tree-sha1 = "0d1245a357cc61c8cd61934c07447aa569ff22e6"
+uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
+version = "1.1.0"
+
+[[Printf]]
+deps = ["Unicode"]
+uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+
+[[ProgressMeter]]
+deps = ["Distributed", "Printf"]
+git-tree-sha1 = "afadeba63d90ff223a6a48d2009434ecee2ec9e8"
+uuid = "92933f4c-e287-5a05-a399-4b506db050ca"
+version = "1.7.1"
+
+[[QuadGK]]
+deps = ["DataStructures", "LinearAlgebra"]
+git-tree-sha1 = "12fbe86da16df6679be7521dfb39fbc861e1dc7b"
+uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
+version = "2.4.1"
+
+[[Quadmath]]
+deps = ["Printf", "Random", "Requires"]
+git-tree-sha1 = "5a8f74af8eae654086a1d058b4ec94ff192e3de0"
+uuid = "be4d8f0f-7fa4-5f49-b795-2f01399ab2dd"
+version = "0.5.5"
+
+[[REPL]]
+deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
+uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
+
+[[Random]]
+deps = ["Serialization"]
+uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+
+[[Ratios]]
+git-tree-sha1 = "37d210f612d70f3f7d57d488cb3b6eff56ad4e41"
+uuid = "c84ed2f1-dad5-54f0-aa8e-dbefe2724439"
+version = "0.4.0"
+
+[[RecipesBase]]
+git-tree-sha1 = "b3fb709f3c97bfc6e948be68beeecb55a0b340ae"
+uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
+version = "1.1.1"
+
+[[Reexport]]
+git-tree-sha1 = "5f6c21241f0f655da3952fd60aa18477cf96c220"
+uuid = "189a3867-3050-52da-a836-e630ba90ab69"
+version = "1.1.0"
+
+[[Requires]]
+deps = ["UUIDs"]
+git-tree-sha1 = "4036a3bd08ac7e968e27c203d45f5fff15020621"
+uuid = "ae029012-a4dd-5104-9daa-d747884805df"
+version = "1.1.3"
+
+[[Rmath]]
+deps = ["Random", "Rmath_jll"]
+git-tree-sha1 = "bf3188feca147ce108c76ad82c2792c57abe7b1f"
+uuid = "79098fc4-a85e-5d69-aa6a-4863f24498fa"
+version = "0.7.0"
+
+[[Rmath_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "68db32dff12bb6127bac73c209881191bf0efbb7"
+uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f"
+version = "0.3.0+0"
+
+[[SHA]]
+uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
+
+[[Serialization]]
+uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
+
+[[Setfield]]
+deps = ["ConstructionBase", "Future", "MacroTools", "Requires"]
+git-tree-sha1 = "d5640fc570fb1b6c54512f0bd3853866bd298b3e"
+uuid = "efcf1570-3423-57d1-acb7-fd33fddbac46"
+version = "0.7.0"
+
+[[SharedArrays]]
+deps = ["Distributed", "Mmap", "Random", "Serialization"]
+uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"
+
+[[SimpleTraits]]
+deps = ["InteractiveUtils", "MacroTools"]
+git-tree-sha1 = "5d7e3f4e11935503d3ecaf7186eac40602e7d231"
+uuid = "699a6c99-e7fa-54fc-8d76-47d257e15c1d"
+version = "0.9.4"
+
+[[Sobol]]
+deps = ["DelimitedFiles", "Random"]
+git-tree-sha1 = "5a74ac22a9daef23705f010f72c81d6925b19df8"
+uuid = "ed01d8cd-4d21-5b2a-85b4-cc3bdc58bad4"
+version = "1.5.0"
+
+[[Sockets]]
+uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
+
+[[SortingAlgorithms]]
+deps = ["DataStructures"]
+git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
+uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
+version = "1.0.1"
+
+[[SparseArrays]]
+deps = ["LinearAlgebra", "Random"]
+uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+
+[[SpecialFunctions]]
+deps = ["ChainRulesCore", "LogExpFunctions", "OpenSpecFun_jll"]
+git-tree-sha1 = "a50550fa3164a8c46747e62063b4d774ac1bcf49"
+uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
+version = "1.5.1"
+
+[[StaticArrays]]
+deps = ["LinearAlgebra", "Random", "Statistics"]
+git-tree-sha1 = "a43a7b58a6e7dc933b2fa2e0ca653ccf8bb8fd0e"
+uuid = "90137ffa-7385-5640-81b9-e52037218182"
+version = "1.2.6"
+
+[[Statistics]]
+deps = ["LinearAlgebra", "SparseArrays"]
+uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+
+[[StatsAPI]]
+git-tree-sha1 = "1958272568dc176a1d881acb797beb909c785510"
+uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
+version = "1.0.0"
+
+[[StatsBase]]
+deps = ["DataAPI", "DataStructures", "LinearAlgebra", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"]
+git-tree-sha1 = "2f6792d523d7448bbe2fec99eca9218f06cc746d"
+uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
+version = "0.33.8"
+
+[[StatsFuns]]
+deps = ["LogExpFunctions", "Rmath", "SpecialFunctions"]
+git-tree-sha1 = "30cd8c360c54081f806b1ee14d2eecbef3c04c49"
+uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
+version = "0.9.8"
+
+[[StringBuilders]]
+git-tree-sha1 = "4bb4facc29d6b2275043151f71c330edb647696b"
+uuid = "db12335b-fddc-5e1b-b0ee-42071d21ae50"
+version = "0.2.3"
+
+[[SuiteSparse]]
+deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
+uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"
+
+[[TOML]]
+deps = ["Dates"]
+uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
+
+[[TableShowUtils]]
+deps = ["DataValues", "Dates", "JSON", "Markdown", "Test"]
+git-tree-sha1 = "14c54e1e96431fb87f0d2f5983f090f1b9d06457"
+uuid = "5e66a065-1f0a-5976-b372-e0b8c017ca10"
+version = "0.2.5"
+
+[[TableTraits]]
+deps = ["IteratorInterfaceExtensions"]
+git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
+uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
+version = "1.0.1"
+
+[[TableTraitsUtils]]
+deps = ["DataValues", "IteratorInterfaceExtensions", "Missings", "TableTraits"]
+git-tree-sha1 = "8fc12ae66deac83e44454e61b02c37b326493233"
+uuid = "382cd787-c1b6-5bf2-a167-d5b971a19bda"
+version = "1.0.1"
+
+[[Tables]]
+deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"]
+git-tree-sha1 = "8ed4a3ea724dac32670b062be3ef1c1de6773ae8"
+uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
+version = "1.4.4"
+
+[[Tar]]
+deps = ["ArgTools", "SHA"]
+uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"
+
+[[Test]]
+deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
+uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[[TextParse]]
+deps = ["CodecZlib", "DataStructures", "Dates", "DoubleFloats", "Mmap", "Nullables", "WeakRefStrings"]
+git-tree-sha1 = "af728c38c839aee693637e15e244074a02f16c68"
+uuid = "e0df1984-e451-5cb5-8b61-797a481e67e3"
+version = "1.0.1"
+
+[[TimeZones]]
+deps = ["Dates", "Future", "LazyArtifacts", "Mocking", "Pkg", "Printf", "RecipesBase", "Serialization", "Unicode"]
+git-tree-sha1 = "81753f400872e5074768c9a77d4c44e70d409ef0"
+uuid = "f269a46b-ccf7-5d73-abea-4c690281aa53"
+version = "1.5.6"
+
+[[TranscodingStreams]]
+deps = ["Random", "Test"]
+git-tree-sha1 = "7c53c35547de1c5b9d46a4797cf6d8253807108c"
+uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
+version = "0.9.5"
+
+[[URIParser]]
+deps = ["Unicode"]
+git-tree-sha1 = "53a9f49546b8d2dd2e688d216421d050c9a31d0d"
+uuid = "30578b45-9adc-5946-b283-645ec420af67"
+version = "0.4.1"
+
+[[URIs]]
+git-tree-sha1 = "97bbe755a53fe859669cd907f2d96aee8d2c1355"
+uuid = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4"
+version = "1.3.0"
+
+[[UUIDs]]
+deps = ["Random", "SHA"]
+uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
+
+[[Unicode]]
+uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
+
+[[Vega]]
+deps = ["DataStructures", "DataValues", "Dates", "FileIO", "FilePaths", "IteratorInterfaceExtensions", "JSON", "JSONSchema", "MacroTools", "NodeJS", "Pkg", "REPL", "Random", "Setfield", "TableTraits", "TableTraitsUtils", "URIParser"]
+git-tree-sha1 = "43f83d3119a868874d18da6bca0f4b5b6aae53f7"
+uuid = "239c3e63-733f-47ad-beb7-a12fde22c578"
+version = "2.3.0"
+
+[[VegaLite]]
+deps = ["Base64", "DataStructures", "DataValues", "Dates", "FileIO", "FilePaths", "IteratorInterfaceExtensions", "JSON", "MacroTools", "NodeJS", "Pkg", "REPL", "Random", "TableTraits", "TableTraitsUtils", "URIParser", "Vega"]
+git-tree-sha1 = "3e23f28af36da21bfb4acef08b144f92ad205660"
+uuid = "112f6efa-9a02-5b7d-90c0-432ed331239a"
+version = "2.6.0"
+
+[[WeakRefStrings]]
+deps = ["DataAPI", "Random", "Test"]
+git-tree-sha1 = "28807f85197eaad3cbd2330386fac1dcb9e7e11d"
+uuid = "ea10d353-3f73-51f8-a26c-33c1cb351aa5"
+version = "0.6.2"
+
+[[WoodburyMatrices]]
+deps = ["LinearAlgebra", "SparseArrays"]
+git-tree-sha1 = "59e2ad8fd1591ea019a5259bd012d7aee15f995c"
+uuid = "efce3f68-66dc-5838-9240-27a6d6f5f9b6"
+version = "0.5.3"
+
+[[ZipFile]]
+deps = ["Libdl", "Printf", "Zlib_jll"]
+git-tree-sha1 = "c3a5637e27e914a7a445b8d0ad063d701931e9f7"
+uuid = "a5390f91-8eb1-5f08-bee0-b1d1ffed6cea"
+version = "0.9.3"
+
+[[Zlib_jll]]
+deps = ["Libdl"]
+uuid = "83775a58-1f1d-513f-b197-d71354ab007a"
+
+[[nghttp2_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"
+
+[[p7zip_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
diff --git a/README.md b/README.md
index 5085e54..7a75566 100644
--- a/README.md
+++ b/README.md
@@ -1,17 +1,15 @@
 # MimiFAIRv2.jl
 
-This is a work-in-progress respository for a Julia-Mimi implementation of the FAIRv2.0 simple climate model. The model description paper can be found at [FaIRv2.0.0: A Generalized Impulse Response Model for Climate Uncertainty and Future Scenario Exploration](https://gmd.copernicus.org/articles/14/3007/2021/gmd-14-3007-2021.html). 
+This is a work-in-progress repository for a Julia-Mimi implementation of the FAIRv2.0 simple climate model. The model description paper can be found at [FaIRv2.0.0: A Generalized Impulse Response Model for Climate Uncertainty and Future Scenario Exploration](https://gmd.copernicus.org/articles/14/3007/2021/gmd-14-3007-2021.html). 
 
 ## Preparing the Software Environment
 
-To add the package without necessarily downloading the full repository, run the following command at the julia package REPL:
+To add the package to your current environment, run the following command at the julia package REPL:
 
 ```julia
-pkg> add https://github.com/FrankErrickson/MimiFAIRv2.jl.git # add MimiFAIRv2 module on main branch
+pkg> add https://github.com/FrankErrickson/MimiFAIRv2.jl.git
 ```
 
-To work locally assuming you have downloaded the full repository and want to work with that version, likely only if you want to actively develope the MimiFAIRv2 package itself, you will later run the following command at the julia REPL, but do not need to call `add`.
-
 ```julia
 include("src/MimiFAIRv2.jl") # load the MimiFAIRv2 module 
 ```
@@ -21,17 +19,14 @@ You probably also want to install the Mimi package into your julia environment,
 ```julia
 pkg> add Mimi
 ```
+
 ## Running the Model
 
 The model uses the Mimi framework and it is highly recommended to read the Mimi  documentation first to understand the code structure. The basic way to access a copy of the default MimiFAIRv2 model and explore the resuts is the following:
 
 ```julia
 using Mimi 
-using MimiFAIRv2  # if you used the `pkg> add` syntax above and thus have the module in your environment
-
-# note that instead of the line above one could call the following line, which
-# may be particularly helpful for package development
-# include("src/MimiFAIRv2.jl")
+using MimiFAIRv2
 
 # Create an instance of MimiFAIRv2.
 m = MimiFAIRv2.get_model() 
diff --git a/test/runtests.jl b/test/runtests.jl
index e510133..da693d4 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,8 +1,9 @@
 using DataFrames
 using CSVFiles
 using Test
+using MimiFAIRv2
+using Mimi
 
-include(joinpath(@__DIR__, "..", "src", "MimiFAIRv2.jl"))
 using MimiFAIRv2: get_model # load `get_model` function to avoid need for `MimiFAIRv2.` prefix
 
 @testset "UI" begin
@@ -81,26 +82,103 @@ using MimiFAIRv2: get_model # load `get_model` function to avoid need for `MimiF
 
 end
 
-@testset "Python Comparison" begin 
+@testset "Python Comparison - With Python Replication Data" begin 
 
     # run model for each possible SSP and compare variable values between julia
     # and python versons
 
-    for SSP in ["ssp119", "ssp126", "ssp245", "ssp370", "ssp585"]
+    for ssp in ["ssp119", "ssp126", "ssp245", "ssp370", "ssp585"]
+
+        #Load initial conditions (just need parameter names) and extract gas names.
+        init_gas_vals     = DataFrame(load(joinpath(@__DIR__, "..", "data", "fair_initial_gas_cycle_conditions_1750.csv"), skiplines_begin=7))
+        montreal_init     = filter(:gas_group => ==("montreal"), init_gas_vals)
+        flourinated_init  = filter(:gas_group => ==("flourinated"), init_gas_vals)
+        aerosol_plus_init = filter(:gas_group => ==("aerosol_plus"), init_gas_vals)
+
+        # Sort arrays of initial conditions for multiple gases so they are listed alphabetically.
+        sort!(montreal_init,     :gas_name)
+        sort!(flourinated_init,  :gas_name)
+        sort!(aerosol_plus_init, :gas_name)
+
+        # Load replication emissions and forcing from Python.
+        python_emiss   = DataFrame(load(joinpath(@__DIR__, "..", "data", "python_replication_data", ssp*"_emissions.csv")))
+        python_exog_RF = DataFrame(load(joinpath(@__DIR__, "..", "data", "python_replication_data", ssp*"_forcing.csv")))[:,:External]
+
+        # Extract emissions arrays for multi-gas groupings.
+        montreal_emissions     = python_emiss[:, Symbol.(montreal_init.gas_name)]
+        flourinated_emissions  = python_emiss[:, Symbol.(flourinated_init.gas_name)]
+        aerosol_plus_emissions = python_emiss[:, Symbol.(aerosol_plus_init.gas_name)]
+
+        # Get a model instance for that SSP.
+        m = MimiFAIRv2.get_model(emissions_forcing_scenario = ssp, start_year=1750, end_year=2100)
+
+        # Update exogenous forcing and emissions data.
+        update_param!(m, :radiative_forcing, :exogenous_RF, python_exog_RF)
+        update_param!(m, :co2_cycle, :E_co2, python_emiss.carbon_dioxide)
+        update_param!(m, :ch4_cycle, :E_ch4, python_emiss.methane)
+        update_param!(m, :n2o_cycle, :E_n2o, python_emiss.nitrous_oxide)
+        update_param!(m, :flourinated_cycles, :E_flourinated, Array(flourinated_emissions))
+        update_param!(m, :montreal_cycles, :E_montreal, Array(montreal_emissions))
+        update_param!(m, :aerosol_plus_cycles, :E_aerosol_plus, Array(aerosol_plus_emissions))
 
-        m = get_model(;emissions_forcing_scenario = SSP, end_year = 2100)
         run(m)
 
-        # compare temperatures (TODO - this is failing)
-        julia_temps = m[:temperature, :T]
+        # Put CO₂ concentration, total radiative forcing, and global temperature anomaly into a dataframe.
+        julia_results = DataFrame(co2 = m[:co2_cycle, :co2], rf = m[:radiative_forcing, :total_RF], temp = m[:temperature, :T])
+
+        # compare temperatures
+        python_results = load(joinpath(@__DIR__, "..", "data", "python_replication_data", string(ssp, "_temperature.csv"))) |> DataFrame
+        @test maximum(abs, julia_results[!, :temp] .- python_results[!, :default]) ≈ 0. atol = 1e-3
+       
+        # compare radiative forcing
+        python_results = load(joinpath(@__DIR__, "..", "data", "python_replication_data", string(ssp, "_forcing.csv"))) |> DataFrame
+        select!(python_results, "Forcing component", "Total")
+        @test maximum(abs, julia_results[!, :rf] .- python_results[!, :Total]) ≈ 0. atol = 1e-3
 
-        data = load(joinpath(@__DIR__, "..", "data", "python_replication_data", string(SSP, "_temperature.csv"))) |> DataFrame
-        rename!(data, [:year, :T])
-        python_temps = data[!, :T]
+        # compare emissions
+        conv = select!(DataFrame(load(joinpath(@__DIR__, "..", "data", "default_gas_cycle_parameters.csv"), skiplines_begin=6)), :gas_name, :emis2conc)
+        filter(:gas_name  => ==("carbon_dioxide"), conv)
 
-        @test maximum(abs, julia_temps .- python_temps) ≈ 0. atol = 1e-3
+        python_results = load(joinpath(@__DIR__, "..", "data", "python_replication_data", string(ssp, "_concentrations.csv"))) |> DataFrame
+        select!(python_results, "Gas name", "carbon_dioxide")
+        @test maximum(abs, julia_results[!, :co2] .- python_results[!, :carbon_dioxide]) ≈ 0. atol = 1e-3
+
+    end
+end
+
+@testset "Python Comparison - With Default Model" begin 
+
+    # run model for each possible SSP and compare variable values between julia
+    # and python versons
+
+    for ssp in ["ssp119", "ssp126", "ssp245", "ssp370", "ssp585"]
+
+        # Get a model instance for that SSP.
+        m = MimiFAIRv2.get_model(emissions_forcing_scenario = ssp, start_year=1750, end_year=2100)
+        run(m)
 
-        # TODO - what other variables do we want to compare?
+        # Put CO₂ concentration, total radiative forcing, and global temperature anomaly into a dataframe.
+        julia_results = DataFrame(co2 = m[:co2_cycle, :co2], rf = m[:radiative_forcing, :total_RF], temp = m[:temperature, :T])
         
+        python_results = DataFrame(:year => collect(1750:2100))
+
+        # add temperature
+        res = load(joinpath(@__DIR__, "..", "data", "python_replication_data", string(ssp, "_temperature.csv"))) |> DataFrame
+        insertcols!(python_results, 2, :temp => res[!, :default])
+
+        # compare radiative forcing
+        res = load(joinpath(@__DIR__, "..", "data", "python_replication_data", string(ssp, "_forcing.csv"))) |> DataFrame
+        insertcols!(python_results, 2, :rf => res[!, :Total])
+
+        # compare emissions
+        res = load(joinpath(@__DIR__, "..", "data", "python_replication_data", string(ssp, "_concentrations.csv"))) |> DataFrame
+        insertcols!(python_results, 2, :co2 => res[!, :carbon_dioxide])
+
+        for name in [:co2, :rf, :temp]
+            maxdiff = maximum(abs, julia_results[!, name] .- python_results[!, name])
+            if maxdiff > 1e-3
+                @warn("Maximum absolute difference for $(ssp)'s $name is $maxdiff")
+            end
+        end
     end
 end
diff --git a/test/test_plot_comparison_code.jl b/test/test_plot_comparison_code.jl
new file mode 100644
index 0000000..3e79841
--- /dev/null
+++ b/test/test_plot_comparison_code.jl
@@ -0,0 +1,152 @@
+using DataFrames
+using CSVFiles
+using Mimi
+
+include(joinpath(@__DIR__, "..", "src/MimiFAIRv2.jl"))
+
+#Load initial conditions (just need parameter names) and extract gas names.
+init_gas_vals     = DataFrame(load(joinpath(@__DIR__, "..", "data", "fair_initial_gas_cycle_conditions_1750.csv"), skiplines_begin=7))
+montreal_init     = filter(:gas_group => ==("montreal"), init_gas_vals)
+flourinated_init  = filter(:gas_group => ==("flourinated"), init_gas_vals)
+aerosol_plus_init = filter(:gas_group => ==("aerosol_plus"), init_gas_vals)
+
+# Sort arrays of initial conditions for multiple gases so they are listed alphabetically.
+sort!(montreal_init,     :gas_name)
+sort!(flourinated_init,  :gas_name)
+sort!(aerosol_plus_init, :gas_name)
+
+# Set SSP (options = ['ssp119','ssp126','ssp245','ssp370','ssp585'])
+ssp = "ssp370"
+
+# Load replication emissions and forcing from Python.
+python_emiss   = DataFrame(load(joinpath(@__DIR__, "..", "data", "python_replication_data", ssp*"_emissions.csv")))
+python_exog_RF = DataFrame(load(joinpath(@__DIR__, "..", "data", "python_replication_data", ssp*"_forcing.csv")))[:,:External]
+
+# Extract emissions arrays for multi-gas groupings.
+montreal_emissions     = python_emiss[:, Symbol.(montreal_init.gas_name)]
+flourinated_emissions  = python_emiss[:, Symbol.(flourinated_init.gas_name)]
+aerosol_plus_emissions = python_emiss[:, Symbol.(aerosol_plus_init.gas_name)]
+
+# Get a model instance for that SSP.
+m = MimiFAIRv2.get_model(emissions_forcing_scenario=ssp, start_year=1750, end_year=2100)
+
+# Update exogenous forcing and emissions data.
+update_param!(m, :radiative_forcing, :exogenous_RF, python_exog_RF)
+update_param!(m, :co2_cycle, :E_co2, python_emiss.carbon_dioxide)
+update_param!(m, :ch4_cycle, :E_ch4, python_emiss.methane)
+update_param!(m, :n2o_cycle, :E_n2o, python_emiss.nitrous_oxide)
+update_param!(m, :flourinated_cycles, :E_flourinated, Array(flourinated_emissions))
+update_param!(m, :montreal_cycles, :E_montreal, Array(montreal_emissions))
+update_param!(m, :aerosol_plus_cycles, :E_aerosol_plus, Array(aerosol_plus_emissions))
+
+run(m)
+
+# Put CO₂ concentration, total radiative forcing, and global temperature anomaly into a dataframe.
+results = DataFrame(co2 = m[:co2_cycle, :co2], rf = m[:radiative_forcing, :total_RF], temp = m[:temperature, :T])
+
+# Save results.
+save("Mimi_FAIR_replication_"*ssp*".csv", results)
+
+
+
+#------------------------------------------------------
+#------------------------------------------------------
+#------------------------------------------------------
+# PLOTTING CODE IN R
+#------------------------------------------------------
+#------------------------------------------------------
+#------------------------------------------------------
+
+library(ggplot2)
+
+setwd("C:/Users/fce/Desktop/fair2.0/MimiFAIRv2.0.jl")
+
+# Load Mimi results
+mimi_ssp119 = read.csv("C:/Users/fce/Desktop/fair2.0/MimiFAIRv2.0.jl/DONT PUT ON GIT - FAIR REPLICATION OF PYTHON RESULTS/Mimi_FAIR_replication_ssp119.csv")
+mimi_ssp126 = read.csv("C:/Users/fce/Desktop/fair2.0/MimiFAIRv2.0.jl/DONT PUT ON GIT - FAIR REPLICATION OF PYTHON RESULTS/Mimi_FAIR_replication_ssp126.csv")
+mimi_ssp245 = read.csv("C:/Users/fce/Desktop/fair2.0/MimiFAIRv2.0.jl/DONT PUT ON GIT - FAIR REPLICATION OF PYTHON RESULTS/Mimi_FAIR_replication_ssp245.csv")
+mimi_ssp370 = read.csv("C:/Users/fce/Desktop/fair2.0/MimiFAIRv2.0.jl/DONT PUT ON GIT - FAIR REPLICATION OF PYTHON RESULTS/Mimi_FAIR_replication_ssp370.csv")
+mimi_ssp585 = read.csv("C:/Users/fce/Desktop/fair2.0/MimiFAIRv2.0.jl/DONT PUT ON GIT - FAIR REPLICATION OF PYTHON RESULTS/Mimi_FAIR_replication_ssp585.csv")
+
+# Load Python temperature results
+python_ssp119 = read.csv("C:/Users/fce/Desktop/fair2.0/MimiFAIRv2.0.jl/data/python_replication_data/ssp119_temperature.csv")
+python_ssp126 = read.csv("C:/Users/fce/Desktop/fair2.0/MimiFAIRv2.0.jl/data/python_replication_data/ssp126_temperature.csv")
+python_ssp245 = read.csv("C:/Users/fce/Desktop/fair2.0/MimiFAIRv2.0.jl/data/python_replication_data/ssp245_temperature.csv")
+python_ssp370 = read.csv("C:/Users/fce/Desktop/fair2.0/MimiFAIRv2.0.jl/data/python_replication_data/ssp370_temperature.csv")
+python_ssp585 = read.csv("C:/Users/fce/Desktop/fair2.0/MimiFAIRv2.0.jl/data/python_replication_data/ssp585_temperature.csv")
+
+# Put data into a dataframe.
+mimi_data   = data.frame(year=1750:2100, ssp119=mimi_ssp119[,3], ssp126=mimi_ssp126[,3], ssp245=mimi_ssp245[,3], ssp370=mimi_ssp370[,3], ssp585=mimi_ssp585[,3])
+python_data = data.frame(year=1750:2100, ssp119=python_ssp119[,2], ssp126=python_ssp126[,2], ssp245=python_ssp245[,2], ssp370=python_ssp370[,2], ssp585=python_ssp585[,2])
+
+# Crop it to all start in the year 2000 (then do historical in black).
+index_2000 = which(1750:2100 == 2000)
+
+short_mimi_data = mimi_data[(index_2000+1):351, ]
+short_python_data = python_data[(index_2000+1):351, ]
+
+# Get historicla data
+historic_mimi_data = mimi_data[1:(index_2000+1), ]
+historic_python_data = python_data[1:(index_2000+1), ]
+
+# Thin python data to every 5-years and to shortened time frame
+thin_indices = seq(1,length(short_python_data[,1]), by=3)
+thin_indices_historic = seq(1,length(historic_python_data[,1]), by=1)
+
+python_data_thinned = short_python_data[thin_indices, ]
+historic_python_data_thinned = historic_python_data[thin_indices_historic, ]
+
+ssp119_col = "#6ace85"
+ssp126_col = "#7f9be7"
+ssp245_col = "#9c4ed2"
+ssp370_col = "darkorange"
+#ssp370_col = "#9c4ed2"
+ssp585_col = "#ea3675"
+
+
+p = ggplot()
+
+p = p + geom_hline(yintercept=0, color="red", linetype="22")
+
+p = p + geom_line(data=historic_mimi_data, aes(x=year, y=ssp119), color="black")
+p = p + geom_point(data=historic_python_data_thinned, aes(x=year, y=ssp119), shape=21, size=1)
+
+p = p + geom_line(data=short_mimi_data, aes(x=year, y=ssp119), color=ssp119_col)
+p = p + geom_point(data=python_data_thinned, aes(x=year, y=ssp119), fill=ssp119_col, shape=21, size=1)
+
+p = p + geom_line(data=short_mimi_data, aes(x=year, y=ssp126), color=ssp126_col)
+p = p + geom_point(data=python_data_thinned, aes(x=year, y=ssp126), fill=ssp126_col, shape=21, size=1)
+
+p = p + geom_line(data=short_mimi_data, aes(x=year, y=ssp245), color=ssp245_col)
+p = p + geom_point(data=python_data_thinned, aes(x=year, y=ssp245), fill=ssp245_col, shape=21, size=1)
+
+p = p + geom_line(data=short_mimi_data, aes(x=year, y=ssp370), color=ssp370_col)
+p = p + geom_point(data=python_data_thinned, aes(x=year, y=ssp370), fill=ssp370_col, shape=21, size=1)
+
+p = p + geom_line(data=short_mimi_data, aes(x=year, y=ssp585), color=ssp585_col)
+p = p + geom_point(data=python_data_thinned, aes(x=year, y=ssp585), fill=ssp585_col, shape=21, size=1)
+
+p = p + xlab("Year")
+p = p + ylab("Global Surface Temperature Anomaly (K)")
+#p = p + labs(title = "Python vs. Julia-Mimi versions of FaIR v2.0")
+p = p + labs(title = "Python vs. Julia-Mimi versions of FaIR v2.0", subtitle = "RCMIP Emission & Forcing Scenarios (1750-2100)")
+
+p = p + theme(panel.background = element_rect(fill = "transparent"),
+                  plot.background = element_rect(fill = "transparent", color="NA"),
+                  panel.grid.minor = element_blank(),
+                  panel.grid.major = element_line(color="gray90", size=0.25),
+                  #axis.line.y = element_blank(),
+                  axis.line = element_line(color="black", size=0.25),
+                  axis.ticks = element_line(color="black", size=0.25),
+                  axis.text = element_text(size=9, colour="black"),
+                  axis.title = element_text(size=9, colour="black"),
+                  legend.position="none",
+                  #plot.title = element_blank(),
+                  #axis.title.y = element_blank(),
+                  axis.ticks.length=unit(.1, "cm"))
+                  #axis.text.y = element_blank(),
+                  #axis.ticks.y = element_blank())
+
+ggsave(p, file="Python_Mimi_FAIR2_temperature_comparison.jpg", device="jpeg", type="cairo", width=200, height=130, unit="mm", dpi=200)
+ggsave(p, file="Python_Mimi_FAIR2_temperature_comparison.pdf", device="pdf", width=200, height=130, unit="mm", dpi=200, useDingbats=FALSE)
+