function [ Vout ] = Diode_Clipper_AA_2( Vin )
% Diode Clipper Antialiased using the 2nd-order antiderivative form

% Output Vector
Vout = zeros(size(Vin));

% Components and Parameters
R = 2.2e3;
Is = 2.52e-9;
VT = 25.85e-3;
N = 1.752;

% Repeated constants
NVT = N*VT;
NVT_sq = NVT*NVT;
NVT_cub = NVT_sq*NVT;
one_NVT = 1/(NVT);
one_sixth = 1/6;

% State variables
Fn1 = 0;
Gn1 = 0;
xn1 = 0; xn2 = 0;

% Error threshold for ill-conditioning
thresh = 10e-10;

% Compute alfa
alfa = (R*Is)*one_NVT;


for n=1:length(Vin)
    
    xn = Vin(n);

    lambda = sign(xn);
    
    beta = lambda*one_NVT;

    % Compute Lambert First
    Lu = Lambert_W_Fritsch(alfa*exp(beta*xn));
    
    % Compute second antiderivative
    Fn = (1/6)*( (xn*xn*xn) - lambda*0.5*NVT_cub*Lu*( Lu * ( 2*Lu + 9) + 12 ) );

    % Central difference
    xn_xn1 = xn - xn1;
    
    % First round of antialiasing
    if (abs(xn_xn1)<thresh)
        
        % First-order escape rule
        x_avg = 0.5*(xn + xn1);
        % Compute FIRST antiderivative (at x_avg)
        lambda = sign(x_avg);
        beta = lambda*one_NVT;
        Lu = Lambert_W_Fritsch(alfa*exp(beta*x_avg));
        Gn = 0.5*( (x_avg*x_avg) - NVT_sq*Lu*(Lu + 2) );
        
    else
        
        % Antiderivative method
        Gn = (Fn - Fn1) / xn_xn1;
        
    end

    % Open difference
    xn_xn2 = xn-xn2;
    
    if (abs(xn_xn2)<thresh)
                                
        % Second-order escape rules
        delta = 0.5 * (xn - 2*xn1 + xn2);
        if abs(delta) < thresh
                        
            % Compute original function
            x_avg = 0.25*(xn + 2*xn1 + xn2);
            lambda = sign(x_avg);
            beta = lambda*one_NVT;
            Lu = Lambert_W_Fritsch(alfa*exp(beta*x_avg));
            Vout(n) = x_avg - lambda*NVT*Lu; 
            
        else
            
            x_avg = 0.5*(xn + xn2);
            
            lambda = sign(x_avg);
            beta = lambda*one_NVT;
            Lu = Lambert_W_Fritsch(alfa*exp(beta*x_avg));
         
            % Compute FIRST antiderivative
            temp1 = 0.5*( (x_avg*x_avg) - NVT_sq*Lu*(Lu + 2) );
            
            % Compute SECOND antiderivative
            temp2 = one_sixth*( (xn*xn*xn) - lambda*0.5*NVT_cub*Lu*( Lu * ( 2*Lu + 9) + 12 ) );
            
            % Output
            delta = 1/delta;
            Vout(n) = 2*delta*(temp1 + delta*(Fn1 - temp2));
           
        end

    else
        
        % Antiderivative Form
        Vout(n) = 2 * (Gn - Gn1)/xn_xn2;
        
    end
    
    % Update state variables
    Fn1 = Fn;
    Gn1 = Gn;
    xn2 = xn1;
    xn1 = xn;
    
end  

end

