%% 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