%+++++++++++++++++++++++++++++++++++
%            Diode Clipper
%         p = 1 N-IT scheme
%        Dr Michele Ducceschi
%       University of Bologna
%           2 April 2021
%+++++++++++++++++++++++++++++++++++

clear all
close all

%+++++++++++++++++++++++++++++++++++
% custom parameters
base_fs  = 44100;    %-- base fs
T        = 0.01 ;    %-- time
OFNIT    = 1 ;       %-- oversampling 
Is       = 2.52e-9 ; %-- sat curr
vt       = 26e-3 ;   %-- th volt
C        = 33e-9 ;   %-- cap
R        = 1e3 ;     %-- res
AmpV     = 1;        %-- input amp
sinF     = 200 ;     %-- input freq
q        = 0.5  ;    %-- free param.
%+++++++++++++++++++++++++++++++++++
% derived parameters
fs       = base_fs*OFNIT ;
k        = 1 / fs ;
Ts       = floor(T*fs) ;
tv       = (0:Ts-1).'/fs ;
vi       = AmpV*sin(sinF*2*pi*tv) ;
c1       = 1/R/C;
c2       = 2*Is/C;
c3       = c2/vt;
c4       = 1/vt ;
%+++++++++++++++++++++++++++++++++++
% init
outNIT   = zeros(Ts,1) ;
x        = 0.00 ;
%+++++++++++++++++++++++++++++++++++
% main loop
for n = 2 : Ts
    
    vin      = c1*vi(n) ;
    xvt      = c4*x ;
    sh       = sinh(xvt) ;
    ch       = sqrt(1 + sh * sh) ;
    f        = x*c1 + c2*sh ;
    fx       = c1 + c2*sh/x ;   
    if abs(x) < 1e-13
        fx = c1 + c3 ; 
    end
    fp        = c3*ch + c1;
    sigma     = q*k*fp   ;
    x         = ((1+sigma)*x-0.5*k*f+k*vin)/(1+sigma+0.5*k*fx) ;
    outNIT(n-1) = x ;
    
end
%+++++++++++++++++++++++++++++++++++
% plot
plot(tv,outNIT);
grid on;
xlabel('t (s)') ; ylabel('x (V)')
%+++++++++++++++++++++++++++++++++++