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