summaryrefslogtreecommitdiff
path: root/src/panorama.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/panorama.c')
-rw-r--r--src/panorama.c369
1 files changed, 369 insertions, 0 deletions
diff --git a/src/panorama.c b/src/panorama.c
new file mode 100644
index 0000000..16dfc9b
--- /dev/null
+++ b/src/panorama.c
@@ -0,0 +1,369 @@
+/*
+ * "$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 <stdio.h>
+#include <unistd.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+
+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 <phi>,<lambda>\n"
+ " -b <phi>,<lambda>\n");
+}