function yp = arco_intero_funzione()

global q0 curva variazione tipologia Restraint n_1 flag_sx flag_dx percentuale

%% %% Arch Analysis

flag_ext = 0;
percentuale = 0.5;
while  flag_ext == 0

    flag = 0;
    contatore = 1 ;
    q0 = -100;                        % Uniformly distributed load [kN/m]
    molla_in = 0;
    molla_fin = 0;
    % tengo eventalmente conto della presenza delle molle alla base 
    flag_sx = 0;
    flag_dx = 0;
    
    percentuale_in = percentuale;
    percentuale_fin = percentuale;
    
    
     while flag == 0
        %% Input parameters
        Arch.f = 30;                  % Height of the arch [m]
        Arch.L = 100;                 % Span length [m]
        n = n_1;                      % Number of intervals of the mesh [/]
        
        NP = n+1;                     % Number of points of the mesh [/]
    
       
    
        %% Material 
        Material.gamma = 78.5;        % Specific weight [kN/m3]
        Material.E = 2.1*10^8;        % Young modulus [kN/m2]
        Material_E = Material.E;
        Material.G = 76903069;        % Shear modulus [kN/m2]
        Material.fyd = 355;           % Design strength [MPa]
        
        %% Transversal section (Circular)
        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.curve = curva;
        Section.variation = variazione;
        Section.type = tipologia;
    
%     if contatore == 1
%         Arch.curve = 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 = \n Fully restrained arch = 1\n Hinged arch = 2\n");
%     end
%         
           
        %% Curve
        
        syms X
        
        if Arch.curve == 1
            % Parabola
            y = (-4*Arch.f/(Arch.L^2))*(X^2) + Arch.f;              % Parabola equation
        elseif Arch.curve ==2
            % Circumference
            y = sqrt(-X^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
        elseif Arch.curve ==3
            m = 4.0*Arch.f/Arch.L;
            y = (-4*m/(Arch.L^3) + 16*Arch.f/(Arch.L^4))*X^4 + (m/Arch.L - 8*Arch.f/(Arch.L^2))*X^2 + Arch.f;                          % 4-th order polynomial
        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)
        
        % Calculation of the incremental intervals
        x_inter = [-Arch.L/2 -Arch.L/2 + linspace(1,n,n).*Arch.L/n];
        y_inter = [0 double(subs(y,X,x_inter(2:n+1)))];
        dy_inter = y_inter(2:n+1) - y_inter(1:n);
        
        % Allocation of memory space for variables (fastening the calculation)
        s_x = zeros(1,n+1);           % Curvilinear coordinates
        
        % Calculation of s(x)
        for i = 1:n
            s_x(i+1) = s_x(i) + sqrt((Arch.L/n)^2 + dy_inter(i)^2);
        end
        
        %% Equispaced nodes along s
        ds = s_x(end)/n;
        s = ds.*(linspace(0,n,n+1));
        s_exp = [-3*ds -2*ds -1*ds s s(end)+ds s(end)+2*ds s(end)+3*ds];
        x_s = interp1(s_x,x_inter,s,'pchip','extrap');
        
        %% Curvature radius and its derivatives
        
        % Mathematical form
        R_x_math = -((1+(dy)^2)^(3/2))/d2y;          % R(x)
        
        R_s = double(subs(R_x_math,X,x_s)); % R[s(x)] [m]
        
        % Calculation of the curvature radius and its derivatives(FDM)             
        dR_s = interp1(s(1:n)+ds/2,diff(R_s)./ds,s,'pchip','extrap');
        d2R_s = interp1(s(1:n)+ds/2,diff(dR_s)./ds,s,'pchip','extrap');
        d3R_s = interp1(s(1:n)+ds/2,diff(d2R_s)./ds,s,'pchip','extrap');
        d4R_s = interp1(s(1:n)+ds/2,diff(d3R_s)./ds,s,'pchip','extrap');
        d5R_s = interp1(s(1:n)+ds/2,diff(d4R_s)./ds,s,'pchip','extrap');
        
        
        %% Transversal section properties along s(x)
        
        % Section dimensions
        if Section.variation == 1
            r1(1:NP) = Section.r_base;      % Sections radius [m]
        elseif Section.variation == 2
            r1 = ((Section.r_base - Section.r)./(s(end)./2)).*abs(s - s(end)./2) + Section.r;   % Sections radius [m] 
        elseif Section.variation == 3
            r1 = Section.r_base+(4/s(end))*(Section.r-Section.r_base).*s-(4/(s(end)^2))*(Section.r-Section.r_base).*(s.^2);   % Sections radius [m]
        end
        
        if Section.variation == 1
            r2(1:NP) = Section.r_base2;      % Sections radius [m]
        else
            r2 = ratio*r1;
        end
        
        t = Section.t;          % Thickness of the hollowed section [m]
        
        % Geometrical parameters
        if Section.type == 1        % Circular
            J = (pi.*(r1.^4))./4;                       % Inertia [m4]
            A = pi.*(r1.^2);                            % Section area [m2]
            W = pi.*(r1.^4)./(4.*r1);                   % Resistant modulus [m3]
        elseif Section.type == 2    % Ellipsoidal
            J = (pi.*(r2.*r1.^3))./4;                   % Inertia [m4]
            A = pi.*(r2.*r1);                           % Section area [m2]
            W = pi.*(r2.*r1.^3)./(4.*r2);               % Resistant modulus [m3]
        elseif Section.type == 3    % Square
            J = (r1.^4)./12;                            % Inertia [m4]
            A = r1.^2;                                  % Section area [m2]
            W = (r1.^3)./6;                             % Resistant modulus [m3]
        elseif Section.type == 4    % Rectangular
            J = r2.*(r1.^3)./12;                        % Inertia [m4]
            A = r2.*r1;                                 % Section area [m2]
            W = r2.*(r1.^2)./6;                         % Resistant modulus [m3]
        elseif Section.type == 5    % Hollowed circular
            J = (pi.*(r1.^4-(r1-t).^4))./4;             % Inertia [m4]
            A = pi.*(r1.^2-(r1-t).^2);                  % Section area [m2]
            W = (pi.*(r1.^4-(r1-t).^4))./(4.*r1);       % Resistant modulus [m3]
        elseif Section.type == 6    % Hollowed square
            J =  (r1.^4 - (r1-2*t).^4)./12;               % Inertia [m4]
            A =  r1.^2 - ((r1-2*t).^2);                     % Section area [m2]
            W =  (r1.^4 - (r1-2*t).^4)./(6.*r1);          % Resistant modulus [m3]
        elseif Section.type == 7    % Hollowed rectangle
            J =  (r2.*(r1.^3) - (r2-2*t).*((r1-2*t).^3))./12;         % Inertia [m4]
            A =  r2.*r1 - ((r2-2*t).*(r1-2*t));                           % Section area [m2]
            W =  (r2.*(r1.^3) - (r2-2*t).*((r1-2*t).^3))./(6.*r1);    % Resistant modulus [m3]
        end
        
        % Smoothing of the inertia in case of linear variation of r1
        if Section.variation == 2
            J_fit = polyfit(s,J,4);
            J = polyval(J_fit,s);
        end
        
        % Inertia and its derivatives - mathematical form
        dJ = interp1(s(1:n)+ds/2,diff(J)./ds,s,'linear','extrap');
        d2J = interp1(s(1:n)+ds/2,diff(dJ)./ds,s,'linear','extrap');
        d3J = interp1(s(1:n)+ds/2,diff(d2J)./ds,s,'linear','extrap');
        
        %% Calculation of applied loads
        alpha = atan(double(subs(dy,X,x_s)));     % Angle of inclination of each mesh segment [rad]
        
        q_sw = -Material.gamma.*A;                          % 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 = q0.*ones(1,NP);
        % 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 = (q_sw + q_ol).*cos(alpha);                     % Load normal to the arch [kN/m]
        Pt = (q_sw + q_ol).*sin(alpha);                     % Load tangent to the arch [kN/m]
        m_load = 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 = interp1(s(1:n)+ds/2,diff(Pn)./ds,s,'pchip','extrap');   % 1st derivative of normal load
        
        %% Calculation of the coefficients of the Equation of the elastic line
        
        % Partial coefficients
        a4 = R_s;
        b3 = R_s;
        b4 = 5.*dR_s;
        c2 = R_s;
        c3 = 4.*dR_s;
        c4 = 10.*d2R_s + 1./R_s;
        d1 = R_s;
        d2 = 3.*dR_s;
        d3 = 6.*d2R_s + 1./R_s;
        d4 = 10.*d3R_s - 4.*dR_s.*(1./((R_s).^2));
        e1 = 2.*dR_s;
        e2 = 3.*d2R_s + 1./R_s;
        e3 = 4.*d3R_s - 3.*dR_s.*(1./((R_s).^2));
        e4 = 5.*d4R_s - 6.*d2R_s.*(1./((R_s).^2)) + 12.*((dR_s).^2).*(1./((R_s).^3));
        f1 = d2R_s + 1./R_s;
        f2 = d3R_s - 2.*dR_s.*(1./((R_s).^2));
        f3 = d4R_s - 3.*d2R_s.*(1./((R_s).^2)) + 6.*((dR_s).^2).*(1./((R_s).^3));
        f4 = d5R_s - 4.*d3R_s.*(1./((R_s).^2)) + 24.*dR_s.*d2R_s.*(1./((R_s).^3)) - 24.*((dR_s).^3).*(1./((R_s).^4));
        g1 = -dR_s.*(1./((R_s).^2));
        g2 = -d2R_s.*(1./((R_s).^2)) + 2.*((dR_s).^2).*(1./((R_s).^3));
        g3 = -d3R_s.*(1./((R_s).^2)) + 6.*dR_s.*d2R_s.*(1./((R_s).^3)) - 6.*((dR_s).^3).*(1./((R_s).^4));
        g4 = -d4R_s.*(1./((R_s).^2)) + 8.*dR_s.*d3R_s.*(1./((R_s).^3)) + 6.*((d2R_s).^2).*(1./((R_s).^3)) -36.*((dR_s).^2).*d2R_s.*(1./((R_s).^4)) + 24.*((dR_s).^4).*(1./((R_s).^5));
        
        % Total coefficients
        a = (R_s.*J).*a4;
        b = (3.*R_s.*dJ + dR_s.*J).*b3 + (R_s.*J).*b4;
        c = (3.*R_s.*d2J + 2.*dR_s.*dJ + J./R_s).*c2 + (3.*R_s.*dJ + dR_s.*J).*c3 + (R_s.*J).*c4;
        d = (R_s.*d3J + dR_s.*d2J + dJ./R_s).*d1 + (3.*R_s.*d2J + 2.*dR_s.*dJ + J./R_s).*d2 + (3.*R_s.*dJ + dR_s.*J).*d3 + (R_s.*J).*d4;
        e = (R_s.*d3J + dR_s.*d2J + dJ./R_s).*e1 + (3.*R_s.*d2J + 2.*dR_s.*dJ + J./R_s).*e2 + (3.*R_s.*dJ + dR_s.*J).*e3 + (R_s.*J).*e4;
        f = (R_s.*d3J + dR_s.*d2J + dJ./R_s).*f1 + (3.*R_s.*d2J + 2.*dR_s.*dJ + J./R_s).*f2 + (3.*R_s.*dJ + dR_s.*J).*f3 + (R_s.*J).*f4;
        g = (R_s.*d3J + dR_s.*d2J + dJ./R_s).*g1 + (3.*R_s.*d2J + 2.*dR_s.*dJ + J./R_s).*g2 + (3.*R_s.*dJ + dR_s.*J).*g3 + (R_s.*J).*g4;
        h = -(Pt + dR_s.*Pn + R_s.*dPn + m_load./R_s)./Material.E;
    
        %Coefficients for M,N,T to solver conditions
    
    % Bending moment coefficients
    aM = R_s;
    bM = 2.*dR_s;
    cM = (d2R_s + 1./R_s);
    dM = -dR_s./(R_s).^2;
    
    % Shear force coefficients
    aT = R_s.*J;
    bT = R_s.*dJ + 3.*dR_s.*J;
    cT = 2.*dR_s.*dJ + 3.*d2R_s.*J + J./R_s;
    dT = dJ.*d2R_s + dJ./R_s + d3R_s.*J - 2.*J.*dR_s./(R_s).^2;
    eT = -(dJ.*(dR_s./(R_s).^2) + J.*(d2R_s./(R_s).^2) -2.*J.*((dR_s.^2)./(R_s).^3));
    
    % Normal force coefficients
    aN = R_s.*J;
    bN = 2.*R_s.*dJ + 4.*dR_s.*J;
    cN = d2J.*R_s + 6.*dJ.*dR_s + 6.*J.*d2R_s + J./R_s;
    dN = 2.*d2J.*dR_s + 6.*dJ.*d2R_s + 2.*dJ./R_s + 4.*J.*d3R_s - 3.*J.*dR_s./(R_s).^2;
    eN = d2J.*d2R_s + 2.*dJ.*d3R_s + J.*d4R_s + d2J./R_s - 4.*dJ.*dR_s./(R_s.^2) - 3.*J.*d2R_s./(R_s.^2) + 6.*J.*(dR_s.^2)./(R_s.^3);
    fN = -d2J.*dR_s./(R_s.^2) - 2.*dJ.*d2R_s./(R_s.^2) - J.*d3R_s./(R_s.^2) + 8.*dJ.*(dR_s.^2)./(R_s.^3) + 2.*J.*dR_s.*d2R_s./(R_s.^3) - 6.*J.*(dR_s.^3)./(R_s.^4);
        
        
        %% Displacements - Finite Differences solution
        ce1 = zeros(1,NP);  ce2 = zeros(1,NP);  ce3 = zeros(1,NP);
        ce4 = zeros(1,NP);  ce5 = zeros(1,NP);  ce6 = zeros(1,NP);
        ce7 = zeros(1,NP);
        
        matA = zeros(NP+6,NP+6);    B = zeros(1,NP+6);
        
        for i = 1:NP
        
        ce1(i) = -(1/60)*f(i)/ds + (1/90)*e(i)/ds^2 + (1/8)*d(i)/ds^3 - (1/6)*c(i)/ds^4 - (1/2)*b(i)/ds^5 + 1*a(i)/ds^6;
        ce7(i) = +(1/60)*f(i)/ds + (1/90)*e(i)/ds^2 - (1/8)*d(i)/ds^3 - (1/6)*c(i)/ds^4 + (1/2)*b(i)/ds^5 + 1*a(i)/ds^6;
        ce2(i) = +(3/20)*f(i)/ds - (3/20)*e(i)/ds^2 - 1*d(i)/ds^3 + 2*c(i)/ds^4 + 2*b(i)/ds^5 - 6*a(i)/ds^6;
        ce6(i) = -(3/20)*f(i)/ds - (3/20)*e(i)/ds^2 + 1*d(i)/ds^3 + 2*c(i)/ds^4 - 2*b(i)/ds^5 - 6*a(i)/ds^6;
        ce3(i) = -(3/4)*f(i)/ds + (3/2)*e(i)/ds^2 + (13/8)*d(i)/ds^3 - (13/2)*c(i)/ds^4 - (5/2)*b(i)/ds^5 + 15*a(i)/ds^6;
        ce5(i) = +(3/4)*f(i)/ds + (3/2)*e(i)/ds^2 - (13/8)*d(i)/ds^3 - (13/2)*c(i)/ds^4 + (5/2)*b(i)/ds^5 + 15*a(i)/ds^6;
        ce4(i) = -(49/18)*e(i)/ds^2 + (28/3)*c(i)/ds^4 - 20*a(i)/ds^6 + 1*g(i);
        
            matA(1,1:NP+6) = [ce1(1) ce2(1) ce3(1) ce4(1) ce5(1) ce6(1) ce7(1) zeros(1,(NP-1))];
        
            if i ~= 1
            matA(i,:) = [zeros(1,i-1) ce1(i) ce2(i) ce3(i) ce4(i) ce5(i) ce6(i) ce7(i) zeros(1,(NP-i))];
            end
        
            B(i) = -h(i);
        
        end
        
        %% Incastro1
        %spostamento u
        cc1_1 = 1;  %i
        %spostamento v
        cc2_1 = [(1/60)*R_s(1)/ds (-3/20)*R_s(1)/ds (3/4)*R_s(1)/ds 0 (-3/4)*R_s(1)/ds (3/20)*R_s(1)/ds (-1/60)*R_s(1)/ds];
        %rotazione phi
        if molla_in ~= 1
            cc3_1 = [(((-1/90)*R_s(1)/ds^2) + ((1/60)*dR_s(1))/ds) (((3/20)*R_s(1)/ds^2) - ((3/20)*dR_s(1))/ds) (((-3/2)*R_s(1)/ds^2) + ((3/4)*dR_s(1))/ds)...
                ((49/18)*R_s(1)/ds^2 - 1/R_s(1)) (((-3/2)*R_s(1)/ds^2) - ((3/4)*dR_s(1))/ds) (((3/20)*R_s(1)/ds^2) + ((3/20)*dR_s(1))/ds) (((-1/90)*R_s(1)/ds^2) - ((1/60)*dR_s(1))/ds)];
        elseif molla_in == 1
            M_in = +Material_E.*[J(1).*((1/8)*aM(1)/ds(1)^3 + (1/90)*bM(1)/ds(1)^2 + (-1/60)*cM(1)/ds(1))...
                                    J(1).*(-1*aM(1)/ds(1)^3 + (-3/20)*bM(1)/ds(1)^2 + (3/20)*cM(1)/ds(1))...
                                    J(1).*((13/8)*aM(1)/ds(1)^3 + (3/2)*bM(1)/ds(1)^2 + (-3/4)*cM(1)/ds(1))...
                                    J(1).*((-49/18)*bM(1)/ds(1)^2 + dM(1))...
                                    J(1).*((-13/8)*aM(1)/ds(1)^3 + (3/2)*bM(1)/ds(1)^2 + (3/4)*cM(1)/ds(1))...
                                    J(1).*(1*aM(1)/ds(1)^3 + (-3/20)*bM(1)/ds(1)^2 + (-3/20)*cM(1)/ds(1))...
                                    J(1).*((-1/8)*aM(1)/ds(1)^3 + (1/90)*bM(1)/ds(1)^2 + (1/60)*cM(1)/ds(1))]; %u(i+3)
        
                phi_in = [((-1/90)*R_s(1)/ds(1)^2 + (1/60)*dR_s(1)/ds(1))...
                        ((3/20)*R_s(1)/ds(1)^2 - (3/20)*dR_s(1)/ds(1))...
                        ((-3/2)*R_s(1)/ds(1)^2 + (3/4)*dR_s(1)/ds(1))...
                        ((49/18)*R_s(1)/ds(1)^2 - 1/R_s(1))...
                        ((-3/2)*R_s(1)/ds(1)^2 - (3/4)*dR_s(1)/ds(1))...
                        ((3/20)*R_s(1)/ds(1)^2 + (3/20)*dR_s(1)/ds(1))...
                        ((-1/90)*R_s(1)/ds(1)^2 - (1/60)*dR_s(1)/ds(1))];
                
                cc3_1 = M_in - (-phi_in.*k_in);
        end
        
        
        matA(NP+1,:) = [zeros(1,3) cc1_1 0 0 zeros(1,NP)];
        matA(NP+2,:) = [cc2_1 zeros(1,NP-1)];
        matA(NP+3,:) = [cc3_1 zeros(1,NP-1)];
        
        B(NP+1) = 0;
        B(NP+2) = 0;
        B(NP+3) = 0;
        
        %% Incastro2
        %spostamento u
        cc1_2 = 1;  %i
        %spostamento v
        cc2_2 = [(1/60)*R_s(NP)/ds (-3/20)*R_s(NP)/ds (3/4)*R_s(NP)/ds 0 (-3/4)*R_s(NP)/ds (3/20)*R_s(NP)/ds (-1/60)*R_s(NP)/ds];
        if molla_fin ~= 1
        %rotazione phi
            cc3_2 = [(((-1/90)*R_s(NP)/ds^2) + ((1/60)*dR_s(NP))/ds) (((3/20)*R_s(NP)/ds^2) - ((3/20)*dR_s(NP))/ds) (((-3/2)*R_s(NP)/ds^2) + ((3/4)*dR_s(NP))/ds)...
                ((49/18)*R_s(NP)/ds^2 - 1/R_s(NP)) (((-3/2)*R_s(NP)/ds^2) - ((3/4)*dR_s(NP))/ds) (((3/20)*R_s(NP)/ds^2) + ((3/20)*dR_s(NP))/ds) (((-1/90)*R_s(NP)/ds^2) - ((1/60)*dR_s(NP))/ds)];
        elseif molla_fin == 1
             M_fin = +Material_E.*[J(end).*((1/8)*aM(end)/ds(end)^3 + (1/90)*bM(end)/ds(end)^2 + (-1/60)*cM(end)/ds(end))...
                                    J(end).*(-1*aM(end)/ds(end)^3 + (-3/20)*bM(end)/ds(end)^2 + (3/20)*cM(end)/ds(end))...
                                    J(end).*((13/8)*aM(end)/ds(end)^3 + (3/2)*bM(end)/ds(end)^2 + (-3/4)*cM(end)/ds(end))...
                                    J(end).*((-49/18)*bM(end)/ds(end)^2 + dM(end))...
                                    J(end).*((-13/8)*aM(end)/ds(end)^3 + (3/2)*bM(end)/ds(end)^2 + (3/4)*cM(end)/ds(end))...
                                    J(end).*(1*aM(end)/ds(end)^3 + (-3/20)*bM(end)/ds(end)^2 + (-3/20)*cM(end)/ds(end))...
                                    J(end).*((-1/8)*aM(end)/ds(end)^3 + (1/90)*bM(end)/ds(end)^2 + (1/60)*cM(end)/ds(end))]; %u(i+3)
        
                phi_fin = [((-1/90)*R_s(end)/ds(end)^2 + (1/60)*dR_s(end)/ds(end))...
                        ((3/20)*R_s(end)/ds(end)^2 - (3/20)*dR_s(end)/ds(end))...
                        ((-3/2)*R_s(end)/ds(end)^2 + (3/4)*dR_s(end)/ds(end))...
                        ((49/18)*R_s(end)/ds(end)^2 - 1/R_s(end))...
                        ((-3/2)*R_s(end)/ds(end)^2 - (3/4)*dR_s(end)/ds(end))...
                        ((3/20)*R_s(end)/ds(end)^2 + (3/20)*dR_s(end)/ds(end))...
                        ((-1/90)*R_s(end)/ds(end)^2 - (1/60)*dR_s(end)/ds(end))];
                
                cc3_2 = M_fin - (phi_fin.*k_fin);
        end
        
        matA(NP+4,:) = [zeros(1,NP) 0 0 cc1_2 zeros(1,3)];
        matA(NP+5,:) = [zeros(1,NP-1) cc2_2];
        matA(NP+6,:) = [zeros(1,NP-1) cc3_2];
        
        B(NP+4) = 0;
        B(NP+5) = 0;
        B(NP+6) = 0;
        
        
        x_res = matA\B';
        u_exp = x_res';
        
        
        % Calculation of the tangential displacements and its derivatives
        u = u_exp(4:NP+3);
        du = ((-1/60).*u_exp(1:NP) + (3/20).*u_exp(2:NP+1) + (-3/4).*u_exp(3:NP+2) + (0).*u_exp(4:NP+3) + (3/4).*u_exp(5:NP+4) + (-3/20).*u_exp(6:NP+5) + (1/60).*u_exp(7:NP+6))./ds;
        d2u = ((1/90).*u_exp(1:NP) + (-3/20).*u_exp(2:NP+1) + (3/2).*u_exp(3:NP+2) + (-49/18).*u_exp(4:NP+3) + (3/2).*u_exp(5:NP+4) + (-3/20).*u_exp(6:NP+5) + (1/90).*u_exp(7:NP+6))./(ds^2);
        d3u = ((1/8).*u_exp(1:NP) + (-1).*u_exp(2:NP+1) + (13/8).*u_exp(3:NP+2) + (0).*u_exp(4:NP+3) + (-13/8).*u_exp(5:NP+4) + (1).*u_exp(6:NP+5) + (-1/8).*u_exp(7:NP+6))./(ds^3);
        d4u = ((-1/6).*u_exp(1:NP) + (2).*u_exp(2:NP+1) + (-13/2).*u_exp(3:NP+2) + (28/3).*u_exp(4:NP+3) + (-13/2).*u_exp(5:NP+4) + (2).*u_exp(6:NP+5) + (-1/6).*u_exp(7:NP+6))./(ds^4);
        d5u = ((-1/2).*u_exp(1:NP) + (2).*u_exp(2:NP+1) + (-5/2).*u_exp(3:NP+2) + (0).*u_exp(4:NP+3) + (5/2).*u_exp(5:NP+4) + (-2).*u_exp(6:NP+5) + (1/2).*u_exp(7:NP+6))./(ds^5);
        
        % Perpendicular displacements [m]
        v = -du.*R_s;
        
        % Rotations [rad]
        phi = -(d2u.*R_s + du.*dR_s + u./R_s);
        
        % Deformation angle and its derivatives
        k = d3u.*d1 + d2u.*e1 + du.*f1 + u.*g1;
        dk = c2.*d4u + d2.*d3u + e2.*d2u + f2.*du + g2.*u;
        d2k = b3.*d5u + c3.*d4u + d3.*d3u + e3.*d2u + f3.*du + g3.*u; 
        
        
        %% Stabilization of d2k in case of discontinuous overload
        for i=1:NP-1
            if q_ol(i) ~= 0
                if q_ol(i+1) == 0
                    d2k(i) = interp1(s(1:i-1), d2k(1:i-1), s(i), 'pchip','extrap');
                    d2k(i+1) = interp1(s(i+2:NP), d2k(i+2:NP), s(i+1), 'pchip','extrap');
                end
            else
                if q_ol(i+1) ~= 0
                    d2k(i) = interp1(s(1:i-1), d2k(1:i-1), s(i), 'pchip','extrap');
                    d2k(i+1) = interp1(s(i+2:NP), d2k(i+2:NP), s(i+1), 'pchip','extrap');
                end
            end            
        end
        
        
        %% Solicitations
        M = Material.E.*J.*k;                                           % Bending Moment [kNm]
        T = Material.E.*-(dJ.*k + J.*dk) + m_load;                           % Shear [kN]
        N = R_s.*Pn - R_s.*Material.E.*(d2J.*k + 2.*dJ.*dk + J.*d2k);       % Axial force [kN]
        
        Dz = -(v.*cos(alpha)+u.*sin(alpha));
        Dx = -(-v.*sin(alpha)+u.*cos(alpha));
        
        Results = [Dx' Dz' phi' M' N' T'];
        z1 = double(subs(y,X,x_s));
        
        %% Stresses
        sigma = (M./W + N./A)./1000;                                    % Normal stress [MPa]
        
        sigma_max = max(abs(sigma));
        
        constraint = sigma_max/Material.fyd - 1;                                 % Strength verification
        
        
        %% Von Mises stress calculation
        z = [-r1; -(2/3).*r1; -(1/3).*r1; zeros(1,length(r1)); (1/3).*r1; (2/3).*r1; r1];
        
        if Section.type == 1                % Circular
            sigma_NM = (z.*M./J + N./A)./1000;
        elseif Section.type == 2            % Ellipsoidal Pay attention, you are using the circular approximation!!!!!!!!!!!!!!!!!!!)
            sigma_NM = (z.*M./J + N./A)./1000;
        elseif Section.type == 3            % Square
            sigma_NM = ((z./2).*M./J + N./A)./1000;
        elseif Section.type == 4            % Rectangular
            sigma_NM = ((z./2).*M./J + N./A)./1000;
        elseif Section.type == 5            % Hollowed circular
            sigma_NM = (z.*M./J + N./A)./1000;
        elseif Section.type == 6            % Hollowed square
            sigma_NM = ((z./2).*M./J + N./A)./1000;
        elseif Section.type == 7            % Hollowed rectangular
            sigma_NM = ((z./2).*M./J + N./A)./1000;
        end
        
        if Section.type == 1                % Circular
            tau_max = ((4/3).*(T./A))./1000;
            tau = -(tau_max./(r1.^2)).*(z.^2) + tau_max;
        elseif Section.type == 2            % Ellipsoidal Pay attention, you are using the circular approximation!!!!!!!!!!!!!!!!!!!)
            tau_max = ((4/3).*(T./A))./1000;
            tau = -(tau_max./(r1.^2)).*(z.^2) + tau_max;
        elseif Section.type == 3            % Square
            tau_max = ((3/2).*(T./A))./1000;
            tau = -(tau_max./((r1/2).^2)).*((z./2).^2) + tau_max;
        elseif Section.type == 4            % Rectangular
            tau_max = ((3/2).*(T./A))./1000;
            tau = -(tau_max./((r1/2).^2)).*((z./2).^2) + tau_max;
        elseif Section.type == 5            % Hollowed circular
            tau_max = (2.*(T./A))./1000;
            tau = -(tau_max./(r1.^2)).*(z.^2) + tau_max;
        elseif Section.type == 6            % Hollowed square
            tau_max = (T./(A - (2.*r1.*t)))./1000;
            tau = tau_max;
        elseif Section.type == 7            % Hollowed rectangular
            tau_max = (T./(A-(2.*r2.*t)))./1000;
            tau = tau_max;   
        end
        
        sigma_VM = ((sigma_NM.^2) + 3.*(tau.^2)).^0.5;
        
        sigma_max_VM = max(max(abs(sigma_VM)));
        
        sigma_VMmax = max(abs(sigma_VM));
        
        AAAA = [s' sigma_VMmax'];
        
        
        
        i = 1;
        yp = [];
        for j = 1:NP
            if sigma_VMmax(j) > Material.fyd
                yp(i) = j;
                i = i+1;
            end
        end
    
        %% Pulizia vettore punti danneggiati
        
        % I punti in esame devono distare almeno due nodi per poter garantire le
        % interpolazioni successive. In questo caso devo eliminare i punti
        % successivi poichè sono comunque inutili in quanto verrebbero
        % eliminati, per il codice successivo.
        xx = length(yp);
        i=1;
        j = 1;
        while i < xx
            if yp(i+1) == yp(i)+1
                if yp(i+1) == NP
                    yp(i) = [];
                    xx = xx-1;
                else
                   yp(i+1) = [];
                   xx = xx-1;
                end   
            end
            i=i+1;
            j=j+1;
            if j == length(yp)
                break
            end   
        end
    
        if isempty(yp) == 1
            flag = 0;
        else
            % caso in cui si ha una molla elastica in una delle due basi o
            % anche entrambe
            if length(yp) == 2 && yp(1) == 1 && yp(end) == NP
                molla_in = 1;
                molla_fin = 1;
                flag = 0;
            elseif length(yp) == 1 && yp(end) == NP
                molla_fin = 1;
                flag = 0;
            elseif length(yp) == 1 && yp(1) == 1
                molla_in = 1;
                flag = 0;
            else
                flag = 1;
            end
        end
        
        % una volta inserite le molle ai piedi si avrà una ridistribuzione
        % delle tensioni. Motivo per il quale in corrsipondenza di tali
        % punti potrei non riavere delle tensioni di snevramento e quindi
        % nel vettore yp non ri ripresentano i valori 1 ed NP --->
        % inserisco un flag che mi dice se le molle in corrispondenza delle
        % reni dell'arco sono state inserite o meno
        if  molla_in == 1
            flag_sx = 1;
        end
        if molla_fin == 1
            flag_dx = 1;
        end
    
        if flag == 0
             q0 = q0-2;
        end
    
        contatore = contatore+1;
        
        %% Determinazione delle rigidezze
        if Restraint == 1
        phi_d = 0.016*(pi);                                                    % [rad]
        elseif Restraint == 2
        phi_d = 0.067*(pi);                                                    % [rad]
        end
    
        L_d_in = phi_d.*(R_s(1));                                              % lunghezza tratto danneggiato [m]
        k_in = Material_E*J(1)/L_d_in;                                         % Rigidezza iniziale [kNm/rad]
        L_d_fin = phi_d.*(R_s(end));                                           % lunghezza tratto danneggiato [m]
        k_fin = Material_E*J(end)/L_d_fin;                                     % Rigidezza iniziale [kNm/rad]
    
    
        J_d_in = J(1).*percentuale_in;
        delta_in = (J(1) - J_d_in)./J(1);
        decremento_in = (1-delta_in)./delta_in;
    
        J_d_fin = J(end).*percentuale_fin;
        delta_fin = (J(end) - J_d_fin)./J(end);
        decremento_fin = (1-delta_fin)./delta_fin;
    
        k_in = k_in.*decremento_in;
        k_fin = k_fin.*decremento_fin;
        
        
        
        %% Plots
        % figure(1)
        % plot(x_s,z1,'o', 'Linewidth',1.5);
        % grid on;
        % title('s(x)', 'fontsize',14);
        % xlabel('x[m]', 'fontsize',12);
        % ylabel('s(x)[m]', 'fontsize',12);
        % %close 2
        
        % figure(1)
        % plot(x_inter, double(subs(y,X,x_inter)),'Linewidth',1.5);
        % grid on;
        % title('Arch shape', 'fontsize',14);
        % xlabel('x[m]', 'fontsize',12);
        % ylabel('y(x)[m]', 'fontsize',12);
        % %close 1
        % 
        % figure(2)
        % plot(x_s,s,'o', 'Linewidth',1.5);
        % grid on;
        % title('s(x)', 'fontsize',14);
        % xlabel('x[m]', 'fontsize',12);
        % ylabel('s(x)[m]', 'fontsize',12);
        % %close 2
        % 
        % figure(3)
        % plot(s,x_s,'o','Linewidth',1.5);
        % grid on;
        % title('x(s)', 'fontsize',14);
        % xlabel('s[m]', 'fontsize',12);
        % ylabel('x(s)[m]', 'fontsize',12);
        % %close 3
        % 
        % figure(4)
        % plot(s,R_s,'Linewidth',1.5);
        % grid on;
        % title('Curvature Radius R(s)', 'fontsize',14);
        % xlabel('s[m]', 'fontsize',12);
        % ylabel('R(s)[m]', 'fontsize',12);
        % %close 4
        % 
        % figure(5)
        % plot(s,J,'Linewidth',1.5);
        % grid on;
        % title('Sectional Inertia I(s)', 'fontsize',14);
        % xlabel('s[m]', 'fontsize',12);
        % ylabel('I(s)[m^4]', 'fontsize',12);
        % %close 5
        % 
        % figure(6)
        % plot(s,Pn,'Linewidth',1.5);
        % hold on
        % plot(s,Pt,'Linewidth',1.5);
        % grid on;
        % title('Applied load', 'fontsize',14);
        % xlabel('s[m]', 'fontsize',12);
        % ylabel('Load[kN]', 'fontsize',12);
        % legend('Pn','Pt','fontsize',10);
        % %close 6
        % 
        % figure(7)
        % plot(s,u,'Linewidth',1.5);
        % grid on;
        % title('Tangential displacements', 'fontsize',14);
        % xlabel('s[m]', 'fontsize',12);
        % ylabel('u[m]', 'fontsize',12);
        % 
        % figure(8)
        % plot(s,v,'Linewidth',1.5);
        % grid on;
        % title('Perpendicular displacements', 'fontsize',14);
        % xlabel('s[m]', 'fontsize',12);
        % ylabel('v[m]', 'fontsize',12);
        % 
        % figure(9)
        % plot(s,phi,'Linewidth',1.5);
        % grid on;
        % title('Rotations', 'fontsize',14);
        % xlabel('s[m]', 'fontsize',12);
        % ylabel('phi[rad]', 'fontsize',12);
        % 
        % figure(10)
        % plot(s,M,'Linewidth',1.5);
        % grid on;
        % title('Bending moment', 'fontsize',14);
        % xlabel('s[m]', 'fontsize',12);
        % ylabel('M[kNm]', 'fontsize',12);
        % 
        % figure(11)
        % plot(s,T,'Linewidth',1.5);
        % grid on;
        % title('Shear', 'fontsize',14);
        % xlabel('s[m]', 'fontsize',12);
        % ylabel('T[kN]', 'fontsize',12);
        % 
        % figure(12)
        % plot(s,N,'Linewidth',1.5);
        % grid on;
        % title('Axial force', 'fontsize',14);
        % xlabel('s[m]', 'fontsize',12);
        % ylabel('N[kN]', 'fontsize',12);
        % 
        % figure(13)
        % plot(s,sigma,'Linewidth',1.5);
        % grid on;
        % title('Maximum normal stress for each section', 'fontsize',14);
        % xlabel('s[m]', 'fontsize',12);
        % ylabel('\sigma[MPa]', 'fontsize',12);
        % 
        % figure(14)
        % plot(s,dJ,'Linewidth',1.5);
        % grid on;
        % title('1st derivative of J', 'fontsize',14);
        % xlabel('s[m]', 'fontsize',12);
        % ylabel('dJ', 'fontsize',12);
        % 
        % figure(15)
        % plot(s,d2J,'Linewidth',1.5);
        % grid on;
        % title('2nd derivative of J', 'fontsize',14);
        % xlabel('s[m]', 'fontsize',12);
        % ylabel('d2J', 'fontsize',12);
        % 
        % figure(16)
        % plot(s,k,'Linewidth',1.5);
        % grid on;
        % title('k', 'fontsize',14);
        % xlabel('s[m]', 'fontsize',12);
        % ylabel('k', 'fontsize',12);
        % 
        % figure(17)
        % plot(s,dk,'Linewidth',1.5);
        % grid on;
        % title('1st derivative of k', 'fontsize',14);
        % xlabel('s[m]', 'fontsize',12);
        % ylabel('dk', 'fontsize',12);
        % 
        % figure(18)
        % plot(s,d2k,'Linewidth',1.5);
        % grid on;
        % title('2nd derivative of k', 'fontsize',14);
        % xlabel('s[m]', 'fontsize',12);
        % ylabel('d2k', 'fontsize',12);
        
        
        % fy_plot_minus = -355*ones(1,length(s));
        % fy_plot_plus = 355*ones(1,length(s));
        % 
        % figure(13)
        % plot(s,r1,'Linewidth',1.5);
        % % hold on
        % % plot(s,fy_plot_plus, 'r', 'Linewidth',1.5);
        % % hold on
        % % plot(s,fy_plot_minus, 'r','Linewidth',1.5);
        % grid on;
        % title('Maximum normal stress for each section', 'fontsize',14);
        % xlabel('s[m]', 'fontsize',12);
        % ylabel('\sigma[MPa]', 'fontsize',12);
        
        %toc
     end

     for i = 1:length(sigma_VMmax)
         if sigma_VMmax(i) > Material.fyd + 5
             flag_ext = 0;
             percentuale = percentuale - 0.1;
             break
         else
             flag_ext = 1;
         end
     end
end
end