function varargout = huckel_lab(varargin)
% HUCKEL_LAB Application M-file for huckel_lab.fig
%    FIG = HUCKEL_LAB launch huckel_lab GUI.
%    HUCKEL_LAB('callback_name', ...) invoke the named callback.

% Last Modified by GUIDE v2.5 17-Jan-2003 15:44:06

if nargin == 0  % LAUNCH GUI

	fig = openfig(mfilename, 'reuse');
	
	% Generate a structure of handles to pass to callbacks, and store it. 
	handles = guihandles(fig);
	handles.LocalCurrentDirectory = pwd;
	handles.MoleculeName = 'Untitled.hcl';
	
	figure_children = get(fig, 'Children');
	
    menu_File_handle = findobj(figure_children, 'flat', 'Tag', 'menu_File');
    menu_File_children = get(menu_File_handle, 'Children');
    menu_File_item_Save_handle = ...
	    findobj(menu_File_children, 'flat', 'Tag', 'menu_File_item_Save');
	handles.menu_File_item_Save = menu_File_item_Save_handle;

    menu_Plot_handle = findobj(figure_children, 'flat', 'Tag', 'menu_Plot');
    menu_Plot_children = get(menu_Plot_handle, 'Children');
    menu_Plot_item_SAN_handle = ...
	    findobj(menu_Plot_children, 'flat', 'Tag', 'menu_Plot_item_SAN');
	handles.menu_Plot_item_SAN = menu_Plot_item_SAN_handle;
	
	guidata(fig, handles);
	
	if nargout > 0
		varargout{1} = fig;
	end

    % Set some figure properties.
    set(fig, 'Color', get(0, 'defaultUicontrolBackgroundColor'));
    set(fig, 'Name', ['H', char(252), 'ckel Lab (', ...
	    handles.MoleculeName, ')']);    
    set(fig, 'DoubleBuffer', 'on');
    set(fig, 'HandleVisibility', 'callback')
    set(fig, 'CloseRequestFcn', ...
        'huckel_lab(''menu_File_Exit_Callback'',gcbo,[],guidata(gcbo))')
    
    % Initialize the axes.
    axesMolecule_Init(handles.axesMolecule, [], handles)
    axesEnergy_Init(handles.axesEnergy, [], handles)

    set(handles.popupmenu1, 'Value', 1)   % "Draw" Mode.

    set(handles.menu_File_item_Save, 'Enable', 'off')
    set(handles.menu_Plot_item_SAN, 'Checked', 'off')
    
    set(handles.pushbutton_down, 'String', 'Dwn')
    
    set(fig, 'Visible', 'on')
    % The GUI is saved with the property 'Visible' set to 'off'
    
elseif ischar(varargin{1}) % INVOKE NAMED SUBFUNCTION OR CALLBACK

	try
		if (nargout)
			[varargout{1:nargout}] = feval(varargin{:}); % FEVAL switchyard
		else
			feval(varargin{:}); % FEVAL switchyard
		end
	catch
		disp(lasterr);
	end

end


%| ABOUT CALLBACKS:
%| GUIDE automatically appends subfunction prototypes to this file, and 
%| sets objects' callback properties to call them through the FEVAL 
%| switchyard above. This comment describes that mechanism.
%|
%| Each callback subfunction declaration has the following form:
%| <SUBFUNCTION_NAME>(H, EVENTDATA, HANDLES, VARARGIN)
%|
%| The subfunction name is composed using the object's Tag and the 
%| callback type separated by '_', e.g. 'slider2_Callback',
%| 'figure1_CloseRequestFcn', 'axis1_ButtondownFcn'.
%|
%| H is the callback object's handle (obtained using GCBO).
%|
%| EVENTDATA is empty, but reserved for future use.
%|
%| HANDLES is a structure containing handles of components in GUI using
%| tags as fieldnames, e.g. handles.figure1, handles.slider2. This
%| structure is created at GUI startup using GUIHANDLES and stored in
%| the figure's application data using GUIDATA. A copy of the structure
%| is passed to each callback.  You can store additional information in
%| this structure at GUI startup, and you can change the structure
%| during callbacks.  Call guidata(h, handles) after changing your
%| copy to replace the stored original so that subsequent callbacks see
%| the updates. Type "help guihandles" and "help guidata" for more
%| information.
%|
%| VARARGIN contains any extra arguments you have passed to the
%| callback. Specify the extra arguments by editing the callback
%| property in the inspector. By default, GUIDE sets the property to:
%| <MFILENAME>('<SUBFUNCTION_NAME>', gcbo, [], guidata(gcbo))
%| Add any extra arguments after the last argument, before the final
%| closing parenthesis.



% --------------------------------------------------------------------
function varargout = pushbuttonCompute_Callback(h, eventdata, handles, varargin)

i = find(handles.Atoms);

if isempty(i)
  error_msg = 'Molecule is not correct.';
  title = ['H', char(252), 'ckel Lab: Error'];
  uiwait(errordlg(error_msg, title, 'modal'))
  return
end

handles.Bonds = full(handles.Bonds);
BondMatrix = handles.Bonds(i, i);
handles.Bonds = sparse(handles.Bonds);
TestVektor = sum(BondMatrix,1);
minbonds = min(TestVektor);
maxbonds = max(TestVektor);

if (minbonds == 0) | (maxbonds > 3)
  error_msg = 'Molecule is not correct.';
  title = ['H', char(252), 'ckel Lab: Error'];
  uiwait(errordlg(error_msg, title, 'modal'))
  return
end

[V, D] = eig(BondMatrix);

% Stoer matrisen lite foer att bryta degenerationer:

BondMatrix(1,1) = BondMatrix(1,1) + 1e-6;

[V, D_pert] = eig(BondMatrix);

x=-diag(D);

% Vi vaender paa vektorn med energier, saa att x(1) aer
% den laegsta och x(length(x)) den hoegsta energin:
x = x(length(x):-1:1);

% Vi goer motsvarande "shift" i matrisen med molekylorbitaler:  
Vikter = V(:,length(x):-1:1);

handles.energies = x;  % Notera att beta < 0.
                       % E = alpha - beta*x
		               % x(1): laegsta energin
		               % x(length(x)): hoegsta energin
handles.orbitals = Vikter;  % Varje kolonn i Vikter motsvarar en molekylorbital.
                            % Molekylorbitalen i kolonn "j"
        		            % har energin E = alpha - beta*x(j).
			                % Med idx = find(handles.Atoms) gäller:
			                % Koefficienten i rad "i" motsvarar
			                % atom(orbital) nummer idx(i).

handles.charge = 0;

set(handles.text4, 'String', ['Charge: ', num2str(handles.charge)])
			    
set(handles.popupmenu1, 'Value', 3)   % "Plot" Mode.  
handles.oldMode = 3;

set([handles.pushbutton_up handles.pushbutton_down], ...
    'Visible', 'on');
set([handles.pushbutton_eminus handles.pushbutton_eplus], ...
    'Visible', 'on');
set([handles.text3 handles.text4], ...
    'Visible', 'on');

empty_indices = find(~handles.Atoms);
empty_handles = handles.AtomHandles(empty_indices);

set(handles.GridHandles, ...
      'Visible', 'off');

set(empty_handles, ...
      'Visible', 'off');  
  
set(handles.AtomHandles, ...
      'ButtonDownFcn', ...
      '');
    
set(handles.BondHandles, ...
    'ButtonDownFcn', ...
    '');

% Radera (eventuell) atomnumrering.
set(handles.AtomLabelHandles, 'Visible', 'off')

txt = get(handles.menu_Plot_item_SAN, 'Checked');
if strcmp(txt, 'on')
  for j=1:length(x)
    set(handles.AtomLabelHandles(i(j)), ...
	'String', num2str(j), ...
	'Visible', 'on')
  end
end

% Radera (eventuellt) aktuellt energinivaadiagram:
% (Ifall anvaendaren trycker paa "Compute" i "Plot" Mode.)

if isfield(handles, 'LevelHandles')
  EnergiSuddning(handles)
  handles = guidata(handles.figure1);   % Laes in de foeraendringar som gjorts
                                        % inuti "EnergiSuddning".
end

% Anropa funktion som ritar energinivådiagram:
  
EnergiRitning(handles)

handles = guidata(handles.figure1);   % Laes in de foeraendringar som gjorts
                                      % inuti "EnergiRitning".
  
% Beraekna kortaste avstaandet mellan tvaa atomer i molekylen:
  
handles.min_bond_length = ComputeMinBondLength(handles);
  
% Anropa funktion som illustrerar molekylorbitalen
% med laegst energi grafiskt:
  
handles.currentMO = 1;

Molekylorbitalritning(handles)
  
guidata(handles.figure1, handles)   % Save the structure handles.



% --------------------------------------------------------------------
function varargout = pushbuttonClear_Callback(h, eventdata, handles, varargin)

set(handles.AtomHandles, ...
      'MarkerSize', handles.markersize_atom, ...
      'MarkerFaceColor', 'none', ...
      'MarkerEdgeColor', 'none', ...
      'ButtonDownFcn', ...
      'huckel_lab(''addatom_Callback'',gcbo,[],guidata(gcbo))', ...
      'Visible', 'on')
 
% Radera (eventuell) atomnumrering.
set(handles.AtomLabelHandles, 'Visible', 'off')
  
noa = length(handles.Atoms);
 
handles.Atoms = zeros(noa, 1);   % Atom-vector.
handles.Bonds = sparse(zeros(noa, noa));   % Bond-matrix.

delete(handles.BondHandles)
handles.BondHandles = [];   % Handles of bonds.

set(handles.popupmenu1, 'Value', 1)   % "Draw" Mode.
handles.oldMode = 1;

set([handles.pushbutton_up handles.pushbutton_down], ...
    'Visible', 'off');
set([handles.pushbutton_eminus handles.pushbutton_eplus], ...
    'Visible', 'off');
set([handles.text3 handles.text4], ...
    'Visible', 'off');

set(handles.GridHandles, ...
      'Visible', 'on');

set(handles.text2, 'String', '')

set(handles.menu_File_item_Save, 'Enable', 'on')
set(handles.figure1, 'Name', ['H', char(252), 'ckel Lab (', ...
	    handles.MoleculeName, '*)']);

if isfield(handles, 'LevelHandles')
  EnergiSuddning(handles)
  handles = guidata(h);   % Laes in de foeraendringar som gjorts
                          % inuti "EnergiSuddning".
end

guidata(h, handles)   % Save the structure handles.



% --------------------------------------------------------------------
function varargout = axesMolecule_Init(h, eventdata, handles, varargin)
% Initialize axesMolecule.
set(h, 'Color', [1 1 1])
set(h, 'XColor', [1 1 1])
set(h, 'YColor', [1 1 1])
set(h, 'XLimMode', 'manual');
set(h, 'YLimMode', 'manual');
set(h, 'XTickLabel', [])
set(h, 'YTickLabel', [])
set(h, 'XTick', [])
set(h, 'YTick', [])
set(h, 'Units', 'points');   % Character units are not square.
Position = get(h, 'Position');
set(h, 'Units', 'characters');   % Change back.
width = Position(3);   % Axes width in points.
height = Position(4);   % Axes height in points.
noxatoms = 15;   % Choose number of atoms in x-direction.
xmin = 0;
xmax = noxatoms+0.5;
xrange = xmax - xmin;
pointsperunit = width/xrange;
handles.pointsperunit = pointsperunit;
yrange = xrange*height/width;
% Compute number of atoms in y-direction.
% The expression for "noyatoms" follows from the demand
% that (the "padding") yendspace >= 1/2.
noyatoms = floor(1 + 2*(yrange-1)/sqrt(3));
yendspace = (yrange - sqrt(3)*(noyatoms-1)/2)/2;
ymin = 0;
ymax = ymin + yrange;
set(h, 'XLim', [xmin xmax], 'YLim', [ymin ymax]);

% Initialize data structures.
noa = noxatoms*noyatoms;   % Total number of atoms.

handles.Atoms = zeros(noa, 1);   % Atom-vector.
handles.Bonds = sparse(zeros(noa, noa));   % Bond-matrix.

handles.AtomHandles = zeros(noa, 1);   % Handles of atoms.
handles.GridHandles = zeros(noa, 1);   % Handles of grid-points.
handles.BondHandles = []; 	       % Handles of bonds.
handles.AtomLabelHandles = zeros(noa, 1); % Handles of atom number labels.

handles.Coords = zeros(noa, 2);   % Coordinate matrix.

% Plot gridpoints.
markersize_grid = ceil(pointsperunit/10);
handles.markersize_grid = markersize_grid;
markersize_atom = ceil(pointsperunit/2);
handles.markersize_atom = markersize_atom;
for i = 1:noyatoms
  y = yendspace + (i-1)*sqrt(3)/2;
  for j = 1:noxatoms
    if mod(i,2)
      x = j;   % i is odd
    else
      x = j-0.5;
    end
    atomnumber = (i-1)*noxatoms + j;
    handles.Coords(atomnumber,1) = x;
    handles.Coords(atomnumber,2) = y;

    handles.GridHandles(atomnumber) = line(x, y, ...
	    'Marker', 'x', ...
	    'MarkerSize', markersize_grid, ...
	    'Parent', h);
    tag = ['Atom:', num2str(atomnumber)];
    handles.AtomHandles(atomnumber) = line(x, y, ...
	    'Marker', 'o', ...
	    'MarkerSize', markersize_atom, ...
        'MarkerFaceColor', 'none', ...
	    'MarkerEdgeColor', 'none', ...
	    'ButtonDownFcn', ...
	    'huckel_lab(''addatom_Callback'',gcbo,[],guidata(gcbo))', ...
	    'Tag', tag, ...
	    'Parent', h);
    handles.AtomLabelHandles(atomnumber) = ...
	text('Position', [x y], ...
	     'Interpreter', 'tex', ...
             'String', '', ...
	     'HorizontalAlignment', 'center', ...
	     'VerticalAlignment', 'middle', ...
	     'Color', [0 0 0], ...
	     'Parent', h, ...
	     'Visible', 'off');
  end
end
guidata(h, handles)   % Save the structure handles.



% --------------------------------------------------------------------
function varargout = axesEnergy_Init(h, eventdata, handles, varargin)
% Initialize axesEnergy.

set(h, 'Color', [0.95 0.95 0.95])
set(h, 'XColor', [0.95 0.95 0.95])
set(h, 'YColor', [0.95 0.95 0.95])
set(h, 'XLimMode', 'manual');
set(h, 'YLimMode', 'manual');
set(h, 'XTickLabel', [])
set(h, 'YTickLabel', [])
set(h, 'XTick', [])
set(h, 'YTick', [])
set(h, 'YAxisLocation', 'right')
set(h, 'FontSize', 8)



% --------------------------------------------------------------------
function varargout = addatom_Callback(h, eventdata, handles, varargin)

tag = get(h, 'Tag');
pos = findstr(tag, ':');
atomnumber = eval(tag(pos+1:length(tag)));
handle = handles.AtomHandles(atomnumber);
set(handle, ...
    'MarkerFaceColor', 'blue', ...
    'MarkerEdgeColor', [0 0 0], ...
    'ButtonDownFcn', ...
    'huckel_lab(''bondplotstart_Callback'',gcbo,[],guidata(gcbo))')
handles.Atoms(atomnumber) = 1;
set(handles.menu_File_item_Save, 'Enable', 'on')
set(handles.figure1, 'Name', ['H', char(252), 'ckel Lab (', ...
	    handles.MoleculeName, '*)']);
guidata(gcbo,handles)   % Save the structure handles.



% --------------------------------------------------------------------
function varargout = bondplotstart_Callback(h, eventdata, handles, varargin)

tag = get(h, 'Tag');
pos = findstr(tag, ':');
start_atom_number = eval(tag(pos+1:length(tag)));

x = handles.Coords(start_atom_number,1);
y = handles.Coords(start_atom_number,2);

bondHandle = line([x x], [y y], ...
    'Color', [0 0 0], ...
    'LineWidth', 2, ...
    'Parent', handles.axesMolecule);

set(gcbf, 'WindowButtonMotionFcn', ...
    ['huckel_lab(''bondplotmove_Callback'',', ...
    num2str(bondHandle, 20), ')'])

set(gcbf, 'WindowButtonUpFcn', ...
    ['huckel_lab(''bondplotstop_Callback'',guidata(gcbo),', ...
    num2str(start_atom_number), ',', ...
    num2str(bondHandle, 20), ')'])



% --------------------------------------------------------------------
function bondplotmove_Callback(bondHandle)
CurrentPoint = get(gca, 'CurrentPoint');   % "gca" may be used since the figure's
                                           % "HandleVisibility" property is set
                                           % to "callback".
XYData = get(bondHandle, {'XData', 'YData'});
% XYData = {XData, YData}
XYData{1}(2) = CurrentPoint(1, 1); % XData(2) = x
XYData{2}(2) = CurrentPoint(1, 2); % YData(2) = y
set(bondHandle, {'XData', 'YData'}, XYData)



% --------------------------------------------------------------------
function bondplotstop_Callback(handles, start_atom_number, bondHandle)

set(gcbf, 'WindowButtonMotionFcn', '')
set(gcbf, 'WindowButtonUpFcn', '')

XYData = get(bondHandle, {'XData', 'YData'});

XData = XYData{1};
YData = XYData{2};

xstart = XData(1);
ystart = YData(1);
x = XData(2);
y = YData(2);

ok = 0;

% Beraekna indexen för alla atomer (besatta gitterpunkter):
atom_indices = find(handles.Atoms);

% Ta bort indexet för bindningens "start-atom":
atom_indices = atom_indices(find(atom_indices ~= start_atom_number));

for j = atom_indices'  % "atom_indices" aer en kolonnvektor
  
  xj = handles.Coords(j,1);
  yj = handles.Coords(j,2);
  distance = norm([x, y] - [xj, yj]);
  
  if (distance < 0.33) & (~handles.Bonds(start_atom_number,j))
    ok = 1;
    
    vektor = [xj - xstart, yj - ystart];
    l = norm(vektor);
    enhetsvektor = vektor/l;
    radius = 0.15;
    XData(1) = xstart + radius*enhetsvektor(1);
    YData(1) = ystart + radius*enhetsvektor(2);
    XData(2) = xj - radius*enhetsvektor(1);
    YData(2) = yj - radius*enhetsvektor(2);
    
    tag = ['Bond:', num2str(min([start_atom_number j])), ':', ...
	    num2str(max([start_atom_number j]))];
    
    set(bondHandle, 'XData', XData, 'YData', YData, 'Tag', tag)
    
    handles.Bonds(start_atom_number, j) = 1;
    handles.Bonds(j, start_atom_number) = 1;
    
    handles.BondHandles = [handles.BondHandles bondHandle];
    
    set(handles.menu_File_item_Save, 'Enable', 'on')
    set(handles.figure1, 'Name', ['H', char(252), 'ckel Lab (', ...
	    handles.MoleculeName, '*)']);
    
    guidata(gcbo,handles)   % Save the structure handles.
    
    break   % Avbryt for-loopen.
    
  end
  
end

if (~ok)
  delete(bondHandle)
end



% --------------------------------------------------------------------
function varargout = popupmenu1_Callback(h, eventdata, handles, varargin)

if isfield(handles, 'oldMode')
  oldMode = handles.oldMode;
else
  oldMode = 1;   % First time!
end

set(handles.text2, 'String', '')

occupied_indices = find(handles.Atoms);
empty_indices = find(~handles.Atoms);
    
occupied_handles = handles.AtomHandles(occupied_indices);
empty_handles = handles.AtomHandles(empty_indices);

val = get(h, 'Value');

handles.oldMode = val;
guidata(handles.figure1,handles)   % Save the structure handles.
  
switch val
    
case 1   % User selected "Draw".

  set([handles.pushbutton_up handles.pushbutton_down], ...
      'Visible', 'off');
  set([handles.pushbutton_eminus handles.pushbutton_eplus], ...
      'Visible', 'off');
  set([handles.text3 handles.text4], ...
      'Visible', 'off');
  
  set(handles.GridHandles, ...
      'Visible', 'on');
  
  set(occupied_handles, ...
        'MarkerSize', handles.markersize_atom, ...
        'MarkerFaceColor', 'blue', ...
	'MarkerEdgeColor', [0 0 0], ...
        'ButtonDownFcn', ...
        'huckel_lab(''bondplotstart_Callback'',gcbo,[],guidata(gcbo))');
  
  set(empty_handles, ...
        'ButtonDownFcn', ...
        'huckel_lab(''addatom_Callback'',gcbo,[],guidata(gcbo))', ...
	'Visible', 'on');

  set(handles.BondHandles, ...
      'ButtonDownFcn', ...
      '');
   
  % Radera (eventuell) atomnumrering.
  set(handles.AtomLabelHandles, 'Visible', 'off')
  
  if isfield(handles, 'LevelHandles')
    EnergiSuddning(handles)
  end
  
case 2   % User selected "Delete".

  set([handles.pushbutton_up handles.pushbutton_down], ...
      'Visible', 'off');
  set([handles.pushbutton_eminus handles.pushbutton_eplus], ...
      'Visible', 'off');
  set([handles.text3 handles.text4], ...
      'Visible', 'off');
  
  set(handles.GridHandles, ...
      'Visible', 'on');

  set(occupied_handles, ...
        'MarkerSize', handles.markersize_atom, ...
        'MarkerFaceColor', 'blue', ...
	'MarkerEdgeColor', [0 0 0], ...
        'ButtonDownFcn', ...
        'huckel_lab(''deleteatom_Callback'',gcbo,[],guidata(gcbo))');
  
  set(empty_handles, ...
        'ButtonDownFcn', ...
        '', ...
	'Visible', 'on');

  set(handles.BondHandles, ...
      'ButtonDownFcn', ...
      'huckel_lab(''deletebond_Callback'',gcbo,[],guidata(gcbo))'); 

  % Radera (eventuell) atomnumrering.
  set(handles.AtomLabelHandles, 'Visible', 'off')
  
  if isfield(handles, 'LevelHandles')
    EnergiSuddning(handles)
  end
  
case 3   % User selected "Plot".  
  
  i = find(handles.Atoms);

  if isempty(i)
  
    error_msg = 'Molecule is not correct!';
    title = ['H', char(252), 'ckel Lab: Error'];
    uiwait(errordlg(error_msg, title, 'modal'))
  
    set(handles.popupmenu1, 'Value', oldMode)   % "Draw" or "Delete".

    handles.oldMode = oldMode;
    guidata(handles.figure1,handles)   % Save the structure handles.
    
    return
  
  end

  handles.Bonds = full(handles.Bonds);
  BondMatrix = handles.Bonds(i, i);
  handles.Bonds = sparse(handles.Bonds);
  TestVektor = sum(BondMatrix,1);
  minbonds = min(TestVektor);
  maxbonds = max(TestVektor);

  if (minbonds == 0) | (maxbonds > 3)

    error_msg = 'Molecule is not correct!';
    title = ['H', char(252), 'ckel Lab: Error'];
    uiwait(errordlg(error_msg, title, 'modal'))
  
    set(handles.popupmenu1, 'Value', oldMode)   % "Draw" or "Delete".
  
    handles.oldMode = oldMode;
    guidata(handles.figure1,handles)   % Save the structure handles.
    
    return
  
  end

  pushbuttonCompute_Callback(h, [], handles)
      
end



% --------------------------------------------------------------------
function varargout = deleteatom_Callback(h, eventdata, handles, varargin)

tag = get(h, 'Tag');

pos = findstr(tag, ':');
atom_number = eval(tag(pos+1:length(tag)));

handle = handles.AtomHandles(atom_number);

set(handle, ...
    'MarkerFaceColor', 'none', ...
    'MarkerEdgeColor', 'none', ...
    'ButtonDownFcn', ...
    '')

handles.Atoms(atom_number) = 0;

% Atoms bonded to the deleted atom:
connected_atoms = find(handles.Bonds(atom_number, :));

for j = connected_atoms
    
    tag = ['Bond:', num2str(min([atom_number j])), ':', ...
	    num2str(max([atom_number j]))];

    handle = findobj(handles.BondHandles, 'flat', 'Tag', tag);
    % "findobj" may be used since the figure's "HandleVisibility"
    % property is set to "callback".
    
    delete(handle)
    
    handles.BondHandles = ...
        handles.BondHandles(find(handles.BondHandles ~= handle));
    
    handles.Bonds(atom_number,j) = 0;
    handles.Bonds(j,atom_number) = 0;
    
end

set(handles.menu_File_item_Save, 'Enable', 'on')
set(handles.figure1, 'Name', ['H', char(252), 'ckel Lab (', ...
	    handles.MoleculeName, '*)']);

guidata(h,handles)   % Save the structure handles.



% --------------------------------------------------------------------
function varargout = deletebond_Callback(h, eventdata, handles, varargin)

tag = get(h, 'Tag');

delete(h)

pos = findstr(tag, ':');
atom_number1 = eval(tag(pos(1)+1:pos(2)-1));
atom_number2 = eval(tag(pos(2)+1:length(tag)));

handles.Bonds(atom_number1,atom_number2) = 0;
handles.Bonds(atom_number2,atom_number1) = 0;

handles.BondHandles = ...
        handles.BondHandles(find(handles.BondHandles ~= h));

set(handles.menu_File_item_Save, 'Enable', 'on')
set(handles.figure1, 'Name', ['H', char(252), 'ckel Lab (', ...
	    handles.MoleculeName, '*)']);    
    
guidata(handles.figure1,handles)   % Save the structure handles.



% --------------------------------------------------------------------
function varargout = menu_File_New_Callback(h, eventdata, handles, varargin)

if strcmp(get(handles.menu_File_item_Save, 'Enable'), 'on')
  title = ['H', char(252), 'ckel Lab: Warning'];
  question = ['Save changes to ', handles.MoleculeName, '?'];
  answer = questdlg(question, title);
  if strcmp(answer, 'Yes')
    Atoms = handles.Atoms;
    Bonds = handles.Bonds;
    save(fullfile(handles.LocalCurrentDirectory, handles.MoleculeName), ...
      'Atoms', 'Bonds', '-mat')   
  elseif strcmp(answer, 'Cancel')
    return
  end
end

huckel_lab('pushbuttonClear_Callback',h,[],handles)
handles = guidata(h); % Laes in de foeraendringar som gjorts
                      % inuti "pushbuttonClear_Callback".
handles.MoleculeName = 'Untitled.hcl';
set(handles.menu_File_item_Save, 'Enable', 'off')
set(handles.figure1, 'Name', ['H', char(252), 'ckel Lab (', ...
        handles.MoleculeName, ')']);
guidata(h, handles)   % Save the structure handles.
  


% --------------------------------------------------------------------
function varargout = menu_File_Open_Callback(h, eventdata, handles, varargin)

if strcmp(get(handles.menu_File_item_Save, 'Enable'), 'on')
  title = ['H', char(252), 'ckel Lab: Warning'];
  question = ['Save changes to ', handles.MoleculeName, '?'];
  answer = questdlg(question, title);
  if strcmp(answer, 'Yes')
    Atoms = handles.Atoms;
    Bonds = handles.Bonds;
    save(fullfile(handles.LocalCurrentDirectory, handles.MoleculeName), ...
      'Atoms', 'Bonds', '-mat')
    set(handles.figure1, 'Name', ['H', char(252), 'ckel Lab (', ...
	    handles.MoleculeName, ')']);
  elseif strcmp(answer, 'Cancel')
    return
  end
end

pos_size = get(handles.figure1, 'Position');
% pos_size = [left, bottom, width, height]
dlg_pos = [pos_size(1)+pos_size(3)/5 pos_size(2)+pos_size(4)/5];
% dlg_pos = [left_dlg, bottom_dlg] = [left+width/5 bottom+height/5]
[handles.MoleculeName, handles.LocalCurrentDirectory, cancel] = ...
    molecule_name_dlg(handles.MoleculeName, ...
                      handles.LocalCurrentDirectory, dlg_pos, 'Open');
if cancel
  return
end

try
  data = load(fullfile(handles.LocalCurrentDirectory, handles.MoleculeName), '-mat');
catch
  title = ['H', char(252), 'ckel Lab: Error'];
  uiwait(errordlg(lasterr, title, 'modal'))
  return
end

% Kompatibilitetstest:

ok = 1;

noa_file = length(data.Atoms);
noa_current_grid = length(handles.Atoms);

if (noa_file < noa_current_grid)   % Molecule on file saved
                                   % on a smaller grid.

  temp_Atoms = zeros(noa_current_grid, 1);
  temp_Atoms(1:noa_file) = data.Atoms;
  data.Atoms = temp_Atoms;
  
  temp_Bonds = zeros(noa_current_grid, noa_current_grid);
  temp_Bonds(1:noa_file, 1:noa_file) = full(data.Bonds);
  data.Bonds = sparse(temp_Bonds);
  
end

if (noa_file > noa_current_grid)   % Molecule on file saved
                                   % on a larger grid.

  overflow = data.Atoms(noa_current_grid+1:noa_file);
  if max(overflow) == 1
    ok = 0;
  else   % No atoms in the "lost" part of the grid.
    data.Atoms = data.Atoms(1:noa_current_grid);
    data.Bonds = full(data.Bonds);
    data.Bonds = data.Bonds(1:noa_current_grid, 1:noa_current_grid);
    data.Bonds = sparse(data.Bonds);
  end

end				   

if ~ok
  title = ['H', char(252), 'ckel Lab: Error'];
  error_text = ...
      ['The molecule ', handles.MoleculeName, ' was saved on a larger grid.'];
  uiwait(errordlg(error_text, title, 'modal'))
  return 
end

huckel_lab('pushbuttonClear_Callback',h,[],handles) % Nb: "handles" sparas inuti
                                                    % "pushbuttonClear_Callback".
handles = guidata(h); % Laes in de foeraendringar som gjorts
                      % inuti "pushbuttonClear_Callback".
  
handles.Atoms = data.Atoms;
handles.Bonds = data.Bonds;
  
occupied_indices = find(handles.Atoms);  
occupied_handles = handles.AtomHandles(occupied_indices);
  
set(occupied_handles, ...
      'MarkerFaceColor', 'blue', ...
      'MarkerEdgeColor', [0 0 0], ...
      'ButtonDownFcn', ...
      'huckel_lab(''bondplotstart_Callback'',gcbo,[],guidata(gcbo))');
  
[i,j] = find(handles.Bonds);  % "i": radindex   "j": kolonnindex
row_indices = i(find(i<j));   % Oevre delen av "handles.Bonds".
column_indices = j(find(i<j));
  
for k = 1:length(row_indices)
    
  i = row_indices(k);
  j = column_indices(k);
    
  xi = handles.Coords(i,1);
  yi = handles.Coords(i,2);
    
  xj = handles.Coords(j,1);
  yj = handles.Coords(j,2);
    
  vektor = [xj - xi, yj - yi];
  l = norm(vektor);
  enhetsvektor = vektor/l;
  radius = 0.15;
        
  XData(1) = xi + radius*enhetsvektor(1);
  YData(1) = yi + radius*enhetsvektor(2);
  XData(2) = xj - radius*enhetsvektor(1);
  YData(2) = yj - radius*enhetsvektor(2);
    
  tag = ['Bond:', num2str(i), ':', num2str(j)]; % Notera att i < j.
    
  handle = line(XData, YData, ...
      'Color', [0 0 0], ...
      'LineWidth', 2, ...
      'Parent', handles.axesMolecule, ...
      'Tag', tag);    
        
  handles.BondHandles = [handles.BondHandles handle];
    
end

set(handles.menu_File_item_Save, 'Enable', 'off')
set(handles.figure1, 'Name', ['H', char(252), 'ckel Lab (', ...
        handles.MoleculeName, ')']);
  
guidata(h,handles)   % Save the structure handles.



% --------------------------------------------------------------------
function varargout = menu_File_Exit_Callback(h, eventdata, handles, varargin)

if strcmp(get(handles.menu_File_item_Save, 'Enable'), 'on')
  title = ['H', char(252), 'ckel Lab: Warning'];
  question = ['Save changes to ', handles.MoleculeName, '?'];
  answer = questdlg(question, title);
  if strcmp(answer, 'Yes')
    Atoms = handles.Atoms;
    Bonds = handles.Bonds;
    save(fullfile(handles.LocalCurrentDirectory, handles.MoleculeName), ...
      'Atoms', 'Bonds', '-mat')   
  elseif strcmp(answer, 'Cancel')
    return
  end
end

closereq



% --------------------------------------------------------------------
function varargout = menu_File_Save_Callback(h, eventdata, handles, varargin)
Atoms = handles.Atoms;
Bonds = handles.Bonds;
save(fullfile(handles.LocalCurrentDirectory, handles.MoleculeName), 'Atoms', 'Bonds', '-mat')
set(handles.menu_File_item_Save, 'Enable', 'off')
set(handles.figure1, 'Name', ['H', char(252), 'ckel Lab (', ...
	    handles.MoleculeName, ')']);



% --------------------------------------------------------------------
function varargout = menu_File_SaveAs_Callback(h, eventdata, handles, varargin)
pos_size = get(handles.figure1, 'Position');
% pos_size = [left, bottom, width, height]
dlg_pos = [pos_size(1)+pos_size(3)/5 pos_size(2)+pos_size(4)/5];
% dlg_pos = [left_dlg, bottom_dlg] = [left+width/5 bottom+height/5]
[handles.MoleculeName, handles.LocalCurrentDirectory, cancel] = ...
    molecule_name_dlg(handles.MoleculeName, ...
                      handles.LocalCurrentDirectory, dlg_pos, 'Save');
if cancel
  return
end

dir_struct = dir(handles.LocalCurrentDirectory);
file_names = {dir_struct.name}';
file_names = file_names(~[dir_struct.isdir]');
file_exists = 'no';
for i = 1:length(file_names)
  if strcmp(deblank(lower(file_names{i})), deblank(lower(handles.MoleculeName)))
    file_exists = 'yes';
    break   % Avbryt for-loopen.
  end
end

if strcmp(file_exists, 'yes')
  title = ['H', char(252), 'ckel Lab: Warning'];
  question = ['The file ', handles.MoleculeName, ' already exists. Overwrite?'];
  answer = questdlg(question, title, 'Cancel');
  if (strcmp(answer, 'No') | strcmp(answer, 'Cancel'))
    return
  end
end
  
Atoms = handles.Atoms;
Bonds = handles.Bonds;
save(fullfile(handles.LocalCurrentDirectory, handles.MoleculeName), 'Atoms', 'Bonds', '-mat')
set(handles.menu_File_item_Save, 'Enable', 'off')
set(handles.figure1, 'Name', ['H', char(252), 'ckel Lab (', ...
	handles.MoleculeName, ')']);
guidata(h, handles)   % Save the structure handles.



% --------------------------------------------------------------------
function varargout = menu_Plot_SAN_Callback(h, eventdata, handles, varargin)
txt = get(handles.menu_Plot_item_SAN, 'Checked');
if strcmp(txt, 'on')
    set(handles.menu_Plot_item_SAN, 'Checked', 'off')
    % Ta bort numrering om vi aer i "Draw Mode":
    mode = get(handles.popupmenu1, 'Value');
    if mode==3
      % Radera atomnumrering.
      set(handles.AtomLabelHandles, 'Visible', 'off')
    end
else
    set(handles.menu_Plot_item_SAN, 'Checked', 'on')
    % Laegg till numrering om vi aer i "Draw Mode":
    mode = get(handles.popupmenu1, 'Value');
    if mode==3
      i = find(handles.Atoms);
      for j=1:length(i)
        set(handles.AtomLabelHandles(i(j)), ...
	  'String', num2str(j), ...
	  'Visible', 'on')
      end
    end
end



% --------------------------------------------------------------------
function varargout = menu_Export_ED_Callback(h, eventdata, handles, varargin)

i = find(handles.Atoms);

if isempty(i)
  error_msg = 'Molecule is not correct.';
  title = ['H', char(252), 'ckel Lab: Error'];
  uiwait(errordlg(error_msg, title, 'modal'))
  return
end

handles.Bonds = full(handles.Bonds);
BondMatrix = handles.Bonds(i, i);
handles.Bonds = sparse(handles.Bonds);
TestVektor = sum(BondMatrix,1);
minbonds = min(TestVektor);
maxbonds = max(TestVektor);

if (minbonds == 0) | (maxbonds > 3)
  error_msg = 'Molecule is not correct.';
  title = ['H', char(252), 'ckel Lab: Error'];
  uiwait(errordlg(error_msg, title, 'modal'))
  return
end

[V, D] = eig(BondMatrix);

% Stoer matrisen lite foer att bryta degenerationer:

BondMatrix(1,1) = BondMatrix(1,1) + 1e-6;

[V, D_pert] = eig(BondMatrix);

BondMatrix(1,1) = BondMatrix(1,1) - 1e-6;

x=-diag(D);

% Vi vaender paa vektorn med energier, saa att x(1) aer
% den laegsta och x(length(x)) den hoegsta energin:
x = x(length(x):-1:1);

% Vi goer motsvarande "shift" i matrisen med molekylorbitaler:  
Vikter = V(:,length(x):-1:1);

assignin('base', 'Huckel_BondMatrix', BondMatrix);
assignin('base', 'Huckel_Energies', -x);   % Exportera egenvärden!
assignin('base', 'Huckel_Orbitals', Vikter);



% --------------------------------------------------------------------
function varargout = ritaMO_Callback(h, eventdata, handles, varargin)
tag = get(h, 'Tag');
pos = findstr(tag, ':');
MO_number = eval(tag(pos+1:length(tag)));
handles.currentMO = MO_number;
guidata(h, handles)   % Save the structure handles.
Molekylorbitalritning(handles)



% --------------------------------------------------------------------
function Molekylorbitalritning(handles)

MO_number = handles.currentMO;

set(handles.LevelHandles, 'LineStyle', '-')
set(handles.LevelHandles(MO_number), 'LineStyle', ':')

weights = handles.orbitals(:, MO_number);
energy = handles.energies(MO_number);

max_abs_weight = max(abs(weights));

idx = find(handles.Atoms);

noa = length(idx);

string = ['MO number: ', num2str(MO_number)];
string = [string, '   '];
string = [string, 'Eigenvalue: ', sprintf('%0.3f', -energy)];

set(handles.text2, 'String', string)

tolerance = 1e-9;

for i = 1:noa

  handle = handles.AtomHandles(idx(i));
  
  size = handles.min_bond_length*handles.pointsperunit;
  size = size*abs(weights(i))/max_abs_weight*0.9;
  
  if abs(weights(i)) > tolerance

    if weights(i)>0
      color = 'green';
    else
      color = 'red';
    end

  else

    color = 'white';
    size = handles.markersize_grid;
    
  end

  set(handle, ...
     'MarkerSize', size, ...
     'MarkerFaceColor', color, ...
     'MarkerEdgeColor', 'none') 
  
end



% --------------------------------------------------------------------
function min_bond_length = ComputeMinBondLength(handles)

[i,j] = find(handles.Bonds);  % "i": radindex   "j": kolonnindex
row_indices = i(find(i<j));   % Oevre delen av "handles.Bonds".
column_indices = j(find(i<j));
  
for k = 1:length(row_indices)
    
  i = row_indices(k);
  j = column_indices(k);
    
  xi = handles.Coords(i,1);
  yi = handles.Coords(i,2);
    
  xj = handles.Coords(j,1);
  yj = handles.Coords(j,2);
    
  vektor = [xj - xi, yj - yi];
  bondlengths(k) = norm(vektor);

end

min_bond_length = min(bondlengths);



% --------------------------------------------------------------------
function EnergiRitning(handles)

energies = handles.energies;

number_of_orbitals = length(energies);

degeneracy = [];
MOs = [];

tolerance = 1e-9;

distinct_energy_level = 1;

degeneracy(distinct_energy_level) = 1;
MOs(distinct_energy_level, :) = [1 1];

for i = 2:number_of_orbitals
  if (energies(i) - energies(i-1)) < tolerance   % Multiple eigen-value.
    degeneracy(distinct_energy_level) = ...
	degeneracy(distinct_energy_level) + 1;
    MOs(distinct_energy_level, 2) = MOs(distinct_energy_level, 2) + 1;
  else
    distinct_energy_level = distinct_energy_level + 1;
    degeneracy(distinct_energy_level) = 1;
    MOs(distinct_energy_level, :) = [i i];
  end
end

distinct_energies = energies(MOs(:,1));

no_of_distinct_energies = length(distinct_energies);

degeneracy = degeneracy';   % Make it a column vector.

% Example: Benzene
%
% energies = [-2; -1; -1; 1; 1; 2]
%
% distinct_energies = [-2; -1; 1; 2]
%
% no_of_distinct_energies = 4
%
% degeneracy = [1; 2; 2; 1]
%
% MOs = [1 1
%        2 3
%        4 5
%        6 6]
%
% (The rows in "degeneracy" and "MOs" correspond
%  to *distinct* energies.
%  First column in "MOs": first MO with this energy.
%  Second column in "MOs": last MO with this energy.)

energy_span = distinct_energies(no_of_distinct_energies) - ...
    distinct_energies(1);

% 0.1 <= y <= 0.9
YTick = 0.1 + ...
    0.8*(distinct_energies - distinct_energies(1))/energy_span;

YTickLabel = cell(no_of_distinct_energies, 1);

for i = 1:no_of_distinct_energies
  YTickLabel{i} = sprintf('%0.3f', -distinct_energies(i));
end

set(handles.axesEnergy, 'YColor', [0 0 0])
set(handles.axesEnergy, 'YTick', YTick)
set(handles.axesEnergy, 'YTickLabel', YTickLabel)
% handles.energy_text1 = text('Position', [0.5 0.985], ...
%     'String', 'E = \alpha+x\mid\beta\mid', ...
%     'Interpreter', 'tex', ...
%     'HorizontalAlignment', 'center', ...
%     'VerticalAlignment', 'cap', ...
%     'Parent', handles.axesEnergy);
handles.energy_text1 = text('Position', [0.5 0.985], ...
    'String', 'E = \beta\lambda', ...
    'Interpreter', 'tex', ...
    'HorizontalAlignment', 'center', ...
    'VerticalAlignment', 'cap', ...
    'Parent', handles.axesEnergy);
% handles.energy_text2 = text('Position', [1.05 0.985], ...
%     'String', 'x', ...
%     'FontSize', 8, ...
%     'HorizontalAlignment', 'left', ...
%     'VerticalAlignment', 'cap', ...
%     'Parent', handles.axesEnergy);
handles.energy_text2 = text('Position', [1.05 0.985], ...
    'String', '\lambda', ...
    'Interpreter', 'tex', ...
    'FontSize', 8, ...
    'HorizontalAlignment', 'left', ...
    'VerticalAlignment', 'cap', ...
    'Parent', handles.axesEnergy);

linewidth = 0.1;

handles.LevelHandles = [];
handles.ElectronHandles = [];

for i = 1:no_of_distinct_energies
  ypos = YTick(i);
  first_electron_level_i = 2*MOs(i,1) - 1;
  if degeneracy(i) == 1
    xpos = 0.5;
  elseif degeneracy(i) == 2
    xpos = [0.3 0.7];
  else
    linewidth = 1/(3*degeneracy(i) + 1);
    xpos = ((1:degeneracy(i))*3 - 1)*linewidth;
  end
  for j = 1:degeneracy(i)

    electron_number_up = first_electron_level_i + j - 1;
    electron_number_down = electron_number_up + degeneracy(i);

    handles.ElectronHandles(electron_number_up) = ...
	text('Position', [xpos(j)-linewidth/3 ypos], ...
	'String', '\uparrow', ...
	'Interpreter', 'tex', ...
	'HorizontalAlignment', 'center', ...
	'VerticalAlignment', 'middle', ...
	'Parent', handles.axesEnergy, ...
	'FontSize', 8, ...
	'Visible', 'off');

    handles.ElectronHandles(electron_number_down) = ...
	text('Position', [xpos(j)+linewidth/3 ypos], ...
	'String', '\downarrow', ...
	'Interpreter', 'tex', ...
	'HorizontalAlignment', 'center', ...
	'VerticalAlignment', 'middle', ...
	'Parent', handles.axesEnergy, ...
	'FontSize', 8, ...
	'Visible', 'off');
    
    MO_number = MOs(i,1) + j - 1;
    
    tag = ['MO:', num2str(MO_number)];
    
    handle = line([xpos(j)-linewidth xpos(j)+linewidth], ...
	[ypos ypos], ...
	'Parent', handles.axesEnergy, ...
        'Tag', tag, ...
	'ButtonDownFcn', ...
	'huckel_lab(''ritaMO_Callback'',gcbo,[],guidata(gcbo))');
    handles.LevelHandles(MO_number) = handle;

  end
end

net_charge = handles.charge;

number_of_electrons = number_of_orbitals - net_charge;

set(handles.ElectronHandles(1:number_of_electrons), 'Visible', 'on')

guidata(handles.figure1, handles)   % Save the structure handles.



% --------------------------------------------------------------------
function EnergiSuddning(handles)

set(handles.axesEnergy, 'YColor', [0.95 0.95 0.95])
set(handles.axesEnergy, 'YTick', [])
set(handles.axesEnergy, 'YTickLabel', [])

delete(handles.energy_text1)
delete(handles.energy_text2)
delete(handles.LevelHandles)
delete(handles.ElectronHandles)

handles = rmfield(handles, {'energy_text1', 'energy_text2', ...
    'LevelHandles', 'ElectronHandles'});

guidata(handles.figure1, handles)   % Save the structure handles.



% --------------------------------------------------------------------
function varargout = pushbutton_up_Callback(h, eventdata, handles, varargin)

number_of_orbitals = length(handles.energies);

handles.currentMO = handles.currentMO + 1;

if handles.currentMO > number_of_orbitals
  handles.currentMO = 1;
end

guidata(h, handles)   % Save the structure handles.

Molekylorbitalritning(handles)



% --------------------------------------------------------------------
function varargout = pushbutton_down_Callback(h, eventdata, handles, varargin)

number_of_orbitals = length(handles.energies);

handles.currentMO = handles.currentMO - 1;

if handles.currentMO < 1
  handles.currentMO = number_of_orbitals;
end

guidata(h, handles)   % Save the structure handles.

Molekylorbitalritning(handles)



% --------------------------------------------------------------------
function varargout = pushbutton_eplus_Callback(h, eventdata, handles, varargin)

number_of_orbitals = length(handles.energies);

handles.charge = handles.charge - 1;

if handles.charge < (-number_of_orbitals)
  handles.charge = -number_of_orbitals;
end

number_of_electrons = number_of_orbitals - handles.charge;

set(handles.ElectronHandles(1:2*number_of_orbitals), 'Visible', 'off')
set(handles.ElectronHandles(1:number_of_electrons), 'Visible', 'on')

if handles.charge > 0
  charge_string = ['+', num2str(handles.charge)];
else
  charge_string = num2str(handles.charge);
end
set(handles.text4, 'String', ['Charge: ', charge_string])

guidata(h, handles)   % Save the structure handles.



% --------------------------------------------------------------------
function varargout = pushbutton_eminus_Callback(h, eventdata, handles, varargin)

number_of_orbitals = length(handles.energies);

handles.charge = handles.charge + 1;

if handles.charge > number_of_orbitals
  handles.charge = number_of_orbitals;
end

number_of_electrons = number_of_orbitals - handles.charge;

set(handles.ElectronHandles(1:2*number_of_orbitals), 'Visible', 'off')
set(handles.ElectronHandles(1:number_of_electrons), 'Visible', 'on')

if handles.charge > 0
  charge_string = ['+', num2str(handles.charge)];
else
  charge_string = num2str(handles.charge);
end
set(handles.text4, 'String', ['Charge: ', charge_string])

guidata(h, handles)   % Save the structure handles.
