Changeset 41387 for branches/eam_branches/relphot.20200519/src/liststats.c

Timestamp:

Jul 7, 2020, 9:30:05 AM (6 years ago)

Author:

eugene

Message:

updates to get tgroups working

File:

• branches/eam_branches/relphot.20200519/src/liststats.c (modified) (2 diffs)

branches/eam_branches/relphot.20200519/src/liststats.c

@@ -172 +172 @@
   for (i = ks; i < ke; i++) {
     M  = SQ (value[i] - Mo);
-    dM = SQ (dvalue[i]);
-    X2 += M / dM;
+    if (dvalue) {
+      dM = SQ (dvalue[i]);
+      X2 += M / dM;
+    }
     dS += M;
   }
-  X2 = X2 / (Nm - 1);
+  X2 = dvalue ? X2 / (Nm - 1) : NAN;
   dS = sqrt (dS / (Nm - 1));
 
@@ -187 +189 @@
   return (TRUE);
 }
+
+int liststats_fit1d (double *value, double *err, double *x, int Npts, StatType *stats, double *dk) {
+
+  double M = 0.0, MZ = 0.0, Z1 = 0.0, Z2 = 0.0, R = 0.0;
+
+  for (int i = 0; i < Npts; i++) {
+    double wt = 1.0 / (err[i]*err[i]);
+    double V = value[i];
+    double Z = x[i];
+    M  += V   * wt;
+    MZ += V*Z * wt;
+    Z1 += Z   * wt;
+    Z2 += Z*Z * wt;
+    R  += 1.0 * wt;
+  }
+
+  double DetInv = (R*Z2 - Z1*Z1);
+  double Det = (fabs(DetInv) < -1e6) ? 0.0 : 1.0 / DetInv; // do not allow unstable solutions
+
+  double zp = (M*Z2 - Z1*MZ)*Det;
+  *dk = (MZ*R - M*Z1)*Det;
+
+  stats->mean = zp;
+  stats->chisq = NAN;
+  return TRUE;
+}
+
 
 // These should probably be tunable: