\documentclass[acus]{tech-lrd}
\usepackage{epsfig}

\begin{document}
\technotenumber{FE-EE-021}
\technoteseries{FE3000}
\technotedate{June 13, 2001}

\title{Low Level RF Modeling for the MEBT Rebunchers}
\author{Lawrence R. Doolittle, LBNL}
\setlength{\titleblockheight}{3.5cm}
\maketitle

{\it
 \small
 This technical note is intended exclusively for internal
use inside the SNS Front-End Systems Group.  The primary purpose of the
technical note series is to promote rapid dissemination of results for
work in progress.  No formal review has taken place to ascertain the
accuracy of its contents, or the consistency/compatibility of the information
with other SNS work.}

\section{Introduction}

The four MEBT rebuncher cavities need control electronics to 
set their amplitude and phase to the values needed for proper
beam propagation.
The raw power to drive them (20\thinspace kW to 40\thinspace kW
at 402.5\thinspace MHz) comes from
commercially packaged grounded-grid triode amplifiers.
The ``Low Level RF'' subsystem discussed in this note needs to:
\begin{itemize}
\item  send an RF waveform to the triode assembly
\item  measure the cavity detune angle (tuner error signal)
\item  use EPICS for command and readback.
\end{itemize}
An architecture for this subsystem was reviewed on May 10, 2001\cite{llrf-review}.
A digital signal processing section (based on an FPGA) is connected
to a low latency measurement of the cavity field vector, and controls
the RF drive vector sent to the triode amplifier.

At its simplest, the digital control unit could turn the RF output
on for 1\thinspace ms after a 60\thinspace Hz trigger,
and provide a readback of the
cavity voltage waveform to the operator.  Previous experience with
the RFQ shows that the decay segment of this waveform can be analyzed 
to determine the cavity detune angle.  The benefit of having a low-latency
measurement/control loop is that digital feedback can be included 
with no additional hardware.

The body of this paper explains a numerical model
of the RF system, and shows the expected behavior of the system.
It shows how feedback and feedforward can be used
to keep the beam loading transient under control.

I also provide documentation for how to work with the numerical
model, which is provided in the form of a {\tt c} program.  The
program, along with this tech note, are available for download\cite{recycle}.

\section{Model setup}

\begin{figure*}[tb]
\centering
\epsfig{file=block.eps, width=0.9\hsize}
\caption{Block diagram of the system as modeled.}
\label{block-model}
\end{figure*}

\begin{figure}[tb]
\centering
\epsfig{file=feedforward.eps, width=0.9\hsize}
\caption{Schematic of feedforward waveform structure.  T2 is adjusted
to match the arrival of the first beam.}
\label{feedforward}
\end{figure}

The block diagram of the RF system model is shown in figure~\ref{block-model}.
An explanation of the feedforward terms
is included in figure~\ref{feedforward}.
The whole loop is modeled at baseband, leaving out all the mixers
and frequency conversion.  This means that all signals are represented
in vector form, and that filter poles are not restricted to 
negative real and conjugate pairs.
In particular, detuning a cavity takes its pole position from a nominal
-1/$\omega_0$ slightly into the imaginary direction.  This shift properly
introduces a coupling between the otherwise linearly independent
vector components of the drive signal.

All the low-pass filters are set up to have low frequency gain of 1.
Thus, the numerical amplitudes used are given in suitably normalized
units.

A unity signal (1,0) coming from the amplifier corresponds to full output,
(20\thinspace kW for cavities~1-3, 40\thinspace kW for cavity~4).
I define the small signal amplifier gain
as unity, which implicitly defines the input voltage scaling.

\begin{figure}[tb]
\centering
\epsfig{file=qei3.eps, width=0.9\hsize}
\caption{Triode amplifier saturation characteristics.}
\label{qei3}
\end{figure}

Triode clipping is described with a harshness parameter $c$, such that the
output signal $v_o$ from the tube based on its input $v_i$ varies according to
$ v_o = v_i  \cdot ( 1+|v_i|^c)^{-1/c} $.
The saturated output amplitude from this equation is 1.
While some phase shift with drive level is observed, this effect
is not yet included in the model.
When $c=5$, as provides a decent fit to measurements shown in
figure~\ref{qei3}, it takes about 1.5 (unitless) drive level from
the control system to reach 97.5\% of this amplitude, 95\% of the power.
This level should correspond to full scale digital output.

% and we
% should put in some kind of clipping on the diagonal, although since we
% are so far into saturation it doesn't make much difference in the final
% output level to the cavity.  It would just be nice to limit both the power
% in the drive stages, and the wind-up in the control algorithms.

\def\Zshunt{Z_{\rm shunt}}
For the cavity, unity is the equilibrium voltage on resonance
given the rated amplifier power.
Rebuncher cavity $\Zshunt T^2/L$ from the Superfish runs\cite{belly2}
is based on a 2.72\thinspace cm length, which is the half-length
of the cell used in the Superfish analysis.
The scaling from current to beam-referred cavity voltage is given by
$\Zshunt T^2$ for the whole cavity (not the half cell), which
involves another multiplication by two.  Furthermore, the computed
value of $\Zshunt$ is based on a cavity $Q$ of~21700.  This needs
to be scaled to the loaded cavity $Q$, estimated (80\% theoretical~$Q_0$,
critical coupling) as~8700.
Apply this impedance
to unity cavity voltage to get unity beam current.

%Let's work with cavity 1, a 30mm bore cavity driven by a single
%triode.
%  Output power level normalized to 20 kW
%  Equilibrium cavity voltage for this power input is 75 keV * sqrt(20/11)
%      = 101 keV, this is "1" for cavity voltage.  Our normalized setpoint is
%      therefore sqrt(11/20) = 0.74.
%  A 38 mA beam induces .038 * 639 = 24.3 keV on the cavity.
%  So to convert beam current to something with the same "units" as the
%      drive amplitude from the triode, divide by 158 mA, so the nominal
%      38 mA converts to 0.24.

The following table shows the key cavity properties
that are used to set up and interpret the model:
\halign to \hsize{\quad#\hfil~&#\hfil~&#\hfil~&#\hfil~&#\hfil~&#\hfil \cr
Cavity          &    1    &    2   &     3   &     4   &       \cr
Bore            &   30    &   36   &    36   &    36   &  mm   \cr
$\Zshunt T^2/L$ &  11.75  &  6.85  &  6.85   &   11.75 &  M$\Omega$/m for $Q_0$\cr
$\Zshunt T^2$   &  256    &  149   &   149   &    256  &  k$\Omega$ for $Q_{\rm L}$\cr
Unity Power     &   20    &   20   &    20   &    40   &  kW   \cr
Design Power    &   11.0  &   6.8  &    8.1  &    28.2 &  kW   \cr
Unity Voltage   &   101   &   77   &    77   &   143   &  keV  \cr
Design Voltage  &   75    &   45   &   49.3  &   120   &  keV  \cr
Setpoint        &   0.74  &   0.58 &    0.64 &   0.84  &       \cr
Unity Current   &   395   &   517  &    517  &   559   &  mA   \cr
}

Beam current {\it du jour} within the minipulse is
38\thinspace mA\cite{par-list},
which converts to normalized values between 0.07 and~0.10,
depending on the cavity.

\begin{figure}[tb]
\centering
\epsfig{file=rise2.eps, width=\hsize}
\caption{Triode amplifier risetime.}
\label{triode-step}
\end{figure}

Figure~\ref{triode-step} shows a measured step response of the amplifier
(at 60\% of full power), along with a calculated curve based on a filter
with poles at -9, -7.5+3j, and -7.5-3j (units are $(\mu {\rm s})^{-1}$),
plus 175\thinspace ns additional delay. 

We will use three pole, 2.5\thinspace MHz bandwidth
(at 50\thinspace MHz carrier)
receive (RX) filters to keep spurious signals out of the receiver.
The manufacturer provided magnitude and delay plots for the customized units.
I get good fits to these plots when I use pole locations 
-18, -9.5-15j, and -9.5+15j, in units of $(\mu {\rm s})^{-1}$.
The manufacturer's plots were faxed in, so I can't show those here.
I do show the results of my three-pole model in figure~\ref{filter-char}.

\begin{figure}[tb]
\centering
\epsfig{file=rxf1.eps, width=0.9\hsize}
\epsfig{file=rxf2.eps, width=0.9\hsize}
\epsfig{file=rxf3.eps, width=0.9\hsize}
\caption{RX Filter amplitude and delay characteristics.}
\label{filter-char}
\end{figure}

So far, the digital section is handled in the same manner as the analog.
The floating point numbers used for the amplitudes are effectively not
quantized.  The integrator accumulates IGRL times
the error signal at each time step.
No clipping as applied; amplifier saturation is
already included as discussed earlier, and the examples shown here do
not send the digital section beyond its planned output range.

The delays around the loop, grouped by physical component, are summarized here:

% 60pt is about centered
\centerline{\vbox{\halign{#\hfil\quad&\hfil#\quad&\hfil#\quad&\hfil#\cr
{\it component} & {\it delay,} ns
                           & DEL1 & DEL2 \cr
Cavity             &    -  &      &      \cr
Cable (32 ft.)     &   35  &      &  35  \cr
RX Filter          &  120  &      &      \cr
ADC (5 cycles)     &  125  &      & 125  \cr
Digital (4 cycles) &  100  &      & 100  \cr
DAC (1 cycle)      &   25  &      &  25  \cr
TX filter          &   50  &   50 &      \cr
Cable (100 ft.)    &  110  &  110 &      \cr
Amplifier          &  450  &  175 &      \cr
Cable (82 ft.)     &   90  &   90 &      \cr
Cavity             &    -  &      &      \cr
{\it Total}        & 1105  &  425 & 285  \cr
}}}
\noindent The DEL1 and DEL2 columns show the contribution to those model elements.
All of the RX filter delay and most of the amplifier delay is generated in the
filter elements, not the explicit delay elements.  Cable delays are based on
a velocity factor of~0.9, a reasonable estimate for low loss foam insulated coaxial cable.

\section{ODE Integration}
While there are many existing solvers for the numeric
Ordinary Differential Equation (ODE) problem,
they tend to be more general than needed, and don't fit well into the
simple architecture chosen for this software.  This architecture is
purposely simple for two reasons: speed (1\thinspace ms simulations at
40\thinspace MHz sample rate involve 40000 sample steps) and embeddability.
Chances are this first round of simulations will be grafted on
to a VHDL or Verilog synthesizable model of the FPGA.  While
packaged mixed-signal simulation software is all the rage these days,
I chose a known-good and {\it libre} framework from which to start,
so that the results can be reproduced without recourse to proprietary
black boxes.  Besides, the mathematics of ODE integration are
well understood.

Start with the first order differential equation for a single-pole
low pass filter.  This is expressed in Laplace form as
$$ {v_o\over v_i} = {1 \over s-p}~~,$$
where $v_i$ and $v_o$ are the input and output signals,
and $p$ is the pole location.
Make the differentiation explicit, and rearrange to get
a form consistent with state-variable numerical ODE integration,
$${dv_o\over dt} = v_i + p v_o~~.$$
The simplest expression for a `next' value $v_o'$ in a
discrete time {\it ansatz} is
$$v_o' = (1+\Delta t \cdot p)v_o + \Delta t \cdot v_i~~.$$
To improve convergence properties in the case where
$\Delta t \cdot p$ is not tiny,
approximate the trajectory of $v_i$ and $v_o$ as linear
within a single time step.
Specifically, assume that $v_o$ changes from $v_o$ to $v_o'$,
and $v_i$ changes from $v_i$ to \hbox{$v_i'$}.
Using what is essentially the Trapezoidal Formula\cite{trapezoid},
our rendition of the discrete time approximation to
the above differential equation becomes
\def\texthalf{{\textstyle {1\over 2}}}
$${v_o' - v_o\over \Delta t} = \texthalf (v_i+v_i') + p \texthalf(v_o+v_o')~~.$$
Solving for $v_o'$, we get
$$v_o' = a v_o + \texthalf b (v_i+v_i')$$
$$a={1+{1\over 2}\Delta t \cdot p \over 1 - {1\over 2} \Delta t \cdot p}~,\quad
  b={            \Delta t         \over 1 - {1\over 2} \Delta t \cdot p}~.$$

This process is coded in c, and tested using a two-pole low pass Butterworth
filter.  Given the transfer function $1/((s+1)^2+1)$, which has poles at
(-1,1) and (-1,-1),
the step response is $1-e^{-t}(\sin x + \cos x)$, for $t>0$.
This analytically known response is plotted in figure~\ref{analytic-check}
against two runs of the numeric ODE solver
with different values for the time step.
The first point after $t=0$ has a relatively large error, effectively due to
differing representations for the input signal near $t=0$.  After this point,
the error is bounded by 1.6\% for $|\Delta t \cdot p| = 0.57$,
and by 0.4\% for $|\Delta t \cdot p| = 0.28$.
The largest $|\Delta t \cdot p|$ used in the real model is 0.45.

\begin{figure}[tb]
\centering
\epsfig{file=analytic.eps, width=\hsize}
\caption{Cross-check of ODE code with analytic solution.}
\label{analytic-check}
\end{figure}

\section{Operating the code}

The {\tt step} program takes its numerical configuration from environment
variables.  The names below have a ``OP\_'' prepended (compile-time
adjustable) in an attempt to avoid namespace collisions.  The table below
summarizes the parameters the user has control of (without recompiling),
along with their default values.  

\noindent\halign{#\hfil\quad&#\hfil\quad&#\hfil\cr
{\it name} & {\it default} &   {\it short description} \cr
DT   & 0.025 &                 digital sample rate  \cr
NPT  & 2900 &                  \# of sample points  \cr
DEL1 & 17 &                    cables plus triode  \cr
TRF1 & -7.5,2.7 -7.5,-2.7 &    triode driver filter  \cr
TRF2 & -9,0 &                  triode output filter  \cr
CAV  & -.075,+.010 &           cavity filter  \cr
DEL2 & 12 &                    cables plus ADC  \cr
RXF  & -18,0 -9.5,-15 -9.5,15 & receive 50 MHz filter  \cr
NOIS & 0.003 &                 rms ADC noise  \cr
BEAM & 0,0.10 &                beam current vector  \cr
CLIP & 5.0 &                   triode clip harshness  \cr
T1   & 750 &                   length of initial ramp  \cr
T2   & 1410 &                  time of feedforward  \cr
DRV1 & 1.50,0.00 &             initial drive vector  \cr
DRV2 & 0.80,0.00 &             second drive vector  \cr
DRV3 & 0.00,-0.05 &            feedforward vector  \cr
SETP & 0.74,0.0 &              feedback setpoint  \cr
GAIN & -6.93,-1.0 &            feedback prop. gain  \cr
IGRL & 0.002 &                 feedback integral term  \cr
}

\noindent{\it Notes:}

Delays must be integer, and are expressed in $\Delta t$ units.

Vector quantities consist of the real and imaginary part,
separated by a comma (but not spaces).

Filters are specified by listing their complex pole locations, in
units of $1/$time.  The default setup uses $\mu$s for this and
\hbox{$\Delta t$}.  Multiple poles must be separated by spaces.

Input and output vectors always use scaled units, as discussed
in the previous section.

When the triode clipping harshness parameter $c$ is~0,
no saturation curve is applied.

The beam current is configurable, but the timing model is not (yet).
The timing model includes the 1\thinspace MHz chopper, and a ramp
from 31\% to 63\% duty factor over 15\thinspace $\mu$s.

The program output (written to {\tt stdout}) is nine columns wide.
The first column is the time,
the others come in pairs, for real and imaginary part.  The four
vectors, in order, are:  FPGA output, triode output, cavity voltage,
and receiver input.

\begin{figure*}[tb]
\centering
\epsfig{file=default.eps, width=0.8\hsize}
\caption{Beam pulse transient as modeled by default parameters.}
\label{default-graph}
\end{figure*}

\section{Performance}

The amplifier and cavity output for default
parameters is shown in figure~\ref{default-graph}.
The analysis runs at the rate of about 16000 time steps per second
on a 400\thinspace MHz Pentium-II computer.
The transient in the imaginary component of the drive signal
at $t$=20\thinspace $\mu$s comes from
the feedback loop correcting for a small amount (8$^\circ$) of
cavity detuning.  A more complex digital section could reduce this
transient by applying some form of feedback during the initial ramp.
Such changes will only involve firmware.

The beam current will induce a voltage on the cavity, an effect known as
beam loading.  This voltage change will not (to first order) effect
the bunching behavior of the cavity, but is phased to decelerate the beam.
For this reason, I show the absolute voltage, instead of the angle relative
to the bunching voltage.
Figure~\ref{imag} shows the scale of the effect for cavity~1, the ``worst''
of the four (higher TTF than cavities~2 and~3, lower stored energy than cavity~4).
The normalized output of the {\tt step} program has been multiplied by~101
to get actual, accelerating mode kV for the vertical scale, which is numerically
the same as the beam energy change in \hbox{keV}.  The sawtooth structure
of the cavity voltage is a result of the 1\thinspace MHz modulation of the
beam current.
Beam modeling results\cite{staples-beam} suggest that SNS will
run properly with induced energy changes in this cavity
as large as about 1.3\thinspace keV
(1$^\circ$ relative to the 75\thinspace keV bunching voltage).

\begin{figure}[tb]
\centering
\epsfig{file=imag.eps, width=\hsize}
\caption{Beam pulse transient for cavity 1.}
\label{imag}
\end{figure}

\section{Conclusions and Future Work}
The digital LLRF control architecture has the flexibility to meet the
needs of the rebuncher cavities.  Including simple proportional/integral
feedback, the transient from beam loading can be kept
smaller than the design tolerance of 1.3\thinspace \hbox{keV}.
Feedforward of a single step can cut the peak of that transient in half.

This model is not truly done, in that more detail can be added: quantization
and digital clipping would probably come next.
As the digital system is virtually assembled in VHDL or Verilog,
this model can be attached so as to provide a complete system simulation.

The feedback parameters for this model were adjusted by hand.
It is of considerable interest to write
software that will adjust the control loop without human intervention.

\begin{thebibliography}{9}

\bibitem{llrf-review}
LLRF Review, A. Ratti, FE-EE-020, June 2001.

\bibitem{recycle}
http://recycle.lbl.gov/$\sim$ldoolitt/llrf-model/

\bibitem{belly2}
Rebuncher Cavity Design Update, Lawrence R. Doolittle, FE-EE-014, March 2000.

\bibitem{par-list}
SNS Parameter List, SNS 100000000-PL0001-R05, May 2001

\bibitem{trapezoid}
Formula 25.5.3 of Handbook of Mathematical Functions, Abramowitz and Stegun,
1964.

\bibitem{staples-beam}
Front End Systems Beam Dynamics, J. Staples, LBNL-45289, January 2000.

\end{thebibliography}

\end{document}
