00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017 #include <grass/Vect.h>
00018
00029 int
00030 Vect_tin_get_z (
00031 struct Map_info *Map,
00032 double tx, double ty,
00033 double *tz, double *angle, double *slope)
00034 {
00035 int i, area, n_points;
00036 struct Plus_head *Plus;
00037 P_AREA *Area;
00038 static struct line_pnts *Points;
00039 static int first_time = 1;
00040 double *x, *y, *z;
00041 double vx1, vx2, vy1, vy2, vz1, vz2;
00042 double a, b, c, d;
00043
00044
00045
00046 Plus = &(Map->plus);
00047 if (first_time == 1) {
00048 Points = Vect_new_line_struct ();
00049 first_time = 0;
00050 }
00051
00052 area = Vect_find_area ( Map, tx, ty );
00053 G_debug (3, "area = %d", area );
00054 if ( area == 0 ) return 0;
00055
00056 Area = Plus->Area[area];
00057 if ( Area->n_isles > 0 ) return -1;
00058
00059 Vect_get_area_points (Map, area, Points);
00060 n_points = Points->n_points;
00061 if ( n_points != 4 ) return -1;
00062
00063 x = Points->x;
00064 y = Points->y;
00065 z = Points->z;
00066 for (i = 0; i < 3; i++) {
00067 G_debug(3, "%d %f %f %f", i, x[i], y[i], z[i] );
00068 }
00069
00070 vx1 = x[1] - x[0];
00071 vy1 = y[1] - y[0];
00072 vz1 = z[1] - z[0];
00073 vx2 = x[2] - x[0];
00074 vy2 = y[2] - y[0];
00075 vz2 = z[2] - z[0];
00076
00077 a = vy1*vz2 - vy2*vz1;
00078 b = vz1*vx2 - vz2*vx1;
00079 c = vx1*vy2 - vx2*vy1;
00080 d = -a*x[0] - b*y[0] - c*z[0];
00081
00082
00083 *tz = -(d + a*tx + b*ty) / c;
00084 G_debug(3, "z = %f", *tz );
00085
00086 return 1;
00087 }
00088
00089