% Usage:
% RecFun_GUI('LKR.sac')
% where LKR is a binary sac file
%
% script to interactively calculate receiver
% functions and compare with observed data
% uses reflectivity code written by S. Day
% code can handle N layers but gui is currently
% restricted to two layers and a halfspace
% Assumes Poisson ratio betwene P and S
% Reads in parameters from a binary sac file
% aligns observed and model on maximum arrival
%
% RJM Feb. 2008
%
function varargout = RecFun_GUI(varargin)
% Last Modified by GUIDE v2.5 20-Feb-2008 09:25:02
% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @RecFun_GUI_OpeningFcn, ...
                   'gui_OutputFcn',  @RecFun_GUI_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin & isstr(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT 
% --- Executes just before RecFun_GUI is made visible.
function RecFun_GUI_OpeningFcn(hObject, eventdata, handles, varargin)
% hObject    handle to figure
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to RecFun_GUI (see VARARGIN)
%
% Choose default command line output for RecFun_GUI
    handles.output = hObject;
    fprintf ('n arg %d \n', size(varargin,2));
    if (size(varargin,2) == 0)
        error('usage: RecFun_GUI(sacfile)');
    end 

    fprintf( 'reading %s\n', varargin{1});

% read sac data file and normalize data
    [handles.obs_data, handles.dt, handles.rayp, handles.npts] = ReadSac(varargin{1});
    handles.obs_data = handles.obs_data ./max(handles.obs_data);
% --------------------------------------------------------

% set default layer thickness for upper and lower crust; 0 is halfspace
% layer P and S velocities, Q and densities
    handles.h = [5 30 0];
    handles.a = [5.0 6.5 8.0];
    handles.rho = [2.8 2.8 3.3];
    
% calculate theta from ray parameter; ReadSac assumes it is in user8
% if not, use 0.05 as default
    if (handles.rayp < 0)
            handles.rayp = 0.05;
    end
    handles.theta = ((180/pi)*asin(handles.rayp*handles.a(2)));
    fprintf('Data: dt %6.3f rayp %6.3f ', handles.dt, handles.rayp);
    fprintf(' npts %d angle %6.2f\n', handles.npts, handles.theta);
    
% run and plot default model; update handle structure
    [handles] = run_model(handles);
    guidata(hObject, handles);

%----------------------------------------------------------
% Now allow for user to change parameters with sliders
%
% Upper crust parameters 
% slider1 sets thickness, slider2 is P velocity

    function slider1_CreateFcn(hObject, eventdata, handles)
    function slider1_Callback(hObject, eventdata, handles)
    handles.h(1) = get(handles.slider1,'Value');
    
    % run and plot model
    [handles] = run_model(handles);
    guidata(hObject, handles);
    
    % slider2 sets P wave velocity
    function slider2_CreateFcn(hObject, eventdata, handles)
    function slider2_Callback(hObject, eventdata, handles)
    handles.a(1) = get(handles.slider2,'Value');
    
    [handles] = run_model(handles);
    guidata(hObject, handles);
    
% end upper crust properties ---------------------------------
    
% Lower crust properties
% slider3 is thickness, slider4 is P wave velocity
% thickness
    function slider3_CreateFcn(hObject, eventdata, handles)
    function slider3_Callback(hObject, eventdata, handles)
    handles.h(2) = get(handles.slider3,'Value');

    % run and plot model
    [handles] = run_model(handles);
    guidata(hObject, handles);

% P wave velocity
    function slider4_CreateFcn(hObject, eventdata, handles)
    function slider4_Callback(hObject, eventdata, handles)

    handles.a(2) = get(handles.slider4,'Value');

    % run and plot model
    [handles] = run_model(handles);
    guidata(hObject, handles);
    
% end lower crust
%--------------------------------------------------------------
% handles output
function varargout = RecFun_GUI_OutputFcn(hObject, eventdata, handles)
varargout{1} = handles.output;
%---------------------------------------------------------------
%
% runs, plots, and prints model; returns modified handles structure
function [handles] = run_model(handles)
% set dependent parameters
    handles.b = (1/sqrt(3)).*handles.a;    
    handles.qs = 100.*handles.b;
    handles.qp = 200.*handles.b;

    handles.model = [handles.h; handles.a; handles.b; handles.rho; handles.qp; handles.qs];
    [handles.time, handles.RF] = RecFun(handles.theta, handles.model, handles.dt, handles.npts);
    hold off
    
    % plot model receiver function (blue)
    plot(handles.time, handles.RF,'b')
    hold on
    N = max(size(handles.time));
    
    % line up model and observed maximum value
    [xmax, ymax1] = max(handles.obs_data);
    [xmax, ymax2] = max(handles.RF);
    diff = ymax1 - ymax2;
    
    % plot observed receiver function (red) and calculate L1 error
    plot(handles.time, handles.obs_data(1+diff:N+diff),'r');
    err = sum(abs(handles.obs_data(1+diff:N+diff) - (handles.RF')));
    
    % print out crustal parameters
    fprintf ('------------------------------\n');
    fprintf('Moho depth: %6.0f km    Error: %4.1f\n', handles.h(1) + handles.h(2), err);
    if (handles.a(1) > handles.a(2)) 
        fprintf('****Warning: low velocity zone!\n');
        fprintf('****Upper crust Vp > lower crust Vp\n');
    end
    for i = 1:3
        fprintf ('Layer %d depth: %4.1f Vp: %6.2f Vs: %6.2f', i, handles.h(i), handles.a(i), handles.b(i));
        fprintf (' rho: %5.2f Qp: %4.0f Qs: %4.0f\n', handles.rho(i), handles.qp(i),handles.qs(i));
    end

%--------------------------------------------------------------

