Index: /trunk/psastro/src/psastroChooseRefstars.c =================================================================== --- /trunk/psastro/src/psastroChooseRefstars.c (revision 9731) +++ /trunk/psastro/src/psastroChooseRefstars.c (revision 9732) @@ -11,6 +11,6 @@ psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, "PSASTRO"); if (!recipe) { - psError(PSASTRO_ERR_CONFIG, true, "Can't find PSASTRO recipe!\n"); - return false; + psError(PSASTRO_ERR_CONFIG, true, "Can't find PSASTRO recipe!\n"); + return false; } @@ -18,9 +18,9 @@ pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT"); if (!input) { - psError(PSASTRO_ERR_CONFIG, true, "Can't find input data!\n"); - return false; + psError(PSASTRO_ERR_CONFIG, true, "Can't find input data!\n"); + return false; } - float fieldExtra = psMetadataLookupS32 (&status, recipe, "PSASTRO.FIELD.EXTRA"); + float fieldExtra = psMetadataLookupS32 (&status, recipe, "PSASTRO.FIELD.EXTRA"); if (!status) fieldExtra = 0.0; @@ -32,56 +32,56 @@ psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process); if (!chip->process || !chip->file_exists) { continue; } - - while ((cell = pmFPAviewNextCell (view, fpa, 1)) != NULL) { + + while ((cell = pmFPAviewNextCell (view, fpa, 1)) != NULL) { psTrace ("psastro", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); if (!cell->process || !cell->file_exists) { continue; } - // process each of the readouts - // XXX there can only be one readout per chip in astrometry, right? - while ((readout = pmFPAviewNextReadout (view, fpa, 1)) != NULL) { - if (! readout->data_exists) { continue; } + // process each of the readouts + // XXX there can only be one readout per chip in astrometry, right? + while ((readout = pmFPAviewNextReadout (view, fpa, 1)) != NULL) { + if (! readout->data_exists) { continue; } - // read WCS data from the corresponding header - pmHDU *hdu = pmFPAviewThisHDU (view, fpa); + // read WCS data from the corresponding header + pmHDU *hdu = pmFPAviewThisHDU (view, fpa); - int Nx = psMetadataLookupS32 (&status, hdu->header, "NAXIS1"); - int Ny = psMetadataLookupS32 (&status, hdu->header, "NAXIS2"); + int Nx = psMetadataLookupS32 (&status, hdu->header, "NAXIS1"); + int Ny = psMetadataLookupS32 (&status, hdu->header, "NAXIS2"); - float minX = -fieldExtra*Nx; - float maxX = (1+fieldExtra)*Nx; - float minY = -fieldExtra*Ny; - float maxY = (1+fieldExtra)*Ny; + float minX = -fieldExtra*Nx; + float maxX = (1+fieldExtra)*Nx; + float minY = -fieldExtra*Ny; + float maxY = (1+fieldExtra)*Ny; - // the refstars is a subset within range of this chip - psArray *refstars = psArrayAlloc (100); + // the refstars is a subset within range of this chip + psArray *refstars = psArrayAllocEmpty (100); - // select the reference objects within range of this readout - // project the reference objects to this chip - for (int i = 0; i < refs->n; i++) { - pmAstromObj *ref = pmAstromObjCopy(refs->data[i]); - - // XXX why is this still a private function? - p_psProject (ref->TP, ref->sky, fpa->projection); - psPlaneDistortApply (ref->FP, fpa->fromTangentPlane, ref->TP, 0.0, 0.0); - psPlaneTransformApply (ref->chip, chip->fromFPA, ref->FP); + // select the reference objects within range of this readout + // project the reference objects to this chip + for (int i = 0; i < refs->n; i++) { + pmAstromObj *ref = pmAstromObjCopy(refs->data[i]); - // limit the X,Y range of the refs to the selected chip - if (ref->chip->x < minX) goto skip; - if (ref->chip->x > maxX) goto skip; - if (ref->chip->y < minY) goto skip; - if (ref->chip->y > maxY) goto skip; + // XXX why is this still a private function? + p_psProject (ref->TP, ref->sky, fpa->projection); + psPlaneDistortApply (ref->FP, fpa->fromTangentPlane, ref->TP, 0.0, 0.0); + psPlaneTransformApply (ref->chip, chip->fromFPA, ref->FP); - psArrayAdd (refstars, 100, ref); - skip: - psFree (ref); - } - psTrace ("psastro", 4, "Added %ld refstars\n", refstars->n); + // limit the X,Y range of the refs to the selected chip + if (ref->chip->x < minX) goto skip; + if (ref->chip->x > maxX) goto skip; + if (ref->chip->y < minY) goto skip; + if (ref->chip->y > maxY) goto skip; - psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.REFSTARS", PS_DATA_ARRAY, "astrometry matches", refstars); - psFree (refstars); + psArrayAdd (refstars, 100, ref); + skip: + psFree (ref); + } + psTrace ("psastro", 4, "Added %ld refstars\n", refstars->n); - psastroRefstarSubset (readout); - } - } + psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.REFSTARS", PS_DATA_ARRAY, "astrometry matches", refstars); + psFree (refstars); + + psastroRefstarSubset (readout); + } + } } Index: /trunk/psastro/src/psastroConvert.c =================================================================== --- /trunk/psastro/src/psastroConvert.c (revision 9731) +++ /trunk/psastro/src/psastroConvert.c (revision 9732) @@ -8,20 +8,20 @@ pmReadout *readout; pmFPAview *view = pmFPAviewAlloc (0); - + while ((chip = pmFPAviewNextChip (view, fpa, 1)) != NULL) { psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process); if (!chip->process || !chip->file_exists) { continue; } - while ((cell = pmFPAviewNextCell (view, fpa, 1)) != NULL) { + while ((cell = pmFPAviewNextCell (view, fpa, 1)) != NULL) { psTrace ("psastro", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); if (!cell->process || !cell->file_exists) { continue; } - // process each of the readouts - while ((readout = pmFPAviewNextReadout (view, fpa, 1)) != NULL) { - if (! readout->data_exists) { continue; } + // process each of the readouts + while ((readout = pmFPAviewNextReadout (view, fpa, 1)) != NULL) { + if (! readout->data_exists) { continue; } - psastroConvertReadout (readout, recipe); - } - } + psastroConvertReadout (readout, recipe); + } + } } psFree (view); @@ -40,5 +40,5 @@ psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.RAWSTARS", PS_DATA_ARRAY, "astrometry objects", rawstars); psFree (rawstars); - + return true; } @@ -47,31 +47,31 @@ psArray *pmSourceToAstromObj (psArray *sources) { - psArray *objects = psArrayAlloc (sources->n); - + psArray *objects = psArrayAllocEmpty (sources->n); + for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - - // only accept the PSF sources? - if (source->type != PM_SOURCE_TYPE_STAR) continue; + pmSource *source = sources->data[i]; - pmModel *model = source->modelPSF; - if (model == NULL) continue; + // only accept the PSF sources? + if (source->type != PM_SOURCE_TYPE_STAR) continue; - psF32 *PAR = model->params->data.F32; + pmModel *model = source->modelPSF; + if (model == NULL) continue; - pmAstromObj *obj = pmAstromObjAlloc (); + psF32 *PAR = model->params->data.F32; - // is the source magnitude calibrated in any sense? - obj->pix->x = PAR[2]; - obj->pix->y = PAR[3]; - obj->Mag = source->psfMag; + pmAstromObj *obj = pmAstromObjAlloc (); - // XXX do we have the information giving the readout and cell offset? - // for the moment, assume chip == cell == readout - *obj->cell = *obj->pix; - *obj->chip = *obj->cell; + // is the source magnitude calibrated in any sense? + obj->pix->x = PAR[2]; + obj->pix->y = PAR[3]; + obj->Mag = source->psfMag; - psArrayAdd (objects, 100, obj); - psFree (obj); + // XXX do we have the information giving the readout and cell offset? + // for the moment, assume chip == cell == readout + *obj->cell = *obj->pix; + *obj->chip = *obj->cell; + + psArrayAdd (objects, 100, obj); + psFree (obj); } return objects; Index: /trunk/psastro/src/psastroLoadRefstars.c =================================================================== --- /trunk/psastro/src/psastroLoadRefstars.c (revision 9731) +++ /trunk/psastro/src/psastroLoadRefstars.c (revision 9732) @@ -24,5 +24,5 @@ char *CATDIR = psMetadataLookupStr(NULL, recipe, "DVO.CATDIR"); if (CATDIR == NULL) { - psLogMsg ("psastro", 2, "warning: missing DVO.CATDIR, using ~/.ptolemyrc definition\n"); + psLogMsg ("psastro", 2, "warning: missing DVO.CATDIR, using ~/.ptolemyrc definition\n"); } @@ -33,5 +33,5 @@ sprintf (tempFile, "/tmp/psastro.XXXXXX"); if ((fd = mkstemp (tempFile)) == -1) { - psError(PSASTRO_ERR_REFSTARS, true, "error creating temp output file\n"); + psError(PSASTRO_ERR_REFSTARS, true, "error creating temp output file\n"); return NULL; } @@ -42,12 +42,12 @@ // use fork to add timeout capability if (ELIXIR_MODE) { - psStringAppend (&catformat, "-D CATFORMAT elixir "); + psStringAppend (&catformat, "-D CATFORMAT elixir "); } else { - psStringAppend (&catformat, "-D CATFORMAT panstarrs "); - } + psStringAppend (&catformat, "-D CATFORMAT panstarrs "); + } if (CATDIR) { - psStringAppend (&catformat, "-D CATDIR %s ", CATDIR); - } + psStringAppend (&catformat, "-D CATDIR %s ", CATDIR); + } psStringAppend (&getstarLine, "getstar %s -D CATMODE mef -maglim %f -region %f %f %f %f -o %s", catformat, MAGmax, RAmin, DECmin, RAmax, DECmax, tempFile); @@ -57,5 +57,5 @@ status = system (getstarLine); if (status) { - psError(PSASTRO_ERR_REFSTARS, true, "error loading reference data\n"); + psError(PSASTRO_ERR_REFSTARS, true, "error loading reference data\n"); return NULL; } @@ -69,7 +69,7 @@ if (ELIXIR_MODE) { - psFitsMoveExtName (fits, "DVO_AVERAGE_ELIXIR"); + psFitsMoveExtName (fits, "DVO_AVERAGE_ELIXIR"); } else { - psFitsMoveExtName (fits, "DVO_AVERAGE_PANSTARRS"); + psFitsMoveExtName (fits, "DVO_AVERAGE_PANSTARRS"); } @@ -84,5 +84,5 @@ // convert the Average table to the pmAstromObj entries psTimerStart ("psastro"); - psArray *refs = psArrayAlloc (table->n); + psArray *refs = psArrayAllocEmpty (table->n); for (int i = 0; i < table->n; i++) { pmAstromObj *ref = pmAstromObjAlloc (); @@ -91,13 +91,13 @@ // DVO tables are stored in degrees - if (ELIXIR_MODE) { - ref->sky->r = RAD_DEG*psMetadataLookupF32 (&status, row, "RA"); - ref->sky->d = RAD_DEG*psMetadataLookupF32 (&status, row, "DEC"); - ref->Mag = 0.001*psMetadataLookupS32 (&status, row, "MAG"); // ELIXIR uses millimags + if (ELIXIR_MODE) { + ref->sky->r = RAD_DEG*psMetadataLookupF32 (&status, row, "RA"); + ref->sky->d = RAD_DEG*psMetadataLookupF32 (&status, row, "DEC"); + ref->Mag = 0.001*psMetadataLookupS32 (&status, row, "MAG"); // ELIXIR uses millimags } else { - ref->sky->r = RAD_DEG*psMetadataLookupF64 (&status, row, "RA"); - ref->sky->d = RAD_DEG*psMetadataLookupF64 (&status, row, "DEC"); - ref->Mag = psMetadataLookupF32 (&status, row, "MAG"); // PANSTARRS uses mags - } + ref->sky->r = RAD_DEG*psMetadataLookupF64 (&status, row, "RA"); + ref->sky->d = RAD_DEG*psMetadataLookupF64 (&status, row, "DEC"); + ref->Mag = psMetadataLookupF32 (&status, row, "MAG"); // PANSTARRS uses mags + } psArrayAdd (refs, 100, ref); Index: /trunk/psastro/src/psastroRefstarSubset.c =================================================================== --- /trunk/psastro/src/psastroRefstarSubset.c (revision 9731) +++ /trunk/psastro/src/psastroRefstarSubset.c (revision 9732) @@ -32,5 +32,5 @@ return NULL; } - + // XXXX test // psFree (rawfunc); @@ -38,5 +38,5 @@ // return true; - // what is the offset between the two lines at the average magnitude? + // what is the offset between the two lines at the average magnitude? double mRef = 0.5*(reffunc->mMin + reffunc->mMax); double logRho = mRef * reffunc->slope + reffunc->offset; @@ -48,5 +48,5 @@ psLogMsg ("psastro", 4, "clipping stars fainter than %f\n", mRefMax); - psArray *subset = psArrayAlloc (100); + psArray *subset = psArrayAllocEmpty (100); for (int i = 0; i < refstars->n; i++) { pmAstromObj *ref = refstars->data[i]; Index: /trunk/psastro/src/psastroUtils.c =================================================================== --- /trunk/psastro/src/psastroUtils.c (revision 9731) +++ /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); }