Ctunes.cpp

Go to the documentation of this file.

/*****************************************************************************/
/*                                                                           */
/* Class TUNE                                                                */
/* ==============                                                            */
/*                                                                           */
/* ASM, September 2001                                                       */
/*****************************************************************************/
#include <algorithm>
#include <memory>
#include <vector>
#include <cstring>

#include "Utility/Inform.h"

#include "Algorithms/Ctunes.h"
#include "Algorithms/lomb.h"

extern Inform *gmsg;

//RANLIB_class rndm(265314159,4);


TUNE_class::TUNE_class():
        ofac(0.0),
        hifac(0.0),
        Qmin(0.0),
        Qmax(0.0)
/*---------------------------------------------------------------------------*
 * constructor
 * ===========
 *
 *---------------------------------------------------------------------------*/
{
}

TUNE_class::~TUNE_class(void)
/*---------------------------------------------------------------------------*
 * destructor
 * ==========
 *
 *---------------------------------------------------------------------------*/
{

    // Do nothing......

}

int TUNE_class::lombAnalysis(std::vector<double> &x, std::vector<double> &y, int /*nhis*/, double Norm)
/*-----------------------------------------------------------------------------
 *  Launch Lomb analysis and plot results
 *  =======================================
 *
 *---------------------------------------------------------------------------*/
{

    int Ndat = x.size();
    int    i, nout, jmax;
    int    pairc;
    int    datcnt = 0;
    int stat = 0;
    double prob, probi;
    double tofac = 0.8;

    LOMB_TYPE tlom;

    CI_lt p, q;

    std::vector<LOMB_TYPE> lodata, lodata2;
    /*---------------------------------------------------------------------------*/

    /*
     * Do Lomb analysis
     * ================
     */

    for(int j = 0; j < Ndat; j++) {
        tlom.x = x[j];
        tlom.y = y[j];
        lodata.push_back(tlom);
    }

    p = lodata.begin();
    q = lodata.end();

    datcnt = (int) count_if(p, q, Lomb_eq(0.));

    if(datcnt > (q - p - 10)) {
        *gmsg << "* Just found " << datcnt << " data points that are == 0!" << endl;
        return(-1);
    }

    // this parameterset works ok in most cases.....
    ofac  = 4.0;
    hifac = 0.8;
    Qmin  = 0.2;
    Qmax  = 0.4;

    std::unique_ptr<LOMB_class> la(new LOMB_class(1));

    stat = 0;
    stat = la->period(&lodata, &lodata2, ofac, hifac, &nout, &jmax, &prob, 0);
    if(stat != 0) {
        *gmsg << "* @C3ERROR: Lomb analysis failed!" << endl;
        return(-1);
    }

    double pairx[nout];
    double pairy[nout];

    pairc = 0;
    for(i = 0; i < nout; i++) {
        if(lodata2[i].y > 2.) {
            pairx[pairc] = lodata2[i].x;
            pairy[pairc] = lodata2[i].y;
            if((pairy[pairc] > pairy[pairc-1]) &&
               (pairy[pairc] > lodata2[i+1].y)) {
                probi = la->signi(&pairy[pairc], &nout, &tofac);
                if(pairy[pairc] > 4.) {
                    *gmsg << std::fixed
                          << std::setw(12) << std::setprecision(8) << pairx[pairc]*Norm << " "
                          << std::setw(8)  << std::setprecision(2) << pairy[pairc] << " "
                          << std::setw(8)  << std::setprecision(3) << probi << " "
                          << i << endl;
                }
            }
            pairc++;
        }
    }

    *gmsg << "* ===> Max: "
          << std::fixed
          << std::setw(12) << std::setprecision(8) << lodata2[jmax].x * Norm << " "
          << std::setw(8)  << std::setprecision(2) << lodata2[jmax].y << " "
          << endl;

    return(0);
}


int TUNE_class::lombAnalysis(double *x, double *y, int Ndat, int /*nhis*/)
/*-----------------------------------------------------------------------------
 *  Launch Lomb analysis and plot results
 *  =======================================
 *
 *---------------------------------------------------------------------------*/
{
    int    i, nout, jmax;
    int    pairc;
    int    datcnt = 0;
    int stat = 0;
    double prob, probi;
    double tofac = 0.8;

    LOMB_TYPE tlom;

    CI_lt p, q;

    std::vector<LOMB_TYPE> lodata, lodata2;
    /*---------------------------------------------------------------------------*/

    *gmsg << "* TUNE_class LombAnalysis requested" << endl;

    /*
     * Do Lomb analysis
     * ================
     */

    for(int j = 0; j < Ndat; j++) {
        tlom.x = x[j];
        tlom.y = y[j];
        lodata.push_back(tlom);
    }

    p = lodata.begin();
    q = lodata.end();

    datcnt = count_if(p, q, Lomb_eq(0.));

    if(datcnt > (q - p - 10)) {
        *gmsg << "* Just found " << datcnt << "data points that are == 0!" << endl;
        return(-1);
    }

    // this parameterset works ok in most cases.....
    ofac  = 4.0;
    hifac = 0.8;
    Qmin  = 0.2;
    Qmax  = 0.4;

    std::unique_ptr<LOMB_class> la(new LOMB_class(1));

    stat = 0;
    stat = la->period(&lodata, &lodata2, ofac, hifac, &nout, &jmax, &prob, 0);
    if(stat != 0) {
        *gmsg << "* @C3ERROR: Lomb analysis failed!" << endl;
        return(-1);
    }

    *gmsg << "* =====> jmax = " << jmax << endl;

    double pairx[nout];
    double pairy[nout];

    *gmsg << "* ********** Peaks in Data:       **************" << endl;

    /*
    ada make histogram
    hbook1(nhis,"Lomb data",nout,
       (float)lodata2[0].x,
       (float)lodata2[nout-1].x);

    */
    pairc = 0;
    for(i = 0; i < nout; i++) {
        /* ada book histogram
        Hf1(nhis,(float)lodata2[i].x,(float)lodata2[i].y);
        */
        if(lodata2[i].y > 2.) {
            pairx[pairc] = lodata2[i].x;
            pairy[pairc] = lodata2[i].y;
            if((pairy[pairc] > pairy[pairc-1]) &&
               (pairy[pairc] > lodata2[i+1].y)) {
                probi = la->signi(&pairy[pairc], &nout, &tofac);
                if(pairy[pairc] > 4.) {
                    *gmsg << std::fixed
                          << std::setw(12) << std::setprecision(8) << pairx[pairc] << " "
                          << std::setw(8)  << std::setprecision(2) << pairy[pairc] << " "
                          << std::setw(8)  << std::setprecision(3) << probi << " "
                          << i << endl;
                }
            }
            pairc++;
        }
    }

    *gmsg << "* ===> Max: "
          << std::fixed
          << std::setw(12) << std::setprecision(8) << lodata2[jmax].x << " "
          << std::setw(8)  << std::setprecision(2) << lodata2[jmax].y << " "
          << endl;

    return(0);

}