% Il codice permette l'individuazione di una configurazione di
% danneggiamento. La rigidezza viene decurtata della stessa aliquota per
% ogni molla ma esse possono essere soggette a valori di decurtazioni
% iniziali differenti

clear all; close all; clc

%% Arch Analysis
%% Input parameters
Arch.f = 30;                                            % Height of the arch [m]
Arch.L = 100;                                           % Span length [m]
Arch.Lhalf = Arch.L/2;
n_1 = 100;                                               % Number of intervals of the mesh [/]
NP_1 = n_1+1;                                           % Number of points of the mesh [/]
q0 = -100;                                              % Uniformly distributed load [kN/m]

%% Material
Material.gamma = 78.5;                                  % Specific weight [kN/m3]
Material.E = 2.1*10^8;                                  % Young modulus [kN/m2]
Material_E = Material.E;
Material.nu = 0.3;                                      % Poisson ratio [-]
Material.G = (Material.E)/(2*(Material.nu+1));          % Shear modulus [kN/m2]
Material.fyd = 355;                                     % Design strength [MPa]

Section.r_base = 0.500;                                 % Sectional 1st dimension at the arc base[m]
Section.r = 0.5*Section.r_base;                         % Reduction factor
Section.r_base2 = 0.5;                                  % Sectional 2nd dimension at the arc base[m]
Section.t = 0.05;                                       % Thickness of hollowed sections [m]
ratio = Section.r_base2/Section.r_base;                 % Ratio between the section dimensions

%% User choices
Arch.curve1 = input("Insert the selected curve type:\n Parabolic = 1\n Circular = 2\n4th order polybomial = 3\n");
Section.variation = input("Insert the type of variation of the section dimension:\n Constant = 1\n Linear = 2\n Quadratic = 3\n");
Section.type = input("Insert the selected section type:\n Circular = 1\n Ellipsoidal = 2\n Squared = 3\n Rectangular = 4\n Hollowed circle = 5\n Hollowed square = 6\n Hollowed rectangle = 7\nInsert the correspondent number:\n");
Restraint = input("Insert the selected restraint condition at the base: \n Fully restrained arch = 1\n Hinged arch = 2\n");

tic

%% Definizione curva iniziale

syms X_1

if Arch.curve1 == 1
    % Parabola
    y_1 = (-4*Arch.f/(Arch.L^2))*(X_1^2) + Arch.f;                                                                  % Parabola equation
elseif Arch.curve1 ==2
    % Circumference
    % y_1 = sqrt(-X_1^2-(-(Arch.L^2)/4)+(((Arch.L^2)/(4*Arch.f)-Arch.f)^2)/4)-((Arch.L^2)/(4*Arch.f)-Arch.f)/2;     % Circumference equation
    coeffA_1 = 0;
    coeffB_1 = ((Arch.L^2)/(4*Arch.f)) + Arch.f;
    coeffC_1 = 0;
    y_1 = (sqrt(-X_1^2 - (coeffA_1*X_1) - coeffC_1 + (coeffB_1^2)/4) - coeffB_1/2) + Arch.f;                        % Circumference equation
elseif Arch.curve1 ==3
    m = 4.0*Arch.f/Arch.L;
    y_1 = (-4*m/(Arch.L^3) + 16*Arch.f/(Arch.L^4))*X_1^4 + (m/Arch.L - 8*Arch.f/(Arch.L^2))*X_1^2 + Arch.f;         % 4-th order polynomial
end


% Derive the curve
dy_1 = diff(y_1,X_1);                 % 1st derivative of the parabola
d2y_1 = diff(dy_1,X_1);               % 2nd derivative of the parabola


% Calculation of the incremental intervals
x_inter_1 = [-Arch.L/2 -Arch.L/2 + linspace(1,n_1,n_1).*Arch.L/n_1];
y_inter_1 = [0 double(subs(y_1,X_1,x_inter_1(2:n_1+1)))];
dy_inter_1 = y_inter_1(2:n_1+1) - y_inter_1(1:n_1);

% Allocation of memory space for variables (fastening the calculation)
s_x_1 = zeros(1,n_1+1); 

% Calculation of s(x)
for i = 1:n_1
    s_x_1(i+1) = s_x_1(i) + sqrt((Arch.L/n_1)^2 + dy_inter_1(i)^2);
end

%% Equispaced nodes along s

ds_1 = s_x_1(end)/n_1;
s_1 = ds_1.*(linspace(0,n_1,n_1+1));
s_exp_1 = [-3*ds_1 -2*ds_1 -1*ds_1 s_1 s_1(end)+ds_1 s_1(end)+2*ds_1 s_1(end)+3*ds_1];
x_s_1 = interp1(s_x_1,x_inter_1,s_1,'pchip','extrap');

%% Vettore punti iniziali

rn = [1 20 51 82 101];
contatore = 1;
i = 1;

% Valore iniziale rigidezza molla alla base a dx 
perc_dx = 0.7;

%% %% Inizio ciclo
while i == 1

    % Ordino eventualmente il vettore
    rn = sort(rn,'ascend');
    
    %% Eliminazione punti doppioni

    % Creo un vettore in cui metto tutte le posizoni dei punti doppioni che
    % quindi vanno incontro a snervamento più volte
    copia = [];
    if length(rn) > 1
        j = 1;
        for i = 1:length(rn)-1
            if rn(i+1) == rn(i)
               copia(j) = i+1;
               j = j+1;
            end
        end
    % Elimino i punti che si ripetono all'interno del vettore in quanto la
    % condizione sulla riduzione della rigidezza avviene tramite
    % memorizzazione di un ulteriore vettore con il quale si riduce la
    % rigidezza 
        if isempty(copia) == 0
            rn(copia) = [];
        end
    end
    
    
    %% Pulizia vettore punti danneggiati
    
    % I punti in esame devono distare almeno due nodi per poter garantire le
    % interpolazioni successive. Con il codice successivo si eliminano dal
    % vettore i punti che distano un nodo dal precedente
    xx = length(rn);
    i=1;
    j = 1;
    while i < xx
        if rn(i+1) == rn(i)+1
           rn(i+1) = [];
           xx = xx-1;
            
        end
        i=i+1;
        j=j+1;
        if j == length(rn)
            break
        end   
    end
    
    % Devo fare in modo che il primo punto e l'ultimo non ricadano all'interno
    % del vettore che definisce la posizione delle molle altrimenti non
    % funziona l'interpolazione.

    % Inizializzazione variabili
    if contatore == 1
        molla_in = 0;
        molla_fin = 0;
        
        if  rn(1) == 1
            rn(1) = [];
            molla_in = 1;
        end
        
        if rn(end) == NP_1
        rn(end) = [];
        molla_fin = NP_1;
        end
    end


    %% Inizio del calcolo per archi separati

    dann = length(rn);

    %Inizizalizzazione vettori
    s_x_point = zeros (1,dann);
    x_s_point = zeros(1,dann);
    y_s_point = zeros(1,dann);
    s_arco = zeros(1,dann);
    n = zeros (1,dann+1);
    Arch.Li = zeros(1,dann);
    Arch.Lfin = zeros(1,dann);
    
    % Definisco il vettore che individua la posizione delle ceriere
    x_point = [x_s_1(rn(1:dann))];

    % Parametro che serve per allocare in memoria i valori all'interno della
    % matrice nella corretta posizione
    posizione = rn;
    
    % Inizializzazione parametri successivi per i diversi cicli. 
    kk = 0;    jj = 0;    ll = 0;    hh = 0;    ff = 0;    zz = 0;    index = 0;
    
    % Ordino il vettore in modo da avere una sequenza ordinata delle lunghezze
    % relative ad ogni tratto di arco. Tuttavia questo passaggio adesso risulta
    % superfluo in quanto il vettore che indica le posizioni in cui si
    % formeranno le cerniere risulta già ordinato fuori dal for
    
    x_point = sort(x_point,'ascend');
    n_a = length(x_point)+1;                                                                %conta numero archi da passare al solver
    
    
    for ind = 1:length(x_point)
        s_x_point(ind) = interp1(x_inter_1,s_x_1,x_point(ind),'pchip','extrap');            %estrapolo la s corrispondente ad un punto x qualsiasi lungo l'arco
        x_s_point(ind) = interp1(s_x_1,x_inter_1,s_x_point(ind),'pchip','extrap');          %estrapolo la x corrispondente alla s estrapolata in precedenza
        y_s_point(ind) = interp1(x_inter_1, y_inter_1, x_s_point(ind),'pchip','extrap');    %definizione della y in corrispondenza della x determinata
    end
    
    % Tecnicamente non avrei bisogno di interpolare x_s_point in quanto
    % x_point cade in corrispondenza di punti facenti parte della mesh e
    % quindi in questo caso i due vettori dovrebbero coincidere
    
    %interpolazione per trovare x(s) cioè la x in corrispondeza dei segmenti della s e quindi trovo i punti
    %x(s) della funzione x_inter nei punti corrispondenti a quelli presenti nel vettore s
            
            
            
    %% Definizione punti per cui passano le parabole
        
    % Definisco le lunghezze dei due tratti di arco riportando il sitema di
    % riferimento all'estremità in basso a sinistra della prima curva
    
    for ind = 1:length(x_point)    
        
        x_s_point(ind) = x_s_point(ind)+Arch.Lhalf;
        x_point(ind) =   x_point(ind)+Arch.Lhalf;
    end
    
    for ind = 1:length(x_point)
    
        if ind == 1 
            Arch.Li(1) = x_point(1);
            Arch.Lfin = Arch.L -x_point(1);
            s_arco(1) = s_x_point(1);
        else
            Arch.Li(ind) = x_point(ind) - x_point(ind-1);
            s_arco(ind) = s_x_point (ind) - s_x_point(ind-1);
        end
    
        % Indica l'altezza di ogni tratto di curva rispetto all'altezza
        % iniziale in chiave della curva completa
        Arch.f2(ind) = Arch.f - y_s_point(ind);
    
    end
    
    % Calcolo fuori dal for l'ultimo valore del vettore nel seguente modo
    Arch.Lfin = Arch.L - sum(Arch.Li);
    s_arch_fin = s_1(end) - sum(s_arco);
    
    % Definisco compiutamente la lunghezza dei tratti di arco
    Arch.Li = [Arch.Li Arch.Lfin];
    s_arco_tot = [s_arco s_arch_fin];
    
    %% Definizione dell'intervallo di discretizzazione per ogni arco
    
    % Allocazione di memoria
    NP = zeros(1,n_a+1);
    HP = zeros(1,n_a);
    
    for ind = 1:length(x_point)+1
    
        %definire n(ind) in modo idoneo e quindi facendo coincidere i putni dei
        %nuovi tratti di arco con i vecchi punti dell'arco completo
        n(1) = 0;
        NP(1) = 0;
        n(ind+1) = round(s_arco_tot(ind)/ds_1);
        NP(ind+1) = n(ind+1) + 1;
    
        n_max = max(n);
        NP_max = max(NP);
    end

    for ind = 1:n_a
        HP(1) = 0;
        if ind ~= 1
            HP(ind) = NP(ind) + HP(ind-1) + 6;
        end
    end
    
    %% 
    % Allocazione di memoria
    x_inter = zeros(length(x_point)+1,NP_max);
    y_inter = zeros(length(x_point)+1,NP_max);
    dy_inter = zeros(length(x_point)+1,NP_max);
    R_s = zeros(length(x_point)+1,NP_max);
    alpha = zeros(length(x_point)+1,NP_max);
    s = zeros(length(x_point)+1,NP_max);
    s_exp = zeros(length(x_point)+1,NP_max+6);
    x_s = zeros(length(x_point)+1,NP_max);
    ds = zeros(1,n_a);
    coeffA = zeros(1,n_a);
    coeffB = zeros(1,n_a);
    coeffC = zeros(1,n_a);
    z1 = zeros(n_a,NP_max);
    
     
    %% Definizione delle diverse curve
    
    % Calcolo dei punti per cui passano gli archi
    x_point_arco = [Arch.L (Arch.L-2.*x_point)];
    Arch.f_arco = [Arch.f Arch.f2];
    
    % Se deve essere inserito un vincolo perfettamente in mezzeria devo
    % usare un'equazione di definizione dela curva differente e quindi in
    % questo modo tengo conto della posizione nel vettore rn per capire
    % quando applicare tale equazione.
    
    intermedio = 0;
    for i = 1:length(x_point_arco)
       if x_point_arco(i) < 0.0001 && x_point_arco(i) > -0.0001
           x_point_arco(i) = 0;
          intermedio = i;
          break
       end
    end
    
    for ind = 1:length(x_point)+1
    
        syms X
        
        if Arch.curve1 == 1
            % Parabola equation
            if ind == intermedio
                y = (-4*Arch.f/(Arch.L^2))*(X^2); 
            else
                y = (-4*Arch.f_arco(ind)/(x_point_arco(ind).^2))*(X^2) + (4*Arch.f_arco(ind)/x_point_arco(ind))*X;        % Parabola equation
            end
        elseif Arch.curve1 == 2
            % Circumference
            if ind == intermedio
                % Circumference    
                coeffA(ind) = 0;
                coeffB(ind) = ((Arch.L^2)/(4*Arch.f)) + Arch.f;
                coeffC(ind) = 0;
                y = sqrt(-X^2 - (coeffA(ind)*X) - coeffC(ind) + (coeffB(ind)^2)/4) - coeffB(ind)/2;                       % Circumference equation
            else
                % Circumference
                coeffA(ind) = -x_point_arco(ind);
                coeffB(ind) = ((x_point_arco(ind).^2)./(4*Arch.f_arco(ind))-Arch.f_arco(ind));
                coeffC(ind) = 0;
        
                y = sqrt(-X^2 -(coeffA(ind)*X)-coeffC(ind)+((coeffB(ind)^2)/4)) - coeffB(ind)/2;                          % Circumference equation
            end
        end
        
        
        % Derive the curve
        dy = diff(y,X);                 % 1st derivative of the parabola
        d2y = diff(dy,X);               % 2nd derivative of the parabola
    
        %% Calculation of s(x)
        
        % Il numero di intervalli deve essere compatibile con quello della
        % prima curva e quindi la somma del numero dei punti delle diverse
        % curve uguale al numero dei punti iniziali della prima curva diminuito
        % del numero di punti corrispondenti al numero di curve
    
        % Calculation of the incremental intervals
        x_inter(ind,1:NP(ind+1)) = [0 linspace(1,n(ind+1),n(ind+1)).*Arch.Li(ind)/n(ind+1)];
        y_inter(ind,1:NP(ind+1)) = [0 double(subs(y,X,x_inter(ind,2:n(ind+1)+1)))];
    
        dy_inter(ind,1:n(ind+1)) = y_inter(ind,2:n(ind+1)+1) - y_inter(ind,1:n(ind+1));
            
        % Allocation of memory space for variables (fastening the calculation)
        % s_x(ind,NP_max) = zeros(length(x_point)+1,NP_max);           
        s_x(ind,1) = 0;                                                             % Curvilinear coordinates
    
        % Calculation of s(x)
        for i = 1:n(ind+1)
            s_x(ind,i+1) = s_x(ind,i) + sqrt((Arch.Li(ind)/n(ind+1))^2 + dy_inter(ind,i)^2);
        end
        
        %% Equispaced nodes along generic "s"
    
        ds(ind) = s_x(ind,NP(ind+1))/n(ind+1);
        s(ind,1:NP(ind+1)) = ds(ind).*(linspace(0,n(ind+1),n(ind+1)+1));
        s_exp(ind,1:NP(ind+1)+6) = [-3*ds(ind) -2*ds(ind) -1*ds(ind) s(ind,1:NP(ind+1)) s(ind,NP(ind+1))+ds(ind) s(ind,NP(ind+1))+2*ds(ind) s(ind,NP(ind+1))+3*ds(ind)];
        x_s(ind,1:NP(ind+1)) = interp1(s_x(ind,1:NP(ind+1)),x_inter(ind,1:NP(ind+1)),s(ind,1:NP(ind+1)),'pchip','extrap');
    
        %% Curvature radius and its derivatives for the generic curve
        
        % Mathematical form
        R_x_math = -((1+(dy)^2)^(3/2))/d2y;                                         % R(x)
        
        R_s(ind,1:NP(ind+1)) = double(subs(R_x_math,X,x_s(ind,1:NP(ind+1))));       % R[s(x)] [m]
        alpha(ind,1:NP(ind+1)) = atan(double(subs(dy,X,x_s(ind,1:NP(ind+1)))));
        z1(ind,1:NP(ind+1)) = double(subs(y,X,x_s(ind,1:NP(ind+1))));               % Coordinata y dei punti della curva [m]

        % Servirà per dopo nella definizione dei carichi ma successivamente avrei dovuto definire
        % nuovamente le curve, per cui è conveniente definirlo adesso.
    end

    
    %% Allocazione di memoria
    J = zeros((length(x_point)+1),NP_max);
    A = zeros((length(x_point)+1),NP_max);
    W = zeros((length(x_point)+1),NP_max);
    J_fit = zeros(length(x_point)+1,5);
    
    r1 = zeros((length(x_point)+1),NP_max);
    r2 = zeros((length(x_point)+1),NP_max);
    t = zeros(1,length(x_point)+1);
    
    dR_s = zeros(length(x_point)+1,NP_max);
    d2R_s = zeros(length(x_point)+1,NP_max);
    d3R_s = zeros(length(x_point)+1,NP_max);
    d4R_s = zeros(length(x_point)+1,NP_max);
    d5R_s = zeros(length(x_point)+1,NP_max);
    
    dJ = zeros((length(x_point)+1),NP_max);
    d2J = zeros((length(x_point)+1),NP_max);
    d3J = zeros((length(x_point)+1),NP_max);
    
    q_sw = zeros((length(x_point)+1),NP_max);
    q_ol = zeros((length(x_point)+1),NP_max);
    Pn = zeros((length(x_point)+1),NP_max);
    Pt = zeros((length(x_point)+1),NP_max);
    dPn = zeros((length(x_point)+1),NP_max);
    m_load = zeros ((length(x_point)+1),NP_max);
    
    a4 = zeros((length(x_point)+1),NP_max);
    b3 = zeros((length(x_point)+1),NP_max);
    b4 = zeros((length(x_point)+1),NP_max);
    c2 = zeros((length(x_point)+1),NP_max);
    c3 = zeros((length(x_point)+1),NP_max);
    c4 = zeros((length(x_point)+1),NP_max);
    d1 = zeros((length(x_point)+1),NP_max);
    d2 = zeros((length(x_point)+1),NP_max);
    d3 = zeros((length(x_point)+1),NP_max);
    d4 = zeros((length(x_point)+1),NP_max);
    e1 = zeros((length(x_point)+1),NP_max);
    e2 = zeros((length(x_point)+1),NP_max);
    e3 = zeros((length(x_point)+1),NP_max);
    e4 = zeros((length(x_point)+1),NP_max);
    f1 = zeros((length(x_point)+1),NP_max);
    f2 = zeros((length(x_point)+1),NP_max);
    f3 = zeros((length(x_point)+1),NP_max);
    f4 = zeros((length(x_point)+1),NP_max);
    g1 = zeros((length(x_point)+1),NP_max);
    g2 = zeros((length(x_point)+1),NP_max);
    g3 = zeros((length(x_point)+1),NP_max);
    g4 = zeros((length(x_point)+1),NP_max);
    a = zeros((length(x_point)+1),NP_max);
    b = zeros((length(x_point)+1),NP_max);
    c = zeros((length(x_point)+1),NP_max);
    d = zeros((length(x_point)+1),NP_max);
    e = zeros((length(x_point)+1),NP_max);
    f = zeros((length(x_point)+1),NP_max);
    g = zeros((length(x_point)+1),NP_max);
    h = zeros((length(x_point)+1),NP_max);
    aM = zeros((length(x_point)+1),NP_max);
    bM = zeros((length(x_point)+1),NP_max);
    cM = zeros((length(x_point)+1),NP_max);
    dM = zeros((length(x_point)+1),NP_max);
    aT = zeros((length(x_point)+1),NP_max);
    bT = zeros((length(x_point)+1),NP_max);
    cT = zeros((length(x_point)+1),NP_max);
    dT = zeros((length(x_point)+1),NP_max);
    eT = zeros((length(x_point)+1),NP_max);
    aN = zeros((length(x_point)+1),NP_max);
    bN = zeros((length(x_point)+1),NP_max);
    cN = zeros((length(x_point)+1),NP_max);
    dN = zeros((length(x_point)+1),NP_max);
    eN = zeros((length(x_point)+1),NP_max);
    fN = zeros((length(x_point)+1),NP_max);
    
    %% 
    for ind = 1:length(x_point)+1    
        % Calculation of the curvature radius and its derivatives(FDM)             
        dR_s(ind,1:NP(ind+1)) =  interp1(s(ind,1:n(ind+1))+ds(ind)/2,diff(R_s(ind,1:NP(ind+1)))./ds(ind),s(ind,1:NP(ind+1)),'pchip','extrap');
        d2R_s(ind,1:NP(ind+1)) = interp1(s(ind,1:n(ind+1))+ds(ind)/2,diff(dR_s(ind,1:NP(ind+1)))./ds(ind),s(ind,1:NP(ind+1)),'pchip','extrap');
        d3R_s(ind,1:NP(ind+1)) = interp1(s(ind,1:n(ind+1))+ds(ind)/2,diff(d2R_s(ind,1:NP(ind+1)))./ds(ind),s(ind,1:NP(ind+1)),'pchip','extrap');
        d4R_s(ind,1:NP(ind+1)) = interp1(s(ind,1:n(ind+1))+ds(ind)/2,diff(d3R_s(ind,1:NP(ind+1)))./ds(ind),s(ind,1:NP(ind+1)),'pchip','extrap');
        d5R_s(ind,1:NP(ind+1)) = interp1(s(ind,1:n(ind+1))+ds(ind)/2,diff(d4R_s(ind,1:NP(ind+1)))./ds(ind),s(ind,1:NP(ind+1)),'pchip','extrap');
        
        %% Transversal section properties along s(x)
    
        % Section dimensions for th entire arch
        if Section.variation == 1
            r1_1(1:NP_1) = Section.r_base;                                                                                                        % Sections radius [m]
        elseif Section.variation == 2
            r1_1 = ((Section.r_base - Section.r)./(s_1(end)./2)).*abs(s_1 - s_1(end)./2) + Section.r;                                 % Sections radius [m] 
        elseif Section.variation == 3
            r1_1 = Section.r_base+(4/s_1(end))*(Section.r-Section.r_base).*s_1-(4/(s_1(end)^2))*(Section.r-Section.r_base).*(s_1.^2); % Sections radius [m]
        end
        
        if Section.variation == 1
            r2_1(1:NP_1) = Section.r_base2;                                % Sections radius [m]
        else
            r2_1 = ratio*r1_1;
        end
        
        % Determino l'andamento dei raggi di curvatura per metà arco,
        % associando i raggi dell'arco intero nelle corrette posizioni
        if ind == 1
            r1(ind,1:NP(ind+1)) = r1_1(1:NP(ind+1));                       % devo lasciare come indice "ind" in quanto parte da 2 e il primo vettore si definisce manualmente
            r2(ind,1:NP(ind+1)) = r2_1(1:NP(ind+1));
            t(ind) = Section.t;                                            % Thickness of the hollowed section [m]
        else
            zz = zz + NP(ind);
            r1(ind,1:NP(ind+1)) = r1_1(zz-ind+2:zz+NP(ind+1)-ind+1);       % devo lasciare come indice "ind" in quanto parte da 2 e il primo vettore si definisce manualmente
            r2(ind,1:NP(ind+1)) = r2_1(zz-ind+2:zz+NP(ind+1)-ind+1);
            t(ind) = Section.t;                                            % Thickness of the hollowed section [m]
        end
        
        
        %% Transversal section properties along s(x)
        
        % Geometrical parameters
        if Section.type == 1        % Circular
            J(ind,1:NP(ind+1)) = (pi.*(r1(ind,1:NP(ind+1)).^4))./4;                                                         % Inertia [m4]
            A(ind,1:NP(ind+1)) = pi.*(r1(ind,1:NP(ind+1)).^2);                                                              % Section area [m2]
            W(ind,1:NP(ind+1)) = pi.*(r1(ind,1:NP(ind+1)).^4)./(4.*r1(ind,1:NP(ind+1)));                                    % Resistant modulus [m3]
        elseif Section.type == 2    % Ellipsoidal
            J(ind,1:NP(ind+1)) = (pi.*(r2(ind,1:NP(ind+1)).*r1(ind,1:NP(ind+1)).^3))./4;                                    % Inertia [m4]
            A(ind,1:NP(ind+1)) = pi.*(r2(ind,1:NP(ind+1)).*r1(ind,1:NP(ind+1)));                                            % Section area [m2]
            W(ind,1:NP(ind+1)) = pi.*(r2(ind,1:NP(ind+1)).*r1(ind,1:NP(ind+1)).^3)./(4.*r2(ind,1:NP(ind+1)));               % Resistant modulus [m3]
        elseif Section.type == 3    % Square
            J(ind,1:NP(ind+1)) = (r1(ind,1:NP(ind+1)).^4)./12;                                                              % Inertia [m4]
            A(ind,1:NP(ind+1)) = r1(ind,1:NP(ind+1)).^2;                                                                    % Section area [m2]
            W(ind,1:NP(ind+1)) = (r1(ind,1:NP(ind+1)).^3)./6;                                                               % Resistant modulus [m3]
        elseif Section.type == 4    % Rectangular
            J(ind,1:NP(ind+1)) = r2(ind,1:NP(ind+1)).*(r1(ind,1:NP(ind+1)).^3)./12;                                         % Inertia [m4]
            A(ind,1:NP(ind+1)) = r2(ind,1:NP(ind+1)).*r1(ind,1:NP(ind+1));                                                  % Section area [m2]
            W(ind,1:NP(ind+1)) = r2(ind,1:NP(ind+1)).*(r1(ind,1:NP(ind+1)).^2)./6;                                          % Resistant modulus [m3]
        elseif Section.type == 5    % Hollowed circular
            J(ind,1:NP(ind+1)) = (pi.*(r1(ind,1:NP(ind+1)).^4-(r1(ind,1:NP(ind+1))-t(ind)).^4))./4;                         % Inertia [m4]
            A(ind,1:NP(ind+1)) = pi.*(r1(ind,1:NP(ind+1)).^2-(r1(ind,1:NP(ind+1))-t(ind)).^2);                              % Section area [m2]
            W(ind,1:NP(ind+1)) = (pi.*(r1(ind,1:NP(ind+1)).^4-(r1(ind,1:NP(ind+1))-t(ind)).^4))./(4.*r1(ind,1:NP(ind+1)));  % Resistant modulus [m3]
        elseif Section.type == 6    % Hollowed square
            J(ind,1:NP(ind+1)) =  (r1(ind,1:NP(ind+1)).^4 - (r1(ind,1:NP(ind+1))-2*t(ind)).^4)./12;                         % Inertia [m4]
            A(ind,1:NP(ind+1)) =  r1(ind,1:NP(ind+1)).^2 - ((r1(ind,1:NP(ind+1))-2*t(ind)).^2);                             % Section area [m2]
            W(ind,1:NP(ind+1)) =  (r1(ind,1:NP(ind+1)).^4 - (r1(ind,1:NP(ind+1))-2*t(ind)).^4)./(6.*r1(ind,1:NP(ind+1)));   % Resistant modulus [m3]
        elseif Section.type == 7    % Hollowed rectangle
            J(ind,1:NP(ind+1)) =  (r2(ind,1:NP(ind+1)).*(r1(ind,1:NP(ind+1)).^3) - (r2(ind,1:NP(ind+1))-2*t(ind)).*((r1(ind,1:NP(ind+1))-2*t(ind)).^3))./12;                        % Inertia [m4]
            A(ind,1:NP(ind+1)) =  r2(ind,1:NP(ind+1)).*r1(ind,1:NP(ind+1)) - ((r2(ind,1:NP(ind+1))-2*t(ind)).*(r1(ind,1:NP(ind+1))-2*t(ind)));                                      % Section area [m2]
            W(ind,1:NP(ind+1)) =  (r2(ind,1:NP(ind+1)).*(r1(ind,1:NP(ind+1)).^3) - (r2(ind,1:NP(ind+1))-2*t(ind)).*((r1(ind,1:NP(ind+1))-2*t(ind)).^3))./(6.*r1(ind,1:NP(ind+1)));  % Resistant modulus [m3]
        end
        
        % Smoothing of the inertia in case of linear variation of r1
        if Section.variation == 2
            J_fit(ind,:) = polyfit(s(ind,1:NP(ind+1)),J(ind,1:NP(ind+1)),4);
            J(ind,1:NP(ind+1)) = polyval(J_fit(ind,:),s(ind,1:NP(ind+1)));
        end
        
        % Inertia and its derivatives - mathematical form
        dJ(ind,1:NP(ind+1))  = interp1(s(ind,1:n(ind+1))+ds(ind)/2,diff(J(ind,1:NP(ind+1)))./ds(ind),s(ind,1:NP(ind+1)),'linear','extrap');
        d2J(ind,1:NP(ind+1)) = interp1(s(ind,1:n(ind+1))+ds(ind)/2,diff(dJ(ind,1:NP(ind+1)))./ds(ind),s(ind,1:NP(ind+1)),'linear','extrap');
        d3J(ind,1:NP(ind+1)) = interp1(s(ind,1:n(ind+1))+ds(ind)/2,diff(d2J(ind,1:NP(ind+1)))./ds(ind),s(ind,1:NP(ind+1)),'linear','extrap');
        
        
        %% Calculation of applied loads 1
    
        %alpha(ind,:) = atan(double(subs(dy,X,x_s(ind,:))));                                                    % Angle of inclination of each mesh segment [rad]
        q_sw(ind,1:NP(ind+1)) = -Material.gamma.*A(ind,1:NP(ind+1));                                            % Self-weight [kN/m]
        % q_ol = interp1(s(1:n)+ds/2,q0.*(x_s(2:n+1)-x_s(1:n))./ds,s,'pchip','extrap');                         % Overload [kN/m]
        % q_sw = 0;
        q_ol(ind,1:NP(ind+1)) = q0.*ones(1,NP(ind+1));
        % q_ol = [ zeros(1, floor(1*NP/3)) ones(1,ceil(1*NP/3)+1)*q0 zeros(1, floor(1*NP/3))];
        % q_ol = [ ones(1,ceil(NP/2))*q0 zeros(1, floor(NP/2))];
        
        
        Pn(ind,1:NP(ind+1)) = (q_sw(ind,1:NP(ind+1)) + q_ol(ind,1:NP(ind+1))).*cos(alpha(ind,1:NP(ind+1)));      % Load normal to the arch [kN/m]
        Pt(ind,1:NP(ind+1)) = (q_sw(ind,1:NP(ind+1)) + q_ol(ind,1:NP(ind+1))).*sin(alpha(ind,1:NP(ind+1)));      % Load tangent to the arch [kN/m]
        m_load(ind,1:NP(ind+1)) = 0;                                                                             % Bending moment [kNm/m]
        
        % Radial load
        % Pn = (q_sw).*cos(alpha) + q_ol;                                               % Load normal to the arch [kN/m]
        % Pt = (q_sw).*sin(alpha);                                                      % Load tangent to the arch [kN/m] 
        
        
        dPn(ind,1:NP(ind+1)) = interp1(s(ind,1:n(ind+1))+ds(ind)/2,diff(Pn(ind,1:NP(ind+1)))./ds(ind),s(ind,1:NP(ind+1)),'pchip','extrap');   % 1st derivative of normal load
        
        %% Calculation of the coefficients of the Equation of the elastic line 1
    
        % Partial coefficients
        a4(ind,1:NP(ind+1)) = R_s(ind,1:NP(ind+1));
        b3(ind,1:NP(ind+1)) = R_s(ind,1:NP(ind+1));
        b4(ind,1:NP(ind+1)) = 5.*dR_s(ind,1:NP(ind+1));
        c2(ind,1:NP(ind+1)) = R_s(ind,1:NP(ind+1));
        c3(ind,1:NP(ind+1)) = 4.*dR_s(ind,1:NP(ind+1));
        c4(ind,1:NP(ind+1)) = 10.*d2R_s(ind,1:NP(ind+1)) + 1./R_s(ind,1:NP(ind+1));
        d1(ind,1:NP(ind+1)) = R_s(ind,1:NP(ind+1));
        d2(ind,1:NP(ind+1)) = 3.*dR_s(ind,1:NP(ind+1));
        d3(ind,1:NP(ind+1)) = 6.*d2R_s(ind,1:NP(ind+1)) + 1./R_s(ind,1:NP(ind+1));
        d4(ind,1:NP(ind+1)) = 10.*d3R_s(ind,1:NP(ind+1)) - 4.*dR_s(ind,1:NP(ind+1)).*(1./((R_s(ind,1:NP(ind+1))).^2));
        e1(ind,1:NP(ind+1)) = 2.*dR_s(ind,1:NP(ind+1));
        e2(ind,1:NP(ind+1)) = 3.*d2R_s(ind,1:NP(ind+1)) + 1./R_s(ind,1:NP(ind+1));
        e3(ind,1:NP(ind+1)) = 4.*d3R_s(ind,1:NP(ind+1)) - 3.*dR_s(ind,1:NP(ind+1)).*(1./((R_s(ind,1:NP(ind+1))).^2));
        e4(ind,1:NP(ind+1)) = 5.*d4R_s(ind,1:NP(ind+1)) - 6.*d2R_s(ind,1:NP(ind+1)).*(1./((R_s(ind,1:NP(ind+1))).^2)) + 12.*((dR_s(ind,1:NP(ind+1))).^2).*(1./((R_s(ind,1:NP(ind+1))).^3));
        f1(ind,1:NP(ind+1)) = d2R_s(ind,1:NP(ind+1)) + 1./R_s(ind,1:NP(ind+1));
        f2(ind,1:NP(ind+1)) = d3R_s(ind,1:NP(ind+1)) - 2.*dR_s(ind,1:NP(ind+1)).*(1./((R_s(ind,1:NP(ind+1))).^2));
        f3(ind,1:NP(ind+1)) = d4R_s(ind,1:NP(ind+1)) - 3.*d2R_s(ind,1:NP(ind+1)).*(1./((R_s(ind,1:NP(ind+1))).^2)) + 6.*((dR_s(ind,1:NP(ind+1))).^2).*(1./((R_s(ind,1:NP(ind+1))).^3));
        f4(ind,1:NP(ind+1)) = d5R_s(ind,1:NP(ind+1)) - 4.*d3R_s(ind,1:NP(ind+1)).*(1./((R_s(ind,1:NP(ind+1))).^2)) + 24.*dR_s(ind,1:NP(ind+1)).*d2R_s(ind,1:NP(ind+1)).*(1./((R_s(ind,1:NP(ind+1))).^3)) - 24.*((dR_s(ind,1:NP(ind+1))).^3).*(1./((R_s(ind,1:NP(ind+1))).^4));
        g1(ind,1:NP(ind+1)) = -dR_s(ind,1:NP(ind+1)).*(1./((R_s(ind,1:NP(ind+1))).^2));
        g2(ind,1:NP(ind+1)) = -d2R_s(ind,1:NP(ind+1)).*(1./((R_s(ind,1:NP(ind+1))).^2)) + 2.*((dR_s(ind,1:NP(ind+1))).^2).*(1./((R_s(ind,1:NP(ind+1))).^3));
        g3(ind,1:NP(ind+1)) = -d3R_s(ind,1:NP(ind+1)).*(1./((R_s(ind,1:NP(ind+1))).^2)) + 6.*dR_s(ind,1:NP(ind+1)).*d2R_s(ind,1:NP(ind+1)).*(1./((R_s(ind,1:NP(ind+1))).^3)) - 6.*((dR_s(ind,1:NP(ind+1))).^3).*(1./((R_s(ind,1:NP(ind+1))).^4));
        g4(ind,1:NP(ind+1)) = -d4R_s(ind,1:NP(ind+1)).*(1./((R_s(ind,1:NP(ind+1))).^2)) + 8.*dR_s(ind,1:NP(ind+1)).*d3R_s(ind,1:NP(ind+1)).*(1./((R_s(ind,1:NP(ind+1))).^3)) + 6.*((d2R_s(ind,1:NP(ind+1))).^2).*(1./((R_s(ind,1:NP(ind+1))).^3)) -36.*((dR_s(ind,1:NP(ind+1))).^2).*d2R_s(ind,1:NP(ind+1)).*(1./((R_s(ind,1:NP(ind+1))).^4)) + 24.*((dR_s(ind,1:NP(ind+1))).^4).*(1./((R_s(ind,1:NP(ind+1))).^5));
        
        % Total coefficients
        a(ind,1:NP(ind+1)) = (R_s(ind,1:NP(ind+1)).*J(ind,1:NP(ind+1))).*a4(ind,1:NP(ind+1));
        b(ind,1:NP(ind+1)) = (3.*R_s(ind,1:NP(ind+1)).*dJ(ind,1:NP(ind+1)) + dR_s(ind,1:NP(ind+1)).*J(ind,1:NP(ind+1))).*b3(ind,1:NP(ind+1)) + (R_s(ind,1:NP(ind+1)).*J(ind,1:NP(ind+1))).*b4(ind,1:NP(ind+1));
        c(ind,1:NP(ind+1)) = (3.*R_s(ind,1:NP(ind+1)).*d2J(ind,1:NP(ind+1)) + 2.*dR_s(ind,1:NP(ind+1)).*dJ(ind,1:NP(ind+1)) + J(ind,1:NP(ind+1))./R_s(ind,1:NP(ind+1))).*c2(ind,1:NP(ind+1)) + (3.*R_s(ind,1:NP(ind+1)).*dJ(ind,1:NP(ind+1)) + dR_s(ind,1:NP(ind+1)).*J(ind,1:NP(ind+1))).*c3(ind,1:NP(ind+1)) + (R_s(ind,1:NP(ind+1)).*J(ind,1:NP(ind+1))).*c4(ind,1:NP(ind+1));
        d(ind,1:NP(ind+1)) = (R_s(ind,1:NP(ind+1)).*d3J(ind,1:NP(ind+1)) + dR_s(ind,1:NP(ind+1)).*d2J(ind,1:NP(ind+1)) + dJ(ind,1:NP(ind+1))./R_s(ind,1:NP(ind+1))).*d1(ind,1:NP(ind+1)) + (3.*R_s(ind,1:NP(ind+1)).*d2J(ind,1:NP(ind+1)) + 2.*dR_s(ind,1:NP(ind+1)).*dJ(ind,1:NP(ind+1)) + J(ind,1:NP(ind+1))./R_s(ind,1:NP(ind+1))).*d2(ind,1:NP(ind+1)) + (3.*R_s(ind,1:NP(ind+1)).*dJ(ind,1:NP(ind+1)) + dR_s(ind,1:NP(ind+1)).*J(ind,1:NP(ind+1))).*d3(ind,1:NP(ind+1)) + (R_s(ind,1:NP(ind+1)).*J(ind,1:NP(ind+1))).*d4(ind,1:NP(ind+1));
        e(ind,1:NP(ind+1)) = (R_s(ind,1:NP(ind+1)).*d3J(ind,1:NP(ind+1)) + dR_s(ind,1:NP(ind+1)).*d2J(ind,1:NP(ind+1)) + dJ(ind,1:NP(ind+1))./R_s(ind,1:NP(ind+1))).*e1(ind,1:NP(ind+1)) + (3.*R_s(ind,1:NP(ind+1)).*d2J(ind,1:NP(ind+1)) + 2.*dR_s(ind,1:NP(ind+1)).*dJ(ind,1:NP(ind+1)) + J(ind,1:NP(ind+1))./R_s(ind,1:NP(ind+1))).*e2(ind,1:NP(ind+1)) + (3.*R_s(ind,1:NP(ind+1)).*dJ(ind,1:NP(ind+1)) + dR_s(ind,1:NP(ind+1)).*J(ind,1:NP(ind+1))).*e3(ind,1:NP(ind+1)) + (R_s(ind,1:NP(ind+1)).*J(ind,1:NP(ind+1))).*e4(ind,1:NP(ind+1));
        f(ind,1:NP(ind+1)) = (R_s(ind,1:NP(ind+1)).*d3J(ind,1:NP(ind+1)) + dR_s(ind,1:NP(ind+1)).*d2J(ind,1:NP(ind+1)) + dJ(ind,1:NP(ind+1))./R_s(ind,1:NP(ind+1))).*f1(ind,1:NP(ind+1)) + (3.*R_s(ind,1:NP(ind+1)).*d2J(ind,1:NP(ind+1)) + 2.*dR_s(ind,1:NP(ind+1)).*dJ(ind,1:NP(ind+1)) + J(ind,1:NP(ind+1))./R_s(ind,1:NP(ind+1))).*f2(ind,1:NP(ind+1)) + (3.*R_s(ind,1:NP(ind+1)).*dJ(ind,1:NP(ind+1)) + dR_s(ind,1:NP(ind+1)).*J(ind,1:NP(ind+1))).*f3(ind,1:NP(ind+1)) + (R_s(ind,1:NP(ind+1)).*J(ind,1:NP(ind+1))).*f4(ind,1:NP(ind+1));
        g(ind,1:NP(ind+1)) = (R_s(ind,1:NP(ind+1)).*d3J(ind,1:NP(ind+1)) + dR_s(ind,1:NP(ind+1)).*d2J(ind,1:NP(ind+1)) + dJ(ind,1:NP(ind+1))./R_s(ind,1:NP(ind+1))).*g1(ind,1:NP(ind+1)) + (3.*R_s(ind,1:NP(ind+1)).*d2J(ind,1:NP(ind+1)) + 2.*dR_s(ind,1:NP(ind+1)).*dJ(ind,1:NP(ind+1)) + J(ind,1:NP(ind+1))./R_s(ind,1:NP(ind+1))).*g2(ind,1:NP(ind+1)) + (3.*R_s(ind,1:NP(ind+1)).*dJ(ind,1:NP(ind+1)) + dR_s(ind,1:NP(ind+1)).*J(ind,1:NP(ind+1))).*g3(ind,1:NP(ind+1)) + (R_s(ind,1:NP(ind+1)).*J(ind,1:NP(ind+1))).*g4(ind,1:NP(ind+1));
        h(ind,1:NP(ind+1)) = -(Pt(ind,1:NP(ind+1)) + dR_s(ind,1:NP(ind+1)).*Pn(ind,1:NP(ind+1)) + R_s(ind,1:NP(ind+1)).*dPn(ind,1:NP(ind+1)) + m_load(ind,1:NP(ind+1))./R_s(ind,1:NP(ind+1)))./Material.E;
        
        %Coefficients for M,N,T to solver conditions
        
        % Bending moment coefficients
        aM(ind,1:NP(ind+1)) = R_s(ind,1:NP(ind+1));
        bM(ind,1:NP(ind+1)) = 2.*dR_s(ind,1:NP(ind+1));
        cM(ind,1:NP(ind+1)) = (d2R_s(ind,1:NP(ind+1)) + 1./R_s(ind,1:NP(ind+1)));
        dM(ind,1:NP(ind+1)) = -dR_s(ind,1:NP(ind+1))./(R_s(ind,1:NP(ind+1))).^2;
        
        % Shear force coefficients
        aT(ind,1:NP(ind+1)) = R_s(ind,1:NP(ind+1)).*J(ind,1:NP(ind+1));
        bT(ind,1:NP(ind+1)) = R_s(ind,1:NP(ind+1)).*dJ(ind,1:NP(ind+1)) + 3.*dR_s(ind,1:NP(ind+1)).*J(ind,1:NP(ind+1));
        cT(ind,1:NP(ind+1)) = 2.*dR_s(ind,1:NP(ind+1)).*dJ(ind,1:NP(ind+1)) + 3.*d2R_s(ind,1:NP(ind+1)).*J(ind,1:NP(ind+1)) + J(ind,1:NP(ind+1))./R_s(ind,1:NP(ind+1));
        dT(ind,1:NP(ind+1)) = dJ(ind,1:NP(ind+1)).*d2R_s(ind,1:NP(ind+1)) + dJ(ind,1:NP(ind+1))./R_s(ind,1:NP(ind+1)) + d3R_s(ind,1:NP(ind+1)).*J(ind,1:NP(ind+1)) - 2.*J(ind,1:NP(ind+1)).*dR_s(ind,1:NP(ind+1))./(R_s(ind,1:NP(ind+1))).^2;
        eT(ind,1:NP(ind+1)) = -(dJ(ind,1:NP(ind+1)).*(dR_s(ind,1:NP(ind+1))./(R_s(ind,1:NP(ind+1))).^2) + J(ind,1:NP(ind+1)).*(d2R_s(ind,1:NP(ind+1))./(R_s(ind,1:NP(ind+1))).^2) -2.*J(ind,1:NP(ind+1)).*((dR_s(ind,1:NP(ind+1)).^2)./(R_s(ind,1:NP(ind+1))).^3));
    
        % Normal force coefficients
        aN(ind,1:NP(ind+1)) = R_s(ind,1:NP(ind+1)).*J(ind,1:NP(ind+1));
        bN(ind,1:NP(ind+1)) = 2.*R_s(ind,1:NP(ind+1)).*dJ(ind,1:NP(ind+1)) + 4.*dR_s(ind,1:NP(ind+1)).*J(ind,1:NP(ind+1));
        cN(ind,1:NP(ind+1)) = d2J(ind,1:NP(ind+1)).*R_s(ind,1:NP(ind+1)) + 6.*dJ(ind,1:NP(ind+1)).*dR_s(ind,1:NP(ind+1)) + 6.*J(ind,1:NP(ind+1)).*d2R_s(ind,1:NP(ind+1)) + J(ind,1:NP(ind+1))./R_s(ind,1:NP(ind+1));
        dN(ind,1:NP(ind+1)) = 2.*d2J(ind,1:NP(ind+1)).*dR_s(ind,1:NP(ind+1)) + 6.*dJ(ind,1:NP(ind+1)).*d2R_s(ind,1:NP(ind+1)) + 2.*dJ(ind,1:NP(ind+1))./R_s(ind,1:NP(ind+1)) + 4.*J(ind,1:NP(ind+1)).*d3R_s(ind,1:NP(ind+1)) - 3.*J(ind,1:NP(ind+1)).*dR_s(ind,1:NP(ind+1))./(R_s(ind,1:NP(ind+1))).^2;
        eN(ind,1:NP(ind+1)) = d2J(ind,1:NP(ind+1)).*d2R_s(ind,1:NP(ind+1)) + 2.*dJ(ind,1:NP(ind+1)).*d3R_s(ind,1:NP(ind+1)) + J(ind,1:NP(ind+1)).*d4R_s(ind,1:NP(ind+1)) + d2J(ind,1:NP(ind+1))./R_s(ind,1:NP(ind+1)) - 4.*dJ(ind,1:NP(ind+1)).*dR_s(ind,1:NP(ind+1))./(R_s(ind,1:NP(ind+1)).^2) - 3.*J(ind,1:NP(ind+1)).*d2R_s(ind,1:NP(ind+1))./(R_s(ind,1:NP(ind+1)).^2) + 6.*J(ind,1:NP(ind+1)).*(dR_s(ind,1:NP(ind+1)).^2)./(R_s(ind,1:NP(ind+1)).^3);
        fN(ind,1:NP(ind+1)) = -d2J(ind,1:NP(ind+1)).*dR_s(ind,1:NP(ind+1))./(R_s(ind,1:NP(ind+1)).^2) - 2.*dJ(ind,1:NP(ind+1)).*d2R_s(ind,1:NP(ind+1))./(R_s(ind,1:NP(ind+1)).^2) - J(ind,1:NP(ind+1)).*d3R_s(ind,1:NP(ind+1))./(R_s(ind,1:NP(ind+1)).^2) + 8.*dJ(ind,1:NP(ind+1)).*(dR_s(ind,1:NP(ind+1)).^2)./(R_s(ind,1:NP(ind+1)).^3) + 2.*J(ind,1:NP(ind+1)).*dR_s(ind,1:NP(ind+1)).*d2R_s(ind,1:NP(ind+1))./(R_s(ind,1:NP(ind+1)).^3) - 6.*J(ind,1:NP(ind+1)).*(dR_s(ind,1:NP(ind+1)).^3)./(R_s(ind,1:NP(ind+1)).^4);
    end
    
    %% Definizione della rigidezza delle molle rotazionali
    % k = EI/L, intendendo L in questo caso come la lunghezza dell'elemento
    % lungo l'ascissa curvilinea.
    
    %Allocazione di memoria
    k_in = zeros(1,length(x_point));
    k_fin = zeros(1,length(x_point));
    L_d_in = zeros(1,length(x_point));
    L_d_fin = zeros(1,length(x_point));
    
    % Tiene conto del fatto che possono esservi delle molle alle estremità
    if molla_in == 1
        k_in = zeros(1,length(x_point)+1);
        L_d_in = zeros(1,length(x_point)+1);
    end
    if molla_fin == NP_1
        k_fin = zeros(1,length(x_point)+1);
        L_d_fin = zeros(1,length(x_point)+1);
    end
    
    if Restraint == 1
    phi_d = 0.016*(pi);                                                    % [rad]
    elseif Restraint == 2
    phi_d = 0.067*(pi);                                                    % [rad]
    end

    % Definisco il vettore delle rigidezze decrementate

    %Presenza di molla alla base
    if molla_in == 1
        L_d_in(1) = phi_d.*(R_s(1,1));                                     % lunghezza tratto danneggiato [m]
        k_in(1) = Material_E*J(1,1)/L_d_in(1);                             % Rigidezza iniziale [kNm/rad]
       for ind = 1:length(x_point)
        L_d_in(ind+1) = phi_d*(R_s(ind+1,1));                              % lunghezza tratto danneggiato [m]
        k_in(ind+1) = Material_E*J(ind+1,1)/L_d_in(ind+1);                 % Rigidezza iniziale [kNm/rad]
       end
    end
    % Presenza di molla in mezzeria
    if molla_fin == NP_1
        L_d_fin(end) = phi_d.*(R_s(end,NP(end)));                          % lunghezza tratto danneggiato [m]
        k_fin(end) = Material_E*J(end,NP(end))/L_d_fin(end);               % Rigidezza finale [kNm/rad]
        for ind = 1:length(x_point)
        L_d_fin(ind) = phi_d*(R_s(ind,NP(ind+1)));                         % lunghezza tratto danneggiato [m]
        k_fin(ind) = Material_E*J(ind,NP(ind+1))/L_d_fin(ind);             % Rigidezza finale [kNm/rad]
       end
    end
    %Assenza di molla in mezzeria
    if molla_fin ~= NP_1
        for ind = 1:length(x_point)
        L_d_fin(ind) = phi_d*(R_s(ind,NP(ind+1)));                         % lunghezza tratto danneggiato [m]
        k_fin(ind) = Material_E*J(ind,NP(ind+1))/L_d_fin(ind);             % Rigidezza finale [kNm/rad]
       end
    end
    %Assenza di molla alla base
    if molla_in ~= 1
        for ind = 1:length(x_point)
        L_d_in(ind) = phi_d*(R_s(ind+1,1));                                % lunghezza tratto danneggiato [m]
        k_in(ind) = Material_E*J(ind+1,1)/L_d_in(ind);                     % Rigidezza iniziale [kNm/rad]
       end
    end
    

    %In questo modo, se il vettore dei punti con cui individuo il
    %danneggiamento non si è aggiornato perchè non vi sono ulteriori
    %danneggiamenti, allora si diminuisce la rigidezza esistente. Nel caso
    %in esame è prorpio il funzionamento del codice che tiene conto già
    %inizialmente di una configurazione di danneggiamento
    
    % Determino il decremento iniziale da cui parte ogni molla. Posso
    % decidere il valore iniziale
    if contatore == 1
        flag = 0;

        percentuale_in = [0.70 0.70 0.70 0.70];

        %Assenza di molla alla base sx
        if molla_in ~= 1
            % Assenza di molla alla base dx
            if molla_fin ~= NP_1
                percentuale_fin = percentuale_in;
            % Presenza di molla alla base dx
            elseif molla_fin == NP_1
                percentuale_fin = [percentuale_in perc_dx];
            end
        %  Presenza di molla alla base sx   
        elseif molla_in == 1
            % Assenza di molla alla base dx
            if molla_fin ~= NP_1
                percentuale_fin = [percentuale_in(2:end)];
            % Presenza di molla alla base dx    
            elseif molla_fin == NP_1
                percentuale_fin = [percentuale_in(2:end) perc_dx];
            end
        end
    
        %Allocazione di memoria
        J_d_in = zeros(1,length(x_point));
        J_d_fin = zeros(1,length(x_point));
        delta_in = zeros(1,length(x_point));
        delta_fin = zeros(1,length(x_point));
        decremento_in = zeros(1,length(x_point));
        decremento_fin = zeros(1,length(x_point));
    
    
        % Tiene conto del fatto che possono esservi delle molle alle estremità
        if molla_in == 1
            J_d_in = zeros(1,length(x_point)+1);
            delta_in = zeros(1,length(x_point)+1);
            decremento_in = zeros(1,length(x_point)+1);
        end
        if molla_fin == NP_1
            J_d_fin = zeros(1,length(x_point)+1);
            delta_fin = zeros(1,length(x_point)+1);
            decremento_fin = zeros(1,length(x_point)+1);
        end
    
        % Presenza di molla alla base
        if molla_in == 1
        J_d_in(1) = J(1,1).*percentuale_in(1);
        delta_in(1) = (J(1,1) - J_d_in(1))./J(1,1);
        decremento_in(1) = (1-delta_in(1))./delta_in(1);
            for ind = 1:length(x_point)
                J_d_in(ind+1) = J(ind+1,1).*percentuale_in(ind+1);
                delta_in(ind+1) = ((J(ind+1,1) - J_d_in(ind+1))./J(ind+1,1));
                decremento_in(ind+1) = (1-delta_in(ind+1))./delta_in(ind+1);
            end
        end
        %Presenza di molla in mezzeria
        if molla_fin == NP_1
            J_d_fin(end) = J(end,NP(end)).*percentuale_fin(end);
            delta_fin(end) = (J(end,NP(end)) - J_d_fin(end))./J(end,NP(end));
            decremento_fin(end) = (1-delta_fin(end))./delta_fin(end);
            for ind = 1:length(x_point)
                J_d_fin(ind) = J(ind,NP(ind+1)).*percentuale_fin(ind);
                delta_fin(ind) = ((J(ind,NP(ind+1)) - J_d_fin(ind))./J(ind,NP(ind+1)));
                decremento_fin(ind) = (1-delta_fin(ind))./delta_fin(ind);
            end
        end
        %Assenza di molla in mezzeria
        if molla_fin ~= NP_1
            for ind = 1:length(x_point)
                J_d_fin(ind) = J(ind,NP(ind+1)).*percentuale_fin(ind);
                delta_fin(ind) = ((J(ind,NP(ind+1)) - J_d_fin(ind))./J(ind,NP(ind+1)));
                decremento_fin(ind) = (1-delta_fin(ind))./delta_fin(ind);
            end
        end
        %Assenza di molla alla base
        if molla_in ~= 1
            for ind = 1:length(x_point)
                J_d_in(ind) = J(ind+1,1).*percentuale_in(ind);
                delta_in(ind) = ((J(ind+1,1) - J_d_in(ind))./J(ind+1,1));
                decremento_in(ind) = (1-delta_in(ind))./delta_in(ind);
            end
        end
    end

    if flag == 1
        k_in = k_in(1,:).*decremento_in(1,:);
        k_fin = k_fin(1,:).*decremento_fin(1,:);
    end
    
    %% Numero di punti progressivo
    
    % Allocazione di memoria
    NP_tot = zeros(1,length(x_point)+2);
    
    for ind = 1:length(x_point)+1
        NP_tot(1) = NP(1);
        NP_tot(ind+1) = NP(ind+1)+NP_tot(ind);
    end
    %% Displacements - Finite Differences solution with system
    
    % Allocazione di memoria
    u_0 = zeros(length(x_point)+1,NP_max+6);
    ce1 = zeros(1,NP_tot(end));
    ce2 = zeros(1,NP_tot(end));
    ce3 = zeros(1,NP_tot(end));
    ce4 = zeros(1,NP_tot(end));
    ce5 = zeros(1,NP_tot(end));
    ce6 = zeros(1,NP_tot(end));
    ce7 = zeros(1,NP_tot(end));
    
    % Il +6 serve per tenere conto della presenza dei punti fittizi per ogni
    % singolo arco e relativamente alle righe per considerare i vincoli esterni
    matA = zeros(NP_tot(end)+(length(x_point)+1)*6,NP_tot(end)+(length(x_point)+1)*6);
    B = zeros(1,NP_tot(end)+(length(x_point)+1)*6);
    
        
    %definisco il numero di equazioni che corrispondono alle equazioni al contorno e quindi 6 per ogni arco
    
    for ind = 1:length(x_point)+1 
        n_cc = ind*6;                       
        n_cc1 = (length(x_point)+1)*6-ind*6;
        
        %% Definizione delle equazioni per tutti i punti
    
        for i = 1:NP(ind+1)
        
            ce1(NP_tot(ind)+i) = -(1/60)*f(ind,i)/ds(ind) + (1/90)*e(ind,i)/ds(ind)^2 + (1/8)*d(ind,i)/ds(ind)^3 - (1/6)*c(ind,i)/ds(ind)^4 - (1/2)*b(ind,i)/ds(ind)^5 + 1*a(ind,i)/ds(ind)^6;
            ce7(NP_tot(ind)+i) = +(1/60)*f(ind,i)/ds(ind) + (1/90)*e(ind,i)/ds(ind)^2 - (1/8)*d(ind,i)/ds(ind)^3 - (1/6)*c(ind,i)/ds(ind)^4 + (1/2)*b(ind,i)/ds(ind)^5 + 1*a(ind,i)/ds(ind)^6;
            ce2(NP_tot(ind)+i) = +(3/20)*f(ind,i)/ds(ind) - (3/20)*e(ind,i)/ds(ind)^2 - 1*d(ind,i)/ds(ind)^3 + 2*c(ind,i)/ds(ind)^4 + 2*b(ind,i)/ds(ind)^5 - 6*a(ind,i)/ds(ind)^6;
            ce6(NP_tot(ind)+i) = -(3/20)*f(ind,i)/ds(ind) - (3/20)*e(ind,i)/ds(ind)^2 + 1*d(ind,i)/ds(ind)^3 + 2*c(ind,i)/ds(ind)^4 - 2*b(ind,i)/ds(ind)^5 - 6*a(ind,i)/ds(ind)^6;
            ce3(NP_tot(ind)+i) = -(3/4)*f(ind,i)/ds(ind) + (3/2)*e(ind,i)/ds(ind)^2 + (13/8)*d(ind,i)/ds(ind)^3 - (13/2)*c(ind,i)/ds(ind)^4 - (5/2)*b(ind,i)/ds(ind)^5 + 15*a(ind,i)/ds(ind)^6;
            ce5(NP_tot(ind)+i) = +(3/4)*f(ind,i)/ds(ind) + (3/2)*e(ind,i)/ds(ind)^2 - (13/8)*d(ind,i)/ds(ind)^3 - (13/2)*c(ind,i)/ds(ind)^4 + (5/2)*b(ind,i)/ds(ind)^5 + 15*a(ind,i)/ds(ind)^6;
            ce4(NP_tot(ind)+i) = -(49/18)*e(ind,i)/ds(ind)^2 + (28/3)*c(ind,i)/ds(ind)^4 - 20*a(ind,i)/ds(ind)^6 + 1*g(ind,i);
    
            %matA(NP_tot(ind)+i,:) = [zeros(1,NP_tot(ind)) zeros(1,n_cc-6) zeros(1,i-1) ce1(NP_tot(ind)+i) ce2(NP_tot(ind)+i) ce3(NP_tot(ind)+i) ce4(NP_tot(ind)+i) ce5(NP_tot(ind)+i) ce6(NP_tot(ind)+i) ce7(NP_tot(ind)+i) zeros(1,(NP(ind+1)-i)) zeros(1,NP_tot(end)-NP_tot(ind+1)) zeros(1,n_cc1)];
            matA(NP_tot(ind)+i,i+n_cc-6+NP_tot(ind):i+6+n_cc-6+NP_tot(ind)) = [ce1(NP_tot(ind)+i) ce2(NP_tot(ind)+i) ce3(NP_tot(ind)+i) ce4(NP_tot(ind)+i) ce5(NP_tot(ind)+i) ce6(NP_tot(ind)+i) ce7(NP_tot(ind)+i)];
        
            B(NP_tot(ind)+i) = -h(ind,i);
        end
    end
    %% Vincolo intermedio interno (incastro)
    
    if contatore == 1
    Restraint_intermediate = input("Inserire la tipologia di vincolo intermedio: \n Incastro arch = 1\n Cerniera = 2\n Molla = 3\n");
    end
    
    for ind = 1:length(x_point)
    
        n_cc = ind*6;
        ind1 = length(x_point)+1-ind;
        n_cc1 = (length(x_point)+1)*6-ind*6;
    
        v_res_1 = [(1/60)*R_s(ind,NP(ind+1))/ds(ind) (-3/20)*R_s(ind,NP(ind+1))/ds(ind) (3/4)*R_s(ind,NP(ind+1))/ds(ind)  0 (-3/4)*R_s(ind,NP(ind+1))/ds(ind) (3/20)*R_s(ind,NP(ind+1))/ds(ind) (-1/60)*R_s(ind,NP(ind+1))/ds(ind)];
        v_res_2 = [(1/60)*R_s(ind+1,1)/ds(ind+1) (-3/20)*R_s(ind+1,1)/ds(ind+1) (3/4)*R_s(ind+1,1)/ds(ind+1) 0 (-3/4)*R_s(ind+1,1)/ds(ind+1) (3/20)*R_s(ind+1,1)/ds(ind+1) (-1/60)*R_s(ind+1,1)/ds(ind+1)];
        
        phi_res_1 = [((-1/90)*R_s(ind,NP(ind+1))/ds(ind)^2 + (1/60)*dR_s(ind,NP(ind+1))/ds(ind))...
                    ((3/20)*R_s(ind,NP(ind+1))/ds(ind)^2 - (3/20)*dR_s(ind,NP(ind+1))/ds(ind))...
                    ((-3/2)*R_s(ind,NP(ind+1))/ds(ind)^2 + (3/4)*dR_s(ind,NP(ind+1))/ds(ind))...
                    ((49/18)*R_s(ind,NP(ind+1))/ds(ind)^2 - 1/R_s(ind,NP(ind+1)))...
                    ((-3/2)*R_s(ind,NP(ind+1))/ds(ind)^2 - (3/4)*dR_s(ind,NP(ind+1))/ds(ind))...
                    ((3/20)*R_s(ind,NP(ind+1))/ds(ind)^2 + (3/20)*dR_s(ind,NP(ind+1))/ds(ind))...
                    ((-1/90)*R_s(ind,NP(ind+1))/ds(ind)^2 - (1/60)*dR_s(ind,NP(ind+1))/ds(ind))];
        
        phi_res_2 = [((-1/90)*R_s(ind+1,1)/ds(ind+1)^2 + (1/60)*dR_s(ind+1,1)/ds(ind+1))...
                    ((3/20)*R_s(ind+1,1)/ds(ind+1)^2 - (3/20)*dR_s(ind+1,1)/ds(ind+1))...
                    ((-3/2)*R_s(ind+1,1)/ds(ind+1)^2 + (3/4)*dR_s(ind+1,1)/ds(ind+1))...
                    ((49/18)*R_s(ind+1,1)/ds(ind+1)^2 - 1/R_s(ind+1,1))...
                    ((-3/2)*R_s(ind+1,1)/ds(ind+1)^2 - (3/4)*dR_s(ind+1,1)/ds(ind+1))...
                    ((3/20)*R_s(ind+1,1)/ds(ind+1)^2 + (3/20)*dR_s(ind+1,1)/ds(ind+1))...
                    ((-1/90)*R_s(ind+1,1)/ds(ind+1)^2 - (1/60)*dR_s(ind+1,1)/ds(ind+1))];
        
        M_res_1 = Material_E.*[J(ind,NP(ind+1)).*((1/8)*aM(ind,NP(ind+1))/ds(ind)^3 + (1/90)*bM(ind,NP(ind+1))/ds(ind)^2 + (-1/60)*cM(ind,NP(ind+1))/ds(ind))...
                                J(ind,NP(ind+1)).*(-1*aM(ind,NP(ind+1))/ds(ind)^3 + (-3/20)*bM(ind,NP(ind+1))/ds(ind)^2 + (3/20)*cM(ind,NP(ind+1))/ds(ind))...
                                J(ind,NP(ind+1)).*((13/8)*aM(ind,NP(ind+1))/ds(ind)^3 + (3/2)*bM(ind,NP(ind+1))/ds(ind)^2 + (-3/4)*cM(ind,NP(ind+1))/ds(ind))...
                                J(ind,NP(ind+1)).*((-49/18)*bM(ind,NP(ind+1))/ds(ind)^2 + dM(ind,NP(ind+1)))...
                                J(ind,NP(ind+1)).*((-13/8)*aM(ind,NP(ind+1))/ds(ind)^3 + (3/2)*bM(ind,NP(ind+1))/ds(ind)^2 + (3/4)*cM(ind,NP(ind+1))/ds(ind))...
                                J(ind,NP(ind+1)).*(1*aM(ind,NP(ind+1))/ds(ind)^3 + (-3/20)*bM(ind,NP(ind+1))/ds(ind)^2 + (-3/20)*cM(ind,NP(ind+1))/ds(ind))...
                                J(ind,NP(ind+1)).*((-1/8)*aM(ind,NP(ind+1))/ds(ind)^3 + (1/90)*bM(ind,NP(ind+1))/ds(ind)^2 + (1/60)*cM(ind,NP(ind+1))/ds(ind))]; %u(i+3)
        
        M_res_2 = Material_E.*[J(ind+1,1).*((1/8)*aM(ind+1,1)/ds(ind+1)^3 + (1/90)*bM(ind+1,1)/ds(ind+1)^2 + (-1/60)*cM(ind+1,1)/ds(ind+1))...
                                J(ind+1,1).*(-1*aM(ind+1,1)/ds(ind+1)^3 + (-3/20)*bM(ind+1,1)/ds(ind+1)^2 + (3/20)*cM(ind+1,1)/ds(ind+1))...
                                J(ind+1,1).*((13/8)*aM(ind+1,1)/ds(ind+1)^3 + (3/2)*bM(ind+1,1)/ds(ind+1)^2 + (-3/4)*cM(ind+1,1)/ds(ind+1))...
                                J(ind+1,1).*((-49/18)*bM(ind+1,1)/ds(ind+1)^2 + dM(ind+1,1))...
                                J(ind+1,1).*((-13/8)*aM(ind+1,1)/ds(ind+1)^3 + (3/2)*bM(ind+1,1)/ds(ind+1)^2 + (3/4)*cM(ind+1,1)/ds(ind+1))...
                                J(ind+1,1).*(1*aM(ind+1,1)/ds(ind+1)^3 + (-3/20)*bM(ind+1,1)/ds(ind+1)^2 + (-3/20)*cM(ind+1,1)/ds(ind+1))...
                                J(ind+1,1).*((-1/8)*aM(ind+1,1)/ds(ind+1)^3 + (1/90)*bM(ind+1,1)/ds(ind+1)^2 + (1/60)*cM(ind+1,1)/ds(ind+1))];
        
        T_res_1 = Material_E.*[((-1/6)*aT(ind,NP(ind+1))/ds(ind)^4 + (1/8)*bT(ind,NP(ind+1))/ds(ind)^3 + (1/90)*cT(ind,NP(ind+1))/ds(ind)^2 + (-1/60)*dT(ind,NP(ind+1))/ds(ind))...
                                ((2)*aT(ind,NP(ind+1))/ds(ind)^4 - 1*bT(ind,NP(ind+1))/ds(ind)^3 + (-3/20)*cT(ind,NP(ind+1))/ds(ind)^2 + (3/20)*dT(ind,NP(ind+1))/ds(ind))...
                                ((-13/2)*aT(ind,NP(ind+1))/ds(ind)^4 + (13/8)*bT(ind,NP(ind+1))/ds(ind)^3 + (3/2)*cT(ind,NP(ind+1))/ds(ind)^2 + (-3/4)*dT(ind,NP(ind+1))/ds(ind))...
                                ((28/3)*aT(ind,NP(ind+1))/ds(ind)^4 + (-49/18)*cT(ind,NP(ind+1))/ds(ind)^2 + eT(ind,NP(ind+1)))...
                                ((-13/2)*aT(ind,NP(ind+1))/ds(ind)^4 + (-13/8)*bT(ind,NP(ind+1))/ds(ind)^3 + (3/2)*cT(ind,NP(ind+1))/ds(ind)^2 + (3/4)*dT(ind,NP(ind+1))/ds(ind))...
                                ((2)*aT(ind,NP(ind+1))/ds(ind)^4 + 1*bT(ind,NP(ind+1))/ds(ind)^3 + (-3/20)*cT(ind,NP(ind+1))/ds(ind)^2 + (-3/20)*dT(ind,NP(ind+1))/ds(ind))...
                                ((-1/6)*aT(ind,NP(ind+1))/ds(ind)^4 + (-1/8)*bT(ind,NP(ind+1))/ds(ind)^3 + (1/90)*cT(ind,NP(ind+1))/ds(ind)^2 + (1/60)*dT(ind,NP(ind+1))/ds(ind))];
        
        T_res_2 = Material_E.*[((-1/6)*aT(ind+1,1)/ds(ind+1)^4 + (1/8)*bT(ind+1,1)/ds(ind+1)^3 + (1/90)*cT(ind+1,1)/ds(ind+1)^2 + (-1/60)*dT(ind+1,1)/ds(ind+1))...
                                ((2)*aT(ind+1,1)/ds(ind+1)^4 - 1*bT(ind+1,1)/ds(ind+1)^3 + (-3/20)*cT(ind+1,1)/ds(ind+1)^2 + (3/20)*dT(ind+1,1)/ds(ind+1))...
                                ((-13/2)*aT(ind+1,1)/ds(ind+1)^4 + (13/8)*bT(ind+1,1)/ds(ind+1)^3 + (3/2)*cT(ind+1,1)/ds(ind+1)^2 + (-3/4)*dT(ind+1,1)/ds(ind+1))...
                                ((28/3)*aT(ind+1,1)/ds(ind+1)^4 + (-49/18)*cT(ind+1,1)/ds(ind+1)^2 + eT(ind+1,1))...
                                ((-13/2)*aT(ind+1,1)/ds(ind+1)^4 + (-13/8)*bT(ind+1,1)/ds(ind+1)^3 + (3/2)*cT(ind+1,1)/ds(ind+1)^2 + (3/4)*dT(ind+1,1)/ds(ind+1))...
                                ((2)*aT(ind+1,1)/ds(ind+1)^4 + 1*bT(ind+1,1)/ds(ind+1)^3 + (-3/20)*cT(ind+1,1)/ds(ind+1)^2 + (-3/20)*dT(ind+1,1)/ds(ind+1))...
                                ((-1/6)*aT(ind+1,1)/ds(ind+1)^4 + (-1/8)*bT(ind+1,1)/ds(ind+1)^3 + (1/90)*cT(ind+1,1)/ds(ind+1)^2 + (1/60)*dT(ind+1,1)/ds(ind+1))];
        
        N_res_1 =  [((R_s(ind,NP(ind+1))*Pn(ind,NP(ind+1)) - Material_E.*R_s(ind,NP(ind+1)).*((-1/2)*aN(ind,NP(ind+1))/ds(ind)^5 + (-1/6)*bN(ind,NP(ind+1))/ds(ind)^4 + (1/8)*cN(ind,NP(ind+1))/ds(ind)^3 + (1/90)*dN(ind,NP(ind+1))/ds(ind)^2 + (-1/60)*eN(ind,NP(ind+1))/ds(ind))))... %u(i-3)
                    ((R_s(ind,NP(ind+1))*Pn(ind,NP(ind+1)) - Material_E.*R_s(ind,NP(ind+1)).*((2)*aN(ind,NP(ind+1))/ds(ind)^5 + (2)*bN(ind,NP(ind+1))/ds(ind)^4 - 1*cN(ind,NP(ind+1))/ds(ind)^3 + (-3/20)*dN(ind,NP(ind+1))/ds(ind)^2 + (3/20)*eN(ind,NP(ind+1))/ds(ind))))... %u(i-2)
                    ((R_s(ind,NP(ind+1))*Pn(ind,NP(ind+1)) - Material_E.*R_s(ind,NP(ind+1)).*((-5/2)*aN(ind,NP(ind+1))/ds(ind)^5 + (-13/2)*bN(ind,NP(ind+1))/ds(ind)^4 + (13/8)*cN(ind,NP(ind+1))/ds(ind)^3 + (3/2)*dN(ind,NP(ind+1))/ds(ind)^2 + (-3/4)*eN(ind,NP(ind+1))/ds(ind))))... %u(i-1)
                    ((R_s(ind,NP(ind+1))*Pn(ind,NP(ind+1)) - Material_E.*R_s(ind,NP(ind+1)).*((28/3)*bN(ind,NP(ind+1))/ds(ind)^4 + (-49/18)*dN(ind,NP(ind+1))/ds(ind)^2 + fN(ind,NP(ind+1)))))... %u(i)
                    ((R_s(ind,NP(ind+1))*Pn(ind,NP(ind+1)) - Material_E.*R_s(ind,NP(ind+1)).*((5/2)*aN(ind,NP(ind+1))/ds(ind)^5 + (-13/2)*bN(ind,NP(ind+1))/ds(ind)^4 + (-13/8)*cN(ind,NP(ind+1))/ds(ind)^3 + (3/2)*dN(ind,NP(ind+1))/ds(ind)^2 + (3/4)*eN(ind,NP(ind+1))/ds(ind))))... %u(i+1)
                    ((R_s(ind,NP(ind+1))*Pn(ind,NP(ind+1)) - Material_E.*R_s(ind,NP(ind+1)).*((-2)*aN(ind,NP(ind+1))/ds(ind)^5 + (2)*bN(ind,NP(ind+1))/ds(ind)^4 + 1*cN(ind,NP(ind+1))/ds(ind)^3 + (-3/20)*dN(ind,NP(ind+1))/ds(ind)^2 + (-3/20)*eN(ind,NP(ind+1))/ds(ind))))... %u(i+2)
                    ((R_s(ind,NP(ind+1))*Pn(ind,NP(ind+1)) - Material_E.*R_s(ind,NP(ind+1)).*((1/2)*aN(ind,NP(ind+1))/ds(ind)^5 + (-1/6)*bN(ind,NP(ind+1))/ds(ind)^4 + (-1/8)*cN(ind,NP(ind+1))/ds(ind)^3 + (1/90)*dN(ind,NP(ind+1))/ds(ind)^2 + (1/60)*eN(ind,NP(ind+1))/ds(ind))))]; %u(i+3)
        
        N_res_2 = [((R_s(ind+1,1)*Pn(ind+1,1) - Material_E.*R_s(ind+1,1).*((-1/2)*aN(ind+1,1)/ds(ind+1)^5 + (-1/6)*bN(ind+1,1)/ds(ind+1)^4 + (1/8)*cN(ind+1,1)/ds(ind+1)^3 + (1/90)*dN(ind+1,1)/ds(ind+1)^2 + (-1/60)*eN(ind+1,1)/ds(ind+1))))...
                   ((R_s(ind+1,1)*Pn(ind+1,1) - Material_E.*R_s(ind+1,1).*((2)*aN(ind+1,1)/ds(ind+1)^5 + (2)*bN(ind+1,1)/ds(ind+1)^4 - 1*cN(ind+1,1)/ds(ind+1)^3 + (-3/20)*dN(ind+1,1)/ds(ind+1)^2 + (3/20)*eN(ind+1,1)/ds(ind+1))))...
                   ((R_s(ind+1,1)*Pn(ind+1,1) - Material_E.*R_s(ind+1,1).*((-5/2)*aN(ind+1,1)/ds(ind+1)^5 + (-13/2)*bN(ind+1,1)/ds(ind+1)^4 + (13/8)*cN(ind+1,1)/ds(ind+1)^3 + (3/2)*dN(ind+1,1)/ds(ind+1)^2 + (-3/4)*eN(ind+1,1)/ds(ind+1))))...
                   ((R_s(ind+1,1)*Pn(ind+1,1) - Material_E.*R_s(ind+1,1).*((28/3)*bN(ind+1,1)/ds(ind+1)^4 + (-49/18)*dN(ind+1,1)/ds(ind+1)^2 + fN(ind+1,1))))...
                   ((R_s(ind+1,1)*Pn(ind+1,1) - Material_E.*R_s(ind+1,1).*((5/2)*aN(ind+1,1)/ds(ind+1)^5 + (-13/2)*bN(ind+1,1)/ds(ind+1)^4 + (-13/8)*cN(ind+1,1)/ds(ind+1)^3 + (3/2)*dN(ind+1,1)/ds(ind+1)^2 + (3/4)*eN(ind+1,1)/ds(ind+1))))...
                   ((R_s(ind+1,1)*Pn(ind+1,1) - Material_E.*R_s(ind+1,1).*((-2)*aN(ind+1,1)/ds(ind+1)^5 + (2)*bN(ind+1,1)/ds(ind+1)^4 + 1*cN(ind+1,1)/ds(ind+1)^3 + (-3/20)*dN(ind+1,1)/ds(ind+1)^2 + (-3/20)*eN(ind+1,1)/ds(ind+1))))...
                   ((R_s(ind+1,1)*Pn(ind+1,1) - Material_E.*R_s(ind+1,1).*((1/2)*aN(ind+1,1)/ds(ind+1)^5 + (-1/6)*bN(ind+1,1)/ds(ind+1)^4 + (-1/8)*cN(ind+1,1)/ds(ind+1)^3 + (1/90)*dN(ind+1,1)/ds(ind+1)^2 + (1/60)*eN(ind+1,1)/ds(ind+1))))];
        
         M_phi_in_1 = k_in(ind).*(phi_res_1);
         M_phi_in_2 = k_in(ind).*(phi_res_2);
         M_phi_fin_1 = k_fin(ind)*(phi_res_1);
         M_phi_fin_2 = k_fin(ind)*(phi_res_2); 
    
    
        if Restraint_intermediate == 1
            %incastro interno
            if ind == 1
            matA(NP_tot(end)+(6*ind+1),posizione(ind):posizione(ind)+13) = [zeros(1,3) 1 zeros(1,6) -1 zeros(1,3)];
            matA(NP_tot(end)+(6*ind+2),posizione(ind):posizione(ind)+13) = [v_res_1 -v_res_2];
            matA(NP_tot(end)+(6*ind+3),posizione(ind):posizione(ind)+13) = [phi_res_1 -phi_res_2];
            matA(NP_tot(end)+(6*ind+4),posizione(ind):posizione(ind)+13) = [M_res_1 -M_res_2];
            matA(NP_tot(end)+(6*ind+5),posizione(ind):posizione(ind)+13) = [T_res_1 -T_res_2];
            matA(NP_tot(end)+(6*ind+6),posizione(ind):posizione(ind)+13) = [N_res_1 -N_res_2];
            else
            matA(NP_tot(end)+(6*ind+1),posizione(ind)+(ind-1)*6+(ind-1):posizione(ind)+13+(ind-1)*6+(ind-1)) = [zeros(1,3) 1 zeros(1,6) -1 zeros(1,3)]; %sistemare
            matA(NP_tot(end)+(6*ind+2),posizione(ind)+(ind-1)*6+(ind-1):posizione(ind)+13+(ind-1)*6+(ind-1)) = [v_res_1 -v_res_2];
            matA(NP_tot(end)+(6*ind+3),posizione(ind)+(ind-1)*6+(ind-1):posizione(ind)+13+(ind-1)*6+(ind-1)) = [phi_res_1 -phi_res_2];
            matA(NP_tot(end)+(6*ind+4),posizione(ind)+(ind-1)*6+(ind-1):posizione(ind)+13+(ind-1)*6+(ind-1)) = [M_res_1 -M_res_2];
            matA(NP_tot(end)+(6*ind+5),posizione(ind)+(ind-1)*6+(ind-1):posizione(ind)+13+(ind-1)*6+(ind-1)) = [T_res_1 -T_res_2];
            matA(NP_tot(end)+(6*ind+6),posizione(ind)+(ind-1)*6+(ind-1):posizione(ind)+13+(ind-1)*6+(ind-1)) = [N_res_1 -N_res_2];
            end
    
    
        elseif Restraint_intermediate == 2
            %cerniera interna
            if ind == 1
            matA(NP_tot(end)+(6*ind+1),posizione(ind):posizione(ind)+13) = [zeros(1,3) 1 zeros(1,6) -1 zeros(1,3)];
            matA(NP_tot(end)+(6*ind+2),posizione(ind):posizione(ind)+13) = [v_res_1 -v_res_2];
            matA(NP_tot(end)+(6*ind+3),posizione(ind):posizione(ind)+13) = [M_res_1 zeros(1,7)];
            matA(NP_tot(end)+(6*ind+4),posizione(ind):posizione(ind)+13) = [zeros(1,7) M_res_2];
            matA(NP_tot(end)+(6*ind+5),posizione(ind):posizione(ind)+13) = [T_res_1 -T_res_2];
            matA(NP_tot(end)+(6*ind+6),posizione(ind):posizione(ind)+13) = [N_res_1 -N_res_2];
            else
            matA(NP_tot(end)+(6*ind+1),posizione(ind)+(ind-1)*6+(ind-1):posizione(ind)+13+(ind-1)*6+(ind-1)) = [zeros(1,3) 1 zeros(1,6) -1 zeros(1,3)]; %sistemare
            matA(NP_tot(end)+(6*ind+2),posizione(ind)+(ind-1)*6+(ind-1):posizione(ind)+13+(ind-1)*6+(ind-1)) = [v_res_1 -v_res_2];
            matA(NP_tot(end)+(6*ind+3),posizione(ind)+(ind-1)*6+(ind-1):posizione(ind)+13+(ind-1)*6+(ind-1)) = [M_res_1 zeros(1,7)];
            matA(NP_tot(end)+(6*ind+4),posizione(ind)+(ind-1)*6+(ind-1):posizione(ind)+13+(ind-1)*6+(ind-1)) = [zeros(1,7) M_res_2];
            matA(NP_tot(end)+(6*ind+5),posizione(ind)+(ind-1)*6+(ind-1):posizione(ind)+13+(ind-1)*6+(ind-1)) = [T_res_1 -T_res_2];
            matA(NP_tot(end)+(6*ind+6),posizione(ind)+(ind-1)*6+(ind-1):posizione(ind)+13+(ind-1)*6+(ind-1)) = [N_res_1 -N_res_2];
            end
    
        elseif Restraint_intermediate == 3
            %molla elastica
            if ind == 1
            matA(NP_tot(end)+(6*ind+1),posizione(ind):posizione(ind)+13) = [zeros(1,3) 1 zeros(1,6) -1 zeros(1,3)];
            matA(NP_tot(end)+(6*ind+2),posizione(ind):posizione(ind)+13) = [v_res_1 -v_res_2];
            matA(NP_tot(end)+(6*ind+3),posizione(ind):posizione(ind)+13) = [M_res_1 -M_res_2];
            matA(NP_tot(end)+(6*ind+4),posizione(ind):posizione(ind)+13) = [-M_phi_fin_1 (M_res_2+M_phi_fin_2)];
            matA(NP_tot(end)+(6*ind+5),posizione(ind):posizione(ind)+13) = [T_res_1 -T_res_2];
            matA(NP_tot(end)+(6*ind+6),posizione(ind):posizione(ind)+13) = [N_res_1 -N_res_2];
            else
            matA(NP_tot(end)+(6*ind+1),posizione(ind)+(ind-1)*6+(ind-1):posizione(ind)+13+(ind-1)*6+(ind-1)) = [zeros(1,3) 1 zeros(1,6) -1 zeros(1,3)]; %sistemare
            matA(NP_tot(end)+(6*ind+2),posizione(ind)+(ind-1)*6+(ind-1):posizione(ind)+13+(ind-1)*6+(ind-1)) = [v_res_1 -v_res_2];
            matA(NP_tot(end)+(6*ind+3),posizione(ind)+(ind-1)*6+(ind-1):posizione(ind)+13+(ind-1)*6+(ind-1)) = [M_res_1 -M_res_2];
            matA(NP_tot(end)+(6*ind+4),posizione(ind)+(ind-1)*6+(ind-1):posizione(ind)+13+(ind-1)*6+(ind-1)) = [-M_phi_fin_1 (M_res_2+M_phi_fin_2)];
            matA(NP_tot(end)+(6*ind+5),posizione(ind)+(ind-1)*6+(ind-1):posizione(ind)+13+(ind-1)*6+(ind-1)) = [T_res_1 -T_res_2];
            matA(NP_tot(end)+(6*ind+6),posizione(ind)+(ind-1)*6+(ind-1):posizione(ind)+13+(ind-1)*6+(ind-1)) = [N_res_1 -N_res_2];
            end
        end
        B(NP_tot(end)+(6*ind+1)) = 0;
        B(NP_tot(end)+(6*ind+2)) = 0;
        B(NP_tot(end)+(6*ind+3)) = 0;
        B(NP_tot(end)+(6*ind+4)) = 0;
        B(NP_tot(end)+(6*ind+5)) = 0;
        B(NP_tot(end)+(6*ind+6)) = 0;
    end
    
    
    %% Definizione delle equazioni relative ai vincoli 

    %% Vincolo 1 (al piede sx)
    
    %spostamento u
    cc1_1 = 1;  %i
    %spostamento v
    cc2_1 = [(1/60)*R_s(1,1)/ds(1) (-3/20)*R_s(1,1)/ds(1) (3/4)*R_s(1,1)/ds(1) 0 (-3/4)*R_s(1,1)/ds(1) (3/20)*R_s(1,1)/ds(1) (-1/60)*R_s(1,1)/ds(1)];
    
    if Restraint == 1          %Incastro
        if molla_in ~= 1
        %rotazione phi
            cc3_1 = [((-1/90)*R_s(1,1)/ds(1)^2 + (1/60)*dR_s(1,1)/ds(1))... %u(i-3)
                    ((3/20)*R_s(1,1)/ds(1)^2 - (3/20)*dR_s(1,1)/ds(1))...   %u(i-2)
                    ((-3/2)*R_s(1,1)/ds(1)^2 + (3/4)*dR_s(1,1)/ds(1))...    %u(i-1)  
                    ((49/18)*R_s(1,1)/ds(1)^2 - 1/R_s(1,1))...              %u(i)  
                    ((-3/2)*R_s(1,1)/ds(1)^2 - (3/4)*dR_s(1,1)/ds(1))...    %u(i+1)  
                    ((3/20)*R_s(1,1)/ds(1)^2 + (3/20)*dR_s(1,1)/ds(1))...   %u(i+2)  
                    ((-1/90)*R_s(1,1)/ds(1)^2 - (1/60)*dR_s(1,1)/ds(1))];   %u(i+3)
    
        else                   %Vincolo molla elastica
            M_in = +Material_E.*[J(1,1).*((1/8)*aM(1,1)/ds(1)^3 + (1/90)*bM(1,1)/ds(1)^2 + (-1/60)*cM(1,1)/ds(1))...
                                J(1,1).*(-1*aM(1,1)/ds(1)^3 + (-3/20)*bM(1,1)/ds(1)^2 + (3/20)*cM(1,1)/ds(1))...
                                J(1,1).*((13/8)*aM(1,1)/ds(1)^3 + (3/2)*bM(1,1)/ds(1)^2 + (-3/4)*cM(1,1)/ds(1))...
                                J(1,1).*((-49/18)*bM(1,1)/ds(1)^2 + dM(1,1))...
                                J(1,1).*((-13/8)*aM(1,1)/ds(1)^3 + (3/2)*bM(1,1)/ds(1)^2 + (3/4)*cM(1,1)/ds(1))...
                                J(1,1).*(1*aM(1,1)/ds(1)^3 + (-3/20)*bM(1,1)/ds(1)^2 + (-3/20)*cM(1,1)/ds(1))...
                                J(1,1).*((-1/8)*aM(1,1)/ds(1)^3 + (1/90)*bM(1,1)/ds(1)^2 + (1/60)*cM(1,1)/ds(1))]; %u(i+3)
    
            phi_in = [((-1/90)*R_s(1,1)/ds(1)^2 + (1/60)*dR_s(1,1)/ds(1))...
                    ((3/20)*R_s(1,1)/ds(1)^2 - (3/20)*dR_s(1,1)/ds(1))...
                    ((-3/2)*R_s(1,1)/ds(1)^2 + (3/4)*dR_s(1,1)/ds(1))...
                    ((49/18)*R_s(1,1)/ds(1)^2 - 1/R_s(1,1))...
                    ((-3/2)*R_s(1,1)/ds(1)^2 - (3/4)*dR_s(1,1)/ds(1))...
                    ((3/20)*R_s(1,1)/ds(1)^2 + (3/20)*dR_s(1,1)/ds(1))...
                    ((-1/90)*R_s(1,1)/ds(1)^2 - (1/60)*dR_s(1,1)/ds(1))];
            
            cc3_1 = M_in - (-phi_in.*k_in(1));
        end
    
    elseif Restraint == 2        % Cerniera
        %Momento
        cc3_1 = Material_E.*[J(1,1).*((1/8)*aM(1,1)/ds(1)^3 + (1/90)*bM(1,1)/ds(1)^2 + (-1/60)*cM(1,1)/ds(1))...   %u(i-3)
                            J(1,1).*(-1*aM(1,1)/ds(1)^3 + (-3/20)*bM(1,1)/ds(1)^2 + (3/20)*cM(1,1)/ds(1))...       %u(i-2)
                            J(1,1).*((13/8)*aM(1,1)/ds(1)^3 + (3/2)*bM(1,1)/ds(1)^2 + (-3/4)*cM(1,1)/ds(1))...     %u(i-1)
                            J(1,1).*((-49/18)*bM(1,1)/ds(1)^2 + dM(1,1))...                                         %u(i)
                            J(1,1).*((-13/8)*aM(1,1)/ds(1)^3 + (3/2)*bM(1,1)/ds(1)^2 + (3/4)*cM(1,1)/ds(1))...     %u(i+1)
                            J(1,1).*(1*aM(1,1)/ds(1)^3 + (-3/20)*bM(1,1)/ds(1)^2 + (-3/20)*cM(1,1)/ds(1))...       %u(i+2)
                            J(1,1).*((-1/8)*aM(1,1)/ds(1)^3 + (1/90)*bM(1,1)/ds(1)^2 + (1/60)*cM(1,1)/ds(1))];     %u(i+3)
    
    end
    
    % n_cc sarà quello finale che rimane alla fine del ciclo for precedente 
    
    matA(NP_tot(end)+1,1:7) = [zeros(1,3) cc1_1 zeros(1,3)];
    matA(NP_tot(end)+2,1:7) = cc2_1;
    matA(NP_tot(end)+3,1:7) = cc3_1;
    
    B(NP_tot(end)+1) = 0;
    B(NP_tot(end)+2) = 0;
    B(NP_tot(end)+3) = 0;
    
    %% Vincolo 2 (al piede dx)
    
    % Mi serve un indice che mi dica quanti archi ho e quindi pari alla
    % lunghezza del vettore x_point +1 ---> n_a
    
     %spostamento u
        cc1_2 = 1;  %i
        %spostamento v
        cc2_2 = [(1/60)*R_s(n_a,NP(end))/ds(n_a) (-3/20)*R_s(n_a,NP(end))/ds(n_a) (3/4)*R_s(n_a,NP(end))/ds(n_a) 0 (-3/4)*R_s(n_a,NP(end))/ds(n_a) (3/20)*R_s(n_a,NP(end))/ds(n_a) (-1/60)*R_s(n_a,NP(end))/ds(n_a)];

    if Restraint == 1
        %Vincolo incastro
        if molla_fin ~= NP_1
            %rotazione phi
            cc3_2 = [((-1/90)*R_s(n_a,NP(end))/ds(n_a)^2 + (1/60)*dR_s(n_a,NP(end))/ds(n_a))...
                    ((3/20)*R_s(n_a,NP(end))/ds(n_a)^2 - (3/20)*dR_s(n_a,NP(end))/ds(n_a))...
                    ((-3/2)*R_s(n_a,NP(end))/ds(n_a)^2 + (3/4)*dR_s(n_a,NP(end))/ds(n_a))...
                    ((49/18)*R_s(n_a,NP(end))/ds(n_a)^2 - 1/R_s(n_a,NP(end)))...
                    ((-3/2)*R_s(n_a,NP(end))/ds(n_a)^2 - (3/4)*dR_s(n_a,NP(end))/ds(n_a))...
                    ((3/20)*R_s(n_a,NP(end))/ds(n_a)^2 + (3/20)*dR_s(n_a,NP(end))/ds(n_a))...
                    ((-1/90)*R_s(n_a,NP(end))/ds(n_a)^2 - (1/60)*dR_s(n_a,NP(end))/ds(n_a))];
    
        else % Molla elastica
            M_fin = +Material_E.*[J(end,NP(end)).*((1/8)*aM(end,NP(end))/ds(end)^3 + (1/90)*bM(end,NP(end))/ds(end)^2 + (-1/60)*cM(end,NP(end))/ds(end))...
                                J(end,NP(end)).*(-1*aM(end,NP(end))/ds(end)^3 + (-3/20)*bM(end,NP(end))/ds(end)^2 + (3/20)*cM(end,NP(end))/ds(end))...
                                J(end,NP(end)).*((13/8)*aM(end,NP(end))/ds(end)^3 + (3/2)*bM(end,NP(end))/ds(end)^2 + (-3/4)*cM(end,NP(end))/ds(end))...
                                J(end,NP(end)).*((-49/18)*bM(end,NP(end))/ds(end)^2 + dM(end,NP(end)))...
                                J(end,NP(end)).*((-13/8)*aM(end,NP(end))/ds(end)^3 + (3/2)*bM(end,NP(end))/ds(end)^2 + (3/4)*cM(end,NP(end))/ds(end))...
                                J(end,NP(end)).*(1*aM(end,NP(end))/ds(end)^3 + (-3/20)*bM(end,NP(end))/ds(end)^2 + (-3/20)*cM(end,NP(end))/ds(end))...
                                J(end,NP(end)).*((-1/8)*aM(end,NP(end))/ds(end)^3 + (1/90)*bM(end,NP(end))/ds(end)^2 + (1/60)*cM(end,NP(end))/ds(end))]; %u(i+3)
    
            phi_fin = [((-1/90)*R_s(end,NP(end))/ds(end)^2 + (1/60)*dR_s(end,NP(end))/ds(end))...
                    ((3/20)*R_s(end,NP(end))/ds(end)^2 - (3/20)*dR_s(end,NP(end))/ds(end))...
                    ((-3/2)*R_s(end,NP(end))/ds(end)^2 + (3/4)*dR_s(end,NP(end))/ds(end))...
                    ((49/18)*R_s(end,NP(end))/ds(end)^2 - 1/R_s(end,NP(end)))...
                    ((-3/2)*R_s(end,NP(end))/ds(end)^2 - (3/4)*dR_s(end,NP(end))/ds(end))...
                    ((3/20)*R_s(end,NP(end))/ds(end)^2 + (3/20)*dR_s(end,NP(end))/ds(end))...
                    ((-1/90)*R_s(end,NP(end))/ds(end)^2 - (1/60)*dR_s(end,NP(end))/ds(end))];
            
            cc3_2 = M_fin - (phi_fin.*k_fin(end));
        end
    
    
    elseif Restraint == 2
        % Vincolo cerniera
            cc3_2 = +Material_E.*[J(end,NP(end)).*((1/8)*aM(end,NP(end))/ds(end)^3 + (1/90)*bM(end,NP(end))/ds(end)^2 + (-1/60)*cM(end,NP(end))/ds(end))...
                                J(end,NP(end)).*(-1*aM(end,NP(end))/ds(end)^3 + (-3/20)*bM(end,NP(end))/ds(end)^2 + (3/20)*cM(end,NP(end))/ds(end))...
                                J(end,NP(end)).*((13/8)*aM(end,NP(end))/ds(end)^3 + (3/2)*bM(end,NP(end))/ds(end)^2 + (-3/4)*cM(end,NP(end))/ds(end))...
                                J(end,NP(end)).*((-49/18)*bM(end,NP(end))/ds(end)^2 + dM(end,NP(end)))...
                                J(end,NP(end)).*((-13/8)*aM(end,NP(end))/ds(end)^3 + (3/2)*bM(end,NP(end))/ds(end)^2 + (3/4)*cM(end,NP(end))/ds(end))...
                                J(end,NP(end)).*(1*aM(end,NP(end))/ds(end)^3 + (-3/20)*bM(end,NP(end))/ds(end)^2 + (-3/20)*cM(end,NP(end))/ds(end))...
                                J(end,NP(end)).*((-1/8)*aM(end,NP(end))/ds(end)^3 + (1/90)*bM(end,NP(end))/ds(end)^2 + (1/60)*cM(end,NP(end))/ds(end))]; %u(i+3)
    end
    
    
    matA(NP_tot(end)+4,NP_tot(end)+6*(n_a-1):NP_tot(end)+6*(n_a-1)+6) = [zeros(1,3) cc1_2 zeros(1,3)];
    matA(NP_tot(end)+5,NP_tot(end)+6*(n_a-1):NP_tot(end)+6*(n_a-1)+6) = cc2_2;
    matA(NP_tot(end)+6,NP_tot(end)+6*(n_a-1):NP_tot(end)+6*(n_a-1)+6) = cc3_2;
    
    B(NP_tot(end)+4) = 0;
    B(NP_tot(end)+5) = 0;
    B(NP_tot(end)+6) = 0;
    
    x_res = matA\B';
    u_sol = x_res';
    
    % Allocazione di memoria
    u_exp_new = zeros(n_a,NP_max+6);
    
    % Una volta risolto il sistema, avrò le soluzioni in un unico vettore e
    % quindi le devo riportare in più vettori, uno per ogni curva definita
    
    for ind = 1:length(x_point)+1
        n_cc = ind*6;
        u_exp_new(ind,1:NP(ind+1)+6) = u_sol(NP_tot(ind)+n_cc-6+1:NP_tot(ind+1)+n_cc);
    end
    
    %% Allocazione di memoria
    
    u = zeros(n_a,NP_max);
    du = zeros(n_a,NP_max);
    d2u = zeros(n_a,NP_max);
    d3u = zeros(n_a,NP_max);
    d4u = zeros(n_a,NP_max);
    d5u = zeros(n_a,NP_max);
    v = zeros(n_a,NP_max);
    phi = zeros(n_a,NP_max);
    k = zeros(n_a,NP_max);
    dk = zeros(n_a,NP_max);
    d2k = zeros(n_a,NP_max);
    M = zeros(n_a,NP_max);
    N = zeros(n_a,NP_max);
    T = zeros(n_a,NP_max);
    Dz = zeros(n_a,NP_max);
    Dx = zeros(n_a,NP_max);
    Results = zeros(NP_1,6);
    Results_sing = zeros(NP_tot(end),6);
    s_plot = zeros(1,NP_1);
    x_plot = zeros(1,NP_1);
    M_plot = zeros(1,NP_1);
    N_plot = zeros(1,NP_1);
    T_plot = zeros(1,NP_1);
    Dx_plot = zeros(1,NP_1);
    Dz_plot = zeros(1,NP_1);
    phi_plot = zeros(1,NP_1);
    sigma = zeros(n_a,NP_max);
    sigma_max = zeros(n_a,1);
    constraint = zeros(n_a,1);
    z = zeros(7*n_a,NP_max);
    sigma_NM = zeros(7*n_a,NP_max);
    tau_max = zeros(n_a,NP_max);
    
    if Section.type == 6 || Section.type == 7
        tau = zeros(n_a,NP_max);
    else
        tau = zeros(7*n_a,NP_max);
    end
    
    sigma_max_VM = zeros(1,n_a);
    sigma_VMmax = zeros(n_a,NP_max);
    sigma_VM = zeros(7*n_a,NP_max);
    s_fin = zeros (1,n_a);
    x_fin = zeros (1,n_a);
    
    %% Calcolo degli spostamenti e le loro derivate
    
    for ind = 1:length(x_point)+1
    
        u(ind,1:NP(ind+1)) = u_exp_new(ind,4:NP(ind+1)+3);
        du(ind,1:NP(ind+1)) = ((-1/60).*u_exp_new(ind,1:NP(ind+1)) + (3/20).*u_exp_new(ind,2:NP(ind+1)+1) + (-3/4).*u_exp_new(ind,3:NP(ind+1)+2) + (0).*u_exp_new(ind,4:NP(ind+1)+3) + (3/4).*u_exp_new(ind,5:NP(ind+1)+4) + (-3/20).*u_exp_new(ind,6:NP(ind+1)+5) + (1/60).*u_exp_new(ind,7:NP(ind+1)+6))./ds(ind);
        d2u(ind,1:NP(ind+1)) = ((1/90).*u_exp_new(ind,1:NP(ind+1)) + (-3/20).*u_exp_new(ind,2:NP(ind+1)+1) + (3/2).*u_exp_new(ind,3:NP(ind+1)+2) + (-49/18).*u_exp_new(ind,4:NP(ind+1)+3) + (3/2).*u_exp_new(ind,5:NP(ind+1)+4) + (-3/20).*u_exp_new(ind,6:NP(ind+1)+5) + (1/90).*u_exp_new(ind,7:NP(ind+1)+6))./(ds(ind)^2);
        d3u(ind,1:NP(ind+1)) = ((1/8).*u_exp_new(ind,1:NP(ind+1)) + (-1).*u_exp_new(ind,2:NP(ind+1)+1) + (13/8).*u_exp_new(ind,3:NP(ind+1)+2) + (0).*u_exp_new(ind,4:NP(ind+1)+3) + (-13/8).*u_exp_new(ind,5:NP(ind+1)+4) + (1).*u_exp_new(ind,6:NP(ind+1)+5) + (-1/8).*u_exp_new(ind,7:NP(ind+1)+6))./(ds(ind)^3);
        d4u(ind,1:NP(ind+1)) = ((-1/6).*u_exp_new(ind,1:NP(ind+1)) + (2).*u_exp_new(ind,2:NP(ind+1)+1) + (-13/2).*u_exp_new(ind,3:NP(ind+1)+2) + (28/3).*u_exp_new(ind,4:NP(ind+1)+3) + (-13/2).*u_exp_new(ind,5:NP(ind+1)+4) + (2).*u_exp_new(ind,6:NP(ind+1)+5) + (-1/6).*u_exp_new(ind,7:NP(ind+1)+6))./(ds(ind)^4);
        d5u(ind,1:NP(ind+1)) = ((-1/2).*u_exp_new(ind,1:NP(ind+1)) + (2).*u_exp_new(ind,2:NP(ind+1)+1) + (-5/2).*u_exp_new(ind,3:NP(ind+1)+2) + (0).*u_exp_new(ind,4:NP(ind+1)+3) + (5/2).*u_exp_new(ind,5:NP(ind+1)+4) + (-2).*u_exp_new(ind,6:NP(ind+1)+5) + (1/2).*u_exp_new(ind,7:NP(ind+1)+6))./(ds(ind)^5);
        
        % Perpendicular displacements [m]
        v(ind,1:NP(ind+1)) = -du(ind,1:NP(ind+1)).*R_s(ind,1:NP(ind+1));
        
        % Rotations [rad]
        phi(ind,1:NP(ind+1)) = -(d2u(ind,1:NP(ind+1)).*R_s(ind,1:NP(ind+1)) + du(ind,1:NP(ind+1)).*dR_s(ind,1:NP(ind+1)) + u(ind,1:NP(ind+1))./R_s(ind,1:NP(ind+1)));
        
        % Deformation angle and its derivatives
        k(ind,1:NP(ind+1)) = d3u(ind,1:NP(ind+1)).*d1(ind,1:NP(ind+1)) + d2u(ind,1:NP(ind+1)).*e1(ind,1:NP(ind+1)) + du(ind,1:NP(ind+1)).*f1(ind,1:NP(ind+1)) + u(ind,1:NP(ind+1)).*g1(ind,1:NP(ind+1));
        dk(ind,1:NP(ind+1)) = c2(ind,1:NP(ind+1)).*d4u(ind,1:NP(ind+1)) + d2(ind,1:NP(ind+1)).*d3u(ind,1:NP(ind+1)) + e2(ind,1:NP(ind+1)).*d2u(ind,1:NP(ind+1)) + f2(ind,1:NP(ind+1)).*du(ind,1:NP(ind+1)) + g2(ind,1:NP(ind+1)).*u(ind,1:NP(ind+1));
        d2k(ind,1:NP(ind+1)) = b3(ind,1:NP(ind+1)).*d5u(ind,1:NP(ind+1)) + c3(ind,1:NP(ind+1)).*d4u(ind,1:NP(ind+1)) + d3(ind,1:NP(ind+1)).*d3u(ind,1:NP(ind+1)) + e3(ind,1:NP(ind+1)).*d2u(ind,1:NP(ind+1)) + f3(ind,1:NP(ind+1)).*du(ind,1:NP(ind+1)) + g3(ind,1:NP(ind+1)).*u(ind,1:NP(ind+1));
    end
    %% Stabilization of d2k in case of discontinuous overload

    % Per poter interpolare fra i punti prima e dopo dvo creare dei vettori
    % interi dei punti dell'ascissa curvilinea, dei carichi q_ol e delle
    % derivate della curvatura
    
    for ind = 1:length(x_point)+1
        kk = kk + NP(ind);
        if ind == 1
            s_fin(1) = 0;
            s_int(kk+1:kk+NP(ind+1)) = s_x(ind,1:NP(ind+1))+s_fin(1);
            q_ol_int(kk+1:kk+NP(ind+1)) = q_ol(ind,1:NP(ind+1));
            d2k_int(kk+1:kk+NP(ind+1)) = d2k(ind,1:NP(ind+1));
        elseif ind  ~= 1
            s_fin(ind) = s_fin(ind-1)+s_x(ind-1,NP(ind));
            s_int(kk+1:kk+NP(ind+1)) = s_x(ind,1:NP(ind+1))+s_fin(ind);
            q_ol_int(kk+1:kk+NP(ind+1)) = q_ol(ind,1:NP(ind+1));
            d2k_int(kk+1:kk+NP(ind+1)) = d2k(ind,1:NP(ind+1));
        end
    end
    kk = 0;

    % for ind = 1:length(x_point)
    %     for i = 1:NP_tot(end-1)
    %         if i == NP_tot(ind+1)
    %            d2k_int(i-1) = interp1(s(ind,1:NP(ind+1)-2), d2k(ind,1:NP(ind+1)-2), s(ind,NP(ind+1)), 'pchip','extrap');
    %            d2k_int(i) = interp1(s(ind+1,2:NP(ind+2)-1), d2k(ind+1,2:NP(ind+2)-1), s(ind+1,1), 'pchip','extrap');
    %         end
    %      end
    % end

    for i=1:NP_tot(end)-3
        if q_ol_int(i) ~= 0
            if q_ol_int(i+1) == 0
                d2k_int(i) = interp1(s_int(1:i-1), d2k_int(1:i-1), s_int(i), 'pchip','extrap');
                d2k_int(i+1) = interp1(s_int(i+2:NP_tot(end)), d2k_int(i+2:NP_tot(end)), s_int(i+1), 'pchip','extrap');
            end
        else
            if q_ol_int(i+1) ~= 0
                d2k_int(i) = interp1(s_int(1:i-1), d2k_int(1:i-1), s_int(i), 'pchip','extrap');
                d2k_int(i+1) = interp1(s_int(i+2:NP_tot(end)), d2k_int(i+2:NP_tot(end)), s_int(i+1), 'pchip','extrap');
            end
        end            
    end
    

   % Devo riportare le curvature in forma matriciale
    for ind = 1:length(x_point)+1
        kk = kk + NP(ind);
        if ind == 1
           d2k(ind,1:NP(ind+1)) = d2k_int(kk+1:kk+NP(ind+1));
        elseif ind ~= 1
           d2k(ind,1:NP(ind+1)) = d2k_int(kk+1:kk+NP(ind+1));
        end
    end
    kk = 0;
    
    %% Calcolo delle sollecitazioni
    for ind = 1:length(x_point)+1
        
        %% Solicitations 1
    
        M(ind,1:NP(ind+1)) = Material.E.*J(ind,1:NP(ind+1)).*k(ind,1:NP(ind+1));                                                                                           % Bending Moment [kNm]
        T(ind,1:NP(ind+1)) = Material.E.*-(dJ(ind,1:NP(ind+1)).*k(ind,1:NP(ind+1)) + J(ind,1:NP(ind+1)).*dk(ind,1:NP(ind+1))) + m_load(ind,1:NP(ind+1));                                                 % Shear [kN]
        N(ind,1:NP(ind+1)) = R_s(ind,1:NP(ind+1)).*Pn(ind,1:NP(ind+1)) - R_s(ind,1:NP(ind+1)).*Material.E.*(d2J(ind,1:NP(ind+1)).*k(ind,1:NP(ind+1)) + 2.*dJ(ind,1:NP(ind+1)).*dk(ind,1:NP(ind+1)) + J(ind,1:NP(ind+1)).*d2k(ind,1:NP(ind+1)));  % Axial force [kN]
        
        %% Displacements
    
        Dz(ind,1:NP(ind+1)) = -(v(ind,1:NP(ind+1)).*cos(alpha(ind,1:NP(ind+1)))+u(ind,1:NP(ind+1)).*sin(alpha(ind,1:NP(ind+1))));
        Dx(ind,1:NP(ind+1)) = -(-v(ind,1:NP(ind+1)).*sin(alpha(ind,1:NP(ind+1)))+u(ind,1:NP(ind+1)).*cos(alpha(ind,1:NP(ind+1))));
    
        %% Plot

        kk = kk + NP(ind);
    
        if ind == 1 && (Restraint_intermediate == 1 || Restraint_intermediate == 2)
            s_fin(1) = 0;
            x_fin(1) = 0;
            s_plot(kk+1:kk+NP(ind+1)) = s_x(ind,1:NP(ind+1))+s_fin(1);
            x_plot(kk+1:kk+NP(ind+1)) = x_s(ind,1:NP(ind+1))+x_fin(1);
            M_plot(kk+1:kk+NP(ind+1)) = M(ind,1:NP(ind+1));
            T_plot(kk+1:kk+NP(ind+1)) = T(ind,1:NP(ind+1));
            N_plot(kk+1:kk+NP(ind+1)) = N(ind,1:NP(ind+1));
            Dx_plot(kk+1:kk+NP(ind+1)) = Dx(ind,1:NP(ind+1));
            Dz_plot(kk+1:kk+NP(ind+1)) = Dz(ind,1:NP(ind+1));
            phi_plot(kk+1:kk+NP(ind+1)) = phi(ind,1:NP(ind+1));
        elseif ind  ~= 1 && (Restraint_intermediate == 1 || Restraint_intermediate ==2)
            s_fin(ind) = s_fin(ind-1)+s_x(ind-1,NP(ind));
            x_fin(ind) = x_fin(ind-1)+x_s(ind-1,NP(ind));
            s_plot(kk-(ind-3):kk+NP(ind+1)-(ind)+1) = s_x(ind,2:NP(ind+1))+s_fin(ind);
            x_plot(kk-(ind-3):kk+NP(ind+1)-(ind)+1) = x_s(ind,2:NP(ind+1))+x_fin(ind);
            M_plot(kk-(ind-3):kk+NP(ind+1)-(ind)+1) = M(ind,2:NP(ind+1));
            T_plot(kk-(ind-3):kk+NP(ind+1)-(ind)+1) = T(ind,2:NP(ind+1));
            N_plot(kk-(ind-3):kk+NP(ind+1)-(ind)+1) = N(ind,2:NP(ind+1));
            Dx_plot(kk-(ind-3):kk+NP(ind+1)-(ind)+1) = Dx(ind,2:NP(ind+1));
            Dz_plot(kk-(ind-3):kk+NP(ind+1)-(ind)+1) = Dz(ind,2:NP(ind+1));
            phi_plot(kk-(ind-3):kk+NP(ind+1)-(ind)+1) = phi(ind,2:NP(ind+1));
        end
    
        if ind == 1 && Restraint_intermediate == 3
            s_fin(1) = 0;
            x_fin(1) = 0;
            s_plot(kk+1:kk+NP(ind+1)) = s_x(ind,1:NP(ind+1))+s_fin(1);
            x_plot(kk+1:kk+NP(ind+1)) = x_s(ind,1:NP(ind+1))+x_fin(1);
            M_plot(kk+1:kk+NP(ind+1)) = M(ind,1:NP(ind+1));
            T_plot(kk+1:kk+NP(ind+1)) = T(ind,1:NP(ind+1));
            N_plot(kk+1:kk+NP(ind+1)) = N(ind,1:NP(ind+1));
            Dx_plot(kk+1:kk+NP(ind+1)) = Dx(ind,1:NP(ind+1));
            Dz_plot(kk+1:kk+NP(ind+1)) = Dz(ind,1:NP(ind+1));
            phi_plot(kk+1:kk+NP(ind+1)) = phi(ind,1:NP(ind+1));
        elseif ind  ~= 1 && Restraint_intermediate == 3
            s_fin(ind) = s_fin(ind-1)+s_x(ind-1,NP(ind));
            x_fin(ind) = x_fin(ind-1)+x_s(ind-1,NP(ind));
            s_plot(kk+1:kk+NP(ind+1)) = s_x(ind,1:NP(ind+1))+s_fin(ind);
            x_plot(kk+1:kk+NP(ind+1)) = x_s(ind,1:NP(ind+1))+x_fin(ind);
            M_plot(kk+1:kk+NP(ind+1)) = M(ind,1:NP(ind+1));
            T_plot(kk+1:kk+NP(ind+1)) = T(ind,1:NP(ind+1));
            N_plot(kk+1:kk+NP(ind+1)) = N(ind,1:NP(ind+1));
            Dx_plot(kk+1:kk+NP(ind+1)) = Dx(ind,1:NP(ind+1));
            Dz_plot(kk+1:kk+NP(ind+1)) = Dz(ind,1:NP(ind+1));
            phi_plot(kk+1:kk+NP(ind+1)) = phi(ind,1:NP(ind+1));
       end
   end
   
   %% Tabella risultati
   for ind = 1:length(x_point)+1

        Results_sing(NP_tot(ind)+1:NP_tot(ind)+NP(ind+1),1:6) = [Dx(ind,1:(NP(ind+1)))' Dz(ind,1:(NP(ind+1)))' phi(ind,1:(NP(ind+1)))' M(ind,1:(NP(ind+1)))' N(ind,1:(NP(ind+1)))' T(ind,1:(NP(ind+1)))'];
        %Results_sing(NP_tot(ind)+1:NP_tot(ind)+NP(ind+1)+1,1:6) = NaN(1,6);
    
        hh = hh + NP(ind+1);
    
        if ind == 1
        Results(1:NP(ind+1)-1,:) = [Dx(1,1:(NP(2)-1))' Dz(1,1:(NP(2)-1))' phi(1,1:(NP(2)-1))' M(1,1:(NP(2)-1))' N(1,1:(NP(2)-1))' T(1,1:(NP(2)-1))'];
        elseif ind~=1
        Results((hh-NP(ind+1)-(ind-2)):(hh-ind),:) = [Dx(ind,1:(NP(ind+1)-1))' Dz(ind,1:(NP(ind+1)-1))' phi(ind,1:(NP(ind+1)-1))' M(ind,1:(NP(ind+1)-1))' N(ind,1:(NP(ind+1)-1))' T(ind,1:(NP(ind+1)-1))'];
        end
        
        if ind == length(x_point)+1
        Results(hh-ind+1,1:6) = [Dx(ind,NP(ind+1)) Dz(ind,NP(ind+1)) phi(ind,NP(ind+1)) M(ind,NP(ind+1)) N(ind,NP(ind+1)) T(ind,NP(ind+1))];
        end

        %% Stresses
        sigma(ind,1:NP(ind+1)) = (M(ind,1:NP(ind+1))./W(ind,1:NP(ind+1)) + N(ind,1:NP(ind+1))./A(ind,1:NP(ind+1)))./1000;      % Normal stress [MPa]
        
        sigma_max(ind,1:NP(ind+1)) = max(abs(sigma(ind,1:NP(ind+1))));
        
        constraint(ind,1:NP(ind+1)) = sigma_max(ind,1:NP(ind+1))/Material.fyd - 1;                                             % Strength verification
        
        %% Von Mises stress calculation
    
        if ind == 1
            z(1:7,1:NP(ind+1)) = [-r1(ind,1:NP(ind+1)); -(2/3).*r1(ind,1:NP(ind+1)); -(1/3).*r1(ind,1:NP(ind+1)); zeros(1,length(r1(ind,1:NP(ind+1)))); (1/3).*r1(ind,1:NP(ind+1)); (2/3).*r1(ind,1:NP(ind+1)); r1(ind,1:NP(ind+1))];
            if Section.type == 1                % Circular
                sigma_NM = (z(1:7,1:NP(ind+1)).*M(ind,1:NP(ind+1))./J(ind,1:NP(ind+1)) + N(ind,1:NP(ind+1))./A(ind,1:NP(ind+1)))./1000;
            elseif Section.type == 2            % Ellipsoidal Pay attention, you are using the circular approximation!!!!!!!!!!!!!!!!!!!)
                sigma_NM = (z(1:7,1:NP(ind+1)).*M(ind,1:NP(ind+1))./J(ind,1:NP(ind+1)) + N(ind,1:NP(ind+1))./A(ind,1:NP(ind+1)))./1000;
            elseif Section.type == 3            % Square
                sigma_NM = ((z(1:7,1:NP(ind+1))./2).*M(ind,1:NP(ind+1))./J(ind,1:NP(ind+1)) + N(ind,1:NP(ind+1))./A(ind,1:NP(ind+1)))./1000;
            elseif Section.type == 4            % Rectangular
                sigma_NM = ((z(1:7,1:NP(ind+1))./2).*M(ind,1:NP(ind+1))./J(ind,1:NP(ind+1)) + N(ind,1:NP(ind+1))./A(ind,1:NP(ind+1)))./1000;
            elseif Section.type == 5            % Hollowed circular
                sigma_NM = (z(1:7,1:NP(ind+1)).*M(ind,1:NP(ind+1))./J(ind,1:NP(ind+1)) + N(ind,1:NP(ind+1))./A(ind,1:NP(ind+1)))./1000;
            elseif Section.type == 6            % Hollowed square
                sigma_NM = ((z(1:7,1:NP(ind+1))./2).*M(ind,1:NP(ind+1))./J(ind,1:NP(ind+1)) + N(ind,1:NP(ind+1))./A(ind,1:NP(ind+1)))./1000;
            elseif Section.type == 7            % Hollowed rectangular
                sigma_NM = ((z(1:7,1:NP(ind+1))./2).*M(ind,1:NP(ind+1))./J(ind,1:NP(ind+1)) + N(ind,1:NP(ind+1))./A(ind,1:NP(ind+1)))./1000;
            end
            
    
            if Section.type == 1                % Circular
                tau_max(ind,1:NP(ind+1)) = ((4/3).*(T(ind,1:NP(ind+1))./A(ind,1:NP(ind+1))))./1000;
                tau(1:7,1:NP(ind+1)) = -(tau_max(ind,1:NP(ind+1))./(r1(ind,1:NP(ind+1)).^2)).*(z(1:7,1:NP(ind+1)).^2) + tau_max(ind,1:NP(ind+1));
            elseif Section.type == 2            % Ellipsoidal Pay attention, you are using the circular approximation!!!!!!!!!!!!!!!!!!!)
                tau_max(ind,1:NP(ind+1)) = ((4/3).*(T(ind,1:NP(ind+1))./A(ind,1:NP(ind+1))))./1000;
                tau(1:7,1:NP(ind+1)) = -(tau_max(ind,1:NP(ind+1))./(r1(ind,1:NP(ind+1)).^2)).*(z(1:7,1:NP(ind+1)).^2) + tau_max(ind,1:NP(ind+1));
            elseif Section.type == 3            % Square
                tau_max(ind,1:NP(ind+1)) = ((3/2).*(T(ind,1:NP(ind+1))./A(ind,1:NP(ind+1))))./1000;
                tau(1:7,1:NP(ind+1)) = -(tau_max(ind,1:NP(ind+1))./((r1(ind,1:NP(ind+1))/2).^2)).*((z(1:7,1:NP(ind+1))./2).^2) + tau_max(ind,1:NP(ind+1));
            elseif Section.type == 4            % Rectangular
                tau_max(ind,1:NP(ind+1)) = ((3/2).*(T(ind,1:NP(ind+1))./A(ind,1:NP(ind+1))))./1000;
                tau(1:7,1:NP(ind+1)) = -(tau_max(ind,1:NP(ind+1))./((r1(ind,1:NP(ind+1))/2).^2)).*((z(1:7,1:NP(ind+1))./2).^2) + tau_max(ind,1:NP(ind+1));
            elseif Section.type == 5            % Hollowed circular
                tau_max(ind,1:NP(ind+1)) = (2.*(T(ind,1:NP(ind+1))./A(ind,1:NP(ind+1))))./1000;
                tau(1:7,1:NP(ind+1)) = -(tau_max(ind,1:NP(ind+1))./(r1(ind,1:NP(ind+1)).^2)).*(z(1:7,1:NP(ind+1)).^2) + tau_max(ind,1:NP(ind+1));
            elseif Section.type == 6            % Hollowed square
                tau_max(ind,1:NP(ind+1)) = (T(ind,1:NP(ind+1))./(A(ind,1:NP(ind+1)) - (2.*r1(ind,1:NP(ind+1)).*t(1))))./1000;
                tau(1,1:NP(ind+1)) = tau_max(1,1:NP(ind+1));
            elseif Section.type == 7            % Hollowed rectangular
                tau_max(ind,1:NP(ind+1)) = (T(ind,1:NP(ind+1))./(A(1,1:NP(ind+1)) - (2.*r2(1,1:NP(ind+1)).*t(1))))./1000;
                tau(1,1:NP(ind+1)) = tau_max(1,1:NP(ind+1));   
            end
            
            if Section.type == 6 || Section.type == 7
                sigma_VM(1:7,1:NP(ind+1)) = ((sigma_NM(1:7,1:NP(ind+1)).^2) + 3.*(tau(1,1:NP(ind+1)).^2)).^0.5;
            else
            sigma_VM(1:7,1:NP(ind+1)) = ((sigma_NM(1:7,1:NP(ind+1)).^2) + 3.*(tau(1:7,1:NP(ind+1)).^2)).^0.5;
            end
        
            sigma_max_VM(ind) = max(max(abs(sigma_VM(1:7,1:NP(ind+1)))));
            
            sigma_VMmax(ind,1:NP(ind+1)) = max(abs(sigma_VM(1:7,1:NP(ind+1))));

            
            %AAAA() = [s(ind,:)' sigma_VMmax(ind,:)'];
    
        else
        
        ff = ff + 7;

        z(ff+1:ff+7,1:NP(ind+1)) = [-r1(ind,1:NP(ind+1)); -(2/3).*r1(ind,1:NP(ind+1)); -(1/3).*r1(ind,1:NP(ind+1)); zeros(1,length(r1(ind,1:NP(ind+1)))); (1/3).*r1(ind,1:NP(ind+1)); (2/3).*r1(ind,1:NP(ind+1)); r1(ind,1:NP(ind+1))];
     
        
            if Section.type == 1                % Circular
                sigma_NM(ff+1:ff+7,1:NP(ind+1)) = (z(ff+1:ff+7,1:NP(ind+1)).*M(ind,1:NP(ind+1))./J(ind,1:NP(ind+1)) + N(ind,1:NP(ind+1))./A(ind,1:NP(ind+1)))./1000;
            elseif Section.type == 2            % Ellipsoidal Pay attention, you are using the circular approximation!!!!!!!!!!!!!!!!!!!)
                sigma_NM(ff+1:ff+7,1:NP(ind+1)) = (z(ff+1:ff+7,1:NP(ind+1)).*M(ind,1:NP(ind+1))./J(ind,1:NP(ind+1)) + N(ind,1:NP(ind+1))./A(ind,1:NP(ind+1)))./1000;
            elseif Section.type == 3            % Square
                sigma_NM(ff+1:ff+7,1:NP(ind+1)) = ((z(ff+1:ff+7,1:NP(ind+1))./2).*M(ind,1:NP(ind+1))./J(ind,1:NP(ind+1)) + N(ind,1:NP(ind+1))./A(ind,1:NP(ind+1)))./1000;
            elseif Section.type == 4            % Rectangular
                sigma_NM(ff+1:ff+7,1:NP(ind+1)) = ((z(ff+1:ff+7,1:NP(ind+1))./2).*M(ind,1:NP(ind+1))./J(ind,1:NP(ind+1)) + N(ind,1:NP(ind+1))./A(ind,1:NP(ind+1)))./1000;
            elseif Section.type == 5            % Hollowed circular
                sigma_NM(ff+1:ff+7,1:NP(ind+1)) = (z(ff+1:ff+7,1:NP(ind+1)).*M(ind,1:NP(ind+1))./J(ind,1:NP(ind+1)) + N(ind,1:NP(ind+1))./A(ind,1:NP(ind+1)))./1000;
            elseif Section.type == 6            % Hollowed square
                sigma_NM(ff+1:ff+7,1:NP(ind+1)) = ((z(ff+1:ff+7,1:NP(ind+1))./2).*M(ind,1:NP(ind+1))./J(ind,1:NP(ind+1)) + N(ind,1:NP(ind+1))./A(ind,1:NP(ind+1)))./1000;
            elseif Section.type == 7            % Hollowed rectangular
                sigma_NM(ff+1:ff+7,1:NP(ind+1)) = ((z(ff+1:ff+7,1:NP(ind+1))./2).*M(ind,1:NP(ind+1))./J(ind,1:NP(ind+1)) + N(ind,1:NP(ind+1))./A(ind,1:NP(ind+1)))./1000;
            end
    
    
            if Section.type == 1                % Circular
                tau_max(ind,1:NP(ind+1)) = ((4/3).*(T(ind,1:NP(ind+1))./A(ind,1:NP(ind+1))))./1000;
                tau(ff+1:ff+7,1:NP(ind+1)) = -(tau_max(ind,1:NP(ind+1))./(r1(ind,1:NP(ind+1)).^2)).*(z(ff+1:ff+7,1:NP(ind+1)).^2) + tau_max(ind,1:NP(ind+1));
            elseif Section.type == 2            % Ellipsoidal Pay attention, you are using the circular approximation!!!!!!!!!!!!!!!!!!!)
                tau_max(ind,1:NP(ind+1)) = ((4/3).*(T(ind,1:NP(ind+1))./A(ind,1:NP(ind+1))))./1000;
                tau(ff+1:ff+7,1:NP(ind+1)) = -(tau_max(ind,1:NP(ind+1))./(r1(ind,1:NP(ind+1)).^2)).*(z(ff+1:ff+7,1:NP(ind+1)).^2) + tau_max(ind,1:NP(ind+1));
            elseif Section.type == 3            % Square
                tau_max(ind,1:NP(ind+1)) = ((3/2).*(T(ind,1:NP(ind+1))./A(ind,1:NP(ind+1))))./1000;
                tau(ff+1:ff+7,1:NP(ind+1)) = -(tau_max(ind,1:NP(ind+1))./((r1(ind,1:NP(ind+1))/2).^2)).*((z(ff+1:ff+7,1:NP(ind+1))./2).^2) + tau_max(ind,1:NP(ind+1));
            elseif Section.type == 4            % Rectangular
                tau_max(ind,1:NP(ind+1)) = ((3/2).*(T(ind,1:NP(ind+1))./A(ind,1:NP(ind+1))))./1000;
                tau(ff+1:ff+7,1:NP(ind+1)) = -(tau_max(ind,1:NP(ind+1))./((r1(ind,1:NP(ind+1))/2).^2)).*((z(ff+1:ff+7,1:NP(ind+1))./2).^2) + tau_max(ind,1:NP(ind+1));
            elseif Section.type == 5            % Hollowed circular
                tau_max(ind,1:NP(ind+1)) = (2.*(T(ind,1:NP(ind+1))./A(ind,1:NP(ind+1))))./1000;
                tau(ff+1:ff+7,1:NP(ind+1)) = -(tau_max(ind,1:NP(ind+1))./(r1(ind,1:NP(ind+1)).^2)).*(z(ff+1:ff+7,1:NP(ind+1)).^2) + tau_max(ind,1:NP(ind+1));
            elseif Section.type == 6            % Hollowed square
                tau_max(ind,1:NP(ind+1)) = (T(ind,1:NP(ind+1))./(A(ind,1:NP(ind+1)) - (2.*r1(ind,1:NP(ind+1)).*t(ind))))./1000;
                tau(ind,1:NP(ind+1)) = tau_max(ind,1:NP(ind+1));
            elseif Section.type == 7            % Hollowed rectangular
                tau_max(ind,1:NP(ind+1)) = (T(ind,1:NP(ind+1))./(A(ind,1:NP(ind+1))-(2.*r2(ind,1:NP(ind+1)).*t(ind))))./1000;
                tau(ind,1:NP(ind+1)) = tau_max(ind,1:NP(ind+1));   
            end
    
            if Section.type == 6 || Section.type == 7
                 sigma_VM(ff+1:ff+7,1:NP(ind+1)) = ((sigma_NM(ff+1:ff+7,1:NP(ind+1)).^2) + 3.*(tau(ind,1:NP(ind+1)).^2)).^0.5;
            else
                 sigma_VM(ff+1:ff+7,1:NP(ind+1)) = ((sigma_NM(ff+1:ff+7,1:NP(ind+1)).^2) + 3.*(tau(ff+1:ff+7,1:NP(ind+1)).^2)).^0.5;
            end
            
            sigma_max_VM(ind) = max(max(abs(sigma_VM(ff+1:ff+7,1:NP(ind+1)))));
            
            sigma_VMmax(ind,1:NP(ind+1)) = max(abs(sigma_VM(ff+1:ff+7,1:NP(ind+1))));
            
            %AAAA = [s(ind,:)' sigma_VMmax(ind,:)'];
        end
   end

   % Creo un vettore che mi raccolga tutte le sigma_VM senza suddividerlo
   % fra i vari archi in modo che ottengo nuovamente i punti dell'arco
   % iniziale

   %Allocazione di memoria
   sigma_VMmax_int = zeros(1,NP_1);
   kk = 0;

   for ind = 1:length(x_point)+1
       if ind == 1
            sigma_VMmax_int(kk+1:kk+NP(ind+1)) = sigma_VMmax(ind,1:NP(ind+1));   
       elseif ind  ~= 1
            kk = kk + NP(ind);
           sigma_VMmax_int(kk-(ind-3):kk+NP(ind+1)-(ind)+1) = sigma_VMmax(ind,2:NP(ind+1));
       end
   end
  %% Individuazione danneggiamento

  i = 1;

  % Determinazione della diminuzione percentuale della rigidezza della
  % sezione. Volendo è anche possibile decurtare la rigidezza in modo
  % differente
  percentuale_in(1,:) = percentuale_in(1,:) - 0.05;
  percentuale_fin(1,:) = percentuale_fin(1,:) - 0.05;


  % Condizione per uscire fuori dal ciclo è che almeno un punto si snervi
  if sigma_VMmax_int(1:NP_1)<Material.fyd
      flag = 1;
  else
      flag = 2;
  end

  %% Definizione di diminuzione della rigidezza

  %Allocazione di memoria
  J_d_in = zeros(1,length(x_point));
  J_d_fin = zeros(1,length(x_point));
  delta_in = zeros(1,length(x_point));
  delta_fin = zeros(1,length(x_point));
  decremento_in = zeros(1,length(x_point));
  decremento_fin = zeros(1,length(x_point));
  % Tiene conto del fatto che possono esservi delle molle alle estremità
  if molla_in == 1
      J_d_in = zeros(1,length(x_point)+1);
      delta_in = zeros(1,length(x_point)+1);
      decremento_in = zeros(1,length(x_point)+1);
  end
  if molla_fin == NP_1
      J_d_fin = zeros(1,length(x_point)+1);
      delta_fin = zeros(1,length(x_point)+1);
      decremento_fin = zeros(1,length(x_point)+1);
  end

  % Determinazione diminuzione di rigidezza nel caso in cui nessun punto
  % risulta snervato, tenendo conto qualora presenti delle molle alle
  % estremità

  if flag == 1

      % Presenza di molla alla base
      if molla_in == 1
      J_d_in(1) = J(1,1).*percentuale_in(1);
      delta_in(1) = (J(1,1) - J_d_in(1))./J(1,1);
      decremento_in(1) = (1-delta_in(1))./delta_in(1);
          for ind = 1:length(x_point)
              J_d_in(ind+1) = J(ind+1,1).*percentuale_in(ind+1);
              delta_in(ind+1) = ((J(ind+1,1) - J_d_in(ind+1))./J(ind+1,1));
              decremento_in(ind+1) = (1-delta_in(ind+1))./delta_in(ind+1);
          end
      end
      %Presenza di molla in mezzeria
      if molla_fin == NP_1
          J_d_fin(end) = J(end,NP(end)).*percentuale_fin(end);
          delta_fin(end) = (J(end,NP(end)) - J_d_fin(end))./J(end,NP(end));
          decremento_fin(end) = (1-delta_fin(end))./delta_fin(end);
          for ind = 1:length(x_point)
              J_d_fin(ind) = J(ind,NP(ind+1)).*percentuale_fin(ind);
              delta_fin(ind) = ((J(ind,NP(ind+1)) - J_d_fin(ind))./J(ind,NP(ind+1)));
              decremento_fin(ind) = (1-delta_fin(ind))./delta_fin(ind);
          end
      end
      %Assenza di molla in mezzeria
      if molla_fin ~= NP_1
          for ind = 1:length(x_point)
              J_d_fin(ind) = J(ind,NP(ind+1)).*percentuale_fin(ind);
              delta_fin(ind) = ((J(ind,NP(ind+1)) - J_d_fin(ind))./J(ind,NP(ind+1)));
              decremento_fin(ind) = (1-delta_fin(ind))./delta_fin(ind);
          end
      end
      %Assenza di molla alla base
      if molla_in ~= 1
          for ind = 1:length(x_point)
              J_d_in(ind) = J(ind+1,1).*percentuale_in(ind);
              delta_in(ind) = ((J(ind+1,1) - J_d_in(ind))./J(ind+1,1));
              decremento_in(ind) = (1-delta_in(ind))./delta_in(ind);
          end
      end
  end


  %% Condizione per uscire fuori dal while
  contatore = contatore + 1;
  if flag == 2 || contatore == 12
      contatore = contatore-1;
   break
  end

  %Imposizione per far continuare il ciclo while
  i=1;
end
toc
 
figure(1)
plot(s_plot,M_plot,'Linewidth',1.5);
grid on;
title('Bending moment', 'fontsize',14);
xlabel('s[m]', 'fontsize',12);
ylabel('M[kNm]', 'fontsize',12);

figure(2)
plot(s_plot,N_plot,'Linewidth',1.5);
grid on;
title('Axial force', 'fontsize',14);
xlabel('s[m]', 'fontsize',12);
ylabel('N[kNm]', 'fontsize',12);

figure(3)
plot(s_plot,T_plot,'Linewidth',1.5);
grid on;
title('Shear force', 'fontsize',14);
xlabel('s[m]', 'fontsize',12);
ylabel('T[kNm]', 'fontsize',12);

figure(4)
plot(s_plot,Dx_plot,'Linewidth',1.5);
grid on;
title('spostamento trasversale', 'fontsize',14);
xlabel('s[m]', 'fontsize',12);
ylabel('Dx[kNm]', 'fontsize',12);

figure(5)
plot(s_plot,Dz_plot,'Linewidth',1.5);
grid on;
title('Spostamento verticale', 'fontsize',14);
xlabel('s[m]', 'fontsize',12);
ylabel('Dy[kNm]', 'fontsize',12);

figure(6)
plot(s_plot,phi_plot,'Linewidth',1.5);
grid on;
title('Rotazione', 'fontsize',14);
xlabel('s[m]', 'fontsize',12);
ylabel('phi[kNm]', 'fontsize',12);


% figure(1)
% plot(x_s_1,z1_1,'o', 'Linewidth',1.5);
% grid on;
% title('s(x)', 'fontsize',14);
% xlabel('x[m]', 'fontsize',12);
% ylabel('s(x)[m]', 'fontsize',12);
% 
% figure(2)
% plot(x_s_2,z1_2,'o', 'Linewidth',1.5);
% grid on;
% title('s(x)', 'fontsize',14);
% xlabel('x[m]', 'fontsize',12);
% ylabel('s(x)[m]', 'fontsize',12);
% 
% z1 = double(subs(y,X,x_s));
% figure(3)
% plot(x_s,z1,'*','Linewidth',1.5);
% grid on;
% title('s(x)', 'fontsize',14);
% xlabel('x[m]', 'fontsize',12);
% ylabel('s(x)[m]', 'fontsize',12);
% hold on
% 
% x_s_tot = [x_s_1-Arch.Lhalf ((x_s_1(end)+x_s_2)-Arch.Lhalf)];
% z1_tot = [z1_1 (z1_1(end)+z1_2)];
% plot(x_s_tot,z1_tot,'o', 'Linewidth',1.5);
% grid on;
% title('s(x)', 'fontsize',14);
% xlabel('x[m]', 'fontsize',12);
% ylabel('s(x)[m]', 'fontsize',12);
% hold off

% figure(1)
% plot(x_1_inter, double(subs(y_1,X_1,x_1_inter)),'Linewidth',1.5);
% grid on;
% title('Arch shape', 'fontsize',14);
% xlabel('x[m]', 'fontsize',12);
% ylabel('y(x)[m]', 'fontsize',12);
 
% figure(2)
% plot(x_3_inter, double(subs(y_3,X_3,x_3_inter)),'Linewidth',1.5);
% grid on;
% title('Arch shape', 'fontsize',14);
% xlabel('x[m]', 'fontsize',12);
% ylabel('y(x)[m]', 'fontsize',12);
% 
% 
% figure(3)
% plot(x_2_inter, double(subs(y_2,X_2,x_2_inter)),'Linewidth',1.5);
% grid on;
% title('Arch shape', 'fontsize',14);
% xlabel('x[m]', 'fontsize',12);
% ylabel('y(x)[m]', 'fontsize',12);

