%function y=particlesim;
clear;

% screen parameters
N = 600;
mbintscalar(N);

screen_x  = 1000; screen_y  = 1000;  screen_margin = 100; 
vx_range  = 200; vy_range  = 200; v_tweak   = .13; v_range = 10;
vx_offset  = 100; vy_offset  = 100; particle_r = 5;  field_div = 100;

mbintscalar(screen_x); mbintscalar(screen_y); mbintscalar(screen_margin); 
mbintscalar(vx_range); mbintscalar(vy_range); mbrealscalar(v_tweak);
mbintscalar(vx_offset); mbintscalar(vy_offset); mbintscalar(particle_r); mbintscalar(field_div);


% polorizer parameters
a_ = .1; ap_ = .9;  loopN = 800; a_margin = .2; charge_mass = 500;
mbrealscalar(a_); mbrealscalar( ap_); mbintscalar( loopN); mbrealscalar( a_margin);
mbrealscalar(charge_mass);

% simulation specific parameters
charge_site_no_move = 0; mbintscalar(charge_site_no_move);
charge_site_just_one_move = 1; mbintscalar(charge_site_just_one_move);

% structure parameters
ST_x = 1; ST_y = 2; ST_vx = 3; ST_vy = 4; ST_type = 5; ST_spin = 6; 
ST_radius = 7; ST_detecttype = 8; ST_length = 8;

mbintscalar(ST_x); mbintscalar( ST_y); mbintscalar( ST_vx); mbintscalar( ST_vy); 
mbintscalar( ST_type); mbintscalar(ST_spin); mbintscalar(ST_radius); 
mbintscalar( ST_detecttype); mbintscalar( ST_length);

% enumerated type values
DT_virtual  = 0; T_particle = 0; T_antiparticle = 1;  
DT_real = 1;     T_positive = 0; T_negative = 1;

mbintscalar(DT_virtual); mbintscalar( T_particle); mbintscalar( T_antiparticle); 
mbintscalar(DT_real); mbintscalar(T_positive); mbintscalar(T_negative);

% set to off boundary values for automatic initialization
particles =   repmat([-screen_x -screen_y  0 0 0 0 0 0],N,1);

% initialize charge sites
particles( 1 , : ) = [ (screen_x/2+50) screen_y/2   0 0 T_positive -1 50 DT_real ];
particles( 2 , : ) = [ (screen_x/2-350) screen_y/2   0 0 T_negative -1 50 DT_real ];

real_particles_list =  find(  particles( : , ST_detecttype) == DT_real);
mbreal(particles); mbintvector(real_particles_list);

SHOW_PLOT = 1; PLOT_STEP = 1; PLOT_vpar = 1; marker_charge_size_ = 50;
PLOT_FRAMEPAUSE = .08;

mbintscalar(SHOW_PLOT); mbintscalar(PLOT_STEP); mbintscalar( PLOT_vpar); 
mbintscalar(marker_charge_size_); mbrealscalar(PLOT_FRAMEPAUSE);

%axis fill;
figure(1);
axis manual;
hold off;

for i = 0:loopN
% isolate virtual particles
virtual_particles = particles(:, ST_detecttype) == DT_virtual; 
mbintvector( virtual_particles);

% calculate the number of them
vpar_N = sum(virtual_particles);
mbintscalar(vpar_N);

% tweak motion of the virtual sea
particles( virtual_particles ,ST_x:ST_y) = particles(  virtual_particles ,ST_x:ST_y) + ...
        particles( virtual_particles ,ST_vx:ST_vy) .* (1 + v_tweak * randn(vpar_N,2));

particles( ~virtual_particles ,ST_x:ST_y) = particles(  ~virtual_particles ,ST_x:ST_y) + ...
        particles( ~virtual_particles ,ST_vx:ST_vy) ;
    
% create new particles for every particle out of bound 
bound_kills = find(  (particles(:,ST_x) <= 0) | (particles(:,ST_x) > screen_x) ...
                   | (particles(:,ST_y) <= 0) | (particles(:,ST_y) > screen_y));
mbintvector(bound_kills);

%O(n^2) collision check
collision =  (repmat( particles( :,ST_x)  , 1 , N)  - repmat( particles(:,ST_x)', N,1)) .^2 + ...
             (repmat( particles( :,ST_y)  , 1 , N)  - repmat( particles(:,ST_y)', N,1)) .^2;
mbreal(collision);
         
%detect them    
an_i   = find( collision < repmat( particles(:,ST_radius).^ 2, 1, N) ); mbintvector(an_i);
an_ix  = mod( (an_i-1), N ) +1;  mbintvector(an_ix);
an_iy  = mod( floor((an_i-1) / N) , N) + 1 ; mbintvector(an_iy);

% minimize comparisions to different particles
remove_ii = (an_ix < an_iy) ; mbintvector(remove_ii);
an_ix = nonzeros( remove_ii.* an_ix); mbintvector(an_ix);
an_iy = nonzeros( remove_ii.* an_iy); mbintvector(an_iy);

clear remove_ii % free unused memory space

% particle/antiparticle collisions with opposite spins only 
collision_kills =  nonzeros( repmat( ...
       (particles(an_ix, ST_detecttype) ==  DT_virtual) & ...
       (particles(an_iy, ST_detecttype) ==  DT_virtual)  & ...
       (particles(an_ix, ST_type ) ~=  particles(an_iy, ST_type)) &  ...
       (particles(an_ix, ST_spin ) ==  -particles(an_iy, ST_spin)) ...
       , 1 , 2 ).*[ an_ix an_iy]  ); 
mbintvector(collision_kills);

% polorizers collisions: virtual/real particle interaction
% APPROACH: Charge Interaction model

polorizers_reactions = (particles(an_ix, ST_detecttype ) ==  DT_real) & ...
                       (particles(an_iy, ST_detecttype ) ==  DT_virtual);
mbintvector( polorizers_reactions);

an_i_size = sum(polorizers_reactions); % allocate as much as we need them
mbintscalar( an_i_size); 

% polorizer mechanism
r_a    =  rand(an_i_size,1) < a_ ;   mbintvector(r_a);
r_ap_  =  rand(an_i_size,1) < ap_;   mbintvector(r_ap_);

% reduce an_ix an_iy to just valid polorization reactions
an_ix_pr =   an_ix(polorizers_reactions);  mbintvector( an_ix_pr);
an_iy_pr =   an_iy(polorizers_reactions);  mbintvector( an_iy_pr);

% calculate for the appropriate reaction scheme
pick_scheme = bitxor( particles( an_iy_pr, ST_spin) > 0, particles( an_ix_pr, ST_type));
keep_p    = (pick_scheme).* r_a   | (~pick_scheme).* r_ap_; mbintvector(keep_p);
convert_p = (pick_scheme).* r_ap_ | (~pick_scheme).* r_a  ; mbintvector(convert_p);  

clear pick_scheme % free memory

% exclusive destruction events  
free_particles =  ((~keep_p) &  (~convert_p)) | (( keep_p) &  (convert_p));
mbintvector( free_particles);

spin_and_type_change =  an_iy_pr( (~keep_p) & (convert_p) ); 
mbintvector( spin_and_type_change);

% polorize the particle(if necessary)
particles( spin_and_type_change, [ST_type ST_spin] ) = ...
           [ ~particles( spin_and_type_change, ST_type ) ...
             -particles( spin_and_type_change, ST_spin ) ];

polorizer_kills =  an_iy_pr(free_particles); % kill them
mbintvector(polorizer_kills);

% reduce them again
an_ix_pr = an_ix_pr( ~free_particles);
an_iy_pr = an_iy_pr( ~free_particles);
  
clear an_ix an_iy spin_and_type_change free_particles keep_p convert_p

% compenstate motion thru the charge site and put the virtual particle 
% away from the charge site
incoming_d_  =      particles( an_iy_pr , ST_x:ST_y ); 
incoming_dd_ =      incoming_d_  - particles( an_ix_pr , ST_x:ST_y ); 
incoming_ddmod2_ =  incoming_dd_(:,1) .^2 + incoming_dd_(:,2) .^2;
outgoing_r   =      particles( an_ix_pr , ST_radius); 
incoming_v_  =      particles( an_iy_pr , ST_vx:ST_vy);
incoming_vmod2_  =  incoming_v_(:,1) .^2  + incoming_v_(:,2)  .^2;
dot_dd_v_    =      dot(incoming_dd_ ,  incoming_v_ , 2);

mbreal( incoming_d_); mbreal( incoming_dd_); mbrealvector( incoming_ddmod2_);
mbrealvector( outgoing_r) ; mbreal( incoming_v_) ; mbrealvector( incoming_vmod2_);

% shift them from repeated reacting
particles( an_iy_pr, ST_x: ST_y ) =  incoming_d_  -  incoming_v_   ...
   .* repmat(  (dot_dd_v_  -  sqrt( incoming_vmod2_ .*  ...
   (  ( ( 1 + a_margin) *outgoing_r ).^2  -    incoming_ddmod2_)  ...
   +  dot_dd_v_ .^2 )) ./ incoming_vmod2_ , 1 , 2) ;

clear incoming_d_ incoming_dd_ incoming_ddmod2_ outgoing_r incoming_v_ incoming_v_  ...
      incoming_vmod2_ dot_dd_v_



if (size(an_ix_pr, 1)> 0) 
% perturb none of the particles or just one particle (depending on the settings)
 an_ix_only =  ((an_ix_pr == 1) & (charge_site_just_one_move)) & ~(charge_site_no_move) ;
 an_ix_pr = an_ix_pr( an_ix_only);
 an_iy_pr = an_iy_pr( an_ix_only);
 clear an_ix_only
end;

if size(an_ix_pr, 1) > 0
 % calculate movement to the charge site from the virtual sea
 interaction_field_ =   sortrows( ...
      [ an_ix_pr  particles( an_iy_pr,ST_spin) .* particles( an_iy_pr,ST_vx)/charge_mass ...
                  particles( an_iy_pr,ST_spin) .* particles( an_iy_pr,ST_vy)/charge_mass ;
            (N+1)   0  0] , [ 1]);

 interaction_fieldmask_     = find( diff(interaction_field_(:,1),1)  );
 interaction_fieldstrength_ = cumsum( interaction_field_(:,2:3),1);
 interaction_fieldstrength_ = diff([  0 0 ; interaction_fieldstrength_(interaction_fieldmask_, :)], 1);
 interaction_fieldorigin_   = interaction_field_( interaction_fieldmask_, 1);

% change velocity accoridingly
 particles( interaction_fieldorigin_ ,  ST_vx:ST_vy ) = ...
     particles( interaction_fieldorigin_,  ST_vx:ST_vy) -  interaction_fieldstrength_;

 clear interaction_field_ interaction_fieldmask_  interaction_fieldstrength_ ...
      interaction_fieldorigin_ an_ix_pr an_iy_pr
end;
% END METOD : Charge interaction model

% destroy wasted particles
redo_particles =  unique([bound_kills ; collision_kills; polorizer_kills]);
mbintvector(redo_particles);

clear polorizers_reactions polorizer_kills collision bound_kills  ...
      collision_kills an_i an_ix an_iy % free memory

% create new particles
new_N = size(redo_particles,1); mbintscalar(new_N);
real_N  = size(real_particles_list,1); mbintscalar(real_N);

particles( redo_particles, :) =  repmat([ 1 1 (1/20) (1/20) 1 2 1 1 ], new_N,1) .* ...
         (floor( rand(new_N,ST_length) .* repmat([ real_N real_N vx_range vy_range 2 2 0 0],new_N,1) ) ...
         + repmat([ 1 1 (-vx_offset) (-vy_offset) 0 (-0.5) particle_r DT_virtual], new_N,1));
     
% spread them randomly and let the virtual-real collision mechanism
% take care of the careful distribution, in the next round

temp_particles_xy =  [ (particles( particles( redo_particles, ST_x),ST_x)  + ... 
                       (1 - a_margin) * rand(new_N,1) .* particles( particles( redo_particles, ST_x),ST_radius) / sqrt(2)) ...
                       (particles( particles( redo_particles, ST_y),ST_y) + ...
                       (1 - a_margin) * rand(new_N,1) .* particles( particles( redo_particles, ST_y) , ST_radius) /sqrt(2))...
                      ]; 
              
mbreal(temp_particles_xy);
particles( redo_particles, ST_x: ST_y )  = temp_particles_xy;

% vx, vy equally distributed in all directions, and in all 'wavelengths' [ v_x, v_y]
particles( redo_particles, ST_vx: ST_vy) =  v_range* particles( redo_particles, ST_vx:ST_vy) .* ...
                                            repmat( rand(new_N,1) ./ sqrt( particles( redo_particles, ST_vx) .^2 + ...
                                                                     particles( redo_particles, ST_vy) .^2) , 1, 2);
clear redo_particles temp_particles

if SHOW_PLOT & (mod( i , PLOT_STEP) == 0)

 % charge sites 
positive_ =  particles(  (particles(:,ST_type) == T_positive) & ~virtual_particles , ST_x:ST_y);
negative_ =  particles(  (particles(:,ST_type) == T_negative) & ~virtual_particles , ST_x:ST_y);

% display them
plot( [ screen_x/2 ] , [ screen_y/2 ], 'ko');
hold on;
axis equal;
axis([ (-screen_margin) (screen_x+screen_margin) (-screen_margin) (screen_y+screen_margin)]);

plot(positive_(:,ST_x), positive_(:,ST_y), 'ro', 'MarkerSize',marker_charge_size_);
plot(negative_(:,ST_x), negative_(:,ST_y), 'bo', 'MarkerSize',marker_charge_size_);


if PLOT_vpar
% seperate particle/anti particles
particle      =  particles(  (particles(:,ST_type) == T_particle   )  & virtual_particles , ST_x:ST_y);
anti_particle =  particles(  (particles(:,ST_type) == T_antiparticle) & virtual_particles , ST_x:ST_y);

% display them
plot(particle(:,ST_x), particle(:,ST_y), 'r.' ,  ...
     anti_particle(:,ST_x), anti_particle(:,ST_y), 'b.' );
axis([ -screen_margin (screen_x+screen_margin) -screen_margin (screen_y+screen_margin)]);

clear particle anti_particle
end;

% TODO: use gauss_maks to show that charge model obeys Gauss's theorem

gauss_mask = virtual_particles ; 
clear virtual_particles positive_ negative_ % free unused memory
mask_N = sum(gauss_mask);
mbintscalar(mask_N);

% calculation of field lines
field_  = sortrows([  floor(particles( gauss_mask,ST_x) / field_div)   ...
                      ceil(particles( gauss_mask,ST_y) / field_div)  ...
                      particles(gauss_mask,ST_spin) .* particles(gauss_mask,ST_vx) ...
                      particles(gauss_mask,ST_spin) .* particles(gauss_mask,ST_vy) ] , [1 2]);
              
fieldmask_     = find(sum( abs(diff([ field_(:,1:2) ] , 1)) ,2) );
fieldstrength_ = cumsum( [ ones(mask_N,1)  field_(:,3:4)],1);
fieldstrength_ = diff([ 0 0 0 ; fieldstrength_( fieldmask_, :)], 1);
fieldstrength_ = fieldstrength_(:,2:3) ; %./  repmat(fieldstrength_(:,1),1,2); 
fieldorigin_   = field_( fieldmask_, 1:2);

clear gauss_region gauss_mask;
               
% show the field lines
quiver( (fieldorigin_(:,1) + .5) * field_div ,(fieldorigin_(:,2) - .5) *field_div , ...
         fieldstrength_(:,1) , fieldstrength_(:,2), 'k');
hold off;

clear field_ fieldmask_ fieldstrength_  fieldorigin_ % free memory

pause(PLOT_FRAMEPAUSE);

end;

end;


