import java.io.*;
import java.lang.UNIXProcess;

class Fish {
  /** Construct a fish with all fields 0. 
   */
  public  Fish () { px = py = vx = vy = mass = 0; }
  /** Construct a fish with the given fields. 
   */
  public  Fish (double x, double y, double v_x, double v_y, double m) 
      { px = x; py = y; vx = v_x; vy = v_y; mass = m; }

  /** Return a copy of this this Fish.
   */
  public Fish copy () {
    return new Fish (px, py, vx, vy, mass);
  }

  /** Return a string version of the Fish.
   */
  public String toString () {
    return ("px = " + px + 
  	    ", py = " + py + 
  	    ", vx = " + vx + 
	    ", vy = " + vy + 
	    ", mass = " + mass);
  }

  /* Data members */
  public double px, py;  // position coordinates
  public double vx, vy;  // velocity 
  public double mass;
}

class FishSim {

  /** Print help message and abort program. 
   */
  private static void usage () {
    System.err.println("usage: FishSim " +
		       "-fish <n> " + 
		       "-dtmin <x> " + 
		       "-dtmax <x> " +
		       "-dxmin <x> " +
		       "-grav <x> " +
		       "-endtime <x> " 
		       );
    System.exit (1);
  }

  /** Report badParam and abort program 
   */
  private static void badParam () {
    System.err.println ("FishSim: bad parameter value\n");
    System.exit (1);
  }

  /** Convert current element of argv array into int and return it. 
    */
  private static int intArg (int current, String single [] single argv) {
    int ret = 0;
    if (current >= argv.length) usage ();
    try { 
      ret = Integer.parseInt(argv[current]); 
    } catch (NumberFormatException x) { 
      usage (); 
    }
    return ret;
  }

  /** Convert current element of argv array into double and return it. 
    */
  private static double floatArg (int current, String single [] single argv) {
    double ret = 0.0;
    if (current >= argv.length) usage ();
    ret = (new Double(argv[current])).doubleValue(); 
    return ret;
  }

  /** Convert current element of argv array into boolean and return it. 
    */
  private static boolean boolArg (int current, String single [] single argv) {
    boolean ret = false;
    if (current >= argv.length) usage ();
    ret = (new Boolean(argv[current])).booleanValue(); 
    return ret;
  }

  /** Check that argument values are reasonable and report their values.
   */
  private static void checkArgs(int myFishCount, double myDtMin, 
			        double myDtMax, double myDxMin, 
                                double myGravConst, double myEndTime) {
    if (myFishCount < 1 || 
        myDtMin > myDtMax ||
        myDxMin < 0 ||
        myEndTime < 0 ) {
      badParam ();
    }

    if (Ti.thisProc() == 0) {
      System.out.println("fish = " + myFishCount);
      System.out.println("min dt = " + myDtMin);
      System.out.println("max dt = " + myDtMax);
      System.out.println("min dx = " + myDxMin);
      System.out.println("gravitational constant = " + myGravConst);
      System.out.println("simulation end = " + myEndTime);
    }
  }

  public static DataOutputStream setupDisplay(int displaySize, int fishCount) {
    DataOutputStream dos = null;
    if (Ti.thisProc() == 0 && displaySize != 0) {
      Process guiProcess = null;
      try {
	guiProcess = Runtime.getRuntime().exec("java FishAWT");
      } catch (IOException e) { System.out.println(e.getMessage()); }
      OutputStream os = guiProcess.getOutputStream();
      dos = new DataOutputStream(new BufferedOutputStream(os));
    }
    // display.writeInt(displaySize);   currently ignored
    try {
      dos.writeInt(fishCount);   
    } catch (Exception e) { System.err.println("Broken display pipe."); }
    return dos;
  }

  public static single void showOcean(Fish [1d] single [1d] allFish,
                                      DataOutputStream display) {
    try {
      foreach (p within allFish.domain()) {
        foreach (fi within allFish[p].domain()) {
          Fish aFish = allFish[p][fi];
          display.writeDouble(aFish.px);   
          display.writeDouble(aFish.py);   
          display.flush();
	}
      }
    } catch (Exception e) { System.err.println("Broken display pipe."); }
  }


  /** processes the input parameters and invokes the simulation. 
   */
  public static single void main (String single [] single argv) {

    /* Initialize some variable common on any parallel program */
    Reduce single reduce = new Reduce();
    int myProc = Ti.thisProc();
    int single numProcs = Ti.numProcs();

    /* Each processor may have a distinct value of the following variables */
    double myMaxSqrSpeed = 0.0d;
    int myFishCount = 0;
    int myFirstFish = 0; 

    /* The following variables have the same value on all processors. */
    int single allFishCount;
    double single allEndTime;
    double single allDtMin;
    double single allDtMax;
    double single allDxMin;
    double single allGravConst;
    double single allMaxSqrSpeed;
    double single allTimestep;
    double single allDt;

    /* The following are used only by processor 0 */
    double myEndTime;
    double myDtMin;
    double myDtMax;
    double myDxMin;
    double myGravConst;
    double myTimestep;
    double myDt;
    int myDisplay;
  
    // argument defaults
    myFishCount = 50; 
    myEndTime = 50.0d;  
    myDtMin = 0.0001d; 
    myDtMax = 0.01d;
    myDxMin = 0.005d;
    myGravConst = 1.0d;
    myDisplay = 200;
  
    // process arguments
    int argc = 0;
    while (argc < argv.length) {
      if (argv[argc].equals("-fish")) {
        myFishCount = intArg (argc+1, argv);
        argc+=2;
      } else if (argv[argc].equals("-dtmin")) {
        myDtMin = floatArg (argc+1, argv);
        argc+=2;
      } else if (argv[argc].equals("-dtmax")) {
        myDtMax = floatArg (argc+1, argv);
        argc+=2;
      } else if (argv[argc].equals("-dxmin")) {
        myDxMin = floatArg (argc+1, argv);
        argc+=2;
      } else if (argv[argc].equals("-grav")) {
        myGravConst = floatArg (argc+1, argv);
        argc+=2;
      } else if (argv[argc].equals("-endtime")) {
        myEndTime = floatArg (argc+1, argv);  
        argc+=2;
      } else if (argv[argc].equals("-display")) {
        myDisplay = intArg (argc+1, argv);  
        argc+=2;
      } else { 
        usage ();
      }
    }

    checkArgs(myFishCount, myDtMin, myDtMax, myDxMin, myGravConst, myEndTime);
    DataOutputStream display = setupDisplay(myDisplay, myFishCount);

    // broadcast basic constants
  
    allFishCount = broadcast myFishCount from 0;
    allDtMin = broadcast myDtMin from 0;
    allDtMax = broadcast myDtMax from 0;
    allDxMin = broadcast myDxMin from 0;
    allGravConst = broadcast myGravConst from 0;
    allEndTime = broadcast myEndTime from 0;
  
    // allocate stuff
  
    myFishCount = allFishCount/numProcs;
    myFirstFish = myFishCount * myProc + 
                  Math.min(allFishCount % numProcs, myProc);
    if (allFishCount % numProcs > myProc) myFishCount++;
  
    RectDomain <1> single allProcs = [0:numProcs-1];
    RectDomain <1> myFishDomain = [0:myFishCount-1];
  
    Fish [1d] single [1d] allFish = new Fish [allProcs][1d];
    double [1d] force_x = new double [myFishDomain];
    double [1d] force_y = new double [myFishDomain];
    Fish [1d] myFish = new Fish [myFishDomain];  
    allFish.exchange(myFish);

    if (myProc == 0) { 
      System.out.println("Number of threads: " + numProcs);
    }
  
    /* initialize fish */
    double z = 2.0d * Math.PI / allFishCount;
    for (int i = 0; i < myFishCount; i++) {
      myFish[i] = new Fish(Math.sin (z * (i + myFirstFish)), 
  				    Math.cos (z * (i + myFirstFish)), 
  				    0,
  				    0,
  				    1.0 / allFishCount);
    }
    Ti.barrier ();
  
    long milliTime = 0;
    /* start timer */
    if (myProc == 0) {
      milliTime = System.currentTimeMillis();
    }
  
    /* loop over the time steps */
    for (allDt = allDtMin, allTimestep = 0; allTimestep < allEndTime; 
      allTimestep += allDt) {

      showOcean(allFish, display);
  
      /* reset force temporaries for next timestep */
      foreach (fi within myFish.domain()) {
        force_x[fi] = 0;
        force_y[fi] = 0;
      }
  
      /* compute forces on my fish */
      foreach (p within allProcs) {
        // move allocation out of loop
        Fish [1d] remoteFish = new Fish [allFish[p].domain()];
        remoteFish.copy(allFish[p]);
        foreach (fj within remoteFish.domain()) {
  	  Fish fishj = (Fish) remoteFish[fj].copy();
  
  	  /* compute global force on fish i */
  	  foreach (fi within myFish.domain()) {
  	    double dx = myFish[fi].px - fishj.px;
  	    double dy = myFish[fi].py - fishj.py;
  
  	    double sqrDist = Math.max(dx*dx + dy*dy, allDxMin*allDxMin);
  	    double gravBase = allGravConst * myFish[fi].mass * 
  	                      fishj.mass / sqrDist;
  
  	    force_x[fi] += dx * gravBase;
  	    force_y[fi] += dy * gravBase;
  	  }
        }
      }
    
      Ti.barrier ();
  
      /* compute new positions and velocities of my fish */
      /* compute maximum velocity of my fish */
      foreach (fi within myFish.domain()) {
        double sqrSpeed;
  
        myFish[fi].px += myFish[fi].vx * allDt;
        myFish[fi].py += myFish[fi].vy * allDt;
        myFish[fi].vx -= force_x[fi] / 
        myFish[fi].mass * allDt;
        myFish[fi].vy -= force_y[fi] / 
        myFish[fi].mass * allDt;
        sqrSpeed = (myFish[fi].vx * myFish[fi].vx) 
  	           + (myFish[fi].vy * myFish[fi].vy);
  
        myMaxSqrSpeed = Math.max (myMaxSqrSpeed, sqrSpeed);
      }  
  
      /* compute global maximum fish velocity */
      allMaxSqrSpeed = reduce.max(myMaxSqrSpeed);
  
      /* recompute dt so that the fastest particle does not move more */
      /* than the minumum distance in the next time step */
      allDt = allDxMin / (double single) Math.sqrt (allMaxSqrSpeed);
      allDt = (double single) Math.min (allDtMax, allDt);
      allDt = (double single) Math.max (allDtMin, allDt);
      if (myProc == 0) {
        System.out.println("Max Square Speed at timestep " 
  	 		   + allTimestep + ": " + allMaxSqrSpeed);
      }
    }
       
    /* stop timer; report time */
    if (myProc == 0) {
      milliTime = System.currentTimeMillis() - milliTime;
      System.out.println("time = " + milliTime/1000 + " sec");
    }
  }
}
  
/** Titanium has a standard library of reduction operations, but
 *  the implementation is not yet part of the official release.
 *  Until it is, this specialized class has a similar interface
 *  with an inefficient implementation. 
 */
class Reduce {
  
  BoxDouble [1d] single dReduce;
  int single numProcs;
  
  single Reduce () {
    numProcs = Ti.numProcs();
    dReduce = new BoxDouble [0:numProcs-1];
    dReduce.exchange(new BoxDouble());
  }
  
  single double single max(double myMax) {
    double single allResult;
 
    Ti.barrier();
    double result;
    int myProc = Ti.thisProc();
  
    dReduce[myProc].value = myMax;
  
    // Could do a tree reduction, but this is centralized for now.
    if (myProc == 0) {
      result = myMax;
      foreach (i within [1:numProcs-1]) {
      result = Math.max(result, dReduce[i].value);
    }
  }
  allResult = broadcast result from 0;
  return allResult;
  }
}

/** Objects to hold double values used in reduction 
 */  
class BoxDouble {
  public double value;
  public BoxDouble () { }
  public BoxDouble (double v) { value = v; }
}