function so = CompStdOrb(so, pe, stTgd)

%CompStdOrb computes polynomial coefficients for all satellites available in precise
%ephemeris peEph
%so = CompStdOrb(so, pe)
%so - Standard orbit (class StdOrb)
%pe - Coordinate ephemeris (class CoordEph)
%stTgd = struct('vDateTime', dtTime, 'cTGD',0) - this structure is created by getTgd() function
%Coordinates for an epoch are computed by method GetSatCoord

%Written by Milan Horemuz, last modified 2005-02-03



x = get(pe,'dX'); %extract x-coordinate for all satellites
y = get(pe,'dY'); 
z = get(pe,'dZ'); 
t = get(pe,'dTime');  %epochs 
cl = get(pe,'dDts');   %clock corrections
prn = get(pe,'iPRN'); %satellie numbers
ne = length(t);     %number of epochs
tb = t(1);          %first epoch, for which the coordinates are given
te = t(ne);  %last epoch
tm = (so.dFitInt - so.dInterval)/2; %lead-in and lead-out interval
Nint = floor((te - tb - 2*tm)/so.dInterval); %number of fit intervals
m = 1;
for i = 1:Nint
    so.dTimeSet(i) = tb + tm + (i-1)*so.dInterval; %time of start of interval
    tg = FindTGD(stTgd, so.dTimeSet(i));  %find value of Tgd closest to time so.dTimeSet(i) for all satellites
    [m, n] = FindInt(t, so.dTimeSet(i)-tm, so.dFitInt, m); %find vector interval (indexes m,n) that corresponds to the fit interval 
    if n < 0     %it was not possible to find data for current interval
        so.dTimeSet(i) = [];
        break;
    end
    tt = (t(m:n) - so.dTimeSet(i));  %reduce time to get better numerical stability
    for j = 1:size(x,1) %loop over satellites
        if prn(j, n) < 1  %no data for satellite j
            coef = zeros(1,size(so.dXCoef,3));
            so.dXCoef(j,i,:) =  coef;   %placeholder for satellite
            so.dYCoef(j,i,:) =  coef;
            so.dZCoef(j,i,:) =  coef;
            continue;
        end
        [coef, S, mu] = polyfit(tt, x(j,m:n), so.iPolyOrder);
        %yy=polyval(coef,tt,[],mu);
        %dd= yy-x(j,m:n);
        %plot(dd);
        so.dXCoef(j,i,:) =  coef;
        so.dXmu(j,i,:) =  mu';
        [coef, S, mu] = polyfit(tt, y(j,m:n), so.iPolyOrder);
        so.dYmu(j,i,:) =  mu';
        so.dYCoef(j,i,:) =  coef;
        [coef, S, mu] = polyfit(tt, z(j,m:n), so.iPolyOrder);
        so.dZCoef(j,i,:) =  coef;
        so.dZmu(j,i,:) =  mu';
        so.dTGD(j,i) = tg(j);  
    end
end
%interval, for which the std. orbit is valid
so.dTbeg = so.dTimeSet(1); %begin
so.dTend = so.dTimeSet(length(so.dTimeSet)) + so.dInterval;  %end
so.iPRN = prn;
so.dtClTime = t;
so.dClErr = cl;