// ********************
// Initialisation
// ********************

clear;
clc;

// ********************
// Fonctions
// ********************

// Saisie des paramètres géométriques

function [A]=param_geo()
    A(1) = input("grand diamètre : D = ");
    A(2) = input("petit diamètre : d = ");
    A(3) = input("rayon : r = ");
endfunction

// Calcul en traction

function [Kt] = Kt_traction(Alpha, Beta)
    drapeau = Beta <= 2;
    select drapeau
    case %t then
        coefs = [0.926, 1.157, -0.099;
            0.012, -3.036, 0.961;
            -0.302, 3.977, -1.744;
            0.365, -2.098, 0.878];
    case %f then
        coefs = [1.200, 0.860, -0.022;
            -1.805, -0.346, -0.038;
            2.198, -0.486, 0.165;
            -0.593, -0.028, -0.106];
    end
    C = coefs(:,1) + coefs(:,2)*sqrt(Beta) + coefs(:,3)*Beta;
    polynome = poly(C', "x", "c") ;
    Kt = horner(polynome, Alpha);
endfunction

// ********************
// Programme principal
// ********************

// Accueil

disp("Facteurs de concentration de contrainte d''un arbre épaulé")

// A(1) = D ; A(2) = d ; A(3) = r

A = param_geo();
t = 0.5*(A(1) - A(2)); // hauteur de filet
Beta = t/A(3); // t/r
Alpha = 2*t/A(1);

Kt_tr = Kt_traction(Alpha, Beta);

if (Beta < 0.1) | (Beta > 20) then
    disp("Attention, valeur extrapolée");
end
disp('Kt traction = '+string(Kt_tr));