clear
clc
close all

%% Cargo el archivo
fid = fopen('StereoCat Frameseries 3D Height Plot Data Saved 2020-07-06T0621Z.txt');
a=textscan(fid,'%q %q %q %q %q %q','Delimiter',',','HeaderLines',2); %Guardo info de cada columna
fclose(fid);

%% Convierto la fecha a un array numerico
for i = 1:3
    for j = 1:length(a{2*i-1})
        if strcmp(a{2*i-1}{j},'') == 0 %Comparo el elemento del array con esapcio vacio
            hora{i}(j) = datenum(a{2*i-1}{j},'yyyy-mm-ddTHH:MM:SSZ');
        else
            break
        end
    end
end

%% Convierto distancias al sol a un array numerico
for i = 1:3
    for j = 1:length(a{2*i})
        if isnan(str2double(strtok(a{2*i}{j},'Rs'))) == 0 %Veo si el valor de la distancia es un numero o no(espacio vacio) (NaN)
            h{i}(j) = str2double(strtok(a{2*i}{j},'Rs'));
        else
            break
        end
    end
end

%% Convierto unidades: distancia a km y tiempo a segundos
rs = 6.957e5; %radio sol en km
dia = 86400; %segundos en un dia
for i = 1:3
    h{i} = h{i}*rs; %R en km
    horanor{i} = (hora{i}-hora{i}(1))*dia; %tiempo normalizado
end

%% Ajusto los datos por una parabola para hallar aceleracion y velocidad inicial
for i = 1:3
    p = polyfit(horanor{i},h{i},2);
    ac(i) = 2*p(1); %km/s2
    vel(i) = p(2); %km/s
    x0(i) = p(3); %km
end

%% Hallo velocidad a 20Rs y hora de eyección

for i = 1:3
    v20(i) = sqrt(2*ac(i)*(20*rs-x0(i))+vel(i)^2);
    t1 = (-vel(i)+sqrt(vel(i)^2-2*ac(i)*(x0(i)-rs)))/(ac(i)); % tiempo1 en la eyeccion
    t2 = (-vel(i)-sqrt(vel(i)^2-2*ac(i)*(x0(i)-rs)))/(ac(i)); % tiempo2 en la eyeccion (este no xq da veloc negativa)
    ey1 = hora{i}(1)+t1/dia;
    ey2 = hora{i}(1)+t2/dia;
    f1=datevec(ey1);
    
    msj = sprintf('Ajuste %d:',i);
    disp(msj)
    fprintf('\n')
    fprintf('Velocidad a 20Rs: %6.2f km/s \n', v20(i))
    fprintf('Fecha de eyeccion: %i-%i-%i %i:%i:%4.2f \n',f1(1),f1(2),f1(3),f1(4),f1(5),f1(6))
    fprintf('\n') 
    fprintf('\n')
end