function stimuli = dots(trial,block,trialMatrix,is_standard,Gparams)

NSEQ        = Gparams.dots.nseq; % number of interleaved sequences
NDOTS       = Gparams.dots.nDots;

if is_standard
    COH = Gparams.dots.standard;
else
    COH = Gparams.dots.standard + trialMatrix.relative_coherence(block,trial);
end

DIRECTION   = trialMatrix.direction(block,trial)*pi/180 - 90*pi/180 ;  % convert direction from degrees to radiants, where 0 degrees is upwards, 90 degrees is rightwards, etc

SPEED       = Gparams.dots.speed;
NFR         = ceil(Gparams.timing.nframes/NSEQ); 
RADIUS      = Gparams.dots.radius;

for var = 1:NSEQ
    
    % generate random starting points within a circular aperture
    rad                     = RADIUS*sqrt(rand(NDOTS,1));
    theta                   = 2*pi*rand(NDOTS,1);
    [pos(:,1), pos(:,2)]    = pol2cart(theta, rad); %convert to cartesian coordinates
    
    for frameNum = 1:NFR
        
        % define the noise dots
        noisedots   = round(NDOTS*(1-abs(COH))); %nr of noisedots based on the coherence
        tmp         = randperm(NDOTS);
        noiseindex  = tmp(1:noisedots); %random subset of dots
        
        % define signal dots
        signalindex = tmp(noisedots+1:end); %the dots that are signal
        
        % move signal dots with a certain speed in the right direction
        pos(signalindex,:)  = [SPEED.*cos(DIRECTION)/Gparams.frameRate+pos(signalindex,1), ...
            SPEED.*sin(DIRECTION)/Gparams.frameRate+pos(signalindex,2)]; %convert to cartesian coordinates
        
        % replot the noisedots somewhere in the aperture
        rad                 = RADIUS*sqrt(rand(noisedots,1));
        theta               = 2*pi*rand(noisedots,1);
        [pos(noiseindex, 1), pos(noiseindex,2)] = pol2cart(theta, rad);
        
        
        %find the dots that have left the aperture
        outindex            = find(sqrt((pos(:,1).^2)+(pos(:,2).^2)) >= RADIUS); %index dots that are outside the aperture
        
        % wrap them around in the direction where they came from
        [theta, rad]        = cart2pol(pos(outindex, 1), pos(outindex,2));
        theta               = theta + pi; %move to the other side of the circle
        rad                 = RADIUS*ones(length(rad),1);
        [pos(outindex, 1), pos(outindex,2)] = pol2cart(theta, rad);
        
        
        % save the positions per variant
        temp(var,frameNum, :, :)  = pos';
        
    end   
    
end

stimuli = squeeze(reshape(temp, 1, NFR*NSEQ, 2, NDOTS));
stimuli = stimuli(1:Gparams.timing.nframes, :, :);

end

