Values not entering IF loop for unknown reason!

Question

0 votes

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

6 commentaires
Afficher 4 commentaires plus anciens Masquer 4 commentaires plus anciens

Phoebe le 23 Fév 2014

Ouvrir dans MATLAB Online

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]