Changeset 6176
- Timestamp:
- Jan 22, 2006, 9:53:37 AM (20 years ago)
- Location:
- trunk/psastro
- Files:
-
- 3 added
- 3 edited
-
doc/notes.txt (modified) (1 diff)
-
src/psastroBuildFPA.c (modified) (1 diff)
-
src/psastroChipAstrom.c (added)
-
src/psastroConvert.c (added)
-
src/psastroUtils.c (modified) (2 diffs)
-
src/psastroWCS.c (added)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psastro/doc/notes.txt
r5510 r6176 10 10 - in the first stage, the per-chip results are applied to the 11 11 chip.toFPA parameter set. 12 12 13 - in the second stage, the chip results are collectively used to 13 14 modify the fpa.toSky terms. -
trunk/psastro/src/psastroBuildFPA.c
r5575 r6176 68 68 return (fpa); 69 69 } 70 71 // interpret header WCS (only handles traditional WCS for the moment)72 bool pmAstromReadWCS (psPlaneTransform **toFPAOut, psProjection **toSkyOut, psMetadata *header) {73 74 psPlaneTransform *toFPA;75 psProjection *toSky;76 psProjectionType type;77 bool status;78 float crval1, crval2, crpix1, crpix2, cdelt1, cdelt2;79 float pc1_1, pc1_2, pc2_1, pc2_2;80 81 // *** interpret header data, convert to crval(i), etc82 char *ctype = psMetadataLookupPtr (&status, header, "CTYPE2");83 if (!status) {84 psLogMsg ("psastro", 2, "warning: no WCS metadata in header\n");85 return false;86 }87 88 // determine projection type89 type = PS_PROJ_NTYPE;90 if (!strcmp (&ctype[4], "-SIN")) type = PS_PROJ_SIN;91 if (!strcmp (&ctype[4], "-TAN")) type = PS_PROJ_TAN;92 if (!strcmp (&ctype[4], "-AIT")) type = PS_PROJ_AIT;93 if (!strcmp (&ctype[4], "-PAR")) type = PS_PROJ_PAR;94 if (type == PS_PROJ_NTYPE) {95 psLogMsg ("psastro", 2, "warning: unknown projection type %s\n", ctype);96 return false;97 }98 99 crval1 = psMetadataLookupF32 (&status, header, "CRVAL1");100 crval2 = psMetadataLookupF32 (&status, header, "CRVAL2");101 crpix1 = psMetadataLookupF32 (&status, header, "CRPIX1");102 crpix2 = psMetadataLookupF32 (&status, header, "CRPIX2");103 104 // test the CDELTi varient105 cdelt1 = psMetadataLookupF32 (&status, header, "CDELT1");106 if (status) {107 cdelt2 = psMetadataLookupF32 (&status, header, "CDELT2");108 109 // test the CROTAi varient:110 double rotate = psMetadataLookupF32 (&status, header, "CROTA2");111 if (status) {112 double Lambda = cdelt2 / cdelt1;113 pc1_1 = cos(rotate*RAD_DEG);114 pc1_2 = -sin(rotate*RAD_DEG) * Lambda;115 pc2_1 = sin(rotate*RAD_DEG) / Lambda;116 pc2_2 = cos(rotate*RAD_DEG);117 goto got_matrix;118 }119 120 // test the PC00i00j varient:121 pc1_1 = psMetadataLookupF32 (&status, header, "PC001001");122 if (status) {123 pc1_2 = psMetadataLookupF32 (&status, header, "PC001002");124 pc2_1 = psMetadataLookupF32 (&status, header, "PC002001");125 pc2_2 = psMetadataLookupF32 (&status, header, "PC002002");126 127 // XXX EAM : add Elixir polynomial terms here eventually128 goto got_matrix;129 }130 psLogMsg ("psastro", 2, "warning: missing rotation matrix?\n");131 return false;132 }133 134 // test the CDi_j varient135 pc1_1 = psMetadataLookupF32 (&status, header, "CD1_1");136 if (status) {137 pc1_2 = psMetadataLookupF32 (&status, header, "CD1_2");138 pc2_1 = psMetadataLookupF32 (&status, header, "CD2_1");139 pc2_2 = psMetadataLookupF32 (&status, header, "CD2_2");140 141 // renormalize to cdelt1, cdelt2, etc142 float scale = hypot (pc1_1, pc1_2);143 cdelt1 = cdelt2 = scale;144 pc1_1 /= scale;145 pc1_2 /= scale;146 pc2_1 /= scale;147 pc2_2 /= scale;148 goto got_matrix;149 }150 psLogMsg ("psastro", 2, "warning: missing rotation matrix?\n");151 return false;152 153 got_matrix:154 155 // XXX EAM : if fpa->toSky and fpa->toTPA are already defined, then the156 // toFPA must be modified to match the crval(i), scale(i) and crpix(i)157 // scale = scale(i)/scale(0) (i == chip #)158 // project crval1(0),crval2(0 using projection159 160 // XXX EAM : psPlaneTransformAlloc uses nTerm not nOrder (bug 581)161 // XXX EAM : I've fixed this in pslib eam_rel8_b2162 toFPA = psPlaneTransformAlloc (1, 1);163 164 toFPA->x->coeff[0][0] = crpix1;165 toFPA->x->coeff[1][0] = pc1_1;166 toFPA->x->coeff[0][1] = pc1_2;167 toFPA->x->mask[1][1] = 1;168 169 toFPA->y->coeff[0][0] = crpix2;170 toFPA->y->coeff[1][0] = pc2_1;171 toFPA->y->coeff[0][1] = pc2_2;172 toFPA->y->mask[1][1] = 1;173 174 // center of projection is (0,0) coordinate of TPA175 toSky = psProjectionAlloc (crval1*RAD_DEG, crval2*RAD_DEG, cdelt1*RAD_DEG, cdelt2*RAD_DEG, type);176 177 *toFPAOut = toFPA;178 *toSkyOut = toSky;179 180 return true;181 }182 183 184 // convert toFPA / toSky components to traditional WCS185 bool pmAstromWriteWCS (psPlaneTransform *toFPA, psProjection *toSky, psMetadata *header) {186 187 switch (toSky->type) {188 case PS_PROJ_SIN:189 psMetadataAdd (header, PS_LIST_TAIL, "CTYPE1", PS_DATA_STRING | PS_META_REPLACE, "", "RA---SIN");190 psMetadataAdd (header, PS_LIST_TAIL, "CTYPE2", PS_DATA_STRING | PS_META_REPLACE, "", "DEC--SIN");191 break;192 case PS_PROJ_TAN:193 psMetadataAdd (header, PS_LIST_TAIL, "CTYPE1", PS_DATA_STRING | PS_META_REPLACE, "", "RA---TAN");194 psMetadataAdd (header, PS_LIST_TAIL, "CTYPE2", PS_DATA_STRING | PS_META_REPLACE, "", "DEC--TAN");195 break;196 case PS_PROJ_AIT:197 psMetadataAdd (header, PS_LIST_TAIL, "CTYPE1", PS_DATA_STRING | PS_META_REPLACE, "", "RA---AIT");198 psMetadataAdd (header, PS_LIST_TAIL, "CTYPE2", PS_DATA_STRING | PS_META_REPLACE, "", "DEC--AIT");199 break;200 case PS_PROJ_PAR:201 psMetadataAdd (header, PS_LIST_TAIL, "CTYPE1", PS_DATA_STRING | PS_META_REPLACE, "", "RA---PAR");202 psMetadataAdd (header, PS_LIST_TAIL, "CTYPE2", PS_DATA_STRING | PS_META_REPLACE, "", "DEC--PAR");203 break;204 default:205 psLogMsg ("psastro", 2, "warning: unknown projection type %s\n", toSky->type);206 return false;207 }208 209 psMetadataAdd (header, PS_LIST_TAIL, "CRVAL1", PS_DATA_F32 | PS_META_REPLACE, "", toSky->R*DEG_RAD);210 psMetadataAdd (header, PS_LIST_TAIL, "CRVAL2", PS_DATA_F32 | PS_META_REPLACE, "", toSky->D*DEG_RAD);211 psMetadataAdd (header, PS_LIST_TAIL, "CRPIX1", PS_DATA_F32 | PS_META_REPLACE, "", toFPA->x->coeff[0][0]);212 psMetadataAdd (header, PS_LIST_TAIL, "CRPIX2", PS_DATA_F32 | PS_META_REPLACE, "", toFPA->y->coeff[0][0]);213 psMetadataAdd (header, PS_LIST_TAIL, "CDELT1", PS_DATA_F32 | PS_META_REPLACE, "", toSky->Xs*DEG_RAD);214 psMetadataAdd (header, PS_LIST_TAIL, "CDELT2", PS_DATA_F32 | PS_META_REPLACE, "", toSky->Ys*DEG_RAD);215 216 psMetadataAdd (header, PS_LIST_TAIL, "PC001001", PS_DATA_F32 | PS_META_REPLACE, "", toFPA->x->coeff[1][0]);217 psMetadataAdd (header, PS_LIST_TAIL, "PC001002", PS_DATA_F32 | PS_META_REPLACE, "", toFPA->x->coeff[0][1]);218 psMetadataAdd (header, PS_LIST_TAIL, "PC002001", PS_DATA_F32 | PS_META_REPLACE, "", toFPA->y->coeff[1][0]);219 psMetadataAdd (header, PS_LIST_TAIL, "PC002002", PS_DATA_F32 | PS_META_REPLACE, "", toFPA->y->coeff[0][1]);220 221 // alternative representations use222 // CD1_1 = PC001001*CDELT1, etc223 // make these representations optional224 225 return true;226 } -
trunk/psastro/src/psastroUtils.c
r5575 r6176 32 32 } 33 33 34 // measure per-chip astrometry terms 35 bool psastroChipAstrom (pmFPA *fpa, psMetadata *config) { 36 37 bool status; 38 psArray *match; 39 pmAstromStats stats; 34 bool psastroNormFPA (pmFPA *fpa, psMetadata *config) { 40 35 41 36 // save the raw astrometry for later reference … … 45 40 for (int i = 0; i < fpa->chips->n; i++) { 46 41 pmChip *chip = fpa->chips->data[i]; 47 48 // cells->n > 1 is not yet well-defined 49 if (chip->cells->n > 1) { 50 psLogMsg ("pmSourcesReadCMP", 3, "undefined behavior for nCells > 1"); 51 return false; 52 } 53 54 for (int j = 0; j < chip->cells->n; j++) { 55 pmCell *cell = chip->cells->data[j]; 56 57 // readouts->n > 1 is not yet well-defined 58 if (cell->readouts->n > 1) { 59 psLogMsg ("pmSourcesReadCMP", 3, "undefined behavior for nReadouts > 1"); 60 return false; 61 } 62 63 // load the corresponding reference data (DVO command) 64 psArray *refstars = psastroLoadReference (cell->header, config); 65 66 for (int k = 0; k < cell->readouts->n; k++) { 67 68 pmReadout *readout = cell->readouts->data[k]; 69 70 // pull out the SUBSET rawstars (a view) 71 psArray *rawstars = psMetadataLookupPtr (&status, readout->analysis, "STARS.SUBSET"); 72 73 // project the rawstars to the current best guess astrometry 74 psastroProjectRawstars (rawstars, readout); 75 76 // use the header & config info to project refstars onto the focal plane 77 psastroProjectRefstars (refstars, readout); 78 79 // testWriteRaw ("ref.inp", refstars); 80 // testWriteRaw ("raw.inp", rawstars); 81 82 // fprintf (stderr, "rawstars:\n"); 83 // psastroDumpStars (rawstars); 84 // fprintf (stderr, "refstars:\n"); 85 // psastroDumpStars (refstars); 86 87 // find initial offset / rotation 88 stats = pmAstromGridMatch (rawstars, refstars, config); 89 90 // adjust the chip.toFPA terms only 91 pmAstromGridApply (chip->toFPA, stats); 92 chip->fromFPA = p_psPlaneTransformLinearInvert (chip->toFPA); 93 94 // use fit result to re-project rawstars 95 psastroProjectRawstars (rawstars, readout); 96 psastroProjectRefstars (refstars, readout); 97 98 // testWriteRaw ("ref.dat", refstars); 99 // testWriteRaw ("raw.dat", rawstars); 100 101 // use small radius to match stars 102 match = pmAstromRadiusMatch (rawstars, refstars, config); 103 104 // improved fit for astrometric terms 105 pmAstromMatchFit (chip->toFPA, rawstars, refstars, match, config); 106 chip->fromFPA = p_psPlaneTransformLinearInvert (chip->toFPA); 107 } 108 pmAstromWriteWCS (chip->toFPA, fpa->toSky, cell->header); 109 } 42 psastroChipAstrom (chip, config); 110 43 } 111 44
Note:
See TracChangeset
for help on using the changeset viewer.
