Index: trunk/psastro/src/psastroUtils.c =================================================================== --- trunk/psastro/src/psastroUtils.c (revision 9373) +++ trunk/psastro/src/psastroUtils.c (revision 9732) @@ -10,14 +10,14 @@ for (int i = 0; i < rawstars->n; i++) { - pmAstromObj *raw = rawstars->data[i]; - - psPlaneTransformApply (raw->FP, chip->toFPA, raw->chip); - psPlaneDistortApply (raw->TP, fpa->toTangentPlane, raw->FP, 0.0, 0.0); - p_psDeproject (raw->sky, raw->TP, fpa->projection); + pmAstromObj *raw = rawstars->data[i]; + + psPlaneTransformApply (raw->FP, chip->toFPA, raw->chip); + psPlaneDistortApply (raw->TP, fpa->toTangentPlane, raw->FP, 0.0, 0.0); + p_psDeproject (raw->sky, raw->TP, fpa->projection); } for (int i = 0; i < refstars->n; i++) { - pmAstromObj *ref = refstars->data[i]; - psPlaneTransformApply (ref->chip, chip->fromFPA, ref->FP); + pmAstromObj *ref = refstars->data[i]; + psPlaneTransformApply (ref->chip, chip->fromFPA, ref->FP); } @@ -35,22 +35,22 @@ for (int i = 0; i < fpa->chips->n; i++) { - pmChip *chip = fpa->chips->data[i]; - for (int j = 0; j < chip->cells->n; j++) { - pmCell *cell = chip->cells->data[j]; - for (int k = 0; k < cell->readouts->n; k++) { - pmReadout *readout = cell->readouts->data[k]; - - psArray *stars = psMetadataLookupPtr (&status, readout->analysis, "STARS.FULLSET"); - stars = psArraySort (stars, pmAstromObjSortByMag); - - int nSubset = PS_MIN (MAX_NSTARS, stars->n); - psArray *subset = psArrayAlloc (nSubset); - - for (int i = 0; i < nSubset; i++) { - subset->data[i] = stars->data[i]; - } - psMetadataAdd (readout->analysis, PS_LIST_TAIL, "STARS.SUBSET", PS_DATA_ARRAY, "stars from analysis", subset); - } - } + pmChip *chip = fpa->chips->data[i]; + for (int j = 0; j < chip->cells->n; j++) { + pmCell *cell = chip->cells->data[j]; + for (int k = 0; k < cell->readouts->n; k++) { + pmReadout *readout = cell->readouts->data[k]; + + psArray *stars = psMetadataLookupPtr (&status, readout->analysis, "STARS.FULLSET"); + stars = psArraySort (stars, pmAstromObjSortByMag); + + int nSubset = PS_MIN (MAX_NSTARS, stars->n); + psArray *subset = psArrayAlloc (nSubset); + + for (int i = 0; i < nSubset; i++) { + subset->data[i] = stars->data[i]; + } + psMetadataAdd (readout->analysis, PS_LIST_TAIL, "STARS.SUBSET", PS_DATA_ARRAY, "stars from analysis", subset); + } + } } return true; @@ -64,107 +64,107 @@ // first pass: measure the per-chip solutions, modify the chip.toFPA terms for (int i = 0; i < fpa->chips->n; i++) { - pmChip *chip = fpa->chips->data[i]; - psastroChipAstrom (chip, config); - } - - // second stage: re-normalize the chip terms, passing the + pmChip *chip = fpa->chips->data[i]; + psastroChipAstrom (chip, config); + } + + // second stage: re-normalize the chip terms, passing the // average rotation and offset values to the fpa.toSky if (RENORM) { - // this code is needed for the mosastro stage, with multiple chip solutions - - double dX, dY, dT, dN; - dX = dY = dT = dN = 0; - - psPlane origin, P1, P2; - origin.x = 0; - origin.y = 0; - - // calculate the average rotation and boresite offset relative to raw - for (int i = 0; i < fpa->chips->n; i++) { - pmChip *iChip = raw->chips->data[i]; - pmChip *oChip = fpa->chips->data[i]; - - // offset of chip - psCoordChipToFPA (&P1, &origin, iChip); - psCoordChipToFPA (&P2, &origin, oChip); - dX += (P2.x - P1.x); - dY += (P2.y - P1.y); - - // get parity-independent rotations for old and new solutions - double T1 = psPlaneTransformGetRotation (iChip->toFPA); - double T2 = psPlaneTransformGetRotation (oChip->toFPA); - dT += T2 - T1; - dN ++; - } - - dT /= dN; - dX /= dN; - dY /= dN; - - // R(T) - double PC1_1 = fpa->toTPA->x->coeff[1][0][0][0]; - double PC1_2 = fpa->toTPA->x->coeff[0][1][0][0]; - double PC2_1 = fpa->toTPA->y->coeff[1][0][0][0]; - double PC2_2 = fpa->toTPA->y->coeff[0][1][0][0]; - - // R(dT) - double dPC1_1 = +cos (dT); - double dPC1_2 = +sin (dT); - double dPC2_1 = -sin (dT); - double dPC2_2 = +cos (dT); - - // R'(T) = R(T) * R(dT) - double pc1_1 = PC1_1*dPC1_1 + PC1_2*dPC2_1; - double pc1_2 = PC1_1*dPC1_2 + PC1_2*dPC2_2; - double pc2_1 = PC2_1*dPC1_1 + PC2_2*dPC2_1; - double pc2_2 = PC2_1*dPC1_2 + PC2_2*dPC2_2; - - double det = 1.0 / (pc1_1*pc2_2 - pc1_2*pc2_1); - - // R'(-T) (matrix inverse, not just rotation inverse -- keeps parity) - double pi1_1 = +pc2_2 * det; - double pi1_2 = -pc1_2 * det; - double pi2_1 = -pc2_1 * det; - double pi2_2 = +pc1_1 * det; - - // apply the new modifcations in rotation and boresite - for (int i = 0; i < fpa->chips->n; i++) { - pmChip *oChip = fpa->chips->data[i]; - - // r(T) - double pr1_1 = oChip->toFPA->x->coeff[1][0]; - double pr1_2 = oChip->toFPA->x->coeff[0][1]; - double pr2_1 = oChip->toFPA->y->coeff[1][0]; - double pr2_2 = oChip->toFPA->y->coeff[0][1]; - - // ri'(T) = R(T) r(t) - double ri1_1 = PC1_1*pr1_1 + PC1_2*pr2_1; - double ri1_2 = PC1_1*pr1_2 + PC1_2*pr2_2; - double ri2_1 = PC2_1*pr1_1 + PC2_2*pr2_1; - double ri2_2 = PC2_1*pr1_2 + PC2_2*pr2_2; - - // r'(T) = R'(-T) ri'(T) - oChip->toFPA->x->coeff[1][0] = pi1_1*ri1_1 + pi1_2*ri2_1; - oChip->toFPA->x->coeff[0][1] = pi1_1*ri1_2 + pi1_2*ri2_2; - oChip->toFPA->y->coeff[1][0] = pi2_1*ri1_1 + pi2_2*ri2_1; - oChip->toFPA->y->coeff[0][1] = pi2_1*ri1_2 + pi2_2*ri2_2; - - double dx = PC1_1*oChip->toFPA->x->coeff[0][0] + PC1_2*oChip->toFPA->y->coeff[0][0] + dX; - double dy = PC2_1*oChip->toFPA->x->coeff[0][0] + PC2_2*oChip->toFPA->y->coeff[0][0] + dY; - - oChip->toFPA->x->coeff[0][0] = pi1_1*dx + pi1_2*dy; - oChip->toFPA->y->coeff[0][0] = pi2_1*dx + pi2_2*dy; - } - - fpa->toTPA->x->coeff[0][0][0][0] -= dX; - fpa->toTPA->y->coeff[0][0][0][0] -= dY; - - fpa->toTPA->x->coeff[1][0][0][0] = pc1_1; - fpa->toTPA->x->coeff[0][1][0][0] = pc1_2; - fpa->toTPA->y->coeff[1][0][0][0] = pc2_1; - fpa->toTPA->y->coeff[0][1][0][0] = pc2_2; - } - return true; + // this code is needed for the mosastro stage, with multiple chip solutions + + double dX, dY, dT, dN; + dX = dY = dT = dN = 0; + + psPlane origin, P1, P2; + origin.x = 0; + origin.y = 0; + + // calculate the average rotation and boresite offset relative to raw + for (int i = 0; i < fpa->chips->n; i++) { + pmChip *iChip = raw->chips->data[i]; + pmChip *oChip = fpa->chips->data[i]; + + // offset of chip + psCoordChipToFPA (&P1, &origin, iChip); + psCoordChipToFPA (&P2, &origin, oChip); + dX += (P2.x - P1.x); + dY += (P2.y - P1.y); + + // get parity-independent rotations for old and new solutions + double T1 = psPlaneTransformGetRotation (iChip->toFPA); + double T2 = psPlaneTransformGetRotation (oChip->toFPA); + dT += T2 - T1; + dN ++; + } + + dT /= dN; + dX /= dN; + dY /= dN; + + // R(T) + double PC1_1 = fpa->toTPA->x->coeff[1][0][0][0]; + double PC1_2 = fpa->toTPA->x->coeff[0][1][0][0]; + double PC2_1 = fpa->toTPA->y->coeff[1][0][0][0]; + double PC2_2 = fpa->toTPA->y->coeff[0][1][0][0]; + + // R(dT) + double dPC1_1 = +cos (dT); + double dPC1_2 = +sin (dT); + double dPC2_1 = -sin (dT); + double dPC2_2 = +cos (dT); + + // R'(T) = R(T) * R(dT) + double pc1_1 = PC1_1*dPC1_1 + PC1_2*dPC2_1; + double pc1_2 = PC1_1*dPC1_2 + PC1_2*dPC2_2; + double pc2_1 = PC2_1*dPC1_1 + PC2_2*dPC2_1; + double pc2_2 = PC2_1*dPC1_2 + PC2_2*dPC2_2; + + double det = 1.0 / (pc1_1*pc2_2 - pc1_2*pc2_1); + + // R'(-T) (matrix inverse, not just rotation inverse -- keeps parity) + double pi1_1 = +pc2_2 * det; + double pi1_2 = -pc1_2 * det; + double pi2_1 = -pc2_1 * det; + double pi2_2 = +pc1_1 * det; + + // apply the new modifcations in rotation and boresite + for (int i = 0; i < fpa->chips->n; i++) { + pmChip *oChip = fpa->chips->data[i]; + + // r(T) + double pr1_1 = oChip->toFPA->x->coeff[1][0]; + double pr1_2 = oChip->toFPA->x->coeff[0][1]; + double pr2_1 = oChip->toFPA->y->coeff[1][0]; + double pr2_2 = oChip->toFPA->y->coeff[0][1]; + + // ri'(T) = R(T) r(t) + double ri1_1 = PC1_1*pr1_1 + PC1_2*pr2_1; + double ri1_2 = PC1_1*pr1_2 + PC1_2*pr2_2; + double ri2_1 = PC2_1*pr1_1 + PC2_2*pr2_1; + double ri2_2 = PC2_1*pr1_2 + PC2_2*pr2_2; + + // r'(T) = R'(-T) ri'(T) + oChip->toFPA->x->coeff[1][0] = pi1_1*ri1_1 + pi1_2*ri2_1; + oChip->toFPA->x->coeff[0][1] = pi1_1*ri1_2 + pi1_2*ri2_2; + oChip->toFPA->y->coeff[1][0] = pi2_1*ri1_1 + pi2_2*ri2_1; + oChip->toFPA->y->coeff[0][1] = pi2_1*ri1_2 + pi2_2*ri2_2; + + double dx = PC1_1*oChip->toFPA->x->coeff[0][0] + PC1_2*oChip->toFPA->y->coeff[0][0] + dX; + double dy = PC2_1*oChip->toFPA->x->coeff[0][0] + PC2_2*oChip->toFPA->y->coeff[0][0] + dY; + + oChip->toFPA->x->coeff[0][0] = pi1_1*dx + pi1_2*dy; + oChip->toFPA->y->coeff[0][0] = pi2_1*dx + pi2_2*dy; + } + + fpa->toTPA->x->coeff[0][0][0][0] -= dX; + fpa->toTPA->y->coeff[0][0][0][0] -= dY; + + fpa->toTPA->x->coeff[1][0][0][0] = pc1_1; + fpa->toTPA->x->coeff[0][1][0][0] = pc1_2; + fpa->toTPA->y->coeff[1][0][0][0] = pc2_1; + fpa->toTPA->y->coeff[0][1][0][0] = pc2_2; + } + return true; } @@ -174,9 +174,9 @@ for (int i = 0; i < input->nX; i++) { - for (int j = 0; j < input->nY; j++) { - output->mask[i][j] = input->mask[i][j]; - output->coeff[i][j] = input->coeff[i][j]; - output->coeffErr[i][j] = input->coeffErr[i][j]; - } + for (int j = 0; j < input->nY; j++) { + output->mask[i][j] = input->mask[i][j]; + output->coeff[i][j] = input->coeff[i][j]; + output->coeffErr[i][j] = input->coeffErr[i][j]; + } } return (output); @@ -188,13 +188,13 @@ for (int i = 0; i < input->nX; i++) { - for (int j = 0; j < input->nY; j++) { - for (int k = 0; k < input->nZ; k++) { - for (int m = 0; m < input->nT; m++) { - output->mask[i][j][k][m] = input->mask[i][j][k][m]; - output->coeff[i][j][k][m] = input->coeff[i][j][k][m]; - output->coeffErr[i][j][k][m] = input->coeffErr[i][j][k][m]; - } - } - } + for (int j = 0; j < input->nY; j++) { + for (int k = 0; k < input->nZ; k++) { + for (int m = 0; m < input->nT; m++) { + output->mask[i][j][k][m] = input->mask[i][j][k][m]; + output->coeff[i][j][k][m] = input->coeff[i][j][k][m]; + output->coeffErr[i][j][k][m] = input->coeffErr[i][j][k][m]; + } + } + } } return (output); @@ -206,18 +206,18 @@ for (int i = 0; i < input->x->nX; i++) { - for (int j = 0; j < input->x->nY; j++) { - for (int k = 0; k < input->x->nZ; k++) { - for (int m = 0; m < input->x->nT; m++) { - // x-terms - output->x->mask[i][j][k][m] = input->x->mask[i][j][k][m]; - output->x->coeff[i][j][k][m] = input->x->coeff[i][j][k][m]; - output->x->coeffErr[i][j][k][m] = input->x->coeffErr[i][j][k][m]; - // y-terms - output->y->mask[i][j][k][m] = input->y->mask[i][j][k][m]; - output->y->coeff[i][j][k][m] = input->y->coeff[i][j][k][m]; - output->y->coeffErr[i][j][k][m] = input->y->coeffErr[i][j][k][m]; - } - } - } + for (int j = 0; j < input->x->nY; j++) { + for (int k = 0; k < input->x->nZ; k++) { + for (int m = 0; m < input->x->nT; m++) { + // x-terms + output->x->mask[i][j][k][m] = input->x->mask[i][j][k][m]; + output->x->coeff[i][j][k][m] = input->x->coeff[i][j][k][m]; + output->x->coeffErr[i][j][k][m] = input->x->coeffErr[i][j][k][m]; + // y-terms + output->y->mask[i][j][k][m] = input->y->mask[i][j][k][m]; + output->y->coeff[i][j][k][m] = input->y->coeff[i][j][k][m]; + output->y->coeffErr[i][j][k][m] = input->y->coeffErr[i][j][k][m]; + } + } + } } return (output); @@ -229,14 +229,14 @@ for (int i = 0; i < input->x->nX; i++) { - for (int j = 0; j < input->x->nY; j++) { - // x-terms - output->x->mask[i][j] = input->x->mask[i][j]; - output->x->coeff[i][j] = input->x->coeff[i][j]; - output->x->coeffErr[i][j] = input->x->coeffErr[i][j]; - // y-terms - output->y->mask[i][j] = input->y->mask[i][j]; - output->y->coeff[i][j] = input->y->coeff[i][j]; - output->y->coeffErr[i][j] = input->y->coeffErr[i][j]; - } + for (int j = 0; j < input->x->nY; j++) { + // x-terms + output->x->mask[i][j] = input->x->mask[i][j]; + output->x->coeff[i][j] = input->x->coeff[i][j]; + output->x->coeffErr[i][j] = input->x->coeffErr[i][j]; + // y-terms + output->y->mask[i][j] = input->y->mask[i][j]; + output->y->coeff[i][j] = input->y->coeff[i][j]; + output->y->coeffErr[i][j] = input->y->coeffErr[i][j]; + } } return (output); @@ -287,5 +287,5 @@ if (map->y->nX < 1) return 0; if (map->y->nY < 1) return 0; - + double pc1_1 = map->x->coeff[1][0]; double pc1_2 = map->x->coeff[0][1]; @@ -299,5 +299,5 @@ double t1 = -atan2 (px*pc1_2, px*pc1_1); double t2 = +atan2 (py*pc2_1, py*pc2_2); - + // careful near -pi,+pi boundary... if (t1 - t2 > M_PI/2) t2 += 2*M_PI; @@ -307,5 +307,5 @@ while (theta < M_PI) theta += 2*M_PI; while (theta > M_PI) theta -= 2*M_PI; - + return (theta); }