IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Dec 13, 2015, 5:50:46 AM (11 years ago)
Author:
eugene
Message:

merged changes from trunk

Location:
branches/eam_branches/ipp-20151113
Files:
13 edited
2 copied

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20151113

  • branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/Makefile

    r39223 r39266  
    2525$(SRC)/cneedles.$(ARCH).o                  \
    2626$(SRC)/cplot.$(ARCH).o             \
     27$(SRC)/crotation.$(ARCH).o         \
    2728$(SRC)/csystem.$(ARCH).o           \
    2829$(SRC)/ctimes.$(ARCH).o    \
    … …  
    3738$(SRC)/fitpm_irls.$(ARCH).o  \
    3839$(SRC)/fitplx_irls.$(ARCH).o  \
     40$(SRC)/astrom_ops.$(ARCH).o        \
    3941$(SRC)/fixwrap.$(ARCH).o           \
    4042$(SRC)/fixcols.$(ARCH).o           \

  • branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/cdensify.c

    r38463 r39266  
    8484  gfits_free_matrix (&bf[0].matrix);
    8585  gfits_free_header (&bf[0].header);
    86   CreateBuffer (bf, Nx, Ny, -32, 0.0, 1.0);
     86  if (!CreateBuffer (bf, Nx, Ny, -32, 0.0, 1.0)) return FALSE;
    8787  strcpy (bf[0].file, "(empty)");
    8888  PutCoords (&newcoords, &bf[0].header);

  • branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/cdhistogram.c

    r38062 r39266  
    101101  gfits_free_matrix (&bf[0].matrix);
    102102  gfits_free_header (&bf[0].header);
    103   CreateBuffer3D (bf, Nx, Ny, Nz, -32, 0.0, 1.0);
     103  if (!CreateBuffer3D (bf, Nx, Ny, Nz, -32, 0.0, 1.0)) return FALSE;
    104104  strcpy (bf[0].file, "(empty)");
    105105  PutCoords (&newcoords, &bf[0].header);

  • branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/fitplx.c

    r38986 r39266  
    11# include "astro.h"
    2 # define J2000 51544.5       /* Modified Julian date at standard epoch J2000 */
    3 
    4 # define ESCAPE(MSG,...) {                      \
    5     gprint (GP_ERR, MSG, __VA_ARGS__);          \
    6     return FALSE; }
    7 
    8 typedef struct {
    9   double *X;
    10   double *Y;
    11   double *t;
    12   double *pX;
    13   double *pY;
    14   double *dX;
    15   double *dY;
    16   int *index;
    17   int Npts;
    18 } PlxFitData;
    19 
    20 typedef struct {
    21   double Ro, dRo;
    22   double Do, dDo;
    23 
    24   double uR, duR;
    25   double uD, duD;
    26 
    27   double p, dp;
    28 
    29   double chisq;
    30   int Nfit;
    31   int getChisq;
    32 } PlxFit;
    33 
    34 int VectorRobustStats (Vector *vector, double *median, double *sigma);
    35 double VectorFractionInterpolate (double *values, float fraction, int Npts);
    36 
    37 int PlxSetMeanEpoch (double *R, double *D, double *T, double *Rmean, double *Dmean, double *Tmean, int *mask, int Ntotal);
    38 int PlxSetEpochPosition (PlxFitData *fitdata, double *R, double *D, double *dR, double *dD, double *T, int *mask, int Ntotal, Coords *coords, double Tmean);
    39 int PlxOutlierClip (PlxFitData *fitdata, int *mask, int Noutlier, float dPsigMax, Vector *dPvec, int VERBOSE);
    40 
    41 int PlxFitDataAlloc (PlxFitData *data, int N);
    42 void PlxFitDataFree (PlxFitData *data);
    43 int PlxBootstrapResample (PlxFitData *src, PlxFitData *tgt);
    44 
    45 int FitPMandPar (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts, int VERBOSE);
    46 int sun_ecliptic (double mjd, double *lambda, double *beta, double *epsilon, double *Radius);
    47 int ParFactor (double *pR, double *pD, double RA, double DEC, double Time);
    482
    493int fitplx (int argc, char **argv) {
    … …  
    10054    gprint (GP_ERR, "  -outlier-tests Nsamples dPsigMax : run Nsample bootstrap-resamples to define the path deviations and reject based on dPsigMax\n");
    10155    gprint (GP_ERR, "  -dPsig vec : save path deviations in vec\n");
     56    gprint (GP_ERR, "  -mask mask : excluded points are marked with a 0 mask value\n");
    10257    return (FALSE);
    10358  }
    … …  
    148103  PlxFit fit; memset (&fit, 0, sizeof(PlxFit));
    149104
    150   // determine dPsig for detections based on Noutlier attempts
     105  // determine dPsig for detections based on Noutlier attempts (mask is updated to mark the bad stars, mask == 0)
    151106  if (Noutlier) {
    152107    int clipRetry = TRUE;
    153108    for (i = 0; clipRetry && (i < 3); i++) {
     109      // XXX NOTE This will segfault is mask is not supplied:
    154110      clipRetry = !PlxOutlierClip (&fitdata, mask, Noutlier, dPsigMax, dPvec, VERBOSE);
    155111
    … …  
    163119  for (i = 0; (VERBOSE == 2) && (i < fitdata.Npts); i++) {
    164120    int n = fitdata.index[i];
    165     int maskValue = mask ? mask[n] : 0;
     121    int maskValue = mask ? mask[n] : 1;
    166122    fprintf (stderr, "%f %f : %f %d : %f %f %f\n", R[n], D[n], T[n], maskValue, fitdata.t[i], fitdata.X[i], fitdata.Y[i]);
    167123  }
    … …  
    401357  fit[0].Nfit = Npts;
    402358  return (TRUE);
    403 }
    404 
    405 /* Low precision formulae for the sun, from Astro. Almanac p. C5 (2012) */
    406 int sun_ecliptic (double mjd, double *lambda, double *beta, double *epsilon, double *Radius) {
    407 
    408   double n = mjd - J2000;                          // day number relative to standard epoch
    409   double L = 280.460 + 0.9856474 * n;              // mean solar longitute (corr. for aberration)
    410   double g = (357.528 + 0.9856003 * n)*RAD_DEG;    // Mean anomaly
    411 
    412   *lambda = L + 1.915 * sin(g) + 0.020 * sin(2*g); // solar longitude in degrees
    413   *beta = 0.0;                                     // approx latitude
    414   *epsilon = (23.439 - 0.0000004 * n);             // obliquity of ecliptic in degrees
    415   *Radius = 1.00014 - 0.01671*cos(g) - 0.00014*cos(2*g); // earth-to-sun dist in AU
    416   return TRUE;
    417 }
    418 
    419 /* given RA, DEC, Time, calculate the parallax factor */
    420 // RA,DEC are decimal degrees
    421 // Time is MJD
    422 int ParFactor (double *pR, double *pD, double RA, double DEC, double Time) {
    423 
    424   double lambda, beta, epsilon, Radius;
    425 
    426   // Time must be mjd
    427   sun_ecliptic (Time, &lambda, &beta, &epsilon, &Radius);
    428 
    429   double lambda_rad = lambda*RAD_DEG;
    430   double epsilon_rad = epsilon*RAD_DEG;
    431   double RA_rad = RA*RAD_DEG;
    432   double DEC_rad = DEC*RAD_DEG;
    433 
    434   double x = Radius*cos(lambda_rad);
    435   double y = Radius*cos(epsilon_rad)*sin(lambda_rad);
    436   double z = Radius*sin(epsilon_rad)*sin(lambda_rad);
    437 
    438   *pR = +(y*cos(RA_rad) - x*sin(RA_rad));
    439   *pD = -(y*sin(RA_rad) + x*cos(RA_rad))*sin(DEC_rad) + z*cos(DEC_rad);
    440 
    441   return TRUE;
    442 }
    443 
    444 // allocate arrays but not the container
    445 int PlxFitDataAlloc (PlxFitData *data, int N) {
    446 
    447   data->Npts = N;
    448   ALLOCATE (data->X, double, N);
    449   ALLOCATE (data->Y, double, N);
    450   ALLOCATE (data->dX, double, N);
    451   ALLOCATE (data->dY, double, N);
    452   ALLOCATE (data->t, double, N);
    453   ALLOCATE (data->pX, double, N);
    454   ALLOCATE (data->pY, double, N);
    455   ALLOCATE (data->index, int, N);
    456   return TRUE;
    457 }
    458 
    459 void PlxFitDataFree (PlxFitData *data) {
    460   FREE (data->X);
    461   FREE (data->Y);
    462   FREE (data->dX);
    463   FREE (data->dY);
    464   FREE (data->t);
    465   FREE (data->pX);
    466   FREE (data->pY);
    467   FREE (data->index);
    468359}
    469360
    … …  
    481372    tgt->pX[i] = src->pX[N];
    482373    tgt->pY[i] = src->pY[N];
     374
     375    // *** make this optional?
     376    tgt->Wx[i] = src->Wx[N];
     377    tgt->Wy[i] = src->Wy[N];
    483378  }
    484379  return TRUE;
    … …  
    539434    pYmin = MIN (pYmin, fitdata->pY[Nsubset]);
    540435    pYmax = MAX (pYmax, fitdata->pY[Nsubset]);
     436
     437    fitdata->Wx[Nsubset] = 1.0;
     438    fitdata->Wy[Nsubset] = 1.0;   
    541439    fitdata->index[Nsubset] = i;
    542440    Nsubset++;

  • branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/fitpm.c

    r37807 r39266  
    11# include "astro.h"
    2 # define J2000 51544.5       /* Modified Julian date at standard epoch J2000 */
    3 
    4 # define ESCAPE(MSG,...) {                      \
    5     gprint (GP_ERR, MSG, __VA_ARGS__);          \
    6     return FALSE; }
    7 
    8 typedef struct {
    9   double Ro, dRo;
    10   double Do, dDo;
    11 
    12   double uR, duR;
    13   double uD, duD;
    14 
    15   double chisq;
    16   int Nfit;
    17 } PMFit;
    18 
    19 int FitPMonly (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, int VERBOSE);
    202
    213int fitpm (int argc, char **argv) {
    … …  
    122104  }
    123105
    124   PMFit fit;
     106  PlxFit fit;
    125107  if (!FitPMonly (&fit, X, dX, Y, dY, t, n, VERBOSE)) {
    126108    return FALSE;
    … …  
    158140
    159141/* do we want an init function which does the alloc and a clear function to free? */
    160 int FitPMonly (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, int VERBOSE) {
     142int FitPMonly (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, int VERBOSE) {
    161143
    162144  int i;
    … …  
    241223  fit[0].duD = sqrt(A[3][3]);
    242224 
     225  fit[0].p   = 0.0;
     226  fit[0].dp  = NAN;
     227
    243228  // add up the chi square for the fit
    244229  chisq = 0.0;

  • branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/fitpm_irls.c

    r39223 r39266  
    11# include "astro.h"
    2 # define J2000 51544.5       /* Modified Julian date at standard epoch J2000 */
    3 
    4 # define ESCAPE(MSG,...) {                      \
    5     gprint (GP_ERR, MSG, __VA_ARGS__);          \
    6     return FALSE; }
    7 
    8 typedef struct {
    9   double Ro, dRo;
    10   double Do, dDo;
    11 
    12   double uR, duR;
    13   double uD, duD;
    14  
    15   double chisq;
    16   int Nfit;
    17 } PMFit_IRLS;
    18 
    19 int FitPMonly_IRLS (PMFit_IRLS *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, int VERBOSE);
    20 int IRLS_converged (PMFit_IRLS *fit);
    21 int weighted_LS (double *T, double *X, double *WX, double *Y, double *WY, int Npts,
    22                  double **A, double **B, int VERBOSE);
    23 double weight_cauchy (double x);
    24 double dpsi_cauchy (double x);
    25 double MAD(double *in, int N);
    26      
    272
    283int fitpm_irls (int argc, char **argv) {
    … …  
    129104  }
    130105
    131   PMFit_IRLS fit;
     106  PlxFit fit;
    132107  if (!FitPMonly_IRLS (&fit, X, dX, Y, dY, t, n, VERBOSE)) {
    133108    return FALSE;
    … …  
    165140
    166141/* do we want an init function which does the alloc and a clear function to free? */
    167 int FitPMonly_IRLS (PMFit_IRLS *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, int VERBOSE) {
     142int FitPMonly_IRLS (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, int VERBOSE) {
    168143
    169144  int i,j;
    … …  
    222197 
    223198  // Solve OLS equation 
    224   if (!weighted_LS(T,X,Wx,Y,Wy,Npts,
     199  if (!weighted_LS_PM(T,X,Wx,Y,Wy,Npts,
    225200                   A,B,VERBOSE)) {
    226201    // Handle fail case
    … …  
    267242
    268243    // Solve
    269     if (!weighted_LS(T,X,Wx,Y,Wy,Npts,
     244    if (!weighted_LS_PM(T,X,Wx,Y,Wy,Npts,
    270245                     A,B,VERBOSE)) {
    271246      // Handle fail case
    … …  
    284259      u[i] = sqrt(SQ(rx[i] / dX[i]) + SQ(ry[i] / dY[i]));
    285260    }
    286     sigma_hat = MAD(u,Npts) / 0.6745;
     261    sigma_hat = MedianAbsDeviation(u,Npts) / 0.6745;
    287262   
    288263    // Check convergence
    … …  
    310285  double Sum_Wx, Sum_Wy;
    311286 
     287  Sum_Wx = 0.0;
     288  Sum_Wy = 0.0;
    312289  ax = 0.0; ay = 0.0;
    313290  bx = 0.0; by = 0.0;
    … …  
    387364}
    388365
    389 
    390 double weight_cauchy (double x) {
    391   double r = x / 2.385;
    392   return (1.0 / (1.0 + SQ(r)));
    393 }
    394 
    395 // dpsi = (d/dx) (x * weight(x))
    396 double dpsi_cauchy (double x) {
    397   double r2 = SQ(x / 2.385);
    398   return ((1.0 - r2) / (SQ(1 + r2)));
    399 }
    400 
    401 
    402 // median absolute deviation
    403 // MAD = median(abs(x - median(x)))
    404 double MAD(double *in, int N) {
    405   double *x;
    406   double median = 0.0;
    407   int i;
    408  
    409   ALLOCATE(x,double,N);
    410   for (i = 0; i < N; i++) {
    411     x[i] = in[i];
    412   }
    413 
    414   dsort(x,N);
    415 
    416   if (N % 2) {
    417     median = 0.5*(x[(int)(0.5*N)] + x[(int)(0.5*N) - 1]);
    418   } else {
    419     median = x[(int)(0.5*N)];
    420   }
    421 
    422   for (i = 0; i < N; i++ ) {
    423     x[i] = fabs(x[i] - median);
    424   }
    425 
    426   dsort(x,N);
    427 
    428   if (N % 2) {
    429     median = 0.5*(x[(int)(0.5*N)] + x[(int)(0.5*N) - 1]);
    430   } else {
    431     median = x[(int)(0.5*N)];
    432   }
    433 
    434   return(median);
    435 }
    436    
    437  
    438  
    439 int weighted_LS (double *T, double *X, double *WX, double *Y, double *WY, int Npts,
    440                  double **A, double **B, int VERBOSE) {
     366int weighted_LS_PM (double *T, double *X, double *WX, double *Y, double *WY, int Npts, double **A, double **B, int VERBOSE) {
    441367
    442368  int i,j;
    … …  
    497423  return TRUE;
    498424}
     425
     426double weight_cauchy (double x) {
     427  double r = x / 2.385;
     428  return (1.0 / (1.0 + SQ(r)));
     429}
     430
     431// dpsi = (d/dx) (x * weight(x))
     432double dpsi_cauchy (double x) {
     433  double r2 = SQ(x / 2.385);
     434  return ((1.0 - r2) / (SQ(1 + r2)));
     435}
     436
     437
     438// median absolute deviation
     439// MAD = median(abs(x - median(x)))
     440double MedianAbsDeviation(double *in, int N) {
     441  double *x;
     442  double median = 0.0;
     443  int i;
     444 
     445  ALLOCATE(x,double,N);
     446  for (i = 0; i < N; i++) {
     447    x[i] = in[i];
     448  }
     449
     450  dsort(x,N);
     451
     452  if (N % 2) {
     453    median = 0.5*(x[(int)(0.5*N)] + x[(int)(0.5*N) - 1]);
     454  } else {
     455    median = x[(int)(0.5*N)];
     456  }
     457
     458  for (i = 0; i < N; i++ ) {
     459    x[i] = fabs(x[i] - median);
     460  }
     461
     462  dsort(x,N);
     463
     464  if (N % 2) {
     465    median = 0.5*(x[(int)(0.5*N)] + x[(int)(0.5*N) - 1]);
     466  } else {
     467    median = x[(int)(0.5*N)];
     468  }
     469
     470  return(median);
     471}
     472 

  • branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/fiximage.c

    r34088 r39266  
    2626  gfits_free_matrix (&mask[0].matrix);
    2727  gfits_free_header (&mask[0].header);
    28   CreateBuffer (mask, Nx, Ny, -32, 0.0, 1.0);
     28  if (!CreateBuffer (mask, Nx, Ny, -32, 0.0, 1.0)) return FALSE;
    2929  strcpy (mask[0].file, "(empty)");
    3030  memset (mask[0].matrix.buffer, 0, Nx*Ny*sizeof(float));

  • branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/getcoords.c

    r37807 r39266  
    7373  gfits_free_matrix (&outC1[0].matrix);
    7474  gfits_free_header (&outC1[0].header);
    75   CreateBuffer (outC1, Nx, Ny, -32, 1.0, 0.0);
     75  if (!CreateBuffer (outC1, Nx, Ny, -32, 1.0, 0.0)) return FALSE;
    7676
    7777  gfits_free_matrix (&outC2[0].matrix);
    7878  gfits_free_header (&outC2[0].header);
    79   CreateBuffer (outC2, Nx, Ny, -32, 1.0, 0.0);
     79  if (!CreateBuffer (outC2, Nx, Ny, -32, 1.0, 0.0)) return FALSE;
    8080
    8181  float *valC1 = (float *) outC1[0].matrix.buffer;

  • branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/imfit.c

    r34088 r39266  
    141141
    142142    strcpy (out[0].file, "(empty)");
    143     CreateBuffer (out, 2*nx, 2*ny, -32, 0.0, 1.0);
     143    if (!CreateBuffer (out, 2*nx, 2*ny, -32, 0.0, 1.0)) return FALSE;
    144144
    145145    /* four panels: 1) raw image. 2) fit  3) raw - fit   4) ?? */

  • branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/init.c

    r39223 r39266  
    99int cneedles                PROTO((int, char **));
    1010int cplot                   PROTO((int, char **));
     11int crotation               PROTO((int, char **));
    1112int csystem                 PROTO((int, char **));
    1213int ctimes                  PROTO((int, char **));
    … …  
    7879  {1, "cneedles",    cneedles,     "plot vectors in sky coordinates"},
    7980  {1, "cplot",       cplot,        "plot vectors in sky coordinates"},
     81  {1, "crotation",   crotation,    "rotate in 3D"},
    8082  {1, "csystem",     csystem,      "convert between coordinate systems"},
    8183  {1, "ctimes",      ctimes,       "convert between time formats"},

  • branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/shimage.c

    r37807 r39266  
    3030  int Ny = src[0].header.Naxis[1];
    3131
    32   gfits_free_matrix (&FrIm[0].matrix); gfits_free_header (&FrIm[0].header); CreateBuffer (FrIm, Nx, Ny, -32, 1.0, 0.0);
    33   gfits_free_matrix (&FiIm[0].matrix); gfits_free_header (&FiIm[0].header); CreateBuffer (FiIm, Nx, Ny, -32, 1.0, 0.0);
     32  gfits_free_matrix (&FrIm[0].matrix); gfits_free_header (&FrIm[0].header); if (!CreateBuffer (FrIm, Nx, Ny, -32, 1.0, 0.0)) return FALSE;
     33  gfits_free_matrix (&FiIm[0].matrix); gfits_free_header (&FiIm[0].header); if (!CreateBuffer (FiIm, Nx, Ny, -32, 1.0, 0.0)) return FALSE;
    3434
    3535  // l=0 allocates space for a single value

  • branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/vshimage.c

    r37807 r39266  
    3232  int Ny = src[0].header.Naxis[1];
    3333
    34   gfits_free_matrix (&Rbim[0].matrix); gfits_free_header (&Rbim[0].header); CreateBuffer (Rbim, Nx, Ny, -32, 1.0, 0.0);
    35   gfits_free_matrix (&Reim[0].matrix); gfits_free_header (&Reim[0].header); CreateBuffer (Reim, Nx, Ny, -32, 1.0, 0.0);
    36   gfits_free_matrix (&Dbim[0].matrix); gfits_free_header (&Dbim[0].header); CreateBuffer (Dbim, Nx, Ny, -32, 1.0, 0.0);
    37   gfits_free_matrix (&Deim[0].matrix); gfits_free_header (&Deim[0].header); CreateBuffer (Deim, Nx, Ny, -32, 1.0, 0.0);
     34  gfits_free_matrix (&Rbim[0].matrix); gfits_free_header (&Rbim[0].header); if (!CreateBuffer (Rbim, Nx, Ny, -32, 1.0, 0.0)) return FALSE;
     35  gfits_free_matrix (&Reim[0].matrix); gfits_free_header (&Reim[0].header); if (!CreateBuffer (Reim, Nx, Ny, -32, 1.0, 0.0)) return FALSE;
     36  gfits_free_matrix (&Dbim[0].matrix); gfits_free_header (&Dbim[0].header); if (!CreateBuffer (Dbim, Nx, Ny, -32, 1.0, 0.0)) return FALSE;
     37  gfits_free_matrix (&Deim[0].matrix); gfits_free_header (&Deim[0].header); if (!CreateBuffer (Deim, Nx, Ny, -32, 1.0, 0.0)) return FALSE;
    3838
    3939  // l=0 allocates space for a single value
Note: See TracChangeset for help on using the changeset viewer.