function [d1,d2,d3] = Span123(phi1,lambda1,phi2,lambda2)
R=6372795; % Радиус в метрах
% все значения широты и долготы должны быть в радианах, а для корректного
% расчёта надо использовать atan для исключения особых точек
% Перевод в радианы pi/180*degrees или встроенной функцией
lambda=lambda2-lambda1;
sigma1=atan(sqrt((cos(phi2).*sin(lambda)).^2+(cos(phi1).*sin(phi2)-sin(phi1).*cos(phi2).*...
    cos(lambda)).^2)./(sin(phi1).*sin(phi2)+cos(phi1).*cos(phi2).*cos(lambda)));
d1=sigma1.*R;
sigma2=acos(sin(phi1).*sin(phi2)+cos(phi1).*cos(phi2).*cos(lambda));
d2=sigma2.*R;
sigma3=2.*asin(sqrt(sin((phi2-phi1)./2).^2+cos(phi1).*cos(phi2).*sin(lambda./2).^2));
d3=sigma3.*R;
end

>> [~,a] = Span123(degtorad(55.0749028368),degtorad(1),degtorad(55.0750255180),degtorad(1))
 
a =
 
  13.645487528257137

function[x,y]=index(Dx,Dy);  % функция вычисления индексов
A=load ('wer.nvg');% начальный файл
R=927; Dx = sqrt(R^2 -A(:,5).^2)% индекс по Х
dR=1.4;n=20;Dy=dR*n% индекс по Y
end

clc; clear all;
A=load ('wer.nvg');% начальный файл
R=927;%минимальное расстояние
global dR
dR=1.4;% разрешение
a = 1; % в строке
b = 300; % в столбце
h = 20;
w = 1024;
L = zeros(h*b,w*a); % пустая матрица
Dy=0;
for i = 0:b
        fname = [num2str(i),'.bmp']; % текущее имя файла
        img = imread(fname); % считываем файл
        Dx = sqrt(R^2 -A(i+1,5).^2);% размер кадра по Х
        [x,y]=index(Dx,Dy);
        !!!!L(x:x+w,y:) = img;
        Dy=Dy+dR*size(img,1);
end
figure,imshow(uint8(L));
title('L')