DataAnalysis.cpp

Go to the documentation of this file.

//DataAnalysis.cpp

/*************************************************************/
//
// GreyLabWindow
//
// Data analysis routines for the GreyLabWindow class
//
//-------------------------------------------------------------
//
// v 1.20 - 2008/03/10 - Matlab export
// v 1.19 - 2007/12/13 - Rounding of z-axis, fixed
// v 1.18 - 2007/12/13 - Rounding of z-axis
// v 1.17 - 2007/11/26 - Fixed auto binning when using functions
// v 1.16 - 2007/06/15 - Updates to cursor labels
// v 1.15 - 2007/06/12 - 2.17.2 improvements
// v 1.14 - 2007/06/06 - Fixed lack of image binding when writing output
// v 1.13 - 2007/05/30 - Zoom back feature
// v 1.12 - 2007/02/17 - Changes for linescan range plotting
// v 1.11 - 2007/02/09 - X target range and better autoranging
// v 1.10 - 2007/02/01 - ysd scatter option
// v 1.9 - 2007/01/06 - Changes for 2.11.0
// v 1.8 - 2006/11/23 - Hires image save etc
// v 1.7 - 2006/11/19 - Autorange Z and equal range Z
// v 1.6 - 2006/11/19 - Fixed zooming
// v 1.5 - 2006/11/18 - Lots of changes for 2.10.0
// v 1.5 - 2006/11/13 - Export all columns options
// v 1.4 - 2006/11/10 - Fixed mouse clicks
// v 1.3 - 2006/09/07 - Acorn export of data
// v 1.2 - 2006/09/06 - Tab on export was the wrong side of the data
// v 1.1 - 2006/08/29 - Fixed a bug crashing the processing when 1st sweep > 0 in a block of sweeps
// v 1.0 - 2006/08/28 - Initial release
//
// See LICENSE.txt for distribution and usage restrictions
// Copyright (c) 2005-2007 Simon Chorley
// www.mylaboratory.co.uk - greylab@mylaboratory.co.uk
//
/*************************************************************/


///Data analysis functions
/**
 * Provides implementation of GreyLabWindow functions to process the data.
 * @file DataAnalysis.cpp
 * @version 1.20 - 2008/03/10
 */


#include "GreyLab.h"


#include <fstream>


FXString strip(FXString s)
{
        unsigned int i;
        for (i = 0; i < s.length(); ++i)
                if (!(((s[i] >= 'a') && (s[i] <= 'z')) || ((s[i] >= 'A') && (s[i] <= 'Z')) || ((s[i] >= '0') && (s[i] <= '9')) || (s[i] == '_')))
                        s[i] = '_';
        return(s);
}


///Rounds a number to the nearest 1, 2 or 5
/**
 * @param step Number to quantise.
 * @return Rounded number.
 */
double GreyLabWindow::quantise(double step)
{
        int i;
        double test = fabs(step)*0.1;
        for (i = 0; i < 500; ++i)
        {
                test *= 10.0;
                if ((test > 0.75) && (test <= 1.5))     //1
                        return(pow(10.0, -i));
                if ((test > 1.5) && (test <= 3.5))      //2
                        return(2.0*pow(10.0, -i));
                if ((test > 3.5) && (test <= 7.5))      //5
                        return(5.0*pow(10, -i));
        }
        test = fabs(step)*10.0;
        for (i = 0; i < 500; ++i)
        {
                test *= 0.1;
                if ((test > 0.75) && (test <= 1.5))     //1
                        return(pow(10.0, i));
                if ((test > 1.5) && (test <= 3.5))      //2
                        return(2.0*pow(10.0, i));
                if ((test > 3.5) && (test <= 7.5))      //5
                        return(5.0*pow(10, i));
        }
        return(1.0);
}


double GreyLabWindow::quantiseupZ(double step)
{
        int i, j;
        for (j = 0; j < 10; ++j)
                if ((fabs(step) > (float)j) && (fabs(step) <= (float)(j+1)))
                        return((float)(j+1));
/*      double test = fabs(step)*0.1;
        for (i = 0; i < 500; ++i)
        {
                test *= 10.0;
                for (j = 1; j < 10; ++j)
                        if ((test > (float)j) && (test <= (float)(j+1)))
                                return((float)(j+1)*pow(10.0, -i));
        }
        test = fabs(step)*10.0;
        for (i = 0; i < 500; ++i)
        {
                test *= 0.1;
                for (j = 1; j < 10; ++j)
                        if ((test > (float)j) && (test <= (float)(j+1)))
                                return((float)(j+1)*pow(10.0, -i));
        }*/
        return(1.0);
}


double GreyLabWindow::quantisedownZ(double step)
{
        int i, j;
        for (j = 0; j < 10; ++j)
                if ((fabs(step) > (float)j) && (fabs(step) <= (float)(j+1)))
                        return((float)(j));
/*      double test = fabs(step)*0.1;
        for (i = 0; i < 500; ++i)
        {
                test *= 10.0;
                for (j = 1; j < 10; ++j)
                        if ((test > (float)j) && (test <= (float)(j+1)))
                                return((float)(j)*pow(10.0, -i));
        }
        test = fabs(step)*10.0;
        for (i = 0; i < 500; ++i)
        {
                test *= 0.1;
                for (j = 1; j < 10; ++j)
                        if ((test > (float)j) && (test <= (float)(j+1)))
                                return((float)(j)*pow(10.0, -i));
        }*/
        return(1.0);
}


void GreyLabWindow::roundZ(double &zlo, double &zhi)
{
        int i = 0;
        double range = zhi-zlo;
        if (range == 0.0)
                return;
        double multiplier = 100/range;
        if (multiplier > 1.0)
        {
                while (multiplier > 1.0)
                {
                        multiplier *= 0.1;
                        ++i;
                }
                multiplier = pow(10, i);
//std::cout << zlo << " -> ";
                zlo *= multiplier;
                if (zlo < 0.0)
                        zlo = (((int)zlo)+(zlo/fabs(zlo))*quantiseupZ((zlo-(int)zlo)))/multiplier;
                else
                        zlo = (((int)zlo)+(zlo/fabs(zlo))*quantisedownZ((zlo-(int)zlo)))/multiplier;
//std::cout << zlo << "\n";
//std::cout << zhi << " -> ";
                zhi *= multiplier;
                if (zhi > 0.0)
                        zhi = (((int)zhi)+(zhi/fabs(zhi))*quantiseupZ((zhi-(int)zhi)))/multiplier;
                else
                        zhi = (((int)zhi)+(zhi/fabs(zhi))*quantisedownZ((zhi-(int)zhi)))/multiplier;
//std::cout << zhi << "\n";
        }
        else if (multiplier <= 1.0)
        {
                while (multiplier < 1.0)
                {
                        multiplier *= 10.0;
                        ++i;
                }
                multiplier = pow(10, 1-i);
//std::cout << zlo << " -> ";
                zlo *= multiplier;
                if (zlo < 0.0)
                        zlo = (((int)zlo)+(zlo/fabs(zlo))*quantiseupZ((zlo-(int)zlo)))/multiplier;
                else
                        zlo = (((int)zlo)+(zlo/fabs(zlo))*quantisedownZ((zlo-(int)zlo)))/multiplier;
//std::cout << zlo << "\n";
//std::cout << zhi << " -> ";
                zhi *= multiplier;
                if (zhi > 0.0)
                        zhi = (((int)zhi)+(zhi/fabs(zhi))*quantiseupZ((zhi-(int)zhi)))/multiplier;
                else
                        zhi = (((int)zhi)+(zhi/fabs(zhi))*quantisedownZ((zhi-(int)zhi)))/multiplier;
//std::cout << zhi << "\n";
        }
}


///Processes the raw data into histogram data.
/**
 * Processes the raw data into histogram data. All paramaters are unused.
 * @remarks Sets GreyLabWindow::pdata#validdata to true on success.
 */
long GreyLabWindow::cmdDAProcess(FXObject*, FXSelector, void*)
{
        if (rawdata == NULL)
                return(1);

        if (gsDCwindow == NULL)
        {
                gsDCwindow = new FXDCWindow((FXDrawable*)datawin->dataImageview);
                gsDCwindow->end();
        }
        if (lsDCwindow == NULL)
        {
                lsDCwindow = new FXDCWindow((FXDrawable*)linewin->lineImageview);
                lsDCwindow->end();
        }

        lstracesIconlist->clearItems();
        lsinfo.newTrace = true;

        //Check functions
        if (onDAParseError(functions[0].InfixToRPN((char*)dafunctionText[0]->getText().text()), 0) == false) return(1);
        if (onDAParseError(functions[1].InfixToRPN((char*)dafunctionText[1]->getText().text()), 1) == false) return(1);
        if (onDAParseError(functions[2].InfixToRPN((char*)dafunctionText[2]->getText().text()), 2) == false) return(1);

        pdata.validdata = false;

        //Scatter?
        bool scatter = dascatterCheckbutton->getCheck();
        
        unsigned long i, j, k;
        //Clear any old data
        if (pdata.sweep != NULL) delete[] pdata.sweep;
        if (pdata.x != NULL) delete[] pdata.x;
        if (pdata.y != NULL) delete[] pdata.y;
        if (pdata.z != NULL)
        {
                for (i = 0; i < pdata.ny; ++i)  //delete z data
                {
                        if (pdata.z[i] != NULL) delete[] pdata.z[i];
                        pdata.z[i] = NULL;
                }
                delete[] pdata.z;               
        }
        if (pdata.scattery != NULL)
        {
                for (i = 0; i < pdata.ny; ++i)  //delete z data
                {
                        if (pdata.scattery[i] != NULL) delete[] pdata.scattery[i];
                        pdata.scattery[i] = NULL;
                }
                delete[] pdata.scattery;
                pdata.scattery = NULL;
        }
        if (pdata.columns != NULL) delete[] pdata.columns;
        pdata.x = NULL;
        pdata.y = NULL;
        pdata.z = NULL;
        pdata.columns = NULL;
        
        //Finalise SD
        sdinfo.stride = FXIntVal(dastrideTextfield->getText());
        if (sdinfo.stride < 1) sdinfo.stride = 1;
        sdinfo.xcol = daxcolListbox->getCurrentItem();
        sdinfo.ycol = daysdListbox->getCurrentItem()-2; //-2: col N, -1: manual, 0.... SD
        sdinfo.zcol = dazcolListbox->getCurrentItem();
        sdinfo.yscattercol = dascatterListbox->getCurrentItem()-1;
        sdinfo.ssweep = FXIntVal(dassweepTextfield->getText());
        sdinfo.esweep = FXIntVal(daesweepTextfield->getText());
        if (sdinfo.ssweep < 0)
        {
                sdinfo.ssweep = 0;
                dassweepTextfield->setText("0");
        }
        if (sdinfo.esweep >= sdinfo.nsweeps)
        {
                sdinfo.esweep = sdinfo.nsweeps-1;
                daesweepTextfield->setText(FXStringVal(sdinfo.esweep));
        }

        if (dasweepselectCheckbutton->getCheck() == false)
        {
//              sdinfo.nselected = (sdinfo.esweep-sdinfo.ssweep+1)/sdinfo.stride;               <-- This doesn't work properly when there are a non-stride multiple of sweeps that may still be valid (in terms of compatibility with the loop below) so replaced with a loop counter.
                sdinfo.nselected = 0;
                for (i = sdinfo.ssweep; i <= sdinfo.esweep; i += sdinfo.stride)
                        ++sdinfo.nselected;
                pdata.sweep = new int[sdinfo.nselected];
                for (i = sdinfo.ssweep; i <= sdinfo.esweep; i += sdinfo.stride)
                        pdata.sweep[(i-sdinfo.ssweep)/sdinfo.stride] = i;
                        
        }
        else
        {
                sdinfo.nselected = 0;
                for (i = sdinfo.ssweep; i <= sdinfo.esweep; i += sdinfo.stride)
                {
                        if (sdIconlist->isItemSelected(i) == true) ++sdinfo.nselected;
                }
                pdata.sweep = new int[sdinfo.nselected];
                int c = 0;
                for (i = sdinfo.ssweep; i <= sdinfo.esweep; i += sdinfo.stride)
                {
                        if (sdIconlist->isItemSelected(i) == true)
                        {
                                pdata.sweep[c] = i;
                                ++c;
                        }
                }
        }

        double val;
        pdata.ny = sdinfo.nselected;
        double lo, hi;
        if (daxautoCheckbutton->getCheck() == true)
        {
                lo = functions[0].Evaluate(FXDoubleVal(daxinminTextfield->getText()), FXDoubleVal(dayminTextfield->getText()), 0.0, 0, NULL);
                hi = functions[0].Evaluate(FXDoubleVal(daxinmaxTextfield->getText()), FXDoubleVal(dayminTextfield->getText()), 0.0, 0, NULL);
                daxtgtminTextfield->setText(FXStringVal(lo));
                daxtgtmaxTextfield->setText(FXStringVal(hi));
        }
        else
        {
                lo = FXDoubleVal(daxtgtminTextfield->getText());
                hi = FXDoubleVal(daxtgtmaxTextfield->getText());
        }
        if (hi < lo)    //Swap?
        {
                val = lo;
                lo = hi;
                hi = val;
        }
        if ((setautobinCheckbutton->getCheck() == true) && (setautoxscaleCheckbutton->getCheck() == true))      //added to redo bins when using functions - SJC 2007/11/26
                daxtgtresTextfield->setText(FXStringVal((hi-lo)/((float)sd[sdinfo.ssweep].npoints/FXFloatVal(setbinsizeTextfield->getText()))));
        double delta = FXDoubleVal(daxtgtresTextfield->getText())/2.0;
        pdata.nx = 1+abs((int)((hi-lo)/(2.0*delta)));

//Allocate arrays
        pdata.x = new double[pdata.nx];
        pdata.y = new double[pdata.ny];
        pdata.z = new point*[pdata.ny];
        if (scatter == true)
                pdata.scattery = new double*[pdata.ny];
        pdata.columns = new double[pdata.nx*pdata.ny*sdinfo.ncols];

//Generate x data
        for (i = 0; i < pdata.nx; ++i)
        {
                pdata.x[i] = lo+((float)i)*(2.0*delta);
//              std::cerr << pdata.x[i] << ',';
        }
        
//      cerr << "ysd" << sdinfo.ycol << '\n';
        if (sdinfo.ycol == -2)                          //generate y data
        {
                for (i = 0; i < pdata.ny; ++i)
                {
                        pdata.y[i] = pdata.sweep[i];
                }
        }
        else if (sdinfo.ycol == -1)
        {
                double s = FXDoubleVal(dayminTextfield->getText());
                double dy = (FXDoubleVal(daymaxTextfield->getText())-s)/((double)pdata.ny);
                for (i = 0; i < pdata.ny; ++i)
                {
                        pdata.y[i] = s+((double)i)*dy;
                }
        }
        else
                for (i = 0; i < pdata.ny; ++i)
                {
                        pdata.y[i] = sd[pdata.sweep[i]].item[sdinfo.ycol];
//                      pdata.y[i] = sd[i*sdinfo.stride+sdinfo.ssweep].item[sdinfo.ycol];
//                      cerr << "sd" << sd[i*sdinfo.stride+sdinfo.ssweep].item[sdinfo.ycol] << 'y' << pdata.y[i] << '\n';
                }

        unsigned long bin;
        unsigned long pt;
        double *row = new double[sdinfo.ncols];
        FXApp *a = getApp();

        statusProgressbar->setTotal(pdata.ny);
        for (i = 0; i < pdata.ny; ++i)  //transform and bin z data
        {
                double *tempyscatter;
                pdata.z[i] = new point[pdata.nx];       //allocate
                for (j = 0; j < pdata.nx; ++j)                          //zero data
                {
                        pdata.z[i][j].data = 0.0;
                        pdata.z[i][j].n = 0;
                }
                if (scatter == true)
                {
                        tempyscatter = new double[sd[pdata.sweep[i]].npoints];
                        pdata.scattery[i] = new double[pdata.nx];       //allocate
                        for (j = 0; j < pdata.nx; ++j)                          //zero data
                                pdata.scattery[i][j] = 0.0;
                        if (sdinfo.yscattercol == -1)
                        {
                                for (j = 0; j < sd[pdata.sweep[i]].npoints; ++j)
                                        tempyscatter[j] = pdata.y[i];   //generate from a ysd
                        }
                        else
                        {
                                for (j = 0; j < sd[pdata.sweep[i]].npoints; ++j)
                                        tempyscatter[j] = rawdata[sdinfo.yscattercol+sdinfo.ncols*pdata.sweep[i]][j];   //from a real column
                        }
                }
                for (j = 0; j < pdata.nx; ++j)          //zero array that stores other columns
                        for (k = 0; k < sdinfo.ncols; ++k)
                                pdata.columns[i*pdata.nx*sdinfo.ncols+j*sdinfo.ncols+k] = 0.0;
                for (j = 0; j < sd[pdata.sweep[i]].npoints; ++j)        //for each point in sweep
                {
                        pt = sdinfo.ncols*pdata.sweep[i];       //current column
//changed in 2.12.0 upgrade
/*                      val = rawdata[sdinfo.xcol+pt][j];                                               //x value
                        bin = (unsigned long)(0.5+(val-lo)/(2.0*delta));        //calculate bin
                        if ((bin >= 0) && (bin < pdata.nx))
                        {
                                for (k = 0; k < sdinfo.ncols; ++k)                              //collect data from row for evaluation
                                {
                                        row[k] = rawdata[k+pt][j];
                                }
                                //Evaluate z point and store
                                pdata.z[i][bin].data += functions[2].Evaluate(val, pdata.y[i], rawdata[sdinfo.zcol+pt][j], sdinfo.ncols, row);
                                if (scatter == true)
                                        pdata.scattery[i][bin] += functions[1].Evaluate(val, rawdata[sdinfo.yscattercol+pt][j], rawdata[sdinfo.zcol+pt][j], sdinfo.ncols, row);
                                for (k = 0; k < sdinfo.ncols; ++k)
                                        pdata.columns[i*pdata.nx*sdinfo.ncols+bin*sdinfo.ncols+k] += rawdata[k+pt][j];
                                ++pdata.z[i][bin].n;
                        }*/
                        for (k = 0; k < sdinfo.ncols; ++k)                              //collect data from row for evaluation
                        {
                                row[k] = rawdata[k+pt][j];
                        }
                        //calculate transformed x of point
                        if (scatter == false)
                                val = functions[0].Evaluate(rawdata[sdinfo.xcol+pt][j], pdata.y[i], rawdata[sdinfo.zcol+pt][j], sdinfo.ncols, row);
                        else
                                val = functions[0].Evaluate(rawdata[sdinfo.xcol+pt][j], tempyscatter[j], rawdata[sdinfo.zcol+pt][j], sdinfo.ncols, row);
                        if (j == 0)                     //if we are not scattering then just do it once
                        {
                                pdata.y[i] = functions[1].Evaluate(val, pdata.y[i], rawdata[sdinfo.zcol+pt][j], sdinfo.ncols, row);
                        }
                        bin = (unsigned long)(0.5+(val-lo)/(2.0*delta));        //calculate bin
                        if ((bin >= 0) && (bin < pdata.nx))
                        {
                                if (scatter == true)
                                {
                                        //transform y
                                        tempyscatter[j] = functions[1].Evaluate(val, tempyscatter[j], rawdata[sdinfo.zcol+pt][j], sdinfo.ncols, row);
                                        pdata.scattery[i][bin] += tempyscatter[j];
                                        //Evaluate z point and store
                                        pdata.z[i][bin].data += functions[2].Evaluate(val, tempyscatter[j], rawdata[sdinfo.zcol+pt][j], sdinfo.ncols, row);
                                }
                                else
                                {
                                        //Evaluate z point and store
                                        pdata.z[i][bin].data += functions[2].Evaluate(val, pdata.y[i], rawdata[sdinfo.zcol+pt][j], sdinfo.ncols, row);
                                }
                                for (k = 0; k < sdinfo.ncols; ++k)
                                        pdata.columns[i*pdata.nx*sdinfo.ncols+bin*sdinfo.ncols+k] += rawdata[k+pt][j];
                                ++pdata.z[i][bin].n;
                        }
                }
                for (j = 0; j < pdata.nx; ++j)                          //average bins
                {
                        pdata.z[i][j].data /= pdata.z[i][j].n;
                        for (k = 0; k < sdinfo.ncols; ++k)      //average array that stores other columns
                                pdata.columns[i*pdata.nx*sdinfo.ncols+j*sdinfo.ncols+k] /= pdata.z[i][j].n;
//                      cerr << 'z' << pdata.z[i][j].data;
                }
                if (scatter == true)
                {
                        for (j = 0; j < pdata.nx; ++j)                          //average bins
                                pdata.scattery[i][j] /= pdata.z[i][j].n;
                        delete[] tempyscatter;
                }

                statusProgressbar->setProgress(i+1);
                a->forceRefresh();
                a->repaint();
        }
        statusProgressbar->setProgress(0);
        delete[] row;

//interpolate
        int m, n;
        int maxinterp = FXIntVal(setinterpTextfield->getText());
        if (setinterpCheckbutton->getCheck() == true)
        {
                if (setinterpTogglebutton->getState() == 1)     //along X
                {
                        for (i = 0; i < pdata.ny; ++i)  //each sweep
                        {
                                for (j = 1; j < pdata.nx-1; ++j)        //each point in sweep
                                {
                                        if (pdata.z[i][j].n == 0)       //if a point is missing
                                        {
                                                if (pdata.z[i][j-1].n != 0)     //check that the one before is there
                                                {
                                                        for (k = j+1; (k < pdata.nx) && (k < j+maxinterp); ++k) //for all points in sweep after the missing one up to the max
                                                        {
                                                                if (pdata.z[i][k].n != 0)       //if a point is there
                                                                {
                                                                        for (m = 0; m < k-j; ++m)       //interpolate all points in the missing range
                                                                        {
                                                                                pdata.z[i][j+m].data = pdata.z[i][j-1].data+(pdata.z[i][k].data-pdata.z[i][j-1].data)*((float)(m+0.5)/(float)(k-j));
                                                                                for (n = 0; n < sdinfo.ncols; ++n)      //interpolate array that stores other columns
                                                                                        pdata.columns[i*pdata.nx*sdinfo.ncols+(j+m)*sdinfo.ncols+n] = pdata.columns[i*pdata.nx*sdinfo.ncols+(j-1)*sdinfo.ncols+n]+(pdata.columns[i*pdata.nx*sdinfo.ncols+k*sdinfo.ncols+n]-pdata.columns[i*pdata.nx*sdinfo.ncols+(j-1)*sdinfo.ncols+n])*((float)(m+0.5)/(float)(k-j));;
                                                                                if (scatter == true)
                                                                                        pdata.scattery[i][j+m] = pdata.scattery[i][j-1]+(pdata.scattery[i][k]-pdata.scattery[i][j-1])*((float)(m+0.5)/(float)(k-j));
                                                                                pdata.z[i][j+m].n = 1;
                                                                        }
                                                                        j += m;
                                                                        break;
                                                                }
                                                        }
                                                }
                                                else    //if one before is missing, no need to check next one as we can't fill the current one
                                                {
                                                        ++j;
                                                }
                                        }
                                }
                        }
                }
                else    //along Y
                {
                        for (j = 0; j < pdata.nx; ++j)  //each point in sweep
                        {
                                for (i = 1; i < pdata.ny-1; ++i)        //each sweep
                                {
                                        if (pdata.z[i][j].n == 0)       //if a point is missing
                                        {
                                                if (pdata.z[i-1][j].n != 0)     //check that the one before is there
                                                {
                                                        for (k = i+1; (k < pdata.ny) && (k < i+maxinterp); ++k) //for all points in sweep after the missing one up to the max
                                                        {
                                                                if (pdata.z[k][j].n != 0)       //if a point is there
                                                                {
                                                                        for (m = 0; m < k-i; ++m)       //interpolate all points in the missing range
                                                                        {
                                                                                pdata.z[i+m][j].data = pdata.z[i-1][j].data+(pdata.z[k][j].data-pdata.z[i-1][j].data)*((float)(m+0.5)/(float)(k-i));
                                                                                for (n = 0; n < sdinfo.ncols; ++n)      //interpolate array that stores other columns
                                                                                        pdata.columns[(i+m)*pdata.nx*sdinfo.ncols+j*sdinfo.ncols+n] = pdata.columns[(i-1)*pdata.nx*sdinfo.ncols+j*sdinfo.ncols+n]+(pdata.columns[k*pdata.nx*sdinfo.ncols+j*sdinfo.ncols+n]-pdata.columns[(i-1)*pdata.nx*sdinfo.ncols+j*sdinfo.ncols+n])*((float)(m+0.5)/(float)(k-j));;
                                                                                if (scatter == true)
                                                                                        pdata.scattery[i+m][j] = pdata.scattery[i-1][j]+(pdata.scattery[k][j]-pdata.scattery[i-1][j])*((float)(m+0.5)/(float)(k-i));
                                                                                pdata.z[i+m][j].n = 1;
                                                                        }
                                                                        i += m;
                                                                        break;
                                                                }
                                                        }
                                                }
                                                else    //if one before is missing, no need to check next one as we can't fill the current one
                                                {
                                                        ++i;
                                                }
                                        }
                                }
                        }
                }
        }

/* 01/02/2007 - moved transformation to before the Z happens and added ability to generate scatter data from non-scatter data for full x/y/z manipulation
        for (i = 0; i < pdata.nx; ++i)  //transform x data
        {
                pdata.x[i] = functions[0].Evaluate(lo+(float)i*2.0*delta, 0.0, 0.0, 0, NULL);
//              cerr << 'x' << pdata.x[i];
        }
        for (i = 0; i < pdata.ny; ++i)  //transform y data
                pdata.y[i] = functions[1].Evaluate(0.0, pdata.y[i], 0.0, 0, NULL);
*/


        int o;
        for (o = 0; o < 2; ++o)         //loop to allow the order to be swapped
        {
                //Differentiate
                if ((dadiffCheckbutton->getCheck() == true)
                        && (((daorderCheckbutton->getCheck() == false) && (o == 0)) || ((daorderCheckbutton->getCheck() == true) && (o == 1))))
                {
                        double *diffed;
//                      if (dadiffTogglebutton->getState() == true)
                        int di = dadiffOptionmenu->getCurrentNo();
                        if (di == 0)
                        {
                                diffed = new double[pdata.nx];
                                if (setaltdifferentialCheckbutton->getCheck() == true)
                                {       //2 point differential
                                        for (i = 0; i < pdata.ny; ++i)
                                        {
                                                for (j = 0; j < pdata.nx-1; ++j)
                                                {
                                                        diffed[j] = (pdata.z[i][j+1].data-pdata.z[i][j].data)/(pdata.x[j+1]-pdata.x[j]);
                                                }
                                                diffed[pdata.nx-1] = diffed[pdata.