% -------- BEGIN GKLMovie.m ---------------
% Plots a movie of the evolution of electrons in their Larmor orbits
% as a function of axial position.  Efficiency evolution is also 
% plotted.  Since the nonlinear electron evolution equation,
%
%  dP/d{\zeta} = i({\Delta} - 1 + |P|^2)P + i F exp{-2({\zeta}/{\mu})^2 }
%
% is solved, a Gaussian electric field profle is assumed, so the equations 
% are true for a gyrotron or the output cavity of a gyroklystron.
%   The plot shows the electric field with its temporal dependance removed, 
% hence the motion of the electrons shows the relative phase between the 
% electrons in the E-field - effectively, the coordinate system is rotating
% at omega*t.  Hence electrons that appear to move in reverse (-theta 
% direction) after the interaction are actually falling behind in phase 
% according to n*phi = theta + omega*t, where n is the harmonic index, phi 
% is the fast time scale phase of the electron, omega is the frequency of 
% the applied RF field, and t is time.
%   This script opens three figures (1,2,3) and overlays figs 2 and 3 into 
% fig 1.  The time it takes to convert to a movie format and also the final 
% file size depends on the size of figure 1.  Generally, the resulting 
% movie is about 25 seconds long and 4 MB.  MATLAB cannot write to the 
% specified file if it is already open in a movie viewer.
%
% NOTE: requires GKLceta2.m, WaitbarRandColor.m
%
% Last edit by Colin Joye, 10/25/04M
%

%------------------- BEGIN User inputs --------------------

%-- No prebunching (gyrotron: max eff 70%) -----
delta = 0.538;   % detuning
q     = 0.0;    % bunching parameter; q=0 is uniform initial distribution.
psi   = 0.84*pi; % relative phase
F     = 0.145;   % normalized field amplitude
mu    = 15.5;    % normalized length

%-- Optimal parameters (max eff 91%) -----
delta = 0.538;   % detuning
q     = 3.17;    % bunching parameter; q=0 is uniform initial distribution.
psi   = 0.84*pi; % relative phase
F     = 0.145;   % normalized field amplitude
mu    = 15.5;    % normalized length


% Number of initial theta_c values to generate (also the number of 
% particles): O(n^0.5)
tsteps = 16;   

% number of axial zeta steps to take: O(n^0.3)
zsteps = 25;   

makemovie = 0; % makes a movie of the evolution (takes a while)
%Mname     = 'C:\MATLAB6p1\work\research\GKLcontour2\CRM_uniform.avi';
Mname     = 'C:\MATLAB6p1\work\research\GKLcontour2\CRM_optimal.avi';

% Set input and output zeta conditions (=3 is from Tran's paper).
% Values >3 show overbunching of a too-long structure.
BunchingFactor = 3;3.32; 2.87;
%------------------- END User inputs --------------------

h  = waitbar(0,['GKLMovie.m : Starting']);
WaitbarRandColor(h);

u = linspace(0,2*pi*(1-1/200),200); % points for unit circle.
 
% ------- BEGIN: Solve the electron momentum eqn ---------------
% define input vector (initial conditions to Diff Eq):
dqp = [delta, q, psi];

% define phases:
theta_c = linspace(0,2*pi*(1-1/tsteps),tsteps);

zeta_a = -sqrt(BunchingFactor)*mu/2; % starting point for solving DE eqn.
zeta_b = -zeta_a;          % ending point for solving DE eqn.
% define zeta vector
zeta = linspace(zeta_a,zeta_b,zsteps);

% solve pendulum equations and return result:
[ceta,P] = GKLceta2(dqp, zeta, theta_c, F, mu);

% calulate efficiency percentage (function of length):
eff = 100*(1 - ceta);
% ------- END: Solve the electron momentum eqn ---------------

% Now step through and make a movie of the evolution:
p = P;
for k=1:zsteps, 
  P = [1;0]*p(k,:) + [0;1]*5*ones(size(p(1,:))); 
  
% -- Plot electron evolution --
figure(2)
clf(2)
set(2,'Color',[1 1 1]);
hold on;

% Put Gaussian profiled electric field arrows in:
Enorm = exp(-(2*zeta(1:k)/mu).^2);
arx   = 0.2*Enorm(end)*[-2, 1, 1,2,1,1,-2];
ary   = 0.02*[-1,-1,-3,0,3,1, 1];
%arxs  = ones(5,1)*arx;
%arys  = [-1:0.5:1]'*ones(1,length(ary))+ ones(5,1)*ary;
%tcolor(1,1:5,1:3) = ones(5,1)*[0.7 0.95 0.85];
arxs  = ones(3,1)*arx;
arys  = [-1:1]'*ones(1,length(ary))+ ones(3,1)*ary;
tcolor(1,1:3,1:3) = ones(3,1)*[0.7 0.95 0.85];
patch(arxs'-1,arys',tcolor,'EdgeColor','none');
%patch(arxs'+1,arys',tcolor,'EdgeColor','none');
patch(arxs(1:2,:)'+1,arys(1:2,:)',tcolor(:,1:2,:),'EdgeColor','none');
text(-1,1.1,'E-field strength','HorizontalAlignment','center','Color',[0.5 0.75 0.75],'FontSize',12);

% Draw electrons 
plot(real(P),imag(P),'.','MarkerSize',18)
plot(sin(u),cos(u),'k.','MarkerSize',3); %unit circle
plot([-1.5 1.5],[0 0],'k:',[0 0],[-1.5 1.5],'k:');
hold off;

axis([-1.5 1.5 -1.5 1.5]);
axis square
%xlabel('Real Part of Norm. Momentum')
%ylabel('Imag Part of Norm. Momentum')
xlabel('x-direction [in Larmor radii]')
ylabel('y-direction [in Larmor radii]')
%title('Uniform Initial Distribution');
%title('Optimal Prebunching Value');

% label plot with initial parameters:
s =           [  'F_0 = ' num2str(F)    ] ;
s = strvcat(s,['\mu_0 = ' num2str(mu)   ]);
s = strvcat(s,[  'q_0 = ' num2str(q)    ]);
s = strvcat(s,['\Delta_0 = ' num2str(delta)]);
s = strvcat(s,['\psi_0 = ' num2str(psi/pi,2) '\pi']);
text(-1.45,-1.05,s);


% -- Plot electron efficiency evolution --  
figure(3)
clf(3)
set(3,'Color',[1 1 1]);
plot(zeta(1:k),eff(1:k),'b',zeta(1:k),100*Enorm,'g--')
axis([zeta(1) zeta(end) 0 100]);
title(['\eta = ' num2str(eff(k),3) '%']);
xlabel('length, \zeta');
%ylabel('eff');


% -- Now overlay efficiency (fig 2) onto electron plot (fig 1)
% (plot overlay code from Amit Kesar).
figure(1)
clf(1)
set(1,'Color',[1 1 1]);
ax = zeros(2,1);
for fi = 1:2
  ax(fi)=subplot(2,1,fi);
end
% Now copy contents of each figure over to destination figure
% Modify position of each axes as it is transferred 
figure(2)
h2 = get(gcf,'Children');
newh2 = copyobj(h2,1); % copy figure 2 into 1
delete(ax(2));

figure(3)
h3 = get(gcf,'Children');
newh3 = copyobj(h3,1); % copy figure 3 into 1
% re-position the new pasted subfigure [xpos ypos xsize ysize]
set(newh3(1),'Position',[0.64 0.8 0.24 0.1]); 
delete(ax(1));


% -- Get the movie frames:
figure(1)
MM(:,k) = getframe(1); 

waitbar(k/zsteps, h,['Generating frames: ' num2str(k/zsteps*100,3) '% done' ]);
end % for k
close(h)

% play the movie in figure:
% movie(MM)

% convert the movie matrix to a movie format:
if makemovie==1
h  = waitbar(0.5,['GKLMovie.m : Converting to movie']);
WaitbarRandColor(h);
movie2avi(MM,Mname,'compression','Cinepak','quality',50,'fps',10);
close(h)
end

% -------- END GKLMovie.m ---------------