-- 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". <----------

-- fish1.sa:  Fish are subject to external current
-- 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:FMUTEX;  -- mutual exclusion locks
   
   main(args:ARRAY{STR}) is
      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;
      nthreads::=clusters*cluster_size;
    
      display := #(256.0);
      SYSTEM::thr_setconcurrency(cluster_size);
      
      mutex1 := #FMUTEX;  --create mutex locks
      mutex2 := #FMUTEX;

      parloop
	 i::=0.upto!(nthreads-1);
      do
	 -- initialize current
	 current::=#CURRENT;
	 -- each thread initializes its fish
	 my_fish:FLIST{FISH} := breed_my_fish(i, nfish_per_thread);
	 dt ::= 0.2; -- initial time step
	 loop 
	    max_acc:FLT := 0.0;
	    max_speed:FLT := 0.0;
	    sum_speed_sq:FLT := 0.0;
	    
	    t::=0.upto!(tsteps);
	    move_my_fish(my_fish, current, dt);
	    compute_my_stats(my_fish, current, inout max_acc, 
			     inout max_speed, inout sum_speed_sq);
	    sync;  -- need to synchronize before accumulating 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.1*global_max_speed/global_max_acc).min(0.05);
	    
	    if t.mod(psteps)=0 then
	       lock mutex2 then
		  -- have to protect by a lock since X calls aren't thread safe
		  draw_my_fish(display, my_fish);
	       end;
	    end;
	 end;
      end;
   end;
   
   -- Initialize fish managed by a single thread
   breed_my_fish(thread_n:INT, n:INT):FLIST{FISH} is
      res::= #FLIST{FISH};
      loop
	 k::=1.upto!(n);
	 res := res.push(#FISH(thread_n*n+k));
      end;
      return res;
   end;
   
   -- compute new position for fish
   move_my_fish(my_fish:FLIST{FISH}, current:CURRENT, dt:FLT) is
      loop my_fish.elt!.compute_pos(current, 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}, current:CURRENT, 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;
	 acc:VEC := current.force(fish.pos);
	 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
   
   compute_pos(current:CURRENT, dt:FLT) is
      pos[0] := pos[0] + vel[0] * dt;
      pos[1] := pos[1] + vel[1] * dt;      
      vec ::= current.force(pos);
      vec.scale_by(dt/mass);
      vel.add(vec);
   end;
   
   create(k:INT):SAME is
      res::=new;
      res.pos := #(2);
      res.vel := #(2);
      res.pos[0] := k.flt*2.0/MAIN::nfish.flt - 1.0;
      res.pos[1] := k.flt*2.0/MAIN::nfish.flt - 1.0; 
      res.vel[0] := -res.pos[1];
      res.vel[1] := res.pos[0];
      res.mass := 1.0 + k.flt / MAIN::nfish.flt;
      return res;
   end;
end;

---------------------------------------------------------------
class CURRENT is
   -- returns current force at a particular point
   -- This particular function describes a whirpool-like current
   attr current_force:VEC; 

   create:SAME is
      res ::= new;
      res.current_force:=#(2);
      return res;
   end;
      
   force(pos:VEC):VEC is
      current_force[0] := -3.0 * pos[1] / pos.length.max(0.01);
      current_force[1] :=  3.0 * pos[0] / pos.length.max(0.01);
      return current_force;
   end;
end;