visualcompute.m

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

% visualcompute
% perform the actual computations given the GUI settings
% called by VISUAL.M
%
% See also VISUAL

figNumber = watchon;

lc  =  get(findobj(gcf,'Tag','locchoice'),'Value');
if lc == 2 && get(findobj(gcf,'Tag','point'),'Value')

    % select area by pointing in map
    watchoff(figNumber);
    uiwait(msgbox ('Click and drag your mouse to select the desired area','','modal'));
    figure(figNumber); zoom off;
    waitforbuttonpress;
    point1 = get(gca,'CurrentPoint');    % button down detected
    finalRect = rbbox;                   % return figure units
    point2 = get(gca,'CurrentPoint');    % button up detected
    x      = [point1(1,1);point2(1,1)];
    y      = [point1(1,2);point2(1,2)];
    figNumber = watchon;

    xmin = round(min(x));
    xmax = round(max(x));
    ymin = round(min(y));
    ymax = round(max(y));

    reso = str2double(get (findobj(gcf,'Tag','resol'),'String'));
    grint = reso;

    pdeg = ymax:-grint:ymin;
    ldeg = xmin:grint:xmax;
    prad = pdeg*pi/180;
    lrad = ldeg*pi/180;

    set (findobj(gcf,'Tag','long1'),'String',num2str(xmin));
    set (findobj(gcf,'Tag','long2'),'String',num2str(xmax));
    set (findobj(gcf,'Tag','lat1'),'String',num2str(ymin));
    set (findobj(gcf,'Tag','lat2'),'String',num2str(ymax));

elseif lc == 1 && get(findobj(gcf,'Tag','point'),'Value')

    % select one point by pointing in map
    watchoff(figNumber);
    uiwait(msgbox ('Click to select the desired location','','modal'));
    figure(figNumber); zoom off;
    [xmin,ymin]= ginput(1); figNumber = watchon;
    prad = ymin * pi /180;
    lrad = xmin * pi /180;

    set (findobj(gcf,'Tag','long1'),'String',num2str(xmin));
    set (findobj(gcf,'Tag','long2'),'String','');
    set (findobj(gcf,'Tag','lat1'),'String',num2str(ymin));
    set (findobj(gcf,'Tag','lat2'),'String','');

elseif get(findobj(gcf,'Tag','type'),'Value') && lc ==1

    % give desired latitude(s) and longitude(s) by typing

    ymin = str2double(get (findobj(gcf,'Tag','lat1'),'String'));
    ymax = str2double(get (findobj(gcf,'Tag','lat2'),'String'));
    xmin = str2double(get (findobj(gcf,'Tag','long1'),'String'));
    xmax = str2double(get (findobj(gcf,'Tag','long2'),'String'));
    prad = ymin*pi/180;
    lrad = xmin*pi/180;


elseif get(findobj(gcf,'Tag','type'),'Value')&& lc == 2

    % area defined typing

    y=[str2double(get (findobj(gcf,'Tag','lat1'),'String'));...
        str2double(get (findobj(gcf,'Tag','lat2'),'String'))];
    x=[str2double(get (findobj(gcf,'Tag','long1'),'String'));...
        str2double(get (findobj(gcf,'Tag','long2'),'String'))];
    ymin = min(y);
    ymax = max(y);
    xmin = min(x);
    xmax = max(x);

    reso = str2double(get (findobj(gcf,'Tag','resol'),'String'));
    grint = reso;
    pdeg  = ymax:-grint:ymin;
    ldeg  = xmin:grint:xmax;
    prad  = pdeg*pi/180;
    lrad  = ldeg*pi/180;


elseif lc == 3

    % computations for entire world with defined resolution

    reso = str2double(get (findobj(gcf,'Tag','resol'),'String'));
    xmin = -180; xmax = 180;
    ymin = -90; ymax = 90;
    grint = reso * pi/180;
    prad = pi/2 : -grint : -pi/2;
    lrad = -pi : grint : pi;

end

height = str2double(get(findobj(gcf,'Tag','height'),'String'));

% find out which output parameter must be computed;
% if satellite tracks are chosen --> plot
figure(figNumber);
out = get(findobj(gcf,'Tag','output'),'Value');
if out == 10
    visual('initialFigure','begin')
    watchoff(figNumber);
    return
end

% frequencies:
% frequencies depend on system choice
SYSTEM = get(findobj(gcf,'Tag','system'),'Value');
rectype = 'non-cc';


if SYSTEM == 2 % GLONASS

    frq = [1602e6;1246e6];

else     % GPS / Galileo

    frq = [1575.42e6;1227.60e6;1176.45e6;1227.60e6;...
        1575.42e6; 1278.75e6;1207.14e6;1176.45e6];

    % receiver type: cross-correlating when cross-correlated code
    % observation is chosen
    if (get(findobj(gcf,'Tag','pL2L1'),'Value'))
        rectype = 'cc';
    end

end

% measurement scenario:
% - single baseline, stationary or roving receiver (sr of rr) (user defined)
% - geometry free
% - single point

if (get(findobj(gcf,'Tag','receiver'),'Value'))==1
    scenario = 'sr';
else
    scenario = 'rr';
end
if (get(findobj(gcf,'Tag','scenario'),'Value'))==3
    scenario = 'gf';
    model = 'baseline';
    sd    = 2;
elseif (get(findobj(gcf,'Tag','scenario'),'Value'))==1
    model = 'single';
    sd    = 1;
else
    model = 'baseline';
    sd    = 2;
end

% standard deviations of observations (user defined)
freqsel  = [get(findobj(gcf,'Tag','pL1'),'Value');...
    get(findobj(gcf,'Tag','pL2'),'Value');...
    get(findobj(gcf,'Tag','pL5'),'Value');...
    get(findobj(gcf,'Tag','pL2L1'),'Value');...
    get(findobj(gcf,'Tag','pE1'),'Value');...
    get(findobj(gcf,'Tag','pE6'),'Value')
    get(findobj(gcf,'Tag','pE5b'),'Value');...
    get(findobj(gcf,'Tag','pE5a'),'Value')];
sigcode  = [str2double(get(findobj(gcf,'Tag','sigc1'),'String'));...
    str2double(get(findobj(gcf,'Tag','sigc2'),'String'));...
    str2double(get(findobj(gcf,'Tag','sigc3'),'String'));...
    str2double(get(findobj(gcf,'Tag','sigc4'),'String'));...
    str2double(get(findobj(gcf,'Tag','sigc1b'),'String'));...
    str2double(get(findobj(gcf,'Tag','sigc4b'),'String'));...
    str2double(get(findobj(gcf,'Tag','sigc3b'),'String'));...
    str2double(get(findobj(gcf,'Tag','sigc2b'),'String'))];
sigphase = [str2double(get(findobj(gcf,'Tag','sigp1'),'String'));...
    str2double(get(findobj(gcf,'Tag','sigp2'),'String'));...
    str2double(get(findobj(gcf,'Tag','sigp3'),'String'));...
    str2double(get(findobj(gcf,'Tag','sigp4'),'String'));...
    str2double(get(findobj(gcf,'Tag','sigp1b'),'String'));...
    str2double(get(findobj(gcf,'Tag','sigp4b'),'String'));...
    str2double(get(findobj(gcf,'Tag','sigp3b'),'String'));...
    str2double(get(findobj(gcf,'Tag','sigp2b'),'String'))];
acode  = [str2double(get(findobj(gcf,'Tag','ac1'),'String'));...
    str2double(get(findobj(gcf,'Tag','ac2'),'String'));...
    str2double(get(findobj(gcf,'Tag','ac3'),'String'));...
    str2double(get(findobj(gcf,'Tag','ac4'),'String'));...
    str2double(get(findobj(gcf,'Tag','ac1b'),'String'));...
    str2double(get(findobj(gcf,'Tag','ac4b'),'String'));...
    str2double(get(findobj(gcf,'Tag','ac3b'),'String'));...
    str2double(get(findobj(gcf,'Tag','ac2b'),'String'))];
aphase = [str2double(get(findobj(gcf,'Tag','ap1'),'String'));...
    str2double(get(findobj(gcf,'Tag','ap2'),'String'));...
    str2double(get(findobj(gcf,'Tag','ap3'),'String'));...
    str2double(get(findobj(gcf,'Tag','ap4'),'String'));...
    str2double(get(findobj(gcf,'Tag','ap1b'),'String'));...
    str2double(get(findobj(gcf,'Tag','ap4b'),'String'));...
    str2double(get(findobj(gcf,'Tag','ap3b'),'String'));...
    str2double(get(findobj(gcf,'Tag','ap2b'),'String'))];
e0code  = [str2double(get(findobj(gcf,'Tag','ec1'),'String'));...
    str2double(get(findobj(gcf,'Tag','ec2'),'String'));...
    str2double(get(findobj(gcf,'Tag','ec3'),'String'));...
    str2double(get(findobj(gcf,'Tag','ec4'),'String'));...
    str2double(get(findobj(gcf,'Tag','ec1b'),'String'));...
    str2double(get(findobj(gcf,'Tag','ec4b'),'String'));...
    str2double(get(findobj(gcf,'Tag','ec3b'),'String'));...
    str2double(get(findobj(gcf,'Tag','ec2b'),'String'))];
e0phase = [str2double(get(findobj(gcf,'Tag','ep1'),'String'));...
    str2double(get(findobj(gcf,'Tag','ep2'),'String'));...
    str2double(get(findobj(gcf,'Tag','ep3'),'String'));...
    str2double(get(findobj(gcf,'Tag','ep4'),'String'));...
    str2double(get(findobj(gcf,'Tag','ep1b'),'String'));...
    str2double(get(findobj(gcf,'Tag','ep4b'),'String'));...
    str2double(get(findobj(gcf,'Tag','ep3b'),'String'));...
    str2double(get(findobj(gcf,'Tag','ep2b'),'String'))];

ifreq    = find(freqsel==1);

freq  = frq(ifreq);
if SYSTEM == 4
    nfreq = length(find(ifreq<5)); % nfreq = no. of GPS frequencies
else
    nfreq = length(ifreq);
end

if get(findobj(gcf,'Tag','codesel'),'Value')
    ncode  = nfreq;
    ncode2 = length(ifreq)-ncode;
else
    ncode  = 0;
    ncode2 = 0;
end
if get(findobj(gcf,'Tag','phasesel'),'Value')
    nphase  = nfreq;
    nphase2 = length(ifreq)-nphase;
else
    nphase = 0;
    nphase2= 0;
end

% wavelengths
lambda = 299792458 ./ freq;
if strcmp(rectype,'cc')
    lambda(2) = 299792458 / (freq(1)-freq(2));
end
LAMBDA = diag(lambda);

sigcode = sigcode(ifreq);
sigphase= sigphase(ifreq);
acode   = acode(ifreq);
aphase  = aphase(ifreq);
e0code  = e0code(ifreq);
e0phase = e0phase(ifreq);

sig  = [sigcode(1:nfreq);sigphase(1:nfreq);sigcode(nfreq+1:end);sigphase(nfreq+1:end)];
a    = [acode(1:nfreq);aphase(1:nfreq);acode(nfreq+1:end);aphase(nfreq+1:end)];
e0   = [e0code(1:nfreq);e0phase(1:nfreq);e0code(nfreq+1:end);e0phase(nfreq+1:end)];

% get ionospheric standard deviation
if out<7
    ic = get(findobj(gcf,'Tag','ioncase'),'Value');
    if ic == 3 % ionosphere weighted

        if get(findobj(gcf,'Tag','sigion'),'Value')
            sdion = str2double(get(findobj(gcf,'Tag','sdionweight'),'String'));
        elseif get(findobj(gcf,'Tag','baseline'),'Value')
            sdion = 0.68 *(frq(1)/freq(1))* ...
                str2double(get(findobj(gcf,'Tag','blionweight'),'String')) / 1000;
        end

    elseif ic == 1 % ionosphere fixed

        sdion = 0;

    else % ionosphere float

        sdion = 99999;

    end

    % single difference ionospheric vc-matrix

    mu    = ([lambda(1:nfreq)]./ lambda(1)).^2;
    mu    = [mu(1:ncode);-mu(1:nphase)];
    s2    = sd*sdion^2;
    
    CI    = s2 * mu * mu';

    if SYSTEM ==4
        mu2   = (lambda(nfreq+1:end)./ lambda(nfreq+1)).^2  ;
        mu2   = [mu2(1:ncode2);-mu2(1:nphase2)];        
        CI2   = s2 * mu2 * mu2';        
    end


    % troposphere model: float or fixed
    % if tropo = 'float' : mapping function = mapfun

    if (get(findobj(gcf,'Tag','tropocase'),'Value'))==1
        tropo = 'Tfixed';
    else
        tropo = 'Tfloat';
        if get(findobj(gcf,'Tag','mapfun'),'Value')==1
            mapfun = 's';  % 1/cos(z)
        else
            mapfun = 'I';  % Ifadis
        end
    end

end
% get ephemerides and compute satellite positions

almanac       = get(findobj(gcf,'Tag','almanac'),'String');
starttime  = [get(findobj(gcf,'Tag','startdate'),'String') ' ' ...
   get(findobj(gcf,'Tag','starttime'),'String')] ;
endtime    = [get(findobj(gcf,'Tag','enddate'),'String') ' ' ...
   get(findobj(gcf,'Tag','endtime'),'String')] ;
no_epochs     = str2double(get(findobj(gcf,'Tag','nepochs'),'String'));

if length(prad) ==1
    interval   = str2double(get(findobj(gcf,'Tag','intv'),'String'));
    tsat       = mktsat ( starttime,endtime,interval);
else
    tsat       = mktsat ( starttime,starttime,1);
    tsat       = tsat(1);
end
if SYSTEM ==4
    eph        = rdyuma(almanac);
    almanac2   = get(findobj(gcf,'Tag','almanac2'),'String');
    eph2       = rdyuma(almanac2);
    i          = find(eph(:,27)==0); % find healthy satellites
    eph        = eph(i,:);
    i          = find(eph2(:,27)==0); % find healthy satellites GALILEO
    eph2       = eph2(i,:);
    [xs,ys,zs] = cpaziele (tsat,eph);
    [xs2,ys2,zs2] = cpaziele (tsat,eph2);
else
    if strcmp(almanac(1:4),'yuma')
        eph = rdyuma(almanac);
    else
        ephtmp =rxnav(almanac);
        ephtmp = sortrows(ephtmp);
        eph(1,:)=ephtmp(1,:);