clear;
%close all;

aa=cputime;             % Used to compute length of time to run simulation (beginning time)
% Initial State Space Representation for Cavity
A=[-2*pi*150 0; -0 -2*pi*150];
B=[2*pi*720 0; 0 2*pi*720];
C=eye(2);
D=[0 0;0 0];

% ZOOM PLOTS AXIS DEFS
ttmin=17.5;
ttmax=19;
Ibmin=-100; 
Ibmax=100;
BPHmin=-5;
BPHmax=5;
ampmin=2.49;
ampmax=2.51;
phasemin=-1;
phasemax=1;
IgImin=0.1;
IgImax=0.5;
IgQmin=-0.5;
IgQmax=0.5;


% BEGIN SL21 VALUES
f=1.5e9;                % Resonance Frequency of cavity in Hz
w0=2*pi*f;              % Resonance Frequency of cavity in rad/sec
cavity_length = 0.5 ;   % length of cavity in meters
% END SL21 VALUES

%OFFSET=2*pi*1.5e9*tan(-30*pi/180)/(2*6.6e6);
OFFSET=2*pi*0;        % Initial cavity (static) detuning in Hz df=f0*tan(psi)/(2*Ql)

% BEGIN SL21 VALUES
V0=5e6*cavity_length;  % Cavity Voltage = Energy gain in cavity in MV
% END SL21 VALUES

% BEGIN LIA'S DOIT_NEW VALUES
%V0=3e6;
% END LIA'S DOIT_NEW VALUES

PHI0=0;                 % Phase of beam wrt the RF crest in degrees (on-crest = 0)

% BEGIN SL21 VALUES
RL=480;                 % Cavity Shunt Impedance in ohms
QL=1.8e7;                 % Loaded Q (unitless)
wL=w0/(2*QL);           % w0/(2*QL) in rad/sec
% END SL21 VALUES

% BEGIN LIA'S DOIT_NEW VALUES
%RL=960/2;
%QL=6.6e6;
%wL=2*pi*1.5e9/(2*QL);
% END LIA'S DOIT_NEW VALUES

Ib=-65e-6;               % Beam Amplitude in Amperes (Added to generator.  Assume that beam takes
                        %      energy out of the cavity (convention: minus sign to take energy out))
Tb=3e-3;                % Beam period in seconds
dtb=0.5e-3;               % Length of macropulse for pulsed beam in seconds (Used to compute pulse
                        %      width 100*(dtb)/(Tb))
t0=6e-3;                % Time of beam turn on in seconds relative to beginning of simulation
tau1=1e-3;              % Decay Time of Mechanical Mode
Km=60/(2.5e6*2.5e6);    % Lorentz Coefficient (Hz/((MV/m)^2)) Amount of Frequency Shift due to
                        %      Mechanical Mode
tau2=0.2;               % NOT USED
wm=2*pi*65;             % NOT USED Mechanical Frequency in rad/sec
w_mfs=2*pi*50;          % Amplitude of microphonic noise frequency effect in rad/sec
w_mfs_m=2*pi*70;        % Frequency of microphonic noise frequency effect in rad/sec
LFS=0;                  % Lorentz Force Switch (Turn On(1)/Off(0) Lorentz detuning effect)
MFS=0;                  % Microphonics Switch (Turn On(1)/Off(0) Microphonics detuning effect)
BS=1;                   % Beam Switch (Turn On(1)/Off(0) Beam
BPH=0;                  % Beam Phase
LS=1;                   % Open (0) or Closed(1) Loop Switch
PBIS=5;                 % Phase Bias
Igmax=4.5e-3;           % NOT USED
mix_co=5.73;            % Ask Curt if this has changed
%phase mixer convertion factor 0.1 V per one degree at 5 MV/m.
GLBG=60;               % Gradient Loop Broadband gain in dB
GLLG=30;                % Gradient Loop Low Frequency gain in dB
GT1=1/(2*pi*150);
GT2=1/(2*pi*1e6);
GT3=1/(2*pi*3e6);
PLBG=100;               % Phase Loop Broadband gain in dB
PLLG=30;                % Phase Loop Low Frequency Gain in dB
PT1=1/(2*pi*150);
PT2=1/(2*pi*1e6);
PT3=1/(2*pi*3e6);
%
% RUN SIMULATION
%
% ORIGINAL CALL FROM LIA'S DOIT_NEW SCRIPT.
%[t,x,y]=gear('cebaf_fdbk_new',2.0e-2,[],[1e-4 1e-4 1e-4 0 0 2]);
%
% gear is the integrator method from an older version of Simulink
% It can be called many ways.  The definition of the arguments
% for the call used here is:
% [time points, state trajectories, output trajectories] =
%            gear(
%              Simulink model name,
%              final time of simulation (starts at 0),
%              initial conditions -- not used here,
%              options [
%                relative error,
%                min step size,
%                max step size,
%                for gear if 1 automatic selection of algorithm depending on stiffness,
%                display parameter (3 = all, 0 = off),
%                plot parameter (1 = plot on, 0 = plot off)
%              ]
%            )
%[t,x,y]=gear('cebaf_fdbk_new_revised',2.0e-2,[],[1e-4 1e-4 1e-4 0 0 2]);
%
% Replacing gear call with sim call.  The definition of the arguments
% for the call used here is:
% [time points, state trajectories, output trajectories] =
%           sim(
%             Simulink model name,
%             final time of simulation (starts at 0),
%             options set using simset (relative error,fixed step solver)
%             external input values -- not used here
%             )
simOptions = simset('Solver','ode15s');
simOptions = simset(simOptions,'RelTol',1.0000e-06);
simOptions = simset(simOptions,'MaxStep',1.0000e-06,...
                    'FixedStep',1.0000e-06);
[t,x,y]=sim('USPAS2005',20.0e-3,simOptions,[]);
clear t x y;
%
% PLOT RESULTS
%
figure(5);
% Beam Amplitude (mA) vs. Time (ms)
subplot(3,2,1);
plot(tt*1e3,Ib*1e6,'r');
xlabel('time (ms)');
ylabel('Beam Amp (uA)');
hold on;
% Beam Phase (degrees) vs. Time (ms)
subplot(3,2,2);
plot(tt*1e3,BPH*180/pi,'r');
xlabel('time (ms)');
ylabel('Beam phase (degree)');
hold on;
% Cavity Amplitude (MV) vs. Time (ms)
subplot(3,2,3);
plot(tt*1e3,amp*1e-6,'r');
xlabel('time (ms)');
ylabel('Cavity amp (MV)');
hold on;
% Cavity Phase (degrees) vs. Time (ms)
subplot(3,2,4);
plot(tt*1e3,phase*180/pi,'r');
xlabel('time (ms)');
ylabel('Cavity phase (degree)');
hold on;
% Reference Amplitude (MV) vs. Time (ms)
%subplot(4,2,5);
%plot(tt*1e3,am_ref*1e-6,'r');
%xlabel('time (ms)');
%ylabel('ref amp (MV)');
%hold on;
% Reference Phase (degrees) vs. Time (ms)
%subplot(4,2,6);
%plot(tt*1e3, ph_ref*180/pi,'r');
%xlabel('time (ms)');
%ylabel('ref phase (degree)');
%hold on;
% Real Generator Current (mA) vs. Time (ms)
subplot(3,2,5);
plot(tt*1e3,Ig_I*1e3,'r');
xlabel('time (ms)');
ylabel('Ig real (mA)');
hold on;
% Imaginary Generator Current (mA) vs. Time (ms)
subplot(3,2,6);
plot(tt*1e3,Ig_Q*1e3,'r');
xlabel('time (ms)');
ylabel('Ig imag (mA)');
hold on;

% Zoomed figures
figure(6)
% Beam Amplitude (uA) vs. Time (ms)
subplot(3,2,1);
plot(tt*1e3,Ib*1e6,'r');
axis([ttmin ttmax Ibmin Ibmax]);
xlabel('time (ms)');
ylabel('Beam Amp (uA)');
hold on;
% Beam Phase (degrees) vs. Time (ms)
subplot(3,2,2);
plot(tt*1e3,BPH*180/pi,'r');
axis([ttmin ttmax BPHmin BPHmax]);
xlabel('time (ms)');
ylabel('Beam phase (degree)');
hold on;
% Cavity Amplitude (MV) vs. Time (ms)
subplot(3,2,3);
plot(tt*1e3,amp*1e-6,'r');
axis([ttmin ttmax ampmin ampmax]);
xlabel('time (ms)');
ylabel('Cavity amp (MV)');
hold on;
% Cavity Phase (degrees) vs. Time (ms)
subplot(3,2,4);
plot(tt*1e3,phase*180/pi,'r');
axis([ttmin ttmax phasemin phasemax]);
xlabel('time (ms)');
ylabel('Cavity phase (degree)');
hold on;
% Reference Amplitude (MV) vs. Time (ms)
%subplot(4,2,5);
%plot(tt*1e3,am_ref*1e-6,'r');
%xlabel('time (ms)');
%ylabel('ref amp (MV)');
%hold on;
% Reference Phase (degrees) vs. Time (ms)
%subplot(4,2,6);
%plot(tt*1e3, ph_ref*180/pi,'r');
%xlabel('time (ms)');
%ylabel('ref phase (degree)');
%hold on;
% Real Generator Current (mA) vs. Time (ms)
subplot(3,2,5);
plot(tt*1e3,Ig_I*1e3,'r');
axis([ttmin ttmax IgImin IgImax]);
xlabel('time (ms)');
ylabel('Ig real (mA)');
hold on;
% Imaginary Generator Current (mA) vs. Time (ms)
subplot(3,2,6);
plot(tt*1e3,Ig_Q*1e3,'r');
axis([ttmin ttmax IgQmin IgQmax]);
xlabel('time (ms)');
ylabel('Ig imag (mA)');
hold on;

bb=cputime;                % Used to compute length of time to run simulation (ending time)
bb-aa                      % length of time to run simulation (ending time - beginning time)