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

Timestamp:

Sep 22, 2006, 11:00:53 AM (20 years ago)

Author:

magnier

Message:

nearly done with pmShutterCorrection

File:

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

trunk/psModules/src/detrend/pmShutterCorrection.c

@@ -51 +51 @@
 */
 
+#if HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdio.h>
+#include <math.h>
+#include <strings.h>
+#include <pslib.h>
+
+#include "pmShutterCorrection.h"
+
 static void pmShutterCorrParsFree (pmShutterCorrPars *pars)
 {
@@ -133 +144 @@
     return (pars);
 }
+
+// linear fit to the counts and exptime, given a value for offref
+pmShutterCorrPars *pmShutterCorrectionLinFit (psVector *exptime, psVector *counts, float offref)
+{
+
+    // this step is identical for all pixels: do it once and save?
+    psVector *x = psVectorAlloc (exptime->n, PS_TYPE_F32);
+    psVector *y = psVectorAlloc (exptime->n, PS_TYPE_F32);
+
+    for (int i = 0; i < exptime->n; i++) {
+        value = 1.0 / (exptime->data.F32[i] + offref);
+        x->data.F32[i] = exptime->data.F32[i] * value;
+        y->data.F32[i] = value;
+    }
+
+    psPolynomial2D *line = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, 1, 1);
+
+    // mask out the terms we will not fit
+    line->mask[0][0] = 1;
+    line->mask[1][1] = 1;
+    line->coeff[0][0] = 0;
+    line->coeff[1][1] = 0;
+
+    psVectorFitPolynomial2D (line, NULL, 0, counts, NULL, x, y);
+
+    pmShutterCorrPars *pars = pmShutterCorrParsAlloc ();
+
+    pars->offref = guess->offref;
+    pars->scale  = line->coeff[1][0];
+    pars->offset = line->coeff[0][1] / line->coeff[1][0];
+
+    return (pars);
+}
+
+static psF32 pmShutterCorrectionModel(psVector *deriv,
+                                      const psVector *params,
+                                      const psVector *x)
+{
+    psF32 A  = params->data.F32[0];
+    psF32 p  = x->data.F32[0] + params->data.F32[1];
+    psF32 q  = 1.0 / (x->data.F32[0] + params->data.F32[2]);
+    psF32 f  = A*p*q;
+
+    if (deriv != NULL) {
+        deriv->data.F32[0] = p*q;
+        deriv->data.F32[1] = A*q;
+        deriv->data.F32[2] = -A*p*q*q;
+    }
+    return(f);
+}
+
+// non-linear fit to the counts and exptime, given a guess for the three parameters
+pmShutterCorrPars *pmShutterCorrectionFullFit (psVector *exptime, psVector *counts, pmShutterCorrPars *guess)
+{
+
+    psMinimization *myMin = psMinimizationAlloc(15, 0.1);
+
+    psVector *params = psVectorAlloc (3, PS_TYPE_F32);
+    params->data.F32[0] = guess->scale;
+    params->data.F32[1] = guess->offset;
+    params->data.F32[2] = guess->offref;
+    params->n = 3;
+
+    // XXX for the moment, don't set any constraints
+    // psMinConstrain *constrain = psMinConstrainAlloc();
+    // constrain->paramDelta = psVectorAlloc (3, PS_TYPE_F32);
+    // constrain->paramMin   = psVectorAlloc (3, PS_TYPE_F32);
+    // constrain->paramMax   = psVectorAlloc (3, PS_TYPE_F32);
+    // constrain->paramMask  = NULL;
+    psMinConstrain *constrain = NULL;
+
+    // XXX ignore covariance matrix for the moment
+    // psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F64);
+    psImage *covar = NULL;
+
+    // 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;
+
+    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]);
+    }
+
+    fitStatus = psMinimizeLMChi2(myMin, covar, params, constrain, x, y, yErr, pmShutterCorrectionModel);
+
+    pmShutterCorrPars *pars = pmShutterCorrParsAlloc ();
+
+    pars->scale  = params->data.F32[0];
+    pars->offset = params->data.F32[1];
+    pars->offref = params->data.F32[2];
+
+    return (pars);
+}