% load rstt C library
rstt_root   = fullfile(pwd, '..');  % '/path/to/rstt'
slbmc       = 'libslbmCshell';  % name of library file
slbmc_dll   = fullfile(rstt_root, 'lib','libslbmCshell.dylib');  % filename
slbmc_hfile = fullfile(rstt_root, 'SLBM_C_shell','include','slbm_C_shell.h');  % h file
addpath(rstt_root);
setenv('RSTT_ROOT', rstt_root);
loadlibrary(slbmc_dll, slbmc_hfile);

% make a struct to hold pointers
p = struct;

% declare path to the model file
modelPath = fullfile(rstt_root, 'models', 'pdu202009Du.geotess');

% choose an arbitrary phase (Pn, Pg, Sn, or Lg)
phase = 'Pn';

% arbitrary source (Meteor Crater, AZ)
srcLatDeg =   35.0274;  % latitude (degrees)
srcLonDeg = -111.0228;  % longitude (degrees)
srcDepKm  =    1.2345;  % depth (km)

% arbitrary receiver (Albuquerque, NM)
rcvLatDeg =   35.1053;  % latitude (degrees)
rcvLonDeg = -106.6294;  % longitude (degrees)
rcvDepKm  =   -1.6000;  % depth (km)

% convert lat/lon from degrees to radians
DEG_PER_RAD = 180 / pi;
RAD_PER_DEG = pi / 180;
srcLatRad = srcLatDeg * RAD_PER_DEG;
srcLonRad = srcLonDeg * RAD_PER_DEG;
rcvLatRad = rcvLatDeg * RAD_PER_DEG;
rcvLonRad = rcvLonDeg * RAD_PER_DEG;

% create a structure to hold pointers
p = struct( ...
    'modelPath',             libpointer('cstring',   modelPath), ...
    'phase',                 libpointer('cstring',   phase),     ...
    'srcLatDeg',             libpointer('doublePtr', srcLatDeg), ...
    'srcLonDeg',             libpointer('doublePtr', srcLonDeg), ...
    'srcDepKm',              libpointer('doublePtr', srcDepKm),  ...
    'rcvLatDeg',             libpointer('doublePtr', rcvLatDeg), ...
    'rcvLonDeg',             libpointer('doublePtr', rcvLonDeg), ...
    'rcvDepKm',              libpointer('doublePtr', rcvDepKm),  ...
    'srcLatRad',             libpointer('doublePtr', srcLatRad), ...
    'srcLonRad',             libpointer('doublePtr', srcLonRad), ...
    'rcvLatRad',             libpointer('doublePtr', rcvLatRad), ...
    'rcvLonRad',             libpointer('doublePtr', rcvLonRad), ...
    'distRad',               libpointer('doublePtr', nan),       ...
    'travelTimeSec',         libpointer('doublePtr', nan),       ...
    'travelTimeUncertSec',   libpointer('doublePtr', nan),       ...
    'travelTime1DUncertSec', libpointer('doublePtr', nan)        ...
    );


% instantiate an RSTT object
ret = calllib(slbmc, 'slbm_shell_create');
% [ret, errorMsg] = calllib(slbmc, 'slbm_shell_getErrorMessage', ' ');  % example call to find out error message

% load the velocity model
ret = calllib(slbmc, 'slbm_shell_loadVelocityModel', p.modelPath);

% create a great circle from source to the receiver
ret = calllib(slbmc, 'slbm_shell_createGreatCircle', p.phase, ...
    p.srcLatRad, p.srcLonRad, p.srcDepKm, ...
    p.rcvLatRad, p.rcvLonRad, p.rcvDepKm  );

% get the distance and travel time from source --> receiver
ret = calllib(slbmc, 'slbm_shell_getDistance',   p.distRad);        % compute distance (rad) and store its value
ret = calllib(slbmc, 'slbm_shell_getTravelTime', p.travelTimeSec);  % compute travel time (sec) and store its value
distDeg = p.distRad.Value * DEG_PER_RAD;  % convert radians --> degrees

% get the travel time uncertainty
ret = calllib(slbmc, 'slbm_shell_getTTUncertainty',   p.travelTimeUncertSec);    % compute uncertainty
ret = calllib(slbmc, 'slbm_shell_getTTUncertainty1D', p.travelTime1DUncertSec);  % compute 1D uncertainty

% print our results
fprintf(['\n'                            ...
    'Matlab Example:\n'                  ...
    '--------------------------------\n' ...
    'Distance:          %9.4f deg\n'     ...
    'Travel Time:       %9.4f sec\n'     ...
    'Uncertainty:       %9.4f sec\n'     ...
    'Uncertainty (1D):  %9.4f sec\n'     ...
    '\n'], distDeg, p.travelTimeSec.Value, p.travelTimeUncertSec.Value, p.travelTime1DUncertSec.Value);
