Index: /trunk/psLib/test/astronomy/Makefile =================================================================== --- /trunk/psLib/test/astronomy/Makefile (revision 2095) +++ /trunk/psLib/test/astronomy/Makefile (revision 2096) @@ -3,6 +3,6 @@ ## Makefile: test/astronomy ## -## $Revision: 1.17 $ $Name: not supported by cvs2svn $ -## $Date: 2004-10-12 23:22:05 $ +## $Revision: 1.18 $ $Name: not supported by cvs2svn $ +## $Date: 2004-10-13 23:58:20 $ ## ## Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -28,4 +28,5 @@ tst_psMetadata_05 \ tst_psMetadata_06 \ + tst_psCoord \ tst_psAstrometry \ tst_psAstrometry01 \ Index: /trunk/psLib/test/astronomy/tst_psAstrometry.c =================================================================== --- /trunk/psLib/test/astronomy/tst_psAstrometry.c (revision 2095) +++ /trunk/psLib/test/astronomy/tst_psAstrometry.c (revision 2096) @@ -5,6 +5,6 @@ * @author Robert DeSonia, MHPCC * - * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-10-13 00:10:50 $ + * @version $Revision: 1.8 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-10-13 23:58:20 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -724,5 +724,4 @@ psFree(readout); psFree(readout2); - psFree(image); return 0; Index: /trunk/psLib/test/astronomy/tst_psAstrometry01.c =================================================================== --- /trunk/psLib/test/astronomy/tst_psAstrometry01.c (revision 2095) +++ /trunk/psLib/test/astronomy/tst_psAstrometry01.c (revision 2096) @@ -5,6 +5,6 @@ * @author GLG, MHPCC * -* @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ -* @date $Date: 2004-10-13 21:08:11 $ +* @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ +* @date $Date: 2004-10-13 23:58:20 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -40,71 +40,147 @@ #define HEIGHT 3055.0 #define LAPSERATE 20.0 -#define NUM_CHIPS 2 -#define NUM_CELLS 4 + #define NUM_READOUTS 2 #define READOUT_NUM_ROWS 8 #define READOUT_NUM_COLS 8 + +#define CELL_WIDTH READOUT_NUM_COLS +#define CELL_HEIGHT READOUT_NUM_ROWS + +#define CELL_GAP 2 + +#define NUM_CELLS 4 +#define CELL_MIN_X 0 +#define CELL_MAX_X CELL_WIDTH +#define CELL_MIN_Y 0 +#define CELL_MAX_Y CELL_HEIGHT + +#define NUM_CHIPS 2 +#define CHIP_MIN_X 0 +#define CHIP_MAX_X ((NUM_CELLS * CELL_WIDTH) + ((NUM_CELLS) * CELL_GAP)) +#define CHIP_MIN_Y 0 +#define CHIP_MAX_Y CELL_HEIGHT + +#define CHIP_WIDTH CHIP_MAX_X +#define CHIP_HEIGHT CHIP_MAX_Y +#define CHIP_GAP 2 + +#define FPA_MIN_X 0 +#define FPA_MAX_X ((NUM_CHIPS * CHIP_WIDTH) + ((NUM_CHIPS) * CHIP_GAP)) +#define FPA_MIN_Y 0 +#define FPA_MAX_Y CHIP_HEIGHT + #define ROW0 0 #define COL0 0 -int genSystem() +static void psPlaneDistortFree(psPlaneDistort *pd) +{ + psFree(pd->x); + psFree(pd->y); +} + +static void psPlaneTransformFree(psPlaneTransform *pd) +{ + psFree(pd->x); + psFree(pd->y); +} + +#define PS_CREATE_2D_IDENTITY_PLANE_TRANSFORM(NAME) \ +{\ + NAME = (psPlaneTransform *) psAlloc(sizeof(psPlaneTransform)); \ + p_psMemSetDeallocator(NAME, (psFreeFcn) psPlaneTransformFree); \ + NAME->x = psDPolynomial2DAlloc(2, 2, PS_POLYNOMIAL_ORD); \ + NAME->y = psDPolynomial2DAlloc(2, 2, PS_POLYNOMIAL_ORD); \ + NAME->x->coeff[1][0] = 1.0; \ + NAME->y->coeff[0][1] = 1.0; \ +} + +#define PS_CREATE_4D_IDENTITY_PLANE_DISTORT(NAME) \ +{\ + NAME = (psPlaneDistort *) psAlloc(sizeof(psPlaneDistort)); \ + p_psMemSetDeallocator(NAME, (psFreeFcn) psPlaneDistortFree); \ + NAME->x = psDPolynomial4DAlloc(2, 2, 2, 2, PS_POLYNOMIAL_ORD); \ + NAME->y = psDPolynomial4DAlloc(2, 2, 2, 2, PS_POLYNOMIAL_ORD); \ + NAME->x->coeff[1][0][0][0] = 1.0; \ + NAME->y->coeff[0][1][0][0] = 1.0; \ +} + +/* + This routine will create a system of FPAs/Chips/Cells/Readouts. For + simplicity, an FPA is defined as a linear array of chips, and a chip is + defined as a linear array of cells, both in the x direction. The transforms + between the various layers take into account the cell/chip and the boundaries + between each. +*/ +psFPA *genSystem() { int i; int j; - psTime* now = psTimeGetTime(PS_TIME_UTC); - psObservatory* obs = psObservatoryAlloc(NAME, LATITUDE, LONGITUDE, HEIGHT, LAPSERATE); - psExposure* exp = psExposureAlloc(RA, DEC, HA, ZD, AZ, - now, ROT_ANGLE, TEMPERATURE, PRESSURE, HUMIDITY, - EXP_TIME, WAVELENGTH, obs); + int k; + psImage *tmpImage; + // psTime* now = psTimeGetTime(PS_TIME_UTC); + // psObservatory* obs = psObservatoryAlloc(NAME, LATITUDE, LONGITUDE, HEIGHT, LAPSERATE); + // psExposure* exp = psExposureAlloc(RA, DEC, HA, ZD, AZ, + // now, ROT_ANGLE, TEMPERATURE, PRESSURE, HUMIDITY, + // EXP_TIME, WAVELENGTH, obs); // psGrommit *grom = psGrommitAlloc(exp); - psFPA *myFPA = psFPAAlloc(NUM_CHIPS, exp); + psFPA *myFPA = psFPAAlloc(NUM_CHIPS, NULL); psChip **chips; psCell **cells; - - myFPA->fromTangentPlane = NULL; - myFPA->toTangentPlane = NULL; + psReadout **readouts; + + // We create the transforms between the FPA and tangent plane. The + // transform is a simple identity transform. + PS_CREATE_4D_IDENTITY_PLANE_DISTORT(myFPA->fromTangentPlane); + PS_CREATE_4D_IDENTITY_PLANE_DISTORT(myFPA->toTangentPlane); chips = (psChip **) myFPA->chips->data; for (i=0;itoFPA = NULL; - chips[i]->fromFPA = NULL; + // We create the transforms between the chip and FPA. The + // transform is a simple identity transform. + PS_CREATE_2D_IDENTITY_PLANE_TRANSFORM(chips[i]->toFPA); + chips[i]->toFPA->x->coeff[0][0] = (double) (i * (CHIP_WIDTH + CHIP_GAP)); + chips[i]->toFPA->x->coeff[1][1] = 0.0; + PS_CREATE_2D_IDENTITY_PLANE_TRANSFORM(chips[i]->fromFPA); + chips[i]->fromFPA->x->coeff[0][0] = (double) (-i * (CHIP_WIDTH + CHIP_GAP)); + chips[i]->fromFPA->x->coeff[1][1] = 0.0; + cells = (psCell **) chips[i]->cells->data; for (j=0;jtoChip = NULL; - cells[j]->fromChip = NULL; - cells[j]->toFPA = NULL; - cells[j]->toTP = NULL; - /* - readouts = (psReadout **) cells[j]->readouts->data; - for (k=0;kchips->data; - for (i=0;icells->data; - for (j=0;jreadouts->data; - for (k=0;ktoChip); + cells[j]->toChip->x->coeff[0][0] = (double) (j * (CELL_WIDTH + CELL_GAP)); + cells[j]->toChip->x->coeff[1][1] = 0.0; + PS_CREATE_2D_IDENTITY_PLANE_TRANSFORM(cells[j]->fromChip); + cells[j]->fromChip->x->coeff[0][0] = (double) (-j * (CELL_WIDTH + CELL_GAP)); + cells[j]->fromChip->x->coeff[1][1] = 0.0; + PS_CREATE_2D_IDENTITY_PLANE_TRANSFORM(cells[j]->toFPA); + cells[j]->toFPA->x->coeff[0][0] = (double) ((i * (CHIP_WIDTH + CHIP_GAP)) + + (j * (CELL_WIDTH + CELL_GAP))); + cells[j]->toFPA->x->coeff[1][1] = 0.0; + PS_CREATE_2D_IDENTITY_PLANE_TRANSFORM(cells[j]->toTP); + cells[j]->toTP->x->coeff[0][0] = cells[j]->toFPA->x->coeff[0][0]; + cells[j]->toTP->x->coeff[1][1] = 0.0; + + readouts = (psReadout **) cells[j]->readouts->data; + for (k=0;kdata.F32[row][col] = (float) ((i * (CHIP_WIDTH + CHIP_GAP)) + + (j * (CELL_WIDTH + CELL_GAP))); + } + } + readouts[k] = psReadoutAlloc((j * (CELL_WIDTH + CELL_GAP)), 0, tmpImage); + + } + } + } + + return(myFPA); } @@ -295,7 +371,194 @@ } +int test3a( void ) +{ + int x; + int y; + psPlane fpaCoord; + psFPA *myFPA = genSystem(); + psCell *myCell = NULL; + psPlane chipCoord; + psPlane cellCoord; + psPlane testCoord; + + int xReadout = 0; + int xFPA = 0; + for (int chip=0;chipreadouts->data; + psReadout *readout = readouts[0]; + + printf("-------------------------------\n"); + for (int i = 0; i < readout->image->numRows ; i++) { + for (int j = 0; j < readout->image->numCols ; j++) { + printf("(%.1f) ", readout->image->data.F32[i][j]); + } + printf("\n"); + } + return(0); +} + + +/****************************************************************************** +This routine tests many Astrometry functions: + psChipInFPA() + psCellInFPA() + psCoordFPAToChip() + psCoordChipToCell() +This routine loops through all valid cell coordinates and maintains a running +set of (x,y) coordinates in the FPA plane. It calls the first two functions +with each FPA coordinate and determines the chip/cell that that coordinate +corresponds to. It then calls the last two functions to determine the actual +chip and cell coordinates. The cell coordinate is verified for correctness. + +It then calls + psCoordCellToChip() +to determine if the chip coordinates that were generated agree with the chip +coords returned above. + +Then, it calls + psCellInChip() +and flags an ERROR if the cell it returns is different from that of +psCellInFPA(). + +Then, it calls + psCoordChipToFPA +to determine if the FPA coordinates are correct. + + + *****************************************************************************/ int test3( void ) { - genSystem(); + int x; + int y; + psPlane fpaCoord; + psFPA *myFPA = genSystem(); + psCell *myCell = NULL; + psPlane chipCoord; + psPlane cellCoord; + psPlane testCoord; + + int xReadout = 0; + int xFPA = 0; + for (int chip=0;chipsinPhi)) { + printf("ERROR: myST->sinPhi is %f, should be %f\n", myST->sinPhi, tmpF32); + testStatus = 1; + } + tmpF32 = cos(NPLAT); + if (FLT_EPSILON < fabs(tmpF32 - myST->cosPhi)) { + printf("ERROR: myST->cosPhi is %f, should be %f\n", myST->cosPhi, tmpF32); + testStatus = 2; + } + + if (FLT_EPSILON < fabs(X0_BIG - myST->Xo)) { + printf("ERROR: myST->Xo is %f, should be %f\n", myST->Xo, X0_BIG); + testStatus = 3; + } + + if (FLT_EPSILON < fabs(X0_SMALL - myST->xo)) { + printf("ERROR: myST->xo is %f, should be %f\n", myST->xo, X0_SMALL); + testStatus = 4; + } + + psFree(myST); + return(testStatus); +} + +#define N 10 +// We do a simple identity transformation on a few x,y pairs. int test2( void ) { - // no need to check for memory leaks, as the runTestSuite does that for you. - - // here is where one implements the tests for generally a single testpoint. - - return 0; // the value that indicates success is part of the testDescription -} - + int i; + int testStatus = 0; + psPlane in; + psPlane out; + psPlaneTransform pt; + + pt.x = psDPolynomial2DAlloc(2, 2, PS_POLYNOMIAL_ORD); + pt.y = psDPolynomial2DAlloc(2, 2, PS_POLYNOMIAL_ORD); + pt.x->coeff[1][0] = 1.0; + pt.y->coeff[0][1] = 1.0; + for (i=0;i (%f, %f)\n", in.x, in.y, out.x, out.y); + } + + psFree(pt.x); + psFree(pt.y); + return(testStatus); +} + +#define COLOR 1.0 +#define MAGNITUDE 1.0 +// We do a simple identity transformation on a few x,y pairs. For x and y, +// we add in the COLOR and MAGNITUDE. int test3( void ) { - // no need to check for memory leaks, as the runTestSuite does that for you. - - // here is where one implements the tests for generally a single testpoint. - - return 0; // the value that indicates success is part of the testDescription -} - + int i; + int testStatus = 0; + psPlane in; + psPlane out; + psPlaneDistort pt; + + pt.x = psDPolynomial4DAlloc(2, 2, 2, 2, PS_POLYNOMIAL_ORD); + pt.y = psDPolynomial4DAlloc(2, 2, 2, 2, PS_POLYNOMIAL_ORD); + pt.x->coeff[1][0][0][0] = 1.0; + pt.x->coeff[0][0][1][0] = 1.0; + pt.x->coeff[0][0][0][1] = 1.0; + pt.y->coeff[0][1][0][0] = 1.0; + pt.y->coeff[0][0][1][0] = 1.0; + pt.y->coeff[0][0][0][1] = 1.0; + for (i=0;i (%f, %f)\n", in.x, in.y, out.x, out.y); + } + + psFree(pt.x); + psFree(pt.y); + return(testStatus); +} + +#define DEG_INC 30.0 +// We do a simple identity transformation on a few RA, DEC pairs. +int test4( void ) +{ + int i; + int testStatus = 0; + psSphere in; + psSphere out; + psSphereTransform *myST = psSphereTransformAlloc(0.0, 0.0, 0.0); + + in.r = 45.0 * DEG_TO_RAD; + in.d = 30.0 * DEG_TO_RAD; + in.rErr = 0.0; + in.dErr = 0.0; + + for (float r=0.0;r<180.0;r+=DEG_INC) { + for (float d=0.0;d<90.0;d+=DEG_INC) { + in.r = r * DEG_TO_RAD; + in.d = d * DEG_TO_RAD; + in.rErr = 0.0; + in.dErr = 0.0; + + psSphereTransformApply(&out, myST, &in); + + if (FLT_EPSILON < fabs(out.r - in.r)) { + printf("ERROR: out.r is %f, should be %f\n", out.r, in.r); + testStatus = 4; + } + + if (FLT_EPSILON < fabs(out.d - in.d)) { + printf("ERROR: out.d is %f, should be %f\n", out.d, in.d); + testStatus = 5; + } + // printf("psSphereTransformApply (%f, %f) -> (%f, %f)\n", out.r, out.d, in.r, in.d); + } + } + + psFree(myST); + return(testStatus); +} + +int test5( void ) +{ + int testStatus = 0; + psTime* now = psTimeGetTime(PS_TIME_UTC); + struct tm *tm_time = psTimeToTM(now); + + // XXX: This test code is simply a copy of the original source code. + double year = (double)(1900 + tm_time->tm_year); + double T = year / 100.0; + double phi = -23.452294 + 0.013013 * T + 0.000001639 * T * T - 0.000000503 * T * T * T; + psSphereTransform *myST = psSphereTransformICRStoEcliptic(*now); + + if (FLT_EPSILON < fabs(sin(phi) - myST->sinPhi)) { + printf("ERROR: myST->sinPhi is %f, should be %f\n", myST->sinPhi, sin(phi)); + testStatus = 1; + } + + if (FLT_EPSILON < fabs(cos(phi) - myST->cosPhi)) { + printf("ERROR: myST->cosPhi is %f, should be %f\n", myST->cosPhi, cos(phi)); + testStatus = 2; + } + + if (FLT_EPSILON < fabs(0.0 - myST->Xo)) { + printf("ERROR: myST->Xo is %f, should be %f\n", myST->Xo, X0_BIG); + testStatus = 3; + } + + if (FLT_EPSILON < fabs(0.0 - myST->xo)) { + printf("ERROR: myST->xo is %f, should be %f\n", myST->xo, X0_SMALL); + testStatus = 4; + } + + psFree(myST); + psFree(tm_time); + psFree(now); + return(testStatus); +} + + +int test6( void ) +{ + int testStatus = 0; + psTime* now = psTimeGetTime(PS_TIME_UTC); + struct tm *tm_time = psTimeToTM(now); + + // XXX: This test code is simply a copy of the original source code. + double year = (double)(1900 + tm_time->tm_year); + double T = year / 100.0; + double phi = +23.452294 - 0.013013 * T - 0.000001639 * T * T + 0.000000503 * T * T * T; + psSphereTransform *myST = psSphereTransformEcliptictoICRS(*now); + + if (FLT_EPSILON < fabs(sin(phi) - myST->sinPhi)) { + printf("ERROR: myST->sinPhi is %f, should be %f\n", myST->sinPhi, sin(phi)); + testStatus = 1; + } + + if (FLT_EPSILON < fabs(cos(phi) - myST->cosPhi)) { + printf("ERROR: myST->cosPhi is %f, should be %f\n", myST->cosPhi, cos(phi)); + testStatus = 2; + } + + if (FLT_EPSILON < fabs(0.0 - myST->Xo)) { + printf("ERROR: myST->Xo is %f, should be %f\n", myST->Xo, X0_BIG); + testStatus = 3; + } + + if (FLT_EPSILON < fabs(0.0 - myST->xo)) { + printf("ERROR: myST->xo is %f, should be %f\n", myST->xo, X0_SMALL); + testStatus = 4; + } + + psFree(myST); + psFree(tm_time); + psFree(now); + return(testStatus); +} + + +int test7( void ) +{ + int testStatus = 0; + // XXX: This test code is simply a copy of the original source code. + double phi = 62.6; + double Xo = 282.25; + double xo = 33.0; + psSphereTransform *myST = psSphereTransformICRStoGalatic(); + + if (FLT_EPSILON < fabs(sin(phi) - myST->sinPhi)) { + printf("ERROR: myST->sinPhi is %f, should be %f\n", myST->sinPhi, sin(phi)); + testStatus = 1; + } + + if (FLT_EPSILON < fabs(cos(phi) - myST->cosPhi)) { + printf("ERROR: myST->cosPhi is %f, should be %f\n", myST->cosPhi, cos(phi)); + testStatus = 2; + } + + if (FLT_EPSILON < fabs(Xo - myST->Xo)) { + printf("ERROR: myST->Xo is %f, should be %f\n", myST->Xo, Xo); + testStatus = 3; + } + + if (FLT_EPSILON < fabs(xo - myST->xo)) { + printf("ERROR: myST->xo is %f, should be %f\n", myST->xo, xo); + testStatus = 4; + } + + psFree(myST); + return(testStatus); +} + +int test8( void ) +{ + int testStatus = 0; + // XXX: This test code is simply a copy of the original source code. + double phi = -62.6; + double Xo = 33.0; + double xo = 282.25; + psSphereTransform *myST = psSphereTransformGalatictoICRS(); + + if (FLT_EPSILON < fabs(sin(phi) - myST->sinPhi)) { + printf("ERROR: myST->sinPhi is %f, should be %f\n", myST->sinPhi, sin(phi)); + testStatus = 1; + } + + if (FLT_EPSILON < fabs(cos(phi) - myST->cosPhi)) { + printf("ERROR: myST->cosPhi is %f, should be %f\n", myST->cosPhi, cos(phi)); + testStatus = 2; + } + + if (FLT_EPSILON < fabs(Xo - myST->Xo)) { + printf("ERROR: myST->Xo is %f, should be %f\n", myST->Xo, Xo); + testStatus = 3; + } + + if (FLT_EPSILON < fabs(xo - myST->xo)) { + printf("ERROR: myST->xo is %f, should be %f\n", myST->xo, xo); + testStatus = 4; + } + + psFree(myST); + return(testStatus); +} + +int test20( void ) +{ + int i; + int testStatus = 0; + psSphere in; + double R; + double expectX; + double expectY; + psPlane *out; + psProjection myProjection; + myProjection.R = 20.0 * DEG_TO_RAD; + myProjection.D = 10.0 * DEG_TO_RAD; + myProjection.Xs = 1.0; + myProjection.Ys = 1.0; + + for (float r=0.0;r<180.0;r+=DEG_INC) { + for (float d=0.0;d<90.0;d+=DEG_INC) { + in.r = r * DEG_TO_RAD; + in.d = d * DEG_TO_RAD; + in.rErr = 0.0; + in.dErr = 0.0; + + /****************************************************************** + Tangent Plane Projection + *****************************************************************/ + myProjection.type = PS_PROJ_TAN; + out = psProject(&in, &myProjection); + + R = PS_COT(in.r) * (180.0 / M_PI); + expectX = R * sin(in.d); + expectY = R * cos(in.d); + PS_COMPARE_TINY_THEN_PRINT_ERROR(out->x, expectX, 1); + PS_COMPARE_TINY_THEN_PRINT_ERROR(out->y, expectY, 2); + psFree(out); + + /****************************************************************** + PS_PROJ_SIN Projection + *****************************************************************/ + myProjection.type = PS_PROJ_SIN; + out = psProject(&in, &myProjection); + + R = cos(in.r) * (180.0 / M_PI); + expectX = R * sin(in.d); + expectY = R * cos(in.d); + PS_COMPARE_TINY_THEN_PRINT_ERROR(out->x, expectX, 3); + PS_COMPARE_TINY_THEN_PRINT_ERROR(out->y, expectY, 4); + psFree(out); + + /****************************************************************** + PS_PROJ_CAR Projection + *****************************************************************/ + myProjection.type = PS_PROJ_CAR; + out = psProject(&in, &myProjection); + + expectX = in.d; + expectY = in.r; + PS_COMPARE_TINY_THEN_PRINT_ERROR(out->x, expectX, 5); + PS_COMPARE_TINY_THEN_PRINT_ERROR(out->y, expectY, 6); + psFree(out); + + /****************************************************************** + PS_PROJ_MER Projection + *****************************************************************/ + myProjection.type = PS_PROJ_MER; + out = psProject(&in, &myProjection); + + expectX = in.d; + expectY = log(tan(45.0 + (0.5 * in.r))) * 180.0 / M_PI; + + PS_COMPARE_TINY_THEN_PRINT_ERROR(out->x, expectX, 7); + PS_COMPARE_TINY_THEN_PRINT_ERROR(out->y, expectY, 8); + psFree(out); + + /****************************************************************** + PS_PROJ_AIT Projection + *****************************************************************/ + myProjection.type = PS_PROJ_AIT; + out = psProject(&in, &myProjection); + + double alpha = 1.0 / ((180.0 / M_PI) * sqrt(1.0 + (cos(in.r) * + cos(0.5 * in.d) * 0.5))); + expectX = 2.0 * alpha * cos(in.r) * sin(0.5 * in.d); + expectY = alpha * sin(in.d); + + PS_COMPARE_TINY_THEN_PRINT_ERROR(out->x, expectX, 7); + PS_COMPARE_TINY_THEN_PRINT_ERROR(out->y, expectY, 8); + psFree(out); + + + } + } + return(testStatus); +} + +#define SPACE_INC 10.0 +int test21( void ) +{ + int testStatus = 0; + psPlane in; + double R; + double expectD; + double expectR; + psSphere *out; + psProjection myProjection; + myProjection.R = 20.0 * DEG_TO_RAD; + myProjection.D = 10.0 * DEG_TO_RAD; + myProjection.Xs = 1.0; + myProjection.Ys = 1.0; + + in.xErr = 0.0; + in.yErr = 0.0; + for (in.x=0.0; in.x<100.0; in.x+=SPACE_INC) { + for (in.y=0.0;in.y<100.0;in.y+=SPACE_INC) { + + /****************************************************************** + Tangent Plane Projection + *****************************************************************/ + myProjection.type = PS_PROJ_TAN; + out = psDeproject(&in, &myProjection); + + R = sqrt((in.x * in.x) + (in.y * in.y)); + expectD = p_psArg(-in.y, in.x); + expectR = atan(180.0 / (R * M_PI)); + + PS_COMPARE_TINY_THEN_PRINT_ERROR(out->r, expectR, 1); + PS_COMPARE_TINY_THEN_PRINT_ERROR(out->d, expectD, 2); + psFree(out); + + /****************************************************************** + PS_PROJ_SIN Projection + *****************************************************************/ + myProjection.type = PS_PROJ_SIN; + out = psDeproject(&in, &myProjection); + + R = sqrt((in.x * in.x) + (in.y * in.y)); + expectD = p_psArg(-in.y, in.x); + expectR = acos((R * M_PI) / 180.0); + + PS_COMPARE_TINY_THEN_PRINT_ERROR(out->r, expectR, 3); + PS_COMPARE_TINY_THEN_PRINT_ERROR(out->d, expectD, 4); + psFree(out); + + /****************************************************************** + PS_PROJ_CAR Projection + *****************************************************************/ + myProjection.type = PS_PROJ_CAR; + out = psDeproject(&in, &myProjection); + + expectD = in.x; + expectR = in.y; + + PS_COMPARE_TINY_THEN_PRINT_ERROR(out->r, expectR, 3); + PS_COMPARE_TINY_THEN_PRINT_ERROR(out->d, expectD, 4); + psFree(out); + + /****************************************************************** + PS_PROJ_MER Projection + *****************************************************************/ + myProjection.type = PS_PROJ_MER; + out = psDeproject(&in, &myProjection); + + expectD = in.x; + expectR = (2.0 * atan(exp((in.y * M_PI / 180.0)))) - 180.0; + + PS_COMPARE_TINY_THEN_PRINT_ERROR(out->r, expectR, 3); + PS_COMPARE_TINY_THEN_PRINT_ERROR(out->d, expectD, 4); + psFree(out); + + + /****************************************************************** + PS_PROJ_AIT Projection + *****************************************************************/ + myProjection.type = PS_PROJ_AIT; + out = psDeproject(&in, &myProjection); + + float chu1 = (in.x * M_PI) / 720.0; + chu1 *= chu1; + float chu2 = (in.y * M_PI) / 360.0; + chu2 *= chu2; + float chu = sqrt(1.0 - chu1 - chu2); + expectD = 2.0 * p_psArg((2.0 * chu * chu) - 1.0, (in.x * chu * M_PI) / 360.0); + expectR = asin((in.y * chu * M_PI) / 180.0); + + PS_COMPARE_TINY_THEN_PRINT_ERROR(out->r, expectR, 3); + PS_COMPARE_TINY_THEN_PRINT_ERROR(out->d, expectD, 4); + psFree(out); + } + } + return(testStatus); +} + +int test40( void ) +{ + int testStatus = 0; + psSphere position1; + psSphere position2; + psSphere *offset = NULL; + + position1.r = 90.0 * DEG_TO_RAD; + position1.d = 45.0 * DEG_TO_RAD; + position1.rErr = 0.0; + position1.dErr = 0.0; + + for (float r=0.0;r<180.0;r+=DEG_INC) { + for (float d=0.0;d<90.0;d+=DEG_INC) { + position2.r = r * DEG_TO_RAD; + position2.d = d * DEG_TO_RAD; + position2.rErr = 0.0; + position2.dErr = 0.0; + + offset = psSphereGetOffset(&position1, &position2, + PS_SPHERICAL, PS_RADIAN); + PS_COMPARE_TINY_THEN_PRINT_ERROR(offset->r, (position2.r - position1.r), 1); + PS_COMPARE_TINY_THEN_PRINT_ERROR(offset->d, (position2.d - position1.d), 1); + psFree(offset); + + offset = psSphereGetOffset(&position1, &position2, + PS_SPHERICAL, PS_DEGREE); + offset->r = offset->r * DEG_TO_RAD; + offset->d = offset->d * DEG_TO_RAD; + PS_COMPARE_TINY_THEN_PRINT_ERROR(offset->r, (position2.r - position1.r), 2); + PS_COMPARE_TINY_THEN_PRINT_ERROR(offset->d, (position2.d - position1.d), 3); + psFree(offset); + + offset = psSphereGetOffset(&position1, &position2, + PS_SPHERICAL, PS_ARCMIN); + offset->r = offset->r * ARCMIN_TO_RAD; + offset->d = offset->d * ARCMIN_TO_RAD; + PS_COMPARE_TINY_THEN_PRINT_ERROR(offset->r, (position2.r - position1.r), 2); + PS_COMPARE_TINY_THEN_PRINT_ERROR(offset->d, (position2.d - position1.d), 3); + psFree(offset); + + offset = psSphereGetOffset(&position1, &position2, + PS_SPHERICAL, PS_ARCSEC); + offset->r = offset->r * ARCSEC_TO_RAD; + offset->d = offset->d * ARCSEC_TO_RAD; + PS_COMPARE_TINY_THEN_PRINT_ERROR(offset->r, (position2.r - position1.r), 2); + PS_COMPARE_TINY_THEN_PRINT_ERROR(offset->d, (position2.d - position1.d), 3); + psFree(offset); + + /* XXX: This code does not work correctly. + offset = psSphereGetOffset(&position1, &position2, + PS_LINEAR, 0); + printf("--------------- (%f, %f) to (%f, %f) is (%f, %f) ---------------\n", + position1.r, position1.d, position2.r, position2.d, + offset->r, offset->d); + PS_COMPARE_TINY_THEN_PRINT_ERROR(offset->r, (position2.r - position1.r), 1); + PS_COMPARE_TINY_THEN_PRINT_ERROR(offset->d, (position2.d - position1.d), 1); + psFree(offset); + */ + } + } + + return(testStatus); +} + +int test41( void ) +{ + int testStatus = 0; + psSphere position1; + psSphere *position2; + psSphere offset; + psSphere tmpOffset; + + position1.r = 90.0 * DEG_TO_RAD; + position1.d = 45.0 * DEG_TO_RAD; + position1.rErr = 0.0; + position1.dErr = 0.0; + + for (float r=0.0;r<180.0;r+=DEG_INC) { + for (float d=0.0;d<90.0;d+=DEG_INC) { + offset.r = r * DEG_TO_RAD; + offset.d = d * DEG_TO_RAD; + offset.rErr = 0.0; + offset.dErr = 0.0; + + position2 = psSphereSetOffset(&position1, &offset, + PS_SPHERICAL, PS_RADIAN); + PS_COMPARE_TINY_THEN_PRINT_ERROR(position2->r, (position1.r + offset.r), 1); + PS_COMPARE_TINY_THEN_PRINT_ERROR(position2->d, (position1.d + offset.d), 1); + psFree(position2); + + tmpOffset.r = offset.r * RAD_TO_DEG; + tmpOffset.d = offset.d * RAD_TO_DEG; + tmpOffset.rErr = 0.0; + tmpOffset.dErr = 0.0; + position2 = psSphereSetOffset(&position1, &tmpOffset, + PS_SPHERICAL, PS_DEGREE); + PS_COMPARE_TINY_THEN_PRINT_ERROR(position2->r, (position1.r + offset.r), 1); + PS_COMPARE_TINY_THEN_PRINT_ERROR(position2->d, (position1.d + offset.d), 1); + psFree(position2); + + tmpOffset.r = offset.r * RAD_TO_ARCMIN; + tmpOffset.d = offset.d * RAD_TO_ARCMIN; + tmpOffset.rErr = 0.0; + tmpOffset.dErr = 0.0; + position2 = psSphereSetOffset(&position1, &tmpOffset, + PS_SPHERICAL, PS_ARCMIN); + PS_COMPARE_TINY_THEN_PRINT_ERROR(position2->r, (position1.r + offset.r), 1); + PS_COMPARE_TINY_THEN_PRINT_ERROR(position2->d, (position1.d + offset.d), 1); + psFree(position2); + + tmpOffset.r = offset.r * RAD_TO_ARCSEC; + tmpOffset.d = offset.d * RAD_TO_ARCSEC; + tmpOffset.rErr = 0.0; + tmpOffset.dErr = 0.0; + position2 = psSphereSetOffset(&position1, &tmpOffset, + PS_SPHERICAL, PS_ARCSEC); + PS_COMPARE_TINY_THEN_PRINT_ERROR(position2->r, (position1.r + offset.r), 1); + PS_COMPARE_TINY_THEN_PRINT_ERROR(position2->d, (position1.d + offset.d), 1); + psFree(position2); + + /* XXX: This code does not work correctly. + position2 = psSphereSetOffset(&position1, &offset, + PS_LINEAR, 0); + printf("--------------- (%f, %f) and (%f, %f) is (%f, %f) ---------------\n", + position1.r, position1.d, offset.r, offset.d, + position2->r, position2->d); + PS_COMPARE_TINY_THEN_PRINT_ERROR(position2->r, (position1.r + offset.r), 1); + PS_COMPARE_TINY_THEN_PRINT_ERROR(position2->d, (position1.d + offset.d), 1); + psFree(position2); + */ + } + } + + return(testStatus); +}