calc_satposition_ecef(acttime, sat, ephem)

function satposition_ecef = calc_satposition_ecef(acttime, sat, ephem)

% Erdparameter G*Me aus WGS-84
GM = 3986004.418e8;
% Winkelgeschwindigkeit der Erde
earth_rot = 7292115.1467e-11;



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

ref_time=datevec(ref_time);
time_in_epoch=etime(acttime,ref_time);

% Mean anomaly berechnen
mean_anomaly = (M_0+sqrt(GM/(a^3))*time_in_epoch); 

% Newton Verfahren zur Berechnung der eccentric anomaly ( E = M + e * sin(E) )
% Initialisierung des Newton Verfahrens
eccentric_anomaly = mean_anomaly + e*sin(mean_anomaly);

for j = 1:7
    % Funktionswert berechnen
    q =  mean_anomaly + e*sin(eccentric_anomaly)-eccentric_anomaly; 
    % Ableitung berechnen
    qdiff =  e*cos(eccentric_anomaly)-1;
    % Newton Korrekturschritt
    eccentric_anomaly = eccentric_anomaly - q/qdiff;
end

% Satposition in der Ellipsenebene berechnen [x; y; 0]
sat_pos_ell =  a*[cos(eccentric_anomaly)-e; sqrt(1-e^2)*sin(eccentric_anomaly);0];


l =  Om_0-siderealtime(acttime);

Rot_l = [cos(-l), sin(-l), 0;
         -sin(-l), cos(-l), 0;
        0, 0, 1];

Rot_i = [1, 0, 0;
         0, cos(-i), sin(-i);
         0, -sin(-i), cos(-i)];
     
Rot_w = [cos(-om), sin(-om), 0;
         -sin(-om), cos(-om), 0;
         0, 0, 1];     
         
         
% Position in ECEF-Koordinaten berechnen
satposition_ecef = Rot_l*Rot_i*Rot_w*sat_pos_ell;