detscan.cpp

Go to the documentation of this file.

#include "MTTensor.h"
#include <fstream>
#include <iostream>
#include <complex>
#include "CUniformRNG.h"
#include "netcdfcpp.h"

using namespace std;
using namespace gplib;

int main()
  {
    /*const complex<double> xx(0.0645,0.0439);
     const complex<double> xy(0.0555,0.0402);
     const complex<double> yx(-0.0736,-0.0695);
     const complex<double> yy(-0.0702,-0.0409);*/

    /*const complex<double> xx(0.5,0.4);
     const complex<double> xy(0.5,0.4);
     const complex<double> yx(-0.5,-0.4);
     const complex<double> yy(-0.6,-0.4);*/

    const complex<double> xx(0.00047, 0.00014);
    const complex<double> xy(0.11653, 0.08244);
    const complex<double> yx(-0.00407, -0.01294);
    const complex<double> yy(0.00309, -0.00417);

    const MTTensor TestData(xx, xy, yx, yy, 0.001);
    const size_t ntestcases = 100;

    cerr << "Det Phi: " << TestData.GetdetPhi() << endl;
    const double imagmax = 2. * yx.imag();
    const double realmax = 2. * yx.real();
    const double imagmin = -2.0 * yx.imag();
    const double realmin = -2.0 * yx.real();
    const double realmin = -2.0 * yx.real();
    const double imagstep = (imagmax - imagmin) / ntestcases;
    const double realstep = (realmax - realmin) / ntestcases;
    cout << "App. res.: " << TestData.GetRhoxx() << " " << TestData.GetRhoxy()
        << " " << TestData.GetRhoyx() << " " << TestData.GetRhoyy() << endl;
    cout << "Phase: " << TestData.GetPhixx() << " " << TestData.GetPhixy()
        << " " << TestData.GetPhiyx() << " " << TestData.GetPhiyy() << endl;
    cout << "Phase tensor: " << TestData.GetPhi11() << " "
        << TestData.GetPhi12() << " " << TestData.GetPhi21() << " "
        << TestData.GetPhi22() << endl;
    cout << "GetPi1() " << TestData.GetPi1() << endl;
    cout << "GetPi2() " << TestData.GetPi2() << endl;
    cout << "GetPhiStrike() " << TestData.GetPhiStrike() << endl;
    cout << "GetPhiMaxNew() " << TestData.GetPhiMaxNew() << endl;
    cout << "GetPhiMinNew() " << TestData.GetPhiMinNew() << endl;
    cout << "GettrPhi() " << TestData.GettrPhi() << endl;
    cout << "GetskPhi() " << TestData.GetskPhi() << endl;
    cout << "GetdetPhi() " << TestData.GetdetPhi() << endl;
    cout << "GetPhi1() " << TestData.GetPhi1() << endl;
    cout << "GetPhi2() " << TestData.GetPhi2() << endl;
    cout << "GetDetreal() " << TestData.GetDetreal() << endl;

    cout << endl << endl;
    //init netcdf
    NcFile combrescdf("det.nc", NcFile::Replace);
    NcDim* xd = combrescdf.add_dim("real", ntestcases);
    NcDim* yd = combrescdf.add_dim("imag", ntestcases);
    NcVar* x = combrescdf.add_var("real", ncFloat, xd);
    NcVar* y = combrescdf.add_var("imag", ncFloat, yd);
    NcVar* phimin = combrescdf.add_var("PhiMin", ncFloat, xd, yd);
    NcVar* phimax = combrescdf.add_var("PhiMax", ncFloat, xd, yd);
    NcVar* alpha = combrescdf.add_var("Alpha", ncFloat, xd, yd);
    NcVar* beta = combrescdf.add_var("Beta", ncFloat, xd, yd);
    NcVar* eta = combrescdf.add_var("Eta", ncFloat, xd, yd);
    float *xvals = new float[xd->size()];
    float *yvals = new float[yd->size()];
    float *phiminvals = new float[xd->size() * yd->size()];
    float *phimaxvals = new float[xd->size() * yd->size()];
    float *alphavals = new float[xd->size() * yd->size()];
    float *betavals = new float[xd->size() * yd->size()];
    float *etavals = new float[xd->size() * yd->size()];
    /*for (double value = 0; value > maxvalue; value += step)
     {
     complex<double> newyx(yx.real(),value);
     MTTensor ScanData(xx,xy,newyx,yy);
     cout << ScanData.GetPhiyx() << " " <<  ScanData.GetdetPhi() << " " << ScanData.GetEta() << " "
     << ScanData.GetBeta_phi() *180.0/cos(-1.0) << " " <<  ScanData.GetAlpha_phi() *180.0/cos(-1.0) << endl;
     }*/
    ofstream iffile("det.txt");
    iffile << ntestcases << " " << ntestcases << " " << 1 << endl;
    for (int i = 0; i < ntestcases; ++i)
      {
        double real = realmin + i * realstep;
        xvals[i] = real;
        for (int j = 0; j < ntestcases; ++j)
          {
            double imag = imagmin + j * imagstep;
            yvals[j] = imag;
            complex<double> newyx(real, imag);
            MTTensor ScanData(xx, xy, newyx, yy, 0.001);
            phiminvals[i * ntestcases + j] = ScanData.GetPhiMin();
            phimaxvals[i * ntestcases + j] = ScanData.GetPhiMax();
            alphavals[i * ntestcases + j] = ScanData.GetAlpha_phi();
            betavals[i * ntestcases + j] = ScanData.GetBeta_phi();
            etavals[i * ntestcases + j] = ScanData.GetEta();
            iffile << ScanData.GetPhiMin() << endl;
            //cout << real << " " << imag <<  " " << ScanData.GetPhiyx() << " " << ScanData.GetRhoyx() << " "<<ScanData.GetdetPhi() << endl;
          }
        x->put(xvals, xd->size());
        y->put(yvals, yd->size());
        phimin->put(phiminvals, phimin->edges());
        phimax->put(phimaxvals, phimax->edges());
        alpha->put(alphavals, alpha->edges());
        beta->put(betavals, beta->edges());
        eta->put(etavals, eta->edges());
      }
    delete[] xvals;
    delete[] yvals;
    delete[] phiminvals;
    delete[] phimaxvals;
  }