%% Synchrotronstrahlung: beta senkrecht Teilchenbeschleunigung

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

% Vers. 1.0 11/2014
% Vers. 3.1 03/2015

%% Winkelkoordinaten
theta = linspace(0,pi,50);%0:pi/40:pi;                     % Polar-Winkel 
phi = linspace(0,2*pi,50);%0:pi/20:2*pi;                   % Azimuth- Winkel

[Phi,Theta] = meshgrid(phi,theta);                         % Abbildungsgitter

r = 26658.883;                                             % m Radius LHC
c = 2.99792458e08;                                         % m/s Lichtgeschwindigkeit
q = 1.60217653e-19;                                            % As  p-Ladung
epsilon0 = 8.854187817e-12;                                % As/Vm 

vorf = q^2/(4*pi*c) * 1/(4*pi*epsilon0) * (c/r)^2;

%% Auswertung der Gleichung dP/dOmega
beta = 0.1;                                                % Teilchengeschwindigkeit
gamma2 = 1./(1 - beta^2);
cos2b = 1./(1 - beta*cos(Theta)).^2;
dP = beta^4 .* cos2b.^2.*(1 - (sin(Theta) .* cos(Phi)).^2 .* cos2b./gamma2); % Laenge in Richtung theta, phi plotten
dP(1)
max(abs(dP(:)))
skal01 = max(abs(dP(:)))*vorf
dP = dP ./ max(abs(dP(:)));                                % Skalieren auf 1
dP = abs(dP);
max(dP(:))

rho = dP .* sin(Theta);                                    % Euklidische Koordinaten
X = rho .* cos(Phi);
Y = rho .* sin(Phi);
Z = dP .* cos(Theta);
%
figure1=figure('renderer','ZBuffer');% Strahlungscharakterisitk plotten
%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);                            % Teilchenbahn und Koordinaten visualisieren 
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.25],'k','linewidth',1.5),plot3([-0.7,0],[0,0],[0,0],'k','linewidth',1.5)
plot3([0.,0],[0,1],[0,0],'k','linewidth',1.5)
xlabel('x'), ylabel('y'),zlabel('z')
axis off
view(-130.5,16)                                            % Abbildungen ausrichten
colormap(gray)                                             % Farbwahl
brighten(figure1,1)                                        % Oberflaechengestaltung
light
lighting gouraud

%% beta 0.6
ah2=axes('position', [0.1500    0.1100    0.7750    0.8150])
beta = 0.6;
gamma2 = 1./(1 - beta^2);
cos2b = 1./(1 - beta*cos(Theta)).^2;
dP = beta^4 .* cos2b.^2.*(1 - (sin(Theta) .* cos(Phi)).^2 .* cos2b./gamma2); % l??nge in Richtung theta, phi plotten
dP(1)
max(abs(dP(:)))
skal06 = max(abs(dP(:)))*vorf
dP = dP ./ max(abs(dP(:)));
dP = abs(dP);
max(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.15],'k','linewidth',1.5),plot3([-0.3,0],[0,0],[0,0],'k','linewidth',1.5)
plot3([0.,0],[0,0.3],[0,0],'k','linewidth',1.5)
xlabel('x'), ylabel('y'),zlabel('z')
view(-130.5,16)
colormap(ah2,gray)
brighten(ah2,1)
light
lighting gouraud
axis off


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

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

beta = 0.9;
gamma2 = 1./(1 - beta^2);
cos2b = 1./(1 - beta*cos(Theta)).^2;
dP = beta^4 .* cos2b.^2.*(1 - (sin(Theta) .* cos(Phi)).^2 .* cos2b./gamma2); % l??nge in Richtung theta, phi plotten
dP(1)
max(abs(dP(:)))
skal09 = max(abs(dP(:)))*vorf
dP = dP ./ max(abs(dP(:)));
dP = abs(dP);
max(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(-pi/12,pi/3);
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.15],'k','linewidth',1.5),plot3([-0.2,0],[0,0],[0,0],'k','linewidth',1.5)
plot3([0.,0],[0,0.1],[0,0],'k','linewidth',1.5)
xlabel('x'), ylabel('y'),zlabel('z')
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)