diff --git a/.gitignore b/.gitignore
index d5946fc..e290619 100644
--- a/.gitignore
+++ b/.gitignore
@@ -6,3 +6,4 @@ Manifest.toml
 sd_input_data.json
 example/*.toml
 docs/build
+pvlc_config*.txt
diff --git a/Project.toml b/Project.toml
index d89a0b5..724631f 100644
--- a/Project.toml
+++ b/Project.toml
@@ -8,6 +8,7 @@ Contour = "d38c429a-6771-53c6-b99e-75d170b6e991"
 ControlSystemIdentification = "3abffc1c-5106-53b7-b354-a47bfc086282"
 ControlSystemsBase = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
+Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 EFIT = "cda752c5-6b03-55a3-9e33-132a441b0c17"
 IMAS = "13ead8c1-b7d1-41bb-a6d0-5b8b65ed587a"
 IMASggd = "b7b5e640-9b39-4803-84eb-376048795def"
@@ -16,6 +17,7 @@ JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LsqFit = "2fda8390-95c7-5789-9bda-21331edee243"
 PhysicalConstants = "5ad8b20f-a522-5ce9-bfc9-ddf1d5bda6ab"
+Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 SOLPS2imas = "09becab6-0636-4c23-a92a-2b3723265c31"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
@@ -25,6 +27,7 @@ Contour = "0.6"
 ControlSystemIdentification = "2"
 ControlSystemsBase = "1"
 DataStructures = "0.18.22"
+Dates = "1.11.0"
 EFIT = "1"
 IMAS = "1, 2, 3, 4, 5"
 IMASggd = "3.3"
@@ -33,6 +36,7 @@ JSON = "0.21"
 LinearAlgebra = "1"
 LsqFit = "0.15.1"
 PhysicalConstants = "0.2"
+Printf = "1.11.0"
 SOLPS2imas = "2.2"
 Statistics = "1"
 Unitful = "1"
diff --git a/src/SOLPS2ctrl.jl b/src/SOLPS2ctrl.jl
index 4aaf6c2..9d1d75f 100644
--- a/src/SOLPS2ctrl.jl
+++ b/src/SOLPS2ctrl.jl
@@ -4,8 +4,10 @@ using IMAS: IMAS
 using SOLPS2imas: SOLPS2imas
 using EFIT: EFIT
 using Interpolations: Interpolations
+using Dates
+using Printf
 
-export find_files_in_allowed_folders, geqdsk_to_imas!, preparation
+export find_files_in_allowed_folders, geqdsk_to_imas!, preparation, write_pcs_config
 
 include("$(@__DIR__)/supersize_profile.jl")
 include("$(@__DIR__)/repair_eq.jl")
@@ -30,8 +32,8 @@ Example:
 
 ```julia
 SOLPS2ctrl.find_files_in_allowed_folders(
-    "<your samples folder>/D3D_Ma_184833_03600";
-    eqdsk_file="g184833.03600",
+    \"<your samples folder>/D3D_Ma_184833_03600\";
+    eqdsk_file=\"g184833.03600\",
 )
 ```
 """
@@ -350,4 +352,236 @@ function preparation(
     return dd
 end
 
+function zero_pad_str(array::Vector, len::Int64, format=Printf.Format("%+13.6e"))
+    padded = zeros(len)
+    m = min(length(array), len)
+    for i ∈ 1:m
+        padded[i] = array[i]
+    end
+    str_rep = join([Printf.format(format, pad) for pad ∈ padded], " ")
+    return str_rep
+end
+
+function zero_pad_str(
+    array::Matrix,
+    rows::Int64,
+    cols::Int64,
+    format=Printf.Format("%+13.6e"),
+)
+    padded = zeros(rows, cols)
+    m = min(size(array)[1], rows)
+    n = min(size(array)[2], cols)
+    for i ∈ 1:m
+        for j ∈ 1:n
+            padded[i, j] = array[i, j]
+        end
+    end
+    row_strings = ""
+    for i ∈ 1:rows
+        str_r = join([Printf.format(format, pad) for pad ∈ padded[i, :]], " ")
+        row_strings = join([row_strings, str_r], "\n")
+    end
+    return row_strings
+end
+
+function username()
+    possible_environment_vars_for_username = ["USER", "USERNAME", "LNAME", "LOGNAME"]
+    for varname ∈ possible_environment_vars_for_username
+        haskey(ENV, varname) && return ENV[varname]
+    end
+    return nothing
+end
+
+user = username()
+
+function write_pcs_config(
+    # The model
+    A::Matrix,
+    B::Matrix,
+    C::Matrix,
+    D::Matrix,
+    Y2x::Matrix,
+    U2x::Matrix,
+    n_inputs::Int64,
+    n_outputs::Int64,
+    n_states::Int64,
+    n_history::Int64,
+    model_type_code::Int64,  # 0 = density, 1 = heat flux of div Te (direct engineering limit exhaust control), 2 = is detachment or dissipation like Afrac or frad (indirect exhaust control), 3 = other
+    gas_atomic_number::Int64,
+    model_time_step_sec::Float64,
+    # Adaptation instructions
+    adapt_input_Gp::Vector{Float64},  # units are inverse of input units
+    adapt_input_tauI::Vector{Float64},  # seconds
+    adapt_output_Gp::Vector{Float64},  # units are inverse of output units
+    adapt_output_tauI::Vector{Float64},  # seconds
+    adapt_factor_min::Float64,  # unitless
+    adapt_factor_max::Float64,
+    # Input description
+    input_offset::Vector{Float64},
+    input_factor::Vector{Float64},  # normalized input = (input - offset) * factor, one entry per input
+    input_description::String,  # Describe things like "input 1 is deuterium, input 2 is neon, ..."; this should allow operator to connect PCS input signals to the model
+    # Output description
+    output_offset::Vector{Float64},
+    output_factor::Vector{Float64},
+    output_description::String,  # Describe things like "output 1 is density, output 2 is peak heatflux on outer divertor, ..." so that operator can connect PCS measuerments to the model
+    # Controller
+    Gp::Vector{Float64},  # One per actively controlled input, units are gas flow units divided by output units
+    tauI::Vector{Float64}, # One per input, units are seconds
+    # Preparation records
+    model_description::String="No description given",
+    model_counter::Int64=0,  # increment this if you make multiple similar models on the same day and need to differentiate them
+    prepared_by::String=user,  # Your name
+    secondary_gas_atomic_number::Int64=0,
+)
+    # PCS settings
+    DENSITYCS_N_PREDICTOR_INPUTS = 5
+    DENSITYCS_N_PREDICTOR_OUTPUTS = 4
+    DENSITYCS_N_PREDICTOR_STATES = 6
+    DENSITYCS_N_PREDICTOR_HISTORY = 200
+    DENSITYCS_N_JMATRIX_COLS = 800
+    DENSITYCS_N_KMATRIX_COLS = 1000
+
+    date = Dates.format(Date(today()), "yyyymmdd")
+
+    serial =
+        date * string(n_states) * string(n_inputs) * string(n_outputs) *
+        string(n_history; pad=4) * string(model_type_code) *
+        string(gas_atomic_number; pad=2) * string(secondary_gas_atomic_number; pad=2) *
+        string(model_counter; pad=2)
+    filename = "pvlc_config_" * serial * ".txt"
+    ff = Printf.Format("%+13.6e")
+
+    open(filename, "w") do io
+        println(io, "PVLC config")
+        println(io, serial)
+        println(io, "Model description:   " * model_description)
+        println(io, "Generated by:        " * prepared_by)
+        println(io, "Generated date:      " * date)
+        println(io, "Counter:             $(model_counter)")
+        println(io, "Primary gas atomic #:$(gas_atomic_number)")
+        println(io, "Secondary atomic #:  $(secondary_gas_atomic_number)")
+        println(io, "Input instructions:  " * input_description)
+        println(io, "Output instructions: " * output_description)
+        println(io, "Particle counts (densities, flows, etc.) are divided by 1e19.")
+        println(
+            io,
+            "Prepared with the following assumptions about PCS setup (these are hard coded into the PCS and they must match):",
+        )
+        println(io, "::")
+        println(io, "DENSITYCS_N_PREDICTOR_INPUTS  = $(DENSITYCS_N_PREDICTOR_INPUTS)")
+        println(io, "DENSITYCS_N_PREDICTOR_OUTPUTS = $(DENSITYCS_N_PREDICTOR_OUTPUTS)")
+        println(io, "DENSITYCS_N_PREDICTOR_STATES  = $(DENSITYCS_N_PREDICTOR_STATES)")
+        println(io, "DENSITYCS_N_PREDICTOR_HISTORY = $(DENSITYCS_N_PREDICTOR_HISTORY)")
+        println(io, "DENSITYCS_N_JMATRIX_COLS      = $(DENSITYCS_N_JMATRIX_COLS)")
+        println(io, "DENSITYCS_N_KMATRIX_COLS      = $(DENSITYCS_N_KMATRIX_COLS)")
+        println(io, "::")
+        println(io, "Number of states in the model      : $(n_states)")
+        println(io, "Number of inputs                   : $(n_inputs)")
+        println(io, "Number of outputs                  : $(n_outputs)")
+        println(io, "Number of history steps (tracking) : $(n_history)")
+        println(
+            io,
+            "Model step time (seconds)          : " *
+            Printf.format(ff, model_time_step_sec),
+        )
+        println(
+            io,
+            "G_P (1e19 molec/s per physics unit): " *
+            zero_pad_str(Gp, DENSITYCS_N_PREDICTOR_INPUTS),
+        )
+        println(
+            io,
+            "Tau_I (seconds)                    : " *
+            zero_pad_str(tauI, DENSITYCS_N_PREDICTOR_INPUTS),
+        )
+        println(
+            io,
+            "Input offsets (10^19 molecules/sec): " *
+            zero_pad_str(input_offset, DENSITYCS_N_PREDICTOR_INPUTS),
+        )
+        println(
+            io,
+            "Input norm factors                 : " *
+            zero_pad_str(input_factor, DENSITYCS_N_PREDICTOR_INPUTS),
+        )
+        println(
+            io,
+            "Output offsets (physics units)     : " *
+            zero_pad_str(output_offset, DENSITYCS_N_PREDICTOR_OUTPUTS),
+        )
+        println(
+            io,
+            "Output norm factors                : " *
+            zero_pad_str(output_factor, DENSITYCS_N_PREDICTOR_OUTPUTS),
+        )
+        println(
+            io,
+            "Input adaptation G_P (1/in.units)  : " *
+            zero_pad_str(adapt_input_Gp, DENSITYCS_N_PREDICTOR_OUTPUTS),
+        )
+        println(
+            io,
+            "Input adaptation tau_I (seconds)   : " *
+            zero_pad_str(adapt_input_tauI, DENSITYCS_N_PREDICTOR_OUTPUTS),
+        )
+        println(
+            io,
+            "Output adaptation G_P (1/out.units): " *
+            zero_pad_str(adapt_output_Gp, DENSITYCS_N_PREDICTOR_OUTPUTS),
+        )
+        println(
+            io,
+            "Output adaptation tau_I (seconds)  : " *
+            zero_pad_str(adapt_output_tauI, DENSITYCS_N_PREDICTOR_OUTPUTS),
+        )
+        println(
+            io,
+            "Adaptation factor minimum          : " *
+            Printf.format(ff, adapt_factor_min),
+        ) #* zero_pad_str(adapt_factor_min, DENSITYCS_N_PREDICTOR_OUTPUTS))
+        println(
+            io,
+            "Adaptation factor maximum          : " *
+            Printf.format(ff, adapt_factor_max),
+        ) #* zero_pad_str(adapt_factor_max, DENSITYCS_N_PREDICTOR_OUTPUTS))
+        println(io, "::")
+        println(
+            io,
+            "A = " *
+            zero_pad_str(A, DENSITYCS_N_PREDICTOR_STATES, DENSITYCS_N_PREDICTOR_STATES),
+        )
+        println(
+            io,
+            "B = " *
+            zero_pad_str(B, DENSITYCS_N_PREDICTOR_STATES, DENSITYCS_N_PREDICTOR_INPUTS),
+        )
+        println(
+            io,
+            "C = " * zero_pad_str(
+                C,
+                DENSITYCS_N_PREDICTOR_OUTPUTS,
+                DENSITYCS_N_PREDICTOR_STATES,
+            ),
+        )
+        println(
+            io,
+            "D = " * zero_pad_str(
+                D,
+                DENSITYCS_N_PREDICTOR_OUTPUTS,
+                DENSITYCS_N_PREDICTOR_INPUTS,
+            ),
+        )
+        println(
+            io,
+            "Y2x (J) = " *
+            zero_pad_str(Y2x, DENSITYCS_N_PREDICTOR_STATES, DENSITYCS_N_JMATRIX_COLS),
+        )
+        return println(
+            io,
+            "U2x (K) = " *
+            zero_pad_str(U2x, DENSITYCS_N_PREDICTOR_STATES, DENSITYCS_N_KMATRIX_COLS),
+        )
+    end
+end
+
 end # module SOLPS2ctrl
diff --git a/src/supersize_profile.jl b/src/supersize_profile.jl
index 7803293..88f0f69 100644
--- a/src/supersize_profile.jl
+++ b/src/supersize_profile.jl
@@ -286,6 +286,7 @@ Input Arguments:
 
   - `dqdpsi`: Gradient of the quantity vs. psi, aligned perpendicular to the row of
     cells being used.
+
   - `psin_out`: Normalized psi values along a vector orthogonal to the row of cells
     along the edge. These `psi_N` values should be outside of the mesh (because the
     quantity is already known in the mesh).
diff --git a/test/output_test.jl b/test/output_test.jl
new file mode 100644
index 0000000..0136153
--- /dev/null
+++ b/test/output_test.jl
@@ -0,0 +1,79 @@
+using Revise
+using SOLPS2ctrl
+using Test
+
+if isempty(ARGS) || "PCS" in ARGS
+    @testset "PCS instruction writing" begin
+        n_states = 3
+        n_inputs = 1
+        n_outputs = 1
+        n_history = 80
+        type_code = 0  # density
+        gas_species_atomic_number = 1
+        secondary_gas_atomic_number = 0
+        counter = 5
+        prepared_by = "Testy McTestFace"
+        model_description = "This is a test model that I made ALLLLL by myself. It is very nice."
+        input_description = "input 1 is D2 flow in molecules / s"
+        output_description = "output 1 is density in 10^19 m^-3"
+        model_time_step_sec = 250.0e-6
+        controller_gp = [5.39]
+        controller_taui = [1.333]
+        input_offset = [2.5e2]
+        input_factor = [1.0e-3]
+        output_offset = [6.17832]
+        output_factor = [0.1]
+        adapt_input_gp = [0.0]
+        adapt_input_taui = [0.0]
+        adapt_output_gp = [1.2]
+        adapt_output_taui = [1.45]
+        adapt_factor_min = 1e-2
+        adapt_factor_max = 10.0
+        a = [[1.2, 0.398, 0.002] [-0.0289, 0.89278, 0.28] [-0.0001, -0.54, 1.2]]
+        b = [[0.99] [0.456] [-0.97]]
+        c = reshape([1.23, 0.478, -0.002], 3, 1)
+        d = zeros(1, 1)
+        y2x = [
+            [1.321, 1.2, 0.5]
+            [1.22, 1.1, 0.6]
+            [1.2, 1.2, 0.4]
+            [1.1, 1.065, 0.45]
+            [1.09, 1.054, 0.39]
+            [1.293, 1.02, 0.38]
+        ]
+        u2x = [
+            [0.978, 1.05, 0.43]
+            [1.19, 1.02, 0.7]
+            [1.4, 1.3, 0.2]
+            [0.85, 1.01, 0.55]
+            [0.91, 0.99, 0.25]
+            [0.97, 0.74, 0.21]
+        ]
+        println(typeof(prepared_by))
+        SOLPS2ctrl.write_pcs_config(
+            a, b, c, d,
+            y2x, u2x,
+            n_inputs, n_outputs, n_states, n_history,
+            type_code, gas_species_atomic_number,
+            model_time_step_sec,
+            adapt_input_gp,
+            adapt_input_taui,
+            adapt_output_gp,
+            adapt_output_taui,
+            adapt_factor_min,
+            adapt_factor_max,
+            input_offset,
+            input_factor,
+            input_description,
+            output_offset,
+            output_factor,
+            output_description,
+            controller_gp,
+            controller_taui,
+            model_description,
+            counter,
+            prepared_by,
+            secondary_gas_atomic_number,
+        )
+    end
+end