-- Copyright (C) International Computer Science Institute, 1996.  COPYRIGHT  --
-- NOTICE: This code is provided "AS IS" WITHOUT ANY WARRANTY and is subject --
-- to the terms of the SATHER LIBRARY GENERAL PUBLIC LICENSE contained in    --
-- the file "Doc/License" of the Sather distribution.  The license is also   --
-- available from ICSI, 1947 Center St., Suite 600, Berkeley CA 94704, USA.  --
--------> Please email comments to "sather-bugs@icsi.berkeley.edu". <----------

-- fish2.sa:  Fish is subject to gravitational attraction
-- author:    Boris Vaysman
-- date:      1-17-96
-------------------------------------------------------------------
class MAIN is
   shared nfish:INT;  -- number of fish in simulation
   shared tsteps:INT; -- simulation duration
   shared psteps:INT; -- display interval
   shared nfish_per_cluster:INT;
   shared nfish_per_thread:INT;
   
   shared display:DISPLAY;
   shared mutex1, mutex2:MUTEX;  --mutual exclusion locks
   
   main(args:ARRAY{STR}) is
      global_fish:ARRAY{FLIST{FISH}};  -- all fish. Each thread manages an
                                       -- an element in the array
      
      global_max_acc:FLT;
      global_max_speed:FLT;
      global_sum_speed_sq:FLT;
      
      display:DISPLAY;
      if args.size < 4 then
	 #OUT + "Usage " + args[0] + "#fish duration display_interval\n";
	 return;
      end;
      nfish := #(args[1]);      tsteps := #(args[2]);
      psteps := #(args[3]);
      
      nfish_per_cluster := nfish / clusters;
      nfish_per_thread := nfish_per_cluster / cluster_size;
      
      display := #(256.0);
      nthreads::=clusters*cluster_size;
      global_fish := #(nthreads);
      SYSTEM::thr_setconcurrency(cluster_size);
      
      -- do initialization: each thread creates its fish
      parloop
	 i::=0.upto!(nthreads-1);
      do
	 global_fish[i] := breed_my_fish(i, nfish_per_thread);
      end;
      
      mutex1 := #MUTEX;  -- create mutex locks
      mutex2 := #MUTEX;      
      
      parloop
	 i::=0.upto!(nthreads-1);
	 my_fish ::= global_fish[i]; -- fish managed by a thread
	 
      do
	 loop 
	    max_acc:FLT := 0.0;
	    max_speed:FLT := 0.0;
	    sum_speed_sq:FLT := 0.0;
	    dt ::= 0.01;  -- initial time step
	    
	    t::=0.upto!(tsteps);
	    -- compute forces acting on my fish
	    loop
	       fish ::= my_fish.elt!;
	       fish.clear_force;
	       -- compute contributions of all other fish (local & global)
	       loop
		  fishBlock ::= global_fish.elt!;
		  loop
		     fish.compute_force(fishBlock.elt!);
		  end;
	       end;
	    end;
	    
	    sync; -- need a barrier before position of fish could be changed

	    move_my_fish(my_fish, dt);
	    compute_my_stats(my_fish, inout max_acc, 
			     inout max_speed, inout sum_speed_sq);
	    sync; -- synchronization is necessary before computing global stats
	    
	    -- compute global stats
	    lock mutex1 then
	       global_max_acc := global_max_acc.max(max_acc);
	       global_max_speed := global_max_speed.max(max_speed);
	       global_sum_speed_sq := global_sum_speed_sq + sum_speed_sq;
	    end;
	    dt := (0.05*global_max_speed/global_max_acc).min(0.1);
	    
	    if t.mod(psteps)=0 then
	       lock mutex2 then  
		  -- need to lock since X calls aren't thread safe
		  draw_my_fish(display, my_fish);
	       end;
	    end;
	 end;
      end;
   end;
   
   -- each thread executes this to initialize its fish
   breed_my_fish(thread_n:INT, n:INT):FLIST{FISH} is
      res::= #FLIST{FISH}(n);
      -- Initialize 
      loop
	 k::=1.upto!(n);
	 res := res.push(#FISH(thread_n*n+k));
      end;
      return res;
   end;
   
   move_my_fish(my_fish:FLIST{FISH}, dt:FLT) is
      loop my_fish.elt!.compute_pos(dt); end;
   end;
   
   draw_my_fish(display:DISPLAY,my_fish:FLIST{FISH}) is
      loop 
	 display.draw_fish(my_fish.elt!);
      end;
   end;   
   
   compute_my_stats(my_fish:FLIST{FISH}, inout max_acc:FLT, 
		    inout max_speed:FLT, inout sum_speed_sq:FLT) is
      -- statistics per thread
      loop
	 fish ::= my_fish.elt!;
	 max_speed := max_speed.max(fish.vel.length);
	 sum_speed_sq := sum_speed_sq + fish.vel.square_length;
	 -- compute acceleration
	 acc ::= fish.force;
	 acc.scale_by(1.0/fish.mass);
	 max_acc := max_acc.max(acc.length);
      end;
   end;
end;

---------------------------------------------------------------
class FISH is
   attr pos: VEC;  -- fish position
   attr vel: VEC;  -- fish velocity
   attr mass: FLT; -- fish mass
   attr force:VEC; -- a force currently acting on the fish
   shared pi:FLT:=3.1415927;
   shared radius:FLT := 4.0;
   attr vec:VEC;   -- a temporary buffer need for some computation to avoid GC
   compute_pos(dt:FLT) is
      -- computes a new position and velocity given the force
      -- This will coded in a more pleasant way when we have a parallel gc
      pos[0] := pos[0] + vel[0] * dt;
      pos[1] := pos[1] + vel[1] * dt;      
      vel[0] := vel[0] + (force[0]/mass)*dt;
      vel[1] := vel[1] + (force[1]/mass)*dt;      
   end;
   
   compute_force(fish:FISH) is
      -- compute force on self due to an argument and add it to 
      -- the force acting on the fish
      vec[0] := (fish.pos[0] - pos[0]);
      vec[1] := (fish.pos[1] - pos[1]);      
      dist_sq ::= vec.square_length.max(0.01);
      -- Use the 2-dimensional gravity rule: F = d * (GMm/d^2)
      grav_base ::= (mass * fish.mass)/dist_sq;
      vec.scale_by(grav_base);
      force.add(vec);
   end;

   clear_force is
      force[0] := 0.0;
      force[1] := 0.0;
   end;
   
   create(k:INT):SAME is
      -- put things in a circle
      res::=new;
      res.pos := #(2);
      res.vel := #(2);
      res.force:=#(2);
      res.vec := #(2);
      delta ::= 2.0*pi/MAIN::nfish.flt;
      
      res.pos[0] := radius*(k.flt*delta).cos;
      res.pos[1] := radius*(k.flt*delta).sin;
      res.vel[0] := 0.0;
      res.vel[1] := 0.0;
      res.mass := 1.0;
      return res;
   end;
end;