UltraScan III
us_math2.cpp
Go to the documentation of this file.
1 
3 #include <stdlib.h>
4 #include <math.h>
5 
6 #include "us_math2.h"
7 #include "us_constants.h"
8 #include "us_dataIO.h"
9 #include "us_matrix.h"
10 #include "us_settings.h"
11 
12 /* The function implements the Box-Muller algorithm for generating
13  * pairs of independent standard normally distributed (zero expectation,
14  * unit variance) random numbers, given a source of uniformly distributed
15  * random numbers.
16 
17  * The function returns the input mean value modified by the standard
18  * deviation weighted by a random normally distributed increment.
19 
20 */
21 
22 double US_Math2::box_muller( double m, double s )
23 {
24  static bool use_last = false;
25 
26  double x1;
27  double x2;
28  double w;
29  double y1;
30  static double y2;
31 
32  if ( use_last ) // Use value from previous call
33  {
34  y1 = y2;
35  use_last = false;
36  }
37  else
38  {
39  do
40  {
41  x1 = 2.0 * ranf() - 1.0;
42  x2 = 2.0 * ranf() - 1.0;
43  w = sq( x1 ) + sq( x2 );
44  } while ( w >= 1.0 );
45 
46  w = sqrt( ( -2.0 * log( w ) ) / w );
47  y1 = x1 * w;
48  y2 = x2 * w;
49  use_last = true;
50  }
51 
52  return m + y1 * s;
53 }
54 
55 // This function returns a randomly distributed double value R
56 // such that 0.0 <= R < 1.0
57 double US_Math2::ranf( void )
58 {
59  return (double)qrand() / ( (double)RAND_MAX + 1.0 );
60 }
61 
62 
63 double US_Math2::linefit( double** x , double** y , double* slope,
64  double* intercept, double* sigma, double* correlation,
65  int arraysize )
66 {
67  double sumx = 0.0;
68  double sumy = 0.0;
69  double sumxy = 0.0;
70  double sumx_sq = 0.0;
71  double sumy_sq = 0.0;
72  double sumchi_sq = 0.0;
73  double average;
74 
75  for ( int i = 0; i < arraysize; i++ ) sumy += (*y)[ i ];
76 
77  average = sumy / arraysize;
78 
79  for ( int i = 0; i < arraysize; i++ )
80  {
81  sumx += (*x)[ i ];
82  sumxy += (*x)[ i ] * (*y)[ i ];
83  sumx_sq += (*x)[ i ] * (*x)[ i ];
84  sumy_sq += (*y)[ i ] * (*y)[ i ];
85 
86  sumchi_sq += ( (*y)[ i ] - average ) * ( (*y)[ i ] - average );
87  }
88 
89  *slope = ( arraysize * sumxy - sumx * sumy ) /
90  ( arraysize * sumx_sq - sumx * sumx );
91 
92  *intercept = ( sumy - *slope * sumx ) / arraysize;
93 
94  *correlation = ( arraysize * sumxy - sumx * sumy ) /
95  ( sqrt( arraysize * sumx_sq - sumx * sumx ) *
96  sqrt( arraysize * sumy_sq - sumy * sumy )
97  );
98 
99  *sigma = sqrt( sumchi_sq ) / arraysize;
100 
101  return average;
102 }
103 /*
104 float US_Math2::linefit( float** x , float** y , float* slope,
105  float* intercept, float* sigma, float* correlation,
106  int arraysize )
107 {
108  float sumx = 0.0;
109  float sumy = 0.0;
110  float sumxy = 0.0;
111  float sumx_sq = 0.0;
112  float sumy_sq = 0.0;
113  float sumchi_sq = 0.0;
114  float average;
115 
116  for ( int i = 0; i < arraysize; i++ ) sumy += (*y)[ i ];
117 
118  average = sumy / arraysize;
119 
120  for ( int i = 0; i < arraysize; i++ )
121  {
122  sumx += (*x)[ i ];
123  sumxy += (*x)[ i ] * (*y)[ i ];
124  sumx_sq += (*x)[ i ] * (*x)[ i ];
125  sumy_sq += (*y)[ i ] * (*y)[ i ];
126 
127  sumchi_sq += ( (*y)[ i ] - average ) * ( (*y)[ i ] - average );
128  }
129 
130  *slope = ( arraysize * sumxy - sumx * sumy ) /
131  ( arraysize * sumx_sq - sumx * sumx );
132 
133  *intercept = ( sumy - *slope * sumx ) / arraysize;
134 
135  *correlation = ( arraysize * sumxy - sumx * sumy ) /
136  ( sqrt( arraysize * sumx_sq - sumx * sumx ) *
137  sqrt( arraysize * sumy_sq - sumy * sumy ) );
138 
139  *sigma = sqrt( sumchi_sq ) / arraysize;
140 
141  return average;
142 }
143 */
144 
145 // This function determines the point on a curve that is nearest to
146 // a given point. The returned point may be an exact curve point or one
147 // interpolated between the two nearest curve points. An optional value
148 // may be returned which is the exact or interpolated value in a third
149 // dimension or associated array.
150 
151 int US_Math2::nearest_curve_point( double* xs, double* ys, const int npoints,
152  bool interp, double& xgiven, double& ygiven,
153  double* xnear, double* ynear, double* zs, double* znear )
154 {
155  int jmin = npoints / 2;
156  int jend = npoints - 1;
157  double xnorm = qAbs( xs[ 0 ] - xs[ jend ] ); // X normalizing factor
158  double ynorm = qAbs( ys[ 0 ] - ys[ jend ] ); // Y normalizing factor
159  xnorm = ( xnorm == 0.0 ) ? 1.0 : ( 1.0 / xnorm );
160  ynorm = ( ynorm == 0.0 ) ? 1.0 : ( 1.0 / ynorm );
161  double xdif = ( xs[ 0 ] - xgiven ) * xnorm;
162  double ydif = ( ys[ 0 ] - ygiven ) * ynorm;
163  double dmin = sqrt( sq( xdif ) + sq( ydif ) );
164 
165  // Find the nearest point to the given point
166  for ( int jj = 0; jj < npoints; jj++ )
167  {
168  xdif = ( xs[ jj ] - xgiven ) * xnorm;
169  ydif = ( ys[ jj ] - ygiven ) * ynorm;
170  double dval = sqrt( sq( xdif ) + sq( ydif ) );
171 
172  if ( dval < dmin )
173  {
174  dmin = dval;
175  jmin = jj;
176  }
177  }
178 
179  if ( interp )
180  { // Interpolate a curve point between nearest point and an adjacent one
181  int jmin1 = jmin - 1; // Index to point preceding the nearest
182  int jmin2 = jmin + 1; // Index to point following the nearest
183  double dmin1 = dmin;
184  double dmin2 = dmin;
185 
186  if ( jmin == 0 )
187  { // If nearest is first on curve, use the following point
188  xdif = ( xs[ jmin2 ] - xgiven ) * xnorm;
189  ydif = ( ys[ jmin2 ] - ygiven ) * ynorm;
190  dmin2 = sqrt( sq( xdif ) + sq( ydif ) );
191  }
192 
193  else if ( jmin2 == npoints )
194  { // If nearest is last on curve, use the preceding point
195  jmin2 = jmin1;
196  xdif = ( xs[ jmin2 ] - xgiven ) * xnorm;
197  ydif = ( ys[ jmin2 ] - ygiven ) * ynorm;
198  dmin2 = sqrt( sq( xdif ) + sq( ydif ) );
199  }
200 
201  else
202  { // If the nearest point is somewhere in the middle of the curve,
203  // determine whether to pick the adjacent point from the one
204  // following or the one preceding. Do this by finding which of the
205  // two gives the closest path to the curve.
206  xdif = ( xs[ jmin1 ] - xgiven ) * xnorm;
207  ydif = ( ys[ jmin1 ] - ygiven ) * ynorm;
208  dmin1 = sqrt( sq( xdif ) + sq( ydif ) ); // Dist of preceding
209  xdif = ( xs[ jmin2 ] - xgiven ) * xnorm;
210  ydif = ( ys[ jmin2 ] - ygiven ) * ynorm;
211  dmin2 = sqrt( sq( xdif ) + sq( ydif ) ); // Dist of following
212  double ratio = dmin1 / ( dmin + dmin1 );
213  double ratio2 = 1.0 - ratio;
214  double x1 = xs[ jmin ] * ratio + xs[ jmin1 ] * ratio2;
215  double y1 = ys[ jmin ] * ratio + ys[ jmin1 ] * ratio2;
216  ratio = dmin2 / ( dmin + dmin2 );
217  ratio2 = 1.0 - ratio;
218  double x2 = xs[ jmin ] * ratio + xs[ jmin2 ] * ratio2;
219  double y2 = ys[ jmin ] * ratio + ys[ jmin2 ] * ratio2;
220  xdif = ( x1 - xgiven ) * xnorm;
221  ydif = ( y1 - ygiven ) * ynorm;
222  double dmin1g = sqrt( sq( xdif ) + sq( ydif ) );
223  xdif = ( x2 - xgiven ) * xnorm;
224  ydif = ( y2 - ygiven ) * ynorm;
225  double dmin2g = sqrt( sq( xdif ) + sq( ydif ) );
226 
227  if ( dmin2g > dmin1g )
228  { // Where the test perpendicular towards the preceding point is
229  // closer to the given point, choose that preceding point to
230  // use as adjacent. Otherwise, default to the point following.
231  dmin2 = dmin1;
232  jmin2 = jmin1;
233  }
234  }
235 
236  // Determine ratio values based on relative distances of the nearest
237  // curve point and a curve point adjacent to it.
238  double ratio = dmin2 / ( dmin + dmin2 );
239  double ratio2 = 1.0 - ratio;
240 
241  if ( xnear != NULL && ynear != NULL )
242  { // Interpolate the curve point nearest to the given point
243  *xnear = xs[ jmin ] * ratio + xs[ jmin2 ] * ratio2;
244  *ynear = ys[ jmin ] * ratio + ys[ jmin2 ] * ratio2;
245  }
246 
247  if ( zs != NULL && znear != NULL )
248  { // Interpolate the additional value associated with the nearest point
249  *znear = zs[ jmin ] * ratio + zs[ jmin2 ] * ratio2;
250  }
251  }
252 
253  else
254  { // Return the exact point found to be the nearest
255  if ( xnear != NULL && ynear != NULL )
256  { // Return the curve point nearest to the given point
257  *xnear = xs[ jmin ];
258  *ynear = ys[ jmin ];
259  }
260 
261  if ( zs != NULL && znear != NULL )
262  { // If desired, return the additional value associated with the point
263  *znear = zs[ jmin ];
264  }
265  }
266 
267  return jmin;
268 }
269 
270 // This function determines the intersection point of two lines.
271 // This version takes as input the slopes and intercepts of the two lines.
272 bool US_Math2::intersect( double& slope1, double& intcp1,
273  double& slope2, double& intcp2, double* xisec, double* yisec )
274 {
275  bool isect = ( slope1 != slope2 ); // Parallel lines do not intersect!
276 
277  if ( isect )
278  { // Return intersection point of non-parallel lines
279  *xisec = ( intcp2 - intcp1 ) / ( slope1 - slope2 );
280  *yisec = *xisec * slope1 + intcp1;
281  }
282 
283  return isect;
284 }
285 
286 // This function determines the intersection point of two fitted lines.
287 // This version takes as input the x and y arrays of two curves
288 // and returns the intersection of straight lines fitted to them.
289 bool US_Math2::intersect( double* x1s, double* y1s, int npoint1,
290  double* x2s, double* y2s, int npoint2, double* xisec, double* yisec )
291 {
292  double slope1;
293  double intcp1;
294  double slope2;
295  double intcp2;
296  double sigma;
297  double corre;
298 
299  // Compute the slope and intercept of a line fitted to the first curve
300  US_Math2::linefit( &x1s, &y1s, &slope1, &intcp1, &sigma, &corre, npoint1 );
301 
302  // Compute the slope and intercept of a line fitted to the second curve
303  US_Math2::linefit( &x2s, &y2s, &slope2, &intcp2, &sigma, &corre, npoint2 );
304 
305  // Use slopes and intercepts to compute and return the intersection point
306  bool isect = intersect( slope1, intcp1, slope2, intcp2, xisec, yisec );
307 
308  return isect;
309 }
310 
311 void US_Math2::calc_vbar( Peptide& pep, const QString& sequence,
312  double temperature )
313 {
314  pep.vbar_sum = 0.0;
315  pep.mw = 0.0;
316  pep.weight = 0.0;
317  pep.e280 = 0.0;
318  pep.vbar20 = 0.0;
319  pep.vbar = 0.0;
320 
321  pep.a = sequence.count( "a", Qt::CaseInsensitive );
322  pep.b = sequence.count( "b", Qt::CaseInsensitive );
323  pep.c = sequence.count( "c", Qt::CaseInsensitive );
324  pep.d = sequence.count( "d", Qt::CaseInsensitive );
325  pep.e = sequence.count( "e", Qt::CaseInsensitive );
326  pep.f = sequence.count( "f", Qt::CaseInsensitive );
327  pep.g = sequence.count( "g", Qt::CaseInsensitive );
328  pep.h = sequence.count( "h", Qt::CaseInsensitive );
329  pep.i = sequence.count( "i", Qt::CaseInsensitive );
330  pep.j = sequence.count( "j", Qt::CaseInsensitive );
331  pep.k = sequence.count( "k", Qt::CaseInsensitive );
332  pep.l = sequence.count( "l", Qt::CaseInsensitive );
333  pep.m = sequence.count( "m", Qt::CaseInsensitive );
334  pep.n = sequence.count( "n", Qt::CaseInsensitive );
335  pep.o = sequence.count( "o", Qt::CaseInsensitive );
336  pep.p = sequence.count( "p", Qt::CaseInsensitive );
337  pep.q = sequence.count( "q", Qt::CaseInsensitive );
338  pep.r = sequence.count( "r", Qt::CaseInsensitive );
339  pep.s = sequence.count( "s", Qt::CaseInsensitive );
340  pep.t = sequence.count( "t", Qt::CaseInsensitive );
341 
342  pep.v = sequence.count( "v", Qt::CaseInsensitive );
343  pep.w = sequence.count( "w", Qt::CaseInsensitive );
344  pep.x = sequence.count( "x", Qt::CaseInsensitive );
345  pep.x += sequence.count( "?" );
346  pep.y = sequence.count( "y", Qt::CaseInsensitive );
347  pep.z = sequence.count( "z", Qt::CaseInsensitive );
348  pep.dab = sequence.count( "+" );
349  pep.dpr = sequence.count( "@" );
350 
351  pep.residues = pep.a + pep.b + pep.c + pep.d + pep.e + pep.f +
352  pep.g + pep.h + pep.i + pep.j + pep.k + pep.l +
353  pep.m + pep.n + pep.o + pep.p + pep.q + pep.r +
354  pep.s + pep.t + pep.v + pep.w + pep.x +
355  pep.y + pep.z + pep.dab + pep.dpr;
356 
357  pep.mw += pep.a * ALANINE_MW;
358  pep.weight += pep.a * ALANINE_MW * ALANINE_VBAR;
359  pep.vbar_sum += pep.a * ALANINE_VBAR;
360 
361  pep.mw += pep.b * M_OR_D_MW;
362  pep.weight += pep.b * M_OR_D_MW * M_OR_D_VBAR;
363  pep.vbar_sum += pep.b * M_OR_D_VBAR;
364 
365  pep.mw += pep.c * CYSTEINE_MW;
366  pep.weight += pep.c * CYSTEINE_MW * CYSTEINE_VBAR;
367  pep.vbar_sum += pep.c * CYSTEINE_VBAR;
368  pep.e280 += pep.c * CYSTEINE_E280;
369 
370  pep.mw += pep.d * ASPARTIC_ACID_MW;
371  pep.weight += pep.d * ASPARTIC_ACID_MW * ASPARTIC_ACID_VBAR;
372  pep.vbar_sum += pep.d * ASPARTIC_ACID_VBAR;
373 
374  pep.mw += pep.e * GLUTAMIC_ACID_MW;
375  pep.weight += pep.e * GLUTAMIC_ACID_MW * GLUTAMIC_ACID_VBAr;
376  pep.vbar_sum += pep.e * GLUTAMIC_ACID_VBAr;
377 
378  pep.mw += pep.f * PHENYLALANINE_MW;
379  pep.weight += pep.f * PHENYLALANINE_MW * PHENYLALANINE_VBAR;
380  pep.vbar_sum += pep.f * PHENYLALANINE_VBAR;
381 
382  pep.mw += pep.g * GLYCINE_MW;
383  pep.weight += pep.g * GLYCINE_MW * GLYCINE_VBAR;
384  pep.vbar_sum += pep.g * GLYCINE_VBAR;
385 
386  pep.mw += pep.h * HISTIDINE_MW;
387  pep.weight += pep.h * HISTIDINE_MW * HISTIDINE_VBAR;
388  pep.vbar_sum += pep.h * HISTIDINE_VBAR;
389 
390  pep.mw += pep.i * ISOLEUCINE_MW;
391  pep.weight += pep.i * ISOLEUCINE_MW * ISOLEUCINE_VBAR;
392  pep.vbar_sum += pep.i * ISOLEUCINE_VBAR;
393 
394  pep.mw += pep.j * HAO_MW; // James Nowick: Hao
395  pep.weight += pep.j * HAO_MW * HAO_VBAR;
396  pep.vbar_sum += pep.j * HAO_VBAR;
397  pep.e280 += pep.j * HAO_E280;
398 
399  // John Kulp: diaminobutyric acid
400  pep.mw += pep.dab * DIAMINOBUTYRIC_ACID_MW;
401  pep.weight += pep.dab * DIAMINOBUTYRIC_ACID_MW * DIAMINOBUTYRIC_ACID_VBAR;
402  pep.vbar_sum += pep.dab * DIAMINOBUTYRIC_ACID_VBAR;
403 
404  // John Kulp: diaminopropanoic acid
405  pep.mw += pep.dpr * DIAMINOPROPANOIC_ACID_MW;
406  pep.weight += pep.dpr * DIAMINOPROPANOIC_ACID_MW * DIAMINOPROPANOIC_ACID_VBAR;
407  pep.vbar_sum += pep.dpr * DIAMINOPROPANOIC_ACID_VBAR;
408 
409  pep.mw += pep.k * LYSINE_MW;
410  pep.weight += pep.k * LYSINE_MW * LYSINE_VBAR;
411  pep.vbar_sum += pep.k * LYSINE_VBAR;
412 
413  pep.mw += pep.l * LEUCINE_MW;
414  pep.weight += pep.l * LEUCINE_MW * LEUCINE_VBAR;
415  pep.vbar_sum += pep.l * LEUCINE_VBAR;
416 
417  pep.mw += pep.m * METHIONINE_MW;
418  pep.weight += pep.m * METHIONINE_MW * METHIONINE_VBAR;
419  pep.vbar_sum += pep.m * METHIONINE_VBAR;
420 
421  pep.mw += pep.n * ASPARAGINE_MW;
422  pep.weight += pep.n * ASPARAGINE_MW * ASPARAGINE_VBAR;
423  pep.vbar_sum += pep.n * ASPARAGINE_VBAR;
424 
425  // James Nowick: delta-linked Ornithine
426  pep.mw += pep.o * ORNITHIN_MW;
427  pep.weight += pep.o * ORNITHIN_MW * ORNITHIN_VBAR;
428  pep.vbar_sum += pep.o * ORNITHIN_VBAR;
429  pep.e280 += pep.o * ORNITHIN_E280;
430 
431  pep.mw += pep.p * PROLINE_MW;
432  pep.weight += pep.p * PROLINE_MW * PROLINE_VBAR;
433  pep.vbar_sum += pep.p * PROLINE_VBAR;
434 
435  pep.mw += pep.q * GLUTAMINE_MW;
436  pep.weight += pep.q * GLUTAMINE_MW * GLUTAMINE_VBAR;
437  pep.vbar_sum += pep.q * GLUTAMINE_VBAR;
438 
439  pep.mw += pep.r * ARGININE_MW;
440  pep.weight += pep.r * ARGININE_MW * ARGININE_VBAR;
441  pep.vbar_sum += pep.r * ARGININE_VBAR;
442 
443  pep.mw += pep.s * SERINE_MW;
444  pep.weight += pep.s * SERINE_MW * SERINE_VBAR;
445  pep.vbar_sum += pep.s * SERINE_VBAR;
446 
447  pep.mw += pep.t * THREONINE_MW;
448  pep.weight += pep.t * THREONINE_MW * THREONINE_VBAR;
449  pep.vbar_sum += pep.t * THREONINE_VBAR;
450 
451  pep.mw += pep.v * VALINE_MW;
452  pep.weight += pep.v * VALINE_MW * VALINE_VBAR;
453  pep.vbar_sum += pep.v * VALINE_VBAR;
454 
455  pep.mw += pep.w * TRYPTOPHAN_MW;
456  pep.weight += pep.w * TRYPTOPHAN_MW * TRYPTOPHAN_VBAR;
457  pep.vbar_sum += pep.w * TRYPTOPHAN_VBAR;
458  pep.e280 += pep.w * TRYPTOPHAN_E280;
459 
460  pep.mw += pep.x * UNK_MW; // Using an average
461  pep.weight += pep.x * UNK_MW * UNK_VBAR;
462  pep.vbar_sum += pep.x * UNK_VBAR;
463 
464  pep.mw += pep.y * TYROSINE_MW;
465  pep.weight += pep.y * TYROSINE_MW * TYROSINE_VBAR;
466  pep.vbar_sum += pep.y * TYROSINE_VBAR;
467  pep.e280 += pep.y * TYROSINE_E280;
468 
469  pep.mw += pep.z * E_OR_Q_MW;
470  pep.weight += pep.z * E_OR_Q_MW * E_OR_Q_VBAR;
471  pep.vbar_sum += pep.z * E_OR_Q_VBAR;
472 
473  // These equations are empirically derived. The values above are for
474  // a temperature of 25C. Adjust the values to 20C from the current
475  // temperature.
476 
477  if ( pep.mw > 0.0 )
478  {
479  pep.vbar20 = ( pep.weight / pep.mw ) - 0.002125;
480  pep.vbar = ( ( pep.weight / pep.mw ) + 4.25e-4 * ( temperature - 25 ) );
481 
482  // Add one water molecule for the end of the chain
483  pep.mw += WATER_MW;
484  }
485 }
486 
487 void US_Math2::data_correction( double t, SolutionData& d )
488 {
489  double xi_max = 1.000028e-3 ;
490  double c0 = 999.83952 ;
491  double c1 = 16.945176 ;
492  double c2 = -7.9870401e-3 ;
493  double c3 = -46.170461e-6 ;
494  double c4 = 105.56302e-9 ;
495  double c5 = -280.54253e-12 ;
496  double b = 16.879850e-3 ;
497 
498  double t2 = t * t;
499  double t3 = t * t2;
500  double t4 = t * t3;
501  double t5 = t * t4;
502 //if ( qAbs(d.vbar-0.72) > 0.001 )
503 //qDebug() << "M2:dacor: t" << t << "dens visc" << d.density << d.viscosity
504 // << "vbar20 vbartb" << d.vbar20 << d.vbar;
505 
526  d.density_wt =
527  xi_max * ( c0 + c1 * t + c2 * t2 + c3 * t3 + c4 * t4 + c5 * t5 ) /
528  ( 1.0 + b * t );
529 
563  double exponent;
564  double t20;
565 
566  if ( t < 20.0 )
567  {
568  c0 = 1301.0;
569  c1 = 998.333;
570  c2 = 8.1855;
571  c3 = 0.00585;
572  c4 = 3.30233;
573 
574  t20 = t - 20.0;
575 
576  exponent = ( c0 / ( c1 + c2 * t20 + c3 * sq( t20 ) ) ) - c4;
577 
578  d.viscosity_wt = 100.0 * pow( 10.0, exponent );
579  }
580  else
581  {
582  c1 = 1.3272;
583  c2 = 1.053e-3;
584  c3 = 105.0;
585 
586  t20 = 20.0 - t;
587 
588  exponent = ( c1 * t20 - c2 * sq( t20 ) ) / ( t + c3 );
589 
590  d.viscosity_wt = VISC_20W * pow( 10.0, exponent );
591  }
592 
593  if ( d.manual )
594  {
595  d.density_tb = d.density;
596  d.viscosity_tb = d.viscosity;
597  }
598  else
599  {
600  d.density_tb = d.density * d.density_wt / DENS_20W;
601  d.viscosity_tb = d.viscosity * d.viscosity_wt / VISC_20W;
602  }
603  d.buoyancyb = 1.0 - d.vbar * d.density_tb;
604  d.buoyancyw = 1.0 - d.vbar20 * DENS_20W;
605  d.s20w_correction = ( d.buoyancyw / d.buoyancyb ) * ( d.viscosity_tb / VISC_20W );
606 
607  double K = t + K0;
608 
609  d.D20w_correction = ( K20 / K ) * ( d.viscosity / VISC_20W );
610 //if ( qAbs(d.vbar-0.72) > 0.001 ) {
611 //qDebug() << "M2:dacor: denstb denswt" << d.density_tb << d.density_wt << "manual" << d.manual;
612 //qDebug() << "M2:dacor: visctb viscwt" << d.viscosity_tb << d.viscosity_wt;
613 //qDebug() << "M2:dacor: buoyb buoyw" << d.buoyancyb << d.buoyancyw;
614 //qDebug() << "M2:dacor: vbtb vb20" << d.vbar << d.vbar20;
615 //qDebug() << "M2:dacor: dentb den20 den" << d.density_tb << DENS_20W << d.density;
616 //qDebug() << "M2:dacor: vistb vis20 vis" << d.viscosity_tb << VISC_20W << d.viscosity;
617 //qDebug() << "M2:dacor: scorr dcorr" << d.s20w_correction << d.D20w_correction; }
618 
619 }
620 
621 double US_Math2::normal_distribution( double sigma, double mean, double x )
622 {
623  double exponent = -sq( ( x - mean ) / sigma ) / 2.0;
624  return exp( exponent ) / sqrt( 2.0 * M_PI * sq( sigma ) );
625 }
626 
627 double US_Math2::time_correction(
628  const QVector< US_DataIO::EditedData >& dataList )
629 {
630  int size = dataList[ 0 ].scanData.size();
631 
632  for ( int ii = 1; ii < dataList.size(); ii++ )
633  size += dataList[ ii ].scanData.size();
634 
635  int count = 0;
636 
637  QVector< double > vecx( size );
638  QVector< double > vecy( size );
639  double* x = vecx.data();
640  double* y = vecy.data();
641 
642  double c[ 2 ]; // Looking for a linear fit
643 
644  for ( int i = 0; i < dataList.size(); i++ )
645  {
646  const US_DataIO::EditedData* e = &dataList[ i ];
647 
648  for ( int j = 0; j < e->scanData.size(); j++ )
649  {
650  x[ count ] = e->scanData[ j ].omega2t;
651  y[ count ] = e->scanData[ j ].seconds;
652  count++;
653  }
654  }
655 
656  US_Matrix::lsfit( c, x, y, count, 2 );
657 
658  return c[ 0 ]; // Return the time value corresponding to zero omega2t
659 }
660 
661 uint US_Math2::randomize( void )
662 {
663  QTime now = QTime::currentTime();
664 
665  uint seed = now.msec()
666  + now.second() * 1000
667  + now.minute() * 60000
668  + now.hour() * 3600000;
669 
670 #ifdef UNIX
671  seed -= getpid();
672 #endif
673 
674  qsrand( seed );
675  return seed;
676 }
677 
678 uint US_Math2::randomize( uint seed )
679 {
680  if ( seed == 0 )
681  seed = randomize();
682  else
683  qsrand( seed );
684 
685  return seed;
686 }
687 
688 double US_Math2::erfc( double x )
689 {
690  // error function erfc(x) with fractional error everywhere less
691  // than 1.2 × 10 7.
692 
693  double z = fabs( x );
694  double t = 1.0 / ( 1.0 + 0.5 * z );
695 
696  double ans = t * exp( -z *
697  z -1.26551223 + t *
698  ( 1.00002368 + t *
699  ( 0.37409196 + t *
700  ( 0.09678418 + t *
701  ( -0.18628806 + t *
702  ( 0.27886807 + t *
703  ( -1.13520398 + t *
704  ( 1.48851587 + t *
705  ( -0.82215223 + t * 0.17087277
706  )
707  )
708  )
709  )
710  )
711  )
712  )
713  )
714  );
715 
716  return x >= 0.0 ? ans : 2.0 - ans;
717 }
718 
719 
720 int US_Math2::nnls( double* a, int a_dim1, int m, int n,
721  double* b,
722  double* x,
723  double* rnorm,
724  double* wp,
725  double* zzp,
726  int* indexp
727  )
728 {
729  int pfeas, ret = 0, iz, jz;
730  double d1, d2, sm, up, ss;
731  int k, j=0, l, izmax=0, ii, jj=0, ip;
732  double temp, wmax, t, alpha, asave, dummy, unorm, ztest, cc;
733 
734  /* Check the parameters and data */
735  if ( m <= 0 || n <= 0 || a == NULL || b == NULL || x == NULL ) return 2;
736 
737  /* Allocate memory for working space, if required */
738  QVector< double > wVec;
739  QVector< double > zVec;
740  QVector< int > iVec;
741 
742  double* w;
743  double* zz;
744  int* index;
745 
746  if ( wp != NULL )
747  w = wp;
748  else
749  {
750  wVec.fill( 0.0, n );
751  w = wVec.data();
752  }
753 
754  if ( zzp != NULL )
755  zz = zzp;
756  else
757  {
758  zVec.fill( 0.0, m );
759  zz = zVec.data();
760  }
761 
762  if ( indexp != NULL )
763  index = indexp;
764  else
765  {
766  iVec.fill( 0, n );
767  index = iVec.data();
768  }
769 
770  if ( w == NULL || zz == NULL || index == NULL ) return 2;
771 
772  /* Initialize the arrays index[] and x[] */
773  for( int k = 0; k < n; k++ )
774  {
775  x [ k ] = 0.0;
776  index[ k ] = k;
777  }
778 
779  int iz2 = n - 1;
780  int iz1 = 0;
781  int nsetp = 0;
782  int npp1 = 0;
783 
784  /* Main loop; quit if all coeffs are already in the solution or */
785  /* if M cols of A have been triangularized */
786 
787  int iter = 0;
788  int itmax = n * 3;
789  int ja_dim1;
790 
791  while ( iz1 <= iz2 && nsetp < m )
792  {
793  /* Compute components of the dual (negative gradient) vector W[] */
794  for ( iz = iz1; iz <= iz2; iz++ )
795  {
796  j = index[ iz ];
797  ja_dim1 = j * a_dim1;
798  sm = 0.0;
799 
800  for( l = npp1; l < m; l++ )
801  sm += a[ l + ja_dim1 ] * b[ l ];
802 
803  w[ j ] = sm;
804  }
805 
806  while ( true )
807  {
808  /* Find largest positive W[j] */
809  for ( wmax = 0.0, iz = iz1; iz <= iz2; iz++ )
810  {
811  j = index[ iz ];
812  if ( w[ j ] > wmax )
813  {
814  wmax = w[ j ];
815  izmax = iz;
816  }
817  }
818 
819  /* Terminate if wmax <= 0.0; */
820  /* it indicates satisfaction of the Kuhn-Tucker conditions */
821  if ( wmax <= 0.0 ) break;
822 
823  iz = izmax;
824  j = index[ iz ];
825 
826  /* The sign of W[ j ] is ok for j to be moved to set P. */
827  /* Begin the transformation and check new diagonal element to avoid */
828  /* near linear dependence. */
829 
830  ja_dim1 = j * a_dim1;
831  asave = a[ npp1 + ja_dim1 ];
832  _nnls_h12( 1, npp1, npp1 + 1, m, &a[ ja_dim1 ],
833  1, &up, &dummy, 1, 1, 0 );
834 
835  unorm = 0.0;
836 
837  if ( nsetp != 0 )
838  for ( l = 0; l < nsetp; l++ )
839  {
840  d1 = a[ l + ja_dim1 ];
841  unorm += d1 * d1;
842  }
843 
844  unorm = sqrt( unorm );
845 
846  d2 = unorm + ( d1 = a[ npp1 + ja_dim1 ], fabs( d1 ) ) * 0.01;
847 
848  if ( ( d2 - unorm ) > 0.0 )
849  {
850  /* Col j is sufficiently independent. Copy B into ZZ, update ZZ
851  and solve for ztest ( = proposed new value for X[ j ] ) */
852 
853  for ( l = 0; l < m; l++ ) zz[ l ] =b[ l ];
854 
855  _nnls_h12( 2, npp1, npp1 + 1, m, &a[ ja_dim1 ],
856  1, &up, zz, 1, 1, 1 );
857 
858  ztest = zz[ npp1 ] / a[ npp1 + ja_dim1 ];
859 
860  /* See if ztest is positive */
861  if ( ztest > 0.0 ) break;
862  }
863 
864  /* Reject j as a candidate to be moved from set Z to set P. Restore */
865  /* A[npp1,j], set W[j]=0., and loop back to test dual coeffs again */
866  a[ npp1 + ja_dim1 ] = asave;
867  w[ j ] = 0.0;
868 
869  } // while( true )
870 
871  if ( wmax <= 0.0 ) break;
872 
873  /* Index j = index[ iz ] has been selected to be moved from set Z to set P.
874  Update B and indices, apply householder transformations to cols in
875  new set Z, zero subdiagonal elts in col j, set W[j]=0. */
876 
877  for ( l = 0; l < m; ++l ) b[ l ] =zz[ l ];
878 
879  index[ iz ] = index[ iz1 ];
880  index[ iz1 ] = j;
881  iz1++;
882  nsetp = npp1 + 1;
883  npp1++;
884 
885  if ( iz1 <= iz2 )
886  for ( jz = iz1; jz <= iz2; jz++ )
887  {
888  jj = index[ jz ];
889  _nnls_h12( 2, nsetp - 1, npp1, m, &a[ ja_dim1 ],
890  1, &up, &a[ jj * a_dim1 ], 1, a_dim1, 1 );
891  }
892 
893  if ( nsetp !=m )
894  for ( l = npp1; l < m; l++ )
895  a[ l + ja_dim1 ] = 0.0;
896 
897  w[ j ]= 0.0;
898 
899  /* Solve the triangular system; store the solution temporarily in Z[] */
900  for ( l = 0; l < nsetp; l++ )
901  {
902  ip = nsetp - ( l + 1 );
903 
904  if ( l != 0 )
905  for ( ii = 0; ii <= ip; ii++ )
906  zz[ ii ] -= a[ ii + jj * a_dim1 ] * zz[ ip + 1 ];
907 
908  jj = index[ ip ];
909  zz[ ip ] /= a[ ip + jj *a_dim1 ];
910  }
911 
912  /* Secondary loop begins here */
913  while ( ++iter < itmax )
914  {
915  /* See if all new constrained coeffs are feasible;
916  if not, compute alpha */
917 
918  for ( alpha = 2.0, ip = 0; ip < nsetp; ip++ )
919  {
920  l = index[ ip ];
921  if ( zz[ ip ] <= 0.0 )
922  {
923  t = -x[ l ] / ( zz[ ip ] - x[ l ]);
924  if ( alpha > t )
925  {
926  alpha = t;
927  jj = ip - 1;
928  }
929  }
930  }
931 
932  /* If all new constrained coeffs are feasible then still alpha==2. */
933  /* If so, then exit from the secondary loop to main loop */
934 
935  if ( alpha == 2.0 ) break;
936 
937  /* Use alpha ( 0.0 < alpha < 1.0 ) to interpolate between old
938  * X and new ZZ */
939 
940  for ( ip = 0; ip < nsetp; ip++ )
941  {
942  l = index[ ip ];
943  x[ l ] += alpha * ( zz[ ip ] - x[ l ] );
944  }
945 
946  /* Modify A and B and the INDEX arrays to move coefficient i */
947  /* from set P to set Z. */
948 
949  k = index[ jj + 1 ];
950  pfeas = 1;
951 
952  do
953  {
954  x[ k ] = 0.0;
955  if ( jj != ( nsetp - 1 ) )
956  {
957  jj++;
958  for ( j = jj + 1; j < nsetp; j++ )
959  {
960  ii = index[ j ];
961  int iia_dim1 = ii * a_dim1;
962  index[ j - 1 ] = ii;
963 
964  _nnls_g1( a [ j - 1 + iia_dim1 ], a[ j + iia_dim1 ],
965  &cc, &ss, &a[ j - 1 + iia_dim1 ] );
966 
967  for ( a[ j + iia_dim1 ] = 0.0, l = 0; l < n; l++ )
968  if ( l != ii )
969  {
970  int jla_dim1 = j + l * a_dim1;
971  /* Apply procedure G2 (CC,SS,A(J-1,L),A(J,L)) */
972  temp = a[ jla_dim1 - 1 ];
973  a[ jla_dim1 - 1 ] = cc * temp + ss * a[ jla_dim1 ];
974  a[ jla_dim1 ] = -ss * temp + cc * a[ jla_dim1 ];
975  }
976 
977  /* Apply procedure G2 (CC,SS,B(J-1),B(J)) */
978  temp = b[ j - 1 ];
979  b[ j - 1 ] = cc * temp + ss * b[ j ];
980  b[ j ] = -ss * temp + cc * b[ j ];
981  }
982  }
983 
984  npp1 = nsetp - 1;
985  nsetp--;
986  iz1--;
987  index[ iz1 ] = k;
988 
989  /* See if the remaining coeffs in set P are feasible; they should */
990  /* be because of the way alpha was determined. If any are */
991  /* infeasible it is due to round-off error. Any that are */
992  /* nonpositive will be set to zero and moved from set P to set Z */
993  for( jj = 0; jj < nsetp; jj++ )
994  {
995  k = index[ jj ];
996  if ( x[ k ] <= 0.0 )
997  {
998  pfeas = 0;
999  break;
1000  }
1001  }
1002  } while ( pfeas == 0 );
1003 
1004  /* Copy b[] into zz[], then solve again and loop back */
1005  for ( k = 0; k < m; k++ )
1006  zz[ k ] = b[ k ];
1007 
1008  for ( l = 0; l < nsetp; l++ )
1009  {
1010  ip = nsetp - ( l + 1 );
1011  if ( l != 0 )
1012  for ( ii = 0; ii <= ip; ii++ )
1013  zz[ ii ] -= a[ ii + jj * a_dim1 ] * zz[ ip + 1 ];
1014 
1015  jj = index[ ip ];
1016  zz[ ip ] /= a[ ip + jj * a_dim1 ];
1017  }
1018  } /* end of secondary loop */
1019 
1020  if ( iter > itmax )
1021  {
1022  ret = 1;
1023  break;
1024  }
1025 
1026  for ( ip = 0; ip < nsetp; ip++ )
1027  {
1028  k = index[ ip ];
1029  x[ k ] = zz [ ip ];
1030  }
1031  } /* end of main loop */
1032 
1033  /* Compute the norm of the final residual vector */
1034 
1035  sm = 0.0;
1036 
1037  if ( npp1 < m )
1038  for ( k = npp1; k < m; k++ )
1039  sm += ( b[ k ] * b[ k ] );
1040  else
1041  for( j = 0; j < n; j++ )
1042  w[ j ] = 0.0;
1043 
1044  if ( rnorm != NULL ) *rnorm = sqrt( sm );
1045 
1046  return ret;
1047 }
1048 
1049 /*****************************************************************************
1050  *
1051  * Compute orthogonal rotation matrix:
1052  * (C, S) so that (C, S)(A) = ( sqrt( A**2 + B**2 ) )
1053  * (-S,C) (-S,C)(B) ( 0 )
1054  * Compute sig = sqrt(A**2+B**2):
1055  * sig is computed last to allow for the possibility that sig may be in
1056  * the same location as A or B.
1057  */
1058 void US_Math2::_nnls_g1( double a, double b, double* cterm,
1059  double* sterm, double* sig )
1060 {
1061  double d1;
1062  double xr;
1063  double yr;
1064 
1065  if ( fabs( a ) > fabs( b ) )
1066  {
1067  xr = b / a;
1068  d1 = xr;
1069  yr = sqrt( d1 * d1 + 1.0 );
1070  d1 = 1.0 / yr;
1071  *cterm = ( a >= 0.0 ? fabs( d1 ) : -fabs( d1 ) );
1072  *sterm = ( *cterm ) * xr;
1073  *sig = fabs( a ) * yr;
1074  }
1075  else if ( b != 0.0 )
1076  {
1077  xr = a / b;
1078  d1 = xr;
1079  yr = sqrt( d1*d1 + 1.0 );
1080  d1 = 1.0 / yr;
1081  *sterm = ( b >= 0.0 ? fabs( d1 ) : -fabs( d1 ) );
1082  *cterm = ( *sterm ) * xr;
1083  *sig = fabs( b ) * yr;
1084  }
1085  else
1086  {
1087  *sig = 0.0;
1088  *cterm = 0.0;
1089  *sterm = 1.0;
1090  }
1091 } /* _nnls_g1 */
1092 
1093 
1094 /*****************************************************************************
1095  *
1096  * Construction and/or application of a single Householder transformation:
1097  * Q = I + U*(U**T)/B
1098  *
1099  * Function returns 0 if succesful, or >0 in case of erroneous parameters.
1100  *
1101  */
1102 int US_Math2::_nnls_h12 (
1103  int mode,
1104  /* mode=1 to construct and apply a Householder transformation, or
1105  mode=2 to apply a previously constructed transformation */
1106  int lpivot, /* Index of the pivot element */
1107  int l1,
1108  int m,
1109  /* Transformation is constructed to zero elements indexed from l1 to M */
1110  double* u,
1111  int u_dim1,
1112  double* up,
1113  /* With mode=1: On entry, u[] must contain the pivot vector.
1114  On exit, u[] and up contain quantities defining the vector u[] of
1115  the Householder transformation. */
1116  /* With mode=2: On entry, u[] and up should contain quantities previously
1117  computed with mode=1. These will not be modified. */
1118  /* u_dim1 is the storage increment between elements. */
1119  double* cm,
1120  /* On entry, cm[] must contain the matrix (set of vectors) to which the
1121  Householder transformation is to be applied. On exit, cm[] will contain
1122  the set of transformed vectors */
1123  int ice, /* Storage increment between elements of vectors in cm[] */
1124  int icv, /* Storage increment between vectors in cm[] */
1125  int ncv /* Nr of vectors in cm[] to be transformed;
1126  if ncv<=0, then no operations will be done on cm[] */
1127  )
1128 {
1129  double d1, d2, b, clinv, cl, sm;
1130  int incr, k, j, i2, i3, i4;
1131 
1132  /* Check parameters */
1133  if ( mode != 1 && mode != 2 ) return 1;
1134 
1135  if ( m < 1 || u == NULL || u_dim1 < 1 || cm == NULL ) return 2;
1136 
1137  if ( lpivot < 0 || lpivot >= l1 || l1 >= m ) return 0;
1138 
1139  /* Function Body */
1140  cl = ( d1 = u[ lpivot * u_dim1 ], fabs( d1 ) );
1141 
1142  if ( mode == 2 ) // Apply transformation I+U*(U**T)/B to cm[]
1143  {
1144  if ( cl <= 0.0 ) return 0;
1145  }
1146  else // Construct the transformation
1147  {
1148  for ( j = l1; j < m; j++ )
1149  { /* Computing MAX */
1150  d2 = ( d1 = u[ j * u_dim1 ], fabs( d1 ) );
1151  if ( d2 > cl ) cl=d2;
1152  }
1153 
1154  if ( cl <= 0.0 ) return 0;
1155 
1156  clinv = 1.0 / cl;
1157 
1158  /* Computing 2nd power */
1159  d1 = u[ lpivot * u_dim1 ] * clinv;
1160  sm = d1 * d1;
1161 
1162  for( j = l1; j < m; j++ )
1163  {
1164  d1 = u[ j * u_dim1 ] * clinv;
1165  sm += d1 * d1;
1166  }
1167 
1168  cl *= sqrt( sm );
1169  if ( u[ lpivot * u_dim1 ] > 0.0 ) cl = -cl;
1170 
1171  *up = u[ lpivot * u_dim1 ] - cl;
1172  u[ lpivot * u_dim1 ] = cl;
1173  }
1174 
1175  if ( ncv <= 0 ) return 0;
1176 
1177  b = (*up) * u[ lpivot * u_dim1 ];
1178 
1179  /* b must be nonpositive here; if b>=0., then return */
1180  if ( b >= 0.0 ) return 0;
1181 
1182  b = 1.0 / b;
1183  i2 = 1 -icv + ice * lpivot;
1184  incr = ice * ( l1 - lpivot );
1185 
1186  for ( j = 0; j < ncv; j++ )
1187  {
1188  i2 += icv;
1189  i3 = i2 + incr;
1190  i4 = i3;
1191  sm = cm[ i2 - 1 ] * (*up);
1192 
1193  for ( k = l1; k < m; k++ )
1194  {
1195  sm += cm[ i3 - 1 ] * u[ k * u_dim1 ];
1196  i3 += ice;
1197  }
1198 
1199  if ( sm != 0.0 )
1200  {
1201  sm *= b;
1202  cm[ i2 - 1 ] += sm * (*up);
1203 
1204  for ( k = l1; k < m; k++ )
1205  {
1206  cm[ i4 - 1 ] += sm * u[ k * u_dim1 ];
1207  i4 += ice;
1208  }
1209  }
1210  }
1211 
1212  return 0;
1213 } /* _nnls_h12 */
1214 
1215 void US_Math2::gaussian_smoothing( QVector< double >& array, int smooth )
1216 {
1217  if ( smooth <= 1 ) return;
1218 
1219  // Apply a normalized Gaussian smoothing kernel that goes out
1220  // to 2 standard deviations
1221 
1222  int points = array.size();
1223 
1224  QVector< double > temp_array( points );
1225  QVector< double > y_weights( smooth );
1226  double x_weight;
1227 
1228  temp_array = array;
1229 
1230  // Standard deviation = 1.0, Mean = 0;
1231  x_weight = 0.0;
1232  y_weights[ 0 ] = normal_distribution( 1.0, 0.0, x_weight );
1233  double increment = 2.0 / (double)smooth;
1234 
1235  // Only calculate half a Gaussian, since the other side is symmetric
1236  for ( int i = 1; i < smooth; i++ )
1237  {
1238  x_weight += increment;
1239  y_weights[ i ] = normal_distribution( 1.0, 0.0, x_weight );
1240  }
1241 
1242  // First, take care of the left border, using an "appearing frame" algorithm,
1243  // starting with half a frame visible:
1244 
1245  for ( int j = 0; j < smooth; j++ ) // Loop over all border point centers
1246  {
1247  double sum = 0.0;
1248  double sum_y = 0.0;
1249  int position = 0;
1250 
1251  // Sum all applicable points on the left of center
1252  for ( int k = j - 1; k >= 0; k-- )
1253  {
1254  position++;
1255  sum += y_weights[ position ] * temp_array[ k ];
1256  sum_y += y_weights[ position ];
1257  }
1258 
1259  position = 0;
1260 
1261  // Sum the weighted points right of center, including center
1262  for ( int k = j; k < j + smooth; k++ )
1263  {
1264  sum += y_weights[ position ] * temp_array[ k ];
1265  sum_y += y_weights[ position ];
1266  position++;
1267  }
1268 
1269  // Normalize by the sum of all weights that were used
1270  array[ j ] = sum / sum_y;
1271  }
1272 
1273  // Now deal with all non-border points:
1274  for ( int j = smooth; j < points - smooth; j++ )
1275  {
1276  double sum = 0.0;
1277  double sum_y = 0.0;
1278  int position = 0;
1279 
1280  // Sum all applicable points on the left of center
1281  for ( int k = j - 1; k >= j - smooth + 1; k-- )
1282  {
1283  position++;
1284  sum += y_weights[ position ] * temp_array[ k ];
1285  sum_y += y_weights[ position ];
1286  }
1287 
1288  position = 0;
1289 
1290  // Sum the weighted points right of center, including center
1291  for ( int k = j; k < j + smooth; k++ )
1292  {
1293  sum += y_weights[ position ] * temp_array[ k ];
1294  sum_y += y_weights[ position ];
1295  position++;
1296  }
1297 
1298  // Normalize by the sum of all weights that were used
1299  array[ j ] = sum / sum_y;
1300  }
1301 
1302  // Now deal with all points from the right border, using a "disappearing
1303  // frame" algorithm, starting with a full frame minus 1 point visible:
1304 
1305  // Loop over all right-border points
1306  for ( int j = points - smooth; j < points; j++ )
1307  {
1308  double sum = 0.0;
1309  double sum_y = 0.0;
1310  int position = 0;
1311 
1312  // Sum all points on the left of center
1313  for ( int k = j - 1; k >= j - smooth + 1; k-- )
1314  {
1315  position++;
1316  sum += y_weights[ position ] * temp_array[ k ];
1317  sum_y += y_weights[ position ];
1318  }
1319 
1320  position = 0;
1321 
1322  // Right of center, including center
1323  for ( int k = j; k < points; k++ )
1324  {
1325  sum += y_weights[ position ] * temp_array[ k ];
1326  sum_y += y_weights[ position ];
1327  position++;
1328  }
1329 
1330  // normalize by the sum of all weights that were used
1331  array[ j ] = sum / sum_y;
1332  }
1333 }
1334 
1335 // Calculate the weighted average of temperature-corrected vbars
1336 double US_Math2::calcCommonVbar( US_Solution& solution, double temperature )
1337 {
1338  double cvbar = TYPICAL_VBAR;
1339  double vbsum = 0.0;
1340  double wtsum = 0.0;
1341 
1342  for ( int ii = 0; ii < solution.analyteInfo.size(); ii++ )
1343  {
1344  double vb20 = solution.analyteInfo[ ii ].analyte.vbar20;
1345  // Use adjusted vbar if PROTEIN
1346  double vbar = solution.analyteInfo[ ii ].analyte.type == US_Analyte::PROTEIN ?
1347  US_Math2::adjust_vbar20( vb20, temperature ) :
1348  vb20;
1349  double wt = solution.analyteInfo[ ii ].analyte.mw * solution.analyteInfo[ ii ].amount;
1350  vbsum += ( vbar * wt );
1351  wtsum += wt;
1352  }
1353 
1354  if ( wtsum != 0.0 )
1355  cvbar = vbsum / wtsum; // Common vbar is the weighted average
1356 
1357  return cvbar;
1358 }
1359 
1360 // Compute the best number of uniform grid repetitions for 2DSA,
1361 // given specified initial total grid points in the 's' and 'k'
1362 // (f/f0 or vbar) dimensions. Modify grid points slightly to insure
1363 // they are multiples of grid repetitions and are with reasonable limits.
1364 int US_Math2::best_grid_reps( int& ngrid_s, int& ngrid_k )
1365 {
1366  const int min_grid = 10; // Grid points at least 10
1367  const int min_reps = 1; // Repetitions at least 1
1368 #if 0
1369  const int max_grid = 200; // Grid points at most 200
1370  const int min_subg = 40; // Sub-grid size at least 40
1371  const int max_subg = 200; // Sub-grid size at most 200
1372  const int max_reps = 40; // Repetitions at most 40
1373 #endif
1374 #if 1
1375  const int max_grid = 800; // Grid points at most 800
1376  const int min_subg = 10; // Sub-grid size at least 10
1377  const int max_subg = 800; // Sub-grid size at most 800
1378  const int max_reps = 160; // Repetitions at most 160
1379 #endif
1380 
1381  // Insure grid points are within reasonable limits
1382  ngrid_s = qMin( max_grid, qMax( min_grid, ngrid_s ) );
1383  ngrid_k = qMin( max_grid, qMax( min_grid, ngrid_k ) );
1384 
1385  // Compute the initial best grid-repetitions value
1386  int nreps_g = qRound( pow( (double)( ngrid_s * ngrid_k ), 0.25 ) );
1387 #if 1
1388 //*DEBUG*
1389  double grfact = 1.0;
1390  QStringList tgrfact=US_Settings::debug_text();
1391 
1392  for ( int ii = 0; ii < tgrfact.count(); ii++ )
1393  {
1394  if ( tgrfact[ ii ].startsWith( "grfact=" ) )
1395  grfact = QString( tgrfact[ ii ]) .section( "=", 1, 1 ).toDouble();
1396  }
1397  nreps_g = qRound( grfact * nreps_g );
1398 //*DEBUG*
1399 #endif
1400 
1401  // Compute the initial sub-grid point counts in each dimension,
1402  // insuring the next highest integers are used.
1403  int nsubg_s = ( ngrid_s + nreps_g - 1 ) / nreps_g;
1404  int nsubg_k = ( ngrid_k + nreps_g - 1 ) / nreps_g;
1405 
1406  // Adjust values until the product yields no more than 200 sub-grid points
1407  while ( ( nsubg_s * nsubg_k ) > max_subg && nreps_g < max_reps )
1408  { // Increase grid-reps and recompute sub-grid points
1409  nreps_g++;
1410  nsubg_s = ( ngrid_s + nreps_g - 1 ) / nreps_g;
1411  nsubg_k = ( ngrid_k + nreps_g - 1 ) / nreps_g;
1412  }
1413 
1414  // Adjust values until the product yields no less than 40 sub-grid points
1415  while ( ( nsubg_s * nsubg_k ) < min_subg && nreps_g > min_reps )
1416  { // Decrease grid-reps and recompute sub-grid points
1417  nreps_g--;
1418  nsubg_s = ( ngrid_s + nreps_g - 1 ) / nreps_g;
1419  nsubg_k = ( ngrid_k + nreps_g - 1 ) / nreps_g;
1420  }
1421 
1422  // Recalculate the total grid points in each dimension
1423  // to be multiples of the grid repetitions
1424  ngrid_s = nsubg_s * nreps_g;
1425  ngrid_k = nsubg_k * nreps_g;
1426 
1427  // Return the computed number of grid repetitions
1428  return nreps_g;
1429 }
1430 
Generated on Sun Jun 21 2015 11:40:42 for UltraScan III by   doxygen 1.8.3.1-20130324