function results
% Function RESULTS is part of program HAVOC, a (free net) 
% least squares traverse adjustment program.  
% Function RESULTS writes data from the adjustment to a 
% nominated ASCII file.

%============================================================================
% Function:  results
%
% Author:
%  Rod Deakin, 
%  School of Mathematical & Geospatial Science, 
%  RMIT University,
%  GPO Box 2476V, MELBOURNE VIC 3001
%  AUSTRALIA
%  email: rod.deakin@rmit.edu.au
%
% Date:
%  Version 1.0   3 August   2006
%
% Functions required:
%  deg2dms
%  join
%
% Remarks:
%  This function prints results of the adjustment to an ASCII text file.
%
% Variables:
%  aik,bik   direction coefficients of the line AT(i) to TO(k)
%  beta      angle (+tve anticlockwise) from x-axis to major axis of
%            ellipse
%  brg       bearing of major axis of error ellipse
%  cik,dik   distance  coefficients of the line AT(i) to TO(k)
%  dz        correction to approximate orientation constant
%  d2r       degree to radian conversion factor = 57.29577951...
%  east      east coordinate
%  ei,ni     east and north coord of station AT(i)
%  ek,nk     east and north coord of station TO(k)
%  flag3,flag4,flag5,  flags for printing headings
%  flag6,flag7,
%  flag9     flag for printing orientation constant
%  i_list    location in list of fixed + floating stations of station AT
%  k_list    location in list of fixed + floating stations of station TO
%  i_serial  serial number of stations AT(i)
%  k_serial  serial number of stations TO(k)
%  index     index for orientation const/coeff in vector Z/Line
%  j,m       loop counters
%  last      serial number of last station AT
%  north     north coordinate
%  phi       bearing of line (degrees)
%  pion2     pi/2
%  serial    serial number of station
%  s         distance
%  se,sn     st.dev's east and north
%  stdev     standard cdeviation
%  theta     "observed" bearing = observed direction + orientation
%  two-beta  2(beta)
%  v         residual
%  VFACT     variance factor
%  X         u,1 vector of parameters in adjustment.  X(1,1) to
%            X(2*floating,1) are corrections to approx. coords.
%            X(2*floating+1,1) to X(u,1) are corrections to
%            approximate orientation constants.
%  z         orientation constant
%  Z         vector of approximate orientation constants
%
%  fixed     no. of fixed stations
%  floating  no. of floating stations
%  dirns     no. of observed directions
%  dists     no. of observed distances
%  cbrgs     no. of constrained bearings
%  cdists    no. of constrained distances
%  cangles   no. of constrained angles
%  orients   no. of orientation constants
%  n         total number of observations (n = dirns + dists)
%  u         no. of unknowns (u = 2*floating + orients)
%  c         no. of constraints
%============================================================================

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;

infilepath  = JOB(1).infile;
outfilepath = JOB(1).outfile;
fidout  = fopen(outfilepath,'wt');

fprintf(fidout,'Program HAVOC: Horizontal Adjustment by Variation Of Coordinates\n');
fprintf(fidout,'\nThis file is   %s\ngenerated from %s',outfilepath,infilepath);
%-------------------------------
% write out adjusted coordinates
%-------------------------------
if fixed == 0
    fprintf(fidout,'\n\nFREE NET ADJUSTMENT - NO FIXED STATIONS');    
else    
    fprintf(fidout,'\n\nFixed Stations:');
    fprintf(fidout,'\n\nSerial     East            North');
    for j = 1:fixed
        serial = STATION(j).serial;
        east   = STATION(j).east;
        north  = STATION(j).north;
        fprintf(fidout,'\n%6d   %13.6f   %13.6f',serial,east,north);
    end  
end
fprintf(fidout,'\n\nAdjusted Stations:');
fprintf(fidout,'\n\nSerial     East            North');
for j = fixed+1:fixed+floating
  serial = STATION(j).serial;
  east   = STATION(j).east;
  north  = STATION(j).north;
  fprintf(fidout,'\n%6d   %13.6f   %14.6f',serial,east,north);
end  
%-----------------------------------------------------------
% write out observed directions and distances with residuals
%-----------------------------------------------------------
flag3 = 0;
flag4 = 0;
index = 0;
last  = 0;
for m = 1:dirns+dists
  %
  % find the stations AT and TO in the list of stations and
  % the east and north coordinates
  for j = 1:fixed+floating
    if (OBS(m).at == STATION(j).serial)
      i_serial = STATION(j).serial;  
      ei = STATION(j).east;  
      ni = STATION(j).north;  
    end  
    if (OBS(m).to == STATION(j).serial)
      k_serial = STATION(j).serial;  
      ek = STATION(j).east;  
      nk = STATION(j).north;  
    end  
  end
  %
  % calculate the bearing and distance
  [phi,s,aik,bik,cik,dik]=join(ei,ni,ek,nk);
  %--------------------------
  % print observed directions
  %--------------------------
  if (OBS(m).code == 3)
    flag9 = 0;
    if (OBS(m).at ~= last)
      index = index+1;
      flag9 = 1;
      dz    = X(2*floating+index,1)/3600;
      z     = Z(index,1) + dz;
      last  = OBS(m).at;
      fprintf(fidout,'\n');
    end  
    theta = OBS(m).val + z;
    if (theta > 360.0)
      theta = theta - 360.0;
    end  
    if (theta > 359.999997222)
      theta = 0.0;  
    end  
    v = (phi - theta)*3600.0;
    stdev = OBS(m).sigma;
    
    if (flag3 == 0)
      fprintf(fidout,'\n\nObserved Directions:\n');
      fprintf(fidout,'\n     station    observation    st.dev     resid');
      fprintf(fidout,'  orient const  plane brg');
      fprintf(fidout,'\n    at      to  deg min sec     sec        sec');
      fprintf(fidout,'   deg min sec   deg min sec     plane dist\n');
      flag3 = 1;
    end  
    % print stations AT(i) and TO(k)
    fprintf(fidout,'\n%6d  %6d',i_serial,k_serial);
    % print observation and residual
    [deg,min,sec]=deg2dms(OBS(m).val);
    fprintf(fidout,'  %3d %2d %5.2f %7.2f  %9.2f',deg,min,sec,stdev,v);
    if (flag9 == 1)
      % print orientation constant
      [deg,min,sec]=deg2dms(z);
      fprintf(fidout,'  %3d %2d %5.2f',deg,min,sec);
      % print computed bearing and distance
      [deg,min,sec]=deg2dms(phi);
      fprintf(fidout,'  %3d %2d %5.2f %12.4f',deg,min,sec,s);
    else
      % print computed bearing and distance
      [deg,min,sec]=deg2dms(phi);
      fprintf(fidout,'%19d %2d %5.2f %12.4f',deg,min,sec,s);
    end  
  end  % end loop for printing directions  
  %-------------------------
  % print observed distances
  %-------------------------
  if (OBS(m).code == 4)
    if(OBS(m).at ~= last)
      last = OBS(m).at;
      fprintf(fidout,'\n');
    end  
    % calculate residual
    v = s - OBS(m).val;
    stdev = OBS(m).sigma;
    if (flag4 == 0)
      fprintf(fidout,'\n\nObserved Distances:\n');
      fprintf(fidout,'\n     station      observation   st.dev      resid');
      fprintf(fidout,'                plane brg');
      fprintf(fidout,'\n    at      to       metres     metres     metres');
      fprintf(fidout,'               deg min sec     plane dist\n');
      flag4 = 1;      
    end  
    % print stations AT(i) and TO(k)
    fprintf(fidout,'\n%6d  %6d',i_serial,k_serial);
    % print observation and residual
    fprintf(fidout,'%13.4f %10.4f %10.4f',OBS(m).val,stdev,v);
    % print computed bearing and distance
    [deg,min,sec]=deg2dms(phi);
    fprintf(fidout,'%18d %2d %5.2f %12.4f',deg,min,sec,s);
  end  % end loop for printing distances
end  % end loop for printing directions and distances

if fixed ~= 0
    %
    % test to see if there are any constraints
    flag5 = 0;
    flag6 = 0;
    flag7 = 0;
    if (c > 0)
        for j = dirns+dists:dirns+dists+c
            if (OBS(j).code == 5)
                if (flag5 == 0)
                    fprintf(fidout,'\n\n\nConstrained Bearings:\n');  
                    fprintf(fidout,'\n     station      bearing');
                    fprintf(fidout,'\n    at      to  deg min sec\n');
                    flag5 = 1;
                end  
                % print stations AT and TO
                % fprintf(fidout,'\n%6d  %6d',i,k);  old line commented out 10 Aug 2004
                fprintf(fidout,'\n%6d  %6d',OBS(j).at,OBS(j).to);  % new line added 10 Aug 2004
                % print observation
                % [deg,min,sec]=deg2dms(OBS(m).val); old line commented out 10 Aug 2004
                [deg,min,sec]=deg2dms(OBS(j).val);  % new line added 10 Aug 2004
                fprintf(fidout,'  %3d %2d %5.2f',deg,min,sec);
            end  
            if (OBS(j).code == 6)
                if (flag6 == 0)
                    fprintf(fidout,'\n\n\nConstrained Distances:\n');  
                    fprintf(fidout,'\n     station        distance');
                    fprintf(fidout,'\n    at      to       metres\n');
                    flag6 = 1;
                end  
                % print stations AT and TO
                % fprintf(fidout,'\n%6d  %6d',i,k); old line commented out 10 Aug 2004
                fprintf(fidout,'\n%6d  %6d',OBS(j).at,OBS(j).to);  % new line added 10 Aug 2004
                % print observation
                fprintf(fidout,'  %10.4f',OBS(j).val);
            end  
            if (OBS(j).code == 7)
                if (flag7 == 0)
                    fprintf(fidout,'\n\n\nConstrained Angles:\n');  
                    fprintf(fidout,'\n        stations           angle');
                    fprintf(fidout,'\n    at    left   right  deg min sec\n');
                    flag5 = 1;
                end  
                % print stations AT and LEFT and RIGHT
                fprintf(fidout,'\n%6d  %6d  %6d',OBS(j).at,OBS(j).left,OBS(j).right);
                % print observation
                [deg,min,sec]=deg2dms(OBS(j).val);
                fprintf(fidout,'  %3d %2d %5.2f',deg,min,sec);
            end  
        end  
    end  % end loop for printing constraints
end
%
%------------------------------------
% print out error ellipse information
%------------------------------------

%                                                    Standard Error Ellipse
%                                     St. Devs       semi-  semi-    bearing
%Serial      East          North      E       N      major  minor    d  m  s
%  2076   325567.889   5880123.456   0.015   0.022   0.022  0.024  156 33 24

fprintf(fidout,'\n\nStandard Deviations and Standard Error Ellipse:');
fprintf(fidout,'\n%62sStandard Error Ellipse',' ');
fprintf(fidout,'\n%40sSt. Devs          semi-      semi-    bearing',' ');
fprintf(fidout,'\nSerial      East          North       E         N         major      minor    d  m  s\n');
%                      1     1030.000       997.000  1.91e-005  2.82e-005  2.83e-005  1.89e-005    7 19 42
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
	  	    %----------------------------------------------
			% Calculate the parameters of the error ellipse
			%----------------------------------------------
			% a = semi-major axis
			% b = semi-minor axis
			% beta = rotation (+tve anticlockwise from x-axis)
            Ep = STATION(j).east;  
            Np = STATION(j).north;  
            i   = (2*k)-1;
			s2e = VFACT*Qxx(i,i)/10000;
			sen = VFACT*Qxx(i,i+1)/10000;
			s2n = VFACT*Qxx(i+1,i+1)/10000;
            se  = sqrt(s2e);
            sn = sqrt(s2n);
			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
                brg = 450-beta;
            else
                brg = 450-beta-90
            end    
            if brg > 360
                brg = brg-360;
            end
            % print out the data
            serial = STATION(j).serial;
            east   = STATION(j).east;
            north  = STATION(j).north;
            fprintf(fidout,'\n%6d  %12.4f   %12.4f',serial,east,north);
            [deg,min,sec]=deg2dms(brg);
            if(a <0.005 || b <0.005)
                fprintf(fidout,'  %8.2e  %8.2e  %8.2e  %8.2e  %3d %2d %2.0f',se,sn,a,b,deg,min,sec);
            else
                fprintf(fidout,'  %7.3f    %7.3f   %7.3f    %7.3f     %3d %2d %2.0f',se,sn,a,b,deg,min,sec);
            end
        end
    end
end
%-------------------------------------------------------
% print out variances of adjusted bearings and distances
%-------------------------------------------------------

fprintf(fidout,'\n\nAdjusted Bearings and Distances:\n');
fprintf(fidout,'\n                plane                     plane');
fprintf(fidout,'\n     station    bearing        st.dev     distance  st.dev');
fprintf(fidout,'\n    at      to  deg min sec     sec          m        m\n');
for m = 1:dirns+dists
    %
    % find the stations AT and TO in the list of stations and
    % the east and north coordinates
    for j = 1:fixed+floating
        if (OBS(m).at == STATION(j).serial)
            i_serial = STATION(j).serial;  
            ei = STATION(j).east;  
            ni = STATION(j).north;  
            i_list = j;
        end  
        if (OBS(m).to == STATION(j).serial)
            k_serial = STATION(j).serial;  
            ek = STATION(j).east;  
            nk = STATION(j).north;
            k_list = j;
        end  
    end
    if (OBS(m).status == 1)  % station AT is floating, station TO is floating
        % variances and covariances of station AT(i)
        i = (2*(i_list-fixed))-1;
        s2ei  = VFACT*Qxx(i,i);
        seini = VFACT*Qxx(i,i+1);
        s2ni  = VFACT*Qxx(i+1,i+1);
        % variances and covariances of station TO(k)
        k = (2*(k_list-fixed))-1;
        s2ek  = VFACT*Qxx(k,k);
        seknk = VFACT*Qxx(k,k+1);
        s2nk  = VFACT*Qxx(k+1,k+1);
        % covariances between stations AT(i) and TO(k)
        seiek = VFACT*Qxx(i,k);
        snink = VFACT*Qxx(i+1,k+1);
        seink = VFACT*Qxx(i,k+1);
        sekni = VFACT*Qxx(k,i+1);
    end    
    if (OBS(m).status == 2)  % station AT(i) is fixed, station TO(k) is floating
        % variances and covariances of station AT(i)
        s2ei  = 0;
        seini = 0;
        s2ni  = 0;
        % variances and covariances of station TO(k)
        k = (2*(k_list-fixed))-1;
        s2ek  = VFACT*Qxx(k,k);
        seknk = VFACT*Qxx(k,k+1);
        s2nk  = VFACT*Qxx(k+1,k+1);
        % covariances between stations AT(i) and TO(k)
        seiek = 0;
        snink = 0;
        seink = 0;
        sekni = 0;
    end
    if (OBS(m).status == 3)  % station AT(i) is floating, station TO(k) is fixed
        % variances and covariances of station AT(i)
        i = (2*(i_list-fixed))-1;
        s2ei  = VFACT*Qxx(i,i);
        seini = VFACT*Qxx(i,i+1);
        s2ni  = VFACT*Qxx(i+1,i+1);
        % variances and covariances of station TO(k)
        s2ek  = 0;
        seknk = 0;
        s2nk  = 0;
        % covariances between stations AT(i) and TO(k)
        seiek = 0;
        snink = 0;
        seink = 0;
        sekni = 0;
    end
    if (OBS(m).status == 4)  % station AT(i) is fixed, station TO(k) is fixed
        % variances and covariances of station AT(i)
        s2ei  = 0;
        seini = 0;
        s2ni  = 0;
        % variances and covariances of station TO(k)
        s2ek  = 0;
        seknk = 0;
        s2nk  = 0;
        % covariances between stations AT(i) and TO(k)
        seiek = 0;
        snink = 0;
        seink = 0;
        sekni = 0;
    end    
    % calculate the adjusted bearing and distance
    [phi,s,aik,bik,cik,dik]=join(ei,ni,ek,nk);
    % calculate variances of adjusted bearings and distances
    s2brg  = bik*bik*(s2ei + s2ek - 2*seiek) + ...
             aik*aik*(s2ni + s2nk - 2*snink) + ...
             2*aik*bik*(seini + seknk - seink - sekni);
    s2dist = dik*dik*(s2ei + s2ek - 2*seiek) + ...
             cik*cik*(s2ni + s2nk - 2*snink) + ...
             2*cik*dik*(seini + seknk - seink - sekni);
    % calculate st.dev's of adjusted bearings and distances
    sbrg  = sqrt(s2brg);
    sdist = sqrt(s2dist)/100;

    % print stations AT and TO
    fprintf(fidout,'\n%6d  %6d',i_serial,k_serial);
    % print computed bearing and st.dev.
    [deg,min,sec]=deg2dms(phi);
    if (OBS(m).status == 4)
        % no standard deviations
        % print computed bearing
        [deg,min,sec]=deg2dms(phi);
        fprintf(fidout,'  %3d %2d %5.2f',deg,min,sec);
        % print computed distance
        fprintf(fidout,'%21.3f',s);
    else
        % print computed bearing and st.dev.
        [deg,min,sec]=deg2dms(phi);
        fprintf(fidout,'  %3d %2d %5.2f %7.2f',deg,min,sec,sbrg);
        % print computed distance and st.dev.
        fprintf(fidout,' %13.4f %8.4f',s,sdist);
    end
end



%--------------------------------------------------
% print out the Variance Factor and Cofactor matrix
%--------------------------------------------------
fprintf(fidout,'\n\n\nVariance Factor = %12.4e',VFACT);
%
fprintf(fidout,'\n\n\nCofactor matrix Qxx (upper triangular portion)');
fprintf(fidout,'\nUnits of metre squared, printed in same order as adjusted stations.');
fprintf(fidout,'\nThe * symbol indicates new row of cofactor matrix beginning at');
fprintf(fidout,'\ndiagonal element, east variance followed by covariances');
fprintf(fidout,'\nthen north variance followed by covariances.');
for i = 1:2*floating
  k = 1;
  fprintf(fidout,'\n* ');
  for j = i:2*floating
    if (k > 6)
      fprintf(fidout,'\n  ');
      k = 1;
    end  
    % note that the elements of Qxx in the program are cm's-squared
    fprintf(fidout,'% 12.4e ',Qxx(i,j)/10000);
    k = k+1;
  end  
end  
%------------------------------
% print out the adjustment info
%------------------------------
if fixed == 0
    fprintf(fidout,'\n\n\nFree Net Adjustment Data:\n');
    fprintf(fidout,'\n%2d fixed stations',fixed);
else
    fprintf(fidout,'\n\n\nAdjustment Data:\n');
    fprintf(fidout,'\n%2d fixed stations',fixed);
end    
fprintf(fidout,'\n%2d floating stations',floating);
fprintf(fidout,'\n%2d observed directions',dirns);
fprintf(fidout,'\n%2d observed distances',dists);
if fixed == 0
   if dists == 0
       c = 4;
       fprintf(fidout,'\n 4 Inner Constraints');
       fprintf(fidout,'\n  (2 position, 1 rotation, 1 scale)');
   else
       c = 3;
       fprintf(fidout,'\n 3 Inner Constraints');
       fprintf(fidout,'\n  (2 position, 1 rotation)');
   end    
else    
    fprintf(fidout,'\n%2d constrained bearings',cbrgs);
    fprintf(fidout,'\n%2d constrained distances',cdists);
    fprintf(fidout,'\n%2d constrained angles',cangles);
end    
fprintf(fidout,'\n%2d orientation constants',orients);
fprintf(fidout,'\n\ndegrees of freedom = observations');
fprintf(fidout,'\n                   - unknowns');
fprintf(fidout,'\n                   + constraints');
fprintf(fidout,'\n                   = %2d',n-u+c);
fprintf(fidout,'\n%2d iterations required for solution',JOB(1).iter);
fprintf(fidout,'\n\n');

fclose(fidout);
