%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%                Space-Time Convergence Stiff String Model
%                   Dr Michele Ducceschi
%                  University of Bologna
%                        29-10-2021
%+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


clear all
close all
clc

%--------------------
% custom parameters

zf      = linspace(0.01,0.0001,20) ;
hVec    = zf ;


T0         = 50 ;
rho        = 8000 ;
L          = 1 ;
r          = 2e-4 ;
E          = 2e11 ;
T          = 0.001 ;

A       = pi * r  * r ;
I       = 0.25 * pi * r  * r * r  * r ;
rA      = rho*A ;
EI      = E*I ;

for sTp = 1 : 2
    %--------------------
    
    out = zeros(length(hVec),1) ;
    Err = zeros(length(hVec),1) ;
    hvec = zeros(length(hVec),1) ;
    kvec = zeros(length(hVec),1) ;
    
    for nFs = 1 : length(hVec)
        
        nFs
        
        h      = hVec(nFs) ;
        %--------------------
        % derived parameters
        
        
        
        M       = round(L / h) ;
        if mod(M,2) ~= 0
            M = M - 1 ;
        end
        h            = L / M ;
        
           
        if sTp == 1
            k = h * sqrt(rA/T0) ;
        end
        
        
        if sTp == 2
            k       = h^2 * sqrt(rA/(T0*h^2+4*EI)) ;
        end
     
        
        Ts = round(T / k) ;
        
        hvec(nFs)    = h ;
        kvec(nFs)    = k ;
        
        %-- matrices
        Id          = speye(M-1,M-1);
        Dv          = ones(M,1);
        D2          = 1/h^2*spdiags([Dv, -2*Dv, Dv],-1:1,M-1,M-1);
        D4          = D2^2;
        Dxp         = 1/h*spdiags([-Dv, Dv],0:1,M-1,M);
        Dxm         = -Dxp.';
        
        if sTp == 1
            Cp = rA*Id/k^2 + 0.5*EI*D4  ;
            C0 = 2*rA/k^2*Id + T0*D2 ;
            Cm = -rA*Id/k^2 - 0.5*EI*D4 ;
        end
        
        if sTp == 2
            Cp = rA/k^2*Id   ;
            C0 = 2*rA/k^2*Id + T0*D2 - EI*D4 ;
            Cm = -rA/k^2*Id  ;
        end
        
        
        disp0   = zeros(M-1,1) ;
        
        for m = 1 : M - 1
            if m * h >= 0.25 * L && m * h <= 0.75 * L
                disp0(m) = 0.5 * ( 1 - cos( 4*pi*(m*h-0.25*L) / L ) );
            end
        end
        

        
        xm      = disp0 ;
        x0      = xm + 0.5 * k^2 * (T0/rA * D2 - EI / rA * D4) * xm ;
        
        for n = 1 : Ts
            
            xp = Cp \ (C0*x0 + Cm*xm) ;
            

            
            if n == Ts
                exsol = 0 ;
                for m = 1 : 1e6
                    Om = sqrt(T0 / rA * (m * pi / L)^2 + EI / rA * (m * pi / L)^4) ;
                    Cm = (Coeffm(m,0.75*L,L) - Coeffm(m,0.25*L,L)) ;
                    Xmat = sin( m * pi * 0.5 ) ;
                    exsol = exsol + Cm * cos(Om*(n+1)*k) * Xmat ;
                end
            end
            
            xm = x0 ;
            x0 = xp ;
            
        end
        

        Err(nFs) =  ( xp(M/2 + 1) - exsol ) ;
        
    end
    
    %Err = abs(out-exsol) ;
    subplot(2,1,1)
    if sTp == 1
    loglog(kvec,abs(Err),'k','marker','x') ; grid on ; hold on ;
    else
        loglog(kvec,abs(Err),'k','marker','o') ; grid on ; hold on ;
    end
    xlabel('k (s)') ; ylabel('|E|') ; legend('Implicit','Explicit')
    subplot(2,1,2)
     if sTp == 1
    loglog(hvec,abs(Err),'k','marker','x') ; grid on ; hold on ;
    else
        loglog(hvec,abs(Err),'k','marker','o') ; grid on ; hold on ;
    end
    xlabel('h (m)') ; ylabel('|E|') ; legend('Implicit','Explicit')
    
end

function y = Coeffm(m,x,L)

if m == 4
    y = 0 ;
else
    y = -1/m/pi*cos(m*pi*x/L) - 1/2/(m+4)/pi*cos((m+4)*pi*x/L) - 1/2/(m-4)/pi*cos((m-4)*pi*x/L) ;
end
end