Для выходных данных цикла одни и те же значения ev после каждой итерации

Question

это длинный сценарий, но я чувствую, что представление всей информации важно для того, чтобы понять, что происходит.Основная проблема возникает, когда я пытаюсь вычислить «ev», которая представляет собой энергию электрона в эВ, и «cs», которая представляет собой поперечную функцию площади, основанную на эВ, учитывая входные данные компонентов скорости.Я сделал это так, чтобы компоненты скорости произвольно менялись на каждом временном шаге, однако, независимо от того, какая скорость принята системой, матрицы ev и cs остаются одинаковыми.Я не могу понять, почему, потому что когда я тестировал каждую функцию без цикла for, они работают правильно и генерируют разные значения.

Функция "ev", на которую ссылается основной скрипт:

function [ev] = evfun(vx,vy,vz)

m_e = 9.11E-31; %Mass of electron
jtoev = 6.24E+18; %Joules to eV conversion

ev = jtoev*(0.5*(m_e)*(vx.^2+vy.^2+vz.^2));
%ev = log(nonfitev);

end

Функция "cs", на которую ссылается основной сценарий "

function cs = crosecarc(ev)

nonfitcs = 34*exp(-((ev-500)/350)^2); %In units of 1E-16 cm^2
cs = nonfitcs*(1E-4); %Converts to 10^-16 m^2

end

Основной сценарий:

filename = 'evtestsheet'; %Specify file name to print to
evmatsheet = 1; %Specify sheet in Excel File
evmatxlRange = 'G3'; %Specify evmat matrix start cell in Excel File
inputvelsheet = 1; %Specify sheet in Excel File
inputvelxlRange = 'D3'; %Specify inputvel matrix start cell in Excel File
inputpossheet = 1; %Specify sheet in Excel File
inputposxlRange = 'A3'; %Specify inputpos matrix start cell in Excel File

%Defining Space and Particle Density based on Pressure PV = nkT
mTorr = 500; %Define pressure in mTorr (From metter reading)
P = (mTorr/1000)*133.322; %mTorr to Pascal conversion to define pressure
k = 1.38E-23; %Boltzman Constant
T = 290; %Temperature in Kelvin. Converted from C  and @ room temp
N = 1; %Total number of particles
ChamTotalV = 1;


%Particle sweep volume is defined later because it depends on electron energy
tstep = 0.05E-8; %Defining time step
%modstep = 2*tstep; %Used when calculating volume swept since the DOE solver splits itself up into 3 times within each time step
tfin = 1.1E-8; %Defining final time
intspan = [0:tstep:tfin]; %Total time span
intspancol = intspan(:); %For Excel plotting purposes
[introw,intcol] = size(intspan);

%This section is just for generating matrices that will be populated later
tindex = [0:1:intcol-2];
tindexcol = tindex(:); %Purely for Excel plotting purposes
[tinrow,tincol] = size(tindex);
PCPmat = zeros(tincol,1);
choicemat = zeros(tincol,1);
colnewv = zeros(intcol,3);
Vsweepmat = zeros(intcol,1);
vmagmat = zeros(intcol,1);
evmat = zeros(intcol,1);
crossareamat = zeros(intcol,1);
inputpos = zeros(intcol,3);
inputvel = zeros(intcol,3);

%Test Initial Conditions
spart = [0.1 0.1 0.];
vpart = [-1.0585E+7 -1.0585E+7 0.5585E+7];

%for p = 1:1:size(s0parts) %particle count
    for t = 0:1:intcol-2; %complete time interval and time step

        [trow,tcol] = size(t);
        %spart = s0parts(p:p,1:3);
        %vpart = v0parts(p:p,1:3);

        %Defining relative position and velocity
        x = spart(1);
        y = spart(2);
        z = spart(3);
        vx = vpart(1);
        vy = vpart(2);
        vz = vpart(3);

        inputvel(t+1,1:3) = [vx vy vz];
        inputpos(t+1,1:3) = [x y z];

        r = sqrt(x.^2 + y.^2 + z.^2);
        vmag = sqrt(vx.^2 + vy.^2 + vz.^2);

        ev = evfun(vx,vy,vz); %Plugs in velocity components of electron into the "evfun" function defined in seperate file to find energy in eV
        cs = crosecarc(ev); %Plugs in energy of electron into the "crosecarc" function defined in seperate file to find crossectional area
        crossarea = cs*4;
        Vsweep = crossarea*vmag*tstep; %Volume of space swept out by particle Vsweep = (area cross section)*(vmag)*(time step)
        %nden = N/ChamTotalV; %Number of particles per unit volume
        nden = 1/Vsweep; %Test line to force collision at every time step
        PCP = nden*Vsweep; %Probability of collision
        choice = rand(1); %Generate random number "choice" between 0 and 1

        %Populate the previously generated matrices with current values
        evmat(t+1,1) = ev;
        crossareamat(t+1,1) = crossarea;
        vmagmat(t+1,1) = vmag;
        Vsweepmat(t+1,1) = Vsweep;
        PCPmat(t+1,1) = PCP;
        choicemat(t+1,1) = choice;

        if choice <= PCP %If the choice by system is accepted within probability amplitude, then collision occurs
            %Assign new, velocity components based on energy prior to collision
            vxrand = 2*rand()-1; %Random # generation for component x between -1 and 1
            vyrand = 2*rand()-1; %Random # generation for component y between -1 and 1
            vzrand = 2*rand()-1; %Random # generation for component z between -1 and 1

            vmagrand = sqrt(vxrand.^2 + vyrand.^2 + vzrand.^2); %Normalize the randomly generated components
            vx = (vxrand/vmagrand)*vmag; %Multiply normalized, random x component by original velocity magnitude
            vy = (vyrand/vmagrand)*vmag; %Multiply normalized, random y component by original velocity magnitude
            vz = (vzrand/vmagrand)*vmag; %Multiply normalized, random z component by original velocity magnitude
        else
            newv = [0 0 0];
        end

        %Redfine the velocity and position components to reference on next if-loop run
        spart(1) = x;
        spart(2) = y;
        spart(3) = z;
        vpart(1) = vx;
        vpart(2) = vy;
        vpart(3) = vz;
        ev = [];
         %Clears ev matrix
    end

    xlswrite(filename,inputvel,inputvelsheet,inputvelxlRange)
xlswrite(filename,inputpos,inputpossheet,inputposxlRange)
    xlswrite(filename,evmat,evmatsheet,evmatxlRange)