function [sys,x0,str,ts] = sfunmoment(t, x, u, flag,...
				      parameters),
%       An S-function which computes the Mirror Moment
%
%       The input is in the Form:
%       u = [V_M(t), phi(t), V_L(t), V_R(t)]
%
%       The parameter structure has the following entries:
%       d0: gap width
%       lm: mirror lenght
%       bm: mirror width
%
%       Heiko Wolfram 2005.
%
%       $Id: sfunmoment.m,v 0.1 2005/08/12 12:57:04 hwol Exp $

%   What is returned by SFUNC at a given point in time, T, depends on the
%   value of the FLAG, the current state vector, X, and the current
%   input vector, U.
%
%   FLAG   RESULT             DESCRIPTION
%   -----  ------             --------------------------------------------
%   0      [SIZES,X0,STR,TS]  Initialization, return system sizes in SYS,
%                             initial state in X0, state ordering strings
%                             in STR, and sample times in TS.
%   1      DX                 Return continuous state derivatives in SYS.
%   2      DS                 Update discrete states SYS = X(n+1)
%   3      Y                  Return outputs in SYS.
%   4      TNEXT              Return next time hit for variable step sample
%                             time in SYS.
%   5                         Reserved for future (root finding).
%   9      []                 Termination, perform any cleanup SYS=[].


  switch flag,
    
    %%%%%%%%%%%%%%%%%%
    % Initialization %
    %%%%%%%%%%%%%%%%%%
   case 0,
    [sys,x0,str,ts]=mdlInitializeSizes;
    
    %%%%%%%%%%%
    % Outputs %
    %%%%%%%%%%%
   case 3,
    sys=mdlOutputs(t, x, u, parameters);
    
    %%%%%%%%%%%%%%%%
    % Unused flags %
    %%%%%%%%%%%%%%%%
   case { 1, 2, 4, 9 },
    sys = [];
    
   case 'moment',
    sys = LocalMoment(parameters, u, x);
    
    %%%%%%%%%%%%%%%%%%%%
    % Unexpected flags %
    %%%%%%%%%%%%%%%%%%%%
   otherwise
    error(['Unhandled flag = ',num2str(flag)]);
    
  end % case flag
  
% end  sfun_moment
    
%
%=============================================================================
% mdlInitializeSizes
% Return the sizes, initial conditions, and sample times for the S-function.
%=============================================================================
%
function [sys,x0,str,ts] = mdlInitializeSizes 
  
  % simsizes for a sizes structure
  sizes = simsizes;
  
  sizes.NumContStates  = 0; % number of continuous states
  sizes.NumDiscStates  = 0; % number of discrete states
  sizes.NumOutputs     = 1; % number of system outputs (length of output y)
  sizes.NumInputs      = 4; % number of system inputs (length of input u)
  sizes.DirFeedthrough = 1; % direct feedthrough flag; used to detect algebraic loops.
  sizes.NumSampleTimes = 0; % at least one sample time is needed
    
  sys = simsizes(sizes);
  
  % initialize the initial conditions
  x0  = []; 
  
  % str is always an empty matrix
  str = [];
  
  % initialize the array of sample times
  ts = [];
% end mdlInitializeSizes


%
%=============================================================================
% mdlOutputs
% Return the block outputs.
%=============================================================================
%
function sys = mdlOutputs(t,x,u,params)

  sys = LocalMoment(params, u(1) + u(4), u(2)) -...
	LocalMoment(params, u(1) + u(3), -u(2));

% end mdlOutputs


%
%=============================================================================
% LocalMoment
% Calculate Moment for one Plate.
%=============================================================================
%
function yout = LocalMoment(params, u, phi)
  
  epsilon0 = 8.8541978E-12;
  
  if abs(phi) < sqrt(eps),
      
    % Taylor Approximation
    yout = (params.bm^2*params.lm*epsilon0*...
	    (-3*params.d0 + 2*params.bm*phi)*u^2)/(48*params.d0^3);
    
  else;
    
    yout = params.lm*epsilon0*u^2*cos(phi)/sin(phi)^2*...
	   ( (params.bm*sin(phi))/(4*params.d0 + 2*params.bm*sin(phi)) - ...
	     atanh((params.bm*sin(phi))/(4*params.d0 + params.bm*sin(phi))) ...
	     );
  end;
  
% end LocalMoment


%% End of File