function varargout = Mays_Stretching_Folding_GUI(varargin)

% MAYS_STRETCHING_FOLDING_GUI MATLAB code for Mays_Stretching_Folding_GUI.fig
%
%      MAYS_STRETCHING_FOLDING_GUI, by itself, creates a new
%      MAYS_STRETCHING_FOLDING_GUI or raises the existing singleton*.
%
%      H = MAYS_STRETCHING_FOLDING_GUI returns the handle to a new
%      MAYS_STRETCHING_FOLDING_GUI or the handle to the existing
%      singleton*.
%
%      MAYS_STRETCHING_FOLDING_GUI('CALLBACK',hObject,eventData,handles,...)
%      calls the local function named CALLBACK in
%      MAYS_STRETCHING_FOLDING_GUI.M with the given input arguments.
%
%      MAYS_STRETCHING_FOLDING_GUI('Property','Value',...) creates a new
%      MAYS_STRETCHING_FOLDING_GUI or raises the existing singleton*.
%      Starting from the left, property value pairs are applied to the GUI
%      before Mays_Stretching_Folding_GUI_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property
%      application stop.  All inputs are passed to
%      Mays_Stretching_Folding_GUI_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
%
% Version 1.0, dcm 3/28/2015, Matlab 8.4.0.150421 (R2014b)
% Version 1.1, dcm 10/9/2015, Matlab 8.4.0.150421 (R2014b)

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

% Last Modified by GUIDE v2.5 11-Oct-2015 21:08:36

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
	'gui_Singleton',  gui_Singleton, ...
	'gui_OpeningFcn', @Mays_Stretching_Folding_GUI_OpeningFcn, ...
	'gui_OutputFcn',  @Mays_Stretching_Folding_GUI_OutputFcn, ...
	'gui_LayoutFcn',  [] , ...
	'gui_Callback',   []);
if nargin && ischar(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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% A

% --- Executes just before Mays_Stretching_Folding_GUI is made visible.
function Mays_Stretching_Folding_GUI_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 Mays_Stretching_Folding_GUI (see VARARGIN)

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

% Update handles structure
guidata(hObject, handles);

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% B

% --- Outputs from this function are returned to the command line.
function varargout = Mays_Stretching_Folding_GUI_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;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1.1

function Input_N_Wells_Callback(hObject, eventdata, handles)
% hObject    handle to Input_N_Wells (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 Input_N_Wells as text
%        str2double(get(hObject,'String')) returns contents of Input_N_Wells as a double

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1.2

% --- Executes during object creation, after setting all properties.
function Input_N_Wells_CreateFcn(hObject, eventdata, handles)
% hObject    handle to Input_N_Wells (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 && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
	set(hObject,'BackgroundColor','white');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 2

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

cla reset

% Get user data from GUI
N_Wells = str2double(get(handles.Input_N_Wells,'String'));
r = str2double(get(handles.Input_Plume_Radius,'String'));

% (from GUTS.M) calculate coordinates of wells on unit circle
points = 2*floor(N_Wells/2);
theta_wells = (2*pi/360)*linspace(0,360-360/points,points);   % radians
r_wells = 1;
x_wells = -r_wells*cos(theta_wells);
y_wells = r_wells*sin(theta_wells);
if mod(N_Wells,2)
	x_wells = [x_wells 0];
	y_wells = [y_wells 0];
end
plot(handles.Display_Axes,x_wells,y_wells,'ko','MarkerFaceColor','white')
for i = 1:points
	text(1.25*x_wells(i),1.25*y_wells(i),num2str(i),'HorizontalAlignment','center','BackgroundColor','w')
end
if mod(N_Wells,2)
	text(-0.25,0,num2str(N_Wells),'HorizontalAlignment','center','BackgroundColor','w')
end
hold on
grid on
axis equal
axis([-2 2 -2 2])
xlabel('Position West-East','fontsize',9)
ylabel('Position South-North','fontsize',9)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 3.1

function Input_Plume_Radius_Callback(hObject, eventdata, handles)
% hObject    handle to Input_Plume_Radius (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 Input_Plume_Radius as text
%        str2double(get(hObject,'String')) returns contents of Input_Plume_Radius as a double

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 3.2

% --- Executes during object creation, after setting all properties.
function Input_Plume_Radius_CreateFcn(hObject, eventdata, handles)
% hObject    handle to Input_Plume_Radius (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 && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
	set(hObject,'BackgroundColor','white');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 4

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

cla reset

% Get user data from GUI
N_Wells = str2double(get(handles.Input_N_Wells,'String'));
r = str2double(get(handles.Input_Plume_Radius,'String'));

% (from GUTS.M) calculate coordinates of wells on unit circle
points = 2*floor(N_Wells/2);
theta_wells = (2*pi/360)*linspace(0,360-360/points,points);   % radians
r_wells = 1;
x_wells = -r_wells*cos(theta_wells);
y_wells = r_wells*sin(theta_wells);
if mod(N_Wells,2)
	x_wells = [x_wells 0];
	y_wells = [y_wells 0];
end

% calculate coordinates of initial plume
theta = (2*pi/360)*linspace(0,360,3601);
x = -r*cos(theta);
y = r*sin(theta);

% plot initial circle with wells
plot(handles.Display_Axes,x_wells,y_wells,'ko','MarkerFaceColor','white')
for i = 1:points
	text(1.25*x_wells(i),1.25*y_wells(i),num2str(i),'HorizontalAlignment','center','BackgroundColor','w')
end
if mod(N_Wells,2)
	text(-0.25,0,num2str(N_Wells),'HorizontalAlignment','center','BackgroundColor','w')
end
hold on
fill(x,y,'b');
grid on
axis equal
axis([-2 2 -2 2])
xlabel('Position West-East','fontsize',9)
ylabel('Position South-North','fontsize',9)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 5.1

function Input_Delay_Callback(hObject, eventdata, handles)
% hObject    handle to Input_Delay (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 Input_Delay as text
%        str2double(get(hObject,'String')) returns contents of Input_Delay as a double

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 5.2

% --- Executes during object creation, after setting all properties.
function Input_Delay_CreateFcn(hObject, eventdata, handles)
% hObject    handle to Input_Delay (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 && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
	set(hObject,'BackgroundColor','white');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 6

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

% clear previous results
delete working_previous.mat

% Get user data from GUI
N_Wells = str2double(get(handles.Input_N_Wells,'String'));

% Pass inputs to subfunction Pump_Scheme
Input_Table_Contents = get(handles.Input_Table);
Input_Table_Data = Input_Table_Contents.Data;
wells = str2double(Input_Table_Data(:,1));		% wells entered as strings
index = find(wells>0)';		% must be row vector to execute one-at-a-time
current = index;
if max(wells) > N_Wells
    set(handles.Display_Message,'String','Error! At least one well number is too large. Please correct this and try again.','ForegroundColor',[1 0 0])
else
    Pump_Scheme(hObject,eventdata,handles,current);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 7.1

function Input_Current_Step_Callback(hObject, eventdata, handles)
% hObject    handle to Input_Current_Step (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 Input_Current_Step as text
%        str2double(get(hObject,'String')) returns contents of Input_Current_Step as a double

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 7.2

% --- Executes during object creation, after setting all properties.
function Input_Current_Step_CreateFcn(hObject, eventdata, handles)
% hObject    handle to Input_Current_Step (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 && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
	set(hObject,'BackgroundColor','white');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 8

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

% Get user data from GUI
current = str2double(get(handles.Input_Current_Step,'String'));
Input_Table_Contents = get(handles.Input_Table);
Input_Table_Data = Input_Table_Contents.Data;
wells = str2double(Input_Table_Data(:,1));		% wells entered as strings

% Run single step of current simulation
if current > length(find(isfinite(wells)))
    set(handles.Display_Message,'String','Error! Current step too large. Please correct and try again.','ForegroundColor',[1 0 0])
else
    Pump_Scheme(hObject,eventdata,handles,current)
end

% Update display after running current step
set(handles.Input_Current_Step,'String',num2str(current+1))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 9

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

% set(handles.Input_Table,'Data',[])
% set(handles.Input_Table,'Data',[])

% code from http://www.mathworks.com/matlabcentral/answers/22845-reset-uitable-data
% set(handles.Input_Table, 'Data', cell(size(get(handles.Input_Table,'Data'))));
% newd = get(handles.Input_Table,'Data');
% charlocs = cellfun(@ischar, newd);
% newd(charlocs) = {''};
% newd(~charlocs) = {[]};
% set(handles.Input_Table, 'Data', newd);

% code from http://www.mathworks.com/matlabcentral/newsreader/view_thread/279207
close(gcbf)
delete working_previous.mat
Mays_Stretching_Folding_GUI

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 10

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

% get user data from GUI
N_Wells = str2double(get(handles.Input_N_Wells,'String'));
current = str2double(get(handles.Input_Current_Step,'String')) - 1;

% load results from previous step
if strcmp(handles.Display_Message.String,'Simulation reset to position before current step shown above.') || ...
        strcmp(handles.Display_Message.String,'Undo no longer available. Please restart from step 1.')
    set(handles.Display_Message,'String','Undo no longer available. Please restart from step 1.','ForegroundColor','k')
elseif current == 1
    set(handles.Display_Message,'String','Undo not available for step 1. Instead, please restart with step 1.','ForegroundColor','k')
elseif exist('working_previous.mat') == 0
    set(handles.Display_Message,'String','Undo only available with Run Current Step. Please use that button starting with step 1.','ForegroundColor','k')
elseif current > 1
    cla reset
    load working_previous.mat
    
    % (from GUTS.M) calculate coordinates of wells on unit circle
    points = 2*floor(N_Wells/2);
    theta_wells = (2*pi/360)*linspace(0,360-360/points,points);   % radians
    r_wells = 1;
    x_wells = -r_wells*cos(theta_wells);
    y_wells = r_wells*sin(theta_wells);
    if mod(N_Wells,2)
        x_wells = [x_wells 0];
        y_wells = [y_wells 0];
    end
    
    % (from GUTS.M) plot wells
    plot(handles.Display_Axes,x_wells,y_wells,'ko','MarkerFaceColor','white')
    for i = 1:points
        text(1.25*x_wells(i),1.25*y_wells(i),num2str(i),'HorizontalAlignment','center','BackgroundColor','w')
    end
    if mod(N_Wells,2)
        text(-0.25,0,num2str(N_Wells),'HorizontalAlignment','center','BackgroundColor','w')
    end
    hold on
    grid on
    axis equal
    axis([-2 2 -2 2])
    xlabel('Position West-East','fontsize',9)
    ylabel('Position South-North','fontsize',9)
    PlotColor = 'b';
    
    % plot new positions (red if ever extracted)
    delete(hPlotData)
    hPlotData = fill(x,y,PlotColor);
    hTextStep = text(-1.9,1.8,['\bf{after step ' num2str(current-1) '}'],'fontsize',9);
    hTextInterfaceL = text(-1.9,1.6,['circumference = ' num2str(InterfaceL,'%.1f')],'fontsize',9);
    hTextOuter = text(-1.9,1.4,['distance = ' num2str(MaxOuter,'%.2f')],'fontsize',9);
    hTextSumV = text(-1.9,1.2,['volume = ' num2str(SumV,'%.1f')],'fontsize',9);      % could be +/-
    hTextScore = text(-1.9,1.0,['\bf{' num2str(score,'%.1f') ' points}'],'fontsize',9);
    
    % Update display after running current step
    set(handles.Input_Current_Step,'String',num2str(current))
    set(handles.Display_Message,'String','Simulation reset to position before current step shown above.','ForegroundColor','k')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 11.1

function Input_Your_Name_Callback(hObject, eventdata, handles)
% hObject    handle to Input_Your_Name (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 Input_Your_Name as text
%        str2double(get(hObject,'String')) returns contents of Input_Your_Name as a double

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 11.2

% --- Executes during object creation, after setting all properties.
function Input_Your_Name_CreateFcn(hObject, eventdata, handles)
% hObject    handle to Input_Your_Name (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 && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
	set(hObject,'BackgroundColor','white');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 12

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

% import data from last step of simulation
load working.mat
Output_Name = get(handles.Input_Your_Name,'String');
Output_Name = strrep(Output_Name, ' ','_');
Output_Time = datestr(now,31);
Output_Time = strrep(Output_Time, ' ','_');
Output_Time = strrep(Output_Time, '-','_');
Output_Time = strrep(Output_Time, ':','_');
Output_N_Wells = get(handles.Input_N_Wells,'String');
Output_Circumference = num2str(InterfaceL);
Output_Outer = num2str(MaxOuter);
Output_SumV = num2str(SumV);
Output_SumAbsV = num2str(SumAbsV);
Output_Score = num2str(score);
Output_Wells = num2cell(wells);
Output_Inject = num2cell(inject);
Output_Volume = num2cell(volume);
Output_Dipole = num2cell(dipole);
Output_Partner = num2cell(partner);
C = {'Name',Output_Name,'','',''; ...
	'N_Wells',Output_N_Wells,'','',''; ...
	'Circumference',Output_Circumference,'','',''; ...
	'Outer',Output_Outer,'','',''; ...
	'SumV',Output_SumV,'','',''; ...
	'SumAbsV',Output_SumAbsV,'','',''; ...
	'Score',Output_Score,'','',''; ...
	'Main','Action','Volume','Dipole','Partner'};
D = [C; Output_Wells Output_Inject Output_Volume Output_Dipole Output_Partner];
T = cell2table(D);
FileName = [Output_Name '_RUN_' Output_Time '.txt'];
writetable(T,FileName,'Delimiter','\t','WriteVariableNames',0)
set(handles.Display_Message,'String',['Success! Your pumping scheme has been saved to ' FileName '. Thanks for playing!'],'ForegroundColor',[0 0 0])

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 13.1

function Display_Message_Callback(hObject, eventdata, handles)
% hObject    handle to Display_Message (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 Display_Message as text
%        str2double(get(hObject,'String')) returns contents of Display_Message as a double

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 13.2

% --- Executes during object creation, after setting all properties.
function Display_Message_CreateFcn(hObject, eventdata, handles)
% hObject    handle to Display_Message (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 && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
	set(hObject,'BackgroundColor','white');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CUSTOM-1

function Pump_Scheme(hObject,eventdata,handles,current)
%
% Function to advance plume using steps specified by "current". To run the
% whole sequence, "current" is a vector of all valid main pump numbers. To
% run only the current step, "current" is the value of the current step.
% This function takes the following inputs (first three c/o GUIDE):
%
% hObject    handle to Button_Run_Sequence (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% current	vector specifying which steps to execute

% if starting with step 1, whole sequence or step-by-step, then start over
if min(current) == 1
	cla reset
	score = 0;
	OK = 1;
% if not starting with step 1, assume step-by-step, load previous
else
	cla reset
	load working.mat
    % save previous results in case user wants to undo
    save('working_previous.mat','x','y','np','Lo','OK','hPlotData', ...	% for next step
        'InterfaceL','MaxOuter','SumV','SumAbsV','score', ...	% for header
        'wells','inject','volume','dipole','partner')			% for table
end

% display running message
set(handles.Display_Message,'String','(running)','ForegroundColor','k')
pause(0.1)

% Get user data from GUI
N_Wells = str2double(get(handles.Input_N_Wells,'String'));
Delay = str2double(get(handles.Input_Delay,'String'));
Input_Table_Contents = get(handles.Input_Table);
Input_Table_Data = Input_Table_Contents.Data;
wells = str2double(Input_Table_Data(:,1));		% wells entered as strings
action = Input_Table_Data(:,2);					% each cell has string
pump = Input_Table_Data(:,3);					% each cell has matrix,
index = cellfun(@isempty,pump);					% so find empty cells,
pump(index) = num2cell(0);						% convert to zero, and
pump = cell2mat(pump);							% convert to double array
dipole = cell2mat(Input_Table_Data(:,4));		% logical + # = double
partner = str2double(Input_Table_Data(:,5));	% wells entered as strings

% Initialize other variables
inject = zeros(length(wells),1);
volume = zeros(length(wells),1);

% Loop through steps to initialize "inject" and "volume"
for i = 1:max(current)
	if strcmp(action(i),'inject')
		inject(i) = 1;
	else
		inject(i) = 0;
	end
	if inject(i)
		volume(i) = pump(i);
	else
		volume(i) = -pump(i);
	end
end

% (from GUTS.M) calculate coordinates of wells on unit circle
points = 2*floor(N_Wells/2);
theta_wells = (2*pi/360)*linspace(0,360-360/points,points);   % radians
r_wells = 1;
x_wells = -r_wells*cos(theta_wells);
y_wells = r_wells*sin(theta_wells);
if mod(N_Wells,2)
	x_wells = [x_wells 0];
	y_wells = [y_wells 0];
end

% (from GUTS.M) calculate coordinates of initial plume r=1/4 (per WRR-1)
r = str2double(get(handles.Input_Plume_Radius,'String'));
if min(current) == 1
	theta = (2*pi/360)*linspace(0,360,3601);    % hard-coded 3601 particles
	x = -r*cos(theta);
	y = r*sin(theta);
end
np = length(x);
SumV = 0;                       % no pumping yet (we are at the IC)
SumAbsV = 0;				
D = separation(x,y);
E = D(2:end,1:end-1);
InterfaceL = trace(E);
if min(current) == 1
	Lo = InterfaceL;	% only assign on first pass (before pumping)
end
outer = r;
MaxOuter = max(outer);

% (from GUTS.M) plot initial position of wells and plume
plot(handles.Display_Axes,x_wells,y_wells,'ko','MarkerFaceColor','white')
for i = 1:points
	text(1.25*x_wells(i),1.25*y_wells(i),num2str(i),'HorizontalAlignment','center','BackgroundColor','w')
end
if mod(N_Wells,2)
	text(-0.25,0,num2str(N_Wells),'HorizontalAlignment','center','BackgroundColor','w')
end
hold on
grid on
axis equal
axis([-2 2 -2 2])
xlabel('Position West-East','fontsize',9)
ylabel('Position South-North','fontsize',9)
PlotColor = 'b';
if min(current) == 1
	hPlotData = fill(x,y,PlotColor);
	hTextStep = text(-1.9,1.8,'\bf{initial circle}','fontsize',9);
	hTextInterfaceL = text(-1.9,1.6,['circumference = ' num2str(InterfaceL,'%.1f')],'fontsize',9);
	hTextOuter = text(-1.9,1.4,['distance = ' num2str(MaxOuter,'%.2f')],'fontsize',9);
	hTextSumV = text(-1.9,1.2,['volume = ' num2str(SumV,'%.1f')],'fontsize',9);    % could be +/-
	hTextScore = text(-1.9,1.0,'\bf{zero points}','fontsize',9);
	pause(Delay)
end

% (from GUTS.M) implement user-specified pumping scheme
for i = current
	x_previous = x;
    y_previous = y;
    clear x_extract
	clear y_extract
	if dipole(i) == 0
		% not dipole, so move plume with PUSH.M
        [x,y] = push(x,y,x_wells(wells(i)),y_wells(wells(i)),volume(i));
		if volume(i) < 0
			x_extract = x_wells(wells(i));
			y_extract = y_wells(wells(i));
		end
	elseif dipole(i) == 1
		% dipole, so move plume with DIPOLEMAP.M
		if wells(i) == partner(i)
			fprintf('Skipping step %.0f, because main and partner wells must be different.\n',i)
			set(handles.Display_Message,'String',['Skipping step ' num2str(i) ', because main and partner wells must be different.'],'ForegroundColor',[1 0 0])
		elseif isnan(partner(i)) || partner(i)<1 || partner(i)>N_Wells
			fprintf('Skipping step %.0f, because partner well not valid.\n',i)
			set(handles.Display_Message,'String',['Skipping step ' num2str(i) ', because partner well not valid.'],'ForegroundColor',[1 0 0])
		elseif volume(i) > 0
            [x,y] = dipolemap(x,y,[x_wells(wells(i)),y_wells(wells(i))],[x_wells(partner(i)),y_wells(partner(i))],volume(i),1,0);
			x_extract = x_wells(partner(i));
			y_extract = y_wells(partner(i));
		elseif volume(i) < 0
			% transpose "partner" and "wells" and use -volume(i)
            [x,y] = dipolemap(x,y,[x_wells(partner(i)),y_wells(partner(i))],[x_wells(wells(i)),y_wells(wells(i))],-volume(i),1,0);
            x_extract = x_wells(wells(i));
			y_extract = y_wells(wells(i));
		end
		np = length(x);
	end
	
	% check and trim in case of active extraction well
	if exist('x_extract','var')
		x_check = [x_extract x];
		y_check = [y_extract y];
		D_check = separation(x_check,y_check);
		check = D_check(1,2:end);
		index = find(check > 1e-6);
		if length(index) < np           % this means extraction!
			PlotColor = 'r';
			OK = 0;
			score = 0;
			x = x(index);
			y = y(index);
			np = length(x);
            percent_remaining = 100*np/3601;
			if np > 0
				fprintf('Removed extracted particle(s). %.0f%% of initial particles remaining.\n',percent_remaining);
				set(handles.Display_Message,'String','Error! Part of the plume was removed. Zero score. Try again!','ForegroundColor',[1 0 0])
			else
				fprintf('Sorry, all particles have now been extracted. GAME OVER!\n')
				set(handles.Display_Message,'String','Sorry, all particles have now been extracted. Zero score. GAME OVER','ForegroundColor',[1 0 0])
				break
			end
		end
	end
	
	% check for out of bounds plume
	if max(abs([x y])) > 2
		if OK
			PlotColor = 'g';
		end
		OK = 0;
		score = 0;
		fprintf('Plume out of bounds.\n');
		set(handles.Display_Message,'String','Error! Plume out of bounds. Zero score. Try again!','ForegroundColor',[0 1 0])
	end
	
	% calculate metrics for new plume
    SumV = sum( (1-dipole(1:i)).*(volume(1:i)));              % count if dipole = 0, but not if dipole = 1
	SumAbsV = sum( (dipole(1:i)+1 ).*abs(volume(1:i)));		% if dipole = 1, count ABSOLUTE pumping volume twice
	D = separation(x,y);
	E = D(2:end,1:end-1);
	InterfaceL = trace(E);
	x_outer = [mean(x) x];
	y_outer = [mean(y) y];
	D_outer = separation(x_outer,y_outer);				% distance from plume centroid
	outer = D_outer(1,2:end);
	MaxOuter = max(outer);
	A = pi*max(outer)^2;
	if OK && SumAbsV>0
		score = 1000*(InterfaceL - Lo)^2/SumAbsV/A;
	end
	
	% plot new positions (red if ever extracted)
	delete(hPlotData)
	if min(current) == 1
		delete(hTextStep)
		delete(hTextInterfaceL)
		delete(hTextOuter)
		delete(hTextSumV)
		delete(hTextScore)
	end
	hPlotData = fill(x,y,PlotColor);
	hTextStep = text(-1.9,1.8,['\bf{after step ' num2str(i) '}'],'fontsize',9);
	hTextInterfaceL = text(-1.9,1.6,['circumference = ' num2str(InterfaceL,'%.1f')],'fontsize',9);
	hTextOuter = text(-1.9,1.4,['distance = ' num2str(MaxOuter,'%.2f')],'fontsize',9);
	hTextSumV = text(-1.9,1.2,['volume = ' num2str(SumV,'%.1f')],'fontsize',9);      % could be +/-
	hTextScore = text(-1.9,1.0,['\bf{' num2str(score,'%.1f') ' points}'],'fontsize',9);
	pause(Delay)
	save('working.mat','x','y','np','Lo','OK','hPlotData', ...	% for next step
		'InterfaceL','MaxOuter','SumV','SumAbsV','score', ...	% for header
		'wells','inject','volume','dipole','partner')			% for table
end

% remove running message (unless previously removed by OK = 0)
if OK
    set(handles.Display_Message,'String','','ForegroundColor','k')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CUSTOM-2

function [x_out,y_out] = push(x_in,y_in,x_well,y_well,V)

% function [x_out,y_out] = push(x_in,y_in,x_well,y_well,V)
%
% This function maps the material line specified by vectors (x_in,y_in) to
% its final position (x_out,y_out) resulting from the well at
% (x_well,y_well) injecting at Q for duration T. Assumes a homogeneous 2D
% infinite aquifer with a fully penetrating well. If particles are
% extracted, they are assumed to remain at the extraction well.
%
% Inputs:
%   x_in, vector of initial particle x-positions
%   y_in, vector of initial particle y-positions
%   x_well, scalar well x-position
%   y_well, scalar well y-position
%   V, injection volume per aquifer thickness (V<0 for extraction)
%
% Output:
%   x_out, vector of final particle x-positions
%   y_out, vector of final particle y-positions
%
% Rev. 0, dcm 10/19/2009, Matlab 7.7.0.471 (R2008b)
% Rev. 1, dcm 10/28/2009, Matlab 7.7.0.471 (R2008b)

% convert to local coordinates with well at (0,0)
x = x_in - x_well;
y = y_in - y_well;

% convert to radial coordinates
r = sqrt(x.^2 + y.^2);
theta = atan2(y,x);

% compute new radius, where T = (t - t_o)
r = real(sqrt(V/pi + r.^2));

% convert to Cartesian coordinates
x = r.*cos(theta);
y = r.*sin(theta);

% convert to global coordinates with well at (x_well,y_well)
x_out = x + x_well;
y_out = y + y_well;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CUSTOM-3

function [x_out,y_out] = dipolemap(x,y,xy_in,xy_out,Q,dt,animate)

% function [x_out,y_out] = dipolemap(x,y,xy_in,xy_out,Q,dt,animate)
%
% This function turns the analytical solution for the pulsed dipole given
% in the appendix of Jones and Aref (1988) into a map, which moves the
% plume forward dt time.
%
% To review Jones and Aref's notation: When alpha = pi/2, the particle is
% injected at "noon" on the injection well at (-a,0) and extracted at
% "noon" from the extraction well at (a,0). When alpha = 3pi/4, the
% particle is injected at 1:30 and extracted at 10:30. When alpha = pi, the
% particle is injected at 3:00 and extracted at 9:00, corresponding to a
% straight path from injection to extraction.
%
%   x and y are vectors (not matrices)
%   xy_in = [x_injection y_injection]
%   xy_out = [x_extraction y_extraction]
%   Q = "Q/b" or "flow per aquifer thickness"
%
% WARNING, deletes any particles on the x-axis between wells
%
% Rev. 0, dcm 10/29/2012, Matlab 7.7.0.471 (R2008b)
% Rev. 1, dcm 11/19/2012, Matlab 7.7.0.471 (R2008b)
% Rev. 2, dcm 3/23/2015, Matlab 8.2.0.701 (R2013b)

% disable graphics by default
if nargin < 7
	animate = 0;
end

% create x_well and y_well to match code from 10/29/2012
x_well = [xy_in(1) xy_out(1)];
y_well = [xy_in(2) xy_out(2)];

% count number of particles
np = length(x);

% half separation distance between wells (same in all coordinates)
a = sqrt( (x_well(1)-x_well(2) )^2 + (y_well(1)-y_well(2) )^2 )/2;

% make x and y into row vectors (not column ) for matrix multiplication
s = size(x);
if s(1) > 1
	x = x';
	y = y';
end

% PLOT #1: initial material line in (x,y) coordinates
if animate == 1
	plot(x,y,'k+')
	hold on
	plot(x_well,y_well,'bo:')
	plot(x_well(1),y_well(1),'bo:','MarkerFaceColor','blue')
	axis equal
	xlabel('x')
	ylabel('y')
	title('DIPOLEMAP.M')
	pause
	close
end

% shift coordinates such that injection at (-a,0), extraction at (a,0)
%   1. shift from (x,y) to (xs, ys) for "x shift" and "y shift"
%   2. rotate axes using equation (10.3.1) in Boas (1983) to find (xp, yp)

% 1. shift axes such that origin of (xs, ys) at midpoint of two wells
x_mid = mean(x_well);
y_mid = mean(y_well);
xs = x - x_mid;
ys = y - y_mid;
xs_well = x_well - x_mid;
ys_well = y_well - y_mid;

% 2. rotate axes by phi to find (xp, yp) for "x prime" and "y prime"
if ys_well(2) >= 0
	% 3/23/2015 c/o http://undocumentedmatlab.com/blog/matlab-numerical-gotchas
	% force argument of inverse trigonometric function to be -1 <= x <= 1
	% phi = acos(xs_well(2)/a);
	argument = xs_well(2)/a;
	argument = min(max(real(argument),-1),1);
	phi = acos(argument);
else
	% 3/23/2015 c/o http://undocumentedmatlab.com/blog/matlab-numerical-gotchas
	% force argument of inverse trigonometric function to be -1 <= x <= 1
	% phi = -acos(xs_well(2)/a);
	argument = xs_well(2)/a;
	argument = min(max(real(argument),-1),1);
	phi = -acos(argument);
end
M = [cos(phi) sin(phi); -sin(phi) cos(phi)];
% matrix multiplication gives only one output, so define XP = [xp; yp]
XP = M*[xs; ys];
xp = XP(1,:);
yp = XP(2,:);
XP_well = M*[xs_well; ys_well];
xp_well = XP_well(1,:);
yp_well = XP_well(2,:);
% 3/23/2015 remove particles on x-axis (not covered in Appendix of JA-1988)
index = find(abs(yp)>1e-6);
reduction = np - length(index);
np = np - reduction;
if reduction > 0
    fprintf('Removed %.0f particles from x-axis.\n',reduction);
end
if np == 0
	error('Sorry, all particles have now been extracted. GAME OVER!')
end
xp = xp(index);
yp = yp(index);

% PLOT #2: material line after shift and rotation in (x',y') coorinates
if animate == 1
	plot(xp,yp,'k+')
	hold on
	plot(xp_well,yp_well,'bo:')
	plot(xp_well(1),yp_well(1),'bo:','MarkerFaceColor','blue')
	axis equal
	xlabel('x''')
	ylabel('y''')
	title('DIPOLEMAP.M')
	pause
end

% mirror around x-axis defined by two wells
index_positive_yp = find(yp >= 0);
index_negative_yp = find(yp < 0);
yp(index_negative_yp) = -yp(index_negative_yp);

% PLOT #2 CONTINUED: plot material line after mirror
if animate == 1
	plot(xp(index_negative_yp),yp(index_negative_yp),'r+')
	pause
end

% calculate initial (alpha, theta)
[alpha, theta] = XYtoAlphaTheta(a,xp,yp);

% use Jones and Aref's (1988) equation (A3) for associated times
t = (2*pi*a^2)*(theta.*cos(alpha) - sin(alpha + theta) + sin(alpha))./(Q*sin(alpha).^3);

% calculate exit times for initial vector of alpha (equation A4)
T = (4*pi*a^2/Q)*((pi - alpha).*cot(alpha) + 1)./(sin(alpha)).^2;

% map particles forward in time
t = t + dt;
guess = 2*(pi - alpha);
for j = 1:np
	% test for extraction
	if t(j) > T(j)
		% put particle at extraction well
		theta(j) = guess(j);
	else
		% calculate theta = f(a,Q,alpha,t)
		% anonymous function "f" feeds extra parameters to FZERO
		% maximum possible range of theta is [0 guess(j)]
		f = @(theta)ja1988a3(theta,a,Q,alpha(j),t(j));
		theta(j) = fzero(f,[0 guess(j)]);
	end
end

% translate (alpha,theta) into (x,y)
[xp_out,yp_out] = AlphaThetatoXY(a,alpha,theta);

% PLOT #2 CONTINUED: material line after dipole but before mirror
if animate == 1
	plot(xp_out(index_positive_yp),yp_out(index_positive_yp),'ko', ...
		xp_out(index_negative_yp),yp_out(index_negative_yp),'ro')
	pause
end

% undo mirror about x-axis defined by two wells
yp(index_negative_yp) = -yp(index_negative_yp);
yp_out(index_negative_yp) = -yp_out(index_negative_yp);

% shift coordinates back to original (x,y)
%   3. rotate axes using equation (10.3.1) in Boas (1983) to find (xs, ys)
%   4. shift from (xs,ys) to (x,y)

% PLOT #2 CONTINUED: material line after dipole but before rotate and shift
if animate == 1
	plot(xp_out,yp_out,'ko')
	pause
	close
end

% 3. rotate axes by -phi to find (xp, yp) for "x prime" and "y prime"
M = [cos(-phi) sin(-phi); -sin(-phi) cos(-phi)];
% matrix multiplication gives only one output, so define XS_out = [xs_out; ys_out]
XS_out = M*[xp_out; yp_out];
xs_out = XS_out(1,:);
ys_out = XS_out(2,:);

% 4. shift from (xs, ys) to (x,y)
x_out = xs_out + x_mid;
y_out = ys_out + y_mid;

% PLOT #3: final material line
if animate == 1
	plot(x,y,'k+',x_out,y_out,'ko')
	hold on
	plot(x_well,y_well,'bo:')
	plot(x_well(1),y_well(1),'bo:','MarkerFaceColor','blue')
	legend('initial','final')
	axis equal
	xlabel('x')
	ylabel('y')
	title('DIPOLEMAP.M')
	pause
	close
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CUSTOM-4

function LHS = ja1988a3(theta,a,Q,alpha,t)

% This function is equation (A3) from Jones and Aref (1988), an analytical
% solution for particle tracking in the standard dipole flow. For physical
% parameters (a,Q) and particle path "alpha" at time "t", (A3) determines
% the particle position "theta".
%
% However, because the equation is implicit in "theta", it must be solved
% using numerical methods. Accordingly, this function puts (A3) into a form
% such that the desired "theta" will make the left hand side (LHS) zero.
% This function will then be solved with Matlab's FZERO function.
%
% Rev. 0, dcm 9/23/2009, Matlab 7.7.0.471 (R2008b)

LHS = theta*cos(alpha) - sin(alpha + theta) + sin(alpha) - Q*sin(alpha)^3*t/2/pi/a^2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CUSTOM-5

function [alpha,theta] = XYtoAlphaTheta(a,x,y)

% function [alpha,theta] = XYtoAlphaTheta(a,x,y)
%
% This function converts coordinate vectors (x,y) into (alpha,theta) for
% dipole wells centered at the origin and separated by distance 2a,
% following the notation of the Appendix to Jones and Aref (1988).
%
% See also AlphaThetatoXY.m
%
% Rev. 0, dcm 8/30/2012, Matlab 7.7.0.471 (R2008b)
% Rev. 1, dcm 11/19/2012, Matlab 7.7.0.471 (R2008b)
% Rev. 2, dcm 3/23/2015, Matlab 8.2.0.701 (R2013b)

% calculate alpha for each point (x,y)
l_alpha = 2*a;
l_in = sqrt((a-x).^2 + y.^2);
l_out = sqrt((a+x).^2 + y.^2);
% 3/23/2015 c/o http://undocumentedmatlab.com/blog/matlab-numerical-gotchas
% force argument of inverse trigonometric function to be -1 <= x <= 1
% alpha = acos((l_alpha.^2 - l_in.^2 - l_out.^2)./(-2*l_in.*l_out))
argument = (l_alpha.^2 - l_in.^2 - l_out.^2)./(-2*l_in.*l_out);
argument = min(max(real(argument),-1),1);
alpha = acos(argument);

% calculate magnitude of y-coordinate of circumcenter
R = a./sin(alpha);
y_D = sqrt(R.^2 - a^2);

% make y_D negative for points where alpha exceeds pi
index = find(alpha > pi/2);
y_D(index) = -y_D(index);

% calculate theta for each point (x,y)
l_theta = sqrt((x+a).^2 + y.^2);
l_mesa = sqrt(x.^2 + (y - y_D).^2);
l_point = sqrt(a^2 + y_D.^2);
% 3/23/2015 c/o http://undocumentedmatlab.com/blog/matlab-numerical-gotchas
% force argument of inverse trigonometric function to be -1 <= x <= 1
% theta = acos((l_theta.^2 - l_mesa.^2 - l_point.^2)./(-2*l_mesa.*l_point))
argument = (l_theta.^2 - l_mesa.^2 - l_point.^2)./(-2*l_mesa.*l_point);
argument = min(max(real(argument),-1),1);
theta = acos(argument);

% correction for x > a (added 11/19/2012)
index = find(x > a);
theta(index) = 2*pi - theta(index);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CUSTOM-6

function [x,y] = AlphaThetatoXY(a,alpha,theta)

% function [x,y] = AlphaThetatoXY(a,alpha,theta)
%
% This function converts coordinate vectors (alpha,theta) into (x,y) for
% dipole wells centered at the origin and separated by distance 2a,
% following the notation of the Appendix to Jones and Aref (1988).
%
% See equations 2-3 on page 2 of calculations 9/2/2009
% See also XYtoAlphaTheta.m
%
% Rev. 0, dcm 8/30/2012, Matlab 7.7.0.471 (R2008b)

% calculate x for each point (alpha,theta)
x = (-a./sin(alpha)).*cos(alpha - pi/2 + theta);

% calculate y for each point (alpha,theta)
y = (a./sin(alpha)).*sin(alpha - pi/2 + theta) - a./sin(alpha) + a*tan(pi/2 - alpha/2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CUSTOM-7

function D = separation(x,y)

% function D = separation(x,y)
%
% This function calculates the separation distance D for a series of N
% points specified by (x,y). Specifically, it inputs two N-vectors "x"
% and "y" and outputs NxN vector D, whose primary diagonal is zeros, and
% which is symmetric. See calculation 9/18/1999. Inspired by Ruben Juanes.
%
% Rev. 0, dcm 10/9/2009, Matlab 7.7.0.471 (R2008b)

% make inputs Nx1 vectors (if necessary)
sizex = size(x);
if sizex(2) > sizex(1)
    x = x';
end
sizey = size(y);
if sizey(2) > sizey(1)
    y = y';
end

% preliminary calculations
N = max(sizex);
I = ones(1,N);

% main calculations
Xa = x*I;
Xb = Xa';
Xc = Xa - Xb;
Ya = y*I;
Yb = Ya';
Yc = Ya - Yb;
D = (Xc.^2 + Yc.^2).^(1/2);
