package cylinder1;

/**
 * Title: cylinder1.java
 * Description: This program performs a finite element analysis
 *   on the interface of an ideal fluid spinning at a single
 *   rotational velocity in a cylinder spinning about its axis.
 *   The x,y coordinates of the interface are saved to a text file,
 *   the coordinates being take with respect to x=0 along the axis of rotation
 *   and y=0 at the bottom of the cylinder.
 *
 *
 * @author Rotating cylinder group, MATH 2413A, Fall 2002
 * @version 1.0
 */
import java.io.*;
public class simulation {

  public simulation() {
  }
  public static void main(String[] args) {
    simulation simulation1 = new simulation();

   //The radius of the circumference of a circle to its diameter
   //in Euclidean space
   double pi=3.141592;
  //unused, can be implemented to model surface tension
  double sigma=.15432358;
  //Density of water in lbm/in^2
  double rho=0.036127292;
  //acceleration of gravity in in^2/s
  double g=386.08858;
  //angular velocity of the cylinder
  double omega=6.86;
  //radius of the cylinder in inches
  double r=1.909859317;
  //height of the water level at rest in inches
  double h0=3.3750;
  //number of times to try to optimize the interface
  int iterations=1000;
  //resolution of the simulation
  int res = 100;
  //width and length of each block of water
  double xscale=2.0/res;
  //height of each block of water
  double yscale=h0/res;

  double xs2=xscale*xscale;
  double ys2=yscale*yscale;

  //amount by which to adjust each stack of water each iteration
  double increment=0.01;

  double energy=0;
  double oldenergy;
  double penergy=-1;

  //the water is modeled by an array, the index being the x coordinate
  //of the stack of water divided by xscale, and the corresponding element
  //of the array being the y coordinate of the top of the stack of water
  //divided by yscale
  double [] cols = new double[res];
  for (int i=0; i<res; i++)
    cols[i]=res;
  for (int iter=0; iter<iterations;iter++)
  {
    //if, after shifting all the stacks, the energy of the system is the
    //same, we have reached an equilibrium state
    if (penergy==energy)
      break;
    penergy=energy;
    if (iter%10==0)
      System.out.println("Iteration #"+iter);
    oldenergy=0;
    for (int i=0; i<res; i++)

    {
      double Icm;
      double mr2;
      double ke;
      double ge;
      //calculate the energy of the stack of water
      // = rotational kinetic energy + gravitational potential energy
      Icm = (1.0/6.0)*xs2*xs2*yscale*cols[i]*rho;
      mr2 = xs2*yscale*cols[i]*rho*i*i*xs2;
      ge=rho*g*xs2*ys2*cols[i]*cols[i]/2;
      ke=8*omega*omega*(Icm+mr2)/2;
      oldenergy+=-ke+ge;
    }
    //This loop modifies the heights of the stacks in hopes of
    //achiving a lower energy state.
    for (int i=0; i<res-1;i++)
    {
      if (cols[i]>=0)
      {
      for (int k=1; k<res-i; k++)
      {
      cols[i]-=increment;
      cols[i+k]+=increment;
      energy=0;
      for (int j=0; j<res; j++)

    {
      double Icm;
      double mr2;
      double ke;
      double ge;
      Icm = (1.0/6.0)*xs2*yscale*cols[j]*rho*(xs2);
      mr2 = xs2*yscale*cols[j]*rho*j*j*xs2;
      ge=rho*g*xs2*ys2*cols[j]*cols[j]/2;
      ke=8*omega*omega*(Icm+mr2)/2;
      energy+=-ke+ge;
     }
      //If the energy of the system after the adjusting is more than it was
      //before the adjusting, return to the orginal state.
      if (oldenergy<=energy)
      {
        cols[i]+=increment;
        cols[i+k]-=increment;
        }
        }
      }
      else break;
    }
  }
  //write the coordinates of the top of the stacks to a file
  PrintWriter fileOut;

  try
  {
    fileOut = new PrintWriter( new FileWriter( "feofea.txt" ) );
    for (int i=0; i<res; i++)
      fileOut.println(new Double(i*xscale).toString()
      + ","+new Double(cols[i]*yscale).toString());

    fileOut.close();
  }
  catch ( Exception e )
  {
    System.err.println(
    "Error when attempting to open, write to, or close feofea.txt" );
  }
 System.out.println("Done.");

}
}