/*********************************************************
  Instances of the Configuration class are complete
  descriptions of the system at a given time.
*********************************************************/

import java.util.Random;

public class Configuration {
  // Size of triangulation
  int xsize,ysize;
  // Data to store the configuration
  Vec3[][] x;
  Vec3[][] na;
  Vec3[][] nb;
  // Random number generator
  Random random;
  // Workspace for localBending
  Vec3[][]  localna;
  Vec3[][]  localnb;
  // One fixed point to remove zeroth order mode
  int xfixed,yfixed;
  // Used to draw
  private int idraw,jdraw,kdraw;
  private int[][] facesx, facesy;
  // ** this bit will only work with square triangulations ** 
  // ** used to work out hexagon to display **
  private int halfSize;
  // Constants
  static float sqrt2=(float)Math.sqrt(2.0);
  static float invSqrt2=(float)(1.0/Math.sqrt(2.0));
  static float sqrt3over2=(float)(Math.sqrt(3.0)/2.0);
  static float invSqrt3over2=(float)(2.0/Math.sqrt(3.0));

  Configuration(int xxsize, int yysize) {
    int i,j;
    xsize=xxsize;
    ysize=yysize;
    xfixed=xsize/2;
    yfixed=ysize/2;
    // create a new configuration object
    x=new Vec3[xsize][ysize];
    na=new Vec3[xsize][ysize];	  // could just be size-1 big
    nb=new Vec3[xsize][ysize];	  // could just be size-1 big
    for (i=0;i<xsize;i++) {
      for (j=0;j<ysize;j++) {
        x[i][j] = new Vec3((float)((float)i-(float)j*0.5),
                           (float)j*invSqrt3over2,(float)0.0);
        na[i][j] = new Vec3();
        nb[i][j] = new Vec3();
      }
    }
    // create a new Random object
    random = new Random();
    // create workspace vectors for localBending
    localna=new Vec3[2][2];
    localnb=new Vec3[2][2];
    for (i=0;i<2;i++) {
      for (j=0;j<2;j++) {
        localna[i][j]=new Vec3();
        localnb[i][j]=new Vec3();
      }
    }
    // for drawing routines
    facesx=new int[xsize][ysize];
    facesy=new int[xsize][ysize];
  }


  public void init(float temp) {
    /* it is intended that temp will be used to set the initial
       roughness in some way */
    int i,j;
    // set initial state (flat)
    for (i=0;i<xsize;i++) {
      for (j=0;j<ysize;j++) {
        x[i][j].setValue((float)((float)i-(float)j*0.5),
                        (float)j*invSqrt2,(float)0.0);
      }
    }
    // calculate normals
    this.normals(); 
  }


  private void normals() {
    int i,j;
    for (i=0;i<(xsize-1);i++) {
      for (j=0;j<(ysize-1);j++) {
        na[i][j].normal(x[i][j],x[i+1][j],x[i+1][j+1]);
        nb[i][j].normal(x[i][j],x[i+1][j+1],x[i][j+1]);
      }
    }
  }


  public int update(float sigma, float kappa, int numTrials) {
    int i,j,k,accepted;
    float energy;
    Vec3 trial=new Vec3();
    // Set accepted conter to zero
    accepted=0;
    for (k=0;k<numTrials;k++) {
      // pick a random site
      // this is rather messy as the following lines can actually give
      // xsize and ysize as their max values which we screen out with
      // the if statement
      i=(int)(random.nextFloat()*xsize);
      j=(int)(random.nextFloat()*ysize);
      // remove zero order mode by fixing one site
      if (((i!=xfixed) || (j!=yfixed)) && ((i<xsize) && (j<ysize))) {
        // pick a random trial vector in box size (2sigma)^2
        trial.setValue((((random.nextFloat()*(float)2.0)-(float)1.0)*sigma),
                       (((random.nextFloat()*(float)2.0)-(float)1.0)*sigma),
                       (((random.nextFloat()*(float)2.0)-(float)1.0)*sigma));
        trial.accAdd(x[i][j]);
        // find local energy without change 
        energy=localGaussian(i,j,x[i][j])+kappa*localBending(i,j,x[i][j]);  
        // find local energy with change     
        energy=energy-(localGaussian(i,j,trial)+kappa*localBending(i,j,trial));  
        // accept move?
        if ( Math.exp((double)energy) > (double)random.nextFloat() ) {
          localUpdate(i,j,trial);
          accepted++;
        }
      }
    }
    measure();
    return(accepted);
  }


  private void localUpdate(int i, int j, Vec3 trial) {
    /* This doesn't worry about i=xsize or j=ysize as we expect
       the na and nb arrays to be oversize */
    int ii, jj;
    x[i][j].setValue(trial);
    if ((ii=i-1)>=0) {
      if ((jj=j-1)>=0) {
        na[ii][jj]=localna[0][0];
        nb[ii][jj]=localnb[0][0];
      }
      na[ii][j]=localna[0][1];
    }
    if ((jj=j-1)>=0) {
      nb[i][jj]=localnb[1][0];
    }
    na[i][j]=localna[1][1];
    nb[i][j]=localnb[1][1];
  }
  
  
  /** Calculates Gaussian stretching energy of the form
      <tt>|x-x'|^2</tt> for all nearest neighbours 
      <tt>x'</tt> around <tt>x</tt> at <tt>(i,j)</tt>.
      Free boundary conditions, ie. no boundary energy. */
  private float localGaussian(int i, int j, Vec3 centrePoint) {
    int ii,jj;
    float energy;
    energy=(float)0.0;
    if ((ii=i-1)>=0) {
      energy=Vec3.distSqd(centrePoint,x[ii][j]);
      if ((jj=j-1)>=0) {
        energy=energy+Vec3.distSqd(centrePoint,x[ii][jj]);
      }
    }
    if ((jj=j-1)>=0) {
      energy=energy+Vec3.distSqd(centrePoint,x[i][jj]);
    }
    if ((jj=j+1)<ysize) {
      energy=energy+Vec3.distSqd(centrePoint,x[i][jj]);
    }
    if ((ii=i+1)<xsize) {
      energy=energy+Vec3.distSqd(centrePoint,x[ii][j]);
      if ((jj=j+1)<ysize) {
        energy=energy+Vec3.distSqd(centrePoint,x[ii][jj]);
      }
    }
    return(energy); 
  }
  
  
  private float localBending(int i, int j, Vec3 centrePoint) {
    int ii,jj;
    float energy;
    energy=(float)0.0;
    /* first calculate normals in local arrays: localna, localnb
       localna[0][0] is just left,below centrePoint
       localna[1][1] is just right,above centrePoint */
    if ((ii=i-1)>=0) {
      if ((jj=j-1)>=0) {
        localnb[0][0].normal(centrePoint,x[ii][j],x[ii][jj]);
        localna[0][0].normal(centrePoint,x[ii][jj],x[i][jj]);
      } else {
        localnb[0][0].setZero();
        localna[0][0].setZero();
      }
      if ((jj=j+1)<ysize) {
        localna[0][1].normal(centrePoint,x[i][jj],x[ii][j]);
      } else {
        localna[0][1].setZero();
      }
    }
    if ((ii=i+1)<xsize) {
      if ((jj=j+1)<ysize) {
        localnb[1][1].normal(centrePoint,x[ii][jj],x[i][jj]);
        localna[1][1].normal(centrePoint,x[ii][j],x[ii][jj]);
      } else {
        localnb[1][1].setZero();
        localna[1][1].setZero();
      }
      if ((jj=j-1)>=0) {
        localnb[1][0].normal(centrePoint,x[i][jj],x[ii][j]);
      } else {  
        localnb[1][0].setZero();
      }
    }
    /* Now use normals to calculate bending energy */
    /* Stage 1: between six triangles around centrePoint. All 
       normals that don't exist are set to zero so we can just 
       calculate blindly. */
    energy=Vec3.dot(localna[0][0],localnb[0][0])+
           Vec3.dot(localnb[0][0],localna[0][1])+
           Vec3.dot(localna[0][1],localnb[1][1])+
           Vec3.dot(localnb[1][1],localna[1][1])+
           Vec3.dot(localna[1][1],localnb[1][0])+
           Vec3.dot(localnb[1][0],localna[0][0]);
    /* Stage 2: outward from six triangles around centrePoint. Here we
       must think about which of the far normals exist. However, if the
       far normal exists then the one for the triangle touching centrePoint 
       must also exist. */
    if (((ii=i-2)>=0) && ((jj=j-1)>=0)) {
      energy=energy+Vec3.dot(na[ii][jj],localnb[0][0]);
    }
    if (((ii=i-1)>=0) && (j<ysize)) {
      energy=energy+Vec3.dot(nb[ii][j],localna[0][1]);
    }
    if ((i<xsize) && ((jj=j+1)<ysize)) {
      energy=energy+Vec3.dot(na[i][jj],localnb[1][1]);
    }
    if (((ii=i+1)<xsize) && (j<ysize)) {
      energy=energy+Vec3.dot(nb[ii][j],localna[1][1]);
    }
    if ((i<xsize) && ((jj=j-1)>=0)) {
      energy=energy+Vec3.dot(na[i][jj],localnb[1][0]);
    }
    if (((ii=i-1)>=0) && ((jj=j-2)>=0)) {
      energy=energy+Vec3.dot(nb[ii][jj],localna[0][0]);
    }   
    return(energy);   
  }
  
  
  public void measure() {
     float energy;
     energy=globalGaussian();
     System.out.println("  "+String.valueOf(energy));
  }

  /** Calculates Gaussian stretching energy of the form
      <tt>|x-x'|^2</tt> for all nearest neighbour pairs 
      <tt>x</tt> and <tt>x'</tt>.
      Free boundary conditions, ie. no boundary energy. */
  private float globalGaussian() {
    int i,j,ii,jj;
    float energy;
    energy=(float)0.0;
    for (i=0;i<xsize;i++) { 
      for (j=0;j<xsize;j++) { 
        if ((ii=i+1)<xsize) {
          energy=energy+Vec3.distSqd(x[i][j],x[ii][j]);
          if ((jj=j+1)<ysize) {
            energy=energy+Vec3.distSqd(x[i][j],x[ii][jj]);
          }
        }
        if ((jj=j+1)<ysize) {
          energy=energy+Vec3.distSqd(x[i][j],x[i][jj]);
        }
      }
    }
    return(energy); 
  }
  
  /** This method initialises a new drawing sequence by calculating
  * screen coordinates for all of the points in the mesh and resetting
  * the count of faces drawn. The coordintates are scaled to fit in
  * the display area <tt>width</tt> pixels by <tt>height</tt> pixels.
  * The method <tt>nextFace()</tt> is used to get each face in turn.
  */
  public void newPaint(int width, int height) {
    int i,j,xoffset, yoffset;
    float a,b,scale;
    halfSize=ysize/(int)2;
    // reset drawing counters
    idraw=xsize-2; jdraw=ysize-2; kdraw=-1;
    // work out how to scale picture
    a=(float)width/((float)(xsize+ysize/2));
    b=(float)height/(float)ysize;
    if (a<b) { scale=a*(float)0.7; } else { scale=b*(float)0.7; }
    xoffset=(width-(int)(scale*(float)(xsize-ysize/2)))/2;
    yoffset=(height+(int)(scale*(float)ysize*(float)0.8))/2;
    // calculate screen coords for all triangulation points
    for (i=0;i<xsize;i++) {
      for (j=0;j<ysize;j++) {
        facesx[i][j]=(int)(scale*(x[i][j].x1))+xoffset;
        facesy[i][j]=yoffset-(int)(scale*(x[i][j].x2*0.8+x[i][j].x3));
      }
    }
  }

  /** After calling the <tt>newPaint(...)</tt> method this method
  is used to update the <tt>face[6]</tt> array with the coordinates
  of the next face to be drawn (three x and three y ordinates). The return value
  is <tt>true</tt> if coordinates are returned, <tt>false</tt>
  otherwise indicating that all faces have been read.
  */
  public boolean nextFace(int[] facex, int[] facey) {
    do {
      kdraw++;
      if (kdraw==2) {
        idraw--;
        if (idraw==0) {
          idraw=xsize-2;
          jdraw--;
          if (jdraw==0) return(false);
        }
        kdraw=0;
      }
    } while ( ((idraw-jdraw)<-halfSize) ||
              (((idraw-jdraw)==-halfSize) && (kdraw==1)) ||
              ((idraw-jdraw)>halfSize) ||
              (((idraw-jdraw)==halfSize) && (kdraw==0)) );
    facex[0]=facesx[idraw][jdraw];
    facey[0]=facesy[idraw][jdraw];
    facex[3]=facex[0];
    facey[3]=facey[0];
    facex[1]=facesx[idraw+1][jdraw+1];
    facey[1]=facesy[idraw+1][jdraw+1];
    if (kdraw==0) {
      // faces b
      facex[2]=facesx[idraw+1][jdraw];
      facey[2]=facesy[idraw+1][jdraw];
    } else {
      // faces a
      facex[2]=facesx[idraw][jdraw+1];
      facey[2]=facesy[idraw][jdraw+1];
    }
    return(true);
  }
}