function compute_potential_Matlab

% Constants
G=6.67428*10^-8; % cm^3 g^-1 s^-2
kpc=3.085678*10^21; % 1kpc = 3.08 10^21 cm
mass=10^7*2.0*10^33; % 10^7 solar mass = 2 10^33 g with ~10^4 stars
% Factor to remain into float range with C Language (in BmyOGG code)
factor_fft = 2*10^40;
% Loading Positions
data=load('Positions_Disk.dat');
% Number of star coordinates : 10240 with Positions_Disk.dat
numStars=numel(data(:,1));
% Size of the edges of mesh
numPoints=256;
% Min and Max for adapting mesh 
x0=linspace(min(data(:,1)),max(data(:,1)),numPoints);
y0=linspace(min(data(:,2)),max(data(:,2)),numPoints);
z0=linspace(min(data(:,3)),max(data(:,3)),numPoints);
% Generating mesh for direct computation
[x y z]=meshgrid(x0,y0,z0);
% Initialization of potential Phi
Phi(1:numPoints,1:numPoints,1:numPoints)=0.0;
% Start time for direct method
tic;
% Direct computing of gravitational potential
for j=1:numPoints
   for k=1:numPoints
      for l=1:numPoints
         Phi(j,k,l)=0.0;
         for i=1:numStars
            Phi(j,k,l)=Phi(j,k,l)-G*mass/((sqrt((x(j,k,l)-data(i,1))^2+(y(j,k,l)-data(i,2))^2+(z(j,k,l)-data(i,3))^2))*kpc);
         end
         Value(j,k,l)=abs(Phi(j,k,l));
      end
   end
end
% indexes for finding minimum values
index_x=find(abs(x0)==min(abs(x0)));
index_y=find(abs(y0)==min(abs(y0)));
index_z=find(abs(z0)==min(abs(z0)));
% Plot potential computed with direct method
[xslice yslice zslice]=meshgrid(1:numPoints,1:numPoints,1:numPoints);
hFig=figure(1);
set(hFig, 'Position', [400 400 750 600]);
slice(xslice,yslice,zslice,Value,[index_x numPoints], [index_y numPoints], [1 index_z]);
shading interp;
view([-42,22]);
hc=colorbar;
set(hc,'position',[0.932 0.3 0.02 0.6]);
caxis([10^14 10^15]);
xlim([0 numPoints]);
ylim([0 numPoints]);
zlim([0 numPoints]);
xlabel('x domain');
ylabel('y domain');
zlabel('z domain');
colormap(jet(1000));
title('Gravitational potential with direct computation');
% End time for direct method
toc;
% Average value and standard deviation
disp('Direct computation : average value and standard deviation');
disp(mean(Value(:)));
disp(std(Value(:)));
% Start time for FFT Matlab method
tic;
% Compute the spacing of all cells
Max_x=max(data(:,1));
Min_x=min(data(:,1));
Max_y=max(data(:,2));
Min_y=min(data(:,2));
Max_z=max(data(:,3));
Min_z=min(data(:,3));
if (Max_x>abs(Min_x))
  space_x=Max_x/(numPoints/2);
  add_correction_x=0;
else
  space_x=abs(Min_x)/(numPoints/2);
  add_correction_x=1;
end
if (Max_y>abs(Min_y))
  space_y=Max_y/(numPoints/2);
  add_correction_y=0;
else
  space_y=abs(Min_y)/(numPoints/2);
  add_correction_y=1;
end
if (Max_z>abs(Min_z))
  space_z=Max_z/(numPoints/2);
  add_correction_z=0;
else
  space_z=abs(Min_z)/(numPoints/2);
  add_correction_z=1;
end
% Initialization of Green and density
green_grid(1:numPoints,1:numPoints,1:numPoints)=0.0;
density(1:numPoints,1:numPoints,1:numPoints)=0.0;
% Computing the density for each cell with CIC method
for i=1:numStars
   % Computing node positions on mesh
   node_x = floor(data(i,1)/space_x) + (numPoints/2 + add_correction_x);
   node_y = floor(data(i,2)/space_y) + (numPoints/2 + add_correction_y);
   node_z = floor(data(i,3)/space_z) + (numPoints/2 + add_correction_z);
   % Set the CIC size fractions using fmod 
   dx = mod(abs(data(i,1)),space_x)/space_x;
   dy = mod(abs(data(i,2)),space_y)/space_y;
   dz = mod(abs(data(i,3)),space_z)/space_z;
   tx = 1.0 - dx;
   ty = 1.0 - dy;
   tz = 1.0 - dz;
   % Adding for setting CIC size fractions
   % Version with (i,j) centered on cells
   if ((dx < 0.5) && (node_x == 1))
     frac_x = dx;
     neighboor_x = 0;
     epsilon_x = 0;
   elseif ((dx > 0.5) && (node_x == numPoints))
     frac_x = dx - 0.5;
     neighboor_x = 0;
     epsilon_x = 0;
   elseif (dx < 0.5) 
     frac_x = dx + 0.5;
     neighboor_x = -1;
     epsilon_x = 1;
   else 
     frac_x = tx + 0.5;
     neighboor_x = 1;
     epsilon_x = 1;
   end
   if ((dy < 0.5) && (node_y == 1))
     frac_y = dy;
     neighboor_y = 0;
     epsilon_y = 0;
   elseif ((dy > 0.5) && (node_y == numPoints))
     frac_y = dy - 0.5;
     neighboor_y = 0;
     epsilon_y = 0;
   elseif (dy < 0.5) 
     frac_y = dy + 0.5;
     neighboor_y = -1;
     epsilon_y = 1;
   else 
     frac_y = ty + 0.5;
     neighboor_y = 1;
     epsilon_y = 1;
   end
   if ((dz < 0.5) && (node_z == 1))
     frac_z = dz;
     neighboor_z = 0;
     epsilon_z = 0;
   elseif ((dz > 0.5) && (node_z == numPoints))
	frac_z = dz - 0.5;
	neighboor_z = 0;
	epsilon_z = 0;
   elseif (dz < 0.5) 
	frac_z = dz + 0.5;
	neighboor_z = -1;
	epsilon_z = 1;
   else 
	frac_z = tz + 0.5;
	neighboor_z = 1;
	epsilon_z = 1;
   end
   % Density with CIC fractions
   density(node_x,node_y,node_z) = density(node_x,node_y,node_z) + mass*(frac_x*frac_y*frac_z)/(space_x*space_y*space_z);
   density(node_x+neighboor_x,node_y,node_z) = density(node_x+neighboor_x,node_y,node_z) + mass*(epsilon_x*(1.0-frac_x)*frac_y*frac_z)/(space_x*space_y*space_z);
   density(node_x,node_y+neighboor_y,node_z) = density(node_x,node_y+neighboor_y,node_z) + mass*(frac_x*epsilon_y*(1.0-frac_y)*frac_z)/(space_x*space_y*space_z);
   density(node_x,node_y,node_z+neighboor_z) = density(node_x,node_y,node_z+neighboor_z) + mass*(frac_x*frac_y*epsilon_z*(1.0-frac_z))/(space_x*space_y*space_z);
   density(node_x,node_y+neighboor_y,node_z+neighboor_z) = density(node_x,node_y+neighboor_y,node_z+neighboor_z) + mass*(frac_x*epsilon_y*(1.0-frac_y)*epsilon_z*...
   (1.0-frac_z))/(space_x*space_y*space_z);
   density(node_x+neighboor_x,node_y+neighboor_y,node_z) = density(node_x+neighboor_x,node_y+neighboor_y,node_z) + mass*(epsilon_x*(1.0-frac_x)*epsilon_y*...
   (1.0-frac_y)*frac_z)/(space_x*space_y*space_z);
   density(node_x+neighboor_x,node_y,node_z+neighboor_z) = density(node_x+neighboor_x,node_y,node_z+neighboor_z) + mass*(epsilon_x*(1.0-frac_x)*frac_y*epsilon_z*...
   (1.0-frac_z))/(space_x*space_y*space_z);
   density(node_x+neighboor_x,node_y+neighboor_y,node_z+neighboor_z) = density(node_x+neighboor_x,node_y+neighboor_y,node_z+neighboor_z) + mass*(epsilon_x*(1.0-frac_x)*...
   epsilon_y*(1.0-frac_y)*epsilon_z*(1.0-frac_z))/(space_x*space_y*space_z);
end
% Important for flip (i,j) = (line,row) => (x,y) axes
p = [2 1 3];
density_flip = permute(density,p);
% Computing discretized Green's function
for i=1:numPoints
   for j=1:numPoints
      for k=1:numPoints
	 dx = sqrt((i-numPoints/2-add_correction_x)*(i-numPoints/2-add_correction_x))*space_x;
	 dy = sqrt((j-numPoints/2-add_correction_y)*(j-numPoints/2-add_correction_y))*space_y;
	 dz = sqrt((k-numPoints/2-add_correction_z)*(k-numPoints/2-add_correction_z))*space_z;
	 % Discretized Green's function  
	 green_grid(i,j,k) = -1.0/(4.0*pi*sqrt(dx*dx + dy*dy + dz*dz));
      end
   end
end
% For avoiding divergence at numPoints/2
green_grid(numPoints/2+add_correction_x,numPoints/2+add_correction_y,numPoints/2+add_correction_z) = 1.0;
% FFT of Green 
fft_green=fftn(green_grid);
% FFT of density 
fft_density=fftn(density_flip);
% Product in Fourier space 
product_fft=4*pi*G*fft_green.*fft_density;
% Backward of product above
fft_potential=real(ifftn(product_fft))*space_x*space_y*space_z/kpc;
% Shift final solution
fft_potential_shift=fftshift(fft_potential);
% Plot slice of potential
hFig=figure(2);
set(hFig, 'Position', [400 400 750 600]);
slice(xslice,yslice,zslice,abs(fft_potential_shift),[numPoints/2+add_correction_x numPoints], [numPoints/2+add_correction_y, numPoints], [1 numPoints/2+add_correction_z]);
shading interp;
view([-42,22]);
hc=colorbar;
set(hc,'position',[0.932 0.3 0.02 0.6]);
caxis([10^14 10^15]);
xlim([0 numPoints]);
ylim([0 numPoints]);
zlim([0 numPoints]);
xlabel('x domain');
ylabel('y domain');
zlabel('z domain');
colormap(jet(1000));
title('Gravitational potential with FFT/Matlab computation');
% End time for FFT Matlab method
toc;
% Average value and standard deviation
disp('FFT/Matlab computation : average value and standard deviation');
disp(mean(abs(fft_potential_shift(:))));
disp(std(abs(fft_potential_shift(:))));
end