import java.io.*;
import java.util.*;

public class Fish {
    public int       proc;
    public double    m;
    public Cart2d    f;
    public Cart2d    pos;
    public Cart2d    vel;
    public Fish() { 
	pos = new Cart2d(); 
	vel = new Cart2d(); 
	f = new Cart2d();
    }
    public Fish(double m, double x, double y) {
	this.m = m;
	pos = new Cart2d(x, y);
	vel = new Cart2d();
        f = new Cart2d();
	proc = 0;
    }

    public void timestep(double dt) {
	pos.x += dt*vel.x;
	pos.y += dt*vel.y;
	vel.x += dt*f.x/m;
	vel.y += dt*f.y/m;
	Aquarium.bounce(this);
    }
    
    public String toString() { return "("+m+","+f+","+vel+","+pos+")"; }
}

public class Cart2d {
    public double x, y;
	 
    public Cart2d() { x = 0.0; y = 0.0; }
    
    public Cart2d(double x, double y) {
	this.x = x;
	this.y = y;
    }
    
    public void subtract(Cart2d a, Cart2d b) {
	x = a.x-b.x;
	y = a.y-b.y;
    }
    
    public double magnitude() { return (x*x+y*y); }
    
    public double distance(Cart2d a) {
	return Math.sqrt((x-a.x)*(x-a.x)+(y-a.y)*(y-a.y));
    }
    
    public String toString() { return "("+x+","+y+")"; }
}

public class Aquarium {

    public static int n_fish = 10;

    public static double max_norm_change = 0.01;
    public static double max_dt = 0.1;
    public static double min_dt = 0.0;
    public static double end_time = 0.25;

    public static int max_steps = 1000;

    public static double min_r = 0.1;
    public static double min_r2 = min_r * min_r;
    public static double G = 1.0;

    public static double output_interval = 0.05;
    public static String output_file = "stdout";

    public static int outputp = 1;
    public static int circlep = 1;

    public static double FISH_MASS = 1.0;
    public static double LEFT_WALL = 1.0;
    public static double RIGHT_WALL = 2.0;
    public static double WALL_SEP = RIGHT_WALL-LEFT_WALL;

    public static double max(double a, double b) {
	if(a>b) return a;
	return b;
    }

    public static double min(double a, double b) {
	if(a<b) return a;
	return b;
    }

    public static double scale_coord(double c) {
	return (c - LEFT_WALL) / WALL_SEP;
    }

    public static double unscale_coord(double c) {
	return (c * WALL_SEP) + LEFT_WALL;
    }

    public static void printUsage() {
	if(Ti.thisProc() == 0) {
	    System.out.println("USAGE: Aquarium [-h|-?] [-u|-c] [-n #] [-t #.##] [...]");
	    System.out.println("\t-h or -?: Display this message.");
	    System.out.println("\t-u: Uniformly distribute the fish.");
	    System.out.println("\t-c: Place the fish evenly along a circle. (default).");
	    System.out.println("\t-n #: Use # fish. (default: "+n_fish+").");
	    System.out.println("\t-t #.##: Set end-time to #.##. (default: "+end_time+")");
	    System.out.println("\t-d #.##: Clamp the minimum time change to #.##. (default: "+min_dt+")");
	    System.out.println("\t-D #.##: Clamp the maximum time change to #.##. (default: "+max_dt+")");
	    System.out.println("\t-x #.##: Try to keep the max speed around #.##. (default: "+max_norm_change+")");
	    System.out.println("\t-s #: Clamp max. number of steps to #. (default: "+max_steps+")");
	    System.out.println("\t-G #.##: Set gravitational constant to #.##. (default: "+G+")");
	    System.out.println("\t-r #.##: Clamp min. dist. between fish to #.##. (default: "+min_r+")");
	    System.out.println("\t-O #.##: Set min. output interval to #.##. (default: "+output_interval+")");
	    System.out.println("\t-o filename: Output locations of fish. (default: "+output_file+")");
	    System.exit(1);
	}
    }

    public static void get_options(String single [] single args) {
	int i;
	for(i = 0; i < args.length; i++) {
	    if(args[i].indexOf("-h") != -1 || args[i].indexOf("-?") != -1) { printUsage(); }
	    else if(args[i].indexOf("-u") != -1) { circlep = 0; }
	    else if(args[i].indexOf("-c") != -1) { circlep = 1; }
	    else if(args[i].indexOf("-n") != -1) { n_fish = (int)Double.valueOf(args[i+1]).doubleValue(); i = i+1; }
	    else if(args[i].indexOf("-t") != -1) { end_time = Double.valueOf(args[i+1]).doubleValue(); i = i+1; }
	    else if(args[i].indexOf("-d") != -1) { min_dt = Double.valueOf(args[i+1]).doubleValue(); i = i+1; }
	    else if(args[i].indexOf("-D") != -1) { max_dt = Double.valueOf(args[i+1]).doubleValue(); i = i+1; }
	    else if(args[i].indexOf("-x") != -1) { max_norm_change = Double.valueOf(args[i+1]).doubleValue(); i = i+1; }
	    else if(args[i].indexOf("-s") != -1) { max_steps = Integer.valueOf(args[i+1]).intValue(); i = i+1; }
	    else if(args[i].indexOf("-G") != -1) { G = Double.valueOf(args[i+1]).doubleValue(); i = i+1; }
	    else if(args[i].indexOf("-r") != -1) { min_r = (int)Double.valueOf(args[i+1]).doubleValue(); min_r2 = min_r*min_r; i = i+1; }
	    else if(args[i].indexOf("-O") != -1) { output_interval = Double.valueOf(args[i+1]).doubleValue(); outputp = 1; i = i+1; }
	    else if(args[i].indexOf("-o") != -1) { output_file = ""+args[i+1]; outputp = 1; i = i+1; }
	    else { printUsage(); }
	}
    }

    public static void bounce(Fish F) {
	while(F.pos.x < LEFT_WALL || F.pos.x > RIGHT_WALL) {
	    if(F.pos.x < LEFT_WALL) {
		F.pos.x = 2.0*LEFT_WALL-F.pos.x;
		F.vel.x = -F.vel.x;
	    }
	    if(F.pos.x > RIGHT_WALL) {
		F.pos.x = 2.0*RIGHT_WALL-F.pos.x;
		F.vel.x = -F.vel.x;
	    }
	}
	while(F.pos.y < LEFT_WALL || F.pos.y > RIGHT_WALL) {
	    if(F.pos.y < LEFT_WALL) {
		F.pos.y = 2.0*LEFT_WALL-F.pos.y;
		F.vel.y = -F.vel.y;
	    }
	    if(F.pos.y > RIGHT_WALL) {
		F.pos.y = 2.0*RIGHT_WALL-F.pos.y;
		F.vel.y = -F.vel.y;
	    }
	}
    }


    public static single void main(String single [] single args) {
	get_options(args);                         // first check arguments

	int    single N = broadcast n_fish from 0;
	int    single MAXSTEPS = broadcast max_steps from 0;
	double single ENDTIME = broadcast end_time from 0;
	int    single i;                           // Synchronous index
	int    j, k;                               // Asynchronous indices
	double single curr_time = 0.0;
	double output_time = 0.0;
	double single dt = 0.0;
	double single max_norm;
	int    single steps;
	PrintStream FILE;
	Timer total_time, comm_time;
	int single NperProc = ((int)Math.floor(((double)N) / ((double)(Ti.numProcs()))));

        Double.setPrecision(4); /* DOB: ensure deterministic output */

	if(Ti.thisProc() == 0) {
	    if(outputp == 1) {
		try { 
                    if (output_file.equals("stdout")) FILE = System.out;
		    else  FILE = new PrintStream(new FileOutputStream(output_file)); 
		    FILE.println("n_fish: "+N);
		}
		catch(IOException fileerr) { System.out.println("IO Error: "+fileerr); System.exit(1); }
	    }
	}

	total_time = new Timer();
	comm_time = new Timer();

	int startIndex = NperProc*Ti.thisProc();
	int endIndex = NperProc*(Ti.thisProc()+1);
	if(Ti.thisProc() == Ti.numProcs()-1) endIndex = N;

	Fish [] single particles = new Fish[N];
	Fish [] myparticles = new Fish[endIndex-startIndex];

	total_time.start();

	comm_time.start();
	Ti.barrier();
	comm_time.stop();

	/* This initialization code you can try to improve (what is it really doing now?). */
	Random R = new Random(100); /* DOB: ensure deterministic output */
	for(i = 0; i < N; i++) {
	    int single proc = Math.min(Ti.numProcs()-1, ((int)Math.floor(((double)i)/((double)(NperProc)))));
	    
	    if(Ti.thisProc() == proc) {
		double mass = 1.0, x, y;
		if(circlep == 1) {
		    double angle = i * (2.0 * Math.PI / N);
		    x = unscale_coord(0.5*Math.cos(angle)+0.5);
		    y = unscale_coord(0.5*Math.sin(angle)+0.5);
		} else {
		    x = unscale_coord(R.nextDouble());
		    y = unscale_coord(R.nextDouble());
		}
		particles[i] = new Fish(mass, x, y);
		particles[i].proc = proc;
	    }

	    comm_time.start();
	    particles[i] = broadcast particles[i] from proc;
	    comm_time.stop();
	}

	for(output_time = 0.0, curr_time = 0.0, steps = 0; 
	    curr_time <= ENDTIME && steps < MAXSTEPS; 
	    curr_time += dt, ++steps) {

	    /* Wait for everyone to catch up. Is this barrier necessary? */
	    comm_time.start();
	    Ti.barrier();
	    comm_time.stop();

	    /* Output results to a file if desired. */
	    if(Ti.thisProc() == 0 && outputp == 1 && curr_time >= output_time) {
		double cg_x = 0.0;
		double cg_y = 0.0;
		FILE.println(""+curr_time+" ("+dt+")");
		for(k = 0; k < N; k++) {
		    cg_x += particles[k].pos.x;
		    cg_y += particles[k].pos.y;
		    FILE.println("  "+k+":"+" ("+scale_coord(particles[k].pos.x)+", "+scale_coord(particles[k].pos.y)+")");
		}
		cg_x /= N;
		cg_y /= N;
		FILE.println("  cg: ("+scale_coord(cg_x)+", "+scale_coord(cg_y)+")");

		output_time = curr_time + output_interval;
	    }

	    /* Copy fish locally. */
	    for(j = startIndex; j < endIndex; j++) {
		myparticles[j-startIndex] = particles[j];
		particles[j].proc = Ti.thisProc();
	    }

	    /* Modify forces on local fish from all the fish at once. */
	    for(j = 0; j < myparticles.length; j++) {
		myparticles[j].f.x = 0.0;
		myparticles[j].f.y = 0.0;
		for(k = 0; k < N; k++) {
		    double dx = particles[k].pos.x - myparticles[j].pos.x;
		    double dy = particles[k].pos.y - myparticles[j].pos.y;
		    if(dx != 0.0 || dy != 0.0) {
			double r2 = Math.max(dx*dx+dy*dy, min_r2);
			double r = Math.sqrt(r2);
			double a = G * FISH_MASS / r2;
			myparticles[j].f.x += a * dx / r;
			myparticles[j].f.y += a * dy / r;
		    }
		}
	    }

	    /* Compute the norm to determine timestep. */
	    double local_max_norm = 0.0;
	    for(j = 0; j < myparticles.length; j++) {
		local_max_norm = Math.max(local_max_norm, Math.abs(myparticles[j].vel.x));
		local_max_norm = Math.max(local_max_norm, Math.abs(myparticles[j].vel.y));
		local_max_norm = Math.max(local_max_norm, Math.abs(myparticles[j].f.x));
		local_max_norm = Math.max(local_max_norm, Math.abs(myparticles[j].f.y));
	    }

	    comm_time.start();
	    max_norm = Reduce.max(local_max_norm);
	    comm_time.stop();

	    dt = broadcast Math.min(Math.max(max_norm_change / max_norm, min_dt), max_dt) from 0;

	    /* Update the positions of the local particles based on timestep. */
	    for(j = 0; j < myparticles.length; j++) myparticles[j].timestep(dt);

	    /* Copy the particles to the other processors. Is this approach the fastest? */
	    comm_time.start();
	    for(i = 0; i < N; i++) {
		int single proc = broadcast particles[i].proc from 0;
		particles[i] = broadcast particles[i] from proc;
	    }
	    comm_time.stop();
	}
	total_time.stop();

	double maxcommtime = Reduce.max(comm_time.millis(), 0);
	double totalcommtime = Reduce.add(comm_time.millis(), 0);
	double maxtotaltime = Reduce.max(total_time.millis(), 0);
	double totaltotaltime = Reduce.add(total_time.millis(), 0);

        Ti.barrier();
	if(Ti.thisProc() == 0) {
	    System.out.println("done."); /* DOB: ensure deterministic output */
/*
	    System.out.println("Number of processors: "+Ti.numProcs());
	    System.out.println("Total combined time taken: "+totaltotaltime);
	    System.out.println("Maximum time taken for a process: "+maxtotaltime);
	    System.out.println("Total combined time in communication: "+totalcommtime);
	    System.out.println("Maximum time in communication for a process: "+maxcommtime);
*/
	    if(outputp == 1 && FILE != System.out) FILE.close();
	}
    }
}