%% Linearbeschleuniger: beta parallel Beschleunigung

% W. Schweizer
% Simulation physikalischer Systeme - Computational Physics mit MATLAB 
% 2016

% Vers. 1.1 03/2015

theta = linspace(0,pi,50);%0:pi/40:pi;                   % polar angle
phi = linspace(0,2*pi,50);%0:pi/20:2*pi;                   % azimuth angle

[Phi,Theta] = meshgrid(phi,theta);    % define the grid


%%
beta = 0.1;
gamma2 = 1 - beta^2;
cos6b = 1./(1 - beta*cos(Theta)).^6;
dP =sin(Theta).^2 .* cos6b; % l??nge in Richtung theta, phi plotten
max(abs(dP(:)))
dP = dP ./ max(abs(dP(:)));
max(abs(dP(:)))
dP = abs(dP);
rho = dP .* sin(Theta);
X = rho .* cos(Phi);
Y = rho .* sin(Phi);
Z = dP .* cos(Theta);

figure1=figure;
axes('position', [0.0500    0.1100    0.7750    0.8150])
hold on
surf(X,Y,Z), %axis equal  
shg
hold on
% % psi = linspace(-pi/4,11*pi/24);
% % x3 = max(dP(:)) .* (cos(psi) - 1);
% % z3 = max(dP(:)) .* sin(psi);
% % y3 = zeros(size(x3));
%plot3(x3,y3,z3,'k --','linewidth',2.)
plot3([0,0],[0,0],[0,1.05],'k','linewidth',1.5),plot3([-1.2,0],[0,0],[0,0],'k','linewidth',1.5)
plot3([0.,0],[0,1.2],[0,0],'k','linewidth',1.5)
text([0,0],[0,0.1],[0,0.85],'z','FontSize',14,'FontWeight','bold')
% xlabel('x'), ylabel('y'),zlabel('z')
view(-130.5,16)
colormap(gray)
brighten(figure1,1)
light
lighting gouraud
axis off

%%
ah2=axes('position', [0.1500    0.1100    0.7750    0.8150])
beta = 0.6;
gamma2 = 1 - beta^2;
cos6b = 1./(1 - beta*cos(Theta)).^6;
dP =sin(Theta).^2 .* cos6b; % l??nge in Richtung theta, phi plotten
max(abs(dP(:)))
dP = dP ./ max(abs(dP(:)));
dP = abs(dP);
rho = dP .* sin(Theta);
X = rho .* cos(Phi);
Y = rho .* sin(Phi);
Z = dP .* cos(Theta);

%figure2=figure;
surf(X,Y,Z), %axis equal  
shg
hold on
% psi = linspace(-2*pi/24,7*pi/24);
% x3 = max(dP(:)) .* (cos(psi) - 1);
% z3 = max(dP(:)) .* sin(psi);
% y3 = zeros(size(x3));
% plot3(x3,y3,z3,'k --','linewidth',2.)
plot3([0,0],[0,0],[0,1.05],'k','linewidth',1.5),plot3([-1.2,0],[0,0],[0,0],'k','linewidth',1.5)
plot3([0.,0],[0,1.2],[0,0],'k','linewidth',1.5)
text([0,0],[0,0.1],[0,0.85],'z','FontSize',14,'FontWeight','bold')
% xlabel('x'), ylabel('y'),zlabel('z')
view(-130.5,16)
colormap(ah2,gray)
brighten(ah2,1)
light
lighting gouraud
axis off


% annotation(figure1,'arrow',[0. 0.],...
%     [0. 1.]);
% figure
% surf(phi,theta,dP),shg
%figure, plot(sin(theta),cos(theta),'*')

%%
ah3=axes('position', [0.2500    0.1100    0.7750    0.8150])
theta = linspace(0,pi,100);%0:pi/40:pi;                   % polar angle
phi = linspace(0,2*pi,25);%0:pi/20:2*pi;                   % azimuth angle

[Phi,Theta] = meshgrid(phi,theta);    % define the grid

beta = 0.9;
gamma2 = 1 - beta^2;
cos6b = 1./(1 - beta*cos(Theta)).^6;
dP =sin(Theta).^2 .* cos6b; % l??nge in Richtung theta, phi plotten
max(abs(dP(:)))
dP = dP ./ max(abs(dP(:)));
dP = abs(dP);
rho = dP .* sin(Theta);
X = rho .* cos(Phi);
Y = rho .* sin(Phi);
Z = dP .* cos(Theta);

%figure3=figure;
surf(X,Y,Z), %axis equal  
shg
hold on
% psi = linspace(-6*pi/24,16*pi/24);
% x3 = max(dP(:)) .* (cos(psi) - 1);
% z3 = max(dP(:)) .* sin(psi);
% y3 = zeros(size(x3));
% plot3(x3,y3,z3,'k --','linewidth',2.)
plot3([0,0],[0,0],[0,1.05],'k','linewidth',1.5),plot3([-1.2,0],[0,0],[0,0],'k','linewidth',1.5)
plot3([0.,0],[0,1.2],[0,0],'k','linewidth',1.5)
text([0,0],[0,0.1],[0,1.05],'z','FontSize',14,'FontWeight','bold')
view(-130.5,16)
colormap(ah3,gray)
brighten(ah3,1)
light
lighting gouraud
axis off


%%  siehe Erlaeuterung Buch zum Nutzen der Camera Toolbar (Figure unter View)