Changeset d0462cf in sasmodels

Timestamp:

Nov 2, 2017 8:50:27 PM (7 years ago)

Author:

GitHub <noreply@…>

Branches:

master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests

Children:

c11d09f

Parents:

f475f90 (diff), 17db833 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.

git-author:

Paul Kienzle <pkienzle@…> (11/02/17 20:50:27)

git-committer:

GitHub <noreply@…> (11/02/17 20:50:27)

Message:

Merge branch 'ticket-776-orientation' into ticket-786

Location:

sasmodels

Files:

• details.py (modified) (2 diffs)
• kernel_iq.c (modified) (12 diffs)
• modelinfo.py (modified) (1 diff)
• models/core_shell_bicelle_elliptical_belt_rough.c (added)
• models/core_shell_bicelle_elliptical_belt_rough.py (added)
• models/img/core_shell_bicelle_belt_geometry.png (added)
• models/img/core_shell_bicelle_belt_parameters.png (added)
• product.py (modified) (2 diffs)
• sasview_model.py (modified) (7 diffs)
• models/core_shell_parallelepiped.c (modified) (4 diffs)

sasmodels/details.py

@@ -193 +193 @@
     #print("off:",offset)
 
-    # Check that we arn't using too many polydispersity loops
+    # Check that we aren't using too many polydispersity loops
     num_active = np.sum(length > 1)
     max_pd = model_info.parameters.max_pd
@@ -318 +318 @@
     parameter set in the vector.
     """
-    value, dispersity, weight = zip(*mesh)
+    _, dispersity, weight = zip(*mesh)
     #weight = [w if len(w)>0 else [1.] for w in weight]
     weight = np.vstack([v.flatten() for v in meshgrid(*weight)])

sasmodels/kernel_iq.c

@@ -13 +13 @@
 // NOTE: the following macros are defined in generate.py:
 //
-//  MAX_PD : the maximum number of dispersity loops allowed
+//  MAX_PD : the maximum number of dispersity loops allowed for this model,
+//      which will be at most modelinfo.MAX_PD.
 //  NUM_PARS : the number of parameters in the parameter table
 //  NUM_VALUES : the number of values to skip at the start of the
@@ -282 +283 @@
 #endif
 
-  // Storage for the current parameter values.  These will be updated as we
-  // walk the dispersity mesh.
+  // ** Fill in the local values table **
+  // Storage for the current parameter values.
+  // These will be updated as we walk the dispersity mesh.
   ParameterBlock local_values;
-
+  //   values[0] is scale
+  //   values[1] is background
+  #ifdef USE_OPENMP
+  #pragma omp parallel for
+  #endif
+  for (int i=0; i < NUM_PARS; i++) {
+    local_values.vector[i] = values[2+i];
+    //if (q_index==0) printf("p%d = %g\n",i, local_values.vector[i]);
+  }
+  //if (q_index==0) printf("NUM_VALUES:%d  NUM_PARS:%d  MAX_PD:%d\n", NUM_VALUES, NUM_PARS, MAX_PD);
+  //if (q_index==0) printf("start:%d  stop:%d\n", pd_start, pd_stop);
+
+  // ** Precompute magnatism values **
 #if defined(MAGNETIC) && NUM_MAGNETIC>0
   // Location of the sld parameters in the parameter vector.
@@ -302 +316 @@
 #endif // MAGNETIC
 
-  // Fill in the initial variables
-  //   values[0] is scale
-  //   values[1] is background
-  #ifdef USE_OPENMP
-  #pragma omp parallel for
-  #endif
-  for (int i=0; i < NUM_PARS; i++) {
-    local_values.vector[i] = values[2+i];
-//if (q_index==0) printf("p%d = %g\n",i, local_values.vector[i]);
-  }
-//if (q_index==0) printf("NUM_VALUES:%d  NUM_PARS:%d  MAX_PD:%d\n", NUM_VALUES, NUM_PARS, MAX_PD);
-//if (q_index==0) printf("start:%d  stop:%d\n", pd_start, pd_stop);
-
+  // ** Fill in the initial results **
   // If pd_start is zero that means that we are starting a new calculation,
   // and must initialize the result to zero.  Otherwise, we are restarting
@@ -336 +338 @@
       for (int q_index=0; q_index < nq; q_index++) result[q_index] = 0.0;
     }
-//if (q_index==0) printf("start %d %g %g\n", pd_start, pd_norm, result[0]);
+    //if (q_index==0) printf("start %d %g %g\n", pd_start, pd_norm, result[0]);
 #endif // !USE_OPENCL
 
@@ -390 +392 @@
       PD_OPEN(0,1)
 
-      ... main loop body ...
+      // ... main loop body ...
+      APPLY_PROJECTION    // convert jitter values to spherical coords
+      BUILD_ROTATION      // construct the rotation matrix qxy => qabc
+      for each q
+          FETCH_Q         // set qx,qy from the q input vector
+          APPLY_ROTATION  // convert qx,qy to qa,qb,qc
+          CALL_KERNEL     // scattering = Iqxy(qa, qb, qc, p1, p2, ...)
 
       ++step;  // increment counter representing position in dispersity mesh
@@ -412 +420 @@
 */
 
-// Define looping variables
-#define PD_INIT(_LOOP) \
-  const int n##_LOOP = details->pd_length[_LOOP]; \
-  const int p##_LOOP = details->pd_par[_LOOP]; \
-  global const double *v##_LOOP = pd_value + details->pd_offset[_LOOP]; \
-  global const double *w##_LOOP = pd_weight + details->pd_offset[_LOOP]; \
-  int i##_LOOP = (pd_start/details->pd_stride[_LOOP])%n##_LOOP;
-
-// Jump into the middle of the dispersity loop
-#define PD_OPEN(_LOOP,_OUTER) \
-  while (i##_LOOP < n##_LOOP) { \
-    local_values.vector[p##_LOOP] = v##_LOOP[i##_LOOP]; \
-    const double weight##_LOOP = w##_LOOP[i##_LOOP] * weight##_OUTER;
-
-// create the variable "weight#=1.0" where # is the outermost level+1 (=MAX_PD).
-#define _PD_OUTERMOST_WEIGHT(_n) const double weight##_n = 1.0;
-#define PD_OUTERMOST_WEIGHT(_n) _PD_OUTERMOST_WEIGHT(_n)
-
-// Close out the loop
-#define PD_CLOSE(_LOOP) \
-    if (step >= pd_stop) break; \
-    ++i##_LOOP; \
-  } \
-  i##_LOOP = 0;
-
-// The main loop body is essentially:
-//
-//    BUILD_ROTATION      // construct the rotation matrix qxy => qabc
-//    for each q
-//        FETCH_Q         // set qx,qy from the q input vector
-//        APPLY_ROTATION  // convert qx,qy to qa,qb,qc
-//        CALL_KERNEL     // scattering = Iqxy(qa, qb, qc, p1, p2, ...)
-//
+
+// ** prepare inner loops **
+
 // Depending on the shape type (radial, axial, triaxial), the variables
-// and calling parameters will be slightly different.  These macros
-// capture the differences in one spot so the rest of the code is easier
-// to read.
+// and calling parameters in the loop body will be slightly different.
+// Macros capture the differences in one spot so the rest of the code
+// is easier to read. The code below both declares variables for the
+// inner loop and defines the macros that use them.
 
 #if defined(CALL_IQ)
@@ -486 +465 @@
   // psi and dpsi are only for IQ_ABC, so they are processed here.
   const double psi = values[details->theta_par+4];
+  local_values.table.psi = 0.;
   #define BUILD_ROTATION() qabc_rotation(&rotation, theta, phi, psi, dtheta, dphi, local_values.table.psi)
   #define APPLY_ROTATION() qabc_apply(rotation, qx, qy, &qa, &qb, &qc)
@@ -501 +481 @@
   const double theta = values[details->theta_par+2];
   const double phi = values[details->theta_par+3];
+  // Make sure jitter angle defaults to zero if there is no jitter distribution
+  local_values.table.theta = 0.;
+  local_values.table.phi = 0.;
   // The "jitter" angles (dtheta, dphi, dpsi) are stored with the
   // dispersity values and copied to the local parameter table as
@@ -523 +506 @@
 #endif // !spherosymmetric projection
 
+// ** define looping macros **
+
+// Define looping variables
+#define PD_INIT(_LOOP) \
+  const int n##_LOOP = details->pd_length[_LOOP]; \
+  const int p##_LOOP = details->pd_par[_LOOP]; \
+  global const double *v##_LOOP = pd_value + details->pd_offset[_LOOP]; \
+  global const double *w##_LOOP = pd_weight + details->pd_offset[_LOOP]; \
+  int i##_LOOP = (pd_start/details->pd_stride[_LOOP])%n##_LOOP;
+
+// Jump into the middle of the dispersity loop
+#define PD_OPEN(_LOOP,_OUTER) \
+  while (i##_LOOP < n##_LOOP) { \
+    local_values.vector[p##_LOOP] = v##_LOOP[i##_LOOP]; \
+    const double weight##_LOOP = w##_LOOP[i##_LOOP] * weight##_OUTER;
+
+// create the variable "weight#=1.0" where # is the outermost level+1 (=MAX_PD).
+#define _PD_OUTERMOST_WEIGHT(_n) const double weight##_n = 1.0;
+#define PD_OUTERMOST_WEIGHT(_n) _PD_OUTERMOST_WEIGHT(_n)
+
+// Close out the loop
+#define PD_CLOSE(_LOOP) \
+    if (step >= pd_stop) break; \
+    ++i##_LOOP; \
+  } \
+  i##_LOOP = 0;
 
 // ====== construct the loops =======
@@ -537 +546 @@
 int step = pd_start;
 
+// *** define loops for each of 0, 1, 2, ..., modelinfo.MAX_PD-1 ***
+
 // define looping variables
 #if MAX_PD>4
@@ -571 +582 @@
   PD_OPEN(0,1)
 #endif
-
 
 //if (q_index==0) {printf("step:%d of %d, pars:",step,pd_stop); for (int i=0; i < NUM_PARS; i++) printf("p%d=%g ",i, local_values.vector[i]); printf("\n");}
@@ -672 +682 @@
 #endif // !USE_OPENCL
 
-// clear the macros in preparation for the next kernel
+// ** clear the macros in preparation for the next kernel **
 #undef PD_INIT
 #undef PD_OPEN

sasmodels/modelinfo.py

@@ -30 +30 @@
     TestCondition = Tuple[ParameterSetUser, TestInput, TestValue]
 
-MAX_PD = 4 #: Maximum number of simultaneously polydisperse parameters
+# If MAX_PD changes, need to change the loop macros in kernel_iq.c
+MAX_PD = 5 #: Maximum number of simultaneously polydisperse parameters
 
 # assumptions about common parameters exist throughout the code, such as:

sasmodels/product.py

@@ -101 +101 @@
     model_info.composition = ('product', [p_info, s_info])
     model_info.control = p_info.control
+    model_info.hidden = p_info.hidden
+    if getattr(p_info, 'profile', None) is not None:
+        profile_pars = set(p.id for p in p_info.parameters.kernel_parameters)
+        def profile(**kwargs):
+            # extract the profile args
+            kwargs = dict((k, v) for k, v in kwargs.items() if k in profile_pars)
+            return p_info.profile(**kwargs)
+    else:
+        profile = None
+    model_info.profile = profile
+    model_info.profile_axes = p_info.profile_axes
+
     # TODO: delegate random to p_info, s_info
     #model_info.random = lambda: {}
@@ -267 +279 @@
     weight = values[nvalues + call_details.num_weights: nvalues + 2*call_details.num_weights]
     npars = model_info.parameters.npars
-    pairs = [(value[offset:offset+length], weight[offset:offset+length])
+    # Note: changed from pairs ([v], [w]) to triples (p, [v], [w]), but the
+    # dispersion mesh code doesn't actually care about the nominal parameter
+    # value, p, so set it to None.
+    pairs = [(None, value[offset:offset+length], weight[offset:offset+length])
              for p, offset, length
              in zip(model_info.parameters.call_parameters[2:2+npars],

sasmodels/sasview_model.py

@@ -205 +205 @@
                                            structure_factor._model_info)
     ConstructedModel = make_model_from_info(model_info)
-    return ConstructedModel(form_factor.multiplicity)    
+    return ConstructedModel(form_factor.multiplicity)
 
 
@@ -673 +673 @@
         #call_details.show()
         #print("pairs", pairs)
+        #for k, p in enumerate(self._model_info.parameters.call_parameters):
+        #    print(k, p.name, *pairs[k])
         #print("params", self.params)
         #print("values", values)
@@ -678 +680 @@
         result = calculator(call_details, values, cutoff=self.cutoff,
                             magnetic=is_magnetic)
+        #print("result", result)
         self._intermediate_results = getattr(calculator, 'results', None)
         calculator.release()
@@ -761 +764 @@
                 return self.multiplicity, [self.multiplicity], [1.0]
             else:
-                # For hidden parameters use the default value.
+                # For hidden parameters use default values.  This sets
+                # scale=1 and background=0 for structure factors
                 default = self._model_info.parameters.defaults.get(par.name, np.NaN)
-                return [default], [1.0]
+                return default, [default], [1.0]
         elif par.polydisperse:
             value = self.params[par.name]
@@ -776 +780 @@
         else:
             value = self.params[par.name]
-            return value, [value if par.relative_pd else 0.0], [1.0]
+            return value, [value], [1.0]
 
 def test_cylinder():
@@ -794 +798 @@
     Model = _make_standard_model('hardsphere')
     model = Model()
+    value2d = model.evalDistribution([0.1, 0.1])
+    value1d = model.evalDistribution(np.array([0.1*np.sqrt(2)]))
+    #print("hardsphere", value1d, value2d)
+    if np.isnan(value1d) or np.isnan(value2d):
+        raise ValueError("hardsphere returns nan")
+
+def test_product():
+    # type: () -> float
+    """
+    Test that 2-D hardsphere model runs and doesn't produce NaN.
+    """
+    S = _make_standard_model('hayter_msa')()
+    P = _make_standard_model('cylinder')()
+    model = MultiplicationModel(P, S)
     value = model.evalDistribution([0.1, 0.1])
     if np.isnan(value):
-        raise ValueError("hardsphere returns null")
+        raise ValueError("cylinder*hatyer_msa returns null")
 
 def test_rpa():
@@ -850 +868 @@
 if __name__ == "__main__":
     print("cylinder(0.1,0.1)=%g"%test_cylinder())
+    #test_product()
+    #test_structure_factor()
+    #print("rpa:", test_rpa())
     #test_empty_distribution()

sasmodels/models/core_shell_parallelepiped.c

@@ -3 +3 @@
     double thick_rim_a, double thick_rim_b, double thick_rim_c)
 {
-    //return length_a * length_b * length_c;
     return length_a * length_b * length_c +
            2.0 * thick_rim_a * length_b * length_c +
@@ -24 +23 @@
     double thick_rim_c)
 {
-    // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c_scaled
+    // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c
     // Did not understand the code completely, it should be rechecked (Miguel Gonzalez)
     //Code is rewritten,the code is compliant with Diva Singhs thesis now (Dirk Honecker)
@@ -30 +29 @@
     const double mu = 0.5 * q * length_b;
 
-    //calculate volume before rescaling (in original code, but not used)
-    //double vol = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c);
-    //double vol = length_a * length_b * length_c +
-    //       2.0 * thick_rim_a * length_b * length_c +
-    //       2.0 * thick_rim_b * length_a * length_c +
-    //       2.0 * thick_rim_c * length_a * length_b;
+    // Scale sides by B
+    const double a_over_b = length_a / length_b;
+    const double c_over_b = length_c / length_b;
 
-    // Scale sides by B
-    const double a_scaled = length_a / length_b;
-    const double c_scaled = length_c / length_b;
-
-    double ta = a_scaled + 2.0*thick_rim_a/length_b; // incorrect ta = (a_scaled + 2.0*thick_rim_a)/length_b;
-    double tb = 1+ 2.0*thick_rim_b/length_b; // incorrect tb = (a_scaled + 2.0*thick_rim_b)/length_b;
-    double tc = c_scaled + 2.0*thick_rim_c/length_b; //not present
+    double tA_over_b = a_over_b + 2.0*thick_rim_a/length_b;
+    double tB_over_b = 1+ 2.0*thick_rim_b/length_b;
+    double tC_over_b = c_over_b + 2.0*thick_rim_c/length_b;
 
     double Vin = length_a * length_b * length_c;
-    //double Vot = (length_a * length_b * length_c +
-    //            2.0 * thick_rim_a * length_b * length_c +
-    //            2.0 * length_a * thick_rim_b * length_c +
-    //            2.0 * length_a * length_b * thick_rim_c);
-    double V1 = (2.0 * thick_rim_a * length_b * length_c);    // incorrect V1 (aa*bb*cc+2*ta*bb*cc)
-    double V2 = (2.0 * length_a * thick_rim_b * length_c);    // incorrect V2(aa*bb*cc+2*aa*tb*cc)
-    double V3 = (2.0 * length_a * length_b * thick_rim_c);    //not present
+    double VtA = (2.0 * thick_rim_a * length_b * length_c);
+    double VtB = (2.0 * length_a * thick_rim_b * length_c);
+    double VtC = (2.0 * length_a * length_b * thick_rim_c);
 
     // Scale factors (note that drC is not used later)
-    const double drho0 = (core_sld-solvent_sld);
-    const double drhoA = (arim_sld-solvent_sld);
-    const double drhoB = (brim_sld-solvent_sld);
-    const double drhoC = (crim_sld-solvent_sld);  // incorrect const double drC_Vot = (crim_sld-solvent_sld)*Vot;
-
+    const double dr0 = (core_sld-solvent_sld);
+    const double drA = (arim_sld-solvent_sld);
+    const double drB = (brim_sld-solvent_sld);
+    const double drC = (crim_sld-solvent_sld);
 
     // Precompute scale factors for combining cross terms from the shape
-    const double scale23 = drhoA*V1;
-    const double scale14 = drhoB*V2;
-    const double scale24 = drhoC*V3;
-    const double scale11 = drho0*Vin;
-    const double scale12 = drho0*Vin - scale23 - scale14 - scale24;
+    const double dr0_Vin = dr0*Vin;
+    const double drA_VtA = drA*VtA;
+    const double drB_VtB = drB*VtB;
+    const double drC_VtC = drC*VtC;
+    const double drV_delta = dr0_Vin - drA_VtA - drB_VtB - drC_VtC;
+
+    /*  *************** algorithm description ******************
+
+    // Rewrite f as x*siC + y*siCt to move the siC/siCt calculation out
+    // of the inner loop.  That is:
+
+    f = (di-da-db-dc) sa sb sc + da sa' sb sc + db sa sb' sc + dc sa sb sc'
+      =  [ (di-da-db-dc) sa sb + da sa' sb + db sa sb' ] sc  + [dc sa sb] sc'
+      = x sc + y sc'
+
+    // where:
+    di = delta rho_core V_core
+    da = delta rho_rimA V_rimA
+    db = delta rho_rimB V_rimB
+    dc = delta rho_rimC V_rimC
+    sa = j_0 (q_a a/2)    // siA the code
+    sb = j_0 (q_b b/2)
+    sc = j_0 (q_c c/2)
+    sa' = j_0(q_a a_rim/2)  // siAt the code
+    sb' = j_0(q_b b_rim/2)
+    sc' = j_0(q_c c_rim/2)
+
+    // qa, qb, and qc are generated using polar coordinates, with the
+    // outer loop integrating over [0,1] after the u-substitution
+    //    sigma = cos(theta), sqrt(1-sigma^2) = sin(theta)
+    // and inner loop integrating over [0, pi/2] as
+    //    uu = phi
+
+    ************************************************************  */
 
     // outer integral (with gauss points), integration limits = 0, 1
-    double outer_total = 0; //initialize integral
-
+    double outer_sum = 0; //initialize integral
     for( int i=0; i<76; i++) {
         double sigma = 0.5 * ( Gauss76Z[i] + 1.0 );
         double mu_proj = mu * sqrt(1.0-sigma*sigma);
 
-        // inner integral (with gauss points), integration limits = 0, 1
-        double inner_total = 0.0;
-        double inner_total_crim = 0.0;
+        // inner integral (with gauss points), integration limits = 0, pi/2
+        const double siC = sas_sinx_x(mu * sigma * c_over_b);
+        const double siCt = sas_sinx_x(mu * sigma * tC_over_b);
+        double inner_sum = 0.0;
         for(int j=0; j<76; j++) {
             const double uu = 0.5 * ( Gauss76Z[j] + 1.0 );
             double sin_uu, cos_uu;
             SINCOS(M_PI_2*uu, sin_uu, cos_uu);
-            const double si1 = sas_sinx_x(mu_proj * sin_uu * a_scaled);
-            const double si2 = sas_sinx_x(mu_proj * cos_uu );
-            const double si3 = sas_sinx_x(mu_proj * sin_uu * ta);
-            const double si4 = sas_sinx_x(mu_proj * cos_uu * tb);
+            const double siA = sas_sinx_x(mu_proj * sin_uu * a_over_b);
+            const double siB = sas_sinx_x(mu_proj * cos_uu );
+            const double siAt = sas_sinx_x(mu_proj * sin_uu * tA_over_b);
+            const double siBt = sas_sinx_x(mu_proj * cos_uu * tB_over_b);
 
-            // Expression in libCylinder.c (neither drC nor Vot are used)
-            const double form = scale12*si1*si2 + scale23*si2*si3 + scale14*si1*si4;
-            const double form_crim = scale11*si1*si2;
+            const double x = drV_delta*siA*siB + drA_VtA*siB*siAt + drB_VtB*siA*siBt;
+            const double form = x*siC + drC_VtC*siA*siB*siCt;
 
-            //  correct FF : sum of square of phase factors
-            inner_total += Gauss76Wt[j] * form * form;
-            inner_total_crim += Gauss76Wt[j] * form_crim * form_crim;
+            inner_sum += Gauss76Wt[j] * form * form;
         }
-        inner_total *= 0.5;
-        inner_total_crim *= 0.5;
+        inner_sum *= 0.5;
         // now sum up the outer integral
-        const double si = sas_sinx_x(mu * c_scaled * sigma);
-        const double si_crim = sas_sinx_x(mu * tc * sigma);
-        outer_total += Gauss76Wt[i] * (inner_total * si * si + inner_total_crim * si_crim * si_crim);
+        outer_sum += Gauss76Wt[i] * inner_sum;
     }
-    outer_total *= 0.5;
+    outer_sum *= 0.5;
 
     //convert from [1e-12 A-1] to [cm-1]
-    return 1.0e-4 * outer_total;
+    return 1.0e-4 * outer_sum;
 }
 
@@ -129 +139 @@
 
     double Vin = length_a * length_b * length_c;
-    double V1 = 2.0 * thick_rim_a * length_b * length_c;    // incorrect V1(aa*bb*cc+2*ta*bb*cc)
-    double V2 = 2.0 * length_a * thick_rim_b * length_c;    // incorrect V2(aa*bb*cc+2*aa*tb*cc)
-    double V3 = 2.0 * length_a * length_b * thick_rim_c;
-    // As for the 1D case, Vot is not used
-    //double Vot = (length_a * length_b * length_c +
-    //              2.0 * thick_rim_a * length_b * length_c +
-    //              2.0 * length_a * thick_rim_b * length_c +
-    //              2.0 * length_a * length_b * thick_rim_c);
+    double VtA = 2.0 * thick_rim_a * length_b * length_c;
+    double VtB = 2.0 * length_a * thick_rim_b * length_c;
+    double VtC = 2.0 * length_a * length_b * thick_rim_c;
 
     // The definitions of ta, tb, tc are not the same as in the 1D case because there is no
     // the scaling by B.
-    double ta = length_a + 2.0*thick_rim_a;
-    double tb = length_b + 2.0*thick_rim_b;
-    double tc = length_c + 2.0*thick_rim_c;
+    double tA = length_a + 2.0*thick_rim_a;
+    double tB = length_b + 2.0*thick_rim_b;
+    double tC = length_c + 2.0*thick_rim_c;
     //handle arg=0 separately, as sin(t)/t -> 1 as t->0
     double siA = sas_sinx_x(0.5*length_a*qa);
     double siB = sas_sinx_x(0.5*length_b*qb);
     double siC = sas_sinx_x(0.5*length_c*qc);
-    double siAt = sas_sinx_x(0.5*ta*qa);
-    double siBt = sas_sinx_x(0.5*tb*qb);
-    double siCt = sas_sinx_x(0.5*tc*qc);
+    double siAt = sas_sinx_x(0.5*tA*qa);
+    double siBt = sas_sinx_x(0.5*tB*qb);
+    double siCt = sas_sinx_x(0.5*tC*qc);
 
 
     // f uses Vin, V1, V2, and V3 and it seems to have more sense than the value computed
     // in the 1D code, but should be checked!
-    double f = ( dr0*siA*siB*siC*Vin
-               + drA*(siAt-siA)*siB*siC*V1
-               + drB*siA*(siBt-siB)*siC*V2
-               + drC*siA*siB*(siCt-siC)*V3);
+    double f = ( dr0*Vin*siA*siB*siC
+               + drA*VtA*(siAt-siA)*siB*siC
+               + drB*VtB*siA*(siBt-siB)*siC
+               + drC*VtC*siA*siB*(siCt-siC));
 
     return 1.0e-4 * f * f;