/*++++++++++++++ .IDENTIFICATION htm1.c .LANGUAGE C .AUTHOR Francois Ochsenbein [CDS] .ENVIRONMENT .KEYWORDS .VERSION 1.0 20-Jul-2001 .COMMENTS Direct implementation of HTM -- small, fast ! Max HTM level (with 32 bits integers): 14 ---------------*/ #include #include #ifndef HTM_LEVEL #define HTM_LEVEL 10 #endif #define acosd(x) (180./M_PI)*acos(x) #define asind(x) (180./M_PI)*asin(x) #define atand(x) (180./M_PI)*atan(x) #define atan2d(x,y) (180./M_PI)*atan2(x,y) #define cosd(x) cos((x)*(M_PI/180.)) #define sind(x) sin((x)*(M_PI/180.)) /* Computation uncertainties: to check limits, accept some extra */ #define EPSILON 5.e-16 /* The 6 root corners of the octaedron numbers 0 .. 5: z x y -x -y -z */ static double HTM6[6][3] = { /* z x y */ 0,0,1, 1,0,0, 0,1,0, -1,0,0, 0,-1,0, 0,0,-1 /* -x -y -z */ }; /* The basic points, at level 0 */ static double *HTM0[] = { HTM6[1], HTM6[5], HTM6[2], /* S0 */ HTM6[2], HTM6[5], HTM6[3], /* S1 */ HTM6[3], HTM6[5], HTM6[4], /* S2 */ HTM6[4], HTM6[5], HTM6[1], /* S3 */ HTM6[1], HTM6[0], HTM6[4], /* N0 */ HTM6[4], HTM6[0], HTM6[3], /* N1 */ HTM6[3], HTM6[0], HTM6[2], /* N2 */ HTM6[2], HTM6[0], HTM6[1], /* N3 */ } ; /* Variables set and retrieved by the actions (callback routines) */ static double NaN ; static int theCount; /* Number of terminal triangles */ static int theNum; static int theLevel ; static int theLevelMask = 0xc << (HTM_LEVEL<<1) ; static int (*theAction)() ; static double *theCenter ; /* x, y, z of Center */ static double *theLimsRA ; /* cos(RA1) sin(RA1) cos(RA2) sin(RA2) */ static double *theLimsDec ; /* cos(DE1) sin(DE1) cos(DE2) sin(DE2) */ static double *theLimsRad ; /* cos(rho) sin(rho) sin2(r/2) */ static double *theTriangleValues ; /*======================================================================= Basic vector operators *=======================================================================*/ #define dotprod(v1, v2) (v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2]) #define norm2(v) dotprod(v,v) static void vecprod (double *v1, double *v2, double *r) /*++++++++++++++++ .PURPOSE Compute vectorial product ( v1 ^ v2) .RETURNS --- .REMARKS Result in r -----------------*/ { r[0] = v1[1]*v2[2] - v1[2]*v2[1]; r[1] = v1[2]*v2[0] - v1[0]*v2[2]; r[2] = v1[0]*v2[1] - v1[1]*v2[0]; } static double det (double **v) /*++++++++++++++++ .PURPOSE Compute determinant .RETURNS The value ( = (v[0] ^ v[1]) . v[2]) -----------------*/ { #if 1 double vec[3] ; vecprod(v[0], v[1], vec) ; return(dotprod(vec, v[2])); #else return(v[0][0]*v[1][1]*v[2][2] + v[0][1]*v[1][2]*v[2][0] + v[0][2]*v[1][0]*v[2][1] - v[0][2]*v[1][1]*v[2][0] - v[0][1]*v[1][0]*v[2][2] - v[0][0]*v[1][2]*v[2][1] ) ; #endif } /*======================================================================= Basic HTM functions *=======================================================================*/ int htm_level(int htm) /*++++++++++++++++ .PURPOSE Find out the level from an HTM node .RETURNS -2/ -1(N or S) / 0 (N0 --> S3) /etc .REMARKS Leftmost bit set to 1 -----------------*/ { int n, level ; if ((htm&(~7)) == 0) { if ((htm == 2) || (htm == 3)) level = -1 ; else level = -2 ; } else { n = (htm>>4) & 0xfffffff ; for (level=0; n; n >>= 2) level++ ; } return(level) ; } int a2htm(char *text) /*++++++++++++++++ .PURPOSE Interpret the HTM string .RETURNS The value / 1 = invalid HTM spec .REMARKS Empty string valid (= emptyset) -----------------*/ { int n,c ; char *p ; n = toupper(*text) ; if (n == 'S') n = 2 ; else if (n == 'N') n = 3; else if (!n) return(0) ; /* Emptyset */ else return(1) ; /* Invalid */ for (p=text+1; *p; p++) { c = *p - '0' ; if (c&(~3)) return(1) ; n = (n<<2) | c ; } return(n) ; } char *htm2a(int htm) /*++++++++++++++++ .PURPOSE Edit (in internal static buffer) the htm number .RETURNS The internal buffer .REMARKS First byte of result is 'N' or 'S' -- other indicate an error -----------------*/ { static char hemisphere[] = {' ', '?', 'S', 'N' } ; static char buf[20] ; char *p ; int n ; /* Edit groups of 2 bits from right to left. */ p = buf + sizeof(buf) ; *--p = 0 ; for (n=htm; n; n >>= 2) *--p = '0' + (n&3) ; /* The leftmost 2-bits are 2 = South, 3 = North */ if (*p) *p = hemisphere[*p - '0'] ; return(p) ; } int htm_get() /*++++++++++++++++ .PURPOSE Get the deepest level .RETURNS The currently defined cell level -----------------*/ { return(htm_level(theLevelMask)); } int htm_set(int level) /*++++++++++++++++ .PURPOSE Set the deepest level .RETURNS The previously defined cell level / -1 -----------------*/ { int old_level, mask ; old_level = htm_level(theLevelMask) ; mask = 0xc << (level << 1) ; if (!mask) { fprintf(stderr, "****Invalid HTM_level: %s (ignored)\n", level) ; return(-1) ; } theLevelMask = mask ; return(old_level) ; } double htm_triangle(double *triangle[3], double u_center[3]) /*++++++++++++++++ .PURPOSE Find the Center and cos(Radius) corresponding to a HTM triangle .RETURNS cos(Radius) .REMARKS The center of the triangle (u1 u2 u3) is u1^u2 + u2^u3 + u3^u1 and the radius is: cos(r) = det(u1 u2 u3) / ||u1^u2 + u2^u3 + u3^u1|| -----------------*/ { double u[3], n ; /* The center is the vector v1^v2 + v2^v3 + v3^v1 */ vecprod(triangle[0], triangle[1], u_center) ; vecprod(triangle[1], triangle[2], u) ; u_center[0] += u[0] ; u_center[1] += u[1] ; u_center[2] += u[2] ; vecprod(triangle[2], triangle[0], u) ; u_center[0] += u[0] ; u_center[1] += u[1] ; u_center[2] += u[2] ; n = sqrt(norm2(u_center)) ; u_center[0] /= n; u_center[1] /= n; u_center[2] /= n; return(det(triangle)/n) ; } /*======================================================================= Internal Routines *=======================================================================*/ static int descent(int htm, double *triangle[3], int (*action)()) /*++++++++++++++++ .PURPOSE Recursive function to go down in the HTM tree. It works as follows: Given an HTM triangle corresponding to the htm number, the action (htm, triangle, 1) is called (examples are check_* below); the action decides (by returning 1 or 0) whether the triangle is worth to be examined (then move down in the HTM tree), or not. At the deepest level (specified by htm_set, 10 by default), theAction is called. .RETURNS 0 .REMARKS theLevelMask (0xc<<(2*level)) must be set !! -----------------*/ { double mid[3][3] ; /* Mid-points */ double *tr4[12] ; /* New triangles */ double n, **a4 ; /* 4-triangles */ int i ; #ifdef DEBUG static int indent ; for (i=indent; --i >= 0; ) putchar(' ') ; printf("....descent(htm=%x)\n", htm) ; #endif if (htm & theLevelMask) { /* At Deepest Level ... */ if ((*action)(htm, triangle, 1)) { if (theAction) (*theAction)(htm, triangle, 1) ; theCount++ ; } return(0) ; } #ifdef DEBUG indent++ ; #endif if ((htm&(~3)) == 0) { /* Initial (htm<4), use HTM0 */ if (htm == 2) /* South */ a4 = HTM0 ; else if (htm == 3) /* Norh */ a4 = HTM0 + 12 ; else { /* Both! */ descent(2, triangle, action) ; descent(3, triangle, action) ; return(0) ; } } else { /* Compute first the halves */ for (i=0; i<3; i++) { mid[0][i] = 0.5*(triangle[1][i] + triangle[2][i]) ; mid[1][i] = 0.5*(triangle[2][i] + triangle[0][i]) ; mid[2][i] = 0.5*(triangle[0][i] + triangle[1][i]) ; } /* ... norm the new vectors */ for (i=0; i<3; i++) { n = sqrt(norm2(mid[i])) ; mid[i][0] /= n ; mid[i][1] /= n ; mid[i][2] /= n ; } /* ... and save the 4-triangles (in trigo order) */ tr4[0] = triangle[0]; tr4[1] = mid[2]; tr4[2] = mid[1]; tr4[3] = triangle[1]; tr4[4] = mid[0]; tr4[5] = mid[2]; tr4[6] = triangle[2]; tr4[7] = mid[1]; tr4[8] = mid[0]; tr4[9] = mid[0]; tr4[10]= mid[1]; tr4[11]= mid[2]; a4 = tr4 ; } htm <<= 2 ; /* Examine the 4 triangles ... */ for (i=4; --i >= 0; htm++, a4 += 3) { if ((*action)(htm, a4, 1)) descent(htm, a4, action) ; } #ifdef DEBUG indent-- ; #endif return(0) ; } static int check_pos(int htm, double *triangle[3], int status) /*++++++++++++++++ .PURPOSE Action Routine to check theCenter is within triangle .RETURNS 0 (no) / 1 (yes) .REMARKS Point in triangle (correctly oriented) when all determinants of a triangle where one edge is replaced by the center are non-negative. -----------------*/ { double *apos ; int i ; if (!status) return(status) ; if (theCount) return(0); /* Just 1 point! */ status = 1 ; for (i=0; (i<3) && status; i++) { apos = triangle[i] ; /* Save vector */ triangle[i] = theCenter ; if (det(triangle) < 0) status = 0 ; /* Point outside */ triangle[i] = apos ; /* Restore vector*/ } if (status) theNum = htm; /* Matching HTM */ return(status) ; } static int check_htm(int htm, double *triangle[3], int status) /*++++++++++++++++ .PURPOSE Action Routine to check whether the HTM number corresponds to theNum .RETURNS 0 (no) / 1 (yes) -----------------*/ { double *aval ; int i ; if (!status) return(status) ; if (theCount) return(0); /* Just 1 point! */ if (htm == theNum) { /* The tree END */ if (theTriangleValues) { memcpy(theTriangleValues, triangle[0], 3*sizeof(double)) ; memcpy(theTriangleValues+3, triangle[1], 3*sizeof(double)) ; memcpy(theTriangleValues+6, triangle[2], 3*sizeof(double)) ; } return(1) ; } i = theLevel - htm_level(htm) ; return ((theNum >> (i<<1)) == htm) ; } static int check_ra(int htm, double *triangle[3], int status) /*++++++++++++++++ .PURPOSE Action Routine to check between RA limits .RETURNS 0 (no) / 1 (yes) .REMARKS Condition: either any of the 3 vertices of the triangle satisfy to the RA condition or at least one of the edges of the triangle intersects one of the RA -----------------*/ { int i, j, n ; double x ; if (!status) return(status) ; #ifdef DEBUG printf("....%-12s", htm2a(htm)); printf("%20.16f%20.16f%20.16f\n", triangle[0][0], triangle[1][0], triangle[2][0]) ; printf(" %-12s", " "); printf("%20.16f%20.16f%20.16f\n", triangle[0][1], triangle[1][1], triangle[2][1]) ; printf(" %-12s", " "); printf("%20.16f%20.16f%20.16f\n", triangle[0][2], triangle[1][2], triangle[2][2]) ; #endif /* Condition 1: is any of the triangle vertices inside the domain ? */ for (i=n=0; (i<3) && (n==0); i++) { n = theLimsRA[0]*triangle[i][1]-theLimsRA[1]*triangle[i][0] >= -EPSILON; if (n) n = theLimsRA[2]*triangle[i][1]-theLimsRA[3]*triangle[i][0] <= EPSILON; } if (n) return(n) ; /* Condition 2: does any of triangle edges intersect the RA lines ? */ for (i=0; (i<3) && (n==0); i++) { j = (i+1)%3 ; /* The next vertex */ x = theLimsRA[0]*triangle[i][1]-theLimsRA[1]*triangle[i][0] ; n = x*(triangle[j][0]*triangle[i][1]-triangle[j][1]*triangle[i][0]) >= -EPSILON ; if (n) n = x*(theLimsRA[0]*triangle[j][1]-theLimsRA[1]*triangle[j][0]) <= EPSILON; } return(n) ; } static int check_dec(int htm, double *triangle[3], int status) /*++++++++++++++++ .PURPOSE Action Routine to check between Dec Limits .RETURNS 0 (no) / 1 (yes) .REMARKS Condition: projections onto z axis -----------------*/ { if (!status) return(status) ; /* z-Projections below minimum ? */ if ( (triangle[0][2] < theLimsDec[1] - EPSILON) && (triangle[1][2] < theLimsDec[1] - EPSILON) && (triangle[2][2] < theLimsDec[1] - EPSILON)) return(0) ; /* z-Projections above maximum ? */ if ( (triangle[0][2] > theLimsDec[1] + EPSILON) && (triangle[1][2] > theLimsDec[1] + EPSILON) && (triangle[2][2] > theLimsDec[1] + EPSILON)) return(0) ; return(1) ; } static int check_radec(int htm, double *triangle[3], int status) /*++++++++++++++++ .PURPOSE Action Routine to check between Dec Limits AND RA lims .RETURNS 0 (no) / 1 (yes) .REMARKS Condition: Distance to center of circle sourronding the triangle. -----------------*/ { if (!status) return(status) ; #ifdef DEBUG printf("....%-12s", htm2a(htm)); printf("%20.16f%20.16f%20.16f\n", triangle[0][0], triangle[1][0], triangle[2][0]) ; printf(" %-12s", " "); printf("%20.16f%20.16f%20.16f\n", triangle[0][1], triangle[1][1], triangle[2][1]) ; printf(" %-12s", " "); printf("%20.16f%20.16f%20.16f\n", triangle[0][2], triangle[1][2], triangle[2][2]) ; #endif if (status = check_ra(htm, triangle, 1)) status = check_dec(htm, triangle, 1); #ifdef DEBUG printf("....%-12s => %d\n", " ", status) ; #endif return(status) ; } static int check_circle(int htm, double *triangle[3], int status) /*++++++++++++++++ .PURPOSE Action Routine to check within a circle .RETURNS 0 (no) / 1 (yes) .REMARKS Condition: Check the triangle is totally external to the circle. The closest distance between the center (C), and the great circle (AB) is: Let C.A = cos(A) C.B = cos(B) A.B = cos(c) there is a point between A and B closer to the center if: cos(c) < (cos(B)/cos(A)) or (cos(A)/cos(B)) and the closest distance is cos(d) = sqrt(cos^2(A) + cos^2(B) - 2 cos(A)cos(B)cos(c))/sin(c) -----------------*/ { double Ccos[3], Tcos[3]; /* cos(r) from Center, in Triangle */ double cos2r; int i0, i1, i2 ; if (!status) return(status) ; /* OK when the center of the circle belongs to triangle... */ if (check_pos(htm, triangle, 1)) return(1) ; Ccos[0] = dotprod(triangle[0], theCenter) ; Ccos[1] = dotprod(triangle[1], theCenter) ; Ccos[2] = dotprod(triangle[2], theCenter) ; if (Ccos[0] >= theLimsRad[0] - EPSILON) return(1) ; if (Ccos[1] >= theLimsRad[0] - EPSILON) return(1) ; if (Ccos[2] >= theLimsRad[0] - EPSILON) return(1) ; /* All vertices outside the circle */ Tcos[0] = dotprod(triangle[1], triangle[2]) ; Tcos[1] = dotprod(triangle[2], triangle[0]) ; Tcos[2] = dotprod(triangle[0], triangle[1]) ; cos2r = theLimsRad[0] * theLimsRad[0] ; for (i0 = 0; i0 < 3; i0++) { i1 = (i0+1)%3; i2 = (i0+2)%3; if ((Ccos[i1]-Ccos[i2]*Tcos[i0])*(Ccos[i2]-Ccos[i1]*Tcos[i0]) <= 0) continue ; /* No closer point between A and B */ if ( (Ccos[i1]*Ccos[i1] + Ccos[i2]*Ccos[i2] - 2.*Ccos[i1]*Ccos[i2]*Tcos[i0]) >= (cos2r * (1. - Tcos[i0]*Tcos[i0]) - EPSILON)) return(1) ; } return(0) ; } /*======================================================================= Position <--> HTM number *=======================================================================*/ int htm_u(double u[3]) /*++++++++++++++++ .PURPOSE Find the HTM cell number from a unit vector .RETURNS The HTM number containing the position. -----------------*/ { double *saved_center ; int saved_num, htm ; /* Save the static values used by action ... */ saved_num = theNum ; saved_center = theCenter ; theNum = 0 ; theCenter = u ; /* Compute directly at level 0 (South / North) */ if (u[2]<=0) htm = 2; else htm = 3; theCount = 0 ; descent(htm, (double **)0, check_pos) ; /* ... Restore static values used by action */ htm = theNum; theNum = saved_num ; theCenter = saved_center ; return(htm) ; } double htm_info(int htm, double triangle[3][3], double u_center[3]) /*++++++++++++++++ .PURPOSE Find the characteristics of a HTM cell .RETURNS cos(Radius) .REMARKS The center of the triangle (u1 u2 u3) is u1^u2 + u2^u3 + u3^u1 and the radius is: cos(r) = det(u1 u2 u3) / ||u1^u2 + u2^u3 + u3^u1|| -----------------*/ { double u[3], radius, *saved_triangles, *atriangle[3]; int i, saved_num, saved_level_mask, saved_level ; /* Save the static values used by the action routine */ saved_level_mask = theLevelMask; saved_level = theLevel ; saved_num = theNum ; saved_triangles = theTriangleValues ; theNum = htm ; theLevel = htm_level(htm) ; atriangle[0] = triangle[0]; atriangle[1] = triangle[1]; atriangle[2] = triangle[2]; theTriangleValues = atriangle[0] ; if (htm_set(theLevel) < 0) { *(int *)(&NaN) = -1 ; /* The NaN value */ return(NaN) ; } /* The initial HTM (at level -1) is 2 or 3. Use as action the function check_htm which verifies the descent. */ i = (htm>>2) & 0xfffffff; i >>= (theLevel<<1) ; theCount = 0 ; descent(i, (double **)0, check_htm) ; /* From the saved triangle, compute the Center and Radius */ /* The center is the vector v1^v2 + v2^v3 + v3^v1 */ radius = htm_triangle(atriangle, u_center) ; /* ... Restore the used static parameters to their original values */ theLevelMask = saved_level_mask ; theNum = saved_num ; theLevel = saved_level ; theTriangleValues = saved_triangles ; return(radius) ; } int htm_circle(double u_center[3], double radius, int (*action)()) /*++++++++++++++++ .PURPOSE Find the triangles affected by a cicular target .RETURNS Number of leaf HTM triangles .REMARKS radius in DEGREES. -----------------*/ { int (*saved_action)() ; double *saved_center, *saved_lims, cosinr[2] ; saved_action = theAction; saved_center = theCenter; saved_lims = theLimsRad; /* Install the static variables... */ theAction = action ; theCenter = u_center ; cosinr[0] = cosd(radius); cosinr[1] = sind(radius) ; theLimsRad = cosinr; theCount = 0 ; descent(0, (double **)0, check_circle) ; /* Restore the original static values */ theAction = saved_action ; theCenter = saved_center ; theLimsRad = saved_lims ; return(theCount) ; } int htm_radec(double alim[2], double dlim[2], int (*action)()) /*++++++++++++++++ .PURPOSE Find the triangles affected by a condition on RA + Dec .RETURNS .REMARKS radius in DEGREES. -----------------*/ { int (*saved_action)(), i ; double *saved_lims[2], cosina[4], cosind[4] ; saved_action = theAction; saved_lims[0]= theLimsRA; saved_lims[1]= theLimsDec; /* Install the static variables... */ theAction = action ; cosina[0] = cosd(alim[0]) ; cosina[1] = sind(alim[0]) ; cosina[2] = cosd(alim[1]) ; cosina[3] = sind(alim[1]) ; i = 0; if (dlim[0] > dlim[1]) i = 1 ; cosind[0] = cosd(dlim[i]) ; cosind[1] = sind(dlim[i]) ; i ^= 1 ; cosind[2] = cosd(dlim[i]) ; cosind[3] = sind(dlim[i]) ; theLimsRA = cosina; theLimsDec= cosind; theCount = 0 ; descent(0, (double **)0, check_radec) ; /* Restore the original static values */ theAction = saved_action ; theLimsRA = saved_lims[0] ; theLimsDec= saved_lims[1] ; return(theCount) ; } /*======================================================================= Test Routines *=======================================================================*/ #ifdef TEST int tr_uo (double u[3], double o[2]) /*++++++++++++++++ .PURPOSE Convert dir. cos. into angles .RETURNS 1 / 0 (x=y=z=0) -----------------*/ { double r2 ; r2 = u[0]*u[0] + u[1]*u[1]; if (r2 == 0.0) { /* In case of poles... */ if (u[2] == 0.0) return(0); o[0] = 0 ; o[1] = u[2]>0.0 ? 90. : -90. ; return(1) ; } o[1] = atand ( u[2] / sqrt(r2)); o[0] = atan2d (u[1] , u[0] ); if (o[0] < 0.0) o[0] += 360.0; return(1) ; } #if 0 FILE *frad ; static int compute_radius(int htm, double *triangle[3], int status) /*++++++++++++++++ .PURPOSE Action Routine to Compute radius & area .RETURNS 1 (yes) -----------------*/ { static int n0 ; double u[3], u_center[3], radius, surf_u3() ; if (!status) return(status) ; /* Only once at level 0 */ if((htm&(~15)) == 0) { if (n0) return(0) ; n0++ ; } /* The center is the vector v1^v2 + v2^v3 + v3^v1 */ radius = acosd(htm_triangle(triangle, u_center)) ; fprintf(frad, " %-12s", htm2a(htm)) ; fprintf(frad, " %12.8f", radius) ; fprintf(frad, " %14.8f\n", surf_u3(triangle[0], triangle[1], triangle[2])) ; return (1) ; } #endif int print_htm(int htm, double *triangle[3], int status) /*++++++++++++++++ .PURPOSE Action used in test program: just print out the HTM. .RETURNS 1 -----------------*/ { double o[2], radius, u_center[3] ; radius = acosd(htm_triangle(triangle, u_center)) ; tr_uo(u_center, o) ; printf("====%-12s st=%d c=%012.8f%+012.8f r=%.8f\n", htm2a(htm), status, o[0], o[1], radius) ; tr_uo(triangle[0], o) ; printf(" %-12s Edge#0=%012.8f%+012.8f # cosr=%16.14f\n", "", o[0], o[1], dotprod(u_center, triangle[0])) ; tr_uo(triangle[1], o) ; printf(" %-12s Edge#1=%012.8f%+012.8f # cosr=%16.14f\n", "", o[0], o[1], dotprod(u_center, triangle[1])) ; tr_uo(triangle[2], o) ; printf(" %-12s Edge#2=%012.8f%+012.8f # cosr=%16.14f\n", "", o[0], o[1], dotprod(u_center, triangle[2])) ; return(1) ; } /* TEST ROUTINE */ main(int argc, char **argv) { char buf[BUFSIZ], *p ; double triangle[3][3], u[3], o[2], ra[2], de[2], atof(), radius ; int htm ; printf("----Set the level as the main option\n") ; ++argv ; if (*argv && isdigit(argv[0][0])) htm_set(atoi(*argv)); #if 0 /* Compute all radiuses */ frad = fopen("/more/francois/gsc2/HTM", "w") ; /* frad = stdout ; */ descent(0, (double **)0, compute_radius) ; fclose(frad) ; #endif while(1) { printf("Give a position: ") ; if (!gets(buf)) break ; for (p=buf; isspace(*p); p++) ; o[0] = atof(buf) ; while (isdigit(*p) || (*p == '.')) p++ ; o[1] = atof(p) ; u[0] = u[1] = cosd(o[1]) ; u[2] = sind(o[1]) ; u[0] *= cosd(o[0]) ; u[1] *= sind(o[0]) ; htm = htm_u(u) ; printf("%010.6f%+010.6f --> %s\n", o[0], o[1], htm2a(htm)) ; /* Convert back center, radius */ radius = acosd(htm_info(htm, triangle, u)) ; tr_uo(u, o) ; printf("%010.6f%+010.6f <-- center, r=%10.6f\n", o[0], o[1], radius) ; tr_uo(triangle[0], o) ; printf("%010.6f%+010.6f <-- Edge #0\n", o[0], o[1]) ; tr_uo(triangle[1], o) ; printf("%010.6f%+010.6f <-- Edge #1\n", o[0], o[1]) ; tr_uo(triangle[2], o) ; printf("%010.6f%+010.6f <-- Edge #2\n", o[0], o[1]) ; /* Check within a small circle */ printf("Give a radius(d) ") ; if (!gets(buf)) break ; radius = atof(buf) ; htm_circle(u, radius, print_htm) ; /* Check now in (RA, DEC) limits */ printf("\nGive Limits in RA (deg): ") ; if (!gets(buf)) break ; for (p=buf; isspace(*p); p++) ; ra[0] = atof(buf) ; while (isdigit(*p) || (*p == '.')) p++ ; ra[1] = atof(p) ; printf("Give Limits in Dec(deg): ") ; if (!gets(buf)) break ; for (p=buf; isspace(*p); p++) ; de[0] = atof(buf) ; if ((*p == '+') || (*p == '-')) p++ ; while (isdigit(*p) || (*p == '.')) p++ ; de[1] = atof(p) ; htm_radec(ra, de, print_htm) ; } } #endif