function  [out1] = RAPTOR_ip_controller(varargin)

%%% Ip feedback contol parameters

ip_controller_params = struct(......
    'method', 'VOH', ...  % control parameters used in some case % can be 'VOH' or 'PsiOH' Used if PsiOH boundary, then give VOH as command
    'kip', 3e-4, ... % proportional feedback for Ip control using VOH
    'VOH0', 3 ... % initial VOH
);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Input processing

if nargin == 0,
    out1 = ip_controller_params;
    return %empty call, probably to get default structures
elseif nargin==5
    [x,u_ff,it,model,params] = deal(varargin{:});
    ip_controller_params_in = params.controllers.ip_controller.params;
    ip_controller_params = local_procinput( ip_controller_params,{ip_controller_params_in} ); %convert list to structure, replace default by user-defined values.
else
    error('must call with 0 or 5 inputs');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Do the actual work

% TBD for this function

% more can be added
if strcmp(params.BC,'Ip');
    % standard Ip prescribed through boundary condition directly, no control
    u(1) = u_ff(1,it);
elseif strcmp(params.BC,'PsiOH') % this is not an issue unless PsiOH boundary condition set
    switch params.control.method
        case 'PsiOH'
            u(1) = u_ff(1);
            % Default, standard case
            % do nothing with u, already contains PsiOH at the right place
            
            % sensitivity follows directly from control vector
            % parametrization, no change
        case 'VOH'
            error('VOH FB control not correctly implemented yet');
            % We are actually controlling VOH instead of PsiOH
            % The B.C. is still in terms of PsiOH, so we integrate VOH numerically
            % in the end the rest boils down to a standard PsiOH case
            
            % Determine VOH
            if u(1,it) < 0; % Ip feedback case
                % negative u(1,it) means VOH proportional feedback on Ip
                %Ipref = -u(1,it); % in MA by choice to get sim magnitude as VOH
                iperror = (-1e6*u(1,it) - model.out.CIp(end,:)*Psihat(:,it));
                voh = params.control.VOH0 + params.control.kip*iperror;
                
                duk_dp(1,:) = 0; % not handled yet, so set to zero
                % handling this would involve feeding back dpsik/dt
                % since the input depends on x as well now.
                % this significantly complicates everything since it
                % provides an additional path for u and the history of u
                % to influence state evolution.
            else
                % VOH is given directly
                voh = u(1,it);
                % new du_dp for this time
                
                if it==1
                    duk_dp(1,:) = 0; % u(1,:) is independent of p
                    duk_dp_kmin1 = 0;
                else
                    %%% Handle sensitivity
                    % u = psiOH_k where
                    % psiOH_k = psiOH_k-1 + dt*VOH
                    % and VOH(k) = u(1,k) given as input
                    %
                    % dpsioh_k/dp = dpsioh_k-1/dp + dt*dvoh/dp
                    
                    % first row relates to B.C. input (PsiOH)
                    duk_dp(1,:) = duk_dp_kmin1 + dt*duk_dp(1,:);
                    duk_dp_kmin1 = duk_dp(1,:); % store previous value
                end
            end
            
            % Determine PsiOH
            if it == 1
                u(1,it) = Psihat(end,1); % start at right flux state
                duk_dp_kmin1 = 0;
            else
                % Euler integrate ohmic voltage to assign to PsiOH
                u(1,it) = u(1,it-1) + voh*dt;
                % so for this modified u we must give du_dp
            end
            
        otherwise
            error('unknown setting for params.control.method');
    end
end
out1=u;