APPENDIX.tex

% Similar equations are applied for all the other meshes of \textit{ElastoDyn} and \textit{BeamDyn} module.
% %     \text{BD}\%u\%F = \text{BD}\%y\%F -  \epsilon_C\  \text{BD}\%y\%v
% \question{We need to make sure we do not damp some of the wanted motions (rotor's speed).}
% 
% in \tt{ED\_InputSolve}
% \begin{lstlisting}
% !p_FAST%CompAero = Module_None
% DO K = 1,SIZE(u_ED%BladePtLoads,1) ! Loop through all blades (p_ED%NumBl)
%   CALL Transfer_Line2_to_Point( y_AD%BladeLoad(k), u_ED%BladePtLoads(k), MeshMapData%AD_L_2_BDED_B(k), ErrStat2, ErrMsg2, u_AD%BladeMotion(k), y_ED%BladeLn2Mesh(k) )
%   ! OR
%   u_ED%BladePtLoads(K)%Force  = 0.0_ReKi
%   u_ED%BladePtLoads(K)%Moment = 0.0_ReKi
% END DO         
% \end{lstlisting}
% 
% % \begin{lstlisting}
% % ! connect AD mesh with BeamDyn
% % DO K=1,NumBl         
% % CALL MeshMapCreate( BD%y(k)%BldMotion, AD%Input(1)%BladeMotion(K), MeshMapData%BDED_L_2_AD_L_B(K), ErrStat2, ErrMsg2 )
% %    CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName//':BD_2_AD_BladeMotion('//TRIM(Num2LStr(K))//')' )
% % CALL MeshMapCreate( AD%y%BladeLoad(K), BD%Input(1,k)%DistrLoad,  MeshMapData%AD_L_2_BDED_B(K), ErrStat2, ErrMsg2 )
% %    CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName//':AD_2_BD_BladeLoad('//TRIM(Num2LStr(K))//')' )
% % END DO
% % \end{lstlisting}
% % 
% 
% 
% ElastoDyn Inputs Meshes
% \begin{lstlisting}
% BladePtLoads{:} "A mesh on each blade, containing aero loads"
% PlatformPtMesh  "A mesh at the platform reference (point Z), containing loads"
% TowerPtLoads    "Tower line2 mesh with tower loads"
% HubPtLoad       "A mesh at the teeter pin, containing loads on the hub"
% NacelleLoads    "From ServoDyn/TMD: loads on the nacelle."
% \end{lstlisting}
% 
% ElastoDyn Outputs Meshes:
% \begin{lstlisting}
% BladeLn2Mesh	{:}		"A mesh on each blade, kinematics of blade nodes"
% PlatformPtMesh	-		"Platform reference point kinematics"
% TowerLn2Mesh	-		"Tower line2 mesh kinematics"
% HubPtMotion14	-		"For AeroDyn v14: motions of the hub"
% HubPtMotion	-	-	"For AeroDyn and Lidar(InflowWind): motions of the hub"
% BladeRootMotion14		-	"For AeroDyn v14: kinematics the blade roots"
% BladeRootMotion	{:}		"For AeroDyn/BeamDyn:  kinematics the blade roots"
% RotorFurlMotion14		-	"For AeroDyn14: motions of the rotor furl point."
% NacelleMotion	-		"For AeroDyn14 & ServoDyn/TMD: motions of the nacelle."
% TowerBaseMotion14		-	"For AeroDyn 14: motions of the tower base"
% \end{lstlisting}
% 
% 
% 
% 
% 
% BeamDyn Outputs Meshes:
% \begin{lstlisting}
% RootMotion  "contains motion" -
% PointLoad   "Applied point forces along beam axis" -
% DistrLoad   "Applied distributed forces along beam axis" -
% HubMotion   "motion (orientation) at the hub" -
% \end{lstlisting}
% 
% BeamDyn Outputs Meshes:
% \begin{lstlisting}
% ReactionForce  "contains force and moments" -
% BldMotion      "Motion (disp,rot,vel, acc) along beam axis" -
% \end{lstlisting}
% 




% --------------------------------------------------------------------------------}
% ---  
% --------------------------------------------------------------------------------{



% \paragraph{Step 1 - Store outputs for one rotor revolution at \tt{NAzimStep} azimuthal positions}
% \begin{itemize}\tightlist
% 
%     \item Set the vector of azimuthal positions $\v{\psi}_0$ of dimension $1\times \tt{NAzimStep}$
%         \begin{align}
%              j=1..\tt{NAzimStep}
%              ,\qquad
%          \v{\psi}_0[j] = \psi_\text{init} +  (j-1) \Delta\psi + \Delta\psi/2, \qquad \Delta\psi\eqdef \frac{2\pi}{\tt{NAzimStep}}
%         \end{align}
%     \item Allocate the matrix $\m{Y}_0$ of dimension $n_{y}\times \tt{NAzimStep}$ to store the output vector from all the modules at each target azimuthal position $\v{\psi}_0$, where $n_{y}$ is the number of outputs from all the modules
% 
%     \item Initialize the following variables:
%         \begin{align}
%          j = 1
%          , \qquad
%          \psi_p = \psi_\text{init}
%          , \qquad
%          \v{y}_p =\operatorname{NaN}(1,n_y)
%               %\label{eq:}
%         \end{align}
%       The subscript $p$ refers to the previous time step, the subscript $c$ will be used for the current time step.
%   \item Perform time stepping for one rotor revolution. The output vector is stored for each target azimuthal position $\v{\psi}_0$. For each time step:
%         \begin{itemize}\tightlist
%             \item Call \tt{FAST\_Solution\_T}
%             \item If $t=0$, store the current rotor speed as the reference one: $\Omega_0$
%             \item Store the current azimuthal angle $\psi_{c}$ 
%             \item Store the current output vector from all the modules $\v{y}_\text{c}$
%             \item If $(\psi_c -\psi_p)>\Delta\psi/2$, return an error: the rotor is spinning too fast, the time step is too large
%             \item If $\psi_c \geq \psi_0[j]$, interpolate the output vector to the target azimuthal position $\psi_0[j]$ using the current output values and the previous ones:
%                 \begin{align}
%                     \v{y}_0 = \v{y}_p + (\v{y}_c-\v{y}_p)\frac{\psi_0[j]-\psi_p}{\psi_c-\psi_p}  %\label{eq:}
%                 \end{align}
%                 Store the interpolated value
%                 \begin{align}
%                     \m{Y}_0[:,j] = \v{y}_0
%                 \end{align}
%                 Increment $j$ %, $j\leftarrow j+1$
% %             \begin{align}
% %                 \v{y}_0[n] = \v{y}(\psi_0[t_n])
% %                   %\label{eq:}
% %             \end{align}
%             \item Set the current values as previous values for the next time step
%                 \begin{align}
%                    \psi_p \leftarrow \psi_c
%                    ,\qquad
%                    \v{y}_p \leftarrow \v{y}_c
%                 \end{align}
%             \item if $j>\tt{NAzimStep}$ exit the loop
%     \end{itemize}
% \end{itemize}
% 
% \paragraph{Step 2 - iterate to a constant rotational speed}
% \begin{itemize}
%     \item From this point, the glue code will enforce a controller offset $p_\text{off}$ to the trimmed controller parameter, and an external damping force will be added to help the convergence of the system.
% \item Initialize the controller offset $p_\text{off}$ to zero.
%     \item Start an iteration loop
%         \begin{itemize} \tightlist
%             \item Compute a reference value for each of the index of the output vector, based on the maximum and minimum values taken over one rotor revolution.
%                 \begin{align}
%                     y_{\text{ref}}[i] = \max(Y_0[:,i]) - \min(Y_0[:,i])
%                 \quad\text{if }\
%                 y_\text{ref}[i]>10^{-6}
%                 ,\quad\text{else }\
%                 y_\text{ref}[i]=1
%                 \end{align}
% 
%             \item Initialize the following variables:
%                 \begin{align}
%                  j = 1
%                  , \qquad
%                  \psi_p = \psi_\text{init}
%                  , \qquad
%                  \v{y}_p =\operatorname{NaN}(1,n_y)
%                       %\label{eq:}
%                 \end{align}
% 
%             \item Perform time stepping for one rotor revolution, where for each time step:
%                 \begin{itemize} \tightlist
%                     \item Call \tt{FAST\_Solution\_T} (Note: the damping is increased via the variable $\epsilon_C$ and the trimmed controller value will receive the offset $p_\text{off}$)
%                     \item Store $\Omega$, $\psi_c$, $\v{y}_c$
% %                     \item If $(\psi_c -\psi_p)>\Delta\psi/2$, return an error: the rotor is spinning too fast, the time step is too large
%                     \item If $\psi_c \geq \psi_0[j]$, interpolate the output vector to the target azimuthal position $\psi_0[j]$ using the current output values and the previous ones:
%                         \begin{align}
%                             \v{y}_0 = \v{y}_p + (\v{y}_c-\v{y}_p)\frac{\psi_0[j]-\psi_p}{\psi_c-\psi_p}  %\label{eq:}
%                         \end{align}
% 
% %                     \item Interpolate the previous output vector $y_0$ to the current value of $\psi$: $\v{y}_{p}=\operatorname{interp}(\psi_0,\v{y}_0,\psi)$ 
%                         Compute the L2-norm of the relative error of the output vector at the azimuthal position $\psi_0[j]$ between the current revolution and the previous one:
%                 \begin{align}
%                     \epsilon_y^2[j] %=\rVert\boldsymbol{\epsilon_y}\lVert^2
%                     = \sum_i \left(\frac{y_0[i]-Y_0[j,i]}{y_\text{ref}[i]}  \right)^2
% %                     ,\qquad \epsilon_{\dot{q}}^2 = \rVert (\boldsymbol{\dot{q}}-\boldsymbol{\dot{q}_0}). \boldsymbol{w}^t\lVert^2
%                 \end{align}
%                 Note: if the variable $y_0[i]$ is in radian or degree, the difference of the variable should be taken modulo $2\pi$. \\
% %             \item 
% %             \item Update damping values (\(C\)), using a stiffness proportional
% %                 damping scaled by the ,
% %                 \[ C_\text{new} = C_\text{start} + k_C \epsilon_{\dot{q}}^2 K\]
% %             \end{itemize}
% 
%             Store the interpolated value
%             \begin{align}
%                 \m{Y}_0[:,j] = \v{y}_0
%             \end{align}
%             Increment $j$%, $j\leftarrow j+1$
% 
%         \item     Compute the error in rotor speed \(\epsilon_\Omega\): the difference between the target speed and the current rotor speed
%                 \begin{align}
%                      \epsilon_\Omega ={\Omega - \Omega_0}
%                 \end{align}
%                 This error is used at each time step to correct the controller parameter
% %             \item 
% %                 The offset on the controller variable is incremented using
% %                 the speed error and a proportional gain $k_p$.
% %                 \begin{align}
% %                     p_\text{off} = p_\text{off} + k_p \epsilon_\Omega %\label{eq:}
% %                 \end{align}
% %                 This offset will be passed to the routine \tt{SD2ED} on the next time step
% 
% 
%         \item Set the current values as previous values for the next time step
%             \begin{align}
%                \psi_p \leftarrow \psi_c
%                ,\qquad
%                \v{y}_p \leftarrow \v{y}_c
%             \end{align}
% 
%         \item If $j>\tt{NAzimStep}$:  \\
%             Exit the iteration loop if $\epsilon^2[j]<\tt{TrimTol}$ for all $j$ (i.e. the output vector has converged for all azimuthal positions). Otherwise, set $j=1$ and continue the time stepping.
% 
% %             \item If $\epsilon_y^2$  is below a tolerance value \tt{Tol} exit the iteration loop
%         \end{itemize}
% %     \item Replace $\v{y}_0$ with the values of $\v{y}_1$, and $\psi_0$ with $\psi_1$: 
% %         \begin{align}
% %         \v{y}_0\leftarrow \v{y}_1
% %         ,\qquad
% %         \psi_0\leftarrow \psi_1
%               %\label{eq:}
% %         \end{align}
%     \item Continue to a next iteration unless the maximum number of iteration is exceeded
% \end{itemize}
% \end{itemize}
% 
% \paragraph{Step 2 - perform linearization}
% \begin{itemize}
% \item The values \tt{Twr\_Kdmp} and \tt{Bld\_Kdmp} are set to 0
% \item The controller off set is kept to its value $p_\text{off}$
% \item The linearization takes place
% \end{itemize}
% 
% 
% 
% 



\newpage
\appendix

\section{Input/output workflow of the controller parameters when linearization is on}
\begin{itemize}
    \item Yaw: ServoDyn returns the Yaw moment \tt{y\%YawMom}. 
        When linearization is enabled, \tt{YCMode=0}, and the routine \tt{CalcStandardYaw} returns the parameters \tt{p\%YawNeut}, hence \tt{YawPosCom}=\tt{YawNeut}. 
        The current nacelle yaw is provided as input (\tt{u\%Yaw}). It is first initialized to \tt{YawNeut}. It is used to set the operating point and in the moment calculation:
        \begin{align}
            \tt{YawMom} = k_\text{yaw} ( \tt{u\%Yaw} - \tt{YawPosCom} )
        \end{align}
        The input \tt{u\%YawAngle} is not used.
        The operating point \tt{y\_op} is set using \tt{y\%YawMom}.
        The operating point \tt{u\_op} is set using \tt{u\%Yaw}.
    \item Pitch: ServoDyn returns the pitch command \tt{y\%BlPitchCom}.
        When linearization is on, \tt{PCMode=0}, and then \tt{y\%BlPitchCom} is first set to \tt{p\%BlPitchInit}, but then the pitch maneuver overrides it. 
        The pitch maneuver is initialized with \tt{u\%BlPitch}.
        Apart from this initialization the input \tt{u\%BlPitch} is not used.
        The operating point \tt{y\_op} is set using the pitch command \tt{y\%BlPitchCom}.
    \item Generator Torque: ServoDyn returns \tt{y\%GenTrq}. 
        When linearization is on, \tt{VSCtrl} is 0 or 1, and the \tt{GenModel} is 1 or 2. 
        A non-zero generator torque is returned (based e.g. on the HSS speed, or the controller region) as long as the generator is online.
        The operating point is set using \tt{y\%GenTrq}.
\end{itemize}





\section{Workflow in FAST7}\label{old-workflow}
The computation is done in \tt{Fast\_lin.f90} within the routine \tt{CalcSteady}.
The workflow is as follows:

\begin{itemize}\tightlist
    \item Assign weights \(\boldsymbol{w}\) (\tt{QWeight},
        \tt{QDWeight}) for each DOF to ensure that each of them have a
        comparable contribution to the norm. For now, same weights for
        positions and speeds. \textbf{New: evaluate this on first iteration
        based on range of values - (max-min)}
    \item Store damping values\\
        \textbf{New: get the part of the system that has damping?}
    \item Save wind speed, as a safety measure to detect if windspeed is unsteady (in that case abort)
    \item Iteration loop
        \begin{itemize} \tightlist
            \item Save DOF position and velocities: \(\boldsymbol{q}\)
                (\tt{QStart}), \(\boldsymbol{\dot{q}}\) (\tt{QDStart})
            \item Start time stepping for one rotor revolution \textbf{Note: currently
                    \tt{NStep} is a constant, so the number of time steps is always
                    the same, and the loop is not based on azimuthal position. Probably
                fair and convenient also when \(\Omega=0\)}

                \begin{itemize} \tightlist
                    \item Call \tt{Solver}
                    \item Increase time step
                \end{itemize}
            \item Compute L2-norm in positions and speeds, \tt{AbsQNorm} and
                \tt{AbsQDNorm}, scaled by weights:
                \[\epsilon_q^2 = \rVert (\boldsymbol{q}-\boldsymbol{q_0}). \boldsymbol{w}^t\lVert^2 = \sum_i \left[(q(i)-q_0(i)) w(i) \right]^2
                    ,\qquad \epsilon_{\dot{q}}^2 = \rVert (\boldsymbol{\dot{q}}-\boldsymbol{\dot{q}_0}). \boldsymbol{w}^t\lVert^2
                \] \textbf{Note: special handling required for azimuth due to
                2pi-modulo}
            \item Compute the error in rotor speed \(\epsilon_\Omega\)
                (\tt{SpeedErr}): the difference between the target speed and the
                average rotor speed over the simulated period (computed during the
                loop above)
                \begin{align}
                    \epsilon_\Omega = \frac{\overline{\Omega} - \Omega_0}{\Omega_0}
                \end{align}
            \item Based on the \tt{TrimCase}, adjust the yaw DOF (either
                \tt{YawNeut} or the proper DOF in \tt{Q}) or the generator
                torque \tt{GenTrq} or the blade pitch \tt{BlPitch}, using
                the speed error and a proportional gain:
                \[ p_\text{new} = p_\text{old} + k_p \epsilon_\Omega\]
            \item Update damping values (\(C\)), using a stiffness proportional
                damping scaled by the ,
                \[ C_\text{new} = C_\text{start} + k_C \epsilon_{\dot{q}}^2 K\]
        \end{itemize}
\end{itemize}

\subsection{FAST 7 inputs}
\begin{itemize}\tightlist
    \item \tt{CalcStdy} This flag determines whether a periodic steady state solution is computed before linearizing the
model. If False, the next three inputs are ignored and the linearization occurs about the initial
conditions specified in FAST’s primary input file. That is, when CalcStdy is False, the operating
point is set to the condition in which all displacements, velocities, and accelerations are zero,
except those specified with nonzero initial conditions (for instance, the azimuth DOF will
increment at a constant rate if and when the rotor is spinning). If CalcStdy is True and
RotSpeed is nonzero, FAST integrates in time until a periodic steady state solution is reached.
The method of solution is determined by the next input, TrimCase. FAST is then linearized
about this periodic operating point. If CalcStdy is True and RotSpeed is zero, FAST will
disable GenDOF (if previously enabled) and integrate in time until a static equilibrium position is
found. FAST is then linearized about this position. The accuracy of the steady state solution is
determined through input convergence tolerances DispTol and VelTol (see below). This input is
not used in the FAST-to-ADAMS preprocessor. (flag)

\item \tt{TrimCase} This switch determines, for a variable speed machine, which control input to trim in order to reach
the desired azimuth-averaged rotor speed indicated through input RotSpeed (which is also the
initial rotor speed). 
 Setting it to 1 causes FAST to trim nacelle yaw command (demand) angle,
while maintaining constant rotor collective blade pitch (indicated by inputs BlPitch i ), to reach the
desired azimuth-averaged rotor speed. With yaw DOF enabled (YawDOF = True), the nacelle
yaw command is the neutral yaw angle, YawNeut, which is passed through FAST’s built-in,
second-order actuator model. With yaw DOF disabled (YawDOF = False), the nacelle yaw
command is the actual nacelle yaw angle. 

Setting TrimCase to 2 causes FAST to trim electrical
generator torque, while maintaining constant rotor collective blade pitch (indicated by inputs
BlPitch i ), to reach the desired azimuth-averaged rotor speed (i.e., Region 2 trim). 

Setting TrimCase to 3 causes FAST to trim rotor collective blade pitch to reach the desired azimuth-
averaged rotor speed (i.e., Region 3 trim). In this case, the initial “guess” blade pitch angles are
given by BlPitch i and the electrical generator torque is determined by the torque-speed
relationship indicated by inputs VSContrl or GenModel. For typical Region 3 trim, collective
pitch can be trimmed while maintaining a constant generator torque by setting TrimCase to 3,
VSContrl to 1, VS\_RtTq to the desired constant generator torque, and VS\_RtGnSp,
VS\_Rgn2K, and VS\_SlPc to 9999.9E-9 (very small don’t cares > 0.0). Input parameter
TrimCase is ignored when either CalcStdy or GenDOF is False. For a constant speed machine,
GenDOF should be set to False when linearizing FAST, in which case, input TrimCase is
ignored. Using values other than 1, 2, or 3 will cause FAST to abort. This input is not used in the
FAST-to-ADAMS preprocessor. (switch)


\item \tt{DispTol} This is the convergence tolerance for the 2-norm of angular displacements in the calculation of
periodic steady state solution. The steady state solution is found when this tolerance and VelTol
are both met. The smaller the number, the tighter the tolerance is. This input is ignored if
CalcStdy is False. This input is not used in the FAST-to-ADAMS preprocessor. (rad)

\item \tt{VelTol} This is the convergence tolerance for the 2-norm of angular velocities in the calculation of the
periodic steady state solution. The steady state solution is found when this tolerance and DispTol
are both met. The smaller the number, the tighter the tolerance is. This input is ignored if
CalcStdy is False. This input is not used in the FAST-to-ADAMS preprocessor. (rad/s)
FAST Users Guide 91 Last updated on August 12, 2005 for version 6.0
Table 15. Linearization Control-Input-File Parameters (concluded).


% Model Linearization

\item \tt{NAzimStep} This is the number of equally spaced rotor azimuth steps in the output periodic linearized model.
The first rotor azimuth location is always the initial azimuth position indicated by inputs Azimuth
and AzimB1Up. The subsequent azimuth steps increment in the direction of rotation. If
RotSpeed is zero, FAST will override NAzimStep and only linearize the model about the initial
azimuth position (as if NAzimStep was set to 1). This input is not used in the FAST-to-ADAMS
preprocessor. (-)

\item \tt{MdlOrder} This is the order of the output linearized model. A setting of 1 causes FAST to output the first-
order representation of the linearized model. A setting of 2 causes FAST to output the second-
order representation of the linearized model. Using values other than 1 or 2 will cause FAST to
abort. This input is not used in the FAST-to-ADAMS preprocessor. (-)
\end{itemize}

% 

% \clearpage








\section{Relevant portions of code}
\fortran
\begin{itemize}
    \item Time stepping: \tt{Fast\_Prog.f90}, around line 100 
\begin{lstlisting}
DO n_t_global = Restart_step, Turbine(1)%p_FAST%n_TMax_m1 
  ! this takes data from n_t_global and gets values at n_t_global + 1
  DO i_turb = 1,NumTurbines
     CALL FAST_Solution_T( t_initial, n_t_global, Turbine(i_turb), ErrStat, ErrMsg )
 \end{lstlisting}

\item Linearization inputs are read in \tt{Fast\_Subs.f90} routine \tt{FAST\_ReadPrimaryFile} about line 2273

\begin{lstlisting}
! Linearize - Linearization analysis (flag)
CALL ReadVar( UnIn, InputFile, p%Linearize, "Linearize", "Linearization analysis (flag)", ErrStat2, ErrMsg2, UnEc)
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
if ( ErrStat >= AbortErrLev ) then
   call cleanup()
   RETURN        
end if      
\end{lstlisting}

\item Copy of states, \tt{Fast\_Subs.f90}, routine \tt{FAST\_Solution}, around line 4400
\begin{lstlisting}
!! STATE_CURR values of x, xd, z, and OtherSt contain values at m_FAST%t_global;
!! STATE_PRED values contain values at t_global_next.
! ElastoDyn: copy final predictions to actual states
CALL ED_CopyContState   (ED%x( STATE_PRED), ED%x( STATE_CURR), MESH_UPDATECOPY, Errstat2, ErrMsg2)
  CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
CALL ED_CopyDiscState   (ED%xd(STATE_PRED), ED%xd(STATE_CURR), MESH_UPDATECOPY, Errstat2, ErrMsg2)  
  CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
CALL ED_CopyConstrState (ED%z( STATE_PRED), ED%z( STATE_CURR), MESH_UPDATECOPY, Errstat2, ErrMsg2)      
  CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
CALL ED_CopyOtherState (ED%OtherSt( STATE_PRED), ED%OtherSt( STATE_CURR), MESH_UPDATECOPY, Errstat2, ErrMsg2)      
  CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
 \end{lstlisting}
      
\item  Elastodyn states and DOFs, e.g. \tt{ElastoDyn.f90}, around line 3600; DOF numbers are defined in \tt{ElastoDyn\_IO.f90}
\begin{lstlisting}
   x%QT (DOF_GeAz) = REAL(InputFileData%Azimuth, R8Ki) - p%AzimB1Up - REAL(Piby2_D, R8Ki)
   x%QDT(DOF_GeAz) = p%RotSpeed                                               ! Rotor speed in rad/sec.
   x%QT (DOF_Yaw) = InputFileData%NacYaw
   x%QDT(DOF_Yaw) = 0.0
 \end{lstlisting}
\item  Elastodyn outputs, e.g. \tt{ElastoDyn.f90}, around line 1727, routine \tt{ED\_CalcOutput}
\begin{lstlisting}
   y%Yaw      = x%QT( DOF_Yaw)
   y%YawRate  = x%QDT(DOF_Yaw)
   y%YawAngle = x%QT( DOF_Yaw) + x%QT(DOF_Y)  !crude approximation for yaw error... (without subtracting it from the wind direction)   
   y%BlPitch  = u%BlPitchCom !OtherState%BlPitch
   y%LSS_Spd  = x%QDT(DOF_GeAz)
   y%HSS_Spd  = ABS(p%GBRatio)*x%QDT(DOF_GeAz)
   y%RotSpeed = x%QDT(DOF_GeAz) + x%QDT(DOF_DrTr)
 \end{lstlisting}
 \item Damping ElastoDyn , routine \tt{FillAugMat}
\begin{lstlisting}
- p%CBF(K,1,1)*x%QDT(DOF_BF(K,1)) &
- p%CBF(K,1,2)*x%QDT(DOF_BF(K,2))
- p%CBF(K,2,1)*x%QDT(DOF_BF(K,1)) &
- p%CBF(K,2,2)*x%QDT(DOF_BF(K,2))
- p%CBE(K,1,1)*x%QDT(DOF_BE(K,1))
- p%CTFA(1,1)*x%QDT(DOF_TFA1) - p%CTFA(1,2)*x%QDT(DOF_TFA2)
- p%CTSS(1,1)*x%QDT(DOF_TSS1) - p%CTSS(1,2)*x%QDT(DOF_TSS2)
- p%CTFA(2,1)*x%QDT(DOF_TFA1) - p%CTFA(2,2)*x%QDT(DOF_TFA2)
- p%CTSS(2,1)*x%QDT(DOF_TSS1) - p%CTSS(2,2)*x%QDT(DOF_TSS2)
\end{lstlisting}
\item Damping for \tt{SubDyn}: \tt{p\%N2OmegaMJDamp}, for \tt{ExtPtfm}:  \tt{p\%Damp(I,J)}, for \tt{MoorDyn}: \tt{Line\%BA} \tt{Line\%Td(J,I)} \tt{Line\%B(3,I)}, for \tt{BeamDyn}: \tt{m\%DampG}, \tt{m\%DampG\_fd}


\item InputFile
\begin{lstlisting}
---------------------- LINEARIZATION -------------------------------------------
False Linearize - {Linearization analysis (flag)}
1   NLinTimes   - {Number of times to linearize (-)}
30, LinTimes    - {List of times at which to linearize (s)}
1   LinInputs   - {Inputs included in linearization (switch) {0=none; 1=standard; 2=all module inputs (debug)}}
1   LinOutputs  - {Outputs included in linearization (switch) {0=none; 1=from OutList(s); 2=all module outputs (debug)}}
False LinOutJac - {Include full Jacobians in linearization output (for debug) (flag) [unused if Linearize=False; used only if LinInputs=LinOutputs=2]}
False LinOutMod - {Write module-level linearization output files}
\end{lstlisting}

\end{itemize}







% The routine \tt{SD2ED} performs the following step:
% 
% This routine should now take \tt{TrimCase} and an offset $p_\text{off}$ in argument.
\begin{align}
    \text{ED}\%u\%p &= \text{SD}\%y\%p + p_\text{off}\\
    \text{ED}\%u\%p &= \text{SD}\%y\%p + k_\Omega \epsilon_\Omega
\end{align}


In \tt{ED\_InputSolve}:
\fortran
\begin{lstlisting}
u_ED%GenTrq     = y_SrvD%GenTrq      +  "__HACK HERE ? __"
u_ED%BlPitchCom = y_SrvD%BlPitchCom  +  "__HACK HERE ? __"
u_ED%YawMom     = y_SrvD%YawMom
\end{lstlisting}

In \tt{SrvD\_InputSolve}:
\begin{lstlisting}
u_SrvD%YawAngle  = y_ED%YawAngle !nacelle yaw plus platform yaw
u_SrvD%Yaw       = y_ED%Yaw   + "__HACK HERE ?__" !nacelle yaw
u_SrvD%BlPitch   = y_ED%BlPitch
u_SrvD%LSS_Spd   = y_ED%LSS_Spd
u_SrvD%HSS_Spd   = y_ED%HSS_Spd
\end{lstlisting}

\paragraph{Controller outputs}
Yaw calc output: 
\begin{lstlisting}
if (TrimCase==1) then
    p%YawNeut = p%YawNeut   + CtrlOffset ) 
    ! OR:
    p%YawMom = - p%YawSpr *( u%Yaw     - YawPosCom  + CtrlOffset ) &   
               - p%YawDamp*( u%YawRate - YawRateCom )                
endif
\end{lstlisting}
CalculateTorque, before each elecpwr calculation
\begin{lstlisting}
if (TrimCase==2) then
    GenTrq = GenTrq + CtrlOffset
endif
\end{lstlisting}
At the end of \tt{Pitch\_CalcOutput}:
\begin{lstlisting}
if (TrimCase==3) then
    y%BlPitchCom  = y%BlPitchCom  + CtrlOffset
endif
\end{lstlisting}
%