UltraScan III
us_vhw_enhanced.cpp
Go to the documentation of this file.
1 
3 #include <QApplication>
4 #include <QtSvg>
5 
6 #include "us_vhw_enhanced.h"
7 #include "us_license_t.h"
8 #include "us_license.h"
9 #include "us_settings.h"
10 #include "us_gui_settings.h"
11 #include "us_matrix.h"
12 #include "us_passwd.h"
13 #include "us_constants.h"
14 #include "us_astfem_rsa.h"
15 
16 #define SEDC_NOVAL -9999.0
17 
18 // main program
19 int main( int argc, char* argv[] )
20 {
21  QApplication application( argc, argv );
22 
23  #include "main1.inc"
24 
25  // License is OK. Start up.
26 
27  US_vHW_Enhanced w;
28  w.show();
29  return application.exec();
30 }
31 
32 // US_vHW_Enhanced class constructor
33 US_vHW_Enhanced::US_vHW_Enhanced() : US_AnalysisBase2()
34 {
35  // set up the GUI (mostly handled in US_AnalysisBase)
36 
37  setWindowTitle( tr( "Enhanced van Holde - Weischet Analysis:" ) );
38 
39  dbg_level = US_Settings::us_debug();
40  pb_dstrpl = us_pushbutton( tr( "Distribution Plot" ) );
41  pb_dstrpl->setEnabled( false );
42  pb_selegr = us_pushbutton( tr( "Select Groups" ) );
43  groupSel = false;
44  pb_selegr->setEnabled( false );
45  us_checkbox( tr( "Plateaus from 2DSA" ), ck_modelpl, true );
46  pb_replot = us_pushbutton( tr( "Refresh Plot" ) );
47  pb_replot->setEnabled( false );
48  us_checkbox( tr( "Manual-only Replot" ), ck_manrepl, false );
49  us_checkbox( tr( "Use Enhanced vHW" ), ck_vhw_enh, true );
50  us_checkbox( tr( "Use FE Data" ), ck_use_fed, false );
51 
52  connect( pb_dstrpl, SIGNAL( clicked() ),
53  this, SLOT( distr_plot() ) );
54  connect( pb_selegr, SIGNAL( clicked() ),
55  this, SLOT( sel_groups() ) );
56  connect( pb_replot, SIGNAL( clicked() ),
57  this, SLOT( plot_refresh() ) );
58  connect( ck_modelpl, SIGNAL( toggled( bool ) ),
59  this, SLOT( data_plot() ) );
60  connect( ck_vhw_enh, SIGNAL( toggled( bool ) ),
61  this, SLOT( data_plot() ) );
62  connect( ck_use_fed, SIGNAL( toggled( bool ) ),
63  this, SLOT( data_plot() ) );
64  connect( pb_save, SIGNAL( clicked() ),
65  this, SLOT( save_data() ) );
66  connect( pb_view, SIGNAL( clicked() ),
67  this, SLOT( view_report() ) );
68 
69  int jr = 2;
70  parameterLayout->addWidget( pb_dstrpl, jr, 0, 1, 2 );
71  parameterLayout->addWidget( pb_replot, jr, 2, 1, 1 );
72  parameterLayout->addWidget( pb_selegr, jr++, 3, 1, 1 );
73  parameterLayout->addWidget( ck_use_fed, jr, 0, 1, 2 );
74  parameterLayout->addWidget( ck_manrepl, jr++, 2, 1, 2 );
75  parameterLayout->addWidget( ck_modelpl, jr, 0, 1, 2 );
76  parameterLayout->addWidget( ck_vhw_enh, jr++, 2, 1, 2 );
77 
78  QLabel* lb_analysis = us_banner( tr( "Analysis Controls" ) );
79  QLabel* lb_scan = us_banner( tr( "Scan Control" ) );
80  QLabel* lb_smoothing = us_label ( tr( "Data Smoothing:" ) );
81  QLabel* lb_boundPercent = us_label ( tr( "% of Boundary:" ) );
82  QLabel* lb_boundPos = us_label ( tr( "Boundary Position (%):" ) );
83 
84  QLabel* lb_from = us_label ( tr( "Scan focus from:" ) );
85  QLabel* lb_to = us_label ( tr( "to:" ) );
86 
87  lb_tolerance = us_label( tr( "Back Diffusion Tolerance:" ) );
88  lb_tolerance->setAlignment( Qt::AlignVCenter | Qt::AlignLeft );
89 
90  ct_tolerance = us_counter( 3, 0.001, 1.0, 0.001 );
91  bdtoler = 0.001;
92  ct_tolerance->setStep( 0.001 );
93  ct_tolerance->setEnabled( true );
94  connect( ct_tolerance, SIGNAL( valueChanged( double ) ),
95  this, SLOT( update_bdtoler( double ) ) );
96 
97  lb_division = us_label( tr( "Divisions:" ) );
98  lb_division->setAlignment( Qt::AlignVCenter | Qt::AlignLeft );
99 
100  ct_division = us_counter( 3, 0.0, 1000.0, 50.0 );
101  ct_division->setStep( 1 );
102  ct_division->setEnabled( true );
103  connect( ct_division, SIGNAL( valueChanged( double ) ),
104  this, SLOT( update_divis( double ) ) );
105 
106  jr = 0;
107  controlsLayout->addWidget( lb_analysis , jr++, 0, 1, 4 );
108  controlsLayout->addWidget( lb_tolerance , jr, 0, 1, 2 );
109  controlsLayout->addWidget( ct_tolerance , jr++, 2, 1, 2 );
110  controlsLayout->addWidget( lb_division , jr, 0, 1, 2 );
111  controlsLayout->addWidget( ct_division , jr++, 2, 1, 2 );
112  controlsLayout->addWidget( lb_smoothing , jr, 0, 1, 2 );
113  controlsLayout->addWidget( ct_smoothing , jr++, 2, 1, 2 );
114  controlsLayout->addWidget( lb_boundPercent , jr, 0, 1, 2 );
115  controlsLayout->addWidget( ct_boundaryPercent, jr++, 2, 1, 2 );
116  controlsLayout->addWidget( lb_boundPos , jr, 0, 1, 2 );
117  controlsLayout->addWidget( ct_boundaryPos , jr++, 2, 1, 2 );
118  controlsLayout->addWidget( lb_scan , jr++, 0, 1, 4 );
119  controlsLayout->addWidget( lb_from , jr, 0, 1, 2 );
120  controlsLayout->addWidget( ct_from , jr++, 2, 1, 2 );
121  controlsLayout->addWidget( lb_to , jr, 0, 1, 2 );
122  controlsLayout->addWidget( ct_to , jr++, 2, 1, 2 );
123  controlsLayout->addWidget( pb_exclude , jr, 0, 1, 2 );
124  controlsLayout->addWidget( pb_reset_exclude , jr++, 2, 1, 2 );
125 
126  connect( pb_help, SIGNAL( clicked() ),
127  this, SLOT( help() ) );
128 
129  dataLoaded = false;
130  haveZone = false;
131  forcePlot = false;
132  skipPlot = false;
133 
134  rightLayout->setStretchFactor( plotLayout1, 3 );
135  rightLayout->setStretchFactor( plotLayout2, 2 );
136 
137  setMaximumSize( qApp->desktop()->size() - QSize( 80, 80 ) );
138  //resize( 100, 100 );
139  //resize( 900, 100 );
140  adjustSize();
141 }
142 
143 // load data
144 void US_vHW_Enhanced::load( void )
145 {
146  skipPlot = true;
147  US_AnalysisBase2::load();
148  skipPlot = false;
149 
150  if ( ! dataLoaded )
151  return;
152 
153  if ( dataList[ 0 ].expType == "Equilibrium" )
154  {
155  QMessageBox::warning( this, tr( "Wrong Data Type" ),
156  tr( "You have selected Equilibrium data, which is not\n"
157  "appropriate for van Holde - Weischet analysis." ) );
158 
159  dataLoaded = false;
160  lw_triples->disconnect();
161  le_id ->disconnect();
162  le_temp ->disconnect();
163  te_desc ->disconnect();
164 
165  lw_triples->clear();
166  le_id ->clear();
167  le_temp ->clear();
168  te_desc ->clear();
169 
170  data_plot1->clear();
171  data_plot2->clear();
172 
173  dataList.clear();
174  rawList .clear();
175  triples .clear();
176  pb_exclude->setEnabled( false );
177  return;
178  }
179 
180  data_plot1->setCanvasBackground( Qt::black );
181  data_plot2->setCanvasBackground( Qt::black );
182  int bord = height() - data_plot1->height() - data_plot2->height() + 12;
183  int mxht = qApp->desktop()->height() - bord;
184  int p1ht = ( mxht * 400 ) / 650;
185  int p2ht = mxht - p1ht;
186  p1ht = qMin( p1ht, 400 );
187  p2ht = qMin( p2ht, 250 );
188 DbgLv(1) << "vhw: mxht p1ht p2ht" << mxht << p1ht << p2ht;
189  data_plot1->setMinimumSize( 600, p1ht );
190  data_plot2->setMinimumSize( 600, p2ht );
191  data_plot2->setAxisAutoScale( QwtPlot::yLeft );
192  data_plot2->setAxisAutoScale( QwtPlot::xBottom );
193 
194  gpick = new US_PlotPicker( data_plot1 );
195  gpick->setSelectionFlags( QwtPicker::PointSelection
196  | QwtPicker::ClickSelection );
197  connect( gpick, SIGNAL( mouseDown( const QwtDoublePoint& ) ),
198  this, SLOT( groupClick( const QwtDoublePoint& ) ) );
199  groupstep = NONE;
200 
201  pb_selegr->setEnabled( true );
202  pb_dstrpl->setEnabled( true );
203  pb_replot->setEnabled( true );
204  haveZone = false;
205 
206  // Set saved flags; initialize simulations lists
207  saved .clear();
208  have_sims.clear();
209  dsimList .clear();
210  for ( int ii = 0; ii < triples.size(); ii++ )
211  {
212  saved << false;
213  have_sims << false;
214  US_DataIO::RawData simraw;
215  US_DataIO::EditedData simdat;
216  US_Model model;
217 
218  simdat = dataList[ ii ];
219 DbgLv(1) << "vhw: init: ii" << ii;
220  US_AstfemMath::initSimData( simraw, simdat, 0.0 );
221 DbgLv(1) << "vhw: init: init'd simraw (s x p)"
222  << simraw.scanCount() << simraw.pointCount();
223 DbgLv(1) << "vhw: init: copy'd simdat (s x p)"
224  << simdat.scanCount() << simdat.pointCount();
225 
226  for ( int jj = 0; jj < simdat.scanCount(); jj++ )
227  {
228  for ( int kk = 0; kk < simdat.pointCount(); kk++ )
229  {
230  simdat.setValue( jj, kk, simraw.value( jj, kk ) );
231  }
232  }
233 
234  dsimList << simdat;
235  modlList << model;
236  }
237 
238  ck_use_fed->disconnect();
239  ck_use_fed->setChecked( false );
240  connect( ck_use_fed, SIGNAL( toggled( bool ) ),
241  this, SLOT( data_plot() ) );
242 
243  update( 0 );
244 }
245 
246 // distribution plot
247 void US_vHW_Enhanced::distr_plot( void )
248 {
249  QVector< double > bfracs;
250  double bterm = 100.0 * boundPct / (double)divsCount;
251  double pterm = 100.0 * positPct + bterm;
252  int npoints = dseds.size();
253 
254  bfracs.resize( npoints );
255 
256  for ( int jj = 0; jj < npoints; jj++ )
257  {
258  bfracs[ jj ] = pterm + bterm * (double)( jj );
259  }
260 
261  US_DistribPlot* dialog = new US_DistribPlot( bfracs, dseds, total_conc );
262  dialog->move( this->pos() + QPoint( 100, 100 ) );
263  dialog->exec();
264  delete dialog;
265  row = lw_triples->currentRow();
266  saved[ row ] = true;
267 }
268 
269 // data plot
270 void US_vHW_Enhanced::data_plot( void )
271 {
272  double xmax = 0.0;
273  double ymax = 0.0;
274  int count = 0;
275  int totalCount;
276 
277 DbgLv(1) << " data_plot: dataLoaded" << dataLoaded << "vbar" << vbar;
278  if ( !dataLoaded || vbar <= 0.0 || skipPlot )
279  return;
280 
281  if ( !forcePlot && ck_manrepl->isChecked() )
282  return;
283 
284  QApplication::setOverrideCursor( QCursor( Qt::WaitCursor ) );
285  row = lw_triples->currentRow();
286  edata = &dataList[ row ];
287  scanCount = edata->scanCount();
288  valueCount = edata->pointCount();
289 
290  if ( have_model() )
291  {
292  if ( ! have_sims[ row ] )
293  {
294  create_simulation();
295  }
296 
297  edata = ck_use_fed->isChecked() ? &dsimList[ row ] : edata;
298  }
299 
300  data_plot2->detachItems();
301  data_plot2->setAxisAutoScale( QwtPlot::yLeft );
302  data_plot2->setAxisAutoScale( QwtPlot::xBottom );
303  boundPct = ct_boundaryPercent->value() / 100.0;
304  positPct = ct_boundaryPos ->value() / 100.0;
305  baseline = calc_baseline();
306 
307  divsCount = qRound( ct_division->value() );
308  totalCount = scanCount * divsCount;
309  QVector< double > pxvec( scanCount );
310  QVector< double > pyvec( totalCount );
311  double* ptx = pxvec.data(); // t,log(p) points
312  double* pty = pyvec.data();
313  divfac = 1.0 / (double)divsCount;
314  correc = solution.s20w_correction * 1.0e13;
315  omega = edata->scanData[ 0 ].rpm * M_PI / 30.0;
316 //*TIMING
317 kcalls[0]+=1;QDateTime sttime=QDateTime::currentDateTime();
318 //*TIMING
319 
320 
321  // Get live scans and original plateaus
322  live_scans();
323 DbgLv(1) << " valueCount totalCount" << valueCount << totalCount;
324 DbgLv(1) << " scanCount divsCount" << scanCount << divsCount;
325 DbgLv(1) << " lscnCount" << lscnCount;
326 
327  mo_plats = ck_modelpl->isChecked() && have_model();
328  vhw_enh = ck_vhw_enh->isChecked();
329 
330  if ( mo_plats )
331  { // Calculate plateaus from a model
332  model_plateaus();
333  }
334 
335  else
336  { // Calculate plateaus from fitting to specified values
337  fitted_plateaus();
338  }
339 
340  // Do the lower plot
341  plot_data2();
342 
343  // Then set up to handle the upper (vHW Extrapolation) plot.
344  // Calculate the division-1 sedimentation coefficient intercept and,
345  // from that, the back diffusion coefficient
346 
347  bdiff_sedc = sedcoeff_intercept();
348 
349  bdiff_coef = back_diff_coeff( bdiff_sedc );
350 
351  init_partials();
352 
353  if ( vhw_enh )
354  { // Calculate all plot points using enhanced method
355  vhw_calcs_enhanced( ptx, pty );
356  }
357 
358  else
359  { // Calculate all plot points using standard method
360  vhw_calcs_standard( ptx, pty );
361  }
362 
363  // Remove points that are singles in a division
364  for ( int jj = 0; jj < divsCount; jj++ )
365  { // Examine each division
366  int count = 0;
367  for ( int ii = 0; ii < lscnCount; ii++ )
368  { // Count unexcluded scan points in a division
369  int kk = divsCount * ii + jj;
370  double xv = ptx[ ii ];
371  double yv = pty[ kk ];
372  if ( xv > 0.0 && yv > 0.0 ) count++;
373  }
374 
375  if ( count == 1 )
376  { // Mark a division with a single point as having none
377  int kk = jj;
378  for ( int ii = 0; ii < lscnCount; ii++, kk += divsCount )
379  pty[ kk ] = SEDC_NOVAL;
380  }
381  }
382 
383  // Draw the vHW extrapolation plot
384  data_plot1->clear();
385  us_grid( data_plot1 );
386 
387  QString ttrip = edata->cell + "/" + edata->channel
388  + "/" + edata->wavelength;
389 
390  data_plot1->setTitle( edata->runID + " (" + ttrip + ")\n"
391  + tr( "vHW Extrapolation Plot" ) );
392 
393  data_plot1->setAxisTitle( QwtPlot::xBottom, tr( "(Time)^-0.5" ) );
394  data_plot1->setAxisTitle( QwtPlot::yLeft ,
395  tr( "Corrected Sed. Coeff. (1e-13 s)" ) );
396 
397  int nxy = ( lscnCount > divsCount ) ? lscnCount : divsCount;
398  QVector< double > xvec( nxy );
399  QVector< double > yvec( nxy );
400  double* x = xvec.data();
401  double* y = yvec.data();
402 
403  QwtPlotCurve* curve;
404  QwtSymbol sym;
405  sym.setStyle( QwtSymbol::Ellipse );
406  sym.setPen ( QPen( Qt::blue ) );
407  sym.setBrush( QBrush( Qt::white ) );
408  sym.setSize ( 8 );
409 
410  int kk = 0; // Index to sed. coeff. values
411 
412  // Set points for each division of each scan
413 
414  for ( int ii = 0; ii < lscnCount; ii++ )
415  {
416  count = 0;
417  double xv = ptx[ ii ]; // Reciprocal square root of time value
418 
419  for ( int jj = 0; jj < divsCount; jj++ )
420  {
421  double yv = pty[ kk++ ]; // Sed.coeff. values for divs in scan
422  if ( xv >= 0.0 && yv >= 0.0 )
423  { // Points in a scan
424  x[ count ] = xv;
425  y[ count ] = yv;
426  xmax = ( xmax > xv ) ? xmax : xv;
427  count++;
428  }
429  }
430 
431  if ( count > 0 )
432  { // Plot the points in a scan
433  curve = us_curve( data_plot1,
434  tr( "Sed Coeff Points, scan %1" ).arg( ii+1 ) );
435 
436  curve->setStyle ( QwtPlotCurve::NoCurve );
437  curve->setSymbol( sym );
438  curve->setData ( x, y, count );
439  }
440  }
441 
442  double slope = 0.0;
443  double intcept = 0.0;
444 
445  // Fit lines for each division to all scan points
446 
447  for ( int jj = 0; jj < divsCount; jj++ )
448  { // Walk thru divisions, fitting line to points from all scans
449  count = 0;
450 
451  for ( int ii = 0; ii < lscnCount; ii++ )
452  {
453  kk = ii * divsCount + jj; // Sed. coeff. index
454 
455  if ( ptx[ ii ] > 0.0 && pty[ kk ] > 0.0 )
456  { // Points for scans in a division
457  x[ count ] = ptx[ ii ];
458  y[ count ] = pty[ kk ];
459  count++;
460  int js = liveScans[ ii ];
461  omega = edata->scanData[ js ].rpm * M_PI / 30.0;
462  }
463  }
464 
465  if ( count > 1 )
466  { // Fit a line to the scan points in a division
467 if(jj<3||jj>(divsCount-9))
468 DbgLv(1) << "plot2 jj count" << jj << count << " y0 yn" << y[0] << y[count-1];
469  double sigma = 0.0;
470  double correl;
471 
472  US_Math2::linefit( &x, &y, &slope, &intcept, &sigma, &correl, count );
473 
474  x[ 0 ] = 0.0; // X from 0.0 to max
475  x[ 1 ] = xmax + 0.001;
476  y[ 0 ] = intcept; // Y from intercept to y at x-max
477  y[ 1 ] = y[ 0 ] + x[ 1 ] * slope;
478 
479  if ( y[ 1 ] < 0.0 )
480  {
481  y[ 1 ] = 0.0;
482  x[ 1 ] = ( slope != 0.0 ? ( -y[ 0 ] / slope ) : xmax );
483  }
484 if(jj<3||jj>(divsCount-9))
485 DbgLv(1) << "plot2 jj" << jj << "x0 y0 x1 y1" << x[0] << y[0] << x[1] << y[1];
486 
487  curve = us_curve( data_plot1, tr( "Fitted Line %1" ).arg( jj ) );
488  curve->setPen( QPen( Qt::yellow ) );
489  curve->setData( x, y, 2 );
490  }
491  }
492 
493  // Set scales, then plot the points and lines
494  xmax *= 1.05;
495  xmax = (double)qRound( ( xmax + 0.0009 ) / 0.001 ) * 0.001;
496  ymax = (double)qRound( ( ymax + 0.3900 ) / 0.400 ) * 0.400;
497 // data_plot1->setAxisScale( QwtPlot::xBottom, 0.0, xmax, 0.005 );
498  data_plot1->setAxisAutoScale( QwtPlot::xBottom);
499  data_plot1->setAxisAutoScale( QwtPlot::yLeft );
500 
501  count = 0;
502 
503  for ( int ii = 0; ii < lscnCount; ii++ )
504  { // Accumulate the points of the back-diffusion cutoff line
505  x[ count ] = bdrads.at( count );
506  y[ count ] = bdcons.at( count );
507 DbgLv(2) << " bd x y k " << x[count] << y[count] << count+1;
508  count++;
509  }
510 
511  // Save all sedcoeff values for report files and plots
512  aseds.clear();
513  for ( int ii = 0; ii < totalCount; ii++ )
514  aseds << pty[ ii ];
515 
516  // Plot the red back-diffusion cutoff line
517  dcurve = us_curve( data_plot2, tr( "Fitted Line BD" ) );
518  dcurve->setPen( QPen( QBrush( Qt::red ), 3.0 ) );
519  dcurve->setData( x, y, count );
520 DbgLv(1) << " DP: xr0 yr0 " << x[0] << y[0];
521 DbgLv(1) << " DP: xrN yrN " << x[count-1] << y[count-1] << count;
522  data_plot2->replot();
523 
524  // Plot any upper plot vertical excluded-scan lines
525  vert_exclude_lines();
526  data_plot1->replot();
527 
528 
529  QApplication::restoreOverrideCursor();
530 //*TIMING
531 QStringList cd;
532 cd << "plot" << "divseds" << "sedc_int" << "fitplat" << "modplat"
533  << "clcstd" << "clcenh" << "avgplat" << "livscns" << "inipart"
534  << "updmid" << "divsed1" << "divsed2" << "divsed3" << "divsed4"
535  << "divsed5" << "sedcoeff" << "findroot" << "none" << "none";
536 kmsecs[0]+=sttime.msecsTo(QDateTime::currentDateTime());
537 for(int mm=0;mm<20;mm++) DbgLv(1) << "TIMINGS: mm" << mm
538  << "kcalls kmsecs" << kcalls[mm] << kmsecs[mm] << cd[mm];
539 //*TIMING
540 }
541 
542 // Write the main report HTML to a stream
543 void US_vHW_Enhanced::write_report( QTextStream& ts )
544 {
545  edata = &dataList[ lw_triples->currentRow() ];
546  ts << html_header( "US_vHW_Enhanced",
547  tr( "van Holde - Weischet Analysis" ), edata );
548  ts << analysis( "" );
549 
550  ts << "\n" + indent( 4 ) + tr( "<h3>Selected Groups:</h3>\n" )
551  + indent( 4 ) + "<table>\n";
552  int ngrp = groupdat.size();
553 
554  if ( ngrp == 0 )
555  ts << table_row( tr( "Groups Selected:" ), tr( "NONE" ) );
556 
557  else
558  {
559  ts << table_row( tr( "Group:" ), tr( "Average S:" ),
560  tr( "Relative Amount:" ) );
561 
562  double tsed = 0.0;
563 
564  for ( int jj = 0; jj < ngrp; jj++ )
565  {
566  ts << table_row(
567  QString().sprintf( "%3d:", jj + 1 ),
568  QString().sprintf( "%6.2f", groupdat.at( jj ).sed ),
569  QString().sprintf( "(%5.2f %%)", groupdat.at( jj ).percent ) );
570  tsed += groupdat.at( jj ).sed * groupdat.at( jj ).percent;
571  }
572 
573  tsed = tsed / ( (double)ngrp * 100.0 );
574  Swavg = ( Swavg > 0.0 ) ? Swavg : ( tsed * 1.0e-13 );
575  }
576 
577  ts << indent( 4 ) + "</table>\n\n";
578 
579  ts << indent( 4 ) + "<br/><table>\n";
580  ts << table_row( tr( "Average S:" ),
581  QString::number( Swavg * 1.0e13 ) );
582  ts << table_row( tr( "Initial Concentration from plateau fit:" ),
583  QString::number( C0 ) + tr( " OD/fringes" ) );
584  ts << indent( 4 ) + "</table>\n";
585 
586  double sl;
587  double ci;
588  double sig;
589  double cor;
590  QVector< double > xvec( lscnCount );
591  QVector< double > yvec( lscnCount );
592  double* x = xvec.data();
593  double* y = yvec.data();
594  QString tscn;
595  QString tpla;
596 
597  for ( int ii = 0; ii < lscnCount; ii++ )
598  { // accumulate time,plateau pairs for line fit
599  int js = liveScans[ ii ];
600  dscan = &edata->scanData[ js ];
601  x[ ii ] = dscan->seconds;
602  y[ ii ] = scPlats[ ii ];
603  }
604 
605  US_Math2::linefit( &x, &y, &sl, &ci, &sig, &cor, lscnCount );
606  C0 = ( C0 == 0.0 ) ? ci : C0;
607 
608  ts << "\n" + indent( 4 ) + "<br/><table>\n";
609  ts << table_row( tr( "Initial Concentration:" ),
610  QString::number( ci ) + " OD" );
611  ts << table_row( tr( "Correlation Coefficient:" ),
612  QString::number( cor ) );
613  ts << table_row( tr( "Standard Deviation:" ),
614  QString::number( sig ) );
615  ts << table_row( tr( "Initial Concentration from exponential fit:" ),
616  QString::number( C0 ) + " OD" );
617  ts << indent( 4 ) + "</table>\n";
618 
619  ts << " </body>\n</html>\n";
620 }
621 
622 // save the enhanced data
623 void US_vHW_Enhanced::save_data( void )
624 {
625  QApplication::setOverrideCursor( QCursor( Qt::WaitCursor ) );
626 QDateTime time0=QDateTime::currentDateTime();
627  row = lw_triples->currentRow();
628  edata = &dataList[ row ];
629  QString tripl = QString( triples.at( row ) ).replace( " / ", "" );
630  QString basernam = US_Settings::resultDir() + "/" + edata->runID
631  + "/vHW." + tripl + ".";
632  QStringList files;
633  QStringList repfiles;
634 
635  // Write results files
636  write_vhw();
637  write_dis();
638  QString data0File = basernam + "extrap.csv";
639  QString data1File = basernam + "s-c-distrib.csv";
640  QString data2File = basernam + "s-c-envelope.csv";
641  files << data0File;
642  files << data1File;
643  files << data2File;
644 
645  if ( groupdat.size() > 0 )
646  {
647  write_model();
648  files << basernam + "fef_model.rpt";
649  }
650 
651  // Set up to write reports files
652  QString basename = US_Settings::reportDir() + "/" + edata->runID
653  + "/vHW." + tripl + ".";
654  QString htmlFile = basename + "report.html";
655  QString plot1File = basename + "velocity.svgz";
656  QString plot2File = basename + "extrap.svgz";
657  QString plot3File = basename + "s-c-distrib.svgz";
658  QString plot4File = basename + "s-c-histo.svgz";
659  QString dsinfFile = QString( basename ).replace( "/vHW.", "/dsinfo." )
660  + "dataset_info.html";
661 
662  // Write a general dataset information file
663  write_dset_report( dsinfFile );
664 
665  // Write the main report file
666  QFile rpt_f( htmlFile );
667 
668  if ( !rpt_f.open( QIODevice::WriteOnly | QIODevice::Text ) )
669  return;
670 
671  QTextStream ts( &rpt_f );
672  write_report( ts );
673  rpt_f.close();
674 
675  // Write the velocity and extrapolation plots
676  write_plot( plot1File, data_plot2 );
677  write_plot( plot2File, data_plot1 );
678 
679  // Save distribution and histogram plots
680  QVector< double > bfracs;
681  double bterm = 100.0 * boundPct / (double)divsCount;
682  double bfrac = 100.0 * positPct;
683  int npoints = dseds.size();
684  bfracs.resize( npoints );
685 
686  for ( int jj = 0; jj < npoints; jj++ )
687  {
688  bfrac += bterm;
689  bfracs[ jj ] = bfrac;
690  }
691 
692  if ( saved[ row ] == false )
693  { // If no plot dialog was opened for this triple, get files now
694  US_DistribPlot* dialog = new US_DistribPlot( bfracs, dseds, total_conc );
695 DbgLv(1) << "(P)PLOT ENV save: plot3File" << plot3File;
696  dialog->save_plots( plot3File, plot4File );
697 
698  delete dialog;
699  }
700 
701  else
702  { // Otherwise, copy the temporary files created earlier
703 DbgLv(1) << "(T)PLOT ENV save: plot3File" << plot3File;
704  copy_data_files( plot3File, plot4File, data2File );
705  }
706 
707  files << " ";
708  files << htmlFile;
709  files << plot1File;
710  files << plot2File;
711  files << plot3File;
712  files << plot4File;
713  files << dsinfFile;
714  update_filelist( repfiles, htmlFile );
715  update_filelist( repfiles, plot1File );
716  update_filelist( repfiles, plot2File );
717  update_filelist( repfiles, plot3File );
718  update_filelist( repfiles, plot4File );
719  update_filelist( repfiles, dsinfFile );
720  update_filelist( repfiles, data0File );
721  update_filelist( repfiles, data1File );
722  update_filelist( repfiles, data2File );
723 
724  // Report files created to the user
725  QString wmsg = tr( "In directories\n%1,\n%2;\nwrote:\n\n" )
726  .arg( files[ 0 ].left( files[ 0 ].lastIndexOf( "/" ) ) )
727  .arg( plot4File.left( plot4File.lastIndexOf( "/" ) ) );
728 
729  for ( int ii = 0; ii < files.size(); ii++ )
730  {
731  QString fname = files[ ii ];
732  fname = fname.mid( fname.lastIndexOf( "/" ) + 1 );
733  wmsg += " " + fname + "\n";
734  }
735 
736 QDateTime time1=QDateTime::currentDateTime();
737  if ( disk_controls->db() )
738  { // Also save report files to the database
739  reportFilesToDB( repfiles );
740 
741  wmsg += tr( "\nFiles were also saved to the database." );
742  }
743 QDateTime time2=QDateTime::currentDateTime();
744 int etim1=time0.msecsTo(time1);
745 int etim2=time1.msecsTo(time2);
746 int etimt=etim1+etim2;
747 int et1pc=(etim1*100)/etimt;
748 int et2pc=(etim2*100)/etimt;
749 DbgLv(1) << "SAVE-FILES: local ms" << etim1 << "=" << et1pc << "%";
750 DbgLv(1) << "SAVE-FILES: DB ms" << etim2 << "=" << et2pc << "%";
751  QApplication::restoreOverrideCursor();
752 
753  QMessageBox::information( this, tr( "Successfully Written" ), wmsg );
754 }
755 
756 // generate result report, pop up dialog and display the report
757 void US_vHW_Enhanced::view_report( void )
758 {
759  // Generate the main report html text
760  QString rtext;
761  QTextStream ts( &rtext );
762  write_report( ts );
763 
764  // Create the report dialog and display the text
765  US_Editor* edit = new US_Editor( US_Editor::LOAD, true );
766  edit->setWindowTitle( "Results: van Holde - Weischet Analysis" );
767  edit->move( this->pos() + QPoint( 100, 100 ) );
768  edit->resize( 700, 600 );
769  edit->e->setFont( QFont( US_GuiSettings::fontFamily(),
770  US_GuiSettings::fontSize() ) );
771  edit->e->setText( rtext );
772  edit->show();
773 }
774 
775 // alternate between selecting set-up and clearing vH-W groups
776 void US_vHW_Enhanced::sel_groups( void )
777 {
778  if ( groupSel )
779  { // had been Select, now doing Clear
780  groupstep = NONE;
781  pb_selegr->setText( tr( "Select Groups" ) );
782 
783  data_plot1->detachItems( QwtPlotItem::Rtti_PlotMarker );
784  data_plot1->replot();
785  groupdat.clear();
786  groupxy.clear();
787  }
788 
789  else
790  { // had been Clear, now doing new Select
791  groupstep = START;
792  pb_selegr->setText( tr( "Clear Groups" ) );
793  groupxy.clear();
794 
795  QMessageBox::information( this, tr( "Group Selection" ),
796  tr( "Please click first ABOVE, then BELOW the intercepts on the\n"
797  "Y-axis of the van Holde - Weischet extrapolation plot to\n"
798  "define groups and to average the S-values within a group.\n\n"
799  "PLEASE NOTE:\n"
800  "This algorithm provides for multiple and also overlapping\n"
801  "groups, such that the total percentage may exceed 100%" ) );
802  }
803 
804  groupSel = !groupSel; // reverse Select/Clear mode
805 }
806 
807 // update density
808 void US_vHW_Enhanced::update_density( double dval )
809 {
810  density = dval;
811 }
812 
813 // update viscosity
814 void US_vHW_Enhanced::update_viscosity( double dval )
815 {
816  viscosity = dval;
817 }
818 
819 // update vbar
820 void US_vHW_Enhanced::update_vbar( double dval )
821 {
822  vbar = dval;
823 }
824 
825 
826 void US_vHW_Enhanced::update_bdtoler( double dval )
827 {
828  bdtoler = dval;
829 
830  data_plot();
831 }
832 
833 void US_vHW_Enhanced::update_divis( double dval )
834 {
835  divsCount = qRound( dval );
836 
837  data_plot();
838 }
839 
840 // Index to first readings value greater than or equal to given concentration
841 int US_vHW_Enhanced::first_gteq( double concenv,
842  QVector< double >& rvalues, int valueCount, int defndx )
843 {
844  int index = defndx; // Set return index to default
845 
846  for ( int jj = 0; jj < valueCount; jj++ )
847  { // Find index where readings value equals or exceeds given concentration
848  if ( rvalues[ jj ] >= concenv )
849  {
850  index = jj; // Set index to found point
851  break; // And return with it
852  }
853  }
854  return index;
855 }
856 
857 // Index to first greater-or-equal value (with default of -1)
858 int US_vHW_Enhanced::first_gteq( double concenv,
859  QVector< double >& rvalues, int valueCount )
860 {
861  return first_gteq( concenv, rvalues, valueCount, -1 );
862 }
863 
864 // Get average scan plateau value for 41 points around user-specified value
865 double US_vHW_Enhanced::avg_plateau( )
866 {
867 //*TIMING
868 kcalls[7]+=1;QDateTime sttime=QDateTime::currentDateTime();
869 //*TIMING
870  double plato = 0.0;
871  int j2 = edata->xindex( edata->plateau ); // Index plateau radius
872  int j1 = max( 0, ( j2 - PA_POINTS ) ); // Point to 20 points before
873  j2 = min( valueCount, ( j2 + PA_POINTS + 1 ) );
874  // and 20 points after
875  for ( int jj = j1; jj < j2; jj++ ) // Sum 41 points centered at the
876  plato += dscan->rvalues[ jj ]; // scan plateau radial position
877 
878  plato /= (double)( j2 - j1 ); // Find and return the average
879 
880 //*TIMING
881 kmsecs[7]+=sttime.msecsTo(QDateTime::currentDateTime());
882 //*TIMING
883  return plato;
884 }
885 
886 // Get sedimentation coefficient for a given concentration
887 double US_vHW_Enhanced::sed_coeff( double cconc, double oterm,
888  double* radP, int* ndxP )
889 {
890 //*TIMING
891 kcalls[16]+=1;QDateTime sttime=QDateTime::currentDateTime();
892 //*TIMING
893  int j2 = first_gteq( cconc, dscan->rvalues, valueCount );
894  double rv0 = -1.0; // Mark radius excluded
895  double sedc = SEDC_NOVAL;
896 
897  if ( j2 >= 0 && oterm >= 0.0 )
898  { // Likely need to interpolate radius from two input values
899  int j1 = j2 - 1;
900 
901  if ( j2 > 0 )
902  { // Interpolate radius value
903  double av1 = dscan->rvalues[ j1 ];
904  double av2 = dscan->rvalues[ j2 ];
905  double rv1 = edata->xvalues[ j1 ];
906  double rv2 = edata->xvalues[ j2 ];
907  double rra = av2 - av1;
908  rra = ( rra == 0.0 ) ? 0.0 : ( ( rv2 - rv1 ) / rra );
909  rv0 = rv1 + ( cconc - av1 ) * rra;
910  j2 = ( av2 - cconc ) > ( cconc - av1 ) ? j1 : j2;
911  }
912 
913  else
914  {
915  rv0 = -1.0;
916  j2 = -1;
917  }
918  }
919 
920  if ( rv0 > 0.0 )
921  { // Use radius and other terms to get corrected sedimentation coeff. value
922  sedc = correc * log( rv0 / edata->meniscus ) / oterm;
923 //DbgLv(2) << "sed_coeff: rv0" << rv0 << " sedc" << sedc << " oterm" << oterm;
924  }
925 
926  if ( radP != NULL ) *radP = rv0;
927  if ( ndxP != NULL ) *ndxP = j2;
928 //*TIMING
929 kmsecs[16]+=sttime.msecsTo(QDateTime::currentDateTime());
930 //*TIMING
931  return sedc;
932 }
933 
934 // Calculate division sedimentation coefficient values (fitted line intercepts)
935 void US_vHW_Enhanced::div_seds( )
936 {
937 //*TIMING
938 kcalls[1]+=1;QDateTime sttime=QDateTime::currentDateTime();
939 //*TIMING
940  QVector< double > xxv( lscnCount );
941  QVector< double > yyv( lscnCount );
942  QVector< double > xrv( lscnCount );
943  QVector< double > yrv( lscnCount );
944  double* xx = xxv.data();
945  double* yy = yyv.data();
946  double* xr = xrv.data();
947  double* yr = yrv.data();
948  int nscnu = 0; // Number used (non-excluded) scans
949  int kscnu = 0; // Count of scans of div not affected by diffusion
950  double bdifcsqr = sqrt( bdiff_coef ); // Sqrt( diff_coeff ) used below
951  double pconc;
952  double mconc;
953  bdtoler = ct_tolerance->value();
954  valueCount = edata->pointCount();
955  double bottom = edata->radius( valueCount - 1 );
956 
957  // Do division-1 determination of base
958 
959 //*TIMING
960 kcalls[11]+=1;QDateTime sttim1=QDateTime::currentDateTime();
961 //*TIMING
962  for ( int ii = 0; ii < lscnCount; ii++ )
963  {
964  double timecor; // Time (corrected)
965  double timesqr; // Square root of corrected time
966  double bdleft; // Back-diffusion left value
967  double xbdleft; // Find-root of bdleft
968  double radD; // Radius of back-diffusion point
969  double omegasq; // Omega squared
970 
971  int js = liveScans[ ii ];
972 DbgLv(1) << " ii nscnu" << ii << nscnu << " js lscnCount" << js << lscnCount;
973  dscan = &edata->scanData[ js ];
974  omega = dscan->rpm * M_PI / 30.0;
975  omegasq = omega * omega;
976  timecor = dscan->seconds - time_correction;
977  timesqr = sqrt( timecor );
978  pconc = baseline + ( scPlats[ ii ] - baseline ) * positPct;
979 
980  xx[ ii ] = 1.0 / timesqr; // Save X (reciprocal sqrt(time))
981 
982  // Accumulate limits based on back diffusion
983 
984 //left=tolerance*pow(diff,0.5)/(2*intercept[0]*omega_s*
985 // (bottom+run_inf.meniscus[selected_cell])/2
986 // *pow(run_inf.time[selected_cell][selected_lambda][i],0.5);
987 //radD=bottom-(2*find_root(left)
988 // *pow((diff*run_inf.time[selected_cell][selected_lambda][i]),0.5));
989 
990  // left = tolerance * sqrt( diff )
991  // / ( 2 * intercept[0] * omega^2
992  // * ( bottom + meniscus ) / 2 * sqrt( time ) )
993 
994  bdleft = bdtoler * bdifcsqr
995  / ( bdiff_sedc * omegasq * ( bottom + edata->meniscus ) * timesqr );
996  xbdleft = find_root( bdleft );
997 
998  // radD = bottom - ( 2 * find_root(left) * sqrt( diff * time ) )
999 
1000  radD = bottom - ( 2.0 * xbdleft * bdifcsqr * timesqr );
1001  radD = qMax( edata->xvalues[ 0 ], qMin( bottom, radD ) );
1002 
1003  int mm = edata->xindex( radD ); // Radius's index
1004 
1005  // Accumulate for this scan of this division
1006  // the back diffusion limit radius and corresponding concentration
1007  xr[ nscnu ] = radD; // BD Radius
1008  yr[ nscnu ] = dscan->rvalues[ mm ]; // BD Concentration
1009 DbgLv(1) << " bottom meniscus bdleft" << bottom << edata->meniscus << bdleft;
1010 DbgLv(1) << " bdifsedc find_root toler" << bdiff_sedc << xbdleft << bdtoler;
1011 DbgLv(1) << "BD x,y " << nscnu+1 << radD << yr[nscnu] << " mm" << mm;
1012 DbgLv(1) << " ii nscnu" << ii << nscnu+1 << " bdsed radD yr"
1013  << bdiff_sedc << radD << yr[nscnu];
1014 
1015  nscnu++; // Bump scans-used count
1016  }
1017 
1018  dseds.clear();
1019  dslos.clear();
1020  dsigs.clear();
1021  dcors.clear();
1022  dpnts.clear();
1023 //*TIMING
1024 kmsecs[11]+=sttim1.msecsTo(QDateTime::currentDateTime());
1025 //*TIMING
1026 //*TIMING
1027 kcalls[12]+=1;QDateTime sttim2=QDateTime::currentDateTime();
1028 //*TIMING
1029 
1030  for ( int jj = 0; jj < divsCount; jj++ )
1031  { // Loop to fit points across scans in a division
1032 //*TIMING
1033 kcalls[13]+=1;QDateTime sttim3=QDateTime::currentDateTime();
1034 //*TIMING
1035  double dsed;
1036  double slope;
1037  double sigma;
1038  double corre;
1039  int ii;
1040  double oterm;
1041  double radC;
1042  double radD;
1043  double range;
1044  double rngFact = boundPct * divfac;
1045  double rngjFact = rngFact * jj;
1046  double conjFact = positPct + rngjFact;
1047 //DbgLv(1) << "div_sed div " << jj+1;
1048 
1049  kscnu = nscnu;
1050 //DbgLv(1) << "DVS jj" << jj << "nscnu" << nscnu;
1051 
1052  // Accumulate y values for this division, across used scans
1053 
1054  for ( int kk = 0; kk < nscnu; kk++ )
1055  { // Accumulate concentration, sed.coeff. for all scans, this div
1056  ii = liveScans[ kk ]; // Scan index
1057  dscan = &edata->scanData[ ii ]; // Scan pointer
1058 //DbgLv(1) << "DVS kk ii" << kk << ii;
1059 
1060  if ( vhw_enh )
1061  {
1062 //DbgLv(1) << "DVS CPijs sz" << CPijs.size();
1063 //DbgLv(1) << "DVS CPijs[0] sz" << CPijs[0].size();
1064 //DbgLv(1) << "DVS mconcs sz" << mconcs.size();
1065 //DbgLv(1) << "DVS mconcs[0] sz" << mconcs[0].size();
1066  cpij = CPijs [ kk ][ jj ]; // Partial concen (increment)
1067  mconc = mconcs[ kk ][ jj ]; // Mid-div concentration
1068  }
1069  else
1070  {
1071  range = scPlats[ kk ] - baseline; // Scan's range
1072  pconc = baseline + range * conjFact; // Base conc. of Div
1073  cpij = range * rngFact; // Partial concentration
1074  mconc = pconc + cpij * 0.5; // Mid-div concentration
1075  }
1076 
1077 //DbgLv(1) << "DVS kk jj" << kk << jj << "mconc cpij" << mconc << cpij;
1078  omega = dscan->rpm * M_PI / 30.0; // Omega
1079  oterm = ( dscan->seconds - time_correction ) * omega * omega;
1080 
1081  yy[ kk ] = sed_coeff( mconc, oterm, &radC, NULL ); // Sed.coeff.
1082  radD = xr[ kk ];
1083 
1084  if ( radC > radD )
1085  { // Gone beyond back-diffusion cutoff: exit loop with truncated list
1086  kscnu = kk;
1087  yy[ kk ] = SEDC_NOVAL;
1088 //DbgLv(1) << " div kscnu" << jj+1 << kscnu
1089 // << " radC radD" << radC << radD << " mconc sedcc" << mconc << yy[kk];
1090  break;
1091  }
1092  }
1093 //*TIMING
1094 kmsecs[13]+=sttim3.msecsTo(QDateTime::currentDateTime());
1095 //*TIMING
1096 //*TIMING
1097 kcalls[14]+=1;QDateTime sttim4=QDateTime::currentDateTime();
1098 //*TIMING
1099 
1100  int kscsv = kscnu;
1101  QVector< double > xtvec;
1102  QVector< double > ytvec;
1103  ii = 0;
1104 
1105  for ( int kk = 0; kk < kscnu; kk++ )
1106  { // Remove any leading points below meniscus
1107  if ( yy[ kk ] != SEDC_NOVAL )
1108  { // Sed coeff value not from below meniscus
1109  xtvec << xx[ kk ];
1110  ytvec << yy[ kk ];
1111  ii++;
1112  }
1113  else
1114  DbgLv(1) << "KS +++ SED<=0.0 divjj" << jj << "scnkk ii" << kk << ii;
1115  }
1116  kscnu = ii;
1117  double* xt = xtvec.data();
1118  double* yt = ytvec.data();
1119 DbgLv(1) << " KS" << kscnu << kscsv << "size xtv ytv" << xtvec.size()
1120  << ytvec.size();
1121 //*TIMING
1122 kmsecs[14]+=sttim4.msecsTo(QDateTime::currentDateTime());
1123 //*TIMING
1124 //*TIMING
1125 kcalls[15]+=1;QDateTime sttim5=QDateTime::currentDateTime();
1126 //*TIMING
1127 
1128  if ( kscnu > 1 )
1129  {
1130  // Calculate and save the division sedcoeff and fitted line slope
1131 
1132  US_Math2::linefit( &xt, &yt, &slope, &dsed, &sigma, &corre, kscnu );
1133 DbgLv(1) << " KS" << kscnu << kscsv << "jj" << jj << "xx0 xxn yy0 yyn" << xt[0]
1134  << xt[kscnu-1] << yt[0] << yt[kscnu-1] << "slo sed" << slope << dsed;
1135  }
1136 
1137  else
1138  {
1139  dsed = yyv[ 0 ];
1140  slope = 0.0;
1141  sigma = 0.0;
1142  corre = 0.0;
1143  slope = 0.0;
1144 ii=(kscnu<1)?1:kscnu;
1145 DbgLv(1) << " kscnu" << kscnu << "yy0 yyn " << yy[0] << yy[ii-1];
1146  continue;
1147  }
1148 
1149  dseds << dsed; // Save fitted line intercept (Sed.Coeff.)
1150  dslos << slope; // Save slope and other fitting values
1151  dsigs << sigma;
1152  dcors << corre;
1153  dpnts << kscnu;
1154 DbgLv(1) << "JJ" << jj << "DSED" << dsed << "dseds size" << dseds.size();
1155 //*TIMING
1156 kmsecs[15]+=sttim5.msecsTo(QDateTime::currentDateTime());
1157 //*TIMING
1158 
1159  }
1160  int kdivs = dseds.size();
1161 int k3=qMax(0,kdivs-3);
1162 int k2=qMax(0,kdivs-2);
1163 int k1=qMax(0,kdivs-1);
1164 DbgLv(1) << " dsed[0] " << dseds[0] << "kdivs" << kdivs;
1165 DbgLv(1) << " dsed[k3] " << dseds[k3];
1166 DbgLv(1) << " dsed[k2] " << dseds[k2];
1167 DbgLv(1) << " dsed[L] " << dseds[k1];
1168 DbgLv(1) << " D_S: xr0 yr0 " << xr[0] << yr[0];
1169 DbgLv(1) << " D_S: xrN yrN " << xr[nscnu-1] << yr[nscnu-1] << nscnu;
1170 
1171  bdrads.clear();
1172  bdcons.clear();
1173 
1174  for ( int kk = 0; kk < nscnu; kk++ )
1175  { // Save the Radius,Concentration points for back-diffusion cutoff curve
1176  bdrads << xr[ kk ];
1177  bdcons << yr[ kk ];
1178  }
1179 
1180 //*TIMING
1181 kmsecs[12]+=sttim2.msecsTo(QDateTime::currentDateTime());
1182 //*TIMING
1183 //*TIMING
1184 kmsecs[1]+=sttime.msecsTo(QDateTime::currentDateTime());
1185 //*TIMING
1186  return;
1187 }
1188 
1189 // Find root X where evaluated Y is virtually equal to a goal, using a
1190 // calculation including the inverse complementary error function (erfc).
1191 double US_vHW_Enhanced::find_root( double goal )
1192 {
1193 //*TIMING
1194 kcalls[17]+=1;QDateTime sttime=QDateTime::currentDateTime();
1195 //*TIMING
1196 #ifdef WIN32
1197 #define erfc(x) US_Math2::erfc(x)
1198 #endif
1199 
1200 #define _FR_MXKNT 100 // Max find-root iteration count
1201  double tolerance = 1.0e-7; // Min difference tolerance
1202  double x1 = 0.0;
1203  double x2 = 10.0;
1204  double xv = 5.0;
1205  double xdiff = 2.5;
1206  double xsqr = xv * xv;
1207  double rsqr_pi = 1.0 / sqrt( M_PI );
1208  double test = exp( -xsqr ) * rsqr_pi - ( xv * erfc( xv ) );
1209  test = ( goal != 0.0 ) ? test : 0.0;
1210 DbgLv(2) << " find_root: goal test" << goal << test << " xv" << xv;
1211 //DbgLv(2) << " erfc(x)" << erfc(xv);
1212  int count = 0;
1213 
1214  // Iterate until the difference between subsequent x value evaluations
1215  // is too small to be relevant (or max count reached);
1216 
1217  while ( fabs( test - goal ) > tolerance )
1218  {
1219  xdiff = ( x2 - x1 ) / 2.0;
1220 
1221  if ( test < goal )
1222  { // At less than goal, adjust top (x2) limit
1223  x2 = xv;
1224  xv -= xdiff;
1225  }
1226 
1227  else
1228  { // At greater than goal, adjust bottom (x1) limit
1229  x1 = xv;
1230  xv += xdiff;
1231  }
1232 
1233  // Then update the test y-value
1234  xsqr = xv * xv;
1235  test = ( 1.0 + 2.0 * xsqr ) * erfc( xv )
1236  - ( 2.0 * xv * exp( -xsqr ) ) * rsqr_pi;
1237 //DbgLv(2) << " find_root: goal test" << goal << test << " x" << xv;
1238 
1239  if ( (++count) > _FR_MXKNT )
1240  break;
1241  }
1242 DbgLv(2) << " find_root: goal test" << goal << test
1243  << " xv" << xv << " count" << count;
1244 
1245 //*TIMING
1246 kmsecs[17]+=sttime.msecsTo(QDateTime::currentDateTime());
1247 //*TIMING
1248  return xv;
1249 }
1250 
1251 // Calculate back diffusion coefficient
1252 double US_vHW_Enhanced::back_diff_coeff( double sedc )
1253 {
1254  double tempera = le_temp->text().section( " ", 0, 0 ).toDouble();
1255  double RT = R * ( K0 + tempera );
1256  double D1 = AVOGADRO * 0.06 * M_PI * viscosity;
1257  double D2 = 0.045 * sedc * vbar * viscosity;
1258  double D3 = 1.0 - vbar * density;
1259 
1260  double bdcoef = RT / ( D1 * sqrt( D2 / D3 ) );
1261 
1262 DbgLv(1) << "BackDiffusion:";
1263 DbgLv(1) << " RT " << RT << " R K0 tempera " << R << K0 << tempera;
1264 DbgLv(1) << " D1 " << D1 << " viscosity AVO " << viscosity << AVOGADRO;
1265 DbgLv(1) << " D2 " << D2 << " sedc vbar " << sedc << vbar;
1266 DbgLv(1) << " D3 " << D3 << " density " << density;
1267 DbgLv(1) << " bdiff_coef" << bdcoef << " = RT/(D1*sqrt(D2/D3))";
1268  return bdcoef;
1269 }
1270 
1271 // Handle mouse clicks for selecting vH-W groups
1272 void US_vHW_Enhanced::groupClick( const QwtDoublePoint& p )
1273 {
1274  QwtPlotMarker* marker;
1275  QwtText label;
1276  GrpInfo cgrdata;
1277  QString gbanner;
1278  int ngroup;
1279 DbgLv(1) << "groupClick: step" << groupstep
1280  << "x y" << p.x() << p.y();
1281 
1282  switch( groupstep )
1283  {
1284  case NONE :
1285  break;
1286 
1287  case START :
1288  groupstep = END;
1289  groupxy << p.x() << p.y(); // add start x,y to list
1290  break;
1291 
1292  case END :
1293  groupstep = START;
1294  groupxy << p.x() << p.y(); // add end x,y to list
1295  ngroup = groupxy.size() / 4;
1296 
1297  add_group_info( );
1298 
1299  marker = new QwtPlotMarker;
1300  cgrdata = groupdat.at( ngroup-1 );
1301 
1302  gbanner = tr( "Group %1: %2 (%3%)" )
1303  .arg( ngroup ).arg( cgrdata.sed ).arg( cgrdata.percent );
1304  label.setFont( QFont( US_GuiSettings::fontFamily(),
1305  US_GuiSettings::fontSize() + 2, QFont::Bold ) );
1306  label.setColor( Qt::darkRed );
1307  label.setBackgroundBrush( QBrush( QColor( 255, 255, 255, 208 ) ) );
1308  label.setText( gbanner );
1309 
1310  marker->setValue( 0.0, cgrdata.sed );
1311  marker->setLabel( label );
1312  marker->setLabelAlignment( Qt::AlignRight | Qt::AlignVCenter );
1313  marker->attach( data_plot1 );
1314 
1315  data_plot1->replot();
1316  break;
1317 
1318  default :
1319  break;
1320  }
1321 DbgLv(1) << "groupClick: nxy val" << groupxy.size();
1322 }
1323 
1324 // add to selected group information
1325 void US_vHW_Enhanced::add_group_info( )
1326 {
1327  int ngroup = groupxy.size() / 4;
1328  int jg = ( ngroup - 1 ) * 4;
1329  GrpInfo grdat;
1330 
1331  grdat.x1 = groupxy.at( jg ); // save pick coordinates
1332  grdat.y1 = groupxy.at( jg + 1 );
1333  grdat.x2 = groupxy.at( jg + 2 );
1334  grdat.y2 = groupxy.at( jg + 3 );
1335  grdat.sed = 0.0; // initialize other variables
1336  grdat.percent = 0.0;
1337  grdat.ndivs = 0;
1338  grdat.idivs.clear();
1339  int divsUsed = dseds.size();
1340 
1341  for ( int jj = 0; jj < divsUsed; jj++ )
1342  { // walk thru all division lines to see if within clicked range
1343  double sed = dseds[ jj ]; // intercept sed coeff
1344  double slope = dslos[ jj ]; // div line slope
1345  double yh = sed + slope * grdat.x1; // line intercept w high x
1346  double yl = sed + slope * grdat.x2; // line intercept w low x
1347 
1348  if ( yh <= grdat.y1 && yl >= grdat.y2 )
1349  { // line is in group: add division to group
1350  grdat.idivs.append( jj ); // add to div index list
1351  grdat.sed += sed; // accumulate sedcoeff sum for average
1352  grdat.ndivs++; // bump included divs count
1353  }
1354  }
1355 
1356  // finish average-sed-coeff and percent calculations; add to groups list
1357  grdat.ndivs = grdat.ndivs > 0 ? grdat.ndivs : 1;
1358  grdat.sed /= (double)grdat.ndivs;
1359  grdat.percent = ( (double)grdat.ndivs / (double)divsCount ) * 100.0;
1360 
1361  groupdat.append( grdat );
1362 }
1363 
1364 // write a file of vHW extrapolation data
1365 void US_vHW_Enhanced::write_vhw()
1366 {
1367  QString dirres = US_Settings::resultDir();
1368  QString dirrep = US_Settings::reportDir();
1369  QString dirname = dirres + "/" + edata->runID;
1370  QString dirrept = dirrep + "/" + edata->runID;
1371 
1372  QDir dirr( dirres );
1373  QDir dirp( dirrep );
1374 
1375  if ( ! dirr.exists( dirname ) )
1376  dirr.mkpath( dirname );
1377 
1378  if ( ! dirp.exists( dirrept ) )
1379  dirp.mkpath( dirrept );
1380 
1381  QString filename = dirname + "/vHW."
1382  + QString( triples.at( row ) ).replace( " / ", "" ) + ".extrap.csv";
1383 
1384  QFile res_f( filename );
1385  double sedc;
1386  int lastDiv = divsCount - 1;
1387  int kk = 0;
1388  const QString fsep( "," );
1389  const QString eoln( "\n" );
1390  const QString blnk( "\"\"" );
1391  //QString control = "\t";
1392  QString control = fsep;
1393 
1394  if ( !res_f.open( QIODevice::WriteOnly | QIODevice::Text ) )
1395  {
1396  return;
1397  }
1398 DbgLv(1) << "WV: filename " << filename;
1399 
1400  QTextStream ts( &res_f );
1401 
1402  // Output header line
1403  ts << "\"1/sqrt(t)\"" << fsep;
1404 
1405  for ( int jj = 0; jj < divsCount; jj++ )
1406  {
1407  QString line = QString().sprintf( "\"D%03dSedCoef\"", jj + 1 );
1408  ts << line << ( jj < lastDiv ? fsep : eoln );
1409  }
1410 
1411  // Output data
1412  for ( int ii = 0; ii < lscnCount; ii++ )
1413  {
1414  // Each output line begins with reciprocal square root of scan time
1415  int js = liveScans[ ii ];
1416  control = fsep;
1417  QString dat = QString().sprintf( "\"%11.8f\"",
1418  ( 1.0 / sqrt( edata->scanData[ js ].seconds - time_correction ) ) );
1419  dat.replace( " ", "" );
1420  ts << dat << control;
1421 
1422  // Balance of line is a list of sedimentation coefficient values for
1423  // the divisions in the scan
1424  for ( int jj = 0; jj < divsCount; jj++ )
1425  {
1426  sedc = aseds[ kk++ ];
1427  if ( sedc > 0 )
1428  {
1429  dat = QString().sprintf( "\"%8.5f\"", sedc );
1430  dat.replace( " ", "" );
1431  }
1432  else
1433  dat = blnk;
1434 
1435  ts << dat << ( jj < lastDiv ? control : eoln );
1436  }
1437  }
1438 
1439  res_f.close();
1440 }
1441 
1442 // write a file of vHW division distribution values
1443 void US_vHW_Enhanced::write_dis()
1444 {
1445  QString filename = US_Settings::resultDir() + "/" + edata->runID
1446  + "/vHW." + QString( triples.at( row ) ).replace( " / ", "" )
1447  + ".s-c-distrib.csv";
1448  QFile res_f( filename );
1449  double bterm = 100.0 * boundPct / (double)divsCount;
1450  double pterm = 100.0 * positPct + bterm;
1451  double bfrac;
1452  QString dline;
1453 
1454  if ( !res_f.open( QIODevice::WriteOnly | QIODevice::Text ) )
1455  {
1456  return;
1457  }
1458 DbgLv(1) << "WD: filename " << filename;
1459 
1460  // write the line-fit variables for each division
1461  QTextStream ts( &res_f );
1462  ts << edata->description << "\n";
1463  ts << tr( "\"%Boundary\",\"Points\",\"Slope\",\"Intercept\","
1464  "\"Sigma\",\"Correlation\"\n" );
1465 
1466  for ( int jj = 0; jj < dseds.size(); jj++ )
1467  {
1468  bfrac = pterm + bterm * (double)( jj );
1469 
1470  dline.sprintf(
1471  "\"%9.2f\",\"%7d\",\"%12.6f\",\"%12.6f\",\"%12.6f\",\"%12.6f\"\n",
1472  bfrac, dpnts[ jj ], dslos[ jj ], dseds[ jj ],
1473  dsigs[ jj ], dcors[ jj ] );
1474  dline.replace( " ", "" );
1475  ts << dline;
1476  }
1477 
1478  res_f.close();
1479 }
1480 
1481 // write a file of vHW detailed division group model data
1482 void US_vHW_Enhanced::write_model()
1483 {
1484  QString filename = US_Settings::resultDir() + "/" + edata->runID
1485  + "/vHW." + QString( triples.at( row ) ).replace( " / ", "" )
1486  + ".fef_model.rpt";
1487  QFile res_f( filename );
1488  int groups = groupdat.size();
1489 
1490  if ( !res_f.open( QIODevice::WriteOnly | QIODevice::Text ) )
1491  {
1492  return;
1493  }
1494 DbgLv(1) << "WM: filename " << filename;
1495 
1496  QTextStream ts( &res_f );
1497 
1498  ts << "*************************************\n";
1499  ts << tr( "* Please do not edit this file! *\n" );
1500  ts << "*************************************\n\n\n";
1501  ts << "3 " << tr( "\t# Fixed Molecular Weight Distribution\n" );
1502  ts << groups << " " << tr( "\t# Number of Components\n" );
1503  ts << edata->meniscus << " " << tr( "\t# Meniscus in cm\n" );
1504  ts << "0.01 " << tr( "\t# Meniscus range in cm\n" );
1505  ts << "0 " << tr( "\t# Meniscus fitting control\n" );
1506  ts << baseline << " " << tr( "\t# Baseline in OD\n" );
1507  ts << "0.01 " << tr( "\t# Baseline range in OD\n" );
1508  ts << "1 " << tr( "\t# Baseline fitting control\n" );
1509  ts << "0.0000 " << tr( "\t# Slope(r) Correction in OD\n" );
1510  ts << "0.0000 " << tr( "\t# Slope(r) Correction range in OD\n" );
1511  ts << "0 " << tr( "\t# Slope(r) Correction fitting control\n" );
1512  ts << "0.0000 " << tr( "\t# Slope(c,r) Correction in OD\n" );
1513  ts << "0.0000 " << tr( "\t# Slope(c,r) Correction range in OD\n" );
1514  ts << "0 " << tr( "\t# Slope(c,r) Correction fitting control\n");
1515  ts << "20 " << tr( "\t# Delta_t in seconds\n" );
1516  ts << "0.001 " << tr( "\t# Delta_t in cm\n" );
1517 
1518  QString line04 = "1 \t" +
1519  tr( "# Sedimentation Coefficient fitting control\n" );
1520  QString line05 = "7.00e-07 \t" +
1521  tr( "# Diffusion Coefficient in D\n" );
1522  QString line06 = "7.00e-08 \t" +
1523  tr( "# Diffusion Coefficient range in D\n" );
1524  QString line07 = "1 \t" +
1525  tr( "# Diffusion Coefficient fitting control\n" );
1526  QString line10 = "0 \t" +
1527  tr( "# Partial Concentration fitting control\n" );
1528  QString line11 = "0.0000 \t" + tr( "# Sigma\n" );
1529  QString line12 = "0.0000 \t" + tr( "# Sigma range\n" );
1530  QString line13 = "0 \t" + tr( "# Sigma fitting control\n" );
1531  QString line14 = "0.0000 \t" + tr( "# Delta\n" );
1532  QString line15 = "0.0000 \t" + tr( "# Delta range\n" );
1533  QString line16 = "0 \t" + tr( "# Delta fitting control\n" );
1534  QString line17 = "1 \t" +
1535  tr( "# Molecular Weight fitting control\n" );
1536  QString line18 = "0 \t" +
1537  tr( "# Part. Spec. Volume fitting control\n" );
1538  double cterm = ( C0 - baseline ) / 100.0;
1539 
1540  for ( int ii = 0; ii < groups; ii++ )
1541  {
1542  double gsed = groupdat.at( ii ).sed * 1.0e-13;
1543  double gcon = groupdat.at( ii ).percent * cterm;
1544  ts << "\n";
1545  ts << tr( "Parameters for Component " ) << ( ii + 1 ) << ":\n\n";
1546  ts << gsed << " "
1547  << tr( "\t# Sedimentation Coefficient in s\n" );
1548  ts << gsed / 10.0 << " "
1549  << tr( "\t# Sedimentation Coefficient range in s\n" );
1550  ts << line04 << line05 << line06 << line07;
1551  ts << gcon << " "
1552  << tr( "\t# Partial Concentration in OD\n" );
1553  ts << gcon / 10.0 << " "
1554  << tr( "\t# Partial Concentration range in OD\n" );
1555  ts << line10 << line11 << line12 << line13;
1556  ts << line14 << line15 << line16 << line17 << line18;
1557  }
1558 
1559  res_f.close();
1560 }
1561 
1562 // text of minutes,seconds or hours,minutes for a given total seconds value
1563 QString US_vHW_Enhanced::text_time( double seconds, int type )
1564 {
1565  int mins = (int)( seconds / 60.0 );
1566  int secs = (int)( seconds - (double)mins * 60.0 );
1567 
1568  if ( type == 0 )
1569  { // fixed-field mins,secs text
1570  QString tmin = QString().sprintf( "%4d", mins );
1571  QString tsec = QString().sprintf( "%3d", secs );
1572  return tr( "%1 min %2 sec" ).arg( tmin ).arg( tsec );
1573  }
1574 
1575  else if ( type == 1 )
1576  { // minutes,seconds text
1577  return tr( "%1 minutes(s) %2 second(s)" ).arg( mins ).arg( secs );
1578  }
1579 
1580  else
1581  { // hours,minutes text
1582  int hrs = (int)( seconds / 3600.0 );
1583  mins = qRound( ( seconds - (double)hrs * 3600.0 ) / 60.0 );
1584  return tr( "%1 hour(s) %2 minute(s)" ).arg( hrs ).arg( mins );
1585  }
1586 }
1587 
1588 // text of minutes and seconds for a given total seconds value
1589 QString US_vHW_Enhanced::text_time( double seconds )
1590 { // default mins,secs text
1591  return text_time( seconds, 0 );
1592 }
1593 
1594 
1595 // pare down files list by including only the last-edit versions
1596 QStringList US_vHW_Enhanced::last_edit_files( QStringList files )
1597 {
1598  QStringList ofiles;
1599  QStringList part;
1600  QString file;
1601  QString test;
1602  QString pfile;
1603  QString ptest;
1604  int nfi = files.size();
1605 
1606  // if only one in list, we need do no more
1607  if ( nfi < 2 )
1608  {
1609  return files;
1610  }
1611 
1612  // make sure files list is in ascending alphabetical order
1613  files.sort();
1614 
1615  // get first file name and its non-editID parts
1616  file = files[ 0 ];
1617  part = file.split( "." );
1618  test = part[ 0 ] + part[ 3 ] + part[ 4 ] + part[ 5 ];
1619 
1620  // skip all but last of any duplicates (differ only in editID)
1621  for ( int ii = 1; ii < nfi; ii++ )
1622  {
1623  pfile = file;
1624  ptest = test;
1625  file = files[ ii ];
1626  part = file.split( "." );
1627  test = part[ 0 ] + part[ 3 ] + part[ 4 ] + part[ 5 ];
1628 
1629  if ( QString::compare( test, ptest ) != 0 )
1630  { // differs by more than just edit, so output previous
1631  ofiles.append( pfile );
1632  }
1633  }
1634 
1635  // output the final
1636  ofiles.append( file );
1637 
1638  return ofiles;
1639 }
1640 
1641 void US_vHW_Enhanced::new_triple( int row )
1642 {
1643  US_AnalysisBase2::new_triple( row );
1644  haveZone = false;
1645 
1646  US_AnalysisBase2::smoothing( ct_smoothing->value() );
1647 }
1648 
1649 void US_vHW_Enhanced::update( int row )
1650 {
1651  haveZone = false;
1652  edata = &dataList[ row ];
1653 //*TIMING
1654 for(int mm=0;mm<20;mm++) { kcalls[mm]=0;kmsecs[mm]=0; }
1655 //*TIMING
1656 
1657  // Do some calculations handled in AnalysisBase, but needed here
1658  // so that model plateau calculations work
1659  time_correction = US_Math2::time_correction( dataList );
1660  solution.density = le_density ->text().toDouble();
1661  solution.viscosity = le_viscosity->text().toDouble();
1662  solution.vbar20 = le_vbar ->text().toDouble();
1663  double avgTemp = edata->average_temperature();
1664  solution.vbar = US_Math2::calcCommonVbar( solution_rec, avgTemp );
1665  US_Math2::data_correction( avgTemp, solution );
1666 
1667  // Do normal analysis triple updating, but suppress plotting for now
1668  skipPlot = true;
1669  US_AnalysisBase2::update( row );
1670  skipPlot = false;
1671 DbgLv(1) << " update: vbar" << vbar;
1672 
1673  if ( vbar <= 0.0 )
1674  {
1675  QMessageBox::warning( this, tr( "Bad Solution Values" ),
1676  tr( "The Vbar for this data is not a positive value (vbar=%1).\n"
1677  "The current data set needs editing of its Solution before\n"
1678  "van Holde - Weischet analysis can proceed." ).arg( vbar ) );
1679  return;
1680  }
1681 
1682  // After triple update has completed, we can proceed with plotting
1683  forcePlot = true;
1684  data_plot();
1685  forcePlot = false;
1686 
1687  // Report on whether a model is available
1688  QString tripl = QString( triples.at( row ) ).replace( " / ", "" );
1689  QString dmsg = te_desc->toPlainText() + "\n";
1690 
1691  if ( have_model() )
1692  {
1693  dmsg += tr( "A model IS implied for %1" ).arg( tripl );
1694  ck_modelpl->setEnabled( true );
1695  ck_use_fed->setEnabled( true );
1696  ck_modelpl->setChecked( true );
1697  }
1698  else
1699  {
1700  dmsg += tr( "NO model is implied for %1" ).arg( tripl );
1701  ck_modelpl->setEnabled( false );
1702  ck_use_fed->setEnabled( false );
1703  ck_modelpl->setChecked( false );
1704  }
1705 
1706  te_desc->setText( dmsg );
1707 
1708  if ( groupSel ) // Turn off Group Select if left from previous triple
1709  sel_groups();
1710 }
1711 
1712 // Calculate the sedimentation coefficient intercept for division 1
1713 double US_vHW_Enhanced::sedcoeff_intercept()
1714 {
1715 //*TIMING
1716 kcalls[2]+=1;QDateTime sttime=QDateTime::currentDateTime();
1717 //*TIMING
1718  QVector< double > xrvec( lscnCount );
1719  QVector< double > yrvec( lscnCount );
1720  double* xr = xrvec.data();
1721  double* yr = yrvec.data();
1722  double sedc = 0.0;
1723  int nscnu = 0;
1724 
1725  for ( int ii = 0; ii < lscnCount; ii++ )
1726  { // Accumulate x,y values: 1/sqrt(time), sed_coeff
1727  int js = liveScans[ ii ];
1728  dscan = &edata->scanData[ js ];
1729 
1730  double range = scPlats[ ii ] - baseline;
1731  double basecut = baseline + range * positPct;
1732  double platcut = basecut + range * boundPct;
1733  double cinc = ( platcut - basecut ) * divfac;
1734  double mconc = basecut + cinc * 0.5;
1735  double omega = dscan->rpm * M_PI / 30.0;
1736  double timecor = dscan->seconds - time_correction;
1737  double oterm = timecor * omega * omega;
1738 
1739  // Get sedimentation coefficient for concentration
1740  sedc = sed_coeff( mconc, oterm );
1741 
1742  if ( sedc > 0.0 )
1743  {
1744  xr[ nscnu ] = 1.0 / sqrt( timecor ); // X = inverse sqrt of time
1745  yr[ nscnu ] = sedc; // Y = sedimentation coeff.
1746  nscnu++; // bump count of used scans
1747  }
1748 DbgLv(2) << " s-i: range baseline basecut platcut mconc" << range << baseline
1749  << basecut << platcut << mconc;
1750  }
1751 
1752  double slope;
1753  double sigma = 0.0;
1754  double corre;
1755 
1756  // Fit a line. Use the intercept as the sedimentation coefficient
1757  US_Math2::linefit( &xr, &yr, &slope, &sedc, &sigma, &corre, nscnu );
1758 
1759  if ( sedc <= 0.0 )
1760  { // If it dropped below zero, back up line until it crosses to positive
1761 DbgLv(1) << "sedco-intcp **SEDC**" << sedc << "x0 y0 xn yn n"
1762  << xr[0] << yr[0] << xr[nscnu-1] << yr[nscnu-1] << nscnu;
1763 DbgLv(1) << "sedco-intcp slope sigma" << slope << sigma;
1764  while ( (--nscnu) > 0 && sedc <= 0.0 )
1765  {
1766  sedc = yr[ nscnu ];
1767  }
1768  }
1769 
1770 DbgLv(1) << "sedco-intcp ++SEDC++" << sedc << nscnu;
1771  sedc *= 1.0e-13;
1772 
1773 //*TIMING
1774 kmsecs[2]+=sttime.msecsTo(QDateTime::currentDateTime());
1775 //*TIMING
1776  return sedc;
1777 }
1778 
1779 // Copy temporary files saved by plot dialog to report/results
1780 void US_vHW_Enhanced::copy_data_files( QString plot1File,
1781  QString plot2File, QString data2File )
1782 {
1783  QString tempbase = US_Settings::tmpDir() + "/vHW.temp.";
1784  QString tplot1File = tempbase + "s-c-distrib.svgz";
1785  QString tplot2File = tempbase + "s-c-histo.svgz";
1786  QString tdata2File = tempbase + "s-c-envelope.csv";
1787  QString tplot3File = tempbase + "s-c-distrib.png";
1788  QString tplot4File = tempbase + "s-c-histo.png";
1789  QString plot3File = QString( plot1File ).section( ".", 0, -2 ) + ".png";
1790  QString plot4File = QString( plot2File ).section( ".", 0, -2 ) + ".png";
1791 
1792  QFile tfp1( tplot1File );
1793  QFile tfp2( tplot2File );
1794  QFile tfd2( tdata2File );
1795  QFile tfp3( tplot3File );
1796  QFile tfp4( tplot4File );
1797 
1798  if ( tfp1.exists() )
1799  {
1800  if ( QFile( plot1File ).remove() )
1801  DbgLv(1) << "CDF: removed:" << plot1File;
1802  if ( tfp1.copy( plot1File ) )
1803  DbgLv(1) << "CDF: copied:" << tplot1File;
1804  }
1805 
1806  if ( tfp2.exists() )
1807  {
1808  if ( QFile( plot2File ).remove() )
1809  DbgLv(1) << "CDF: removed:" << plot2File;
1810  if ( tfp2.copy( plot2File ) )
1811  DbgLv(1) << "CDF: copied:" << tplot2File;
1812  }
1813 
1814  if ( tfd2.exists() )
1815  {
1816  if ( QFile( data2File ).remove() )
1817  DbgLv(1) << "CDF: removed:" << data2File;
1818  if ( tfd2.copy( data2File ) )
1819  DbgLv(1) << "CDF: copied:" << tdata2File;
1820  }
1821 
1822  if ( tfp3.exists() )
1823  {
1824  if ( QFile( plot3File ).remove() )
1825  DbgLv(1) << "CDF: removed:" << plot3File;
1826  if ( tfp3.copy( plot3File ) )
1827  DbgLv(1) << "CDF: copied:" << tplot3File;
1828  }
1829 
1830  if ( tfp4.exists() )
1831  {
1832  if ( QFile( plot4File ).remove() )
1833  DbgLv(1) << "CDF: removed:" << plot4File;
1834  if ( tfp4.copy( plot4File ) )
1835  DbgLv(1) << "CDF: copied:" << tplot4File;
1836  }
1837 }
1838 
1839 // Calculate scan plateau concentrations by fitting initial values
1840 bool US_vHW_Enhanced::fitted_plateaus()
1841 {
1842 //*TIMING
1843 kcalls[3]+=1;QDateTime sttime=QDateTime::currentDateTime();
1844 //*TIMING
1845  int totalCount;
1846 DbgLv(1) << "FITTED_PLATEAUS";
1847  if ( !dataLoaded || vbar <= 0.0 )
1848  return false;
1849 
1850  valueCount = edata->pointCount();
1851  scanCount = edata->scanCount();
1852  divsCount = qRound( ct_division->value() );
1853  totalCount = scanCount * divsCount;
1854  divfac = 1.0 / (double)divsCount;
1855  boundPct = ct_boundaryPercent->value() / 100.0;
1856  positPct = ct_boundaryPos ->value() / 100.0;
1857  baseline = calc_baseline();
1858  correc = solution.s20w_correction * 1.0e13;
1859  C0 = 0.0;
1860  Swavg = 0.0;
1861  omega = edata->scanData[ 0 ].rpm * M_PI / 30.0;
1862 
1863  // Do first experimental plateau calcs based on horizontal zones
1864 
1865  int nrelp = 0;
1866  QVector< double > pxvec( lscnCount );
1867  QVector< double > pyvec( totalCount );
1868  QVector< double > sxvec( lscnCount );
1869  QVector< double > syvec( lscnCount );
1870  QVector< double > fyvec( lscnCount );
1871  double* ptx = pxvec.data(); // t,log(p) points
1872  double* pty = pyvec.data();
1873  double* csx = sxvec.data(); // count,slope points
1874  double* csy = syvec.data();
1875  double* fsy = fyvec.data(); // fitted slope
1876  double plateau;
1877 
1878  QList< double > plats;
1879  QList< int > isrel;
1880  QList< int > isunr;
1881 
1882  // accumulate reliable,unreliable scans and plateaus
1883 
1884  scPlats.fill( 0.0, lscnCount );
1885  nrelp = 0;
1886 
1887 DbgLv(1) << " Initial reliable (non-excluded) scan t,log(avg-plat):";
1888  // Initially reliable plateaus are all average non-excluded
1889  for ( int ii = 0; ii < lscnCount; ii++ )
1890  {
1891  int js = liveScans[ ii ];
1892  dscan = &edata->scanData[ js ];
1893  plateau = avg_plateau();
1894  scPlats[ ii ] = plateau;
1895  ptx[ nrelp ] = dscan->seconds - time_correction;
1896  pty[ nrelp ] = log( plateau );
1897 
1898  isrel.append( ii );
1899 DbgLv(1) << ptx[nrelp] << pty[nrelp];
1900  nrelp++;
1901  }
1902 
1903  double slope;
1904  double intcp;
1905  double sigma;
1906  double corre;
1907  int krelp = nrelp;
1908  int jrelp = nrelp;
1909  int kf = 0;
1910 
1911  // Accumulate count,slope points for line fits to the t,log(p) points
1912 
1913  for ( int jj = 6; jj < nrelp; jj++ )
1914  {
1915  US_Math2::linefit( &ptx, &pty, &slope, &intcp, &sigma, &corre, jj );
1916 
1917  csx[ kf ] = (double)jj;
1918  csy[ kf++ ] = slope;
1919 DbgLv(1) << "k,slope point: " << kf << jj << slope;
1920  }
1921 
1922  // Create a 3rd-order polynomial fit to the count,slope curve
1923 
1924  double coefs[ 4 ];
1925 
1926  US_Matrix::lsfit( coefs, csx, csy, kf, 4 );
1927 DbgLv(1) << "3rd-ord poly fit coefs:" << coefs[0] << coefs[1] << coefs[2]
1928  << coefs[3];
1929 
1930  for ( int jj = 0; jj < kf; jj++ )
1931  {
1932  double xx = csx[ jj ];
1933  double yy = coefs[ 0 ] + coefs[ 1 ] * xx + coefs[ 2 ] * xx * xx
1934  + coefs[ 3 ] * xx * xx * xx;
1935  fsy[ jj ] = yy;
1936 DbgLv(1) << " k,fit-slope point: " << jj << xx << yy << " raw slope" << csy[jj];
1937  }
1938 
1939  // Now find the spot where the derivative of the fit crosses zero
1940  double prevy = fsy[ 1 ];
1941  double dfac = 1.0 / qAbs( fsy[ kf - 1 ] - prevy ); // norm. deriv. factor
1942  double prevd = ( prevy - fsy[ 0 ] ) * dfac; // normalized derivative
1943  jrelp = -1;
1944  int kdecr = 0;
1945 
1946 DbgLv(1) << "Fitted slopes derivative points:";
1947  for ( int jj = 2; jj < kf; jj++ )
1948  {
1949  double currx = csx[ jj ];
1950  double curry = fsy[ jj ];
1951  double currd = ( curry - prevy ) * dfac; // normalized derivative
1952 //DbgLv(1) << " k cx cy cd pd" << jj << currx << curry << currd << prevd;
1953 DbgLv(1) << " k" << jj << " x fslo" << currx << curry << " deriv" << currd;
1954 
1955  if ( currd == 0.0 ||
1956  ( currd > 0.0 && prevd < 0.0 ) ||
1957  ( currd < 0.0 && prevd > 0.0 ) )
1958  { // Zero point or zero crossing
1959  jrelp = jj - 1;
1960 DbgLv(1) << " Z-CROSS jrelp" << jrelp;
1961 if (dbg_level>=1 ) { for ( int mm = jj + 1; mm < kf; mm++ ) {
1962  prevy = curry; prevd = currd;
1963  currx = csx[mm]; curry = fsy[mm];
1964  currd = ( curry - prevy ) * dfac;
1965  DbgLv(1) << " k cx cy cd pd" << mm << currx << curry << currd << prevd;
1966 }}
1967  break;
1968  }
1969 
1970  if ( prevy < 0.0 && curry < prevy )
1971  { // Count consecutive decreases in negative slopes
1972  kdecr++;
1973  if ( kdecr > 2 )
1974  { // 3rd consecutive decrease in negative slope, break at this point
1975  jrelp = jj;
1976 DbgLv(1) << " KDECR>2 jrelp prevy curry" << jrelp << prevy << curry;
1977  break;
1978  }
1979  }
1980 
1981  else
1982  { // Reset count of consecutive decreasing negative slopes
1983  kdecr = 0;
1984  }
1985 
1986  prevy = curry;
1987  prevd = currd;
1988  }
1989 
1990 DbgLv(1) << " JRELP" << jrelp << "KF" << kf;
1991  if ( jrelp < 0 || jrelp > ( kf - 2 ) )
1992  {
1993  jrelp = min( 2, ( kf - 1 ) );
1994 DbgLv(1) << " NO z-cross/k-decr, so default jrelp" << jrelp;
1995  }
1996 
1997 
1998  // Do final fit to the determined number of leading points
1999  krelp = (int)csx[ jrelp ];
2000 DbgLv(1) << " jrelp krelp" << jrelp << krelp;
2001 
2002  US_Math2::linefit( &ptx, &pty, &slope, &intcp, &sigma, &corre, krelp );
2003 DbgLv(1) << " KRELP" << krelp << " slope intcp sigma"
2004  << slope << intcp << sigma;
2005 
2006  Swavg = slope / ( -2.0 * omega * omega ); // Swavg func of slope
2007  C0 = exp( intcp ); // C0 func of intercept
2008  total_conc = C0;
2009 DbgLv(1) << "Swavg(c): " << Swavg*correc << " C0: " << C0 ;
2010 
2011  // Determine Cp for all the scans based on fitted line:
2012  // y = ax + b, using "a" and "b" determined above.
2013  // Since y = log( Cp ), we get Cp by exponentiating
2014  // the left-hand term.
2015 
2016  for ( int ii = 0; ii < lscnCount; ii++ )
2017  { // each extrapolated plateau is exp of Y for X of corrected time
2018  int js = liveScans[ ii ];
2019  dscan = &edata->scanData[ js ];
2020  double tc = dscan->seconds - time_correction;
2021  dscan->plateau = exp( tc * slope + intcp );
2022  scPlats[ ii ] = dscan->plateau;
2023 DbgLv(1) << " jj scan plateau " << ii << ii+1 << scPlats[ii];
2024  }
2025 //*TIMING
2026 kmsecs[3]+=sttime.msecsTo(QDateTime::currentDateTime());
2027 //*TIMING
2028 DbgLv(1) << "total_conc(fitted): " << total_conc;
2029  return true;
2030 }
2031 
2032 // Calculate scan plateau concentrations from a 2DSA model
2033 bool US_vHW_Enhanced::model_plateaus()
2034 {
2035 //*TIMING
2036 kcalls[4]+=1;QDateTime sttime=QDateTime::currentDateTime();
2037 //*TIMING
2038  get_model();
2039 
2040  int ncomp = model.components.size();
2041  double scorr = -2.0 / solution.s20w_correction;
2042 DbgLv(1) << "Mpl:ncmp" << ncomp << "scorr" << scorr << solution.s20w_correction;
2043  expda = &dataList[ row ];
2044  simda = &dsimList[ row ];
2045  edata = expda;
2046  valueCount = expda->pointCount();
2047  scanCount = expda->scanCount();
2048  scPlats.fill( 0.0, lscnCount );
2049  C0 = 0.0;
2050  Swavg = 0.0;
2051  total_conc = 0.0;
2052 
2053  for ( int ii = 0; ii < lscnCount; ii++ )
2054  {
2055  int js = liveScans[ ii ];
2056  expsc = &expda->scanData[ js ];
2057  simsc = &simda->scanData[ js ];
2058  double omega = expsc->rpm * M_PI / 30.0;
2059  double oterm = ( expsc->seconds - time_correction ) * omega * omega;
2060  oterm *= scorr;
2061  double cplat = 0.0;
2062 
2063  for ( int jj = 0; jj < ncomp; jj++ )
2064  {
2065  US_Model::SimulationComponent* sc = &model.components[ jj ];
2066  double conc = sc->signal_concentration;
2067  double sval = sc->s;
2068  cplat += ( conc * exp( oterm * sval ) );
2069 
2070  if ( ii == 0 )
2071  total_conc += conc;
2072  }
2073 //DbgLv(1) << "Mpl: scan" << ii << "cplat" << cplat;
2074 
2075  scPlats[ ii ] = cplat;
2076  expsc->plateau = cplat;
2077  simsc->plateau = cplat;
2078  C0 += cplat;
2079  }
2080 
2081  C0 = ( lscnCount < 1 ) ? 0.0 : ( C0 / (double)lscnCount );
2082 //*TIMING
2083 kmsecs[4]+=sttime.msecsTo(QDateTime::currentDateTime());
2084 //*TIMING
2085 DbgLv(1) << "total_conc(model): " << total_conc;
2086  return true;
2087 }
2088 
2089 // Do a manually forced re-plot of the data
2090 void US_vHW_Enhanced::plot_refresh()
2091 {
2092  forcePlot = true;
2093  data_plot();
2094  forcePlot = false;
2095 }
2096 
2097 // Calculate X,Y extrapolation points using standard method
2098 void US_vHW_Enhanced::vhw_calcs_standard( double* ptx, double* pty )
2099 {
2100 //*TIMING
2101 kcalls[5]+=1;QDateTime sttime=QDateTime::currentDateTime();
2102 //*TIMING
2103  int totalCount;
2104  double cconc;
2105  double pconc;
2106  double mconc;
2107  double cinc;
2108  double cinch;
2109  double oterm;
2110  divsCount = qRound( ct_division->value() );
2111  totalCount = scanCount * divsCount;
2112  divfac = 1.0 / (double)divsCount;
2113  C0 = 0.0;
2114  Swavg = 0.0;
2115  edata = ( ck_use_fed->isChecked() && have_model() )
2116  ? &dsimList[ row ] : &dataList[ row ];
2117 
2118  div_seds();
2119 
2120  for ( int ii = 0; ii < totalCount; ii++ )
2121  pty[ ii ] = SEDC_NOVAL;
2122 
2123  int kk = 0; // Index to sed. coeff. values
2124  int kl = 0; // Index/count of live scans
2125  valueCount = edata->pointCount();
2126 
2127  // Calculate the corrected sedimentation coefficients
2128 
2129  for ( int ii = 0; ii < lscnCount; ii++ )
2130  {
2131  range = scPlats[ ii ] - baseline;
2132  cinc = range * boundPct * divfac;
2133  int js = liveScans[ ii ];
2134  dscan = &edata->scanData[ js ];
2135  double timev = dscan->seconds - time_correction;
2136  double timex = 1.0 / sqrt( timev );
2137  double bdrad = bdrads.at( kl ); // Back-diffus cutoff radius for scan
2138  double bdcon = bdcons.at( kl++ ); // Back-diffus cutoff concentration
2139  double divrad = 0.0; // Division radius value
2140 DbgLv(1) << "scn liv" << ii+1 << kl
2141  << " radius concen time" << bdrad << bdcon << timev;
2142 
2143  ptx[ ii ] = timex; // Save corrected time and accum max
2144 
2145  range = scPlats[ ii ] - baseline;
2146  cconc = baseline + range * positPct; // Initial conc for span
2147  cinc = range * boundPct * divfac;
2148  cinch = cinc * 0.5;
2149  omega = dscan->rpm * M_PI / 30.0;
2150  oterm = ( timev > 0.0 ) ? ( timev * omega * omega ) : -1.0;
2151 
2152  for ( int jj = 0; jj < divsCount; jj++ )
2153  { // walk through division points; get sed. coeff. by place in readings
2154  pconc = cconc; // Div base
2155  cconc = pconc + cinc; // Absolute concentration
2156  mconc = pconc + cinch; // Mid div concentration
2157 
2158  sedc = sed_coeff( mconc, oterm, &divrad, NULL );
2159 
2160  if ( divrad > bdrad )
2161  { // Mark a point to be excluded by back-diffusion
2162  sedc = SEDC_NOVAL;
2163 DbgLv(1) << "CPm: *excl* div" << jj+1 << " drad dcon " << divrad << mconc;
2164  }
2165 DbgLv(1) << "CPm: ii jj" << ii << jj << "kk iSedc" << kk << sedc
2166  << "mconc divrad bdrad" << mconc << divrad << bdrad;
2167 
2168  // Y value of point is sedcoeff; accumulate y max
2169  pty[ kk++ ] = sedc;
2170  }
2171  }
2172 //*TIMING
2173 kmsecs[5]+=sttime.msecsTo(QDateTime::currentDateTime());
2174 //*TIMING
2175 }
2176 
2177 // Calculate X,Y extrapolation points using enhanced method
2178 void US_vHW_Enhanced::vhw_calcs_enhanced( double* ptx, double* pty )
2179 {
2180 //*TIMING
2181 kcalls[6]+=1;QDateTime sttime=QDateTime::currentDateTime();
2182 //*TIMING
2183  int count = 0;
2184  int totalCount;
2185  double mconc;
2186  double eterm;
2187  double oterm;
2188 
2189  edata = ( ck_use_fed->isChecked() && have_model() )
2190  ? &dsimList[ row ] : &dataList[ row ];
2191  scanCount = edata->scanCount();
2192  valueCount = edata->pointCount();
2193  boundPct = ct_boundaryPercent->value() / 100.0;
2194  positPct = ct_boundaryPos ->value() / 100.0;
2195  baseline = calc_baseline();
2196 
2197  divsCount = qRound( ct_division->value() );
2198  totalCount = scanCount * divsCount;
2199  divfac = 1.0 / (double)divsCount;
2200  correc = solution.s20w_correction * 1.0e13;
2201  omega = edata->scanData[ 0 ].rpm * M_PI / 30.0;
2202 
2203 
2204  // Iterate to adjust plateaus until none needed or max iters reached
2205 
2206  int iter = 1;
2207  int mxiter = 3; // maximum iterations
2208  double avdthr = 2.0e-5; // threshold cp-absavg-diff
2209 
2210  while( iter <= mxiter )
2211  {
2212  double avgdif = 0.0;
2213  count = 0;
2214 DbgLv(1) << "iter mxiter " << iter << mxiter;
2215 
2216  // Get division sedimentation coefficient values (intercepts)
2217 
2218  div_seds();
2219 
2220  int divsUsed = dseds.size();
2221 
2222  // Reset division plateaus
2223 
2224  for ( int ii = 0; ii < lscnCount; ii++ )
2225  {
2226  int js = liveScans[ ii ];
2227  dscan = &edata->scanData[ js ];
2228  omega = dscan->rpm * M_PI / 30.0;
2229  oterm = ( dscan->seconds - time_correction ) * omega * omega;
2230  eterm = -2.0 * oterm / correc;
2231  c0term = ( C0 - baseline ) * boundPct * divfac;
2232  span = ( scPlats[ ii ] - baseline ) * boundPct;
2233  sumcpij = 0.0;
2234 
2235  // Split the difference between divisions
2236 
2237  for ( int jj = 0; jj < divsCount; jj++ )
2238  { // Recalculate partial concentrations based on sedcoeff intercepts
2239  sedc = ( jj < divsUsed ) ? dseds[ jj ] : SEDC_NOVAL;
2240  cpij = ( sedc != SEDC_NOVAL ) ?
2241  ( c0term * exp( sedc * eterm ) ) :
2242  CPijs[ ii ][ jj ];
2243  CPijs[ ii ][ jj ] = cpij;
2244  sumcpij += cpij;
2245 //DbgLv(1) << " div " << jj+1 << " tcdps cpij " << tcpds.at(jj) << cpij;
2246  }
2247 
2248  // Set to split span-sum difference over each division
2249  sdiff = ( span - sumcpij ) * divfac;
2250 
2251  for ( int jj = 0; jj < divsCount; jj++ )
2252  { // Spread the difference to each partial plateau concentration
2253  CPijs[ ii ][ jj ] += sdiff;
2254  }
2255 
2256  avgdif += qAbs( sdiff ); // Sum of difference magnitudes
2257  count++;
2258 DbgLv(1) << " iter scn " << iter << ii+1 << " sumcpij span "
2259  << sumcpij << span << " sdiff sumabsdif" << sdiff << avgdif;
2260  }
2261 
2262  // Insure we have mid-division concentrations for newest partials
2263  update_mid_concs();
2264 
2265  avgdif /= (double)count; // Average of difference magnitudes
2266 DbgLv(1) << " iter" << iter << " avg(abs(sdiff))" << avgdif;
2267 
2268  if ( avgdif < avdthr ) // If differences are small, we're done
2269  {
2270 DbgLv(1) << " +++ avgdif < avdthr (" << avgdif << avdthr << ") +++";
2271  break;
2272  }
2273 
2274  iter++;
2275  }
2276 
2277  for ( int ii = 0; ii < totalCount; ii++ )
2278  pty[ ii ] = SEDC_NOVAL;
2279 
2280  int kk = 0; // Index to sed. coeff. values
2281  int kl = 0; // Index/count of live scans
2282  valueCount = edata->pointCount();
2283 
2284  // Calculate the corrected sedimentation coefficients
2285 
2286  for ( int ii = 0; ii < lscnCount; ii++ )
2287  {
2288  int js = liveScans[ ii ];
2289  dscan = &edata->scanData[ js ];
2290 
2291  double timev = dscan->seconds - time_correction;
2292  double timex = 1.0 / sqrt( timev );
2293  double bdrad = bdrads.at( kl ); // Back-diffus cutoff radius for scan
2294  double bdcon = bdcons.at( kl++ ); // Back-diffus cutoff concentration
2295  double divrad = 0.0; // Division radius value
2296 DbgLv(1) << "scn liv" << ii+1 << kl
2297  << " radius concen time" << bdrad << bdcon << timev;
2298 
2299  ptx[ ii ] = timex; // Save corrected time and accum max
2300 
2301  omega = dscan->rpm * M_PI / 30.0;
2302  oterm = ( timev > 0.0 ) ? ( timev * omega * omega ) : -1.0;
2303 
2304  for ( int jj = 0; jj < divsCount; jj++ )
2305  { // walk through division points; get sed. coeff. by place in readings
2306  mconc = mconcs[ ii ][ jj ]; // Mid div concentration
2307 
2308  sedc = sed_coeff( mconc, oterm, &divrad, NULL );
2309 
2310  if ( divrad > bdrad )
2311  { // Mark a point to be excluded by back-diffusion
2312  sedc = SEDC_NOVAL;
2313 DbgLv(1) << " *excl* div" << jj+1 << " drad dcon " << divrad << mconc;
2314  }
2315 
2316  // Y value of point is sedcoeff
2317  pty[ kk++ ] = sedc;
2318  }
2319  }
2320 //*TIMING
2321 kmsecs[6]+=sttime.msecsTo(QDateTime::currentDateTime());
2322 //*TIMING
2323 
2324 }
2325 
2326 // Flag whether we have a model to use for finite-element plateau determination
2327 bool US_vHW_Enhanced::have_model()
2328 {
2329  row = lw_triples->currentRow();
2330  ti_noise = tinoises[ row ];
2331  ri_noise = rinoises[ row ];
2332  int n_ti_noi = ti_noise.values.size();
2333  int n_ri_noi = ri_noise.values.size();
2334  QString modelGUID;
2335 
2336  if ( n_ti_noi > 0 )
2337  modelGUID = ti_noise.modelGUID;
2338  else if ( n_ri_noi > 0 )
2339  modelGUID = ri_noise.modelGUID;
2340 
2341  return ( ! modelGUID.isEmpty() && modelGUID.length() == 36 );
2342 }
2343 
2344 // Plot any upper plot vertical excluded-scan lines
2345 void US_vHW_Enhanced::vert_exclude_lines()
2346 {
2347  double ymax;
2348  double xx[ 2 ];
2349  double yy[ 2 ];
2350  int frsc = qRound( ct_from->value() );
2351  int tosc = qRound( ct_to ->value() );
2352 
2353  if ( tosc <= 0 ) return;
2354 
2355  // First remove previously drawn red vertical lines
2356  QwtPlotItemList list = data_plot1->itemList();
2357  for ( int ii = 0; ii < list.size(); ii++ )
2358  {
2359  QwtPlotItem* curve = list[ ii ];
2360  if ( curve->title().text().contains( "Exclude Marker" ) )
2361  {
2362  curve->detach();
2363  }
2364  }
2365 
2366  ymax = 0.0;
2367  for ( int ii = 0; ii < aseds.size(); ii++ )
2368  ymax = qMax( ymax, aseds[ ii ] ); // Max sed.coeff. value
2369  frsc = ( frsc < 1 ) ? 1 : frsc;
2370  yy[ 0 ] = 0.5;
2371  yy[ 1 ] = ymax - 0.1;
2372 
2373  for ( int ii = 0; ii < lscnCount; ii++ )
2374  {
2375  int ixsc = ii + 1;
2376  if ( ixsc < frsc ) continue;
2377  if ( ixsc > tosc ) break;
2378 
2379  int js = liveScans[ ii ];
2380  xx[ 0 ] = 1.0 / sqrt( edata->scanData[ js ].seconds - time_correction );
2381  xx[ 1 ] = xx[ 0 ];
2382  curve = us_curve( data_plot1,
2383  tr( "Scan %1 Exclude Marker" ).arg( ixsc ) );
2384  curve->setPen( QPen( QBrush( Qt::red ), 1.0 ) );
2385  curve->setData( xx, yy, 2 );
2386  }
2387 }
2388 
2389 // Slot to handle exclude-from change
2390 void US_vHW_Enhanced::exclude_from( double from )
2391 {
2392  double to = ct_to->value();
2393 
2394  if ( to < from )
2395  { // Adjust exclude-to if need be
2396  ct_to->disconnect();
2397  ct_to->setValue( from );
2398 
2399  connect( ct_to, SIGNAL( valueChanged( double ) ),
2400  SLOT ( exclude_to ( double ) ) );
2401  }
2402 
2403  // Mark upper plot excluded scans, then lower plot ones
2404  skipPlot = true;
2405  vert_exclude_lines();
2406  data_plot1->replot();
2407 
2408  plot_data2();
2409  skipPlot = false;
2410 }
2411 
2412 // Slot to handle exclude-to change
2413 void US_vHW_Enhanced::exclude_to( double to )
2414 {
2415 DbgLv(1) << "(1)TO=" << to;
2416  double from = ct_from->value();
2417 
2418  if ( from > to || from == 0 )
2419  { // Adjust exclude-from if need be
2420 DbgLv(1) << "(2)TO=" << to;
2421  ct_from->disconnect();
2422  if ( from > to )
2423  ct_from->setValue( to );
2424  else if ( to > 0.0 )
2425  ct_from->setValue( 1.0 );
2426 
2427  connect( ct_from, SIGNAL( valueChanged( double ) ),
2428  SLOT ( exclude_from( double ) ) );
2429 DbgLv(1) << "(3)TO=" << to;
2430  }
2431 
2432  // Mark upper plot excluded scans, then lower plot ones
2433  if ( to > 0.0 )
2434  {
2435 DbgLv(1) << "(4)TO=" << to;
2436  skipPlot = true;
2437  vert_exclude_lines();
2438  data_plot1->replot();
2439 
2440  plot_data2();
2441  skipPlot = false;
2442 DbgLv(1) << "(5)TO=" << to;
2443  }
2444 
2445  else // Special case of to=0 after exclude clicked
2446  {
2447 DbgLv(1) << "(8)TO=" << to;
2448  data_plot();
2449 DbgLv(1) << "(9)TO=" << to;
2450  }
2451 }
2452 
2453 // Build vectors of live scan information
2454 void US_vHW_Enhanced::live_scans()
2455 {
2456 //*TIMING
2457 kcalls[8]+=1;QDateTime sttime=QDateTime::currentDateTime();
2458 //*TIMING
2459  row = lw_triples->currentRow();
2460  edata = &dataList[ row ];
2461  scanCount = edata->scanCount();
2462  lscnCount = 0;
2463  liveScans.clear();
2464 
2465  for ( int ii = 0; ii < scanCount; ii++ )
2466  { // For non-excluded scans, save scan index and plateau concentration
2467  if ( ! excludedScans.contains( ii ) )
2468  {
2469  liveScans << ii; // Save original scan index
2470  lscnCount++; // Bump count of live scans
2471  }
2472  }
2473 //*TIMING
2474 kmsecs[8]+=sttime.msecsTo(QDateTime::currentDateTime());
2475 //*TIMING
2476 }
2477 
2478 // Build initial versions of partial-conc and mid-div conc vectors
2479 void US_vHW_Enhanced::init_partials()
2480 {
2481 //*TIMING
2482 kcalls[9]+=1;QDateTime sttime=QDateTime::currentDateTime();
2483 //*TIMING
2484  double mconc; // Mid-division concentration
2485  double cinc; // Constant division concentration increment
2486  divsCount = qRound( ct_division->value() );
2487  divfac = 1.0 / (double)divsCount;
2488  CPijs .clear();
2489  mconcs.clear();
2490 
2491  for ( int ii = 0; ii < lscnCount; ii++ )
2492  { // Populate partial concentration lists
2493  QVector< double > div_pconcs;
2494  QVector< double > div_mconcs;
2495  range = scPlats[ ii ] - baseline; // Range for this scan
2496  cinc = range * boundPct * divfac; // Conc. incr. this scan
2497  mconc = baseline + range * positPct + cinc * 0.5; // 1st Mid-div conc
2498  div_pconcs.clear();
2499  div_mconcs.clear();
2500 
2501  for ( int jj = 0; jj < divsCount; jj++ )
2502  {
2503  div_pconcs << cinc; // Save partial concentration
2504  div_mconcs << mconc; // Save mid-div concentration
2505  mconc += cinc; // Bump to next division
2506  }
2507 
2508  CPijs << div_pconcs; // Save a vector for each scan
2509  mconcs << div_mconcs;
2510  }
2511 //*TIMING
2512 kmsecs[9]+=sttime.msecsTo(QDateTime::currentDateTime());
2513 //*TIMING
2514 }
2515 
2516 // Update mid-division concentrations
2517 void US_vHW_Enhanced::update_mid_concs()
2518 {
2519 //*TIMING
2520 kcalls[10]+=1;QDateTime sttime=QDateTime::currentDateTime();
2521 //*TIMING
2522  double bconc; // Base concentration of a division
2523  double mconc; // Mid concentration of a division
2524  double cpij; // Partial concentration of a division
2525 
2526  for ( int ii = 0; ii < lscnCount; ii++ )
2527  { // Populate mid-division concentration vector of each scan
2528  range = scPlats[ ii ] - baseline;
2529  bconc = baseline + range * positPct; // Division 1 base conc.
2530 
2531  for ( int jj = 0; jj < divsCount; jj++ )
2532  { // Determine new concentrations in divisions of the scan
2533  cpij = CPijs[ ii ][ jj ]; // Partial concentration of div.
2534  mconc = bconc + cpij * 0.5; // Mid-division concentration
2535  bconc += cpij; // Next division base
2536  mconcs[ ii ][ jj ] = mconc; // Save new mid-div concentration
2537 //DbgLv(1) << "UPD: mconc cpij" << mconc << cpij << "ii jj" << ii << jj;
2538  }
2539  }
2540 //*TIMING
2541 kmsecs[10]+=sttime.msecsTo(QDateTime::currentDateTime());
2542 //*TIMING
2543 }
2544 
2545 // Create simulation for a data set (if possible)
2546 void US_vHW_Enhanced::create_simulation()
2547 {
2548  row = lw_triples->currentRow();
2549  expda = &dataList[ row ];
2550  simda = &dsimList[ row ];
2551  edata = simda;
2552  US_DataIO::RawData simraw;
2553  US_SimulationParameters simparams;
2554  US_Passwd pw;
2555  US_DB2* dbP = disk_controls->db() ? new US_DB2( pw.getPasswd() ) : 0;
2556 DbgLv(1) << " cr.sim.: row dbP" << row << dbP;
2557 
2558  // Initialize simulation parameters and simulation data
2559  simparams.initFromData( dbP, *expda, true );
2560 simparams.debug();
2561  US_AstfemMath::initSimData( simraw, *expda, 0.0 );
2562 
2563  // Insure we have a model for the current triple
2564  get_model();
2565 DbgLv(1) << " cr.sim.: model:" << model.description;
2566 
2567  US_Model amodel = model;
2568  solution.density = le_density ->text().toDouble();
2569  solution.viscosity = le_viscosity->text().toDouble();
2570  double avgTemp = expda->average_temperature();
2571  bool cnst_vbar = model.constant_vbar();
2572  bool vari_vbar = ! cnst_vbar;
2573 
2574  if ( cnst_vbar )
2575  { // Compute data corrections for constant vbar
2576  solution.vbar20 = le_vbar ->text().toDouble();
2577  solution.vbar = US_Math2::calcCommonVbar( solution_rec, avgTemp );
2578  US_Math2::data_correction( avgTemp, solution );
2579  }
2580 
2581  // Convert the model to experiment space
2582 
2583  for ( int ii = 0; ii < amodel.components.size(); ii++ )
2584  {
2585  if ( vari_vbar )
2586  { // Compute data corrections for varying vbar
2587  solution.vbar20 = amodel.components[ ii ].vbar20;
2588  solution.vbar = US_Math2::adjust_vbar( solution.vbar20, avgTemp );
2589  US_Math2::data_correction( avgTemp, solution );
2590  }
2591 
2592  amodel.components[ ii ].s /= solution.s20w_correction;
2593  amodel.components[ ii ].D /= solution.D20w_correction;
2594  }
2595 
2596  // Compute simulation data for the current triple and its model
2597 DbgLv(1) << " cr.sim.: new astfem";
2598  US_Astfem_RSA* astfem_rsa = new US_Astfem_RSA( amodel, simparams );
2599 DbgLv(1) << " cr.sim.: astfem calc";
2600  astfem_rsa->calculate( simraw );
2601 DbgLv(1) << " cr.sim.: astfem calc RTN";
2602 int ll=edata->scanCount()-1;
2603 int nn=edata->pointCount()-1;
2604 int mm=nn/2;
2605 DbgLv(1) << " cr.sim.: sd 00 0m 0n l0 lm ln"
2606  << simraw.value( 0, 0)
2607  << simraw.value( 0,mm)
2608  << simraw.value( 0,nn)
2609  << simraw.value(ll, 0)
2610  << simraw.value(ll,mm)
2611  << simraw.value(ll,nn);
2612 
2613  // Set values in the editedData simulation
2614  for ( int ii = 0; ii < simda->scanCount(); ii++ )
2615  for ( int jj = 0; jj < simda->pointCount(); jj++ )
2616  simda->setValue( ii, jj, simraw.value( ii, jj ) );
2617 
2618  // Flag that the current triple has a simulation
2619  have_sims[ row ] = true;
2620 }
2621 
2622 // Show lower data plot depending on "Use FE Data" setting
2623 void US_vHW_Enhanced::plot_data2()
2624 {
2625  if ( ! ck_use_fed->isChecked() || ! have_model() )
2626  { // Normal experimental data plot
2627  US_AnalysisBase2::data_plot();
2628  return;
2629  }
2630 
2631  // Plot of simulation data when Use FE Data is checked
2632  int row = lw_triples->currentRow();
2633  simda = &dsimList[ row ];
2634 
2635  QString dataType = tr( "Absorbance" );
2636  if ( simda->dataType == "RI" ) dataType = tr( "Intensity" );
2637  if ( simda->dataType == "WI" ) dataType = tr( "Intensity" );
2638  if ( simda->dataType == "IP" ) dataType = tr( "Interference" );
2639  if ( simda->dataType == "FI" ) dataType = tr( "Fluorescence" );
2640 
2641  QString header = tr( "Simulation Velocity Data for\n ") + simda->runID
2642  + " (" + simda->cell + "/" + simda->channel
2643  + "/" + simda->wavelength + ")";
2644  data_plot2->setTitle( header );
2645 
2646  header = dataType + tr( " at " ) + simda->wavelength + tr( " nm" );
2647  data_plot2->setAxisTitle( QwtPlot::yLeft, header );
2648 
2649  header = tr( "Radius (cm) " );
2650  data_plot2->setAxisTitle( QwtPlot::xBottom, header );
2651 
2652  data_plot2->clear();
2653  us_grid( data_plot2 );
2654 
2655  int scan_number = 0;
2656  int from = (int)ct_from->value();
2657  int to = (int)ct_to ->value();
2658 
2659  int scanCount = simda->scanCount();
2660  int points = simda->pointCount();
2661  double boundaryPct = ct_boundaryPercent->value() / 100.0;
2662  boundaryPct = ct_boundaryPercent->isEnabled() ? boundaryPct : 9.0;
2663  double positionPct = ct_boundaryPos ->value() / 100.0;
2664  double baseline = calc_baseline();
2665 
2666  QVector< double > rvec( points );
2667  QVector< double > vvec( points );
2668  double* rr = rvec.data();
2669  double* vv = vvec.data();
2670 
2671  // Calculate basic parameters for other functions
2672  time_correction = US_Math2::time_correction( dsimList );
2673 
2674  solution.density = le_density ->text().toDouble();
2675  solution.viscosity = le_viscosity->text().toDouble();
2676  solution.vbar20 = le_vbar ->text().toDouble();
2677  solution.manual = manual;
2678  double avgTemp = simda->average_temperature();
2679  solution.vbar = US_Math2::calcCommonVbar( solution_rec, avgTemp );
2680 
2681  US_Math2::data_correction( avgTemp, solution );
2682 
2683  // Draw curves
2684  for ( int ii = 0; ii < scanCount; ii++ )
2685  {
2686  if ( excludedScans.contains( ii ) ) continue;
2687 
2688  scan_number++;
2689  bool highlight = ( scan_number >= from && scan_number <= to );
2690  simsc = &simda->scanData[ ii ];
2691  double range = simsc->plateau - baseline;
2692  double lower_limit = baseline + range * positionPct;
2693  double upper_limit = lower_limit + range * boundaryPct;
2694  QString curvname = tr( "Curve %1 " ).arg( ii );
2695 
2696  int jj = 0;
2697  int count = 0;
2698 
2699  // Plot each scan in (up to) three segments: below, in, and above
2700  // the specified boundaries
2701  while ( jj < points && simsc->rvalues[ jj ] < lower_limit )
2702  {
2703  rr[ count ] = simda->xvalues[ jj ];
2704  vv[ count ] = simsc->rvalues[ jj ];
2705  jj++;
2706  count++;
2707  }
2708 
2709  QString title;
2710  QwtPlotCurve* cc;
2711 
2712  if ( count > 1 )
2713  {
2714  title = curvname + tr( "below range" );
2715  cc = us_curve( data_plot2, title );
2716 
2717  if ( highlight )
2718  cc->setPen( QPen( Qt::red ) );
2719  else
2720  cc->setPen( QPen( Qt::cyan ) );
2721 
2722  cc->setData( rr, vv, count );
2723  }
2724 
2725  count = 0;
2726 
2727  while ( jj < points && simsc->rvalues[ jj ] < upper_limit )
2728  {
2729  rr[ count ] = simda->xvalues[ jj ];
2730  vv[ count ] = simsc->rvalues[ jj ];
2731  jj++;
2732  count++;
2733  }
2734 
2735  if ( count > 1 )
2736  {
2737  title = curvname + tr( "in range" );
2738  cc = us_curve( data_plot2, title );
2739 
2740  if ( highlight )
2741  cc->setPen( QPen( Qt::red ) );
2742  else
2743  cc->setPen( QPen( US_GuiSettings::plotCurve() ) );
2744 
2745  cc->setData( rr, vv, count );
2746  }
2747 
2748  count = 0;
2749 
2750  while ( jj < points )
2751  {
2752  rr[ count ] = simda->xvalues[ jj ];
2753  vv[ count ] = simsc->rvalues[ jj ];
2754  jj++;
2755  count++;
2756  }
2757 
2758  if ( count > 1 )
2759  {
2760  title = curvname + tr( "above range" );
2761  cc = us_curve( data_plot2, title );
2762 
2763  if ( highlight )
2764  cc->setPen( QPen( Qt::red ) );
2765  else
2766  cc->setPen( QPen( Qt::cyan ) );
2767 
2768  cc->setData( rr, vv, count );
2769  }
2770  }
2771 
2772  data_plot2->replot();
2773 
2774  return;
2775 }
2776 
2777 // Insure we get the model associated with the current noises
2778 void US_vHW_Enhanced::get_model()
2779 {
2780  row = lw_triples->currentRow();
2781  int n_ti_noi = ti_noise.values.size();
2782  int n_ri_noi = ri_noise.values.size();
2783  model = modlList[ row ];
2784  QString mmodlGUID = model.modelGUID;
2785  QString modelGUID;
2786 
2787  modelGUID = ( n_ti_noi > 0 ) ? ti_noise.modelGUID : modelGUID;
2788  modelGUID = ( n_ri_noi > 0 ) ? ri_noise.modelGUID : modelGUID;
2789 
2790 DbgLv(1) << "gMo: modelGUID" << modelGUID;
2791  if ( modelGUID.isEmpty() || modelGUID.length() != 36 ||
2792  ( ! mmodlGUID.isEmpty() && mmodlGUID == modelGUID ) )
2793  { // No new model load if none exists or already loaded matches noise
2794 DbgLv(1) << "gMo: MODEL empty/matches IDln" << modelGUID.length()
2795  << "modID" << modelGUID << "mmoID" << mmodlGUID;
2796  return;
2797  }
2798 
2799  US_Passwd pw;
2800  bool isDB = disk_controls->db();
2801  US_DB2* db = isDB ? new US_DB2( pw.getPasswd() ) : 0;
2802 
2803  model.load( isDB, modelGUID, db );
2804 
2805  mmodlGUID = model.modelGUID;
2806  modlList[ row ] = model;
2807 DbgLv(1) << "gMo: loaded model: modelGUID" << modelGUID << mmodlGUID;
2808 }
2809 
Generated on Sun Jun 21 2015 11:40:43 for UltraScan III by   doxygen 1.8.3.1-20130324