twolines2kepler(file, sat)

function ephem = twolines2kepler(file, sat)

% twolines2kepler transforms a "two lines" file in Kepler elements
% see http://www.movingsatellites.com/tle-fmt.htm

% Right ascension = Om_0
% Inclinatoin = i
% Argument of perigee = om
% mean semi-major axis = a
% eccentricity = e
% reference epoch = ref_time
% mean anomaly = M_0 

galileodata=fopen(file);

if strcmp(sat,'giovea') %getting ephemerides for Giove-A (lines 2-3)
    fgetl(galileodata);
    data1=fgetl(galileodata);
    data2=fgetl(galileodata);
else                    %getting ephemerides for Giove-B (lines 5-6)
    fgetl(galileodata); fgetl(galileodata); fgetl(galileodata);fgetl(galileodata);
    data1=fgetl(galileodata);
    data2=fgetl(galileodata);
end;

fclose(galileodata);
clear galileodata;

format('long');

% read content
epochyear=data1(1,19:20);
epochday=data1(1,21:32);
i=data2(1,9:16);
Om_0=data2(1,18:25);
e=data2(1,27:33);
om=data2(1,35:42);
M_0=data2(1,44:51);
mean_motion=data2(1,53:63);
d_mean_motion=data1(1,34:43);
dd_mean_motion=data1(1,46:50);
dd_mean_motion_power=data1(1,52);

% add a decimal point
e=[46,e];
dd_mean_motion=[46,dd_mean_motion];

% transform to numeric values
epochyear=sscanf(epochyear,'%f');
epochday=sscanf(epochday,'%f');
i=sscanf(i,'%f');
Om_0=sscanf(Om_0,'%f');
e=sscanf(e,'%f');
om=sscanf(om,'%f');
M_0=sscanf(M_0,'%f');
mean_motion=sscanf(mean_motion,'%f');
d_mean_motion=sscanf(d_mean_motion,'%f');
dd_mean_motion=sscanf(dd_mean_motion,'%f');
dd_mean_motion_power=sscanf(dd_mean_motion_power,'%f');
dd_mean_motion=dd_mean_motion*10^dd_mean_motion_power;

dOm_0=0;    %"Rate of change in the measurement of the angle of right ascension as defined in the Right Ascension mnemonic" not used in two line element sets

% transform to rad
i=i/180*pi;
Om_0=Om_0/180*pi;
om=om/180*pi;
M_0=M_0/180*pi;

% Calculates missing kepler elements
GM=3986004.418e8;
n=(2*pi*mean_motion)/(86400);
a=(GM/(n^2))^(1/3);

% Calculates the reference epoch in matlab format 
%(year month day hour minute second)
epochyear=2000+epochyear;
ref_time=datenum(epochyear,1,1)+epochday-1;   


ephem.Om_0 = Om_0; 
ephem.i = i;
ephem.om = om;
ephem.a = a;
ephem.e = e;
ephem.ref_time = ref_time ;
ephem.M_0 = M_0;
ephem.dOm_0 = dOm_0;