% Dirac Talbot with octave


%

NT = 1024; % Number of nodes in time and angle

% Momentum list
momL = [1,2,6,13];
momL = [momL,-momL];
% q^2
Q = 160/9

% Momentum list
momL = [2,18,26];
momL = [momL,-momL];
% q^2
Q = 805/9

% Momentum list
momL = [3,5,15];
momL = [momL,-momL];
% q^2
Q = 189/4

trev = 2*pi/sqrt(momL(1)^2+Q)*15/1
%trev=6*pi


th_r = linspace(-pi,pi,NT);
t_r = linspace(0,1,NT);
[T,TH] = meshgrid(t_r,th_r);


S1 = 0;
S2 = 0;
for m = momL
m
	wm = sqrt(m^2+Q);
	% normalization coefficient
	normc = sqrt(1+ sqrt(Q)/wm)/sqrt(2);
	normc2 = normc*i*m/(sqrt(Q)+wm);
	S1 = S1+ normc*exp(i*m*TH-i*wm*T*trev);
	S2 = S2+ normc2*exp(i*m*TH-i*wm*T*trev);
end


%S2=0

R = sqrt(abs(S1).^2 + abs(S2).^2);

%R = cos(3*TH)*cos(15*T*trev/2) + cos(5*TH)*cos(17*T*trev/2) + cos(15*TH)*cos(33*T*trev/2);
%R = abs(exp(i*3*TH)*exp(-i*15*T*trev/2) + cos(5*TH)*exp(-i*17*T*trev/2) + cos(15*TH)*exp(-i*33*T*trev/2));
%R = abs(cos(3*TH)*exp(-i*15*T*trev/2) + cos(5*TH)*exp(-i*17*T*trev/2) + cos(15*TH)*exp(-i*33*T*trev/2));

figure(1)
surf(T,TH,R)
colormap(jet(256));
%colormap(gray);
shading interp;
view(0,90)
colorbar
axis([0 1 -pi pi])
grid off
%saveas(gcf,'../../images/carpet1_5ts.png')
saveas(gcf,'./borri.png')


