Changeset 39266 for branches/eam_branches/ipp-20151113/Ohana/src/opihi/cmd.astro/fitpm_irls.c

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:

• . (modified) (1 prop)
• Ohana/src/opihi/cmd.astro/fitpm_irls.c (modified) (9 diffs)

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

@@ -1 +1 @@
 # include "astro.h"
-# define J2000 51544.5       /* Modified Julian date at standard epoch J2000 */
-
-# define ESCAPE(MSG,...) {                      \
-    gprint (GP_ERR, MSG, __VA_ARGS__);          \
-    return FALSE; }
-
-typedef struct {
-  double Ro, dRo;
-  double Do, dDo;
-
-  double uR, duR;
-  double uD, duD;
-  
-  double chisq;
-  int Nfit;
-} PMFit_IRLS;
-
-int FitPMonly_IRLS (PMFit_IRLS *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, int VERBOSE);
-int IRLS_converged (PMFit_IRLS *fit);
-int weighted_LS (double *T, double *X, double *WX, double *Y, double *WY, int Npts,
-                 double **A, double **B, int VERBOSE);
-double weight_cauchy (double x);
-double dpsi_cauchy (double x);
-double MAD(double *in, int N);
-      
 
 int fitpm_irls (int argc, char **argv) {
@@ -129 +104 @@
   }
 
-  PMFit_IRLS fit;
+  PlxFit fit;
   if (!FitPMonly_IRLS (&fit, X, dX, Y, dY, t, n, VERBOSE)) {
     return FALSE;
@@ -165 +140 @@
 
 /* do we want an init function which does the alloc and a clear function to free? */
-int FitPMonly_IRLS (PMFit_IRLS *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, int VERBOSE) {
+int FitPMonly_IRLS (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, int VERBOSE) {
 
   int i,j;
@@ -222 +197 @@
   
   // Solve OLS equation  
-  if (!weighted_LS(T,X,Wx,Y,Wy,Npts,
+  if (!weighted_LS_PM(T,X,Wx,Y,Wy,Npts,
                    A,B,VERBOSE)) {
     // Handle fail case
@@ -267 +242 @@
 
     // Solve
-    if (!weighted_LS(T,X,Wx,Y,Wy,Npts,
+    if (!weighted_LS_PM(T,X,Wx,Y,Wy,Npts,
                      A,B,VERBOSE)) {
       // Handle fail case
@@ -284 +259 @@
       u[i] = sqrt(SQ(rx[i] / dX[i]) + SQ(ry[i] / dY[i]));
     }
-    sigma_hat = MAD(u,Npts) / 0.6745;
+    sigma_hat = MedianAbsDeviation(u,Npts) / 0.6745;
     
     // Check convergence
@@ -310 +285 @@
   double Sum_Wx, Sum_Wy;
   
+  Sum_Wx = 0.0;
+  Sum_Wy = 0.0;
   ax = 0.0; ay = 0.0;
   bx = 0.0; by = 0.0;
@@ -387 +364 @@
 }
 
-
-double weight_cauchy (double x) {
-  double r = x / 2.385;
-  return (1.0 / (1.0 + SQ(r)));
-}
-
-// dpsi = (d/dx) (x * weight(x))
-double dpsi_cauchy (double x) {
-  double r2 = SQ(x / 2.385);
-  return ((1.0 - r2) / (SQ(1 + r2)));
-}
-
-
-// median absolute deviation
-// MAD = median(abs(x - median(x)))
-double MAD(double *in, int N) {
-  double *x;
-  double median = 0.0;
-  int i;
-  
-  ALLOCATE(x,double,N);
-  for (i = 0; i < N; i++) {
-    x[i] = in[i];
-  }
-
-  dsort(x,N);
-
-  if (N % 2) {
-    median = 0.5*(x[(int)(0.5*N)] + x[(int)(0.5*N) - 1]);
-  } else {
-    median = x[(int)(0.5*N)];
-  }
-
-  for (i = 0; i < N; i++ ) {
-    x[i] = fabs(x[i] - median);
-  }
-
-  dsort(x,N);
-
-  if (N % 2) {
-    median = 0.5*(x[(int)(0.5*N)] + x[(int)(0.5*N) - 1]);
-  } else {
-    median = x[(int)(0.5*N)];
-  }
-
-  return(median);
-}
-    
-  
-  
-int weighted_LS (double *T, double *X, double *WX, double *Y, double *WY, int Npts,
-                 double **A, double **B, int VERBOSE) {
+int weighted_LS_PM (double *T, double *X, double *WX, double *Y, double *WY, int Npts, double **A, double **B, int VERBOSE) {
 
   int i,j;
@@ -497 +423 @@
   return TRUE;
 }
+
+double weight_cauchy (double x) {
+  double r = x / 2.385;
+  return (1.0 / (1.0 + SQ(r)));
+}
+
+// dpsi = (d/dx) (x * weight(x))
+double dpsi_cauchy (double x) {
+  double r2 = SQ(x / 2.385);
+  return ((1.0 - r2) / (SQ(1 + r2)));
+}
+
+
+// median absolute deviation
+// MAD = median(abs(x - median(x)))
+double MedianAbsDeviation(double *in, int N) {
+  double *x;
+  double median = 0.0;
+  int i;
+  
+  ALLOCATE(x,double,N);
+  for (i = 0; i < N; i++) {
+    x[i] = in[i];
+  }
+
+  dsort(x,N);
+
+  if (N % 2) {
+    median = 0.5*(x[(int)(0.5*N)] + x[(int)(0.5*N) - 1]);
+  } else {
+    median = x[(int)(0.5*N)];
+  }
+
+  for (i = 0; i < N; i++ ) {
+    x[i] = fabs(x[i] - median);
+  }
+
+  dsort(x,N);
+
+  if (N % 2) {
+    median = 0.5*(x[(int)(0.5*N)] + x[(int)(0.5*N) - 1]);
+  } else {
+    median = x[(int)(0.5*N)];
+  }
+
+  return(median);
+}
+