visualcompute.m

In the planning stage of a GNSS-measurement campaign it is useful to compute several design para-met

        j=1;
        for i=2:length(ephtmp)
            if ephtmp(i,1)==ephtmp(i-1,1),
                continue
            else
                j=j+1;
                eph(j,:)=ephtmp(i,:);
            end
        end
    end
    i          = find(eph(:,27)==0); % find healthy satellites
    eph        = eph(i,:);
    [xs,ys,zs] = cpaziele (tsat,eph);
end

% cutoff elevation and number of epochs
cutoff     = str2double(get(findobj(gcf,'Tag','cutoff'),'String'));

% testing parameter (default)
lam0 = 17.0747;


% MDB/BNR/MDE
if ismember(out,[1 2 4 5 6])
    % observation for which to compute output:
    % carrier slip or code outlier
    if (get(findobj(gcf,'Tag','errtype'),'Value'))==1
        output = 'slip';
    else
        output = 'outlier';
    end
    outfreq = (get(findobj(gcf,'Tag','outputfreq'),'Value'));
    if (SYSTEM==4) && (outfreq > 3)
        outputfreq = find(freq==frq(outfreq+1));
        if length(outputfreq)>1, outputfreq=outputfreq(2); end
    elseif SYSTEM==3
        outputfreq = find(freq==frq(outfreq+4));
    else
        outputfreq = find(freq==frq(outfreq));
        outputfreq = outputfreq(1);
    end

    if (outputfreq > nfreq) && (SYSTEM ~= 4)
        h2 = msgbox ('Error: output cannot be computed with the number of frequencies as given');
        watchoff(figNumber);
    elseif isempty(outputfreq)
        h2 = msgbox ('Error: wrong frequency choice for outlier / cycle slip');
        watchoff(figNumber);
    elseif (outputfreq > nfreq) && (SYSTEM == 4)
        sys        = 2;
        ml         = lambda(outputfreq);
        outputfreq = outputfreq - nfreq;
        nc         = ncode2;
        np         = nphase2;
    else
        sys        = 1;
        ml         = lambda(outputfreq);
        nc         = ncode;
        np         = nphase;
    end

    startepoch = str2double(get(findobj(gcf,'Tag','startepoch'),'String'));
    if startepoch > no_epochs
        h2 = msgbox ('Error: epoch specified in output frame > number of epochs');
        watchoff(figNumber);
        return
    end
    % definition of canonical vectors that select observation for
    % which to compute output c_{i} or c_{i+f} in eq.(3.11),(3.12)
    if strcmp(output,'slip')
        if np == 0
            h2 = msgbox ('Error: carrier slip not possible with code-only');
            watchoff(figNumber);
            return
        end
        Ifreq    = eye(np);
        c1       = Ifreq(:,outputfreq);
        c2       = zeros(nc,1);
        c        = [c2 ; c1];
        slip     = no_epochs - startepoch + 1;
        cl       = [zeros(no_epochs-slip,1);ones(slip,1)];  % vector sl in eq.(3.12)
    elseif strcmp(output,'outlier')
        if nc == 0
            h2 = msgbox ('Error: code outlier not possible with phase-only');
            watchoff(figNumber);
            return
        end
        Ifreq    = eye(nc);
        c1       = Ifreq(:,outputfreq);
        c2       = zeros(np,1);
        if strcmp(rectype,'cc')&&(outputfreq==1), c1 = [1;1]; end
        c        =  [c1 ; c2];
        slip     = 1;
        cl       = eye(no_epochs);
        cl       = cl(:,slip);  
    end

end

% definition of vector with epochs and size of slip window
if length(prad)==1
    time = 1:size(tsat,2);
else
    time = 1;
end

% initialize output matrix
OUT = zeros(length(prad),length(lrad));
OUT2= OUT;
% re-order vectors with latitudes and longitudes and compute corresponding XYZ-coordinates
lp = length(prad);
ll = length(lrad);
i=1;
j = ll;
for p = 1: lp
    plh(i:j,:) = [prad(p)*ones(ll,1) lrad' height*ones(ll,1)];  % vector [latitudes longitudes heights]
    i = j+1;
    j = i + ll -1;
end
xyz      = plh2xyzwgs(plh);

% set parameters for and create waitbar
wend     = size(tsat,2) + lp * ll - 1;
widx     = 0;
watchoff(figNumber);
wh       = waitbar(0,'Busy...');
widx2    = 1;
% loop over latitudes and longitudes, when only one location:
% loop over all epochs
for t = time             % loop over all start times for option 'one point'

    for p = 1:lp*ll       % loop over all grid points for option 'area' or 'world'

        widx = widx + 1;

        el       = cpelev1(xs(:,t),ys(:,t),zs(:,t),plh(p,:),xyz(p,:));
        idx1     = find (el > cutoff);
        no_sv    = length(idx1);
        el       = el(idx1);
        refsat   = find(el==max(el));
        xsv      = xs(idx1,t);
        ysv      = ys(idx1,t);
        zsv      = zs(idx1,t);
        zen      = deg2rad(90 - el)';
        no_sv2= 0;
        if SYSTEM == 4
            el2       = cpelev1(xs2(:,t),ys2(:,t),zs2(:,t),plh(p,:),xyz(p,:));
            idx2      = find (el2 > cutoff);
            el2       = el2(idx2);
            no_sv2    = length(idx2);
            refsat2   = find(el2==max(el2));
            xsv2      = xs2(idx2,t);
            ysv2      = ys2(idx2,t);
            zsv2      = zs2(idx2,t);
            zen2      = deg2rad(90 - el2)';
        end
        if out   == 9

            % --------------------
            % number of satellites
            % --------------------

            OUT(p) = no_sv + no_sv2;
            titel    = ' number of satellites ';

            % ____________________

        else
            OUT2(p) = no_sv + no_sv2;
            LOS = zeros(no_sv,3);
            if ~strcmp(scenario,'gf')
                LOS = cplos([xsv ysv zsv],xyz(p,:));   % LOS-vectors
            end
            LOS2 = [];
            if SYSTEM == 4
                LOS2 = cplos([xsv2 ysv2 zsv2],xyz(p,:));
                m2  = no_sv2 - 1;
                if strcmp(model,'single')
                    m2 = no_sv2;
                end
            end
            m  = no_sv - 1;
            if strcmp(model,'single')
                m = no_sv;
            end

            if out < 7  % for internal or external reliability or success rates

                if no_sv<4, OUT(p)=NaN; continue; end

                Im   = eye(m);
                em   = ones(m,1);
                if strcmp(model,'single')
                    D = Im;
                else
                    D = [Im(:,1:refsat-1) -em Im(:,refsat:m)];
                end
                if strcmp(tropo,'Tfixed')
                    Mz = [];
                else
                    Mz   = tropomap(zen,plh(3),mapfun);    % tropo. wet delay mapping function
                end
                if ~strcmp(scenario,'gf')

                    if ~strcmp(model,'single')
                        M1 = D * [LOS Mz];           % [\bar{G} \Psi] as in eq.(2.22)
                    else
                        M1 = [LOS em Mz];            % [G em \Psi] as in eq.(2.5)
                    end
                else
                    M1 = Im;
                end
                M = []; % see eq.(19)
                if strcmp(rectype,'non-cc')
                    n1 = nphase + ncode;
                    n2 = 1;
                    Qy = [];
                    while (n1 > 0 )
                        M = [M ; M1];
                        n1 = n1 - 1;
                        %elevation dependent weighting as in eq.(2.13) and (2.15)
                        Wi = sd * (sig(n2) * diag(1 + a(n2) * exp(-el/e0(n2)))).^2;
                        Qy = blkdiag(Qy,D*Wi*D'); %first term after =-sign in eq.(2.20)/(2.21)
                        n2 = n2 + 1;
                    end
                else
                    M = [M1;zeros(size(M1))];
                    W1 = sd * (sig(1) * diag(1 + a(1) * exp(-el/e0(1)))).^2;
                    W2 = sd * (sig(2) * diag(1 + a(2) * exp(-el/e0(2)))).^2;
                    Qy = blkdiag(D*W1*D',D*W2*D');
                    if nphase > 0
                        M = [M;M1;M1];
                        W1 = sd * (sig(3) * diag(1 + a(3) * exp(-el/e0(3)))).^2;
                        W2 = sd * (sig(4) * diag(1 + a(4) * exp(-el/e0(4)))).^2;
                        Qy = blkdiag(Qy,D*W1*D',D*W2*D');
                    end
                end
                if sdion > 0
                    QI   = kron(CI,D*D');       % 2nd term after =-sign in eq.(2.20)/(2.21)
                else
                    QI   = zeros(size(Qy));
                end
                iQy = inv(Qy+QI);               % inverse of matrix in eq.(2.20)/(2.21)
                N =[];  % see eq.(19)
                if nphase>0
                    N1 = zeros(ncode*m);
                    N2 = kron(LAMBDA(1:nphase,1:nphase),Im);
                    N = [N1;N2];
                end
                NQN  = N' * iQy * N;
                MQN  = M' * iQy * N;
                MQM  = M' * iQy * M;
                iMQM = inv(MQM);
                
                % also for GALILEO if system=GPS + GALILEO
                if SYSTEM == 4
                    Im2   = eye(m2);
                    em2   = ones(m2,1);
                    if strcmp(model,'single')
                        D2 = Im2;
                    else
                        D2 = [Im2(:,1:refsat2-1) -em2 Im2(:,refsat2:m2)];
                    end
                    if strcmp(tropo,'Tfixed')
                        Mz2 = [];
                    else
                        Mz2   = tropomap(zen2,plh(3),mapfun);    % tropo. wet delay mapping function
                    end
                    if ~strcmp(scenario,'gf')
                        if ~strcmp(model,'single')
                            M21 = D2 * [LOS2 Mz2];       % [\bar{G} \Psi] as in eq.(2.22)
                        else
                            M21 = [LOS2 em2 Mz2];        % [G em \Psi] as in eq.(2.5)
                        end
                    else
                        M21 = Im2;
                    end
                    M2 = [];   % see eq.(2.22)
                    n1 = nphase2 + ncode2;
                    Qy2 = [];
                    while (n1 > 0 )
                        M2 = [M2 ; M21];
                        n1 = n1 - 1;
                        %elevation dependent weighting as in eq.(2.13)and (2.15)
                        Wi = sd * (sig(n2) * diag(1 + a(n2) * exp(-el2/e0(n2)))).^2;

                        Qy2 = blkdiag(Qy2,D2*Wi*D2');
                        n2 = n2 + 1;
                    end
                    if sdion > 0
                        QI2   = kron(CI2,D2*D2'); % 2nd term after =-sign in eq.(2.20)/(2.21)
                    else
                        QI2   = zeros(size(Qy2));
                    end
                    iQy2 = inv(Qy2+QI2);
                    N2 =[];   % see eq.(2.22)
                    if nphase2>0
                        N21 = zeros(ncode2*m2);
                        N22 = kron(LAMBDA(nphase+1:end,nphase+1:end),Im2);
                        N2 = [N21;N22];
                    end

                    NQN2 = N2' * iQy2 * N2;
                    MQN2 = M2' * iQy2 * N2;
                    MQM2 = M2' * iQy2 * M2;
                    iMQM2= inv(MQM2);
                    if strcmp(scenario,'gf')
                        NQN    = blkdiag(NQN,NQN2);
                        MQN    = blkdiag(MQN,MQN2);
                        iMQM   = blkdiag(iMQM,iMQM2);
                    end
                end

                if SYSTEM == 4 && ~strcmp(scenario,'gf')
                    iMQM  = inv(MQM + MQM2);
                    ins1  = NQN - MQN' * iMQM * MQN;
                    ins12 = -MQN' * iMQM * MQN2;
                    ins2  = NQN2 - MQN2' * iMQM * MQN2;
                    ins = [ins1 ins12;ins12' ins2];
                else                    
                    ins   = NQN - MQN' * iMQM * MQN;
                end
                Q1    = inv(ins);
                Qahat = (1/no_epochs)*Q1 ;% eq.(3.23)
                % ---------------------
                % --- success-rates ---
                % ---------------------
                if ismember(out,[4 5 6])
                    if strcmp(model,'single')
                        h2 = msgbox ('Error: you should choose the baseline scenario');
                        watchoff(figNumber);
                        return
                    end
                    [Qz,Z,LQz,DQz]= decorrel(Qahat);
                    iL = inv(LQz');
                    OUT(p) = prod ( 2 * normcdf(1./(2*sqrt(DQz))) -1 );   % eq.(3.31)

                    if out == 6
                        SR(p) = OUT(p);
                    end

                    titel   = ' success rates ';
                end