function  RDC_spec = make_RDC_spec(Rparams,th,N_coher,N_total)

%MAKE_RDC_SPEC -- Frame-by-frame X-Y spec of a random-dot cinematogram
%   
%  RDC_spec = make_RDC_spec(Rparams)
%  RDC_spec = make_RDC_spec(Rparams,th,N_coher,N_total)
%
%  Rparams is a structure produced by RDC_PARAMS.
%  TH specifies the direction of coherent motion [deg. from vertical].
%  N_COHER specifies the number of coherently moving dots.
%  N_TOTAL specifies the total number of dots within the aperture.
%  If not supplied, TH, N_COHER, and N_TOTAL default to 'direction',
%  'N_coher', and 'N_dots' parameter fields, respectively.
%
%  The result RDC_SPEC is an Nx3xT matrix with tuples (x,y,age)
%  for each of N dots across T frames. The coherent dots are on top.
%  Use RENDER_RDC or RENDER_RDC_DELUXE to render.
%
%  Algorithm (Britten et al, 1992):
%    Dots are generated in a larger square and then clipped within a
%    circular aperture with diameter Rparams.aperture.  Each frame is
%    guaranteed to contain exactly N_total dots within the aperture.
%    On each successive frame, N_coher dots are sampled at random and
%    updated with a global motion vector relative to their respective
%    locations on the previous frame. The remaining (N_total-N_coher)
%    dots are replaced with newly generated dots at random locations.
%    Dots that move outside the aperture are dealt with depending on
%    the Rparams.wraparound parameter:
%     * wraparound=0: Just compensate attrition with new random dots
%     * wraparound=1=default:  Identify the dots that would go out of
%          bounds if moved. Instead of moving them, replace them with
%          dots diametrically opposed w.r.t their location in the circle.
%     * wraparound=2: Move the dots first, then identify those who had
%          actually gone out, then flip them w.r.t their new location.
%          This mode is flawed, particularly at high motion speeds.
%          It thends to concentrate the dots near the midline of motion.
%
%  Example:
%    Rparams = RDC_params ;
%    RDC_spec = make_RDC_spec(Rparams) ; % 100 incoherent dots
%    M = RDC_movie(Rparams,RDC_spec) ;   % diagnostic plot
%    movie(M,3)                          % replay movie 3 times
%
%  Example:
%    RDC_movie(Rparams,make_RDC_spec(Rparams,-90,20))    % 20% coh leftward
%    RDC_movie(Rparams,make_RDC_spec(Rparams,-90,20,20)) % 100% coherence
%
%  Typical usage:
%    RDC_spec = make_RDC_spec(Rparams,45,50) ;   % 50% coh diagonal
%    RDC = render_RDC(Rparams,RDC_spec) ;
%    M = RDC_movie(Rparams,RDC) ;
% 
%  Example -- Shows the flaw in wraparound=2
%    Rparams.wraparound=2; Rparams.dot_speed=50; Rparams.duration=150;
%    Rparams = RDC_params(Rparams)     % update all derivative fields
%    RDC_spec = make_RDC_spec(Rparams,45,100) ;
%    M = RDC_movie(Rparams,RDC_spec) ;
%
%  Reference:
%    Britten, K., Shadlen, M., Newsome, W., & Movshon, J. (1992).
%    The analysis of visual motion: A comparison of neuronal and
%    psychophysical performance. _J of Neuroscience_, 12 (12), 4745-4765.
%
% See also RDC_PARAMS, RENDER_RDC, RDC_MOVIE, RENDER_RDC_DELUXE.

% Original coding by Alex Petrov, Ohio State University
% $Revision: 2.0 $  $Date: 2007-12-05 $
%
% Part of the utils toolbox version 1.2 for MATLAB version 5 and up.
% http://alexpetrov.com/softw/utils/
% Copyright (c) Alexander Petrov 1999-2008, http://alexpetrov.com
% Please read the LICENSE and NO WARRANTY statement in ../utils_license.m

% 2.0     2007-12-05  ap -- Switch to the Newsome algorithm
% 1.0     2006-06-04  ap -- Individ.dots moved coherently for many frames


%%-- Handle input args and decipher params
if (nargin<4) ; N_total = Rparams.N_dots ;  end
if (nargin<3) ; N_coher = Rparams.N_coher ;  end
if (nargin<2) ; th = Rparams.direction ; end

T = Rparams.duration ;         % total movie length, in frames
Rmax = Rparams.aperture/2 ;
Rmax2 = Rmax^2 ;
dr = Rparams.displacement ;    % travel distance per frame, deg.vis.ang
th1 = th.*(pi/180) ;           % motion direction in radians
dxy = repmat(dr.*[sin(th1) cos(th1)],N_coher,1) ;    % 0==upward
redundancy = 1.4 ;             % 1.1 * 4/pi
wraparound = Rparams.wraparound ;

RDC_spec = zeros(N_total,3,T) ;
if (N_total==0) ; return ; end   % return an empty matrix

%%-- Main loop
for t = 1:T
    %- Coherent dots generated first
    if ((t==1) | (N_coher<=0))   % | used instead of || for OS9 compat.
        % Fully random field
        n2go = N_total ;
        beg = 1 ;              % random dots begin at index 1
    else  %% Move N_COHER dots from the previous frame
        % A freshly sampled set of dots move on every frame
        if (N_coher==N_total)
            dots2move = RDC_spec(:,:,t-1) ;
        else
            idx = randperm(N_total) ;
            dots2move = RDC_spec(idx(1:N_coher),:,t-1) ;
        end
        % Apply a coherent motion vector
        xy = dots2move(:,1:2) + dxy ;
        % Deal with the dots that moved outside the aperture
        R2 = sum(xy.^2,2) ;
        if (wraparound==0)
            % Just compensate attrition with fresh random dots
            idx = find(R2 < Rmax2) ; % keep only those inside
            n = length(idx) ;
            RDC_spec(1:n,1:2,t) = xy(idx,:) ;
            RDC_spec(1:n,3,t) = dots2move(idx,3) + 1 ;
        else  % Systematically flip the dots that exited on the far side
            idx = find(R2 > Rmax2) ;  % flip those outside
            if (wraparound==1)
                % Flip the dots from their original location
                xy(idx,:) = -dots2move(idx,1:2) ;
            else  % (wraparound==2)
                % Move and then flip
                R = sqrt(R2(idx)) ;
                wrapping_factor = (R-2*Rmax)./R ;    % b/n -1 and -.9
                xy(idx,1) = xy(idx,1) .* wrapping_factor ;
                xy(idx,2) = xy(idx,2) .* wrapping_factor ;
            end
            dots2move(idx,3) = -1 ;  % incremented below
            RDC_spec(1:N_coher,1:2,t) = xy ;
            RDC_spec(1:N_coher,3,t) = dots2move(:,3) + 1 ;
            n = N_coher ;
        end
        % N2GO dots remain to be generated from scratch at random
        n2go = N_total - n ;
        beg = n + 1 ;
    end
    %- Fill the remainder with N2GO randomly generated dots
    while (n2go>0)
        xy = (2*Rmax).*rand(ceil(redundancy*n2go),2) - Rmax ;
        idx = find(sum(xy.^2,2) < Rmax2) ;
        n = min(length(idx),N_total-beg+1) ;
        RDC_spec(beg:beg+n-1,1:2,t) = xy(idx(1:n),:) ;
        % RDC_spec(_,3,t) is initialized to 0
        beg = beg + n ;
        n2go = n2go - n ;
    end
end


%--- Return RDC_SPEC
%%%%% End of file MAKE_RDC_SPEC.M
