Changeset 8926 for trunk/psModules/src/detrend/pmShutterCorrection.c

Timestamp:

Sep 24, 2006, 3:06:28 PM (20 years ago)

Author:

magnier

Message:

finished, tested pmShutterCorrection

File:

• trunk/psModules/src/detrend/pmShutterCorrection.c (modified) (9 diffs)

trunk/psModules/src/detrend/pmShutterCorrection.c

@@ -101 +101 @@
     // XXX we could examine the top 2 or 3 values and decide if we
     // extended exptime enough or median clip.
-    pars->scale = counts->data.F64[N-1];
+    pars->scale = counts->data.F32[N-1];
 
     // fit a line to the lowest three points and extrapolate to 0.0
-    tmpX = psVectorAlloc (2, PS_TYPE_F64);
-    tmpY = psVectorAlloc (2, PS_TYPE_F64);
-
-    tmpX->data.F64[0] = exptime->data.F64[0];
-    tmpX->data.F64[1] = exptime->data.F64[1];
-    tmpY->data.F64[0] = counts->data.F64[0];
-    tmpY->data.F64[1] = counts->data.F64[1];
+    tmpX = psVectorAlloc (2, PS_TYPE_F32);
+    tmpY = psVectorAlloc (2, PS_TYPE_F32);
+    tmpX->n = tmpY->n = 2;
+
+    tmpX->data.F32[0] = exptime->data.F32[0];
+    tmpX->data.F32[1] = exptime->data.F32[1];
+    tmpY->data.F32[0] = counts->data.F32[0];
+    tmpY->data.F32[1] = counts->data.F32[1];
 
     // fit a line and extrapolate the fit to 0.0
@@ -125 +126 @@
     float value = pars->scale * (1 + ratio) / 2.0;
 
-    int Nm = psVectorBracket (exptime, value, (ratio < 1.0));
-    int Np = (Nm == N - 1) ? Nm - 1 : Nm + 1;
-
-    tmpX->data.F64[0] = counts->data.F64[Nm];
-    tmpX->data.F64[1] = counts->data.F64[Np];
-    tmpY->data.F64[0] = exptime->data.F64[Nm];
-    tmpY->data.F64[1] = exptime->data.F64[Np];
+    int Np;
+    if (ratio < 1.0) {
+        Np = psVectorBracket (counts, value, true);
+    } else {
+        Np = psVectorBracketDescend (counts, value, true);
+    }
+    int Nm = (Np == 0) ? 1 : Np - 1;
+
+    tmpX->data.F32[0] = counts->data.F32[Nm];
+    tmpX->data.F32[1] = counts->data.F32[Np];
+    tmpY->data.F32[0] = exptime->data.F32[Nm];
+    tmpY->data.F32[1] = exptime->data.F32[Np];
 
     // fit a line and extrapolate the fit to counts = A (1 + dTk/dTo) : exptime = dTo
@@ -146 +152 @@
 
 // linear fit to the counts and exptime, given a value for offref
-pmShutterCorrPars *pmShutterCorrectionLinFit (psVector *exptime, psVector *counts, float offref)
+pmShutterCorrPars *pmShutterCorrectionLinFit (psVector *exptime, psVector *counts, psVector *cntError, float offref)
 {
 
@@ -152 +158 @@
     psVector *x = psVectorAlloc (exptime->n, PS_TYPE_F32);
     psVector *y = psVectorAlloc (exptime->n, PS_TYPE_F32);
+    x->n = y->n = exptime->n;
 
     for (int i = 0; i < exptime->n; i++) {
@@ -167 +174 @@
     line->coeff[1][1] = 0;
 
-    psVectorFitPolynomial2D (line, NULL, 0, counts, NULL, x, y);
+    psVectorFitPolynomial2D (line, NULL, 0, counts, cntError, x, y);
 
     pmShutterCorrPars *pars = pmShutterCorrParsAlloc ();
@@ -174 +181 @@
     pars->scale  = line->coeff[1][0];
     pars->offset = line->coeff[0][1] / line->coeff[1][0];
+
+    psFree (x);
+    psFree (y);
+    psFree (line);
 
     return (pars);
@@ -196 +207 @@
 
 // non-linear fit to the counts and exptime, given a guess for the three parameters
-pmShutterCorrPars *pmShutterCorrectionFullFit (psVector *exptime, psVector *counts, pmShutterCorrPars *guess)
+pmShutterCorrPars *pmShutterCorrectionFullFit (psVector *exptime, psVector *counts, psVector *cntError, pmShutterCorrPars *guess)
 {
 
@@ -221 +232 @@
     // construct the coordinate and value entries (y is counts)
     psArray *x = psArrayAlloc(exptime->n);
-    psVector *y = counts;
-    psVector *yErr = psVectorAlloc(exptime->n, PS_TYPE_F32);
-    x->n = yErr->n = exptime->n;
+    x->n = exptime->n;
 
     for (int i = 0; i < exptime->n; i++) {
         psVector *coord = psVectorAlloc(1, PS_TYPE_F32);
         coord->n = 1;
-
         coord->data.F32[0] = exptime->data.F32[i];
-
         x->data[i] = coord;
-        yErr->data.F32[i] = sqrt(y->data.F32[i]);
-    }
-
-    psMinimizeLMChi2(myMin, covar, params, constrain, x, y, yErr, pmShutterCorrectionModel);
+    }
+
+    psMinimizeLMChi2(myMin, covar, params, constrain, x, counts, cntError, pmShutterCorrectionModel);
 
     pmShutterCorrPars *pars = pmShutterCorrParsAlloc ();
@@ -243 +249 @@
     pars->offref = params->data.F32[2];
 
-    return (pars);
-}
+    psFree (myMin);
+    psFree (params);
+    psFree (x);
+
+    return (pars);
+}