00001 /* 00002 * cartography.c ~rick/src/util 00003 * 00004 * Functions for mapping between plate, heliographic, and various map 00005 * coordinate systems, including scaling. 00006 * 00007 * Contents: 00008 * img2sphere Map from plate location to heliographic 00009 * coordinates 00010 * plane2sphere Map from map location to heliographic or 00011 * geographical coordinates 00012 * sphere2img Map from heliographic coordinates to plate 00013 * location 00014 * sphere2plane Map from heliographic/geographic coordinates 00015 * to map location 00016 * 00017 * Responsible: Rick Bogart RBogart@solar.Stanford.EDU 00018 * 00019 * Usage: 00020 * int img2sphere (double x, double y, double ang_r, double latc, 00021 * double lonc, double pa, double *rho, double *lat, double *lon, 00022 * double *sinlat, double *coslat, double *sig, double *mu, double *chi); 00023 * int plane2sphere (double x, double y, double latc, double lonc, 00024 * double *lat, double *lon, int projection); 00025 * int sphere2img (double lat, double lon, double latc, double lonc, 00026 * double *x, double *y, double xcenter, double ycenter, 00027 * double rsun, double peff, double ecc, double chi, 00028 * int xinvrt, int yinvrt); 00029 * int sphere2plane (double lat, double lon, double latc, double lonc, 00030 * double *x, double *y, int projection); 00031 * 00032 * Bugs: 00033 * It is assumed that the function atan2() returns the correct value 00034 * for the quadrant; not all libraries may support this feature. 00035 * sphere2img uses a constant for the correction due to the finite 00036 * distance to the sun. 00037 * sphere2plane and plane2sphere are not true inverses for the following 00038 * projections: (Lambert) cylindrical equal area, sinusoidal equal area 00039 * (Sanson- 00040 * Flamsteed), and Mercator; for these projections, sphere2plane is 00041 * implemented as the normal projection, while plane2sphere is implemented 00042 * as the oblique projection tangent at the normal to the central meridian 00043 * plane2sphere doesn't return 1 if the x coordinate would map to a point 00044 * off the surface. 00045 * 00046 * Planned updates: 00047 * Provide appropriate oblique versions for cylindrical and 00048 * pseudo-cylindrical projections in sphere2plane 00049 * 00050 * Revision history is at end of file. 00051 */ 00052 #include <math.h> 00053 00054 #define RECTANGULAR (0) 00055 #define CASSINI (1) 00056 #define MERCATOR (2) 00057 #define CYLEQA (3) 00058 #define SINEQA (4) 00059 #define GNOMONIC (5) 00060 #define POSTEL (6) 00061 #define STEREOGRAPHIC (7) 00062 #define ORTHOGRAPHIC (8) 00063 #define LAMBERT (9) 00064 00065 static double arc_distance (double lat, double lon, double latc, double lonc) { 00066 double cosa = sin (lat) * sin (latc) + 00067 cos (lat) * cos (latc) * cos (lon - lonc); 00068 return acos (cosa); 00069 } 00070 00071 int img2sphere (double x, double y, double ang_r, double latc, double lonc, 00072 double pa, double *rho, double *lat, double *lon, double *sinlat, 00073 double *coslat, double *sig, double *mu, double *chi) { 00074 /* 00075 * Map projected coordinates (x, y) to (lon, lat) and (rho | sig, chi) 00076 * 00077 * Arguments: 00078 * x } Plate locations, in units of the image radius, relative 00079 * y } to the image center 00080 * ang_r Apparent semi-diameter of sun (angular radius of sun at 00081 * the observer's tangent line) 00082 * latc Latitude of disc center, uncorrected for light travel time 00083 * lonc Longitude of disc center 00084 * pa Position angle of solar north on image, measured eastward 00085 * from north (sky coordinates) 00086 * Return values: 00087 * rho Angle point:sun_center:observer 00088 * lon Heliographic longitude 00089 * lat Heliographic latitude 00090 * sinlat sine of heliographic latitude 00091 * coslat cosine of heliographic latitude 00092 * sig Angle point:observer:sun_center 00093 * mu cosine of angle between the point:observer line and the 00094 * local normal 00095 * chi Position angle on image measured westward from solar 00096 * north 00097 * 00098 * All angles are in radians. 00099 * Return value is 1 if point is outside solar radius (in which case the 00100 * heliographic coordinates and mu are meaningless), 0 otherwise. 00101 * It is assumed that the image is direct; the x or y coordinates require a 00102 * sign change if the image is inverted. 00103 * 00104 */ 00105 static double ang_r0 = 0.0, sinang_r = 0.0, tanang_r = 0.0; 00106 static double latc0 = 0.0, coslatc = 1.0, sinlatc = 0.0; 00107 double cosr, sinr, sinlon, sinsig; 00108 00109 if (ang_r != ang_r0) { 00110 sinang_r = sin (ang_r); 00111 tanang_r = tan (ang_r); 00112 ang_r0 = ang_r; 00113 } 00114 if (latc != latc0) { 00115 sinlatc = sin (latc); 00116 coslatc = cos (latc); 00117 latc0 = latc; 00118 } 00119 *chi = atan2 (x, y) + pa; 00120 while (*chi > 2 * M_PI) *chi -= 2 * M_PI; 00121 while (*chi < 0.0) *chi += 2 * M_PI; 00122 /* Camera curvature correction, no small angle approximations */ 00123 *sig = atan (hypot (x, y) * tanang_r); 00124 sinsig = sin (*sig); 00125 *rho = asin (sinsig / sinang_r) - *sig; 00126 if (*sig > ang_r) return (-1); 00127 *mu = cos (*rho + *sig); 00128 sinr = sin (*rho); 00129 cosr = cos (*rho); 00130 00131 *sinlat = sinlatc * cos (*rho) + coslatc * sinr * cos (*chi); 00132 *coslat = sqrt (1.0 - *sinlat * *sinlat); 00133 *lat = asin (*sinlat); 00134 sinlon = sinr * sin (*chi) / *coslat; 00135 *lon = asin (sinlon); 00136 if (cosr < (*sinlat * sinlatc)) *lon = M_PI - *lon; 00137 *lon += lonc; 00138 while (*lon < 0.0) *lon += 2 * M_PI; 00139 while (*lon >= 2 * M_PI) *lon -= 2 * M_PI; 00140 return (0); 00141 } 00142 00143 int plane2sphere (double x, double y, double latc, double lonc, 00144 double *lat, double *lon, int projection) { 00145 /* 00146 * Perform the inverse mapping from rectangular coordinates x, y on a map 00147 * in a particular projection to heliographic (or geographic) coordinates 00148 * latitude and longitude (in radians). 00149 * The map coordinates are first transformed into arc and azimuth coordinates 00150 * relative to the center of the map according to the appropriate inverse 00151 * transformation for the projection, and thence to latitude and longitude 00152 * from the known heliographic coordinates of the map center (in radians). 00153 * The scale of the map coordinates is assumed to be in units of radians at 00154 * the map center (or other appropriate location of minimum distortion). 00155 * 00156 * Arguments: 00157 * x } Map coordinates, in units of radians at the scale 00158 * y } appropriate to the map center 00159 * latc Latitude of the map center (in radians) 00160 * lonc Longitude of the map center (in radians) 00161 * *lat Returned latitude (in radians) 00162 * *lon Returned longitude (in radians) 00163 * projection A code specifying the map projection to be used: see below 00164 * 00165 * The following projections are supported: 00166 * RECTANGULAR A "rectangular" mapping of x and y directly to 00167 * longitude and latitude, respectively; it is the 00168 * normal cylindrical equidistant projection (plate 00169 * carrée) tangent at the equator and equidistant 00170 * along meridians. Central latitudes off the equator 00171 * merely result in displacement of the map in y 00172 * Also known as CYLEQD 00173 * CASSINI The transverse cylindrical equidistant projection 00174 * (Cassini-Soldner) equidistant along great circles 00175 * perpendicular to the central meridian 00176 * MERCATOR Mercator's conformal projection, in which paths of 00177 * constant bearing are straight lines 00178 * CYLEQA Lambert's normal equal cylindrical (equal-area) 00179 * projection, in which evenly-spaced meridians are 00180 * evenly spaced in x and evenly-spaced parallels are 00181 * separated by the cosine of the latitude 00182 * SINEQA The Sanson-Flamsteed sinusoidal equal-area projection, 00183 * in which evenly-spaced parallels are evenly spaced in 00184 * y and meridians are sinusoidal curves 00185 * GNOMONIC The gnomonic, or central, projection, in which all 00186 * straight lines correspond to geodesics; projection 00187 * from the center of the sphere onto a tangent plane 00188 * POSTEL Postel's azimuthal equidistant projection, in which 00189 * straight lines through the center of the map are 00190 * geodesics with a uniform scale 00191 * STEREOGRAPHIC The stereographic projection, mapping from the 00192 * antipode of the map center onto a tangent plane 00193 * ORTHOGRAPHIC The orthographic projection, mapping from infinity 00194 * onto a tangent plane 00195 * LAMBERT Lambert's azimuthal equivalent projection 00196 * 00197 * The function returns -1 if the requested point on the map does not project 00198 * back to the sphere or is not a principal value, 1 if it projects to a 00199 * point on a hidden hemisphere (if that makes sense), 0 otherwise 00200 */ 00201 static double latc0 = 0.0, sinlatc = 0.0, coslatc = 1.0 ; 00202 double r, rm, test; 00203 double cosr, sinr, cosp, sinp, coslat, sinlat, sinlon; 00204 double sinphi, cosphi, phicom; 00205 int status = 0; 00206 00207 if (latc != latc0) { 00208 coslatc = cos (latc); 00209 sinlatc = sin (latc); 00210 } 00211 latc0 = latc; 00212 00213 switch (projection) { 00214 case (RECTANGULAR): 00215 00216 // printf(" RECTANGULAR projection=%d\n", projection); 00217 00218 *lon = lonc + x; 00219 *lat = latc + y; 00220 if (arc_distance (*lat, *lon, latc, lonc) > M_PI_2) status = 1; 00221 if (fabs (x) > M_PI || fabs (y) > M_PI_2) status = -1; 00222 return status; 00223 case (CASSINI): { 00224 00225 // printf(" CASSINI projection=%d\n", projection); 00226 00227 double sinx = sin (x); 00228 double cosy = cos (y + latc); 00229 double siny = sin (y + latc); 00230 *lat = acos (sqrt (cosy * cosy + siny * siny * sinx * sinx)); 00231 if (y < -latc) *lat *= -1; 00232 *lon = (fabs (*lat) < M_PI_2) ? lonc + asin (sinx / cos (*lat)) : lonc; 00233 if (y > (M_PI_2 - latc) || y < (-M_PI_2 - latc)) 00234 *lon = 2 * lonc + M_PI - *lon; 00235 if (*lon < -M_PI) *lon += 2* M_PI; 00236 if (*lon > M_PI) *lon -= 2 * M_PI; 00237 if (arc_distance (*lat, *lon, latc, lonc) > M_PI_2) status = 1; 00238 if (fabs (x) > M_PI || fabs (y) > M_PI_2) status = -1; 00239 return status; 00240 } 00241 case (CYLEQA): 00242 00243 // printf(" CYLEQA projection=%d\n", projection); 00244 00245 if (fabs (y) > 1.0) { 00246 y = copysign (1.0, y); 00247 status = -1; 00248 } 00249 cosphi = sqrt (1.0 - y*y); 00250 *lat = asin ((y * coslatc) + (cosphi * cos (x) * sinlatc)); 00251 test = (cos (*lat) == 0.0) ? 0.0 : cosphi * sin (x) / cos (*lat); 00252 *lon = asin (test) + lonc; 00253 if (fabs (x) > M_PI_2) { 00254 status = 1; 00255 while (x > M_PI_2) { 00256 *lon = M_PI - *lon; 00257 x -= M_PI; 00258 } 00259 while (x < -M_PI_2) { 00260 *lon = -M_PI - *lon; 00261 x += M_PI; 00262 } 00263 } 00264 if (arc_distance (*lat, *lon, latc, lonc) > M_PI_2) status = 1; 00265 return status; 00266 case (SINEQA): 00267 00268 // printf(" SINEQA projection=%d\n", projection); 00269 00270 cosphi = cos (y); 00271 if (cosphi <= 0.0) { 00272 *lat = y; 00273 *lon = lonc; 00274 if (cosphi < 0.0) status = -1; 00275 return status; 00276 } 00277 *lat = asin ((sin (y) * coslatc) + (cosphi * cos (x/cosphi) * sinlatc)); 00278 coslat = cos (*lat); 00279 if (coslat <= 0.0) { 00280 *lon = lonc; 00281 if (coslat < 0.0) status = 1; 00282 return status; 00283 } 00284 test = cosphi * sin (x/cosphi) / coslat; 00285 *lon = asin (test) + lonc; 00286 if (fabs (x) > M_PI * cosphi) return (-1); 00287 if (fabs (x) > M_PI_2) { 00288 status = 1; 00289 while (x > M_PI_2) { 00290 *lon = M_PI - *lon; 00291 x -= M_PI; 00292 } 00293 while (x < -M_PI_2) { 00294 *lon = -M_PI - *lon; 00295 x += M_PI; 00296 } 00297 } 00298 /* 00299 if (arc_distance (*lat, *lon, latc, lonc) > M_PI_2) status = 1; 00300 */ 00301 return status; 00302 case (MERCATOR): 00303 00304 // printf(" MERCATOR projection=%d\n", projection); 00305 00306 /* 00307 phicom = 2.0 * atan (exp (y)); 00308 sinphi = -cos (phicom); 00309 cosphi = sin (phicom); 00310 *lat = asin ((sinphi * coslatc) + (cosphi * cos (x) * sinlatc)); 00311 *lon = asin (cosphi * sin (x) / cos (*lat)) + lonc; 00312 */ 00313 // modified for CGEM, normal mercator 00314 00315 *lon = x + lonc; 00316 double yc_merc = log (tan (M_PI_4 + 0.5 * latc)); 00317 *lat = atan(exp(y + yc_merc)) * 2 - M_PI_2; 00318 if (arc_distance (*lat, *lon, latc, lonc) > M_PI_2) status = 1; 00319 if (fabs (x) > M_PI_2) status = -1; 00320 00321 return status; 00322 } 00323 /* Convert to polar coordinates */ 00324 r = hypot (x, y); 00325 cosp = (r == 0.0) ? 1.0 : x / r; 00326 sinp = (r == 0.0) ? 0.0 : y / r; 00327 /* Convert to arc */ 00328 switch (projection) { 00329 case (POSTEL): 00330 00331 // printf(" POSTEL projection=%d\n", projection); 00332 00333 rm = r; 00334 if (rm > M_PI_2) status = 1; 00335 break; 00336 case (GNOMONIC): 00337 00338 // printf(" GNOMONIC projection=%d\n", projection); 00339 00340 rm = atan (r); 00341 break; 00342 case (STEREOGRAPHIC): 00343 00344 // printf(" STEREOGRAPHIC projection=%d\n", projection); 00345 00346 rm = 2 * atan (0.5 * r); 00347 if (rm > M_PI_2) status = 1; 00348 break; 00349 case (ORTHOGRAPHIC): 00350 00351 // printf(" ORTHOGRAPHIC projection=%d\n", projection); 00352 00353 if ( r > 1.0 ) { 00354 r = 1.0; 00355 status = -1; 00356 } 00357 rm = asin (r); 00358 break; 00359 case (LAMBERT): 00360 if ( r > 2.0 ) { 00361 r = 2.0; 00362 status = -1; 00363 } 00364 rm = 2 * asin (0.5 * r); 00365 if (rm > M_PI_2 && status == 0) status = 1; 00366 00367 // printf(" LAMBERT projection=%d\n", projection); 00368 00369 break; 00370 } 00371 cosr = cos (rm); 00372 sinr = sin (rm); 00373 /* Convert to latitude-longitude */ 00374 sinlat = sinlatc * cosr + coslatc * sinr * sinp; 00375 *lat = asin (sinlat); 00376 coslat = cos (*lat); 00377 sinlon = (coslat == 0.0) ? 0.0 : sinr * cosp / coslat; 00378 /* This should never happen except for roundoff errors, but just in case: */ 00379 /* 00380 if (sinlon + 1.0 <= 0.0) *lon = -M_PI_2; 00381 else if (sinlon - 1.0 >= 0.0) *lon = M_PI_2; 00382 else 00383 */ 00384 *lon = asin (sinlon); 00385 /* Correction suggested by Dick Shine */ 00386 if (cosr < (sinlat * sinlatc)) *lon = M_PI - *lon; 00387 *lon += lonc; 00388 return status; 00389 } 00390 00391 int sphere2img (double lat, double lon, double latc, double lonc, 00392 double *x, double *y, double xcenter, double ycenter, 00393 double rsun, double peff, double ecc, double chi, 00394 int xinvrt, int yinvrt) { 00395 /* 00396 * Perform a mapping from heliographic coordinates latitude and longitude 00397 * (in radians) to plate location on an image of the sun. The plate 00398 * location is in units of the image radius and is given relative to 00399 * the image center. The function returns 1 if the requested point is 00400 * on the far side (>90 deg from disc center), 0 otherwise. 00401 * 00402 * Arguments: 00403 * lat Latitude (in radians) 00404 * lon Longitude (in radians) 00405 * latc Heliographic latitude of the disc center (in radians) 00406 * lonc Heliographic longitude of the disc center (in radians) 00407 * *x } Plate locations, in units of the image radius, relative 00408 * *y } to the image center 00409 * xcenter } Plate locations of the image center, in units of the 00410 * ycenter } image radius, and measured from an arbitrary origin 00411 * (presumably the plate center or a corner) 00412 * rsun Apparent semi-diameter of the solar disc, in plate 00413 * coordinates 00414 * peff Position angle of the heliographic pole, measured 00415 * eastward from north, relative to the north direction 00416 * on the plate, in radians 00417 * ecc Eccentricity of the fit ellipse presumed due to image 00418 * distortion (no distortion in direction of major axis) 00419 * chi Position angle of the major axis of the fit ellipse, 00420 * measure eastward from north, relative to the north 00421 * direction on the plate, in radians (ignored if ecc = 0) 00422 * xinvrt} Flag parameters: if not equal to 0, the respective 00423 * yinvrt} coordinates on the image x and y are inverted 00424 * 00425 * The heliographic coordinates are first mapped into the polar coordinates 00426 * in an orthographic projection centered at the appropriate location and 00427 * oriented with north in direction of increasing y and west in direction 00428 * of increasing x. The radial coordinate is corrected for foreshortening 00429 * due to the finite distance to the Sun. If the eccentricity of the fit 00430 * ellipse is non-zero the coordinate of the mapped point is proportionately 00431 * reduced in the direction parallel to the minor axis. 00432 * 00433 * Bugs: 00434 * The finite distance correction uses a fixed apparent semi-diameter 00435 * of 16'01'' appropriate to 1.0 AU. In principle the plate radius could 00436 * be used, but this would require the plate scale to be supplied and the 00437 * correction would probably be erroneous and in any case negligible. 00438 * 00439 * The ellipsoidal correction has not been tested very thoroughly. 00440 * 00441 * The return value is based on a test which does not take foreshortening 00442 * into account. 00443 */ 00444 static double sin_asd = 0.004660, cos_asd = 0.99998914; 00445 /* appropriate to 1 AU */ 00446 static double last_latc = 0.0, cos_latc = 1.0, sin_latc = 0.0; 00447 double r, cos_cang, xr, yr; 00448 double sin_lat, cos_lat, cos_lat_lon, cospa, sinpa; 00449 double squash, cchi, schi, c2chi, s2chi, xp, yp; 00450 int hemisphere; 00451 00452 if (latc != last_latc) { 00453 sin_latc = sin (latc); 00454 cos_latc = cos (latc); 00455 last_latc = latc; 00456 } 00457 sin_lat = sin (lat); 00458 cos_lat = cos (lat); 00459 cos_lat_lon = cos_lat * cos (lon - lonc); 00460 00461 cos_cang = sin_lat * sin_latc + cos_latc * cos_lat_lon; 00462 hemisphere = (cos_cang < 0.0) ? 1 : 0; 00463 r = rsun * cos_asd / (1.0 - cos_cang * sin_asd); 00464 xr = r * cos_lat * sin (lon - lonc); 00465 yr = r * (sin_lat * cos_latc - sin_latc * cos_lat_lon); 00466 /* Change sign for inverted images */ 00467 if (xinvrt) xr *= -1.0; 00468 if (yinvrt) yr *= -1.0; 00469 /* Correction for ellipsoidal squashing of image */ 00470 if (ecc > 0.0 && ecc < 1.0) { 00471 squash = sqrt (1.0 - ecc * ecc); 00472 cchi = cos (chi); 00473 schi = sin (chi); 00474 s2chi = schi * schi; 00475 c2chi = 1.0 - s2chi; 00476 xp = xr * (s2chi + squash * c2chi) - yr * (1.0 - squash) * schi * cchi; 00477 yp = yr * (c2chi + squash * s2chi) - xr * (1.0 - squash) * schi * cchi; 00478 xr = xp; 00479 yr = yp; 00480 } 00481 00482 cospa = cos (peff); 00483 sinpa = sin (peff); 00484 *x = xr * cospa - yr * sinpa; 00485 *y = xr * sinpa + yr * cospa; 00486 00487 *y += ycenter; 00488 *x += xcenter; 00489 00490 return (hemisphere); 00491 } 00492 00493 int sphere2plane (double lat, double lon, double latc, double lonc, 00494 double *x, double *y, int projection) { 00495 /* 00496 * Perform a mapping from heliographic (or geographic or celestial) 00497 * coordinates latitude and longitude (in radians) to map location in 00498 * the given projection. The function returns 1 if the requested point is 00499 * on the far side (>90 deg from disc center), 0 otherwise. 00500 * 00501 * Arguments: 00502 * lat Latitude (in radians) 00503 * lon Longitude (in radians) 00504 * latc Heliographic latitude of the disc center (in radians) 00505 * lonc Heliographic longitude of the disc center (in radians) 00506 * *x } Plate locations, in units of the image radius, relative 00507 * *y } to the image center 00508 * projection code specifying the map projection to be used: 00509 * see plane2sphere 00510 */ 00511 static double last_latc = 0.0, cos_latc = 1.0, sin_latc = 0.0, yc_merc; 00512 double r, rm, cos_cang; 00513 double sin_lat, cos_lat, cos_lat_lon; 00514 double cos_phi, sin_phi; 00515 int hemisphere; 00516 00517 if (latc != last_latc) { 00518 sin_latc = sin (latc); 00519 cos_latc = cos (latc); 00520 last_latc = latc; 00521 yc_merc = log (tan (M_PI_4 + 0.5 * latc)); 00522 } 00523 sin_lat = sin (lat); 00524 cos_lat = cos (lat); 00525 cos_lat_lon = cos_lat * cos (lon - lonc); 00526 cos_cang = sin_lat * sin_latc + cos_latc * cos_lat_lon; 00527 hemisphere = (cos_cang < 0.0) ? 1 : 0; 00528 /* Cylindrical projections */ 00529 switch (projection) { 00530 /* Normal cylindrical equidistant */ 00531 case (RECTANGULAR): 00532 *x = lon - lonc; 00533 *y = lat - latc; 00534 return hemisphere; 00535 /* Transverse cylindrical equidistant */ 00536 case (CASSINI): 00537 *x = asin (cos_lat * sin (lon - lonc)); 00538 *y = atan2 (tan (lat), cos (lon - lonc)) - latc; 00539 return hemisphere; 00540 /* Normal cylindrical equivalent - differs from sphere2plane */ 00541 case (CYLEQA): 00542 *x = lon - lonc; 00543 *y = sin_lat - sin_latc; 00544 return hemisphere; 00545 /* Normal sinusoidal equivalent - differs from sphere2plane */ 00546 case (SINEQA): 00547 *x = cos_lat * (lon - lonc); 00548 *y = lat - latc; 00549 return hemisphere; 00550 /* Normal Mercator - differs from sphere2plane */ 00551 case (MERCATOR): 00552 *x = lon - lonc; 00553 *y = log (tan (M_PI_4 + 0.5 * lat)) - yc_merc; 00554 return hemisphere; 00555 } 00556 /* Azimuthal projections */ 00557 rm = acos (cos_cang); 00558 00559 switch (projection) { 00560 case (POSTEL): 00561 r = rm; 00562 break; 00563 case (GNOMONIC): 00564 r = tan (rm); 00565 break; 00566 case (STEREOGRAPHIC): 00567 r = 2.0 * tan (0.5 * rm); 00568 break; 00569 case (ORTHOGRAPHIC): 00570 r = sin (rm); 00571 break; 00572 case (LAMBERT): 00573 r = 2.0 * sin (0.5 * rm); 00574 break; 00575 default: 00576 return -1; 00577 } 00578 if (rm != 0.) { 00579 *x = r * cos_lat * sin (lon - lonc) / sin (rm); 00580 *y = r * (sin_lat * cos_latc - sin_latc * cos_lat_lon) / sin(rm); 00581 } else { 00582 *x = 0.; 00583 *y = 0.; 00584 } 00585 return hemisphere; 00586 } 00587 00588 /* 00589 * Revision History 00590 * v 0.9 94.08.03 Rick Bogart created this file 00591 * 94.10.27 R Bogart reorganized for inclusion in libM 00592 * 94.11.24 R Bogart fixed bug in sinusoidal equal-area 00593 * mapping and minor bug in cylindrical equal-area mapoing 00594 * (latc != 0); added Mercator projection; more comments 00595 * 94.12.12 R Bogart & Luiz Sa fixed non-azimuthal projections 00596 * in plane2sphere to support center of map at other than L0L0; 00597 * removed unnecessary (and incorrect) scaling in azimuthal 00598 * mappings 00599 * v 1.0 95.08.18 R Bogart & L Sa implemented optimum_scale(); 00600 * added arguments and code to sphere2img() to deal with elliptic 00601 * images and image inversion; additional comments; 00602 * minor fix in SINEQA option of plane2sphere(); changed sign on 00603 * on position angle used in sphere2img() to conform with standard 00604 * usage (positive eastward from north in sky) 00605 * 96.05.17 R Bogart fixed bug in longitude calculation of 00606 * points "over the pole" for azimuthal projections in plane2sphere() 00607 * and corrected bug in sphere2img(). 00608 * 96.06.07 R Bogart return status of sphere2img indicates 00609 * whether mapped point is on far side of globe. 00610 * 96.12.18 R Bogart removed debugging prints from function 00611 * optimum_scale() 00612 * 97.01.21 R Bogart plane2sphere now returns status 00613 * 98.09.02 R Bogart added img2sphere 00614 * 99.02.25 R Bogart fixed calculated chi in img2sphere to 00615 * conform to man page description; plane2sphere now returns 1 if 00616 * requested point is not principal value in rectangular, cylindrical 00617 * equal-area, sinusoidal equal-area, and Mercator projections 00618 * 99.12.13 R Bogart added sphere2plane 00619 * 00.04.24 R Bogart modifications to plane2sphere to avoid 00620 * occasional roundoff errors leading to undefined results at the 00621 * limb 00622 * 05.01.07 R Bogart fixed bug in longitude calculation in 00623 * plane2sphere() for points more than 90deg from target longitude 00624 * of map center 00625 * 06.11.10 R Bogart fixed return value of sphere2plane for 00626 * cylindrical projection(s); added support for Cassini's projection 00627 * (transverse cylindrical equidistant); added support for normal 00628 * cylindrical and pseudocylindrical projections to sphere2plane 00629 * 09.07.02 R Bogart copied to private location for inclusion in 00630 * DRMS modules and standalone code; removed optimum_scale(), pow2scale(), 00631 * name2proj(), proj2name() 00632 */
1.5.7.1