function varargout = five_membered(varargin)

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @five_membered_OpeningFcn, ...
                   'gui_OutputFcn',  @five_membered_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


% --- Executes just before five_membered is made visible.
function five_membered_OpeningFcn(hObject, eventdata, handles, varargin)

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

% Update handles structure
warning off all;




handles.dir = 'C:\pseudo_V1.00\';





handles.karplus = 1;
handles.output_format = '-djpeg100';
handles.output_ext = '.jpg';
handles.output_dpi = '-r300';
handles.sums = zeros(0,25);
handles.atoms_selected_value = 0;
handles.atoms_selected = [0 0 0 0 0 0 0 0 0 0 0 0];
handles.parameters = [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1];
handles.electronegativities = [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];

handles.all_couplings = [1 2      0 0 0 0 0 0 0   8 16 13  7  1 1 0 ; 
                         1 8   -120 0 0 0 0 0 0  16  2 13  7  1 1 0 ; 
                         7 8      0 0 0 0 0 0 0  16  2  1 13  1 1 0 ; 
                         2 7    120 0 0 0 0 0 0   8 16  1 13  1 1 0 ; 
                         2 3      0 0 0 0 0 0 0   9 18 15  8  1 1 0 ; 
                         2 9   -120 0 0 0 0 0 0  18  3 15  8  1 1 0 ; 
                         8 9      0 0 0 0 0 0 0  18  3  2 15  1 1 0 ; 
                         3 8    120 0 0 0 0 0 0   9 18  2 15  1 1 0 ; 
                         3 4      0 0 0 0 0 0 0  10 20 17  9  1 1 0 ; 
                         3 10  -120 0 0 0 0 0 0  20  4 17  9  1 1 0 ; 
                         9 10     0 0 0 0 0 0 0  20  4  3 17  1 1 0 ; 
                         4 9    120 0 0 0 0 0 0  10 20  3 17  1 1 0 ; 
                         4 5      0 0 0 0 0 0 0  11 22 19 10  1 1 0 ; 
                         4 11  -120 0 0 0 0 0 0  22  5 19 10  1 1 0 ; 
                         10 11    0 0 0 0 0 0 0  22  5  4 19  1 1 0 ; 
                         5 10   120 0 0 0 0 0 0  11 22  4 19  1 1 0 ; 
                         1 5      0 0 0 0 0 0 0  7  14 21 11  1 1 0 ; 
                         5 7   -120 0 0 0 0 0 0  14  1 21 11  1 1 0 ; 
                         7 11     0 0 0 0 0 0 0  14  1  5 21  1 1 0 ; 
                         1 11   120 0 0 0 0 0 0  7  14  5 21  1 1 0]; 

guidata(hObject, handles);

% --- Outputs from this function are returned to the command line.
function varargout = five_membered_OutputFcn(hObject, eventdata, handles) 
varargout{1} = handles.output;


% --------------------------------------------------------------------
function couplings_Callback(hObject, eventdata, handles)

update_constants(hObject, eventdata, handles);
update_en(hObject, eventdata, handles);

contents = get(handles.couplings, 'Value' );
num = size(contents);
num = num(2);

set( handles.t1,'value' , 0);
set( handles.t2,'value' , 0);
set( handles.t3,'value' , 0);
set( handles.t4,'value' , 0);
set( handles.t5,'value' , 0);
set( handles.b1,'value' , 0);
set( handles.b2,'value' , 0);
set( handles.b3,'value' , 0);
set( handles.b4,'value' , 0);
set( handles.b5,'value' , 0);

for iii = 1:20
 handles.atoms_selected(iii) = 0;
end
handles.atoms_selected_value = 0;

ii = 0;
if num == 1
 for i = 1:20
     if handles.all_couplings(i,4) == 1
         ii = ii + 1;
     end
     if ii == contents
       ii = 21;
       
       switch i
           case 1
          set( handles.t1,'value' , 1);
          set( handles.t2,'value' , 1);
          handles.atoms_selected(handles.all_couplings(i,1)) = 1;
          handles.atoms_selected(handles.all_couplings(i,2)) = 1;
           case 2
          set( handles.t1,'value' , 1);
          set( handles.b2,'value' , 1);
          handles.atoms_selected(handles.all_couplings(i,1)) = 1;
          handles.atoms_selected(handles.all_couplings(i,2)) = 1;
           case 3
          set( handles.b1,'value' , 1);
          set( handles.b2,'value' , 1);
          handles.atoms_selected(handles.all_couplings(i,1)) = 1;
          handles.atoms_selected(handles.all_couplings(i,2)) = 1;
           case 4
          set( handles.b1,'value' , 1);
          set( handles.t2,'value' , 1);
          handles.atoms_selected(handles.all_couplings(i,1)) = 1;
          handles.atoms_selected(handles.all_couplings(i,2)) = 1;
           case 5
          set( handles.t2,'value' , 1);
          set( handles.t3,'value' , 1);
          handles.atoms_selected(handles.all_couplings(i,1)) = 1;
          handles.atoms_selected(handles.all_couplings(i,2)) = 1;
           case 6
          set( handles.t2,'value' , 1);
          set( handles.b3,'value' , 1);
          handles.atoms_selected(handles.all_couplings(i,1)) = 1;
          handles.atoms_selected(handles.all_couplings(i,2)) = 1;
           case 7
          set( handles.b2,'value' , 1);
          set( handles.b3,'value' , 1);
          handles.atoms_selected(handles.all_couplings(i,1)) = 1;
          handles.atoms_selected(handles.all_couplings(i,2)) = 1;
           case 8
          set( handles.b2,'value' , 1);
          set( handles.t3,'value' , 1);
          handles.atoms_selected(handles.all_couplings(i,1)) = 1;
          handles.atoms_selected(handles.all_couplings(i,2)) = 1;
           case 9
          set( handles.t3,'value' , 1);
          set( handles.t4,'value' , 1);
          handles.atoms_selected(handles.all_couplings(i,1)) = 1;
          handles.atoms_selected(handles.all_couplings(i,2)) = 1;
           case 10
          set( handles.t3,'value' , 1);
          set( handles.b4,'value' , 1);
          handles.atoms_selected(handles.all_couplings(i,1)) = 1;
          handles.atoms_selected(handles.all_couplings(i,2)) = 1;
           case 11
          set( handles.b3,'value' , 1);
          set( handles.b4,'value' , 1);
          handles.atoms_selected(handles.all_couplings(i,1)) = 1;
          handles.atoms_selected(handles.all_couplings(i,2)) = 1;
           case 12
          set( handles.b3,'value' , 1);
          set( handles.t4,'value' , 1);
          handles.atoms_selected(handles.all_couplings(i,1)) = 1;
          handles.atoms_selected(handles.all_couplings(i,2)) = 1;
           case 13
          set( handles.t4,'value' , 1);
          set( handles.t5,'value' , 1);
          handles.atoms_selected(handles.all_couplings(i,1)) = 1;
          handles.atoms_selected(handles.all_couplings(i,2)) = 1;
           case 14
          set( handles.t4,'value' , 1);
          set( handles.b5,'value' , 1);
          handles.atoms_selected(handles.all_couplings(i,1)) = 1;
          handles.atoms_selected(handles.all_couplings(i,2)) = 1;
           case 15
          set( handles.b4,'value' , 1);
          set( handles.b5,'value' , 1);
          handles.atoms_selected(handles.all_couplings(i,1)) = 1;
          handles.atoms_selected(handles.all_couplings(i,2)) = 1;
           case 16
          set( handles.b4,'value' , 1);
          set( handles.t5,'value' , 1);
          handles.atoms_selected(handles.all_couplings(i,1)) = 1;
          handles.atoms_selected(handles.all_couplings(i,2)) = 1;
           case 17
          set( handles.t5,'value' , 1);
          set( handles.t1,'value' , 1);
          handles.atoms_selected(handles.all_couplings(i,1)) = 1;
          handles.atoms_selected(handles.all_couplings(i,2)) = 1;
           case 18
          set( handles.t5,'value' , 1);
          set( handles.b1,'value' , 1);
          handles.atoms_selected(handles.all_couplings(i,1)) = 1;
          handles.atoms_selected(handles.all_couplings(i,2)) = 1;
           case 19
          set( handles.b5,'value' , 1);
          set( handles.b1,'value' , 1);
          handles.atoms_selected(handles.all_couplings(i,1)) = 1;
          handles.atoms_selected(handles.all_couplings(i,2)) = 1;
           case 20
          set( handles.b5,'value' , 1);
          set( handles.t1,'value' , 1);
          handles.atoms_selected(handles.all_couplings(i,1)) = 1;
          handles.atoms_selected(handles.all_couplings(i,2)) = 1;
       end
       
     handles.atoms_selected_value = 2;
     end
 end
end

guidata(hObject, handles);

% --------------------------------------------------------------------
function update_constants(hObject, eventdata, handles)

contents = get(handles.couplings, 'Value' );
num = size(contents);
num = num(2);

if num == 1

iii = 0;
aaa = 0;
for ii = 1:20
 if handles.all_couplings(ii,4) == 1
    iii = iii +1 ;
     if iii == contents
       aaa = ii;
    end
 end
end
if aaa > 0
ss = [];
for ii = 5:9
sss = num2str(handles.all_couplings(aaa,ii), '%05.2f\n') ;
ss = vertcat( ss , sss);
end

set(handles.constants, 'String' , ss );
set(handles.update_en, 'Visible' , 'on');
set(handles.update_ot, 'Visible' , 'on');
set(handles.en, 'Visible' , 'on');
set(handles.ot, 'Visible' , 'on');
set(handles.electron, 'Visible' , 'on');
set(handles.other, 'Visible' , 'on');
set(handles.text7, 'Visible' , 'on');
set(handles.text8, 'Visible' , 'on');
set(handles.text9, 'Visible' , 'on');
set(handles.constant, 'Visible' , 'on');
set(handles.constants, 'Visible' , 'on');
set(handles.update_constant, 'Visible' , 'on');
set(handles.text10, 'Visible' , 'off');
else 
%single selected 

sumnum = size(handles.sums);
sumnum = sumnum(1);
sums = contents - iii;

if sums <= sumnum

ss = [];
for ii = 21:25
sss = num2str(handles.sums(sums,ii), '%05.2f\n') ;
ss = vertcat( ss , sss);    
end

set(handles.constants, 'String' , ss );
set(handles.constants, 'Value' , 1 );

ss =[ 'This sum of couplings includes:'];
for ii = 1:20
   if  handles.sums(sums,ii) == 1
     if ii < 10
         iii = '0';
     else
         iii = '';
     end
     ss = vertcat( ss , [' - Coupling ' iii num2str(ii) '.                ']);
   end
end
set(handles.text10, 'String' , ss );

set(handles.update_en, 'Visible' , 'off');
set(handles.update_ot, 'Visible' , 'off');
set(handles.en, 'Visible' , 'off');
set(handles.ot, 'Visible' , 'off');
set(handles.electron, 'Visible' , 'off');
set(handles.other, 'Visible' , 'off');
set(handles.text7, 'Visible' , 'on');
set(handles.text8, 'Visible' , 'off');
set(handles.text9, 'Visible' , 'off');
set(handles.constant, 'Visible' , 'on');
set(handles.constants, 'Visible' , 'on');
set(handles.update_constant, 'Visible' , 'on');
set(handles.text10, 'Visible' , 'on');
else

set(handles.constants, 'Value' , 1 );
set(handles.constants, 'String' , [] );   

set(handles.update_en, 'Visible' , 'off');
set(handles.update_ot, 'Visible' , 'off');
set(handles.en, 'Visible' , 'off');
set(handles.ot, 'Visible' , 'off');
set(handles.electron, 'Visible' , 'off');
set(handles.other, 'Visible' , 'off');
set(handles.text7, 'Visible' , 'off');
set(handles.text8, 'Visible' , 'off');
set(handles.text9, 'Visible' , 'off');
set(handles.constant, 'Visible' , 'off');
set(handles.constants, 'Visible' , 'off');
set(handles.update_constant, 'Visible' , 'off');
set(handles.text10, 'Visible' , 'off');
end
        
end

guidata(hObject, handles);

else

% multiple selected    
set(handles.constants, 'Value' , 1 );
set(handles.constants, 'String' , [] );   

set(handles.update_en, 'Visible' , 'off');
set(handles.update_ot, 'Visible' , 'off');
set(handles.en, 'Visible' , 'off');
set(handles.ot, 'Visible' , 'off');
set(handles.electron, 'Visible' , 'off');
set(handles.other, 'Visible' , 'off');
set(handles.text7, 'Visible' , 'off');
set(handles.text8, 'Visible' , 'off');
set(handles.text9, 'Visible' , 'off');
set(handles.constant, 'Visible' , 'off');
set(handles.constants, 'Visible' , 'off');
set(handles.update_constant, 'Visible' , 'off');
set(handles.text10, 'Visible' , 'off');

end
% --------------------------------------------------------------------
function update_en(hObject, eventdata, handles)

contents = get(handles.couplings, 'Value' );
iii = 0;
aaa = 0;
for ii = 1:20
 if handles.all_couplings(ii,4) == 1
    iii = iii +1 ;
     if iii == contents
       aaa = ii;
    end
 end
end
if aaa > 0
ss = [];
for ii = 10:13
enn = handles.all_couplings(aaa,ii);

sss = num2str(handles.electronegativities(enn), '%05.2f\n') ;
ss = vertcat( ss , sss);
end

set(handles.electron, 'String' , ss );

else 

set(handles.electron, 'Value' , 1 );
set(handles.electron, 'String' , [] );   

end

if aaa > 0
ss = num2str(handles.all_couplings(aaa,14), '%07.2f\n') ;
sss = num2str(handles.all_couplings(aaa,3), '%07.2f\n') ;
ss = vertcat( ss , sss);
sss = num2str(handles.all_couplings(aaa,15), '%07.2f\n') ;
ss = vertcat( ss , sss);
sss = num2str(handles.all_couplings(aaa,16), '%07.2f\n') ;
ss = vertcat( ss , sss);

set(handles.other, 'String' , ss );

else 

set(handles.other, 'Value' , 1 );
set(handles.other, 'String' , [] );   

end

guidata(hObject, handles);

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

% Hint: listbox 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


% --- Executes on selection change in constants.
function constants_Callback(hObject, eventdata, handles)
% hObject    handle to constants (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = get(hObject,'String') returns constants contents as cell array
%        contents{get(hObject,'Value')} returns selected item from constants

contents = get(hObject, 'Value' );
if contents >= 1 
else   
 set(hObject, 'Value', 1 )   ;
end

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

% Hint: listbox 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

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

if handles.atoms_selected_value == 2
a = [ 0 0 ];
aa = 1;
for ii = 1:12
    if handles.atoms_selected(ii) == 1
      a(aa) = ii;
      aa = aa + 1;
    end
end

aa = 0;
for ii = 1:20
   if a == handles.all_couplings(ii,1:2)
       aa = ii;
   end
end

if aa ~= 0
 handles.all_couplings(aa,4) = 1;
 guidata(hObject, handles);
 
 
end

else
if handles.atoms_selected_value == 0
    
 contents = get(handles.couplings, 'Value' );
 num = size(contents);
 num = num(2);
 if num > 1
 
     newentry = zeros(1,25);
     for j = 1:num
     q = 0;
     for i = 1:20
       if handles.all_couplings(i,4) == 1
          q = q + 1;
          if q == contents(j)
              newentry(i) = 1;
          end
       end
     end
     end
 handles.sums = vertcat( handles.sums , newentry );
 end
end
end

recalc_cpl(hObject, eventdata, handles)
guidata(hObject, handles);
update_constants(hObject, eventdata, handles);
update_en(hObject, eventdata, handles);

function recalc_cpl(hObject, eventdata, handles)

 ss = [];
 for ii = 1:20
   if handles.all_couplings(ii,4) == 1
       sss = '';
       if ii < 10
          sss = [sss  '0']; 
       end
       sss = [sss  num2str(ii) '. Coupling '];
       if handles.all_couplings(ii,1) < 7
          sss = [ sss 'T' num2str(handles.all_couplings(ii,1)) ' - '];
       else
          sss = [ sss 'B' num2str(handles.all_couplings(ii,1)-6) ' - ']; 
       end
       if handles.all_couplings(ii,2) < 7
          sss = [ sss 'T' num2str(handles.all_couplings(ii,2)) ];
       else
          sss = [ sss 'B' num2str(handles.all_couplings(ii,2)-6) ]; 
       end       
       ss = vertcat(ss, sss);
   end
 end

 sums = size(handles.sums);
 sums = sums(1);
 
 for ii = 1:sums
     num = 20 + ii;
     sss = [  num2str(num) '. Sum of couplings'];
     ss = vertcat(ss, sss);     
 end
 
 set(handles.couplings, 'Value' , 1);
 set(handles.couplings, 'String' , ss );
 set(handles.constants, 'Value' , 1);
 set(handles.constants, 'String' , [] );

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

contents = get(handles.couplings, 'Value' );
num = size(contents);
num = num(2);

for j = num:-1:1
iii = 0;
del = 0;
for ii = 1:20
 if handles.all_couplings(ii,4) == 1
    iii = iii +1 ;
    if iii == contents(j)
       del = size(handles.sums);
       if del(1) ~= 1
        del = sum(handles.sums);
       else
         del = handles.sums;
       end
       if del(ii) == 0
         handles.all_couplings(ii,4)  = 0;
       else
         fprintf('Cannot remove coupling %2i : Involved in Sum \n',ii)
       end
       del = 1;
    end
 end
end

if del == 0
    del = contents(j) - iii ;
    numsum = size(handles.sums);
    numsum = numsum(1);
    if del <= numsum
     handles.sums(del,:) = [];
    end
    
end
end

guidata(hObject, handles);
recalc_cpl(hObject, eventdata, handles);
update_constants(hObject, eventdata, handles);
update_en(hObject, eventdata, handles);

function atom_select(hObject, eventdata, handles, iii)

selected = get(hObject,'Value');
 if (selected == 1 )
     handles.atoms_selected_value = handles.atoms_selected_value + 1;
     handles.atoms_selected(iii) = 1;    
  else
     handles.atoms_selected_value = handles.atoms_selected_value - 1;
     handles.atoms_selected(iii) = 0;    
 end

guidata(hObject, handles);


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


% --- Executes during object creation, after setting all properties.
function constant_CreateFcn(hObject, eventdata, handles)
% hObject    handle to constant (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


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

coupling     = get(handles.couplings, 'Value');
temperature  = get(handles.constants, 'Value');
cnst         = get(handles.constant , 'String');
iii = 0;

try
    upd = 0;
    cnst = str2num(cnst);
    for ii = 1:20
    if handles.all_couplings(ii,4) == 1 
      iii = iii + 1;
      if iii == coupling
          upd = 1;
          handles.all_couplings(ii,4+temperature) = cnst;
      end
    end
    end
    if upd == 0 
        handles.sums(coupling - iii , 20 + temperature ) = cnst;
    end 
        
catch
end

guidata(hObject, handles);
update_constants(hObject, eventdata, handles);

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

try
 handles.parameters(1) = str2num( get(handles.Pa_value , 'String') ) ;
catch  
end
try
 handles.parameters(2) = str2num( get(handles.phia_value , 'String')) ;
catch  
end
try
 handles.parameters(4) = str2num( get(handles.Pb_value , 'String')) ;
catch  
end
try
 handles.parameters(5) = str2num( get(handles.phib_value , 'String')) ;
catch  
end
try
 handles.parameters(7) = str2num( get(handles.partition_value , 'String')) ;
catch  
end
handles.parameters(8)  = get(handles.Pa_check , 'Value' ) ;
handles.parameters(9)  = get(handles.phia_check , 'Value' );
handles.parameters(11) = get(handles.Pb_check , 'Value') ;
handles.parameters(12) = get(handles.phib_check , 'Value') ;
handles.parameters(14) = get(handles.partition_check , 'Value') ;
handles.parameters(15) = get(handles.temp , 'Value') ;

guidata(hObject, handles);

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


set(handles.Pa_value     , 'String' , num2str( handles.parameters(1) , '%5.2d')  ) ;
set(handles.phia_value     , 'String' , num2str( handles.parameters(2) , '%5.2d')  ) ;
set(handles.Pb_value     , 'String' , num2str( handles.parameters(4) , '%5.2d')  ) ;
set(handles.phib_value     , 'String' , num2str( handles.parameters(5) , '%5.2d')  ) ;
set(handles.partition_value , 'String' , num2str( handles.parameters(7) , '%5.2d')  ) ;
set(handles.Pa_check     , 'Value' , handles.parameters(8)  ) ;
set(handles.phia_check     , 'Value' , handles.parameters(9)  );
set(handles.Pb_check     , 'Value' , handles.parameters(11) ) ;
set(handles.phib_check     , 'Value' , handles.parameters(12) ) ;
set(handles.partition_check , 'Value' , handles.parameters(14) ) ;
set(handles.temp            , 'Value' , handles.parameters(15) ) ;

guidata(hObject, handles);


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


% --- Executes on button press in Pa_check.
function Pa_check_Callback(hObject, eventdata, handles)
% hObject    handle to Pa_check (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 Pa_check


% --- Executes on button press in phia_check.
function phia_check_Callback(hObject, eventdata, handles)
% hObject    handle to phia_check (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 phia_check


% --- Executes on button press in pa_check.
function pa_check_Callback(hObject, eventdata, handles)
% hObject    handle to pa_check (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 pa_check


% --- Executes on button press in partition_check.
function partition_check_Callback(hObject, eventdata, handles)
% hObject    handle to partition_check (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 partition_check



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


% --- Executes during object creation, after setting all properties.
function Pa_value_CreateFcn(hObject, eventdata, handles)
% hObject    handle to Pa_value (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



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


% --- Executes during object creation, after setting all properties.
function phia_value_CreateFcn(hObject, eventdata, handles)
% hObject    handle to phia_value (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



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


% --- Executes during object creation, after setting all properties.
function pa_value_CreateFcn(hObject, eventdata, handles)
% hObject    handle to pa_value (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


% --- Executes on button press in Pb_check.
function Pb_check_Callback(hObject, eventdata, handles)
% hObject    handle to Pb_check (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 Pb_check


% --- Executes on button press in phib_check.
function phib_check_Callback(hObject, eventdata, handles)
% hObject    handle to phib_check (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 phib_check


% --- Executes on button press in pb_check.
function pb_check_Callback(hObject, eventdata, handles)
% hObject    handle to pb_check (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 pb_check



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


% --- Executes during object creation, after setting all properties.
function Pb_value_CreateFcn(hObject, eventdata, handles)
% hObject    handle to Pb_value (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



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


% --- Executes during object creation, after setting all properties.
function phib_value_CreateFcn(hObject, eventdata, handles)
% hObject    handle to phib_value (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



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


% --- Executes during object creation, after setting all properties.
function pb_value_CreateFcn(hObject, eventdata, handles)
% hObject    handle to pb_value (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



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


% --- Executes during object creation, after setting all properties.
function partition_value_CreateFcn(hObject, eventdata, handles)
% hObject    handle to partition_value (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



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


% --- Executes during object creation, after setting all properties.
function file_path_CreateFcn(hObject, eventdata, handles)
% hObject    handle to file_path (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


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

direct = [num2str(date()) '_' num2str(datenum(rem(now(),1))*24)];
mkdir(handles.dir, direct);
direct = [handles.dir direct '\'];

ppath = [ handles.dir get( handles.file_path , 'String' ) '.par'];
csvwrite(ppath, handles.parameters);
ppath = [ handles.dir get( handles.file_path , 'String' ) '.en'];
csvwrite(ppath, handles.electronegativities);
ppath = [ handles.dir get( handles.file_path , 'String' ) '.cpl'];
csvwrite(ppath, handles.all_couplings);
ppath = [ handles.dir get( handles.file_path , 'String' ) '.sum'];
csvwrite(ppath, handles.sums);

ppath = [ direct get( handles.file_path , 'String' ) '.par'];
csvwrite(ppath, handles.parameters);
ppath = [ direct get( handles.file_path , 'String' ) '.en'];
csvwrite(ppath, handles.electronegativities);
ppath = [ direct get( handles.file_path , 'String' ) '.cpl'];
csvwrite(ppath, handles.all_couplings);
ppath = [ direct get( handles.file_path , 'String' ) '.sum'];
csvwrite(ppath, handles.sums);

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

try
ppath = [ handles.dir get( handles.file_path , 'String' ) '.par'];
handles.parameters = csvread(ppath);
ppath = [ handles.dir get( handles.file_path , 'String' ) '.en'];
handles.electronegativities =  csvread(ppath);
ppath = [ handles.dir get( handles.file_path , 'String' ) '.cpl'];
handles.all_couplings = csvread(ppath);
catch

    fprintf('\n\n-->  ERROR: files not found / Load failed\n');

end

handles.sums = zeros(0,25);
try
ppath = [ handles.dir get( handles.file_path , 'String' ) '.sum'];
handles.sums = csvread(ppath);

catch

    fprintf('\n\n-->  WARNING: .Sum file not found\n');

end

if sum(size(handles.sums)) == 1
   handles.sums = zeros(0,25);
end

guidata(hObject, handles);
recalc_cpl(hObject, eventdata, handles);
update_constants(hObject, eventdata, handles);
update_en(hObject, eventdata, handles);
set_parameters(hObject, eventdata, handles);


% --- Executes on selection change in electron.
function electron_Callback(hObject, eventdata, handles)
% hObject    handle to electron (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = get(hObject,'String') returns electron contents as cell array
%        contents{get(hObject,'Value')} returns selected item from electron

contents = get(hObject, 'Value' );
if contents >= 1 
else    
 set(hObject, 'Value', 1 )   ;
end


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

% Hint: listbox 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



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


% --- Executes during object creation, after setting all properties.
function edit10_CreateFcn(hObject, eventdata, handles)
% hObject    handle to edit10 (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


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

coupling     = get(handles.couplings, 'Value');
elen         = get(handles.electron, 'Value');
cnst         = get(handles.en , 'String');
iii = 0;

try
    cnst = str2num(cnst);
    for ii = 1:20
    if handles.all_couplings(ii,4) == 1 
      iii = iii + 1;
      if iii == coupling
          enn = handles.all_couplings(ii,9+elen) ;
          handles.electronegativities(enn) = cnst;
      end
    end
    end
catch
    
end

guidata(hObject, handles);
update_en(hObject, eventdata, handles);


% --- Executes on selection change in en_select.
function en_select_Callback(hObject, eventdata, handles)
% hObject    handle to en_select (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = get(hObject,'String') returns en_select contents as cell array
%        contents{get(hObject,'Value')} returns selected item from en_select


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

% Hint: popupmenu 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



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


% --- Executes during object creation, after setting all properties.
function en_CreateFcn(hObject, eventdata, handles)
% hObject    handle to en (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




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


% --------------------------------------------------------------------
function calc1_Callback(hObject, eventdata, handles)
  fprintf ( '\n\n////////////////////////////////////////////////////////////// \n'); 
  fprintf ( '////////////////////////////////////////////////////////////// \n\n');  
  fprintf ( '  START Pseudorotational calculation \n');
  fprintf ( '\n////////////////////////////////////////////////////////////// \n');
  fprintf ( '////////////////////////////////////////////////////////////// \n\n');

try
 handles.parameters(1) = str2num( get(handles.Pa_value , 'String') ) ;
catch  
end
try
 handles.parameters(2) = str2num( get(handles.phia_value , 'String')) ;
catch  
end
try
 handles.parameters(4) = str2num( get(handles.Pb_value , 'String')) ;
catch  
end
try
 handles.parameters(5) = str2num( get(handles.phib_value , 'String')) ;
catch  
end
try
 handles.parameters(7) = str2num( get(handles.partition_value , 'String')) ;
catch  
end
handles.parameters(8)  = get(handles.Pa_check , 'Value' ) ;
handles.parameters(9)  = get(handles.phia_check , 'Value' );
handles.parameters(11) = get(handles.Pb_check , 'Value') ;
handles.parameters(12) = get(handles.phib_check , 'Value') ;
handles.parameters(14) = get(handles.partition_check , 'Value') ;
handles.parameters(15) = get(handles.temp , 'Value') ;

guidata(hObject, handles);

calc2_Callback(hObject, eventdata, handles);


% --------------------------------------------------------------------
function calc2_Callback(hObject, eventdata, handles)

elapsed_time = clock;

textout = [];

xxx= [];
yyy= [];
matrix = zeros(7);
parameters = horzcat( handles.parameters(1:6) * pi() / 180 , handles.parameters(7));  
counter  = 1;
counter2 = 1;  
numtemp = handles.parameters(15);

for ii = 1:7
    if handles.parameters(7+ii) == 1
        if ii ~= 7
          xxx(counter)  = parameters(ii);
          counter = counter + 1;
        else
          for iii = 1:numtemp
            xxx(counter)  = parameters(ii)/100;
            counter = counter + 1;
          end
        end
    else
        if ii ~= 7
          yyy(counter2)  = parameters(ii);
          counter2 = counter2 + 1;
        else
          for iii = 1:numtemp
            yyy(counter2)  = parameters(ii)/100;
            counter2 = counter2 + 1;
          end
        end
    end
end

counter2 = 1;   
for ii = 1:7
    if handles.parameters(7+ii) == 1
        if ii ~= 7
          matrix(counter2,ii) = 1;
          counter2 = counter2 + 1;    
        else
          for iii = 1:numtemp
            matrix(counter2,ii + iii - 1) = 1;
            counter2 = counter2 + 1; 
          end
        end
    else
        if ii ~= 7
          matrix(counter,ii) = 1;
          counter = counter + 1;
        else
          for iii = 1:numtemp 
            matrix(counter,ii + iii - 1) = 1;
            counter = counter + 1;              
          end
        end
    end
end

counter = 1;
cpl = [];
for ii = 1:20
   if handles.all_couplings(ii,4) == 1     
      cpl = vertcat ( cpl , handles.all_couplings(ii,:) );  
      cpl ( counter , 10 ) = handles.electronegativities( cpl ( counter , 10 ) );
      cpl ( counter , 11 ) = handles.electronegativities( cpl ( counter , 11 ) );
      cpl ( counter , 12 ) = handles.electronegativities( cpl ( counter , 12 ) );
      cpl ( counter , 13 ) = handles.electronegativities( cpl ( counter , 13 ) );    
      cpl ( counter ,  4 ) = (ii - 1 - mod(ii - 1,4) + 4 ) /4 ;
      counter = counter + 1 ;
   end 
end

if size(cpl) > 0
    
cpl = horzcat ( cpl(:,5:5+numtemp-1) , cpl(:,10:13) , cpl(:,3)/180*pi() , cpl(:,4) , cpl(:,14:15), cpl(:,16)/180*pi() );
sums = handles.sums;
numsum = size(sums);
numsum = numsum(1);

if numsum > 0
for i = 20:-1:1
   if handles.all_couplings(i,4) == 0
       sums(:,i) = [];
   end
end
end

fff = @(var) search_minimum(var,yyy,matrix,cpl , numtemp , 0 , handles.karplus, sums);
options = optimset('MaxIter', 250, 'GradObj', 'off');


  
if size(xxx) > 0
  qq = size(xxx);
  qq = qq(2);
if handles.parameters(14) == 1  
  
  A = zeros(numtemp*2 , qq);
  B = zeros(numtemp*2 , 1);
  counter=1;
  for i = 1:numtemp
   A( counter     , qq - numtemp + i ) = 1;
   A( counter + 1 , qq - numtemp + i ) = -1;
   B( counter     , 1 ) = 1;
   B( counter + 1 , 1 ) = 0;
   counter = counter + 2 ;
  end

else
    A = zeros(1, qq);
    B = [0];
end

QQ = [];
QQQ = [];
jj = 0;
for ii = 1:6
    if handles.parameters(7+ii) == 1
         jj = jj + 1;
         CC = zeros(2,qq);
         CC(1:2,jj) = [1 , -1];
         DD = [2*pi() , 2*pi()]' ;
         if ii == 2
             DD = [pi()/3 , pi()/3]' ;
         end
         if ii == 5
             DD = [pi()/3 , pi()/3]' ;
         end
         QQ  = vertcat( QQ  , CC );
         QQQ = vertcat( QQQ , DD );
    end
end

A  = vertcat( QQ   , A );
B  = vertcat( QQQ  , B );

  [optimised_parameters, total_error] = fmincon(fff,xxx, A,B, [],[],[],[],[],options);
  
  opt_parameters = horzcat( optimised_parameters , yyy ) * matrix;
  opt_parameters(1:6) = opt_parameters(1:6) * 180 / pi();
  opt_parameters(7:6+numtemp) = opt_parameters(7:6+numtemp) * 100;

  if opt_parameters(2) < 0
     opt_parameters(2) = - opt_parameters(2);
     opt_parameters(1) = opt_parameters(1) + 180;
  end
  if opt_parameters(5) < 0
     opt_parameters(5) = - opt_parameters(5);
     opt_parameters(4) = opt_parameters(4) + 180;
  end
  textout = horzcat(textout, sprintf ( '  Optimized parameters \n'));
  textout = horzcat(textout, sprintf ( '-------------------------------------------------------------- \n\n'));
  textout = horzcat(textout, sprintf ( '   Conformation 1      \n') ); 
  textout = horzcat(textout, sprintf ( '    P     : %8.3f      \n' ,opt_parameters(1)) ); 
  textout = horzcat(textout, sprintf ( '    Phi   : %8.3f      \n\n' ,opt_parameters(2)));
  textout = horzcat(textout, sprintf ( '   Conformation 2      \n'));
  textout = horzcat(textout, sprintf ( '    P     : %8.3f      \n' ,opt_parameters(4)));
  textout = horzcat(textout, sprintf ( '    Phi   : %8.3f      \n\n' ,opt_parameters(5)));
  textout = horzcat(textout, sprintf ( '   Temperature Coefficients \n'));  
    
  fprintf ( '\n\n-------------------------------------------------------------- \n');  
  fprintf ( '  Optimized parameters \n');
  fprintf ( '-------------------------------------------------------------- \n\n');
  fprintf ( '   Conformation 1      \n');
  fprintf ( '    P     : %8.3f      \n' ,opt_parameters(1));
  fprintf ( '    Phi_m : %8.3f      \n\n' ,opt_parameters(2));
  fprintf ( '   Conformation 2      \n');
  fprintf ( '    P     : %8.3f      \n' ,opt_parameters(4));
  fprintf ( '    Phi_m : %8.3f      \n\n' ,opt_parameters(5));
  fprintf ( '   Temperature Coefficients \n');  
  
  for i = 7:6+numtemp
    fprintf ( '    %%Conformation1  : %8.3f\n' ,opt_parameters(i));
    textout = horzcat(textout, sprintf ( '    %%Conformation1  : %8.3f\n' ,opt_parameters(i)));  
  end
  
  
  total_error = search_minimum(optimised_parameters,yyy,matrix,cpl , numtemp , 1 , handles.karplus, sums) ;

  opt_parameters(3) = [];
  opt_parameters(5) = [];
  
  assignin('base','optimized_parameters', opt_parameters)  ; 
%     % start graph 

      options = optimset('MaxIter', 250, 'Display', 'off', 'GradObj', 'off');

      parnum = size(xxx);
      parnum = parnum(2);
      parnum2 = parnum;

      if handles.parameters(14) == 1
       parnum2 = parnum - numtemp;
      end
      ncpl = size(cpl);  
      ncpl = ncpl(1);
      fin = [];
      for i = 1:parnum
         unknown = optimised_parameters;

         conA = A;
         conB = B;


         if i <= parnum2
            unknown(i) = unknown(i) + 3 / 180 * pi();
            dum = 3 / 180 * pi();
         else
            unknown(i) = unknown(i) + 0.05;
            dum = 0.05;
         end

         conAA = zeros(2,parnum);
         conAA(1,i) = 1;
         conAA(2,i) = -1;

         conA = vertcat(conA , conAA );
         conB = vertcat(conB , [ unknown(i) -unknown(i) ]' );

         [opar, te] = fmincon(fff,unknown, conA, conB, [],[],[],[],[],options);

         final = horzcat( (opar - optimised_parameters)/ dum , te - total_error);
         fin = vertcat(fin, final);   
      end

    fprintf ( '\n-------------------------------------------------------------- \n');
    fprintf ( '  ERROR ANALYSIS\n');
    fprintf ( '-------------------------------------------------------------- \n\n');

    fprintf ( '  i    ');
    for i = 1 : parnum 
        fprintf ( '   dP%1i/dPi  ', i);
    end    
    fprintf ( '  dRMSD/dPi ', i);
    fprintf ( '\n');
    for i = 1 : parnum
        fprintf ( ' %2i  ' , i );
        fprintf ( ' %10.3f ', fin(i,:));
        fprintf ( '\n');
    end

  assignin('base','err_analysis_parameters', fin)  ; 

% begin error analysis couplings  
  
    cpl_org = cpl;
    fin = [];
    
    for i = 1 : ncpl
        
     cpl = cpl_org;
     
     for ii = 1 : numtemp
         if cpl(i,ii) > 0
           cpl(i,ii) = cpl(i,ii) + 0.5;
         end
     end
     
     fff = @(var) search_minimum(var,yyy,matrix,cpl , numtemp , 0 , handles.karplus, sums);
     [opar, te] = fmincon(fff,xxx, A,B, [],[],[],[],[],options);
     
     final = horzcat( (opar - optimised_parameters)/ 0.5 , te - total_error);

     for ii = 1:parnum
         if ii <= parnum2
            final(ii) = final(ii) * 180 / pi();    
         else
            final(ii) = final(ii) * 100 ;
         end
     end
     
     fin = vertcat(fin, final);       
    
    end

    fprintf ( '\n  i    ');
    for i = 1 : parnum 
        fprintf ( '   dP%1i/dCi  ', i);
    end    
    fprintf ( '  dRMSD/dCi ', i);
    fprintf ( '\n');
    for i = 1 : ncpl
        fprintf ( ' %2i  ' , i );
        fprintf ( ' %10.3f ', fin(i,:));
        fprintf ( '\n');
    end

    assignin('base','err_analysis_couplings', fin)  ;
    cpl = cpl_org;    
 
 % end error analysis couplings  
 
 % begin error analysis sums of couplings  
  
    sums_org = sums;
    fin = [];
    if numsum > 0
    
    for i = 1 : numsum
        
     sums = sums_org;
     
     for ii = 1 : numtemp
         if sums(i,ii + ncpl) > 0
           sums(i,ii + ncpl) = sums(i,ii + ncpl) + 0.5;
         end
     end
     
     fff = @(var) search_minimum(var,yyy,matrix,cpl , numtemp , 0 , handles.karplus, sums);
     [opar, te] = fmincon(fff,xxx, A,B, [],[],[],[],[],options);
     
     final = horzcat( (opar - optimised_parameters)/ 0.5 , te - total_error);

     for ii = 1:parnum
         if ii <= parnum2
            final(ii) = final(ii) * 180 / pi();    
         else
            final(ii) = final(ii) * 100 ;
         end
     end
     
     fin = vertcat(fin, final);       
    
    end

    fprintf ( '\n  i    ');
    for i = 1 : parnum 
        fprintf ( '   dP%1i/dSi  ', i);
    end    
    fprintf ( '  dRMSD/dSi ', i);
    fprintf ( '\n');
    for i = 1 : numsum
        fprintf ( ' %2i  ' , i );
        fprintf ( ' %10.3f ', fin(i,:));
        fprintf ( '\n');
    end

    assignin('base','err_analysis_sums', fin)  ;
    sums = sums_org;    
 
    end
 % end error analysis sums of couplings 
 
else
  total_error = search_minimum(xxx,yyy,matrix,cpl , numtemp , 1 , handles.karplus, sums) ;
  
  yyy = yyy * 100;
  yyy(1:6) = yyy(1:6) * 1.8 / pi();
    
  textout = horzcat(textout, sprintf ( '  Input parameters \n'));
  textout = horzcat(textout, sprintf ( '-------------------------------------------------------------- \n\n'));
  textout = horzcat(textout, sprintf ( '   Conformation 1      \n') ); 
  textout = horzcat(textout, sprintf ( '    P     : %8.3f      \n' ,yyy(1)) ); 
  textout = horzcat(textout, sprintf ( '    Phi   : %8.3f      \n\n' ,yyy(2)));
  textout = horzcat(textout, sprintf ( '   Conformation 2      \n'));
  textout = horzcat(textout, sprintf ( '    P     : %8.3f      \n' ,yyy(4)));
  textout = horzcat(textout, sprintf ( '    Phi   : %8.3f      \n\n' ,yyy(5)));
  textout = horzcat(textout, sprintf ( '   Temperature Coefficients \n'));  
    
  fprintf ( '\n\n-------------------------------------------------------------- \n');  
  fprintf ( '  Input parameters \n');
  fprintf ( '-------------------------------------------------------------- \n\n');
  fprintf ( '   Conformation 1      \n');
  fprintf ( '    P     : %8.3f      \n' ,yyy(1));
  fprintf ( '    Phi_m : %8.3f      \n\n' ,yyy(2));
  fprintf ( '   Conformation 2      \n');
  fprintf ( '    P     : %8.3f      \n' ,yyy(4));
  fprintf ( '    Phi_m : %8.3f      \n\n' ,yyy(5));
  fprintf ( '   Temperature Coefficients \n');  
  
  for i = 7:6+numtemp
    fprintf ( '    %%Conformation1  : %8.3f\n' ,yyy(i));
    textout = horzcat(textout, sprintf ( '    %%Conformation1  : %8.3f\n' ,yyy(i)));  
  end


end

numsum = size(sums);
numsum = numsum(1);
ncpl = size(cpl);  
ncpl = ncpl(1);

    fprintf ( '\n-------------------------------------------------------------- \n');
    fprintf ( '  RMSD                : %10.2f Hz\n' , total_error);
    textout = horzcat(textout, sprintf ( '\n  RMSD                : %10.2f Hz\n' , total_error));
    
elapsed_time = - etime(elapsed_time , clock);
  fprintf ( '  Total time          : %10.2f s\n' , elapsed_time);

    fprintf ( '-------------------------------------------------------------- \n\n');
assignin('base','total_error', total_error)  ; 
else
   fprintf ( '\n-------------------------------------------------------------- \n');
   fprintf ( '  No couplings !!!\n');
   fprintf ( '-------------------------------------------------------------- \n\n\n'); 
end    
    
assignin('base','textout1',textout);
guidata(hObject, handles);

% --------------------------------------------------------------------
function error = search_minimum ( unk1 , kno1 , link , cpl_data , nt , print , keq, sums)

ommit = [ 0 0 0 0 0 ];

parameters = horzcat( unk1 , kno1 ) * link;
ncpl = size(cpl_data);  
cols = ncpl(2);
ncpl = ncpl(1);

for i = 1:ncpl
   lambda =  cpl_data(i, nt+1:nt+4);
   
   for ii = 1:5 
   phi(ii,1) = parameters(2) * cpl_data(i,cols - 1)* cos( parameters(1) + cpl_data(i,cols) + 4/5 * pi() * (ii-3) ) ;
   phi(ii,2) = parameters(5) * cpl_data(i,cols - 1)* cos( parameters(4) + cpl_data(i,cols) + 4/5 * pi() * (ii-3) ) ;
   end
   
   switch keq
       case 1 
         c0 = 6.97 - 0.58 * sum(lambda) - 0.24 * ( lambda(1)*lambda(2) + lambda(3)*lambda(4) ) ;
         c1 = -1.06 ;
         c2 = 6.55 - 0.82 * sum(lambda) + 0.20 * ( lambda(1)*lambda(4) + lambda(3)*lambda(2) ) ;
         c3 = -0.54 ;
         s1 = 0;
         s2 = 0.68 * ( lambda(1)^2 - lambda(2)^2 + lambda(3)^2 - lambda(4)^2 ) ;
       case 2
         c0 =  6.675 - 0.376 * sum(lambda) - 0.133 * ( lambda(1)*lambda(2) + lambda(3)*lambda(4) ) ;
         c1 = -0.976 + 0.096 * sum(lambda) ;
         c2 =  6.278 - 0.572 * sum(lambda) + 0.131 * ( lambda(1)*lambda(4) + lambda(3)*lambda(2) );
         c3 = -0.147 ;
         s1 = -0.073 * ( lambda(1) - lambda(2) + lambda(3) - lambda(4) );
         s2 =  0.869 * ( lambda(1) - lambda(2) + lambda(3) - lambda(4) ) - 0.148 * ( lambda(1)*lambda(3) - lambda(2)*lambda(4) ) ;      
       case 3
         c0 =  6.736 - 0.422 * sum(lambda) - 0.133 * ( lambda(1)*lambda(2) + lambda(3)*lambda(4) ) + 0.038 * ( lambda(1)*lambda(3) + lambda(2)*lambda(4) ) + 0.040 * ( lambda(1)*lambda(4) + lambda(3)*lambda(2) );
         c1 = -0.893 + 0.039 * sum(lambda) + 0.047 * ( lambda(1)*lambda(2) + lambda(3)*lambda(4) ) + 0.027 * ( lambda(1)*lambda(3) + lambda(2)*lambda(4) ) - 0.003 * ( lambda(1)*lambda(4) + lambda(3)*lambda(2) );
         c2 =  6.144 - 0.481 * sum(lambda) - 0.070 * ( lambda(1)*lambda(2) + lambda(3)*lambda(4) ) - 0.062 * ( lambda(1)*lambda(3) + lambda(2)*lambda(4) ) + 0.150 * ( lambda(1)*lambda(4) + lambda(3)*lambda(2) );
         c3 = -0.147 ;
         s1 = -0.077 * ( lambda(1) - lambda(2) + lambda(3) - lambda(4) ) + 0.006 * ( lambda(1)*lambda(3) - lambda(2)*lambda(4) ) ;
         s2 =  0.865 * ( lambda(1) - lambda(2) + lambda(3) - lambda(4) ) - 0.146 * ( lambda(1)*lambda(3) - lambda(2)*lambda(4) ) ;      
       case 4 
         c0 = 7.01 - 0.58 * ( lambda(1)+lambda(2)+lambda(3)+lambda(4) ) - 0.24 * ( (lambda(1)*lambda(2)) + (lambda(3)*lambda(4)) ) ;
         c1 = -1.08 ;
         c2 = 6.54 - 0.82 * ( lambda(1)+lambda(2)+lambda(3)+lambda(4) ) + 0.20 * ( (lambda(1)*lambda(4)) + (lambda(3)*lambda(2)) ) ;
         c3 = -0.49 ;
         s1 = 0;
         s2 = 0.68 * ( lambda(1)^2 - lambda(2)^2 + lambda(3)^2 - lambda(4)^2 ) ;
   end

   phi1 = cpl_data(i,cols - 2) * phi( cpl_data(i,cols - 3) , 1) + cpl_data(i,cols - 4);
   phi2 = cpl_data(i,cols - 2) * phi( cpl_data(i,cols - 3) , 2) + cpl_data(i,cols - 4);
   cpl(i,1) = c0 + c1 * cos( phi1 ) + c2 * cos( 2 * phi1 ) + c3 * cos ( 3 * phi1 ) + s1 * sin ( phi1 ) + s2 * sin ( 2 * phi1 );
   cpl(i,2) = c0 + c1 * cos( phi2 ) + c2 * cos( 2 * phi2 ) + c3 * cos ( 3 * phi2 ) + s1 * sin ( phi2 ) + s2 * sin ( 2 * phi2 ); 
end

for i = 1:nt
   weighted(:,i) = parameters(6+i) *  cpl(:,1) + (1-parameters(6+i)) *  cpl(:,2);
end

sw = weighted;

for i = 1:ncpl
    for ii = 1:nt
        if cpl_data(i,ii) <= 0
            ommit(ii) = ommit(ii) + 1;
            weighted(i,ii) = cpl_data(i,ii);
        end
    end
end

numsum = size(sums);
numsum = numsum(1);

sumerror = 0;
if numsum > 0
   sumerror = zeros(numsum,nt);
   for i = 1:numsum
       for ii = 1:nt       
         if sums(i, ncpl + ii) <=0
             ommit(ii) = ommit(ii) + 1;
             sumerror(i,ii) = 0;
         else
             dummy = sum( sums(i,1:ncpl) .* sw(:,ii)' );
             sumerror(i,ii) = dummy - sums(i,ii+ncpl) ;
         end
       end
   end
end


error = weighted - cpl_data(:,1:nt) ;

if print == 1

textout = [];
    
textout = horzcat(textout,sprintf ( ' \n'));    
textout = horzcat(textout,sprintf ( '  Endocyclic torsion angles\n' ));    
textout = horzcat(textout,sprintf ( '-------------------------------------------------------------- \n\n\n'));
textout = horzcat(textout,sprintf ( '            Conf.1       Conf.2 \n' ));
textout = horzcat(textout,sprintf ( '   Phi0: %10.3f   %10.3f\n' , phi(3,1) * 180 / pi(), phi(3,2)* 180 / pi() ));
textout = horzcat(textout,sprintf ( '   Phi1: %10.3f   %10.3f\n' , phi(4,1) * 180 / pi(), phi(4,2)* 180 / pi() ));
textout = horzcat(textout,sprintf ( '   Phi2: %10.3f   %10.3f\n' , phi(5,1) * 180 / pi(), phi(5,2)* 180 / pi() ));
textout = horzcat(textout,sprintf ( '   Phi3: %10.3f   %10.3f\n' , phi(1,1) * 180 / pi(), phi(1,2)* 180 / pi() ));
textout = horzcat(textout,sprintf ( '   Phi4: %10.3f   %10.3f\n\n\n' , phi(2,1) * 180 / pi(), phi(2,2)* 180 / pi() ));
    
fprintf ( ' \n');    
fprintf ( '-------------------------------------------------------------- \n');
fprintf ( '  Endocyclic torsion angles\n' );    
fprintf ( '-------------------------------------------------------------- \n\n\n');
fprintf ( '            Conf.1       Conf.2 \n' );
fprintf ( '   Phi0: %10.3f   %10.3f\n' , phi(3,1) * 180 / pi(), phi(3,2)* 180 / pi() );
fprintf ( '   Phi1: %10.3f   %10.3f\n' , phi(4,1) * 180 / pi(), phi(4,2)* 180 / pi() );
fprintf ( '   Phi2: %10.3f   %10.3f\n' , phi(5,1) * 180 / pi(), phi(5,2)* 180 / pi() );
fprintf ( '   Phi3: %10.3f   %10.3f\n' , phi(1,1) * 180 / pi(), phi(1,2)* 180 / pi() );
fprintf ( '   Phi4: %10.3f   %10.3f\n\n\n' , phi(2,1) * 180 / pi(), phi(2,2)* 180 / pi() );

  assignin('base','endo_torsion_angles', vertcat(phi(3:5,:) * 180 / pi(),phi(1:2,:) * 180 / pi()))  ; 

textout = horzcat(textout,sprintf ( '  Final couplings\n' ));    
textout = horzcat(textout,sprintf ( '-------------------------------------------------------------- \n\n'));
  
fprintf ( '-------------------------------------------------------------- \n');
fprintf ( '  Final couplings\n' );    
fprintf ( '-------------------------------------------------------------- \n\n');

for i = 1:nt
   rmsd = 0;
   
   if i == 1
   textout = horzcat(textout,sprintf ( '   Temperature %1i:\n\n' , i ));
   textout = horzcat(textout,sprintf ( '    ------------------------------------\n'));
   textout = horzcat(textout,sprintf ( '    Conf1  Conf2    Avg.   Exp.  Diff.\n'));
   textout = horzcat(textout,sprintf ( '    ------------------------------------\n'));
   for ii = 1:ncpl
     textout = horzcat(textout, sprintf ( '    %5.2f  %5.2f  %5.2f  %5.2f  %6.2f \n' , cpl(ii,1) , cpl(ii,2) , weighted(ii,i) , cpl_data(ii,i) ,  weighted(ii,i) - cpl_data(ii,i)));
   end
   if numsum > 0
     textout = horzcat(textout, sprintf ( '    ------------------------------------\n'));
      for ii = 1:numsum
          textout = horzcat(textout, sprintf('     Sum %2i :     %5.2f  %5.2f  %6.2f\n',ii, sum( sums(ii,1:ncpl) .* sw(:,i)' ), sums(ii, i+ncpl), sumerror(ii,i) ));
      end
   end
   end
   
   fprintf ( '   Temperature %1i:\n\n' , i );
   fprintf ( '    ------------------------------------\n');
   fprintf ( '    Conf1  Conf2    Avg.   Exp.  Diff.\n');
   fprintf ( '    ------------------------------------\n');
   
   for ii = 1:ncpl
       
      fprintf ( '    %5.2f  %5.2f  %5.2f  %5.2f  %6.2f \n' , cpl(ii,1) , cpl(ii,2) , weighted(ii,i) , cpl_data(ii,i) ,  weighted(ii,i) - cpl_data(ii,i));
      rmsd = rmsd + ( weighted(ii,i) - cpl_data(ii,i) )^2; 
   end
   
   if numsum > 0
      fprintf ( '    ------------------------------------\n');
      for ii = 1:numsum
            fprintf('     Sum %2i :     %5.2f  %5.2f  %6.2f\n',ii, sum( sums(ii,1:ncpl) .* sw(:,i)' ), sums(ii, i+ncpl), sumerror(ii,i) );
            rmsd = rmsd + sumerror(ii,i)^2;
      end
   end
   
   rmsd = sqrt(rmsd/(ncpl+numsum-ommit(i)));
   
   if i == 1
   textout = horzcat(textout, sprintf ( '    ------------------------------------\n'));
   textout = horzcat(textout, sprintf ( '    RMSD : %8.2f Hz\n\n' , rmsd));  
   end
   
   fprintf ( '    ------------------------------------\n');
   fprintf ( '    RMSD : %8.2f Hz\n\n' , rmsd);    

   assignin('base','textout2', textout);
   
end

weighted_couplings =  weighted;
error_NMR_calc = error  ;

end

totommit = sum(ommit);

sumerror = sumerror .* sumerror;
sumerror = sum(sum(sumerror));
error = error .* error;
error = sum(sum(error));

error = error + sumerror;
error = error / ( nt*(ncpl+numsum) - totommit);
error = sqrt(error);

% --- Executes on selection change in other.
function other_Callback(hObject, eventdata, handles)
% hObject    handle to other (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = get(hObject,'String') returns other contents as cell array
%        contents{get(hObject,'Value')} returns selected item from other

contents = get(hObject, 'Value' );
if contents >= 1 
else    
 set(hObject, 'Value', 1 )   ;
end


set(handles.ot,'String',''); 
if contents == 1
   set(handles.ot,'String','endo-exo A');   
end
if contents == 2
   set(handles.ot,'String','endo-exo B');   
end
if contents == 3
   set(handles.ot,'String','Diez a');   
end
if contents == 4
   set(handles.ot,'String','Diez e');   
end



guidata(hObject, handles);

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

% Hint: listbox 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


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

coupling     = get(handles.couplings, 'Value');
oth          = get(handles.other, 'Value');
cnst         = get(handles.ot , 'String');
iii = 0;

try

col = 0;    
if oth == 1
  col = 14;  
end
if oth == 2
  col = 3;  
end
if col > 0
    cnst = str2num(cnst);
    for ii = 1:20
    if handles.all_couplings(ii,4) == 1 
      iii = iii + 1;
      if iii == coupling
          handles.all_couplings(ii,col) = cnst ;
      end
    end
    end
end

if oth == 3
  col = 15;  
end
if oth == 4
  col = 16;  
end    
if col > 14
    cnst = str2num(cnst);
    for ii = 1:20
    if handles.all_couplings(ii,4) == 1 
      iii = iii + 1;
      if iii == coupling
          bondnr = floor((ii - 1)/4) ; 
          for loopnr = 1:4
              handles.all_couplings(bondnr*4+loopnr,col) = cnst ;
          end
      end
    end
    end       
end
catch
    
end

guidata(hObject, handles);
update_en(hObject, eventdata, handles);


function ot_Callback(hObject, eventdata, handles)
function ot_CreateFcn(hObject, eventdata, handles)

if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --------------------------------------------------------------------
function get_blank_Callback(hObject, eventdata, handles)

ppath =  [handles.dir 'templates\blank.par'];
handles.parameters = csvread(ppath);
ppath =  [handles.dir 'templates\blank.en'];
handles.electronegativities =  csvread(ppath);
ppath =  [handles.dir 'templates\blank.cpl'];
handles.all_couplings = csvread(ppath);
handles.sums = zeros(0,25);
    
guidata(hObject, handles);
recalc_cpl(hObject, eventdata, handles);
update_constants(hObject, eventdata, handles);
update_en(hObject, eventdata, handles);
set_parameters(hObject, eventdata, handles);

% --------------------------------------------------------------------
function get_ribo_Callback(hObject, eventdata, handles)

ppath =  [handles.dir 'templates\ribo.par'];
handles.parameters = csvread(ppath);
ppath =  [handles.dir 'templates\ribo.en'];
handles.electronegativities =  csvread(ppath);
ppath =  [handles.dir 'templates\ribo.cpl'];
handles.all_couplings = csvread(ppath);
handles.sums = zeros(0,25);
    
guidata(hObject, handles);
recalc_cpl(hObject, eventdata, handles);
update_constants(hObject, eventdata, handles);
update_en(hObject, eventdata, handles);
set_parameters(hObject, eventdata, handles);

% --------------------------------------------------------------------
function get_deox_Callback(hObject, eventdata, handles)

ppath =  [handles.dir 'templates\deox.par'];
handles.parameters = csvread(ppath);
ppath =  [handles.dir 'templates\deox.en'];
handles.electronegativities =  csvread(ppath);
ppath =  [handles.dir 'templates\deox.cpl'];
handles.all_couplings = csvread(ppath);
handles.sums = zeros(0,25);
 
guidata(hObject, handles);
recalc_cpl(hObject, eventdata, handles);
update_constants(hObject, eventdata, handles);
update_en(hObject, eventdata, handles);
set_parameters(hObject, eventdata, handles);

function GetPars_Callback(hObject, eventdata, handles)

% --------------------------------------------------------------------
function get_ex1_Callback(hObject, eventdata, handles)

ppath =  [handles.dir 'examples\1\ex.par'] ;
handles.parameters = csvread(ppath);
ppath =  [handles.dir 'examples\1\ex.en'] ;
handles.electronegativities =  csvread(ppath);
ppath =  [handles.dir 'examples\1\ex.cpl'] ;
handles.all_couplings = csvread(ppath);
handles.sums = zeros(0,25);

    
guidata(hObject, handles);
recalc_cpl(hObject, eventdata, handles);
update_constants(hObject, eventdata, handles);
update_en(hObject, eventdata, handles);
set_parameters(hObject, eventdata, handles);


% --------------------------------------------------------------------
function get_ex2_Callback(hObject, eventdata, handles)

ppath =  [handles.dir 'examples\2\ex.par'] ;
handles.parameters = csvread(ppath);
ppath =  [handles.dir 'examples\2\ex.en'] ;
handles.electronegativities =  csvread(ppath);
ppath =  [handles.dir 'examples\2\ex.cpl'] ;
handles.all_couplings = csvread(ppath);
handles.sums = zeros(0,25);

    
guidata(hObject, handles);
recalc_cpl(hObject, eventdata, handles);
update_constants(hObject, eventdata, handles);
update_en(hObject, eventdata, handles);
set_parameters(hObject, eventdata, handles);

% --------------------------------------------------------------------
function Maxiter_Callback(hObject, eventdata, handles)

% --------------------------------------------------------------------
function edit_en_Callback(hObject, eventdata, handles)

assignin('base','elec',handles.electronegativities);
electroneg;


% --------------------------------------------------------------------
function prefs_Callback(hObject, eventdata, handles)
% --------------------------------------------------------------------

function scan_Callback(hObject, eventdata, handles)
elapsed_time = clock;

fprintf ( '\n\n////////////////////////////////////////////////////////////// \n');
fprintf ( '////////////////////////////////////////////////////////////// \n\n');
fprintf ( '  START Pseudorotational scan \n');
fprintf ( '\n////////////////////////////////////////////////////////////// \n');
fprintf ( '////////////////////////////////////////////////////////////// \n\n');


loop_p = 20;
loop_phi = 40;

handles.parameters(1) = loop_p;
handles.parameters(2) = loop_phi;
handles.parameters(3) = 0;
handles.parameters(4) = 0;
handles.parameters(5) = 0;
handles.parameters(6) = 0;
handles.parameters(7) = 50;
handles.parameters(8) = 0;
handles.parameters(9) = 0;
handles.parameters(10) = 0;
handles.parameters(11) = get(handles.Pb_check , 'Value');
handles.parameters(12) = get(handles.phib_check , 'Value');
handles.parameters(13) = 0;
handles.parameters(14) = get(handles.partition_check , 'Value') ;
handles.parameters(15) = get(handles.temp , 'Value') ;

try
    if handles.parameters(14) == 0
        handles.parameters(7) = str2num( get(handles.partition_value , 'String')) ;
    end
catch
end
try
    handles.parameters(4) = str2num( get(handles.Pb_value , 'String')) ;
catch
end
try
    handles.parameters(5) = str2num( get(handles.phib_value , 'String')) ;
catch
end



xxx= [];
yyy= [];
matrix = zeros(7);
parameters = horzcat( handles.parameters(1:6) * pi() / 180 , handles.parameters(7));
counter  = 1;
counter2 = 1;
numtemp = handles.parameters(15);

for ii = 1:7
    if handles.parameters(7+ii) == 1
        if ii ~= 7
            xxx(counter)  = parameters(ii);
            counter = counter + 1;
        else
            for iii = 1:numtemp
                xxx(counter)  = parameters(ii)/100;
                counter = counter + 1;
            end
        end
    else
        if ii ~= 7
            yyy(counter2)  = parameters(ii);
            counter2 = counter2 + 1;
        else
            for iii = 1:numtemp
                yyy(counter2)  = parameters(ii)/100;
                counter2 = counter2 + 1;
            end
        end
    end
end

counter2 = 1;
for ii = 1:7
    if handles.parameters(7+ii) == 1
        if ii ~= 7
            matrix(counter2,ii) = 1;
            counter2 = counter2 + 1;
        else
            for iii = 1:numtemp
                matrix(counter2,ii + iii - 1) = 1;
                counter2 = counter2 + 1;
            end
        end
    else
        if ii ~= 7
            matrix(counter,ii) = 1;
            counter = counter + 1;
        else
            for iii = 1:numtemp
                matrix(counter,ii + iii - 1) = 1;
                counter = counter + 1;
            end
        end
    end
end


counter = 1;
cpl = [];
for ii = 1:20
    if handles.all_couplings(ii,4) == 1
        cpl = vertcat ( cpl , handles.all_couplings(ii,:) );
        cpl ( counter , 10 ) = handles.electronegativities( cpl ( counter , 10 ) );
        cpl ( counter , 11 ) = handles.electronegativities( cpl ( counter , 11 ) );
        cpl ( counter , 12 ) = handles.electronegativities( cpl ( counter , 12 ) );
        cpl ( counter , 13 ) = handles.electronegativities( cpl ( counter , 13 ) );
        cpl ( counter ,  4 ) = (ii - 1 - mod(ii - 1,4) + 4 ) /4 ;
        counter = counter + 1 ;
    end
end

if size(cpl) > 0

    cpl = horzcat ( cpl(:,5:5+numtemp-1) , cpl(:,10:13) , cpl(:,3)/180*pi() , cpl(:,4) , cpl(:,14:15), cpl(:,16)/180*pi() );

    ncpl = size(cpl);
    ncpl = ncpl(1);
    sums = handles.sums;
    numsum = size(sums);
    numsum = numsum(1);

    if numsum > 0
        for i = 20:-1:1
            if handles.all_couplings(i,4) == 0
                sums(:,i) = [];
            end
        end
    end


    options = optimset('MaxIter', 250, 'Display', 'off', 'GradObj', 'off');
    elapsed_time = clock;

     
        final  = [];
        final2 = [];
        fprintf('\n\nProgress: [');
        
        grd = 60; 
        grd1 = grd + 1; 
                    
        if size(xxx) > 0

                    qq = size(xxx);
        qq = qq(2);
        if handles.parameters(14) == 1

            A = zeros(numtemp*2 , qq);
            B = zeros(numtemp*2 , 1);
            counter=1;
            for i = 1:numtemp
                A( counter     , qq - numtemp + i ) = 1;
                A( counter + 1 , qq - numtemp + i ) = -1;
                B( counter     , 1 ) = 1;
                B( counter + 1 , 1 ) = 0;
                counter = counter + 2 ;
            end

        else
            A = zeros(1, qq);
            B = [0];
        end

        QQ = [];
        QQQ = [];
        jj = 0;
        for ii = 1:6
            if handles.parameters(7+ii) == 1
                jj = jj + 1;
                CC = zeros(2,qq);
                CC(1:2,jj) = [1 , -1];
                DD = [4*pi() , 4*pi()]' ;
                if ii == 2
                    DD = [pi()/3 , pi()/3]' ;
                end
                if ii == 5
                    DD = [pi()/3 , pi()/3]' ;
                end
                QQ  = vertcat( QQ  , CC );
                QQQ = vertcat( QQQ , DD );
            end
        end

        A  = vertcat( QQ   , A );
        B  = vertcat( QQQ  , B );


        
              
        index1 =  handles.parameters(11) ;
        index2 =  handles.parameters(12) * (1 + index1) ;
        
              if index2 ~= 0
                xxx(index2) = 40 / 180 * pi();
              end
        ii = 0;    
        for i = 1:grd1                                                                                 %loop for phi_max

            klo(1) =-pi()/3 + (i-1) / grd * pi() * 2/3 ;
            ii = ii + 1;
            if ii >= grd1/21 
                ii = 0;
                fprintf('%2.0f ',floor(i/grd1*100))               
            end
            for j = 1:grd1                                                                              % loop for P
                
                yyy(2) = -pi()/3 + (j-1) / grd * pi() * 2/3 ;
                [ yyy(1) , yyy(2) ] = cart2pol( klo(1) , yyy(2) );         
                
                if yyy(2) > 1.02*pi()/3
                    final = vertcat(final, horzcat( yyy(1:2) , NaN, xxx ));
                else
                    
                if index1 ~= 0
                    xxx(index1) = mod(yyy(1) + pi(),2*pi());
                end

                fff = @(var) search_minimum(var,yyy,matrix,cpl , numtemp , 0 , handles.karplus, sums);
                [optimised_parameters, total_error] = fmincon(fff,xxx, A,B, [],[],[],[],[],options);

                if ( yyy(2) < pi()/12)
                if index1 ~= 0
                    xxx(index1) =yyy(1) ;
                    fff = @(var) search_minimum(var,yyy,matrix,cpl , numtemp , 0 , handles.karplus, sums);
                    [optimised_parameters2, total_error2] = fmincon(fff,xxx, A,B, [],[],[],[],[],options);
                    total_error = min(total_error, total_error2);
                    if total_error == total_error2
                      optimised_parameters = optimised_parameters2; 
                    end
                end                 
                end
                
                final = vertcat(final, horzcat( yyy(1:2) , total_error, optimised_parameters));
                end
            end
        end
        
        else
        
        ii = 0;
        for i = 1:grd1                                                                                
            klo(1) = -pi()/3 + (i-1) / grd * pi() * 2/3 ;  
            ii = ii + 1;
            if ii >= grd1/21 
                ii = 0;
                fprintf('%2.0f ',floor(i/grd1*100))               
            end
            for j = 1:grd1                                                                             
                yyy(2) = -pi()/3 + (j-1) / grd * pi() * 2/3 ; 
                [ yyy(1) , yyy(2) ] = cart2pol( klo(1) , yyy(2) );
                if yyy(2) > 1.02*pi()/3
                    final = vertcat(final, horzcat( yyy(1:2) , NaN, xxx ));
                else
                total_error = search_minimum(xxx,yyy,matrix,cpl , numtemp , 0 , handles.karplus, sums);
                final = vertcat(final, horzcat( yyy(1:2) , total_error, xxx ));
                end
            end
        end            
                      
        end
        
        fprintf( ']\n\n');

        [te,mini] = min(final(:,3));

        A = linspace(-1,1,grd1);    
        C = reshape(final(:,3),grd1,grd1);
        
        assignin('base','final',final);
        
h = figure(999);   
[hh,h] =contour(A,A,C,50); %,'EdgeColor','none')

LL = linspace(min(min(C)),1.5*min(min(C))+0.1,20);
set(h,'LevelList',LL);
hold;


for i = 1 : 361
 H(i) = cos(i/180*3.141592);
 I(i) = sin(i/180*3.141592);
 J(i) = 5;
end

ola = plot3(H*1.04,I*1.04,J,'o','Marker','o','Markersize' , 10.0  );
set(ola, 'MarkerEdgeColor', [0.99 0.99 0.99]);
set(ola, 'MarkerFaceColor', [0.99 0.99 0.99]);
set(ola, 'Color', [0.99 0.99 0.99]);
set(ola, 'LineWidth', 40);

ola = mmpolar([0 0],[1 1],'x', 'FontName', 'Palatino Linotype','FontSize', 8, 'RTickOffset' , 0.08,'TTickOffset' , 0.10, 'RLimit', [0 60], 'TTickValue', [0 30 60 90 120 150 180 210 240 270 300 330], 'RTickValue', [0 15 30 45 60], 'RTickLabel', {'' '15\circ' '30\circ' '45\circ' '60\circ'}, 'RTickAngle', [75] );
view([90 -90]);
set (ola, 'Visible' , 'Off' );

colorbar;

print( '-f999', handles.output_dpi ,  handles.output_format , [ handles.dir get( handles.file_path , 'String' ) handles.output_ext] );

end


% --------------------------------------------------------------------
% -----------   MENU  ITEMS  -----------------------------------------
% --------------------------------------------------------------------

% --------------------------------------------------------------------
function Format_Callback(hObject, eventdata, handles)
% --------------------------------------------------------------------
function Resolution_Callback(hObject, eventdata, handles)
% --------------------------------------------------------------------
function out_jpg_Callback(hObject, eventdata, handles)
handles.output_format = '-djpeg100';
handles.output_ext = '.jpg';
set(handles.out_jpg,'checked','off');
set(handles.out_eps,'checked','off');
set(handles.out_pdf,'checked','off');
set(handles.out_tiff,'checked','off');
set(handles.out_fig,'checked','off');
set(hObject,'checked','on');
guidata(hObject, handles);
% --------------------------------------------------------------------
function out_eps_Callback(hObject, eventdata, handles)
handles.output_format = '-depsc';
handles.output_ext = '.eps';
set(handles.out_jpg,'checked','off');
set(handles.out_eps,'checked','off');
set(handles.out_pdf,'checked','off');
set(handles.out_tiff,'checked','off');
set(handles.out_fig,'checked','off');
set(hObject,'checked','on');
guidata(hObject, handles);
% --------------------------------------------------------------------
function out_pdf_Callback(hObject, eventdata, handles)
handles.output_format = '-dpdf';
handles.output_ext = '.pdf';
set(handles.out_jpg,'checked','off');
set(handles.out_eps,'checked','off');
set(handles.out_pdf,'checked','off');
set(handles.out_tiff,'checked','off');
set(handles.out_fig,'checked','off');
set(hObject,'checked','on');
guidata(hObject, handles);
% --------------------------------------------------------------------
function out_tiff_Callback(hObject, eventdata, handles)
handles.output_format = '-dtiff';
handles.output_ext = '.tiff';
set(handles.out_jpg,'checked','off');
set(handles.out_eps,'checked','off');
set(handles.out_pdf,'checked','off');
set(handles.out_tiff,'checked','off');
set(handles.out_fig,'checked','off');
set(hObject,'checked','on');
guidata(hObject, handles);
% --------------------------------------------------------------------
function out_fig_Callback(hObject, eventdata, handles)
handles.output_format = 'fig';
handles.output_ext = '.fig';
set(handles.out_jpg,'checked','off');
set(handles.out_eps,'checked','off');
set(handles.out_pdf,'checked','off');
set(handles.out_tiff,'checked','off');
set(handles.out_fig,'checked','off');
set(hObject,'checked','on');
guidata(hObject, handles);
% --------------------------------------------------------------------
function out_75_Callback(hObject, eventdata, handles)
handles.output_dpi = '-r75';
set(handles.out_75,'checked','off');
set(handles.out_150,'checked','off');
set(handles.out_300,'checked','off');
set(handles.out_600,'checked','off');
set(hObject,'checked','on');
guidata(hObject, handles);
% --------------------------------------------------------------------
function out_150_Callback(hObject, eventdata, handles)
handles.output_dpi = '-r150';
set(handles.out_75,'checked','off');
set(handles.out_150,'checked','off');
set(handles.out_300,'checked','off');
set(handles.out_600,'checked','off');
set(hObject,'checked','on');
guidata(hObject, handles);
% --------------------------------------------------------------------
function out_300_Callback(hObject, eventdata, handles)
handles.output_dpi = '-r300';
set(handles.out_75,'checked','off');
set(handles.out_150,'checked','off');
set(handles.out_300,'checked','off');
set(handles.out_600,'checked','off');
set(hObject,'checked','on');
guidata(hObject, handles);
% --------------------------------------------------------------------
function out_600_Callback(hObject, eventdata, handles)
handles.output_dpi = '-r600';
set(handles.out_75,'checked','off');
set(handles.out_150,'checked','off');
set(handles.out_300,'checked','off');
set(handles.out_600,'checked','off');
set(hObject,'checked','on');
guidata(hObject, handles);
% --------------------------------------------------------------------
function Karplus_Callback(hObject, eventdata, handles)
% --------------------------------------------------------------------
function par09_Callback(hObject, eventdata, handles)
handles.karplus = 1;
set(handles.par09,'checked','off');
set(handles.par09b,'checked','off');
set(handles.par12,'checked','off');
set(handles.par20,'checked','off');
set(hObject,'checked','on');
guidata(hObject, handles);
% --------------------------------------------------------------------
function par12_Callback(hObject, eventdata, handles)
handles.karplus = 2;
set(handles.par09,'checked','off');
set(handles.par09b,'checked','off');
set(handles.par12,'checked','off');
set(handles.par20,'checked','off');
set(hObject,'checked','on');
guidata(hObject, handles);
% --------------------------------------------------------------------
function par20_Callback(hObject, eventdata, handles)
handles.karplus = 3;
set(handles.par09,'checked','off');
set(handles.par09b,'checked','off');
set(handles.par12,'checked','off');
set(handles.par20,'checked','off');
set(hObject,'checked','on');
guidata(hObject, handles);
% --------------------------------------------------------------------
function par09b_Callback(hObject, eventdata, handles)
handles.karplus = 4;
set(handles.par09,'checked','off');
set(handles.par09b,'checked','off');
set(handles.par12,'checked','off');
set(handles.par20,'checked','off');
set(hObject,'checked','on');
guidata(hObject, handles);


% --- Executes on selection change in tempe.
function temp_Callback(hObject, eventdata, handles)
% hObject    handle to tempe (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = get(hObject,'String') returns tempe contents as cell array
%        contents{get(hObject,'Value')} returns selected item from tempe


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

% Hint: popupmenu 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


% --- Executes on selection change in tempe.
function tempe_Callback(hObject, eventdata, handles)
% hObject    handle to tempe (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = get(hObject,'String') returns tempe contents as cell array
%        contents{get(hObject,'Value')} returns selected item from tempe


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

% Hint: popupmenu 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

function out=mmpolar(varargin)
%MMPOLAR Polar Plot with Settable Properties.
% MMPOLAR(Theta,Rho) creates a polar coordinate plot using the angle Theta
% in RADIANS and radius in Rho. Rho can contain negative values.
% MMPOLAR(Theta,Rho,S) creates the plot using the line spec given by S. See
% the function PLOT for information about S.
% MMPOLAR(Theta1,Rho1,S1,Theta2,Rho2,S2,...) plots all the defined curves.
%
% MMPOLAR(Theta1,Rho1,S1,...,'PName',PValue,...) plots all defined curves,
% and sets plot property names to the corresponding property values.
% MMPOLAR(Theta1,Rho1,S1,...,P) plots all the defined curves, and uses the
% structure P having fieldnames equal to plot property names to set
% corresponding property values contained in the associated fields.
%
% H=MMPOLAR(Theta,Rho,...) returns handles to lines or lineseries objects.
% For example, set(H,'LineWidth',2) sets all linewidths to 2 points.
% Note: 'LineWidth' is NOT a property that can be set with MMPOLAR. It must
% be set as shown above by using the SET function on the line handles H.
%
% MMPOLAR('PName',PValue,...) sets the property names to the corresponding
% property values. See below for property name/value pairs. Just as with
% the function SET 'PName' is case insensitive and need only be unique.
% MMPOLAR with no input argument returns a structure with fieldnames equal
% to property names each containing the associated property values.
% MMPOLAR(P) sets property values using the structure P as described above.
% MMPOLAR('PName') returns the property value associated with 'PName'.
% MMPOLAR({'PName1','PName2',...}) returns multiple property values in a
% cell array.
% MMPOLAR(Hax,...) uses the axes having handle Hax.
%
% Examples: MMPOLAR(Theta,Rho,S,'Style','compass') creates a polar plot with
% theta=0 pointing North and theta increasing in the clockwise direction.
%
% MMPOLAR(Theta,Rho,S) creates a cartesian polar plot where theta=0 is along
% the x-axis and theta increases in the counterclockwise direction.
%
% MMPOLAR works with HOLD, XLABEL, YLABEL, TITLE, ZOOM, SUBPLOT
% but does not work with AXIS, GRID (Use MMPOLAR properties to set these)
%
% See also POLAR, PLOT, HOLD
%
% PROPERTY          VALUE {Default}  DESCRIPTION
% Style             {cartesian} | compass  shortcut to two common polar
%                     styles. Cartesian: theta=0 points east and increases
%                     going north. Compass: theta=0 points north and
%                     increases going east. See TDirection and TZeroDirection.
% Axis              {on} | off  shortcut for grids, ticks, border,
%                     backgroundcolor, visibility
% Border            {on} | off  shortcut for axis border, tick mark visibility.
% Grid              {on} | off  shortcut for visibility of rho and theta grids
% RLimit            [Rmin Rmax] rho axis limits, may be negative values
% TLimit            [Tmin Tmax] theta axis limits in RADIANS
% RTickUnits        {''} string added to last rho tick label to denote units
% TTickScale        {degrees} | radians  theta axis tick label scaling
% TDirection        cw | {ccw} direction of increasing theta
% TZeroDirection    North | {East} | South | West  theta=0 axis direction
%
% BackgroundColor   {w}  colorspec for axis background color
% BorderColor       {k} colorspec for axis border and tick mark colors
% FontName          string  font name for tick labels
% FontSize          scalar  font size for tick labels
% FontWeight        {normal} | bold  font weight for tick labels
% TickLength        {.02} normalized length of rho and theta axis tick marks
%
% RGridColor        {k} colorspec for rho axis grid color
% RGridLineStyle    - | -- | {:} | -.  rho axis grid line style
% RGridLineWidth    {0.5}  rho axis grid line width in points
% RGridVisible      {on} | off  rho axis grid visibility
% RTickAngle        [scalar]  angular position of rho axis tick labels in
%                             TTickScale units
% RTickOffset       {.04} Normalized radial offset for rho tick labels
% RTickLabel        string cell array containing rho axis tick labels
% RTickLabelVisible {on} | off  visibility of rho axis tick labels
% RTickLabelHalign  {center} | left | right  horizontal
%                             alignment of rho axis tick labels
% RTickLabelValign  {middle} | top | cap | baseline | bottom  vertical
%                             alignment of rho axis tick labels
% RTickValue        [vector]  vector containing rho axis tick positions
% RTickVisible      {on} | off  rho axis tick visibility
%
% TGridColor        colorspec for theta axis grid color
% TGridLineStyle    - | -- | {:} | -.  theta axis grid line style
% TGridLineWidth    {0.5}  theta axis grid line width in points
% TGridVisible      {on} | off  theta axis grid visibility
% TTickDelta        theta axis tick spacing in TTickScale units
%                   {15 degrees or pi/12 radians}
% TTickDirection    {in} | out  direction of theta tick marks
% TTickOffset       {.08} normalized radial offset of theta tick labels
% TTickLabel        string cell array containing theta axis tick labels
% TTickLabelVisible {on} | off  visiblity of theta axis tick labels
% TTickSign         {+-} | + sign of theta tick labels
% TTickValue        [vector]  vector of theta ticks in TTickScale units
% TTickVisible      {on} | off  theta axis tick visibility

% D.C. Hanselman, University of Maine, Orono, ME 04469
% MasteringMatlab@yahoo.com
% Mastering MATLAB 7
% 2005-04-25, 2006-01-18, 2006-04-06, 2006-05-17, 2006-05-18
% 2006-10-03, 2007-03-04, 2008-03-18

%--------------------------------------------------------------------------
% Parse Inputs                                                 Parse Inputs
%--------------------------------------------------------------------------
% Find MMPOLAR axes if it exists
nargi=nargin;
% find MMPOLAR axes if it is supplied or if it is the current axes
if nargi>0 && isscalar(varargin{1}) && ishandle(varargin{1})
   HAxes=varargin{1}; % see if first argument is an MMPOLAR axes
   if strcmp(get(HAxes,'Tag'),'MMPOLAR_Axes')
      HFig=ancestor(HAxes,'figure');
      HoldIsON=strcmp(get(HAxes,'nextplot'),'add')...
            && strcmp(get(HFig,'nextplot'),'add');
      P=getappdata(HAxes,'MMPOLAR_Properties');
      Pfn=fieldnames(P);
      varargin(1)=[]; % strip initial axes handle off varargin
      nargi=nargi-1;  % varargin now contains rest of input arguments
   else
      local_error('First Argument is Not a Valid MMPOLAR Axes Handle.')
   end
else % see if MMPOLAR axes is current axes
   HFig=get(0,'CurrentFigure');
   if isempty(HFig)
      HAxes=[];
      Pfn=fieldnames(local_getDefaults);
      HoldIsON=false;
   else
      HAxes=get(HFig,'CurrentAxes');
      if isempty(HAxes)
         Pfn=fieldnames(local_getDefaults);
         HoldIsON=false;
      else
         if strcmp(get(HAxes,'Tag'),'MMPOLAR_Axes')
            HoldIsON=strcmp(get(HAxes,'nextplot'),'add')...
                  && strcmp(get(HFig,'nextplot'),'add');
            P=getappdata(HAxes,'MMPOLAR_Properties');
            Pfn=fieldnames(P);
         else % no MMPOLAR axes exists
            HAxes=[];
            Pfn=fieldnames(local_getDefaults);
            HoldIsON=false;
            set(HAxes,'NextPlot','replace') % hold off
         end
      end
   end
end
%--------------------------------------------------------------------------
% Consider input arguments                         Consider input arguments
%--------------------------------------------------------------------------
if nargi==0   % MMPOLAR() MMPOLAR() MMPOLAR() MMPOLAR() MMPOLAR() MMPOLAR()
   if ~isempty(HAxes)
      out=P; % return property structure if it exists
      return
   else
      local_error('No MMPOLAR Axes exists or is not Current Axes.')
   end
end
if nargi==1   % Consider SET and GET Requests Consider SET and GET Requests
   if ~isempty(HAxes)
      arg=varargin{1};
      if ischar(arg)   % MMPOLAR('Pname') MMPOLAR('Pname') MMPOLAR('Pname')
         [fn,errmsg]=local_isfield(Pfn,arg);
         error(errmsg)
         out=P.(fn);
         return
      elseif iscellstr(arg)              % MMPOLAR({'PName1','PName2',...})
         nc=length(arg);
         out=cell(1,nc);
         for k=1:nc
            [fn,errmsg]=local_isfield(Pfn,arg{k});
            error(errmsg)
            out{k}=P.(fn);
         end
         return
      elseif isstruct(arg)    % MMPOLAR(S) MMPOLAR(S) MMPOLAR(S) MMPOLAR(S)
         Sfn=fieldnames(arg);
         for k=1:length(Sfn)
            [fn,errmsg]=local_isfield(Pfn,Sfn{k});
            error(errmsg)
            S.(fn)=arg.(Sfn{k});
         end
         local_updatePlot(HAxes,S);
         return
      else
         local_error('Unknown Input Argument.')
      end
   else
      local_error('No MMPOLAR exists or is not Current Axes.')
   end
end
%           MMPOLAR('PName1',PValue1,'PName2',PValue2,'PName3',PValue3,...)
if rem(nargi,2)==0 && ischar(varargin{1}) && ~isempty(HAxes)
   for k=1:2:nargi-1
      PName=varargin{k};
      if ischar(PName)
         [fn,errmsg]=local_isfield(Pfn,PName);
         error(errmsg)
         S.(fn)=varargin{k+1};
      else
         local_error('String Input Property Name Argument Expected.')
      end
   end
   local_updatePlot(HAxes,S)
   return
elseif ischar(varargin{1})      % Unknown Input Unknown Input Unknown Input
   local_error('Unknown Input Arguments or NO MMPOLAR Axes Exists.')
   
elseif isnumeric(varargin{1})%MMPOLAR(Theta,Rho,...) MMPOLAR(Theta,Rho,...)
   % find out if there are appended 'PName',PValue pairs or a structure P
   last=[];
   k=3; % 'Pname' or P can't appear before 3rd argument
   while k<=nargi 
      vark=varargin{k};
      k=k+1;
      if ischar(vark)
         fn=local_isfield(Pfn,vark);
         if ~isempty(fn)
            if isempty(last)
               last=k-1;
            end
            S.(fn)=varargin{k};
            k=k+1; % skip known PValue
         end
      elseif isstruct(vark) % found appended structure
         if isempty(last)
            last=k-1;
         end
         Sfn=fieldnames(vark);
         for ki=1:length(Sfn)
            [fn,errmsg]=local_isfield(Pfn,Sfn{ki});
            error(errmsg)
            S.(fn)=vark.(Sfn{ki});
         end
      end
   end
   if ~isempty(last)
      varargin(last:end)=[]; % strip properties and values from input
   end
else
   local_error('Unknown Input Arguments.')
end
%--------------------------------------------------------------------------
% Now have valid data for plotting         Now have valid data for plotting
%--------------------------------------------------------------------------
if HoldIsON % a current held plot exists
   
   D=getappdata(HAxes,'MMPOLAR_Data');                    % get stored data
   P=getappdata(HAxes,'MMPOLAR_Properties');
   
   tmpaxes=axes('Position',get(HAxes,'Position'));
   try % the plot function should work with new data
      Hlines=plot(tmpaxes,varargin{:});
   catch
      delete(tmpaxes)
      local_error('Input Arguments Not Understood.')
   end
   D.TData=[D.TData; get(Hlines,{'XData'})];   % add to held data
   D.RData=[D.RData; get(Hlines,{'YData'})];
   D.LineColor=[D.LineColor; get(Hlines,{'Color'})];
   D.LineStyle=[D.LineStyle; get(Hlines,{'LineStyle'})];
   D.Marker=[D.Marker; get(Hlines,{'Marker'})];
   D.NumLines=length(D.TData);
   
   delete(Hlines)                % got the data, lines are no longer needed
   delete(D.HLines)              % delete original lines as well
   
   set(tmpaxes,'NextPlot','add') % hold on
   for k=1:D.NumLines        % plot ALL data to find new RTicks and RLimits
      plot(tmpaxes,D.TData{k},D.RData{k})
   end
   P.RLimit=get(tmpaxes,'YLim');                          % Rho axis limits
   P.RTickValue=get(tmpaxes,'YTick');              % Default Rho axis ticks
   delete(tmpaxes)                        % Temporary axes no longer needed
   [P,D]=local_getRTickValue(HAxes,P,D);              % get rho tick values
   
   D.RDataN=cell(D.NumLines,1);
   for k=1:D.NumLines % normalize rho data for plotting
      D.TData(k)={mod(D.TData{k},2*pi)}; % map theta into [0 2*pi]
      D.RDataN(k)={(D.RData{k}-D.RMin)/D.RLimitDiff};
   end
   
   P.TLimit=[0 2*pi];                                    % plot full circle
   [P,D]=local_getTTickValue(P,D);                  % get theta tick values
   [P,D]=local_placeAxesPatch(HAxes,P,D,1);% draw axes patch, border, ticks
   [P,D]=local_placeRGrid(HAxes,P,D,1);                     % Draw Rho Grid
   [P,D]=local_placeTGrid(HAxes,P,D,1);                   % Draw Theta Grid
   [P,D]=local_placeTTickLabel(HAxes,P,D,1);        % Add Theta Tick Labels
   [P,D]=local_placeRTickLabel(HAxes,P,D,1);         % Add Rho Tick Lablels
      
else % Hold is OFF
   
   try % the plot function should work now
      HAxes=newplot; % create axes
      D.HLines=plot(HAxes,varargin{:});
   catch
      delete(gcf)
      local_error('Input Arguments Not Understood.')
   end
   HFig=ancestor(HAxes,'figure');
   D.NumLines=length(D.HLines);                  % get all data for storage
   D.TData=get(D.HLines,{'XData'});
   D.RData=get(D.HLines,{'YData'});
   D.LineColor=get(D.HLines,{'Color'});
   D.LineStyle=get(D.HLines,{'LineStyle'});
   D.Marker=get(D.HLines,{'Marker'});
   
   P=local_getDefaults;          % get default properties, update as needed
   
   P.RLimit=get(HAxes,'YLim');                            % Rho axis limits
   P.RTickValue=get(HAxes,'YTick');                % Default Rho axis ticks
   [P,D]=local_getRTickValue(HAxes,P,D);              % get rho tick values
   
   D.RDataN=cell(D.NumLines,1);
   for k=1:D.NumLines                              % Condition plotted data
      % wrap angles into first revolution
      D.TData{k}=mod(D.TData{k},2*pi);
      % normalize rho data for plotting
      D.RDataN(k)={(D.RData{k}-D.RMin)/D.RLimitDiff};
   end
   P.TLimit=[0 2*pi];                                    % plot full circle
   [P,D]=local_getTTickValue(P,D);                  % get theta tick values
   delete(D.HLines)         % clear cartesian lines, then create polar axes
   [P,D]=local_placeAxesPatch(HAxes,P,D);  % draw axes patch, border, ticks
   [P,D]=local_placeRGrid(HAxes,P,D);                       % Draw Rho Grid
   [P,D]=local_placeTGrid(HAxes,P,D);                     % Draw Theta Grid
   [P,D]=local_placeTTickLabel(HAxes,P,D);          % Add Theta Tick Labels
   [P,D]=local_placeRTickLabel(HAxes,P,D);           % Add Rho Tick Lablels

end

xylims=[-1 1]*1.08;
% Finalize Axes View                                     Finalize Axes View
set(HAxes,'DataAspectRatio',[1 1 1],....
          'XLimMode','manual','YLimMode','manual',...
          'XLim',xylims,'YLim',xylims,...
          'Visible','Off','Tag','MMPOLAR_Axes')
Hlabels=get(HAxes,{'Xlabel','YLabel', 'Title'});
set([Hlabels{:}],'Visible','on') % make labels visible

% Plot the Data                                               Plot the Data
D.HLines=zeros(D.NumLines,1);      % storage for lineseries handles
set([HFig,HAxes],'NextPlot','add') % hold on
for k=1:D.NumLines                 % plot the normalized data
   tdata=D.TData{k};
   rdata=D.RDataN{k};
   xdata=rdata.*cos(tdata);
   ydata=rdata.*sin(tdata);
   D.HLines(k)=plot(HAxes,xdata,ydata,...
                  'Color',D.LineColor{k},...
                  'LineStyle',D.LineStyle{k},...
                  'Marker',D.Marker{k});
end
if HoldIsON
   set([HFig,HAxes],'NextPlot','add') % hold on
else
   set([HFig,HAxes],'NextPlot','replace') % hold off
end

% Store Data                                                     Store Data
setappdata(HAxes,'MMPOLAR_Properties',P)
setappdata(HAxes,'MMPOLAR_Data',D)

if nargout % output handles if requested
   out=D.HLines;
end

% Update Plot with 'PName' PValue pairs if they exist
if exist('S','var')==1
   local_updatePlot(HAxes,S)
end
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
% Local Functions                                           Local Functions
%--------------------------------------------------------------------------
function local_updatePlot(HAxes,S)              % local_updatePlot(HAxes,S)
% update MMPOLAR plot properties
% S contains known properties

P=getappdata(HAxes,'MMPOLAR_Properties');
D=getappdata(HAxes,'MMPOLAR_Data');
Sfn=fieldnames(S);

for kk=1:length(Sfn)
   switch Sfn{kk}
   case 'Axis'                                                       % Axis
      [istrue,onoff]=local_isonoff(S.Axis);
      if istrue
         set(D.HAPatch,'Visible',onoff)
         set(D.HRGrid,'Visible',onoff)
         set(D.HTGrid,'Visible',onoff)
         set(D.HRTick,'Visible',onoff)
         set(D.HTTick,'Visible',onoff)
         set(D.HRTickLabel,'Visible',onoff)
         set(D.HTTickLabel,'Visible',onoff)
         P.RGridVisible=onoff;
         P.TGridVisible=onoff;
         P.RTickLabelVisible=onoff;
         P.TTickLabelVisible=onoff;         
      else
         local_error('Unknown ''Axis'' Property Value.')
      end
   case 'BackgroundColor'                                 % BackgroundColor
      [istrue,cs]=local_iscolorspec(S.BackgroundColor);
      if istrue
         set(D.HAPatch,'FaceColor',cs)
         P.BackgroundColor=cs;
      else
         local_error('Unknown ''BackgroundColor'' Property Value.')
      end
   case 'Border'                                                   % Border
      [istrue,onoff]=local_isonoff(S.Border);
      if istrue && strcmp(onoff,'on')
         set(D.HAPatch,'EdgeColor',P.BorderColor)
         set(D.HTTick,'Visible','on')
         set(D.HRTick,'Visible','on')
         P.RTickVisible='on';
         P.TTickVisibel='on';
      elseif istrue && strcmp(onoff,'off')
         set(D.HAPatch,'EdgeColor','none')
         set(D.HTTick,'Visible','off')
         set(D.HRTick,'Visible','off')
         P.RTickVisible='off';
         P.TTickVisibel='off';
      else
         local_error('Unknown ''Border'' Property Value.')
      end
   case 'BorderColor'                                         % BorderColor
      [istrue,cs]=local_iscolorspec(S.BorderColor);
      if istrue
         P.BorderColor=cs;
         set(D.HAPatch,'EdgeColor',cs)
         set(D.HRTick,'Color',cs)
         set(D.HTTick,'Color',cs)
      else
         local_error('Unknown ''BorderColor'' Property Value.')
      end
   case 'FontName'                                               % FontName
      if ischar(S.FontName) && any(strcmpi(listfonts,S.FontName))
         set([D.HRTickLabel; D.HTTickLabel],'FontName',S.FontName)
         P.FontName=S.FontName;
      else
         local_error('Unknown ''FontName'' Property Value.')
      end
   case 'FontSize'                                               % FontSize
      if isnumeric(S.FontSize) && isscalar(S.FontSize)
         set([D.HRTickLabel; D.HTTickLabel],'FontSize',S.FontSize)
         P.FontSize=S.FontSize;
      else
         local_error('Unknown ''FontSize'' Property Value.')
      end
   case 'FontWeight'                                           % FontWeight
      if ischar(S.FontWeight) && ...
            (strncmpi(S.FontWeight,'normal',3)...
            ||strncmpi(S.FontWeight,'bold',3))
         set([D.HRTickLabel; D.HTTickLabel],'FontWeight',S.FontWeight)
         P.FontWeight=S.FontWeight;
      else
         local_error('Unknown ''FontWeight'' Property Value.')
      end      
   case 'Grid'                                                       % Grid
      [istrue,onoff]=local_isonoff(S.Grid);
      if istrue
         set(D.HRGrid,'Visible',onoff)
         set(D.HTGrid,'Visible',onoff)
         P.Grid=onoff;
         P.RGridVisible=onoff;
         P.TGridVisible=onoff;
      else
         local_error('Unknown ''Grid'' Property Value.')
      end
   case 'RGridColor'                                           % RGridColor
      [istrue,cs]=local_iscolorspec(S.RGridColor);
      if istrue
         set(D.HRGrid,'Color',cs)
         set(D.HRTickLabel,'Color',cs)
         P.RGridColor=cs;
      else
         local_error('Unknown ''RGridColor'' Property Value.')
      end
   case 'RGridLineStyle'                                   % RGridLineStyle
      if local_islinespec(S.RGridLineStyle)
         set(D.HRGrid,'LineStyle',S.RGridLineStyle)
         P.RGridLineStyle=S.RGridLineStyle;
      else
         local_error('Unknown ''RGridLineStyle'' Property Value.')
      end               
   case 'RGridLineWidth'                                   % RGridLineWidth
      if isnumeric(S.RGridLineWidth) && isscalar(S.RGridLineWidth)
         set(D.HRGrid,'LineWidth',S.RGridLineWidth)
         P.RGridLineWidth=S.RGridLineWidth;
      else
         local_error('Unknown ''RGridLineWidth'' Property Value.')
      end
   case 'RGridVisible'                                       % RGridVisible
      [istrue,onoff]=local_isonoff(S.RGridVisible);
      if istrue
         set(D.HRGrid,'Visible',onoff)
         P.RGridVisible=onoff;
      else
         local_error('Unknown ''RGridVisible'' Property Value.')
      end      
   case 'RLimit'                                                   % RLimit
      if isnumeric(S.RLimit) && numel(S.RLimit)==2
         S.RLimit=[min(S.RLimit) max(S.RLimit)];
         S.RLimit(isinf(S.RLimit))=P.RLimit(isinf(S.RLimit));
         [P,D]=local_getRTickValue(HAxes,P,D,S);
         [P,D]=local_placeRGrid(HAxes,P,D,S);
         [P,D]=local_placeRTickLabel(HAxes,P,D,S);
         % rescale rho data to new limits
         for k=1:length(D.RData)
            D.RDataN(k)={(D.RData{k}-D.RMin)/D.RLimitDiff};
            D.RDataN{k}(D.RDataN{k}>1)=NaN; % hide data outside limits
            D.RDataN{k}(D.RDataN{k}<0)=NaN;
            theta=D.TData{k};
            xdata=D.RDataN{k}.*cos(theta);
            ydata=D.RDataN{k}.*sin(theta);
            set(D.HLines(k),'XData',xdata,'YData',ydata)
         end
      else
         local_error('Unknown ''RLimit'' Property Value.')
      end               
   case 'RTickAngle'                                           % RTickAngle
      if isnumeric(S.RTickAngle) && isscalar(S.RTickAngle)
         rad=S.RTickAngle;
         if strcmp(P.TTickScale,'degrees')
            rad=S.RTickAngle*pi/180;
         end
         if P.TLimit(1)>P.TLimit(2) && (rad<P.TLimit(2) || rad>P.TLimit(1))
            P.RTickAngle=S.RTickAngle;
            D.RTickAngle=rad;
         elseif rad>P.TLimit(1) && rad<P.TLimit(2)
            P.RTickAngle=S.RTickAngle;
            D.RTickAngle=rad;            
         else
            local_error('RTickAngle not within Theta Axis Limits.')
         end
         for k=1:D.RTickLabelN % ignore innermost tick
            xdata=(P.RTickOffset+D.RTickRadius(k))*cos(D.RTickAngle);
            ydata=(P.RTickOffset+D.RTickRadius(k))*sin(D.RTickAngle);
            set(D.HRTickLabel(k),'Position',[xdata ydata])
         end
         phi=asin(P.TickLength./(2*D.RTickRadius));
         tdata=[D.RTickAngle-phi; D.RTickAngle+zeros(size(D.RTickRadius))
                D.RTickAngle+phi; NaN(size(D.RTickRadius))];
         rdata=[D.RTickRadius; D.RTickRadius
                D.RTickRadius; NaN(size(D.RTickRadius))];
         xdata=rdata(:).*cos(tdata(:));
         ydata=rdata(:).*sin(tdata(:));
         set(D.HRTick,'XData',xdata,'YData',ydata) % move rho ticks
      else
         local_error('Unknown ''RTickAngle'' Property Value.')
      end
   case 'RTickLabel'                                           % RTickLabel
      if iscellstr(S.RTickLabel)
         NumS=length(S.RTickLabel);
         for k=1:D.RTickLabelN
            str=S.RTickLabel{rem(k-1,NumS)+1};
            set(D.HRTickLabel(k),'String',str)
         end         
         P.RTickLabel=S.RTickLabel;
      else
         local_error('Unknown ''RTickLabel'' Property Value.')
      end
   case 'RTickLabelHalign'                               % RTickLabelHalign
      fnames={'left' 'center' 'right'};
      out=local_isfield(fnames,S.RTickLabelHalign);
      if ~isempty(out)
         P.RTickLabelHalign=out;
         set(D.HRTickLabel,'HorizontalAlignment',out)
      else
         local_error('Unknown ''RTickLabelHalign'' Property Value.')
      end         
   case 'RTickLabelValign'                               % RTickLabelValign
      fnames={'top' 'cap' 'middle' 'baseline' 'bottom'};
      out=local_isfield(fnames,S.RTickLabelValign);
      if ~isempty(out)
         P.RTickLabelValign=out;
         set(D.HRTickLabel,'VerticalAlignment',out)
      else
         local_error('Unknown ''RTickLabelValign'' Property Value.')
      end
   case 'RTickLabelVisible'                             % RTickLabelVisible
      [istrue,onoff]=local_isonoff(S.RTickLabelVisible);
      if istrue
         set(D.HRTickLabel,'Visible',onoff)
         P.RTickLabelVisible=onoff;
      else
         local_error('Unknown ''RTickLabelVisible'' Property Value.')
      end
   case 'RTickOffset'                                         % RTickOffset
      if isnumeric(S.RTickOffset) && isscalar(S.RTickOffset)
         P.RTickOffset=S.RTickOffset;
         for k=1:D.RTickLabelN
            xdata=(P.RTickOffset+D.RTickRadius(k))*cos(D.RTickAngle);
            ydata=(P.RTickOffset+D.RTickRadius(k))*sin(D.RTickAngle);
            set(D.HRTickLabel(k),'Position',[xdata ydata])
         end
      else
          local_error('Unknown ''RTickOffset'' Property Value.')
      end
   case 'RTickUnits'                                           % RTickUnits
      if ischar(S.RTickUnits)
         tmp=char(get(D.HRTickLabel(end),'String'));
         if ~isempty(P.RTickUnits)
            idx=strfind(tmp,P.RTickUnits);
            tmp=[tmp(1:idx(end)-1) S.RTickUnits];
         else
            tmp=[tmp S.RTickUnits]; %#ok
         end
         set(D.HRTickLabel(end),'String',tmp)
         P.RTickUnits=S.RTickUnits;
      else
         local_error('Unknown ''RTickUnits'' Property Value.')
      end
   case 'RTickValue'                                           % RTickValue
      if isnumeric(S.RTickValue) && numel(S.RTickValue)>0
         S.RTickValue=S.RTickValue(S.RTickValue>=P.RLimit(1)...
            & S.RTickValue<=P.RLimit(2));
         if length(S.RTickValue)>1
            P.RTickValue=S.RTickValue;
            D.RTickLabelN=length(P.RTickValue);
            D.RTickRadius=(P.RTickValue-D.RMin)/D.RLimitDiff;
            [P,D]=local_placeRGrid(HAxes,P,D,S);
            [P,D]=local_placeRTickLabel(HAxes,P,D,S);
         end
      else
         local_error('Unknown ''RTickValue'' Property Value.')
      end
   case 'RTickVisible'                                       % RTickVisible
      [istrue,onoff]=local_isonoff(S.RTickVisible);
      if istrue
         set(D.HRTick,'Visible',onoff)
         P.RTickVisible=onoff;
      else
         local_error('Unknown ''RTickVisible'' Property Value.')
      end
   case 'Style'                                                     % Style
      if strncmpi(S.Style,'cartesian',3)   % Cartesian style
         set(HAxes,'View',[0 90])
         P.TDirection='ccw';
         P.TZeroDirection='east';
         P.Style='cartesian';
      elseif strncmpi(S.Style,'compass',3) % Compass style
         set(HAxes,'View',[90 -90])
         P.TDirection='cw';
         P.TZeroDirection='north';
         P.Style='compass';
      else
         local_error('Unknown ''Style'' Property Value.')
      end
   case 'TDirection'                                           % TDirection
      if ischar(S.TDirection) && strcmpi(S.TDirection,'cw')
         P.TDirection='cw';
         if strcmp(P.TZeroDirection,'north')
            set(HAxes,'View',[90 -90])
         elseif strcmp(P.TZeroDirection,'east')
            set(HAxes,'View',[0 -90])
         elseif strcmp(P.TZeroDirection,'south')
            set(HAxes,'View',[270 -90])
         elseif strcmp(P.TZeroDirection,'west')
            set(HAxes,'View',[180 -90])
         end
      elseif ischar(S.TDirection) && strcmpi(S.TDirection,'ccw')
         P.TDirection='ccw';
         if strcmp(P.TZeroDirection,'north')
            set(HAxes,'View',[270 90])
         elseif strcmp(P.TZeroDirection,'east')
            set(HAxes,'View',[0 90])
         elseif strcmp(P.TZeroDirection,'south')
            set(HAxes,'View',[90 90])
         elseif strcmp(P.TZeroDirection,'west')
            set(HAxes,'View',[180 90])
         end
      else
         local_error('Unknown ''TDirection'' Property Value.')
      end
      P.Style='unknown';
   case 'TGridColor'                                           % TGridColor
      [istrue,cs]=local_iscolorspec(S.TGridColor);
      if istrue
         set(D.HTGrid,'Color',cs)
         set(D.HTTickLabel,'Color',cs)
         P.TGridColor=cs;
      else
         local_error('Unknown ''TGridColor'' Property Value.')
      end
   case 'TGridLineStyle'                                   % TGridLineStyle
      if local_islinespec(S.TGridLineStyle)
         set(D.HTGrid,'LineStyle',S.TGridLineStyle)
         P.TGridLineStyle=S.TGridLineStyle;
      else
         local_error('Unknown ''TGridLineStyle'' Property Value.')
      end
   case 'TGridLineWidth'                                   % TGridLineWidth
      if isnumeric(S.TGridLineWidth) && isscalar(S.TGridLineWidth)
         set(D.HTGrid,'LineWidth',S.TGridLineWidth)
         P.TGridLineWidth=S.TGridLineWidth;
      else
         local_error('Unknown ''TGridLineWidth'' Property Value.')
      end
   case 'TGridVisible'                                       % TGridVisible
      [istrue,onoff]=local_isonoff(S.TGridVisible);
      if istrue
         set(D.HTGrid,'Visible',onoff)
         P.TGridVisible=onoff;
      else
         local_error('Unknown ''TGridVisible'' Property Value.')
      end      
   case 'TickLength'                                           % TickLength
      if isnumeric(S.TickLength) && isscalar(S.TickLength)
         P.TickLength=max(min(abs(S.TickLength),0.1),.001);
         tdir=2*strcmp(P.TTickDirection,'in')-1;
         tdata=[D.TTickValue;D.TTickValue;NaN(1,D.TTickLabelN)];
         rdata=[ones(1,D.TTickLabelN)
                (1-tdir*P.TickLength)+zeros(1,D.TTickLabelN)
                NaN(1,D.TTickLabelN)];
         xdata=rdata(:).*cos(tdata(:));
         ydata=rdata(:).*sin(tdata(:));
         set(D.HTTick,'XData',xdata,'YData',ydata) % theta ticks
         phi=asin(P.TickLength./(2*D.RTickRadius));
         tdata=[D.RTickAngle-phi; D.RTickAngle+zeros(size(D.RTickRadius))
                D.RTickAngle+phi; NaN(size(D.RTickRadius))];
         rdata=[D.RTickRadius; D.RTickRadius
                D.RTickRadius; NaN(size(D.RTickRadius))];
         xdata=rdata(:).*cos(tdata(:));
         ydata=rdata(:).*sin(tdata(:));
         set(D.HRTick,'XData',xdata,'YData',ydata) % rho ticks
      else
         local_error('Unknown ''TickLength'' Property Value.')
      end
   case 'TLimit'                                                   % TLimit
      if isnumeric(S.TLimit) && numel(S.TLimit)==2
         if abs(diff(S.TLimit))>1.9*pi   % make full circle if close
            P.TLimit=[0 2*pi];
         else
            P.TLimit=mod(S.TLimit,2*pi); % move limits to range 0 to 2pi
         end
         [P,D]=local_getTTickValue(P,D,S);
         [P,D]=local_placeAxesPatch(HAxes,P,D,S);
         [P,D]=local_placeRGrid(HAxes,P,D,S);
         [P,D]=local_placeTGrid(HAxes,P,D,S);
         [P,D]=local_placeTTickLabel(HAxes,P,D,S);
         [P,D]=local_placeRTickLabel(HAxes,P,D,S);
         for k=1:length(D.TData) % hide data outside TLimits
            tdata=D.TData{k};
            if P.TLimit(1)>P.TLimit(2)
               tdata(tdata<P.TLimit(1) & tdata>P.TLimit(2))=NaN;
            else            
               tdata(tdata<P.TLimit(1) | tdata>P.TLimit(2))=NaN;
            end
            xdata=D.RDataN{k}.*cos(tdata);
            ydata=D.RDataN{k}.*sin(tdata);
            set(D.HLines(k),'XData',xdata,'YData',ydata)
         end
      else
         local_error('Unknown ''TLimit'' Property Value.')
      end
   case 'TTickDelta'                                           % TTickDelta
      if isnumeric(S.TTickDelta) && isscalar(S.TTickDelta)
         if strcmp(P.TTickScale,'degrees')
            P.TTickDelta=min(max(abs(S.TTickDelta),5),90);
         else
            P.TTickDelta=min(max(abs(S.TTickDelta),pi/36),pi/2);
         end
         [P,D]=local_getTTickValue(P,D,S);
         [P,D]=local_placeTGrid(HAxes,P,D,S);
         [P,D]=local_placeTTickLabel(HAxes,P,D,S);
      else
         local_error('Unknown ''TTickDelta'' Property Value.')
      end
   case 'TTickDirection'                                   % TTickDirection
      if ischar(S.TTickDirection) &&...
        (strcmpi(S.TTickDirection,'out') || strcmpi(S.TTickDirection,'in'))
         P.TTickDirection=S.TTickDirection;
         tdir=2*strcmp(P.TTickDirection,'in')-1;
         tdata=[D.TTickValue;D.TTickValue;NaN(1,D.TTickLabelN)];
         rdata=[ones(1,D.TTickLabelN)
                (1-tdir*P.TickLength)+zeros(1,D.TTickLabelN)
                NaN(1,D.TTickLabelN)];
         xdata=rdata(:).*cos(tdata(:));
         ydata=rdata(:).*sin(tdata(:));
         set(D.HTTick,'XData',xdata,'YData',ydata) % theta ticks
      else
         local_error('Unknown ''TTickDirection'' Property Value.')
      end
   case 'TTickLabel'                                           % TTickLabel
      if iscellstr(S.TTickLabel)
         NumS=length(S.TTickLabel);
         for k=1:D.TTickLabelN
            str=S.TTickLabel{rem(k-1,NumS)+1};
            set(D.HTTickLabel(k),'String',str)
         end         
         P.TTickLabel=S.TTickLabel;
      else
         local_error('Unknown ''TTickLabel'' Property Value.')
      end
   case 'TTickLabelVisible'                             % TTickLabelVisible
      [istrue,onoff]=local_isonoff(S.TTickLabelVisible);
      if istrue
         set(D.HTTickLabel,'Visible',onoff)
         P.TTickLabelVisible=onoff;
      else
         local_error('Unknown ''TTickLabelVisible'' Property Value.')
      end
   case 'TTickOffset'                                         % TTickOffset
      if isnumeric(S.TTickOffset) && isscalar(S.TTickOffset)
         P.TTickOffset=S.TTickOffset;
         for k=1:D.TTickLabelN
            xdata=(1+P.TTickOffset)*cos(D.TTickValue(k));
            ydata=(1+P.TTickOffset)*sin(D.TTickValue(k));
            set(D.HTTickLabel(k),'Position',[xdata ydata])
         end
      else
         local_error('Unknown ''TTickOffset'' Property Value.')
      end
   case 'TTickScale'                                           % TTickScale
      if ischar(S.TTickScale) && strncmpi(S.TTickScale,'degrees',3)...
                              && strcmp(P.TTickScale,'radians')
         P.TTickScale='degrees';
         P.TTickDelta=P.TTickDelta*180/pi;
         P.RTickAngle=P.RTickAngle*180/pi;
         [P,D]=local_getTTickValue(P,D,S);
         [P,D]=local_placeTTickLabel(HAxes,P,D,S);
      elseif ischar(S.TTickScale) && strncmpi(S.TTickScale,'radians',3)...
                                  && strcmp(P.TTickScale,'degrees')
         P.TTickScale='radians';
         P.TTickDelta=P.TTickDelta*pi/180;
         P.RTickAngle=P.RTickAngle*pi/180;
         [P,D]=local_getTTickValue(P,D,S);
         [P,D]=local_placeTTickLabel(HAxes,P,D,S);
      elseif ~ischar(S.TTickScale)
         local_error('Unknown ''TTickScale'' Property Value.')
      end
   case 'TTickSign'                                             % TTickSign
      if ischar(S.TTickSign)
         if strcmp(S.TTickSign,'+')
            P.TTickSign='+';
         else
            P.TTickSign='+-';
         end
         [P,D]=local_getTTickValue(P,D,S);
         [P,D]=local_placeTTickLabel(HAxes,P,D,S);
      else
         local_error('Unknown ''TTickSign'' Property Value.')
      end
   case 'TTickValue'                                           % TTickValue
      if isnumeric(S.TTickValue) && numel(S.TTickValue)>0
         TTick=S.TTickValue(:)';
         if strcmp(P.TTickScale,'degrees')
            TTick=TTick*pi/180;
         end
         if P.TLimit(1)>P.TLimit(2)
            idx=TTick<=P.TLimit(2) | TTick>=P.TLimit(1); % keepers
         else
            idx=TTick>=P.TLimit(1) & TTick<=P.TLimit(2); % keepers
         end
         S.TTickValue=S.TTickValue(idx);
         if length(S.TTickValue)>1
            P.TTickValue=S.TTickValue(:)';
            D.TTickValue=TTick(idx);
            D.TTickLabelN=length(P.TTickValue);
            [P,D]=local_placeTGrid(HAxes,P,D,S);
            [P,D]=local_placeTTickLabel(HAxes,P,D,S);
         end
      else
         local_error('Unknown ''TTickValue'' Property Value.')
      end
   case 'TTickVisible'                                       % TTickVisible
      [istrue,onoff]=local_isonoff(S.TTickVisible);
      if istrue
         set(D.HTTick,'Visible',onoff)
         P.TTickVisible=onoff;
      else
         local_error('Unknown ''RTickVisible'' Property Value.')
      end
   case 'TZeroDirection'                                   % TZeroDirection
      if ischar(S.TZeroDirection) && strncmpi(S.TZeroDirection,'north',1)
         P.TZeroDirection='north';
         if strcmp(P.TDirection,'ccw')
            set(HAxes,'View',[270 90])
         elseif strcmp(P.TDirection,'cw')
            set(HAxes,'View',[90 -90])
         end
      elseif ischar(S.TZeroDirection) && strncmpi(S.TZeroDirection,'east',1)
         P.TZeroDirection='east';
         if strcmp(P.TDirection,'ccw')
            set(HAxes,'View',[0 90])
         elseif strcmp(P.TDirection,'cw')
            set(HAxes,'View',[0 -90])
         end
      elseif ischar(S.TZeroDirection) && strncmpi(S.TZeroDirection,'south',1)
         P.TZeroDirection='south';
         if strcmp(P.TDirection,'ccw')
            set(HAxes,'View',[90 90])
         elseif strcmp(P.TDirection,'cw')
            set(HAxes,'View',[270 -90])
         end
      elseif ischar(S.TZeroDirection) && strncmpi(S.TZeroDirection,'west',1)
         P.TZeroDirection='west';
         if strcmp(P.TDirection,'ccw')
            set(HAxes,'View',[180 90])
         elseif strcmp(P.TDirection,'cw')
            set(HAxes,'View',[180 -90])
         end
      else
         local_error('Unknown ''TZeroDirection'' Property Value.')
      end
      P.Style='unknown';
   end
end
setappdata(HAxes,'MMPOLAR_Properties',P)
setappdata(HAxes,'MMPOLAR_Data',D)
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
function [out,errmsg]=local_isfield(fnames,str)             % local_isfield
% compare str to fnames, if found, return complete fieldname
% otherwise return error and empty string.
% fnames is cell array, str is a string
% outputs are strings

% look for exact match first
idx=find(strcmpi(fnames,str));

if isempty(idx) % no exact match, so look for more general match
   idx=find(strncmpi(str,fnames,max(length(str),2)));
end
if numel(idx)==1 % unique match found
   out=fnames{idx};
   errmsg='';
else             % trouble
   out='';
   errmsg=sprintf('Unknown or Not Unique Property: %s',str);
end
%--------------------------------------------------------------------------
function [istrue,cs]=local_iscolorspec(arg)             % local_iscolorspec
% see if arg is a valid color specification including 'none'
rgb={[1 0 0],[0 1 0],[0 0 1],[0 1 1],[1 0 1],[1 1 0],[0 0 0],[1 1 1],'none'};
cc='rgbcmykwn';
istrue=false;
cs=[];
if ~isempty(arg) && ischar(arg)
   idx=find(arg(1)==cc);
   if ~isempty(idx)
      istrue=true;
      cs=rgb{idx};
   end 
elseif ~isempty(arg) && isnumeric(arg)
   istrue=all(size(arg)==[1 3])...
      && all(arg>=0) && all(arg<=1);
   cs=arg;
else
   istrue=false;
end
%--------------------------------------------------------------------------
function istrue=local_islinespec(arg)                    % local_islinespec
% see if arg is a valid line style specification
str={'-','-.','--',':'};
if isempty(arg)
   istrue=false;
elseif ischar(arg) && length(arg)<3
   istrue=any(strncmp(str,arg,length(arg)));
else
   istrue=false;
end
%--------------------------------------------------------------------------
function [istrue,s]=local_isonoff(arg)                      % local_isonoff
% see if arg is either on or off
istrue=false;
s='';
if ~isempty(arg) && ischar(arg) && strcmpi(arg,'on')
   istrue=true;
   s='on';
elseif ~isempty(arg) && ischar(arg) && strcmpi(arg,'off')
   istrue=true;
   s='off';
end
%--------------------------------------------------------------------------
function local_error(arg)  
% Add message identifier to error message and post error
if ~isempty(arg) && ischar(arg)
   error('MMPOLAR:error',arg)
end
%--------------------------------------------------------------------------
function [P,D]=local_placeAxesPatch(HAxes,P,D,S)  %#ok local_placeAxesPatch
% Draw Axes Border and Patch
tinc=pi/250;
if P.TLimit(1)>P.TLimit(2)
   theta=[P.TLimit(1):tinc:(P.TLimit(2)+2*pi-eps) P.TLimit(2)+2*pi];
else
   theta=[P.TLimit(1):tinc:P.TLimit(2)-eps P.TLimit(2)];
end   
costheta=cos(theta);
sintheta=sin(theta);
if abs(diff(P.TLimit))<2*(1-eps)*pi; % less than 4 quadrant box
   xdata=[0 costheta 0];
   ydata=[0 sintheta 0];
else % four quadrant box
   xdata=costheta;
   ydata=sintheta;
end
% let the axes grow for less than 4 quadrant box
set(HAxes,'Xlim',[min(xdata) max(xdata)],'Ylim',[min(ydata) max(ydata)])
if nargin==3 % new plot
   D.HAPatch=patch('XData',xdata,'YData',ydata,...
      'Parent',HAxes,...
      'LineStyle','-',...
      'Linewidth',2*P.RGridLineWidth,...
      'EdgeColor',P.BorderColor,...
      'FaceColor',P.BackgroundColor,...
      'FaceAlpha',0,...
      'HandleVisibility','off',...
      'HitTest','off');
else % old plot, update data
   set(D.HAPatch,'XData',xdata,'YData',ydata)
end

%      'FaceColor',P.BackgroundColor,...

%--------------------------------------------------------------------------
function [P,D]=local_placeRGrid(HAxes,P,D,S)          %#ok local_placeRGrid
tinc=pi/250;
if P.TLimit(1)>P.TLimit(2)
   theta=[P.TLimit(1):tinc:(P.TLimit(2)+2*pi-eps) P.TLimit(2)+2*pi];
else
   theta=[P.TLimit(1):tinc:P.TLimit(2)-eps P.TLimit(2)];
end   
costheta=cos(theta);
sintheta=sin(theta);
xdata=[];
ydata=[];
% no outer grid if outer tick is at outer RLimit
D.RGridN=length(P.RTickValue(P.RTickValue<P.RLimit(2)));
for k=1:D.RGridN
   xdata=[xdata D.RTickRadius(k)*costheta NaN];  %#ok
   ydata=[ydata D.RTickRadius(k)*sintheta NaN];  %#ok
end
if nargin<4 % new grid
   D.HRGrid=line(xdata,ydata,...
      'Parent',HAxes,...
      'LineStyle',P.RGridLineStyle,...
      'LineWidth',P.RGridLineWidth,...
      'Color',P.RGridColor,...
      'HandleVisibility','off',...
      'HitTest','off');
else
   set(D.HRGrid,'Xdata',xdata,'YData',ydata)
end
% Draw Rho Axis Tick Marks
phi=asin(P.TickLength./(2*D.RTickRadius));
tdata=[D.RTickAngle-phi; D.RTickAngle+zeros(size(D.RTickRadius))
       D.RTickAngle+phi; NaN(size(D.RTickRadius))];
rdata=[D.RTickRadius; D.RTickRadius
       D.RTickRadius; NaN(size(D.RTickRadius))];
xdata=rdata(:).*cos(tdata(:));
ydata=rdata(:).*sin(tdata(:));
if nargin==3 % new plot
   D.HRTick=line(xdata,ydata,...
      'Parent',HAxes,...
      'LineStyle','-',...
      'LineWidth',2*P.RGridLineWidth,...
      'Color',P.BorderColor,...
      'HandleVisibility','off',...
      'HitTest','off');
else % old plot, update data
   set(D.HRTick,'XData',xdata,'YData',ydata)
end
%--------------------------------------------------------------------------
function [P,D]=local_placeRTickLabel(HAxes,P,D,S)%#ok local_placeRTickLabel
% Draw Rho Tick Labels
D.RScale=floor(log10(max(abs(P.RTickValue))));
if abs(D.RScale)<2
   D.RScale=0;
end
P.RTickLabel=cell(D.RTickLabelN,1);
if nargin==4 % delete old labels and create new ones
   delete(D.HRTickLabel)
end
D.HRTickLabel=zeros(D.RTickLabelN,1);
for k=1:D.RTickLabelN
   xdata=(P.RTickOffset+D.RTickRadius(k))*cos(D.RTickAngle);
   ydata=(P.RTickOffset+D.RTickRadius(k))*sin(D.RTickAngle);
   P.RTickLabel{k}=num2str(P.RTickValue(k)*10^(-D.RScale));
   if (k==D.RTickLabelN) && (D.RScale~=0)
      P.RTickLabel{k}=[P.RTickLabel{k} sprintf('\\times10^{%d}',D.RScale)];
   end
   D.HRTickLabel(k)=text(xdata,ydata,P.RTickLabel{k},...
      'Parent',HAxes,...
      'Color',P.TGridColor,...
      'FontName',P.FontName,...
      'FontSize',P.FontSize,...
      'FontWeight',P.FontWeight',...
      'HorizontalAlignment',P.RTickLabelHalign,...
      'VerticalAlignment',P.RTickLabelValign,...
      'Clipping','off',...
      'HandleVisibility','off',...
      'HitTest','off');
end
%--------------------------------------------------------------------------
function [P,D]=local_getRTickValue(HAxes,P,D,S)       % local_getRTickValue
% get RTicks
if nargin==4 % updating current ticks given new rho axis limits
   tmpaxes=axes('Position',get(HAxes,'Position'));
   line([0 1],S.RLimit,'Parent',tmpaxes);
   P.RTickValue=get(tmpaxes,'YTick');
   P.RLimit=S.RLimit;
   delete(tmpaxes) % got ticks, don't need axes anymore
end
%P.RTickValue(end)=P.RLimit(2); % place last tick at outer rho limit
NumRTick=length(P.RTickValue); % reduce ticks if too many
if NumRTick>6
   if rem(NumRTick,2)~=0  % odd number keep alternate ones
      P.RTickValue=P.RTickValue(1:2:end);
   else % even number, add one tick, then keep alternate ones
      if P.RTickValue(1)==0	% keep lowest tick if zero, add one at outside
         P.RTickValue=[P.RTickValue(1:2:end-1) ...
                        2*P.RTickValue(end)-P.RTickValue(end-1)];
         P.RTickValue(P.RTickValue>P.RLimit(2))=[]; % no label past limit
      else                    % add one tick at inside
         P.RTickValue=[2*P.RTickValue(1)-P.RTickValue(2)...
                        P.RTickValue(2:2:end)];
      end
      P.RLimit(1)=P.RTickValue(1); % make first tick lower axis limit
   end
end
if NumRTick<3
   m=sum(P.RTickValue)/2;
   P.RTickValue(3)=P.RTickValue(2);
   P.RTickValue(2)=m;
end
D.RLimitDiff=diff(P.RLimit);
D.RMin=P.RLimit(1);
if abs(P.RTickValue(1)-P.RLimit(1)) < abs(D.RLimitDiff)/100
   P.RTickValue(1)=[]; % throw out inner tick if at inner axis limit
end
D.RTickRadius=(P.RTickValue-D.RMin)/D.RLimitDiff;
D.RTickLabelN=length(P.RTickValue);
%--------------------------------------------------------------------------
function [P,D]=local_placeTGrid(HAxes,P,D,S)          %#ok local_placeTGrid
xdata=[];
ydata=[];
costheta=cos(D.TTickValue);
sintheta=sin(D.TTickValue);
% no grid on first or last ticks are axis limits if less than 4 quadrants
ki=1;
ke=length(D.TTickValue);
if abs(diff(P.TLimit))<2*(1-eps)*pi; % less than 4 quadrant box
   ki=1+(D.TTickValue(1)==P.TLimit(1));
   ke=length(D.TTickValue)-(D.TTickValue(end)==P.TLimit(2));
end
D.TGridN=ke-ki+1;
for k=ki:ke
   xdata=[xdata 0 costheta(k) NaN];  %#ok
   ydata=[ydata 0 sintheta(k) NaN];  %#ok
end
if nargin<4 % new grid
   D.HTGrid=line(xdata,ydata,...
      'Parent',HAxes,...
      'LineStyle',P.TGridLineStyle,...
      'LineWidth',P.TGridLineWidth,...
      'Color',P.TGridColor,...
      'HandleVisibility','off',...
      'HitTest','off');
else
   set(D.HTGrid,'Xdata',xdata,'YData',ydata)
end
% Draw Theta Axis Tick Marks
tdir=2*strcmp(P.TTickDirection,'in')-1;
tdata=[D.TTickValue; D.TTickValue; NaN(1,D.TTickLabelN)];
rdata=[ones(1,D.TTickLabelN)
       (1-tdir*P.TickLength)+zeros(1,D.TTickLabelN)
       NaN(1,D.TTickLabelN)];
xdata=rdata(:).*cos(tdata(:));
ydata=rdata(:).*sin(tdata(:));
if nargin==3 % new plot
   D.HTTick=line(xdata,ydata,...
      'Parent',HAxes,...
      'LineStyle','-',...
      'LineWidth',2*P.RGridLineWidth,...
      'Color',P.BorderColor,...
      'Clipping','off',...
      'HandleVisibility','off',...
      'HitTest','off');
else % old plot, update data
   set(D.HTTick,'XData',xdata,'YData',ydata)
end
%--------------------------------------------------------------------------
function [P,D]=local_placeTTickLabel(HAxes,P,D,S)%#ok local_placeTTickLabel
% Draw Theta Ticks
if nargin==4 % delete old labels and create new ones
   delete(D.HTTickLabel)
end
P.TTickLabel=cell(D.TTickLabelN,1);
D.HTTickLabel=zeros(D.TTickLabelN,1);
for k=1:D.TTickLabelN % label ticks at theta axis limits
   xdata=(1+P.TTickOffset)*cos(D.TTickValue(k));
   ydata=(1+P.TTickOffset)*sin(D.TTickValue(k));
   if strcmp(P.TTickScale,'radians')
      [n,d]=rat(P.TTickValue(k)/pi); % ticks as fractions
      if n==0
         Tstr='0';
      elseif n==1 && d==1
         Tstr='\pi';
      elseif n==-1 && d==1
         Tstr='-\pi';
      elseif n==1
         Tstr=['\pi/' num2str(d)];
      elseif n==-1
         Tstr=['-\pi/' num2str(d)];
      elseif d==1
         Tstr=[num2str(n) '\pi'];
      else
         Tstr=[num2str(n) '\pi/' num2str(d)];
      end
      P.TTickLabel{k}=Tstr;
   else % degrees
      P.TTickLabel{k}=[num2str(P.TTickValue(k)) '\circ'];
      if P.TTickValue(k)==-180
         P.TTickLabel{k}=['\pm' P.TTickLabel{k}(2:end)];
      end
   end
   D.HTTickLabel(k)=text(xdata,ydata,P.TTickLabel{k},...
      'Parent',HAxes,...
      'Color',P.TGridColor,...
      'FontName',P.FontName,...
      'FontSize',P.FontSize,...
      'FontWeight',P.FontWeight',...
      'HorizontalAlignment','center',...
      'VerticalAlignment','middle',...
      'Clipping','off',...
      'HandleVisibility','off',...
      'HitTest','off');
end
%--------------------------------------------------------------------------
function [P,D]=local_getTTickValue(P,D,S)          %#ok local_getTTickValue
% Get Theta Ticks
if strcmp(P.TTickScale,'degrees') % ticks are in degrees
   TTick=0:P.TTickDelta:360; % possible ticks
   if P.TLimit(1)>P.TLimit(2)
      idx=TTick<=P.TLimit(2)*180/pi | TTick>=P.TLimit(1)*180/pi; % keepers
   else
      idx=TTick>=P.TLimit(1)*180/pi & TTick<=P.TLimit(2)*180/pi; % keepers
   end
   P.TTickValue=TTick(idx);
   P.TTickValue=unique(rem(P.TTickValue,360)); % get unique ticks
   D.TTickValue=P.TTickValue*pi/180; % store ticks in radians
   if strcmp(P.TTickSign,'+-')
      tmp=P.TTickValue>=180;
      P.TTickValue(tmp)=P.TTickValue(tmp)-360;
   end
else                          % ticks are in radians
   TTick=(0:P.TTickDelta*180/pi:360)*pi/180; % possible ticks
   if P.TLimit(1)>P.TLimit(2)
      idx=TTick<=P.TLimit(2) | TTick>=P.TLimit(1); % keepers
   else
      idx=TTick>=P.TLimit(1) & TTick<=P.TLimit(2); % keepers
   end
   P.TTickValue=TTick(idx);
   P.TTickValue=unique(rem(P.TTickValue,2*pi)); % get unique ticks
   D.TTickValue=P.TTickValue; % store ticks in radians for plotting
   if strcmp(P.TTickSign,'+-')
      tmp=P.TTickValue>=pi;
      P.TTickValue(tmp)=P.TTickValue(tmp)-2*pi;
   end
end
D.TTickLabelN=length(D.TTickValue);
kmid=max(1,floor(length(P.TTickValue)/2));
P.RTickAngle=(D.TTickValue(kmid)+D.TTickValue(kmid+1))/2;
D.RTickAngle=P.RTickAngle;
if strcmp(P.TTickScale,'degrees')
   P.RTickAngle=P.RTickAngle*180/pi;
end
%--------------------------------------------------------------------------
function out=local_getDefaults                          % local_getDefaults
out.Style='cartesian';
out.Axis='on';
out.Border='on';
out.Grid='on';
out.RLimit=[0 1];
out.TLimit=[0 2*pi];
out.RTickUnits='';
out.TTickScale='degrees';
out.TDirection='ccw';
out.TZeroDirection='east';

out.BackgroundColor=get(0,'defaultaxescolor');
out.BorderColor=[0 0 0];
out.FontName=get(0,'defaultaxesfontname');
out.FontSize=get(0,'defaultaxesfontsize');
out.FontWeight=get(0,'defaultaxesfontweight');
out.TickLength=0.02;

out.RGridColor=get(0,'defaultaxesycolor');
out.RGridLineStyle=get(0,'defaultaxesgridlinestyle');
out.RGridLineWidth=get(0,'defaultaxeslinewidth');
out.RGridVisible='on';
out.RTickAngle=0;
out.RTickOffset=0.04;
out.RTickLabel='0|0.5|1.0';
out.RTickLabelHalign='center';
out.RTickLabelValign='middle';
out.RTickLabelVisible='on';
out.RTickValue=[0 .5 1];
out.RTickVisible='on';

out.TGridColor=get(0,'defaultaxesxcolor');
out.TGridVisible='on';
out.TGridLineStyle=get(0,'defaultaxesgridlinestyle');
out.TGridLineWidth=get(0,'defaultaxeslinewidth');
out.TTickOffset=0.08;
out.TTickDelta=15;
out.TTickDirection='in';
out.TTickLabel='';
out.TTickLabelVisible='on';
out.TTickSign='+-';
out.TTickValue=0:15:359;
out.TTickVisible='on';
%--------------------------------------------------------------------------
% Description of Plot Data Stored           Description of Plot Data Stored
%--------------------------------------------------------------------------
%
% VARIABLE        TYPE        DESCRIPTION
% D.TData         Cell        Raw input theta data
% D.RData         Cell        Raw input rho data
% D.RDataN        Cell        Normalized rho data
% D.LineColor     Cell        Line colors of plotted data
% D.LineStyle     Cell        Line styles of plotted data
% D.Marker        Cell        Markers of plotted data
% D.RLimitDiff    Double      Rho axis limit difference
% D.RMin          Double      Rho axis minimum
% D.RScale        Double      Power of ten scaling rho axis ticks
% D.RTickAngle    Double      Angle of rho tick labels in radians
% D.RTickRadius   Double      Radial position of rho tick labels
% D.RGridN        Double      Number of rho axis grid lines
% D.RTickLabelN   Double      Number of rho axis tick labels
% D.TGridN        Double      Number of theta axis grid lines      
% D.TTickLabelN   Double      Number of theta axis tick labels
% D.TTickValue    Double      Theta tick values in RADIANS
%--------------------------------------------------------------------------
% Description of Stored Handles               Description of Stored Handles
%--------------------------------------------------------------------------
%
% HANDLE          TYPE        DESCRIPTION
% D.HAPatch       Patch       Axis background and border
% D.HTTick        Line        Theta tick marks
% D.HRTick        Line        Rho tick marks
% D.HRGrid        Line        Rho grid
% D.HTTickLabel   Text        Theta Tick Labels
% D.HRTIckLabel   Text        Rho Tick Labels
% D.HLines        Line        Plotted Data
%
%--------------------------------------------------------------------------
