First commit (yeah, its a mess)

.gitignore

@@ -0,0 +1,79 @@
+# vscode .gitignore
+.vscode/*
+!.vscode/settings.json
+!.vscode/tasks.json
+!.vscode/launch.json
+!.vscode/extensions.json
+!.vscode/*.code-snippets
+
+# Local History for Visual Studio Code
+.history/
+
+# Built Visual Studio Code Extensions
+*.vsix
+
+# julia .gitignore
+# Files generated by invoking Julia with --code-coverage
+*.jl.cov
+*.jl.*.cov
+
+# Files generated by invoking Julia with --track-allocation
+*.jl.mem
+
+# System-specific files and directories generated by the BinaryProvider and BinDeps packages
+# They contain absolute paths specific to the host computer, and so should not be committed
+deps/deps.jl
+deps/build.log
+deps/downloads/
+deps/usr/
+deps/src/
+
+# Build artifacts for creating documentation generated by the Documenter package
+docs/build/
+docs/site/
+
+# File generated by Pkg, the package manager, based on a corresponding Project.toml
+# It records a fixed state of all packages used by the project. As such, it should not be
+# committed for packages, but should be committed for applications that require a static
+# environment.
+Manifest.toml
+
+# fortran .gitignore
+# Prerequisites
+*.d
+
+# Compiled Object files
+*.slo
+*.lo
+*.o
+*.obj
+
+# Precompiled Headers
+*.gch
+*.pch
+
+# Compiled Dynamic libraries
+*.so
+*.dylib
+*.dll
+
+# Fortran module files
+*.mod
+*.smod
+
+# Compiled Static libraries
+*.lai
+*.la
+*.a
+*.lib
+
+# Executables
+*.exe
+*.out
+*.app
+
+# Project specific gitignores
+# Original notes from when developing
+notes/*
+# My terrible, newbie attempt at version control (yes, shouldve jsut learned git)
+backups/*

0-100abs.jl

@@ -0,0 +1,24 @@
+include("nanoconc.jl")
+using PyPlot
+
+particledata = ([2.5 1.],
+                [5. 1.],
+                [10. 1.],
+                [20. 1.],
+                [30. 1.],
+                [40. 1.],
+                [50. 1.],
+                [60. 1.],
+                [70. 1.],
+                [80. 1.],
+                [90. 1.],
+                [100. 1.])
+
+mat = nanoconc.loadmaterial("Gold",disp=false)
+predict(x) = nanoconc.abspredict(1.333,250.,900.,0x00001000,x,mat,1000.,1.)
+spectra = predict.(particledata)
+
+for spectrum in spectra
+    plot(spectrum[:,1],spectrum[:,2])
+end
+show()

0-100test.jl

@@ -0,0 +1,39 @@
+println("\nImporting modules...")
+
+include("nanoconc.jl")
+using PyPlot
+using JLD2
+
+println("Loading parameters...")
+
+particledata = ([2.5 1.],
+                [5. 1.],
+                [10. 1.],
+                [20. 1.],
+                [30. 1.],
+                [40. 1.],
+                [50. 1.],
+                [60. 1.],
+                [70. 1.],
+                [80. 1.],
+                [90. 1.],
+                [100. 1.])
+
+mat = nanoconc.loadmaterial("Gold",disp=false)
+predict(x) = nanoconc.qpredict(1.333,250.,900.,0x00001000,x,mat)
+println("Calculating spectra...")
+print("Calculation time: ")
+
+@time spectra = predict.(particledata)
+
+println("Saving output...")
+
+println("Plotting spectra...")
+
+for spectrum in spectra
+    plot(spectrum[:,1],spectrum[:,2])
+end
+println("Displaying spectra...")
+println("Script complete!")
+show()
+println()

20nm.jl

@@ -0,0 +1,28 @@
+println("\nImporting modules...")
+
+include("nanoconc.jl")
+using PyPlot
+
+println("Loading parameters...")
+
+particledata = [20. 1.]
+
+mat = nanoconc.loadmaterial("Gold",disp=false)
+predict(x) = nanoconc.abspredict(1.333,450.,650.,0x00001000,x,mat,1000.,1.)
+
+println("Calculating spectra...")
+print("Calculation time: ")
+
+@time spectrum = predict(particledata)
+
+println("Finding lambda max...")
+#println("\nDebug data:\nfindmax = $(findmax(spectrum[:,2]))\nmaxind = $(findmax(spectrum[:,2])[2])")
+lmax = spectrum[findmax(spectrum[:,2])[2],1]
+println("lambda max = $(lmax)")
+
+println("Plotting spectra...")
+
+plot(spectrum[:,1],spectrum[:,2])
+println("Displaying spectra...")
+show()
+println("Script complete!")

Project.toml

@@ -0,0 +1,9 @@
+[deps]
+CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
+Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
+Memoize = "c03570c3-d221-55d1-a50c-7939bbd78826"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
+XLSX = "fdbf4ff8-1666-58a4-91e7-1b58723a45e0"

abstest.jl

@@ -0,0 +1,37 @@
+println("\nImporting modules...")
+
+include("src/nanoconc.jl")
+using CSV
+using DataFrames
+using Profile, ProfileVega
+
+println("Loading parameters...")
+
+refmed = 1.333 # refractive index of the medium
+wavel1, wavel2 = 250., 900. # wavelengths to calculate between
+numval = 0x00001000 # number of values to calculate in spectrum
+particledata = [2.5 1.;
+                 5. 1.;
+                 10. 1.;
+                 20. 1.;
+                 30. 1.;
+                 40. 1.;
+                 50. 1.;
+                 60. 1.;
+                 70. 1.;
+                 80. 1.;
+                 90. 1.;
+                 100. 1.]
+materialdata = nanoconc.loadmaterial("Gold",disp=false)
+ppml = 10000. # particles per ml nanoconc.quantumcalc.numtomol(ppml) -> conc
+d0 = 1. # path length
+params = (refmed,wavel1,wavel2,numval,particledata,materialdata,ppml,d0)
+spectrum = nanoconc.abspredict(params...)
+
+println("Calculating absorbance spectrum...")
+print("Calculation time: ")
+
+@time spectrum = nanoconc.abspredict(params...)
+CSV.write("abs.csv", DataFrame(spectrum, :auto), writeheader=false)
+
+# @profview nanoconc.abspredict(params...)

src/miemfp.jl

@@ -0,0 +1,390 @@
+"""
+A module for calculating the expected extinction coefficients and backscattering
+parameters for coated/uncoated particles of a given size/material
+"""
+@fastmath module miemfp
+
+#############################################################################
+# STATUS NOTES:                                                             #
+# Synchronous parallelization implemented for bare calculations, notable    #
+# performance increase. Need to do same for coated                          #
+#############################################################################
+
+using Base.Threads
+using Distributed
+using Memoize
+using Interpolations
+
+const PI_5_996E3::Float64 = pi * 5.996e3
+const TWO_PI::Float64 = 2 * pi
+const BHCOAT_DEL::Float64 = 1e-8
+
+"Helper function for calculating om"
+@memoize function _om_calc(wavel::Float64)::Float64
+    PI_5_996E3 / wavel
+end
+
+"Helper function for calculating om0"
+@memoize function _om0r_calc(om0::Float64, fv::Float64, radcor::Float64)::Float64
+    om0 + (fv / (radcor * 1e7))
+end
+
+"Helper function for calculating om^2"
+@memoize function _om_sq_calc(om::Float64)::Float64
+    om^2
+end
+
+"Helper function for calculating om0^2"
+@memoize function _om0_sq_calc(om0::Float64)::Float64
+    om0^2
+end
+
+"Helper function for calculating omp^2"
+@memoize function _omp_sq_calc(omp::Float64)::Float64
+    omp^2
+end
+
+"Helper function for calculating om0r^2"
+@memoize function _om0r_sq_calc(om0r::Float64)::Float64
+    om0r^2
+end
+
+" function that performs manipulations to om necessary for efficient mfp"
+function _om_manipulations(wavel::Float64, fv::Float64, radcor::Float64,
+    omp::Float64, om0::Float64)::Tuple{Float64,Float64,Float64,Float64,Float64,Float64,Float64}
+
+    om = _om_calc(wavel)
+    om0r = _om0r_calc(om0, fv, radcor)
+    om_sq = _om_sq_calc(om)
+    om0_sq = _om0_sq_calc(om0)
+    omp_sq = _omp_sq_calc(omp)
+    om0r_sq = _om0r_sq_calc(om0r)
+    om_sq_plus_om0_sq = om_sq + om0_sq
+    return (om, om0r, om_sq, om0_sq, omp_sq, om0r_sq, om_sq_plus_om0_sq)
+end
+
+"""
+Corrects refractive index and attenuation coefficient values to account for
+the mean free-path effect on the free conductance electrons in the particles
+(translated from Haiss et als FORTRAN implementation)
+"""
+ function mfp(fv::Float64, wavel::Float64, radcor::Float64, omp::Float64,
+    om0::Float64, rn::Float64, rk::Float64)::Tuple{Float64,Float64,Float64}
+
+    om, om0r, om_sq, _, omp_sq, om0r_sq, om_sq_plus_om0_sq =
+        _om_manipulations(wavel, fv, radcor, omp, om0)
+
+    b1 = rn^2 - rk^2 - (1 - (omp_sq / om_sq_plus_om0_sq))
+    b2 = 2 * rn * rk - (omp_sq * om0 / (om * om_sq_plus_om0_sq))
+    a1r = 1 - (omp_sq / (om_sq + om0r_sq))
+    a2r = omp_sq * om0r / (om * (om_sq + om0r_sq))
+
+    r_const = sqrt((a1r / 2 + b1 / 2)^2 + (a2r / 2 + b2 / 2)^2)
+    rnr = sqrt((a1r + b1) / 2 + r_const)
+    rkr = sqrt(-(a1r + b1) / 2 + r_const)
+    return (rnr, rkr, om0r)
+end
+
+"""
+Solves for Qext accounting for the effect of particle coatings on extinction
+and scattering coefficients of particles (translated from Bohren/Huffman's
+FORTRAN77 implementation for coated particles)
+"""
+function bhcoat(x::Float64, y::Float64, rfrel1::ComplexF64, rfrel2::ComplexF64
+)::Tuple{Float64,Float64,Float64}
+    x1::ComplexF64 = rfrel1 * x
+    x2::ComplexF64 = rfrel2 * x
+    y2::ComplexF64 = rfrel2 * y
+    nstop::UInt64 = UInt64(round(y + 4 * cbrt(y) + 1))
+    refrel::ComplexF64 = rfrel2 / rfrel1
+
+    d0x1::ComplexF64 = cot(x1)
+    d0x2::ComplexF64 = cot(x2)
+    d0y2::ComplexF64 = cot(y2)
+    psi0y::Float64, psi1y::Float64 = cos(y), sin(y)
+    chi0y::Float64, chi1y::Float64 = -sin(y), cos(y)
+    # xi0y::ComplexF64 = ComplexF64(psi0y, -chi0y)
+    xi1y::ComplexF64 = ComplexF64(psi1y, -chi1y)
+    chi0y2::ComplexF64, chi1y2::ComplexF64 = -sin(y2), cos(y2)
+    chi0x2::ComplexF64, chi1x2::ComplexF64 = -sin(x2), cos(x2)
+
+    qsca::Float64, qext::Float64, xback::ComplexF64 = 0.0, 0.0, ComplexF64(0.0, 0.0)
+    iflag::Bool = false
+
+    @inbounds for n in (1:nstop)::UnitRange{Int64}
+        two_n_minus_1::Float64 = 2 * n - 1
+        psiy::Float64 = two_n_minus_1 * psi1y / y - psi0y
+        chiy::Float64 = two_n_minus_1 * chi1y / y - chi0y
+        xiy::ComplexF64 = ComplexF64(psiy, -chiy)
+        d1y2::ComplexF64 = 1.0 / (n / y2 - d0y2) - n / y2
+
+        if iflag == false
+            d1x1 = 1.0 / (n / x1 - d0x1) - n / x1
+            d1x2 = 1.0 / (n / x2 - d0x2) - n / x2
+            chix2 = two_n_minus_1 * chi1x2 / x2 - chi0x2
+            chiy2 = two_n_minus_1 * chi1y2 / y2 - chi0y2
+            chipx2::ComplexF64 = chi1x2 - n * chix2 / x2
+            chipy2 = chi1y2 - n * chiy2 / y2
+            rack_denominator = chix2 * d1x2 - chipx2
+            ancap::ComplexF64 = ((refrel * d1x1 - d1x2) / (refrel * d1x1 * chix2 - chipx2)) / rack_denominator
+            brack = ancap * (chiy2 * d1y2 - chipy2)
+            bncap::ComplexF64 = ((refrel * d1x2 - d1x1) / (refrel * chipx2 - d1x1 * chix2)) / rack_denominator
+            crack = bncap * (chiy2 * d1y2 - chipy2)
+            amess1::ComplexF64 = brack * chipy2
+            amess2::ComplexF64 = brack * chiy2
+            amess3::ComplexF64 = crack * chipy2
+            amess4::ComplexF64 = crack * chiy2
+            if !((abs(amess1) > BHCOAT_DEL * abs(d1y2))
+                 || (abs(amess2) > BHCOAT_DEL)
+                 || (abs(amess3) > BHCOAT_DEL * abs(d1y2))
+                 || (abs(amess4) > BHCOAT_DEL))
+                brack, crack = ComplexF64(0.0, 0.0), ComplexF64(0.0, 0.0)
+                iflag = true
+            end
+        end
+        dnbar::ComplexF64 = (d1y2 - brack * chipy2) / (1.0 - brack * chiy2)
+        gnbar::ComplexF64 = (d1y2 - crack * chipy2) / (1.0 - crack * chiy2)
+        n_over_y::Float64 = n / y
+        xiy_minus_xi1y::ComplexF64 = xiy - xi1y
+        psiy_minus_psi1y::Float64 = psiy - psi1y
+        an_mult::ComplexF64 = dnbar / rfrel2 + n_over_y
+        an::ComplexF64 = (an_mult * psiy_minus_psi1y) / (an_mult * xiy_minus_xi1y)
+        bn_mult::ComplexF64 = rfrel2 * gnbar + n_over_y
+        bn::ComplexF64 = (bn_mult * psiy_minus_psi1y) / (bn_mult * xiy_minus_xi1y)
+        two_n_plus_1::Float64 = 2 * n + 1
+        qsca += two_n_plus_1 * (abs(an)^2 + abs(bn)^2)
+        xback_sign::Float64 = if iseven(n)
+            1.0
+        else
+            -1.0
+        end
+        xback += xback_sign * two_n_plus_1 * (an - bn)
+        qext += two_n_plus_1 * (real(an) + real(bn))
+        psi0y, psi1y = psi1y, psiy
+        chi0y, chi1y = chi1y, chiy
+        xi1y = psi1y - chi1y * im
+        chi0x2, chi1x2 = chi1x2, chix2
+        chi0y2, chi1y2 = chi1y2, chiy2
+        d0x1, d0x2, d0y2 = d1x1, d1x2, d1y2
+    end
+    y_sq = y^2
+    two_over_y_sq = 2.0 / y_sq
+    qsca = two_over_y_sq * qsca
+    qext = two_over_y_sq * qext
+    qback = real((1 / y_sq) * (xback * conj(xback)))
+    return (qext, qsca, qback)
+end
+
+"""
+Finds scattering parameters for uncoated particles. Heavily modified, but originally
+based on a combination of the original FORTRAN77 BHMIE algorithm and a Python
+implementation attributed to Herbert Kaiser hosted on ScatterLib (http://scatterlib.wikidot.com/mie)
+"""
+function bhmie(x::Float64, refrel::ComplexF64, nang::UInt32
+)::Tuple{Float64,Float64,Float64,Array{ComplexF64,1},Array{ComplexF64,1}}
+    y::ComplexF64 = x * refrel
+    nstop::UInt32 = UInt32(round(x + 4.0 * cbrt(x) + 2.0))
+    nn::UInt32 = UInt32(round(max(nstop, abs(y)) + 14))
+    amu::Vector{Float64} = cos.((1.570796327 / Float64(nang - 1)).*Float64.(0:nang-1))    
+    
+    d = Vector{ComplexF64}(undef, nn)
+    d[nn] = ComplexF64(0.0, 0.0)
+    @simd for n in nn:-1:2
+        let n_over_y = n / y
+            @views d[n-1] = n_over_y - (1.0 / (d[n] + n_over_y))
+        end
+    end
+
+    psi0 = cos(x)
+    psi1 = sin(x)
+    qsca::Float64 = 0.0
+    # xi0 = ComplexF64(psi0, psi1)
+    xi1 = ComplexF64(psi1, -psi0)
+    chi1 = psi0
+    chi0 = -psi1
+    s1_1 = zeros(ComplexF64, nang)
+    s1_2 = zeros(ComplexF64, nang)
+    s2_1 = zeros(ComplexF64, nang)
+    s2_2 = zeros(ComplexF64, nang)
+    pi0::Vector{Float64} = zeros(nang)
+    pi1::Vector{Float64} = ones(nang)
+    tau::Vector{Float64} = similar(Vector{Float64}, nang)
+    @inbounds @simd for n in (1:nstop)::UnitRange{Int64}
+        psi, chi = let nm1x2 = 2 * n - 1
+            nm1x2 * psi1 / x - psi0,
+            nm1x2 * chi1 / x - chi0
+        end
+        xi = psi - (chi * im)
+        an::ComplexF64, bn::ComplexF64 = let @views dn = d[n]
+            let an_mult::ComplexF64 = dn / refrel + n / x
+                (an_mult * psi - psi1) / (an_mult * xi - xi1)
+            end,
+            let bn_mult::ComplexF64 = refrel * dn + n / x
+                (bn_mult * psi - psi1) / (bn_mult * xi - xi1)
+            end
+        end
+
+        qsca += (2 * n + 1) * ((abs(an)^2) + (abs(bn)^2))
+        pi_ = pi1
+
+        let np1x2 = 2 * n + 1, tau = n * amu .* pi_ - (n + 1) .* pi0
+            let fn = np1x2 / (n * (n + 1)), anpi = an * pi_, antau = an * tau, bnpi = bn * pi_, bntau = bn * tau
+                s1_1::Vector{ComplexF64} .+= fn * (anpi + bntau)
+                s2_1::Vector{ComplexF64} .+= fn * (antau + bnpi)
+                if isodd(n)
+                    s1_2::Vector{ComplexF64} .+= fn * (anpi - bntau)
+                    s2_2::Vector{ComplexF64} .+= fn * (bnpi - antau)
+                else
+                    s1_2::Vector{ComplexF64} .-= fn * (anpi - bntau)
+                    s2_2::Vector{ComplexF64} .-= fn * (bnpi - antau)
+                end
+            end
+            pi1 = (np1x2 * amu .* pi_ - (n + 1) .* pi0) / n
+        end
+        psi0, chi0 = psi1, chi1
+        psi1, chi1 = psi, chi
+        xi1 = ComplexF64(psi1, -chi1)
+        pi0 = pi_
+    end
+    @inbounds s1::Vector{ComplexF64} = vcat(s1_1, s1_2[1:-1:end-1])
+    @inbounds s2::Vector{ComplexF64} = vcat(s2_1, s2_2[1:-1:end-1])
+    x_sq::Float64 = x^2
+    qsca *= (2.0 / (x_sq))
+    qext::Float64, qback::Float64 = let four_over_x_sq::Float64 = 4.0 / x_sq
+        (four_over_x_sq) * real(s1[1]),
+        (four_over_x_sq) * (abs(s1[2*nang-1])^2)
+    end
+    return (qext, qsca, qback, s1, s2)
+end
+
+"Helper function to calculate the apparent cross section"
+@memoize function _calc_apparent_cross_section(radcore::Float64, refmed::Float64)::Float64
+    TWO_PI * radcore * refmed
+end
+
+# The following is a messy setup to allow for caching of the interpolation objects
+# The reason this has been implemented this way is because Interpolations.jl has
+# very confusing behavior when it comes to typing and caching. This is a workaround.
+struct InterpolationPair
+    rn::Interpolations.Extrapolation{Float64,1,Interpolations.GriddedInterpolation{Float64,1,Vector{Float64},Interpolations.Gridded{Interpolations.Linear{Interpolations.Throw{Interpolations.OnGrid}}},Tuple{Vector{Float64}}},Interpolations.Gridded{Interpolations.Linear{Interpolations.Throw{Interpolations.OnGrid}}},Interpolations.Throw{Nothing}}
+    rk::Interpolations.Extrapolation{Float64,1,Interpolations.GriddedInterpolation{Float64,1,Vector{Float64},Interpolations.Gridded{Interpolations.Linear{Interpolations.Throw{Interpolations.OnGrid}}},Tuple{Vector{Float64}}},Interpolations.Gridded{Interpolations.Linear{Interpolations.Throw{Interpolations.OnGrid}}},Interpolations.Throw{Nothing}}
+end
+
+const _cache = Dict{UInt64,InterpolationPair}()
+
+"Helper function to get the interpolation objects for rn and rk"
+ function _get_rnrk_interp_objects(refcore::Array{Float64,2})::InterpolationPair
+    refcore_hash = hash(refcore)
+    if !haskey(_cache, refcore_hash)
+        _cache[refcore_hash] = InterpolationPair(LinearInterpolation(refcore[:, 1], refcore[:, 2]), LinearInterpolation(refcore[:, 1], refcore[:, 3]))
+    end
+    return _cache[refcore_hash]
+end
+
+"Helper function to calculate the delta wavelength for a given wavelength range and number of values"
+@memoize function _calc_delta_wl(wavel1::Float64, wavel2::Float64, numval::UInt32)::Float64
+    (wavel2 - wavel1) / (numval - 1)
+end
+
+"Helper function to calculate the size parameter"
+@memoize function _calc_size_parameter(apparent_cross_section::Float64, wavel::Float64)::Float64
+    apparent_cross_section / wavel
+end
+
+"Wrapper for partial dispatch of the bhmie function"
+struct _bhmie
+    nang::UInt32
+end
+
+function (p::_bhmie)(x::Float64, refrel::ComplexF64)::Float64
+    first(bhmie(x, refrel, p.nang))::Float64
+end
+
+"Wrapper for partial dispatch of the mfp function"
+struct _mfp
+    fv::Float64
+    radcore::Float64
+    omp::Float64
+    om0::Float64
+end
+
+function(p::_mfp)(wavel::Float64, rn::Float64, rk::Float64)::ComplexF64
+    refre1, refim1, _ = mfp(p.fv, wavel, p.radcore, p.omp, p.om0, rn, rk)
+    return ComplexF64(refre1, refim1)::ComplexF64
+end
+
+"""
+Calculates the expected extinction coefficient spectrum for particles of a given
+size for uncoated particles (modified/translated from Haiss et als FORTRAN implementation)
+"""
+function qbare(wavel1::Float64, wavel2::Float64, numval::UInt32, scangles::UInt32,
+    refmed::Float64, radcore::Float64, omp::Float64, om0::Float64,
+    fv::Float64, refcore::Array{Float64,2})::Matrix{Float64}
+
+    # calculate the new wavelengths
+    qarray::Matrix{Float64} = let wavelengths::StepRangeLen{Float64} = wavel1:_calc_delta_wl(wavel1, wavel2, numval):wavel2
+        hcat(wavelengths, Vector{Float64}(undef, numval))
+    end
+    interp_pair::InterpolationPair = _get_rnrk_interp_objects(refcore)
+    let interp_ref_rn = interp_pair.rn, interp_ref_rk = interp_pair.rk, map_fn = _bhmie(scangles), mfp_fn = _mfp(fv, radcore, omp, om0)
+        @inbounds @threads for i in 0x00000001:numval
+            wl = qarray[i, 1]
+            qarray[i, 2] = map_fn(
+                _calc_apparent_cross_section(radcore, refmed) / wl,
+                mfp_fn(wl, interp_ref_rn(wl), interp_ref_rk(wl)) / refmed
+            )
+        end
+    end
+
+    return qarray
+end
+
+"""
+Calculates the expected extinction coefficient spectrum for particles of a given
+size for coated particles (translated from Haiss et als FORTRAN implementation)
+"""
+function qcoat(wavel1::Float64, wavel2::Float64, numval::Int64, refmed::Float64,
+    radcor::Float64, refre1::Float64, refim1::Float64,
+    radcot::Float64, refre2::Float64, refim2::Float64,
+    omp::Float64, om0::Float64, fv::Float64, refcore::Array{Float64,2})
+    # TODO: Add static output!
+    # refcore is an array containing the spectrum of refractive index vs particle size
+    tau::Float64 = (1 / om0) * 1e-14
+    fmpinf::Float64 = fv * tau
+    delta::Float64 = (wavel2 - wavel1) / (numval - 1)
+
+    wavrnrk::Vector{Tuple{Float64,Float64,Float64}} = Vector{Tuple{Float64,Float64,Float64}}()
+    k::UInt32 = 1
+    @simd for i in 1:numval
+        wavel::Float64 = wavel1 + (i - 1) * delta
+        # find values for refractive index of particles using linear interpolation:
+        if (wavel != refcore[1, 1])
+            while (refcore[k, 1] < wavel)
+                k += 1
+            end
+            push!(wavrnrk, (wavel,
+                (refcore[k-1, 2] + (wavel - refcore[k-1, 1]) / (refcore[k, 1] - refcore[k-1, 1]) * (refcore[k, 2] - refcore[k-1, 2])),
+                (refcore[k-1, 3] + (wavel - refcore[k-1, 1]) / (refcore[k, 1] - refcore[k-1, 1]) * (refcore[k, 3] - refcore[k-1, 3])))
+            )
+        end
+    end
+    qarray::Array{Float64} = Array{Float64}(0, 4)
+    @inbounds @simd for vals in wavrnrk
+        wavel::Float64, rn::Float64, rk::Float64 = vals
+        # adjust for "mean free path effect" on free electrons in metal core
+        refre1, refim1, om0r = mfp(fv, wavel, radcor, omp, om0, rn, rk)
+        # calculate relative size parameters x and y for core and coating respectively
+        two_pi_times_refractive_ratio::Float64 = TWO_PI * (refmed / refre1)
+        x = radcor * two_pi_times_refractive_ratio
+        y = radcot * two_pi_times_refractive_ratio
+        # calculate ComplexF64 refractive indices:
+        rfrel1 = refre1 + (refim1 * im) / refmed
+        rfrel2 = refre2 + (refim2 * im) / refmed
+        # appends data to qarray as a row in the format [wavel qext qsca qback]
+        # bhcoat outputs are adjusted to be relative to core instead of coating
+        qarray = vcat(qarray, [wavel (bhcoat(x, y, rfrel1, rfrel2) .* ((radcot / radcor)^2))...])
+    end
+    return qarray
+end
+
+end

src/nanoconc.jl

@@ -0,0 +1,324 @@
+"""
+A module that uses the functions contained in miemfp and quantumcalc to find the
+concentration of a nanoparticle colloid given enough information
+"""
+@fastmath module nanoconc
+
+export fetchinfo, dispinfo, addmaterial, loadmaterial, listmaterials, delmaterial, editmaterial, qpredict, abspredict
+
+include("miemfp.jl")
+include("quantumcalc.jl")
+
+using CSV
+using DataFrames
+using HDF5
+
+const matfile = "data/materials.h5" # this is the location of the material datafile
+
+"""
+Fetches info string for a saved material
+"""
+function fetchinfo(material::String)::String
+    return string("\nMaterial:\t\t", material,
+        "\n- Valence:\t\t", h5readattr(matfile, material)["z"],
+        "\n- Atomic Mass:\t\t", h5readattr(matfile, material)["am"],
+        "amu\n- Density:\t\t", h5readattr(matfile, material)["rho"],
+        "g/cm^3\n- Resistivity:\t\t", h5readattr(matfile, material)["res"],
+        "g/cm^3\nDerived values:\n- Plasma frequency:\t", h5readattr(matfile, material)["omp"],
+        "e14Hz\n- Collision frequency:\t", h5readattr(matfile, material)["om0"],
+        "e14Hz\n- Fermi Velocity:\t", h5readattr(matfile, material)["fv"], "cm/s\n",
+        h5readattr(matfile, material)["info"],
+        "\nDescription:\n\"", h5readattr(matfile, material)["description"],
+        "\"\n")
+end
+
+"""
+Displays info for saved material
+"""
+function dispinfo(material::String)
+    print(fetchinfo(material))
+end
+
+"""
+stores material data for future use in the materials HDF5 data file with the given
+material name and description. Requires valence (z), atomic mass (am), density (rho)(g/cm^3)
+resistivity (res)(ohm*m) and a complex refractive index spectrum from a file at filepath.
+Refractive index data is taken from an appropriately formatted csv file (formatted with
+columns entitled "w","n" and "k" for wavelength in nm, n and k respectively)
+"""
+function addmaterial(z::Float64, am::Float64, rho::Float64, res::Float64,
+    filepath::String, material::String, description::String;
+    disp::Bool=true)
+    try
+        flag = h5open(matfile, "r") do file
+            has(file, material)
+        end
+    catch
+        flag = false
+    end
+
+    if !flag
+        df::DataFrames.DataFrame = CSV.read(filepath)
+        data::Array{Float64,2} = (convert(Array{Float64,2}, hcat(df[:w], df[:n], df[:k])))
+        h5write(matfile, material, data)
+        omp, om0, fv = quantumcalc.mieparams(z, am, rho, res)
+        println(omp)
+        println(om0)
+        println(fv)
+        h5writeattr(matfile, material, Dict(
+            "z" => z,
+            "am" => am,
+            "rho" => rho,
+            "res" => res,
+            "omp" => omp,
+            "om0" => om0,
+            "fv" => fv,
+            "description" => description,
+            "info" => string("Wavelength Range:\t", minimum(data[:, 1]),
+                "nm-", maximum(data[:, 1]),
+                "nm\nNumber of datapoints:\t",
+                size(data)[1])))
+        if disp == true
+            println("\nAdded material:")
+            dispinfo(material)
+        end
+    elseif disp == true
+        println("\nMaterial ", material, " Already exists!\n\nExisting material:")
+        dispinfo(material)
+    end
+end
+
+"""
+An alternative version of the addmaterial function for materials with known mie
+    parameters. Requires input of plasma freq (omp) in Hz/1e14, collision freq (om0) in Hz/1e14 and
+    Fermi velocity (fv) in cm/s
+"""
+function addmaterial(omp::Float64, om0::Float64, fv::Float64,
+    filepath::String, material::String, description::String;
+    disp::Bool=true)
+    try
+        flag = h5open(matfile, "r") do file
+            has(file, material)
+        end
+    catch
+        flag = false
+    end
+
+    if !flag
+        df::DataFrames.DataFrame = CSV.read(filepath)
+        data::Array{Float64,2} = (convert(Array{Float64,2}, hcat(df[:w], df[:n], df[:k])))
+        h5write(matfile, material, data)
+        println(omp)
+        println(om0)
+        println(fv)
+        h5writeattr(matfile, material, Dict(
+            "omp" => omp,
+            "om0" => om0,
+            "fv" => fv,
+            "description" => description,
+            "info" => string("Wavelength Range:\t", minimum(data[:, 1]),
+                "nm-", maximum(data[:, 1]),
+                "nm\nNumber of datapoints:\t",
+                size(data)[1])))
+        if disp == true
+            println("\nAdded material:")
+            dispinfo(material)
+        end
+    elseif disp == true
+        println("\nMaterial ", material, " Already exists!\n\nExisting material:")
+        dispinfo(material)
+    end
+end
+
+"""
+This function loads the stored miemfp data for a saved material. omp, om0, fv and
+refractive index array are returned as a static tuple that can be readily splatted
+into the last few args of the miemfp.qbare and miemfp.qcoat functions
+"""
+function loadmaterial(material::String; disp::Bool=true
+)::Tuple{Float64,Float64,Float64,Array{Float64,2}}
+    local output::Tuple{Float64,Float64,Float64,Array{Float64,2}}
+    h5open(matfile, "r") do file
+        material_data = file[material]
+        output = (
+            read(material_data["omp"]),
+            read(material_data["om0"]),
+            read(material_data["fv"]),
+            read(material_data)
+        )
+    end
+    if disp == true
+        println("\nLoading material:")
+        dispinfo(material)
+        println()
+    end
+    return output
+end
+
+"""
+This function lists all materials for which refractive index data is stored
+"""
+function listmaterials()
+    h5open(matfile, "r") do file
+        flag::Bool = false
+        for i in file
+            flag = true
+            break
+        end
+        if flag == true
+            println("\nSaved materials:")
+            for material in file
+                println(name(material))
+            end
+        else
+            println("\nNo materials saved!")
+        end
+    end
+end
+
+"""
+This function deletes a material from the materials datafile
+"""
+function delmaterial(material::String; disp::Bool=true)
+    h5open(matfile, "r+") do file
+        if exists(file, material)
+            if disp == true
+                println("\nDeleting material:")
+                dispinfo(material)
+                println()
+            end
+            o_delete(file, material)
+        else
+            println("\nMaterial \"", material, "\" not found!")
+        end
+    end
+end
+
+"""
+This function allows editing of specific parameters for a stored material
+"""
+function editmaterial(material::String, param::String, value; disp::Bool=true)
+    h5open(matfile, "r+") do file
+        if !exists(file, material)
+            if disp
+                println("\nMaterial \"", material, "\" not found!")
+            end
+            return
+        end
+    end
+    matattr = h5readattr(matfile, material)
+    matdata = h5read(matfile, material)
+    if param == "refind"
+        val0 = matdata
+        value::typeof(matdata)
+        matdata = value
+    else
+        val0 = matattr[param]
+        value::typeof(val0)
+        matattr[param] = value
+    end
+    h5open(matfile, "r+") do file
+        o_delete(file, material)
+    end
+    h5write(matfile, material, matdata)
+    h5writeattr(matfile, material, matattr)
+    if disp
+        println(material, " Parameter ", param, " changed from ", val0, " to ", value)
+    end
+end
+
+"""
+A function that predicts the expected average per particle extinction efficiency
+spectrum for a colloid within a wavelength range (wavel1 - wavel2) with a set
+number of points (numval) per spectrum. Also requires refractive index of the
+medium, a particledata array containing the diameters of the particles in solution
+(in nm) and their relative amounts, and the saved material data array
+"""
+function qpredict(refmed::Float64, wavel1::Float64, wavel2::Float64,
+    numval::UInt32, particledata::Array{Float64,2},
+    materialdata::Tuple{Float64,Float64,Float64,Array{Float64,2}}
+)::Array{Float64,2}
+    # normalise relative amounts of particles to a total of 1
+    particledata[:, 2] = particledata[:, 2] ./ sum(particledata[:, 2])
+    # convert particle diameters in particledata to radii
+    particledata[:, 1] = particledata[:, 1] ./ 2
+    # define a wrapper function for later vectorised calculations
+    spec(x) = miemfp.qbare(wavel1, wavel2, numval, 0x00000002, refmed, x, materialdata...)
+    # note: the UInt32 above for "scangles" must be at least 2 for accurate Qext. Higher "scangles"
+    ##	does not increase Qext precision but will be important if modding to also account for Qsca
+    # perform vectorised spectrum prediction on particledata array
+    spectra = spec.(particledata[:, 1])
+    # for each spectrum, multiply all values by their relative amounts from the particledata array
+    @inbounds @simd for i in eachindex(spectra)
+        spectra[i][:, 2] = particledata[i, 2] .* spectra[i][:, 2]
+    end
+    # prep an array for final predicted spectrum
+    data::Array{Float64,2} = Array{Float64,2}(undef, numval, 2)
+    data[:, 1] = spectra[1][:, 1]
+    data[:, 2] = zeros(numval)
+    # for each weighted spectrum, add it to the total spectrum in the "data" array
+    @inbounds @simd for x in spectra
+        data[:, 2] = data[:, 2] .+ x[:, 2]
+    end
+    return data
+end
+
+struct q_to_sigma
+    geometric_cross_section::Vector{Float64}
+
+    function q_to_sigma(diameter::Vector{Float64})::q_to_sigma
+        new(π .* ((diameter ./ 2).^2))
+    end 
+end
+
+function (p::q_to_sigma)(q::Matrix{Float64})::Array{Float64, 3}
+    # Calculate the extinction cross-section
+    hcat([q .* x for x in p.geometric_cross_section])
+end
+
+function sigmapredict(refmed::Float64, wavel1::Float64, wavel2::Float64,
+    numval::UInt32, particledata::Array{Float64,2},
+    materialdata::Tuple{Float64,Float64,Float64,Array{Float64,2}}
+)::Array{Float64,2}
+    # normalise relative amounts of particles to a total of 1
+    particledata[:, 2] = particledata[:, 2] ./ sum(particledata[:, 2])
+    # convert particle diameters in particledata to radii
+    particledata[:, 1] = particledata[:, 1] ./ 2
+    # define a wrapper function for later vectorised calculations
+    spec(x) = miemfp.qbare(wavel1, wavel2, numval, 0x00000002, refmed, x, materialdata...)
+    # note: the UInt32 above for "scangles" must be at least 2 for accurate Qext. Higher "scangles"
+    ##	does not increase Qext precision but will be important if modding to also account for Qsca
+    # perform vectorised spectrum prediction on particledata array
+    spectra = spec.(particledata[:, 1])
+    println(size(spectra))
+    # for each spectrum, multiply all values by their relative amounts from the particledata array
+    @inbounds @simd for i in eachindex(spectra)
+        spectra[i][:, 2] = particledata[i, 2] .* spectra[i][:, 2]
+    end
+    # prep an array for final predicted spectrum
+    data::Array{Float64,2} = Array{Float64,2}(undef, numval, 2)
+    data[:, 1] = spectra[1][:, 1]
+    data[:, 2] = zeros(numval)
+    # for each weighted spectrum, add it to the total spectrum in the "data" array
+    @inbounds @simd for x in spectra
+        data[:, 2] = data[:, 2] .+ x[:, 2]
+    end
+    return data
+end
+
+function abspredict(refmed::Float64, wavel1::Float64, wavel2::Float64,
+    numval::UInt32, particledata::Array{Float64,2},
+    materialdata::Tuple{Float64,Float64,Float64,Array{Float64,2}},
+    ppml::Float64, d0::Float64)::Array{Float64,2}
+    data = qpredict(refmed, wavel1, wavel2, numval, particledata, materialdata)
+    predict(qext) = quantumcalc.predictabs(
+        qext,
+        sum([particledata[i, 1] * particledata[i, 2] for i in 1:size(particledata)[1]]) / sum(particledata[:, 2]),
+        ppml,
+        d0
+    )
+    data[:, 2] = predict.(data[:, 2])
+    return data
+end
+
+end

src/quantumcalc.jl

@@ -0,0 +1,196 @@
+# NOTE: Future functionality to add
+# calculation of plasma frequency/collision frequency for non-metals
+#       can be done based on electron's "effective mass" and charge
+#       carrier density instead of e- density?
+
+"""
+A module containing functions that apply various quantum physics formulae
+for use with "miemfp.jl" in spectrum prediction/concentration calculation
+"""
+@fastmath module quantumcalc
+
+using Memoize
+
+# Planck's constant in JS according to the CODATA 2014 definition
+const h::Float64 = 6.626070040e-34
+# Reduced Planck's Constant
+const hbar::Float64 = h / (2 * pi)
+# Fine structure constant based on the CODATA 2014 definition
+const alpha::Float64 = 7.2973525664e-3
+# Speed of light in vacuum by the CODATA 2014 definition
+const c::Float64 = 299752458
+# Elementary charge in C according to the CODATA 2014 definition
+const ec::Float64 = 1.6021766208e-19
+# Boltzmann Constant in J/K according to CODATA 2014 definition
+const kb::Float64 = 1.38064852e-23
+# Electron mass in kg according to the 2014 CODATA definition of 1amu and ref: DOI:10.1038/nature13026
+const me::Float64 = 9.109383555654034e-31
+# Bohr radius in m
+const a0::Float64 = hbar / (me * c * alpha)
+# Avogadro's constant in mol^-1 according to the CODATA 2014 definition
+const A::Float64 = 6.022140857e23
+# Precomputed constants for optimization of functions
+const A_E6::Float64 = A * 1e6
+const PI_SQUARED::Float64 = pi^2
+const THREE_PI_SQUARED::Float64 = 3 * PI_SQUARED
+const TWO_PI::Float64 = 2 * pi
+const FOUR_PI::Float64 = 2 * TWO_PI
+const TWO_PI_TIMES_11_44E15::Float64 = TWO_PI * 11.44e15
+const EC_SQUARED::Float64 = ec^2
+const A_PI_LOG10_EULER_OVER_1000::Float64 = (A * pi * log10(ℯ)) / 1000
+const NUM_TO_MOL_CONV::Float64 = 1000 / A
+
+
+"""
+Calculates the free electron density (ne) in m^-3 of an elemental material given
+its valence (z), atomic mass in amu (am) and mass density (rho) in g/cm^2
+Source: ISBN 0-03-083993-9 (page 4)
+"""
+function metalfed(z::Float64, rho::Float64, am::Float64)::Float64
+    (A_E6 * ((z * rho) / am))
+end
+
+"""
+Calculates the Fermi wave vector (kf) of electrons in a material given
+electron density (ne)
+Source: ISBN 0-03-083993-9 (page 36)
+"""
+function fermivec(ne::Float64)::Float64
+    cbrt(ne * THREE_PI_SQUARED)
+end
+
+"""
+Calculates the Fermi velocity (fv) in m/s of electrons in a material given Fermi
+wave vector (kf)
+Source: ISBN 0-03-083993-9 (page 36)
+"""
+function fermivel(kf::Float64; mstar::Float64=me)::Float64
+    (hbar / mstar) * kf
+end
+
+"""
+Calculates the electron sphere radius (rs) given its free electron density (ne)
+Source: ISBN 0-03-083993-9 (page 4)
+"""
+function spherad(ne::Float64)::Float64
+    cbrt(3 / (FOUR_PI * ne))
+end
+
+"""
+Calculates the plasma frequency (omp) of a material in Hz give its electron
+sphere radius (rs)
+Source: ISBN 0-03-083993-9 (page 758)
+"""
+function plasmafreq(rs::Float64)::Float64
+    TWO_PI_TIMES_11_44E15 * ((rs / a0)^(-1.5))
+end
+
+"""
+Calculates the collision frequency (om0) in Hz for a material given its free
+electron density (ne) and its resistivity (res)
+Source: ISBN 0-03-083993-9 (page 8)
+"""
+function collfreq(ne::Float64, res::Float64; mstar::Float64=me)::Float64
+    r((EC_SQUARED) * ne * res) / mstar
+end
+
+"""
+Calculates the plasma freq (omp) in Hz/1e14, collision freq (om0) in Hz/1e14 and
+Fermi velocity (fv) in cm/s of an elemental material given its valence (z),
+atomic mass (am) in amu, mass density (rho) in g/cm^2 and its resisitivity (res)
+in Ohm*m as a tuple of (omp,om0,fv) for use in Mie model algorithms contained in
+the "miemfp" module
+"""
+function mieparams(z::Float64, am::Float64, rho::Float64, res::Float64)
+    ne::Float64 = metalfed(z, rho, am)::Tuple{Float64,Float64,Float64}
+    @sync begin
+        @async omp = plasmafreq(spherad(ne)) * 1e-14
+        @async om0 = collfreq(ne, res) * 1e-14
+        @async fv = fermivel(fermivec(ne)) * 1e2
+    end
+
+    (omp, om0, fv)
+end
+
+"""
+Converts Qext to the molar attenuation coefficient (κ) based on the average radius
+of the particles present. Formula sourced from ISSN 1061-933X
+formula: κ = π * R^2 * Qext * Na * log(e) / 1000
+"""
+function attencoeff(qext::Float64, r::Float64)::Float64
+    (r^2) * qext * A_PI_LOG10_EULER_OVER_1000
+end
+
+"""
+Converts units per millilitre to moles per litre
+"""
+@memoize function numtomol(n::Float64)::Float64
+    n * NUM_TO_MOL_CONV
+end
+
+"""
+Finds absorbance based on molar attenutaion coefficient (κ), concentration (molar)
+(c) and path length (in cm) (d0) based on the Beer/Lambert law
+"""
+function attentoabs(kappa::Float64, c::Float64, d0::Float64)::Float64
+    kappa * c * d0
+end
+
+"""
+Predicts the abs value as displayed on a spectrum from the extinction efficiency
+(qext), average particle radius (avgr), particles per ml (ppml) and path length (l)
+"""
+function predictabs(qext::Float64, avgr::Float64, ppml::Float64, l::Float64)::Float64
+    log10(attentoabs(attencoeff(qext, avgr), numtomol(ppml), l))
+end
+
+# Below is an old, failed version of the predictabs code above, kept in case useful in later debugging
+#
+# @doc """
+# This function takes an absorbance (Abs), extinction coefficient (qext), average
+# particle radius in nm (r) and path length (d0) and returns a count of particles
+# per ml
+# """ ->
+# function partperml(Abs::Float64,qext::Float64,r::Float64,d0::Float64)
+#     return (log(10) * Abs) / (pi * (r ^ 2) *qext * d0) ::Float64
+# end
+#
+# @doc """
+# This function takes number of particles per ml, extinction coefficient (qext),
+# average particle radius in nm (r) and path length (d0) and returns a predicted
+# absorbance value (from DOI: 10.1016/j.saa.2017.10.047)
+# """ ->
+# function abspredict(N::Float64,qext::Float64,r::Float64,d0::Float64)
+#     # note: log10 here is not part of formula but it is required to convert to AU (i think)
+#     return log10((pi * (r^2) *qext * d0 * N) / log(10)) ::Float64
+# end
+
+############## FUNCTIONS BELOW THIS POINT ARE PROTOTYPES #######################
+
+# NOTE: if successful change mieparams to miemetal
+function XXmieionic(z::Float64, rho::Float64, mm::Float64, shc::Float64, res::Float64)::Tuple{Float64,Float64,Float64}
+    # mm is molar mass in g/mol
+    ne::Float64 = metalfed(z, rho, mm)
+    mstar::Float64 = (shc * hbar^2 * (3 * PI_SQUARED * ne)^(2 / 3)) / (PI_SQUARED * kb^2 * A * z)
+
+    return (
+        plasmafreq(spherad(ne)) * 1e-14,
+        collfreq(ne, res, mstar=mstar) * 1e-14,
+        fermivel(fermivec(ne), mstar=mstar) * 1e2
+    )
+end
+
+function mieionic(z::Float64, rho::Float64, mm::Float64, shc::Float64, res::Float64)::Tuple{Float64,Float64,Float64}
+    # mm is molar mass in g/mol
+    ne::Float64 = metalfed(z, rho, mm)
+    rs::Float64 = spherad(ne)
+    mstar::Float64 = shc * me / 0.169 * z * ((rs / a0)^2) * 1e-4
+
+    return (
+        plasmafreq(rs) * 1e-14,
+        collfreq(ne, res, mstar=mstar) * 1e-14,
+        fermivel(fermivec(ne), mstar=mstar) * 1e2
+    )
+end
+
+end

test/miemfp_tests.jl

@@ -0,0 +1,31 @@
+using Test
+
+if !@isdefined TestUtils
+    include("testutils.jl")
+end
+if !@isdefined nanoconc
+    include("../src/nanoconc.jl")
+end
+
+@testset "miemfp" begin
+    @testset "miemfp.bhmie" begin
+        TestUtils.test_from_serialized(
+            nanoconc.miemfp.bhmie,
+            "test/data/Main.nanoconc.miemfp.bhmie.ser"
+        )
+    end
+
+    @testset "miemfp.mfp" begin
+        TestUtils.test_from_serialized(
+            nanoconc.miemfp.mfp,
+            "test/data/Main.nanoconc.miemfp.mfp.ser"
+        )
+    end
+
+    @testset "miemfp.qbare" begin
+        TestUtils.test_from_serialized(
+            nanoconc.miemfp.qbare,
+            "test/data/Main.nanoconc.miemfp.qbare.ser"
+        )
+    end
+end

test/nanoconc_tests.jl

@@ -0,0 +1,24 @@
+using Test
+
+if !@isdefined TestUtils
+    include("testutils.jl")
+end
+if !@isdefined nanoconc
+    include("../src/nanoconc.jl")
+end
+
+@testset "nanoconc" begin
+    @testset "qpredict" begin
+        TestUtils.test_from_serialized(
+            nanoconc.qpredict,
+            "test/data/Main.nanoconc.qpredict.ser"
+        )
+    end
+
+    @testset "nanoconc" begin
+        TestUtils.test_from_serialized(
+            nanoconc.abspredict,
+            "test/data/Main.nanoconc.abspredict.ser"
+        )
+    end
+end

test/quantumcalc_tests.jl

@@ -0,0 +1,38 @@
+using Test
+
+if !@isdefined TestUtils
+    include("testutils.jl")
+end
+if !@isdefined nanoconc
+    include("../src/nanoconc.jl")
+end
+
+@testset "quantumcalc" begin
+    @testset "quantumcalc.attencoeff" begin
+        TestUtils.test_from_serialized(
+            nanoconc.quantumcalc.attencoeff,
+            "test/data/Main.nanoconc.quantumcalc.attencoeff.ser"
+        )
+    end
+
+    @testset "quantumcalc.attentoabs" begin
+        TestUtils.test_from_serialized(
+            nanoconc.quantumcalc.attentoabs,
+            "test/data/Main.nanoconc.quantumcalc.attentoabs.ser"
+        )
+    end
+
+    @testset "quantumcalc.numtomol" begin
+        TestUtils.test_from_serialized(
+            nanoconc.quantumcalc.numtomol,
+            "test/data/Main.nanoconc.quantumcalc.numtomol.ser"
+        )
+    end
+
+    @testset "quantumcalc.predictabs" begin
+        TestUtils.test_from_serialized(
+            nanoconc.quantumcalc.predictabs,
+            "test/data/Main.nanoconc.quantumcalc.predictabs.ser"
+        )
+    end
+end

test/runtests.jl

@@ -0,0 +1,8 @@
+using Test
+
+include("testutils.jl")
+include("../src/nanoconc.jl")
+
+include("nanoconc_tests.jl")
+include("miemfp_tests.jl")
+include("quantumcalc_tests.jl")

test/testutils.jl

@@ -0,0 +1,25 @@
+module TestUtils
+using Serialization
+
+using Test
+
+function fieldvalues(obj)
+    [getfield(obj, f) for f in fieldnames(typeof(obj))]
+end
+
+deep_compare(a, b; rtol::Real=sqrt(eps()), atol::Real=0.) = a == b # base case: use ==
+deep_compare(a::Union{AbstractFloat, Complex}, b::Union{AbstractFloat, Complex}; rtol::Real=sqrt(eps()), atol::Real=0.) = isapprox(a, b) # for floats, use isapprox
+deep_compare(a::Union{Array{AbstractFloat}, Array{Complex}}, b::Union{Array{AbstractFloat}, Array{Complex}}; rtol::Real=sqrt(eps()), atol::Real=0.) = isapprox(a, b) # for arrays of floats, use isapprox element-wise
+deep_compare(a::AbstractArray, b::AbstractArray; rtol::Real=sqrt(eps()), atol::Real=0.) = all(deep_compare.(a, b; rtol, atol)) # for arrays of other types, recurse
+deep_compare(a::Tuple, b::Tuple; rtol::Real=sqrt(eps()), atol::Real=0.) = all(deep_compare.(a, b; rtol, atol)) # for tuples, recurse
+deep_compare(a::T, b::T; rtol::Real=sqrt(eps()), atol::Real=0.) where {T <: Any} = deep_compare(fieldvalues(a), fieldvalues(b); rtol, atol) # for composite types, recurse
+
+function test_from_serialized(fn::Function, filename::String)
+    argskwargs, out = open(filename, "r") do f
+        deserialize(f)
+    end
+
+    @test deep_compare([fn(a...; kw...) for (a, kw) in argskwargs], out)
+end
+
+end # module TestUtils

val_data/uv-vis/.~lock.UV VIS_Rasodyne _1.xlsx#

@@ -0,0 +1 @@
+Cian Hughes,cianh,cian-worklaptop,04.07.2023 11:49,file:///home/cianh/.config/libreoffice/4;

val_data/uv-vis/Ag_17092021/.~lock.AG1.xlsx#

@@ -0,0 +1 @@
+Cian Hughes,cianh,cian-worklaptop,04.07.2023 14:52,file:///home/cianh/.config/libreoffice/4;