![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to lox.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [Research Unix] / researchv10no / cmd / map / lox|
Return to lox.c CVS log | Up to [Research Unix] / researchv10no / cmd / map / lox |
1.1 root 1: #include <stdio.h>
2: #include <math.h>
3: /* development of a loxodrome from a unit sphere onto a plane
4: the loxodrome passes through 0,0 on the sphere,
5: inclined a (alpha) from the equator
6:
7: latitude p (phi) positive north
8: longitude t (theta) positive east - opposite to "map"
9: all angles in radians
10:
11: x,y are coords of developed curve
12: c (chi) is angle of inclination of developed curve from x axis
13:
14: dx/dp = cos(c)/sin(a)
15: dy/dp = sin(c)/sin(a)
16: dc/dp = tan(p)/tan(a)
17: dt/dp = 1/(cos(p)*tan(a))
18:
19: cc lox.c -lport -lF77 -lI77
20: a.out alpha here alpha is in degrees
21: */
22: extern char *realloc();
23: static double renorm();
24: static bsearch(), getxy(), calcf(), solve();
25: #define DEG (180/PI)
26:
27: /* curlicue is an ordered sequence of points
28: along the midline of the spiral, in both geographic
29: and developed coordinates */
30: struct curl {
31: double phi,x,y,theta,chi;
32: } *curlicue; /* the midline */
33: int nc; /* number of points on the midline */
34:
35: double a = PI/10;/* 45 degrees */
36: double sa,ca,ta; /* sin(a), cos(a), tan(a) */
37: static signa = 1;
38: static double z[4]; /* x==z[0],y==z[1],t=z[2],c==z[3] */
39:
40: /* inputs to dodes, port library integration routine */
41:
42: static int n = 4;
43: static double z0[4]; /* initial values = 0,0,0,a */
44: static double p0; /* initial value of p = 0 */
45: static double p1 = 86*PI/180; /* final value */
46: static double dt; /* current time step */
47: static double deg = PI/180; /* 1 degree */
48: static float errpar[2] = {1.e-7,1.e-7};
49:
50: static
51: f(p,z,n,r)
52: double *p, *z, *r;
53: int *n;
54: {
55: r[0] = cos(z[3])/sa;
56: r[1] = sin(z[3])/sa;
57: r[2] = 1/(cos(*p)*ta);
58: r[3] = sin(*p)*r[2];
59: }
60: static
61: handle(t0,z0,t1,z1,n,dt,ts,e)
62: double *t0,*z0,*t1,*z1,*dt,*ts;
63: float *e;
64: int *n;
65: {
66: if(*dt > deg) *dt = deg;
67: if(fabs(z1[3]-z0[3])>.1 && *dt >.00001)
68: *dt *= .1/fabs(z1[3]-z0[3]);
69: if(*t1 + *dt >= PI/2)
70: *dt = (PI/2 - *t1)/2;
71: if(*t1 + *dt > *ts)
72: *dt = *ts - *t1 + .00001;
73: curlicue = (struct curl*)realloc((char*)curlicue,
74: (nc+1)*sizeof(*curlicue));
75: curlicue[nc].phi = *t1;
76: curlicue[nc].x = z1[0];
77: curlicue[nc].y = z1[1];
78: curlicue[nc].theta = z1[2];
79: curlicue[nc].chi = z1[3];
80: nc++;
81: printf("%f %f %f %f %f\n",*t1*DEG,
82: z1[0]*DEG,z1[1]*DEG,z1[2]*DEG,z1[3]*DEG);
83: fflush(stdout);
84: if(*t0==*t1) abort();
85: }
86: main(argc,argv)
87: char **argv;
88: {
89: double u,v,x,y;
90: a = atof(argv[1])/DEG;
91: if(a < 0)
92: signa = -1;
93: sa = sin(a);
94: ca = cos(a);
95: ta = sa/ca;
96: z0[3] = a;
97: dt = deg;
98: handle(°,z0,&p0,z0,&n,°,&p1,errpar);
99: dodes_(f,z0,&n,&p0,&p1,&dt,errpar,handle);
100: while(scanf("%lf %lf",&u,&v)==2) {
101: double u1 = u/DEG, v1=v/DEG;
102: double d,p,t;
103: v1 = renorm(u1,bsearch(u1,0),v1);
104: solve(u1,v1,&d,&p,&t);
105: printf("%f %f %f %f %f %f\n",
106: u,v,v1*DEG,d*DEG,p*DEG,t*DEG);
107: fflush(stdout);
108: getxy(d,p,&x,&y);
109: printf("%f %f\n",x*DEG,y*DEG);
110: fflush(stdout);
111: }
112: }
113: /* find the index for latitude (i=0) or longitude (i=1)
114: in array curlicue. In case of longitude, u must have
115: been normalized by renorm() */
116: static
117: bsearch(u,i)
118: double u;
119: {
120: int top = nc;
121: int bot = 0;
122: int mid;
123: while(bot<top-1) {
124: mid = (bot+top)/2;
125: if(i==0&&curlicue[mid].phi<=u||
126: i==1&&(signa>0&&curlicue[mid].theta<=u||
127: signa<0&&curlicue[mid].theta>=u))
128: bot = mid;
129: else
130: top = mid;
131: }
132: return bot;
133: }
134: /* wind v thru multiples of 2*PI until it lies within PI
135: of the (linearly interpolated) midline */
136: static double
137: renorm(u,n,v)
138: double u,v;
139: {
140: double t,fract;
141: fract = (u-curlicue[n].phi)/
142: (curlicue[n+1].phi-curlicue[n].phi);
143: t = curlicue[n].theta*(1-fract) +
144: curlicue[n+1].theta*fract;
145: if(signa > 0) while(v<t-PI)
146: v += 2*PI;
147: else while(v>t+PI)
148: v -= 2*PI;
149: return v;
150: }
151: /* given the height (d) and foot (p,t) of the perpendicular projection
152: of a point onto the midline, return the x-y coordinates of the point*/
153: static
154: getxy(d,p,xp,yp)
155: double d,p,*xp,*yp;
156: {
157: int n, sign = 1;
158: double fract, x0, y0, c0;
159: if(p < 0) {
160: p = -p;
161: sign = -1;
162: }
163: n = bsearch(p,0);
164: fract = (p-curlicue[n].phi)/
165: (curlicue[n+1].phi-curlicue[n].phi);
166: /*printf("n,fract %d %f\n",n,fract);*/
167: x0 = sign*(curlicue[n].x*(1-fract) +
168: curlicue[n+1].x*fract);
169: y0 = sign*(curlicue[n].y*(1-fract) +
170: curlicue[n+1].y*fract);
171: c0 = curlicue[n].chi*(1-fract) +
172: curlicue[n+1].chi*fract;
173: *xp = x0 - d*sin(c0);
174: *yp = y0 + d*cos(c0);
175: }
176:
177: /*
178: solve simultaneously for d,p,t (dist to midline,
179: phi lat on midline,
180: theta lon on midline)
181: given a, u, v (rhumb angle from parallels, lat, east lon)
182:
183: sin d sin a = cos u sin(t-v) law of sines
184: sin u = cos d sin p + sin d cos p cos a law of cosines
185: t tan a = log tan (p/2 + PI/4) eqn of rhumb line
186:
187: initial guess: d=u-p, t=v, p from table lookup in curlicue
188: */
189: static double su,cu,v0;
190: static int n3 = 3, jmax = 50;
191: static double acc = 1e-6;
192: static
193: calcf(n,x,nf,g)
194: int *n, *nf;
195: double *x,*g;
196: {
197: double f[3], sd;
198: if(x[1]>=PI/2) {
199: *nf = 0;
200: return;
201: }
202: sd = sin(x[0]);
203: f[0] = sd*sa-cu*sin(x[2]-v0);
204: f[1] = su - cos(x[0])*sin(x[1]) - sd*cos(x[1])*ca;
205: f[2] = x[2]*ta - log(tan(x[1]/2+PI/4));
206: /*printf("calcf: d p t %f %f %f",
207: x[0]*DEG,x[1]*DEG,x[2]*DEG);
208: printf(" resid %f %f %f\n",f[0],f[1],f[2]);
209: fflush(stdout);*/
210: *g = f[0]*f[0]+f[1]*f[1]+f[2]*f[2];
211: }
212: /* given renormed lat-lon (u,v) find the height (*dp) and foot
213: (*pp,*tp) of a perpendicular to the midline. precondition: u>=0 */
214: static
215: solve(u,v,dp,pp,tp)
216: double u,v,*dp,*pp,*tp;
217: {
218: double z[3]; /* d, p, t */
219: int n;
220: if(u < 0) {
221: solve(-u,-v,dp,pp,tp);
222: *dp = -*dp;
223: *pp = -*pp;
224: *tp = -*tp;
225: return;
226: }
227: n = bsearch(signa*fabs(v),1);
228: su = sin(u);
229: cu = cos(u);
230: z[1] = curlicue[n].phi; /* initial guess */
231: if(signa*v < 0) z[1] = -z[1];
232: z[0] = u - z[1];
233: z[2] = v0 = v;
234: printf("inputs (deg) %f %f\n",u*DEG,v*DEG);
235: printf("interval %d guess (deg) %f %f %f\n",n,
236: z[0]*DEG,z[1]*DEG,z[2]*DEG);
237: fflush(stdout);
238: dsmnf_(&n3,z,calcf,&jmax,&acc);
239: *dp = z[0];
240: *pp = z[1];
241: *tp = z[2];
242: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.