function [Lat_m,Lon_m]=geogr2geomag(Lat_g,Lon_g,Lat_p,Lon_p)
% Вычисление геомагнитных координат по географическим и координатам
% геомагнитного полюса
% Вызов
%        [Lat_m,Lon_m]=geogr2geomag(Lat_g,Lon_g,Lat_p,Lon_p);
%
% Lat_g,Lon_g - географические широта и долгота пункта;
% Lat_p,Lon_p - географические широта и долгота геомагнитного полюса;  
%               (д.б. вектором размерности Lat_g или скалярной величиной)
% Lat_m,Lon_m - геомагнитные широта и долгота пункта
%
% Долгота - в диапазоне -180...+180 град. !!!
% Все величины - в градусах !!!!

% 27.08.2013 - дополнено до варианта с векторами
% 31.01.2022 - использую ф-ции преобразования из Map Toolbox

% параметры преобразования из старого полюса в новый
origin = newpole(Lat_p,Lon_p,'degrees');

[Lat_m,Lon_m]=rotatem(Lat_g,Lon_g,origin,'forward','degrees');
Lon_m=Lon_m+180;

% замена 31.01.2022
% if length(Lat_p)==1,    % задан скаляр
%     Lat_p=ones(length(Lat_g),1)*Lat_p;
%     Lon_p=ones(length(Lat_g),1)*Lon_p;
% end
% 
% Lat_m=asin(sin(d2r(Lat_g)).*sin(d2r(Lat_p))+...
%            cos(d2r(Lat_g)).*cos(d2r(Lat_p)).*cos(d2r(Lon_g-Lon_p)));
% Lat_m=r2d(Lat_m);
% 
% Lon_m=asin(cos(d2r(Lat_g)).*sin(d2r(Lon_g-Lon_p))./cos(d2r(Lat_m)));
% Lon_m=r2d(Lon_m);

% % Дополнение широты (90-широта) и долгота -> в радианах
% % Tgr=d2r(90-Lat_g);
% % Lgr=d2r(Lon_g);
% % Tpr=d2r(90-Lat_p);
% % Lpr=d2r(Lon_p);
% 
% % Tmr=acos(cos(Tpr).*cos(Tgr)+sin(Tpr).*sin(Tgr).*cos(Lgr-Lpr));
% % Lat_m=90-r2d(Tmr);
% % 
% % Lmr=asin(sin(Tgr)./sin(Tmr).*sin(Lgr-Lpr));
% % Lon_m=r2d(Lmr);

