\documentclass[a4paper]{article}

\usepackage{amsfonts}
\usepackage{amsmath}
\usepackage[english]{babel}
\usepackage[squaren]{SIunits}
\usepackage[latin1]{inputenc}
\usepackage{graphicx}
\usepackage{hyperref}
\usepackage{epsfig}
\usepackage{subfig}
\usepackage{color}
\usepackage[]{algorithm2e}
\usepackage{verbatim}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Renew commands 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\newcommand{\de}{\partial}
\newcommand{\dd}{\mathrm{d}}

\newcommand{\er}{\mathbf{e_r}}
\newcommand{\ez}{\mathbf{e_z}}
\newcommand{\et}{\mathbf{e_\theta}}
\newcommand{\ep}{\mathbf{e_\phi}}
\newcommand{\B}{\mathbf{B}}
\newcommand{\Bp}{\mathbf{B_p}}
\newcommand{\dl}{\mathbf{dl_p}}

\newcommand{\E}{\mathbf{E}}
\renewcommand{\j}{\mathbf{j}}
\newcommand{\psih}{\hat{\psi}}

\newcommand{\be}{\begin{equation}}
\newcommand{\ee}{\end{equation}}

\newcommand{\ba}{\begin{align}}
\newcommand{\ea}{\end{align}}
\newcommand{\R}{\mathbb{R}}
\newcommand{\didi}[2]{\ensuremath{\frac{\dd{#1}}{\dd{#2}}}}
\newcommand{\dede}[2]{\ensuremath{\frac{\partial{#1}}{\partial{#2}}}}
\newcommand{\ddedde}[2]{\ensuremath{\frac{\partial^2{#1}}{\partial{#2}^2}}}
\newcommand{\psihat}{\hat{\psi}}
\newcommand{\Tehat}{\hat{T_e}}
\newcommand{\rhohat}{\hat{\rho}}
\newcommand{\psihb}{\hat{\boldsymbol{\psi}}}

\newcommand{\cmd}[1]{\textsf{#1}}

\newcommand{\Simulink}{\textsc{Simulink}}
\newcommand{\Matlab}{\textsc{Matlab}}

% units
\newcommand{\eV}{\electronvolt}
\newcommand{\keV}{\kilo\electronvolt}

% fonts
\newcommand{\bs}[1]{\boldsymbol{#1}}


%%%%%%%%%%%%%%%%%%%%%%

\title{Notes on coupling an equilibrium code to the RAPTOR suite of codes}
\author{F. Felici, federico.felici@epfl.ch}


\bibliographystyle{unsrt}

\begin{document}
\maketitle
Version 2.2

\tableofcontents

\section{Introduction}
This document describes the quantities that depend on the plasma equilibrium, that are required as input for the real-time suite of codes linked to RAPTOR.

This suite is broadly defined as a set of algorithms to perform real-time physics based estimation and prediction of plasma profiles. It may include RAPTOR and other observers such as neObs, and source calculation codes such as RABBIT and TORBEAM.

It is assumed that a real-time equlibrium code returns $\psi(R,Z)$ on a known $(R,Z)$ grid. We describe here the post-processing calculations that should be done based on this $\psi(R,Z)$ map. These can be done either as post-processing steps within the equilibrium code itself, or in separate codes/modules.

The outputs should all be given in the code's own COCOS convention \cite{Sauter2012}\footnote{Known COCOS values:, LIUQE (COCOS=17), JANET (COCOS=17), CHEASE (COCOS=2), EQUINOX (COCOS=7). E-mail the author with your COCOS number to add it to the list.}, and should be consistent with the experimental signs (e.g. the code's output $F=RB_\phi$ must flip sign with the experimental $B_{0}$, and $q$ must flip sign with both $B_0$ and $I_p$). 


\section{Flux surface integrals}
We propose a specific set of geometry profiles to be calculated in the equilibrium code, for specific values of normalized $\psi$. These are chosen to avoid infinite quantities at the edge or badly defined quantities at the axis.

One has to choose on which values of $\psi_N = (\psi-\psi_a)/(\psi_b - \psi_a)$, where $\psi_a$ and $\psi_b$ are the flux at the axis and boundary respectively, to compute the contour integrals and shape parameters. This set of psiN points on which the geometry quantities are evaluated should ideally be given as an input parameter vector to the code. A default choice is described below:

\subsection{Choice of default $\psi_N$ grid}

Near the axis, one is limited by numerical errors due to grid resolution in computing the contour integrals. Close to the boundary for diverted plasmas, however, one would like to have more contours since the value of some contour integrals rapidly diverges to zero. It is proposed to quadratically distribute the points as follows:

Starting from an equidistant grid $x \in [0,1]$ with number of points $N_p$, the values of $\psi_N$ for which the contours are calculated $\psi_N = 1-(1-x)^2$.
For 16 points, this gives
\be
\psi_N = 
\left[\begin{array}{c|c|c|c|c|c|c|c}
 0    & 0.1289   &  0.2489   & 0.3600  &  0.4622  &  0.5556  &  0.6400  &  0.7156  \\  \hline
 0.7822 &  0.8400  &  0.8889  &  0.9289  &  0.9500  &  0.9822   &  0.9956   & 1.0000
\end{array}\right]
\label{eqn:psiN}
\ee
Note that the point at 0.95 should be at 0.96 according to the formula, but has been moved in the interest of calculating $q_{95}$ at the same time as the other contour integrals.


\subsection{Flux-surface quantities required as output from the equilibrium code}
\label{sec:FSA} 
We define in Table~\ref{tab:fsa_inputs} the profiles to be calculated as outputs of the equilibrium code. We also specify their (analytical) values at the magnetic axis and for a diverted last closed flux surface.

\begin{table}[h!!]
\begin{tabular}{|p{2.5cm}|p{1.0cm}|p{1.5cm}|p{2.cm}|p{1.5cm}|p{3.5cm}|c|}\hline
Quantity & Suggested name & Value on axis &  Value on flux surfaces without $\times$-point & Value on separatrix with $\times$-point & Relation to RAPTOR internal variable (COCOS=11)  \\\hline
$\left\langle {1}/{R} \right\rangle $& \verb|iR| & $1/R_{a}$   & $C_0/C_1$ & $1/R_{\times}$ & -- \\
$\left\langle {1}/{R^2} \right\rangle$& \verb|iR2| & $1/R_{a}^2$ & $C_2/C_1$ & $1/R_{\times}^2$ & \verb|iR2|\\
$\left\langle B_p^2 \right\rangle$ & \verb|Bp2|& 0 & $C_3/C_1$ & 0 & \verb|dpsi2oR2| $= (4\pi^2)(\cdot)$ \\
$\left\langle (RB_p)^2 \right\rangle$  & \verb|R2Bp2|& 0 & $C_4/C_1$ & 0  & \verb|dpsi2| $= (4\pi^2)(\cdot)$ \\
$\left\langle RB_p \right\rangle$  & \verb|RBp|& 0 & $C_5/C_1$ & 0 & \verb|dpsi| $= (2\pi)(\cdot)$  \\
$1/(\oint \frac{\mathbf{B_p\!\cdot\! d\ell}}{B_p^2})$ & \verb|H|& $ \frac{F(0)}{2\pi R_{ax}^2 q_0}$ &  $1/C_1$ & 0 & \verb|dpsidV| $=H$  \\
$(1/q)$ & \verb|iota| & $1/q_0^\dagger$ &  $2\pi/(FC_2)$ & 0 & \verb|iota|\\ 
$F=RB_\phi$ & \verb|F| & $F(0)$ & $F(\psi)$ & $F(\psi_b)$ & \verb|F| \\ 
$\langle \frac{j_\phi}{R} \rangle$ & \verb|jphioR| & $p'+\frac{1}{\mu_0}\langle\frac{1}{R^2}\rangle FF'$ & $p'+\frac{1}{\mu_0}\langle\frac{1}{R^2}\rangle FF'$ & $p'+\frac{1}{\mu_0}\langle\frac{1}{R^2}\rangle FF'$  & - \\  \hline
\end{tabular}\caption{Geometric profiles to be returned by equilibrium code and expressions for correct evaluation including axis and LCFS.\\
$*$Note the relation
$2\pi/\oint \frac{d\ell}{B_p} = 2\pi H = F \left\langle \frac{1}{R^2} \right\rangle \frac{1}{q}$. Hence only three out of the four quantities $H$, $F$, $1/q$, $\langle1/R^2 \rangle$ are required, and the third can always be recovered.
$*$Note also that the surface area of the 3D flux surface (not to be confused with the poloidal cross-sectional area) is equal to $S_\psi = 2\pi C_5$.\\\
$^\dagger$Note that $q_0$ can be computed accurately from a local expansion of the equilibrium $\psi(R,Z)$ around the magnetic axis using  equation (88) in \cite{Moret2015}.
}
\label{tab:fsa_inputs}
\end{table}


\subsubsection{Contour integral definition with cosistent signs}
Quantities in Table~\ref{tab:fsa_inputs} are defined as ratios of contour integrals on the surfaces of constant $\psi$ on the poloidal cross-section. Following CHEASE \cite{Lutjens1996}, one defines the contour integrals $C_0 \ldots C_5$ of the form:

\be
C_i = \oint |\nabla \psi|^\mu R^\nu \dd \ell
\label{eqn:Ci}
\ee
where $\dd\ell$ is the element of the poloidal flux surface contour. To be independent of the definition of $\psi$ in each code, they are written here in terms of the poloidal field $B_p$, and rewritten in generic that works for any COCOS choice of a particular code. 

To be fully consistent, the contour integrals should be evaluated with the proper sign of the contour integration, with the line integration element $\dl$ in the direction of positive $\theta$ for the code's coordinate system. Alternatively, absolute values can be used in the calculations, but the results must be corrected depending on the sign of plasma current and the sign convention of the code (see \eqref{eqn:C0}--\eqref{eqn:C5}).

Using the COCOS convention we may write
\be
|B_p| = \frac{1}{(2\pi)^{e_{B_p}}} \frac{|\nabla \psi|}{R}
\quad \hbox{and} \quad 
\oint \Bp \cdot \dl = \sigma_{I_p}\sigma_{\rho\theta\phi}\oint |B_p| dl_p
\ee

This gives the following definitions of the contour integrals
\be
C_{0}= \oint \frac{\Bp\cdot \dl}{B_p^2R} = \sigma_{I_p}\sigma_{\rho\theta\phi}\oint \frac{1}{R} \frac{\dd l_p}{|B_p|} =  \sigma_{I_p}\sigma_{\rho\theta\phi}(2 \pi)^{e_{Bp}} \oint \frac{\dd l_p}{|\nabla \psi|}  \label{eqn:C0}
\ee
\be
C_{1}= \oint \frac{\Bp\cdot \dl}{B_p^2}  = \sigma_{I_p}\sigma_{\rho\theta\phi} \oint \frac{\dd l_p}{|B_p|} = \sigma_{I_p}\sigma_{\rho\theta\phi}(2 \pi)^{e_{Bp}} \oint \frac{R\dd l_p}{|\nabla \psi|} 
\ee
\be
C_{2}=  \oint \frac{\Bp\cdot \dl}{R^2B_p^2}  = \sigma_{I_p}\sigma_{\rho\theta\phi} \oint \frac{1}{R^2} \frac{\dd l_p}{|B_p|} = \sigma_{I_p}\sigma_{\rho\theta\phi}(2 \pi)^{e_{Bp}} \oint \frac{\dd l_p}{R|\nabla \psi|} 
\ee
\be
C_{3} = \oint  \Bp\cdot \dl  = \sigma_{\rho\theta\phi} \mu_0 I_p(\psi) = \frac{\sigma_{I_p}\sigma_{\rho\theta\phi}}{(2 \pi)^{e_{Bp}}} \oint
\frac{|\nabla \psi|}{R} \; \dd l_p
\ee
\be
C_{4}= \oint  R^2\Bp\cdot \dl  = \sigma_{I_p}\sigma_{\rho\theta\phi}\oint R^2 B_p^2 \frac{\dd l_p}{|B_p|} = \frac{\sigma_{I_p}\sigma_{\rho\theta\phi}}{(2 \pi)^{e_{Bp}}}
\oint R |\nabla \psi|  \; \dd l_p
\ee
\be
C_{5}= \oint R \; \dd  l_p \label{eqn:C5}
\ee

Here, $\sigma_{Ip}$ is the sign of the plasma current in the code's COCOS convention. $\sigma_{\rho\theta\phi}$, $\sigma_{Bp}$ and $e_{Bp}$ are the values corresponding to the code's COCOS choice \cite{Sauter2012}.
 
\subsubsection{Note on values at the axis or LCFS}
For closed flux surfaces excluding the axis or diverted last closed flux surface, the contour integrals can be computed by numerical integration around the contour. The accuracy of the contour integration is determined by the spatial grid resolution and numerical contour integration method, and may be compromised close to the axis (due to very few grid points) and close to a diverted LCFS (due to strong curvature of the flux surfaces close to the $\times$-point.

Consider cases with $\mu<0$. In these cases both $\dd \ell$ and $|\nabla \psi|^\mu$ go to infinity at the axis, and the result is ill-defined. At a LCFS with $\times$-point, $|\nabla\psi|^\mu$ goes to infinity and the integral goes to infinity. This is why many flux surface quantities need to be evaluated \emph{only} at closed, or limited flux surfaces to compute the expressions of Table~\ref{tab:fsa_inputs}. Values at the axis and diverted LCFS should never be computed by contour integration, nor should they be extrapolated. Instead, their exact values are prescribed in Table~\ref{tab:fsa_inputs}.

\subsubsection{Outline of numerical algorithm to compute $C_i$}
\begin{algorithm}[H]
% \KwData{this text}
% \KwResult{how to write algorithm with \LaTeX2e }
\textsc{ Initialization}: load $({\psi_{N,grid}},R_{grid},Z_{grid})$\;
\textsc{ Step}:
 Get data for time $k$: $(\Psi_N),R_{ax},Z_{ax}$ \;
 \For{Value $i$ on $\psi_{N,grid}$}{
  $\psi_{N,i} = \psi_{N,grid}[i]$\;
  \uIf{$\psi_{N,i}=0$ OR $(\psi_{N,i}=1$ AND PlasmaIsDiverted)}
  {Do not compute contour integrals, treat later using special rules.}
\uElse(Closed flux surface){
  $(R_c,Z_c)$ = ContourFinder$(\Psi_N,R_{grid},Z_{grid},\psi_{N,grid},R_{ax},Z_{ax})$ \footnote{CountourFinder() should use known axis position to find good starting points for each contour}\;
  $L_c$ = ComputeSegmentLength$(R_c,Z_c)$\;
  $|\nabla \psi|_c$ = InterpolatingDerivative$(\Psi_N,R_{grid},Z_{grid},R_c,Z_c)$\;
  \For{$j=0:5$}{
  $G_c = |\nabla \psi|_c^\mu R_c^{\nu}$ \quad \hbox{(Integrand from eq\eqref{eqn:Ci})}\;
  $C_j$ = TrapezoidalInterpolation$(L_c,G_c)$
   }
   }
 }
 \caption{Outline of algorithm for computing contour integrals $C_0,\ldots,C_5$}
\end{algorithm}


\section{Shape descriptors}
Some codes require, e.g. for neoclassical calculations and for some advanced transport models, profiles of the shaping information for each flux surface. Expressions for these quantities are given in Table~\ref{tab:shape_profiles}. The values should be returned for each flux surface on the same $\psi_N$ grid as used in Section~\ref{sec:FSA}. 
\begin{table}[h!!]
\begin{tabular}{|c|p{1cm}|c|p{2cm}|}\hline
Quantity & Value on axis &  Flux surfaces with or without $\times$-point & suggested variable name \\\hline
$R_{geom}$ & $R_{ax}$  & $(R_{max}+R_{min})/2$ & \verb|Rgeom|\\ \hline
$\epsilon$ & 0 &  $a / R_{geom} = (R_{max}-R_{min})/(R_{max}+R_{min})$  & \verb|epsilon|\\ \hline
$\kappa $ & 0 & $(Z_{upper} - Z_{lower})/(R_{max}-R_{min})$  & \verb|kappa| \\ \hline
$\delta $ & 0 & $(2R_{geom}-R_{upper} - R_{lower})/(R_{max}-R_{min})$ & \verb|delta|  \\ \hline
\end{tabular}
\caption{Flux surface shape quantities}
\label{tab:shape_profiles}
\end{table}

\section{Surface/Volume integrated quantities}
In addition to the flux quantities of Section~\ref{sec:FSA}, some integral quantities are needed. These should be returned on the same $\psi_N$ grid as the flux quantities of Section~\ref{sec:FSA}.

Depending on the equilibrium reconstruction code, it may or may not be easily possible to perform 2D surface integration over flux surfaces. If accurate integration is possible, then this is preferable, otherwise,  contour integrals can be integrated over $\psi$. Both options are listed in Table~\ref{tab:integral_quantities}.

\begin{table}[h!!]
\begin{tabular}{|p{1.6cm}|p{0.5cm}|p{1cm}|p{1cm}|p{3.5cm}|p{3cm}|}\hline
Quantity & unit &Suggested name &Value on axis & Value as surface integrals in (R,Z) plane&  Value as integrals over $\psi$  \\\hline
$V(\psi)$ (Volume ) &[m$^3$] & \verb|Vol| & 0 & $2\pi\int^{\Omega_p} R dS_{\phi}$ &$2\pi c_* \int_{\psi_a}^{\psi} C_1 d\psi \;(\star)$  \\
$S(\psi)$ (Pol. area) & [m$^2$] & \verb|Area| & 0 & $\int^{\Omega_p} dS_{\phi}$  & $c_* \int_{\psi_a}^{\psi} C_0 d\psi\;(\star)$  \\
$\Phi(\psi)$ (Tor. flux) & [Wb] & \verb|Phi| & 0 &  $\int^{\Omega_p} \frac{F}{R} dS_{\phi}$  &$c_* \int_{\psi_a}^{\psi} q d\psi\;(\star)$  \\
$\rho_{tor_N}(\psi)$ & [-] & \verb|rhotor| & 0 & - &$\sqrt{\Phi(\psi)/\Phi(\psi_b)}$ \\
$\mu_0I_p(\psi)$ (current) & [A] & \verb|mu0Ip| & 0 &  $\mu_0\int^{\psi_N<1} j_\phi dS_{\phi}$  & $\sigma_{\rho\theta\phi}C_3$  \\
\hline
\end{tabular}\caption{Quantities related to surface integrals to be computed by equilibrium code. Beware that integrands marked with $(\star)$ diverge (are infinity at the LCFS), so trapezoidal integration is discouraged (see below).
$c_*=2\pi^{-e_{B_p}}\sigma_{\rho\theta\phi}\sigma_{B_p}$ }
\label{tab:integral_quantities}
\end{table}

\paragraph{Alternative integration scheme for diverging LCFS quantities}
In case one is integrating over $\psi$, some profiles take infinite value on the LCFS in case of an $\times$-point (as discussed previously). Integrating these quantities numerically up to the edge can then give inaccurate results. A typical example is the integration of $q$ to obtain the toroidal flux. Using COCOS=11 sign conventions for brevity:
 \be	 
 \Phi(\psi) = \int_{\psi_a}^{\psi} q(\psi) \dd \psi	 
 \label{eqn:Phi}	 
 \ee	 
 If $q$ is equal to infinity for diverted LCFS, and grows to large values at $\psi_N$ close to 1, it is easy to make an error in the value of $\Phi_b = \Phi(\psi_b)$, which is needed for the definition of $\rho_{tor}$, for example. A strategy to compute $\Phi(\psi)$ more accurately is as follows:	 
 	 
 \begin{itemize}	 
 \item Given $1/q$ on some grid of known $\psi$ values (with grid index $i$), and recalling $\frac{1}{q} = \dede{\psi}{\Phi}$
 \item Write $\int_{\Phi_i}^{\Phi_{i+1}} 1/q \; \dd \Phi = \psi_{i+1} - \psi_i$	 
 \item Formulate a numerical approximation to this integral. For trapeze integration, we get:	 
 \be	 
 \frac{1}{2} \left(\frac{1}{q_i}+\frac{1}{q_{i+1}}\right)(\Phi_{i+1}-\Phi_i) = (\psi_{i+1} - \psi_i )	 
 \ee	 
 yielding a recursive equation for $\Phi_{i+1}$	 
 \be	 
 \Phi_{i+1} = \Phi_i + 2\left(\frac{1}{q_i}+\frac{1}{q_{i+1}}\right)^{-1} (\psi_{i+1} - \psi_i )	 
 \ee	 
 Note that the great advantage of this scheme is that $1/q_{N} = 0$ for a diverted LCFS, thus entirely avoiding the singularity.
 \end{itemize}	

\section{Scalar quantities}
Finally, the following scalar quantities are desired outputs of the code.
See Table \ref{tab:scalar_quantities}
\begin{table}[h!!]
\begin{tabular}{|p{2cm}|p{1.5cm}|p{1cm}|p{3.5cm}|p{3cm}|}\hline
Quantity & unit &Sugg. name & Value as surface integrals in (R,Z) plane&  Value as integrals over $\psi$  \\\hline
$W_k$ (Kin. energy) & [J] & \verb|Wk| &  $2\pi\frac{3}{2}\int^{\Omega_p}  p(\psi) RdS_{\phi}$ & $2\pi c_* \int_{\psi_a}^{\psi_b} p C_1 d\psi \;(\star)$   \\
$W_{pol}$ (Pol. magnetic energy) & [J] & \verb|Wpol|  &  $2\pi\int^{\Omega_p} \frac{ B_p^2}{2\mu_0} RdS_\phi$ & $\pi c_* \sigma_{I_p} \int_{\psi_a}^{\psi_b} I_p(\psi) d\psi$  \\ 
$W_{tor}$ (Plasma contribution to toroidal magnetic energy) & [J] & \verb|Wtor| &  $2\pi\int^{\Omega_p} \frac{ B_\phi^2-B^2_{\phi,vac}}{2\mu_0} RdS_\phi$ & $\frac{2\pi c_*}{2\mu_0}   \int_{\psi_a}^{\psi_b} (F^2-R_0^2B_0^2) C_2 d\psi \;(\star)$  \\ 
$\Phi_{p}$ (Plasma contribution to toroidal flux) & [Wb] & \verb|Phip| &  $\int^{\Omega_p}  (B_\phi-B_{\phi,vac}) dS_\phi$ & $\Phi_b -c_* R_0B_0  \int_{\psi_a}^{\psi_b} C_2 d\psi \;(\star)$  \\ 
$\psi_{a}$  & COCOS & \verb|psia| & &  \\
$\psi_{b}$  & COCOS & \verb|psib| & &  \\
$R_{ax}$  & [m] & \verb|Rax| & &  \\
$Z_{ax}$  & [m] & \verb|Zax| & &  \\
$I_p$  & [A] & \verb|Iptot| & &  \\
$B_0$  & [T] & \verb|B0| & &  \\
\hline
\end{tabular}
\caption{Scalar quantities to be computed by equilibrium code. `COCOS' units indicates that the unit depends on the chosen COCOS convention. $(\star)$ indicates again that the integrand diverges at a diverted LCFS.}
\label{tab:scalar_quantities}
\end{table}

\newpage
\section{Calculations for diagnostics interpretation}
These only apply when the equilibrium code is to be interfaced to RAPTOR-observer or another code that uses real-time kinetic diagnostic data. The equilibrium code should then return some extra information on the mapping between equilibrium and a diagnostic. Depending on the preference of the equilibrium code developer, this is done either as part of the equilibrium code itself, or as a post-processing step in a separate module that receives the equilibrium map information from the equilibrium code.

\subsection{Point-wise diagnostics}

Inputs:
\begin{itemize}
\item Two vectors of $(R,Z)$ [m] coordinates of point-wise diagnostic measurements (e.g. Thomson Scattering)
\end{itemize}
\noindent
Outputs:
\begin{itemize}
\item $\psi_N$ and $\rho_{tor}$ values inside the plasma corresponding to these (R,Z) points, with a $-1$ value indicating the the point is outside the plasma region.
\end{itemize}


\subsubsection{Outline of fast algorithm for evaluation of $\Psi_N$ and $\rho_{tor}$ on fixed $R,Z$ points}
\begin{algorithm}[H]
 \textsc{Initialization:} Given a fixed 2D grid $(R_{grid},Z_{grid})$ and fixed evaluation points $(R_c,Z_c)$:\\
$W_c$ = InterpolationWeights$(R_{grid},Z_{grid},R_c,Z_c)$\footnote{Compute weight matrix $W_c$ that interpolates quantities from $(R,Z)$ grid points onto the points $R_c,Z_c$, such that $\bar{\Psi}_{N,c} = W_c\bar{\Psi}_N$ where $\bar{\Psi}_{N,c}$ is the vector of $\Psi$ values on $R_c,Z_c$ and $\bar{\Psi}_N$ is the vector of $\Psi_N$ values on the $(R_{grid},Z_{grid})$ grid.}\\
\textsc{Step}: {Get data for time $k$: ${\bar\Psi}_N$, $\psi_{N,grid}$, $\rho_{tor,N}$ on $\psi_{N,grid}$ (from Table~\ref{tab:integral_quantities}) \\
$\bar{\Psi}_{N,c} = W_j\bar{\Psi}_N$ \quad \hbox{(Compute $\Psi_N$ values on points along chord $j$)} \;
$\rho_{tor,N,c}$ = Interpolate1D($\psi_{N,grid},\rho_{tor,N},\bar{\Psi}_{N,c}$)\\
}
 \caption{Outline of algorithm for computing integrals of basis functions on diagnostic chords.}
 \label{alg:fixRZ}
\end{algorithm}


\subsection{Line-integrated diagnostics}

Inputs:
\begin{itemize}
\item Vector of $n_l$ launch points $(R,Z)$ (in [m]) and orientation vector $[k_r,k_z,k_\phi]$ (such that $\mathbf{k} = k_r\mathbf{e_r} + k_z\mathbf{e_z} + k_\phi\mathbf{e_\phi}$ is the direction of the chord towards the plasma, of any line integrating diagnostics.
\item $m \times n_b$ matrix containing: Values of $m$ basis functions $b_i(\bar{\psi}_N)$ for $i=\{1,\ldots,m\}$, evaluated on $\bar{\psi}_N$ grid (vector of $n$ points). This may be a different grid than the input grid $\psi_N$ defined in \eqref{eqn:psiN}.
\item $n_b$-sized vector containing $\bar{\psi}_N$ grid on which basis functions (above) are evaluated. above.
\item Alternatively to the two points above, analytical expressions for the $m$ basis functions $b_i(\psi_N)$ can be prescribed, then only $m$ should be an input parameter to the code. For this option, B-splines are recommended due to the fact that they have a local support.
\end{itemize}
\noindent
Outputs:
\begin{itemize}
\item Matrix ($m\times n_l$) of line integrals over the chord for each basis function, each entry $i,j$ being 
\be 
  M_{ij}=\int b_i(\psi_N(R,Z)) ds_j
\ee
With $ds_j$ the line element along the $j$th diagnostic chord. In case of double-pass chords, the full path length should be considered.
\end{itemize}

\subsubsection{Outline of algorithm for line integration of basis functions a on a fixed chord}

\begin{algorithm}[H]
% \KwData{this text}
% \KwResult{how to write algorithm with \LaTeX2e }
 \textsc{Initialization:}\;
 Given a fixed 2D grid $(R_{grid},Z_{grid})$ grid and fixed launch point and launch direction:\\
\For{Each chord $j$}{
 Compute a set of points $(R_j,Z_j)$ that lies along the chord.\\
 Compute the path lengths $L_j$ for each line segment.\\
$W_j$ = InterpolationWeights$(R_{grid},Z_{grid},R_j,Z_j)$ (Similar to Algorithm~\ref{alg:fixRZ})\\
}
\textsc{Step}: \\
Get data for time $k$: ${\bar\Psi}_N$ (2D psi map)\\
\For{Each chord $j$}{
$\bar{\Psi}_{N,j} = W_j\bar{\Psi}_N$ \quad {(Compute $\Psi_N$ values on points along chord $j$)} \;
\For{Each basis function $i$}
{
$V_{i,j} = b_i({\bar\Psi}_{N,j})$, {(Value of basis functions on each chord point)}\\
$M_{ij} =$ TrapezoidalIntegration$(L_j,b_{i,j})$ \quad (Line integral of basis function on chord)\\
}
 }
 \caption{Outline of algorithm for computing integrals of basis functions on diagnostic chords.}
\end{algorithm}


\subsection{ECE (cold approximation, X mode)}
Inputs:
\begin{itemize}
\item Vector ($n_{ECE}\times 1$) with frequency of each ECE channel [GHz].
\item A vector of integers ($n_{ECE}\times 1$) indicating the harmonic used (1,2,3..)
\item Two matrices ($n_{ECE}\times 3$) of launch points and receiver vectors of the antenna (pointing into the vessel) (see definition above).
\item Plasma electron density values on a $\psi_N$ grid [$\meter^{-3}$]
\item $psiN$ grid used to specify electron density profile.
\end{itemize}
\noindent
Outputs:
\begin{itemize}
\item Vector of $n_{ECE}$ with values of $\psi_N$ where $B_{tot}(R,Z) = \sqrt{B_\phi^2 + B_p^2}$ matches for the first time the ECE cold resonance condition $f[\hertz]=q_eB/2\pi m_e$, along the antenna path.
\item Vector of booleans indicating whether an ECE channel has a resonance inside the plasma (1=yes, 0=no).
\item Vector of $n_{ECE}$ values $f_{cutoff,max}$ indicating the maximum fraction of the ECE cutoff density  encountered during the path between the launcher and the cold absorption point. 

For X-mode, there is only an upper density limit: $n_{cutoff,X} = l(l-1)B^2 \times 9.7197\times 10^{18}$
hence 
\be
f_{cutoff,max} = \max_{s}(f_{cutoff}(s)) = \max_{s}(n_e(s)/n_{cutoff,X}(s))
\ee 
where $s$ is the distance along the antenna launch vector, evaluated between the launch point and the cold resonance point, and \be
n_{cutoff,X} = l(l-1)B^2 \times 9.7197\times 10^{18}
\ee 
Here, $l$ is the harmonic number and $B$ is the total field in [T].

\end{itemize}

\section*{Changelist}
\begin{itemize}
\item v2.2: Remove $2\pi$ from H definition
\item v2.1: Add algorithm outlines
\item v2.0: Complete rewrite with details of diagnostics interfaces, integral quantities etc (for JET, July 2018)
\end{itemize}

\bibliography{/Users/ffelici/Documents/SPC/References/Library.bib}


  \end{document}


