excel - 将蒙特卡罗计算从 VBA 代码更改为 MATLAB

Question

我正在尝试将此子句因子蒙特卡罗计算代码（用 Visual Basic 编写）：https://descanso.jpl.nasa.gov/SciTechBook/series1/Goebel_AppG_Clausing.pdf 转换为 MATLAB 代码。

当我尝试在 Excel 中运行此代码时，它总是在第一次运行时返回一致的值。但是当我在 MATLAB 中实现时，它会随机返回子句因子。您能否帮我编写 MATLAB 代码，以便它可以像在 VBA 中那样做。

我在 VBA 代码中谈论的值是Range("C11") = (rBottom ^ 2) * iescape / npart

这是我解决 clausingFactor 的 MATLAB 代码

    thickScreen = 0.5;
    thickAccel = 1;
    rScreen = 0.8;
    rAccel = 0.5;
    gridSpace = 1;
    npart = 313;
    Pi = 3.14159265358979;

%assumes rTop = 1
    rBottom = rScreen/rAccel;
    lenBottom = (thickScreen + gridSpace)/rAccel;
    lenTop = thickAccel / rAccel;
    Length = lenTop + lenBottom;  
    iescape = 0;
    maxcount = 0;
    icount = 0;
    nlost = 0;
    vztot = 0;
    vz0tot = 0;
    for ipart = 1: npart
%launch form bottom
        notgone = true;
        r0 = rBottom * sqrt(rand());
        z0 = 0;
        costheta = sqrt(1 - rand());
        if (costheta > 0.99999)
            costheta = 0.99999;
        end
        Phi = 2 * Pi * rand();
        sintheta = sqrt(1 - costheta);
        vx = cos(Phi) * sintheta;
        vy = sin(Phi) * sintheta;
        vz = costheta;
        rf = rBottom;
        t = (vx * r0 + sqrt((vx^2 + vy^2) * rf^2 - (vy * r0)^2)) / (vx^2 + vy^2);
        z = z0 + vz * t;
        vz0tot = vz0tot + vz;
        icount = 0;
        while notgone
            icount = icount +1;
            if (z < lenBottom)
%hit wall of bottom cylinder and is re-emitted
                r0 = rBottom;
                z0 = z;
                costheta = sqrt(1-rand());
                if (costheta > 0.99999)
                    costheta = 0.99999;
                end
                Phi = 2 * Pi * rand();
                sintheta = sqrt(1 - costheta^2);
                vz = cos(Phi) * sintheta;
                vy = sin(Phi) * sintheta;
                vx = costheta;
                rf = rBottom;
                t = (vx * r0 + sqrt((vx^2 + vy^2) * rf^2 - (vy * r0)^2)) / (vx^2 + vy^2);
                z = z0 + t * vz;
            end %bottom cylinder re-emission
            if ((z >= lenBottom) && (z0 < lenBottom))
                %emitted below but going up
                %find radius at lenBottom
                t = (lenBottom - z0) / vz;
                r = sqrt((r0 - vx * t)^2 + (vy * t)^2);
                if (r <= 1)
                    %continuing upward
                    rf = 1;
                    t = (vx * r0 + sqrt((vx^2 + vy^2) * rf^2 - (vy * r0)^2)) / (vx^2 + vy^2);
                    z = z0 + vz*t;
                else
                    % hit the upstram side of the accel grid and is
                    % re-emitted downward
                    r0 = r;
                    z0 = lenBottom;
                    costheta = sqrt(1 - rand());
                    if (costheta > 0.99999)
                        costheta = 0.99999;
                    end
                    Phi = 2 * Pi * rand();
                    sintheta = sqrt(1 - costheta^2);
                    vz = cos(Phi) * sintheta;
                    vy = sin(Phi) * sintheta;
                    vx = costheta;
                    rf = rBottom;
                    t = (vx * r0 + sqrt((vx^2 + vy^2) * rf^2 - (vy * r0)^2)) / (vx^2 + vy^2);
                    z = z0 + t * vz;
                end
            end %end upward
            if ((z >= lenBottom) && (z<= Length))
                %hit the upper cylinder wall and is re-emitted
                r0 = 1;
                z0 = z;
                costheta = sqrt(1 - rand());
                if (costheta > 0.99999)
                    costheta = 0.99999;
                end
                Phi = 2 * Pi * rand();
                sintheta = sqrt(1 - costheta^2);
                vz = cos(Phi) * sintheta;
                vy = sin(Phi) * sintheta;
                vx = costheta;
                rf = 1;
                t = (vx * r0 + sqrt((vx^2 + vy^2) * rf^2 - (vy * r0)^2)) / (vx^2 + vy^2);
                z = z0 + t * vz;
                if (z < lenBottom)
                    %find z when particle hits the bottom cylinder
                    rf = rBottom;
                    if ((vx^2 + vy^2) * rf^2 - (vy * r0)^2 < 0 )
                        t = (vx * r0) / (vx^2 + vy^2);
                        %if sqrt argument is less than 0 then set sqr term
                        %to 0 12 May 2004
                    else
                        t = (vx * r0 + sqrt((vx^2 + vy^2) * rf^2 - (vy * r0)^2)) / (vx^2 + vy^2);
                    end
                    z = z0 + vz * t;
                end
            end %end upper cylinder emission
            if (z < 0)
                notgone = false;
            end
            if (z > Length)
                iescape = iescape + 1;
                vztot = vztot + vz;
                notgone = false;
            end
            if (icount > 1000)
                notgone = false;
                icount = 0;
                nlost = nlost + 1;
            end
        end
        if (maxcount < icount)
            maxcount = icount;
        end
    end
    clausingFactor = (rBottom^2) * iescape / npart;
    vz0av = vz0tot / npart;
    vzav = vztot / iescape;
    DenCor = vz0av / vzav;

谢谢你。

标签： excelvbamatlabmontecarlo