function grav_gui
close all;clear all
global harr
fh = figure;
uicontrol('Style','text',...
        'Position',[20 390 80 20],...
        'String','Trayecoria')
uicontrol('Style', 'popup',...
    'String','Manual|Circular|Eliptica|Parabolica|Hiperbolica',...
    'Position', [20 340 80 50],'Callback', @inputdata);
uicontrol('Style', 'pushbutton', 'String', 'Limpiar',...
    'Position', [20 20 50 20],'Callback', @limpiar);
ah1 = axes('Parent',fh,'units','pixels','Position',[120 30 340 340]);
scatter(0,0,100,[0 0.8 0],'filled')
xlim([-20 20]);ylim([-20 20]);grid on
    function inputdata(hObj,eventdata)
        val = get(hObj,'Value');
        if val==1
            [x,y] = ginput(2);
            xi=x(1);
            yi=y(1);
            vxi=0.1*(x(2)-x(1));
            vyi=0.1*(y(2)-y(1));
            % Transform from data space to figure space
            [arrowx,arrowy] = dsxy2figxy(ah1, x, y);
            harr=annotation('arrow',arrowx,arrowy);
            tf=200;
        elseif val==2
            xi=0;
            yi=10;
            vxi=sqrt(1/yi);
            vyi=0;
            tf=200;
        elseif val==3
            xi=0;
            yi=10;
            vxi=sqrt(1/yi)-0.1;
            vyi=0;
            tf=150;
        elseif val==4
            xi=-20;
            yi=20;
            alpha=pi/8;
            v=sqrt(2/(sqrt(xi^2+yi^2)));
            vxi=sin(alpha)*v;
            vyi=-cos(alpha)*v;
            tf=200;
        elseif val==5
            xi=-20;
            yi=20;
            alpha=pi/6;
            v=sqrt(2/(sqrt(xi^2+yi^2)))+0.05;
            vxi=sin(alpha)*v;
            vyi=-cos(alpha)*v;
            tf=200;
        end
        grav_k1(xi,yi,vxi,vyi,tf);
    end
    function limpiar(hObj,eventdata)
        cla;
        if exist('harr','var')
        delete(harr)%clf(harr)
        end
        ah1 = axes('Parent',fh,'units','pixels','Position',[120 30 340 340]);
        scatter(0,0,100,[0 0.8 0],'filled')
        xlim([-20 20]);ylim([-20 20]);grid on
    end
end
function grav_k1(xi,yi,vxi,vyi,tf)
%Problema 6 de la guía de teoremas de conservación

global G dt M
% Valores iniciales
G = 1;     %Constante de gravitación
M = 1 ;    %Masa terrestre
dt = 0.1; %Paso de integración [dt] = s

% Vectores
t = 0:dt:tf;         %Vector tiempo [t] = s
x = nan(length(t),4); %Vector de posición y velocidad

x(1,:) = [xi yi vxi vyi];

hold on
h1=scatter(x(1,1),x(1,2),50,[0 0 0],'filled','XDataSource','x(i,1)','YDataSource','x(i,2)');
h2=plot(x(:,1),x(:,2),'k--','XDataSource','x(:,1)','YDataSource','x(:,2)');
xlim([-20 20]);ylim([-20 20]);
axis square;grid on;
for i = 2:length(t)
    x(i,:)=RK(x(i-1,:));
    if mod(i,10)==1
        refreshdata(h1,'caller');
        refreshdata(h2,'caller');
        drawnow
    end
end
end
function a=a(x)
global G M
a(3:4)=-G*M*x(1:2)/(x(1)^2+x(2)^2)^(3/2);
a(1:2)=x(3:4);
end
function x=RK(x)
% Método de Runge-Kutta de 4º orden
global dt
k1=a(x);
x_estrella1=x+0.5*dt*k1;
k2=a(x_estrella1);
x_estrella2=x+0.5*dt*k2;
k3=a(x_estrella2);
x_estrella3=x+k3*dt;
k4=a(x_estrella3);
x= x+ dt*(k1 + 2*k2 +2*k3+k4) / 6;
end

