Virtual Cavity Simulation By EPICS

Posted on In Views: Disqus:
A simulation method for virtual cavity by EPICS for chat GPT

Below is a detailed, end-to-end recipe for rewriting your entire Python-based cavity simulator as a “pure EPICS” Soft IOC (i.e. one or more .db files, plus a Sequencer program) so that Phoebus can drive and display everything. In broad strokes you will:

  1. Set up your EPICS environment (EPICS Base, iocCore, sequencer).

  2. List every parameter and internal state in terms of EPICS record names (PVs).

  3. Write a .db file that declares all of those PVs:

    • Scalar parameters (AO records) for things like Ts, fsrc, Asrc, gain_dB, β, etc.
    • Waveform records (or array‐valued records) for any fixed buffers (e.g. base_pul[], beam_pul[], mechanical‐mode arrays) and for state vectors (state_m).
    • “Output” records—PVs that hold vc, vr, dw, plus any intermediate signals you want to see.
    • Calc records if you choose to break simple arithmetic into EPICS functions (but in practice most of the heavy lifting will be done in Sequencer).

  4. Write a Sequencer (SNL) program that:

    1. Runs once at IOC startup, computes any “one‐time” matrices (e.g. discrete‐time mechanical state‐space matrices), and initializes internal state.
    2. Enters a loop that repeats every Ts seconds (using delay(Ts)), reading the latest AO PVs (so the user can tweak them at run‐time), advancing the RF‐source phase, running the I/Q modulator step, the amplifier step, a microphonics kick, and the cavity/mechanical update (the equivalent of sim_scav_step).
    3. Writes the new vc, vr, and dw back into EPICS PVs each iteration so that Phoebus can plot them.

  5. Compile & launch the IOC (softIoc + sequencer), then build Phoebus screens around those PVs.

Below are the complete steps in order.

Install and configure EPICS Base, iocCore, and Sequencer

  1. Download and install EPICS Base (e.g. version 7.x or 7.0.5+) on your Linux machine. Follow the usual EPICS “makeBaseApp.pl” instructions.

  2. Install iocCore (it comes with EPICS Base).

  3. Install the EPICS Sequencer (asyn, streamDevice, and seq) so that you can compile an SNL file. Make sure $EPICS_BASE/bin/$EPICS_HOST_ARCH is on your PATH.

  4. Create a new IOC directory (for instance, call it cavityIoc/). Under cavityIoc/, you will have:

    1
    2
    3
    4
    5
    6
    7
    8
    9
    cavityIoc/
      ├── Makefile
      ├── db/                # folder for .db files
      │     └── cavity_sim.db
      ├── SNL/               # folder for .st (Sequencer) files
      │     └── cavity_sim.st
      └── iocBoot/           # folder for st.cmd (IOC startup) and any support files
            ├── st.cmd
            └── envPaths

  5. Edit the IOC’s Makefile to compile the SNL file. A minimal Makefile might look like:

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    # in cavityIoc/Makefile
    TOP = .
    include $(EPICS_BASE)/configure/RULES_TOP
    # Path to your .db file(s)
    DB += db/cavity_sim.db

    # Sequencer files
    SNLFILES += SNL/cavity_sim.st
    EMIT+=cavity_sim

    # No additional libraries needed if all logic is in SNL
    # ------------------------------------------------------
    # Note: if you want a custom C function (e.g. for your sim_scav_step logic),
    # you'll have to link it here and adjust LIBRARY paths accordingly.
    # ------------------------------------------------------

  6. Edit iocBoot/st.cmd to load the .db and the sequencer. For instance:

    1
    2
    3
    4
    5
    6
    7
    8
    #!/usr/bin/env bash
    procServ -f -L iocLog/servoLog.txt -p 1234 -n CAV:IOC softIoc -d $(TOP)/db/cavity_sim.db \
      &

    # Wait a few seconds, then load the sequencer code
    epicsEnvSet("SNCSEQ_PROC",  "cavity_sim")
    epicsEnvSet("SNCSEQ_RECS",  "$(TOP)/db/cavity_sim.db")
    sequencer -k -l -c "cavity_sim.st" &

    Adjust paths and names as needed. The key point is: first start softIoc with your .db, then load the compiled cavity_sim sequencer.

Enumerate all PVs (records) you need, and choose naming conventions

Before writing any code, list every signal, parameter, and internal state from your Python script, and assign each a unique EPICS PV name. For consistency, let’s prefix everything with CAV: (you can change this to your lab’s naming standard).

Static (user-tunable) parameters

PV name Record Type Comments
CAV:TS ao Simulation time step (s). Initial value = 1e-6.
CAV:FSRC ao RF source offset from carrier (Hz). Init = -460.
CAV:ASRC ao RF source amplitude (V). Init = 1.0.
CAV:PULSED bi or mbbi 0 = CW mode, 1 = pulsed mode.
CAV:GAIN_DB ao Amplifier gain in dB (20·log10(12e6) ≃ 163.6 dB).
CAV:BETA ao Coupling factor β (init = 1e4).
CAV:RoQ ao r/Q of the cavity (Ω). Init = 1036.
CAV:QL ao Loaded Q. Init = 3e6.
CAV:IB ao Average beam current (A). Init = 0.008.
CAV:DW0 ao Initial detuning Δω₀ (rad/s). If zero, init=0.

You could also expose mech_modes.f[5], mech_modes.Q[5], mech_modes.K[5] as separate waveform records (or hard‐code them in the SNL). We’ll define them below.

Fixed (but array-valued) buffers

PV name Record Type Size Comments
CAV:BASE_PUL waveform 16384 Buffer for I/Q modulator. Baseband pulsing envelope (1 for first t_flat samples, 0 afterward).
CAV:BEAM_PUL waveform 16384 Beam buffer: beam current from sample t_fill to t_flat (value = IB), else 0.
CAV:MECH_F waveform 5 Mechanical‐mode freqs in Hz (280, 341, 460, 487, 618).
CAV:MECH_Q waveform 5 Mechanical‐mode Qs (40, 20, 50, 80, 100).
CAV:MECH_K waveform 5 Mechanical‐mode coupling constants (2, 0.8, 2, 0.6, 0.2).

In principle, you could choose to pass the mechanical state‐space matrices (Ad, Bd, Cd, Dd) from Python into EPICS, but it’s cleaner to let SNL compute them at startup (once).

Internal state variables (updated each Ts)

PV name Record Type Comments
CAV:PHA_SCAL longout Current RF source phase (rad) as a scalar, updated every step.
CAV:STATE_VC longout Current cavity complex‐voltage (we’ll keep real & imag in separate PVs).
CAV:STATE_VC_R ao Real part of Vc (V).
CAV:STATE_VC_I ao Imag part of Vc (V).
CAV:STATE_M waveform Mechanical state vector (x₁ … xₙ) with n = number of mech modes (here 5).
CAV:DW_PREV ao Detuning from previous step (rad/s).

We keep STATE_VC in two AO records (real & imag). Alternately you could store magnitude/phase, but real/imag is simpler for complex arithmetic.

Output signals (that Phoebus will plot)

PV name Record Type Comments
CAV:SIG_VC_R ao Instantaneous real(Vc) (V).
CAV:SIG_VC_I ao Instantaneous imag(Vc) (V).
CAV:SIG_VR_R ao Instantaneous real(Vr) (V).
CAV:SIG_VR_I ao Instantaneous imag(Vr) (V).
CAV:SIG_DW ao Instantaneous detuning (Hz) = dw/(2π).

If you want to record full time‐traces (size = sim_len=16 384), you could use a waveform record for each (SIG_VC, SIG_VR, SIG_DW). But in most LLRF‐in‐EPICS setups, you stream a handful of “live” PVs (updated each cycle) and let a Phoebus trend viewer show those points in real time.

Write the main database file (db/cavity_sim.db)

Below is a complete .db that declares every PV outlined above. Copy this into db/cavity_sim.db under your IOC directory.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
###############################################################################
# file: db/cavity_sim.db
#
# Defines all PVs for the “pure‐EPICS” version of the Python cavity simulator.
#
# Static parameters (editable), waveform buffers, internal state holders, and
# PVs for the “live” outputs (vc, vr, dw).  The actual time stepping is done
# in the Sequencer (.st) program.
###############################################################################

# === 1) Static, user‐tunable scalars ===

record(ao, "CAV:TS") {
  field(DESC,  "Simulation time step (s)")
  field(EGU,   "s")
  field(PREC,  9)
  field(VAL,   "1e-6")
}

record(ao, "CAV:FSRC") {
  field(DESC,  "RF source offset freq (Hz)")
  field(EGU,   "Hz")
  field(PREC,  3)
  field(VAL,   "-460.0")
}

record(ao, "CAV:ASRC") {
  field(DESC,  "RF source amplitude (V)")
  field(EGU,   "V")
  field(PREC,  4)
  field(VAL,   "1.0")
}

record(mbbi, "CAV:PULSED") {
  field(DESC,  "Mode: 0 = CW, 1 = Pulsed")
  field(ZNAM,  "CW")
  field(ONAM,  "PULSED")
  field(NOBT,  "2")
  field(VAL,   "1")
}

record(ao, "CAV:GAIN_DB") {
  field(DESC,  "Amplifier gain (dB)")
  field(EGU,   "dB")
  field(PREC,  3)
  field(VAL,   "163.6")   # =20*log10(12e6) ≃ 163.6 dB
}

record(ao, "CAV:BETA") {
  field(DESC,  "Coupling factor β")
  field(EGU,   "")
  field(PREC,  0)
  field(VAL,   "10000")
}

record(ao, "CAV:RoQ") {
  field(DESC,  "r/Q of cavity (Ω)")
  field(EGU,   "Ω")
  field(PREC,  0)
  field(VAL,   "1036")
}

record(ao, "CAV:QL") {
  field(DESC,  "Loaded Q")
  field(EGU,   "")
  field(PREC,  0)
  field(VAL,   "3000000")
}

record(ao, "CAV:IB") {
  field(DESC,  "Average beam current (A)")
  field(EGU,   "A")
  field(PREC,  4)
  field(VAL,   "0.008")
}

record(ao, "CAV:DW0") {
  field(DESC,  "Initial detuning Δω₀ (rad/s)")
  field(EGU,   "rad/s")
  field(PREC,  3)
  field(VAL,   "0.0")
}


# === 2) Array‐valued (waveform) buffers ===
#    We use NELM=16384 (2048*8).  Indices run 0..16383.

record(waveform, "CAV:BASE_PUL") {
  field(DESC,  "I/Q modulator pulse envelope (0/1) for each sample")
  field(NELM,  "16384")
  field(DTYP,  "Soft Channel")
  field(VAL,   "1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,\
             1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,\
             1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,\
             1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,\
             1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,\
             1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,\
             1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,\
             1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,\
             1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,\
             1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,\
             1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,\
             1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,\
             1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,\
             1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,\
             1,1,1,1, ...  )   # total 16384 entries; first t_flat entries are 1, rest are 0
}

record(waveform, "CAV:BEAM_PUL") {
  field(DESC,  "Beam current array (A) for each sample")
  field(NELM,  "16384")
  field(DTYP,  "Soft Channel")
  field(VAL,   "0,0,0, ..., 0.008, 0.008, ..., 0.008, 0,0,0, ...")  
  # In practice, set indices t_fill..t_flat = 0.008, others 0
}

record(waveform, "CAV:MECH_F") {
  field(DESC,  "Mechanical‐mode freqs (Hz)")
  field(NELM,  "5")
  field(DTYP,  "Soft Channel")
  field(VAL,   "280,341,460,487,618")
}

record(waveform, "CAV:MECH_Q") {
  field(DESC,  "Mechanical‐mode Qs")
  field(NELM,  "5")
  field(DTYP,  "Soft Channel")
  field(VAL,   "40,20,50,80,100")
}

record(waveform, "CAV:MECH_K") {
  field(DESC,  "Mechanical‐mode couplings")
  field(NELM,  "5")
  field(DTYP,  "Soft Channel")
  field(VAL,   "2,0.8,2,0.6,0.2")
}


# === 3) Internal state variables ===

record(longout, "CAV:PHA_SCAL") {
  field(DESC,  "RF source phase (rad)")
  field(PREC,  6)
  field(VAL,   "0.0")
}

record(ao, "CAV:STATE_VC_R") {
  field(DESC,  "Real part of Vc (V)")
  field(PREC,  6)
  field(VAL,   "0.0")
}

record(ao, "CAV:STATE_VC_I") {
  field(DESC,  "Imag part of Vc (V)")
  field(PREC,  6)
  field(VAL,   "0.0")
}

record(waveform, "CAV:STATE_M") {
  field(DESC,  "Mechanical state vector [x₁…x₅]")
  field(NELM,  "5")
  field(DTYP,  "Soft Channel")
  field(VAL,   "0,0,0,0,0")
}

record(ao, "CAV:DW_PREV") {
  field(DESC,  "Detuning from previous step (rad/s)")
  field(EGU,   "rad/s")
  field(PREC,  3)
  field(VAL,   "0.0")
}


# === 4) Output signals (written each TS by SNL) ===

record(ao, "CAV:SIG_VC_R") {
  field(DESC,  "Current real(Vc) (V)")
  field(EGU,   "V")
  field(PREC,  6)
  field(VAL,   "0.0")
}

record(ao, "CAV:SIG_VC_I") {
  field(DESC,  "Current imag(Vc) (V)")
  field(EGU,   "V")
  field(PREC,  6)
  field(VAL,   "0.0")
}

record(ao, "CAV:SIG_VR_R") {
  field(DESC,  "Current real(Vr) (V)")
  field(EGU,   "V")
  field(PREC,  6)
  field(VAL,   "0.0")
}

record(ao, "CAV:SIG_VR_I") {
  field(DESC,  "Current imag(Vr) (V)")
  field(EGU,   "V")
  field(PREC,  6)
  field(VAL,   "0.0")
}

record(ao, "CAV:SIG_DW") {
  field(DESC,  "Current detuning (Hz)")
  field(EGU,   "Hz")
  field(PREC,  3)
  field(VAL,   "0.0")
}

Notes on the .db file above

  • The VAL lines for CAV:BASE_PUL and CAV:BEAM_PUL must be filled with exactly 16 384 comma-separated numbers. In practice, you can generate these arrays offline (e.g. with Python), then paste them into the .db file.
  • All “AO” records with VAL set to “0.0” or some initial number can be edited at run-time by Phoebus.
  • We chose to store STATE_M as a length-5 waveform because there are 5 mechanical modes. If you prefer, you can split STATE_M into 5 separate AO records (CAV:STATE_M1, …, CAV:STATE_M5), but a single waveform is tidier.
  • CAV:PHA_SCAL is a longout so it can hold the continuously‐increasing phase (rad). We kept that as a scalar in SNL.
  • The “output” PVs CAV:SIG_VC_R/I, CAV:SIG_VR_R/I, and CAV:SIG_DW are AO records that will be written on every Ts by the sequencer.

Write the Sequencer program (SNL/cavity_sim.st)

You now need an SNL file that:

  1. Runs once at IOC startup:

    • Reads in the mechanical‐mode arrays (MECH_F, MECH_Q, MECH_K).
    • Computes the continuous-time state-space (Aₘ, Bₘ, Cₘ, Dₘ) for the mechanical part (just like cav_ss_mech(mech_modes) in Python).
    • Discretizes it ((A_d, B_d, C_d, D_d) = ss_discrete(Aₘ, Bₘ, Cₘ, Dₘ, Ts)).
    • Computes constants like wh = π·f₀/QL and RL = 0.5·(r/Q)·QL.
    • Initializes internal state (pha = 0, state_vc = 0 + j·0, state_m = [0 … 0]^T, dw_prev = DW0).

  2. Steps at each “time step” (every Ts):

    • Reads fresh parameter PVs: TS, FSRC, ASRC, GAIN_DB, BETA, RoQ, QL, IB, DW0, PULSED.

    • Advances the RF phase:

      1
      pha = pha + 2·π·fsrc·Ts

    • Computes the complex RF source S0 = Asrc · exp(i·pha), i.e.

      1
      2
      S0_R = Asrc * cos(pha)
      S0_I = Asrc * sin(pha)

    • Chooses “pulsed” or “CW” I/Q modulator output:

      1
      2
      3
      4
      5
      6
      7
      idx = buf_id mod 16384
      if (PULSED == 1) then
        base_coeff = BASE_PUL[idx]        # 0 or 1
      else
        base_coeff = 1.0
      S1_R = S0_R * base_coeff
      S1_I = S0_I * base_coeff

    • Amplifier:

      1
      2
      3
      gain_lin = 10^(GAIN_DB / 20)
      S2_R = S1_R * gain_lin
      S2_I = S1_I * gain_lin

    • Microphonic detuning:

      1
      2
      3
      4
      5
      # SNL’s seq function rand( ) returns uniform rand in [0,1). For Gaussian, you can approximate
      # using Box–Muller if you want. For simplicity, we’ll do small uniform random:
      rnd = rand(1) - 0.5        ; # uniform in [–0.5, +0.5)
      dw_micr = 2·π·rnd·20       ; # scale so sigma ≃ 10 Hz
      dw_tot  = dw_prev + dw_micr  ; # total detuning (rad/s)

      If you need a true Gaussian, you can code a Box–Muller inside SNL (two calls to rand()).

    • Cavity + mechanical update: call a custom subroutine that implements your discrete‐time cavity ODE + mechanical ODE for one step. In Python you used:

      1
      2
      3
      4
      status, vc, vr, dw, state_vc, state_m = sim_scav_step(
          wh, dw_prev, dw0 + dw_micr, S2, state_vc, Ts,
          β = beta, state_m0 = state_m, Am = A_d, Bm = B_d,
          Cm = C_d, Dm = D_d, mech_exe = True)

      In EPICS/SNL you cannot call a Python function. Instead you have two options:

      1. Write a C subroutine (linked into Sequencer) that replicates sim_scav_step in its entirety (i.e. the discrete‐time cavity/mechanical update). Let’s name it scav_step(…). You compile it as part of your IOC, then call it from SNL:

        1
        2
        3
        4
        5
        6
        7
        8
        9
        10
        11
        12
        13
        14
        15
        16
        17
        18
        19
        20
        21
        22
        23
        24
        25
        26
        27
        28
        /* SNL_C_subr.c */
        #include <math.h>
        #include "cav_ss.h"   /* header for your discrete matrices, maybe generated */
        void scav_step(double wh, double dw_prev, double dw0_total,
                       double VF_R, double VF_I,
                       double SVc_R, double SVc_I, double Ts,
                       double beta,
                       double stateM[], /* length = 5 */
                       double A_d[][5], double B_d[][5], double C_d[][5], double D_d[][5],
                       double *newVc_R, double *newVc_I,
                       double *newVr_R, double *newVr_I,
                       double *newDw,   /* rad/s */
                       double newStateM[])
        {
            /* 
              Implement exactly what sim_scav_step does:
              - Compute cavity differential equation:
                  dVc/dt = –(ωh + j·(dw_prev+dw0_total))·Vc + (2·ωh)·S2 – (2·ωh/β)·Vc – (2·ωh/β)·0
                Then discretize by Ts or use the A_d,B_d etc. outputs for state_m update.
              - Update mechanical state:  x[k+1] = A_d * x[k] + B_d * detuning
              - Compute detuning (dw) from mechanical motion:  dw_mech = C_d * x[k] + D_d * detuning_in
              - Compute new Vc, Vr from formulas in sim_scav_step
              - Return newVc_R, newVc_I, newVr_R, newVr_I, newDw, newStateM[]

              In practice you can literally copy the code from llrflibs/… 
              and adapt it to pure C.                                                
            */
        }

        Then, in your SNL cavity_sim.st, you declare:

        1
        2
        3
        4
        5
        6
        7
        8
        9
        10
        11
        12
        13
        14
        15
        16
        17
        18
        19
        20
        21
        22
        23
        24
        25
        26
        27
        28
        29
        30
        31
        32
        33
        34
        35
        36
        37
        38
        39
        40
        41
        42
        43
        44
        45
        46
        47
        48
        49
        50
        51
        52
        53
        54
        55
        56
        57
        58
        59
        60
        61
        62
        63
        64
        65
        66
        67
        68
        69
        70
        71
        72
        73
        74
        75
        76
        77
        78
        79
        80
        81
        82
        83
        84
        85
        86
        87
        88
        89
        90
        91
        92
        93
        94
        95
        96
        97
        98
        99
        100
        101
        102
        103
        104
        105
        106
        107
        108
        109
        110
        111
        112
        113
        114
        115
        116
        117
        118
        119
        120
        121
        122
        123
        124
        125
        126
        127
        128
        129
        130
        131
        132
        133
        134
        135
        136
        137
        138
        139
        140
        141
        142
        143
        144
        145
        146
        147
        148
        149
        150
        151
        152
        153
        154
        155
        156
        157
        158
        159
        160
        161
        162
        163
        164
        165
        166
        167
        168
        169
        170
        171
        172
        #include "epicsEvent.h"
        #include "cav_ss.h"    /* where you declared the subroutine prototype */
        #include "cadef.h"

        epicsShareFunc
        void scav_step(double wh, double dw_prev, double dw0_total,
                       double VF_R,  double VF_I,
                       double SVc_R, double SVc_I,
                       double Ts, double beta,
                       double stateM[], /* length = 5 */
                       double A_d[][5], double B_d[][5], double C_d[][5], double D_d[][5],
                       double *newVc_R, double *newVc_I,
                       double *newVr_R, double *newVr_I,
                       double *newDw,   /* rad/s */
                       double newStateM[]);

        program cavity_sim
        {
          double pha, Ts, fsrc, Asrc, gain_dB, beta, rOQ, QL, IB, dw0, wh, RL;
          int pulsed;
          double base_val, beam_val;
          int buf_id;
          double S0_R, S0_I, S1_R, S1_I, S2_R, S2_I;
          double dw_prev, dw_tot, dw_micr;
          double VC_R, VC_I, VR_R, VR_I, newVc_R, newVc_I, newVr_R, newVr_I, newDw;
          double state_vc_R, state_vc_I;
          double stateM[5], newStateM[5];
          double A_d[5][5], B_d[5][5], C_d[5][5], D_d[5][5];
          double mech_f[5], mech_Q[5], mech_K[5];
          int i;
          
          double tmpArray[16384];   /* for reading waveform values */
          
          /* --------------- Initialization --------------- */
          /* 1) Read constant PVs once (will refresh inside loop too) */
          dvGet("CAV:TS",       &Ts);
          dvGet("CAV:FSRC",     &fsrc);
          dvGet("CAV:ASRC",     &Asrc);
          dvGet("CAV:GAIN_DB",  &gain_dB);
          dvGet("CAV:BETA",     &beta);
          dvGet("CAV:RoQ",      &rOQ);
          dvGet("CAV:QL",       &QL);
          dvGet("CAV:IB",       &IB);
          dvGet("CAV:DW0",      &dw0);
          dvGet("CAV:PULSED",   &pulsed);

          /* 2) Read mechanical‐mode arrays */
          dvGet("CAV:MECH_F", mech_f);
          dvGet("CAV:MECH_Q", mech_Q);
          dvGet("CAV:MECH_K", mech_K);

          /* 3) Compute RL and wh */
          RL = 0.5 * rOQ * QL;
          wh = 3.141592653589793 * (1.3e9) / QL;  # f0=1.3GHz
          wh = M_PI * (1.3e9) / QL;  # f0=1.3GHz
          /* 4) Build continuous‐time A_m, B_m, C_m, D_m for mech modes, then discretize.
           *    You must write a helper C function in 'cav_ss.h/.c' that does exactly:
           *      cav_ss_mech(...) → returns (A_m, B_m, C_m, D_m), 
           *      then ss_discrete((A_m,B_m,C_m,D_m), Ts) → (A_d, B_d, C_d, D_d).
           *    For brevity, assume we call:
           */
          compute_ss_discrete(mech_f, mech_Q, mech_K, 5, Ts, A_d, B_d, C_d, D_d);

          /* 5) Initialize internal state */
          pha        = 0.0;
          state_vc_R = 0.0;
          state_vc_I = 0.0;
          for(i=0; i<5; i++) {
              stateM[i] = 0.0;
          }
          dw_prev = dw0;   /* initial detuning (rad/s)  */

          buf_id = 0;      /* start at sample 0 */

          /* 6) Populate BASE_PUL buffer: (we could read it from PV each time; instead:
           *    BASE_PUL[i] = 1 for i<t_flat (1300), else 0. We do that once:
           */
          dvGetWaveform("CAV:BASE_PUL", tmpArray, 16384);
          /* But since BASE_PUL is static in the .db file, we actually do not need to re‐fill it;
           * the .db value is already set. */

          /* 7) Similarly, BEAM_PUL is loaded from .db at IOC startup. */

          /* --------------- Real‐time TS loop --------------- */
          while (TRUE) {
            /***** 1) Re‐read dynamic PVs so user can adjust on the fly *****/
            dvGet("CAV:TS",       &Ts);
            dvGet("CAV:FSRC",     &fsrc);
            dvGet("CAV:ASRC",     &Asrc);
            dvGet("CAV:GAIN_DB",  &gain_dB);
            dvGet("CAV:BETA",     &beta);
            dvGet("CAV:RoQ",      &rOQ);
            dvGet("CAV:QL",       &QL);
            dvGet("CAV:IB",       &IB);
            dvGet("CAV:DW0",      &dw0);
            dvGet("CAV:PULSED",   &pulsed);

            /* 2) Advance RF phase */
            pha = pha + 2.0 * 3.141592653589793 * fsrc * Ts;
            dvPut("CAV:PHA_SCAL", pha);

            /* 3) Compute S0 = Asrc * exp(i·pha) */
            S0_R = Asrc * cos(pha);
            S0_I = Asrc * sin(pha);

            /* 4) I/Q modulator: if pulsed, multiply by BASE_PUL[buf_id] */
            /*    We read one element from the BASE_PUL waveform PV: */
            dvGetWaveformElement("CAV:BASE_PUL", buf_id, &base_val);
            if (pulsed == 1) {
              S1_R = S0_R * base_val;
              S1_I = S0_I * base_val;
            } else {
              S1_R = S0_R;
              S1_I = S0_I;
            }

            /* 5) Amplifier: gain_lin = 10^(gain_dB/20) */
            double gain_lin = pow(10.0, gain_dB / 20.0);
            S2_R = S1_R * gain_lin;
            S2_I = S1_I * gain_lin;

            /* 6) Microphonics: uniform random approx (for Gaussian, use Box-Muller) */
            double rnd1 = rand() / (RAND_MAX + 1.0);
            double rnd2 = rand() / (RAND_MAX + 1.0);
            /* Box–Muller for zero‐mean unit‐std dev: */
            double gm = sqrt(-2.0 * log(rnd1)) * cos(2.0 * 3.141592653589793 * rnd2);
            dw_micr = 2.0 * 3.141592653589793 * gm * 10.0;  # σ=10 Hz
            dw_tot = dw_prev + dw_micr + dw0;

            /* 7) Fetch current cavity state */
            dvGet("CAV:STATE_VC_R", &state_vc_R);
            dvGet("CAV:STATE_VC_I", &state_vc_I);
            dvGet("CAV:DW_PREV",    &dw_prev);
            dvGetWaveform("CAV:STATE_M", stateM, 5);

            /* 8) Call custom C routine that does one discrete‐step: */
            scav_step(wh, dw_prev, dw_tot,
                      S2_R, S2_I,
                      state_vc_R, state_vc_I,
                      Ts, beta,
                      stateM, A_d, B_d, C_d, D_d,
                      &newVc_R, &newVc_I,
                      &newVr_R, &newVr_I,
                      &newDw,     /* rad/s */
                      newStateM);

            /* 9) Write back new state_vc, new state_m, new dw */
            dvPut("CAV:STATE_VC_R", newVc_R);
            dvPut("CAV:STATE_VC_I", newVc_I);

            /*    We only store dw_prev = newDw, so next iteration uses it: */
            dvPut("CAV:DW_PREV", newDw);

            /*    Write out new mechanical state vector: */
            dvPutWaveform("CAV:STATE_M", newStateM, 5);

            /* 10) Write “live” output PVs so Phoebus can trend them: */
            dvPut("CAV:SIG_VC_R", newVc_R);
            dvPut("CAV:SIG_VC_I", newVc_I);
            dvPut("CAV:SIG_VR_R", newVr_R);
            dvPut("CAV:SIG_VR_I", newVr_I);
            dvPut("CAV:SIG_DW",   (newDw / (2.0 * 3.141592653589793)));  # convert to Hz

            /* 11) Increment buffer index, wrap around at 16384 */
            buf_id++;
            if (buf_id >= 16384) buf_id = 0;

            /* 12) Delay for Ts seconds */
            delay(Ts);
          } /* end while(TRUE) */

        } /* end program cavity_sim */

Explanation of key parts:

  • compute_ss_discrete(…)
    You must write a helper C function (e.g. in cav_ss.c) that takes the mechanical mode arrays (mech_f[], mech_Q[], mech_K[]) and returns the continuous‐time state‐space (Aₘ, Bₘ, Cₘ, Dₘ) for those five modes, then discretizes them at the current Ts (using a zero‐order hold). In Python you used cav_ss_mech(...) followed by ss_discrete(...). In C, you can port that code (there are simple textbook formulas for a second‐order mechanical resonance) or literally translate the Python logic from llrflibs.rf_sim.

  • scav_step(…)
    This is where the heart of sim_scav_step lives. It must implement exactly what Python’s sim_scav_step does:

    1. Compute the new mechanical detuning (dw_mech) from C_d * x[k] + D_d * detuning_in.
    2. Solve the discrete‐time cavity ODE (with half‐bandwidth wh, coupling β, RL, etc.) to get vc[k+1] from vc[k] and the new drive vf_step (S2).
    3. Compute vr (reflected voltage) from vc and vf, typically vr = (2·ωh/β)·Vc – Vf or however your convention is.
    4. Return newVc_R, newVc_I, newVr_R, newVr_I, newDw (rad/s), and the new mechanical state newStateM[].

    The exact formulas appear inside llrflibs/rf_sim.py, so you must port them to C. We assume you create <cav_ss.h> and <cav_ss.c> that supply both compute_ss_discrete(...) and scav_step(...). Then compile them as part of your IOC (adjust the Makefile to add these sources).

  • dvGet / dvPut / dvGetWaveform / dvPutWaveform
    These are SNL built-in calls that read or write EPICS PVs. For example, dvGet("CAV:TS", &Ts) fetches the current value of CAV:TS into the local SNL variable Ts.

  • delay(Ts)
    Causes the SNL program to sleep for Ts seconds. If Ts = 1e-6, the sequencer will loop at 1 MHz. In practice, runs at 1 MHz are very tight; watch CPU load. If you do not need a full 1 MHz scan (for example, if you can tolerate 100 kHz), set CAV:TS = 1e-5 (10 µs) to lighten the CPU.

Write or generate the C helper code (cav_ss.c and cav_ss.h)

You need two pieces of C code:

  1. compute_ss_discrete(…)

    • Input:

      • mech_f[5] (Hz),
      • mech_Q[5],
      • mech_K[5],
      • num_modes = 5,
      • Ts (s).

    • Output (passed by pointer):

      • A_d[5][5], B_d[5][5], C_d[5][5], D_d[5][5]

    These are the discretized mechanical‐mode state matrices, just like your Python did with ss_discrete(...). The continuous system for each mode m is:

    If you bundle all 5 modes into a 10×10 A_m, B_m, C_m, D_m system, you can build it in C, then do the matrix‐exponential (zero‐order hold) to get A_d = exp(A_m·Ts), and so on. A simpler approach is to use a Taylor expansion if Ts is small. Whichever route you choose, make sure your C code returns exactly the same A_d, B_d, C_d, D_d that Python’s ss_discrete would.

  2. scav_step(…)

    • Input: all the arguments listed in the SNL prototype above.
    • Behavior: compute one discrete‐time update of the cavity+mechanical system.
    • Output: new Vc_R, Vc_I, Vr_R, Vr_I, dw, newStateM[].

    The formulas come directly from llrflibs.rf_sim:sim_scav_step. You must translate them line by line into C.

Compile both cav_ss.c and scav_step.c (or keep them together) into an object library and link them with your IOC. In your IOC’s Makefile, add:

1
2
3
4
5
6
7
8
9
PROD_IOC = cavityIoc
# List the C source files you need:
cavityIoc_SRCS += \
    cav_ss.c \
    scav_step.c

# Sequencer files:
SNCFLAGS += +r
cavityIoc_LIBS += seq pv

Adjust as needed so that when you run make, you get an IOC executable that contains both the SNL program and your C routines.

Compile and launch the Soft IOC

  1. At the top of cavityIoc/Makefile, ensure you have:

    1
    2
    TOP = .
    include $(EPICS_BASE)/configure/RULES_TOP

  2. From inside cavityIoc/, run:

    1
    make

    This should compile your SNL (cavity_sim.st → cavity_sim.c → cavity_sim.o) and link it with cav_ss.o and scav_step.o into an IOC executable (e.g. cavityIoc).

  3. Start the IOC by running:

    1
    2
    cd iocBoot
    ./st.cmd

    You should see output indicating EPICS is loading cavity_sim.db, then printing “Sequencer loaded: cavity_sim” and “Program cavity_sim started”.

  4. Check that the PVs exist. In another terminal, use:

    1
    caget CAV:TS CAV:FSRC CAV:ASRC CAV:SIG_VC_R CAV:SIG_DW

    You should get back their initial values (e.g. 1e-6, -460.0, 1.0, 0.0, 0.0). If Euler’s algorithm or the discrete integration is running, after a few milliseconds you’ll see CAV:SIG_VC_R jump off zero (once the loop executes).

Build a Phoebus display

Once the IOC is running and streaming the “live” PVs (updated each Ts), open Phoebus and:

  1. Create a new Screen.

  2. Add a “Stripchart” widget (Trend) and bind it to CAV:SIG_VC_R (y-axis) vs. time.

  3. Add another Trend for CAV:SIG_VC_I. Or combine both real/imag into one “Complex” plot if you like.

  4. Add a third Trend for CAV:SIG_DW (in Hz).

  5. Add numeric input fields (Text‐Entry or Slider) bound to:

    • CAV:TS (so you can increase or decrease step size)
    • CAV:FSRC (offset frequency)
    • CAV:ASRC (source amplitude)
    • CAV:GAIN_DB (amplifier gain)
    • CAV:IB (beam current)
    • CAV:PULSED (you can switch CW vs. pulsed)
    • CAV:RoQ, CAV:QL, CAV:BETA if you want to change cavity loading at run‐time.

  6. Optionally add read-only display fields (Read Only Text) showing CAV:PHA_SCAL, so you can watch the RF phase drift.

Now, as you hit “Apply” or move the slider, the Sequencer loop will read the new PV value on its next iteration and immediately incorporate it into the next scav_step. Your Phoebus trends will show how vc, vr, and dw evolve over time.

Recap of “Whole Steps”

To summarize, here is the step-by-step checklist you must follow:

  1. Install EPICS Base, iocCore, and Sequencer on your development machine.

  2. Create an IOC directory structure (with db/ and SNL/ subfolders).

  3. Write db/cavity_sim.db containing:

    • AO records for static parameters (Ts, fsrc, Asrc, gain_dB, β, r/Q, QL, IB, DW0, PULSED).
    • Waveform records for the 16 384-point arrays (BASE_PUL, BEAM_PUL) and the 5-element mechanical‐mode arrays (MECH_F, MECH_Q, MECH_K).
    • Records for internal state: CAV:PHA_SCAL (longout), CAV:STATE_VC_R/I, CAV:STATE_M[5] (waveform), CAV:DW_PREV (AO).
    • Output PVs (AO) for CAV:SIG_VC_R/I, CAV:SIG_VR_R/I, CAV:SIG_DW.

  4. Write C helper code (cav_ss.c/.h and scav_step.c/.h):

    • compute_ss_discrete(...) to build and discretize the mechanical‐mode state matrices.
    • scav_step(...) to implement one discrete‐time cavity + mechanical update (port from Python’s sim_scav_step).

  5. Write the Sequencer program (SNL/cavity_sim.st) that:

    1. Imports your C routines.

    2. In initialize{} or at the top of program cavity_sim, reads the mechanical arrays, computes A_d, B_d, …, initializes state (pha, state_vc, state_m, dw_prev, buf_id).

    3. Enters a while(TRUE) loop that runs once per Ts:

      • dvGet(…) for all parameter PVs.
      • Advances the RF phase: pha += 2π·fsrc·Ts.
      • Computes S0 = Asrc·exp(i·pha), “I/Q modulate” by reading CAV:BASE_PUL[buf_id] if pulsed, then amplify by 10^(gain_dB/20).
      • Creates a small microphonic dw_micr (e.g. Box–Muller for Gaussian noise).
      • Calls your C routine scav_step(...) to update vc, vr, dw and state_m.
      • dvPut(...) the new state back into CAV:STATE_VC_R/I, CAV:STATE_M, CAV:DW_PREV.
      • dvPut(...) the “live” outputs into CAV:SIG_VC_R/I, CAV:SIG_VR_R/I, CAV:SIG_DW.
      • Increment buf_id modulo 16 384.
      • delay(Ts).

  6. Edit your IOC’s Makefile so that it:

    • Includes your C source files (cav_ss.c, scav_step.c, plus the SNL output).
    • Lists db/cavity_sim.db under DB +=
    • Lists SNL/cavity_sim.st under SNLFILES += and sets EMIT += cavity_sim.
    • Ensures that on make you end up with a working IOC named, say, cavityIoc.

  7. Write iocBoot/st.cmd (or equivalent) that:

    • Launches softIoc -d db/cavity_sim.db.
    • Loads and starts the Sequencer (sequencer -k -l -c "cavity_sim.st").

  8. Run make; ./iocBoot/st.cmd. Confirm the IOC is up by using cainfo or caget on a few PVs (CAV:TS, CAV:SIG_VC_R, etc.).

  9. Open Phoebus, create a new screen, and add:

    • Slider or Numeric Entry widgets bound to your user-tunable AO PVs (CAV:TS, CAV:FSRC, CAV:ASRC, CAV:GAIN_DB, CAV:IB, etc.).
    • Stripchart/Trend widgets bound to CAV:SIG_VC_R, CAV:SIG_VC_I, and CAV:SIG_DW.
    • Status displays (Text) for things like CAV:PHA_SCAL (RF phase) or CAV:STATE_M (mechanical states) if desired.

  10. Test interactivity: adjust CAV:FSRC or CAV:GAIN_DB on the Phoebus panel and watch CAV:SIG_VC_R respond in real time. Switch CAV:PULSED between 0/1 to see pulsed vs. CW behavior.

Tips and Caveats

  • Timing resolution: A sequencer loop at Ts = 1 µs is very tight. Real-time OS jitter can affect accuracy. If you don’t truly need 1 MHz updates, consider Ts = 10 µs (100 kHz) to lighten CPU load.
  • Waveform data sizes: If you actually want to store the entire sig_vc[16384] as a waveform, declare record(waveform, "CAV:SIG_VC") { NELM = 16384 } and write each new sample into SIG_VC[buf_id] via dvPutWaveformElement("CAV:SIG_VC", buf_id, newVc_R + i·newVc_I). But this can fill buffers quickly. For simple live plots, updating a single “latest point” PV is easier.
  • Gaussian vs. Uniform noise: Python’s np.random.randn() yields true Gaussian. In SNL, you must implement your own Box–Muller if you need Gaussian microphonics. Otherwise, a small uniform random (e.g. (rand()–0.5)*2π*10) may suffice.
  • Discretization of mechanical modes: If you do not want to code a full compute_ss_discrete in C, you can precompute the 5×5 discrete matrices offline (once) for your nominal Ts = 1 µs, then paste them into the SNL as constants. However, if the user changes CAV:TS at run-time, those matrices must be recomputed. In that case, it’s best to keep the C routine that recomputes A_d,B_d,C_d,D_d whenever Ts changes.
  • Complex arithmetic: We keep real/imag parts in separate PVs always. Never try to concatenate them into one.
  • Memory layout: SNL’s double stateM[5] is a local array in the while(TRUE) loop. We fetch its contents from CAV:STATE_M (waveform) at the start of each iteration, then write the new stateM back at the end. That ensures mechanical state is persistent.
  • Buffer indices: We keep a local int buf_id that tracks which sample from the waveform buffers (BASE_PUL, BEAM_PUL) to pick this iteration. Wrap it at 16 384 so you cycle continuously.

In summary

  • You have now mapped every single variable from the Python simulation into EPICS PVs and records:

    1. Static parameters (AO, MBBI).
    2. Buffers (WAVEFORM).
    3. Internal state (LONGOUT, AO, WAVEFORM).
    4. Outputs (AO).

  • You wrote an SNL program that runs exactly the same logic as your Python loop, except each operation is either:

    1. A direct C call in scav_step(...), or
    2. Arithmetic inside SNL (phase update, amplifier gain, pulsed flags).

  • At each Ts, you read from EPICS PVs (dvGet), update internal variables, call scav_step, then dvPut back to PVs.

  • Phoebus can subscribe to those PVs and display live‐updating cavity voltage, phase, and detuning.

Once this is in place, you have a fully EPICS‐native Soft IOC that reproduces your original Python simulation, but now every parameter and every dynamic output is exposed as a PV—and can be driven or displayed by Phoebus (or CSS, EDM, CSS-BOY, etc.). Maintenance is harder than a Python‐centric IOC, but if your lab policy requires “no Python inside the IOC,” this is the standard EPICS way to do it.