% GKLcontour2:
% Following T.M. Tran and B.G. Danly in their paper, "Optimization
% of gyroklystron efficiency," Phys. Fluids, 4/1986, pp1274-1281.  
% Equations marked "EQN #" are referred from this paper.
%   The purpose is to generate plots of the optimal efficiency 
% "eta_perp," detuning parameter "delta," bunching parameter "q" 
% and relative phase parameter "psi" as functions of the normalized 
% field strength "F" and the normalized length "mu" (plots on page 1279).
%   This is accomplished by use of the matlab 'fminsearch' command 
% operating on the 'ode45' differential equation solver.
%
% Note: requires the function GKLceta2.m, which calculates the complementary
% effociency using ODE45 and a nonlinear differential equation.  Also 
% requires the function sec2time.m and WaitbarRandColor.m
%   This program may take a long time to run.  For rsteps = 20, tsteps = 150, 
% zsteps = 300, it took ~45 minutes on my 3GHz machine (half-processor).
%   The plots at the end often require a lot of attention to obtain useful
% results.  Often this means putting the contour labels on by hand ('manual' 
% mode) and choosing which contours to show, e.g. 'con=[0.4:0.1:1.2]'.
%
% Last edit by Colin Joye, 8/20/03

clear all;
tic;

% ------------ USER INPUT ------------------------------------------
% 0 = start value, 1 = end value, mesh controlled by 'rsteps' value.

% optimal efficiency region.
F0     = 0.12;  % normalized field amplitude 
F1     = 0.18;  % normalized field amplitude 
mu0    = 14;    % normalized length 
mu1    = 18;    % normalized length 

F0     = 0.05;  % normalized field amplitude 
F1     = 0.16;  % normalized field amplitude 
mu0    = 3.7;   % normalized length 
mu1    = 4.7;   % normalized length 

F0     = 0.01;  % normalized field amplitude 
F1     = 0.50;  % normalized field amplitude 
mu0    = 0.5;   % normalized length 
mu1    = 18;    % normalized length 

rsteps = 10;    % (>1) number of steps to take for F & mu resolution: O(n^2)
tsteps = 10;   % number initial theta_c values to generate: O(n^0.5)
zsteps = 25;   % number of zeta steps to take: O(n^0.3)

% method of placing numbers on curves = 'manual'/[];  If 'manual' is chosen, 
% you must click the points on the curve for 'clabel'.
mode = []; 

% '=1' saves contour variables as *.mat files, '=0' does not.
keep = 1;       
% ------------ END USER INPUT --------------------------------------

% Initial guess values (guess high)
delta0 = 0.54;    % detuning factor 
q0     = 3.17;    % bunching parameter
psi0   = 0.84*pi;   % relative phase between cavities
dqp0   = [delta0, q0, psi0 ]; % vector containing guesses

% generate F and mu arrays
inc   = (F1 - F0) / rsteps;
F     = [F0:inc:(F1-inc)];
inc   = (mu1 - mu0) / rsteps;
mu    = [mu0:inc:(mu1-inc)];

% declare contour matrices as funcions of F and mu ("fm").
eta_fm   = zeros(rsteps);
q_fm     = eta_fm;
delta_fm = eta_fm;
psi_fm   = eta_fm;

% define theta vector
theta_c = [0:1/(tsteps-1):1]*2*pi;

disp([' ']);
h = waitbar(0,'Computing variations');
WaitbarRandColor(h);

for f = 1:rsteps, % master loop for F values.
for m = 1:rsteps, % master loop for mu values.

% Usual starting and ending points taken for DE solution (max efficiency).
zeta_a = -sqrt(3)*mu(m)/2; % starting point for solving DE eqn.
zeta_b = -zeta_a;          % ending point for solving DE eqn.
% define zeta vector
zinc = (zeta_b - zeta_a) / (zsteps-1);
zeta = [zeta_a: zinc :zeta_b];

% use matlab command to optimize for delta, q and psi. MUCH faster than 
% writing one using for-loops!!  'fminbnd' would be better, since it lets 
% you set constraints on the values, but it can only optimize 1-D functions.
[dqp_min,ceta_min] = fminsearch( @GKLceta2, dqp0, [] , zeta, theta_c, F(f), mu(m) );

% assign minimized values to contour variables
eta_fm(f,m)   = 1 - ceta_min;
delta_fm(f,m) = dqp_min(1);
q_fm(f,m)     = dqp_min(2);
psi_fm(f,m)   = dqp_min(3);

% set initial guesses to optimal values (cuts time down noticeably), but 
% double them.  I think the q/psi go negative if the initial guess is too small.
%dqp0 = 2*abs(dqp_min);

% show waitbar and compute estimated time remaining
pcnt = (f-1)/rsteps + m/rsteps^2 ;
waitbar( pcnt ,h,['Computing: ' num2str(pcnt*100,3) '% done']);
pcnt = rsteps^2 / (rsteps * (f-1) + m) - 1;
time_est = sec2time( ceil(toc*pcnt) );
disp( ['Time Remaining: ' time_est] );
end
end
close(h)

% save data after a long run!
if(keep==1),
save eta eta_fm;
save delta delta_fm;
save q q_fm;
save psi psi_fm;
end

% find optimal parameters
eta_opt = max(max(eta_fm));
[fopt,mopt] = ind2sub( size(eta_fm), find( eta_fm == eta_opt ) );
fopt = fopt(1); mopt = mopt(1); % incase multiple maxes occur.

% record optimal parameters
Fopt      = F(fopt);
muopt     = mu(mopt);
delta_opt = delta_fm(fopt,mopt);
q_opt     = q_fm(fopt,mopt);
psi_opt   = psi_fm(fopt,mopt);

% display actual computation time
disp(['Actual run time: ' sec2time(toc)]);

% -------------- Figures ----------------------
figure(1)
clf(1)
set(1,'color',[1 1 1]); % set figure background color to white for printing
A = mu(1); B = mu(end); C = F(1); D = F(end);
hx = [A B];
hy = Fopt*[1 1];
vx = muopt*[1 1];
vy = [C D];

subplot(221)
con = 0.2:0.1:0.9; % contour values
hold on;
plot(hx,hy,'m:',vx,vy,'m:');
[c,h] = contour(mu, F, eta_fm);%, con);
hold off;
clabel(c,h,mode);
axis([A B C D]);
title(['Max eff: ' num2str(100*eta_opt,3) '%, F_{opt} = ' num2str(Fopt,4) ...
    ', \mu_{opt} = ' num2str(muopt,3)]);
xlabel('\mu');
ylabel('F');

subplot(222)
con = floor(min(min(q_fm))):1:ceil(max(max(q_fm))); % contour values
hold on;
plot(hx,hy,'m:',vx,vy,'m:');
[c,h] = contour(mu, F, q_fm);%, con);
hold off;
clabel(c,h,mode);
axis([A B C D]);
title(['Optimal q: ' num2str(q_opt,3)]);
xlabel('\mu');
ylabel('F');

subplot(223)
%con = floor(min(min(delta_fm))):0.1:ceil(max(max(delta_fm))); % contour values
hold on;
plot(hx,hy,'m:',vx,vy,'m:');
[c,h] = contour(mu, F, delta_fm);%, con);
hold off;
clabel(c,h,mode);
axis([A B C D]);
title(['Optimal \Delta: ' num2str(delta_opt,3)]);
xlabel('\mu');
ylabel('F');

subplot(224)
con = 0.1:0.2:1.1;
hold on;
plot(hx,hy,'m:',vx,vy,'m:');
[c,h] = contour(mu, F, psi_fm/pi);%, con);
hold off;
clabel(c,h,mode);
axis([A B C D]);
title(['Optimal \psi: ' num2str(psi_opt/pi,3) '\pi']);
xlabel('\mu');
ylabel('F');

% --------------- End figures ----------------------


% ---------- END GKLcontour.m -------------------------