function [lon,lat]=arc_p(lon0,lat0,s,teta1,teta2,n)
%Программа вычисл?ет вектора долгот и широт дуги круга на сфере
%
%Вызов:  [lon,lat]=arc_p(lon0,lat0,s,[teta1,teta2,[n]])
%
%  lon0,lat0 - долгота и широта центра окружности (в град.)
%          s - радиус окружности (в град.)
%teta1,teta2 - начальный и конечный углы, определ?ющие положение
%              дуги на окружности (0 - направление на север;
%              в град.). По умолчанию - 0 и 360.
%          n - количество точек на дуге (по умолчанию - 180);
%
%    lon,lat - векторы долгот и широт размерности n, определ?ющие
%              заданную дугу (в град.)
%
%Долгота измен?етс? от 0 до 360 град., широта - от -90(S) до +90(N).
%Требуема? дуга затем может быть нарисована.
%

if nargin<4
   teta1=0; teta2=360;n=180;
elseif nargin<6
   n=180;
end

% ....... вектор углов поворота в рад.
TETA=(d2r(teta1):d2r((teta2-teta1)/n):d2r(teta2))';
lr = d2r(lon0);
fr = d2r(lat0);
sr = d2r(s);

if lat0==+90 | lat0==-90      % центр - на полюсах
   lon = r2d(TETA);
   lat = zeros(length(TETA),1)+(90-s);
else
   lat = asin( sin(fr)*cos(sr)+cos(fr)*sin(sr)*cos(TETA) ); % текуща? широта
   k = lat==pi/2 | lat==-pi/2;
   cl = ( cos(sr)-sin(lat)*sin(fr) )./( cos(lat)*cos(fr) ); % cos долготы
   sl = sin(sr)*sin(TETA)./cos(lat);                        % sin долготы
   lat = r2d(lat);
   lon=r2d(lr+atan2(sl,cl));
end
