% Dirac Talbot with octave


%

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

% Momentum list
momL = [3,19,39,71];
momL = [momL,-momL];
% q^2
Q = 791




trev = 2*pi/sqrt(momL(1)^2+Q)*10
%trev = 3.6276

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


S1 = 0;
S2 = 0;
disp('--')
for m = momL
	wm = sqrt(m^2+Q);
(wm/sqrt(momL(1)^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);
%	S1 = S1+ normc*exp(i*m*TH-i*wm*T*trev)/sqrt(abs(m));
%	S2 = S2+ normc2*exp(i*m*TH-i*wm*T*trev)/sqrt(abs(m));
	S1 = S1+ normc*exp(i*m*TH-i*wm*T*trev)/(abs(m)).^0.25;
	S2 = S2+ normc2*exp(i*m*TH-i*wm*T*trev)/(abs(m)).^0.25;
end


R = sqrt(abs(S1).^2 + abs(S2).^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,'borri.png')


