Why does my function that should modify cells from an input return exactly the input every time?

I am attempting to model diffusion for an MRI project. In doing so, I created a function called diffuseCells() that takes several inputs and should modify the contents of the cells according to the rules in the functions. When I work through an individual example in the Command Window, everything seems to work. However, when I run the function it returns a cell structure that is exactly the same as the input structure. Any help would be greatly appreciated - I have worked on this single bug for hours and can't find the error. At this point I'm more or less banging my head against the wall. The function is below, and the entire script is below that.

The offending function:

%Inputs: cell array of protons, phantom, Gx(t),Gy(t), b0, dt,t
%Output: cell array of protons for t+dt
function t_cell = diffuseCells(p_cell,p, G1, G2,b0,dt,t)
%Calculate gradients 
Gxc = G1(t)*linspace(round(-1*size(p_cell,2)/2), round(size(p_cell,2)/2), size(p_cell,2));
Gyc = G2(t)*linspace(round(-1*size(p_cell,1)/2), round(size(p_cell,1)/2), size(p_cell,1));

%Initialize temp cell storage for randomwalk calc
t_cell = cell(size(p_cell));  % Initialize 
for i = 1:numel(p_cell)
    t_cell{i} = p_cell{i};  % Copy each element 
end
%t_cell = p_cell;
%at each timepoint iterate through cells and update spins and phases,
%then update position according to random walk
for i=1:size(p_cell,2) % x coord
    for j=size(p_cell,1) % y coord
        a = cell2mat(p_cell(i,j));
        %update phase1 = mod(1,phase0 + spin (Hz) * dt (uS)),1), spin=4257 Hz/G * (b0 + Gx +Gy) (G) =
        %gamma*field strength
        a(1,:) = 4257*(b0+Gxc(i)+Gyc(j));
        a(2,:) = mod(a(2,:)+a(1,:)*dt/(10^6),1);
        
        %Storage for a after random walk
        a_r=a;
        %update position according to random walk
        for k=1:size(a,2)
            current = a(:,k);
            [dx,dy] = randomwalk(p(i,j)); %use diffusion coefficients stored in phantom as n for random walk
            if(p(i+dx,j+dy) == 0) %If moving out of object, stop
                dx=0;
                dy=0;
            end
            
            %move proton if it needs to be moved
            if(dx ~= 0 || dy ~= 0)
                %place in new cell
                %t_cell{i+dx,j+dy} = [];
                if(i+dx <= size(p_cell,2) && i+dx >= 1 && y+dy <= size(p_cell,1) && y+dy >=1)
                    t_cell{i+dx,j+dy} = [cell2mat(p_cell(i+dx,j+dy)),current];
                end
                %remove from old cell
                a_r(:,k) = [];
            end
        end
        %update cell with new proton entries
        %t_cell{i,j} = [];
        t_cell{i,j} = a_r;
    end
end
end

Entire script:

%%
%define phantom, convert to cells
pt = phantom('Modified Shepp-Logan', 256);
%imagesc(p)
%Process phantom to have pixel values as diffusion coefficients
%unique(pt) = -0.0000,0,0.1,0.2,0.3,0.4,1 so multiply by ten to get
%diffusion coefficients 
%NOTE: divide output by ten when displaying phantom to get good contrast
p = round(pt*10);
imagesc(pt/10)
p_c = makeCells(p);
%Define gradients and dimensions for k space traversal
FOVx = 25.6;FOVy = 19.2;dx = .4; dy = 2.4;
[G1,G2,G3,dt,Npts,Gx,Gy,kx,ky] = makeEPI(FOVx,FOVy,dx,dy);
b0 = 3*10000; %b0 = 3T * 10,000G/T
%Loop through time and simulate diffusion
%for t=1:size(Gx,1)
%    nextCell = diffuseCells(p_c,p,Gx,Gy,b0,dt,t);
%    p_c = nextCell;
%end
%%
%troubleshooting area
%[dx, dy] = randomwalk(0)
nextCell = diffuseCells(p_c,p,Gx,Gy,b0,dt,1);
%[dx,dy] = randomwalk(p(i,j)); 
%%
%Inputs: cell array of protons, phantom, Gx(t),Gy(t), b0, dt,t
%Output: cell array of protons for t+dt
function t_cell = diffuseCells(p_cell,p, G1, G2,b0,dt,t)
%Calculate gradients 
Gxc = G1(t)*linspace(round(-1*size(p_cell,2)/2), round(size(p_cell,2)/2), size(p_cell,2));
Gyc = G2(t)*linspace(round(-1*size(p_cell,1)/2), round(size(p_cell,1)/2), size(p_cell,1));

%Initialize temp cell storage for randomwalk calc
t_cell = cell(size(p_cell));  % Initialize 
for i = 1:numel(p_cell)
    t_cell{i} = p_cell{i};  % Copy each element 
end
%t_cell = p_cell;
%at each timepoint iterate through cells and update spins and phases,
%then update position according to random walk
for i=1:size(p_cell,2) % x coord
    for j=size(p_cell,1) % y coord
        a = cell2mat(p_cell(i,j));
        %update phase1 = mod(1,phase0 + spin (Hz) * dt (uS)),1), spin=4257 Hz/G * (b0 + Gx +Gy) (G) =
        %gamma*field strength
        a(1,:) = 4257*(b0+Gxc(i)+Gyc(j));
        a(2,:) = mod(a(2,:)+a(1,:)*dt/(10^6),1);
        
        %Storage for a after random walk
        a_r=a;
        %update position according to random walk
        for k=1:size(a,2)
            current = a(:,k);
            [dx,dy] = randomwalk(p(i,j)); %use diffusion coefficients stored in phantom as n for random walk
            if(p(i+dx,j+dy) == 0) %If moving out of object, stop
                dx=0;
                dy=0;
            end
            
            %move proton if it needs to be moved
            if(dx ~= 0 || dy ~= 0)
                %place in new cell
                %t_cell{i+dx,j+dy} = [];
                if(i+dx <= size(p_cell,2) && i+dx >= 1 && y+dy <= size(p_cell,1) && y+dy >=1)
                    t_cell{i+dx,j+dy} = [cell2mat(p_cell(i+dx,j+dy)),current];
                end
                %remove from old cell
                a_r(:,k) = [];
            end
        end
        %update cell with new proton entries
        %t_cell{i,j} = [];
        t_cell{i,j} = a_r;
    end
end
end
function [G1,G2,G3,dt,Npts,gx,gy,kx,ky] = makeEPI(FOVx,FOVy,dx,dy)
% Inputs
% FOVx = x-axis field of view in cm
% FOVy = y-axis field of view in cm
% dx = x-axis resolution in cm
% dy = y-axis resolution in cm
% 
%
% Outputs (note: all amplitudes are positive)
% G1:  G1 gradient amplitude in G/cm;
% G2:  G2 gradient amplitude in G/cm;
% G3:  G3 gradient amplitude in G/cm;
% dt:  readout ADC sampling rate in microseconds;
% Npts: number of points in the gradient vector (sampled at dTS = 10 usec); 
% gx: Npts x 1 column vector; x-gradient sampled at dTS; 
% gy: Npts x 1 column vector; y-gradient sampled at dTS; 
% kx = (Npts+1) x 1  column vector with kx values; sampled at dTS; 
% ky = (Npts+1) x 1  column vector with ky values; sampled at dTS; 
%
gam = 4257; %gamma/(2*pi) in Hz/G
G1dephase = 320;  % width of readout dephaser in microseconds.
G1readout = 640;  % width of each readout gradient in microseconds.
G2time = 10; % width of each phase encode blip in microseconds.
G3time = 320; % width of phase encode dephaser in microseconds.
dTS = G2time;  % sampling rate for waveforms; 
% Inputs
% FOVx = x-axis field of view in cm
% FOVy = y-axis field of view in cm
% dx = x-axis resolution in cm
% dy = y-axis resolution in cm
% 
%
% Outputs (note: all amplitudes are positive)
% G1:  G1 gradient amplitude in G/cm;
% G2:  G2 gradient amplitude in G/cm;
% G3:  G3 gradient amplitude in G/cm;
% dt:  readout ADC sampling rate in microseconds;
% Npts: number of points in the gradient vector (sampled at dTS = 10 usec); 
% gx: Npts x 1 column vector; x-gradient sampled at dTS; 
% gy: Npts x 1 column vector; y-gradient sampled at dTS; 
% kx = (Npts+1) x 1  column vector with kx values; sampled at dTS; 
% ky = (Npts+1) x 1  column vector with ky values; sampled at dTS; 
%
% Notes:
% 1) assume the gradients begin with t = 0 at the start of the dephasers.
% 2) The first element of each gradient vector describes the gradient for time interval [0 dTS]
% 3) the first element of each k-space vector should be zero -- i.e. start at
% the center of k-space.  
% 4) You may assume that the gradients can change instantaneously.  
gam = 4257; %gamma/(2*pi) in Hz/G
G1dephase = 320;  % width of readout dephaser in microseconds.
G1readout = 640;  % width of each readout gradient in microseconds.
G2time = 10; % width of each phase encode blip in microseconds.
G3time = 320; % width of phase encode dephaser in microseconds.
dTS = G2time;  % sampling rate for waveforms; 
% PUT YOUR CODE HERE; replace NaN  or 1 values with correct values.
% Calculating gradient amplitudes (in G/cm)
G1 = 1/(dx*gam*G1readout*10^-6); % Gx = 1/(FOVx*gam*G1readout*10^-6)
G2 = 1/(FOVy*gam*G2time*10^-6); % G2: phase encode gradient amplitude, typically for phase encoding (y-direction)
G3 = 1/(2*dy*gam*G3time*10^-6);
%1/(dy*gam*G3time*10^-6); %depends on dy
%dy/(gam*G3time*10^-6);
%gam*2*pi/(FOVy*dy);
%(pi/2)/(gam*G3time);
%(pi/2)/(gam*dy*2*pi*G3time*10^-6);
%(pi/2)/(gam*dy*2*pi*G3time*10^-6); %90 degree shift
%1/(FOVy*gam*G3time*10^-6); % G3: phase encode dephaser, how to make? shift to pi/2?
% Time intervals and Npts calculation
%wx = 1/dx;
%wy = 1/dy;
%dt = G2time;
dt = 1/(gam*G1*FOVx)*10^6;
%1/(2*max(wx,wy))*10^6; % us
%Npts = round((G1dephase+(FOVy/dy)*(G1readout+G2time-1))/dt);
Nx = FOVx/dx;
Ny = FOVy/dy;
Npts = round((G1dephase+Ny*G1readout+(Ny-1)*G2time)/dTS);
%FOVx/dx
%(FOVx/dx)*(FOVy/dy); 
%Npts = (4*(G1dephase+G1readout)+G1dephase)/dTs


gx = zeros(Npts, 1);  % x-gradient
gy = zeros(Npts, 1);  % y-gradient
kx = zeros(Npts + 1, 1);  % k-space values for x-direction
ky = zeros(Npts + 1, 1);  % k-space values for y-direction



Nd = round(G1dephase/dTS); %dephase
Nr = round(G1readout/dTS); %readout
Nc = Nd+1; %current
Nb = round(G2time/dTS); %blip

%Initialize dephase
gx(1:Nd) = -G1;
gy(1:round(G3time/dTS))=-G3;

for i=1:Ny
    j = mod(i,2);
    if(j==1)
        gx(Nc:Nc+Nr-1)=G1; 
    else
        gx(Nc:Nc+Nr-1)=-G1;
    end
    
    if(i<Ny)
        gy(Nc+Nr:Nc+Nr+Nb-1)=G2;
    end
    
    Nc = Nc + Nr + Nb;
end
gx(G1dephase/dt+1)=G1;
% Calculate k-space vectors (kx, ky) based on the gradient waveforms
for i = 2:Npts + 1
    kx(i) = kx(i-1) + gam * gx(i-1) * (dt * 1e-6);  % kx calculation in Hz/cm
    ky(i) = ky(i-1) + gam * gy(i-1) * (dt * 1e-6);  % ky calculation in Hz/cm
end
figure; clf;
nexttile;
dTS = 10;
totalTime = Npts*dTS
taxis = [0:dTS:totalTime];
plot(taxis(2:end),gx);
title("Gx vs Time")
xlabel("Time (uS)")
ylabel("Gradient amplitude (G/cm)")
legend("Gx")
nexttile;grid;
plot(taxis(2:end),gy);
title("Gy vs Time")
xlabel("Time (uS)")
ylabel("Gradient amplitude (G/cm)")
legend("Gx","Gy")
nexttile;grid;
plot(kx,ky,'-o');
title("Kx vs Ky")
xlabel("Kx (1/cm)")
ylabel("Ky (1/cm)")
grid;
end
%Inputs: phantom with pixel values as diffusion coefficients (must be
%integers)
%Output: cell array with 3 protons per pixel
function p_cell = makeCells(p)
p_cell = cell(size(p));
for i=1:size(p,1)
    for j=1:size(p,2)
        %initialize spins as 0, phases as 0 
        p_cell{i,j} = [1,0,3;0,0,0];
    end
end
end
%Inputs: # of steps to take
%Output: change in x and change in y position
function [dx, dy] = randomwalk(n)
dx = 0;
dy = 0;
for i=1:n
    %Uniform random variate [0 1]
    var = rand();
    %choose one of 8 directions randomly
    if(0 <= var && var <= 1/8)
        %d='UL';
        dx = dx-1;
        dy = dy+1;
    elseif(1/8 < var && var <= 1/4)
        %d='U';
        dy = dy+1;
    elseif(1/4 < var && var <= 3/8)
        %d='UR';
        dy = dy+1;
        dx = dx+1;
    elseif(3/8 < var && var <= 1/2)
        %d='R';
        dx = dx+1;
    elseif(1/2 < var && var <= 5/8)
        %d='DR';
        dy = dy-1;
        dx = dx+1;
    elseif(5/8 < var && var <= 3/4)
        %d='D';
        dy = dy-1;
    elseif(3/4< var && var <= 7/8)
        %d='DL';
        dy = dy-1;
        dx = dx-1;
    else
        %d='UL';
        dy = dy+1;
        dx = dx-1;
    end
end
end

Answer 1

[Stephen23]

This line:

for j=size(p_cell,1) % y coord

should be

for j=1:size(p_cell,1) % y coord
%     ^^

After that you will need to fix some y's which should be j's. And some indexing.

[Skyler]

Thanks! saved me a lot of trouble.

[Stephen23]

@Skyler: if it helped you then please remember to click the accept button!