/* * "$Id: panorama.c,v 1.7 2005/04/11 14:40:18 hofmann Exp $" * * flpsed program. * * Copyright 2004,2005 by Johannes Hofmann * * This program is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public * License as published by the Free Software Foundation; either * version 2 of the License, or (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public * License along with this library; if not, write to the Free Software * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 * USA. */ #include #include #include #include #include double pi, deg2rad; double distance(double phi_a, double lam_a, double phi_b, double lam_b); double sin_alpha(double lam_a, double lam_b, double phi_b, double c); double cos_alpha(double phi_a, double phi_b, double c); double alpha(double phi_a, double lam_a, double phi_b, double lam_b); double center_angle(double alph_a, double alph_b, double d1, double d2); int read_file(double phi_a, double lam_a, char *name, int dist, double h_d_ratio, double mm_per_deg, double a_center); int main(int argc, char **argv) { char c; char *sep, *tmp, *data_file = NULL; int errflg = 0; int dist = 100000; double phi_a, lam_a, phi_b, lam_b; double h_d_ratio = 0.1, mm_per_deg = 10.0, dist_c_m1 = 0.0; double dist_c_m2 = 0.0, a_center = 0.0; char *pos = NULL, *m1 = NULL, *m2 = NULL, *m_c = NULL; pi = asin(1.0) * 2.0; deg2rad = pi / 180.0; while ((c = getopt(argc, argv, "a:b:f:d:r:p:m:n:s:t:c:")) != EOF) { switch (c) { case 'a': tmp = strdup(optarg); if ((sep = strchr(tmp, ',')) == NULL) { errflg++; break; } *sep = '\0'; phi_a = atof(tmp) * deg2rad; lam_a = atof(sep+1) * deg2rad; free(tmp); break; case 'b': tmp = strdup(optarg); if ((sep = strchr(tmp, ',')) == NULL) { usage(); } *sep = '\0'; phi_b = atof(tmp) * deg2rad; lam_b = atof(sep+1) * deg2rad; free(tmp); break; case 'f': data_file = optarg; break; case 'd': dist = atoi(optarg); break; case 'r': h_d_ratio = atof(optarg); break; case 's': dist_c_m1 = atof(optarg); break; case 't': dist_c_m2 = atof(optarg); break; case 'p': pos = optarg; break; case 'm': m1 = optarg; break; case 'n': m2 = optarg; break; case 'c': m_c = optarg; break; } } if (errflg) { usage(); exit(1); } if (!data_file) { usage(); exit(1); } if (pos){ int n = get_pos(data_file, pos, &phi_a, &lam_a); if (n == 0) { fprintf(stderr, "No matching entry for %s in datafile.\n", pos); exit(1); } else if (n > 1) { fprintf(stderr, "More than one matching entry for %s in datafile.\n", pos); exit(1); } } if (dist_c_m1 > 0 && dist_c_m2 && m1 && m2) { double phi, lam; double a1, a2, a; int n; n = get_pos(data_file, m1, &phi, &lam); if (n == 0) { fprintf(stderr, "No matching entry for %s in datafile.\n", m1); exit(1); } else if (n > 1) { fprintf(stderr, "More than one matching entry for %s in datafile.\n", m1); exit(1); } a1 = alpha(phi_a, lam_a, phi, lam); n = get_pos(data_file, m2, &phi, &lam); if (n == 0) { fprintf(stderr, "No matching entry for %s in datafile.\n", m2); exit(1); } else if (n > 1) { fprintf(stderr, "More than one matching entry for %s in datafile.\n", m2); exit(1); } a2 = alpha(phi_a, lam_a, phi, lam); a_center = center_angle(a1, a2, dist_c_m1, dist_c_m2); fprintf(stderr, "center_angle = %f\n", a_center / deg2rad); } read_file(phi_a, lam_a, data_file, dist, h_d_ratio, mm_per_deg, a_center); } #define PS_HEADER "%!\n \ /mountain { % expects name, height, dist, deg on stack\n \ newpath % Start a new path\n \ gsave % Keep rotations temporary\n \ rotate % Rotate by degrees on stack\n \ dup \n \ 0 moveto \n \ add % End line at 20 + dist + height\n \ dup % Duplicate end of line pos\n \ 0 lineto\n \ 10 add % Start text 10 from end of line \n \ 0 moveto\n \ show\n \ stroke\n \ grestore % Get back the unrotated state\n \ } def \n \ /Times-Roman findfont \n \ 8 scalefont \n \ setfont \n \ 300 344 translate\n \ 5 0 moveto\n \ 0 0 5 0 360 arc\n \ stroke %11 0 moveto\n" #define PS_HEADER_PHOTO "%!\n \ /mountain { % expects name, height, dist, deg on stack\n \ newpath\n \ a_center sub\n \ mm_per_deg mul\n \ 0.352778 div % convert from mm to PS points\n \ dup\n \ dup\n \ 200 exch moveto\n \ 500 exch lineto\n \ 500 exch moveto\n \ pop \n \ pop \n \ show\n \ stroke\n \ } def\n \ \n \ /Times-Roman findfont \n \ 6 scalefont \n \ setfont \n 0 420 translate \n" int get_pos(char *filename, char *name, double *phi, double *lam) { FILE *fp; char buf[1024]; char *vals[10]; char **ap, *bp; double height; int found = 0; fp = fopen(filename, "r"); if (!fp) { perror("fopen"); return 1; } while (fgets(buf, sizeof(buf), fp)) { bp = buf; for (ap = vals; (*ap = strsep(&bp, ",")) != NULL;) if (++ap >= &vals[10]) break; if (strstr(vals[1], name)) { *phi = atof(vals[3]) * deg2rad; *lam = atof(vals[4]) * deg2rad; height = atof(vals[5]); fprintf(stderr, "Found matching entry: %s (%fm)\n", vals[1], height); found++; } } fclose(fp); return found; } int read_file(double phi_a, double lam_a, char *name, int dist, double h_d_ratio, double mm_per_deg, double a_center) { FILE *fp; char buf[1024]; char *vals[10]; char **ap, *bp; double phi_b, lam_b, height; fp = fopen(name, "r"); if (!fp) { perror("fopen"); return 1; } printf("%s", PS_HEADER_PHOTO); printf("/mm_per_deg %f def\n", mm_per_deg); printf("/a_center %f def\n", a_center); while (fgets(buf, sizeof(buf), fp)) { bp = buf; for (ap = vals; (*ap = strsep(&bp, ",")) != NULL;) if (++ap >= &vals[10]) break; phi_b = atof(vals[3]) * deg2rad; lam_b = atof(vals[4]) * deg2rad; if (phi_a == phi_b && lam_a == lam_b) { continue; } height = atof(vals[5]); if (distance(phi_a, lam_a, phi_b, lam_b) * 6368000.0 > dist) { continue; } if (height / (distance(phi_a, lam_a, phi_b, lam_b)* 6368000) < h_d_ratio) { continue; } printf("(%s) %f %f %f mountain\n", vals[1], height / 20.0, distance(phi_a, lam_a, phi_b, lam_b) * 25000.0, alpha(phi_a, lam_a, phi_b, lam_b)); #if 0 printf("%f %f %s %fm\n", alpha(phi_a, lam_a, phi_b, lam_b), distance(phi_a, lam_a, phi_b, lam_b) * 6368000, vals[1], height); #endif } printf("showpage\n"); fclose(fp); return 0; } double distance(double phi_a, double lam_a, double phi_b, double lam_b) { return acos(sin(phi_a)*sin(phi_b) + cos(phi_a)*cos(phi_b)*cos(lam_a - lam_b)); } double sin_alpha(double lam_a, double lam_b, double phi_b, double c) { return sin(lam_b - lam_a) * cos(phi_b) / sin(c); } double cos_alpha(double phi_a, double phi_b, double c) { return (sin(phi_b) - sin(phi_a) * cos(c)) / (cos(phi_a) * sin(c)); } double alpha(double phi_a, double lam_a, double phi_b, double lam_b) { double dist, sin_alph, cos_alph, alph; dist = distance(phi_a, lam_a, phi_b, lam_b); sin_alph = sin_alpha(lam_a, lam_b, phi_b, dist); cos_alph = cos_alpha(phi_a, phi_b, dist); if (sin_alph > 0) { alph = acos(cos_alph); } else { alph = 2.0 * pi - acos(cos_alph); } return alph; } double center_angle(double alph_a, double alph_b, double d1, double d2) { double tan_a, tan_b; fprintf(stderr, "a=%f, b=%f d1=%f d2=%f\n", alph_a / deg2rad, alph_b / deg2rad, d1, d2); tan_a = tan(alph_a - alph_b); tan_b = (d2 - d1 - ((sqrt((d2*(d2 - (2*d1*(1 + (2*(tan_a * tan_a)))))) + (d1*d1))))) / (2*d2*tan_a); return alph_a + atan(tan_b); } void usage() { fprintf(stderr, "usage: panorama -a ,\n" " -b ,\n"); }