function out = ArenaPathsSim(numSims,simSteps,headingDrift,bias,visualizationFlag)

    bodyLength = 2.5; % convert units to body length
    R=25.4/bodyLength;

    out.X=zeros(simSteps,numSims);
    out.Y=zeros(simSteps,numSims);
    out.inROI=zeros(simSteps,numSims);
    out.tmd=zeros(simSteps,numSims);

    speed=1/bodyLength;
    reflectEpsilon=0.05*2*pi;

    for i=1:numSims

        r=Inf;
        while r>R
            x=(rand()*2*R-R);
            y=(rand()*2*R-R);
            [~, r]=cart2pol(x,y);
        end

        heading=rand()*2*pi;

        for t=1:simSteps
            [th, ~]=cart2pol(x,y);
            th=mod(th,2*pi);
            deltaHeading=randn()*2*pi*headingDrift+bias;
            heading=mod(heading+deltaHeading,2*pi);

            newPosX=x+speed*cos(heading);
            newPosY=y+speed*sin(heading);
            [~, newr]=cart2pol(newPosX,newPosY);
            if newr>R %wall colision
                deltaThetas=linspace(-pi,pi,200);
                deltaThetaPDF=normpdf(deltaThetas,0,headingDrift);
                thetaBool=abs(angdiff(deltaThetas,th-heading))<pi/2;

                deltaThetaPDF(thetaBool)=0;

                if sum(deltaThetaPDF) > 10^-100;
                    deltaThetaCDF=cumsum(deltaThetaPDF);

                    deltaThetaIndex=sum(rand > [0 , deltaThetaCDF(1:end-1)]);
                    heading=heading+deltaThetas(deltaThetaIndex);
                elseif rand()<abs(angdiff(th+pi/2,heading))/pi
                    heading=th-(pi/2+reflectEpsilon);
                else
                     heading=th+(pi/2+reflectEpsilon);
                end
                x=x+speed*cos(heading);
                y=y+speed*sin(heading);
            else
                x=newPosX;
                y=newPosY;
            end

            out.X(t,i)=x;
            out.Y(t,i)=y;
            [th, r]=cart2pol(x,y);
            if r<0.8*R
                out.inROI(t,i)=1;
            end
            out.tmd(t,i)=angdiff(th-heading);

        end
    end

    if numSims==1
        binCenters=0.00001-pi:2*pi/25:2*pi-pi;
            % without masking paths close to edge
        h=hist(out.tmd,binCenters);
        h = h/sum(h);
        out.binCenters = binCenters;
        out.histogramFrequency = h;
        out.mu = sum(h.*sin(binCenters));
            % with masking paths close to edge
        out.tmdMasked = out.tmd.*out.inROI;
        out.tmdMasked = out.tmdMasked(out.tmdMasked~=0);
        hMasked=hist(out.tmdMasked,binCenters);
        hMasked = hMasked/sum(hMasked);
        out.histogramFrequencyMasked = hMasked;
        out.muMasked = sum(hMasked.*sin(binCenters));
        
        if visualizationFlag == 1
            close all
            figure
            hold on
            scatter(out.X(out.tmdMasked>0),out.Y(out.tmdMasked>0),'r.');
            scatter(out.X(out.tmdMasked<0),out.Y(out.tmdMasked<0),'b.');


            figure
            bar(out.binCenters,out.histogramFrequencyMasked)
        end
    end
end