IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 9732


Ignore:
Timestamp:
Oct 24, 2006, 1:58:13 PM (20 years ago)
Author:
Paul Price
Message:

Changed definition of psArrayAlloc.

Location:
trunk/psastro/src
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/psastro/src/psastroChooseRefstars.c

    r9627 r9732  
    1111    psMetadata *recipe  = psMetadataLookupPtr (NULL, config->recipes, "PSASTRO");
    1212    if (!recipe) {
    13         psError(PSASTRO_ERR_CONFIG, true, "Can't find PSASTRO recipe!\n");
    14         return false;
     13        psError(PSASTRO_ERR_CONFIG, true, "Can't find PSASTRO recipe!\n");
     14        return false;
    1515    }
    1616
    … …  
    1818    pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT");
    1919    if (!input) {
    20         psError(PSASTRO_ERR_CONFIG, true, "Can't find input data!\n");
    21         return false;
     20        psError(PSASTRO_ERR_CONFIG, true, "Can't find input data!\n");
     21        return false;
    2222    }
    2323
    24     float fieldExtra = psMetadataLookupS32 (&status, recipe, "PSASTRO.FIELD.EXTRA"); 
     24    float fieldExtra = psMetadataLookupS32 (&status, recipe, "PSASTRO.FIELD.EXTRA");
    2525    if (!status) fieldExtra = 0.0;
    2626
    … …  
    3232        psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
    3333        if (!chip->process || !chip->file_exists) { continue; }
    34        
    35         while ((cell = pmFPAviewNextCell (view, fpa, 1)) != NULL) {
     34
     35        while ((cell = pmFPAviewNextCell (view, fpa, 1)) != NULL) {
    3636            psTrace ("psastro", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
    3737            if (!cell->process || !cell->file_exists) { continue; }
    3838
    39             // process each of the readouts
    40             // XXX there can only be one readout per chip in astrometry, right?
    41             while ((readout = pmFPAviewNextReadout (view, fpa, 1)) != NULL) {
    42                 if (! readout->data_exists) { continue; }
     39            // process each of the readouts
     40            // XXX there can only be one readout per chip in astrometry, right?
     41            while ((readout = pmFPAviewNextReadout (view, fpa, 1)) != NULL) {
     42                if (! readout->data_exists) { continue; }
    4343
    44                 // read WCS data from the corresponding header
    45                 pmHDU *hdu = pmFPAviewThisHDU (view, fpa);
     44                // read WCS data from the corresponding header
     45                pmHDU *hdu = pmFPAviewThisHDU (view, fpa);
    4646
    47                 int Nx = psMetadataLookupS32 (&status, hdu->header, "NAXIS1");
    48                 int Ny = psMetadataLookupS32 (&status, hdu->header, "NAXIS2");
     47                int Nx = psMetadataLookupS32 (&status, hdu->header, "NAXIS1");
     48                int Ny = psMetadataLookupS32 (&status, hdu->header, "NAXIS2");
    4949
    50                 float minX = -fieldExtra*Nx;
    51                 float maxX = (1+fieldExtra)*Nx;
    52                 float minY = -fieldExtra*Ny;
    53                 float maxY = (1+fieldExtra)*Ny;
     50                float minX = -fieldExtra*Nx;
     51                float maxX = (1+fieldExtra)*Nx;
     52                float minY = -fieldExtra*Ny;
     53                float maxY = (1+fieldExtra)*Ny;
    5454
    55                 // the refstars is a subset within range of this chip
    56                 psArray *refstars = psArrayAlloc (100);
     55                // the refstars is a subset within range of this chip
     56                psArray *refstars = psArrayAllocEmpty (100);
    5757
    58                 // select the reference objects within range of this readout
    59                 // project the reference objects to this chip
    60                 for (int i = 0; i < refs->n; i++) {
    61                     pmAstromObj *ref = pmAstromObjCopy(refs->data[i]);
    62        
    63                     // XXX why is this still a private function?
    64                     p_psProject (ref->TP, ref->sky, fpa->projection);
    65                     psPlaneDistortApply (ref->FP, fpa->fromTangentPlane, ref->TP, 0.0, 0.0);
    66                     psPlaneTransformApply (ref->chip, chip->fromFPA, ref->FP);
     58                // select the reference objects within range of this readout
     59                // project the reference objects to this chip
     60                for (int i = 0; i < refs->n; i++) {
     61                    pmAstromObj *ref = pmAstromObjCopy(refs->data[i]);
    6762
    68                     // limit the X,Y range of the refs to the selected chip
    69                     if (ref->chip->x < minX) goto skip;
    70                     if (ref->chip->x > maxX) goto skip;
    71                     if (ref->chip->y < minY) goto skip;
    72                     if (ref->chip->y > maxY) goto skip;
     63                    // XXX why is this still a private function?
     64                    p_psProject (ref->TP, ref->sky, fpa->projection);
     65                    psPlaneDistortApply (ref->FP, fpa->fromTangentPlane, ref->TP, 0.0, 0.0);
     66                    psPlaneTransformApply (ref->chip, chip->fromFPA, ref->FP);
    7367
    74                     psArrayAdd (refstars, 100, ref);
    75                 skip:
    76                     psFree (ref);
    77                 }
    78                 psTrace ("psastro", 4, "Added %ld refstars\n", refstars->n);
     68                    // limit the X,Y range of the refs to the selected chip
     69                    if (ref->chip->x < minX) goto skip;
     70                    if (ref->chip->x > maxX) goto skip;
     71                    if (ref->chip->y < minY) goto skip;
     72                    if (ref->chip->y > maxY) goto skip;
    7973
    80                 psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.REFSTARS", PS_DATA_ARRAY, "astrometry matches", refstars);
    81                 psFree (refstars);
     74                    psArrayAdd (refstars, 100, ref);
     75                skip:
     76                    psFree (ref);
     77                }
     78                psTrace ("psastro", 4, "Added %ld refstars\n", refstars->n);
    8279
    83                 psastroRefstarSubset (readout);
    84             }
    85         }
     80                psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.REFSTARS", PS_DATA_ARRAY, "astrometry matches", refstars);
     81                psFree (refstars);
     82
     83                psastroRefstarSubset (readout);
     84            }
     85        }
    8686    }
    8787

  • trunk/psastro/src/psastroConvert.c

    r9574 r9732  
    88    pmReadout *readout;
    99    pmFPAview *view = pmFPAviewAlloc (0);
    10  
     10
    1111    while ((chip = pmFPAviewNextChip (view, fpa, 1)) != NULL) {
    1212        psTrace ("psastro", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
    1313        if (!chip->process || !chip->file_exists) { continue; }
    1414
    15         while ((cell = pmFPAviewNextCell (view, fpa, 1)) != NULL) {
     15        while ((cell = pmFPAviewNextCell (view, fpa, 1)) != NULL) {
    1616            psTrace ("psastro", 4, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
    1717            if (!cell->process || !cell->file_exists) { continue; }
    1818
    19             // process each of the readouts
    20             while ((readout = pmFPAviewNextReadout (view, fpa, 1)) != NULL) {
    21                 if (! readout->data_exists) { continue; }
     19            // process each of the readouts
     20            while ((readout = pmFPAviewNextReadout (view, fpa, 1)) != NULL) {
     21                if (! readout->data_exists) { continue; }
    2222
    23                 psastroConvertReadout (readout, recipe);
    24             }
    25         }
     23                psastroConvertReadout (readout, recipe);
     24            }
     25        }
    2626    }
    2727    psFree (view);
    … …  
    4040    psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.RAWSTARS", PS_DATA_ARRAY, "astrometry objects", rawstars);
    4141    psFree (rawstars);
    42  
     42
    4343    return true;
    4444}
    … …  
    4747psArray *pmSourceToAstromObj (psArray *sources) {
    4848
    49     psArray *objects = psArrayAlloc (sources->n);
    50    
     49    psArray *objects = psArrayAllocEmpty (sources->n);
     50
    5151    for (int i = 0; i < sources->n; i++) {
    52         pmSource *source = sources->data[i];
    53        
    54         // only accept the PSF sources?
    55         if (source->type != PM_SOURCE_TYPE_STAR) continue;
     52        pmSource *source = sources->data[i];
    5653
    57         pmModel *model = source->modelPSF;
    58         if (model == NULL) continue;
     54        // only accept the PSF sources?
     55        if (source->type != PM_SOURCE_TYPE_STAR) continue;
    5956
    60         psF32 *PAR = model->params->data.F32;
     57        pmModel *model = source->modelPSF;
     58        if (model == NULL) continue;
    6159
    62         pmAstromObj *obj = pmAstromObjAlloc ();
     60        psF32 *PAR = model->params->data.F32;
    6361
    64         // is the source magnitude calibrated in any sense?
    65         obj->pix->x = PAR[2];
    66         obj->pix->y = PAR[3];
    67         obj->Mag = source->psfMag;
     62        pmAstromObj *obj = pmAstromObjAlloc ();
    6863
    69         // XXX do we have the information giving the readout and cell offset?
    70         // for the moment, assume chip == cell == readout
    71         *obj->cell = *obj->pix;
    72         *obj->chip = *obj->cell;
     64        // is the source magnitude calibrated in any sense?
     65        obj->pix->x = PAR[2];
     66        obj->pix->y = PAR[3];
     67        obj->Mag = source->psfMag;
    7368
    74         psArrayAdd (objects, 100, obj);
    75         psFree (obj);
     69        // XXX do we have the information giving the readout and cell offset?
     70        // for the moment, assume chip == cell == readout
     71        *obj->cell = *obj->pix;
     72        *obj->chip = *obj->cell;
     73
     74        psArrayAdd (objects, 100, obj);
     75        psFree (obj);
    7676    }
    7777    return objects;

  • trunk/psastro/src/psastroLoadRefstars.c

    r9645 r9732  
    2424    char *CATDIR = psMetadataLookupStr(NULL, recipe, "DVO.CATDIR");
    2525    if (CATDIR == NULL) {
    26         psLogMsg ("psastro", 2, "warning: missing DVO.CATDIR, using ~/.ptolemyrc definition\n");
     26        psLogMsg ("psastro", 2, "warning: missing DVO.CATDIR, using ~/.ptolemyrc definition\n");
    2727    }
    2828
    … …  
    3333    sprintf (tempFile, "/tmp/psastro.XXXXXX");
    3434    if ((fd = mkstemp (tempFile)) == -1) {
    35         psError(PSASTRO_ERR_REFSTARS, true, "error creating temp output file\n");
     35        psError(PSASTRO_ERR_REFSTARS, true, "error creating temp output file\n");
    3636        return NULL;
    3737    }
    … …  
    4242    // use fork to add timeout capability
    4343    if (ELIXIR_MODE) {
    44         psStringAppend (&catformat, "-D CATFORMAT elixir ");
     44        psStringAppend (&catformat, "-D CATFORMAT elixir ");
    4545    } else {
    46         psStringAppend (&catformat, "-D CATFORMAT panstarrs ");
    47     }   
     46        psStringAppend (&catformat, "-D CATFORMAT panstarrs ");
     47    }
    4848
    4949    if (CATDIR) {
    50         psStringAppend (&catformat, "-D CATDIR %s ", CATDIR);
    51     } 
     50        psStringAppend (&catformat, "-D CATDIR %s ", CATDIR);
     51    }
    5252
    5353    psStringAppend (&getstarLine, "getstar %s -D CATMODE mef -maglim %f -region %f %f %f %f -o %s", catformat, MAGmax, RAmin, DECmin, RAmax, DECmax, tempFile);
    … …  
    5757    status = system (getstarLine);
    5858    if (status) {
    59         psError(PSASTRO_ERR_REFSTARS, true, "error loading reference data\n");
     59        psError(PSASTRO_ERR_REFSTARS, true, "error loading reference data\n");
    6060        return NULL;
    6161    }
    … …  
    6969
    7070    if (ELIXIR_MODE) {
    71         psFitsMoveExtName (fits, "DVO_AVERAGE_ELIXIR");
     71        psFitsMoveExtName (fits, "DVO_AVERAGE_ELIXIR");
    7272    } else {
    73         psFitsMoveExtName (fits, "DVO_AVERAGE_PANSTARRS");
     73        psFitsMoveExtName (fits, "DVO_AVERAGE_PANSTARRS");
    7474    }
    7575
    … …  
    8484    // convert the Average table to the pmAstromObj entries
    8585    psTimerStart ("psastro");
    86     psArray *refs = psArrayAlloc (table->n);
     86    psArray *refs = psArrayAllocEmpty (table->n);
    8787    for (int i = 0; i < table->n; i++) {
    8888        pmAstromObj *ref = pmAstromObjAlloc ();
    … …  
    9191
    9292        // DVO tables are stored in degrees
    93         if (ELIXIR_MODE) {
    94             ref->sky->r   = RAD_DEG*psMetadataLookupF32 (&status, row, "RA");
    95             ref->sky->d   = RAD_DEG*psMetadataLookupF32 (&status, row, "DEC");
    96             ref->Mag      = 0.001*psMetadataLookupS32 (&status, row, "MAG");  // ELIXIR uses millimags
     93        if (ELIXIR_MODE) {
     94            ref->sky->r   = RAD_DEG*psMetadataLookupF32 (&status, row, "RA");
     95            ref->sky->d   = RAD_DEG*psMetadataLookupF32 (&status, row, "DEC");
     96            ref->Mag      = 0.001*psMetadataLookupS32 (&status, row, "MAG");  // ELIXIR uses millimags
    9797        } else {
    98             ref->sky->r   = RAD_DEG*psMetadataLookupF64 (&status, row, "RA");
    99             ref->sky->d   = RAD_DEG*psMetadataLookupF64 (&status, row, "DEC");
    100             ref->Mag      = psMetadataLookupF32 (&status, row, "MAG"); // PANSTARRS uses mags
    101         }
     98            ref->sky->r   = RAD_DEG*psMetadataLookupF64 (&status, row, "RA");
     99            ref->sky->d   = RAD_DEG*psMetadataLookupF64 (&status, row, "DEC");
     100            ref->Mag      = psMetadataLookupF32 (&status, row, "MAG"); // PANSTARRS uses mags
     101        }
    102102
    103103        psArrayAdd (refs, 100, ref);

  • trunk/psastro/src/psastroRefstarSubset.c

    r9649 r9732  
    3232    return NULL;
    3333  }
    34  
     34
    3535  // XXXX test
    3636  // psFree (rawfunc);
    … …  
    3838  // return true;
    3939
    40   // what is the offset between the two lines at the average magnitude? 
     40  // what is the offset between the two lines at the average magnitude?
    4141  double mRef = 0.5*(reffunc->mMin + reffunc->mMax);
    4242  double logRho = mRef * reffunc->slope + reffunc->offset;
    … …  
    4848  psLogMsg ("psastro", 4, "clipping stars fainter than %f\n", mRefMax);
    4949
    50   psArray *subset = psArrayAlloc (100);
     50  psArray *subset = psArrayAllocEmpty (100);
    5151  for (int i = 0; i < refstars->n; i++) {
    5252    pmAstromObj *ref = refstars->data[i];

  • trunk/psastro/src/psastroUtils.c

    r9373 r9732  
    1010
    1111    for (int i = 0; i < rawstars->n; i++) {
    12         pmAstromObj *raw = rawstars->data[i];
    13        
    14         psPlaneTransformApply (raw->FP, chip->toFPA, raw->chip);
    15         psPlaneDistortApply (raw->TP, fpa->toTangentPlane, raw->FP, 0.0, 0.0);
    16         p_psDeproject (raw->sky, raw->TP, fpa->projection);
     12        pmAstromObj *raw = rawstars->data[i];
     13
     14        psPlaneTransformApply (raw->FP, chip->toFPA, raw->chip);
     15        psPlaneDistortApply (raw->TP, fpa->toTangentPlane, raw->FP, 0.0, 0.0);
     16        p_psDeproject (raw->sky, raw->TP, fpa->projection);
    1717    }
    1818
    1919    for (int i = 0; i < refstars->n; i++) {
    20         pmAstromObj *ref = refstars->data[i];
    21         psPlaneTransformApply (ref->chip, chip->fromFPA, ref->FP);
     20        pmAstromObj *ref = refstars->data[i];
     21        psPlaneTransformApply (ref->chip, chip->fromFPA, ref->FP);
    2222    }
    2323
    … …  
    3535
    3636    for (int i = 0; i < fpa->chips->n; i++) {
    37         pmChip *chip = fpa->chips->data[i];
    38         for (int j = 0; j < chip->cells->n; j++) {
    39             pmCell *cell = chip->cells->data[j];
    40             for (int k = 0; k < cell->readouts->n; k++) {
    41                 pmReadout *readout = cell->readouts->data[k];
    42 
    43                 psArray *stars = psMetadataLookupPtr (&status, readout->analysis, "STARS.FULLSET");
    44                 stars = psArraySort (stars, pmAstromObjSortByMag);
    45 
    46                 int nSubset = PS_MIN (MAX_NSTARS, stars->n);
    47                 psArray *subset = psArrayAlloc (nSubset);
    48 
    49                 for (int i = 0; i < nSubset; i++) {
    50                     subset->data[i] = stars->data[i];
    51                 }
    52                 psMetadataAdd (readout->analysis, PS_LIST_TAIL, "STARS.SUBSET", PS_DATA_ARRAY, "stars from analysis", subset);
    53             }
    54         }
     37        pmChip *chip = fpa->chips->data[i];
     38        for (int j = 0; j < chip->cells->n; j++) {
     39            pmCell *cell = chip->cells->data[j];
     40            for (int k = 0; k < cell->readouts->n; k++) {
     41                pmReadout *readout = cell->readouts->data[k];
     42
     43                psArray *stars = psMetadataLookupPtr (&status, readout->analysis, "STARS.FULLSET");
     44                stars = psArraySort (stars, pmAstromObjSortByMag);
     45
     46                int nSubset = PS_MIN (MAX_NSTARS, stars->n);
     47                psArray *subset = psArrayAlloc (nSubset);
     48
     49                for (int i = 0; i < nSubset; i++) {
     50                    subset->data[i] = stars->data[i];
     51                }
     52                psMetadataAdd (readout->analysis, PS_LIST_TAIL, "STARS.SUBSET", PS_DATA_ARRAY, "stars from analysis", subset);
     53            }
     54        }
    5555    }
    5656    return true;
    … …  
    6464    // first pass: measure the per-chip solutions, modify the chip.toFPA terms
    6565    for (int i = 0; i < fpa->chips->n; i++) {
    66         pmChip *chip = fpa->chips->data[i];
    67         psastroChipAstrom (chip, config);
    68     }
    69 
    70     // second stage: re-normalize the chip terms, passing the 
     66        pmChip *chip = fpa->chips->data[i];
     67        psastroChipAstrom (chip, config);
     68    }
     69
     70    // second stage: re-normalize the chip terms, passing the
    7171    // average rotation and offset values to the fpa.toSky
    7272    if (RENORM) {
    7373
    74         // this code is needed for the mosastro stage, with multiple chip solutions
    75 
    76         double dX, dY, dT, dN;
    77         dX = dY = dT = dN = 0;
    78 
    79         psPlane origin, P1, P2;
    80         origin.x = 0;
    81         origin.y = 0;
    82 
    83         // calculate the average rotation and boresite offset relative to raw
    84         for (int i = 0; i < fpa->chips->n; i++) {
    85             pmChip *iChip = raw->chips->data[i];
    86             pmChip *oChip = fpa->chips->data[i];
    87 
    88             // offset of chip
    89             psCoordChipToFPA (&P1, &origin, iChip);
    90             psCoordChipToFPA (&P2, &origin, oChip);
    91             dX += (P2.x - P1.x);
    92             dY += (P2.y - P1.y);
    93 
    94             // get parity-independent rotations for old and new solutions
    95             double T1 = psPlaneTransformGetRotation (iChip->toFPA);
    96             double T2 = psPlaneTransformGetRotation (oChip->toFPA);
    97             dT += T2 - T1;
    98             dN ++;
    99         }
    100 
    101         dT /= dN;
    102         dX /= dN;
    103         dY /= dN;
    104 
    105         // R(T)
    106         double PC1_1 = fpa->toTPA->x->coeff[1][0][0][0];
    107         double PC1_2 = fpa->toTPA->x->coeff[0][1][0][0];
    108         double PC2_1 = fpa->toTPA->y->coeff[1][0][0][0];
    109         double PC2_2 = fpa->toTPA->y->coeff[0][1][0][0];
    110 
    111         // R(dT)
    112         double dPC1_1 = +cos (dT);
    113         double dPC1_2 = +sin (dT);
    114         double dPC2_1 = -sin (dT);
    115         double dPC2_2 = +cos (dT);
    116 
    117         // R'(T) = R(T) * R(dT)
    118         double pc1_1 = PC1_1*dPC1_1 + PC1_2*dPC2_1;
    119         double pc1_2 = PC1_1*dPC1_2 + PC1_2*dPC2_2;
    120         double pc2_1 = PC2_1*dPC1_1 + PC2_2*dPC2_1;
    121         double pc2_2 = PC2_1*dPC1_2 + PC2_2*dPC2_2;
    122 
    123         double det = 1.0 / (pc1_1*pc2_2 - pc1_2*pc2_1);
    124 
    125         // R'(-T)  (matrix inverse, not just rotation inverse -- keeps parity)
    126         double pi1_1 = +pc2_2 * det;
    127         double pi1_2 = -pc1_2 * det;
    128         double pi2_1 = -pc2_1 * det;
    129         double pi2_2 = +pc1_1 * det;
    130    
    131         // apply the new modifcations in rotation and boresite
    132         for (int i = 0; i < fpa->chips->n; i++) {
    133             pmChip *oChip = fpa->chips->data[i];
    134 
    135             // r(T)
    136             double pr1_1 = oChip->toFPA->x->coeff[1][0];
    137             double pr1_2 = oChip->toFPA->x->coeff[0][1];
    138             double pr2_1 = oChip->toFPA->y->coeff[1][0];
    139             double pr2_2 = oChip->toFPA->y->coeff[0][1];
    140 
    141             // ri'(T) = R(T) r(t)
    142             double ri1_1 = PC1_1*pr1_1 + PC1_2*pr2_1;
    143             double ri1_2 = PC1_1*pr1_2 + PC1_2*pr2_2;
    144             double ri2_1 = PC2_1*pr1_1 + PC2_2*pr2_1;
    145             double ri2_2 = PC2_1*pr1_2 + PC2_2*pr2_2;
    146 
    147             // r'(T) = R'(-T) ri'(T)
    148             oChip->toFPA->x->coeff[1][0] = pi1_1*ri1_1 + pi1_2*ri2_1;
    149             oChip->toFPA->x->coeff[0][1] = pi1_1*ri1_2 + pi1_2*ri2_2;
    150             oChip->toFPA->y->coeff[1][0] = pi2_1*ri1_1 + pi2_2*ri2_1;
    151             oChip->toFPA->y->coeff[0][1] = pi2_1*ri1_2 + pi2_2*ri2_2;
    152 
    153             double dx = PC1_1*oChip->toFPA->x->coeff[0][0] + PC1_2*oChip->toFPA->y->coeff[0][0] + dX;
    154             double dy = PC2_1*oChip->toFPA->x->coeff[0][0] + PC2_2*oChip->toFPA->y->coeff[0][0] + dY;
    155 
    156             oChip->toFPA->x->coeff[0][0] = pi1_1*dx + pi1_2*dy;
    157             oChip->toFPA->y->coeff[0][0] = pi2_1*dx + pi2_2*dy;
    158         }
    159 
    160         fpa->toTPA->x->coeff[0][0][0][0] -= dX;
    161         fpa->toTPA->y->coeff[0][0][0][0] -= dY;
    162 
    163         fpa->toTPA->x->coeff[1][0][0][0] = pc1_1;
    164         fpa->toTPA->x->coeff[0][1][0][0] = pc1_2;
    165         fpa->toTPA->y->coeff[1][0][0][0] = pc2_1;
    166         fpa->toTPA->y->coeff[0][1][0][0] = pc2_2;
    167     }
    168     return true; 
     74        // this code is needed for the mosastro stage, with multiple chip solutions
     75
     76        double dX, dY, dT, dN;
     77        dX = dY = dT = dN = 0;
     78
     79        psPlane origin, P1, P2;
     80        origin.x = 0;
     81        origin.y = 0;
     82
     83        // calculate the average rotation and boresite offset relative to raw
     84        for (int i = 0; i < fpa->chips->n; i++) {
     85            pmChip *iChip = raw->chips->data[i];
     86            pmChip *oChip = fpa->chips->data[i];
     87
     88            // offset of chip
     89            psCoordChipToFPA (&P1, &origin, iChip);
     90            psCoordChipToFPA (&P2, &origin, oChip);
     91            dX += (P2.x - P1.x);
     92            dY += (P2.y - P1.y);
     93
     94            // get parity-independent rotations for old and new solutions
     95            double T1 = psPlaneTransformGetRotation (iChip->toFPA);
     96            double T2 = psPlaneTransformGetRotation (oChip->toFPA);
     97            dT += T2 - T1;
     98            dN ++;
     99        }
     100
     101        dT /= dN;
     102        dX /= dN;
     103        dY /= dN;
     104
     105        // R(T)
     106        double PC1_1 = fpa->toTPA->x->coeff[1][0][0][0];
     107        double PC1_2 = fpa->toTPA->x->coeff[0][1][0][0];
     108        double PC2_1 = fpa->toTPA->y->coeff[1][0][0][0];
     109        double PC2_2 = fpa->toTPA->y->coeff[0][1][0][0];
     110
     111        // R(dT)
     112        double dPC1_1 = +cos (dT);
     113        double dPC1_2 = +sin (dT);
     114        double dPC2_1 = -sin (dT);
     115        double dPC2_2 = +cos (dT);
     116
     117        // R'(T) = R(T) * R(dT)
     118        double pc1_1 = PC1_1*dPC1_1 + PC1_2*dPC2_1;
     119        double pc1_2 = PC1_1*dPC1_2 + PC1_2*dPC2_2;
     120        double pc2_1 = PC2_1*dPC1_1 + PC2_2*dPC2_1;
     121        double pc2_2 = PC2_1*dPC1_2 + PC2_2*dPC2_2;
     122
     123        double det = 1.0 / (pc1_1*pc2_2 - pc1_2*pc2_1);
     124
     125        // R'(-T)  (matrix inverse, not just rotation inverse -- keeps parity)
     126        double pi1_1 = +pc2_2 * det;
     127        double pi1_2 = -pc1_2 * det;
     128        double pi2_1 = -pc2_1 * det;
     129        double pi2_2 = +pc1_1 * det;
     130
     131        // apply the new modifcations in rotation and boresite
     132        for (int i = 0; i < fpa->chips->n; i++) {
     133            pmChip *oChip = fpa->chips->data[i];
     134
     135            // r(T)
     136            double pr1_1 = oChip->toFPA->x->coeff[1][0];
     137            double pr1_2 = oChip->toFPA->x->coeff[0][1];
     138            double pr2_1 = oChip->toFPA->y->coeff[1][0];
     139            double pr2_2 = oChip->toFPA->y->coeff[0][1];
     140
     141            // ri'(T) = R(T) r(t)
     142            double ri1_1 = PC1_1*pr1_1 + PC1_2*pr2_1;
     143            double ri1_2 = PC1_1*pr1_2 + PC1_2*pr2_2;
     144            double ri2_1 = PC2_1*pr1_1 + PC2_2*pr2_1;
     145            double ri2_2 = PC2_1*pr1_2 + PC2_2*pr2_2;
     146
     147            // r'(T) = R'(-T) ri'(T)
     148            oChip->toFPA->x->coeff[1][0] = pi1_1*ri1_1 + pi1_2*ri2_1;
     149            oChip->toFPA->x->coeff[0][1] = pi1_1*ri1_2 + pi1_2*ri2_2;
     150            oChip->toFPA->y->coeff[1][0] = pi2_1*ri1_1 + pi2_2*ri2_1;
     151            oChip->toFPA->y->coeff[0][1] = pi2_1*ri1_2 + pi2_2*ri2_2;
     152
     153            double dx = PC1_1*oChip->toFPA->x->coeff[0][0] + PC1_2*oChip->toFPA->y->coeff[0][0] + dX;
     154            double dy = PC2_1*oChip->toFPA->x->coeff[0][0] + PC2_2*oChip->toFPA->y->coeff[0][0] + dY;
     155
     156            oChip->toFPA->x->coeff[0][0] = pi1_1*dx + pi1_2*dy;
     157            oChip->toFPA->y->coeff[0][0] = pi2_1*dx + pi2_2*dy;
     158        }
     159
     160        fpa->toTPA->x->coeff[0][0][0][0] -= dX;
     161        fpa->toTPA->y->coeff[0][0][0][0] -= dY;
     162
     163        fpa->toTPA->x->coeff[1][0][0][0] = pc1_1;
     164        fpa->toTPA->x->coeff[0][1][0][0] = pc1_2;
     165        fpa->toTPA->y->coeff[1][0][0][0] = pc2_1;
     166        fpa->toTPA->y->coeff[0][1][0][0] = pc2_2;
     167    }
     168    return true;
    169169}
    170170
    … …  
    174174
    175175    for (int i = 0; i < input->nX; i++) {
    176         for (int j = 0; j < input->nY; j++) {
    177             output->mask[i][j]     = input->mask[i][j];
    178             output->coeff[i][j]    = input->coeff[i][j];
    179             output->coeffErr[i][j] = input->coeffErr[i][j];
    180         }
     176        for (int j = 0; j < input->nY; j++) {
     177            output->mask[i][j]     = input->mask[i][j];
     178            output->coeff[i][j]    = input->coeff[i][j];
     179            output->coeffErr[i][j] = input->coeffErr[i][j];
     180        }
    181181    }
    182182    return (output);
    … …  
    188188
    189189    for (int i = 0; i < input->nX; i++) {
    190         for (int j = 0; j < input->nY; j++) {
    191             for (int k = 0; k < input->nZ; k++) {
    192                 for (int m = 0; m < input->nT; m++) {
    193                     output->mask[i][j][k][m]     = input->mask[i][j][k][m];
    194                     output->coeff[i][j][k][m]    = input->coeff[i][j][k][m];
    195                     output->coeffErr[i][j][k][m] = input->coeffErr[i][j][k][m];
    196                 }
    197             }
    198         }
     190        for (int j = 0; j < input->nY; j++) {
     191            for (int k = 0; k < input->nZ; k++) {
     192                for (int m = 0; m < input->nT; m++) {
     193                    output->mask[i][j][k][m]     = input->mask[i][j][k][m];
     194                    output->coeff[i][j][k][m]    = input->coeff[i][j][k][m];
     195                    output->coeffErr[i][j][k][m] = input->coeffErr[i][j][k][m];
     196                }
     197            }
     198        }
    199199    }
    200200    return (output);
    … …  
    206206
    207207    for (int i = 0; i < input->x->nX; i++) {
    208         for (int j = 0; j < input->x->nY; j++) {
    209             for (int k = 0; k < input->x->nZ; k++) {
    210                 for (int m = 0; m < input->x->nT; m++) {
    211                     // x-terms
    212                     output->x->mask[i][j][k][m]     = input->x->mask[i][j][k][m];
    213                     output->x->coeff[i][j][k][m]    = input->x->coeff[i][j][k][m];
    214                     output->x->coeffErr[i][j][k][m] = input->x->coeffErr[i][j][k][m];
    215                     // y-terms
    216                     output->y->mask[i][j][k][m]     = input->y->mask[i][j][k][m];
    217                     output->y->coeff[i][j][k][m]    = input->y->coeff[i][j][k][m];
    218                     output->y->coeffErr[i][j][k][m] = input->y->coeffErr[i][j][k][m];
    219                 }
    220             }
    221         }
     208        for (int j = 0; j < input->x->nY; j++) {
     209            for (int k = 0; k < input->x->nZ; k++) {
     210                for (int m = 0; m < input->x->nT; m++) {
     211                    // x-terms
     212                    output->x->mask[i][j][k][m]     = input->x->mask[i][j][k][m];
     213                    output->x->coeff[i][j][k][m]    = input->x->coeff[i][j][k][m];
     214                    output->x->coeffErr[i][j][k][m] = input->x->coeffErr[i][j][k][m];
     215                    // y-terms
     216                    output->y->mask[i][j][k][m]     = input->y->mask[i][j][k][m];
     217                    output->y->coeff[i][j][k][m]    = input->y->coeff[i][j][k][m];
     218                    output->y->coeffErr[i][j][k][m] = input->y->coeffErr[i][j][k][m];
     219                }
     220            }
     221        }
    222222    }
    223223    return (output);
    … …  
    229229
    230230    for (int i = 0; i < input->x->nX; i++) {
    231         for (int j = 0; j < input->x->nY; j++) {
    232             // x-terms
    233             output->x->mask[i][j]     = input->x->mask[i][j];
    234             output->x->coeff[i][j]    = input->x->coeff[i][j];
    235             output->x->coeffErr[i][j] = input->x->coeffErr[i][j];
    236             // y-terms
    237             output->y->mask[i][j]     = input->y->mask[i][j];
    238             output->y->coeff[i][j]    = input->y->coeff[i][j];
    239             output->y->coeffErr[i][j] = input->y->coeffErr[i][j];
    240         }
     231        for (int j = 0; j < input->x->nY; j++) {
     232            // x-terms
     233            output->x->mask[i][j]     = input->x->mask[i][j];
     234            output->x->coeff[i][j]    = input->x->coeff[i][j];
     235            output->x->coeffErr[i][j] = input->x->coeffErr[i][j];
     236            // y-terms
     237            output->y->mask[i][j]     = input->y->mask[i][j];
     238            output->y->coeff[i][j]    = input->y->coeff[i][j];
     239            output->y->coeffErr[i][j] = input->y->coeffErr[i][j];
     240        }
    241241    }
    242242    return (output);
    … …  
    287287    if (map->y->nX < 1) return 0;
    288288    if (map->y->nY < 1) return 0;
    289    
     289
    290290    double pc1_1 = map->x->coeff[1][0];
    291291    double pc1_2 = map->x->coeff[0][1];
    … …  
    299299    double t1 = -atan2 (px*pc1_2, px*pc1_1);
    300300    double t2 = +atan2 (py*pc2_1, py*pc2_2);
    301    
     301
    302302    // careful near -pi,+pi boundary...
    303303    if (t1 - t2 > M_PI/2) t2 += 2*M_PI;
    … …  
    307307    while (theta < M_PI) theta += 2*M_PI;
    308308    while (theta > M_PI) theta -= 2*M_PI;
    309    
     309
    310310    return (theta);
    311311}
Note: See TracChangeset for help on using the changeset viewer.