import java.awt.*;
import java.awt.image.*;
// 1D Lagrangian hydro; perfect gas, shockwave, no Q.
// Should be symmetrical!
class Main {
  static final int ms = 20;
  static final double dt = .003, dur = 12;
  static final boolean se = true;
  static class zone{double x, v, T;}
  static final int X = 800, Y = 700;
  // m[j].T is the temperature of the material zone between m[j].x and m[j+1].x
  // Each zone has mass 1; pv=NRt
  static zone[] m = new zone[ms+1];
  static double co = se?0:ms/2;
  static BufferedImage I = new BufferedImage(X, Y, BufferedImage.TYPE_INT_RGB);
  static WritableRaster wr = I.getRaster();
  static int[] mk = new int[] {255, 255, 255};
  static {for(int j = 0; j<ms+1; ++j) {
            m[j] = new zone(); m[j].x = j-co; m[j].v = 0; m[j].T = 1.;}}
  static double xs = X/dur, ys = Y/(ms+1);
  static void pd(int j, double t) {
     try {wr.setPixel((int)(xs*t), (int)(ys*(m[j].x+co)), mk);
      if(false) if(j< ms) wr.setPixel((int)(xs*t), (int)(ys*((m[j].x+m[j+1].x)/2+co)),
            new int[] {(int)m[j].T*20, 255, 255});}
        catch (ArrayIndexOutOfBoundsException d) {System.out.print(j);
            System.out.print("www j="); System.out.print(j);
            System.out.print(" y="); System.out.print(m[j].x);
            System.out.print(" "); System.out.println(ys*m[j].x);}}
  static void pnv(String n, Double v){
     System.out.print(n); System.out.print(v);}
  static class MainCanvas extends Canvas {
     public void paint(Graphics g){
       g.drawImage(Main.I, 0, 0, Color.red, null);}}
 
    static public void main( String[] args ){
      Frame f = new Frame( "1D Hydro" );
       Canvas can = new MainCanvas();
       f.add(can, "Center");
       f.setSize(new Dimension(X, Y+15));
       for(double t = 0; t<dur; t += dt){
          pd(0, t);
          m[0].v = 1; // Boundary condition for shock wave.
          m[ms].v = se?0:-m[0].v;
          wr.setPixel((int)(xs*t), (int)(ys*(m[0].x+co)), mk);
          double x0 = m[0].x;
          m[0].x += m[0].v*dt;
          for(int j=0; j<ms; ++j) {
             double V0 = m[j+1].x - x0; x0 = m[j+1].x;
             m[j+1].x += m[j+1].v*dt;
             double V = m[j+1].x - m[j].x, dV = V0  - V; // positive for compression
             double P = m[j].T/V;
             if(false) {
               pnv("\nm[j].T=", m[j].T);
               pnv(" P=", P);
           //    pnv(" V0=", V0);
             //  pnv(" V=", V);
               }
             m[j].T += dV*P;
             m[j].v -= P*dt; m[j+1].v += P*dt;
             pd(j+1, t); 
     }}
   f.setVisible(true);
   for(int j = 0; j<ms+1; ++j) {System.out.print(m[j].x);
      System.out.print(", "); System.out.println(m[j].v);}
}}