%% script to read shear wave splitting database text files:
% this is only an example of how to handle the file format and how to
% generate some simple output in Matlab. If you come up with a more
% sophisticated script, and want to share it with others, don't hesitate 
% to contact the database Administrators. We might be putting it online ;-)

% A. Wuestefeld, March 2007



% file locations:__________________________________________________________
ref_file = 'Q:\PhD\SplitlabHP\DB\admin\data\splittingDB_ref.txt';
DB_file = 'Q:\PhD\SplitlabHP\DB\admin\data\splittingDB.txt';


% read references__________________________________________________________
[RefNo, Tag, Authors, Year, Title, Journal, No, Pages, DOI, Link, Comment]=...
    textread(ref_file, '%u%s%s%f%s%s%s%s%s%s%s',... 
                       'headerlines', 1,...
                       'whitespace',  '',...
                       'delimiter',  '|');


% read splitting data______________________________________________________
% We have to read "phi" into a temporary STRING to account for NULL stations
[id, Station, Lat, Long, tmp, dt, refID, Phase, pmax, pmin, dtmax, dtmin, remark]=...
    textread(DB_file, '%u%s%f%f%s%f%u%s%f%f%f%f%s',...
                      'headerlines', 1,...
                      'whitespace',  '',...
                      'delimiter',   '|',...
                      'emptyvalue',  nan);

%% now convert temp-phi into double-precission:
NULL    = find(strcmp(lower(tmp), 'null')); %find isotropic stations
NoValue = find(strcmp(lower(tmp), ' '));    
empty   = union(NoValue, NULL);

for k = 1:length(tmp)
    if ismember(k, empty);
        phi(k,1) = nan;
    else
        phi(k,1) = sscanf(tmp{k} ,'%f');
    end
end
clear tmp k 






%%  ___________plotting examples:__________________________________________
% 
%% XXXXXXXXXXXX FIGURE 1 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
%% Plot coastline (requires Mapping Toolbox)
figure(1)
load coast
plot(long, lat, 'k-')
hold on


%% plot Null stations as red circles;
plot(Long(empty), Lat(empty), 'ro');


%% plot other (non-Null) stations as blue crosses:
other = setdiff(1:length(Lat), empty);
plot(Long(other), Lat(other), 'bx');
hold off


%beautify
axis equal
axis([-180 180 -90 90])




%% XXXXXXXXXXXX FIGURE 2 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
%% make simple fast axis plot
figure(2)
plot(long, lat, 'k-')
hold on


% set length of marker: factor = 1 means that one second delay is displayed
% as one degree on globe
factor = 2; 
x = sin(phi(other)*pi/180) .* dt(other) * factor/2; %YES, division by 2!!!
y = cos(phi(other)*pi/180) .* dt(other) * factor/2; 
XX = [Long(other)-x   Long(other)+x   nan(size(x))]';
YY = [Lat(other)-y    Lat(other)+y    nan(size(y))]';

plot(XX, YY, 'm-', Long(empty), Lat(empty), 'b.');

hold off


%beautify
axis equal
axis([-180 180 -90 90])


%% XXXXXXXXXXXX FIGURE 3 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
% histogram plot; similar to that in Silver, 1996
figure (3)
subplot(2,1,1)
  hist(dt,40)
  title('Global delay time distribution')
  
subplot(2,1,2) 
  phi2 = mod(phi+90, 180)-90;
  hist(phi2, 45) % somewhat pointless, but for the fun of it...
  xlim([-90 90])
  title('Global fast axis azimuth distribution')