Changeset eca05c8 in sasview for pr_inversion/c_extensions/invertor.c

Timestamp:

Apr 30, 2008 4:29:36 PM (16 years ago)

Author:

Mathieu Doucet <doucetm@…>

Branches:

master, ESS_GUI, ESS_GUI_Docs, ESS_GUI_batch_fitting, ESS_GUI_bumps_abstraction, ESS_GUI_iss1116, ESS_GUI_iss879, ESS_GUI_iss959, ESS_GUI_opencl, ESS_GUI_ordering, ESS_GUI_sync_sascalc, costrafo411, magnetic_scatt, release-4.1.1, release-4.1.2, release-4.2.2, release_4.0.1, ticket-1009, ticket-1094-headless, ticket-1242-2d-resolution, ticket-1243, ticket-1249, ticket885, unittest-saveload

Children:

34ae302

Parents:

47f695c9

Message:

Added P(r) fit test

File:

• pr_inversion/c_extensions/invertor.c (modified) (3 diffs)

pr_inversion/c_extensions/invertor.c

@@ -40 +40 @@
 double ortho(double d_max, int n, double r) {
         return 2.0*r*sin(pi*n*r/d_max);
-        
 }
 
@@ -76 +75 @@
  */
 double pr(double *pars, double d_max, int n_c, double r) {
-    double sum = 0.0; 
+    double sum = 0.0;  
         int i;
     for (i=0; i<n_c; i++) {
@@ -84 +83 @@
 }
 
+/**
+ * P(r) calculated from the expansion, with errors
+ */
+void pr_err(double *pars, double *err, double d_max, int n_c, 
+                double r, double *pr_value, double *pr_value_err) {
+    double sum = 0.0; 
+    double sum_err = 0.0;
+    double func_value;
+        int i;
+    for (i=0; i<n_c; i++) {
+        func_value = ortho(d_max, i+1, r);
+        sum += pars[i] * func_value;
+        sum_err += err[i]*err[i]*func_value*func_value;
+    }
+    *pr_value = sum;
+    if (sum_err>0) {
+        *pr_value_err = sqrt(sum_err);
+    } else {
+        *pr_value_err = sum;
+    }
+} 
 
+/**
+ * dP(r)/dr calculated from the expansion.
+ */
+double dprdr(double *pars, double d_max, int n_c, double r) {
+    double sum = 0.0; 
+        int i;
+    for (i=0; i<n_c; i++) {
+        sum += pars[i] * 2.0*(sin(pi*(i+1)*r/d_max) + pi*(i+1)*r/d_max * cos(pi*(i+1)*r/d_max));
+    }
+    return sum;
+}
+
+/**
+ * regularization term calculated from the expansion.
+ */
+double reg_term(double *pars, double d_max, int n_c) {
+    double sum = 0.0; 
+    double r;
+    double deriv;
+        int i;
+    for (i=0; i<25; i++) {
+        r = d_max/25.0*i;
+        deriv = dprdr(pars, d_max, n_c, r);
+        sum += deriv*deriv;
+    }
+    return sum/25.0*d_max;
+}
+
+