function [M_i, M_e] = calcola_M(Ai_At,Ae_At)

gamma = 1.4;
a = gamma + 1;
b = gamma - 1;
c = 0.5*(gamma - 1);
d = sqrt(gamma);

F_1 = d*(1/(1+c))^(a/(2*b));

fun = @(M_i) F_1/( d*M_i*(1/(1+c*M_i^2))^(a/(2*b)) ) - Ai_At;
fun2 = @(M_e) F_1/( d*M_e*(1/(1+c*M_e^2))^(a/(2*b)) ) - Ae_At;

%% M_I

M_i0 = 0.12; % valore iniziale M_i
options = optimset('Display','iter','Tolfun',1e-12,'TolX',1e-12);

M_i = fsolve(fun, M_i0, options);

%% M_e

M_e0 = 5; % valore iniziale M_e
options = optimset('Display','iter','Tolfun',1e-12,'TolX',1e-12);
M_e = fsolve(fun2, M_e0, options);



