% Want to just read in the psis and qs calculated previously and use them to
% calculate g(r) and S(k) for given values of rho

% This is the version of the code used by Adda-Bedia et al., J Chem. Phys.,
% 129, 144506 (2008) to study higher EVEN dimensions. Here it is set for
% the case of two dimensions.

m=0; % this works in d=2*(m+1) dimensions so here d=2

maxits=10; % The Maximum number of iterations
N=maxits; % the number of piecewise bits that we actually want to take

k=1001; % The number of points within each piecewise section

dt=1/(k-1); % The gap between points

psidata=load('2D_DX1000psi.txt');
qdata=load('2D_DX1000q.txt');
psi=psidata(1:(N*(k-1)+1),2:(N+1));
Q=qdata(1:k,2:(N+1)); % Need to be very careful about the difference between q and Q !

% Here is where we define the values of rho at which the various quantities are
% calculated
rho1=0.1;
rhos1=(2*pi)^m*rho1 % remember that the series is lambda, not rho
g1=ones((N-1)*(k-1)+1,1);

rho2=0.3;
rhos2=(2*pi)^m*rho2
g2=ones((N-1)*(k-1)+1,1);

rho3=0.5;
rhos3=(2*pi)^m*rho3
g3=ones((N-1)*(k-1)+1,1);

disp('I am now calculating g(r)')

for a=1:(maxits-1)
	for i=0:(a*(k-1))
		r=1+i*dt;
		temp=-0.5*dt*(a+1)*psi(i+1,a+1)/sqrt((a+1)^2-r^2)+psi(i+1,a+1)*sqrt((a+1)^2-r^2);
		for j=(i+1):(a*(k-1)-1)
			temp=temp+dt*(1+j*dt)/sqrt((1+j*dt)^2-r^2)*(psi(j+1,a+1)-psi(i+1,a+1));
		end
		g1(i+1)=g1(i+1)+rhos1^a*temp/pi;
		g2(i+1)=g2(i+1)+rhos2^a*temp/pi;
		g3(i+1)=g3(i+1)+rhos3^a*temp/pi;
    end
	disp(strcat('I am on iteration no. ',num2str(a)))
end

x=1+(0:((N-1)*(k-1)))*dt;
figure(1)
plot(x,g1,'r-','linewidth',2)
hold on
plot(x,g2,'g-','linewidth',2)
q=plot(x,g3,'b-','linewidth',2);
axis([0.75 5 0.5 3])

set(get(q,'parent'),'FontSize',30);
xlabel('r','FontSize',40)
ylabel('g(r)','FontSize',40)
text(3.25,1.9,['D=2'],'FontSize',40)
set(get(q,'parent'),'linewidth',3)
set(get(q,'parent'),'ticklength',[0.02 0.05])
title('The pair correlation function')
hold off

disp('I am now calculating S(k)')

dq=0.25;
maxq=20;
npts=maxq/dq;
q=(dq:dq:maxq)';

Qc=zeros(npts,N+1);   % The cosine transform of the components of Q
Qs=zeros(npts,N+1);   % The sine transform of the components of Q
for n=0:(N-1)
    Qc(:,n+1)=0.5*dt*Q(1,n+1);
    for j=1:(k-2)
        Qc(:,n+1)=Qc(:,n+1)+dt*Q(j+1,n+1)*cos(j*dt*q);
        Qs(:,n+1)=Qs(:,n+1)+dt*Q(j+1,n+1)*sin(j*dt*q);
    end
end

qcval1=ones(npts,1);
qsval1=zeros(npts,1);
qcval2=ones(npts,1);
qsval2=zeros(npts,1);
qcval3=ones(npts,1);
qsval3=zeros(npts,1);


for n=0:(N-1)
    qcval1=qcval1-rhos1^(n+1)*Qc(:,n+1);
    qsval1=qsval1-rhos1^(n+1)*Qs(:,n+1);
    qcval2=qcval2-rhos2^(n+1)*Qc(:,n+1);
    qsval2=qsval2-rhos2^(n+1)*Qs(:,n+1);
    qcval3=qcval3-rhos3^(n+1)*Qc(:,n+1);
    qsval3=qsval3-rhos3^(n+1)*Qs(:,n+1);
end

qval1=qsval1.^2+qcval1.^2;
s1=1./qval1;
qval2=qsval2.^2+qcval2.^2;
s2=1./qval2;
qval3=qsval3.^2+qcval3.^2;
s3=1./qval3;


figure(2)
plot(q,s1,'r-','linewidth',2)
hold on
plot(q,s2,'g-','linewidth',2)
r=plot(q,s3,'b-','linewidth',2);
set(get(r,'parent'),'FontSize',30);
text(3.25,1.9,['D=2'],'FontSize',40)
xlabel('k','FontSize',40)
ylabel('S(k)','FontSize',40)
set(get(r,'parent'),'linewidth',3)
set(get(r,'parent'),'ticklength',[0.02 0.05])
title('Fourier transform of the structure function')
hold off



