| Triangle | \n" + " | Side 1 | \n" + " | Side 2 | \n" + " | Side 3 | \n" + " | Area\n\n"; public static void main(String args[]) throws IOException { String bd = "/home/rob/projects/projects_int/landshare/purchase/Itamar/survey/triangulator.htm"; earthRadius(); //compute the Earth's radius + mean height at the site sidesAndAreas(); //compute lengths of sides of all triangles tf = new FileWriter(bd); //open output file double A = 0; tf.write(TH); for(int i = 0; i < triangles.length; i++) { //print side lengths and area of each triangle A += areas[i][3]; tf.write(" |
|---|---|---|---|---|---|---|---|---|
| " + L[i] + " | \n" + " | " + twoDec(areas[i][0]) + " | \n" + " | " + twoDec(areas[i][1]) + " | \n" + " | " + twoDec(areas[i][2]) + " | \n" + " | " + twoDec(areas[i][3]) + "\n\n"); } tf.write(" |
Total Area = " + twoDec(A) + " square metres.\n");
double b = A - 5000.30; // error in square metres
String s = "+"; if(b < 0) {s = "-"; b = -b;}
tf.write("
Error: " + s + twoDec(b) + " square metres.\n");
tf.write("
Equivalent to a square of side " + twoDec(Math.sqrt(b)) + " metres.\n");
tf.write("
Earth radius at mean site latitude = " + twoDec(R) + " metres.\n");
tf.write("
This radius includes the mean elevation\n");
tf.write("
of the site above sea level: " + twoDec(H) + " metres.\n\n");
tf.write("
Within this site and its viscinity:");
tf.write("
0.01 second of latitude = " + twoDec(Δφ * 100) + " cm.\n");
tf.write("
0.01 second of longitude = " + twoDec(Δλ * 100) + " cm.\n\n");
tf.close();
}
/* computes the side-lengths and area of each triangle of the triangulated plot */
private static void sidesAndAreas(){ //side-lengths and areas of triangles
for(int i = 0; i < 6; i++) { //for each of the triangles of plot
toRads(triangles[i][0]);
double φ1 = φ; //latitude of first point in radians
double λ1 = λ; //longitude of first point in radians
toRads(triangles[i][1]);
double φ2 = φ; //latitude of second point in radians
double λ2 = λ; //longitude of second point in radians
toRads(triangles[i][2]);
double φ3 = φ; //latitude of third point in radians
double λ3 = λ; //longitude of third point in radians
double a = side(φ1, λ1, φ2, λ2); //length of first side
double b = side(φ2, λ2, φ3, λ3); //length of second side
double c = side(φ3, λ3, φ1, λ1); //length of third side
double A = heron(a, b, c); //area of the triangle
areas[i][0] = a; //store side lengths and area of this triangle in array
areas[i][1] = b;
areas[i][2] = c;
areas[i][3] = A;
}
}
/* accesses the lattitude and longitude of a point n in the latlng[][] array
(in degrees:minutes:seconds), converts them to radians, stores
them respectively in the doubles φ (latitude) and λ (longitude) */
private static void toRads(int n) { //n = index number of point
φ = Math.toRadians(latlng[n][0]
+ latlng[n][1] / 60 //latitude of point n in radians
+ latlng[n][2] / 3600);
λ = Math.toRadians(latlng[n][3]
+ latlng[n][4] / 60 //longitude of point n in radians
+ latlng[n][5] / 3600);
}
/* This formula is only an approximation when applied to the Earth, because the
Earth is not a perfect sphere. There are small errors, typically of the order
of 0.1% (assuming the geometric mean radius is used, because of this slight
ellipticity of the planet. A more accurate method, which takes into account
the Earth's ellipticity, is given by Vincenty's formulae. */
private static double side(double φ1, double λ1, double φ2, double λ2){
return R * invhav(havsin(φ1 - φ2)
+ Math.cos(φ1) * Math.cos(φ2) * havsin(λ1 - λ2));
}
private static double havsin(double θ){ //computes the haversine of θ
double x = Math.sin(θ/2);
return x * x;
}
private static double invhav(double θ){ // computes the inverse haversine of θ
return 2 * Math.asin(Math.sqrt(θ));
}
/* The method invented by Heron (or Hero) of Alexandria to compute the
area of a general plane triangle from the lengths of its sides */
private static double heron(double a, double b, double c){
double s = (a + b + c) / 2; //semi-circumference of the triangle
return Math.sqrt(s * (s - a) * (s - b) * (s - c));
}
/* called only once at the beginning to establish the radius
of the Earth at the mid-latitude of the plot */
private static void earthRadius(){
double re = 6378137.000, //equitorial radius using GRS80 / WGS84 (NAD83) ellipsoid
rp = 6356752.314; //polar radius 6356752.31414034724399572718485825958879
φm = 0;
for(int i = 0; i < latlng.length; i++) {
toRads(i);
φm += φ; //sum the latitudes of all the points that mark the site
}
φm /= latlng.length; //mean latitude of the site
double sd = Math.sin(φm), //sine of the mean latitude of the site
re2 = re * re,
rp2 = rp * rp, //squared values
sd2 = sd * sd;
for(int i = 0; i < GH.length; i++)
H += GH[i]; //sum the heights of all the points that mark the site
H /= GH.length; //mean height of the site
R = H + re * Math.sqrt(1 + sd2 * (rp2 - re2) / re2);
// Measurement precision at this site:
Δφ = 2* π * R / (360 * 60 * 60 * 100); //0.01 second of latitude in metres
Δλ = Δφ * Math.cos(φm); //0.01 second of longitude in metres
}
/* accepts a double and returns it as a string
truncated or padded to two decimal places */
private static String twoDec(double d){
d += .005; //round to 2 decimal places
String s = "" + d;
int x = s.indexOf('.');
int y = s.length();
if(x == -1)
s += ".00";
else if(x == (y - 1))
s += "00";
else if(x == (y - 2))
s += "0";
else
s = s.substring(0, x + 3);
return s;
}
}