UltraScan III
us_eqmath.cpp
Go to the documentation of this file.
1 
3 #include "us_eqmath.h"
4 #include "us_settings.h"
5 #include "us_math2.h"
6 #include "us_constants.h"
7 #include <cfloat>
8 
9 const double dflt_min = (double)FLT_MIN;
10 const double dflt_max = (double)FLT_MAX;
11 
12 // Main constructor: pass references to EditedData, ScanEdit, EqScanFit,
13 // and EqRunFit objects needed by methods of this EqMath object
14 US_EqMath::US_EqMath(
15  QVector< US_DataIO::EditedData >& dataList,
16  QVector< ScanEdit >& scedits,
17  QVector< EqScanFit >& scanfits,
18  EqRunFit& runfit )
19  : QObject(),
20  dataList ( dataList ),
21  scedits ( scedits ),
22  scanfits ( scanfits ),
23  runfit ( runfit )
24 {
25 
26  dbg_level = US_Settings::us_debug();
27 }
28 
29 // Initialize parameters
30 void US_EqMath::init_params( int modx, bool update_mw,
31  QList< double >& ds_vbar20s, QList< double >&aud_pars )
32 {
33  modelx = modx;
34 
35  // Find the index to data corresponding to the first fitted scan
36  int jdx = -1;
37 
38  for ( int ii = 0; ii < scanfits.size(); ii++ )
39  {
40  if ( scanfits[ ii ].scanFit && jdx < 0 )
41  {
42  jdx = scedits[ ii ].dsindex;
43  break;
44  }
45  }
46 
47 DbgLv(1) << "EM:IP: jdx" << jdx << "modelx" << modelx << "update_mw" << update_mw;
48  if ( jdx < 0 ) return;
49 
50  double molecwt = runfit.mw_vals[ 0 ];
51  double portion = molecwt;
52  double mwll;
53  double mwul;
54  double uvbar;
55  double dnumc;
56  double mwinc;
57  double tempa;
58  double total;
59 
60  // Set runfit and scanfits parameters for the selected model
61  switch( modelx )
62  {
63  case 0: // 0: "1-Component, Ideal"
64  runfit.vbar_vals[ 0 ] = ds_vbar20s[ jdx ];
65  runfit.vbar_rngs[ 0 ] = runfit.vbar_vals[ 0 ] * 0.2;
66  break;
67  case 1: // 1: "2-Component, Ideal, Noninteracting"
68  if ( update_mw )
69  {
70  runfit.mw_vals[ 0 ] = portion - ( portion * 0.2 );
71  runfit.mw_vals[ 1 ] = portion + ( portion * 0.2 );
72  }
73 
74  for ( int ii = 0; ii < runfit.nbr_comps; ii++ )
75  {
76  runfit.mw_rngs [ ii ] = runfit.mw_vals[ ii ] * 0.2;
77  runfit.vbar_vals[ ii ] = ds_vbar20s[ jdx ];
78  runfit.vbar_rngs[ ii ] = runfit.vbar_vals[ ii ] * 0.2;
79  }
80 
81  if ( (int)( 100.0 * runfit.vbar_vals[ 1 ] ) == 72 )
82  {
83  runfit.vbar_vals[ 1 ] = runfit.vbar_vals[ 0 ];
84  runfit.vbar_rngs[ 1 ] = runfit.vbar_vals[ 1 ] * 0.2;
85  }
86 
87  for ( int ii = 0; ii < scanfits.size(); ii++ )
88  {
89  if ( ! scanfits[ ii ].scanFit ) continue;
90 
91  portion = exp( scanfits[ ii ].amp_vals[ 0 ] ) * 0.5;
92  scanfits[ ii ].amp_vals[ 0 ] = log( portion * 0.7 );
93  scanfits[ ii ].amp_vals[ 1 ] = log( portion * 0.3 );
94  scanfits[ ii ].amp_rngs[ 0 ] = scanfits[ ii ].amp_vals[ 0 ] * 0.2;
95  scanfits[ ii ].amp_rngs[ 1 ] = scanfits[ ii ].amp_rngs[ 0 ];
96  }
97  break;
98  case 2: // 2: "3-Component, Ideal, Noninteracting"
99  if ( update_mw )
100  {
101  runfit.mw_vals[ 0 ] = molecwt - ( molecwt * 0.2 );
102  runfit.mw_vals[ 1 ] = molecwt;
103  runfit.mw_vals[ 2 ] = molecwt + ( molecwt * 0.2 );
104  }
105 
106 DbgLv(1) << "EM:IP: C2: nbrcomps" << runfit.nbr_comps;
107  for ( int ii = 0; ii < runfit.nbr_comps; ii++ )
108  {
109  runfit.mw_rngs [ ii ] = runfit.mw_vals[ ii ] * 0.2;
110  runfit.vbar_vals[ ii ] = ds_vbar20s[ jdx ];
111  runfit.vbar_rngs[ ii ] = runfit.vbar_vals[ ii ] * 0.2;
112  }
113 
114  if ( (int)( 100.0 * runfit.vbar_vals[ 1 ] ) == 72 )
115  {
116  runfit.vbar_vals[ 1 ] = runfit.vbar_vals[ 0 ];
117  runfit.vbar_rngs[ 1 ] = runfit.vbar_vals[ 1 ] * 0.2;
118  }
119  if ( (int)( 100.0 * runfit.vbar_vals[ 2 ] ) == 72 )
120  {
121  runfit.vbar_vals[ 2 ] = runfit.vbar_vals[ 0 ];
122  runfit.vbar_rngs[ 2 ] = runfit.vbar_vals[ 2 ] * 0.2;
123  }
124 
125  total = runfit.mw_vals[ 0 ] + runfit.mw_vals[ 1 ]
126  + runfit.mw_vals[ 2 ];
127 
128 DbgLv(1) << "EM:IP: C2: total" << total;
129  for ( int ii = 0; ii < scanfits.size(); ii++ )
130  {
131  if ( ! scanfits[ ii ].scanFit ) continue;
132 
133  portion = exp( scanfits[ ii ].amp_vals[ 0 ] ) / 3.0;
134  scanfits[ ii ].amp_vals[ 0 ] = log( portion * 0.6 );
135  scanfits[ ii ].amp_vals[ 1 ] = log( portion * 0.3 );
136  scanfits[ ii ].amp_vals[ 2 ] = log( portion * 0.1 );
137  scanfits[ ii ].amp_rngs[ 0 ] = scanfits[ ii ].amp_vals[ 0 ] * 0.2;
138  scanfits[ ii ].amp_rngs[ 1 ] = scanfits[ ii ].amp_rngs[ 0 ];
139  scanfits[ ii ].amp_rngs[ 2 ] = scanfits[ ii ].amp_rngs[ 0 ];
140  }
141  break;
142  case 3: // 3: "Fixed Molecular Weight Distribution"
143  runfit.nbr_comps = aud_pars[ 0 ];
144  mwll = aud_pars[ 1 ];
145  mwul = aud_pars[ 2 ];
146  uvbar = aud_pars[ 3 ];
147  uvbar = uvbar > 0.0 ? uvbar : ds_vbar20s[ jdx ];
148  dnumc = (double)runfit.nbr_comps;
149  mwinc = ( mwul - mwll ) / ( dnumc - 1.0 );
150  tempa = log( 1.0e-7 / dnumc );
151 
152  for ( int ii = 0; ii < scanfits.size(); ii++ )
153  {
154  if ( ! scanfits[ ii ].scanFit ) continue;
155 
156  for ( int jj = 0; jj < runfit.nbr_comps; jj++ )
157  {
158  scanfits[ ii ].amp_vals[ jj ] = tempa;
159  scanfits[ ii ].amp_rngs[ jj ] =
160  scanfits[ ii ].amp_vals[ jj ] * 0.2;
161  scanfits[ ii ].baseline = 0.0;
162  scanfits[ ii ].baseln_rng = 0.05;
163  runfit.mw_vals[ jj ] = mwll + ( (double)jj * mwinc );
164  runfit.mw_rngs[ jj ] = runfit.mw_vals[ jj ] * 0.2;
165  runfit.vbar_vals[ jj ] = uvbar;
166  runfit.vbar_rngs[ jj ] = uvbar * 0.2;
167  }
168  }
169  break;
170  case 4: // 4: "Monomer-Dimer Equilibrium"
171  case 5: // 5: "Monomer-Trimer Equilibrium"
172  case 6: // 6: "Monomer-Tetramer Equilibrium"
173  case 7: // 7: "Monomer-Pentamer Equilibrium"
174  case 8: // 8: "Monomer-Hexamer Equilibrium"
175  case 9: // 9: "Monomer-Heptamer Equilibrium"
176  case 10: // 10: "User-Defined Monomer-Nmer Equilibrium"
177  case 11: // 11: "Monomer-Dimer-Trimer Equilibrium"
178  case 12: // 12: "Monomer-Dimer-Tetramer Equilibrium"
179  case 13: // 13: "User-Defined Monomer - N-mer - M-mer Equilibrium"
180  runfit.eq_vals[ 0 ] = -1.0;
181  runfit.eq_vals[ 1 ] = -1.0e4;
182  runfit.eq_rngs[ 0 ] = 5.0;
183  runfit.eq_rngs[ 1 ] = 5.0;
184  runfit.vbar_vals[ 0 ] = ds_vbar20s[ jdx ];
185  runfit.vbar_rngs[ 0 ] = runfit.vbar_vals[ 0 ] * 0.2;
186  break;
187  case 14: // 14: "2-Component Hetero-Association: A + B <=> AB"
188  if ( update_mw )
189  {
190  runfit.mw_vals[ 0 ] = molecwt - ( molecwt * 0.2 );
191  runfit.mw_vals[ 1 ] = molecwt + ( molecwt * 0.2 );
192  }
193 
194  runfit.mw_rngs [ 0 ] = runfit.mw_vals[ 0 ] * 0.2;
195  runfit.mw_rngs [ 1 ] = runfit.mw_vals[ 1 ] * 0.2;
196  runfit.vbar_vals[ 0 ] = ds_vbar20s[ jdx ];
197  runfit.vbar_rngs[ 0 ] = runfit.vbar_vals[ 0 ] * 0.2;
198  runfit.vbar_vals[ 1 ] = ds_vbar20s[ jdx ];
199  runfit.vbar_rngs[ 1 ] = runfit.vbar_vals[ 1 ] * 0.2;
200 
201  if ( (int)( 100.0 * runfit.vbar_vals[ 1 ] ) == 72 )
202  {
203  runfit.vbar_vals[ 1 ] = runfit.vbar_vals[ 0 ];
204  runfit.vbar_rngs[ 1 ] = runfit.vbar_vals[ 1 ] * 0.2;
205  }
206 
207  runfit.eq_vals[ 0 ] = -1.0e4;
208  runfit.eq_rngs[ 0 ] = 5.0e2;
209 
210  for ( int ii = 0; ii < scanfits.size(); ii++ )
211  {
212  if ( ! scanfits[ ii ].scanFit ) continue;
213 
214  portion = exp( scanfits[ ii ].amp_vals[ 0 ] ) / 2.0;
215  scanfits[ ii ].amp_vals[ 0 ] = log( portion );
216  scanfits[ ii ].amp_vals[ 1 ] = scanfits[ ii ].amp_vals[ 0 ];
217  scanfits[ ii ].amp_rngs[ 0 ] = scanfits[ ii ].amp_vals[ 0 ] * 0.2;
218  scanfits[ ii ].amp_rngs[ 1 ] = scanfits[ ii ].amp_rngs[ 0 ];
219  }
220  break;
221  case 15: // 15: "U-Defined self/Hetero-Assoc.: A + B <=> AB, nA <=> An"
222  if ( update_mw )
223  {
224  runfit.mw_vals[ 0 ] = molecwt - ( molecwt * 0.2 );
225  runfit.mw_vals[ 1 ] = molecwt + ( molecwt * 0.2 );
226  }
227 
228  runfit.mw_rngs [ 0 ] = runfit.mw_vals[ 0 ] * 0.2;
229  runfit.mw_rngs [ 1 ] = runfit.mw_vals[ 1 ] * 0.2;
230  runfit.vbar_vals[ 0 ] = ds_vbar20s[ jdx ];
231  runfit.vbar_rngs[ 0 ] = runfit.vbar_vals[ 0 ] * 0.2;
232  runfit.vbar_vals[ 1 ] = ds_vbar20s[ jdx ];
233  runfit.vbar_rngs[ 1 ] = runfit.vbar_vals[ 1 ] * 0.2;
234 
235  if ( (int)( 100.0 * runfit.vbar_vals[ 1 ] ) == 72 )
236  {
237  runfit.vbar_vals[ 1 ] = runfit.vbar_vals[ 0 ];
238  runfit.vbar_rngs[ 1 ] = runfit.vbar_vals[ 1 ] * 0.2;
239  }
240 
241  runfit.eq_vals[ 0 ] = -1.0e4;
242  runfit.eq_rngs[ 0 ] = 5.0e2;
243  runfit.eq_vals[ 1 ] = -1.0e4;
244  runfit.eq_rngs[ 1 ] = 5.0e2;
245 
246  for ( int ii = 0; ii < scanfits.size(); ii++ )
247  {
248  if ( ! scanfits[ ii ].scanFit ) continue;
249 
250  portion = exp( scanfits[ ii ].amp_vals[ 0 ] ) / 2.0;
251  scanfits[ ii ].amp_vals[ 0 ] = log( portion );
252  scanfits[ ii ].amp_vals[ 1 ] = scanfits[ ii ].amp_vals[ 0 ];
253  scanfits[ ii ].amp_rngs[ 0 ] = scanfits[ ii ].amp_vals[ 0 ] * 0.2;
254  scanfits[ ii ].amp_rngs[ 1 ] = scanfits[ ii ].amp_rngs[ 0 ];
255  }
256  break;
257  case 16: // 16: "U-Defined Monomer-Nmer, some monomer is incompetent"
258  runfit.mw_vals [ 1 ] = runfit.mw_vals[ 0 ];
259  runfit.mw_rngs [ 0 ] = runfit.mw_vals[ 0 ] * 0.2;
260  runfit.mw_rngs [ 1 ] = runfit.mw_vals[ 1 ] * 0.2;
261  runfit.vbar_vals[ 0 ] = ds_vbar20s[ jdx ];
262  runfit.vbar_rngs[ 0 ] = runfit.vbar_vals[ 0 ] * 0.2;
263  runfit.vbar_vals[ 1 ] = ds_vbar20s[ jdx ];
264  runfit.vbar_rngs[ 1 ] = runfit.vbar_vals[ 1 ] * 0.2;
265  runfit.eq_vals[ 0 ] = -0.1;
266  runfit.eq_rngs[ 0 ] = 5.0;
267 
268  for ( int ii = 0; ii < scanfits.size(); ii++ )
269  {
270  if ( ! scanfits[ ii ].scanFit ) continue;
271 
272  scanfits[ ii ].amp_vals[ 1 ] = -1.0e4;
273  scanfits[ ii ].amp_rngs[ 1 ] = 1.0e-3;
274  }
275  break;
276  case 17: // 17: "User-Defined Monomer-Nmer, some Nmer is incompetent"
277  runfit.mw_vals [ 1 ] = runfit.stoichs[ 0 ] * runfit.mw_vals[ 0 ];
278  runfit.mw_rngs [ 0 ] = runfit.mw_vals[ 0 ] * 0.2;
279  runfit.mw_rngs [ 1 ] = runfit.mw_vals[ 1 ] * 0.2;
280  runfit.vbar_vals[ 0 ] = ds_vbar20s[ jdx ];
281  runfit.vbar_rngs[ 0 ] = runfit.vbar_vals[ 0 ] * 0.2;
282  runfit.vbar_vals[ 1 ] = ds_vbar20s[ jdx ];
283  runfit.vbar_rngs[ 1 ] = runfit.vbar_vals[ 1 ] * 0.2;
284  runfit.eq_vals[ 0 ] = -0.1;
285  runfit.eq_rngs[ 0 ] = 5.0;
286 
287  for ( int ii = 0; ii < scanfits.size(); ii++ )
288  {
289  if ( ! scanfits[ ii ].scanFit ) continue;
290 
291  scanfits[ ii ].amp_vals[ 1 ] = -1.0e4;
292  scanfits[ ii ].amp_rngs[ 1 ] = 1.0e-3;
293  }
294  break;
295  case 18: // 18: "User-Defined irreversible Monomer-Nmer"
296  runfit.mw_vals [ 1 ] = runfit.stoichs[ 0 ] * runfit.mw_vals[ 0 ];
297  runfit.mw_rngs [ 0 ] = runfit.mw_vals[ 0 ] * 0.2;
298  runfit.mw_rngs [ 1 ] = runfit.mw_vals[ 1 ] * 0.2;
299  runfit.vbar_vals[ 0 ] = ds_vbar20s[ jdx ];
300  runfit.vbar_rngs[ 0 ] = runfit.vbar_vals[ 0 ] * 0.2;
301  runfit.vbar_vals[ 1 ] = ds_vbar20s[ jdx ];
302  runfit.vbar_rngs[ 1 ] = runfit.vbar_vals[ 1 ] * 0.2;
303 
304  for ( int ii = 0; ii < scanfits.size(); ii++ )
305  {
306  if ( ! scanfits[ ii ].scanFit ) continue;
307 
308  portion = exp( scanfits[ ii ].amp_vals[ 0 ] ) / 2.0;
309  scanfits[ ii ].amp_vals[ 0 ] = log( portion * 0.7 );
310  scanfits[ ii ].amp_vals[ 1 ] = log( portion * 0.3 );
311  scanfits[ ii ].amp_rngs[ 0 ] = scanfits[ ii ].amp_vals[ 0 ] * 0.2;
312  scanfits[ ii ].amp_rngs[ 1 ] = scanfits[ ii ].amp_rngs[ 0 ];
313  }
314  break;
315  case 19: // 19: "User-Defined Monomer-Nmer plus contaminant"
316  runfit.mw_vals [ 1 ] = runfit.mw_vals[ 0 ];
317  runfit.mw_rngs [ 0 ] = runfit.mw_vals[ 0 ] * 0.2;
318  runfit.mw_rngs [ 1 ] = runfit.mw_vals[ 1 ] * 0.2;
319  runfit.vbar_vals[ 0 ] = ds_vbar20s[ jdx ];
320  runfit.vbar_rngs[ 0 ] = runfit.vbar_vals[ 0 ] * 0.2;
321  runfit.vbar_vals[ 1 ] = ds_vbar20s[ jdx ];
322  runfit.vbar_rngs[ 1 ] = runfit.vbar_vals[ 1 ] * 0.2;
323  runfit.eq_vals[ 0 ] = -1.0e4;
324  runfit.eq_rngs[ 0 ] = 5.0;
325 
326  for ( int ii = 0; ii < scanfits.size(); ii++ )
327  {
328  if ( ! scanfits[ ii ].scanFit ) continue;
329 
330  scanfits[ ii ].amp_vals[ 1 ] = -1.0e4;
331  scanfits[ ii ].amp_rngs[ 1 ] = 1.0e-3;
332  }
333  break;
334  }
335 #if 0
336 double aaa[ 12 ]; double ppp[ 16 ];
337 for ( int ii=0;ii<12;ii++ ) aaa[ii]=(double)( ii+1 );
338 for ( int ii=0;ii<16;ii++ ) ppp[ii]=0.0;
339 double* aa[3]; double* pp[4];
340 aa[0]=aaa;aa[1]=aaa+4;aa[2]=aaa+8;
341 pp[0]=ppp;pp[1]=ppp+4;pp[2]=ppp+8;pp[3]=ppp+12;
342 DbgLv(2) << "AA: " << aaa[0] << aaa[1] << aaa[2] << aaa[3];
343 DbgLv(2) << "AA: " << aaa[4] << aaa[5] << aaa[6] << aaa[7];
344 DbgLv(2) << "AA: " << aaa[8] << aaa[9] << aaa[10] << aaa[11];
345 US_Matrix::tmm( aa, pp, 3, 4, false );
346 DbgLv(2) << "==US_Matrix::tmm( aa, pp, 4, 3, false )==";
347 DbgLv(2) << " PP: " << ppp[0] << ppp[1] << ppp[2] << ppp[3];
348 DbgLv(2) << " PP: " << ppp[4] << ppp[5] << ppp[6] << ppp[7];
349 DbgLv(2) << " PP: " << ppp[8] << ppp[9] << ppp[10] << ppp[11];
350 DbgLv(2) << " PP: " << ppp[12] << ppp[13] << ppp[14] << ppp[15];
351 US_Matrix::tmm( aa, pp, 3, 4, true );
352 DbgLv(2) << "==US_Matrix::tmm( aa, pp, 4, 3, true )==";
353 DbgLv(2) << " PP: " << ppp[0] << ppp[1] << ppp[2] << ppp[3];
354 DbgLv(2) << " PP: " << ppp[4] << ppp[5] << ppp[6] << ppp[7];
355 DbgLv(2) << " PP: " << ppp[8] << ppp[9] << ppp[10] << ppp[11];
356 DbgLv(2) << " PP: " << ppp[12] << ppp[13] << ppp[14] << ppp[15];
357 #endif
358 
359 
360 }
361 
362 // Initialize for fit
363 void US_EqMath::init_fit( int modx, int methx, FitCtrlPar& fitpars )
364 {
365  modelx = modx; // Model type index
366  nlsmeth = methx; // NLS method index
367 DbgLv(1) << "EM:IF: modelx nlsmeth nc" << modelx << nlsmeth << runfit.nbr_comps;
368 
369  // Find the index to the first fitted scan and count data sets, points
370  ffitx = -1;
371  ntpts = 0;
372  ndsets = 0;
373  nfpars = 0;
374  nspts = 0;
375  nslpts = 0;
376  v_setpts .clear();
377  v_setlpts.clear();
378 
379  for ( int ii = 0; ii < scanfits.size(); ii++ )
380  {
381  if ( ! scanfits[ ii ].scanFit ) continue;
382 
383  if ( ffitx < 0 ) ffitx = ii;
384 
385  nspts = scanfits[ ii ].stop_ndx - scanfits[ ii ].start_ndx + 1;
386  v_setpts << nspts;
387  v_setlpts << nspts;
388  ndsets++;
389  ntpts += nspts;
390  }
391 DbgLv(1) << "EM:IF: scan 1 ntpts ndsets nfpars" << ntpts << ndsets << nfpars;
392 
393  if ( ffitx < 0 || ntpts == 0 ) return;
394 
395  // Count the number of fitting parameters
396  for ( int jj = 0; jj < runfit.nbr_comps; jj++ )
397  {
398  if ( runfit.mw_fits[ jj ] == true )
399  nfpars++;
400 
401  if ( runfit.vbar_fits[ jj ] == true )
402  nfpars++;
403 
404  if ( runfit.viri_fits[ jj ] == true )
405  nfpars++;
406  }
407 DbgLv(1) << "EM:IF: scan 2 ntpts ndsets nfpars" << ntpts << ndsets << nfpars;
408 
409  for ( int ii = 0; ii < scanfits.size(); ii++ )
410  {
411  if ( ! scanfits[ ii ].scanFit ) continue;
412 
413  EqScanFit* scnf = &scanfits[ ii ];
414 
415  for ( int jj = 0; jj < runfit.nbr_comps; jj++ )
416  if ( scnf->amp_fits[ jj ] )
417  nfpars++;
418 
419  if ( scnf->baseln_fit )
420  nfpars++;
421  }
422 DbgLv(1) << "EM:IF: scan 3 ntpts ndsets nfpars" << ntpts << ndsets << nfpars;
423 
424  for ( int jj = 0; jj < runfit.nbr_assocs; jj++ )
425  if ( runfit.eq_fits[ jj ] )
426  nfpars++;
427 DbgLv(1) << "EM:IF: scan 4 ntpts ndsets nfpars" << ntpts << ndsets << nfpars;
428 
429  // Allocate vectors and matrices
430  v_yraw .fill( 0.0, ntpts );
431  v_yguess.fill( 0.0, ntpts );
432  v_ydelta.fill( 0.0, ntpts );
433  v_BB .fill( 0.0, nfpars );
434  v_guess .fill( 0.0, nfpars );
435  v_tguess.fill( 0.0, nfpars );
436 
437  m_jacobi.fill( 0, ntpts );
438  m_info .fill( 0, nfpars );
439  m_LLtrns.fill( 0, nfpars );
440  m_dcr2 .fill( 0, ndsets );
441  m_dlncr2.fill( 0, ndsets );
442  m_lncr2 .fill( 0, ndsets );
443 
444  int njacs = ntpts * nfpars;
445  int nfpsq = nfpars * nfpars;
446  v_jacobi.fill( 0.0, njacs );
447  v_info .fill( 0.0, nfpsq );
448  v_LLtrns.fill( 0.0, nfpsq );
449  v_lncr2 .fill( 0.0, ntpts );
450  v_dlncr2.fill( 0.0, ntpts );
451  v_dcr2 .fill( 0.0, ntpts );
452 
453  int dsx = 0;
454  int ptx = 0;
455  d_lncr2 = v_lncr2 .data();
456  d_dlncr2 = v_dlncr2.data();
457  d_dcr2 = v_dcr2 .data();
458 
459  for ( int ii = 0; ii < scanfits.size(); ii++ )
460  {
461  if ( ! scanfits[ ii ].scanFit ) continue;
462 
463  EqScanFit* scnf = &scanfits[ ii ];
464 
465  m_dlncr2[ dsx ] = d_dlncr2;
466  m_lncr2 [ dsx ] = d_lncr2;
467  m_dcr2 [ dsx ] = d_dcr2;
468 
469  nspts = v_setpts[ dsx++ ];
470  d_dlncr2 += nspts;
471  d_lncr2 += nspts;
472  d_dcr2 += nspts;
473  int jy = scnf->start_ndx;
474 
475  for ( int jj = 0; jj < nspts; jj++ )
476  v_yraw[ ptx++ ] = scnf->yvs[ jy++ ];
477  }
478 DbgLv(1) << "EM:IF: matr fill 1 complete";
479 
480  // Initialize parameter guess
481  guess_mapForward( v_guess.data() );
482 DbgLv(1) << "EM:IF: g map Forw complete";
483 
484  // Set up jacobian, info, and LLtranspose matrices
485  d_jacobi = v_jacobi.data();
486  d_info = v_info .data();
487  d_LLtrns = v_LLtrns.data();
488 
489  for ( int ii = 0; ii < ntpts; ii++ )
490  {
491  m_jacobi[ ii ] = d_jacobi;
492  d_jacobi += nfpars;
493  }
494 DbgLv(1) << "EM:IF: jacobi fill complete";
495 
496  for ( int ii = 0; ii < nfpars; ii++ )
497  {
498  m_info [ ii ] = d_info;
499  m_LLtrns[ ii ] = d_LLtrns;
500  d_info += nfpars;
501  d_LLtrns += nfpars;
502  }
503 DbgLv(1) << "EM:IF: inf/trns fill complete";
504  setpts = v_setpts .data();
505  setlpts = v_setlpts.data();
506  y_raw = v_yraw .data();
507  y_guess = v_yguess .data();
508  y_delta = v_ydelta .data();
509  BB = v_BB .data();
510  guess = v_guess .data();
511  tguess = v_tguess .data();
512  jacobian = m_jacobi .data();
513  info = m_info .data();
514  LLtr = m_LLtrns .data();
515  dcr2 = m_dcr2 .data();
516  dlncr2 = m_dlncr2 .data();
517  lncr2 = m_lncr2 .data();
518 
519 DbgLv(1) << "EM:FI: ffitx" << ffitx << "modelx" << modelx;
520  // Return counts and other fit parameters to caller
521  fitpars.nlsmeth = nlsmeth;
522  fitpars.modelx = modelx;
523  fitpars.ntpts = ntpts;
524  fitpars.ndsets = ndsets;
525  fitpars.nfpars = nfpars;
526  fitpars.setpts = setpts;
527  fitpars.setlpts = setlpts;
528  fitpars.y_raw = y_raw;
529  fitpars.y_guess = y_guess;
530  fitpars.y_delta = y_delta;
531  fitpars.BB = BB;
532  fitpars.guess = guess;
533  fitpars.tguess = tguess;
534  fitpars.jacobian = jacobian;
535  fitpars.info = info;
536  fitpars.LLtr = LLtr;
537  fitpars.dcr2 = dcr2;
538  fitpars.dlncr2 = dlncr2;
539  fitpars.lncr2 = lncr2;
540 DbgLv(1) << "EM:FI: ktpts kdsets kfpars" << ntpts << ndsets << nfpars;
541 }
542 
543 // Fill in guess parameter vector
544 // Parameter order (needs to be maintained so Jacobian columns match):
545 // for each component k:
546 // 1. Molecular Weight ( k )
547 // 2. Vbar ( k )
548 // 3. Virial Coefficient ( k )
549 //
550 // for each scan:
551 // for each component:
552 // 4. Amplitude
553 // 5. Baseline
554 //
555 // for each association constant:
556 // 6. Association constant
557 //
558 void US_EqMath::guess_mapForward( double* vguess )
559 {
560  int jpx = 0;
561 DbgLv(1) << "EM:gmF: ncomps" << runfit.nbr_comps;
562 
563  for ( int jj = 0; jj < runfit.nbr_comps; jj++ )
564  {
565 //DbgLv(1) << "EM:gmF: jj jpx" << jj << jpx;
566  if ( runfit.mw_fits[ jj ] )
567  {
568  vguess[ jpx ] = runfit.mw_vals[ jj ];
569  runfit.mw_ndxs[ jj ] = jpx++;
570  }
571 
572  if ( runfit.vbar_fits[ jj ] )
573  {
574  vguess[ jpx ] = runfit.vbar_vals[ jj ];
575  runfit.vbar_ndxs[ jj ] = jpx++;
576  }
577 
578  if ( runfit.viri_fits[ jj ] )
579  {
580  vguess[ jpx ] = runfit.viri_vals[ jj ];
581  runfit.viri_ndxs[ jj ] = jpx++;
582  }
583  }
584 
585  for ( int ii = 0; ii < scanfits.size(); ii++ )
586  {
587  if ( ! scanfits[ ii ].scanFit ) continue;
588 
589 //DbgLv(1) << "EM:gmF: ii jpx" << ii << jpx;
590  EqScanFit* scnf = &scanfits[ ii ];
591 
592  for ( int jj = 0; jj < runfit.nbr_comps; jj++ )
593  {
594  if ( scnf->amp_fits[ jj ] )
595  {
596  vguess[ jpx ] = scnf->amp_vals[ jj ];
597  scnf->amp_ndxs[ jj ] = jpx++;
598  }
599  }
600 
601  if ( scnf->baseln_fit )
602  {
603  vguess[ jpx ] = scnf->baseline;
604  scnf->baseln_ndx = jpx++;
605  }
606  }
607 
608  for ( int jj = 0; jj < runfit.nbr_assocs; jj++ )
609  {
610 //DbgLv(1) << "EM:gmF: as_jj jpx" << jj << jpx;
611  if ( runfit.eq_fits[ jj ] )
612  {
613  vguess[ jpx ] = runfit.eq_vals[ jj ];
614  runfit.eq_ndxs[ jj ] = jpx++;
615  }
616  }
617 DbgLv(1) << "EM:gmF: jpx" << jpx;
618 }
619 
620 // Map parameters back from guesses
621 void US_EqMath::parameter_mapBackward( double* vguess )
622 {
623  int jpx = 0;
624 
625  for ( int jj = 0; jj < runfit.nbr_comps; jj++ )
626  {
627  if ( runfit.mw_fits[ jj ] )
628  runfit.mw_vals[ jj ] = vguess[ jpx++ ];
629 
630  if ( runfit.vbar_fits[ jj ] )
631  runfit.vbar_vals[ jj ] = vguess[ jpx++ ];
632 
633  if ( runfit.viri_fits[ jj ] )
634  runfit.viri_vals[ jj ] = vguess[ jpx++ ];
635  }
636 
637  for ( int ii = 0; ii < scanfits.size(); ii++ )
638  {
639  if ( ! scanfits[ ii ].scanFit ) continue;
640 
641  EqScanFit* scnf = &scanfits[ ii ];
642 
643  for ( int jj = 0; jj < runfit.nbr_comps; jj++ )
644  {
645  if ( scnf->amp_fits[ jj ] )
646  scnf->amp_vals[ jj ] = vguess[ jpx++ ];
647 
648  if ( scnf->baseln_fit )
649  scnf->baseline = vguess[ jpx++ ];
650  }
651  }
652 
653  for ( int jj = 0; jj < runfit.nbr_assocs; jj++ )
654  if ( runfit.eq_fits[ jj ] )
655  runfit.eq_vals[ jj ] = vguess[ jpx++ ];
656 
657  if ( modelx == 3 )
658  {
659  EqScanFit* sc1f = &scanfits[ ffitx ];
660 
661  for ( int ii = ffitx + 1; ii < scanfits.size(); ii++ )
662  {
663  EqScanFit* scnf = &scanfits[ ii ];
664 
665  if ( ! scnf->scanFit ) continue;
666 
667  for ( int jj = 0; jj < runfit.nbr_comps; jj++ )
668  scnf->amp_vals[ jj ] = sc1f->amp_vals[ jj ];
669  }
670  }
671 
672 }
673 
674 // Calculate the jacobian matrix
675 int US_EqMath::calc_jacobian()
676 {
677  int stat = 0;
678  int ncomps = runfit.nbr_comps;
679  int mcomp = max( ncomps, 4 );
680  int jpx = 0;
681  int jdx = 0;
682  QVector< double > v_ufunc( mcomp );
683  QVector< double > v_vbar ( mcomp );
684  QVector< double > v_buoy ( mcomp );
685 
686  v_jacobi.fill( 0.0, ntpts * nfpars );
687 
688  switch( modelx )
689  {
690  case 0: // 0: "1-Component, Ideal"
691  case 1: // 1: "2-Component, Ideal, Noninteracting"
692  case 2: // 2: "3-Component, Ideal, Noninteracting"
693  case 3: // 3: "Fixed Molecular Weight Distribution"
694  for ( int ii = 0; ii < scanfits.size(); ii++ )
695  {
696  EqScanFit* scnf = &scanfits[ ii ];
697 
698  if ( ! scnf->scanFit ) continue;
699 
700  int jstx = scnf->start_ndx;
701  double xm_sq = sq( scnf->xvs[ jstx ] );
702  double omega_sq = sq( M_PI * scnf->rpm / 30.0 );
703  double tempera = scnf->tempera;
704  double density = scnf->density;
705  double dconst = omega_sq / ( 2.0 * R * ( K0 + tempera ) );
706 
707  for ( int kk = 0; kk < ncomps; kk++ )
708  {
709  v_vbar[ kk ] = US_Math2::adjust_vbar20( runfit.vbar_vals[ kk ],
710  tempera );
711  v_buoy[ kk ] = ( 1.0 - v_vbar[ kk ] * density );
712  }
713 
714  for ( int jj = jstx; jj < jstx + v_setpts[ jdx ]; jj++ )
715  {
716  double xv = sq( scnf->xvs[ jj ] ) - xm_sq;
717 
718  for ( int kk = 0; kk < ncomps; kk++ )
719  {
720  double buoy = v_buoy[ kk ];
721  double mwv = runfit.mw_vals[ kk ];
722  double ampv = scnf->amp_vals[ kk ];
723  double ufunc = exp( ampv + dconst * mwv * buoy * xv );
724  v_ufunc[ kk ] = ufunc;
725 
726  if ( runfit.mw_fits[ kk ] )
727  m_jacobi[ jpx ][ runfit.mw_ndxs[ kk ] ]
728  = dconst * xv * buoy * ufunc;
729 
730  if ( runfit.vbar_fits[ kk ] )
731  m_jacobi[ jpx ][ runfit.vbar_ndxs[ kk ] ]
732  = (-1.0 ) * dconst * mwv * xv * ufunc * density;
733 
734  if ( scnf->amp_fits[ kk ] )
735  m_jacobi[ jpx ][ scnf->amp_ndxs[ kk ] ] = ufunc;
736  }
737 
738  if ( scnf->baseln_fit )
739  m_jacobi[ jpx ][ scnf->baseln_ndx ] = 1.0;
740 
741  jpx++;
742  }
743 
744  jdx++;
745  }
746  break;
747  case 4: // 4: "Monomer-Dimer Equilibrium"
748  case 5: // 5: "Monomer-Trimer Equilibrium"
749  case 6: // 6: "Monomer-Tetramer Equilibrium"
750  case 7: // 7: "Monomer-Pentamer Equilibrium"
751  case 8: // 8: "Monomer-Hexamer Equilibrium"
752  case 9: // 9: "Monomer-Heptamer Equilibrium"
753  case 10: // 10: "User-Defined Monomer-Nmer Equilibrium"
754  {
755  double stoich1 = (double)( modelx - 2 );
756  double stoiexp = stoich1 - 1.0;
757  runfit.stoichs[ 0 ] = stoich1;
758 
759  for ( int ii = 0; ii < scanfits.size(); ii++ )
760  {
761  EqScanFit* scnf = &scanfits[ ii ];
762 
763  if ( ! scnf->scanFit ) continue;
764 
765  int jstx = scnf->start_ndx;
766  double xm_sq = sq( scnf->xvs[ jstx ] );
767  double omega_sq = sq( M_PI * scnf->rpm / 30.0 );
768  double tempera = scnf->tempera;
769  double density = scnf->density;
770  double dconst = omega_sq / ( 2.0 * R * ( K0 + tempera ) );
771  double ODcorrec = log( stoich1 / pow( scnf->extincts[ 0 ]
772  * scnf->pathlen , stoiexp ) );
773  v_vbar[ 0 ] = US_Math2::adjust_vbar20( runfit.vbar_vals[ 0 ],
774  tempera );
775  v_buoy[ 0 ] = ( 1.0 - v_vbar[ 0 ] * density );
776  double buoy = v_buoy[ 0 ];
777 
778  for ( int jj = jstx; jj < jstx + v_setpts[ jdx ]; jj++ )
779  {
780  double xv = sq( scnf->xvs[ jj ] ) - xm_sq;
781  double mwv = runfit.mw_vals[ 0 ];
782  double ampv = scnf->amp_vals[ 0 ];
783  double ufunc0 = exp( ampv + dconst * mwv * buoy * xv );
784  double ufunc1 = exp( stoich1 * ampv + runfit.eq_vals[ 0 ] +
785  ODcorrec + dconst * stoich1 * mwv * buoy * ufunc0 );
786  v_ufunc[ 0 ] = ufunc0;
787  v_ufunc[ 1 ] = ufunc1;
788  double dcoeff = dconst * xv * buoy;
789 
790  if ( runfit.mw_fits[ 0 ] )
791  m_jacobi[ jpx ][ runfit.mw_ndxs[ 0 ] ]
792  = dcoeff * ufunc0 + dcoeff * ufunc1 * stoich1;
793 
794  dcoeff = dconst * mwv * xv * density;
795 
796  if ( runfit.vbar_fits[ 0 ] )
797  m_jacobi[ jpx ][ runfit.vbar_ndxs[ 0 ] ]
798  = (-1.0 ) * dcoeff * ufunc0 - dcoeff * ufunc1 * stoich1;
799 
800  if ( scnf->amp_fits[ 0 ] )
801  m_jacobi[ jpx ][ scnf->amp_ndxs[ 0 ] ]
802  = ufunc0 + ufunc1 * stoich1;
803 
804  if ( runfit.eq_fits[ 0 ] )
805  m_jacobi[ jpx ][ runfit.eq_ndxs[ 0 ] ] = ufunc1;
806 
807  if ( scnf->baseln_fit )
808  m_jacobi[ jpx ][ scnf->baseln_ndx ] = 1.0;
809 
810  jpx++;
811  }
812 
813  jdx++;
814  }
815  break;
816  }
817  case 11: // 11: "Monomer-Dimer-Trimer Equilibrium"
818  case 12: // 12: "Monomer-Dimer-Tetramer Equilibrium"
819  case 13: // 13: "User-Defined Monomer - N-mer - M-mer Equilibrium"
820  {
821  double stoich1 = 2.0;
822  double stoich2 = (double)( modelx - 8 );
823  runfit.stoichs[ 0 ] = stoich1;
824  runfit.stoichs[ 1 ] = stoich2;
825  double stoiexp = stoich1 - 1.0;
826  double stoiex2 = stoich2 - 1.0;
827 
828  for ( int ii = 0; ii < scanfits.size(); ii++ )
829  {
830  EqScanFit* scnf = &scanfits[ ii ];
831 
832  if ( ! scnf->scanFit ) continue;
833 
834  int jstx = scnf->start_ndx;
835  double xm_sq = sq( scnf->xvs[ jstx ] );
836  double omega_sq = sq( M_PI * scnf->rpm / 30.0 );
837  double tempera = scnf->tempera;
838  double density = scnf->density;
839  double dconst = omega_sq / ( 2.0 * R * ( K0 + tempera ) );
840  double ODcorec1 = log( stoich1 / pow( scnf->extincts[ 0 ]
841  * scnf->pathlen , stoiexp ) );
842  double ODcorec2 = log( stoich2 / pow( scnf->extincts[ 0 ]
843  * scnf->pathlen , stoiex2 ) );
844  v_vbar[ 0 ] = US_Math2::adjust_vbar20( runfit.vbar_vals[ 0 ],
845  tempera );
846  v_buoy[ 0 ] = ( 1.0 - v_vbar[ 0 ] * density );
847  double buoy = v_buoy[ 0 ];
848 
849  for ( int jj = jstx; jj < jstx + v_setpts[ jdx ]; jj++ )
850  {
851  double xv = sq( scnf->xvs[ jj ] ) - xm_sq;
852  double mwv = runfit.mw_vals[ 0 ];
853  double ampv = scnf->amp_vals[ 0 ];
854  double ufunc0 = exp( ampv + dconst * mwv * buoy * xv );
855  double ufunc1 = exp( stoich1 * ampv + runfit.eq_vals[ 0 ] +
856  ODcorec1 + dconst * stoich1 * mwv * buoy * xv );
857  double ufunc2 = exp( stoich2 * ampv + runfit.eq_vals[ 1 ] +
858  ODcorec2 + dconst * stoich2 * mwv * buoy * xv );
859  v_ufunc[ 0 ] = ufunc0;
860  v_ufunc[ 1 ] = ufunc1;
861  v_ufunc[ 2 ] = ufunc2;
862  double dcoeff = dconst * xv * buoy;
863 
864  if ( runfit.mw_fits[ 0 ] )
865  m_jacobi[ jpx ][ runfit.mw_ndxs[ 0 ] ]
866  = dcoeff * ufunc0 * + dcoeff * ufunc1 * stoich1
867  + dcoeff * ufunc2 * stoich2;
868 
869  dcoeff = dconst * mwv * xv * density;
870 
871  if ( runfit.vbar_fits[ 0 ] )
872  m_jacobi[ jpx ][ runfit.vbar_ndxs[ 0 ] ]
873  = (-1.0 ) * dcoeff * ufunc0 - dcoeff * ufunc1 * stoich1
874  - dcoeff * ufunc2 * stoich2;
875 
876  if ( scnf->amp_fits[ 0 ] )
877  m_jacobi[ jpx ][ scnf->amp_ndxs[ 0 ] ]
878  = ufunc0 + ufunc1 * stoich1 + ufunc2 * stoich2;
879 
880  if ( runfit.eq_fits[ 0 ] )
881  m_jacobi[ jpx ][ runfit.eq_ndxs[ 0 ] ] = ufunc1;
882 
883  if ( runfit.eq_fits[ 1 ] )
884  m_jacobi[ jpx ][ runfit.eq_ndxs[ 1 ] ] = ufunc2;
885 
886  if ( scnf->baseln_fit )
887  m_jacobi[ jpx ][ scnf->baseln_ndx ] = 1.0;
888 
889  jpx++;
890  }
891 
892  jdx++;
893  }
894  break;
895  }
896  case 14: // 14: "2-Component Hetero-Association: A + B <=> AB"
897  {
898  double mwv0 = runfit.mw_vals[ 0 ];
899  double mwv1 = runfit.mw_vals[ 1 ];
900  double mw_ab = mwv0 + mwv1;
901 
902  for ( int ii = 0; ii < scanfits.size(); ii++ )
903  {
904  EqScanFit* scnf = &scanfits[ ii ];
905 
906  if ( ! scnf->scanFit ) continue;
907 
908  int jstx = scnf->start_ndx;
909  double xm_sq = sq( scnf->xvs[ jstx ] );
910  double omega_sq = sq( M_PI * scnf->rpm / 30.0 );
911  double tempera = scnf->tempera;
912  double density = scnf->density;
913  double dconst = omega_sq / ( 2.0 * R * ( K0 + tempera ) );
914  double extinc0 = scnf->extincts[ 0 ];
915  double extinc1 = scnf->extincts[ 1 ];
916  double ODcorrec = log( ( extinc0 + extinc1 ) /
917  ( scnf->pathlen * extinc0 * extinc1 ) );
918  v_vbar[ 0 ] = US_Math2::adjust_vbar20( runfit.vbar_vals[ 0 ],
919  tempera );
920  v_vbar[ 1 ] = US_Math2::adjust_vbar20( runfit.vbar_vals[ 1 ],
921  tempera );
922  v_vbar[ 2 ] = ( v_vbar[ 0 ] * mwv0 + v_vbar[ 1 ] * mwv1 )
923  / mw_ab;
924  v_buoy[ 0 ] = ( 1.0 - v_vbar[ 0 ] * density );
925  v_buoy[ 1 ] = ( 1.0 - v_vbar[ 1 ] * density );
926  v_buoy[ 2 ] = ( 1.0 - v_vbar[ 2 ] * density );
927  double buoy0 = v_buoy[ 0 ];
928  double buoy1 = v_buoy[ 1 ];
929  double buoy2 = v_buoy[ 2 ];
930 
931  for ( int jj = jstx; jj < jstx + v_setpts[ jdx ]; jj++ )
932  {
933  double xv = sq( scnf->xvs[ jj ] ) - xm_sq;
934  double ampv0 = scnf->amp_vals[ 0 ];
935  double ampv1 = scnf->amp_vals[ 1 ];
936  double constx = dconst * xv;
937  double ufunc0 = exp( ampv0 + constx * mwv0 * buoy0 );
938  double ufunc1 = exp( ampv1 + constx * mwv1 * buoy1 );
939  double ufunc2 = exp( ampv0 + ampv1 + runfit.eq_vals[ 0 ] +
940  ODcorrec + constx * mw_ab * buoy2 );
941  v_ufunc[ 0 ] = ufunc0;
942  v_ufunc[ 1 ] = ufunc1;
943  v_ufunc[ 2 ] = ufunc2;
944 
945  if ( runfit.mw_fits[ 0 ] )
946  m_jacobi[ jpx ][ runfit.mw_ndxs[ 0 ] ]
947  = constx * buoy0 * ufunc0 +
948  constx * buoy2 +
949  constx * ( v_vbar[ 2 ] * density -
950  v_vbar[ 0 ] * density ) + ufunc2;
951 
952  if ( runfit.vbar_fits[ 0 ] )
953  m_jacobi[ jpx ][ runfit.vbar_ndxs[ 0 ] ]
954  = (-1.0 ) * constx * mwv1 * density * ufunc1
955  - constx * mwv1 * density * ufunc2;
956 
957  if ( scnf->amp_fits[ 0 ] )
958  m_jacobi[ jpx ][ scnf->amp_ndxs[ 0 ] ] = ufunc0 + ufunc2;
959 
960  if ( runfit.eq_fits[ 0 ] )
961  m_jacobi[ jpx ][ runfit.eq_ndxs[ 0 ] ] = ufunc2;
962 
963  if ( scnf->baseln_fit )
964  m_jacobi[ jpx ][ scnf->baseln_ndx ] = 1.0;
965 
966  jpx++;
967  }
968 
969  jdx++;
970  }
971  break;
972  }
973  case 15: // 15: "U-Defined self/Hetero-Assoc.: A + B <=> AB, nA <=> An"
974  {
975  double mwv0 = runfit.mw_vals[ 0 ];
976  double mwv1 = runfit.mw_vals[ 1 ];
977  double mw_ab = mwv0 + mwv1;
978  double stoich1 = runfit.stoichs[ 0 ];
979  double stoiexp = stoich1 - 1.0;
980 
981  for ( int ii = 0; ii < scanfits.size(); ii++ )
982  {
983  EqScanFit* scnf = &scanfits[ ii ];
984 
985  if ( ! scnf->scanFit ) continue;
986 
987  int jstx = scnf->start_ndx;
988  double xm_sq = sq( scnf->xvs[ jstx ] );
989  double omega_sq = sq( M_PI * scnf->rpm / 30.0 );
990  double tempera = scnf->tempera;
991  double density = scnf->density;
992  double dconst = omega_sq / ( 2.0 * R * ( K0 + tempera ) );
993  double extinc0 = scnf->extincts[ 0 ];
994  double extinc1 = scnf->extincts[ 1 ];
995  double ODcorec1 = log( ( extinc0 + extinc0 ) /
996  ( scnf->pathlen * extinc0 * extinc1 ) );
997  double ODcorec2 = log( stoich1 /
998  pow( scnf->pathlen * extinc0, stoiexp ) );
999  v_vbar[ 0 ] = US_Math2::adjust_vbar20( runfit.vbar_vals[ 0 ],
1000  tempera );
1001  v_vbar[ 1 ] = US_Math2::adjust_vbar20( runfit.vbar_vals[ 1 ],
1002  tempera );
1003  v_vbar[ 2 ] = ( v_vbar[ 0 ] * mwv0 + v_vbar[ 1 ] * mwv1 )
1004  / mw_ab;
1005  v_buoy[ 0 ] = ( 1.0 - v_vbar[ 0 ] * density );
1006  v_buoy[ 1 ] = ( 1.0 - v_vbar[ 1 ] * density );
1007  v_buoy[ 2 ] = ( 1.0 - v_vbar[ 2 ] * density );
1008  double buoy0 = v_buoy[ 0 ];
1009  double buoy1 = v_buoy[ 1 ];
1010  double buoy2 = v_buoy[ 2 ];
1011 
1012  for ( int jj = jstx; jj < jstx + v_setpts[ jdx ]; jj++ )
1013  {
1014  double xv = sq( scnf->xvs[ jj ] ) - xm_sq;
1015  double ampv0 = scnf->amp_vals[ 0 ];
1016  double ampv1 = scnf->amp_vals[ 1 ];
1017  double constx = dconst * xv;
1018  double ufunc0 = exp( ampv0 + constx * mwv0 * buoy0 );
1019  double ufunc1 = exp( ampv1 + constx * mwv1 * buoy1 );
1020  double ufunc2 = exp( ampv0 + ampv1 + runfit.eq_vals[ 0 ] +
1021  ODcorec1 + constx * mw_ab * buoy2 );
1022  double ufunc3 = exp( stoich1 * ampv0 + runfit.eq_vals[ 0 ] +
1023  ODcorec2 + constx * stoich1 * mw_ab * buoy0 );
1024  v_ufunc[ 0 ] = ufunc0;
1025  v_ufunc[ 1 ] = ufunc1;
1026  v_ufunc[ 2 ] = ufunc2;
1027  v_ufunc[ 3 ] = ufunc3;
1028 
1029  if ( runfit.mw_fits[ 0 ] )
1030  m_jacobi[ jpx ][ runfit.mw_ndxs[ 0 ] ]
1031  = constx * buoy0 * ufunc0 +
1032  constx * buoy2 +
1033  constx * ( v_vbar[ 2 ] * density -
1034  v_vbar[ 0 ] * density ) +
1035  constx * buoy0 * ufunc3 * stoich1 + ufunc2;
1036 
1037  if ( runfit.mw_fits[ 1 ] )
1038  m_jacobi[ jpx ][ runfit.mw_ndxs[ 1 ] ]
1039  = constx * buoy1 * ufunc1 +
1040  constx * buoy2 +
1041  constx * ( v_vbar[ 2 ] * density -
1042  v_vbar[ 1 ] * density ) + ufunc2;
1043 
1044  if ( runfit.vbar_fits[ 0 ] )
1045  m_jacobi[ jpx ][ runfit.vbar_ndxs[ 0 ] ]
1046  = (-1.0 ) * constx * mwv0 * density * ufunc0
1047  - constx * mwv0 * density * ufunc3 * stoich1
1048  - constx * mwv0 * density * ufunc2;
1049 
1050  if ( runfit.vbar_fits[ 1 ] )
1051  m_jacobi[ jpx ][ runfit.vbar_ndxs[ 1 ] ]
1052  = (-1.0 ) * constx * mwv1 * density * ufunc1
1053  - constx * mwv1 * density * ufunc2;
1054 
1055  if ( scnf->amp_fits[ 0 ] )
1056  m_jacobi[ jpx ][ scnf->amp_ndxs[ 0 ] ]
1057  = ufunc0 + ufunc3 * stoich1 + ufunc2;
1058 
1059  if ( scnf->amp_fits[ 1 ] )
1060  m_jacobi[ jpx ][ scnf->amp_ndxs[ 1 ] ] = ufunc1 + ufunc2;
1061 
1062  if ( runfit.eq_fits[ 0 ] )
1063  m_jacobi[ jpx ][ runfit.eq_ndxs[ 0 ] ] = ufunc2;
1064 
1065  if ( runfit.eq_fits[ 1 ] )
1066  m_jacobi[ jpx ][ runfit.eq_ndxs[ 1 ] ] = ufunc3;
1067 
1068  if ( scnf->baseln_fit )
1069  m_jacobi[ jpx ][ scnf->baseln_ndx ] = 1.0;
1070 
1071  jpx++;
1072  }
1073 
1074  jdx++;
1075  }
1076  break;
1077  }
1078  case 16: // 16: "U-Defined Monomer-Nmer, some monomer is incompetent"
1079  {
1080  double mwv0 = runfit.mw_vals[ 0 ];
1081  double mwv1 = runfit.mw_vals[ 1 ];
1082  double stoich1 = runfit.stoichs[ 0 ];
1083  double stoiexp = stoich1 - 1.0;
1084 
1085  for ( int ii = 0; ii < scanfits.size(); ii++ )
1086  {
1087  EqScanFit* scnf = &scanfits[ ii ];
1088 
1089  if ( ! scnf->scanFit ) continue;
1090 
1091  int jstx = scnf->start_ndx;
1092  double xm_sq = sq( scnf->xvs[ jstx ] );
1093  double omega_sq = sq( M_PI * scnf->rpm / 30.0 );
1094  double tempera = scnf->tempera;
1095  double density = scnf->density;
1096  double dconst = omega_sq / ( 2.0 * R * ( K0 + tempera ) );
1097  double extinc0 = scnf->extincts[ 0 ];
1098  double ODcorec1 = log( stoich1 /
1099  pow( extinc0 * scnf->pathlen, stoiexp ) );
1100  v_vbar[ 0 ] = US_Math2::adjust_vbar20( runfit.vbar_vals[ 0 ],
1101  tempera );
1102  v_buoy[ 0 ] = ( 1.0 - v_vbar[ 0 ] * density );
1103  double buoy0 = v_buoy[ 0 ];
1104 
1105  for ( int jj = jstx; jj < jstx + v_setpts[ jdx ]; jj++ )
1106  {
1107  double xv = sq( scnf->xvs[ jj ] ) - xm_sq;
1108  double ampv0 = scnf->amp_vals[ 0 ];
1109  double ampv1 = scnf->amp_vals[ 1 ];
1110  double constx = dconst * xv;
1111  double ufunc0 = exp( ampv0 + constx * mwv0 * buoy0 );
1112  double ufunc1 = exp( ampv1 + constx * mwv1 * buoy0 );
1113  double ufunc2 = exp( stoich1 * ampv0 + runfit.eq_vals[ 0 ] +
1114  ODcorec1 + constx * stoich1 * mwv0 * buoy0 );
1115  v_ufunc[ 0 ] = ufunc0;
1116  v_ufunc[ 1 ] = ufunc1;
1117  v_ufunc[ 2 ] = ufunc2;
1118 
1119  if ( runfit.mw_fits[ 0 ] )
1120  m_jacobi[ jpx ][ runfit.mw_ndxs[ 0 ] ]
1121  = constx * buoy0 * ufunc0 +
1122  constx * buoy0 * ufunc1 +
1123  constx * buoy0 * ufunc2 * stoich1;
1124 
1125  if ( runfit.vbar_fits[ 0 ] )
1126  m_jacobi[ jpx ][ runfit.vbar_ndxs[ 0 ] ]
1127  = (-1.0 ) * constx * mwv0 * density * ufunc0
1128  - constx * mwv0 * density * ufunc1
1129  - constx * mwv0 * density * ufunc2 * stoich1;
1130 
1131  if ( scnf->amp_fits[ 0 ] )
1132  m_jacobi[ jpx ][ scnf->amp_ndxs[ 0 ] ]
1133  = ufunc0 + stoich1 * ufunc2;
1134 
1135  if ( scnf->amp_fits[ 1 ] )
1136  m_jacobi[ jpx ][ scnf->amp_ndxs[ 1 ] ] = ufunc1;
1137 
1138  if ( runfit.eq_fits[ 0 ] )
1139  m_jacobi[ jpx ][ runfit.eq_ndxs[ 0 ] ] = ufunc2;
1140 
1141  if ( scnf->baseln_fit )
1142  m_jacobi[ jpx ][ scnf->baseln_ndx ] = 1.0;
1143 
1144  jpx++;
1145  }
1146 
1147  jdx++;
1148  }
1149  break;
1150  }
1151  case 17: // 17: "User-Defined Monomer-Nmer, some Nmer is incompetent"
1152  {
1153  double mwv0 = runfit.mw_vals[ 0 ];
1154  double mwv1 = runfit.mw_vals[ 1 ];
1155  double stoich1 = runfit.stoichs[ 0 ];
1156  double stoiexp = stoich1 - 1.0;
1157 
1158  for ( int ii = 0; ii < scanfits.size(); ii++ )
1159  {
1160  EqScanFit* scnf = &scanfits[ ii ];
1161 
1162  if ( ! scnf->scanFit ) continue;
1163 
1164  int jstx = scnf->start_ndx;
1165  double xm_sq = sq( scnf->xvs[ jstx ] );
1166  double omega_sq = sq( M_PI * scnf->rpm / 30.0 );
1167  double tempera = scnf->tempera;
1168  double density = scnf->density;
1169  double dconst = omega_sq / ( 2.0 * R * ( K0 + tempera ) );
1170  double extinc0 = scnf->extincts[ 0 ];
1171  double ODcorec1 = log( stoich1 /
1172  pow( scnf->pathlen * extinc0, stoiexp ) );
1173  v_vbar[ 0 ] = US_Math2::adjust_vbar20( runfit.vbar_vals[ 0 ],
1174  tempera );
1175  v_buoy[ 0 ] = ( 1.0 - v_vbar[ 0 ] * density );
1176  double buoy0 = v_buoy[ 0 ];
1177 
1178  for ( int jj = jstx; jj < jstx + v_setpts[ jdx ]; jj++ )
1179  {
1180  double xv = sq( scnf->xvs[ jj ] ) - xm_sq;
1181  double ampv0 = scnf->amp_vals[ 0 ];
1182  double ampv1 = scnf->amp_vals[ 1 ];
1183  double constx = dconst * xv;
1184  double ufunc0 = exp( ampv0 + constx * mwv0 * buoy0 );
1185  double ufunc1 = exp( ampv1 + constx * mwv1 * buoy0 );
1186  double ufunc2 = exp( stoich1 * ampv0 + runfit.eq_vals[ 0 ] +
1187  ODcorec1 + constx * stoich1 * mwv0 * buoy0 );
1188  v_ufunc[ 0 ] = ufunc0;
1189  v_ufunc[ 1 ] = ufunc1;
1190  v_ufunc[ 2 ] = ufunc2;
1191 
1192  if ( runfit.mw_fits[ 0 ] )
1193  m_jacobi[ jpx ][ runfit.mw_ndxs[ 0 ] ]
1194  = constx * buoy0 * ufunc0 +
1195  constx * buoy0 * ufunc1 * stoich1 +
1196  constx * buoy0 * ufunc2 * stoich1;
1197 
1198  if ( runfit.vbar_fits[ 0 ] )
1199  m_jacobi[ jpx ][ runfit.vbar_ndxs[ 1 ] ]
1200  = (-1.0 ) * constx * mwv0 * density * ufunc0
1201  - constx * mwv0 * density * ufunc1 * stoich1
1202  - constx * mwv0 * density * ufunc2 * stoich1;
1203 
1204  if ( scnf->amp_fits[ 0 ] )
1205  m_jacobi[ jpx ][ scnf->amp_ndxs[ 0 ] ]
1206  = ufunc0 + ufunc2 * stoich1;
1207 
1208  if ( scnf->amp_fits[ 1 ] )
1209  m_jacobi[ jpx ][ scnf->amp_ndxs[ 1 ] ] = ufunc1;
1210 
1211  if ( runfit.eq_fits[ 0 ] )
1212  m_jacobi[ jpx ][ runfit.eq_ndxs[ 0 ] ] = ufunc2;
1213 
1214  if ( scnf->baseln_fit )
1215  m_jacobi[ jpx ][ scnf->baseln_ndx ] = 1.0;
1216 
1217  jpx++;
1218  }
1219 
1220  jdx++;
1221  }
1222  break;
1223  }
1224  case 18: // 18: "User-Defined irreversible Monomer-Nmer"
1225  {
1226  double mwv0 = runfit.mw_vals[ 0 ];
1227  double stoich1 = runfit.stoichs[ 0 ];
1228  double mwv1 = mwv0 * stoich1;
1229 
1230  for ( int ii = 0; ii < scanfits.size(); ii++ )
1231  {
1232  EqScanFit* scnf = &scanfits[ ii ];
1233 
1234  if ( ! scnf->scanFit ) continue;
1235 
1236  int jstx = scnf->start_ndx;
1237  double xm_sq = sq( scnf->xvs[ jstx ] );
1238  double omega_sq = sq( M_PI * scnf->rpm / 30.0 );
1239  double tempera = scnf->tempera;
1240  double density = scnf->density;
1241  double dconst = omega_sq / ( 2.0 * R * ( K0 + tempera ) );
1242  v_vbar[ 0 ] = US_Math2::adjust_vbar20( runfit.vbar_vals[ 0 ],
1243  tempera );
1244  v_buoy[ 0 ] = ( 1.0 - v_vbar[ 0 ] * density );
1245  double buoy0 = v_buoy[ 0 ];
1246  double ampv0 = scnf->amp_vals[ 0 ];
1247  double ampv1 = scnf->amp_vals[ 1 ];
1248 
1249  for ( int jj = jstx; jj < jstx + v_setpts[ jdx ]; jj++ )
1250  {
1251  double xv = sq( scnf->xvs[ jj ] ) - xm_sq;
1252  double constx = dconst * xv;
1253  double ufunc0 = exp( ampv0 + constx * mwv0 * buoy0 );
1254  double ufunc1 = exp( ampv1 + constx * mwv1 * buoy0 );
1255  v_ufunc[ 0 ] = ufunc0;
1256  v_ufunc[ 1 ] = ufunc1;
1257 
1258  if ( runfit.mw_fits[ 0 ] )
1259  m_jacobi[ jpx ][ runfit.mw_ndxs[ 0 ] ]
1260  = constx * buoy0 * ufunc0 +
1261  constx * buoy0 * ufunc1 * stoich1;
1262 
1263  if ( runfit.vbar_fits[ 0 ] )
1264  m_jacobi[ jpx ][ runfit.vbar_ndxs[ 0 ] ]
1265  = (-1.0 ) * constx * mwv0 * density * ufunc0
1266  - constx * mwv0 * density * ufunc1 * stoich1;
1267 
1268  if ( scnf->amp_fits[ 0 ] )
1269  m_jacobi[ jpx ][ scnf->amp_ndxs[ 0 ] ] = ufunc0;
1270 
1271  if ( scnf->amp_fits[ 1 ] )
1272  m_jacobi[ jpx ][ scnf->amp_ndxs[ 1 ] ] = ufunc1;
1273 
1274  if ( scnf->baseln_fit )
1275  m_jacobi[ jpx ][ scnf->baseln_ndx ] = 1.0;
1276 
1277  jpx++;
1278  }
1279 
1280  jdx++;
1281  }
1282  break;
1283  }
1284  case 19: // 19: "User-Defined Monomer-Nmer plus contaminant"
1285  {
1286  double mwv0 = runfit.mw_vals[ 0 ];
1287  double mwv1 = runfit.mw_vals[ 1 ];
1288  double stoich1 = runfit.stoichs[ 0 ];
1289  double stoiexp = stoich1 - 1.0;
1290 
1291  for ( int ii = 0; ii < scanfits.size(); ii++ )
1292  {
1293  EqScanFit* scnf = &scanfits[ ii ];
1294 
1295  if ( ! scnf->scanFit ) continue;
1296 
1297  int jstx = scnf->start_ndx;
1298  double xm_sq = sq( scnf->xvs[ jstx ] );
1299  double omega_sq = sq( M_PI * scnf->rpm / 30.0 );
1300  double tempera = scnf->tempera;
1301  double density = scnf->density;
1302  double dconst = omega_sq / ( 2.0 * R * ( K0 + tempera ) );
1303  double extinc0 = scnf->extincts[ 0 ];
1304  double ODcorec1 = log( stoich1 /
1305  pow( scnf->pathlen * extinc0, stoiexp ) );
1306  v_vbar[ 0 ] = US_Math2::adjust_vbar20( runfit.vbar_vals[ 0 ],
1307  tempera );
1308  v_vbar[ 1 ] = US_Math2::adjust_vbar20( runfit.vbar_vals[ 1 ],
1309  tempera );
1310  v_buoy[ 0 ] = ( 1.0 - v_vbar[ 0 ] * density );
1311  v_buoy[ 1 ] = ( 1.0 - v_vbar[ 1 ] * density );
1312  double buoy0 = v_buoy[ 0 ];
1313  double buoy1 = v_buoy[ 1 ];
1314  double ampv0 = scnf->amp_vals[ 0 ];
1315  double ampv1 = scnf->amp_vals[ 1 ];
1316 
1317  for ( int jj = jstx; jj < jstx + v_setpts[ jdx ]; jj++ )
1318  {
1319  double xv = sq( scnf->xvs[ jj ] ) - xm_sq;
1320  double constx = dconst * xv;
1321  double ufunc0 = exp( ampv0 + constx * mwv0 * buoy0 );
1322  double ufunc1 = exp( ampv1 + constx * mwv1 * buoy1 );
1323  double ufunc2 = exp( stoich1 * ampv0 + runfit.eq_vals[ 0 ] +
1324  ODcorec1 + constx * stoich1 * mwv0 * buoy0 );
1325  v_ufunc[ 0 ] = ufunc0;
1326  v_ufunc[ 1 ] = ufunc1;
1327  v_ufunc[ 2 ] = ufunc2;
1328 
1329  if ( runfit.mw_fits[ 0 ] )
1330  m_jacobi[ jpx ][ runfit.mw_ndxs[ 0 ] ]
1331  = constx * buoy0 * ufunc0 +
1332  constx * buoy0 * ufunc2 * stoich1;
1333 
1334  if ( runfit.mw_fits[ 1 ] )
1335  m_jacobi[ jpx ][ runfit.mw_ndxs[ 1 ] ]
1336  = constx * buoy1 * ufunc1;
1337 
1338  if ( runfit.vbar_fits[ 0 ] )
1339  m_jacobi[ jpx ][ runfit.vbar_ndxs[ 0 ] ]
1340  = (-1.0 ) * constx * mwv0 * density * ufunc0
1341  - constx * mwv0 * density * ufunc2;
1342 
1343  if ( scnf->amp_fits[ 0 ] )
1344  m_jacobi[ jpx ][ scnf->amp_ndxs[ 0 ] ]
1345  = ufunc0 + ufunc2 * stoich1;
1346 
1347  if ( scnf->amp_fits[ 1 ] )
1348  m_jacobi[ jpx ][ scnf->amp_ndxs[ 1 ] ] = ufunc1;
1349 
1350  if ( runfit.eq_fits[ 0 ] )
1351  m_jacobi[ jpx ][ runfit.eq_ndxs[ 0 ] ] = ufunc2;
1352 
1353  if ( scnf->baseln_fit )
1354  m_jacobi[ jpx ][ scnf->baseln_ndx ] = 1.0;
1355 
1356  jpx++;
1357  }
1358 
1359  jdx++;
1360  }
1361  break;
1362  }
1363  }
1364 
1365  return stat;
1366 }
1367 
1368 // Calculate the chi-squared for the fixed molecular weight estimate
1369 // for a single component model
1370 double US_EqMath::calc_testParameter( double mwval )
1371 {
1372  int points = 0;
1373  double chi_sq = 0.0;
1374  double x_temp = 0.0;
1375  double dconst;
1376  double buoyancy_tb;
1377  double omega_s;
1378  double x_val;
1379  double y_val;
1380 
1381  QVector< double* > mmat; // M matrix (vector of pointers to arrays)
1382  QVector< double > mvec; // Concatenated matrix doubles
1383  QVector< double > yrvec; // Raw Y-values vector
1384  QVector< double > covec; // Coefficients vector
1385  QVector< QPointF > xyvec; // X,Y vector (of points)
1386 
1387  // Determine maximum number of points for vectors
1388  for ( int jes = 0; jes < scanfits.size(); jes++ )
1389  points = max( points, scanfits[ jes ].xvs.size() );
1390 
1391 DbgLv(1) << "ctPar: max points" << points;
1392  mmat .fill( NULL, points ); // Initialize matrix (array pointers)
1393  mvec .fill( 0.0, points * 2 ); // Initialize vector of matrix elements
1394  yrvec.fill( 0.0, points ); // Initialize vector of raw Y's
1395  covec.fill( 0.0, 2 );
1396  xyvec.fill( QPointF( 0.0, 0.0 ), points );
1397 
1398  for ( int jes = 0; jes < scanfits.size(); jes++ )
1399  { // Build up chi-squared for fitted scans
1400  EqScanFit* scnf = &scanfits[ jes ];
1401 
1402  if ( ! scnf->scanFit )
1403  continue; // Ignore any scan that is not fitted
1404 
1405  int strtx = scnf->start_ndx; // Get scan parameters
1406  int stopx = scnf->stop_ndx;
1407  int stopn = stopx + 1;
1408  double fstemp = scnf->tempera;
1409  double fsdens = scnf->density;
1410  double fsvisc = scnf->viscosity;
1411  double fsrpm = scnf->rpm;
1412  double fsvbar = runfit.vbar_vals[ 0 ];
1413 //DbgLv(1) << "ctPar: jes strtx stopx" << jes << strtx << stopx;
1414 
1415  US_Math2::SolutionData solution; // Calculate constant based on
1416  solution.vbar = fsvbar; // buoyancy, based on vbar,density...
1417  solution.vbar20 = fsvbar;
1418  solution.density = fsdens;
1419  solution.viscosity = fsvisc;
1420  US_Math2::data_correction( fstemp, solution );
1421  buoyancy_tb = solution.buoyancyb;
1422 
1423  omega_s = sq( fsrpm * M_PI / 30.0 );
1424  dconst = ( buoyancy_tb * omega_s ) / ( 2.0 * R * ( K0 + fstemp ) );
1425  x_temp = scnf->xvs[ strtx ];
1426  x_temp = sq( x_temp );
1427 //DbgLv(1) << "ctPar: xt x0 xN" << x_temp << scnf->xvs[strtx] << scnf->xvs[stopx];
1428  points = 0;
1429 
1430  for ( int jxy = strtx; jxy < stopn; jxy++ )
1431  { // Accumulate the X,Y values within the scan's start,stop range
1432  x_val = scnf->xvs[ jxy ];
1433  y_val = scnf->yvs[ jxy ];
1434  x_val = dconst * ( sq( x_val ) - x_temp );
1435  xyvec[ points++ ] = QPointF( x_val, y_val );
1436  }
1437 
1438  double** MM = mmat .data(); // Pointer to matrix data
1439  double* yraw = yrvec.data(); // Pointer to raw-Y data
1440  double* mvls = mvec .data(); // Pointer to matrix values
1441  double* coeffs = covec.data(); // Pointer to output coefficients
1442  int kk = 0;
1443 
1444  for ( int ii = 0; ii < points; ii++ )
1445  { // Build matrix and vector needed for LS
1446  x_val = xyvec[ ii ].x(); // Raw X value
1447  y_val = xyvec[ ii ].y(); // Raw Y value
1448 
1449  // The following builds and stores MM[ii][0] and MM[ii][1]
1450  mmat[ ii ] = mvls + kk; // Store pointer to 2-point array
1451  mvec[ kk++ ] = 1.0; // Update matrix column
1452  mvec[ kk++ ] = exp( mwval * x_val );
1453 
1454  yrvec[ ii ] = y_val; // Update Y-Raw vector
1455  }
1456 //DbgLv(1) << "ctPar: yrvec[0]" << yrvec[0] << "points" << points;
1457 //DbgLv(1) << "ctPar: mm00 mm01" << MM[0][0] << MM[0][1];
1458 //int N=points-1;
1459 //DbgLv(1) << "ctPar: mmN0 mN01" << MM[N][0] << MM[N][1];
1460 //DbgLv(1) << "ctPar: mwval x_val" << mwval << x_val;
1461 
1462  // Get coefficients thru general Least Squares, order=2
1463  genLeastSquaresOrd2( MM, points, yraw, &coeffs );
1464 
1465  for ( int ii = 0; ii < points; ii++ )
1466  { // Update the chi-squared value
1467  double chi = yrvec[ ii ] - ( coeffs[ 0 ] +
1468  coeffs[ 1 ] * exp( mwval * xyvec[ ii ].x() ) );
1469  chi_sq += sq( chi );
1470  }
1471 //DbgLv(1) << "ctPar: chi_sq" << chi_sq;
1472 
1473  scnf->baseline = coeffs[ 0 ];
1474  scnf->baseln_rng = 0.05;
1475 
1476  if ( coeffs[ 1 ] < dflt_min || isNan( coeffs[ 1 ] ) )
1477  {
1478 DbgLv(1) << "ctPar: ** jes" << jes << "coeffs1 " << covec[1];
1479  scnf->amp_vals[ 0 ] = dflt_min;
1480  }
1481 
1482  else
1483  {
1484  scnf->amp_vals[ 0 ] = log( coeffs[ 1 ] );
1485  }
1486 
1487  scnf->amp_rngs[ 0 ] = scnf->amp_vals[ 0 ] * 0.2;
1488  }
1489 DbgLv(1) << "ctPar: coeffs 0 1 " << covec[0] << covec[1];
1490 DbgLv(1) << "ctPar:chi_sq" << chi_sq;
1491 
1492  return chi_sq;
1493 }
1494 
1495 // Find the minimum residual.
1496 // Residual values are f0, f1, f2; evaluated at x0, x1, x2.
1497 // x0, x1, x2 are multipliers for incremental change of the parameter.
1498 // Calculate bracket: assume the minimum is between x0=0 and some stepsize
1499 // x2 away from x0. Then find an x1 in the middle between x0 and x2; and
1500 // calculate f1(x1).
1501 double US_EqMath::linesearch()
1502 {
1503  double residm = 0.0;
1504  double old_f0 = 0.0;
1505  double old_f1 = 0.0;
1506  double old_f2 = 0.0;
1507  double x0 = 100.0;
1508  double x1 = 5000.0;
1509  double x2 = 10000.0;
1510  double hh = 0.01;
1511  double toler = 100.0;
1512  double errmx;
1513  double xmin;
1514  double fmin;
1515  int iter = 1;
1516 
1517  double f0 = calc_testParameter( x0 );
1518  double f1 = calc_testParameter( x1 );
1519  double f2 = calc_testParameter( x2 );
1520 DbgLv(1) << "MinRes:iter" << iter << " f0 f1 f2" << f0 << f1 << f2;
1521 DbgLv(1) << "MinRes: dflt_min dflt_max" << dflt_min << dflt_max;
1522 
1523  if ( dataList[ 0 ].dataType == "RA" )
1524  errmx = 1.0e4; // absorbance
1525  else
1526  errmx = 1.0e12; // larger value for interference
1527 
1528 
1529  while ( f0 >= errmx || f0 < 0.0 ||
1530  f1 >= errmx || f1 < 0.0 ||
1531  f2 >= errmx || f2 < 0.0 )
1532  { // Assure f1,f2 are between 0.0 and errmx
1533  x1 *= 0.5;
1534  x2 *= 0.5;
1535 
1536  f1 = calc_testParameter( x1 );
1537  f2 = calc_testParameter( x2 );
1538 
1539  if ( x1 < dflt_min )
1540  return ( -1.0 );
1541  }
1542 DbgLv(1) << " x1 x2" << x1 << x2 << " f1 f2" << f1 << f2 << "errmx" << errmx;
1543 
1544  bool check_flag = true;
1545 
1546  while ( check_flag )
1547  {
1548  if ( ( isNan( f0 ) && isNan( f1 ) ) ||
1549  ( isNan( f1 ) && isNan( f2 ) ) ||
1550  ( isNan( f0 ) && isNan( f1 ) ) )
1551  return ( -1.0 ); // At least two values screwed up
1552 
1553  if ( ( qAbs( f2 - old_f2 ) < dflt_min ) &&
1554  ( qAbs( f1 - old_f1 ) < dflt_min ) &&
1555  ( qAbs( f0 - old_f0 ) < dflt_min ) )
1556  return ( 0.0 ); // Solution converged horizontally
1557 
1558  old_f0 = f0;
1559  old_f1 = f1;
1560  old_f2 = f2;
1561 DbgLv(1) << " old f0 f1 f2" << f0 << f1 << f2;
1562 
1563 DbgLv(1) << " test-0a";
1564 bool t1 = ( qAbs( f2 - f0 ) < dflt_min );
1565 bool t2 = ( qAbs( f1 - f0 ) < dflt_min );
1566 bool t3 = ( f0 > f1 );
1567 bool t4 = ( qAbs( f2 - f1 ) < dflt_min );
1568 bool t12 = t1 && t2;
1569 bool t34 = t3 && t4;
1570 DbgLv(1) << " t-0a t1-4" << t1 << t2 << t3 << t4 << "t12,34" << t12 << t34;
1571  if ( ( ( qAbs( f2 - f0 ) < dflt_min ) &&
1572  ( qAbs( f1 - f0 ) < dflt_min ) ) ||
1573  ( ( f0 > f1 ) &&
1574  ( qAbs( f2 - f1 ) < dflt_min ) ) )
1575  return ( 0.0 );
1576 
1577 DbgLv(1) << " test-0b";
1578  if ( ( qAbs( x0 ) < dflt_min ) &&
1579  ( qAbs( x1 ) < dflt_min ) &&
1580  ( qAbs( x2 ) < dflt_min ) )
1581  return ( 0.0 );
1582 
1583 DbgLv(1) << " test-0c";
1584  if ( ( ( qAbs( f0 - f1 ) < dflt_min ) &&
1585  ( qAbs( f1 - f2 ) < dflt_min ) ) ||
1586  ( ( f2 > f1 ) &&
1587  ( qAbs( f0 - f1 ) < dflt_min ) ) )
1588  return ( 0.0 );
1589 
1590 DbgLv(1) << " test-0 x0 x1 x2" << x0 << x1 << x2;
1591  if ( ( f0 > f1 ) && ( f2 > f1 ) ) // We have a bracket
1592  {
1593 DbgLv(1) << " update-0 f0 f1 f2" << f0 << f1 << f2;
1594  check_flag = false;
1595  break;
1596  }
1597 
1598  if ( ( f2 > f1 && f1 > f0 ) ||
1599  ( f1 > f0 && f1 > f2 ) ||
1600  ( f1 == f2 && f1 > f0 ) )
1601  {
1602  x2 = x1;
1603  f2 = f1;
1604  x1 = ( x2 + x1 ) * 0.5;
1605  f1 = calc_testParameter( x1 );
1606 DbgLv(1) << " update-1 f0 f1 f2" << f0 << f1 << f2;
1607  }
1608 
1609  else if ( ( f0 > f1 ) && ( f1 > f2 ) )
1610  {
1611  x0 = x1;
1612  f0 = f1;
1613  x1 = x2;
1614  f1 = f2;
1615  x2 = x2 + ( ( pow( 2.0, (double)( iter + 2 ) ) ) + hh );
1616  f2 = calc_testParameter( x2 );
1617 DbgLv(1) << " update-2 f0 f1 f2" << f0 << f1 << f2;
1618  }
1619 
1620  iter++;
1621 DbgLv(1) << "MinRes: iter" << iter << " f0 f1 f2" << f0 << f1 << f2;
1622  }
1623 
1624 DbgLv(1) << "MinRes: iter" << iter << " f0 f1 f2" << f0 << f1 << f2;
1625  x1 = ( x0 + x2 ) * 0.5;
1626  hh = x1 - x0;
1627  f1 = calc_testParameter( x1 );
1628 
1629  while ( true )
1630  {
1631  if ( f0 < f1 ) // Shift left
1632  {
1633  x2 = x1;
1634  f2 = f1;
1635  x1 = x0;
1636  f1 = f0;
1637  x0 = x1 - hh;
1638  f0 = calc_testParameter( x0 );
1639  }
1640 
1641  if ( f2 < f1 ) // Shift right
1642  {
1643  x0 = x1;
1644  f0 = f1;
1645  x1 = x2;
1646  f1 = f2;
1647  x2 = x1 + hh;
1648  f2 = calc_testParameter( x2 );
1649  }
1650 
1651  if ( qAbs( f0 - ( f1 * 2.0 ) + f2 ) < dflt_min )
1652  {
1653  residm = 0.0;
1654  break;
1655  }
1656 
1657  xmin = x1 + ( hh * ( f0 - f2 ) ) / ( 2.0 * ( f0 - f1 * 2.0 + f2 ) );
1658  fmin = calc_testParameter( xmin );
1659 
1660  if ( fmin < f1 )
1661  {
1662  x1 = xmin;
1663  f1 = fmin;
1664  }
1665 
1666  hh = hh * 0.5;
1667 
1668  if ( hh < toler )
1669  {
1670  residm = x1; // Shouldn't it be "f" we return?
1671  //residm = f1;
1672  break;
1673  }
1674 
1675  x0 = x1 - hh;
1676  x2 = x1 + hh;
1677  f0 = calc_testParameter( x0 );
1678  f2 = calc_testParameter( x2 );
1679 DbgLv(1) << " update-9 x0 x2 f0 f2" << x0 << x2 << f0 << f2;
1680  }
1681 
1682  return residm;
1683 }
1684 
1685 // Calculate the B matrix ( = Jacobian-tranpose * ydelta )
1686 void US_EqMath::calc_B()
1687 {
1688  // Calculate B = J' * d
1689  US_Matrix::tvv( jacobian, y_delta, BB, ntpts, nfpars );
1690 }
1691 
1692 // Calculate the delta array and the variance value
1693 double US_EqMath::calc_residuals()
1694 {
1695  double residual = 0.0;
1696 
1697  for ( int ii = 0; ii < ntpts; ii++ )
1698  {
1699  double delta = y_raw[ ii ] - y_guess[ ii ];
1700 
1701  y_delta[ ii ] = delta;
1702  residual += sq( delta );
1703  }
1704 
1705  if ( residual > dflt_max )
1706  residual = -1.0;
1707 
1708  else
1709  residual /= (double)ntpts;
1710 
1711 if(residual<0.0)
1712  qDebug() << "CalcResid: resid" << residual << "delta 0,1,m,n"
1713  << y_delta[0] << y_delta[1] << y_delta[ntpts-2] << y_delta[ntpts-1];
1714 
1715  return residual;
1716 }
1717 
1718 // Calculate the model from a set of guesses
1719 int US_EqMath::calc_model( double* guess )
1720 {
1721  int stat = 0;
1722 
1723  parameter_mapBackward( guess );
1724 
1725  int ncomps = runfit.nbr_comps;
1726  int jpx = 0; // Total points count
1727  int jlx = 0; // Log count
1728  int jdx = 0; // Datasets count
1729  int mcomp = max( ncomps, 4 );
1730  QVector< double > v_vbar ( mcomp );
1731  QVector< double > v_buoy ( mcomp );
1732 
1733  v_setlpts.clear();
1734 
1735  switch( modelx )
1736  {
1737  case 0: // 0: "1-Component, Ideal"
1738  case 1: // 1: "2-Component, Ideal, Noninteracting"
1739  case 2: // 2: "3-Component, Ideal, Noninteracting"
1740  case 3: // 3: "Fixed Molecular Weight Distribution"
1741  for ( int ii = 0; ii < scanfits.size(); ii++ )
1742  {
1743  EqScanFit* scnf = &scanfits[ ii ];
1744 
1745  if ( ! scnf->scanFit ) continue;
1746 
1747  int jstx = scnf->start_ndx;
1748  double xm_sq = sq( scnf->xvs[ jstx ] );
1749  double omega_sq = sq( M_PI * scnf->rpm / 30.0 );
1750  double tempera = scnf->tempera;
1751  double density = scnf->density;
1752  double dconst = omega_sq / ( 2.0 * R * ( K0 + tempera ) );
1753 
1754  for ( int kk = 0; kk < ncomps; kk++ )
1755  {
1756  v_vbar[ kk ] = US_Math2::adjust_vbar20( runfit.vbar_vals[ kk ],
1757  tempera );
1758  v_buoy[ kk ] = ( 1.0 - v_vbar[ kk ] * density );
1759  }
1760 
1761  jlx = 0;
1762 
1763  for ( int jj = jstx; jj < jstx + v_setpts[ jdx ]; jj++ )
1764  {
1765  double xv = sq( scnf->xvs[ jj ] ) - xm_sq;
1766  double cguess = 0.0;
1767 
1768  for ( int kk = 0; kk < ncomps; kk++ )
1769  {
1770  double buoy = v_buoy[ kk ];
1771  double mwv = runfit.mw_vals[ kk ];
1772  double ampv = scnf->amp_vals[ kk ];
1773  double ufunc = exp( ampv + dconst * mwv * buoy * xv );
1774  cguess += ufunc;
1775  }
1776 
1777  y_guess[ jpx ] = cguess + scnf->baseline;
1778 
1779  if ( cguess > 0.0 )
1780  {
1781  lncr2[ jdx ][ jlx++ ] = log( cguess );
1782  }
1783 
1784  jpx++;
1785  }
1786 
1787  v_setlpts << jlx;
1788  jdx++;
1789  }
1790  break;
1791  case 4: // 4: "Monomer-Dimer Equilibrium"
1792  case 5: // 5: "Monomer-Trimer Equilibrium"
1793  case 6: // 6: "Monomer-Tetramer Equilibrium"
1794  case 7: // 7: "Monomer-Pentamer Equilibrium"
1795  case 8: // 8: "Monomer-Hexamer Equilibrium"
1796  case 9: // 9: "Monomer-Heptamer Equilibrium"
1797  case 10: // 10: "User-Defined Monomer-Nmer Equilibrium"
1798  {
1799  double stoich1 = (double)( modelx - 2 );
1800  double stoiexp = stoich1 - 1.0;
1801  runfit.stoichs[ 0 ] = stoich1;
1802 
1803  for ( int ii = 0; ii < scanfits.size(); ii++ )
1804  {
1805  EqScanFit* scnf = &scanfits[ ii ];
1806 
1807  if ( ! scnf->scanFit ) continue;
1808 
1809  int jstx = scnf->start_ndx;
1810  double xm_sq = sq( scnf->xvs[ jstx ] );
1811  double omega_sq = sq( M_PI * scnf->rpm / 30.0 );
1812  double tempera = scnf->tempera;
1813  double density = scnf->density;
1814  double dconst = omega_sq / ( 2.0 * R * ( K0 + tempera ) );
1815  double ODcorrec = log( stoich1 / pow( scnf->extincts[ 0 ]
1816  * scnf->pathlen , stoiexp ) );
1817  v_vbar[ 0 ] = US_Math2::adjust_vbar20( runfit.vbar_vals[ 0 ],
1818  tempera );
1819  double buoy = ( 1.0 - v_vbar[ 0 ] * density );
1820  jlx = 0;
1821 
1822  for ( int jj = jstx; jj < jstx + v_setpts[ jdx ]; jj++ )
1823  {
1824  double xv = sq( scnf->xvs[ jj ] ) - xm_sq;
1825  double mwv = runfit.mw_vals[ 0 ];
1826  double ampv = scnf->amp_vals[ 0 ];
1827  double ufunc0 = exp( ampv + dconst * mwv * buoy * xv );
1828  double ufunc1 = exp( stoich1 * ampv + runfit.eq_vals[ 0 ] +
1829  ODcorrec + dconst * stoich1 * mwv * buoy * ufunc0 );
1830  double cguess = ufunc0 + ufunc1;
1831 
1832  y_guess[ jpx ] = cguess + scnf->baseline;
1833 
1834  if ( cguess > 0.0 )
1835  {
1836  lncr2[ jdx ][ jlx++ ] = log( cguess );
1837  }
1838 
1839  jpx++;
1840  }
1841 
1842  v_setlpts << jlx;
1843  jdx++;
1844  }
1845  break;
1846  }
1847  case 11: // 11: "Monomer-Dimer-Trimer Equilibrium"
1848  case 12: // 12: "Monomer-Dimer-Tetramer Equilibrium"
1849  case 13: // 13: "User-Defined Monomer - N-mer - M-mer Equilibrium"
1850  {
1851  double stoich1 = 2.0;
1852  double stoich2 = (double)( modelx - 8 );
1853  runfit.stoichs[ 0 ] = stoich1;
1854  runfit.stoichs[ 1 ] = stoich2;
1855  double stoiexp = stoich1 - 1.0;
1856  double stoiex2 = stoich2 - 1.0;
1857 
1858  for ( int ii = 0; ii < scanfits.size(); ii++ )
1859  {
1860  EqScanFit* scnf = &scanfits[ ii ];
1861 
1862  if ( ! scnf->scanFit ) continue;
1863 
1864  int jstx = scnf->start_ndx;
1865  double xm_sq = sq( scnf->xvs[ jstx ] );
1866  double omega_sq = sq( M_PI * scnf->rpm / 30.0 );
1867  double tempera = scnf->tempera;
1868  double density = scnf->density;
1869  double dconst = omega_sq / ( 2.0 * R * ( K0 + tempera ) );
1870  double ODcorec1 = log( stoich1 / pow( scnf->extincts[ 0 ]
1871  * scnf->pathlen , stoiexp ) );
1872  double ODcorec2 = log( stoich2 / pow( scnf->extincts[ 0 ]
1873  * scnf->pathlen , stoiex2 ) );
1874  v_vbar[ 0 ] = US_Math2::adjust_vbar20( runfit.vbar_vals[ 0 ],
1875  tempera );
1876  double buoy = ( 1.0 - v_vbar[ 0 ] * density );
1877  jlx = 0;
1878 
1879  for ( int jj = jstx; jj < jstx + v_setpts[ jdx ]; jj++ )
1880  {
1881  double xv = sq( scnf->xvs[ jj ] ) - xm_sq;
1882  double mwv = runfit.mw_vals[ 0 ];
1883  double ampv = scnf->amp_vals[ 0 ];
1884  double ufunc0 = exp( ampv + dconst * mwv * buoy * xv );
1885  double ufunc1 = exp( stoich1 * ampv + runfit.eq_vals[ 0 ] +
1886  ODcorec1 + dconst * stoich1 * mwv * buoy * xv );
1887  double ufunc2 = exp( stoich2 * ampv + runfit.eq_vals[ 1 ] +
1888  ODcorec2 + dconst * stoich2 * mwv * buoy * xv );
1889  double cguess = ufunc0 + ufunc1 + ufunc2;
1890 
1891  y_guess[ jpx ] = cguess + scnf->baseline;
1892 
1893  if ( cguess > 0.0 )
1894  {
1895  lncr2[ jdx ][ jlx++ ] = log( cguess );
1896  }
1897 
1898  jpx++;
1899  }
1900 
1901  v_setlpts << jlx;
1902  jdx++;
1903  }
1904  break;
1905  }
1906  case 14: // 14: "2-Component Hetero-Association: A + B <=> AB"
1907  {
1908  double mwv0 = runfit.mw_vals[ 0 ];
1909  double mwv1 = runfit.mw_vals[ 1 ];
1910  double mw_ab = mwv0 + mwv1;
1911 
1912  for ( int ii = 0; ii < scanfits.size(); ii++ )
1913  {
1914  EqScanFit* scnf = &scanfits[ ii ];
1915 
1916  if ( ! scnf->scanFit ) continue;
1917 
1918  int jstx = scnf->start_ndx;
1919  double xm_sq = sq( scnf->xvs[ jstx ] );
1920  double omega_sq = sq( M_PI * scnf->rpm / 30.0 );
1921  double tempera = scnf->tempera;
1922  double density = scnf->density;
1923  double dconst = omega_sq / ( 2.0 * R * ( K0 + tempera ) );
1924  double extinc0 = scnf->extincts[ 0 ];
1925  double extinc1 = scnf->extincts[ 1 ];
1926  double ODcorrec = log( ( extinc0 + extinc1 ) /
1927  ( scnf->pathlen * extinc0 * extinc1 ) );
1928  v_vbar[ 0 ] = US_Math2::adjust_vbar20( runfit.vbar_vals[ 0 ],
1929  tempera );
1930  v_vbar[ 1 ] = US_Math2::adjust_vbar20( runfit.vbar_vals[ 1 ],
1931  tempera );
1932  v_vbar[ 2 ] = ( v_vbar[ 0 ] * mwv0 + v_vbar[ 1 ] * mwv1 )
1933  / mw_ab;
1934  double buoy0 = ( 1.0 - v_vbar[ 0 ] * density );
1935  double buoy1 = ( 1.0 - v_vbar[ 1 ] * density );
1936  double buoy2 = ( 1.0 - v_vbar[ 2 ] * density );
1937  jlx = 0;
1938 
1939  for ( int jj = jstx; jj < jstx + v_setpts[ jdx ]; jj++ )
1940  {
1941  double xv = sq( scnf->xvs[ jj ] ) - xm_sq;
1942  double ampv0 = scnf->amp_vals[ 0 ];
1943  double ampv1 = scnf->amp_vals[ 1 ];
1944  double constx = dconst * xv;
1945  double ufunc0 = exp( ampv0 + constx * mwv0 * buoy0 );
1946  double ufunc1 = exp( ampv1 + constx * mwv1 * buoy1 );
1947  double ufunc2 = exp( ampv0 + ampv1 + runfit.eq_vals[ 0 ] +
1948  ODcorrec + constx * mw_ab * buoy2 );
1949  double cguess = ufunc0 + ufunc1 + ufunc2;
1950 
1951  y_guess[ jpx ] = cguess + scnf->baseline;
1952 
1953  if ( cguess > 0.0 )
1954  {
1955  lncr2[ jdx ][ jlx++ ] = log( cguess );
1956  }
1957 
1958  jpx++;
1959  }
1960 
1961  v_setlpts << jlx;
1962  jdx++;
1963  }
1964  break;
1965  }
1966  case 15: // 15: "U-Defined self/Hetero-Assoc.: A + B <=> AB, nA <=> An"
1967  {
1968  double mwv0 = runfit.mw_vals[ 0 ];
1969  double mwv1 = runfit.mw_vals[ 1 ];
1970  double mw_ab = mwv0 + mwv1;
1971  double stoich1 = runfit.stoichs[ 0 ];
1972  double stoiexp = stoich1 - 1.0;
1973 
1974  for ( int ii = 0; ii < scanfits.size(); ii++ )
1975  {
1976  EqScanFit* scnf = &scanfits[ ii ];
1977 
1978  if ( ! scnf->scanFit ) continue;
1979 
1980  int jstx = scnf->start_ndx;
1981  double xm_sq = sq( scnf->xvs[ jstx ] );
1982  double omega_sq = sq( M_PI * scnf->rpm / 30.0 );
1983  double tempera = scnf->tempera;
1984  double density = scnf->density;
1985  double dconst = omega_sq / ( 2.0 * R * ( K0 + tempera ) );
1986  double extinc0 = scnf->extincts[ 0 ];
1987  double extinc1 = scnf->extincts[ 1 ];
1988  double ODcorec1 = log( ( extinc0 + extinc0 ) /
1989  ( scnf->pathlen * extinc0 * extinc1 ) );
1990  double ODcorec2 = log( stoich1 /
1991  pow( scnf->pathlen * extinc0, stoiexp ) );
1992  v_vbar[ 0 ] = US_Math2::adjust_vbar20( runfit.vbar_vals[ 0 ],
1993  tempera );
1994  v_vbar[ 1 ] = US_Math2::adjust_vbar20( runfit.vbar_vals[ 1 ],
1995  tempera );
1996  v_vbar[ 2 ] = ( v_vbar[ 0 ] * mwv0 + v_vbar[ 1 ] * mwv1 )
1997  / mw_ab;
1998  double buoy0 = ( 1.0 - v_vbar[ 0 ] * density );
1999  double buoy1 = ( 1.0 - v_vbar[ 1 ] * density );
2000  double buoy2 = ( 1.0 - v_vbar[ 2 ] * density );
2001  jlx = 0;
2002 
2003  for ( int jj = jstx; jj < jstx + v_setpts[ jdx ]; jj++ )
2004  {
2005  double xv = sq( scnf->xvs[ jj ] ) - xm_sq;
2006  double ampv0 = scnf->amp_vals[ 0 ];
2007  double ampv1 = scnf->amp_vals[ 1 ];
2008  double constx = dconst * xv;
2009  double ufunc0 = exp( ampv0 + constx * mwv0 * buoy0 );
2010  double ufunc1 = exp( ampv1 + constx * mwv1 * buoy1 );
2011  double ufunc2 = exp( ampv0 + ampv1 + runfit.eq_vals[ 0 ] +
2012  ODcorec1 + constx * mw_ab * buoy2 );
2013  double ufunc3 = exp( stoich1 * ampv0 + runfit.eq_vals[ 0 ] +
2014  ODcorec2 + constx * stoich1 * mw_ab * buoy0 );
2015  double cguess = ufunc0 + ufunc1 + ufunc2 + ufunc3;
2016 
2017  y_guess[ jpx ] = cguess + scnf->baseline;
2018 
2019  if ( cguess > 0.0 )
2020  lncr2[ jdx ][ jlx++ ] = log( cguess );
2021 
2022  jpx++;
2023  }
2024 
2025  v_setlpts << jlx;
2026  jdx++;
2027  }
2028  break;
2029  }
2030  case 16: // 16: "U-Defined Monomer-Nmer, some monomer is incompetent"
2031  {
2032  double mwv0 = runfit.mw_vals[ 0 ];
2033  double mwv1 = runfit.mw_vals[ 1 ];
2034  double stoich1 = runfit.stoichs[ 0 ];
2035  double stoiexp = stoich1 - 1.0;
2036 
2037  for ( int ii = 0; ii < scanfits.size(); ii++ )
2038  {
2039  EqScanFit* scnf = &scanfits[ ii ];
2040 
2041  if ( ! scnf->scanFit ) continue;
2042 
2043  int jstx = scnf->start_ndx;
2044  double xm_sq = sq( scnf->xvs[ jstx ] );
2045  double omega_sq = sq( M_PI * scnf->rpm / 30.0 );
2046  double tempera = scnf->tempera;
2047  double density = scnf->density;
2048  double dconst = omega_sq / ( 2.0 * R * ( K0 + tempera ) );
2049  double extinc0 = scnf->extincts[ 0 ];
2050  double ODcorec1 = log( stoich1 /
2051  pow( extinc0 * scnf->pathlen, stoiexp ) );
2052  v_vbar[ 0 ] = US_Math2::adjust_vbar20( runfit.vbar_vals[ 0 ],
2053  tempera );
2054  double buoy0 = ( 1.0 - v_vbar[ 0 ] * density );
2055  jlx = 0;
2056 
2057  for ( int jj = jstx; jj < jstx + v_setpts[ jdx ]; jj++ )
2058  {
2059  double xv = sq( scnf->xvs[ jj ] ) - xm_sq;
2060  double ampv0 = scnf->amp_vals[ 0 ];
2061  double ampv1 = scnf->amp_vals[ 1 ];
2062  double constx = dconst * xv;
2063  double ufunc0 = exp( ampv0 + constx * mwv0 * buoy0 );
2064  double ufunc1 = exp( ampv1 + constx * mwv1 * buoy0 );
2065  double ufunc2 = exp( stoich1 * ampv0 + runfit.eq_vals[ 0 ] +
2066  ODcorec1 + constx * stoich1 * mwv0 * buoy0 );
2067  double cguess = ufunc0 + ufunc1 + ufunc2;
2068 
2069  y_guess[ jpx ] = cguess + scnf->baseline;
2070 
2071  if ( cguess > 0.0 )
2072  lncr2[ jdx ][ jlx++ ] = log( cguess );
2073 
2074  jpx++;
2075  }
2076 
2077  v_setlpts << jlx;
2078  jdx++;
2079  }
2080  break;
2081  }
2082  case 17: // 17: "User-Defined Monomer-Nmer, some Nmer is incompetent"
2083  {
2084  double mwv0 = runfit.mw_vals[ 0 ];
2085  double mwv1 = runfit.mw_vals[ 1 ];
2086  double stoich1 = runfit.stoichs[ 0 ];
2087  double stoiexp = stoich1 - 1.0;
2088 
2089  for ( int ii = 0; ii < scanfits.size(); ii++ )
2090  {
2091  EqScanFit* scnf = &scanfits[ ii ];
2092 
2093  if ( ! scnf->scanFit ) continue;
2094 
2095  int jstx = scnf->start_ndx;
2096  double xm_sq = sq( scnf->xvs[ jstx ] );
2097  double omega_sq = sq( M_PI * scnf->rpm / 30.0 );
2098  double tempera = scnf->tempera;
2099  double density = scnf->density;
2100  double dconst = omega_sq / ( 2.0 * R * ( K0 + tempera ) );
2101  double extinc0 = scnf->extincts[ 0 ];
2102  double ODcorec1 = log( stoich1 /
2103  pow( scnf->pathlen * extinc0, stoiexp ) );
2104  v_vbar[ 0 ] = US_Math2::adjust_vbar20( runfit.vbar_vals[ 0 ],
2105  tempera );
2106  double buoy0 = ( 1.0 - v_vbar[ 0 ] * density );
2107  jlx = 0;
2108 
2109  for ( int jj = jstx; jj < jstx + v_setpts[ jdx ]; jj++ )
2110  {
2111  double xv = sq( scnf->xvs[ jj ] ) - xm_sq;
2112  double ampv0 = scnf->amp_vals[ 0 ];
2113  double ampv1 = scnf->amp_vals[ 1 ];
2114  double constx = dconst * xv;
2115  double ufunc0 = exp( ampv0 + constx * mwv0 * buoy0 );
2116  double ufunc1 = exp( ampv1 + constx * mwv1 * buoy0 );
2117  double ufunc2 = exp( stoich1 * ampv0 + runfit.eq_vals[ 0 ] +
2118  ODcorec1 + constx * stoich1 * mwv0 * buoy0 );
2119  double cguess = ufunc0 + ufunc1 + ufunc2;
2120 
2121  y_guess[ jpx ] = cguess + scnf->baseline;
2122 
2123  if ( cguess > 0.0 )
2124  lncr2[ jdx ][ jlx++ ] = log( cguess );
2125 
2126  jpx++;
2127  }
2128 
2129  v_setlpts << jlx;
2130  jdx++;
2131  }
2132  break;
2133  }
2134  case 18: // 18: "User-Defined irreversible Monomer-Nmer"
2135  {
2136  double mwv0 = runfit.mw_vals[ 0 ];
2137  double stoich1 = runfit.stoichs[ 0 ];
2138  double mwv1 = mwv0 * stoich1;
2139 
2140  for ( int ii = 0; ii < scanfits.size(); ii++ )
2141  {
2142  EqScanFit* scnf = &scanfits[ ii ];
2143 
2144  if ( ! scnf->scanFit ) continue;
2145 
2146  int jstx = scnf->start_ndx;
2147  double xm_sq = sq( scnf->xvs[ jstx ] );
2148  double omega_sq = sq( M_PI * scnf->rpm / 30.0 );
2149  double tempera = scnf->tempera;
2150  double density = scnf->density;
2151  double dconst = omega_sq / ( 2.0 * R * ( K0 + tempera ) );
2152  v_vbar[ 0 ] = US_Math2::adjust_vbar20( runfit.vbar_vals[ 0 ],
2153  tempera );
2154  double buoy0 = ( 1.0 - v_vbar[ 0 ] * density );
2155  double ampv0 = scnf->amp_vals[ 0 ];
2156  double ampv1 = scnf->amp_vals[ 1 ];
2157  jlx = 0;
2158 
2159  for ( int jj = jstx; jj < jstx + v_setpts[ jdx ]; jj++ )
2160  {
2161  double xv = sq( scnf->xvs[ jj ] ) - xm_sq;
2162  double constx = dconst * xv;
2163  double ufunc0 = exp( ampv0 + constx * mwv0 * buoy0 );
2164  double ufunc1 = exp( ampv1 + constx * mwv1 * buoy0 );
2165  double cguess = ufunc0 + ufunc1;
2166 
2167  y_guess[ jpx ] = cguess + scnf->baseline;
2168 
2169  if ( cguess > 0.0 )
2170  lncr2[ jdx ][ jlx++ ] = log( cguess );
2171 
2172  jpx++;
2173  }
2174 
2175  v_setlpts << jlx;
2176  jdx++;
2177  }
2178  break;
2179  }
2180  case 19: // 19: "User-Defined Monomer-Nmer plus contaminant"
2181  {
2182  double mwv0 = runfit.mw_vals[ 0 ];
2183  double mwv1 = runfit.mw_vals[ 1 ];
2184  double stoich1 = runfit.stoichs[ 0 ];
2185  double stoiexp = stoich1 - 1.0;
2186 
2187  for ( int ii = 0; ii < scanfits.size(); ii++ )
2188  {
2189  EqScanFit* scnf = &scanfits[ ii ];
2190 
2191  if ( ! scnf->scanFit ) continue;
2192 
2193  int jstx = scnf->start_ndx;
2194  double xm_sq = sq( scnf->xvs[ jstx ] );
2195  double omega_sq = sq( M_PI * scnf->rpm / 30.0 );
2196  double tempera = scnf->tempera;
2197  double density = scnf->density;
2198  double dconst = omega_sq / ( 2.0 * R * ( K0 + tempera ) );
2199  double extinc0 = scnf->extincts[ 0 ];
2200  double ODcorec1 = log( stoich1 /
2201  pow( scnf->pathlen * extinc0, stoiexp ) );
2202  v_vbar[ 0 ] = US_Math2::adjust_vbar20( runfit.vbar_vals[ 0 ],
2203  tempera );
2204  v_vbar[ 1 ] = US_Math2::adjust_vbar20( runfit.vbar_vals[ 1 ],
2205  tempera );
2206  double buoy0 = ( 1.0 - v_vbar[ 0 ] * density );
2207  double buoy1 = ( 1.0 - v_vbar[ 1 ] * density );
2208  double ampv0 = scnf->amp_vals[ 0 ];
2209  double ampv1 = scnf->amp_vals[ 1 ];
2210  jlx = 0;
2211 
2212  for ( int jj = jstx; jj < jstx + v_setpts[ jdx ]; jj++ )
2213  {
2214  double xv = sq( scnf->xvs[ jj ] ) - xm_sq;
2215  double constx = dconst * xv;
2216  double ufunc0 = exp( ampv0 + constx * mwv0 * buoy0 );
2217  double ufunc1 = exp( ampv1 + constx * mwv1 * buoy1 );
2218  double ufunc2 = exp( stoich1 * ampv0 + runfit.eq_vals[ 0 ] +
2219  ODcorec1 + constx * stoich1 * mwv0 * buoy0 );
2220  double cguess = ufunc0 + ufunc1 + ufunc2;
2221 
2222  y_guess[ jpx ] = cguess + scnf->baseline;
2223 
2224  if ( cguess > 0.0 )
2225  lncr2[ jdx ][ jlx++ ] = log( cguess );
2226 
2227  jpx++;
2228  }
2229 
2230  v_setlpts << jlx;
2231  jdx++;
2232  }
2233  break;
2234  }
2235  }
2236  return stat;
2237 }
2238 
2239 // Solve a general least squares (order=2)
2240 void US_EqMath::genLeastSquaresOrd2( double** MM, int points,
2241  double* y_raw, double** coeff )
2242 {
2243  double* aa[ 2 ]; // 2 x 2 matrix
2244  double a0[ 2 ]; // matrix column 1
2245  double a1[ 2 ]; // matrix column 2
2246  double bb[ 2 ]; // 2-point vector
2247  aa[ 0 ] = (double*)a0; // finish matrix creation
2248  aa[ 1 ] = (double*)a1;
2249  a0[ 0 ] = 0.0; // initialize work matrix
2250  a0[ 1 ] = 0.0;
2251  a1[ 0 ] = 0.0;
2252  a1[ 1 ] = 0.0;
2253 
2254  for ( int ii = 0; ii < 2; ii++ )
2255  { // Fill in the lower triangle of the 2 x 2 "A" matrix
2256  for ( int jj = 0; jj <= ii; jj++ )
2257  {
2258  double dotp = 0.0;
2259 
2260  for ( int kk = 0; kk < points; kk++ )
2261  dotp += ( MM[ kk ][ ii ] * MM[ kk ][ jj ] );
2262 
2263  aa[ ii ][ jj ] = dotp;
2264  }
2265  }
2266 
2267  for ( int jj = 0; jj < 2; jj++ )
2268  { // Fill in the 2 values of the "B" vector
2269  double dotp = 0.0;
2270 
2271  for ( int kk = 0; kk < points; kk++ )
2272  dotp += ( y_raw[ kk ] * MM[ kk ][ jj ] );
2273 
2274  bb[ jj ] = dotp;
2275  }
2276 
2277  // Do 2nd-order Cholesky decomposition and solve system
2278  Cholesky_DecompOrd2( (double**)aa );
2279  Cholesky_SolveSysOrd2( (double**)aa, bb );
2280 
2281  // Result is the coefficients vector
2282  (*coeff)[ 0 ] = bb[ 0 ];
2283  (*coeff)[ 1 ] = bb[ 1 ];
2284 //DbgLv(1) << "LS2: CSS b0 b1" << bb[0] << bb[1];
2285 }
2286 
2287 // Cholesky Decomposition for order=2
2288 bool US_EqMath::Cholesky_DecompOrd2( double** aa )
2289 {
2290 DbgLv(1) << "LS2:CDec a00 a10 a11" << aa[0][0] << aa[1][0] << aa[1][1];
2291  aa[ 0 ][ 0 ] = sqrt( aa[ 0 ][ 0 ] );
2292  aa[ 1 ][ 0 ] = aa[ 1 ][ 0 ] / aa[ 0 ][ 0 ];
2293 
2294  double sum = sq( aa[ 1 ][ 0 ] );
2295  double diff = aa[ 1 ][ 1 ] - sum;
2296 DbgLv(1) << "LS2:CDec sum diff" << sum << diff;
2297 
2298  if ( diff <= 0.0 ) return false;
2299 
2300  aa[ 1 ][ 1 ] = sqrt( diff );
2301  aa[ 0 ][ 1 ] = 0.0;
2302 DbgLv(1) << "LS2:CDec a00 a10 a11" << aa[0][0] << aa[1][0] << aa[1][1];
2303 
2304  return true;
2305 }
2306 
2307 // Cholesky Solve System for order=2
2308 bool US_EqMath::Cholesky_SolveSysOrd2( double** LL, double* bb )
2309 {
2310  // Forward substitution
2311 
2312  bb[ 0 ] = bb[ 0 ] / LL[ 0 ][ 0 ];
2313  bb[ 1 ] = bb[ 1 ] - LL[ 1 ][ 0 ] * bb[ 0 ];
2314  bb[ 1 ] = bb[ 1 ] / LL[ 1 ][ 1 ];
2315 
2316  // Backward substitution
2317 
2318  bb[ 1 ] = bb[ 1 ] / LL[ 1 ][ 1 ]; // AGAIN????
2319  bb[ 0 ] = bb[ 0 ] - LL[ 1 ][ 0 ] * bb[ 1 ];
2320  bb[ 0 ] = bb[ 0 ] / LL[ 0 ][ 0 ];
2321 
2322  return true;
2323 }
2324 
2325 // Convenience function: is a double value NAN (Not A valid Number)?
2326 bool US_EqMath::isNan( double value )
2327 {
2328  if ( value != value ) // NAN, the one case where value!=value
2329  return true;
2330 
2331  double avalue = qAbs( value ); // Also mark NAN if beyond float min,max
2332 
2333  if ( avalue < dflt_min || avalue > dflt_max )
2334  return true;
2335 
2336  return false;
2337 }
2338 
2339 // Calculate runs counts
2340 void US_EqMath::calc_runs()
2341 {
2342  // Accumulate runs and +/- counts for scans
2343  runfit.nbr_runs = 0;
2344  int jptx = 0;
2345 
2346  for ( int ii = 0; ii < scanfits.size(); ii++ )
2347  {
2348  EqScanFit* scnf = &scanfits[ ii ];
2349 
2350  if ( ! scnf->scanFit ) continue;
2351 
2352  scnf->nbr_posr = 0;
2353  scnf->nbr_negr = 0;
2354  scnf->runs = 0;
2355  bool lastpos = ( y_guess[ jptx ] > scnf->xvs[ scnf->start_ndx ] );
2356 
2357  for ( int jj = scnf->start_ndx; jj < scnf->stop_ndx + 1; jj++ )
2358  {
2359  bool thispos = ( y_guess[ jptx++ ] > scnf->yvs[ jj ] );
2360 
2361  if ( thispos )
2362  {
2363  scnf->nbr_posr++;
2364  if ( ! lastpos ) scnf->runs++;
2365  }
2366 
2367  else
2368  {
2369  scnf->nbr_negr++;
2370  if ( lastpos ) scnf->runs++;
2371  }
2372 
2373  lastpos = thispos;
2374  }
2375 
2376  runfit.nbr_runs += scnf->runs;
2377  }
2378 }
2379 
2380 // Calculate the integral from meniscus to bottom for each exponential term
2381 void US_EqMath::calc_integral()
2382 {
2383  int ncomps = runfit.nbr_comps;
2384 
2385  switch( modelx )
2386  {
2387  case 0: // 0: "1-Component, Ideal"
2388  case 1: // 1: "2-Component, Ideal, Noninteracting"
2389  case 2: // 2: "3-Component, Ideal, Noninteracting"
2390  for ( int ii = 0; ii < scanfits.size(); ii++ )
2391  {
2392  EqScanFit* scnf = &scanfits[ ii ];
2393 
2394  if ( ! scnf->scanFit ) continue;
2395 
2396  double bottom = calc_bottom( scnf->rpm );
2397  scnf->bottom = bottom;
2398  double deltr = ( bottom - scnf->meniscus ) / 5000.0;
2399  double xm_sq = sq( scnf->xvs[ scnf->start_ndx ] );
2400  double omega_sq = sq( M_PI * scnf->rpm / 30.0 );
2401  double tempera = scnf->tempera;
2402  double density = scnf->density;
2403  double dconst = omega_sq / ( 2.0 * R * ( K0 + tempera ) );
2404 
2405  for ( int kk = 0; kk < ncomps; kk++ )
2406  {
2407  double vbar = US_Math2::adjust_vbar20( runfit.vbar_vals[ kk ],
2408  tempera );
2409  double buoy = ( 1.0 - vbar * density );
2410  double xval = scnf->meniscus;
2411  double amplv = scnf->amp_vals[ kk ];
2412  double ampl2 = exp( amplv );
2413  double mwfac = runfit.mw_vals[ kk ] * dconst * buoy;
2414  double sumi = 0.0;
2415 
2416  for ( int mm = 0; mm < 4999; mm++ )
2417  {
2418  xval += deltr;
2419  double xvsq = sq( xval ) - xm_sq;
2420  double ampl1 = ampl2;
2421  ampl2 = exp( amplv + mwfac * xvsq );
2422  sumi += ( deltr * ( ampl1 + ampl2 ) / 2.0 );
2423  }
2424 
2425  scnf->integral[ kk ] = sumi;
2426  }
2427  }
2428  break;
2429  case 3: // 3: "Fixed Molecular Weight Distribution"
2430  for ( int ii = 0; ii < scanfits.size(); ii++ )
2431  {
2432  EqScanFit* scnf = &scanfits[ ii ];
2433 
2434  if ( ! scnf->scanFit ) continue;
2435 
2436  double bottom = calc_bottom( scnf->rpm );
2437  scnf->bottom = bottom;
2438  double deltr = ( bottom - scnf->meniscus ) / 5000.0;
2439  double xm_sq = sq( scnf->xvs[ scnf->start_ndx ] );
2440  double omega_sq = sq( M_PI * scnf->rpm / 30.0 );
2441  double tempera = scnf->tempera;
2442  double density = scnf->density;
2443  double dconst = omega_sq / ( 2.0 * R * ( K0 + tempera ) );
2444 
2445  for ( int kk = 0; kk < ncomps; kk++ )
2446  {
2447  double vbar = US_Math2::adjust_vbar20( runfit.vbar_vals[ kk ],
2448  tempera );
2449  double buoy = ( 1.0 - vbar * density );
2450  double xval = scnf->meniscus;
2451  double amplv = scnf->amp_vals[ kk ];
2452  double ampl2 = exp( amplv );
2453  double mwfac = runfit.mw_vals[ kk ] * dconst * buoy;
2454  double sumi = 0.0;
2455 
2456  for ( int mm = 0; mm < 4999; mm++ )
2457  {
2458  xval += deltr;
2459  double xvsq = sq( xval ) - xm_sq;
2460  double ampl1 = ampl2;
2461  ampl2 = exp( amplv + mwfac * xvsq );
2462  sumi += ( deltr * ( ampl1 + ampl2 ) / 2.0 );
2463  }
2464 
2465  scnf->integral[ kk ] = sumi;
2466  }
2467  }
2468  break;
2469  case 4: // 4: "Monomer-Dimer Equilibrium"
2470  case 5: // 5: "Monomer-Trimer Equilibrium"
2471  case 6: // 6: "Monomer-Tetramer Equilibrium"
2472  case 7: // 7: "Monomer-Pentamer Equilibrium"
2473  case 8: // 8: "Monomer-Hexamer Equilibrium"
2474  case 9: // 9: "Monomer-Heptamer Equilibrium"
2475  case 10: // 10: "User-Defined Monomer-Nmer Equilibrium"
2476  {
2477  double stoich1 = (double)( modelx - 2 );
2478  double stoiexp = stoich1 - 1.0;
2479  runfit.stoichs[ 0 ] = stoich1;
2480  double eqval0 = runfit.eq_vals[ 0 ];
2481 
2482  for ( int ii = 0; ii < scanfits.size(); ii++ )
2483  {
2484  EqScanFit* scnf = &scanfits[ ii ];
2485 
2486  if ( ! scnf->scanFit ) continue;
2487 
2488  double bottom = calc_bottom( scnf->rpm );
2489  scnf->bottom = bottom;
2490  double deltr = ( bottom - scnf->meniscus ) / 5000.0;
2491  double xm_sq = sq( scnf->xvs[ scnf->start_ndx ] );
2492  double omega_sq = sq( M_PI * scnf->rpm / 30.0 );
2493  double tempera = scnf->tempera;
2494  double density = scnf->density;
2495  double ODcorrec = log( stoich1 / pow( scnf->extincts[ 0 ]
2496  * scnf->pathlen , stoiexp ) ) + eqval0;
2497  double dconst = omega_sq / ( 2.0 * R * ( K0 + tempera ) );
2498 
2499  double vbar = US_Math2::adjust_vbar20( runfit.vbar_vals[ 0 ],
2500  tempera );
2501  double buoy = ( 1.0 - vbar * density );
2502  double xval = scnf->meniscus;
2503  double amplv = scnf->amp_vals[ 0 ];
2504  double ampl2 = exp( amplv );
2505  double mwfac = runfit.mw_vals[ 0 ] * dconst * buoy;
2506  double sumi = 0.0;
2507 
2508  for ( int mm = 0; mm < 4999; mm++ )
2509  {
2510  xval += deltr;
2511  double xvsq = sq( xval ) - xm_sq;
2512  double ampl1 = ampl2;
2513  ampl2 = exp( amplv + mwfac * xvsq );
2514  sumi += ( deltr * ( ampl1 + ampl2 ) / 2.0 );
2515  }
2516 
2517  scnf->integral[ 0 ] = sumi;
2518 
2519  xval = scnf->meniscus;
2520  sumi = 0.0;
2521  amplv = scnf->amp_vals[ 0 ] * stoich1;
2522  ampl2 = exp( amplv + ODcorrec );
2523  mwfac *= stoich1;
2524 
2525  for ( int mm = 0; mm < 4999; mm++ )
2526  {
2527  xval += deltr;
2528  double xvsq = sq( xval ) - xm_sq;
2529  double ampl1 = ampl2;
2530  ampl2 = exp( amplv + ODcorrec + mwfac * xvsq );
2531  sumi += ( deltr * ( ampl1 + ampl2 ) / 2.0 );
2532  }
2533 
2534  scnf->integral[ 1 ] = sumi;
2535  }
2536  break;
2537  }
2538  case 11: // 11: "Monomer-Dimer-Trimer Equilibrium"
2539  case 12: // 12: "Monomer-Dimer-Tetramer Equilibrium"
2540  case 13: // 13: "User-Defined Monomer - N-mer - M-mer Equilibrium"
2541  {
2542  double stoich1 = 2.0;
2543  double stoich2 = (double)( modelx - 8 );
2544  runfit.stoichs[ 0 ] = stoich1;
2545  runfit.stoichs[ 1 ] = stoich2;
2546  //double stoiexp = stoich1 - 1.0;
2547  //double stoiex2 = stoich2 - 1.0;
2548  double stoiexp = stoich1;
2549  double stoiex2 = stoich2;
2550  double eqval0 = runfit.eq_vals[ 0 ];
2551  double eqval1 = runfit.eq_vals[ 1 ];
2552 
2553  for ( int ii = 0; ii < scanfits.size(); ii++ )
2554  {
2555  EqScanFit* scnf = &scanfits[ ii ];
2556 
2557  if ( ! scnf->scanFit ) continue;
2558 
2559  double bottom = calc_bottom( scnf->rpm );
2560  scnf->bottom = bottom;
2561  double deltr = ( bottom - scnf->meniscus ) / 5000.0;
2562  double xm_sq = sq( scnf->xvs[ scnf->start_ndx ] );
2563  double omega_sq = sq( M_PI * scnf->rpm / 30.0 );
2564  double tempera = scnf->tempera;
2565  double density = scnf->density;
2566  double ODcorre1 = log( stoich1 / pow( scnf->extincts[ 0 ]
2567  * scnf->pathlen , stoiexp ) ) + eqval0;
2568  double ODcorre2 = log( stoich2 / pow( scnf->extincts[ 0 ]
2569  * scnf->pathlen , stoiex2 ) ) + eqval1;
2570  double dconst = omega_sq / ( 2.0 * R * ( K0 + tempera ) );
2571 
2572  double vbar = US_Math2::adjust_vbar20( runfit.vbar_vals[ 0 ],
2573  tempera );
2574  double buoy = ( 1.0 - vbar * density );
2575  double xval = scnf->meniscus;
2576  double amplv = scnf->amp_vals[ 0 ];
2577  double ampl2 = exp( amplv );
2578  double mwfac = runfit.mw_vals[ 0 ] * dconst * buoy;
2579  double sumi = 0.0;
2580 
2581  for ( int mm = 0; mm < 4999; mm++ )
2582  {
2583  xval += deltr;
2584  double xvsq = sq( xval ) - xm_sq;
2585  double ampl1 = ampl2;
2586  ampl2 = exp( amplv + mwfac * xvsq );
2587  sumi += ( deltr * ( ampl1 + ampl2 ) / 2.0 );
2588  }
2589 
2590  scnf->integral[ 0 ] = sumi;
2591 
2592  xval = scnf->meniscus;
2593  amplv = scnf->amp_vals[ 0 ] * stoich1;
2594  ampl2 = exp( amplv + ODcorre1 );
2595  mwfac = runfit.mw_vals[ 0 ] * dconst * buoy * stoich1;
2596  sumi = 0.0;
2597 
2598  for ( int mm = 0; mm < 4999; mm++ )
2599  {
2600  xval += deltr;
2601  double xvsq = sq( xval ) - xm_sq;
2602  double ampl1 = ampl2;
2603  ampl2 = exp( amplv + ODcorre1 + mwfac * xvsq );
2604  sumi += ( deltr * ( ampl1 + ampl2 ) / 2.0 );
2605  }
2606 
2607  scnf->integral[ 1 ] = sumi;
2608 
2609  xval = scnf->meniscus;
2610  amplv = scnf->amp_vals[ 0 ] * stoich2;
2611  ampl2 = exp( amplv + ODcorre2 );
2612  mwfac = runfit.mw_vals[ 0 ] * dconst * buoy * stoich2;
2613  sumi = 0.0;
2614 
2615  for ( int mm = 0; mm < 4999; mm++ )
2616  {
2617  xval += deltr;
2618  double xvsq = sq( xval ) - xm_sq;
2619  double ampl1 = ampl2;
2620  ampl2 = exp( amplv + ODcorre1 + mwfac * xvsq );
2621  sumi += ( deltr * ( ampl1 + ampl2 ) / 2.0 );
2622  }
2623 
2624  scnf->integral[ 2 ] = sumi;
2625 
2626  }
2627  break;
2628  }
2629  case 14: // 14: "2-Component Hetero-Association: A + B <=> AB"
2630  {
2631  double mwv0 = runfit.mw_vals[ 0 ];
2632  double mwv1 = runfit.mw_vals[ 1 ];
2633  double mwv2 = runfit.mw_vals[ 2 ];
2634  double mw_ab = mwv0 + mwv2;
2635  double eqval0 = runfit.eq_vals[ 0 ];
2636 
2637  for ( int ii = 0; ii < scanfits.size(); ii++ )
2638  {
2639  EqScanFit* scnf = &scanfits[ ii ];
2640 
2641  if ( ! scnf->scanFit ) continue;
2642 
2643  double bottom = calc_bottom( scnf->rpm );
2644  scnf->bottom = bottom;
2645  double deltr = ( bottom - scnf->meniscus ) / 5000.0;
2646  double xm_sq = sq( scnf->xvs[ scnf->start_ndx ] );
2647  double omega_sq = sq( M_PI * scnf->rpm / 30.0 );
2648  double tempera = scnf->tempera;
2649  double density = scnf->density;
2650  double dconst = omega_sq / ( 2.0 * R * ( K0 + tempera ) );
2651 
2652  double vbar0 = US_Math2::adjust_vbar20( runfit.vbar_vals[ 0 ],
2653  tempera );
2654  double vbar1 = US_Math2::adjust_vbar20( runfit.vbar_vals[ 1 ],
2655  tempera );
2656  double vbar2 = ( vbar0 * mwv0 + vbar1 * mwv1 ) / mw_ab;
2657  double buoy0 = ( 1.0 - vbar0 * density );
2658  double buoy1 = ( 1.0 - vbar1 * density );
2659  double buoy2 = ( 1.0 - vbar2 * density );
2660 
2661  double xval = scnf->meniscus;
2662  double amplv = scnf->amp_vals[ 0 ];
2663  double ampl2 = exp( amplv );
2664  double mwfac = mwv0 * dconst * buoy0;
2665  double sumi = 0.0;
2666 
2667  for ( int mm = 0; mm < 4999; mm++ )
2668  {
2669  xval += deltr;
2670  double xvsq = sq( xval ) - xm_sq;
2671  double ampl1 = ampl2;
2672  ampl2 = exp( amplv + mwfac * xvsq );
2673  sumi += ( deltr * ( ampl1 + ampl2 ) / 2.0 );
2674  }
2675 
2676  scnf->integral[ 0 ] = sumi;
2677 
2678  xval = scnf->meniscus;
2679  amplv = scnf->amp_vals[ 1 ];
2680  ampl2 = exp( amplv );
2681  mwfac = mwv1 * dconst * buoy1;
2682  sumi = 0.0;
2683 
2684  for ( int mm = 0; mm < 4999; mm++ )
2685  {
2686  xval += deltr;
2687  double xvsq = sq( xval ) - xm_sq;
2688  double ampl1 = ampl2;
2689  ampl2 = exp( amplv + mwfac * xvsq );
2690  sumi += ( deltr * ( ampl1 + ampl2 ) / 2.0 );
2691  }
2692 
2693  scnf->integral[ 1 ] = sumi;
2694 
2695  double extinc0 = scnf->extincts[ 0 ];
2696  double extinc1 = scnf->extincts[ 1 ];
2697  double extin_ab = extinc0 + extinc1;
2698  double ODcorrec = log( extin_ab
2699  / ( extinc0 * extinc1 * scnf->pathlen ) )
2700  + eqval0;
2701  xval = scnf->meniscus;
2702  amplv = scnf->amp_vals[ 0 ] + scnf->amp_vals[ 1 ];
2703  ampl2 = exp( amplv + ODcorrec );
2704  mwfac = mw_ab * dconst * buoy2;
2705  sumi = 0.0;
2706 
2707  for ( int mm = 0; mm < 4999; mm++ )
2708  {
2709  xval += deltr;
2710  double xvsq = sq( xval ) - xm_sq;
2711  double ampl1 = ampl2;
2712  ampl2 = exp( amplv + ODcorrec + mwfac * xvsq );
2713  sumi += ( deltr * ( ampl1 + ampl2 ) / 2.0 );
2714  }
2715 
2716  scnf->integral[ 2 ] = sumi;
2717 
2718  }
2719  break;
2720  }
2721  case 15: // 15: "U-Defined self/Hetero-Assoc.: A + B <=> AB, nA <=> An"
2722  case 16: // 16: "U-Defined Monomer-Nmer, some monomer is incompetent"
2723  case 17: // 17: "User-Defined Monomer-Nmer, some Nmer is incompetent"
2724  case 18: // 18: "User-Defined irreversible Monomer-Nmer"
2725  case 19: // 19: "User-Defined Monomer-Nmer plus contaminant"
2726  break;
2727  }
2728 }
2729 
2730 // Calculate the bottom for a scan, using rotor information and speed
2731 double US_EqMath::calc_bottom( double rpm )
2732 {
2733  return ( runfit.bottom_pos + runfit.rcoeffs[ 0 ] * rpm
2734  + runfit.rcoeffs[ 1 ] * sq( rpm ) );
2735 }
Generated on Sun Jun 21 2015 11:40:41 for UltraScan III by   doxygen 1.8.3.1-20130324