Changeset 39086 for trunk

Timestamp:

Nov 9, 2015, 2:46:18 PM (11 years ago)

Author:

watersc1

Message:

Working version of IRLS proper motion fitting.

File:

• trunk/Ohana/src/opihi/cmd.astro/fitpm_irls.c (modified) (6 diffs)

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

@@ -19 +19 @@
 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) {
@@ -161 +167 @@
 int FitPMonly_IRLS (PMFit_IRLS *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, int VERBOSE) {
 
-  int i;
+  int i,j;
 
   static double **A, **B;
@@ -171 +177 @@
   double *Beta, *Beta_prev;
   
-  double sigma_robust, sigma_ols, sigma_final, sigma_hat;
+  double sigma_ols, sigma_hat;
   double *Wx, *Wy;
   double *rx, *ry;
-  double *ux, *uy;
+  //  double *ux, *uy;
+  double *u;
   int dof = 2 * Npts - 4;
   int p   = 4;
   int n   = 2 * Npts;
   double tolerance;
-  int convergence;
+  int converged;
   int iterations;
   
@@ -324 +331 @@
   by /= 1.0 * (Npts - p);
   
-  sigma_robust_x = lamba * sqrt(bx) * sigma_hat * 2.385 / ax;
-  sigma_robust_y = lamba * sqrt(by) * sigma_hat * 2.385 / ay;
+  sigma_robust_x = lambda * sqrt(bx) * sigma_hat * 2.385 / ax;
+  sigma_robust_y = lambda * sqrt(by) * sigma_hat * 2.385 / ay;
 
   // This is actually sigma^2, as that's the factor in the covariance (dumouchel 4.1)
-  sigma_final_x  = fmax(SQ(sigma_robust_x), (n * SQ(sigma_robust_x) + SQ(p * sigma_OLS)) / (n + SQ(p)));
-  sigma_final_y  = fmax(SQ(sigma_robust_y), (n * SQ(sigma_robust_y) + SQ(p * sigma_OLS)) / (n + SQ(p)));
+  sigma_final_x  = MAX(SQ(sigma_robust_x), (n * SQ(sigma_robust_x) + SQ(p * sigma_ols)) / (n + SQ(p)));
+  sigma_final_y  = MAX(SQ(sigma_robust_y), (n * SQ(sigma_robust_y) + SQ(p * sigma_ols)) / (n + SQ(p)));
 
   for (i = 0; i < 4; i++) {
@@ -337 +344 @@
         Cov[i][j] *= sigma_final_x;
       }
-      elsif ((i > 1)&&(j > 1)) { // Lower portion
+      else if ((i > 1)&&(j > 1)) { // Lower portion
         Cov[i][j] *= sigma_final_y;
       }
@@ -475 +482 @@
     if (VERBOSE) fprintf (stderr, "error in fit\n");
     if (VERBOSE == 2) {
-      int j;
       for (i = 0; i < 4; i++) {
         for (j = 0; j < 4; j++) {