%================================================================
%  NLC simple example
%  un-ideal tuned machine, two beamlines
%================================================================
%
clear all;
close all;

% electrons after IP
fname = input('Enter file name of beam1 data: ','s');
beam1      = load(fname);
pbeam=250.;
w0_laser=1.e-3; %1000 micron laser w0
yoff=-2.; %nominal y offset is -2cm at chicane
R22=-0.595; %TRANSPORT from IP to mid-chicane
R44=-0.443; %TRANSPORT from IP to mid-chicane
i1=[1 0 0 0 0 0];
i2=[0 1 0 0 0 0];
i3=[0 0 1 0 0 0 ];
i4=[0 0 0 1 0 0];
i5=[0 0 0 0 1 0];
i6=[0 0 0 0 0 1];
e_beam1=i1*beam1';
angx_beam1=i2*beam1'*1e6; %units are micro-radians
angy_beam1=i3*beam1'*1e6; %units are micro-radians
ang_beam1=sqrt(angx_beam1.^2+angy_beam1.^2); %units are micro-radians
DeltaP=(1.-cos((pbeam/0.44065)*ang_beam1*1.e-6));
z_beam1=i4*beam1'*1e-4; %units are cm
x_beam1=i5*beam1'*1e-7; %units are cm
y_beam1=(i6*beam1'-yoff)*1e-7; %units are cm
r2_beam1=1.e-4*(x_beam1.^2+y_beam1.^2); %units are m^2
laserwt=exp(-2*r2_beam1/w0_laser^2);
DeltaP_laserwt=(1./sum(laserwt))*laserwt*DeltaP';
% printout some results on means, stds
mean_e_beam1=(mean(e_beam1)-pbeam)/250*1e6 % in ppm
mean_x_beam1=mean(x_beam1)
mean_y_beam1=mean(y_beam1)
mean_z_beam1=mean(z_beam1)
mean_angx_beam1=mean(angx_beam1)
mean_angy_beam1=mean(angy_beam1)
std_e_beam1=std(e_beam1)/pbeam*1e6 % in ppm
std_x_beam1=std(x_beam1)
std_y_beam1=std(y_beam1)
std_z_beam1=std(z_beam1)
std_angx_beam1=std(angx_beam1)
std_angy_beam1=std(angy_beam1)
mean_DeltaP=mean(DeltaP)
DeltaP_laserwt

% PLOTS
subplot(1,2,1)
hist(abs(e_beam1),100);
grid on
%axis([0 20 -0.025 0.0]);
xlabel('Beam1 Energy (GeV)','fontsize',14);
%set(gca,'yscale','log')

subplot(1,2,2)
hist(abs(e_beam1),2000);
grid on
axis([247.5 252.5 0 500]);
xlabel('Beam1 Energy (GeV)','fontsize',14);
%set(gca,'yscale','log')

figure
subplot(2,2,1)
hist(angx_beam1,100);
grid on
%axis([0 20 -0.025 0.0]);
xlabel('Electron X Angle (urad)','fontsize',14);

subplot(2,2,2)
hist(angx_beam1,1000);
grid on
axis([-1000 1000 0 500]);
xlabel('Electron X Angle (urad)','fontsize',14);

subplot(2,2,3)
hist(angy_beam1,100);
grid on
%axis([0 20 -0.025 0.0]);
xlabel('Electron Y Angle (urad)','fontsize',14);

subplot(2,2,4)
hist(angy_beam1,1000);
grid on
axis([-200 200 0 250]);
xlabel('Electron Y Angle (urad)','fontsize',14);

figure
subplot(1,2,1)
plot(sqrt(r2_beam1)*1e3,laserwt,'r.');
grid on
axis([0 5 0 1]);
xlabel('beam radius (mm)','fontsize',14);
ylabel('laser weight','fontsize',14);

subplot(1,2,2)
hist(DeltaP,100);
grid on
%axis([0 20 -0.025 0.0]);
xlabel('Avg depolarization','fontsize',14);

figure
subplot(2,2,1)
hist(x_beam1,100);
grid on
%axis([0 20 -0.025 0.0]);
xlabel('Electron X (cm)','fontsize',14);

subplot(2,2,2)
hist(x_beam1,5000);
grid on
axis([-0.1 0.1 0 2000]);
xlabel('Electron X (cm)','fontsize',14);

subplot(2,2,3)
hist(y_beam1,100);
grid on
%axis([0 20 -0.025 0.0]);
xlabel('Electron Y (cm)','fontsize',14);

subplot(2,2,4)
hist(y_beam1,5000);
grid on
axis([-0.1 0.1 0 500]);
xlabel('Electron Y (cm)','fontsize',14);

figure
hist(z_beam1,100);
grid on
%axis([0 20 -0.025 0.0]);
xlabel('Electron Z (cm)','fontsize',14);

figure
subplot(2,2,1)
plot(x_beam1,y_beam1,'r.')
grid on
xlabel('Electron X (cm)','fontsize',14);
ylabel('Electron Y (cm)','fontsize',14);

subplot(2,2,2)
plot(x_beam1,y_beam1,'r.')
axis([-0.1 0.1 -0.1 0.1]);
grid on
xlabel('Electron X (cm)','fontsize',14);
ylabel('Electron Y (cm)','fontsize',14);

subplot(2,2,3)
plot(angx_beam1,angy_beam1,'r.')
grid on
xlabel('Electron X Angle (urad)','fontsize',14);
ylabel('Electron Y Angle (urad)','fontsize',14);

subplot(2,2,4)
plot(angx_beam1,angy_beam1,'r.')
axis([-500 500 -200 200]);
grid on
xlabel('Electron X Angle (urad)','fontsize',14);
ylabel('Electron Y Angle (urad)','fontsize',14);

figure
subplot(2,2,1)
plot(x_beam1,angx_beam1,'r.')
grid on
xlabel('Electron X (cm)','fontsize',14);
ylabel('Electron X Angle (urad)','fontsize',14);

subplot(2,2,2)
plot(x_beam1,angx_beam1,'r.')
axis([-0.1 0.1 -500 500]);
grid on
xlabel('Electron X (cm)','fontsize',14);
ylabel('Electron X Angle (urad)','fontsize',14);

subplot(2,2,3)
plot(y_beam1,angy_beam1,'r.')
grid on
xlabel('Electron Y (cm)','fontsize',14);
ylabel('Electron Y Angle (urad)','fontsize',14);

subplot(2,2,4)
plot(y_beam1,angy_beam1,'r.')
axis([-0.1 0.1 -200 200]);
grid on
xlabel('Electron Y (cm)','fontsize',14);
ylabel('Electron Y Angle (urad)','fontsize',14);

figure
subplot(2,2,1)
plot(e_beam1,x_beam1,'r.')
axis([50 250 -4 4])
grid on
xlabel('Electron Energy (GeV)','fontsize',14);
ylabel('Electron X (cm)','fontsize',14);

subplot(2,2,2)
plot(e_beam1,x_beam1,'r.')
axis([247.5 252.5 -0.04 +0.04]);
grid on
xlabel('Electron Energy (GeV)','fontsize',14);
ylabel('Electron X (cm)','fontsize',14);

subplot(2,2,3)
plot(e_beam1,y_beam1,'r.')
axis([50 250 -8 2])
grid on
xlabel('Electron Energy (GeV)','fontsize',14);
ylabel('Electron Y (cm)','fontsize',14);

subplot(2,2,4)
plot(e_beam1,y_beam1,'r.')
axis([247.5 252.5 -0.04 +0.04]);
grid on
xlabel('Electron Energy (GeV)','fontsize',14);
ylabel('Electron Y (cm)','fontsize',14);

figure
subplot(1,2,1)
plot(e_beam1,angx_beam1,'r.')
axis([50 250 -1000 800])
grid on
xlabel('Electron Energy (GeV)','fontsize',14);
ylabel('Electron X Angle (urad)','fontsize',14);

subplot(1,2,2)
plot(e_beam1,angy_beam1,'r.')
axis([50 250 -1000 1000])
grid on
xlabel('Electron Energy (GeV)','fontsize',14);
ylabel('Electron Y Angle (urad)','fontsize',14);

%print -dpdf -f1 beam1_chicane_1.pdf
%print -dpdf -f2 beam1_chicane_2.pdf
%print -dpdf -f3 beam1_chicane_3.pdf
%print -dpdf -f4 beam1_chicane_4.pdf
%print -dpdf -f5 beam1_chicane_5.pdf
%print -dpdf -f6 beam1_chicane_6.pdf
%print -dpdf -f7 beam1_chicane_7.pdf
%print -dpdf -f8 beam1_chicane_8.pdf
%print -dpdf -f9 beam1_chicane_9.pdf

%
%exit