// Уравновешивание замкнутого теодолитного хода
// scilab 6.1 
// version: 2021-03-19  2.10

clear
funcprot(0)

// Вспомогательные функции
function v=to_degree(deg, mins, secs)
    // Преобразование из градусов, минут, секунда 
    // в градусы и десятые градуса для sind/cosd
    v = deg + mins/60 + secs/60/60
endfunction

function [deg, mins, secs]=from_degree(v)
    deg  = floor(v)
    mins_secs = v - floor(v)
    
    mins  = floor(mins_secs * 60)
    secs = round((mins_secs * 60 - mins) * 60)
endfunction    

// Функции для печати значений в правильном формате

function print_dgr(label, valvec)
    // Печатает значение углов в градусах и десятых градуса как градусы, минуты, секунды
    function __print_dgr(label, varargin)
    
      if size(varargin) == 1 then
        
          v = varargin(1)
          [deg, mins, secs] = from_degree(v)
          printf("%s: %d°%d''%d""\n", label, deg, mins, secs)

      elseif size(varargin) == 2 then   
    
          idx = varargin(1)
          v = varargin(2)
          [deg, mins, secs] = from_degree(v)
          printf("%s %d: %d°%d''%d""\n", label, idx, deg, mins, secs)

      else
        assert_checktrue ( false );    
      end
    endfunction

    // Печатает вектор углов
    for i = 1 : size(valvec, "*")
        __print_dgr(label,i,valvec(i))
    end
endfunction    

function print_xy(label, valvec_x, valvec_y)
    // Печатает вектор координат

    function __print_xy(label, varargin)
       // Печатает координаты
    
      if size(varargin) == 2 then
        
          x = varargin(1)
          y = varargin(2)        
          printf("%s: %0.4f,%0.4f\n", label, x, y)
        
      elseif size(varargin) == 3 then
        
          idx = varargin(1)
          x   = varargin(2)
          y   = varargin(3)        
          printf("%s %d: %0.4f,%0.4f\n", label, idx, x, y)

      else
          assert_checktrue ( false );    
      end        
        
   endfunction

    assert_checkequal (size(valvec_x,"*"), size(valvec_y,"*"));

    for i = 1 : size(valvec_x, "*")
        __print_xy(label,i,valvec_x(i), valvec_y(i))
    end
endfunction    

// *****  Дано ******
// Измеренные горизонтальные углы
angles = [ to_degree(105, 24, 0), to_degree(123, 10, 0), to_degree( 79, 24, 0), to_degree(140, 39, 0), to_degree( 91, 21, 0) ]

// Горизонтальное положение: 12 23 34 45 51 ...
d_hors = [ 135.62, 159.82, 142.15,  138.61, 153.71 ]

// Дирекционные углы, 
// первый должен быть задан, остальные вычисляются
dirs = [ to_degree(44, 44, 0) ]

// Координаты точек, 
// координаты первой точки должны быть заданы, остальные вычисляются
pX = [ 250.00 ]
pY = [-780.00 ]

// *********** Начало вычислений *****
// Количество точек
n_points = size(angles, "*")

// Количество горизонтальных положений должно быть равно количеству точек
assert_checkequal (size(d_hors,"*"), n_points);

// Определяют практическую сумму углов Sпр и угловую невязку f, 
//  равную разности практической и теоретической сумм углов в полигоне 

// Теоретическая сумма углов в полигоне
St = 180 * (n_points - 2)

// Практическая сумма углол в полигоне 
Sp = sum(angles)

// Угловая невязка
f = abs(St - Sp)

// Промежуточная печать
print_dgr("Теоретическая сумма уголов", St)
print_dgr("Практическая сумма уголов", Sp)
print_dgr("Угловая невязка", f)

// Угловую невязку f распределяют с обратным знаком поровну на все измеренные углы, с точностью до 0',1.
// Поправка
ppr = f/n_points
print_dgr("Поправка", ppr)

// Углы с поправкой
angles = angles + ppr

// Вычисляют дирекционные углы (азимуты) линий.
// dir_1 - pre-calculated, don't set

for i = 2 : n_points
    dirs(i) = dirs(i - 1) + 180 - angles(i)
end    

print_dgr("Дирекционный угол", dirs)

// Вычисляем румбы, в дальнейших расчетах они не участвуют
for i = 1 : n_points
    if dirs(i) <= 90 then
        rumbs(i) = dirs(i)
    elseif dirs(i) > 90 && dirs(i) <= 180 then   
        rumbs(i) = 180 - dirs(i)
    elseif dirs(i) > 180 && dirs(i) <= 270 then   
        rumbs(i) = dirs(i) - 180
   elseif dirs(i) > 270 then   
        rumbs(i) = 360 - dirs(i)
   end
end   

print_dgr("Румб", rumbs)

// Приращение координат
for i = 1 : n_points
    dX(i) = d_hors(i) * cosd(dirs(i))
    dY(i) = d_hors(i) * sind(dirs(i))
end

print_xy("Приращение координат", dX, dY)

// Невязка равна сумме приращений 
FdX = sum(dX)
FdY = sum(dY)

print_xy("Невязка", FdX, FdY)

// Поправки
Sd = sum(d_hors)

for i = 1 : n_points
  sigma_X(i) = (FdX * d_hors(i))/ Sd 
  sigma_Y(i) = (FdY * d_hors(i))/ Sd 
end

print_xy("Поправки", sigma_X, sigma_Y)

// Приращение координат с поправками 
for i = 1: n_points
    dXc(i) = dX(i) + sigma_X(i) * -1
    dYc(i) = dY(i) + sigma_Y(i) * -1
end

print_xy("Приращение координат с поправками", dXc, dYc)

// Невязка равна сумме приращений 
FdXc = sum(dXc)
FdYc = sum(dYc)

print_xy("Невязка с поправками", FdXc, FdYc)

// Невязка с поправками должна быть равна нулю
assert_checkequal (round(FdXc), 0);
assert_checkequal (round(FdYc), 0);

// Вычисляем координаты
for i = 2 : n_points
   pX(i) = pX(i-1) + dXc(i-1)
   pY(i) = pY(i-1) + dYc(i-1)
end

// Полигон должен быть замкнутым
pX(n_points + 1) = pX(n_points) + dXc(n_points)
pY(n_points + 1) = pY(n_points) + dYc(n_points)

print_xy("Координаты", pX, pY)

assert_checkequal (round(pX(n_points + 1)), pX(1));
assert_checkequal (round(pY(n_points + 1)), pY(1));

plot(pX, pY, "b")

// The End

Details

Published: 22 March 2021

Hits: 1131