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