function shapeopt(controlswitch)
%
% shapeopt v 1.0
%
% This program optimizes the shape of the wire winding for an inductor using 2-D simulation
%
% Requirements: Matlab R11 or later versions.
%					 User must calculate number of turns, gap width required for desired inductance
%					 For more information visit http://thayer.dartmouth.edu/inductor
% Inputs:  The user chooses to either ignore or account for the radius of the inductor centerpost
%
% Ignoring the radius means that the radius is large compared to the window height and window
% width.  
%
% Accounting for the radius considers the three dimensional geometry.  This is a 
%  more accurate simulation for most inductors.
%
%	The user specifies a filename so that the program can save parameters in a 
%	.mat file.
%
%	The GUI window will open up.  The user specifies 11 parameters of the inductor in
%	editable text boxes. For more information on these parameters, click the box 
%	with the parameter name.
%
%Output :		3 plots
%	Figure(5) shows a cross section of the inductor, with the optimal placement of the wire
%	shaded in green.  The centerpost is on the left, and the gap is in the left center.
%
%	Figure(4) shows a contour plot of the magniture of the magnetic field B
%	Figure(3) shows a quiver plot of the magnetic field B
%
%	In the command window, the program will display the following data
%
%	Totalloss
%	Acloss
%	Dcloss
%	Total Area available for winding
%	Area used
%	Fraction of Area used 
%	Average magnitude of B squared
%	Ac Resistance
%	Dc Resistance
%	n = number of parallel strands per turn
%  Units are in milimeters, Hz, tesla, watts, and ohms

% Shapeopt by Jiankun Hu, Charles Sullivan, Rob Jensen, Dan Gestwick, Reid Hutchins
% Created at Thayer School of Engineering - Dartmouth College - Hanover NH 1997 - 2002

% Rob Jensen July 2002



%make sure only one copy of shapeopt is running at once 
h_figs = get(0, 'children');
fig_exists = 0;

if nargin < 1
   
	for fig = h_figs'
   	fig_exists = strcmp(get(fig, 'tag'), 'mainwindow');
   
  	 if fig_exists
    	  figure(fig);
    	  disp('Shapeopt already running');
    	  return;
  	 end
	end
end

if nargin <1, flag = 'startup';  end;
if nargin ==1, flag = 'guicalled'; end;

if nargin == 1
    if controlswitch == 6
        desoptions;     %displays a diagram of inductor dimensions in figure(2)
        return
    end
    
end
    
  switch flag    
   case 'startup'
            
            %this is a patch for Matlab 6 - previously the question boxes would
            %only display partially.  It is a known problem (see Mathworks website).
     
            filename = 0;
                
            
            for q = 1:100000
                dummy = log(q);
            end
    
            %this is a patch for Matlab 6 - previously the question boxes would
            %only display partially.  It is a known problem (see Mathworks website).
           
            % the user will choose to account for or ignore the centerpost radius.  
			quest_ans = questdlg('Do you want to account for the inductor centerpost radius?',...
   				'Choose',...
   				'Ignore',...
   				'Account',...
      			'cancel');
            

  
  			 if strcmp(quest_ans, 'cancel')
  			    return;
 			 elseif strcmp(quest_ans, 'Ignore')
  			    whichprogram = 0;                %means length of a turn must be entered on gui
 			 else
 			    whichprogram = 1;                %means diameter must be entered on gui
             end
             
             filename = initialize(whichprogram); %initializes variables or loads parameters
             gui(whichprogram, filename);         %start the gui
             
    case 'guicalled'
            guicalled(controlswitch);   % a recursive call to shapeopt to do a design
    end                                 % useful because UI callbacks in Matlab execute
                                        % at the scope of the command line, not where they should.
   
function fig = gui(whichprogram, filename)
% This is the machine-generated representation of a Handle Graphics object
% and its children.  Note that handle values may change when these objects
% are re-created. This may cause problems with any callbacks written to
% depend on the value of the handle at the time the object was saved.
% This problem is solved by saving the output as a FIG-file.
%
% To reopen this object, just type the name of the M-file at the MATLAB
% prompt. The M-file and its associated MAT-file must be on your path.
% 
% NOTE: certain newer features in MATLAB may not have been saved in this
% M-file due to limitations of this format, which has been superseded by
% FIG-files.  Figures which have been annotated using the plot editor tools
% are incompatible with the M-file/MAT-file format, and should be saved as
% FIG-files.

% Rob Jensen July 2002



% colormap for figure (machine generated)
mat0(1, :) = [		    0         0    0.5625];
mat0(2, :) = [        0         0    0.6250];
mat0(3, :) = [        0         0    0.6875];
mat0(4, :) = [        0         0    0.7500];
mat0(5, :) = [        0         0    0.8125];
mat0(6, :) = [        0         0    0.8750];
mat0(7, :) = [        0         0    0.9375];
mat0(8, :) = [        0         0    1.0000];
mat0(9, :) = [        0    0.0625    1.0000];
mat0(10, :) = [        0    0.1250    1.0000];
mat0(11, :) = [        0    0.1875    1.0000];
mat0(12, :) = [        0    0.2500    1.0000];
mat0(13, :) = [        0    0.3125    1.0000];
mat0(14, :) = [        0    0.3750    1.0000];
mat0(15, :) = [        0    0.4375    1.0000];
mat0(16, :) = [        0    0.5000    1.0000];
mat0(17, :) = [        0    0.5625    1.0000];
mat0(18, :) = [        0    0.6250    1.0000];
mat0(19, :) = [        0    0.6875    1.0000];
mat0(20, :) = [        0    0.7500    1.0000];
mat0(21, :) = [        0    0.8125    1.0000];
mat0(22, :) = [        0    0.8750    1.0000];
mat0(23, :) = [        0    0.9375    1.0000];
mat0(24, :) = [        0    1.0000    1.0000];
mat0(25, :) = [   0.0625    1.0000    1.0000];
mat0(26, :) = [   0.1250    1.0000    0.9375];
mat0(27, :) = [   0.1875    1.0000    0.8750];
mat0(28, :) = [   0.2500    1.0000    0.8125];
mat0(29, :) = [   0.3125    1.0000    0.7500];
mat0(30, :) = [   0.3750    1.0000    0.6875];
mat0(31, :) = [   0.4375    1.0000    0.6250];
mat0(32, :) = [   0.5000    1.0000    0.5625];
mat0(33, :) = [   0.5625    1.0000    0.5000];
mat0(34, :) = [   0.6250    1.0000    0.4375];
mat0(35, :) = [   0.6875    1.0000    0.3750];
mat0(36, :) = [   0.7500    1.0000    0.3125];
mat0(37, :) = [   0.8125    1.0000    0.2500];
mat0(38, :) = [   0.8750    1.0000    0.1875];
mat0(39, :) = [   0.9375    1.0000    0.1250];
mat0(40, :) = [   1.0000    1.0000    0.0625];
mat0(41, :) = [   1.0000    1.0000         0];
mat0(42, :) = [   1.0000    0.9375         0];
mat0(43, :) = [   1.0000    0.8750         0];
mat0(44, :) = [   1.0000    0.8125         0];
mat0(45, :) = [   1.0000    0.7500         0];
mat0(46, :) = [   1.0000    0.6875         0];
mat0(47, :) = [   1.0000    0.6250         0];
mat0(48, :) = [   1.0000    0.5625         0];
mat0(49, :) = [   1.0000    0.5000         0];
mat0(50, :) = [   1.0000    0.4375         0];
mat0(51, :) = [   1.0000    0.3750         0];
mat0(52, :) = [   1.0000    0.3125         0];
mat0(53, :) = [   1.0000    0.2500         0];
mat0(54, :) = [   1.0000    0.1875         0];
mat0(55, :) = [   1.0000    0.1250         0];
mat0(56, :) = [   1.0000    0.0625         0];
mat0(57, :) = [   1.0000         0         0];
mat0(58, :) = [   0.9375         0         0];
mat0(59, :) = [   0.8750         0         0];
mat0(60, :) = [   0.8125         0         0];
mat0(61, :) = [   0.7500         0         0];
mat0(62, :) = [   0.6875         0         0];
mat0(63, :) = [   0.6250         0         0];
mat0(64, :) = [   0.5625         0         0];


%  gui callbacks and createfunctions.  
mat1 =  '';
mat10 =  ' ';
mat11 =  'filename = get(findobj(''tag'', ''mainwindow''), ''userdata''); load(filename); f_handle = gcbo;  set(gcbo, ''String'', num2str(Settings(10))); Handle(4) = gcbo; save(filename, ''Handle'', ''Settings'',  ''-APPEND''); ';
mat12 =  'filename = get(findobj(''tag'', ''mainwindow''), ''userdata''); load(filename); gw_handle = gcbo;  set(gcbo, ''String'', num2str(Settings(7))); Handle(3) = gcbo; save(filename, ''Handle'', ''Settings'',  ''-APPEND'');';
mat13 =  '                                                                                                                               ';
mat14 =  'filename = get(findobj(''tag'', ''mainwindow''), ''userdata''); load(filename);  set(gcbo, ''String'', num2str(Settings(2))); Handle(2) = gcbo; save(filename, ''Handle'', ''Settings'',  ''-APPEND'');';
mat15 =  '                                                                                                                                                                 ';
mat16 =  'filename = get(findobj(''tag'', ''mainwindow''), ''userdata''); load(filename);  bw_handle = gcbo;  set(gcbo, ''String'', num2str(Settings(1)));  Handle(1) = gcbo;  save(filename, ''Handle'', ''Settings'',  ''-APPEND'');        ';
mat17 =  'helpdlg(''Window Width - Bw. Width of the winding window in milimeters--the dimension parallel to the leg with the gap in it.  For example, with the gapped centerpost, this is the length of the centerpost (not including the parts beyond the window).  '');';
mat18 =  'helpdlg(''Window Height - Hw  Height of the winding window in milimeters--the dimension perpendicular to the leg with the gap in it.  For example, with the gapped centerpost, this is the distance from the centerpost to the inside of the outer core leg.  '');';
mat19 =  'helpdlg(''Gap Width - This is the length of the air gap, normally in the centerpost.   All distances are in milimeters.  For a diagram, push "show dimensions."'');';
mat2 =  'close';
mat20 =  'CreateStruct.WindowStyle=''replace''; CreateStruct.Interpreter=''tex''; H = msgbox(''.                                                                                                                                                               .                                                                   Frequency - The frequency of sinusoidal current in the winding.  For nonsinusoidal waveforms "equivalent frequency" can be used.  f_{eq} = 1/(2pi)*rms{dI/dt}/(I_{rms}), where I_{rms} is total rms current including ac and dc components.  For a sinusoid with dc offset, f_{eq} = f_{sine} * sqrt( I_{ac}^2/(I_{ac}^2 + I_{dc}^2)) where I_{ac} is the rms value of the ac component.  For more detail, see Appendix C in C.R. Sullivan.  "Optimal Choice for Number of Strands in a Litz-Wire Transformer Winding." IEEE Transactions on Power Electronics, 14(2), March 1999, pp. 283-291, available at  http://thayer.dartmouth.edu/inductor/litzj.pdf'', ''Help'', ''help'', CreateStruct); set(H, ''resize'', ''on'')';
mat21 =  'helpdlg(''Strand Diameter - Diameter of one individual strand in milimeters.  The winding is assumed to be wound with litz wire comprising a number of these strands; the program will calculate the optimal number of strands to be used.  AWG 44, 0.05 mm diameter, is often a good compromise between cost and loss.'');';
mat22 =  'helpdlg(''Number of turns in the winding.  This must be calculated based on core loss and saturations characteristics before running this winding optimization program.'');';
mat23 =  'helpdlg(''Packing Factor - The fraction of cross-sectional area that is copper in a space filled with wire is Fp*pi/4.  In other words, Fp is the packing factor relative to perfect square packing of cylinders. It accounts for both insulation space and for imperfect packing.   A typical number for litz wire is Fp = 0.35, resulting in 28% of winding area occupied by copper.'');';mat24 =  'helpdlg(''Total current - This is the current through a turn of wire - Irms in amps '');';
mat25 =  'delete(findobj(''tag'',''mainwindow'')); shapeopt;                               ';
mat26 =  '                                                                    ';
if ~whichprogram
    mat27 =  'filename = get(findobj(''tag'', ''mainwindow''), ''userdata''); load(filename);  r_handle = gcbo; set(gcbo, ''String'', num2str(Settings(17))); Handle(11) = gcbo;  save(filename, ''Handle'', ''Settings'',  ''-APPEND''); ';
    mat28 =  'helpdlg(''This is the length of a turn in milimeters.   For better accuracy, chose the axisymmetric model and specify the centerpost radius (ie run the program taking into account the centerpost radius).'');';
    mat29 =  [...
      'filename = get(findobj(''tag'', ''mainwindow''), ''userdata'');  load(filename);'...
      'bwstring = get(Handle(1), ''string'');'...
      'bw = eval(bwstring)/1000;'...
      'gwstring = get(Handle(3), ''string'');'...
      'gw = eval(gwstring)/1000;'...
      'hwstring = get(Handle(2), ''String'');'...
      'hw = eval(hwstring)/1000;'...
      'N_string = get(Handle(6), ''string'');'...
      'N = eval(N_string);'...
      'I_string = get(Handle(8), ''string'');'...
      'I_total = eval(I_string)*N;'...
      'bbstring = get(Handle(12), ''string'');'...
       'bb = eval(bbstring)/1000;'...
       'bhstring = get(Handle(13), ''string'');'...
       'bh = eval(bhstring)/1000;'...
       'lstring = get(Handle(11), ''string'');'...
       'l = eval(lstring)/1000;'...
      'if((bw ~= Settings(1)) | (l ~= Settings(17)) | (bb ~= Settings(19)) | (bh~= Settings(20)) | (gw ~= Settings(7)) | (I_total ~= Settings(9)) | (N ~= Settings(12))  | Settings(15)),'...
      'controlswitch = 0;'...
   	'elseif  hw ~= Settings(2),'...
      'controlswitch = 1;'...
      'else, '...
      'controlswitch = 2;'...
   	'end,'...
      'shapeopt(controlswitch);'];
else
    mat27 =  'filename = get(findobj(''tag'', ''mainwindow''), ''userdata''); load(filename);  r_handle = gcbo; set(gcbo, ''String'', num2str(Settings(16))); Handle(9) = gcbo;  save(filename, ''Handle'', ''Settings'',  ''-APPEND''); ';
    mat28 =  'helpdlg(''Core centerpost diameter in milimeters.  For a square centerpost, use d = perimeter/pi.'');';
    mat29 =  [...
      ' filename = get(findobj(''tag'', ''mainwindow''), ''userdata''); load(filename);'...
      'bwstring = get(Handle(1), ''string'');'...
      'bw = eval(bwstring)/1000;'...
      'gwstring = get(Handle(3), ''string'');'...
      'gw = eval(gwstring)/1000;'...
      'hwstring = get(Handle(2), ''String'');'...
      'hw = eval(hwstring)/1000;'...
      'N_string = get(Handle(6), ''string'');'...
      'N = eval(N_string);'...
      'I_string = get(Handle(8), ''string'');'...
      'I_total = eval(I_string)*N;'...
      'bbstring = get(Handle(12), ''string'');'...
       'bb = eval(bbstring)/1000;'...
       'bhstring = get(Handle(13), ''string'');'...
       'bh = eval(bhstring)/1000;'...
       'rstring = get(Handle(9), ''string'');'...
       'r = eval(rstring)/(2*1000);'...
      'if((bw ~= Settings(1)) | (r ~= Settings(16)) | (bb ~= Settings(19)) | (bh~= Settings(20)) | (gw ~= Settings(7)) | (I_total ~= Settings(9)) | (N ~= Settings(12))  | Settings(15)),'...
      'controlswitch = 0;'...
   	'elseif  hw ~= Settings(2),'...
      'controlswitch = 1;'...
      'else, '...
      'controlswitch = 2;'...
   	'end,'...
      'shapeopt(controlswitch);'];
end

mat3 =  'filename = get(findobj(''tag'', ''mainwindow''), ''userdata'');  load(filename);  I_handle = gcbo;  set(gcbo, ''String'', num2str(Settings(9))); Handle(8) = gcbo; save(filename, ''Handle'', ''Settings'',  ''-APPEND'');';
mat4 =  '                                 ';
mat5 =  'filename = get(findobj(''tag'', ''mainwindow''), ''userdata'');  load(filename); Fp_handle = gcbo; set(gcbo, ''String'', num2str(Settings(13))); Handle(7) = gcbo;  save(filename, ''Handle'', ''Settings'',  ''-APPEND'');';
mat6 =  '                                 ';
mat7 =  'filename = get(findobj(''tag'', ''mainwindow''), ''userdata'');  load(filename);  N_handle = gcbo;  set(gcbo, ''String'', num2str(Settings(12)));   Handle(6) = gcbo;  save(filename, ''Handle'', ''Settings'',  ''-APPEND'');  ';
mat8 =  '                                 ';
mat9 =  'filename = get(findobj(''tag'', ''mainwindow''), ''userdata''); load(filename);  d_handle = gcbo;  set(gcbo, ''String'', num2str(Settings(11)));  Handle(5) = gcbo;  save(filename, ''Handle'', ''Settings'',  ''-APPEND'');  ';

runbuttonmat = ['filename = get(findobj(''tag'', ''mainwindow''), ''userdata'');  load(filename);'...
          'if doption(1) ==1,'...
            'shapeopt(7);'...
            'elseif doption(1) ==2,'...
            'shapeopt(8);'...
            'else,'...
            mat29...
           'end'];
% definte parts of the user interface


% the figure itself
h0 = figure('Units','points', ...
	'Color',[0.8 0.8 0.8], ...
	'Colormap',mat0, ...
   'CreateFcn',mat1, ...
   'DeleteFcn', 'get(findobj(''tag'',''designoptions''),''userdata''); close; shapeopt(5)',...
	'Position',[200 40 450 490], ...
	'Tag','mainwindow', ...
	'ToolBar','none', ...
   'FileName','', ...
   'Name',filename,...
   'NumberTitle','off',...
   'UserData',filename);

showdimensions

% edittext boxes
bb = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[1 1 1], ...
	'Callback','', ...
	'CreateFcn','filename = get(findobj(''tag'', ''mainwindow''), ''userdata''); load(filename);  set(gcbo, ''String'', num2str(Settings(19))); Handle(12) = gcbo; save(filename, ''Handle'', ''Settings'',  ''-APPEND'');', ...
	'ListboxTop',0, ...
	'Position',[125.25 369.44 45 15], ...
	'String','1', ...
	'Style','edit', ...
   'Tag','EditText10');
bh = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[1 1 1], ...
	'Callback','', ...
	'CreateFcn','filename = get(findobj(''tag'', ''mainwindow''), ''userdata''); load(filename);  set(gcbo, ''String'', num2str(Settings(20))); Handle(13) = gcbo; save(filename, ''Handle'', ''Settings'',  ''-APPEND'');', ...
	'ListboxTop',0, ...
	'Position',[125.25 336.22 45 15], ...
	'String','1', ...
	'Style','edit', ...
   'Tag','EditText11');  

h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[1 1 1], ...
	'Callback',mat15, ...
	'CreateFcn',mat16, ...
	'ListboxTop',0, ...
	'Position',[125.25 303 45 15], ...
	'String','0.01', ...
	'Style','edit', ...
   'Tag','EditText1');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[1 1 1], ...
	'Callback',mat13, ...
	'CreateFcn',mat14, ...
	'ListboxTop',0, ...
	'Position',[125.25 269.53125 45 15], ...
	'String','0.01', ...
	'Style','edit', ...
   'Tag','EditText2');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[1 1 1], ...
	'Callback','                                                             ', ...
	'CreateFcn',mat12, ...
	'ListboxTop',0, ...
	'Position',[125.25 236.0625 45 15], ...
	'String','0.001', ...
	'Style','edit', ...
   'Tag','EditText3');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[1 1 1], ...
	'Callback',mat10, ...
	'CreateFcn',mat11, ...
	'ListboxTop',0, ...
	'Position',[125.25 202.59375 45 15], ...
	'String','0', ...
	'Style','edit', ...
   'Tag','EditText4');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[1 1 1], ...
	'Callback',mat8, ...
	'CreateFcn',mat9, ...
	'ListboxTop',0, ...
	'Position',[125.25 169.125 45 15], ...
	'String','0.0001', ...
	'Style','edit', ...
	'Tag','diameter');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[1 1 1], ...
	'Callback',mat6, ...
	'CreateFcn',mat7, ...
	'ListboxTop',0, ...
	'Position',[125.25 135.65625 45 15], ...
	'String','1', ...
	'Style','edit', ...
   'Tag','EditText6');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[1 1 1], ...
	'Callback',mat4, ...
	'CreateFcn',mat5, ...
	'ListboxTop',0, ...
	'Position',[125.25 102.1875 45 15], ...
	'String','0.85', ...
	'Style','edit', ...
	'Tag','EditText7');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[1 1 1], ...
	'Callback','                                                      ', ...
	'CreateFcn',mat3, ...
	'ListboxTop',0, ...
	'Position',[125.25 68.71875 45 15], ...
	'String','1', ...
	'Style','edit', ...
	'Tag','EditText8');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[1 1 1], ...
	'Callback',mat26, ...
	'CreateFcn',mat27, ...
	'ListboxTop',0, ...
	'Position',[125.25 35.25 45 15], ...
	'String','1', ...
	'Style','edit', ...
   'Tag','EditText9');

%display the units on the gui - text boxes below

h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[.8 .8 .8], ...
	'Position',[175.25 298 25 15], ...
	'String','mm', ...
	'Style','text', ...
	'Tag','Units1');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[.8 .8 .8], ...
	'Position',[175.25 265.53125 25 15], ...
	'String','mm', ...
	'Style','text', ...
	'Tag','Units2');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[.8 .8 .8], ...
	'Position',[175.25 231.0625 25 15], ...
	'String','mm', ...
	'Style','text', ...
	'Tag','Units3');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[.8 .8 .8], ...
	'Position',[175.25 197.59375 25 15], ...
	'String','Hz', ...
	'Style','text', ...
	'Tag','Units4');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[.8 .8 .8], ...
	'Position',[175.25 164.125 25 15], ...
	'String','mm', ...
	'Style','text', ...
	'Tag','Units5');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[.8 .8 .8], ...
	'Position',[175.25 130.65625 25 15], ...
	'String','', ...
	'Style','text', ...
	'Tag','Units6');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[.8 .8 .8], ...
	'Position',[175.25 97.1875 25 15], ...
	'String','', ...
	'Style','text', ...
	'Tag','Units7');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[.8 .8 .8], ...
	'Position',[175.25 63.71875 25 15], ...
	'String','A', ...
	'Style','text', ...
	'Tag','Units8');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[.8 .8 .8], ...
	'Position',[175.25 30.25 25 15], ...
	'String','mm', ...
	'Style','text', ...
	'Tag','Units9');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[.8 .8 .8], ...
	'Position',[175.25 331 25 15], ...
	'String','mm', ...
	'Style','text', ...
	'Tag','Units10');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[.8 .8 .8], ...
	'Position',[175.25 364 25 15], ...
	'String','mm', ...
	'Style','text', ...
	'Tag','Units11');

% static text boxes
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[.8 .8 .8], ...
	'ListboxTop',0, ...
	'Position',[22.5 430 434.25 33.75], ...
   'String','Shapeopt - Optimization of Inductor Shape', ...
   'fontsize', 10,...
   'fontweight', 'bold',...
	'Style','text', ...
	'Tag','StaticText16');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[.8 .8 .8], ...
	'ListboxTop',0, ...
	'Position',[124.5 400 52.5 30.75], ...
	'String','Enter Parameters', ...
	'Style','text', ...
	'Tag','StaticText26');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[.8 .8 .8], ...
	'ListboxTop',0, ...
	'Position',[27.75 400 66 29.25], ...
	'String','Paramameters - Click for Description', ...
	'Style','text', ...
	'Tag','StaticText27');


% pushbuttons
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[.75 .75 .75], ...
	'Callback',mat17, ...
	'ListboxTop',0, ...
	'Position',[24.375 300 80 20], ...
	'String','Window Breadth', ...
	'Tag','Pushbutton3');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[.75 .75 .75], ...
	'Callback',mat18, ...
	'ListboxTop',0, ...
	'Position',[23.25 265.625 80 20], ...
	'String','Window Height', ...
	'Tag','Pushbutton4');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[.75 .75 .75], ...
	'Callback',mat19, ...
	'ListboxTop',0, ...
	'Position',[23.25 233 80 20], ...
	'String','Gap Length', ...
	'Tag','Pushbutton5');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[.75 .75 .75], ...
	'Callback',mat20, ...
	'ListboxTop',0, ...
	'Position',[23.25 198.375 80 20], ...
	'String','Frequency', ...
	'Tag','Pushbutton6');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[.75 .75 .75], ...
	'Callback',mat21, ...
	'ListboxTop',0, ...
	'Position',[22.875 165.75 80 20], ...
	'String','Strand Diameter', ...
	'Tag','Pushbutton7');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[.75 .75 .75], ...
	'Callback','shapeopt(9)', ...
	'ListboxTop',0, ...
	'Position',[210 170 30 15], ...
	'String','AWG', ...
	'Tag','AWG');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[.75 .75 .75], ...
	'Callback',mat22, ...
	'ListboxTop',0, ...
	'Position',[22.875 133.125 80 20], ...
	'String','Number of Turns', ...
	'Tag','Pushbutton8');
h1 = uicontrol('Parent',h0, ...
   'Units','points', ...
   'BackgroundColor',[.75 .75 .75], ...
   'Callback',mat23, ...
	'ListboxTop',0, ...
	'Position',[23.25 97.5 80 20], ...
	'String','Packing Factor', ...
	'Tag','Pushbutton9');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[.75 .75 .75], ...
	'Callback',mat24, ...
	'ListboxTop',0, ...
	'Position',[24.375 65.625 80 20], ...
	'String','Current (RMS)', ...
	'Tag','Pushbutton10');

if ~whichprogram        %length of a turn
    h1 = uicontrol('Parent',h0, ...
	    'Units','points', ...
	    'BackgroundColor',[.75 .75 .75], ...
	    'Callback',mat28, ...
	    'ListboxTop',0, ...
	    'Position',[24.75 33.75 80 20], ...
	    'String','Length of a turn', ...
	    'Tag','Pushbutton10');

else                    %radius
    h1 = uicontrol('Parent',h0, ...
	    'Units','points', ...
	    'BackgroundColor',[.75 .75 .75], ...
	    'Callback',mat28, ...
	    'ListboxTop',0, ...
	    'Position',[22.5 35.25 80 20], ...
	    'String','Centerpost Diameter', ...
	    'Tag','Pushbutton10');
end

h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[.75 .75 .75], ...
	'Callback','helpdlg(''The breadth of the bobbin is ususally slightly smaller than the breadth of the winding window (perpeneidular to the gap). '');', ...
	'ListboxTop',0, ...
	'Position',[24.75 368.75 80 20], ...
	'String','Bobbin breadth', ...
	'Tag','Pushbutton13');

h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[.75 .75 .75], ...
	'Callback','helpdlg(''The height of the winding window is defined as the distance from the outer edge of the winding window to the near edge of the bobbin. '');', ...
	'ListboxTop',0, ...
	'Position',[24.75 334 80 20], ...
	'String','Bobbin height', ...
	'Tag','Pushbutton14');

a = 70; % used to move the boxes on the right up and down for fine-tuning of display

% the key buttons on the right of the gui
% utility buttons
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[.75 .75 .75], ...
	'Callback',mat25, ...
	'ListboxTop',0, ...
	'Position',[328 288.75+a 87.75 29.25], ...
	'String','Load/ New', ...
	'Tag','Pushbutton11');
h1 = uicontrol('Parent',h0, ...
   'Units','points', ...
   'BackgroundColor',[.75 .75 .75], ...
   'Callback',mat2, ...
	'ListboxTop',0, ...
	'Position',[327.25 250.5+a 89.25 26.25], ...
	'String','Close/Save', ...
   'Tag','Pushbutton2');

% design buttons
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[.75 .75 .75], ...
	'Callback','shapeopt(6)', ...
	'ListboxTop',0, ...
	'Position',[220.75 288.75+a 87 27.75], ...
	'String','Options', ...
	'Tag','Pushbutton15');
h1 = uicontrol('Parent',h0, ...
	'Units','points', ...
	'BackgroundColor',[.75 .75 .75], ...
	'Callback',runbuttonmat, ...
	'ListboxTop',0, ...
	'Position',[220.75 249.75+a 88.5 27], ...
	'String','Run', ...
    'fontweight', 'bold',...
	'Tag','Pushbutton1');

%radio buttons

h1 = uicontrol(h0,'style','radio',...
    'Units','points', ...
    'BackgroundColor',[0.8 0.8 0.8], ...
    'callback','set(findobj(''tag'', ''centergaps''), ''value'', 1); set(findobj(''tag'', ''outergaps''), ''value'', 0); set(findobj(''tag'', ''allgapsround''), ''value'', 0); set(findobj(''tag'', ''allgapssquare''), ''value'', 0); ',...
    'string', 'Gapped Centerpost Only',...
    'tag', 'centergaps',...
    'Value',[1],...
    'position', [220.75 280 200 20]);
    
h1 = uicontrol(h0,'style','radio',...
    'Units','points', ...
    'BackgroundColor',[0.8 0.8 0.8], ...
    'callback','set(findobj(''tag'', ''centergaps''), ''value'', 0); set(findobj(''tag'', ''outergaps''), ''value'', 1); set(findobj(''tag'', ''allgapsround''), ''value'', 0); set(findobj(''tag'', ''allgapssquare''), ''value'', 0); ',...
    'string', 'Outer Legs Gapped Only',...
    'tag', 'outergaps',...
    'position', [220.75 260 200 20]);
    
h1 = uicontrol(h0,'style','radio',...
    'Units','points', ...
    'BackgroundColor',[0.8 0.8 0.8], ...
    'callback','set(findobj(''tag'', ''centergaps''), ''value'', 0); set(findobj(''tag'', ''outergaps''), ''value'', 0); set(findobj(''tag'', ''allrapsround''), ''value'', 1); set(findobj(''tag'', ''allgapssquare''), ''value'', 0); ',...
    'string', 'All Legs Gapped (round centerpost)',...
    'tag', 'allgapsround',...
    'position', [220.75 240 200 20]);
h1 = uicontrol(h0,'style','radio',...
    'Units','points', ...
    'BackgroundColor',[0.8 0.8 0.8], ...
    'callback','set(findobj(''tag'', ''centergaps''), ''value'', 0); set(findobj(''tag'', ''outergaps''), ''value'', 0); set(findobj(''tag'', ''allgapsround''), ''value'', 0); set(findobj(''tag'', ''allgapssquare''), ''value'', 1); ',...
    'string', 'All Legs Gapped (square centerpost)',...
    'tag', 'allgapssquare',...
    'position', [220.75 220 200 20]);

%  display the menu labels for changing grid resolution and number of images

h_menu_options = uimenu(h0, 'label', 'Options');

h_menu_setgrid = uimenu(h_menu_options, 'label', 'Change Grid Reslution');

h_menu_grid10 = uimenu('parent', h_menu_setgrid, 'label', '10x10', 'callback', 'filename = get(findobj(''tag'', ''mainwindow''), ''userdata'');  Settings(15) = 1; Settings(5) = 10; Settings(6) = 10; save(filename, ''Settings'', ''-APPEND''); disp(''Grid is now 10x10''); ');
h_menu_grid20 = uimenu('parent', h_menu_setgrid, 'label', '20x20 (default)', 'callback', 'filename = get(findobj(''tag'', ''mainwindow''), ''userdata'');  Settings(15) = 1; Settings(5) = 20; Settings(6) = 20; save(filename, ''Settings'', ''-APPEND'');disp(''Grid is now 20x20''); ');
h_menu_grid30 = uimenu('parent', h_menu_setgrid, 'label', '30x30', 'callback', 'filename = get(findobj(''tag'', ''mainwindow''), ''userdata'');  Settings(15) = 1; Settings(5) = 30; Settings(6) = 30; save(filename, ''Settings'', ''-APPEND'');disp(''Grid is now 30x30''); ');
h_menu_grid40 = uimenu('parent', h_menu_setgrid, 'label', '40x40', 'callback', 'filename = get(findobj(''tag'', ''mainwindow''), ''userdata'');  Settings(15) = 1; Settings(5) = 40; Settings(6) = 40; save(filename, ''Settings'', ''-APPEND'');disp(''Grid is now 40x40''); ');
h_menu_grid50 = uimenu('parent', h_menu_setgrid, 'label', '50x50', 'callback', 'filename = get(findobj(''tag'', ''mainwindow''), ''userdata'');  Settings(15) = 1; Settings(5) = 50; Settings(6) = 50; save(filename, ''Settings'', ''-APPEND'');disp(''Grid is now 50x50''); ');
h_menu_grid60 = uimenu('parent', h_menu_setgrid, 'label', '60x60', 'callback', 'filename = get(findobj(''tag'', ''mainwindow''), ''userdata'');  Settings(15) = 1; Settings(5) = 60; Settings(6) = 60; save(filename, ''Settings'', ''-APPEND'');disp(''Grid is now 60x60''); ');
h_menu_grid70 = uimenu('parent', h_menu_setgrid, 'label', '70x70', 'callback', 'filename = get(findobj(''tag'', ''mainwindow''), ''userdata'');  Settings(15) = 1; Settings(5) = 70; Settings(6) = 70; save(filename, ''Settings'', ''-APPEND'');disp(''Grid is now 70x70''); ');
h_menu_grid80 = uimenu('parent', h_menu_setgrid, 'label', '80x80', 'callback', 'filename = get(findobj(''tag'', ''mainwindow''), ''userdata'');  Settings(15) = 1; Settings(5) = 80; Settings(6) = 80; save(filename, ''Settings'', ''-APPEND'');disp(''Grid is now 80x80''); ');
h_menu_grid90 = uimenu('parent', h_menu_setgrid, 'label', '90x90', 'callback', 'filename = get(findobj(''tag'', ''mainwindow''), ''userdata'');  Settings(15) = 1; Settings(5) = 90; Settings(6) = 90; save(filename, ''Settings'', ''-APPEND'');disp(''Grid is now 90x90''); ')';
h_menu_grid100 = uimenu('parent', h_menu_setgrid, 'label', '100x100', 'callback', 'filename = get(findobj(''tag'', ''mainwindow''), ''userdata'');  Settings(15) = 1; Settings(5) = 100; Settings(6) = 100; save(filename, ''Settings'', ''-APPEND'');disp(''Grid is now 100x100''); ')';

h_menu_setimage = uimenu(h_menu_options, 'label', 'Change Image Matrix');

h_menu_image3 = uimenu('parent', h_menu_setimage, 'label', '3x3', 'callback', 'filename = get(findobj(''tag'', ''mainwindow''), ''userdata'');  Settings(15) = 1; Settings(3) = 1; Settings(4) = 1; save(filename, ''Settings'', ''-APPEND''); disp(''Switched number of image units to 3x3''); ');
h_menu_image5 = uimenu('parent', h_menu_setimage, 'label', '5x5 (default)', 'callback', 'filename = get(findobj(''tag'', ''mainwindow''), ''userdata'');  Settings(15) = 1; Settings(3) = 2; Settings(4) = 2; save(filename, ''Settings'', ''-APPEND''); disp(''Switched number of image units to 5x5''); ');
h_menu_image7 = uimenu('parent', h_menu_setimage, 'label', '7x7', 'callback', 'filename = get(findobj(''tag'', ''mainwindow''), ''userdata'');  Settings(15) = 1; Settings(3) = 3; Settings(4) = 3; save(filename, ''Settings'', ''-APPEND''); disp(''Switched number of image units to 7x7''); ');


if nargout > 0, fig = h0; end


function filename = initialize(whichprogram)
%Called by shapeopt to initialize the variables or load variables from a file

% Rob Jensen July 2002

filename = 0;

Xpmult = 2;
andyswitch =1;  % can turn off the extended flux path calculations
                % and run the old version of shapeopt by setting this to 0

    
            for q = 1:200000
                dummy = log(q);
            end

    %this is a patch for Matlab 6 - previously the question boxes would
    %only display partially.  It is a known problem (see Mathworks website).

	quest_ans = questdlg('Would you like to start a new file or open an existing one?',...
   	'Choose',...
   	'New',...
   	'Open',...
      'cancel');
  

  
%  If a new file is chosen, initialize user-chosen variables to 0
if (strcmp(quest_ans, 'New') | strcmp(quest_ans, 'cancel'))
   
   	  while filename == 0   % force the user to enter a filename for saving parameters
   		[filename, pathname] = uiputfile('*.*', 'Save as');
     	assignin('base', 'filename', filename);
      end
     
      if ischar(filename)
         cd(pathname);

   		bw = 0;
		hw = 0;
		gw = 0;
		f  = 0;
		d  = 0;                        %initialization
		N  = 0;
        Fp = 0;
        I_total = 0;
        r  = 0;
        wirelength = 0;
        bb = 0;
        bh = 0;
    
   
   		Settings(1) = bw;    
		Settings(2) = hw;  
		Settings(7) = gw;      
		Settings(9) = I_total;         %user-defined settings (later)
		Settings(10) = f;                                                                                                                                                                                                                                                                                                                                                                                                                                                  
		Settings(11) = d;                                                                                                                                                                                                                                                                                                                                                                                                                                                  
		Settings(12) = N;                                                                                                                                                                                                                                                                                                                                                                                                                                                  
        Settings(13) = Fp; 
        Settings(16) = r;
        Settings(17) = wirelength;
        Settings(18) = whichprogram;
        Settings(19) = bb;
        Settings(20) = bh;
        
   
   		bw_handle = 0;
		hw_handle = 0;
		gw_handle = 0;
		f_handle  = 0;				%the handles for certain editable text boxes
		d_handle  = 0;
		N_handle  = 0;
		Fp_handle = 0;
        I_handle = 0;
        r_handle = 0;
        prog_handle = 0;
        wireL_handle =0;
        bb_handle = 0;
        bh_handle = 0;
        
         
   		Handle(1) = bw_handle;
		Handle(2) = hw_handle;
		Handle(3) = gw_handle;
		Handle(4) = f_handle;
		Handle(5) = d_handle;       %initialization
		Handle(6) = N_handle;
		Handle(7) = Fp_handle;
        Handle(8) = I_handle;
        Handle(9) = r_handle;
        Handle(10) = prog_handle;
        Handle(11) = wireL_handle;
        Handle(12) = bb_handle;
        Handle(13) = bh_handle;
      
         
		Settings(3) =   2;          %unit of x                            
		Settings(4) =   2;          %unit of y                            
		Settings(5) =   20;         %xdiv                                      
		Settings(6) =   20;         %ydiv                                      
		Settings(8) =   10;         %gapdiv                                         
		Settings(14)=   .1;         %step                                      
		Settings(15)=   1;          %program not run yet (0 if run)  
        Settings(21) = Xpmult;       
        
        windheight = 0;
        Settings(22) = windheight;
   
        doption(1) = 3;
        doption(2) = 1;
        doption(3) = 1;
        doption(4) = 1;
        
   		save(filename, 'Settings', 'doption','Handle', 'andyswitch');
      else          %if the user tries not entering a filename
         disp('Error, a filename must be entered to start the program');
      end
      
   elseif strcmp(quest_ans, 'Open') 
   %open the file and extract the user-defined variables
      
      while filename == 0
     		[filename, pathname] = uigetfile('*.*', 'Choose file to load');       
     		assignin('base', 'filename', filename);
      end     
           
      if ischar(filename)
      	load(filename);
   		cd(pathname);

   		bw_handle = 0;
		hw_handle = 0;
		gw_handle = 0;
		f_handle  = 0;
		d_handle  = 0;
		N_handle  = 0;
		Fp_handle = 0;                          %initialize
        I_handle = 0;
        r_handle = 0;
        prog_handle = 0;
        wireL_handle = 0;
        bb_handle = 0;
        bh_handle = 0;
        
   		Handle(1) = bw_handle;
		Handle(2) = hw_handle;
		Handle(3) = gw_handle;
		Handle(4) = f_handle;
		Handle(5) = d_handle;
		Handle(6) = N_handle;
		Handle(7) = Fp_handle;                  %intialize
        Handle(8) = I_handle;
        Handle(9) = r_handle;
        Handle(10) = prog_handle;
        Handle(11) = wireL_handle;
        Handle(12) = bb_handle;
        Handle(13) = bh_handle;
         
         
        Settings(1) = Settings(1)*1000;         %scale to put in mm for gui
        Settings(2) = Settings(2)*1000;
        Settings(7) = Settings(7)*1000;
        Settings(11) = Settings(11)*1000;
        Settings(16) = Settings(16)*1000*2;     %scale to diameter here
        Settings(17) = Settings(17)*1000;       
        Settings(19) = Settings(19)*1000;
        Settings(20) = Settings(20)*1000;
        Settings(21) = Xpmult;  
        
        windheight = 0;
        Settings(22) = windheight;
           
        Settings(3) =2;                    %unit of x                            
	    Settings(4) =2;                    %unit of y                            
		Settings(5) =20;              		%xdiv                                      
		Settings(6) = 20;             		%ydiv                                      
		Settings(8) =10;         			%gapdiv                                         
		Settings(14) =.1;             		%step                                      
		Settings(15) = 1;                  	%program not run yet     
        
        doption(1) = 3;
        doption(2) = 1;
        doption(3) = 1;
        doption(4) = 1;
   
   		cd(pathname);
      
      Settings(15) = 1;                  	%program not run yet  
      Settings(9) = Settings(9)/Settings(12); %I_total is stored, reset to I(in one turn)
   		Settings(18) = whichprogram;
   		save(filename,'Settings', 'Handle','doption', 'andyswitch','-APPEND');
      else              %if the user tries not entering a filename
         disp('Error, a filename must be entered to start the program');
      end 
   else                 %%if the user tries not entering a filename
      disp('Error, a filename must be entered to start the program');
end
   

function guicalled(control)
%guicalled
%Rob Jensen Aug 7, 2000
%guicalled runs the COOS algorithims as they are needed.  This is called by the Optimize Shape button
%The control flow here greatly reduces the time needed for multiple optimization runs
%because the algorithims are only called as needed.  
%
% Control == 7 filled bobbin design (no optimization)
% Control == 8 design by picking strands and doing a rectangular shape



filename = get(findobj('tag','mainwindow'), 'userdata');                                                                                   
load(filename);

%The files all need to be run once before they can start to be done out of order
%Setting 15 is 1 if they have not all been run yet, and it is 0
%if they have all been run.  This only applies to the optimize shape button

%parameters are divided by 1000 to convert to meters for COOS

Settings(1) = eval(get(Handle(1), 'String'))/1000;      %bw
Settings(2) = eval(get(Handle(2), 'String'))/1000;  	  %hw
Settings(7) = eval(get(Handle(3), 'String'))/1000;   	  %gw 
I_one_turn = eval(get(Handle(8), 'String'));		        %I
I_total = eval(get(Handle(6), 'String')) * I_one_turn;  %I_total = NI
Settings(9) = I_total;  									     %I_total
Settings(10) = eval(get(Handle(4), 'String'));          %f                                                                                                                                                                                                                                                                                                                                                                                                                                                  
Settings(11) = eval(get(Handle(5), 'String'))/1000;     %d                                                                                                                                                                                                                                                                                                                                                                                                                                                  
Settings(12) = eval(get(Handle(6), 'String'));          %N                                                                                                                                                                                                                                                                                                                                                                                                                                                
Settings(13) = eval(get(Handle(7), 'String'));          %Fp
Settings(19) = eval(get(Handle(12), 'String'))/1000;          %bobbin breadth
Settings(20) = eval(get(Handle(13), 'String'))/1000;          %bobbin height

if Settings(18)
Settings(16) = eval(get(Handle(9), 'String'))/(2*1000);     %r

else
Settings(17) = eval(get(Handle(11),'String'))/1000;     %turn length
end

if mod(Settings(12), 1) > 0
    disp('Error, number of turns must be an integer.')
    return
end

if Settings(19) > Settings(1)
    disp('Error, bobbin breadth can not be larger than window breadth.');
    return
end

if Settings(20) > Settings(2)
    disp('Error, bobbin height can not be larger than window height.');
    return
end

save(filename, 'Settings', '-APPEND');

% if filling the bobbin
if control == 7
    AA0;
    AA1;
    AA2;
    displayoutput;
    return
end

% if selecting the number of strands
if control == 8
   answer = inputdlg('Enter the number of strands', 'Enter the number of strands');

   
   if isempty(answer)
       return
   end
       
     n = eval(answer{1});
      
   if mod(n, 1) >0 | n <= 0
       disp('Error, number of strands must be a positive integer.');
       return
   end
   windingheight = n*Settings(12)*Settings(11)^2/(Settings(13)*Settings(19));
   if windingheight > Settings(20)
       disp('Error, too many strands.  Available bobbin space exceeded.');
       return
   end
   
   Settings(22) = windingheight;
   save(filename, 'Settings', '-APPEND');
   
   AA0;
   AA1;
   AA2;
   displayoutput;
   return
end

if control ==9
    awg;
    return
end
   
if control >3	    %close shapeopt
    filename =0;
   return;
end

%Settings(15) = 0;
Settings(15) = 1;

corerad = Settings(16);					%set up the parameters to pass to Newloss
wheight = Settings(1);
userchoice = Settings(18);
wspacing = Settings(5);

if  (control == 0 )
%runs this if certain settings have been changed


save(filename, 'Settings', '-APPEND');
AA0;           
AA1;                                                                                                                                                                                                                                                                                                                                                                                                                                                                
AA2;                                    %the actual COOS algorithim                                                                                                                                                                                                                                                                                                                                                                                                                            
AA4mod;                                                                                                                                                                                                                                                                                                                                                                                                                                                            

elseif  (control == 1)
%runs this if certain settings have been changed
                                                                                                                                                                                                                                                                                                                                                                  
save(filename, 'Settings', '-APPEND');
AA1;                                                                                                                                                                                                                                                                                                                                                                                                                                                                
AA2;                                                                                                                                                                                                                                                                                                                                                                                                                                                                
AA4mod; 
                                                                                                                                                                                                                                                                                                                                                                                                                                                             
else                                                                                                                                                                                                                                                                                                                                                                                                                                                                
                                                                                                                                                                                                                                                                                                                                                                                                                                
save(filename, 'Settings', '-APPEND');
AA2;                                                                                                                                                                                                                                                                                                                                                                                                                                                                
AA4mod;      

end                       


function AA0()
%File name: AA0.m
%Revision: rob jensen summer00- sets up settings now instead of having to use A0 also
%	Settings are now saved in one file of which the user specifies the name
%
%Purpose:
%        To calculate the component of magnetic field due to gap current.
%
%Synopsis:
%        No requirements.
%
%Description:
%        Similar to A7, this program will try to compute the B field due to gap current
%        



format compact

filename = get(findobj('tag','mainwindow'), 'userdata');                                                                             
load(filename);  

%set values for the width, height of the winding window
bw = Settings(1);
hw = Settings(2);

bbreadth = Settings(19);
bheight = Settings(20);

% these are longer needed, because they are done in the gui, but
% are left here for completeness

%Settings(3) = 2;                    %unit of x                            
%Settings(4) = 2;                    %unit of y  

%Settings(5) = 10;              	%xdiv                                      
%Settings(6) = 10;             		%ydiv        

%Settings(8) =10;         			%gapdiv                                         
%Settings(14) =.1;             		%step                                      
 

%select an arbitrary value of images of unit on one side
%the fundamental winding window and its first order image on the left form a unit
unit_of_X = Settings(3);
unit_of_Y = Settings(4); %totally there will be (2*unit_of_X+1)*(2*unit_of_Y+1)=9 units

%set values for division of bw and hw
x_div = Settings(5);
y_div = Settings(6);

gw = Settings(7);
gap_div = Settings(8);

%set value for total current
I_total = Settings(9);

%other parameters
f = Settings(10);
d = Settings(11);
N = Settings(12);
Fp = Settings(13);

Xpmult = Settings(21);
windheight = Settings(22);

%--------------------------------------------------------------------
%compute x_grid,y_grid
x_grid=bheight/x_div;
y_grid=bbreadth/y_div;

%create X,Y vector for quiver plot
%X=linspace(bheight/x_div,bheight,x_div)-1/2*x_grid;
%Y=linspace(bbreadth/y_div,bbreadth,y_div)-(y_div+1)/2*y_grid;
%save(filename, 'X', 'Y', 'Settings', '-APPEND');


%R0_x_component=linspace(bw/x_div,bw,x_div)-1/2*x_grid;
%R0_y_component=linspace(hw/y_div,hw,y_div)-(y_div+1)/2*y_grid;
%R0 is the matrix in which every cell contains the coordinates of the center
%of a grid in fundamental winding window
R0_x_component=linspace(bheight/x_div,bheight,x_div)-1/2*x_grid+(hw-bheight);
R0_y_component=linspace(bbreadth/y_div,bbreadth,y_div)-(y_div+1)/2*y_grid;
R0=R0_x_component'*ones(1,y_div)+j*ones(x_div,1)*R0_y_component;
save(filename, 'R0','-APPEND');
clear R0_x_component R0_y_component;

%--------------------------------------------------------------------
%initial setting of the shape of the winding is a rectangle
%pattern_right=ones(y_div,x_div/2);
%pattern_left=zeros(y_div,x_div/2);
%pattern_new = [pattern_left, pattern_right];

if doption(1) == 2
    %initial setting of the shape of the winding is a rectangle
    pattern_right=zeros(y_div, floor(x_div*(1-(windheight/bheight))));
    pattern_left=ones(y_div, x_div-1-floor(x_div*(1-(windheight/bheight))));
    leftover = 1-mod((x_div*(1-(windheight/bheight))), 1);
    leftovercolumn = ones(y_div, 1)*leftover;

    pattern_new = [pattern_left, leftovercolumn, pattern_right];

else
    pattern_new = ones(x_div, y_div);
    
end
    
%start with a new pattern based on experience and the results we got
%pattern_new(6:15,1:5) = 0;
%pattern_new(7:14,6) = 0;
%pattern_new(8:13,7) = 0;
%pattern_new(9:12,8) = 0;
%disp(pattern_new);
%save Pattern pattern_new;
%save SETTING pattern_new -APPEND;
save(filename, 'pattern_new', '-APPEND');
%--------------------------------------------------------------------

u0=4*pi*1e-07;  %permeability of free space

omega_0=2*pi*f; %angular frequency
rohc=1/(5.7e7); %resistivity of conductor (copper)

pattern_matrix=ones(x_div,y_div);
%pattern_matrix(:,1:12)=0;    %set up initial pattern to match the simulation setup in t3
%pattern_matrix(:,17:20)=0;      %I don't really know why so I'm commenting this out

Bgap=zeros(x_div,y_div);
BgapT=zeros(x_div,y_div);
I_per_gap_div=I_total/gap_div;

%parameters for subroutine
para(1)=unit_of_X;
para(2)=unit_of_Y;
para(3)=bw;
para(4)=hw;
para(5)=x_div;
para(6)=y_div;
para(7)=gw;
para(8)=gap_div;
para(9)=I_per_gap_div;
para(10)=omega_0;
para(11)=rohc;
para(12)=d;
para(13)=N;
para(14)=Fp;

const2=pi*omega_0^2/(128*rohc);
para(15)=const2;

%parameters for function Bfinite1.m
pvec=zeros(1,10);
pvec(1:9)=para(1:9);
clear para;

gap_grid=gw/gap_div;

%create a vector of gap current positions of the base unit
%create a vector of gap current positions of the base unit

if get(findobj('tag', 'centergaps'), 'value')
     Rgbase=linspace(-(gw-gap_grid)/2,(gw-gap_grid)/2,gap_div);
     RgbaseT = Rgbase;
     condition = 0;
elseif get(findobj('tag', 'outergaps'), 'value')
     Rgbase=linspace(-(gw-gap_grid)/2,(gw-gap_grid)/2,gap_div);
     RgbaseT = linspace(-(gw-gap_grid)/2,(gw-gap_grid)/2,gap_div);
     condition = 1;
elseif get(findobj('tag', 'allgapsround'), 'value')
     Rgbase1=linspace(-(gw-gap_grid)/2,(gw-gap_grid)/2,gap_div);
     Rgbase2=linspace(-(gw-gap_grid)/2,(gw-gap_grid)/2,gap_div);
     Rgbase2T=linspace(-(gw-gap_grid)/2,(gw-gap_grid)/2,gap_div);
     Rgbase = [Rgbase1, Rgbase2];
     RgbaseT = [Rgbase1, Rgbase2T];

     k = 1.23;  % see Andy Hoke paper on E-core transformers
     l = ((bw/2)+pi*Settings(16)/(2*sqrt(2)));
     center_perim = 2*pi*Settings(16);
     outer_perim = center_perim/sqrt(2);
     
     Center_reluc_gap = gw/(u0*pi*Settings(16)^2);
     Center_reluc_faces = pi/(center_perim*u0*(1+log(pi*l/(2*gw))));
     Center_reluc_corners = 1/(u0*k*l);
     Center_reluc = 1/((1/Center_reluc_corners)+(1/Center_reluc_faces)+(1/Center_reluc_gap));
     Outer_reluc_corners = 1/(u0*k*l);
     Outer_reluc_faces =  pi/(outer_perim*u0*(1+log(pi*l/(2*gw))));
     Outer_reluc_gap = 2*Center_reluc_gap;
     % ignore the reluctance of the faces for now - they're probably pretty close
     Outer_par_comb = .5/((1/Outer_reluc_corners)+(1/Outer_reluc_faces)+(1/Outer_reluc_gap));
    
     I_per_gap_div_outer = (I_total/gap_div)*Outer_par_comb/(Outer_par_comb+Center_reluc);
     I_per_gap_div_inner = (I_total/gap_div)*Center_reluc/(Outer_par_comb+Center_reluc);

elseif get(findobj('tag', 'allgapssquare'), 'value')
     Rgbase1=linspace(-(gw-gap_grid)/2,(gw-gap_grid)/2,gap_div);
     Rgbase2=linspace(-(gw-gap_grid)/2,(gw-gap_grid)/2,gap_div);
     Rgbase2T=linspace(-(gw-gap_grid)/2,(gw-gap_grid)/2,gap_div);
     Rgbase = [Rgbase1, Rgbase2];
     RgbaseT = [Rgbase1, Rgbase2T];
    
     I_per_gap_div = I_total/(gap_div*2);
     k = 1.23;  % see Andy Hoke paper on E-core transformers
     l = ((bw/2)+pi*Settings(16)/(2*sqrt(2)));
     center_perim = 2*pi*Settings(16);
     outer_perim = center_perim/sqrt(2);
     
     Center_reluc_gap = gw/(u0*(pi*Settings(16)/2)^2);
     Center_reluc_faces = pi/(center_perim*u0*(1+log(pi*l/(2*gw))));
     Center_reluc_corners = 1/(u0*k*l);
     Center_reluc = 1/((1/Center_reluc_corners)+(1/Center_reluc_faces)+(1/Center_reluc_gap));
     Outer_reluc_corners = 1/(u0*k*l);
     Outer_reluc_faces =  pi/(outer_perim*u0*(1+log(pi*l/(2*gw))));
     Outer_reluc_gap = 2*Center_reluc_gap;
     % ignore the reluctance of the faces for now - they're probably pretty close
     Outer_par_comb = .5/((1/Outer_reluc_corners)+(1/Outer_reluc_faces)+(1/Outer_reluc_gap));
     
     I_per_gap_div_outer = (I_total/gap_div)*Outer_par_comb/(Outer_par_comb+Center_reluc);
     I_per_gap_div_inner = (I_total/gap_div)*Center_reluc/(Outer_par_comb+Center_reluc);
else
end
 
Rgbase=Rgbase'*j;
RgbaseT =RgbaseT'*j;

for k=1:length(Rgbase)
    
   if k>length(Rgbase)/2 & (get(findobj('tag', 'allgapsround'), 'value') | get(findobj('tag', 'allgapssquare'), 'value'))
       pvec(9) = I_per_gap_div_outer;
       condition = 1;
   elseif get(findobj('tag', 'allgapsround'), 'value') |  get(findobj('tag', 'allgapssquare'), 'value')
        pvec(9) = I_per_gap_div_inner;
        condition = 0;
   end
    
   pvec(10)=Rgbase(k);
   Btemp=Bfinite(pvec, R0, condition);
   Bgap=Bgap+Btemp;         %parallel path to core
   
   pvec(3) = bw*Xpmult;
   pvec(4) = hw*Xpmult;
   pvec(10) = RgbaseT(k);
   BtempT = Bfinite(pvec, R0, condition);  %extended path
   BgapT = BgapT+BtempT;
   
    pvec(3) = bw;
    pvec(4) = hw;
end


save(filename, 'Bgap','BgapT','-APPEND');

% This figure is not necessary, but could be included if desired
% Fields due to the gap 

%   lat = linspace(0, bheight*1000, x_div);
%	vert = linspace(0, bbreadth*1000, y_div);  %vectors for plots

%figure(6)
%subplot(2,2,3)
%v = linspace(max(max(abs(Bgap))), min(min(abs(Bgap))), 20);
%contourf(lat,vert,abs(Bgap),v);
%title('Field due to gap - parallel section');

%subplot(2,2,4)
%v = linspace(max(max(abs(BgapT))), min(min(abs(BgapT))), 20);
%contourf(lat,vert,abs(BgapT),v);
%title('Perpendicular section')

%subplot(2,2,1)
%quiver(lat,vert,real(Bgap),imag(Bgap));
%title('Field due to gap - parallel section')

%subplot(2,2,2)
%quiver(lat,vert,real(BgapT),imag(BgapT));
%title('Perpendicular section')


function AA1()
%file name: AA1.m        Oct. 30, 1997
%Revision:
%        July, 2002, Rob Jensen - added waitbar, fixed the mixup of bw and hw, etc
%
%Purpose:
%        Try to make intelligent use of the multidimensional array feature in Matlab 5 and use
%        it to store the data of the magnetic field due to current elements. 
%
%Synopsis:
%        No special requirements. Parameters have to be changed in the gui.
%
%Description:
%        Create a multidimensional arry to store the data of magnetic field. Temporary we set
%        it to be three dimensional. Each page is a (ydiv,xdiv) matrix. E.g., page0 of this
%        huge array is the magnetic field due to current filament at position (0,-1/2*hw). 
%        Current in every element is assumed to be unit 1. Therefore in later programs using 
%        this array, the value of each page has to be scaled by a factor of 
%        Itotal/sum(sum(pattern))*pattern(y,x)
%
%Conclusion:
%        It yields the same result as we got before. But this time all the data of magnetic 
%        field is save in a three dimensional matrix. Later computation can simply look up 
%        this huge matrix instead of calling subroutine Bfinite.m and Bfinite1.m again and 
%        again. 



format long g;

u0=4*pi*1e-07;

filename = get(findobj('tag','mainwindow'), 'userdata');                                                                            
load(filename);
%set arbitrary values for the width, height of the winding window
bw = Settings(1);
hw = Settings(2);

bbreadth = Settings(19);
bheight = Settings(20);

%select an arbitrary value of images of unit on one side
%the fundamental winding window and its first order image on the left form a unit
unit_of_X = Settings(3);
unit_of_Y = Settings(4); %totally there will be (2*unit_of_X+1)*(2*unit_of_Y+1)=9 units

%set arbitrary values for division of bw and hw
x_div = Settings(5);
y_div = Settings(6);

gw = Settings(7);
gap_div = Settings(8);

%set arbitrary value for total current
I_total = Settings(9);

f = Settings(10);
d = Settings(11);
N = Settings(12);
Fp = Settings(13);

Xpmult = Settings(21);

%compute x_grid,y_grid
x_grid=bw/x_div;
y_grid=hw/y_div;

omega_0=2*pi*f;
rohc=1/(5.7e7);

%--------------------------------------------------------------------
%create X,Y vector for quiver plot
%X=linspace(bw/_div,bw,x_div)-1/2*x_grid;
%Y=linspace(hw/y_div,hw,y_div)-(y_div+1)/2*y_grid;
%save(filename, 'X', 'Y', '-APPEND');

%R0 is the matrix in which every cell contains the coordinates of the center
%of a grid in fundamental winding window
%R0_x_component=linspace(bw/x_div,bw,x_div)-1/2*x_grid;
%R0_y_component=linspace(hw/y_div,hw,y_div)-(y_div+1)/2*y_grid;
%R0=R0_x_component'*ones(1,y_div)+j*ones(x_div,1)*R0_y_component;
%save R0 R0;
%save SETTING R0 -APPEND;
%save(filename, 'R0', '-APPEND');

%clear R0_x_component R0_y_component

%set up initial conditions
%pattern_matrix=ones(y_div,x_div);
%pattern_matrix=ones(x_div,y_div);

Bdata=zeros(x_div,y_div);
BdataT=zeros(x_div,y_div);

%I_block=I_total/sum(sum(pattern_matrix));
I_block=1;    %That's why the data has to be scaled. 

%parameters for subroutine
para(1)=unit_of_X;
para(2)=unit_of_Y;
para(3)=bw;
para(4)=hw;
para(5)=x_div;
para(6)=y_div;
para(7)=gw;
para(8)=gap_div;
para(9)=I_block;
para(10)=omega_0;
para(11)=rohc;
para(12)=d;
para(13)=N;
para(14)=Fp;

%parameters for function Bfinite1.m
pvec=zeros(1,10);
pvec(1:9)=para(1:9);
clear para;
page=1;

row=0;  % just to be user friendly
h = waitbar(row/y_div, 'Computing Magnetic Field');

for xindex=1:x_div
   for yindex=1:y_div
      R1=R0(xindex,yindex);
      pvec(10)=R1;
      Btemp=Bfinite1(pvec,R0);
      Bdata(:,:,page)=Btemp;
      
      pvec(3) = bw*Xpmult;
      pvec(4) = hw*Xpmult;
      BtempT = Bfinite1(pvec, R0);
      BdataT(:,:,page) = BtempT;
      
      pvec(3) = bw;
      pvec(4) = hw;
      
      page=page+1;
      
   end

row = row + 1;
   waitbar(row/y_div);
   pause(.01)
end

close(h);    %get rid of the waitbar

save(filename, 'Bdata', 'BdataT','-APPEND');


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

%figure
%quiver(X,Y,real(Bconductor),imag(Bconductor));
%title('B field due to all current filaments computed by using superposition and Bfinite.m');
%save bconductor Bconductor;
%save SETTING Bconductor -APPEND;
%save(filename, 'Bconductor','-APPEND');


function AA2()
%File name: AA2.m
%Revision: 
%			Rob Jenensen July 2001, fixed Irms in a calculation.
%
%Purpose:
%			Prepare data for finding the optimum shape of the winding by numerical approach.
%
%Synopsis:
%			No command line arguments. All the source data are in BDATA.mat.
%
%Description:
%			After generating the data of the magnetic field due to a unit 
%			current inside every small block, we start to optimize the 
%			shape of the winding again using numerical approach.
%        The optimization process is not implemented here. In order to make debugging  
%        easier, it's implemented in a seperate file AA3.m.


%A0.m, AA0.m, AA1.m must have been executed before AA2.m runs. 
%disp('AA2 starting')
format long g;


 filename = get(findobj('tag','mainwindow'), 'userdata');                                                                                 
load(filename);
%set arbitrary values for the width, height of the winding window

hw = Settings(1);
bw = Settings(2);

bbreadth = Settings(19);
bheight = Settings(20);

%select an arbitrary value of images of unit on one side
%the fundamental winding window and its first order image on the left form a unit
unit_of_X = Settings(3);
unit_of_Y = Settings(4); %totally there will be (2*unit_of_X+1)*(2*unit_of_Y+1)=9 units

%set arbitrary values for division of bw and hw
x_div = Settings(5);
y_div = Settings(6);

gw = Settings(7);
gap_div = Settings(8);

%set arbitrary value for total current
I_total = Settings(9);

f = Settings(10);
d = Settings(11);
N = Settings(12);
Fp = Settings(13);

corerad = Settings(16);					%set up the parameters to pass to Newloss
wheight = Settings(1);
userchoice = Settings(18);
wspacing = Settings(5);

%compute x_grid,y_grid
%x_grid=bw/x_div;
%y_grid=hw/y_div;

omega_0=2*pi*f;
rohc=1/(5.7e7);
u0=4*pi*1e-07;

%The initial shape of the winding is no more the one defined in A0.m.
%load Pattern;

scale_factor=I_total/sum(sum(pattern_new));

c0=4*rohc*I_total^2*x_div*y_div/(pi*Fp*bbreadth*bheight);
c1=pi*omega_0^2*bbreadth*bheight*Fp*d^2/(64*rohc*x_div*y_div);      %adjusted to take Irms
%c1=pi*omega_0^2*bw*hw*Fp*d^2/(128*rohc*x_div*y_div);
%The above line appears to be wrong - corrected above

LossPara=zeros(1,2);
LossPara(1)=c0;
LossPara(2)=c1;

Xpmult = Settings(21);

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


BC=zeros(x_div,y_div);
BCT=zeros(x_div,y_div);

for xindex=1:x_div
   for yindex=1:y_div
      
      if pattern_new(yindex,xindex)~=0
         pagenum=(xindex-1)*y_div+yindex;
         Btemp=Bdata(:,:,pagenum);
         BC=BC+pattern_new(yindex,xindex).*scale_factor.*Btemp;
         
         BtempT=BdataT(:,:,pagenum);
         BCT=BCT+pattern_new(yindex,xindex).*scale_factor.*BtempT;
         
      end
      
   end
end


Btotal=BC+Bgap;
Bcond = BC;

BtotalT = BCT+BgapT;
BcondT = BCT;

Ppresent = pattern_new;

save(filename, 'BC', 'Bcond', 'Btotal','Ppresent','LossPara','BCT','BcondT','BtotalT','-APPEND');

%load Xvec;
%load Yvec;
%figure(1)
%quiver(X,Y,real(BC),imag(BC));
%title('B field due to all current filaments computed by using superposition and Bfinite.m');
%figure(2)
%quiver(X,Y,real(Bgap),imag(Bgap));
%title('B field due to the air gap');

%[totalloss,acloss,dcloss]=Newloss(LossPara,pattern_new,Btotal, corerad, wheight, userchoice, wspacing)

%P=pattern_new;
%P(5,1)=0;
%P(6,1)=0;
%pagenum1=(1-1)*10+5;
%Btemp1=scale_factor.*Bdata(:,:,pagenum1);
%pagenum2=(1-1)*10+6;
%Btemp2=scale_factor.*Bdata(:,:,pagenum2);
%scale_new=(sum(sum(P))+2)/sum(sum(P));
%Btotal_new=scale_new.*(BC-Btemp1-Btemp2)+Bgap;

%[totalloss_new,acloss_new,dcloss_new]=Newloss(LossPara,P,Btotal_new, corerad, wheight, userchoice, wspacing)
%disp('AA2 ending')






function AA4mod()
%Dan Gestwick and Reid Hutchins modification to Jiankin Hu's code
%This code is stable and is approximately 9 times faster.
%File name: AA4mod.m
%Revision: 1.0
%
%Purpose: 
%        Find the optimum shape of the winding by numerical approach. 
%
%Synopsis:
%        Several matrices are basic necessary data.
%        BC: the B field due to conductor in winding window. It's generated by AA2.m.
%        Bgap: the B field due to air gap in winding window. It's generated by AA1.m.
%        Bdata: a 3-d matrix, where every page is the B field in the winding window due
%               to a unity current in a specific position.
%        Common settings: These settings, like x_div,y_div, they have to be the same as
%        in the other files.
%        Pattern: the present pattern of the winding. Generated by AA2.m.
%
%Purpose:
%        The algorithm here is based on AA3.m. 
%        The extension is: now we start to add 'step' back to elements. 
%			The first stage is applying a very simple and straight forward 
%			algorithm. In this algorithm, the initial condition is that every
%			element of matrix pattern is 1. The algorithm only removes strands
%			and never add them back. 
%
%        First step of optimization: try to remove every element respectively. Save the
%        value of total loss in a matrix, in which an element like LossM(yindex,xindex)
%        stores the value of total loss if Pattern_new(yindex,xindex) is removed. 
%        Suppose that LossM(y',x') is the maximum in this matrix. This element will be 
%        removed in this round. This maximum will also be saved in a matrix, maybe a 
%        vector, and the position (x',y') will also be saved for future tracing.
%

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

% flops(0);
format long g;
%format short;
u0=4*pi*1e-07;

addStep = 1;  %add a 'step' to the current element
removeStep = 0; %remove a 'step' from the current element

%load SETTING;
filename = get(findobj('tag','mainwindow'), 'userdata');                                                                              
load(filename);
%set arbitrary values for the width, height of the winding window
hw = Settings(1);
bw = Settings(2);

bbreadth = Settings(19);
bheight = Settings(20);

%select an arbitrary value of images of unit on one side
%the fundamental winding window and its first order image on the left form a unit
unit_of_X = Settings(3);
unit_of_Y = Settings(4); %totally there will be (2*unit_of_X+1)*(2*unit_of_Y+1)=9 units

%set arbitrary values for division of bw and hw
x_div = Settings(5);
y_div = Settings(6);

gw = Settings(7);
gap_div = Settings(8);

%set arbitrary value for total current
I_total = Settings(9);
%Overwrite default value, might not be a good practice, but can meet curiosity more quickly
%I_total = 1;

f = Settings(10);
d = Settings(11);
N = Settings(12);
Fp = Settings(13);

Xpmult = Settings(21);

corerad = Settings(16);					%set up the parameters to pass to Newloss
wheight = Settings(1);
userchoice = Settings(18);
wspacing = Settings(5);

whichprogram = Settings(15);

%compute x_grid,y_grid
%x_grid=bw/x_div;
%y_grid=hw/y_div;

omega_0=2*pi*f;
rohc=1/(5.7e7);
%--------------------------------------------------------------------
%f = 1e6;

%--------------------------------------------------------------------
%load the necessary data
%load Pattern;

xempty = Settings(6) - round(Settings(6)*bbreadth/bw);
yempty = Settings(5) - round(Settings(5)*bheight/(hw*2));

Ppresent=pattern_new;
clear pattern_new;

%trigger=Settings(14);
trigger=.1;
step=1;  

LossPara=zeros(1,2);
c0=4*rohc*I_total^2*x_div*y_div/(pi*Fp*bbreadth*bheight);
%c1=pi*omega_0^2*bw*hw*Fp*d^2/(128*rohc*x_div*y_div);  This appears to be wrong I correct it on the next line
c1=pi*omega_0^2*bbreadth*bheight*Fp*d^2/(64*rohc*x_div*y_div);	 %Scaled to take Irms
LossPara(1)=c0;
LossPara(2)=c1;
[totalloss,acloss,dcloss]=Newloss(LossPara,Ppresent,Btotal, Settings);
%the totaloss computed here is usually the highest
highestLoss=totalloss;
TotalLowestRecord = [highestLoss];  %keep track of every detail
ACLowestRecord=acloss;
DCLowestRecord=dcloss;
coordinate = [0 0];
%There's actually no coordiate (0,0) for centers of blocks 
AddRemove = [0]; %create a vector tracking the add or remove information

%a flag used to control the flow
%To to remove some conductor at least once
improvement = 1;  
TLossAdd=zeros(y_div/2,x_div);
ALossAdd=zeros(y_div/2,x_div);
DLossAdd=zeros(y_div/2,x_div);
TLossRemove=zeros(y_div/2,x_div);
ALossRemove=zeros(y_div/2,x_div);
DLossRemove=zeros(y_div/2,x_div);

working = 0;   %initialize counters
count = 0;     

if  doption(4) 
    figh = figure;  %animation of opt
end

while step>=trigger,
      
   scale_current=step*I_total/sum(sum(Ppresent));
   %Every page of Bdata stores the B field due to a current equals to '1' in a specific block
   
   scale_remove=sum(sum(Ppresent))/(sum(sum(Ppresent))-2*step);
   scale_add=sum(sum(Ppresent))/(sum(sum(Ppresent))+2*step);
   %after removing or adding 2*step conductor, the total current should be restored to make 
   %it the same
   
   for x1=1:x_div
      for y1=1:1:y_div/2  %the other half will be (x1,10-y1+1)
         %Ideally there're three options
         %1) P = 0: Try adding only.
         %2) P = step -- 1-step: Try both adding and removing.
         %3) P = 1: Try removing only.
         %But due to the rounding problem
         %1) -0.01*step <= P < 0.99*step: Try adding only.
         %2) 0.99*step <= P < 1-0.99*step: Try both adding and removing.
         %3) 1-0.99*step <= P: Try removing only.
         
         Ptemp=Ppresent;
         
         if Ptemp(y1,x1) >= -0.01*step & Ptemp(y1,x1) < 0.99*step %Try adding only
            
            Padd=Ptemp;
            Padd(y1,x1)=Ptemp(y1,x1)+step;
            Padd(y_div-y1+1,x1)=Ptemp(y_div-y1+1,x1)+step;
            
            pagenum1=(x1-1)*y_div+y1;
            pagenum2=(x1-1)*y_div+(y_div-y1+1);
            
            Bunit_add=scale_current.*(Bdata(:,:,pagenum1)+Bdata(:,:,pagenum2));
            Bcond_add=scale_add.*(Bunit_add+Bcond);
            Btotal_add=Bcond_add+Bgap;
            
            Bunit_addT=scale_current.*(BdataT(:,:,pagenum1)+BdataT(:,:,pagenum2));
            Bcond_addT=scale_add.*(Bunit_addT+BcondT);
            Btotal_addT=Bcond_addT+BgapT;
            
            [tladd,aladd,dladd]=Newloss(LossPara,Padd,Btotal_add, Settings);
            %if 'step' is added to Ppresent(y1,x1), the total loss, dc loss, 
            %ac loss will be as following
            
            [tladdT,aladdT,dladdT]=Newloss(LossPara,Padd,Btotal_addT, Settings);
            
            % here, the losses must be scaled by the weighting factor
            % between the parallel and permendicular sections
            
            if whichprogram & andyswitch
                w = weight(corerad, Padd, hw, bheight);
            else
                w=0;
            end
            
            tladd = w*tladdT+(1-w)*tladd;
            aladd = w*aladdT+(1-w)*aladd;
            dladd= w*dladdT+(1-w)*dladd;
            
            TLossAdd(y1,x1)=tladd;  
            ALossAdd(y1,x1)=aladd;
            DLossAdd(y1,x1)=dladd;
           
            TLossRemove(y1,x1)=highestLoss;  %saturate the removing loss matrix
            ALossRemove(y1,x1)=highestLoss;
            DLossRemove(y1,x1)=highestLoss;
            
            
         elseif Ptemp(y1,x1) >= 0.99*step & Ptemp(y1,x1) <= (1-0.99*step)
            %Try both adding and removing
            Padd=Ptemp;
            Premove=Ptemp;
            
            Premove(y1,x1)=Ptemp(y1,x1)-step;
            Premove(y_div-y1+1,x1)=Ptemp(y_div-y1+1,x1)-step;
            
            Padd(y1,x1)=Ptemp(y1,x1)+step;
            Padd(y_div-y1+1,x1)=Ptemp(y_div-y1+1,x1)+step;
            
            pagenum1=(x1-1)*y_div+y1;
            pagenum2=(x1-1)*y_div+(y_div-y1+1);
            
            Bunit_add=scale_current.*(Bdata(:,:,pagenum1)+Bdata(:,:,pagenum2));
            Bunit_remove=-Bunit_add;
            
            Bcond_add=scale_add.*(Bunit_add+Bcond);
            Btotal_add=Bcond_add+Bgap;
            
            Bunit_addT=scale_current.*(BdataT(:,:,pagenum1)+BdataT(:,:,pagenum2));
            Bunit_removeT=-Bunit_addT;
            
            Bcond_addT=scale_add.*(Bunit_addT+BcondT);
            Btotal_addT=Bcond_addT+BgapT;
         
            
            [tladd,aladd,dladd]=Newloss(LossPara,Padd,Btotal_add, Settings);
            
            [tladdT,aladdT,dladdT]=Newloss(LossPara,Padd,Btotal_addT, Settings);
            
            % here, the losses must be scaled by the weighting factor
            % between the parallel and permendicular sections
            
            if whichprogram & andyswitch
                w = weight(corerad, Padd, hw, bheight);
            else
                w=0;
            end
                
                
            tladd = w*tladdT+(1-w)*tladd;
            aladd = w*aladdT+(1-w)*aladd;
            dladd = w*dladdT+(1-w)*dladd;
            
            TLossAdd(y1,x1)=tladd;  
            ALossAdd(y1,x1)=aladd;
            DLossAdd(y1,x1)=dladd;
  
            
            Bcond_remove=scale_remove.*(Bunit_remove+Bcond);
            Btotal_remove=Bcond_remove+Bgap;
            
            Bcond_removeT=scale_remove.*(Bunit_removeT+BcondT);
            Btotal_removeT=Bcond_removeT+BgapT;
            
            [tlremove,alremove,dlremove]=Newloss(LossPara,Premove,Btotal_remove, Settings);
            
            [tlremoveT,alremoveT,dlremoveT]=Newloss(LossPara,Premove,Btotal_removeT, Settings);
            
            % here, the losses must be scaled by the weighting factor
            % between the parallel and permendicular sections
            
            if whichprogram & andyswitch
                w = weight(corerad, Premove, hw, bheight);
            else 
                w=0;
            end
            
            tlremove = w*tlremoveT+(1-w)*tlremove;
            alremove = w*alremoveT+(1-w)*alremove;
            dlremove = w*dlremoveT+(1-w)*dlremove;
            
            TLossRemove(y1,x1)=tlremove;
            ALossRemove(y1,x1)=alremove;
            DLossRemove(y1,x1)=dlremove;
            
            
         elseif Ptemp(y1,x1) > (1-0.99*step) & Ptemp(y1,x1) < (1+0.01*step) %Try removing only 
            Premove=Ptemp;
            
            Premove(y1,x1)=Ptemp(y1,x1)-step;
            Premove(y_div-y1+1,x1)=Ptemp(y_div-y1+1,x1)-step;
            
            pagenum1=(x1-1)*y_div+y1;
            pagenum2=(x1-1)*y_div+(y_div-y1+1);
            
            Bunit_remove=-scale_current.*(Bdata(:,:,pagenum1)+Bdata(:,:,pagenum2));
            
            Bcond_remove=scale_remove.*(Bunit_remove+Bcond);
            Btotal_remove=Bcond_remove+Bgap;
            
            Bunit_removeT=-scale_current.*(BdataT(:,:,pagenum1)+BdataT(:,:,pagenum2));
            
            Bcond_removeT=scale_remove.*(Bunit_removeT+BcondT);
            Btotal_removeT=Bcond_removeT+BgapT;
            
            [tlremove,alremove,dlremove]=Newloss(LossPara,Premove,Btotal_remove, Settings);
            
            [tlremoveT,alremoveT,dlremoveT]=Newloss(LossPara,Premove,Btotal_removeT, Settings);
            
            % here, the losses must be scaled by the weighting factor
            % between the parallel and permendicular sections
            
            if whichprogram & andyswitch
                w = weight(corerad, Premove, hw, bheight);
            else
                w=0;
            end
                
                
            tlremove = w*tlremoveT+(1-w)*tlremove;
            alremove = w*alremoveT+(1-w)*alremove;
            dlremove = w*dlremoveT+(1-w)*dlremove;
            
            TLossRemove(y1,x1)=tlremove;
            ALossRemove(y1,x1)=alremove;
            DLossRemove(y1,x1)=dlremove;

            TLossAdd(y1,x1)=highestLoss;  %saturate the adding loss matrix
            ALossAdd(y1,x1)=highestLoss;
            DLossAdd(y1,x1)=highestLoss;
            
         else
            disp('Problem has ocured. Check the booundary settings.');
            disp(Ppresent);
            disp(x1);
            disp(y1);
            pause;
         end
         
         %It will save execution time to implement adding 'step' to the current element 
         %within this double 'for' loop
         %Bunit doesn't have to be calculated again. 'adding' only changes the sign of 'Bunit'
         
      end
   end
   
   %Let's compare the Loss after adding or removing a 'step'
   AddMin=minimum(TLossAdd);
   RemoveMin=minimum(TLossRemove);
   
   %totalloss is the loss before this optimization step
   %improvement=1, this optimization is reasonable
   %improvement=-1 or 0, give up this optimization step
 
   if AddMin(3) <= RemoveMin(3) & AddMin(3) < totalloss  %Adding is beneficial
      improvement = 1; %Acknowledge this adding
      
      yindex=AddMin(1);
      xindex=AddMin(2);
      
      Ppresent(yindex,xindex)=Ppresent(yindex,xindex)+step;
      Ppresent(y_div-yindex+1,xindex)=Ppresent(y_div-yindex+1,xindex)+step;
      
      coordinate=[coordinate;xindex yindex];
      AddRemove=[AddRemove;addStep];

      pagenum1=(xindex-1)*y_div+yindex;
      pagenum2=(xindex-1)*y_div+(y_div-yindex+1);
      Bunit_add=scale_current.*(Bdata(:,:,pagenum1)+Bdata(:,:,pagenum2));
      Bcond=scale_add.*(Bcond+Bunit_add);

      Bunit_addT=scale_current.*(BdataT(:,:,pagenum1)+BdataT(:,:,pagenum2));
      BcondT=scale_add.*(BcondT+Bunit_addT);
      
      
      [totalTestAdd,acTestAdd,dcTestAdd]=Newloss(LossPara,Ppresent,Bcond+Bgap, Settings);
      
      [totalTestAddT,acTestAddT,dcTestAddT]=Newloss(LossPara,Ppresent,BcondT+BgapT, Settings);
      
      if whichprogram & andyswitch
        w = weight(corerad, Ppresent, hw, bheight);
      else
        w=0;
      end
        
      totalloss= w*totalTestAddT + (1-w)*totalTestAdd;
      
      totalTestAdd = w*totalTestAddT+(1-w)*totalTestAdd;
      acTestAdd = w*acTestAddT+(1-w)*acTestAdd;
      dcTestAdd = w*dcTestAddT+(1-w)*dcTestAdd;
      
      TotalLowestRecord=[TotalLowestRecord;AddMin(3)];
      ACLowestRecord=[ACLowestRecord;acTestAdd];
      DCLowestRecord=[DCLowestRecord;dcTestAdd];
      
      if abs(totalTestAdd-AddMin(3)) > 1e-2
         disp('Error has occured during this adding.');
         pause;  %usually this should not happen
      end
      
      
   elseif RemoveMin(3) < AddMin(3) & RemoveMin(3) < totalloss  %Removing is better
      improvement = 1; %Acknowledge this adding
      
      yindex=RemoveMin(1);
      xindex=RemoveMin(2);
      
      Ppresent(yindex,xindex)=Ppresent(yindex,xindex)-step;
      Ppresent(y_div-yindex+1,xindex)=Ppresent(y_div-yindex+1,xindex)-step;
      
      coordinate=[coordinate;xindex yindex];
      AddRemove=[AddRemove;removeStep];

      pagenum1=(xindex-1)*y_div+yindex;
      pagenum2=(xindex-1)*y_div+(y_div-yindex+1);
      Bunit_remove=-scale_current.*(Bdata(:,:,pagenum1)+Bdata(:,:,pagenum2));
      Bcond=scale_remove.*(Bcond+Bunit_remove);
      
  
      
      Bunit_removeT=-scale_current.*(BdataT(:,:,pagenum1)+BdataT(:,:,pagenum2));
      BcondT=scale_remove.*(BcondT+Bunit_removeT);
      
                  % here, the losses must be scaled by the weighting factor
            % between the parallel and permendicular sections

      [totalTestRemove,acTestRemove,dcTestRemove]=Newloss(LossPara,Ppresent,Bcond+Bgap, Settings);
      
      [totalTestRemoveT,acTestRemoveT,dcTestRemoveT]=Newloss(LossPara,Ppresent,BcondT+BgapT, Settings);
      
      if whichprogram & andyswitch
        w = weight(corerad, Ppresent, hw, bheight);
      else
        w = 0;
      end
      
      totalloss=w*totalTestRemoveT+(1-w)*totalTestRemove;
      
      totalTestRemove = w*totalTestRemoveT+(1-w)*totalTestRemove;
      acTestRemove = w*acTestRemoveT+(1-w)*acTestRemove;
      dcTestRemove = w*dcTestRemoveT+(1-w)*dcTestRemove;
      
      TotalLowestRecord=[TotalLowestRecord;RemoveMin(3)];
      ACLowestRecord=[ACLowestRecord;acTestRemove];
      DCLowestRecord=[DCLowestRecord;dcTestRemove];
      
      if abs(totalTestRemove-RemoveMin(3)) > 1e-2
         disp('Error has occured during this removing.');
      end

   else
      step = step/3;
            
   end
   
   %disp('Ppresent is:');
   %disp(Ppresent);
   
   %count = count+1;   	%to show the algorithim is still going
   %if count - 4 == 1
    %  working = working+1
     % count = 0;
     %end

    lat = linspace(0, bheight*1000, x_div);
	vert = linspace(0, bbreadth*1000, y_div);  %vectors for plots
  
  if  doption(4)  
			figh;
			contourf(lat, vert, Ppresent, [0.5 0.5]);
			axis equal;
			axis([0 bheight*1000 0 bbreadth*1000]);
			title('Optimal shape of the winding (wire placement in green)  - gap at (0, bw/2)');
			xlabel('bobbin height (mm)');
         ylabel('bobbin breadth (mm)');
   
	end
pause(.01)
working = working+1;
end %end of while loop    
 
% figh;
%close
     
 %[totalloss,acloss,dcloss]=Newloss(LossPara,Ppresent,Btotal, corerad, wheight, userchoice, wspacing);
%[totalloss,acloss,dcloss]=Newloss(LossPara,Ppresent,Btotal, Settings)
%Btotal hasn't been updated.

[totalloss,acloss,dcloss]=Newloss(LossPara,Ppresent,Bcond+Bgap, Settings);
Btotal = Bcond+Bgap;
[totallossT,aclossT,dclossT]=Newloss(LossPara,Ppresent,BcondT+BgapT, Settings);
BtotalT = BcondT+BgapT;
 save(filename, 'AddRemove', 'coordinate', 'TotalLowestRecord',...
   'ACLowestRecord', 'DCLowestRecord', 'Ppresent', 'Bcond', 'Bgap', 'Btotal', 'BcondT', 'BgapT', 'BtotalT','-APPEND');

displayoutput;


function  displayoutput()
% Displays the plots and final calculations



% flops(0);
format long g;
%format short;
u0=4*pi*1e-07;

filename = get(findobj('tag','mainwindow'), 'userdata');                                                                               
load(filename);
%set arbitrary values for the width, height of the winding window
hw = Settings(1);
bw = Settings(2);

bbreadth = Settings(19);
bheight = Settings(20);

%select an arbitrary value of images of unit on one side
%the fundamental winding window and its first order image on the left form a unit
unit_of_X = Settings(3);
unit_of_Y = Settings(4); %totally there will be (2*unit_of_X+1)*(2*unit_of_Y+1)=9 units

%set arbitrary values for division of bw and hw
x_div = Settings(5);
y_div = Settings(6);

gw = Settings(7);
gap_div = Settings(8);

%set arbitrary value for total current
I_total = Settings(9);

f = Settings(10);
d = Settings(11);
N = Settings(12);
Fp = Settings(13);

corerad = Settings(16);					%set up the parameters to pass to Newloss
wheight = Settings(1);
userchoice = Settings(18);
wspacing = Settings(5);

whichprogram = Settings(15);

      if whichprogram & andyswitch
        w = weight(corerad, Ppresent, hw, bheight);
      else
        w = 0;
      end


%compute x_grid,y_grid
%x_grid=bw/x_div;
%y_grid=hw/y_div;

omega_0=2*pi*f;
rohc=1/(5.7e7);

                                                                               
load(filename);


   lat = linspace(0, bheight*1000, x_div);
	vert = linspace(0, bbreadth*1000, y_div);  %vectors for plots
  
  if  doption(4)
			figure
	        contourf(lat, vert, Ppresent, [0.5 0.5]); 
            if min(min(Ppresent)) == 1
                fill([0, bheight*1000, bheight*1000, 0],[0, 0, bbreadth*1000, bbreadth*1000],'g');   %to fix the case where the bobbin is fulll
            end
			axis equal;
			axis([0 bheight*1000 0 bbreadth*1000]);
			title('Optimal shape of the winding (wire placement in green)  - gap at (0, bw/2)');
			xlabel('bobbin height (mm)');
         ylabel('bobbin breadth (mm)');
   
	end
pause(.01)

 %end of while loop

 if doption(2) | doption(3)
     figh2 = figure;
 end

if doption(2)
      figh2;

      if doption(3)
         subplot(2,2,1); 
     else
         subplot(1,2,1);
     end
			%quiver(lat,vert,real(Btotal),imag(Btotal));
			quiver(lat,vert,real(Bcond+Bgap),imag(Bcond+Bgap));
			axis equal;
			axis([0 bheight*1000 0 bbreadth*1000]);
			title('Total B field - gap at (0, bw/2), Plane parallel to core');
			xlabel('bobbin height (mm)');
         ylabel('bobbin breadth (mm)');
         
     if doption(3)
         subplot(2,2,2); 
     else
         subplot(1,2,2);
     end
			%quiver(lat,vert,real(Btotal),imag(Btotal));
			quiver(lat,vert,real(BcondT+BgapT),imag(BcondT+BgapT));
			axis equal;
			axis([0 bheight*1000 0 bbreadth*1000]);
			title('Plane perpendicular to core');
			xlabel('bobbin height (mm)');
         ylabel('bobbin breadth (mm)');
         
         
         
         
end
      

if  doption(3)
		figh2;

	v = linspace(max(max(abs(Bcond+Bgap))), min(min(abs(Bcond+Bgap))), 20);

    if doption(2)
        subplot(2,2,3);
    else
        subplot(1,2,1);
    end
	%contourf(lat, vert, abs(Btotal), v);
	contourf(lat, vert, abs(Bcond+Bgap), v);
		axis equal;
		axis([0 bheight*1000 0 bbreadth(1)*1000]);
		title('Magnitude of B  - gap at (0, bw/2), plane parallel to core');
		xlabel('bobbin height (mm)');
		ylabel('bobbin breadth (mm)');
    if doption(2)
        subplot(2,2,4);
    else
        subplot(1,2,2);
    end
   v  = linspace(max(max(abs(BcondT+BgapT))), min(min(abs(BcondT+BgapT))), 20);
     	contourf(lat, vert, abs(BcondT+BgapT), v);
		axis equal;
		axis([0 bheight*1000 0 bbreadth(1)*1000]);
		title('Plane perpendicular to core');
		xlabel('bobbin height (mm)');
		ylabel('bobbin breadth (mm)');  
        
        
        
end

if doption(3) | doption(4)
    set(gcf, 'units','points')
    set(gcf, 'position',[0 0 450 450])
end


   %disp(Ppresent);
%save up the results for future analysis

corerad = Settings(16);					%set up the parameters to pass to Newloss
wheight = Settings(1);
userchoice = Settings(18);
wspacing = Settings(5);
Pattern = Ppresent;




x_div = Settings(5);
y_div = Settings(6);

%The following lines display the results of the optimization in the command window

[totalloss,acloss,dcloss]=Newloss(LossPara,Ppresent,Bcond+Bgap, Settings);
[totallossT,aclossT,dclossT]=Newloss(LossPara,Ppresent,BcondT+BgapT, Settings);

totalloss = w*totallossT+(1-w)*totalloss
acloss = w*aclossT+(1-w)*acloss
dcloss = w*dclossT+(1-w)*dcloss

tot_area = Settings(19)*Settings(20)*1e6  %to put it into mm^2
A_used_fraction = sum(sum(Ppresent))/x_div/y_div
area_used = tot_area*A_used_fraction
B_squared_avg = sum(sum(Ppresent.*(abs(Btotal).^2)))/(x_div*y_div)
Rac = totalloss/((Settings(9))/Settings(12))^2
Rdc = dcloss/((Settings(9))/Settings(12))^2
n = (area_used*1e-6 * Settings(13))/ (Settings(12)*(Settings(11)^2))

disp('Now, the number of strands is rounded, and the loss is adjusted accordingly')
suggested_n = round(n)
new_acloss = acloss*(suggested_n/n)
new_dcloss = dcloss*(n/suggested_n)
new_totalloss = new_acloss+new_dcloss

disp('Units are in Watts, Ohms, Tesla, and Milimeters');




final_shape = Ppresent;

B = B_squared_avg;
   

v = [1:wspacing];
onemat = ones(1, Settings(6))';
vtwodim = onemat*v;
R = (((vtwodim -.5)*wheight)./wspacing) + corerad + Settings(2)-Settings(20);

lavg = (sum(sum(Pattern.*R))*2*pi)/(sum(sum(Pattern)));

intBsqR = sum(sum(Pattern.*(abs(B).^2).*R*2*pi));
%This will be in units of power - not power/meter

% the following code may be used to find the winding height for center gapped des

%if find(sum(Pattern)==0
%windheight = (bheight/x_div)*(min(find(sum(Pattern)==0))-1)*1000
%Settings(22) = windheight;
%end
%flops

save(filename, 'tot_area','w','A_used_fraction', 'area_used',...
   'B_squared_avg', 'Rac', 'lavg', 'Rdc', 'intBsqR', 'n', 'LossPara', 'final_shape', 'totalloss', 'acloss', 'dcloss', '-APPEND');




function y=Bfinite(p_vec, R0, condition)
%File name: Bfinite.m
%
%Purpose:
%        Calculate B field due to one element carry current which flows in the direction of 
%        into the plane (cross).
%
%Synopsis:
%        y=Bfinite(pvec)
%        where pvec is a vector of parameters with all the necessary information for this 
%        function, y is a matrix containing the value of B field at the center of every
%        element. Size of y is (ydiv, xdiv).
%
%Revision:
%        None.
%
%new fast vectorized version



%load SETTING;
%                                                                              
%load(filename);
unit_of_X=p_vec(1);
unit_of_Y=p_vec(2);
bw=p_vec(3);
hw=p_vec(4);
x_div=p_vec(5);
y_div=p_vec(6);
gw=p_vec(7);
gap_div=p_vec(8);
I_unit=p_vec(9);
R1_1=p_vec(10);
R1_2=-real(R1_1)+j*imag(R1_1); %get R1's first order image of left side

%bbreadth = Settings(19);
%bheight = Settings(20);

u0=4*pi*1e-7;

x_grid=hw/x_div;
y_grid=bw/y_div;

%Note the that direction of this current filament is opposite to the direction of current
%flowing in the conductors
const1=u0*(j)*I_unit/(2*pi);

%matrix R0 stores all positions within the fundamental winding window
%this matrix should be generated by the main program
%load R0;
R0_vector=reshape(real(R0)'+j*imag(R0)',1,x_div*y_div);

%create a matrix storing the coordinates of the center of every unit
x_temp=(2*hw).*linspace(-unit_of_X,unit_of_X,2*unit_of_X+1);
y_temp=bw.*linspace(-unit_of_Y,unit_of_Y,2*unit_of_Y+1);
center_matrix=ones(size(y_temp))'*x_temp+j.*y_temp'*ones(size(x_temp));
%center_matrix = center_matrix+.5*hw+.5*bw*j;

if condition        % outer gap;
    center_matrix = center_matrix+hw;
end

R1_1_matrix=R1_1+center_matrix;
size_tmp=size(center_matrix);
R1_2_matrix=R1_2+center_matrix; 
R1_matrix=zeros(size_tmp(1),2*size_tmp(2));
R1_matrix(:,1:2:(2*size_tmp(2)-1))=R1_1_matrix;
R1_matrix(:,2:2:2*size_tmp(2))=R1_2_matrix;
R1_vector=reshape(R1_matrix,1,2*size_tmp(1)*size_tmp(2));

Rp=ones(size(R1_vector))'*R0_vector-R1_vector.'*ones(size(R0_vector));
Rp_abs=abs(Rp);
size_of_Rp=size(Rp);

%for xindex=1:size_of_Rp(1)
 %  for yindex=1:size_of_Rp(2)
  %    if Rp_abs(xindex,yindex)==0
   %      Rp_abs(xindex,yindex)=1;
  %    end
 %  end
%end
     % This line does what above loop does--puts in a dummy value of 1 to avoid divide by zero errors.
Rp_abs = ((Rp_abs == 0)+Rp_abs);
     
Bfinite_temp=const1.*Rp./(Rp_abs.^2);
Bfinite=reshape(sum(Bfinite_temp),x_div,y_div);

y=Bfinite;






function y=Bfinite1(p_vec, R0)
%File name: Bfinite1.m
%
%Purpose:
%        Calculate B field due to one element carry current which flows in the direction of 
%        out of the plane (dot).
%
%Synopsis:
%        y=Bfinite1(pvec)
%        where pvec is a vector of parameters with all the necessary information for this 
%        function, y is a matrix containing the value of B field at the center of every
%        element. Size of y is (ydiv, xdiv).
%
%Revision:
%        vectorized stuff near end speeds it up and replaces commented out nested for loop
%			C Sullivan spring 2000
%
%This function calculates B field due to a current filament (current in a grid) 
%at the center of every grid within the fundamental winding window
%direction of current: out of the plane (dot)

%load SETTING;
%                                                                                  
%load(filename);

unit_of_X=p_vec(1);
unit_of_Y=p_vec(2);
bw=p_vec(3);
hw=p_vec(4);
x_div=p_vec(5);
y_div=p_vec(6);
gw=p_vec(7);
gap_div=p_vec(8);
I_unit=p_vec(9);
R1_1=p_vec(10);
R1_2=-real(R1_1)+j*imag(R1_1);  %get R1's first order image of left side



u0=4*pi*1e-7;

x_grid=hw/x_div;
y_grid=bw/y_div;

%Note the that direction of this current filament is opposite to the direction of current
%flowing in the conductors
const1=u0*(-j)*I_unit/(2*pi);

%matrix R0 stores all positions within the fundamental winding window
%this matrix should be generated by the main program
%load R0;
R0_vector=reshape(R0.',1,x_div*y_div); 

%create a matrix storing the coordinates of the center of every unit
x_temp=(2*hw).*linspace(-unit_of_X,unit_of_X,2*unit_of_X+1);
y_temp=bw.*linspace(-unit_of_Y,unit_of_Y,2*unit_of_Y+1);
center_matrix=ones(size(y_temp))'*x_temp+j.*y_temp'*ones(size(x_temp));

R1_1_matrix=R1_1+center_matrix;
size_tmp=size(center_matrix);
R1_2_matrix=R1_2+center_matrix;
R1_matrix=zeros(size_tmp(1),2*size_tmp(2));
R1_matrix(:,1:2:(2*size_tmp(2)-1))=R1_1_matrix;
R1_matrix(:,2:2:2*size_tmp(2))=R1_2_matrix;
R1_vector=reshape(R1_matrix,1,2*size_tmp(1)*size_tmp(2));

Rp=ones(size(R1_vector))'*R0_vector-R1_vector.'*ones(size(R0_vector));
Rp_abs=abs(Rp);
size_of_Rp=size(Rp);

%for xindex=1:size_of_Rp(1)
 %  for yindex=1:size_of_Rp(2)
  %    if Rp_abs(xindex,yindex)==0
   %      Rp_abs(xindex,yindex)=1;
  %    end
 %  end
%end
     % This line does what above loop does--puts in a dummy value of 1 to avoid divide by zero errors.
Rp_abs = ((Rp_abs == 0)+Rp_abs);


Bfinite1_temp=const1.*(Rp./Rp_abs.^2);
B1=reshape(sum(Bfinite1_temp),y_div,x_div);

y=B1;




function [TotalLoss,ACLoss,DCLoss]=Newloss(LossPara,Pattern,B, Settings)
%Revision: 1.0
%
%Purpose:
%        Compute the loss for a specific configuration.
%
%
%Description:
%        Refer to document 'Thesis12.doc' for definitions of the parameters.
%        Pdc/meter = c0*1/sum(sum(Pattern))
%        Ppe/meter = c1*sum(sum(Pattern.*abs(B)^2));
%        See comments below for notes on what units the if-else outputs
%        outputs will be in units of power (not power/length)


corerad = Settings(16);					%set up the parameters to pass to Newloss
wheight = Settings(2);
userchoice = Settings(18);
wspacing = Settings(5);
bheight = Settings(20);
hw = Settings(2);

if ~userchoice
   
   %loss = loss/unit length* Settings(17) %correction factor for length of wire
   DCLoss=LossPara(1).*(1/sum(sum(Pattern))) * Settings(17); 
   ACLoss=LossPara(2).*sum(sum(Pattern.*(abs(B).^2))) * Settings(17);
   

else
   
%wspacing = Settings(5);
   
v = [1:wspacing];
onemat = ones(1, Settings(6))';
vtwodim = onemat*v;
R = (((vtwodim -.5)*bheight)./wspacing) + corerad + hw-bheight;

lavg = (sum(sum(Pattern.*R))*2*pi)/(sum(sum(Pattern)));

DCLoss=LossPara(1).*(1/sum(sum(Pattern)))*lavg;
ACLoss=LossPara(2).*sum(sum(Pattern.*(abs(B).^2).*R*2*pi));
%This will be in units of power - not power/mmeter

end

TotalLoss=DCLoss+ACLoss;

function y=minimum(M)
%File name: MINIMUM.m
%
%Purpose:
%        Find the index and value of the minimum value of a matrix.
%
%Synopsis:
%        y=minimum(Matrix)
%        where y(1) is the row index of the minimum value, 
%        y(2) is the column index of the minimum value,
%        y(3) is the value of the minimum
%
%Description:
%        None.
%
%Revision:
%        None.
%
[row_vector,row_number]=min(M);
[column_value,column_number]=min(row_vector);
y=[row_number(column_number),column_number(1),column_value];


function showdimensions()
% This function displays a figure which shows the definition of the dimensions

subplot(2,2,4)
axis square
axis off
axis([0 10 0 10])
f = 9;

% the winding window
line([2,2],[2,4.8],'linewidth',3)
line([2,2],[5.2,8],'linewidth', 3)
line([2,8],[8,8],'linewidth', 3)
line([8,8],[8,2],'linewidth', 3)
line([8,2],[2,2],'linewidth', 3)

text(3, 7.5, 'Window Breadth','fontsize',f, 'color', 'r')
%text(0, 6.5, 'Breadth','fontsize',f, 'color', 'r')
text(3,1.5,'Window Height','fontsize',f, 'color', 'r')
text(.5,5,'Gap','fontsize',f)

%arrows
a = 2;
line([0.35+a, 0.35+a], [2,8], 'color', 'r')
line([.35+a, .5+a], [8,7.5], 'color', 'r')
line([.35+a, .2+a], [8,7.5], 'color', 'r')

line([2,8],[1,1], 'color', 'r')
line([8, 7.5],[1,1.5], 'color', 'r')
line([8,7.5], [1,.5], 'color', 'r')

%the bobbin
line([3,8],[3,3],'linewidth', 3, 'color', 'k')
line([8,8],[3,7],'linewidth', 3, 'color', 'k')
line([3,8],[7,7],'linewidth', 3, 'color', 'k')
line([3,3],[3,7],'linewidth', 3, 'color', 'k')

text(3.5,6, 'Bobbin Breadth','fontsize',f)
text(3.5,4,'Bobbin Height','fontsize',f)

% bobbin arrows
line([3.35, 3.35], [3,7], 'color', 'k')
line([3.35, 3.5], [7,6.5], 'color', 'k')
line([3.35, 3.2], [7,6.5], 'color', 'k')

line([3,8],[3.5,3.5], 'color', 'k')
line([8, 7.5],[3.5,3.75], 'color', 'k')
line([8,7.5], [3.5,3.25], 'color', 'k')



function fig = desoptions()
% This is the machine-generated representation of a Handle Graphics object
% and its children.  Note that handle values may change when these objects
% are re-created. This may cause problems with any callbacks written to
% depend on the value of the handle at the time the object was saved.
% This problem is solved by saving the output as a FIG-file.
%
% To reopen this object, just type the name of the M-file at the MATLAB
% prompt. The M-file and its associated MAT-file must be on your path.
% 
% NOTE: certain newer features in MATLAB may not have been saved in this
% M-file due to limitations of this format, which has been superseded by
% FIG-files.  Figures which have been annotated using the plot editor tools
% are incompatible with the M-file/MAT-file format, and should be saved as
% FIG-files.

% Rob Jensen July 2002



% colormap for figure (machine generated)
mat0(1, :) = [		    0         0    0.5625];
mat0(2, :) = [        0         0    0.6250];
mat0(3, :) = [        0         0    0.6875];
mat0(4, :) = [        0         0    0.7500];
mat0(5, :) = [        0         0    0.8125];
mat0(6, :) = [        0         0    0.8750];
mat0(7, :) = [        0         0    0.9375];
mat0(8, :) = [        0         0    1.0000];
mat0(9, :) = [        0    0.0625    1.0000];
mat0(10, :) = [        0    0.1250    1.0000];
mat0(11, :) = [        0    0.1875    1.0000];
mat0(12, :) = [        0    0.2500    1.0000];
mat0(13, :) = [        0    0.3125    1.0000];
mat0(14, :) = [        0    0.3750    1.0000];
mat0(15, :) = [        0    0.4375    1.0000];
mat0(16, :) = [        0    0.5000    1.0000];
mat0(17, :) = [        0    0.5625    1.0000];
mat0(18, :) = [        0    0.6250    1.0000];
mat0(19, :) = [        0    0.6875    1.0000];
mat0(20, :) = [        0    0.7500    1.0000];
mat0(21, :) = [        0    0.8125    1.0000];
mat0(22, :) = [        0    0.8750    1.0000];
mat0(23, :) = [        0    0.9375    1.0000];
mat0(24, :) = [        0    1.0000    1.0000];
mat0(25, :) = [   0.0625    1.0000    1.0000];
mat0(26, :) = [   0.1250    1.0000    0.9375];
mat0(27, :) = [   0.1875    1.0000    0.8750];
mat0(28, :) = [   0.2500    1.0000    0.8125];
mat0(29, :) = [   0.3125    1.0000    0.7500];
mat0(30, :) = [   0.3750    1.0000    0.6875];
mat0(31, :) = [   0.4375    1.0000    0.6250];
mat0(32, :) = [   0.5000    1.0000    0.5625];
mat0(33, :) = [   0.5625    1.0000    0.5000];
mat0(34, :) = [   0.6250    1.0000    0.4375];
mat0(35, :) = [   0.6875    1.0000    0.3750];
mat0(36, :) = [   0.7500    1.0000    0.3125];
mat0(37, :) = [   0.8125    1.0000    0.2500];
mat0(38, :) = [   0.8750    1.0000    0.1875];
mat0(39, :) = [   0.9375    1.0000    0.1250];
mat0(40, :) = [   1.0000    1.0000    0.0625];
mat0(41, :) = [   1.0000    1.0000         0];
mat0(42, :) = [   1.0000    0.9375         0];
mat0(43, :) = [   1.0000    0.8750         0];
mat0(44, :) = [   1.0000    0.8125         0];
mat0(45, :) = [   1.0000    0.7500         0];
mat0(46, :) = [   1.0000    0.6875         0];
mat0(47, :) = [   1.0000    0.6250         0];
mat0(48, :) = [   1.0000    0.5625         0];
mat0(49, :) = [   1.0000    0.5000         0];
mat0(50, :) = [   1.0000    0.4375         0];
mat0(51, :) = [   1.0000    0.3750         0];
mat0(52, :) = [   1.0000    0.3125         0];
mat0(53, :) = [   1.0000    0.2500         0];
mat0(54, :) = [   1.0000    0.1875         0];
mat0(55, :) = [   1.0000    0.1250         0];
mat0(56, :) = [   1.0000    0.0625         0];
mat0(57, :) = [   1.0000         0         0];
mat0(58, :) = [   0.9375         0         0];
mat0(59, :) = [   0.8750         0         0];
mat0(60, :) = [   0.8125         0         0];
mat0(61, :) = [   0.7500         0         0];
mat0(62, :) = [   0.6875         0         0];
mat0(63, :) = [   0.6250         0         0];
mat0(64, :) = [   0.5625         0         0];

filename = get(findobj('tag','mainwindow'), 'userdata');
load(filename);


% the figure itself
h_fig = figure('Units','points', ...
	'Color',[0.8 0.8 0.8], ...
	'Colormap',mat0, ...
   'CreateFcn','', ...
   'DeleteFcn', '',...
	'Position',[250 100 200 250], ...
	'Tag','designoptions', ...
	'ToolBar','none', ...
    'NumberTitle','off',...
   'UserData',doption);

h_radio(1) = uicontrol(h_fig,'style','radio',...
    'BackgroundColor',[0.8 0.8 0.8], ...
    'createfcn','filename = get(findobj(''tag'', ''mainwindow''), ''userdata''); load(filename); , if doption(1) ==1, set(findobj(''tag'', ''h_radio(1)''), ''value'',1); end,',...
    'callback','set(findobj(''tag'', ''h_radio(1)''), ''value'', 1); set(findobj(''tag'', ''h_radio(2)''), ''value'', 0); set(findobj(''tag'', ''h_radio(3)''), ''value'', 0); filename = get(findobj(''tag'', ''mainwindow''),''userdata''); doption(1) = 1; save(filename, ''doption'',  ''-APPEND'')',...
    'string', 'Full Bobbin',...
    'tag', 'h_radio(1)',...
    'position', [50 180 200 20]);
    
h_radio(2) = uicontrol(h_fig,'style','radio',...
    'BackgroundColor',[0.8 0.8 0.8], ...
    'createfcn','filename = get(findobj(''tag'', ''mainwindow''), ''userdata''); load(filename); if doption(1) ==2, set(findobj(''tag'', ''h_radio(2)''), ''value'',1); end,',...
    'callback','set(findobj(''tag'', ''h_radio(1)''), ''value'', 0); set(findobj(''tag'', ''h_radio(2)''), ''value'', 1); set(findobj(''tag'', ''h_radio(3)''), ''value'', 0); filename = get(findobj(''tag'', ''mainwindow''),''userdata''); doption(1) = 2; save(filename, ''doption'',  ''-APPEND'')',...
    'string', 'Fixed Strand Number',...
    'tag', 'h_radio(2)',...
    'position', [50 205 200 20]);
    
h_radio(3) = uicontrol(h_fig,'style','radio',...
    'BackgroundColor',[0.8 0.8 0.8], ...
     'createfcn','filename = get(findobj(''tag'', ''mainwindow''), ''userdata'');  load(filename); if doption(1) ==3, set(findobj(''tag'', ''h_radio(3)''), ''value'',1); end,',...
    'callback','set(findobj(''tag'', ''h_radio(1)''), ''value'', 0); set(findobj(''tag'', ''h_radio(2)''), ''value'', 0); set(findobj(''tag'', ''h_radio(3)''), ''value'', 1); filename = get(findobj(''tag'', ''mainwindow''),''userdata''); doption(1) = 3; save(filename, ''doption'',  ''-APPEND'')',...
    'string', 'Shape Optimization ',...
    'fontweight','bold',...
    'tag', 'h_radio(3)',...
    'position', [50 230 200 20]);

left = 50;
%  checkboxes
quiver_box = uicontrol('Parent', h_fig,...
   'BackgroundColor',[0.8 0.8 0.8], ...
   'Position', [left 125 250 15],...
    'createfcn','filename = get(findobj(''tag'', ''mainwindow''), ''userdata''); load(filename); if doption(2), set(findobj(''tag'', ''quiver''), ''value'',1); end,',...
   'callback','filename = get(findobj(''tag'', ''mainwindow''),''userdata'');  load(filename);  doption(2) = ~doption(2); save(filename, ''doption'',  ''-APPEND'');',...
   'String', 'Plot Magnetic Field',...
   'Style', 'checkbox', ...
   'tag', 'quiver');

contour_box = uicontrol('Parent', h_fig,...
   'BackgroundColor',[0.8 0.8 0.8], ...
   'Position', [left 100 250 15],...
    'createfcn','filename = get(findobj(''tag'', ''mainwindow''), ''userdata''); load(filename); if doption(3), set(findobj(''tag'', ''contour''), ''value'',1); end,',...
   'callback','filename = get(findobj(''tag'', ''mainwindow''),''userdata''); load(filename); doption(3) = ~doption(3); save(filename, ''doption'',  ''-APPEND'');',...
   'String', 'Plot Magnetic Field Magnitude Contours',...
   'Style', 'checkbox', ...
   'tag', 'contour');   

shape_box = uicontrol('Parent', h_fig,...
   'BackgroundColor',[0.8 0.8 0.8], ...
   'Position', [left 75 250 15],...
       'createfcn','filename = get(findobj(''tag'', ''mainwindow''), ''userdata''); load(filename); if doption(4), set(findobj(''tag'', ''shape''), ''value'',1); end,',...
   'callback','filename = get(findobj(''tag'', ''mainwindow''),''userdata''); load(filename);  doption(4) = ~doption(4); save(filename,  ''doption'',  ''-APPEND'');',...
   'String', 'Show Shape',...
   'Style', 'checkbox', ...
   'tag', 'shape');

% pushbuttons

%okmat = ['filename = get(findobj(''tag'', mainwindow),''userdata'');'...
%         'if get(findobj(''tag'', h_radio(1)), ''value''),'...
%         'doption(1) = 1;'...
%         'elseif get(findobj(''tag'', h_radio(2)), ''value''),'...
%         'doption(1) = 2;'...
%         'else,'...
%         'doption(1) = 3;'...
%         'end,'...
%         'if get(findobj(''tag'', quiver), ''value''),'...
%         'doption(2) = 1;'...
%         'end,'...
%'if get(findobj(''tag'', contour), ''value''),'...
%         'doption(3) = 1;'...
%'end,'...
%         'if get(findobj(''tag'', shape), ''value''),'...
%         'doption(4) = 1;'...
%         'end,'...
%         'save(doption,  ''-APPEND'');'...
%         'close;'];

h1 = uicontrol('Parent',h_fig, ...
	'BackgroundColor',[.75 .75 .75], ...
	'Callback','close', ...
	'ListboxTop',0, ...
	'Position',[120 40 50 20], ...
	'String','Ok', ...
    'fontweight', 'bold',...
	'Tag','Pushbutton1');
%static text
h1 = uicontrol('Parent',h_fig, ...
	'BackgroundColor',[.8 .8 .8], ...
	'ListboxTop',0, ...
	'Position',[20 250 200 30.75], ...
	'String','Choose a Design Option', ...
	'Style','text', ...
    'fontweight','bold',...
	'Tag','s');

h1 = uicontrol('Parent',h_fig, ...
	'BackgroundColor',[.8 .8 .8], ...
	'ListboxTop',0, ...
    'fontsize',12,...
	'Position',[5 290 280 30.75], ...
	'String','Options', ...
	'Style','text', ...
    'fontweight','bold',...
	'Tag','s');

h1 = uicontrol('Parent',h_fig, ...
	'BackgroundColor',[.8 .8 .8], ...
	'ListboxTop',0, ...
	'Position',[20 140 200 30.75], ...
	'String','Choose Display Options', ...
	'Style','text', ...
    'fontweight','bold',...
	'Tag','s');

if nargout > 0, fig = h0; end

function w = weight(r, shape, hw, bheight)
% computes the weight that the extended path will be given in the averageing
% this is basically the fraction of wire subject to the ex path versus the parallel path

dimsum = sum(shape);
radii = r+hw-bheight + linspace(0+.5*bheight/length(dimsum), bheight-.5*bheight/length(dimsum), length(dimsum));
angles = 2*(pi-2*acos(r./radii));
arclengths = angles.*radii;
area2 = sum(arclengths.*dimsum);
totarea = sum(dimsum*2*pi.*radii);
w = 1- (area2/totarea);

function awg 

answer = inputdlg('Enter AWG #', 'Enter AWG #');

   
   if isempty(answer)
       return
   end
       
    size = eval(answer{1});
      
   if mod(size, 1) >0 | size < 6 | size > 54 
       disp('Error, AWG # must be a positive integer between 6 and 54.');
       return
   end

   switch size
    case 6
        dia  =4.104;
    case 7
        dia  =3.657;
    case 8
        dia  =3.260;
    case 9
        dia  =2.905;
    case 10
        dia  =2.589;
    case 11
        dia  =2.308;
    case 12
        dia  =2.057;
    case 13
        dia  =1.833;
    case 14
        dia  =1.634;
    case 15
        dia  =1.456;
    case 16
        dia  =1.298;
    case 17
        dia  =1.156;
    case 18
        dia  =1.024;
    case 19
        dia  =0.910;
    case 20
        dia  =0.812;
    case 21
        dia  =0.724;
    case 22
        dia  =0.642;
    case 23
        dia  =0.574;
    case 24
        dia  =0.510;
    case 25
        dia  =0.454;
    case 26
        dia  =0.404;
    case 27
        dia  =0.360;
    case 28
        dia  =0.320;
    case 29
        dia  =0.287;
    case 30
        dia  =0.254;
    case 31
        dia  =0.226;
    case 32
        dia  =0.203;
    case 33
        dia  =0.180;
    case 34
        dia  =0.160;
    case 35
        dia  =0.142;
    case 36
        dia  =0.127;
    case 37
        dia  =0.114;
    case 38
        dia  =0.102;
    case 39
        dia  =0.089;
    case 40
        dia  =0.078;
    case 41
        dia  =0.071;
    case 42
        dia  =0.063;
    case 43
        dia  =0.056;
    case 44
        dia  =0.050;
    case 46
        dia = .040;
    case 48
        dia = .031;
    case 50
        dia  =0.025;
    case 52
        dia = .020;
    case 54
        dia = .015;
    otherwise
        dia = 0;
    end
    
    set(findobj('tag', 'diameter'), 'string', dia);