%% VGP Calculator using array data of Inc and Dec
% RSD ==>    SLat = -21.1665      SLong = 55.5029
% JC  ==>    SLat = -49.7         SLong = 150.5
%
% Virtual Geomagnetic Pole (VGP) Calculator 
% Takes Site Latitude and Longitude as well as
% samples Inclination and Declination and
% calculates VGP Latitude and Longitude

%
clear all; close all; clc;


%
disp('Enter site latitude: ');
SLat = input('');
disp('Now enter site Longitude: ');
SLong = input('');
disp(' The data-file must have [Dec Inc a95] for each flow.');
disp(' Enter data-file name:');
fname = input('','s');
data = load(fname);

%
figure(1)
worldmap world
geoshow('landareas.shp', 'FaceColor', [0.5 1.0 0.5]);
hold on;
a95 = 0;

%
[nrow, ncol, na95] = size(data);
VgpOutData = zeros(nrow,ncol,na95);
dVgpOutData = zeros(nrow,2);

for q = 1:nrow
    
    SamInc = data(q,2) ;    %Sample Inclination
    SamDec = data(q,1) ;    %Sample Declination
    a95 = data(q,3);
    %
    InQ = (2/tand(SamInc));
    p = atand(InQ);
    % Here, I compensated for angular distance p to be 
    % a positive parameter from 0 to 180 (Nate)
    
    if p < 0
        p = p + 180;
    elseif p > 180
        p = p - 180;
    end
    % Here I (Nate) added uncertainties in VGP Lat and Lon.
    dp = a95*(0.5)*(1+3*cosd(p)*cosd(p));
    
    dm = a95*sind(p)/cosd(SamInc);
    
    dVgpOutData(q,1) = dp;
    dVgpOutData(q,2) = dm;
    
    t1=sind(SLat); %fine
    t2=cosd(p);
    t3=cosd(SLat); %fine
    t4=sind(p);
    t5=cosd(SamDec); 
    
    VagInput = (((t1)*(t2)) + ((t3)*(t4)*(t5)));
    
    p1=(t1)*(t2);
    p2= (t3)*(t4)*(t5);
    
    
    VgpLatitude = asind( VagInput );
    
    Beta = asind((sind(p))*(sind(SamDec))/(cosd(VgpLatitude))) ;
    
    
    if cosd(p) >= sind(SLat)*sind(VgpLatitude);
        VgpLongitude = SLong + Beta;
    else
        VgpLongitude = SLong + 180 - Beta;
    end
    % I have to keep 0 < Long < 360 for map
    if VgpLongitude > 360
        VgpLongitude = VgpLongitude - 360;
    end
    %
    VgpOutData(q,1) = VgpLatitude;
    VgpOutData(q,2) = VgpLongitude;
    VgpOutData(q,3) = a95;
    %
    geoshow(VgpOutData(q,1), VgpOutData(q,2), 'DisplayType', 'point', 'Marker', '.' ,'Color', 'k','MarkerSize', 2+q); 
    hold on;
    
    for t = 0:360
        eLat = dm*cosd(t)*cosd(360-SamDec) - dp*sind(t)*sind(360-SamDec) + VgpOutData(q,1);
        eLon = dm*cosd(t)*sind(360-SamDec) + dp*sind(t)*cosd(360-SamDec) + VgpOutData(q,2);
        geoshow(eLat,eLon, 'DisplayType', 'point', 'Marker', '.' ,'Color', 'c','MarkerSize', 1); 
        hold on;
    end
    
end

for kk = 1:nrow-1

    [gclat,gclong] = track2('gc',VgpOutData(kk,1),VgpOutData(kk,2),...
                             VgpOutData(kk+1,1),VgpOutData(kk+1,2));

    plotm(gclat,gclong,'k-'); 

end