//
//  Moon.java
//  Generate Huygens Suite Landscape
//
//  Created by Cara MacNish on 1/04/05.
//  Copyright 2005 CSSE, UWA. All rights reserved.
//
import java.io.*;

public class Moon {
  
  private long IA = 16807l;      //training
  private int maxRecDepth = 12;
  private int bmean = 20;
  private long seed = 1;
  
  // Constructors
  
  public Moon (long IA, int maxRecDepth, int bmean, long seed) {
    this.IA = IA;
    this.maxRecDepth = maxRecDepth;
    this.bmean = bmean;
    this.seed = seed;
  }
  
  public Moon (long IA, int bmean, long seed) {
    this.IA = IA;
    this.bmean = bmean;
    this.seed = seed;
  }
  
  public Moon (long IA, long seed) {
    this.IA = IA;
    this.seed = seed;
  }
  
  public Moon (int bmean, long seed) {
    this.bmean = bmean;
    this.seed = seed;
  }
  
  public Moon (long seed) {
    this.seed = seed;
  }
  
  public Moon () {
  }
  
  // Core methods
      
  double getDepth (double[] pos) {
    // map pos into [0,1), [0,1)
    pos[1] = pos[1] % 1;
    pos[2] = pos[2] % 1;
    if (pos[1] < 0) pos[1] = pos[1]+1;
    if (pos[2] < 0) pos[2] = pos[2]+1;
    return getDepthLocal (pos, 0, seed);
  }
 
  double getDepth2 (double[] pos2) {
    double[] pos = {0, pos2[0], pos2[1]};
    return getDepth(pos);
  }
  
  double getDepth (double x, double y) {
    double[] pos = {0, x, y};
    return getDepth(pos);
  }
  
  double getDepthLocal (double[] pos, int recDepth, long seed) {
    double depth = 0;
    long span = (long) Math.pow(10,recDepth);
    long[] bin = {0, (long) Math.ceil(pos[1]*span), (long) Math.ceil(pos[2]*span)};
    
    long newSeed;
    double[] pos1 = {0,0,0};
    long[] square1 = {0,0,0};
    int[] offset = {0,0,0};

    // get contribution from each of the 9 relevant squares...
    for (int i = -1; i<2; i++) {
      for (int j = -1; j<2; j++) {
        System.arraycopy(pos,0,pos1,0,3);
        System.arraycopy(bin,0,square1,0,3);
        offset[1] = i; offset[2] = j; 
        translate(pos1, square1, offset, span);
        double td = getDepthWRTSquare(pos1, square1, recDepth, seed);
        depth = depth + td;  // getDepthWRTSquare(pos1, square1, recDepth, bmean, seed);
      }
    }

    // now fire recursion to next level down...
    if (recDepth < maxRecDepth) {
      // unique seeds are possible to depth 7
      if (recDepth < 8) newSeed = (long) Math.pow(10,2*recDepth) + seed;
      // otherwise settle for any different (repeat) seed
      else newSeed = seed+1;
      depth = depth + getDepthLocal(pos,recDepth+1,newSeed);
    }
    return depth;
  }

  
  void translate (double[] pos, long[] square, int[] offset, long span) {
    for (int i=1; i<3; i++) {
      square[i] = square[i] + offset[i];
      if (square[i] == 0) {
        square[i] = span;
        pos[i] = pos[i]+1;
      }
      else if (square[i] > span) {
        square[i] = 1;
        pos[i] = pos[i]-1;
      }
    }
  }
  
  
  double getDepthWRTSquare (double[] pos, long[] square, int recDepth, long seed) {
    
    long localSeed = 0;
    double bstdev = bmean/5.0;
    double x = pos[1];
    double y = pos[2];
    long span = (long) (Math.pow(10,recDepth));
    double scale = 1.0/span;
    
    if (recDepth > 7)  localSeed = seed + (square[1]-1)*10 + square[2]-1;
    else localSeed = seed + (square[1]-1)*span + (square[2]-1);
    
    Ran0 ran = new Ran0(IA, localSeed);
    
    double gdev = ran.gauss(bmean, bstdev);
    int numboulders = (int) Math.abs(Math.round(gdev));
    
    double depth = 0;
    for (int i=1; i<=numboulders; i++) {
      double diameter = 1/(10 - 9*ran.uniform()) * scale; // get diameter from skewed distribution
      double radius = diameter/2.0;
      double radsq = radius*radius;
      double[] centre = {0, (square[1]-ran.uniform())*scale, (square[2]-ran.uniform())*scale};
      double dist1 = Math.abs(pos[1] - centre[1]);
      double dist2 = Math.abs(pos[2] - centre[2]);
      double distsq = dist1*dist1 + dist2*dist2;
      double distinc = 0;
      if (distsq < radsq) distinc = Math.sqrt(radsq-distsq);
      depth = depth-distinc;
    }
    return depth;
  }
  
  // Utility methods
  
  /**
    * Get depths for an array of probes.
   */
  double[] getDepths (double[][] probes) {
    double[] result = new double[probes.length];
    for (int i=0; i<probes.length; i++)
      result[i] = getDepth2(probes[i]);
    return result;
  }
  
  // array of grid evaluations
  double[][] getDepths (int probesPerAxis) {
    double[][] depths = new double[probesPerAxis][probesPerAxis];
    double step = 1.0/probesPerAxis;
    for (int i = 0; i < probesPerAxis; i++) {
      for (int j = 0; j < probesPerAxis; j++) {
        double[] pos = {0, i*step, j*step};
        depths[i][j] = getDepth(pos);
      }
    }
    return depths;
  }
  
  // minimum of grid over whole surface
  double getMinimum (int probesPerAxis) {
    double minimum =0;
    double step = 1.0/probesPerAxis;
    for (int i = 0; i < probesPerAxis; i++) {
      for (int j = 0; j < probesPerAxis; j++) {
        double[] pos = {0, i*step, j*step};
        double depth = getDepth(pos);
        if (depth < minimum) minimum=depth;
      }
      if (i%100 == 0) System.out.println(i);
    }
    return minimum;
  }
  
  // write grid of depths for plotting graph
  public void writeDepths (String dataDir, int probesPerAxis) throws IOException {
    //    int maxRecDepth = 8;
    double depth;
    File dir = new File(dataDir);
    int totalProbes=probesPerAxis+1;
    String name =  "moon_"+bmean+"_"+seed+"-"+totalProbes+"x"+totalProbes+".txt";
    File target = new File(dir,name);
    target.createNewFile();
    PrintWriter out = new PrintWriter(new BufferedWriter(new FileWriter(target)));
    double step = 1.0/probesPerAxis;
    for (int y = 0; y <= probesPerAxis; y++) {
      for (int x = 0; x <= probesPerAxis; x++) {
        double[] pos = {0, x*step, y*step};
        depth = getDepth(pos);
        out.print(depth+" ");
      }
      out.println();
    }
    out.flush();
    out.close();
  }
 
  /**
    * Write grid of depths within square at resolution determined by first and last. Useful for
   * demonstrating fractal properties.
   */
  public void writeResolutions (String dataDir, int probesPerAxis, double first, double last) throws IOException {
    //  int maxRecDepth = 12;
    double depth;
    File dir = new File(dataDir);
    String name =  "moon_"+bmean+"_"+seed+"_"+first+"_"+last;
    File target = new File(dir,name);
    target.createNewFile();
    PrintWriter out = new PrintWriter(new BufferedWriter(new FileWriter(target)));
    double step = (last-first)/probesPerAxis;
    for (int y = 0; y <= probesPerAxis; y++) {
      for (int x = 0; x <= probesPerAxis; x++) {
        double[] pos = {0, first+x*step, first+y*step};
        depth = getDepth(pos);
        out.print(depth+" ");
      }
      out.println();
    }
    out.flush();
    out.close();
  }
  
  /**
    * Write sequence of depths along line at resolution determined by first and last. Useful for
   * demonstrating fractal properties.
   */
  public void writeResolutions1D (String dataDir, int probesPerAxis, double first, double last) throws IOException {
    //    int maxRecDepth = 12;
    double depth;
    File dir = new File(dataDir);
    String name =  "moon_"+bmean+"_"+seed+"_"+first+"_"+last;
    File target = new File(dir,name);
    target.createNewFile();
    PrintWriter out = new PrintWriter(new BufferedWriter(new FileWriter(target)));
    double step = (last-first)/probesPerAxis;
    for (int x = 0; x <= probesPerAxis; x++) {
      double[] pos = {0, first+x*step, 0.5};
      depth = getDepth(pos);
      out.print(depth+" ");
    }
    out.println();
    System.out.println(name);
    out.flush();
    out.close();
  }
  
}