/Proj4/geod.c

http://github.com/route-me/route-me · C · 240 lines · 227 code · 9 blank · 4 comment · 74 complexity · f226012d902f38ac71af486e3bc639ef MD5 · raw file

  1. #ifndef lint
  2. static const char SCCSID[]="@(#)geod.c 4.8 95/09/23 GIE REL";
  3. #endif
  4. /* <<<< Geodesic filter program >>>> */
  5. # include "projects.h"
  6. # include "geodesic.h"
  7. # include "emess.h"
  8. # include <ctype.h>
  9. # include <stdio.h>
  10. # include <string.h>
  11. # define MAXLINE 200
  12. # define MAX_PARGS 50
  13. # define TAB putchar('\t')
  14. static int
  15. fullout = 0, /* output full set of geodesic values */
  16. tag = '#', /* beginning of line tag character */
  17. pos_azi = 0, /* output azimuths as positive values */
  18. inverse = 0; /* != 0 then inverse geodesic */
  19. static char
  20. *oform = (char *)0, /* output format for decimal degrees */
  21. *osform = "%.3f", /* output format for S */
  22. pline[50], /* work string */
  23. *usage =
  24. "%s\nusage: %s [ -afFIptTwW [args] ] [ +opts[=arg] ] [ files ]\n";
  25. static void
  26. printLL(double p, double l) {
  27. if (oform) {
  28. (void)printf(oform, p * RAD_TO_DEG); TAB;
  29. (void)printf(oform, l * RAD_TO_DEG);
  30. } else {
  31. (void)fputs(rtodms(pline, p, 'N', 'S'),stdout); TAB;
  32. (void)fputs(rtodms(pline, l, 'E', 'W'),stdout);
  33. }
  34. }
  35. static void
  36. do_arc(void) {
  37. double az;
  38. printLL(phi2, lam2); putchar('\n');
  39. for (az = al12; n_alpha--; ) {
  40. al12 = az = adjlon(az + del_alpha);
  41. geod_pre();
  42. geod_for();
  43. printLL(phi2, lam2); putchar('\n');
  44. }
  45. }
  46. static void /* generate intermediate geodesic coordinates */
  47. do_geod(void) {
  48. double phil, laml, del_S;
  49. phil = phi2;
  50. laml = lam2;
  51. printLL(phi1, lam1); putchar('\n');
  52. for ( geod_S = del_S = geod_S / n_S; --n_S; geod_S += del_S) {
  53. geod_for();
  54. printLL(phi2, lam2); putchar('\n');
  55. }
  56. printLL(phil, laml); putchar('\n');
  57. }
  58. void static /* file processing function */
  59. process(FILE *fid) {
  60. char line[MAXLINE+3], *s;
  61. for (;;) {
  62. ++emess_dat.File_line;
  63. if (!(s = fgets(line, MAXLINE, fid)))
  64. break;
  65. if (!strchr(s, '\n')) { /* overlong line */
  66. int c;
  67. strcat(s, "\n");
  68. /* gobble up to newline */
  69. while ((c = fgetc(fid)) != EOF && c != '\n') ;
  70. }
  71. if (*s == tag) {
  72. fputs(line, stdout);
  73. continue;
  74. }
  75. phi1 = dmstor(s, &s);
  76. lam1 = dmstor(s, &s);
  77. if (inverse) {
  78. phi2 = dmstor(s, &s);
  79. lam2 = dmstor(s, &s);
  80. geod_inv();
  81. } else {
  82. al12 = dmstor(s, &s);
  83. geod_S = strtod(s, &s) * to_meter;
  84. geod_pre();
  85. geod_for();
  86. }
  87. if (!*s && (s > line)) --s; /* assumed we gobbled \n */
  88. if (pos_azi) {
  89. if (al12 < 0.) al12 += TWOPI;
  90. if (al21 < 0.) al21 += TWOPI;
  91. }
  92. if (fullout) {
  93. printLL(phi1, lam1); TAB;
  94. printLL(phi2, lam2); TAB;
  95. if (oform) {
  96. (void)printf(oform, al12 * RAD_TO_DEG); TAB;
  97. (void)printf(oform, al21 * RAD_TO_DEG); TAB;
  98. (void)printf(osform, geod_S * fr_meter);
  99. } else {
  100. (void)fputs(rtodms(pline, al12, 0, 0), stdout); TAB;
  101. (void)fputs(rtodms(pline, al21, 0, 0), stdout); TAB;
  102. (void)printf(osform, geod_S * fr_meter);
  103. }
  104. } else if (inverse)
  105. if (oform) {
  106. (void)printf(oform, al12 * RAD_TO_DEG); TAB;
  107. (void)printf(oform, al21 * RAD_TO_DEG); TAB;
  108. (void)printf(osform, geod_S * fr_meter);
  109. } else {
  110. (void)fputs(rtodms(pline, al12, 0, 0), stdout); TAB;
  111. (void)fputs(rtodms(pline, al21, 0, 0), stdout); TAB;
  112. (void)printf(osform, geod_S * fr_meter);
  113. }
  114. else {
  115. printLL(phi2, lam2); TAB;
  116. if (oform)
  117. (void)printf(oform, al21 * RAD_TO_DEG);
  118. else
  119. (void)fputs(rtodms(pline, al21, 0, 0), stdout);
  120. }
  121. (void)fputs(s, stdout);
  122. }
  123. }
  124. static char *pargv[MAX_PARGS];
  125. static int pargc = 0;
  126. int main(int argc, char **argv) {
  127. char *arg, **eargv = argv, *strnchr();
  128. FILE *fid;
  129. static int eargc = 0, c;
  130. if (emess_dat.Prog_name = strrchr(*argv,'/')) ++emess_dat.Prog_name;
  131. else emess_dat.Prog_name = *argv;
  132. inverse = ! strncmp(emess_dat.Prog_name, "inv", 3);
  133. if (argc <= 1 ) {
  134. (void)fprintf(stderr, usage, pj_get_release(),
  135. emess_dat.Prog_name);
  136. exit (0);
  137. }
  138. /* process run line arguments */
  139. while (--argc > 0) { /* collect run line arguments */
  140. if(**++argv == '-') for(arg = *argv;;) {
  141. switch(*++arg) {
  142. case '\0': /* position of "stdin" */
  143. if (arg[-1] == '-') eargv[eargc++] = "-";
  144. break;
  145. case 'a': /* output full set of values */
  146. fullout = 1;
  147. continue;
  148. case 'I': /* alt. inverse spec. */
  149. inverse = 1;
  150. continue;
  151. case 't': /* set col. one char */
  152. if (arg[1]) tag = *++arg;
  153. else emess(1,"missing -t col. 1 tag");
  154. continue;
  155. case 'W': /* specify seconds precision */
  156. case 'w': /* -W for constant field width */
  157. if ((c = arg[1]) && isdigit(c)) {
  158. set_rtodms(c - '0', *arg == 'W');
  159. ++arg;
  160. } else
  161. emess(1,"-W argument missing or non-digit");
  162. continue;
  163. case 'f': /* alternate output format degrees or xy */
  164. if (--argc <= 0)
  165. noargument: emess(1,"missing argument for -%c",*arg);
  166. oform = *++argv;
  167. continue;
  168. case 'F': /* alternate output format degrees or xy */
  169. if (--argc <= 0) goto noargument;
  170. osform = *++argv;
  171. continue;
  172. case 'l':
  173. if (!arg[1] || arg[1] == 'e') { /* list of ellipsoids */
  174. struct PJ_ELLPS *le;
  175. for (le=pj_get_ellps_ref(); le->id ; ++le)
  176. (void)printf("%9s %-16s %-16s %s\n",
  177. le->id, le->major, le->ell, le->name);
  178. } else if (arg[1] == 'u') { /* list of units */
  179. struct PJ_UNITS *lu;
  180. for (lu = pj_get_units_ref();lu->id ; ++lu)
  181. (void)printf("%12s %-20s %s\n",
  182. lu->id, lu->to_meter, lu->name);
  183. } else
  184. emess(1,"invalid list option: l%c",arg[1]);
  185. exit( 0 );
  186. case 'p': /* output azimuths as positive */
  187. pos_azi = 1;
  188. continue;
  189. default:
  190. emess(1, "invalid option: -%c",*arg);
  191. break;
  192. }
  193. break;
  194. } else if (**argv == '+') /* + argument */
  195. if (pargc < MAX_PARGS)
  196. pargv[pargc++] = *argv + 1;
  197. else
  198. emess(1,"overflowed + argument table");
  199. else /* assumed to be input file name(s) */
  200. eargv[eargc++] = *argv;
  201. }
  202. /* done with parameter and control input */
  203. geod_set(pargc, pargv); /* setup projection */
  204. if ((n_alpha || n_S) && eargc)
  205. emess(1,"files specified for arc/geodesic mode");
  206. if (n_alpha)
  207. do_arc();
  208. else if (n_S)
  209. do_geod();
  210. else { /* process input file list */
  211. if (eargc == 0) /* if no specific files force sysin */
  212. eargv[eargc++] = "-";
  213. for ( ; eargc-- ; ++eargv) {
  214. if (**eargv == '-') {
  215. fid = stdin;
  216. emess_dat.File_name = "<stdin>";
  217. } else {
  218. if ((fid = fopen(*eargv, "r")) == NULL) {
  219. emess(-2, *eargv, "input file");
  220. continue;
  221. }
  222. emess_dat.File_name = *eargv;
  223. }
  224. emess_dat.File_line = 0;
  225. process(fid);
  226. (void)fclose(fid);
  227. emess_dat.File_name = (char *)0;
  228. }
  229. }
  230. exit(0); /* normal completion */
  231. }