function varargout = QDControl(varargin)
% QDCONTROL M-file for QDControl.fig
%      QDCONTROL, by itself, creates a new QDCONTROL or raises the existing
%      singleton*.
%
%      H = QDCONTROL returns the handle to a new QDCONTROL or the handle to
%      the existing singleton*.
%
%      QDCONTROL('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in QDCONTROL.M with the given input arguments.
%
%      QDCONTROL('Property','Value',...) creates a new QDCONTROL or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before QDControl_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to QDControl_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help QDControl

% Last Modified by GUIDE v2.5 22-Feb-2005 12:49:48

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @QDControl_OpeningFcn, ...
                   'gui_OutputFcn',  @QDControl_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 QDControl is made visible.
function QDControl_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to QDControl (see VARARGIN)

% Choose default command line output for QDControl
handles.output = hObject;

% UIWAIT makes QDControl wait for user response (see UIRESUME)
% uiwait(handles.figure1);

handles.n0 = 1.52;
handles.n1 = 1.0;
handles.lambda = 0.57;

nn = 20;
handles.pix = 6.45;
[x,y] = meshgrid(-nn:nn,-nn:nn);
handles.p = angle(x+i*y);
handles.r = sqrt(x.^2+y.^2);

handles.kappa = str2double(get(handles.kappa_edit,'String'));
handles.ratio = str2double(get(handles.ratio_edit,'String'));
handles.omega = str2double(get(handles.omega_edit,'String'));
handles.theta = str2double(get(handles.theta_edit,'String'));
handles.phi = 0;
if nargin 
    pos = find(strcmpi('kappa', varargin));
    if ~isempty(pos)
        handles.kappa = varargin{pos+1};
        set(handles.kappa_edit,'String',num2str(handles.kappa));
    end
    pos = find(strcmpi('ratio', varargin));
    if ~isempty(pos)
        handles.ratio = varargin{pos+1};
        set(handles.ratio_edit,'String',num2str(handles.ratio));
    end
    pos = find(strcmpi('omega', varargin));
    if ~isempty(pos)
        handles.omega = varargin{pos+1};
        set(handles.omega_edit,'String',num2str(handles.omega));
    end
    pos = find(strcmpi('theta', varargin));
    if ~isempty(pos)
        handles.theta = varargin{pos+1};
        set(handles.theta_edit,'String',num2str(handles.theta));
    end
    pos = find(strcmpi('phi', varargin));
    if ~isempty(pos)
        handles.phi = varargin{pos+1};
    end
end

set(handles.theta_slider,'sliderstep',1/90*[1 10],'max',90,'min',0,'Value',handles.theta);
set(handles.omega_slider,'sliderstep',1/180*[1 10],'max',180,'min',0,'Value',handles.omega);
set(handles.kappa_slider,'sliderstep',0.01*[1 5],'max',1,'min',-1,'Value',handles.kappa);
set(handles.ratio_slider,'sliderstep',0.02*[1 5],'max',1,'min',0,'Value',handles.ratio);
set(handles.phi_slider,'sliderstep',1/360*[1 10],'max',360,'min',0,'Value',handles.phi);

handles.mag = str2double(get(handles.mag_edit,'String'));
handles.focus = str2double(get(handles.foc_edit,'String'));
handles.na = str2double(get(handles.na_edit,'String'));
if nargin
    pos = find(strcmpi('mag', varargin));
    if ~isempty(pos)
        handles.mag = varargin{pos+1};
        set(handles.mag_edit,'String',num2str(handles.mag));
    end
    pos = find(strcmpi('focus', varargin));
    if ~isempty(pos)
        handles.focus = varargin{pos+1};
        set(handles.foc_edit,'String',num2str(handles.focus));
    end
    pos = find(strcmpi('na', varargin));
    if ~isempty(pos)
        handles.na = varargin{pos+1};
        set(handles.na_edit,'String',num2str(handles.na));
    end
end

[tmp, tmp, tmp, rho, tmp, fxx0, fxx2, fxz, byx0, byx2, byz] = ...
    SEPDipole([0 2],0,handles.na,handles.n0,handles.n1,handles.n1,[],0,[],handles.lambda,handles.mag,handles.focus);

f00 = fxx0.*conj(byx0);
f01 = -fxx0.*conj(byz);
f02 = -fxx0.*conj(byx2);
f10 = -fxz.*conj(byx0);
f11 = fxz.*conj(byz);
f12 = fxz.*conj(byx2);
f20 = -fxx2.*conj(byx0);
f21 = fxx2.*conj(byz);
f22 = fxx2.*conj(byx2);
handles.rho = rho/handles.pix;
handles.f00 = f00;
handles.f01 = f01;
handles.f02 = f02;
handles.f10 = f10;
handles.f11 = f11;
handles.f12 = f12;
handles.f20 = f20;
handles.f21 = f21;
handles.f22 = f22;

handles.int = int_compute(handles.theta, handles.omega, handles.kappa, handles.ratio, handles.rho, handles.r, handles.p, ...
    handles.f00, handles.f01, handles.f02, handles.f10, handles.f11, handles.f12, handles.f20, handles.f21, handles.f22, handles.phi, get(handles.radiobutton1,'Value'));

if nargin
    pos = find(strcmpi('image',varargin));
    if ~isempty(pos)
        axes(handles.imm);
        mim(varargin{pos+1});
    end
end
axes(handles.pic);
mim(handles.int); title(int2str(round(handles.phi)),'fontsize',12)

% Update handles structure
guidata(hObject, handles);


% --- Outputs from this function are returned to the command line.
function varargout = QDControl_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;



% --- Executes during object creation, after setting all properties.
function na_edit_CreateFcn(hObject, eventdata, handles)
% hObject    handle to na_edit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end

set(hObject, 'String', '1.2');

function na_edit_Callback(hObject, eventdata, handles)
% hObject    handle to na_edit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of na_edit as text
%        str2double(get(hObject,'String')) returns contents of na_edit as a double
handles.na = str2double(get(hObject,'String'));
[tmp, tmp, tmp, rho, tmp, fxx0, fxx2, fxz, byx0, byx2, byz] = ...
    SEPDipole([0 2],0,handles.na,handles.n0,handles.n1,handles.n1,[],0,[],handles.lambda,handles.mag,handles.focus);
f00 = fxx0.*conj(byx0);
f01 = -fxx0.*conj(byz);
f02 = -fxx0.*conj(byx2);
f10 = -fxz.*conj(byx0);
f11 = fxz.*conj(byz);
f12 = fxz.*conj(byx2);
f20 = -fxx2.*conj(byx0);
f21 = fxx2.*conj(byz);
f22 = fxx2.*conj(byx2);
handles.f00 = f00;
handles.f01 = f01;
handles.f02 = f02;
handles.f10 = f10;
handles.f11 = f11;
handles.f12 = f12;
handles.f20 = f20;
handles.f21 = f21;
handles.f22 = f22;
handles.int = int_compute(handles.theta, handles.omega, handles.kappa, handles.ratio, handles.rho, handles.r, handles.p, ...
    handles.f00, handles.f01, handles.f02, handles.f10, handles.f11, handles.f12, handles.f20, handles.f21, handles.f22, handles.phi, get(handles.radiobutton1,'Value'));
mim(handles.int); title(int2str(round(handles.phi)),'fontsize',12)
guidata(hObject, handles);


% --- Executes during object creation, after setting all properties.
function mag_edit_CreateFcn(hObject, eventdata, handles)
% hObject    handle to mag_edit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end
set(hObject, 'String', '110');


function mag_edit_Callback(hObject, eventdata, handles)
% hObject    handle to mag_edit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of mag_edit as text
%        str2double(get(hObject,'String')) returns contents of mag_edit as a double
handles.mag = str2double(get(hObject,'String'));
[tmp, tmp, tmp, rho, tmp, fxx0, fxx2, fxz, byx0, byx2, byz] = ...
    SEPDipole([0 2],0,handles.na,handles.n0,handles.n1,handles.n1,[],0,[],handles.lambda,handles.mag,handles.focus);
f00 = fxx0.*conj(byx0);
f01 = -fxx0.*conj(byz);
f02 = -fxx0.*conj(byx2);
f10 = -fxz.*conj(byx0);
f11 = fxz.*conj(byz);
f12 = fxz.*conj(byx2);
f20 = -fxx2.*conj(byx0);
f21 = fxx2.*conj(byz);
f22 = fxx2.*conj(byx2);
handles.rho = rho/handles.pix;
handles.f00 = f00;
handles.f01 = f01;
handles.f02 = f02;
handles.f10 = f10;
handles.f11 = f11;
handles.f12 = f12;
handles.f20 = f20;
handles.f21 = f21;
handles.f22 = f22;
handles.int = int_compute(handles.theta, handles.omega, handles.kappa, handles.ratio, handles.rho, handles.r, handles.p, ...
    handles.f00, handles.f01, handles.f02, handles.f10, handles.f11, handles.f12, handles.f20, handles.f21, handles.f22, handles.phi, get(handles.radiobutton1,'Value'));
mim(handles.int); title(int2str(round(handles.phi)),'fontsize',12)
guidata(hObject, handles);


% --- Executes during object creation, after setting all properties.
function foc_edit_CreateFcn(hObject, eventdata, handles)
% hObject    handle to foc_edit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end
set(hObject, 'String', '1.2');

function foc_edit_Callback(hObject, eventdata, handles)
% hObject    handle to foc_edit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of foc_edit as text
%        str2double(get(hObject,'String')) returns contents of foc_edit as a double
handles.focus = str2double(get(hObject,'String'));
[tmp, tmp, tmp, rho, tmp, fxx0, fxx2, fxz, byx0, byx2, byz] = ...
    SEPDipole([0 2],0,handles.na,handles.n0,handles.n1,handles.n1,[],0,[],handles.lambda,handles.mag,handles.focus);
f00 = fxx0.*conj(byx0);
f01 = -fxx0.*conj(byz);
f02 = -fxx0.*conj(byx2);
f10 = -fxz.*conj(byx0);
f11 = fxz.*conj(byz);
f12 = fxz.*conj(byx2);
f20 = -fxx2.*conj(byx0);
f21 = fxx2.*conj(byz);
f22 = fxx2.*conj(byx2);
handles.f00 = f00;
handles.f01 = f01;
handles.f02 = f02;
handles.f10 = f10;
handles.f11 = f11;
handles.f12 = f12;
handles.f20 = f20;
handles.f21 = f21;
handles.f22 = f22;
handles.int = int_compute(handles.theta, handles.omega, handles.kappa, handles.ratio, handles.rho, handles.r, handles.p, ...
    handles.f00, handles.f01, handles.f02, handles.f10, handles.f11, handles.f12, handles.f20, handles.f21, handles.f22, handles.phi, get(handles.radiobutton1,'Value'));
mim(handles.int); title(int2str(round(handles.phi)),'fontsize',12)
guidata(hObject, handles);


% --- Executes during object creation, after setting all properties.
function omega_edit_CreateFcn(hObject, eventdata, handles)
% hObject    handle to omega_edit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end
set(hObject, 'String', '0');

function omega_edit_Callback(hObject, eventdata, handles)
% hObject    handle to omega_edit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of omega_edit as text
%        str2double(get(hObject,'String')) returns contents of omega_edit as a double
handles.omega = str2double(get(hObject,'String'));
set(handles.omega_slider,'Value',handles.omega);
handles.int = int_compute(handles.theta, handles.omega, handles.kappa, handles.ratio, handles.rho, handles.r, handles.p, ...
    handles.f00, handles.f01, handles.f02, handles.f10, handles.f11, handles.f12, handles.f20, handles.f21, handles.f22, handles.phi, get(handles.radiobutton1,'Value'));
mim(handles.int); title(int2str(round(handles.phi)),'fontsize',12)
guidata(hObject, handles);


% --- Executes during object creation, after setting all properties.
function theta_edit_CreateFcn(hObject, eventdata, handles)
% hObject    handle to theta_edit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end
set(hObject, 'String', '0');

function theta_edit_Callback(hObject, eventdata, handles)
% hObject    handle to theta_edit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of theta_edit as text
%        str2double(get(hObject,'String')) returns contents of theta_edit as a double
handles.theta = str2double(get(hObject,'String'));
set(handles.theta_slider,'Value',handles.theta);
handles.int = int_compute(handles.theta, handles.omega, handles.kappa, handles.ratio, handles.rho, handles.r, handles.p, ...
    handles.f00, handles.f01, handles.f02, handles.f10, handles.f11, handles.f12, handles.f20, handles.f21, handles.f22, handles.phi, get(handles.radiobutton1,'Value'));
mim(handles.int); title(int2str(round(handles.phi)),'fontsize',12)
guidata(hObject, handles);


% --- Executes during object creation, after setting all properties.
function kappa_edit_CreateFcn(hObject, eventdata, handles)
% hObject    handle to kappa_edit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end
set(hObject, 'String', '0');

function kappa_edit_Callback(hObject, eventdata, handles)
% hObject    handle to kappa_edit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of kappa_edit as text
%        str2double(get(hObject,'String')) returns contents of kappa_edit as a double
handles.kappa = str2double(get(hObject,'String'));
set(handles.kappa_slider,'Value',handles.kappa);
handles.int = int_compute(handles.theta, handles.omega, handles.kappa, handles.ratio, handles.rho, handles.r, handles.p, ...
    handles.f00, handles.f01, handles.f02, handles.f10, handles.f11, handles.f12, handles.f20, handles.f21, handles.f22, handles.phi, get(handles.radiobutton1,'Value'));
mim(handles.int); title(int2str(round(handles.phi)),'fontsize',12)
guidata(hObject, handles);


% --- Executes during object creation, after setting all properties.
function ratio_edit_CreateFcn(hObject, eventdata, handles)
% hObject    handle to ratio_edit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end
set(hObject, 'String', '0');

function ratio_edit_Callback(hObject, eventdata, handles)
% hObject    handle to ratio_edit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of ratio_edit as text
%        str2double(get(hObject,'String')) returns contents of ratio_edit as a double
handles.ratio = str2double(get(hObject,'String'));
set(handles.ratio_slider,'Value',handles.ratio);
handles.int = int_compute(handles.theta, handles.omega, handles.kappa, handles.ratio, handles.rho, handles.r, handles.p, ...
    handles.f00, handles.f01, handles.f02, handles.f10, handles.f11, handles.f12, handles.f20, handles.f21, handles.f22, handles.phi, get(handles.radiobutton1,'Value'));
mim(handles.int); title(int2str(round(handles.phi)),'fontsize',12)
guidata(hObject, handles);


% --- Executes during object creation, after setting all properties.
function pic_CreateFcn(hObject, eventdata, handles)
% hObject    handle to pic (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: place code in OpeningFcn to populate pic


% --- Executes during object creation, after setting all properties.
function theta_slider_CreateFcn(hObject, eventdata, handles)
% hObject    handle to theta_slider (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: slider controls usually have a light gray background, change
%       'usewhitebg' to 0 to use default.  See ISPC and COMPUTER.
usewhitebg = 1;
if usewhitebg
    set(hObject,'BackgroundColor',[.9 .9 .9]);
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end


% --- Executes on slider movement.
function theta_slider_Callback(hObject, eventdata, handles)
% hObject    handle to theta_slider (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider
handles.theta = get(hObject,'Value');
set(handles.theta_edit,'String',num2str(handles.theta));
handles.int = int_compute(handles.theta, handles.omega, handles.kappa, handles.ratio, handles.rho, handles.r, handles.p, ...
    handles.f00, handles.f01, handles.f02, handles.f10, handles.f11, handles.f12, handles.f20, handles.f21, handles.f22, handles.phi, get(handles.radiobutton1,'Value'));
mim(handles.int); title(int2str(round(handles.phi)),'fontsize',12)
guidata(hObject, handles);


% --- Executes during object creation, after setting all properties.
function omega_slider_CreateFcn(hObject, eventdata, handles)
% hObject    handle to omega_slider (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: slider controls usually have a light gray background, change
%       'usewhitebg' to 0 to use default.  See ISPC and COMPUTER.
usewhitebg = 1;
if usewhitebg
    set(hObject,'BackgroundColor',[.9 .9 .9]);
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end


% --- Executes on slider movement.
function omega_slider_Callback(hObject, eventdata, handles)
% hObject    handle to omega_slider (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider
handles.omega = get(hObject,'Value');
set(handles.omega_edit,'String',num2str(handles.omega));
handles.int = int_compute(handles.theta, handles.omega, handles.kappa, handles.ratio, handles.rho, handles.r, handles.p, ...
    handles.f00, handles.f01, handles.f02, handles.f10, handles.f11, handles.f12, handles.f20, handles.f21, handles.f22, handles.phi, get(handles.radiobutton1,'Value'));
mim(handles.int); title(int2str(round(handles.phi)),'fontsize',12)
guidata(hObject, handles);


% --- Executes during object creation, after setting all properties.
function kappa_slider_CreateFcn(hObject, eventdata, handles)
% hObject    handle to kappa_slider (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: slider controls usually have a light gray background, change
%       'usewhitebg' to 0 to use default.  See ISPC and COMPUTER.
usewhitebg = 1;
if usewhitebg
    set(hObject,'BackgroundColor',[.9 .9 .9]);
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end


% --- Executes on slider movement.
function kappa_slider_Callback(hObject, eventdata, handles)
% hObject    handle to kappa_slider (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider
handles.kappa = get(hObject,'Value');
set(handles.kappa_edit,'String',num2str(handles.kappa));
handles.int = int_compute(handles.theta, handles.omega, handles.kappa, handles.ratio, handles.rho, handles.r, handles.p, ...
    handles.f00, handles.f01, handles.f02, handles.f10, handles.f11, handles.f12, handles.f20, handles.f21, handles.f22, handles.phi, get(handles.radiobutton1,'Value'));
mim(handles.int); title(int2str(round(handles.phi)),'fontsize',12)
guidata(hObject, handles);


% --- Executes during object creation, after setting all properties.
function ratio_slider_CreateFcn(hObject, eventdata, handles)
% hObject    handle to ratio_slider (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: slider controls usually have a light gray background, change
%       'usewhitebg' to 0 to use default.  See ISPC and COMPUTER.
usewhitebg = 1;
if usewhitebg
    set(hObject,'BackgroundColor',[.9 .9 .9]);
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end


% --- Executes on slider movement.
function ratio_slider_Callback(hObject, eventdata, handles)
% hObject    handle to ratio_slider (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider

handles.ratio = get(hObject,'Value');
set(handles.ratio_edit,'String',num2str(handles.ratio));
handles.int = int_compute(handles.theta, handles.omega, handles.kappa, handles.ratio, handles.rho, handles.r, handles.p, ...
    handles.f00, handles.f01, handles.f02, handles.f10, handles.f11, handles.f12, handles.f20, handles.f21, handles.f22, handles.phi, get(handles.radiobutton1,'Value'));
mim(handles.int); title(int2str(round(handles.phi)),'fontsize',12)
guidata(hObject, handles);


% --- Executes during object creation, after setting all properties.
function phi_slider_CreateFcn(hObject, eventdata, handles)
% hObject    handle to phi_slider (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: slider controls usually have a light gray background, change
%       'usewhitebg' to 0 to use default.  See ISPC and COMPUTER.
usewhitebg = 1;
if usewhitebg
    set(hObject,'BackgroundColor',[.9 .9 .9]);
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end


% --- Executes on slider movement.
function phi_slider_Callback(hObject, eventdata, handles)
% hObject    handle to phi_slider (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider
handles.phi = get(hObject,'Value');
handles.int = int_compute(handles.theta, handles.omega, handles.kappa, handles.ratio, handles.rho, handles.r, handles.p, ...
    handles.f00, handles.f01, handles.f02, handles.f10, handles.f11, handles.f12, handles.f20, handles.f21, handles.f22, handles.phi, get(handles.radiobutton1,'Value'));
mim(handles.int); title(int2str(round(handles.phi)),'fontsize',12)
guidata(hObject, handles);



% --- Executes during object creation, after setting all properties.
function radiobutton1_CreateFcn(hObject, eventdata, handles)
% hObject    handle to radiobutton1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called
set(hObject,'Value',1);


% --- Executes during object creation, after setting all properties.
function radiobutton2_CreateFcn(hObject, eventdata, handles)
% hObject    handle to radiobutton2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called
set(hObject,'Value',0);


% --- Executes on button press in radiobutton1.
function radiobutton1_Callback(hObject, eventdata, handles)
% hObject    handle to radiobutton1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of radiobutton1
set(handles.radiobutton2,'Value',0);
handles.int = int_compute(handles.theta, handles.omega, handles.kappa, handles.ratio, handles.rho, handles.r, handles.p, ...
    handles.f00, handles.f01, handles.f02, handles.f10, handles.f11, handles.f12, handles.f20, handles.f21, handles.f22, ...
    handles.phi, 1);
mim(handles.int); title(int2str(round(handles.phi)),'fontsize',12)
guidata(hObject, handles);


% --- Executes on button press in radiobutton2.
function radiobutton2_Callback(hObject, eventdata, handles)
% hObject    handle to radiobutton2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of radiobutton2
set(handles.radiobutton1,'Value',0);
handles.int = int_compute(handles.theta, handles.omega, handles.kappa, handles.ratio, handles.rho, handles.r, handles.p, ...
    handles.f00, handles.f01, handles.f02, handles.f10, handles.f11, handles.f12, handles.f20, handles.f21, handles.f22, ...
    handles.phi, 0);
mim(handles.int); title(int2str(round(handles.phi)),'fontsize',12)
guidata(hObject, handles);


% --- Executes during object creation, after setting all properties.
function imm_CreateFcn(hObject, eventdata, handles)
% hObject    handle to imm (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: place code in OpeningFcn to populate imm
return

function int = int_compute(theta, omega, kappa, ratio, rho, r, p, f00, f01, f02, f10, f11, f12, f20, f21, f22, phi, flag);

theta = theta/180*pi;
omega = omega/180*pi;
p = p - phi/180*pi;

if flag==1
    psi = 0;
    intc1(1,:) = 3*f00+3*f22+f11 - f00*cos(2*psi)*(1-cos(2*theta)+kappa*(3+cos(2*theta))*cos(2*omega)) + ...
        (f00-f11+f22)*(cos(2*theta)+2*kappa*cos(2*omega)*sin(theta)^2) + 4*kappa*f00*cos(theta)*sin(2*psi)*sin(2*omega);
    intc1(2,:) = -2*(2*f01-f21+2*f10-f12)*sin(theta)*(cos(psi)*cos(theta)*(1-kappa*cos(2*omega))+kappa*sin(psi)*sin(2*omega));
    intc1(3,:) = -3*f20+f11-3*f02 + (f20+f02)*cos(2*psi)*(kappa*(3+cos(2*theta))*cos(2*omega)+2*sin(theta)^2) - ...
        (f20+f11+f02)*(cos(2*theta)+2*kappa*cos(2*omega)*sin(theta)^2) - 4*(f20+f02)*kappa*cos(theta)*sin(2*psi)*sin(2*omega);
    intc1(4,:) = 2*(f21+f12)*sin(theta)*(cos(psi)*cos(theta)*(1-kappa*cos(2*omega))+kappa*sin(psi)*sin(2*omega));
    intc1(5,:) = -f22*(cos(2*psi)*(kappa*(3+cos(2*theta))*cos(2*omega)+2*sin(theta)^2)-4*kappa*cos(theta)*sin(2*psi)*sin(2*omega));
    ints1(1,:) = 2*(f21+f12)*sin(theta)*(cos(theta)*(1-kappa*cos(2*omega))*sin(psi)-kappa*cos(psi)*sin(2*omega));
    ints1(2,:) = (f20+f02)*(sin(2*psi)*(kappa*(3+cos(2*theta))*cos(2*omega)+2*sin(theta)^2)+4*kappa*cos(2*psi)*cos(theta)*sin(2*omega));
    ints1(3,:) = ints1(1,:);
    ints1(4,:) = -f22*(sin(2*psi)*(kappa*(3+cos(2*theta))*cos(2*omega)-2*sin(theta)^2)+4*kappa*cos(2*psi)*cos(theta)*sin(2*omega));
    
    psi = pi/2;
    intc2(1,:) = 3*f00+3*f22+f11 - f00*cos(2*psi)*(1-cos(2*theta)+kappa*(3+cos(2*theta))*cos(2*omega)) + ...
        (f00-f11+f22)*(cos(2*theta)+2*kappa*cos(2*omega)*sin(theta)^2) + 4*kappa*f00*cos(theta)*sin(2*psi)*sin(2*omega);
    intc2(2,:) = -2*(2*f01-f21+2*f10-f12)*sin(theta)*(cos(psi)*cos(theta)*(1-kappa*cos(2*omega))+kappa*sin(psi)*sin(2*omega));
    intc2(3,:) = -3*f20+f11-3*f02 + (f20+f02)*cos(2*psi)*(kappa*(3+cos(2*theta))*cos(2*omega)+2*sin(theta)^2) - ...
        (f20+f11+f02)*(cos(2*theta)+2*kappa*cos(2*omega)*sin(theta)^2) - 4*(f20+f02)*kappa*cos(theta)*sin(2*psi)*sin(2*omega);
    intc2(4,:) = 2*(f21+f12)*sin(theta)*(cos(psi)*cos(theta)*(1-kappa*cos(2*omega))+kappa*sin(psi)*sin(2*omega));
    intc2(5,:) = -f22*(cos(2*psi)*(kappa*(3+cos(2*theta))*cos(2*omega)+2*sin(theta)^2)-4*kappa*cos(theta)*sin(2*psi)*sin(2*omega));
    ints2(1,:) = 2*(f21+f12)*sin(theta)*(cos(theta)*(1-kappa*cos(2*omega))*sin(psi)-kappa*cos(psi)*sin(2*omega));
    ints2(2,:) = (f20+f02)*(sin(2*psi)*(kappa*(3+cos(2*theta))*cos(2*omega)+2*sin(theta)^2)+4*kappa*cos(2*psi)*cos(theta)*sin(2*omega));
    ints2(3,:) = ints2(1,:);
    ints2(4,:) = -f22*(sin(2*psi)*(kappa*(3+cos(2*theta))*cos(2*omega)-2*sin(theta)^2)+4*kappa*cos(2*psi)*cos(theta)*sin(2*omega));
    
    % int = interp1(rho,real(intc1(1,:)+intc2(1,:)),r,'cubic');
    int = interp1(rho,real(intc1(1,:)),r,'cubic') + interp1(rho,real(intc2(1,:)),r,'cubic');
    for jj=1:4 
        int = int + cos(jj*p).*interp1(rho,real(intc1(jj+1,:)),r,'cubic') + ...
            sin(jj*p).*interp1(rho,real(ints1(jj,:)),r,'cubic'); 
        int = int + cos(jj*(p+pi/2)).*interp1(rho,real(intc2(jj+1,:)),r,'cubic') + ...
            sin(jj*(p+pi/2)).*interp1(rho,real(ints2(jj,:)),r,'cubic');                 
    end
    
    intc1(1,:) = 4*(f00+f22)*sin(theta)^2+4*f11*cos(theta)^2;
    intc1(2,:) = -2*(f01-f21+f10-f12)*sin(2*theta);
    intc1(3,:) = -4*(f20+f02)*sin(theta)^2;
    
    tmp = interp1(rho,real(intc1(1,:)),r,'cubic');
    for jj=1:2 
        tmp = tmp + cos(jj*p).*interp1(rho,real(intc1(jj+1,:)),r,'cubic'); 
    end
else
    intc1(1,:) = (f00+f22)*(cos(omega)^2*cos(theta)^2 + sin(omega)^2) + f11*cos(omega)^2*sin(theta)^2;
    intc1(2,:) = -0.5*(f01-f21+f10-f12)*sin(2*theta)*cos(omega)^2;
    intc1(3,:) = 0.25*(f20+f02)*(1-3*cos(2*omega)-2*cos(omega)^2*cos(2*theta));
    ints1(1,:) = -0.5*(f01-f21+f10-f12)*sin(2*omega)*sin(theta);
    ints1(2,:) = -(f20+f02)*sin(2*omega)*cos(theta);
    
    omega = omega + pi/2;
    intc2(1,:) = (f00+f22)*(cos(omega)^2*cos(theta)^2 + sin(omega)^2) + f11*cos(omega)^2*sin(theta)^2;
    intc2(2,:) = -0.5*(f01-f21+f10-f12)*sin(2*theta)*cos(omega)^2;
    intc2(3,:) = 0.25*(f20+f02)*(1-3*cos(2*omega)-2*cos(omega)^2*cos(2*theta));
    ints2(1,:) = -0.5*(f01-f21+f10-f12)*sin(2*omega)*sin(theta);
    ints2(2,:) = -(f20+f02)*sin(2*omega)*cos(theta);

    int = (1-kappa)/2*interp1(rho,real(intc1(1,:)),r,'cubic') + (1+kappa)/2*interp1(rho,real(intc2(1,:)),r,'cubic');
    for jj=1:2
        int = int + (1-kappa)/2*cos(jj*p).*interp1(rho,real(intc1(jj+1,:)),r,'cubic') + ...
            (1-kappa)/2*sin(jj*p).*interp1(rho,real(ints1(jj,:)),r,'cubic'); 
        int = int + (1+kappa)/2*cos(jj*p).*interp1(rho,real(intc2(jj+1,:)),r,'cubic') + ...
            (1+kappa)/2*sin(jj*p).*interp1(rho,real(ints2(jj,:)),r,'cubic');                 
    end
    
    intc1(1,:) = (f00+f22)*sin(theta)^2+f11*cos(theta)^2;
    intc1(2,:) = -0.5*(f01-f21+f10-f12)*sin(2*theta);
    intc1(3,:) = -(f20+f02)*sin(theta)^2;
    
    tmp = interp1(rho,real(intc1(1,:)),r,'cubic');
    for jj=1:2 
        tmp = tmp + cos(jj*p).*interp1(rho,real(intc1(jj+1,:)),r,'cubic'); 
    end
    
end

int = (1-ratio)*int + ratio*tmp;
return

% --- Executes during object deletion, before destroying properties.
function pic_DeleteFcn(hObject, eventdata, handles)
% hObject    handle to pic (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)


% --- Executes during object deletion, before destroying properties.
function imm_DeleteFcn(hObject, eventdata, handles)
% hObject    handle to imm (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)


% --- Executes during object deletion, before destroying properties.
function figure1_DeleteFcn(hObject, eventdata, handles)
% hObject    handle to figure1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)


function [intx, inty, intz, rho, phi, fxx0, fxx2, fxz, byx0, byx2, byz, s0, dphase, v, pc, ps, psi] = SEPDipole(rho, z, NA, n1, n, n2, d1, d, d2, lambda, mag, focpos, atf, alpha)

% SEPDipole calculates the emission pattern of dipoles while imaged with a camera
%
% [intx, inty, intz, rho, phi, fxx0, fxx2, fxz, byx0, byx2, byzs0, dphase, v, pc, ps, psi] = SEPDipole(rho, z, NA, n1, n, n2, d1, d, d2, lambda, mag, focpos)
%
% input arguments:
% r        - two-element vector containing the minimum and maximum value of the rho coordinate
% z        - distance of molecule from surface
% NA       - numerical aperture
% n1       - refractive index vector (bottom to top) below molecule
% n        - refractive index of molecule's layer
% n2       - refractive index vector (bottom to top) above molecule
% d1       - thickness vector for layers below molecule
% d        - thickness of embedding layer
% d2       - thickness vector for layers above molecule
% lambda   - wavelength
% mag      - magnification
% focpos   - focus position along optical axis
%
% output arguments:
% int*     - intensity distributions
% f**      - electric field component distributions
% b**      - magnetic field component distributions
% rho, phi - 2d coordinate fields
% s0       - aperture variable
% dphase   - wave front deviation
% v        - angular distribution of radiation for vertical dipole
% pc, ps   - cosine and sine part of angular distribution of radiation for parallel dipole 
% psi      - emission angle for angular distribution of radiation 

maxnum = 1e3;
ni = 1.0; % it is assumed that imaging medium is air

ki = 2*pi/lambda*ni; 

if length(rho)==2
    drho = lambda/50;
    rho = rho(1) + (0:(rho(2)-rho(1))/drho)*drho;
elseif length(rho)>2
    drho = diff(rho(1:2));
end
rho = mag*rho; rho = rho(:)';
row = ones(size(rho));
    
chimin = 0;
chimax = asin(NA/mag/ni);
dchi = (chimax-chimin)/maxnum;
chi = chimin + (0.5:maxnum)*dchi;

ci = cos(chi);
si = sin(chi);
s0 = ni/n1(1)*mag*si;
c0 = sqrt(1-s0.^2);
psi = asin(s0);
[v,pc,ps] = DipoleL(psi,2*pi*z/lambda,n1,n,n2,2*pi*d1/lambda,2*pi*d/lambda,2*pi*d2/lambda);
v = v.'; ps = ps.'; pc = pc.';

if nargin>12 & ~isempty(atf)
    if length(atf)==1 % accounting for reflection losses when using water immersion; atf = ref index of cover slide
        [tmp, tms, tmp, tms] = Fresnel(n1(1)*c0,n1(1),atf);
        [mp, ms, mp, ms] = Fresnel(sqrt(atf^2-n1(1)^2*s0.^2),atf,n1(1));
    else % + aberration for water immersion; atf(2) = thickness mismatch
        [tmp, tms, tmp, tms] = Fresnel(n1(1)*c0,[n1(1) atf(1) atf(1)], 2*pi*atf(2)/lambda);        
        [mp, ms, mp, ms] = Fresnel(sqrt(atf(1)^2-n1(1)^2*s0.^2),[atf(1) n1(1) n1(1)], -2*pi*atf(2)/lambda);        
    end
    v = tmp.*mp.*v;
    pc = tmp.*mp.*pc;    
    ps = tms.*ms.*ps;
end

phase = -n1(1).*c0*focpos;
if nargin>13 & ~isempty(alpha)
    dphase = 0;
    for j=1:length(alpha)
        dphase = dphase + alpha(j)*(s0/max(s0)).^(2*j);
    end
    phase = phase + dphase;
    plot(s0/max(s0),dphase);
end
%alternative formulation but same result:
%phase = ni*mag*focpos*s0./c0.*si-n1(1)./c0*focpos; 
ez = exp(i*2*pi/lambda*phase);

fac = si.*sqrt(ci./c0); % aplanatic objective

barg = ki*si'*rho;
j0 = besselj(0,barg); j1 = besselj(1,barg); j2 = besselj(2,barg);

ezi = fac.*ci.*ez;
ezr = fac.*ez;

fxx0 = (ezi.*pc+ezr.*ps)*j0; % cos(0*phi)-component
fxx2 = -(ezi.*pc-ezr.*ps)*j2; % cos(2*phi)-component
fxz =  -2*i*(ezi.*v)*j1; % cos(1*phi)-component

byx0 = (ezr.*pc+ezi.*ps)*j0; % cos(0*phi)-component
byx2 = -(ezr.*pc-ezi.*ps)*j2; % cos(2*phi)-component
byz = -2*i*(ezr.*v)*j1; % cos(1*phi)-component

phi = 0:pi/100:2*pi; phi = phi'; col = ones(size(phi));
intx = real((col*fxx0+cos(2*phi)*fxx2).*conj(col*byx0+cos(2*phi)*byx2)); % int for x-dipole along x-pol
inty = real((sin(2*phi)*fxx2).*conj(sin(2*phi)*byx2)); % int for x-dipole along y-pol
intz = col*real(fxz.*conj(byz)); % int for z-dipole along x-pol

s0 = s0/max(s0);


function [v,pc,ps,tp,ts,tp1,ts1,fac] = DipoleL(theta,z,n1,n,n2,d1,d,d2)

% [v,pc,ps] = DipoleL(theta,z,n,d) calculates the electric field amplitudes of dipole radiation along emission angle theta 
% of a dipole at distance z from an interface within a layer 
% theta - direction of radiation downwards
% z  - molecule's distance from the bottom of its layer
% n1 - vector of refractive indices of the stack below the molecule's layer
% n  - refracive index of the molecule's layer
% n2 - vector of refractive indices of the stack above the the molecule's layer
% d1 - vector of layer thickness values of the stack below the molecule's layer ( length(d1)=length(n1)-1 )
% d  - thickness of molecule's layer
% d2 - vector of layer thickness values of the stack above the molecule's layer ( length(d2)=length(n2)-1 )

z = z(:)';
col = ones(size(z));
theta = abs(theta(:));
ind = theta<pi/2;

v = zeros(length(theta),length(z));
pc = v; ps = v;

if sum(ind)>0
    tmp = theta(ind);
    w = sqrt(n^2-n1(1)^2*sin(tmp).^2);
    [rpu, rsu, tpu, tsu] = Fresnel(w,[n n2],d2);
    [rpd, rsd, tpd, tsd] = Fresnel(w,[n n1(end:-1:1)],d1(end:-1:1));    
    tp = (tpd./(1-rpu.*rpd.*exp(2*i*w.'*d))).';
    ts = (tsd./(1-rsu.*rsd.*exp(2*i*w.'*d))).';
    tp1 = tp.*(rpu.*exp(2*i*w.'*d)).';
    ts1 = ts.*(rsu.*exp(2*i*w.'*d)).';
    fac = sqrt(n1(1))*n1(1)*cos(tmp)./w;
    ez = exp(i*w*z);
    v(ind,:) = ((fac.*n1(1)/n.*sin(tmp).*tp)*col).*ez + ((fac.*n1(1)/n.*sin(tmp).*tp1)*col)./ez;
    pc(ind,:) = ((fac.*w/n.*tp)*col).*ez - ((fac.*w/n.*tp1)*col)./ez;
    ps(ind,:) = ((fac.*ts)*col).*ez + ((fac.*ts1)*col)./ez;
end


function [rp,rs,tp,ts] = Fresnel(w1,n1,n2);

maxval = exp(100);
minval = exp(-100);

w1 = w1(:).'; n1 = n1(:).'; n2 = n2(:).';

if length(n1)==1 & length(n2)==1
    w2 = sqrt(n2^2-n1^2+w1.^2);
    w2(imag(w2)<0) = conj(w2(imag(w2)<0));
    rp = (w1*n2^2-w2*n1^2)./(w1*n2^2+w2*n1^2);
    rs = (w1-w2)./(w1+w2);
    tp = 2*n1*n2*w1./(w1*n2^2+w2*n1^2);
    ts = 2*w1./(w1+w2);
elseif length(n1)==length(n2)+2
    if length(n1)==2
        [rp,rs,tp,ts] = Fresnel(w1,n1(1),n1(2));
    else
        n = n1; d = [0 n2 0];
        w(1,:) = w1;
        for j = 2:length(n)    
            w(j,:) = sqrt(n(j)^2-n(1)^2+w1.^2);
            w(j,imag(w(j,:))<0) = conj(w(j,imag(w(j,:))<0));        
        end
        
        j = length(n);
        M11 = (w(j,:)./w(j-1,:)*n(j-1)/n(j)+n(j)/n(j-1))/2;
        M12 = (-w(j,:)./w(j-1,:)*n(j-1)/n(j)+n(j)/n(j-1))/2;
        M21 = M12;
        M22 = M11;
        for j=length(n)-1:-1:2
            M11 = exp(-i*w(j,:)*d(j)).*M11;
            M12 = exp(-i*w(j,:)*d(j)).*M12;
            M21 = exp(i*w(j,:)*d(j)).*M21;        
            M22 = exp(i*w(j,:)*d(j)).*M22;                
            N11 = (w(j,:)./w(j-1,:)*n(j-1)/n(j)+n(j)/n(j-1))/2;
            N12 = (-w(j,:)./w(j-1,:)*n(j-1)/n(j)+n(j)/n(j-1))/2;
            N21 = N12;
            N22 = N11;
            tmp11 = N11.*M11+N12.*M21;
            tmp12 = N11.*M12+N12.*M22;
            tmp21 = N21.*M11+N22.*M21;
            tmp22 = N21.*M12+N22.*M22;
            M11 = tmp11;
            M12 = tmp12;
            M21 = tmp21;
            M22 = tmp22;
        end
        
        rp = M21./M11;
        tp = 1./M11;
        
        j = length(n);
        M11 = (w(j,:)./w(j-1,:)+1)/2;
        M12 = (-w(j,:)./w(j-1,:)+1)/2;
        M21 = M12;
        M22 = M11;
        for j=length(n)-1:-1:2
            M11 = exp(-i*w(j,:)*d(j)).*M11;
            M12 = exp(-i*w(j,:)*d(j)).*M12;        
            M21 = exp(i*w(j,:)*d(j)).*M21;        
            M22 = exp(i*w(j,:)*d(j)).*M22;                
            N11 = (w(j,:)./w(j-1,:)+1)/2;
            N12 = (-w(j,:)./w(j-1,:)+1)/2;
            N21 = N12;
            N22 = N11;
            tmp11 = N11.*M11+N12.*M21;
            tmp12 = N11.*M12+N12.*M22;
            tmp21 = N21.*M11+N22.*M21;
            tmp22 = N21.*M12+N22.*M22;
            M11 = tmp11;
            M12 = tmp12;
            M21 = tmp21;
            M22 = tmp22;
        end
        rs = M21./M11;
        ts = 1./M11;
    end
end


function mim(x,p1,p2)

if nargin==1
    nd = ndims(x);
    switch nd
        case 2
            imagesc(x);
            axis image
            colormap hot
            axis off
        case 3
            nj = size(x,3);
            for j=1:nj
                subplot(1,nj,j);
                mim(x(:,:,j))
            end
        case 4
            nj = size(x,3);
            nk = size(x,4);
            for j=1:nj
                for k=1:nk
                    subplot(nk,nj,j+(k-1)*nj);
                    mim(x(:,:,j,k))
                end
            end
    end
end
if nargin==2
    if prod(size(p1))==2
        imagesc(x,p1);
        axis image
        colormap hot
        axis off
    elseif p1=='h'
        mim(x);
        colorbar('h');
    elseif p1=='v'
        mim(x);
        colorbar;
    elseif isstr(p1)
        mim(x); 
        [a,b]=size(x);
        text(b*(1-0.025*length(p1)),0.06*a,p1,'FontName','Times','FontSize',16,'Color','w');
    else
        c = zeros(size(x,1),size(x,2),3);
        mm = min(p1(isfinite(p1)));
        tmp = max(p1(isfinite(p1)))-mm;
        if tmp>0
            p1 = (p1-mm)/tmp;
        else
            p1 = zeros(size(p1));
        end
        mm = min(x(isfinite(p1)));
        tmp = max(x(isfinite(p1)))-mm;
        if tmp>0
            x = 1 + (x-mm)/tmp*63;
        else
            x = 64*ones(size(x));
        end
        x(isnan(x)) = 1;
        x(x<1) = 1;
        x(x>64) = 64;
        col = double(jet);
        for j=1:3
            c(:,:,j) = ((floor(x)+1-x).*reshape(col(floor(x),j),size(x)) + (x-floor(x)).*reshape(col(ceil(x),j),size(x))).*p1;
        end
        if tmp>0
            subplot('position',[0.1 0.1 0.7 0.8]);
            imagesc(c);
            axis image
            axis off
            subplot('position',[0.8 0.1 0.1 0.8]);
            pcolor([0 0.1]*tmp,mm+(1:size(x,1))/size(x,1)*tmp,(mm+(1:size(x,1))'/size(x,1)*tmp)*ones(1,2))
            shading flat
            axis image
            set(gca,'xtick',[],'yaxislocation','right')
            colormap(jet)            
        else
            imagesc(c);
            axis image
            axis off
        end
    end
end
if nargin==3
    if size(x,1)==size(p1,1) & size(x,2)==size(p1,2)
        if prod(size(p2))==2
            c = zeros(size(x,1),size(x,2),3);
            mm = min(p1(isfinite(p1)));
            tmp = max(p1(isfinite(p1)))-mm;
            if tmp>0
                p1 = (p1-mm)/tmp;
            else
                p1 = zeros(size(p1));
            end
            mm = p2(1);
            tmp = p2(2) - mm;
            if tmp>0
                x = 1 + (x-mm)/tmp*63;
            else
                x = 64*ones(size(x));
            end
            x(isnan(x)) = 1;
            x(x<1) = 1;
            x(x>64) = 64;
            col = double(jet);
            for j=1:3
                c(:,:,j) = ((floor(x)+1-x).*reshape(col(floor(x),j),size(x)) + (x-floor(x)).*reshape(col(ceil(x),j),size(x))).*p1;
            end
            if tmp>0
                subplot('position',[0.1 0.1 0.7 0.8]);
                imagesc(c);
                axis image
                axis off
                subplot('position',[0.8 0.1 0.1 0.8]);
                pcolor([0 0.1]*tmp,mm+(1:size(x,1))/size(x,1)*tmp,(mm+(1:size(x,1))'/size(x,1)*tmp)*ones(1,2))
                shading flat
                axis image
                set(gca,'xtick',[],'yaxislocation','right')
                colormap(jet)
            else
                imagesc(c);
                axis image
                axis off
            end
        end
    elseif prod(size(p1))==2
        mim(x,p1);
        if p2=='h'
            colorbar('h');
        elseif p2=='v'
            colorbar;
        end
    elseif prod(size(p2))==2
        mim(x,p2);
        if p1=='h'
            colorbar('h');
        elseif p1=='v'
            colorbar;
        end
    end
end
figure(gcf);