Changeset 13839
- Timestamp:
- Jun 14, 2007, 3:03:22 PM (19 years ago)
- Location:
- trunk/psLib/src/astro
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/psLib/src/astro/psCoord.c
r12737 r13839 1 1 /** @file psCoord.c 2 *3 * @brief Contains basic coordinate transformation definitions and operations4 *5 * This file defines the basic types for astronomical coordinate6 * transformation7 *8 * @ingroup CoordinateTransform9 *10 * @author GLG, MHPCC11 *12 * @version $Revision: 1.138$ $Name: not supported by cvs2svn $13 * @date $Date: 2007-04-04 21:10:03$14 *15 * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii16 */2 * 3 * @brief Contains basic coordinate transformation definitions and operations 4 * 5 * This file defines the basic types for astronomical coordinate 6 * transformation 7 * 8 * @ingroup CoordinateTransform 9 * 10 * @author GLG, MHPCC 11 * 12 * @version $Revision: 1.139 $ $Name: not supported by cvs2svn $ 13 * @date $Date: 2007-06-15 01:03:22 $ 14 * 15 * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii 16 */ 17 17 /******************************************************************************/ 18 18 /* INCLUDE FILES */ … … 72 72 XXX: below is the code using the standard matrix representation. note that 73 73 this inversion requires x->nX == 1, y->nY == 1 and x->nY <= 1, y->nX <= 1 74 74 *****************************************************************************/ 75 75 psPlaneTransform *p_psPlaneTransformLinearInvert(psPlaneTransform *transform) 76 76 { … … 79 79 PS_ASSERT_PTR_NON_NULL(transform->y, NULL); 80 80 if ((transform->x->nX < 1) || 81 82 83 81 (transform->x->nY < 1) || 82 (transform->y->nX < 1) || 83 (transform->y->nY < 1)) { 84 84 psLogMsg(__func__, PS_LOG_WARN, "WARNING: input transform is not invertible."); 85 85 return(NULL); … … 150 150 151 151 Why isn't this called p_psIsPlaneTransformLinear()? 152 152 *****************************************************************************/ 153 153 bool p_psIsProjectionLinear(psPlaneTransform *transform) 154 154 { … … 161 161 if (transform->x->coeff[i][j] != 0.0) { 162 162 if (!(((i == 0) && (j == 0)) || 163 164 163 ((i == 0) && (j == 1)) || 164 ((i == 1) && (j == 0)))) { 165 165 return(false); 166 166 } … … 173 173 if (transform->y->coeff[i][j] != 0.0) { 174 174 if (!(((i == 0) && (j == 0)) || 175 176 175 ((i == 0) && (j == 1)) || 176 ((i == 1) && (j == 0)))) { 177 177 return(false); 178 178 } … … 304 304 This transformation takes into account parameters beyond an objects spatial 305 305 coordinates: term3 and term4 (magnitude and color). 306 306 *****************************************************************************/ 307 307 psPlane* psPlaneDistortApply( 308 308 psPlane* out, … … 379 379 psProjectionClass class = PS_PROJECTION_CLASS_NONE; 380 380 switch (projection->type) { 381 case PS_PROJ_LIN:382 case PS_PROJ_PLY:383 case PS_PROJ_WRP:381 case PS_PROJ_LIN: 382 case PS_PROJ_PLY: 383 case PS_PROJ_WRP: 384 384 class = PS_PROJECTION_CLASS_CARTESIAN; 385 385 break; 386 case PS_PROJ_TAN: 387 case PS_PROJ_DIS: 388 case PS_PROJ_SIN: 389 case PS_PROJ_STG: 390 case PS_PROJ_ZEA: 391 case PS_PROJ_ZPL: 386 case PS_PROJ_TAN: 387 case PS_PROJ_TNX: 388 case PS_PROJ_DIS: 389 case PS_PROJ_SIN: 390 case PS_PROJ_STG: 391 case PS_PROJ_ZEA: 392 case PS_PROJ_ZPL: 392 393 class = PS_PROJECTION_CLASS_ZENITHAL; 393 394 break; 394 case PS_PROJ_AIT:395 case PS_PROJ_PAR:396 case PS_PROJ_GLS:397 case PS_PROJ_MER:398 case PS_PROJ_CAR:395 case PS_PROJ_AIT: 396 case PS_PROJ_PAR: 397 case PS_PROJ_GLS: 398 case PS_PROJ_MER: 399 case PS_PROJ_CAR: 399 400 class = PS_PROJECTION_CLASS_PSEUDOCYL; 400 401 break; 401 402 402 default:403 default: 403 404 psAbort("invalid projection"); 404 405 break; … … 414 415 415 416 switch (class) { 416 case PS_PROJECTION_CLASS_CARTESIAN:417 case PS_PROJECTION_CLASS_CARTESIAN: 417 418 out->x = (coord->r - projection->R); 418 419 out->y = (coord->d - projection->D); 419 420 break; 420 421 421 case PS_PROJECTION_CLASS_ZENITHAL:422 case PS_PROJECTION_CLASS_ZENITHAL: 422 423 sinDp = sin(projection->D); 423 424 cosDp = cos(projection->D); … … 427 428 cosDelta = cos(coord->d); 428 429 429 430 # if ELIXIR_CODE 430 431 431 432 sinTheta = cosDelta*cosAlpha*cosDp + sinDelta*sinDp; 432 433 434 433 sinPhiCT = cosDelta*sinAlpha; 434 cosPhiCT = cosDelta*cosAlpha*sinDp - sinDelta*cosDp; 435 # else 435 436 436 437 sinTheta = cosDelta*cosAlpha*cosDp + sinDelta*sinDp; 437 438 439 440 441 switch (projection->type) 442 { 443 case PS_PROJ_TAN: 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 case PS_PROJECTION_CLASS_PSEUDOCYL:438 sinPhiCT = -cosDelta*sinAlpha; 439 cosPhiCT = -cosDelta*cosAlpha*sinDp + sinDelta*cosDp; 440 # endif 441 // Perform the specified projection 442 switch (projection->type) { 443 case PS_PROJ_TAN: 444 case PS_PROJ_TNX: // XXX this is not quite right 445 case PS_PROJ_DIS: 446 // Gnomonic projection 447 out->x = +sinPhiCT / sinTheta; 448 out->y = -cosPhiCT / sinTheta; 449 break; 450 case PS_PROJ_SIN: 451 // Othrographic projection 452 out->x = +sinPhiCT; 453 out->y = -cosPhiCT; 454 break; 455 case PS_PROJ_STG: 456 // Othrographic projection 457 out->x = +2*sinPhiCT / (1 + sinTheta); 458 out->y = -2*cosPhiCT / (1 + sinTheta); 459 break; 460 case PS_PROJ_ZEA: 461 case PS_PROJ_ZPL: 462 // Zenithal Equal-Area 463 zeta = M_SQRT2 / sqrt (1 + sinTheta); 464 out->x = +zeta * sinPhiCT; 465 out->y = -zeta * cosPhiCT; 466 break; 467 default: 468 psAbort("invalid projection"); 469 break; 470 } 471 break; 472 473 case PS_PROJECTION_CLASS_PSEUDOCYL: 473 474 phi = coord->r - projection->R; 474 475 theta = coord->d - projection->D; 475 476 switch (projection->type) 476 477 { 477 478 case PS_PROJ_AIT: 478 479 // Hammer-Aitoff projection 479 480 zeta = 1.0/sqrt(0.5*(1.0+cos(theta)*cos(phi/2.0))); … … 481 482 out->y = zeta*sin(theta); 482 483 break; 483 484 case PS_PROJ_GLS: 484 485 // projection name? 485 486 out->x = phi * cos(theta); 486 487 out->x = theta; 487 488 break; 488 489 case PS_PROJ_PAR: 489 490 // Parabolic projection 490 491 out->x = phi*(2.0*cos(2.0*theta/3.0) - 1.0); 491 492 out->y = M_PI*sin(theta/3.0); 492 493 break; 493 494 494 case PS_PROJ_CAR: 495 case PS_PROJ_MER: 495 496 psAbort("projection not yet implemented"); 496 497 break; 497 498 default: 498 499 psAbort("invalid projection"); 499 500 break; … … 501 502 break; 502 503 503 default:504 default: 504 505 psAbort("invalid projection"); 505 506 break; … … 535 536 psProjectionClass class = PS_PROJECTION_CLASS_NONE; 536 537 switch (projection->type) { 537 case PS_PROJ_LIN:538 case PS_PROJ_PLY:539 case PS_PROJ_WRP:538 case PS_PROJ_LIN: 539 case PS_PROJ_PLY: 540 case PS_PROJ_WRP: 540 541 class = PS_PROJECTION_CLASS_CARTESIAN; 541 542 break; 542 case PS_PROJ_TAN: 543 case PS_PROJ_DIS: 544 case PS_PROJ_SIN: 545 case PS_PROJ_STG: 546 case PS_PROJ_ZEA: 547 case PS_PROJ_ZPL: 543 case PS_PROJ_TAN: 544 case PS_PROJ_TNX: 545 case PS_PROJ_DIS: 546 case PS_PROJ_SIN: 547 case PS_PROJ_STG: 548 case PS_PROJ_ZEA: 549 case PS_PROJ_ZPL: 548 550 class = PS_PROJECTION_CLASS_ZENITHAL; 549 551 break; 550 case PS_PROJ_AIT:551 case PS_PROJ_PAR:552 case PS_PROJ_GLS:553 case PS_PROJ_CAR:554 case PS_PROJ_MER:552 case PS_PROJ_AIT: 553 case PS_PROJ_PAR: 554 case PS_PROJ_GLS: 555 case PS_PROJ_CAR: 556 case PS_PROJ_MER: 555 557 class = PS_PROJECTION_CLASS_PSEUDOCYL; 556 558 break; 557 559 558 default:560 default: 559 561 psAbort("invalid projection"); 560 562 break; … … 575 577 // Perform inverse projection 576 578 switch (class) { 577 case PS_PROJECTION_CLASS_CARTESIAN:579 case PS_PROJECTION_CLASS_CARTESIAN: 578 580 out->r = coord->x*projection->Xs + projection->R; 579 581 out->d = coord->y*projection->Ys + projection->D; 580 582 break; 581 583 582 case PS_PROJECTION_CLASS_ZENITHAL:584 case PS_PROJECTION_CLASS_ZENITHAL: 583 585 // Remove plate scales 584 586 x = coord->x*projection->Xs; … … 588 590 cosPhi = (R == 0) ? 1.0 : -y / R; 589 591 590 switch (projection->type) 591 { 592 case PS_PROJ_TAN: 593 592 switch (projection->type) { 593 case PS_PROJ_TAN: 594 case PS_PROJ_TNX:// XXX this is not quite right 595 case PS_PROJ_DIS: 594 596 // Gnonomic deprojection 595 597 rho = sqrt (1 + R*R); … … 597 599 cosTheta = R / rho; 598 600 break; 599 601 case PS_PROJ_SIN: 600 602 // Orhtographic deprojection 601 603 sinTheta = sqrt (1 - R*R); 602 604 cosTheta = R; 603 605 break; 604 606 case PS_PROJ_STG: 605 607 psAbort("STG not defined"); 606 608 break; 607 608 609 case PS_PROJ_ZEA: 610 case PS_PROJ_ZPL: 609 611 if (R > 2) 610 612 return NULL; … … 612 614 cosTheta = sqrt (1 - PS_SQR(sinTheta)); 613 615 break; 614 616 default: 615 617 psAbort("invalid projection"); 616 618 break; … … 622 624 // Convert from projection spherical coordinates 623 625 // psLib versions: 624 626 # if ELIXIR_CODE 625 627 // XXX the elixir version : does the ADD have a sign error? 626 628 psF64 delta = asin(sinTheta*sinDp - cosTheta*cosPhi*cosDp); 627 629 psF64 sinAlpha = +cosTheta*sinPhi; 628 630 psF64 cosAlpha = +cosTheta*cosPhi*sinDp + sinTheta*cosDp; 629 630 631 631 # else 632 633 psF64 delta = asin(sinTheta*sinDp + cosTheta*cosPhi*cosDp); 632 634 psF64 sinAlpha = -cosTheta*sinPhi; 633 635 psF64 cosAlpha = -cosTheta*cosPhi*sinDp + sinTheta*cosDp; 634 636 # endif 635 637 636 638 out->d = delta; … … 638 640 break; 639 641 640 case PS_PROJECTION_CLASS_PSEUDOCYL:642 case PS_PROJECTION_CLASS_PSEUDOCYL: 641 643 switch (projection->type) 642 644 { 643 645 case PS_PROJ_AIT: 644 646 // Hammer-Aitoff deprojection 645 647 // XXX EAM : need range check on z^2 : must be > 0 … … 652 654 theta = asin(y*rho); 653 655 break; 654 656 case PS_PROJ_PAR: 655 657 // Parabolic deprojection 656 658 rho = y/M_PI; … … 658 660 theta = 3.0*asin(rho); 659 661 break; 660 662 case PS_PROJ_GLS: 661 663 phi = x/cos(y); 662 664 theta = y; 663 665 break; 664 666 default: 665 667 psAbort("invalid projection"); 666 668 break; … … 669 671 out->d = theta + projection->D; 670 672 break; 671 default:673 default: 672 674 psAbort("invalid projection"); 673 675 break; … … 682 684 multiplies them. Basically, for each non-zero coeff in the trans1 coeff[][] 683 685 array, you must multiply by all non-zero coeffs in trans2. 684 686 *****************************************************************************/ 685 687 static psPolynomial2D *multiplyDPoly2D( 686 688 psPolynomial2D *trans1, … … 713 715 /***************************************************************************** 714 716 psPlaneTransformCombine(out, trans1, trans2) 715 717 *****************************************************************************/ 716 718 psPlaneTransform *psPlaneTransformCombine( 717 719 psPlaneTransform *out, … … 757 759 } else { 758 760 if ((out->x->nX == orderX) && 759 760 761 761 (out->x->nY == orderY) && 762 (out->y->nX == orderX) && 763 (out->y->nY == orderY)) { 762 764 myPT = out; 763 765 // … … 885 887 XXX: This code ignores nRejIter and sigmaClip. We must call the ClipFit 886 888 routines instead. 887 889 *****************************************************************************/ 888 890 bool psPlaneTransformFit( 889 891 psPlaneTransform *trans, … … 932 934 psPlaneTransformInvert(out, in, region, nSamples) 933 935 934 936 *****************************************************************************/ 935 937 psPlaneTransform *psPlaneTransformInvert( 936 938 psPlaneTransform *out, … … 972 974 } else { 973 975 if ((out->x->nX == order) && (out->x->nY == order) && 974 976 (out->y->nX == order) && (out->y->nX == order)) { 975 977 myPT = out; 976 978 } else { … … 1030 1032 const psPlaneTransform *transformation, 1031 1033 const psPlane *coord 1032 )1034 ) 1033 1035 { 1034 1036 PS_ASSERT_PTR_NON_NULL(transformation, NULL); … … 1121 1123 fxnVal = psPlaneTransformApply(fxnVal, inToOut, coord); 1122 1124 if (fabs(fxnVal->x - coord->x) <= fabs(deriv->x) && 1123 1125 fabs(fxnVal->y - coord->y) <= fabs(deriv->y)) { 1124 1126 int x = (int)(ceil(fabs(deriv->x))); 1125 1127 int y = (int)(ceil(fabs(deriv->y))); … … 1205 1207 1206 1208 switch (type) { 1207 case PS_PROJ_LIN:1209 case PS_PROJ_LIN: 1208 1210 psStringAppend (&name, "%s-LIN", prefix); 1209 1211 return name; 1210 case PS_PROJ_PLY:1212 case PS_PROJ_PLY: 1211 1213 psStringAppend (&name, "%s-PLY", prefix); 1212 1214 return name; 1213 case PS_PROJ_WRP:1215 case PS_PROJ_WRP: 1214 1216 psStringAppend (&name, "%s-WRP", prefix); 1215 1217 return name; 1216 case PS_PROJ_TAN:1218 case PS_PROJ_TAN: 1217 1219 psStringAppend (&name, "%s-TAN", prefix); 1218 1220 return name; 1219 case PS_PROJ_DIS: 1221 case PS_PROJ_TNX: 1222 psStringAppend (&name, "%s-TNX", prefix); 1223 return name; 1224 case PS_PROJ_DIS: 1220 1225 psStringAppend (&name, "%s-DIS", prefix); 1221 1226 return name; 1222 case PS_PROJ_SIN:1227 case PS_PROJ_SIN: 1223 1228 psStringAppend (&name, "%s-SIN", prefix); 1224 1229 return name; 1225 case PS_PROJ_STG:1230 case PS_PROJ_STG: 1226 1231 psStringAppend (&name, "%s-STG", prefix); 1227 1232 return name; 1228 case PS_PROJ_AIT:1233 case PS_PROJ_AIT: 1229 1234 psStringAppend (&name, "%s-AIT", prefix); 1230 1235 return name; 1231 case PS_PROJ_PAR:1236 case PS_PROJ_PAR: 1232 1237 psStringAppend (&name, "%s-PAR", prefix); 1233 1238 return name; 1234 case PS_PROJ_GLS:1239 case PS_PROJ_GLS: 1235 1240 psStringAppend (&name, "%s-GLS", prefix); 1236 1241 return name; 1237 case PS_PROJ_CAR:1242 case PS_PROJ_CAR: 1238 1243 psStringAppend (&name, "%s-CAR", prefix); 1239 1244 return name; 1240 case PS_PROJ_MER:1245 case PS_PROJ_MER: 1241 1246 psStringAppend (&name, "%s-MER", prefix); 1242 1247 return name; 1243 1248 1244 default:1249 default: 1245 1250 psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unknown projection type: %d\n", type); 1246 1251 return NULL; … … 1257 1262 if (!strcmp (&name[4], "-WRP")) return PS_PROJ_WRP; 1258 1263 if (!strcmp (&name[4], "-TAN")) return PS_PROJ_TAN; 1264 if (!strcmp (&name[4], "-TNX")) return PS_PROJ_TNX; 1259 1265 if (!strcmp (&name[4], "-DIS")) return PS_PROJ_DIS; 1260 1266 if (!strcmp (&name[4], "-SIN")) return PS_PROJ_SIN; -
trunk/psLib/src/astro/psCoord.h
r12737 r13839 8 8 * @author GLG, MHPCC 9 9 * 10 * @version $Revision: 1. 59$ $Name: not supported by cvs2svn $11 * @date $Date: 2007-0 4-04 21:10:03$10 * @version $Revision: 1.60 $ $Name: not supported by cvs2svn $ 11 * @date $Date: 2007-06-15 01:03:22 $ 12 12 * 13 13 * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii … … 127 127 PS_PROJ_SIN, ///< Sine projection 128 128 PS_PROJ_STG, ///< Sine projection 129 PS_PROJ_TNX, ///< Sine projection 129 130 PS_PROJ_ZEA, ///< Sine projection 130 131 PS_PROJ_ZPL, ///< Sine projection
Note:
See TracChangeset
for help on using the changeset viewer.
