bht_step2(t,epsilon,a)

function [x,dx,ddx] = bht_step2(t,epsilon,a)
% BHT_STEP2
% 
% Syntax
% 
%     [x,dtx,dttx] = BHT_STEP2(t,epsilon,h)
% 
% Description
% 
%     This function returns the values x and the values dtx and dttx of
%     the first and second derivatives of the 2*pi-periodic twice
%     continuously step functions, evaluated at the point t.
%     BHT_STEP1 is called by the function BHT_ROUND_RECTANGLE. 
%     The input epsilon is a smooting parameter (it set how much the angles 
%     are rounded and has to be within the interval ]0,1]) and h is the 
%     steps' height. 
%
% See also BHT_ROUND_RECTANGLE

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                         %
%                      BIOHYDRODYNAMICS MATLAB TOOLBOX                    %
%                                                                         %
%                         A. MUNNIER and B. PINCON                        %
%                                                                         %
% [e-mail]        [e-mail]   %
% http://www.iecn.u-nancy.fr/~munnier  http://www.iecn.u-nancy.fr/~pincon % 
%                                                                         %
%                      INSTITUT ELIE CARTAN, NANCY 1                      %
%                       http://www.iecn.u-nancy.fr/                       %
%                                                                         %
%                      INRIA Lorraine, Projet CORIDA                      %
%                    http://www.iecn.u-nancy.fr/~corida/                  %
%                                                                         %  
%                                                                         %  
%                                                                         % 
%                            August 15th 2008                             %
%                                                                         %
%                                               GNU GPL v3 licence        %
%                                                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%##########################################################################
P = [35/128,0,-45/32,0,189/64,0,-105/32,0,315/128,0];
Pp = [315/128,0,-315/32,0,945/64,0,-315/32,0,315/128];
Ppp = [315/16,0,-945/16,0,945/16,0,-315/16,0];

%--------------------------------------------------------------------------
n = length(t);
t = reshape(t,1,n);

%--------------------------------------------------------------------------
I = find(t<=epsilon);
x(I) = -a.*polyval(P,t(I)./epsilon);
dx(I) = -a.*polyval(Pp,t(I)./epsilon)./epsilon;
ddx(I) = -a.*polyval(Ppp,t(I)./epsilon)./epsilon^2;

%--------------------------------------------------------------------------
I = find(t>epsilon & t<=pi-epsilon);
x(I) = -a;
dx(I) = 0;
ddx(I) = 0;

%--------------------------------------------------------------------------
I = find(t>pi-epsilon & t<=pi+epsilon);
Il = length(I);
x(I) = a.*polyval(P,(t(I)-pi.*ones(1,Il))./epsilon);
dx(I) = a.*polyval(Pp,(t(I)-pi.*ones(1,Il))./epsilon)./epsilon;
ddx(I) = a.*polyval(Ppp,(t(I)-pi.*ones(1,Il))./epsilon)./epsilon^2;

%--------------------------------------------------------------------------
I = find(t>pi+epsilon & t<=2*pi-epsilon);
x(I) = a;
dx(I) = 0;
ddx(I) = 0;

%--------------------------------------------------------------------------
I = find(t>2*pi-epsilon);
Il = length(I);
x(I) = -a.*polyval(P,(t(I)-2*pi.*ones(1,Il))./epsilon);
dx(I) = -a.*polyval(Pp,(t(I)-2*pi.*ones(1,Il))./epsilon)./epsilon;
ddx(I) = -a.*polyval(Ppp,(t(I)-2*pi.*ones(1,Il))./epsilon)./epsilon^2;