% 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,47];
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




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

%R = cos(2*TH)*cos(29*T/3) + cos(18*TH)*cos(61*T/3) + cos(26*TH)*cos(83*T/3);

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')


