// conversion of Steve's quick and dirty ruby approximation of Mercator projection from www.openstreetmap.org

public class OSMMercator {

  protected double clat, clon, dlat, dlon, degrees_per_pixel;
  protected double tx, ty, bx, by;
  protected double w, h;
  protected double QUARTER_PI = Math.PI / 4.0;
  protected double PIby360 = Math.PI / 360.0;

  protected double xdiv, ydiv;

  /** init me with your centre lat/lon, the number of degrees per pixel and the size of your image */
  public OSMMercator(double clat, double clon, double degrees_per_pixel, int w, int h) {
  
    this.clat = clat;
    this.clon = clon;
    this.degrees_per_pixel = degrees_per_pixel;
    this.w = w;
    this.h = h;
    dlon = this.w / 2.0 * degrees_per_pixel;
    dlat = this.h / 2.0 * degrees_per_pixel * Math.cos(clat * Math.PI / 180.0);

    tx = xsheet(clon - dlon);
    ty = ysheet(clat - dlat);

    bx = xsheet(clon + dlon);
    by = ysheet(clat + dlat);

    xdiv = 1.0 / (bx - tx) * w;
    ydiv = 1.0 / (by - ty) * h;

  }

  public double kilometerinpixels() {
    return 40008.0  / 360.0 * degrees_per_pixel;
  }

  // the following two functions will give you the x/y on the entire sheet

  public double ysheet(double lat) {
    return Math.log(Math.tan(QUARTER_PI + (lat * PIby360)));
  }
  
  public double xsheet(double lon) {
    return lon;
  }

  // and these two will give you the right points on your image. all the constants can be reduced to speed things up. FIXME

  public double y(double lat) {
    return h - ((ysheet(lat) - ty) * ydiv);
  }

  public double x(double lon) {
    return (xsheet(lon) - tx) * xdiv;
  }

  public boolean projectable(double lat, double lon) {
    // top left?                                     bottom right?
    return lat < clat + dlat && lon > clon - dlon && lat > clat - dlat && lon < clon + dlon;
  }

}