function VGPPair = IncDec2VGP(IncGr,DecGr,GeoLatGrad,GeoLonGrad)
%IncDec2VGP Values  Transform ID-data to VGP data.
%   all input values are scalars
% formulas and notations from  R.T.Merril, M.W.McElhinny, P.L. McFadden 
% "The magnetic field of the Earth" pp94-95
%   author: Andrey Khokhlov. 
%   $Revision: 0.25 $  $Date: 2011/10/05 18:15:33 $
    
         TO_GRAD=180.0/pi;       
         AA=atan(0.5*tan(IncGr/TO_GRAD)); %lambda -- paleolatitude theta = _PI/2 - lambda
         
         BB=asin(sin(GeoLatGrad/TO_GRAD)*sin(AA)+cos(GeoLatGrad/TO_GRAD)*cos(AA)*cos(DecGr/TO_GRAD)); %lambda_p
         VGPPair(1) =  BB*TO_GRAD ; %lattitude
         CC=asin(cos(AA)*sin(DecGr/TO_GRAD)/cos(BB)); %beta
         if sin(AA) < sin(GeoLatGrad/TO_GRAD)*sin(BB) 
             VGPPair(2)=GeoLonGrad +180.0 -CC*TO_GRAD;
         else
              VGPPair(2)=GeoLonGrad +CC*TO_GRAD ;
         end
         
         if VGPPair(2) > 360.0
             VGPPair(2)=VGPPair(2)-360.0;
         end
         
         if VGPPair(2) < 0.0
             VGPPair(2)=VGPPair(2)+360.0;
         end    
         return
         
%          
%          AA=atan(0.5*tan(IncGr/TO_GRAD)); 
%          
%          BB=asin(sin(GeoLatGrad/TO_GRAD)*sin(AA) +cos(GeoLatGrad/TO_GRAD)*cos(AA)*cos(DecGr/TO_GRAD)); 
%          VGPLatGr =  BB*TO_GRAD ; %lattitude
%          CC=asin(cos(AA)*sin(DecGr/TO_GRAD)/cos(BB)); %beta
% 
%          VGPLonGr = sin(AA) < sin(GeoLatGrad/TO_GRAD)*sin(BB)?
%                       GeoLonGrad +180.0 -CC*TO_GRAD : GeoLonGrad +CC*TO_GRAD ;
%          
%          if(VGPLonGr > 360.0) VGPLonGr-= 360.0;
%          if(VGPLonGr < 0.0) VGPLonGr+= 360.0;
