IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 39228 for trunk


Ignore:
Timestamp:
Dec 3, 2015, 5:34:47 PM (11 years ago)
Author:
eugene
Message:

merge common code again (fitpm/fitpm_irls)

Location:
trunk/Ohana/src/opihi/cmd.astro
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/opihi/cmd.astro/astrom_ops.c

    r39227 r39228  
    4545
    4646  data->Npts = N;
    47   ALLOCATE (data->X, double, N);
    48   ALLOCATE (data->Y, double, N);
     47  ALLOCATE (data->X,  double, N);
     48  ALLOCATE (data->Y,  double, N);
    4949  ALLOCATE (data->dX, double, N);
    5050  ALLOCATE (data->dY, double, N);
    51   ALLOCATE (data->t, double, N);
     51  ALLOCATE (data->t,  double, N);
    5252  ALLOCATE (data->pX, double, N);
    5353  ALLOCATE (data->pY, double, N);

  • trunk/Ohana/src/opihi/cmd.astro/fitplx_irls.c

    r39227 r39228  
    340340 
    341341  // Solve OLS equation
    342   if (!Plx_weighted_LS(T,pR,pD,X,Wx,Y,Wy,Npts,
    343                        A,B,VERBOSE)) {
     342  if (!weighted_LS_PLX(T,pR,pD,X,Wx,Y,Wy,Npts,A,B,VERBOSE)) {
    344343    // Handle fail cases gracefully.
    345344    return(FALSE);
    … …  
    386385    // Assign W
    387386    for (i = 0; i < Npts; i++) {
    388       Wx[i] = Plx_weight_cauchy(rx[i] / dX[i]);
    389       Wy[i] = Plx_weight_cauchy(ry[i] / dY[i]);
     387      Wx[i] = weight_cauchy(rx[i] / dX[i]);
     388      Wy[i] = weight_cauchy(ry[i] / dY[i]);
    390389    }
    391390
    392391    // Solve
    393     if (!Plx_weighted_LS(T,pD,pR,X,Wx,Y,Wy,Npts,
    394                          A,B,VERBOSE)) {
     392    if (!weighted_LS_PLX(T,pD,pR,X,Wx,Y,Wy,Npts,A,B,VERBOSE)) {
    395393      // Handle fail case
    396394      return(FALSE);
    … …  
    410408
    411409    // Calculate sigma_hat from distribution of residual magnitudes
    412     sigma_hat = Plx_MAD(u,Npts) / 0.6745;
     410    sigma_hat = MedianAbsDeviation(u,Npts) / 0.6745;
    413411
    414412    // Check convergence
    … …  
    442440  lambda = 0.0;
    443441  for (i = 0; i < Npts; i++) {
    444     Wx[i] = Plx_weight_cauchy(rx[i] / dX[i]);
    445     Wy[i] = Plx_weight_cauchy(ry[i] / dY[i]);
    446 
    447     ax += Plx_dpsi_cauchy(rx[i] / dX[i]);
    448     ay += Plx_dpsi_cauchy(ry[i] / dY[i]);
     442    Wx[i] = weight_cauchy(rx[i] / dX[i]);
     443    Wy[i] = weight_cauchy(ry[i] / dY[i]);
     444
     445    ax += dpsi_cauchy(rx[i] / dX[i]);
     446    ay += dpsi_cauchy(ry[i] / dY[i]);
    449447
    450448    bx += SQ(Wx[i]);
    … …  
    526524}
    527525
    528 double Plx_weight_cauchy (double x) {
    529   double r = x / 2.385;
    530   return (1.0 / (1.0 + SQ(r)));
    531 }
    532 
    533 // dpsi = (d/dx) (x * weight(x))
    534 double Plx_dpsi_cauchy (double x) {
    535   double r2 = SQ(x / 2.385);
    536   return ((1.0 - r2) / (SQ(1 + r2)));
    537 }
    538 
    539 
    540 // median absolute deviation
    541 // MAD = median(abs(x - median(x)))
    542 double Plx_MAD(double *in, int N) {
    543   double *x;
    544   double median = 0.0;
    545   int i;
    546 
    547   ALLOCATE(x,double,N);
    548   for (i = 0; i < N; i++) {
    549     x[i] = in[i];
    550   }
    551 
    552   dsort(x,N);
    553 
    554   if (N % 2) {
    555     median = 0.5*(x[(int)(0.5*N)] + x[(int)(0.5*N) - 1]);
    556   } else {
    557     median = x[(int)(0.5*N)];
    558   }
    559 
    560   for (i = 0; i < N; i++ ) {
    561     x[i] = fabs(x[i] - median);
    562   }
    563 
    564   dsort(x,N);
    565 
    566   if (N % 2) {
    567     median = 0.5*(x[(int)(0.5*N)] + x[(int)(0.5*N) - 1]);
    568   } else {
    569     median = x[(int)(0.5*N)];
    570   }
    571 
    572   return(median);
    573 }
    574 
    575 int Plx_weighted_LS (double *T, double *pR, double *pD, double *X, double *WX, double *Y, double *WY, int Npts,
     526int weighted_LS_PLX (double *T, double *pR, double *pD, double *X, double *WX, double *Y, double *WY, int Npts,
    576527                     double **A, double **B, int VERBOSE) {
    577528

  • trunk/Ohana/src/opihi/cmd.astro/fitpm.c

    r37807 r39228  
    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;

  • trunk/Ohana/src/opihi/cmd.astro/fitpm_irls.c

    r39226 r39228  
    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
    … …  
    389364}
    390365
    391 
    392 double weight_cauchy (double x) {
    393   double r = x / 2.385;
    394   return (1.0 / (1.0 + SQ(r)));
    395 }
    396 
    397 // dpsi = (d/dx) (x * weight(x))
    398 double dpsi_cauchy (double x) {
    399   double r2 = SQ(x / 2.385);
    400   return ((1.0 - r2) / (SQ(1 + r2)));
    401 }
    402 
    403 
    404 // median absolute deviation
    405 // MAD = median(abs(x - median(x)))
    406 double MAD(double *in, int N) {
    407   double *x;
    408   double median = 0.0;
    409   int i;
    410  
    411   ALLOCATE(x,double,N);
    412   for (i = 0; i < N; i++) {
    413     x[i] = in[i];
    414   }
    415 
    416   dsort(x,N);
    417 
    418   if (N % 2) {
    419     median = 0.5*(x[(int)(0.5*N)] + x[(int)(0.5*N) - 1]);
    420   } else {
    421     median = x[(int)(0.5*N)];
    422   }
    423 
    424   for (i = 0; i < N; i++ ) {
    425     x[i] = fabs(x[i] - median);
    426   }
    427 
    428   dsort(x,N);
    429 
    430   if (N % 2) {
    431     median = 0.5*(x[(int)(0.5*N)] + x[(int)(0.5*N) - 1]);
    432   } else {
    433     median = x[(int)(0.5*N)];
    434   }
    435 
    436   return(median);
    437 }
    438    
    439  
    440  
    441 int weighted_LS (double *T, double *X, double *WX, double *Y, double *WY, int Npts,
    442                  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) {
    443367
    444368  int i,j;
    … …  
    499423  return TRUE;
    500424}
     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 
Note: See TracChangeset for help on using the changeset viewer.