Values not entering IF loop for unknown reason!

Hi I am really struggling with the following, I have code that loops over a time interval, within this I have an if statement from which a function file is ran. From my values I know this should work but it seems nothing is entering the if statement? IF anyone could help on the following it would be much appreciated!!

CODE

time=0; %Initialise time

tfinal=40; 

diagstep=10; 

diagcounter=0;

while time<tfinal

     if ((round(time/diagstep) * diagstep) == time) 

         [xcounter] = ParticleCounterFun(x0, particleno)   

         diagcounter=diagcounter+1
     end

%Where the function file ParticleCounterFun is;

function [ xcounter ] = ParticleCounterFun( x0, particleno )

%Initialise counters used in E field diagnostics
counter01=0;

counter12=0;

counter23=0;

counter34=0;

counter45=0;

    for np=1:particleno
         if (0<=x0(np)) & (x0(np)<1)
             counter01=counter01+1;
         end
         if (1<=x0(np)) & (x0(np)<2)
             counter12=counter12+1;
         end
         if (2<=x0(np)) & (x0(np)<3)
             counter23=counter23+1;
         end
         if (3<=x0(np)) & (x0(np)<4)
             counter34=counter34+1;
         end
         if (4<=x0(np)) & (x0(np)<5)
             counter45=counter45+1;
         end

     end
     xcounter=[counter01 counter12 counter23 counter34 counter45]

end

I really need this solved as I am unable to continue with work on my project so any help would be greatly appreciated! Thanks in advance

[the cyclist]

I can't run your code because you don't indicate what x0 and particleno are.

[Mischa Kim]

Seems to be a duplicate of http://www.mathworks.de/matlabcentral/answers/117142-why-do-i-get-error-undefined-function-or-variable.

[Peter Nave]

time does not seem to change inside the loop.

[Phoebe]

this is an extract the full code is as follows;

%Physical constants

q=-1; %Charge of particle

m=1; %Mass of particle

xstart=0;

ystart=0;

dx=1;

dy=1;

particleno=4;

nrows=sqrt(particleno);

ncolumns=nrows;

nrows=nrows+1;

B=[0 0 8]; %Magnetic field strength (T)

time=0; %Initialise time

dt=0.01; %Define time step

h=dt/2;

tfinal=40; %State time for code to stop

v0=[2 2 0];

%Constants used in numerical procedure

b=norm(B,2); %Takes absolute value of vector B

t=(((q*b)/m)*h);

s=(2*t)/(1+(t^2));

con=q*h/m;

%Storage of three consecutive positions

x3=0; x2=0; x1=0;

y3=0; y2=0; y1=0;

%Particle Initialisation

yposition=ystart;

xposition=xstart;

np=1;

for ny=1:nrows

    for nx = 1:ncolumns
     x0(np)=xstart+(nx*dx);
     y0(np)=yposition;
     np=np+1;
     if (np>particleno) break
     end
    end
     if (np>particleno) break
     end
    yposition=ystart+(ny*dy);

end for np = 1 : particleno vx0(np)=v0(1); vy0(np)=v0(2); x3(np)=0; x2(np)=0; x1(np)=0; y3(np)=0; y2(np)=0; y1(np)=0; end

My=max(y0) + 10; %Maximum value of y in Box

my=-My; %Minimum value of y in Box

Mx=max(x0) + 10; %Maximum value of x in Box

mx=-Mx;%Minimum value of x in Box

diagstep=10; %define step size for E field diagnostic to count number of particles at that time

diagcounter=0;

while time<tfinal for np = 1 : particleno %Run Function to calculate new velocity and position v0=[vx0(np) vy0(np) 0]; My=5; Ey=2*cos((2*pi)*(y0(np)/My)); E=[0 Ey 0 ]; [v1]=particlepusherfun( v0, con, E, t, s);

    x=(dt*v1(1))+x0(np); %New position x component
    y=(dt*v1(2))+y0(np); %New position y component

    %Collect and store 3 consecutive x and y values of new position
    x3(np)=x2(np);  x2(np)=x1(np);  x1(np)=x;
    y3(np)=y2(np);  y2(np)=y1(np);  y1(np)=y;

   %Reset position, velocity and time for loop
   x0(np)=x1(np);
   y0(np)=y1(np);
   vx0(np)=v1(1);
   vy0(np)=v1(2);
   time=time+dt;  

   %Run Function to calculate guiding centre for diagnostic test
   [a, b]=guidingcentrefun(x1(np), x2(np), x3(np), y1(np), y2(np), y3(np));

    hleg1=plot(x,y,'m'); %plot x and y values
   hold on; 
   if (x3(np)~= 0)
      hleg2=plot(a,b,'b'); %plot guiding centres
   end

   end %End np=1:particle LOOP

   if ((round(time/diagstep) * diagstep) == time) %If the integer value of time/diagstep is = that time count it 
       [xcounter] = ParticleCounterFun(x0, particleno)   
       diagcounter=diagcounter+1
   end

 %Xcounter(diagcounter)=xcounter;
 %xposition=xposition+1;
end %End while time<tfinal LOOP

%Annotate graphs title('Motion of a particle in an Electric & Magnetic Field'); xlabel('x distance (m)','FontSize',16) ylabel('y distance (m)','FontSize',16) hleg3=legend([hleg1 hleg2],'Particle Positions','Guiding Centre'); set(hleg3,'Location','NorthEastOutside'); %Output final coordinates of guiding centre C=[a b 0]

FUNCTION FILES;

function [ v1] = pariclepusherfun( v0, con, E, t, s)

 %Add half the electric impulse to pre-rotational velocity vector

    vpre=v0+(con*E);
    vxpre=vpre(1); %x component
    vypre=vpre(2); %y component

    %Rotate velocity components (magnetic field contribution)
    vdash=vxpre + (vypre*t);
    vypost=vypre - (vdash*s);
    vxpost=vdash + (vypost*t);
    vpost=[vxpost vypost 0];

    %Add remaining half of electric impulse to post-rotational velocity
    v1=vpost+(con*E); 
end

function [ a b ] = guidingcentrefun(x1, x2, x3, y1, y2, y3)

%Calculate guiding centre with simultaneous eqns

k1=(-2*x1)+(2*x2);
k2=(-2*y1)+(2*y2);
k3=(x2^2) -(x1^2) + (y2^2) - (y1^2);
k4=(-2*x1)+(2*x3);
k5=(-2*y1)+(2*y3);
k6=(x3^2)-(x1^2)+(y3^2)-(y1^2);
b=((k1*k6)-(k3*k4))/((k1*k5)-(k2*k4));
a=(k3-(k2*b))/k1;
end

function [ xcounter ] = ParticleCounterFun( x0, particleno )

%Initialise counters used in E field diagnostics

counter01=0;

counter12=0;

counter23=0;

counter34=0;

counter45=0;

for np=1:particleno
     if ((0<=x0(np)) & (x0(np)<1))
         counter01=counter01+1;
     end
     if ((1<=x0(np)) & (x0(np)<2))
         counter12=counter12+1;
     end
     if ((2<=x0(np)) & (x0(np)<3))
         counter23=counter23+1;
     end
     if ((3<=x0(np)) & (x0(np)<4))
         counter34=counter34+1;
     end
     if ((4<=x0(np)) & (x0(np)<5))
         counter45=counter45+1;
     end

   end

   xcounter=[counter01 counter12 counter23 counter34 counter45]

end

[Image Analyst]

Phoebe: Please read this: http://www.mathworks.com/matlabcentral/answers/13205-tutorial-how-to-format-your-question-with-markup. Then delete all that in your prior comment and attach your m-file with the paper clip icon. Then read this: http://blogs.mathworks.com/videos/2012/07/03/debugging-in-matlab/

[Star Strider]

My guess is that you’re yet another victim of rounding error.

Just after

while time < tfinal

insert:

test = (round(time/diagstep) * diagstep) - time

My guess is that it will be non-zero, so:

(round(time/diagstep) * diagstep) == time

will never evaluate to ‘true’.