1 xudong 1.1 /*
2 * cartography.c ~rick/src/util
3 *
4 * Functions for mapping between plate, heliographic, and various map
5 * coordinate systems, including scaling.
6 *
7 * Contents:
8 * img2sphere Map from plate location to heliographic
9 * coordinates
10 * plane2sphere Map from map location to heliographic or
11 * geographical coordinates
12 * sphere2img Map from heliographic coordinates to plate
13 * location
14 * sphere2plane Map from heliographic/geographic coordinates
15 * to map location
16 *
17 * Responsible: Rick Bogart RBogart@solar.Stanford.EDU
18 *
19 * Usage:
20 * int img2sphere (double x, double y, double ang_r, double latc,
21 * double lonc, double pa, double *rho, double *lat, double *lon,
22 xudong 1.1 * double *sinlat, double *coslat, double *sig, double *mu, double *chi);
23 * int plane2sphere (double x, double y, double latc, double lonc,
24 * double *lat, double *lon, int projection);
25 * int sphere2img (double lat, double lon, double latc, double lonc,
26 * double *x, double *y, double xcenter, double ycenter,
27 * double rsun, double peff, double ecc, double chi,
28 * int xinvrt, int yinvrt);
29 * int sphere2plane (double lat, double lon, double latc, double lonc,
30 * double *x, double *y, int projection);
31 *
32 * Bugs:
33 * It is assumed that the function atan2() returns the correct value
34 * for the quadrant; not all libraries may support this feature.
35 * sphere2img uses a constant for the correction due to the finite
36 * distance to the sun.
37 * sphere2plane and plane2sphere are not true inverses for the following
38 * projections: (Lambert) cylindrical equal area, sinusoidal equal area
39 * (Sanson-
40 * Flamsteed), and Mercator; for these projections, sphere2plane is
41 * implemented as the normal projection, while plane2sphere is implemented
42 * as the oblique projection tangent at the normal to the central meridian
43 xudong 1.1 * plane2sphere doesn't return 1 if the x coordinate would map to a point
44 * off the surface.
45 *
46 * Planned updates:
47 * Provide appropriate oblique versions for cylindrical and
48 * pseudo-cylindrical projections in sphere2plane
49 *
50 * Revision history is at end of file.
51 */
52 #include <math.h>
53
54 #define RECTANGULAR (0)
55 #define CASSINI (1)
56 #define MERCATOR (2)
57 #define CYLEQA (3)
58 #define SINEQA (4)
59 #define GNOMONIC (5)
60 #define POSTEL (6)
61 #define STEREOGRAPHIC (7)
62 #define ORTHOGRAPHIC (8)
63 #define LAMBERT (9)
64 xudong 1.1
65 static double arc_distance (double lat, double lon, double latc, double lonc) {
66 double cosa = sin (lat) * sin (latc) +
67 cos (lat) * cos (latc) * cos (lon - lonc);
68 return acos (cosa);
69 }
70
71 int img2sphere (double x, double y, double ang_r, double latc, double lonc,
72 double pa, double *rho, double *lat, double *lon, double *sinlat,
73 double *coslat, double *sig, double *mu, double *chi) {
74 /*
75 * Map projected coordinates (x, y) to (lon, lat) and (rho | sig, chi)
76 *
77 * Arguments:
78 * x } Plate locations, in units of the image radius, relative
79 * y } to the image center
80 * ang_r Apparent semi-diameter of sun (angular radius of sun at
81 * the observer's tangent line)
82 * latc Latitude of disc center, uncorrected for light travel time
83 * lonc Longitude of disc center
84 * pa Position angle of solar north on image, measured eastward
85 xudong 1.1 * from north (sky coordinates)
86 * Return values:
87 * rho Angle point:sun_center:observer
88 * lon Heliographic longitude
89 * lat Heliographic latitude
90 * sinlat sine of heliographic latitude
91 * coslat cosine of heliographic latitude
92 * sig Angle point:observer:sun_center
93 * mu cosine of angle between the point:observer line and the
94 * local normal
95 * chi Position angle on image measured westward from solar
96 * north
97 *
98 * All angles are in radians.
99 * Return value is 1 if point is outside solar radius (in which case the
100 * heliographic coordinates and mu are meaningless), 0 otherwise.
101 * It is assumed that the image is direct; the x or y coordinates require a
102 * sign change if the image is inverted.
103 *
104 */
105 static double ang_r0 = 0.0, sinang_r = 0.0, tanang_r = 0.0;
106 xudong 1.1 static double latc0 = 0.0, coslatc = 1.0, sinlatc = 0.0;
107 double cosr, sinr, sinlon, sinsig;
108
109 if (ang_r != ang_r0) {
110 sinang_r = sin (ang_r);
111 tanang_r = tan (ang_r);
112 ang_r0 = ang_r;
113 }
114 if (latc != latc0) {
115 sinlatc = sin (latc);
116 coslatc = cos (latc);
117 latc0 = latc;
118 }
119 *chi = atan2 (x, y) + pa;
120 while (*chi > 2 * M_PI) *chi -= 2 * M_PI;
121 while (*chi < 0.0) *chi += 2 * M_PI;
122 /* Camera curvature correction, no small angle approximations */
123 *sig = atan (hypot (x, y) * tanang_r);
124 sinsig = sin (*sig);
125 *rho = asin (sinsig / sinang_r) - *sig;
126 if (*sig > ang_r) return (-1);
127 xudong 1.1 *mu = cos (*rho + *sig);
128 sinr = sin (*rho);
129 cosr = cos (*rho);
130
131 *sinlat = sinlatc * cos (*rho) + coslatc * sinr * cos (*chi);
132 *coslat = sqrt (1.0 - *sinlat * *sinlat);
133 *lat = asin (*sinlat);
134 sinlon = sinr * sin (*chi) / *coslat;
135 *lon = asin (sinlon);
136 if (cosr < (*sinlat * sinlatc)) *lon = M_PI - *lon;
137 *lon += lonc;
138 while (*lon < 0.0) *lon += 2 * M_PI;
139 while (*lon >= 2 * M_PI) *lon -= 2 * M_PI;
140 return (0);
141 }
142
143 int plane2sphere (double x, double y, double latc, double lonc,
144 double *lat, double *lon, int projection) {
145 /*
146 * Perform the inverse mapping from rectangular coordinates x, y on a map
147 * in a particular projection to heliographic (or geographic) coordinates
148 xudong 1.1 * latitude and longitude (in radians).
149 * The map coordinates are first transformed into arc and azimuth coordinates
150 * relative to the center of the map according to the appropriate inverse
151 * transformation for the projection, and thence to latitude and longitude
152 * from the known heliographic coordinates of the map center (in radians).
153 * The scale of the map coordinates is assumed to be in units of radians at
154 * the map center (or other appropriate location of minimum distortion).
155 *
156 * Arguments:
157 * x } Map coordinates, in units of radians at the scale
158 * y } appropriate to the map center
159 * latc Latitude of the map center (in radians)
160 * lonc Longitude of the map center (in radians)
161 * *lat Returned latitude (in radians)
162 * *lon Returned longitude (in radians)
163 * projection A code specifying the map projection to be used: see below
164 *
165 * The following projections are supported:
166 * RECTANGULAR A "rectangular" mapping of x and y directly to
167 * longitude and latitude, respectively; it is the
168 * normal cylindrical equidistant projection (plate
169 xudong 1.1 * carrée) tangent at the equator and equidistant
170 * along meridians. Central latitudes off the equator
171 * merely result in displacement of the map in y
172 * Also known as CYLEQD
173 * CASSINI The transverse cylindrical equidistant projection
174 * (Cassini-Soldner) equidistant along great circles
175 * perpendicular to the central meridian
176 * MERCATOR Mercator's conformal projection, in which paths of
177 * constant bearing are straight lines
178 * CYLEQA Lambert's normal equal cylindrical (equal-area)
179 * projection, in which evenly-spaced meridians are
180 * evenly spaced in x and evenly-spaced parallels are
181 * separated by the cosine of the latitude
182 * SINEQA The Sanson-Flamsteed sinusoidal equal-area projection,
183 * in which evenly-spaced parallels are evenly spaced in
184 * y and meridians are sinusoidal curves
185 * GNOMONIC The gnomonic, or central, projection, in which all
186 * straight lines correspond to geodesics; projection
187 * from the center of the sphere onto a tangent plane
188 * POSTEL Postel's azimuthal equidistant projection, in which
189 * straight lines through the center of the map are
190 xudong 1.1 * geodesics with a uniform scale
191 * STEREOGRAPHIC The stereographic projection, mapping from the
192 * antipode of the map center onto a tangent plane
193 * ORTHOGRAPHIC The orthographic projection, mapping from infinity
194 * onto a tangent plane
195 * LAMBERT Lambert's azimuthal equivalent projection
196 *
197 * The function returns -1 if the requested point on the map does not project
198 * back to the sphere or is not a principal value, 1 if it projects to a
199 * point on a hidden hemisphere (if that makes sense), 0 otherwise
200 */
201 static double latc0 = 0.0, sinlatc = 0.0, coslatc = 1.0 ;
202 double r, rm, test;
203 double cosr, sinr, cosp, sinp, coslat, sinlat, sinlon;
204 double sinphi, cosphi, phicom;
205 int status = 0;
206
207 if (latc != latc0) {
208 coslatc = cos (latc);
209 sinlatc = sin (latc);
210 }
211 xudong 1.1 latc0 = latc;
212
213 switch (projection) {
214 case (RECTANGULAR):
215
216 // printf(" RECTANGULAR projection=%d\n", projection);
217
218 *lon = lonc + x;
219 *lat = latc + y;
220 if (arc_distance (*lat, *lon, latc, lonc) > M_PI_2) status = 1;
221 if (fabs (x) > M_PI || fabs (y) > M_PI_2) status = -1;
222 return status;
223 case (CASSINI): {
224
225 // printf(" CASSINI projection=%d\n", projection);
226
227 double sinx = sin (x);
228 double cosy = cos (y + latc);
229 double siny = sin (y + latc);
230 *lat = acos (sqrt (cosy * cosy + siny * siny * sinx * sinx));
231 if (y < -latc) *lat *= -1;
232 xudong 1.1 *lon = (fabs (*lat) < M_PI_2) ? lonc + asin (sinx / cos (*lat)) : lonc;
233 if (y > (M_PI_2 - latc) || y < (-M_PI_2 - latc))
234 *lon = 2 * lonc + M_PI - *lon;
235 if (*lon < -M_PI) *lon += 2* M_PI;
236 if (*lon > M_PI) *lon -= 2 * M_PI;
237 if (arc_distance (*lat, *lon, latc, lonc) > M_PI_2) status = 1;
238 if (fabs (x) > M_PI || fabs (y) > M_PI_2) status = -1;
239 return status;
240 }
241 case (CYLEQA):
242
243 // printf(" CYLEQA projection=%d\n", projection);
244
245 if (fabs (y) > 1.0) {
246 y = copysign (1.0, y);
247 status = -1;
248 }
249 cosphi = sqrt (1.0 - y*y);
250 *lat = asin ((y * coslatc) + (cosphi * cos (x) * sinlatc));
251 test = (cos (*lat) == 0.0) ? 0.0 : cosphi * sin (x) / cos (*lat);
252 *lon = asin (test) + lonc;
253 xudong 1.1 if (fabs (x) > M_PI_2) {
254 status = 1;
255 while (x > M_PI_2) {
256 *lon = M_PI - *lon;
257 x -= M_PI;
258 }
259 while (x < -M_PI_2) {
260 *lon = -M_PI - *lon;
261 x += M_PI;
262 }
263 }
264 if (arc_distance (*lat, *lon, latc, lonc) > M_PI_2) status = 1;
265 return status;
266 case (SINEQA):
267
268 // printf(" SINEQA projection=%d\n", projection);
269
270 cosphi = cos (y);
271 if (cosphi <= 0.0) {
272 *lat = y;
273 *lon = lonc;
274 xudong 1.1 if (cosphi < 0.0) status = -1;
275 return status;
276 }
277 *lat = asin ((sin (y) * coslatc) + (cosphi * cos (x/cosphi) * sinlatc));
278 coslat = cos (*lat);
279 if (coslat <= 0.0) {
280 *lon = lonc;
281 if (coslat < 0.0) status = 1;
282 return status;
283 }
284 test = cosphi * sin (x/cosphi) / coslat;
285 *lon = asin (test) + lonc;
286 if (fabs (x) > M_PI * cosphi) return (-1);
287 if (fabs (x) > M_PI_2) {
288 status = 1;
289 while (x > M_PI_2) {
290 *lon = M_PI - *lon;
291 x -= M_PI;
292 }
293 while (x < -M_PI_2) {
294 *lon = -M_PI - *lon;
295 xudong 1.1 x += M_PI;
296 }
297 }
298 /*
299 if (arc_distance (*lat, *lon, latc, lonc) > M_PI_2) status = 1;
300 */
301 return status;
302 case (MERCATOR):
303
304 // printf(" MERCATOR projection=%d\n", projection);
305
306 phicom = 2.0 * atan (exp (y));
307 sinphi = -cos (phicom);
308 cosphi = sin (phicom);
309 *lat = asin ((sinphi * coslatc) + (cosphi * cos (x) * sinlatc));
310 *lon = asin (cosphi * sin (x) / cos (*lat)) + lonc;
311 if (arc_distance (*lat, *lon, latc, lonc) > M_PI_2) status = 1;
312 if (fabs (x) > M_PI_2) status = -1;
313
314 return status;
315 }
316 xudong 1.1 /* Convert to polar coordinates */
317 r = hypot (x, y);
318 cosp = (r == 0.0) ? 1.0 : x / r;
319 sinp = (r == 0.0) ? 0.0 : y / r;
320 /* Convert to arc */
321 switch (projection) {
322 case (POSTEL):
323
324 // printf(" POSTEL projection=%d\n", projection);
325
326 rm = r;
327 if (rm > M_PI_2) status = 1;
328 break;
329 case (GNOMONIC):
330
331 // printf(" GNOMONIC projection=%d\n", projection);
332
333 rm = atan (r);
334 break;
335 case (STEREOGRAPHIC):
336
337 xudong 1.1 // printf(" STEREOGRAPHIC projection=%d\n", projection);
338
339 rm = 2 * atan (0.5 * r);
340 if (rm > M_PI_2) status = 1;
341 break;
342 case (ORTHOGRAPHIC):
343
344 // printf(" ORTHOGRAPHIC projection=%d\n", projection);
345
346 if ( r > 1.0 ) {
347 r = 1.0;
348 status = -1;
349 }
350 rm = asin (r);
351 break;
352 case (LAMBERT):
353 if ( r > 2.0 ) {
354 r = 2.0;
355 status = -1;
356 }
357 rm = 2 * asin (0.5 * r);
358 xudong 1.1 if (rm > M_PI_2 && status == 0) status = 1;
359
360 // printf(" LAMBERT projection=%d\n", projection);
361
362 break;
363 }
364 cosr = cos (rm);
365 sinr = sin (rm);
366 /* Convert to latitude-longitude */
367 sinlat = sinlatc * cosr + coslatc * sinr * sinp;
368 *lat = asin (sinlat);
369 coslat = cos (*lat);
370 sinlon = (coslat == 0.0) ? 0.0 : sinr * cosp / coslat;
371 /* This should never happen except for roundoff errors, but just in case: */
372 /*
373 if (sinlon + 1.0 <= 0.0) *lon = -M_PI_2;
374 else if (sinlon - 1.0 >= 0.0) *lon = M_PI_2;
375 else
376 */
377 *lon = asin (sinlon);
378 /* Correction suggested by Dick Shine */
379 xudong 1.1 if (cosr < (sinlat * sinlatc)) *lon = M_PI - *lon;
380 *lon += lonc;
381 return status;
382 }
383
384 int sphere2img (double lat, double lon, double latc, double lonc,
385 double *x, double *y, double xcenter, double ycenter,
386 double rsun, double peff, double ecc, double chi,
387 int xinvrt, int yinvrt) {
388 /*
389 * Perform a mapping from heliographic coordinates latitude and longitude
390 * (in radians) to plate location on an image of the sun. The plate
391 * location is in units of the image radius and is given relative to
392 * the image center. The function returns 1 if the requested point is
393 * on the far side (>90 deg from disc center), 0 otherwise.
394 *
395 * Arguments:
396 * lat Latitude (in radians)
397 * lon Longitude (in radians)
398 * latc Heliographic latitude of the disc center (in radians)
399 * lonc Heliographic longitude of the disc center (in radians)
400 xudong 1.1 * *x } Plate locations, in units of the image radius, relative
401 * *y } to the image center
402 * xcenter } Plate locations of the image center, in units of the
403 * ycenter } image radius, and measured from an arbitrary origin
404 * (presumably the plate center or a corner)
405 * rsun Apparent semi-diameter of the solar disc, in plate
406 * coordinates
407 * peff Position angle of the heliographic pole, measured
408 * eastward from north, relative to the north direction
409 * on the plate, in radians
410 * ecc Eccentricity of the fit ellipse presumed due to image
411 * distortion (no distortion in direction of major axis)
412 * chi Position angle of the major axis of the fit ellipse,
413 * measure eastward from north, relative to the north
414 * direction on the plate, in radians (ignored if ecc = 0)
415 * xinvrt} Flag parameters: if not equal to 0, the respective
416 * yinvrt} coordinates on the image x and y are inverted
417 *
418 * The heliographic coordinates are first mapped into the polar coordinates
419 * in an orthographic projection centered at the appropriate location and
420 * oriented with north in direction of increasing y and west in direction
421 xudong 1.1 * of increasing x. The radial coordinate is corrected for foreshortening
422 * due to the finite distance to the Sun. If the eccentricity of the fit
423 * ellipse is non-zero the coordinate of the mapped point is proportionately
424 * reduced in the direction parallel to the minor axis.
425 *
426 * Bugs:
427 * The finite distance correction uses a fixed apparent semi-diameter
428 * of 16'01'' appropriate to 1.0 AU. In principle the plate radius could
429 * be used, but this would require the plate scale to be supplied and the
430 * correction would probably be erroneous and in any case negligible.
431 *
432 * The ellipsoidal correction has not been tested very thoroughly.
433 *
434 * The return value is based on a test which does not take foreshortening
435 * into account.
436 */
437 static double sin_asd = 0.004660, cos_asd = 0.99998914;
438 /* appropriate to 1 AU */
439 static double last_latc = 0.0, cos_latc = 1.0, sin_latc = 0.0;
440 double r, cos_cang, xr, yr;
441 double sin_lat, cos_lat, cos_lat_lon, cospa, sinpa;
442 xudong 1.1 double squash, cchi, schi, c2chi, s2chi, xp, yp;
443 int hemisphere;
444
445 if (latc != last_latc) {
446 sin_latc = sin (latc);
447 cos_latc = cos (latc);
448 last_latc = latc;
449 }
450 sin_lat = sin (lat);
451 cos_lat = cos (lat);
452 cos_lat_lon = cos_lat * cos (lon - lonc);
453
454 cos_cang = sin_lat * sin_latc + cos_latc * cos_lat_lon;
455 hemisphere = (cos_cang < 0.0) ? 1 : 0;
456 r = rsun * cos_asd / (1.0 - cos_cang * sin_asd);
457 xr = r * cos_lat * sin (lon - lonc);
458 yr = r * (sin_lat * cos_latc - sin_latc * cos_lat_lon);
459 /* Change sign for inverted images */
460 if (xinvrt) xr *= -1.0;
461 if (yinvrt) yr *= -1.0;
462 /* Correction for ellipsoidal squashing of image */
463 xudong 1.1 if (ecc > 0.0 && ecc < 1.0) {
464 squash = sqrt (1.0 - ecc * ecc);
465 cchi = cos (chi);
466 schi = sin (chi);
467 s2chi = schi * schi;
468 c2chi = 1.0 - s2chi;
469 xp = xr * (s2chi + squash * c2chi) - yr * (1.0 - squash) * schi * cchi;
470 yp = yr * (c2chi + squash * s2chi) - xr * (1.0 - squash) * schi * cchi;
471 xr = xp;
472 yr = yp;
473 }
474
475 cospa = cos (peff);
476 sinpa = sin (peff);
477 *x = xr * cospa - yr * sinpa;
478 *y = xr * sinpa + yr * cospa;
479
480 *y += ycenter;
481 *x += xcenter;
482
483 return (hemisphere);
484 xudong 1.1 }
485
486 int sphere2plane (double lat, double lon, double latc, double lonc,
487 double *x, double *y, int projection) {
488 /*
489 * Perform a mapping from heliographic (or geographic or celestial)
490 * coordinates latitude and longitude (in radians) to map location in
491 * the given projection. The function returns 1 if the requested point is
492 * on the far side (>90 deg from disc center), 0 otherwise.
493 *
494 * Arguments:
495 * lat Latitude (in radians)
496 * lon Longitude (in radians)
497 * latc Heliographic latitude of the disc center (in radians)
498 * lonc Heliographic longitude of the disc center (in radians)
499 * *x } Plate locations, in units of the image radius, relative
500 * *y } to the image center
501 * projection code specifying the map projection to be used:
502 * see plane2sphere
503 */
504 static double last_latc = 0.0, cos_latc = 1.0, sin_latc = 0.0, yc_merc;
505 xudong 1.1 double r, rm, cos_cang;
506 double sin_lat, cos_lat, cos_lat_lon;
507 double cos_phi, sin_phi;
508 int hemisphere;
509
510 if (latc != last_latc) {
511 sin_latc = sin (latc);
512 cos_latc = cos (latc);
513 last_latc = latc;
514 yc_merc = log (tan (M_PI_4 + 0.5 * latc));
515 }
516 sin_lat = sin (lat);
517 cos_lat = cos (lat);
518 cos_lat_lon = cos_lat * cos (lon - lonc);
519 cos_cang = sin_lat * sin_latc + cos_latc * cos_lat_lon;
520 hemisphere = (cos_cang < 0.0) ? 1 : 0;
521 /* Cylindrical projections */
522 switch (projection) {
523 /* Normal cylindrical equidistant */
524 case (RECTANGULAR):
525 *x = lon - lonc;
526 xudong 1.1 *y = lat - latc;
527 return hemisphere;
528 /* Transverse cylindrical equidistant */
529 case (CASSINI):
530 *x = asin (cos_lat * sin (lon - lonc));
531 *y = atan2 (tan (lat), cos (lon - lonc)) - latc;
532 return hemisphere;
533 /* Normal cylindrical equivalent - differs from sphere2plane */
534 case (CYLEQA):
535 *x = lon - lonc;
536 *y = sin_lat - sin_latc;
537 return hemisphere;
538 /* Normal sinusoidal equivalent - differs from sphere2plane */
539 case (SINEQA):
540 *x = cos_lat * (lon - lonc);
541 *y = lat - latc;
542 return hemisphere;
543 /* Normal Mercator - differs from sphere2plane */
544 case (MERCATOR):
545 *x = lon - lonc;
546 *y = log (tan (M_PI_4 + 0.5 * lat)) - yc_merc;
547 xudong 1.1 return hemisphere;
548 }
549 /* Azimuthal projections */
550 rm = acos (cos_cang);
551
552 switch (projection) {
553 case (POSTEL):
554 r = rm;
555 break;
556 case (GNOMONIC):
557 r = tan (rm);
558 break;
559 case (STEREOGRAPHIC):
560 r = 2.0 * tan (0.5 * rm);
561 break;
562 case (ORTHOGRAPHIC):
563 r = sin (rm);
564 break;
565 case (LAMBERT):
566 r = 2.0 * sin (0.5 * rm);
567 break;
568 xudong 1.1 default:
569 return -1;
570 }
571 if (rm != 0.) {
572 *x = r * cos_lat * sin (lon - lonc) / sin (rm);
573 *y = r * (sin_lat * cos_latc - sin_latc * cos_lat_lon) / sin(rm);
574 } else {
575 *x = 0.;
576 *y = 0.;
577 }
578 return hemisphere;
579 }
580
581 /*
582 * Revision History
583 * v 0.9 94.08.03 Rick Bogart created this file
584 * 94.10.27 R Bogart reorganized for inclusion in libM
585 * 94.11.24 R Bogart fixed bug in sinusoidal equal-area
586 * mapping and minor bug in cylindrical equal-area mapoing
587 * (latc != 0); added Mercator projection; more comments
588 * 94.12.12 R Bogart & Luiz Sa fixed non-azimuthal projections
589 xudong 1.1 * in plane2sphere to support center of map at other than L0L0;
590 * removed unnecessary (and incorrect) scaling in azimuthal
591 * mappings
592 * v 1.0 95.08.18 R Bogart & L Sa implemented optimum_scale();
593 * added arguments and code to sphere2img() to deal with elliptic
594 * images and image inversion; additional comments;
595 * minor fix in SINEQA option of plane2sphere(); changed sign on
596 * on position angle used in sphere2img() to conform with standard
597 * usage (positive eastward from north in sky)
598 * 96.05.17 R Bogart fixed bug in longitude calculation of
599 * points "over the pole" for azimuthal projections in plane2sphere()
600 * and corrected bug in sphere2img().
601 * 96.06.07 R Bogart return status of sphere2img indicates
602 * whether mapped point is on far side of globe.
603 * 96.12.18 R Bogart removed debugging prints from function
604 * optimum_scale()
605 * 97.01.21 R Bogart plane2sphere now returns status
606 * 98.09.02 R Bogart added img2sphere
607 * 99.02.25 R Bogart fixed calculated chi in img2sphere to
608 * conform to man page description; plane2sphere now returns 1 if
609 * requested point is not principal value in rectangular, cylindrical
610 xudong 1.1 * equal-area, sinusoidal equal-area, and Mercator projections
611 * 99.12.13 R Bogart added sphere2plane
612 * 00.04.24 R Bogart modifications to plane2sphere to avoid
613 * occasional roundoff errors leading to undefined results at the
614 * limb
615 * 05.01.07 R Bogart fixed bug in longitude calculation in
616 * plane2sphere() for points more than 90deg from target longitude
617 * of map center
618 * 06.11.10 R Bogart fixed return value of sphere2plane for
619 * cylindrical projection(s); added support for Cassini's projection
620 * (transverse cylindrical equidistant); added support for normal
621 * cylindrical and pseudocylindrical projections to sphere2plane
622 * 09.07.02 R Bogart copied to private location for inclusion in
623 * DRMS modules and standalone code; removed optimum_scale(), pow2scale(),
624 * name2proj(), proj2name()
625 */
|
Return to
cartography.c
CVS log
Up to [Development] / JSOC / proj / sharp / apps