IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 5560


Ignore:
Timestamp:
Nov 21, 2005, 11:03:53 AM (20 years ago)
Author:
eugene
Message:

various cleanups and fixes

Location:
trunk/psastro
Files:
10 edited

Legend:

Unmodified
Added
Removed

  • trunk/psastro/Makefile

    r5360 r5560  
    1 default: psphot
     1default: psastro
    22help:
    3         @echo "USAGE: make psphot"
     3        @echo "USAGE: make psastro"
    44
    55CC      =       gcc
    … …  
    1212IPSLIB  :=      $(shell pslib-config --cflags)
    1313
    14 INCS    =       $(IPSLIB)
    15 # LIBS  =       -lpsmodule $(LPSLIB)
    16 LIBS    =       $(LPSLIB)
     14LPSMOD  :=      $(shell psmodule-config --libs)
     15IPSMOD  :=      $(shell psmodule-config --cflags)
     16
     17INCS    =       $(IPSLIB) $(IPSMOD)
     18LIBS    =       $(LPSLIB) $(LPSMOD)
    1719CFLAGS  =       $(INCS) -std=c99 -Wall -Werror -g
    1820# CFLAGS        =       $(INCS) -std=c99 -g
    1921LFLAGS  =       $(LIBS)
    2022
    21 PSPHOT = \
    22 $(SRC)/psphot.$(ARCH).o            \
    23 $(SRC)/psphotArguments.$(ARCH).o   \
    24 $(SRC)/psphotSetup.$(ARCH).o       \
    25 $(SRC)/psphotImageStats.$(ARCH).o  \
    26 $(SRC)/psphotSourceStats.$(ARCH).o \
    27 $(SRC)/psphotChoosePSF.$(ARCH).o   \
    28 $(SRC)/psphotApplyPSF.$(ARCH).o    \
    29 $(SRC)/psphotFitGalaxies.$(ARCH).o \
    30 $(SRC)/psphotOutput.$(ARCH).o      \
    31 $(SRC)/psphotMarkPSF.$(ARCH).o     \
    32 $(SRC)/psphotSubtractPSF.$(ARCH).o \
    33 $(SRC)/psphotSortBySN.$(ARCH).o    \
    34 $(SRC)/psphotDefinePixels.$(ARCH).o \
    35 $(SRC)/psphotMagnitudes.$(ARCH).o  \
     23PSASTRO = \
     24$(SRC)/psastro.$(ARCH).o           \
     25$(SRC)/psastroIO.$(ARCH).o         \
     26$(SRC)/psastroBuildFPA.$(ARCH).o   \
     27$(SRC)/psastroArguments.$(ARCH).o  \
     28$(SRC)/psastroUtils.$(ARCH).o      \
    3629$(SRC)/psLibUtils.$(ARCH).o        \
    37 $(SRC)/psLine.$(ARCH).o            \
    38 $(SRC)/psMinimize.$(ARCH).o        \
    39 $(SRC)/psPolynomials.$(ARCH).o     \
    40 $(SRC)/psEllipse.$(ARCH).o         \
    41 $(SRC)/psModulesUtils.$(ARCH).o    \
    42 $(SRC)/pmPeaksSigmaLimit.$(ARCH).o \
    43 $(SRC)/pmPSF.$(ARCH).o             \
    44 $(SRC)/pmPSFtry.$(ARCH).o          \
    45 $(SRC)/psImageData.$(ARCH).o       \
    46 $(SRC)/pmModelGroup.$(ARCH).o       \
    47 $(SRC)/pmObjects_EAM.$(ARCH).o
     30$(SRC)/pmAstrom.$(ARCH).o          \
     31$(SRC)/pmAstromGrid.$(ARCH).o
    4832
    49 MODELS = \
    50 $(SRC)/models/pmModel_GAUSS.c      \
    51 $(SRC)/models/pmModel_PGAUSS.c     \
    52 $(SRC)/models/pmModel_QGAUSS.c     \
    53 $(SRC)/models/pmModel_SGAUSS.c     \
    54 $(SRC)/models/pmModel_TGAUSS.c     
     33PSASTRO-MKTEST = \
     34$(SRC)/psastro-mktest.$(ARCH).o    \
     35$(SRC)/psastroArguments.$(ARCH).o  \
     36$(SRC)/psLibUtils.$(ARCH).o        \
     37$(SRC)/pmAstrom.$(ARCH).o
    5538
    56 $(SRC)/pmModelGroup.$(ARCH).o: $(MODELS)
     39# $(SRC)/psModUtils.$(ARCH).o      \
    5740
    58 # deprecated
     41psastro: $(BIN)/psastro.$(ARCH)
     42$(BIN)/psastro.$(ARCH) : $(PSASTRO)
     43$(PSASTRO): $(SRC)/psastro.h $(SRC)/pmAstrom.h
    5944
    60 FITSOURCE = \
    61 $(SRC)/fitsource.$(ARCH).o \
    62 $(SRC)/onesource.$(ARCH).o \
    63 $(SRC)/fs_args.$(ARCH).o \
    64 $(SRC)/psphot-utils.$(ARCH).o \
    65 $(SRC)/psPolynomials.$(ARCH).o \
    66 $(SRC)/psUtils.$(ARCH).o \
    67 $(SRC)/setup.$(ARCH).o \
    68 $(SRC)/LocalSky.$(ARCH).o
     45psastro-mktest: $(BIN)/psastro-mktest.$(ARCH)
     46$(BIN)/psastro-mktest.$(ARCH) : $(PSASTRO-MKTEST)
     47$(PSASTRO-MKTEST): $(SRC)/psastro.h $(SRC)/pmAstrom.h
    6948
    70 psphot: $(BIN)/psphot.$(ARCH)
    71 $(BIN)/psphot.$(ARCH) : $(PSPHOT)
    72 $(PSPHOT): $(SRC)/psphot.h
    73 
    74 fitsource: $(BIN)/fitsource.$(ARCH)
    75 $(BIN)/fitsource.$(ARCH) : $(FITSOURCE)
    76 $(FITSOURCE): $(SRC)/psphot.h
    77 
    78 INSTALL = psphot
     49INSTALL = psastro psastro-mktest
    7950
    8051# dependancy rules for binary code #########################

  • trunk/psastro/src/pmAstrom.c

    r5510 r5560  
    11# include "pmAstrom.h"
    22
     3// sort by mag (descending)
     4int pmAstromObjSortByFPX (const void **a, const void **b)
     5{
     6    pmAstromObj *A = *(pmAstromObj **)a;
     7    pmAstromObj *B = *(pmAstromObj **)b;
     8
     9    psF32 diff = A->FP.x - B->FP.x;
     10    if (diff > FLT_EPSILON) return (+1);
     11    if (diff < FLT_EPSILON) return (-1);
     12    return (0);
     13}
     14
     15psPlane *psCoordReadoutToCell (psPlane *cellpix, psPlane *readpix, pmReadout *readout) {
     16
     17    if (cellpix == NULL) {
     18        cellpix = psPlaneAlloc ();
     19    }
     20
     21    cellpix->x = readpix->x*readout->colBins + readout->col0;
     22    cellpix->y = readpix->y*readout->rowBins + readout->row0;
     23
     24    return (cellpix);
     25}
     26
     27psPlane *psCoordCellToReadout (psPlane *readpix, psPlane *cellpix, pmReadout *readout) {
     28
     29    if (readpix == NULL) {
     30        readpix = psPlaneAlloc ();
     31    }
     32
     33    readpix->x = (cellpix->x - readout->col0) / readout->colBins;
     34    readpix->y = (cellpix->y - readout->row0) / readout->rowBins;
     35
     36    return (readpix);
     37}
     38
     39psPlane* psCoordChipToCell_EAM(psPlane* cellCoord,
     40                           const psPlane* chipCoord,
     41                           const pmCell* cell)
     42{
     43    PS_ASSERT_PTR_NON_NULL(chipCoord, NULL);
     44    PS_ASSERT_PTR_NON_NULL(cell, NULL);
     45    PS_ASSERT_PTR_NON_NULL(cell->parent, NULL);
     46
     47    // XXX EAM : why was this being done?
     48    // pmCell *tmpCell = pmCellInChip(chipCoord, cell->parent);
     49    PS_ASSERT_PTR_NON_NULL(cell->toChip, NULL);
     50    psPlaneTransform *tmpChipToCell = p_psPlaneTransformLinearInvert(cell->toChip);
     51    PS_ASSERT_PTR_NON_NULL(tmpChipToCell, NULL);
     52    cellCoord = psPlaneTransformApply(cellCoord, tmpChipToCell, chipCoord);
     53    psFree(tmpChipToCell);
     54    return(cellCoord);
     55}
     56
     57bool psastroProjectFPA (pmFPA *fpa, char *starlist, bool toSky) {
     58
     59    bool status;
     60
     61    for (int i = 0; i < fpa->chips->n; i++) {
     62        pmChip *chip = fpa->chips->data[i];
     63        for (int j = 0; j < chip->cells->n; j++) {
     64            pmCell *cell = chip->cells->data[j];
     65            for (int k = 0; k < cell->readouts->n; k++) {
     66                pmReadout *readout = cell->readouts->data[k];
     67                psArray *stars = psMetadataLookupPtr (&status, readout->analysis, starlist);
     68                if (toSky) {
     69                    psastroProjectRawstars (stars, readout);
     70                } else {
     71                    psastroProjectRefstars (stars, readout);
     72                }
     73            }
     74        }
     75    }
     76    return true;
     77}
     78 
     79bool psastroProjectRawstars (psArray *stars, pmReadout *readout) {
     80
     81    pmCell *cell = readout->parent;
     82    pmChip *chip = cell->parent;
     83    pmFPA  *fpa  = chip->parent;
     84
     85    for (int i = 0; i < stars->n; i++) {
     86        pmAstromObj *star = stars->data[i];
     87        psCoordReadoutToCell (&star->cell, &star->pix, readout);
     88        psCoordCellToChip (&star->chip, &star->cell, cell);
     89        psCoordChipToFPA (&star->FP, &star->chip, chip);
     90        psCoordFPAToTP (&star->TP, &star->FP, 0.0, 0.0, fpa);
     91        psCoordTPToSky (&star->sky, &star->TP, fpa->projection);
     92    }
     93    return true;
     94}
     95 
     96bool psastroProjectRefstars (psArray *stars, pmReadout *readout) {
     97
     98    pmCell *cell = readout->parent;
     99    pmChip *chip = cell->parent;
     100    pmFPA  *fpa  = chip->parent;
     101
     102    for (int i = 0; i < stars->n; i++) {
     103        pmAstromObj *star = stars->data[i];
     104        psCoordSkyToTP (&star->TP, &star->sky, fpa->projection);
     105        psCoordTPToFPA (&star->FP, &star->TP, 0.0, 0.0, fpa);
     106        psCoordFPAToChip (&star->chip, &star->FP, chip);
     107        psCoordChipToCell_EAM (&star->cell, &star->chip, cell);
     108        psCoordCellToReadout (&star->pix, &star->cell, readout);
     109    }
     110    return true;
     111}
     112 
    3113psArray *pmAstromRadiusMatch (psArray *st1, psArray *st2, psMetadata *config) {
    4114
    5     // assume the lists are sorted, or sort in place?
    6     // sort st1 by P
    7     // sort st2 by P
    8 
    9     double dX, dY;
     115    bool status;
     116    int jStart;
     117    double dX, dY, dR;
     118
     119    double RADIUS = psMetadataLookupF32 (&status, config, "PSASTRO.MATCH.RADIUS");
     120    double RADIUS_SQR = PS_SQR (RADIUS);
     121
     122    // sort both lists by Focal Plane X coord
     123    st1 = psArraySort (st1, pmAstromObjSortByFPX);
     124    st2 = psArraySort (st2, pmAstromObjSortByFPX);
     125
     126    psArray *matches = psArrayAlloc (100);
     127    matches->n = 0;
    10128
    11129    int i = 0;
    12130    int j = 0;
    13 
    14     double RADIUS = psMetadataLookupR32 (&status, myHeader, "ASTROM_MATCH_RADIUS");
    15     double RADIUS_SQR = PS_SQR (RADIUS);
    16 
    17     psArray *matches = psArrayAlloc (100);
    18     matches->n = 0;
    19 
    20131    while ((i < st1->n) && (j < st2->n)) {
    21132
    … …  
    34145           
    35146            dX = ((pmAstromObj *)st1->data[i])->FP.x - ((pmAstromObj *)st2->data[j])->FP.x;
    36             dQ = ((pmAstromObj *)st1->data[i])->FP.y - ((pmAstromObj *)st2->data[j])->FP.y;
    37             dR = dX*dX + dQ*dQ;
     147            dY = ((pmAstromObj *)st1->data[i])->FP.y - ((pmAstromObj *)st2->data[j])->FP.y;
     148            dR = dX*dX + dY*dY;
    38149
    39150            if (dR > RADIUS_SQR) {
    … …  
    44155            // got a match; add to output list
    45156            pmAstromMatch *match = pmAstromMatchAlloc (i, j);
    46             psArrayAdd (matches, match, 100);
     157            psArrayAdd (matches, 100, match);
    47158
    48159            j++;
    … …  
    51162        i++;
    52163    }
    53     return (match);
     164    return (matches);
    54165}
    55166
    56167// take two matched star lists and fit a psPlaneTransform between them
    57168//
    58 psPlaneTransform *pmAstromMatchedListFit (psPlaneTransform *map, psArray *st1, psArray *st2, psArray *match, psMetadata *config) {
     169psPlaneTransform *pmAstromMatchFit (psPlaneTransform *map, psArray *raw, psArray *ref, psArray *match, psMetadata *config) {
     170
     171    bool status;
     172    pmAstromObj *rawStar, *refStar;
     173    pmAstromMatch *pair;
    59174
    60175    if (map == NULL) {
    61         // int nX = psMetadataLookupS32 (&status, myHeader, "CHIP.NX");
    62         // int nY = psMetadataLookupS32 (&status, myHeader, "CHIP.NY");
    63         map = psPlaneTransform (2, 2);
    64     }
    65 
    66     psStats *stats = psStatsAlloc (CLIPPED_MEAN);
    67     stats->nSigma = psMetadataLookupS32 (&status, myHeader, "CHIP.NITER");
    68     stats->nSigma = psMetadataLookupS32 (&status, myHeader, "CHIP.NSIGMA");
    69 
    70     psVector *X = psVectorAlloc (match->n);
    71     psVector *Y = psVectorAlloc (match->n);
    72     psVector *x = psVectorAlloc (match->n);
    73     psVector *y = psVectorAlloc (match->n);
     176        int nX = psMetadataLookupS32 (&status, config, "PSASTRO.CHIP.NX");
     177        if (!status) nX = 1;
     178        int nY = psMetadataLookupS32 (&status, config, "PSASTRO.CHIP.NY");
     179        if (!status) nY = 1;
     180        map = psPlaneTransformAlloc (nX, nY);
     181        // XXX EAM : not clear that we are allowed to use orders > 1
     182    }
     183
     184    psStats *stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
     185    stats->clipSigma = psMetadataLookupF32 (&status, config, "PSASTRO.CHIP.NSIGMA");
     186    stats->clipIter  = psMetadataLookupS32 (&status, config, "PSASTRO.CHIP.NITER");
     187
     188    // XXX EAM : clip fit seems to only work for F32!
     189    // XXX EAM : clip fit fails to handle NULL error
     190    psVector *X = psVectorAlloc (match->n, PS_TYPE_F32);
     191    psVector *Y = psVectorAlloc (match->n, PS_TYPE_F32);
     192    psVector *x = psVectorAlloc (match->n, PS_TYPE_F32);
     193    psVector *y = psVectorAlloc (match->n, PS_TYPE_F32);
     194    psVector *wt = psVectorAlloc (match->n, PS_TYPE_F32);
    74195 
    75196    // take the matched stars, first fit
    76197    for (int i = 0; i < match->n; i++) {
    77198
    78         pair = match->data[i];
    79         ob1 = st1->data[pair->i1];
    80         ob2 = st2->data[pair->i2];
    81 
    82         X->data.F64[i] = ob1->chip.x;
    83         Y->data.F64[i] = ob1->chip.y;
    84 
    85         x->data.F64[i] = ob2->FP.x;
    86         y->data.F64[i] = ob2->FP.y;
     199        pair    = match->data[i];
     200        rawStar = raw->data[pair->i1];
     201        refStar = ref->data[pair->i2];
     202
     203        X->data.F32[i] = rawStar->chip.x;
     204        Y->data.F32[i] = rawStar->chip.y;
     205
     206        x->data.F32[i] = refStar->FP.x;
     207        y->data.F32[i] = refStar->FP.y;
     208
     209        wt->data.F32[i] = 1.0;
    87210    }
    88211
    89212    // no masking, no errors
    90     psVectorClipFitPolynomial2D (map->x, stats, NULL, 0, X, NULL, x, y);
    91     psVectorClipFitPolynomial2D (map->y, stats, NULL, 0, Y, NULL, x, y);
     213    psVectorClipFitPolynomial2D (map->x, stats, NULL, 0, X, wt, x, y);
     214    psVectorClipFitPolynomial2D (map->y, stats, NULL, 0, Y, wt, x, y);
    92215    psFree (x);
    93216    psFree (y);
    94217    psFree (X);
    95218    psFree (Y);
     219    psFree (wt);
    96220    psFree (stats);
    97221
    … …  
    103227psArray *pmAstromRotateObj (psArray *old, psPlane center, double angle) {
    104228
     229    double X, Y;
    105230    pmAstromObj *newObj;
     231    pmAstromObj *oldObj;
    106232
    107233    psArray *new = psArrayAlloc (old->n);
    108234
    109     cs = cos(angle);
    110     sn = sin(angle);
    111     xCenter = center.x;
    112     yCenter = center.y;
     235    double cs = cos(angle);
     236    double sn = sin(angle);
     237    double xCenter = center.x;
     238    double yCenter = center.y;
    113239   
    114240    for (int i = 0; i < old->n; i++) {
    … …  
    120246        Y = oldObj->FP.y - yCenter;
    121247       
    122         newObj->FP.x = X*cs + Y*sn;
    123         newObj->FP.y = Y*cs - X*sn;
     248        newObj->FP.x = X*cs + Y*sn + xCenter;
     249        newObj->FP.y = Y*cs - X*sn + yCenter;
    124250       
    125251        new->data[i] = newObj;
    … …  
    148274  obj->TP.y = 0;
    149275
    150   obj->sky.x = 0;
    151   obj->sky.y = 0;
    152 
    153   obj->mag = 0;
     276  obj->sky.r = 0;
     277  obj->sky.d = 0;
     278
     279  obj->Mag = 0;
     280  obj->dMag = 0;
    154281
    155282  return (obj);

  • trunk/psastro/src/pmAstrom.h

    r5510 r5560  
    33# include <unistd.h>   // for unlink
    44# include <pslib.h>
    5 
    6 psArray *pmAstromGridMatch (psArray *s1, psArray *s2, pmAstromGridMatchOpt *opt);
     5# include <pmAstrometry.h>
    76
    87typedef struct {
    … …  
    2221
    2322typedef struct {
    24     double minAngle;
    25     double maxAngle;
    26     double delAngle;
    27 } pmAstromOpt;
    28 
    29 typedef struct {
    3023    psPlane center;
    3124    psPlane offset;
    32     double angle;
    33     double minMetric;
    34     double minVar;
    35     int    nMatch;
    36 } pmAstromGridMatchStats;
     25    double  angle;
     26    double  minMetric;
     27    double  minVar;
     28    int     nMatch;
     29} pmAstromStats;
    3730
     31pmAstromObj      *pmAstromObjAlloc (void);
     32pmAstromObj      *pmAstromObjCopy (pmAstromObj *old);
     33pmAstromMatch    *pmAstromMatchAlloc (int i1, int i2);
     34psArray          *pmAstromRadiusMatch (psArray *st1, psArray *st2, psMetadata *config);
     35psArray          *pmAstromRotateObj (psArray *old, psPlane center, double angle);
     36psPlaneTransform *pmAstromMatchFit (psPlaneTransform *map, psArray *st1, psArray *st2, psArray *match, psMetadata *config);
     37
     38int               pmAstromObjSortByFPX (const void **a, const void **b);
     39
     40bool              psastroProjectFPA (pmFPA *fpa, char *starlist, bool toSky);
     41bool              psastroProjectRawstars (psArray *stars, pmReadout *readout);
     42bool              psastroProjectRefstars (psArray *stars, pmReadout *readout);
     43
     44psPlane* psCoordChipToCell_EAM(
     45    psPlane* out,                      ///< a plane struct to recycle. If NULL, a new struct is created
     46    const psPlane* in,                 ///< the Chip coordinate
     47    const pmCell* cell                 ///< the cell of interest
     48);
     49
     50psPlaneTransform *p_psPlaneTransformLinearInvert_EAM(psPlaneTransform *transform);

  • trunk/psastro/src/pmAstromGrid.c

    r5510 r5560  
    11# include "pmAstrom.h"
    22
    3 static double maxOffset;  // maximum allowed offset between lists, in raw pixels
     3static double maxOffpix;  // maximum allowed offset between lists, in raw pixels
    44static double Scale;      // grid pixel scale
    55static double Offset;     // deltas to pixels
    … …  
    77// local function to convert x,y coords to grid bins
    88// it requires the globals defined above
    9 static bool AstromGridBin (int *pOut, int *qOut, double dP, int dQ) {
    10 
    11     if (fabs(dP) > maxOffset) return false;
    12     if (fabs(dQ) > maxOffset) return false;
    13 
    14     *pOut = dP / Scale + Offset;
    15     *qOut = dQ / Scale + Offset;
     9static bool AstromGridBin (int *dx, int *dy, double dX, double dY) {
     10
     11    if (fabs(dX) > maxOffpix) return false;
     12    if (fabs(dY) > maxOffpix) return false;
     13
     14    *dx = dX / Scale + Offset;
     15    *dy = dY / Scale + Offset;
    1616    return true;
    1717}
    1818
    19 // match two star lists
    20 pmAstromGridMatchStats pmAstromGridMatch (psArray *st1, psArray *st2, psMetadata *config) {
    21 
    22     double xMin, xMax, yMin, yMax;
    23     pmAstromObj *ob1, *ob2;
    24 
    25     pmAstromGridMatchStat minStat, newStat;
    26     psPlane center;
    27 
    28     // find center of the st2 field (focal-plane coords)
    29     xMin = yMin = +1e10;
    30     xMax = yMax = -1e10;
    31     for (int i = 0; i < st2->n; i++) {
    32         ob2 = (pmAstromObj *)st2->data[i];
    33         xMin = PS_MIN (ob2->FP.x, xMin);
    34         xMax = PS_MAX (ob2->FP.x, xMax);
    35         yMin = PS_MIN (ob2->FP.y, yMin);
    36         yMax = PS_MAX (ob2->FP.y, yMax);
    37     }
    38     center.x = 0.5*(xMin + xMax);
    39     center.y = 0.5*(yMin + yMax);
    40 
    41     double minAngle = psMetadataLookupF32 (&status, config, "GRID.MIN.ANGLE");
    42     double maxAngle = psMetadataLookupF32 (&status, config, "GRID.MAX.ANGLE");
    43     double delAngle = psMetadataLookupF32 (&status, config, "GRID.DEL.ANGLE");
    44 
    45     minStat.minMetric = 1e10;
    46     for (angle = minAngle; angle <= maxAngle; angle += delAngle) {
    47         st2r = pmAstromRotateObj (st2, center, angle);
    48         newStat = pmAstromGridMatchAngle (st1, st2r, config);
    49         newStat.angle  = angle;
    50         newStat.center = center;
    51         if (newStat.minMetric < minStat.minMetric) {
    52             minStat = newStat;
    53         }
    54         psFree (st2r);
    55     }
    56     return (minStat);
    57 }
    58 
    5919// match the two lists using the binned delta-delta max
    60 pmAstromGridMatchStats pmAstromGridMatchAngle (psArray *st1, psArray *st2, psMetadata *config) {
    61 
     20pmAstromStats pmAstromGridAngle (psArray *raw, psArray *ref, psMetadata *config) {
     21
     22    bool status;
    6223    int nPix;       // size of matching grid
     24    int nPixHalf;   // half-size of matching grid
    6325    double dX, dY;  // offset between a possible matched pair
    6426    int iX, iY;     // corresponding grid bin
    6527
    66     pmAstromObj *ob1, *ob2 // short-cut pointers to the objects
    67     pmAstromGridMatchStats matchStats; // output match statistics
     28    pmAstromObj *ob1, *ob2; // short-cut pointers to the objects
     29    pmAstromStats stats;    // output match statistics
    6830
    6931    // max allowed offset in either X or Y directions
    70     double gridOffset = psMetadataLookupF32 (&status, config, "GRID.OFFSET");
     32    double gridOffset = psMetadataLookupF32 (&status, config, "PSASTRO.GRID.OFFSET");
    7133
    7234    // sampling scale of the grid
    73     double gridScale  = psMetadataLookupF32 (&status, config, "GRID.SCALE");
     35    double gridScale  = psMetadataLookupF32 (&status, config, "PSASTRO.GRID.SCALE");
    7436
    7537    // set the static scaling factors
    76     nPix = (int)(gridOffset / gridScale + 0.5);  // half-grid
    77     nPix = 2*nPix + 1;
     38    nPixHalf = (int)(gridOffset / gridScale + 0.5);  // half-grid
     39    nPix = 2*nPixHalf + 1;                           // full grid width
    7840
    7941    // these are globals used by p_pmAstromGridBin
    80     maxOffpix = gridScale * (nPix + 0.5);
     42    maxOffpix = gridScale * (nPixHalf + 0.5);            // max offset from true center
    8143    Offset    = maxOffpix / gridScale;
    8244    Scale     = gridScale;
    8345
    84     // XXX EAM : can we assume the allocated image is init-ed?
    85     psImage *gridNP = psImageAlloc (nPix, nPix, PS_TYPE_S32);
     46    // images used as accumulators for the loop below
     47    psImage *gridNP = psImageAlloc (nPix, nPix, PS_TYPE_U32);
    8648    psImage *gridDX = psImageAlloc (nPix, nPix, PS_TYPE_F32);
    8749    psImage *gridDY = psImageAlloc (nPix, nPix, PS_TYPE_F32);
    8850    psImage *gridD2 = psImageAlloc (nPix, nPix, PS_TYPE_F32);
    89     psImageInit (gridNP);
    90     psImageInit (gridDX);
    91     psImageInit (gridDY);
    92     psImageInit (gridD2);
     51    psImageInit (gridNP, 0);
     52    psImageInit (gridDX, 0);
     53    psImageInit (gridDY, 0);
     54    psImageInit (gridD2, 0);
    9355
    9456    // short-cut names for grid images
    95     psS32 *NP = gridNP->data.S32;
    96     psF32 *DX = gridDP->data.F32;
    97     psF32 *DY = gridDQ->data.F32;
    98     psF32 *D2 = gridD2->data.F32;
     57    psS32 **NP = gridNP->data.S32;
     58    psF32 **DX = gridDX->data.F32;
     59    psF32 **DY = gridDY->data.F32;
     60    psF32 **D2 = gridD2->data.F32;
    9961
    10062    // accumulate grids for focal plane (L,M) matches
    101     for (int i = 0; i < st1->n; i++) {
    102         ob1 = (pmAstromObj *)st1->data[i];
    103         for (int j = 0; j < st2->n; j++) {
    104             ob2 = (pmAstromObj *)st2->data[i];
     63    for (int i = 0; i < raw->n; i++) {
     64        ob1 = (pmAstromObj *)raw->data[i];
     65        for (int j = 0; j < ref->n; j++) {
     66            ob2 = (pmAstromObj *)ref->data[j];
    10567            dX = ob1->FP.x - ob2->FP.x;
    10668            dY = ob1->FP.y - ob2->FP.y;
    10769
     70            // fprintf (stderr, "dX,dY: %8.2f %8.2f : %8.2f %8.2f : %8.2f %8.2f\n", dX, dY, ob1->FP.x, ob2->FP.x, ob1->FP.y, ob2->FP.y);
    10871            // find bin coordinates for this delta-delta
    109             if (!p_pmAstromGridBin (&iX, &iY, dX, dY)) {
     72            if (!AstromGridBin (&iX, &iY, dX, dY)) {
    11073                continue; // matched pair is too far offset
    111 
    11274            }
    11375
    11476            // accumulate bin stats
    11577            NP[iY][iX] ++;
    116             DP[iY][iX] += dX;
    117             DQ[iY][iX] += dY;
     78            DX[iY][iX] += dX;
     79            DY[iY][iX] += dY;
    11880            D2[iY][iX] += PS_SQR(dX) + PS_SQR(dY);
    11981        }
    … …  
    12486        double minMetric = 1e10;
    12587        double minVar = 1e10;
    126         double minX = -1;
    127         double minY = -1;
     88        int minX = -1;
     89        int minY = -1;
    12890        double metric, var;
    12991
    13092        // find the max pixel
    131         psStats *stats = psStatsAlloc (PS_STAT_MAX);
    132         stats = psImageStats (stats, gridNP);
     93        psStats *imStats = psStatsAlloc (PS_STAT_MAX);
     94        imStats = psImageStats (imStats, gridNP, NULL, 0);
    13395
    13496        // only check bins with at least 1/2 of max bin
    135         int minNpt = 0.5*stats->max;
     97        int minNpts = 0.5*imStats->max;
     98        fprintf (stderr, "minNpts: %d, max: %d\n", minNpts, (int)(imStats->max));
    13699
    137100        // find the 'best' bin
    138         for (int j = 0; j < gridNP->nY; j++) {
    139             for (int i = 0; i < gridNP->nX; i++) {
     101        for (int j = 0; j < gridNP->numRows; j++) {
     102            for (int i = 0; i < gridNP->numCols; i++) {
    140103
    141104                if (NP[j][i] < minNpts) continue;
    … …  
    145108                metric = var / PS_SQR(PS_SQR(NP[j][i]));
    146109           
     110                fprintf (stderr, "try : %f %f (%d pts, %f var, %f met)\n", DX[j][i]/NP[j][i], DY[j][i]/NP[j][i], NP[j][i], var, metric);
     111
    147112                if (metric < minMetric) {
    148113                    minMetric = metric;
    … …  
    155120
    156121        // convert the bin to delta-delta
    157         matchStats.offset.x  = DP[minY][minX] / NP[minY][minX];
    158         matchStats.offset.y  = DQ[minY][minX] / NP[minY][minX];
    159         matchStats.minMetric = minMetric;
    160         matchStats.minVar    = minVar;
    161         matchStats.nMatch    = NP[minY][minX];
    162     }
    163     return (matchStats);
     122        stats.offset.x  = DX[minY][minX] / NP[minY][minX];
     123        stats.offset.y  = DY[minY][minX] / NP[minY][minX];
     124        stats.minMetric = minMetric;
     125        stats.minVar    = minVar;
     126        stats.nMatch    = NP[minY][minX];
     127
     128        // XXX EAM : This routine, and pmAstromGridMatch, need to handle failure cases better
     129    }
     130    return (stats);
     131}
     132
     133// match two star lists
     134pmAstromStats pmAstromGridMatch (psArray *raw, psArray *ref, psMetadata *config) {
     135
     136    bool status;
     137    double xMin, xMax, yMin, yMax;
     138    pmAstromObj *obj;
     139    psArray *rot;
     140
     141    pmAstromStats minStat, newStat;
     142    psPlane center;
     143
     144    // find center of the raw field (focal-plane coords)
     145    xMin = yMin = +1e10;
     146    xMax = yMax = -1e10;
     147    for (int i = 0; i < raw->n; i++) {
     148        obj = (pmAstromObj *)raw->data[i];
     149        xMin = PS_MIN (obj->FP.x, xMin);
     150        xMax = PS_MAX (obj->FP.x, xMax);
     151        yMin = PS_MIN (obj->FP.y, yMin);
     152        yMax = PS_MAX (obj->FP.y, yMax);
     153    }
     154    center.x = 0.5*(xMin + xMax);
     155    center.y = 0.5*(yMin + yMax);
     156
     157    double minAngle = psMetadataLookupF32 (&status, config, "PSASTRO.GRID.MIN.ANGLE");
     158    double maxAngle = psMetadataLookupF32 (&status, config, "PSASTRO.GRID.MAX.ANGLE");
     159    double delAngle = psMetadataLookupF32 (&status, config, "PSASTRO.GRID.DEL.ANGLE");
     160
     161    minStat.minMetric = 1e10;
     162    for (double angle = minAngle; angle <= maxAngle; angle += delAngle) {
     163        rot = pmAstromRotateObj (raw, center, angle);
     164        newStat = pmAstromGridAngle (rot, ref, config);
     165        newStat.angle  = angle;
     166        newStat.center = center;
     167        if (newStat.minMetric < minStat.minMetric) {
     168            minStat = newStat;
     169        }
     170        psFree (rot);
     171    }
     172    fprintf (stderr, "best: %f %f (%d pts, %f var, %f met)\n", newStat.offset.x, newStat.offset.y, newStat.nMatch, newStat.minVar, newStat.minMetric);
     173    return (minStat);
    164174}
    165175
    166176// apply the measured FPA offset and rotation (stat) to the fpa astrom structures
    167 psFPA *pmAstromGridApply (psPlaneTransform *map, pmAstromGridMatchStat stat) {
    168 
    169     // stat.angle, stat.center, stat.offse
    170     // I think i need to know the center reference....
    171 
    172     return (fpa);
     177psPlaneTransform *pmAstromGridApply (psPlaneTransform *map, pmAstromStats stat) {
     178
     179    double cs = cos (stat.angle);
     180    double sn = sin (stat.angle);
     181   
     182    double dx = (map->x->coeff[0][0] - stat.center.x);
     183    double dy = (map->y->coeff[0][0] - stat.center.y);
     184
     185    // new offset
     186    map->x->coeff[0][0] =  cs*dx + sn*dy - stat.offset.x + stat.center.x;
     187    map->y->coeff[0][0] = -sn*dx + cs*dy - stat.offset.y + stat.center.y;
     188
     189    // original rotation matrix
     190    double pc1_1 = map->x->coeff[1][0];
     191    double pc1_2 = map->x->coeff[0][1];
     192    double pc2_1 = map->y->coeff[1][0];
     193    double pc2_2 = map->y->coeff[0][1];
     194
     195    // new rotation matrix
     196    map->x->coeff[1][0] = +cs*pc1_1 + sn*pc2_1;
     197    map->x->coeff[0][1] = +cs*pc1_2 + sn*pc2_2;
     198    map->y->coeff[1][0] = -sn*pc1_1 + cs*pc2_1;
     199    map->y->coeff[0][1] = -sn*pc1_2 + cs*pc2_2;
     200
     201    return (map);
    173202}
    174203

  • trunk/psastro/src/psModUtils.c

    r5506 r5560  
    44
    55// template sample alloc/free pair:
     6# if (0)
    67static void psFooFree (psFoo *foo) {
    78
    … …  
    910  return;
    1011}
    11 
    1212psFoo *psFooAlloc (int i1, int i2) {
    1313
    … …  
    1717  return (foo);
    1818}
     19# endif
    1920
    2021// pmFPA
    … …  
    3132}
    3233
    33 pmFPA *pmFPAAlloc (int i1, int i2) {
     34pmFPA *pmFPAAlloc (void) {
    3435
    3536  pmFPA *fpa = psAlloc (sizeof(pmFPA));
    … …  
    5657}
    5758
    58 pmChip *pmChipAlloc (int i1, int i2) {
     59pmChip *pmChipAlloc (void) {
    5960
    6061  pmChip *chip = psAlloc (sizeof(pmChip));
    … …  
    7980}
    8081
    81 pmCell *pmCellAlloc (int i1, int i2) {
     82pmCell *pmCellAlloc (void) {
    8283
    8384  pmCell *cell = psAlloc (sizeof(pmCell));
    … …  
    9596  if (readout == NULL) return;
    9697
    97   psFree (readout->objects);
     98  psFree (readout->stars);
    9899
    99100  return;
    100101}
    101102
    102 pmReadout *pmReadoutAlloc (int i1, int i2) {
     103pmReadout *pmReadoutAlloc (void) {
    103104
    104105  pmReadout *readout = psAlloc (sizeof(pmReadout));
    … …  
    109110  readout->colBins = 1;
    110111  readout->rowBins = 1;
    111   readout->objects = NULL;
     112  readout->stars = NULL;
    112113
    113114  return (readout);

  • trunk/psastro/src/psastro.c

    r5510 r5560  
    44
    55    psMetadata *header = NULL;
    6     pmAstromGridMatchStat stat;
    7     pmAstromMatch *match;
    86
    97    // load configuration information
    … …  
    119
    1210    // load the input data (cmp file)
    13     psArray *rawstars = pmSourcesReadCMP (&header, argv[1]);
     11    psArray *rawstars = psastroReadCMP (&header, argv[2]);
    1412
    1513    // simple layout interpretation, basic WCS conversion
    16     psFPA *fpa = psastroBuildFPA (header, rawstars);
     14    pmFPA *fpa = psastroBuildFPA (header, rawstars);
    1715
    1816    // limit fit to bright stars only
    19     psFPA *subset = psastroSelectBrightStars (fpa, config);
    20 
    21     // use the header & config info to project rawstars on the focal plane
    22     psastroProjectRawstars (subset);
     17    psastroSelectBrightStars (fpa, config);
    2318
    2419    // fpa and subset point to the same astrometry terms
    25     psastroChipAstrom (subset, config);
     20    psastroChipAstrom (fpa, config);
    2621
    2722    // write out data (cmp file)
    28     psastroWriteCMP (fpa, argv[2]);
     23    psastroWriteCMP (fpa, argv[3]);
     24
     25    psFree (config);
     26    psFree (fpa);
    2927
    3028    exit (0);

  • trunk/psastro/src/psastro.h

    r5509 r5560  
    33# include <unistd.h>   // for unlink
    44# include <pslib.h>
     5# include <pmAstrometry.h>
     6# include "psLibUtils.h"
     7# include "pmAstrom.h"
    58
     9# define fromTPA fromTangentPlane
     10# define toTPA toTangentPlane
     11# define toSky projection
     12
     13# if (0)
    614typedef struct
    715{
    … …  
    1927    psPlaneTransform *fromFPA;          ///< Transformation from FPA to chip coordinates
    2028    psArray *cells;                     ///< The cells
     29    pmFPA *fpa;
    2130}
    2231pmChip;
    … …  
    2837    psArray *readouts;                  ///< The readouts (referred to by number)
    2938    psMetadata *header;                 ///< header always corresponds to a cell
     39    pmChip *chip;
    3040}
    3141pmCell;
    … …  
    3949    int rowBins;                        ///< Amount of binning in y-dimension
    4050    psArray *stars;                     ///< sources detected / measured on readout
     51    pmCell *cell;
    4152}
    4253pmReadout;
     54# endif
    4355
    44 // dummy structure for template code (psFooAlloc, etc)
    45 typedef struct
    46 {
    47     int i;
    48 }
    49 psFoo;
     56pmAstromStats     pmAstromGridAngle (psArray *st1, psArray *st2, psMetadata *config);
     57psPlaneTransform *pmAstromGridApply (psPlaneTransform *map, pmAstromStats stat);
     58pmAstromStats     pmAstromGridMatch (psArray *st1, psArray *st2, psMetadata *config);
     59int               pmAstromObjSortByMag (const void **a, const void **b);
     60bool              pmCellInterpretWCS (pmCell *cell, psMetadata *header) ;
     61pmFPA            *pmFPACopyAstrom (pmFPA *inFPA);
     62psMetadata       *psastroArguments (int *argc, char **argv);
     63pmFPA            *psastroBuildFPA (psMetadata *header, psArray *stars);
     64bool              psastroChipAstrom (pmFPA *fpa, psMetadata *config);
     65psArray          *psastroLoadReference (psMetadata *header, psMetadata *config);
     66psArray          *psastroReadCMP (psMetadata **header, char *filename);
     67bool              psastroSelectBrightStars (pmFPA *fpa, psMetadata *config);
     68bool              psastroWriteCMP (pmFPA *fpa, char *filename);
     69
     70psMetadata *testArguments (int *argc, char **argv);
     71
     72# if (0)
     73pmFPA            *pmFPAAlloc (void);
     74pmChip           *pmChipAlloc (void);
     75pmCell           *pmCellAlloc (void);
     76pmReadout        *pmReadoutAlloc (void);
     77# endif
     78
     79# define SIGN(X)  (((X) == 0) ? 0 : ((fabs((double)(X))) / (X)))
     80
     81psPolynomial2D *psPolynomial2DCopy (psPolynomial2D *input);
     82psPolynomial4D *psPolynomial4DCopy (psPolynomial4D *input);
     83psPlaneDistort *psPlaneDistortCopy (psPlaneDistort *input);
     84psPlaneTransform *psPlaneTransformCopy (psPlaneTransform *input);
     85psProjection *psProjectionCopy (psProjection *input);
     86psPlane *psCoordReadoutToCell (psPlane *cellpix, psPlane *readpix, pmReadout *readout);
     87psPlane *psCoordCellToReadout (psPlane *readpix, psPlane *cellpix, pmReadout *readout);
     88double psPlaneTransformGetRotation (psPlaneTransform *map);
     89
     90psPlaneDistort *psPlaneDistortInvert(psPlaneDistort *distort);
     91psPlaneTransform *psPlaneTransformLinearInvert(psPlaneTransform *transform);

  • trunk/psastro/src/psastroArguments.c

    r5505 r5560  
    1 # include "psphot.h"
     1# include "psastro.h"
    22
    3 static void usage (void) ;
     3static void usage (void);
     4static void usage_test (void);
    45
    5 psMetadata *psphotArguments (int *argc, char **argv) {
     6psMetadata *psastroArguments (int *argc, char **argv) {
    67
    7   int N, Nfail;
    8   int mode = PS_DATA_STRING | PS_META_REPLACE;
     8    unsigned int Nfail;
    99
    10   // basic pslib options
    11   fprintf (stderr, "starting... %s\n", psLibVersion());
    12   psLogSetFormat ("M");
    13   psLogArguments (argc, argv);
    14   psTraceArguments (argc, argv);
     10    // basic pslib options
     11    fprintf (stderr, "starting... %s\n", psLibVersion());
     12    psLogSetFormat ("M");
     13    psLogArguments (argc, argv);
     14    psTraceArguments (argc, argv);
    1515
    16   if (*argc != 2) usage ();
     16    if (*argc != 4) usage ();
    1717
    18   // load config information
    19   psMetadata *config = psMetadataAlloc ();
    20   config = psMetadataConfigParse (config, &Nfail, argv[3], FALSE);
    21   return (config);
     18    // load config information
     19    psMetadata *config = psMetadataAlloc ();
     20    config = psMetadataConfigParse (config, &Nfail, argv[1], FALSE);
     21    return (config);
    2222}
    2323
    2424static void usage (void) {
    25 
    26     fprintf (stderr, "USAGE: psastro (cmpfile)\n");
     25    fprintf (stderr, "USAGE: psastro (config) (input) (output)\n");
    2726    exit (2);
    2827}
     28
     29psMetadata *testArguments (int *argc, char **argv) {
     30
     31    unsigned int Nfail;
     32
     33    // basic pslib options
     34    fprintf (stderr, "starting... %s\n", psLibVersion());
     35    psLogSetFormat ("M");
     36    psLogArguments (argc, argv);
     37    psTraceArguments (argc, argv);
     38
     39    if (*argc != 6) usage_test ();
     40
     41    // load config information
     42    psMetadata *config = psMetadataAlloc ();
     43    config = psMetadataConfigParse (config, &Nfail, argv[1], FALSE);
     44    return (config);
     45}
     46
     47static void usage_test (void) {
     48    fprintf (stderr, "USAGE: psastro-mktest (config) (cmpfile) (reflist) (rawfile) (rawfile)\n");
     49    exit (2);
     50}

  • trunk/psastro/src/psastroBuildFPA.c

    r5510 r5560  
    1818    // allocate the structures
    1919
    20     pmFPA *fpa = pmFPAAlloc ();
     20    pmFPA *fpa = pmFPAAlloc (NULL, NULL);
    2121
    22     fpa->chips = psArrayAlloc (Nchips);
     22    psArrayRealloc (fpa->chips, Nchips);
    2323    for (int i = 0; i < Nchips; i++) {
     24        pmChip *chip = pmChipAlloc (fpa);
    2425
    25         pmChip *chip = pmChipAlloc ();
    26         chip->fpa = fpa; // assign parent fpa (view only; don't free)
     26        psArrayRealloc (chip->cells, Ncells);
     27        for (int j = 0; j < Ncells; j++) {
     28            pmCell *cell = pmCellAlloc (chip);
    2729
    28         chip->cells = psArrayAlloc (Ncells);
    29         for (int j = 0; j < Ncells; j++) {
    30 
    31             pmCell *cell = pmCellAlloc ();
    32             cell->chip = chip;      // assign parent chip (view only; don't free)
    33             cell->header = header;  // XXX EAM : extend this function to take more than header
    34 
     30            // XXX EAM : extend this function to take more than one header
     31            cell->header = header;     
    3532            pmCellInterpretWCS (cell, header);
    3633
    37             cell->readouts = psArrayAlloc (Nreadouts);
     34            psArrayRealloc (cell->readouts, Nreadouts);
    3835            for (int k = 0; k < Nreadouts; k++) {
     36                pmReadout *readout = pmReadoutAlloc (cell);
    3937
    40                 pmReadout *readout = pmReadoutAlloc ();
    41                 readout->stars  = stars;   // XXX EAM : need more than one stars...
     38                // XXX EAM : extend this function to take more than one stars...
     39                psMetadataAdd (readout->analysis, PS_LIST_TAIL, "STARS.FULLSET", PS_DATA_ARRAY, "stars from analysis", stars);
    4240
    43                 cells->readouts->data[j] = readout;
     41                // XXX EAM : this information should be in the header...
     42                readout->col0 = readout->row0 = 0;
     43                readout->colBins = readout->rowBins = 1;
     44
     45                cell->readouts->data[j] = readout;
    4446            }
    45             chips->cells->data[j] = cell;
     47            chip->cells->data[j] = cell;
    4648        }
    4749        fpa->chips->data[i] = chip;
    … …  
    5355bool pmCellInterpretWCS (pmCell *cell, psMetadata *header) {
    5456
     57    pmChip *chip;
     58    pmFPA  *fpa;
     59    bool status;
     60    psProjectionType type;
    5561    float crval1, crval2, crpix1, crpix2, cdelt1, cdelt2;
    5662    float pc1_1, pc1_2, pc2_1, pc2_2;
    5763
    5864    // *** interpret header data, convert to crval(i), etc
    59     char *ctype = pmMetadataLookupPtr (&status, header, "CTYPE2");
     65    char *ctype = psMetadataLookupPtr (&status, header, "CTYPE2");
    6066    if (!status) {
    6167        psLogMsg ("psastro", 2, "warning: no WCS metadata in header\n");
    … …  
    8591
    8692        // test the CROTAi varient:
    87         float rotate = psMetadataLookupF32 (&status, header, "CROTA2");
     93        double rotate = psMetadataLookupF32 (&status, header, "CROTA2");
    8894        if (status) {
    89             Lambda = cdelt2 / cdelt1;
     95            double Lambda = cdelt2 / cdelt1;
    9096            pc1_1 =  cos(rotate*RAD_DEG);
    9197            pc1_2 = -sin(rotate*RAD_DEG) * Lambda;
    … …  
    131137
    132138    // XXX EAM : these must already be set
    133     chip = cell->chip;
    134     fpa = chip->fpa;
     139    chip = cell->parent;
     140    fpa  = chip->parent;
    135141
    136142    // set toChip to identity as default
    137     // XXX EAM : psPlaneTransformAlloc uses nTerm not nOrder (bug 581)
    138     // XXX EAM : define these in the config?
    139     // int nX = psMetadataLookupS32 (&status, myHeader, "CHIP.NX");
    140     // int nY = psMetadataLookupS32 (&status, myHeader, "CHIP.NY");
    141     cell->toChip   = psPlaneTransformAlloc (2, 2);
     143    // XXX EAM : can we actually use higher order?
     144    int nX = psMetadataLookupS32 (&status, header, "PSASTRO.CHIP.NX");
     145    if (!status) nX = 1;
     146    int nY = psMetadataLookupS32 (&status, header, "PSASTRO.CHIP.NY");
     147    if (!status) nY = 1;
     148
     149    cell->toChip   = psPlaneTransformAlloc (1, 1);
    142150    cell->toChip->x->coeff[1][0] = 1;
    143151    cell->toChip->x->mask[1][1]  = 1;
    … …  
    152160
    153161    // XXX EAM : psPlaneTransformAlloc uses nTerm not nOrder (bug 581)
    154     chip->toFPA   = psPlaneTransformAlloc (2, 2);
     162    // XXX EAM : I've fixed this in pslib eam_rel8_b2
     163    chip->toFPA   = psPlaneTransformAlloc (1, 1);
    155164   
    156165    chip->toFPA->x->coeff[0][0] = crpix1;
    … …  
    164173    chip->toFPA->y->mask[1][1]  = 1;
    165174
     175    chip->fromFPA = p_psPlaneTransformLinearInvert (chip->toFPA);
     176
    166177    // set toTPA to identity as default
    167178    // XXX EAM : psPlaneDistortAlloc uses nTerm not nOrder (bug 581)
    168179    if (fpa->toTPA == NULL) {
    169         fpa->toTPA   = psPlaneDistortAlloc (2, 2, 1, 1);
     180        fpa->toTPA   = psPlaneDistortAlloc (1, 1, 0, 0);
    170181        fpa->toTPA->x->coeff[1][0][0][0] = 1;
    171182        fpa->toTPA->x->mask[1][1][0][0]  = 1;
    … …  
    173184        fpa->toTPA->y->coeff[0][1][0][0] = 1;
    174185        fpa->toTPA->y->mask[1][1][0][0]  = 1;
     186        fpa->fromTangentPlane = psPlaneDistortInvert (fpa->toTangentPlane);
    175187    } else {
    176188        psLogMsg ("psastro", 2, "warning: fpa distortion already defined\n");
    … …  
    178190
    179191    // center of projection is (0,0) coordinate of TPA
    180     fpa->toSky = psProjectionAlloc (crval1, crval2, cdelt1, cdelt2, type);
     192    fpa->toSky = psProjectionAlloc (crval1*RAD_DEG, crval2*RAD_DEG, cdelt1*RAD_DEG, cdelt2*RAD_DEG, type);
    181193    return true;
    182194}

  • trunk/psastro/src/psastroUtils.c

    r5510 r5560  
    11# include "psastro.h"
    2 
    3 // crude function to load CMP data (header + ascii)
    4 // XXX EAM : assumes fixed line, expects NSTARS entries
    5 psArrray *pmSourcesReadCMP (psMetadata **header, char *filename) {
    6 
    7     psMetadata *myHeader;
    8     char line[80];
    9 
    10     psFits *fits = psFitsAlloc (filename);
    11     myHeader = psFitsReadHeader (fits);
    12     psFree (fits);
    13 
    14     // how many lines in the header?
    15     // XXX EAM : is this calculation robust?
    16     nLines = myHeader->list->n;
    17     nBytes = nLines * 80;
    18     if (nBytes % 2880) {
    19         nBlock = 1 + (int)(nBytes / 2880);
    20         nBytes = nBlock * 2880;
    21     }
    22 
    23     // re-open, seek to end of header
    24     FILE *f = fopen (filename, "r");
    25     if (f == NULL) {
    26         psLogMsg ("pmSourcesReadCMP", 3, "can't open output file for input %s\n", filename);
    27         return NULL;
    28     }
    29     fseek (f, nBytes, SEEK_SET);
    30 
    31     // prepare array to store data
    32     nStars = psMetadataLookupS32 (&status, myHeader, "NSTARS");
    33     psArray *stars = psArrayAlloc (nStars);
    34     stars->n = 0;
    35 
    36     // we have fixed bytes / line : use that info
    37     for (i = 0; i < nStars; i++) {
    38         fread (line, 1, 67, f);
    39         sscanf (line, "%lf %lf %lf %lf", &X, &Y, &Mag, &dMag);
    40        
    41         pmAstromObj *obj = pmAstromObjAlloc ();
    42         obj->pix.X    = X;
    43         obj->pix.Y    = Y;
    44         obj->pix.Mag  = Mag;
    45         obj->pix.dMag = dMag;
    46         psArrayAdd (stars, obj, 100);
    47     }
    48     fclose (f);
    49 
    50     *header = myHeader;
    51     return (stars);
     2# define RENORM 0
     3
     4bool testWriteRaw (char *filename, psArray *sources);
     5
     6void psastroDumpStars (psArray *sources) {
     7
     8    for (int i = 0; i < sources->n; i++) {
     9        pmAstromObj *star = sources->data[i];
     10
     11        fprintf (stderr, "%8.2f %8.2f   %8.2f %8.2f   %8.2f %8.2f   %8.2f %8.2f   %8.2f %8.2f   %10.6f %10.6f   %8.2f %8.2f\n",
     12                 star->pix.x, star->pix.y,
     13                 star->cell.x, star->cell.y,
     14                 star->chip.x, star->chip.y,
     15                 star->FP.x, star->FP.y,
     16                 star->TP.x, star->TP.y,
     17                 star->sky.r*DEG_RAD, star->sky.d*DEG_RAD,
     18                 star->Mag, star->dMag);
     19    }
    5220}
    5321
    5422// sort by mag (descending)
    55 int psastroSortByMag (const void **a, const void **b)
     23int pmAstromObjSortByMag (const void **a, const void **b)
    5624{
    57     pmSource *A = *(pmSource **)a;
    58     pmSource *B = *(pmSource **)b;
    59 
    60     psF32 fA = (A->moments == NULL) ? 0 : A->moments->SN;
    61     psF32 fB = (B->moments == NULL) ? 0 : B->moments->SN;
    62     if (isnan (fA)) fA = 0;
    63     if (isnan (fB)) fB = 0;
    64 
    65     psF32 diff = fA - fB;
     25    pmAstromObj *A = *(pmAstromObj **)a;
     26    pmAstromObj *B = *(pmAstromObj **)b;
     27
     28    psF32 diff = A->Mag - B->Mag;
    6629    if (diff > FLT_EPSILON) return (-1);
    6730    if (diff < FLT_EPSILON) return (+1);
    … …  
    6932}
    7033
    71 psArray *psastroSelectBrightStars (pmFPA *inFPA, psMetadata *config) {
    72 
    73     pmFPA *sbFPA = pmFPACopy (fpa);
     34bool psastroSelectBrightStars (pmFPA *fpa, psMetadata *config) {
     35
     36    bool status;
    7437
    7538    // add exclusions for objects on some basis?
    76     int MAX_NSTARS = pmMetadataLookupS32 (&status, config, "MAX_NSTARS");
    77 
    78     for (int i = 0; i < inFPA->chips->n; i++) {
    79         pmChip *inChip = inFPA->chips->data[i];
    80         pmChip *sbChip = sbFPA->chips->data[i];
    81         for (int j = 0; j < Ncells; j++) {
    82             pmCell *inCell = inChip->cells->data[j];
    83 o           pmCell *sbCell = sbChip->cells->data[j];
    84             for (int k = 0; k < Nreadouts; k++) {
    85                 pmReadout *inReadout = inCell->readouts->data[k];
    86                 pmReadout *sbReadout = sbCell->readouts->data[k];
    87 
    88                 inReadout->stars = psArraySort (inReadout->stars, psastroSortByMag);
    89 
    90                 nSubset = PS_MIN (MAX_NSTARS, stars->n);
    91 
    92                 psArray *subset = psArrayAlloc (nSubset);
    93 
    94                 for (int i = 0; i < nSubset; i++) {
    95                     subset->data[i] = inReadout->stars->data[i];
    96                 }
    97                 sbReadout->stars = subset;
    98             }
    99         }
    100     }
    101     return (subset);
    102 }
    103 
    104 // fake this by loading from a text file
    105 psArray *psastroLoadReference (char *filename) {
    106 
    107     psMetadata *myHeader;
    108     char line[80];
    109 
    110     psFits *fits = psFitsAlloc (filename);
    111     myHeader = psFitsReadHeader (fits);
    112     psFree (fits);
    113 
    114     // how many lines in the header?
    115     // XXX EAM : is this calculation robust?
    116     nLines = header->list->n;
    117     nBytes = nLines * 80;
    118     if (nBytes % 2880) {
    119         nBlock = 1 + (int)(nBytes / 2880);
    120         nBytes = nBlock * 2880;
    121     }
    122 
    123     // re-open, seek to end of header
    124     FILE *f = fopen (filename, "r");
    125     if (f == NULL) {
    126         psLogMsg ("pmSourcesReadCMP", 3, "can't open output file for input %s\n", filename);
    127         return NULL;
    128     }
    129     fseek (f, nBytes, SEEK_SET);
    130 
    131     // prepare array to store data
    132     nStars = psMetadataLookupS32 (&status, myHeader, "NSTARS");
    133     psArray *stars = psArrayAlloc (nStars);
    134     stars->n = 0;
    135 
    136     // we have fixed bytes / line : use that info
    137     for (i = 0; i < nStars; i++) {
    138         fread (line, 1, 67, f);
    139         sscanf (line, "%lf %lf %lf %lf", &X, &Y, &Mag, &dMag);
    140        
    141         pmAstromObj *obj = pmAstromObjAlloc (X, Y, Mag, dMag);
    142         psArrayAdd (stars, obj, 100);
    143     }
    144     fclose (f);
    145 
    146     *header = myHeader;
    147     return (stars);
    148 }
    149 
    150 pmFPA *pmFPACopy (pmFPA *inFPA) {
    151 
    152     pmFPA *fpa = pmFPAAlloc ();
    153 
    154     fpa->toSky   = psMemCopy (inFPA->toSky);
    155     fpa->toTPA   = psMemCopy (inFPA->toTPA);
    156     fpa->fromTPA = psMemCopy (inFPA->fromTPA);
    157 
    158     fpa->chips = psArrayAlloc (inFPA->chips->n);
    159     for (int i = 0; i < inFPA->chips->n; i++) {
    160 
    161         pmChip *inChip = inFPA->chips->data[i];
    162 
    163         pmChip *chip = pmChipAlloc ();
    164         chip->fpa = fpa; // assign parent fpa (view only; don't free)
    165 
    166         chip->toFPA = psMemCopy (inChip->toFPA);
    167         chip->fromFPA = psMemCopy (inChip->fromFPA);
    168 
    169         chip->cells = psArrayAlloc (inChip->cells->n);
    170         for (int j = 0; j < Ncells; j++) {
    171 
    172             pmCell *inCell = inChip->cells->data[j];
    173 
    174             pmCell *cell = pmCellAlloc ();
    175             cell->chip = chip;      // assign parent chip (view only; don't free)
    176 
    177             cell->header = psMemCopy (inCell->header);
    178             cell->toChip = psMemCopy (inCell->toChip);
    179 
    180             cell->readouts = psArrayAlloc (inCell);
    181             for (int k = 0; k < Nreadouts; k++) {
    182 
    183                 pmReadout *inReadout = inCell->readouts->data[k];
    184 
    185                 pmReadout *readout = pmReadoutAlloc ();
    186 
    187                 *readout = *inReadout;
    188                 readout->stars = NULL;
    189                
    190                 cell->readouts->data[k] = readout;
    191             }
    192             chip->cells->data[j] = cell;
    193         }
    194         fpa->chips->data[i] = chip;
    195     }
    196     return (fpa);
    197 }
    198 
    199 bool psastroProjectFPA (pmFPA *fpa, bool toSky) {
     39    int MAX_NSTARS = psMetadataLookupS32 (&status, config, "PSASTRO.STARS.MAX");
    20040
    20141    for (int i = 0; i < fpa->chips->n; i++) {
    … …  
    20545            for (int k = 0; k < cell->readouts->n; k++) {
    20646                pmReadout *readout = cell->readouts->data[k];
    207                 if (toSky) {
    208                     psastroProjectRawstars (readout->stars, readout);
    209                 } else {
    210                     psastroProjectRefstars (readout->stars, readout);
     47
     48                psArray *stars = psMetadataLookupPtr (&status, readout->analysis, "STARS.FULLSET");
     49                stars = psArraySort (stars, pmAstromObjSortByMag);
     50
     51                int nSubset = PS_MIN (MAX_NSTARS, stars->n);
     52                psArray *subset = psArrayAlloc (nSubset);
     53
     54                for (int i = 0; i < nSubset; i++) {
     55                    subset->data[i] = stars->data[i];
    21156                }
     57                psMetadataAdd (readout->analysis, PS_LIST_TAIL, "STARS.SUBSET", PS_DATA_ARRAY, "stars from analysis", subset);
    21258            }
    21359        }
    … …  
    21561    return true;
    21662}
    217  
    218 bool psastroProjectRawstars (psArray *stars, pmReadout *readout) {
    219 
    220     pmCell *cell = readout->cell;
    221     pmChip *chip = cell->chip;
    222     pmFPA  *fpa  = chip->cell;
    223 
    224     for (int i = 0; i < readout->stars->n; i++) {
    225         pmAstromObj *star = readout->stars->data[i];
    226         psCoordReadoutToCell (&star->cell, &star->pix, readout);
    227         psCoordCellToChip (&star->chip, &star->cell, cell);
    228         psCoordChipToFP (&star->FP, &star->chip, chip);
    229         psCoordFPtoTP (&star->TP, &star->FP, fpa);
    230         psCoordTPtoSky (&star->sky, &star->TP, fpa);
    231     }
    232     return true;
    233 }
    234  
    235 bool psastroProjectRefstars (psArray *stars, pmReadout *readout) {
    236 
    237     pmCell *cell = readout->cell;
    238     pmChip *chip = cell->chip;
    239     pmFPA  *fpa  = chip->cell;
    240 
    241     for (int i = 0; i < stars->n; i++) {
    242         pmAstromObj *star = stars->data[i];
    243 
    244         psCoordSkytoTP (&star->TP, &star->sky, fpa);
    245         psCoordTPtoFP (&star->FP, &star->TP, fpa);
    246         psCoordFPtoChip (&star->chip, &star->FP, chip);
    247         psCoordChipToCell (&star->cell, &star->chip, cell);
    248         psCoordReadoutToCell (&star->pix, &star->cell, readout);
    249     }
    250     return true;
    251 }
    252  
     63
     64pmFPA *pmFPACopyAstrom (pmFPA *inFPA) {
     65
     66    pmFPA *fpa = pmFPAAlloc (NULL, NULL);
     67
     68    // copy FPA astrometry data
     69    fpa->toSky   = psProjectionCopy (inFPA->toSky);
     70    fpa->toTPA   = psPlaneDistortCopy (inFPA->toTPA);
     71    fpa->fromTPA = psPlaneDistortCopy (inFPA->fromTPA);
     72
     73    psArrayRealloc (fpa->chips, inFPA->chips->n);
     74    for (int i = 0; i < inFPA->chips->n; i++) {
     75        pmChip *inChip = inFPA->chips->data[i];
     76        pmChip *chip = pmChipAlloc (fpa);
     77
     78        // copy Chip astrometry data
     79        chip->toFPA = psPlaneTransformCopy (inChip->toFPA);
     80        chip->fromFPA = psPlaneTransformCopy (inChip->fromFPA);
     81
     82        psArrayRealloc (chip->cells, inChip->cells->n);
     83        for (int j = 0; j < inChip->cells->n; j++) {
     84            pmCell *inCell = inChip->cells->data[j];
     85            pmCell *cell = pmCellAlloc (chip);
     86
     87            // cell.header is a view on inCell.header
     88            cell->header = psMemCopy (inCell->header);
     89
     90            // copy Cell astrometry data
     91            cell->toChip = psPlaneTransformCopy (inCell->toChip);
     92
     93            psArrayRealloc (cell->readouts, inCell->readouts->n);
     94            for (int k = 0; k < inCell->readouts->n; k++) {
     95                pmReadout *inReadout = inCell->readouts->data[k];
     96                pmReadout *readout = pmReadoutAlloc (cell);
     97
     98                // copy Readout data
     99                *readout = *inReadout;
     100
     101                cell->readouts->data[k] = readout;
     102            }
     103            chip->cells->data[j] = cell;
     104        }
     105        fpa->chips->data[i] = chip;
     106    }
     107    return (fpa);
     108}
     109
    253110// measure per-chip astrometry terms
    254111bool psastroChipAstrom (pmFPA *fpa, psMetadata *config) {
    255112
    256     pmFPA *raw = pmFPACopy (fpa);
     113    bool status;
     114    psArray *match;
     115    pmAstromStats stats;
     116
     117    // save the raw astrometry for later reference
     118    pmFPA *raw = pmFPACopyAstrom (fpa);
    257119
    258120    // first pass: measure the per-chip solutions, modify the chip.toFPA terms
    … …  
    282144                pmReadout *readout = cell->readouts->data[k];
    283145
    284                 // use the header & config info to project refstars on the focal plane
     146                // pull out the SUBSET rawstars (a view)
     147                psArray *rawstars = psMetadataLookupPtr (&status, readout->analysis, "STARS.SUBSET");
     148
     149                // project the rawstars to the current best guess astrometry
     150                psastroProjectRawstars (rawstars, readout);
     151
     152                // use the header & config info to project refstars onto the focal plane
    285153                psastroProjectRefstars (refstars, readout);
    286154
     155                testWriteRaw ("ref.inp", refstars);
     156                testWriteRaw ("raw.inp", rawstars);
     157
     158                // fprintf (stderr, "rawstars:\n");
     159                // psastroDumpStars (rawstars);
     160                // fprintf (stderr, "refstars:\n");
     161                // psastroDumpStars (refstars);
     162
    287163                // find initial offset / rotation
    288                 stat = pmAstromGridMatch (readout->stars, refstars, config);
     164                stats = pmAstromGridMatch (rawstars, refstars, config);
    289165
    290166                // adjust the chip.toFPA terms only
    291                 pmAstromGridApply (chip->toFPA, stat);
     167                pmAstromGridApply (chip->toFPA, stats);
     168                chip->fromFPA = p_psPlaneTransformLinearInvert (chip->toFPA);
    292169
    293170                // use fit result to re-project rawstars
    294                 psastroProjectRawstars (readout->stars, readout);
    295                 psastroProjectRefstars (refstars,       readout);
     171                psastroProjectRawstars (rawstars, readout);
     172                psastroProjectRefstars (refstars, readout);
     173
     174                testWriteRaw ("ref.dat", refstars);
     175                testWriteRaw ("raw.dat", rawstars);
    296176   
    297177                // use small radius to match stars
    298                 match = pmAstromRadiusMatch (rawstars, refstars, options);
     178                match = pmAstromRadiusMatch (rawstars, refstars, config);
    299179
    300180                // fit astrometric terms
    301                 pmAstromMatchedListFit (chip->toFPA, readout->stars, refstars, match, options);
     181                pmAstromMatchFit (chip->toFPA, rawstars, refstars, match, config);
    302182            }
    303183        }
    … …  
    306186    // second stage: re-normalize the chip terms, passing the
    307187    // average rotation and offset values to the fpa.toSky
    308            
    309 # if (0)
    310     // code disabled for now
    311     // this is not as trivial as it seems at first:
    312     // applying a new rotation to the FPA implies changes to the
    313     // (x,y) reference coordinates for the chips as well as a rotation
    314     // I need to calculate both the x,y offset effect and the rotation
    315     // in one step.  what about the effect of higher order terms in either
    316     // the chip.toFPA or the fpa.toTPA?
    317     // this calculation also needs to account for scale effects between the
    318     // fpa and the chips.
    319 
    320     // calculate the average rotation and boresite offset relative to raw
    321     for (int i = 0; i < fpa->chips->n; i++) {
    322         pmChip *iChip = raw->chips->data[i];
    323         pmChip *oChip = fpa->chips->data[i];
    324 
    325         T1 = psPlaneTransformGetRotation (iChip->toFPA);
    326         T2 = psPlaneTransformGetRotation (oChip->toFPA);
    327         dT += T1 - T2;
    328         dN ++;
    329 
     188    if (RENORM) {
     189
     190        // this code is needed for the mosastro stage, with multiple chip solutions
     191
     192        double dX, dY, dT, dN;
     193        dX = dY = dT = dN = 0;
     194
     195        psPlane origin, P1, P2;
     196        origin.x = 0;
     197        origin.y = 0;
     198
     199        // calculate the average rotation and boresite offset relative to raw
     200        for (int i = 0; i < fpa->chips->n; i++) {
     201            pmChip *iChip = raw->chips->data[i];
     202            pmChip *oChip = fpa->chips->data[i];
     203
     204            // offset of chip
     205            psCoordChipToFPA (&P1, &origin, iChip);
     206            psCoordChipToFPA (&P2, &origin, oChip);
     207            dX += (P2.x - P1.x);
     208            dY += (P2.y - P1.y);
     209
     210            // get parity-independent rotations for old and new solutions
     211            double T1 = psPlaneTransformGetRotation (iChip->toFPA);
     212            double T2 = psPlaneTransformGetRotation (oChip->toFPA);
     213            dT += T2 - T1;
     214            dN ++;
     215        }
     216
     217        dT /= dN;
     218        dX /= dN;
     219        dY /= dN;
     220
     221        // R(T)
     222        double PC1_1 = fpa->toTPA->x->coeff[1][0][0][0];
     223        double PC1_2 = fpa->toTPA->x->coeff[0][1][0][0];
     224        double PC2_1 = fpa->toTPA->y->coeff[1][0][0][0];
     225        double PC2_2 = fpa->toTPA->y->coeff[0][1][0][0];
     226
     227        // R(dT)
     228        double dPC1_1 = +cos (dT);
     229        double dPC1_2 = +sin (dT);
     230        double dPC2_1 = -sin (dT);
     231        double dPC2_2 = +cos (dT);
     232
     233        // R'(T) = R(T) * R(dT)
     234        double pc1_1 = PC1_1*dPC1_1 + PC1_2*dPC2_1;
     235        double pc1_2 = PC1_1*dPC1_2 + PC1_2*dPC2_2;
     236        double pc2_1 = PC2_1*dPC1_1 + PC2_2*dPC2_1;
     237        double pc2_2 = PC2_1*dPC1_2 + PC2_2*dPC2_2;
     238
     239        double det = 1.0 / (pc1_1*pc2_2 - pc1_2*pc2_1);
     240
     241        // R'(-T)  (matrix inverse, not just rotation inverse -- keeps parity)
     242        double pi1_1 = +pc2_2 * det;
     243        double pi1_2 = -pc1_2 * det;
     244        double pi2_1 = -pc2_1 * det;
     245        double pi2_2 = +pc1_1 * det;
     246   
     247        // apply the new modifcations in rotation and boresite
     248        for (int i = 0; i < fpa->chips->n; i++) {
     249            pmChip *oChip = fpa->chips->data[i];
     250
     251            // r(T)
     252            double pr1_1 = oChip->toFPA->x->coeff[1][0];
     253            double pr1_2 = oChip->toFPA->x->coeff[0][1];
     254            double pr2_1 = oChip->toFPA->y->coeff[1][0];
     255            double pr2_2 = oChip->toFPA->y->coeff[0][1];
     256
     257            // ri'(T) = R(T) r(t)
     258            double ri1_1 = PC1_1*pr1_1 + PC1_2*pr2_1;
     259            double ri1_2 = PC1_1*pr1_2 + PC1_2*pr2_2;
     260            double ri2_1 = PC2_1*pr1_1 + PC2_2*pr2_1;
     261            double ri2_2 = PC2_1*pr1_2 + PC2_2*pr2_2;
     262
     263            // r'(T) = R'(-T) ri'(T)
     264            oChip->toFPA->x->coeff[1][0] = pi1_1*ri1_1 + pi1_2*ri2_1;
     265            oChip->toFPA->x->coeff[0][1] = pi1_1*ri1_2 + pi1_2*ri2_2;
     266            oChip->toFPA->y->coeff[1][0] = pi2_1*ri1_1 + pi2_2*ri2_1;
     267            oChip->toFPA->y->coeff[0][1] = pi2_1*ri1_2 + pi2_2*ri2_2;
     268
     269            double dx = PC1_1*oChip->toFPA->x->coeff[0][0] + PC1_2*oChip->toFPA->y->coeff[0][0] + dX;
     270            double dy = PC2_1*oChip->toFPA->x->coeff[0][0] + PC2_2*oChip->toFPA->y->coeff[0][0] + dY;
     271
     272            oChip->toFPA->x->coeff[0][0] = pi1_1*dx + pi1_2*dy;
     273            oChip->toFPA->y->coeff[0][0] = pi2_1*dx + pi2_2*dy;
     274        }
     275
     276        fpa->toTPA->x->coeff[0][0][0][0] -= dX;
     277        fpa->toTPA->y->coeff[0][0][0][0] -= dY;
     278
     279        fpa->toTPA->x->coeff[1][0][0][0] = pc1_1;
     280        fpa->toTPA->x->coeff[0][1][0][0] = pc1_2;
     281        fpa->toTPA->y->coeff[1][0][0][0] = pc2_1;
     282        fpa->toTPA->y->coeff[0][1][0][0] = pc2_2;
     283    }
     284    return true;
     285}
     286
     287// returns the rotation term, forcing positive parity
     288double psPlaneTransformGetRotation (psPlaneTransform *map) {
     289
     290    if (map->x->nX < 1) return 0;
     291    if (map->x->nY < 1) return 0;
     292
     293    if (map->y->nX < 1) return 0;
     294    if (map->y->nY < 1) return 0;
     295   
     296    double pc1_1 = map->x->coeff[1][0];
     297    double pc1_2 = map->x->coeff[0][1];
     298    double pc2_1 = map->y->coeff[1][0];
     299    double pc2_2 = map->y->coeff[0][1];
     300
     301    double px = SIGN (pc1_1);
     302    double py = SIGN (pc2_2);
     303
     304    // both x and y terms imply an angle. take the average
     305    double t1 = -atan2 (px*pc1_2, px*pc1_1);
     306    double t2 = +atan2 (py*pc2_1, py*pc2_2);
     307   
     308    // careful near -pi,+pi boundary...
     309    if (t1 - t2 > M_PI/2) t2 += 2*M_PI;
     310    if (t2 - t1 > M_PI/2) t1 += 2*M_PI;
     311
     312    double theta = 0.5*(t1 + t2);
     313    while (theta < M_PI) theta += 2*M_PI;
     314    while (theta > M_PI) theta -= 2*M_PI;
     315   
     316    return (theta);
     317}
     318
     319// elixir-style pseudo FITS table (header + ascii list)
     320bool testWriteRaw (char *filename, psArray *sources) {
     321
     322    int i;
     323
     324    // re-open, add data to end of file
     325    FILE *f = fopen (filename, "w");
     326    if (f == NULL) {
     327        psLogMsg ("WriteSourceOBJ", 3, "can't open output file for output %s\n", filename);
     328        return false;
     329    }
     330
     331    for (i = 0; i < sources->n; i++) {
    330332       
    331     }
    332 
    333     dT /= dN;
    334 
    335     psPlaneTransformSetRotation (fpa->toTPA, dT);
    336        
    337     // apply the new modifcations in rotation and boresite
    338     for (int i = 0; i < fpa->chips->n; i++) {
    339         pmChip *iChip = raw->chips->data[i];
    340         pmChip *oChip = fpa->chips->data[i];
    341 
    342         psPlaneTransformSetRotation (oChip->toFPA, -dT);
    343     }
    344 # endif
    345 }
     333        pmAstromObj *star = sources->data[i];
     334
     335        fprintf (f, "%8.2f %8.2f   %8.2f %8.2f   %8.2f %8.2f   %8.2f %8.2f   %8.2f %8.2f   %10.6f %10.6f   %8.2f %8.2f\n",
     336                 star->pix.x, star->pix.y,
     337                 star->cell.x, star->cell.y,
     338                 star->chip.x, star->chip.y,
     339                 star->FP.x, star->FP.y,
     340                 star->TP.x, star->TP.y,
     341                 star->sky.r*DEG_RAD, star->sky.d*DEG_RAD,
     342                 star->Mag, star->dMag);
     343    }
     344    fclose (f);
     345    return true;
     346}
Note: See TracChangeset for help on using the changeset viewer.