clc;
h=[0,3,6,9,12,15,18,21,24,27,30,33]; %Altura en km
d=[1.91000000000000,1.21000000000000,0.770000000000000,0.470000000000000,0.310000000000000,0.190000000000000,0.120000000000000,0.0750000000000000,0.0460000000000000,0.0290000000000000,0.0180000000000000,0.0110000000000000]; %Densidad en kg/m3

%plot(l,p,figura,tamaño de la figura,tamaño,color de relleno, color, color del borde, color)
plot(h,d,'o','MarkerSize',4,'MarkerFaceColor','r','MarkerEdgeColor','r')
%suma de todas los elementos de un vector
format long
A=[length(h), sum(h), sum(h.^2); sum(h), sum(h.^2), sum(h.^3);sum(h.^2), sum(h.^3), sum(h.^4)];
B=[sum(d);sum(d.*h);sum(d.*h.^2)];
x=inv(A)*B; %valores a0 a1 y a2
%extraccion de coeficientes del modelo
a_0=x(1);
a_1=x(2);
a_2=x(3);

syms ll
%POLINOMIO DE GRADO 2
P_2=a_2*ll^2+a_1*ll+a_0;
vpa(P_2,5)
%grado 3
%suma de todas los elementos de un vector
format long
A=[length(h), sum(h), sum(h.^2), sum(h.^3); sum(h), sum(h.^2), sum(h.^3), sum(h.^4);sum(h.^2), sum(h.^3), sum(h.^4), sum(h.^5);sum(h.^3), sum(h.^4), sum(h.^5), sum(h.^6)];
B=[sum(d);sum(d.*h);sum(d.*h.^2);sum(d.*h.^3)];
x=inv(A)*B; %valores a0 a1 y a2
%extraccion de coeficientes del modelo
a_0=x(1);
a_1=x(2);
a_2=x(3);
a_3=x(4);
syms ll
%POLINOMIO DE GRADO 3
P_3=a_3*ll^3+a_2*ll^2+a_1*ll+a_0;
vpa(P_3,5)
%GRAFICANDO
hold on
aa=ezplot(P_2,[0,max(h)]);
set(aa,'Color','g');
a2=ezplot(P_3,[0,max(h)]);
set(a2,'Color','b');
legend('Pts. experimentales', 'P_2(ll)', 'P_3(ll)');
% r^2 POLINOMIO DE SEGUNDO GRADO
sr=vpa(sum((d-subs(P_2,h)).^2),5)
st=vpa(sum((d-mean(d)).^2),5)
r=sqrt((st-sr)/st)
sr=vpa(sum((d-subs(P_3,h)).^2),5)
st=vpa(sum((d-mean(d)).^2),5)
r=sqrt((st-sr)/st)
disp('densidad del aire en la Cdmx:')
disp(vpa(subs(P_3,2.25)))
disp('densidad del aire en el Everest:')
disp(vpa(subs(P_3,1.2)))

%el de grado tres numericamente es el mejor
%pero decidimos con respecto a la grafica y a numericamente.