function [retval,iimax] = RadSim(A,B,C,E,X,Y,M,D,a,dr,T,dt,P,RadialDist)
% usage: [events,iimax] = RadSim(A,B,C,D,X,Y,M,D,a,dr,T,dt,P,RadialDist)
% description: returns events from a radial diffusion model simulation
% and optionally the number of space annuli used in the simulation.
% The model is a radially symmetric diffusion problem with 
% letf Dirichlet boundary condition. The PDE of the problem is
%  N_t = (1/r)*(D*r*N_r)_r + f(r,t,N), 0 < r < infinity
% The solution is Nsoln(r,t,M,D,a).  The physical meaning of Nsoln is
% physical "event", e.g., criminal activity, density. So we've hardwired
% the name 'Nsoln' into this routine, but we waffle on meaning of M, D, a.
% We're building in flexibility, but input parameters a,M,D have
% a specific use in the particular implementation of this file,
% where, in fact, we mean to have f(r,t,N) = a*N.
% Here are the meanings of the input parameters:
% A,B: left and right endpoints of the view window (x coords)
% C,E: bottom and top endpoints of the view window (y coords)
% X,Y: coordinates of the center of the diffusion process
% M: instantaneous source of "diffusive activity" placed at center at t=0.
% D: diffusion coefficient
% a: growth/decay term
% dr: size of concentric ring in which activity is placed
% T: number of time steps to use, starting at t=0, ending at t=T*dt
% dt: time increment
% P: probability that an event will be observed. (NB: it might be
%    more realistic to have P be a position dependent probability..later!
% RadialDist: we can distribute the data randomly in various ways.
% 
% There is the question of exactly what a discrete event really is.  Here
% is an idea: think of an event as occupying both space and time. The
% number of events occurring in a area A time time interval dt is the
% volume of the solid with appropriate area times appropriate time times
% the value of the density function at some point in this pseudo-3D box.
% ("Pseudo" because coordinates t and r determine such a point.)
%   For now, we just do this:
%  Parital Circle(0): Events are tallied in annulus of thickness drxdt
%                     at each time step to border of window. If less than one,
%                     skip; else, divide annulus into equal area sectors
%                     of event mass one. Determine if each is observed by
%                     probability P, and if so, randomly place the event
%                     in the appropriate sector.
%  Partial Circle(1): Same as above, but event placed at center of sector.
%                     Note, one can place events randomly in a sector or
%                     regularly.  See rem in switch case 1 below.
%     Full Circle(2): Events are tallied as above, then they are counted
%                     one by one and tested for occurrence on a random
%                     location on entire annulus with given probablilty P.
% In all of above cases, if event location is not in view window, it
% simply isn't shown.  So we ignore barriers of city boundaries, etc.
% Further, we prevent loss of outlying events by the following strategy:
% if the quantity of events in the annulus is less than one, preserve the
% inner radius of this annulus and pass the fraction of an event therein 
% on to the next annulus. When one event is reached, use a random location
% between the initial inner radius and last outer radius.
% The returned value is a Events X 3 matrix whose rows represent an event 
% at the space-time point (x,y,j*dt).  We set a maximum number of events in
% the interest of efficiency, and we might toss this limitation out later.
% 

% Last modified: 7/16/05

% local variables:
% MaxEvents: upper bound on number of events recorded
% radmax: radius enclosing view window
% iimax: number of annuli to loop over
% ndx: index into returned array of events
% eventrem: remainder of event quantity in previous annulus
% ii,jj: index variables
% t: current time
% radbase: inner radius for determination of event location


% initialization
MaxEvents = 5000000;
retval = zeros(MaxEvents,3);
% minimal amount of bulletproofing
if ((X<A) || (X>B) || (Y<C) || (Y>E)) % bad center?
    retval = 0;
    return
end
% find max radius
radmax = max([norm([X,Y]-[A,C]),norm([X,Y]-[A,E]), ...
        norm([X,Y]-[B,C]),norm([X,Y]-[B,E])]);
iimax = ceil(radmax/dr);
ndx = 1; 
% main loop
eventrem = 0;
for jj = 0:T-1 % loop over time increments
    t = jj*dt; % current time
  radbase = 0;
  for ii = 1:iimax % now loop over annuli
    radin = (ii-1)*dr;
    radout = radin+dr;
    radmid = (radin+radout)/2;
    area = pi*(radin+radin+dr)*dr*dt;
    eventstmp = Nsoln(radmid,t+dt/2,M,D,a)*area;
    % derem one pair of the following code lines
    % truncation method .. the default
    %eventstmp = eventstmp + eventrem; % number of events
    %events = floor(eventstmp); % trunc to integer and save remainder
    % rounding method
    eventstmp = max(0,eventstmp + eventrem);% number of events, must be >=0
    events = round(eventstmp); % round to integer and save remainder    
    eventrem = eventstmp-events;
    if (events ~= 0) 
      radii = radbase + rand(1,events)*(radout - radbase); % random radii
      times = t + rand(1,events)*dt; % random times 
      switch RadialDist
        case 0 % Partial Circle, random location in sector
          start = rand*2*pi;          
          theta = linspace(start,start+2*pi,events+1);
          theta = theta(1:events);
          theta = theta + rand(1,events)*2*pi/events;
        case 1 % Partial Circle, center location in sector
          % it is instructive flip rem'ing of next two lines
          start = rand*2*pi;
          %start = 0;      
          theta = linspace(start,start+2*pi,events+1);
          theta = theta(1:events) +pi/events;
        otherwise  % Full Circle and anything else
          theta = rand(1,events)*2*pi;
      end
      % convert to screen coordinates
      xcrime = X + radii.*cos(theta);
      ycrime = Y + radii.*sin(theta);
      revents = rand(events,1);
      eventoccurred = ( revents <= P); %  actual observed events
      for kk = 1:events % loop over each event
        if (eventoccurred(kk) && (ndx<=MaxEvents))
            if ((xcrime(kk)>=A) && (xcrime(kk)<=B) && (ycrime(kk)>=C) ...
                    &&  (ycrime(kk)<=E))
                retval(ndx,:) = [xcrime(kk),ycrime(kk),times(kk)];
                ndx = ndx+1;
            end
        end
      end
      radbase = radout; % set up base radius for next iteration
    end
  end
end
retval = retval(1:ndx-1,:); % trim returned vector
 