function HavocPlot

global STATION OBS JOB X Z VFACT Qxx               

fixed    = JOB(1).fixed;       % number of fixed stations
floating = JOB(1).floating;    % number of floating stations
dirns    = JOB(1).dirns;       % number of directions
dists    = JOB(1).dists;       % number of distances
cbrgs    = JOB(1).cbrgs;       % number of constrained bearings
cdists   = JOB(1).cdists;      % number of constrained distances
cangles  = JOB(1).cangles;     % number of constrained angles
orients  = JOB(1).orients;     % number of orientation constants
n        = JOB(1).n;           % number of observation equations
u        = JOB(1).u;           % number of unknowns
c        = JOB(1).c;           % number of constraint equations

d2r   = 180/pi;
pion2 = pi/2;
%-----------------
% plot the network
%-----------------
for k = 1:fixed+floating
  E(k) = STATION(k).east;
  N(k) = STATION(k).north;
end  
figure(1);
clf(1);
hold on;
%box on;
box off;
plot(E,N,'k^');
%grid on;
grid off;
axis equal;
%axis on;
axis off;
E = [];
N = [];
for m = 1:(dirns+dists)
    for j = 1:fixed+floating
        if (OBS(m).at == STATION(j).serial)
            E(1) = STATION(j).east;  
            N(1) = STATION(j).north;  
        end  
        if (OBS(m).to == STATION(j).serial)
            E(2) = STATION(j).east;  
            N(2) = STATION(j).north;  
        end  
    end
    plot(E,N,'k-');
end
title('Plot of Network');
xlabel('East (m)');
ylabel('North (m)');

%------------------------------------------
% plot error ellipses for floating stations
%------------------------------------------
k = 0;
for j = 1:fixed+floating
    if (STATION(j).code == 2)
        k = k+1;
        flag   = 0;  % flag is false; no error ellipse for this point
        if (c == 0);
            flag = 1;  % flag is true; error ellipse for this point  
        else
            serial = STATION(j).serial;
            for m = 1:(dirns+dists+cbrgs+cdists+cangles)
                if ((OBS(m).code == 3) && (OBS(m).status == 1))
                    if (OBS(m).at == serial || OBS(m).to == serial) 
                        flag = 1;
                    end  
                end  
                if ((OBS(m).code == 4) && (OBS(m).status == 1))
                    if (OBS(m).at == serial || OBS(m).to == serial) 
                        flag = 1;
                    end  
                end  
                if ((OBS(m).code == 5) && (OBS(m).status == 1))
                    if (OBS(m).at == serial || OBS(m).to == serial) 
                        flag = 1;
                    end  
                end  
                if ((OBS(m).code == 6) && (OBS(m).status == 1))
                    if (OBS(m).at == serial || OBS(m).to == serial) 
                        flag = 1;
                    end  
                end  
                if (OBS(m).code == 7)
                    if (OBS(m).status < 4 || OBS(m).status == 5 || OBS(m).status == 7) 
                        if (OBS(m).at == serial || OBS(m).left == serial || OBS(m).right == serial)
                            flag = 1;
                        end  
                    end  
                end  
            end  
        end
    
        if flag == 1
            Ep = STATION(j).east;  
            Np = STATION(j).north;  
	  	    %----------------------------------------------
			% Calculate the parameters of the error ellipse
			%----------------------------------------------
			% a = semi-major axis
			% b = semi-minor axis
			% theta = rotation
            i   = (2*k)-1;
			s2e = Qxx(i,i);
			sen = Qxx(i,i+1);
			s2n = Qxx(i+1,i+1);
			w   = sqrt((s2e-s2n)^2 + (2*sen)^2);
            a   = sqrt(0.5*(s2e + s2n + w));
			b   = sqrt(0.5*(s2e + s2n - w));
			two_beta = atan2((2*sen),(s2e-s2n));
			
            % compute second derivative
            f_dd = 2*(s2n-s2e)*cos(two_beta) - 4*sen*sin(two_beta);
            % test the second derivative to determine max or min
            beta = two_beta/2*d2r;
            if f_dd < 0
                theta = beta;
            else
                theta = beta-90
            end    
            
            %axis(k,1) = a;
			%axis(k,2) = b;
			%axis(k,3) = theta;
			
			%-------------------------------------------------------------------
			% Call function 'ellipse' with parameters a,b,theta and receive back
			% X,Y coordinates whose origin is at the centre of the ellipse
			%-------------------------------------------------------------------
			a = a*120;  % 150 is a scale factor
			b = b*120;
			[X,Y] = ellipse(a,b,theta);  
			X = X + Ep;
			Y = Y + Np;
			
			%-------------------------------------------------------------------
			% Set the X,Y coordinates of the major and minor axes of the ellipse
			%-------------------------------------------------------------------
			aX(1) = X(180);
			aY(1) = Y(180);
			aX(2) = X(360);
			aY(2) = Y(360);
			bX(1) = X(90);
			bY(1) = Y(90);
			bX(2) = X(270);
			bY(2) = Y(270);
			
			%----------------------------------------
			% plot the error ellipse for the location
			%----------------------------------------
			plot(X,Y,'r-',aX,aY,'b-',bX,bY,'b-');
        end  
    end  % end if for floating stations  
end