testMeasuredDBMieModel.cc

Go to the documentation of this file.

#include <cppunit/extensions/HelperMacros.h>

#include <det/Detector.h>
#include <fdet/FDetector.h>
#include <fdet/Eye.h>

#include <fwk/CoordinateSystemRegistry.h>

#include <atm/AerosolDB.h>
#include <atm/AttenuationResult.h>
#include <atm/ScatteringResult.h>
#include <atm/MeasuredDBMieModel.h>

#include <fwk/CentralConfig.h>

#include <utl/TimeStamp.h>
#include <utl/UTCDateTime.h>
#include <utl/ErrorLogger.h>
#include <utl/Point.h>
#include <utl/TabulatedFunctionErrors.h>
#include <utl/ReferenceEllipsoid.h>

#include <tst/Verify.h>

#include <iostream>

using namespace std;
using namespace det;
using namespace fdet;
using namespace atm;
using namespace fwk;
using namespace utl;
using namespace tst;


class testMeasuredDBMieModel : public CppUnit::TestFixture {

  CPPUNIT_TEST_SUITE(testMeasuredDBMieModel);

  CPPUNIT_TEST(testZoneSelection);
  CPPUNIT_TEST(testAttenuation);
  CPPUNIT_TEST(testLambdaDependence);
  CPPUNIT_TEST(testScattering);
  // Disable out-of-bound check (temporarily I home) as corrections
  // for bug report 418 prevent out-of-bound conditions by construction.
  //  CPPUNIT_TEST_EXCEPTION(testOutOfBound, OutOfBoundException);
  CPPUNIT_TEST(testHasData);
  CPPUNIT_TEST_EXCEPTION(testNoData, NoDataForModelException);
  CPPUNIT_TEST(testNoDataMultipleCalls); // checks bug 456
  CPPUNIT_TEST(testUncertainty);

  CPPUNIT_TEST_SUITE_END();

public:
  void
  setUp()
  {
    CentralConfig::GetInstance(BOOTSTRAPFILE);
    ErrorLogger::GetInstance().SetVerbosity(Verbosity::eVerbose);
  }

  void tearDown() { }

  void
  testZoneSelection()
  {
    Detector::GetInstance().Update(UTCDateTime(2005,1,1,0,0,0).GetTimeStamp());
    const Atmosphere& theAtm = Detector::GetInstance().GetAtmosphere();

    vector<double> wlength;
    wlength.push_back(300*nanometer);
    wlength.push_back(355*nanometer);
    wlength.push_back(400*nanometer);

    // Los Leones
    //
    CoordinateSystemPtr
      CSll(Detector::GetInstance().GetFDetector().GetEye("Los Leones").GetLocalCoordinateSystem() );

    Point All(0, 0, 0, CSll);
    Point Bll(0, 0, 9000*m, CSll);

    AttenuationResult AttABll = theAtm.EvaluateMieAttenuation(All, Bll, wlength);
    CPPUNIT_ASSERT(Verify<CloseTo>(AttABll.GetTransmissionFactor().Y(355*nanometer),0.979985, 1e-5));

    // Coihueco
    CoordinateSystemPtr
      CSco(Detector::GetInstance().GetFDetector().GetEye("Coihueco").GetLocalCoordinateSystem() );

    Point Aco(0, 0, 0,      CSco);
    Point Bco(0, 0, 9000*m, CSco);

    AttenuationResult AttABco = theAtm.EvaluateMieAttenuation(Aco, Bco, wlength);
    CPPUNIT_ASSERT(Verify<CloseTo>(AttABco.GetTransmissionFactor().Y(355.*nanometer),0.391086 , 1e-5));
  }

  void
  testAttenuation()
  {
    Detector::GetInstance().Update(UTCDateTime(2005,1,1,0,0,0).GetTimeStamp());

    CoordinateSystemPtr
      CSll(Detector::GetInstance().GetFDetector().GetEye("Los Leones").GetLocalCoordinateSystem() );

    ReferenceEllipsoid wgs84(ReferenceEllipsoid::Get(ReferenceEllipsoid::eWGS84));

    //
    // To test scattering and attenuation we define several
    // points in the sky and compute scattering/attenuation between
    // different pairs.  Points are also defined in >1 coordinate system
    // to ensure that the calculations are not CS dependent.
    //
    // 5 points are defined, arranged as indicated (schematically) below.
    //
    //    C
    //
    //
    //
    //
    //    B    D                                   E
    //
    //    A   -------- z=0 in local CS ----------
    //

    Point All(0,      0, 0,      CSll);
    Point Bll(0,      0, 2600*m, CSll);
    Point Cll(0,      0, 9000*m, CSll);
    Point Dll(0, 2000*m, 2600*m, CSll);
    Point Ell(0, 9000*m, 2600*m, CSll);

    CoordinateSystemPtr CSst( Detector::GetInstance().GetSiteCoordinateSystem() );
    Point Cst(0, 0, 10000*m, CSst);

    vector<double> wlength;
    wlength.push_back(300*nanometer);
    wlength.push_back(355*nanometer);
    wlength.push_back(400*nanometer);

    const Atmosphere& theAtm = Detector::GetInstance().GetAtmosphere();

    AttenuationResult attABll = theAtm.EvaluateMieAttenuation(All,
                                                              Bll,
                                                              wlength);
    CPPUNIT_ASSERT(Verify<CloseTo>(attABll.GetTransmissionFactor().Y(355*nanometer), 0.981456, 1e-3));

    AttenuationResult attACll = theAtm.EvaluateMieAttenuation(All,
                                                              Cll,
                                                              wlength);
    CPPUNIT_ASSERT(Verify<CloseTo>(attACll.GetTransmissionFactor().Y(300*nanometer), 0.979984, 1e-3));

    AttenuationResult attADll = theAtm.EvaluateMieAttenuation(All,
                                                              Dll,
                                                              wlength);
    CPPUNIT_ASSERT(Verify<CloseTo>(attADll.GetTransmissionFactor().Y(300*nanometer), 0.976664, 1e-3));


    // Following two checks need some additional attention.
    // Not clear if MeasuredDBMieModel handles case of
    // points close together in height in the correct way
    //
    AttenuationResult attBDll = theAtm.EvaluateMieAttenuation(Bll,
                                                              Dll,
                                                              wlength);
    CPPUNIT_ASSERT(Verify<CloseTo>(attBDll.GetTransmissionFactor().Y(300*nanometer), 0.997001, 1e-3));

    AttenuationResult attBEll = theAtm.EvaluateMieAttenuation(Bll,
                                                              Ell,
                                                              wlength);
    CPPUNIT_ASSERT(Verify<CloseTo>(attBEll.GetTransmissionFactor().Y(300*nanometer), 0.986609, 1e-3));


//     // Set parameters to check bug 397
//     //
//     CoordinateSystemPtr
//       csLL(Detector::GetInstance().GetFDetector().GetEye("Los Leones").GetEyeCoordinateSystem() );
//     Point Pa(0.*m, 0.*m,  0.*m,   csLL);
//     Point Pb(0.*m, 0.*m, 11000.*m, csLL);
//     AttenuationResult attPaPb = theAtm.EvaluateMieAttenuation(Pa, Pb, wlength);

//     cout << "tansmission fact " << attPaPb.GetTransmissionFactor().Y(300.*nanometer) << endl;

  }

  void
  testLambdaDependence()
  {
    Detector::GetInstance().Update(UTCDateTime(2005,1,1,0,0,0).GetTimeStamp());
    const Atmosphere& theAtm = Detector::GetInstance().GetAtmosphere();

    // Coihueco
    CoordinateSystemPtr
      CSco(Detector::GetInstance().GetFDetector().GetEye("Coihueco").GetLocalCoordinateSystem() );

    Point Aco(0, 0, 0*m,   CSco);
    Point Bco(0, 0, 200*m, CSco);

    vector<double> wlength;
    wlength.push_back(300*nanometer);
    wlength.push_back(355*nanometer);
    wlength.push_back(400*nanometer);

    AttenuationResult att = theAtm.EvaluateMieAttenuation(Aco,
                                                          Bco,
                                                          wlength);

    CPPUNIT_ASSERT(Verify<CloseTo>( log(att.GetTransmissionFactor().Y(375*nanometer)),
                                  log(att.GetTransmissionFactor().Y(355*nanometer)) * pow(375./355.,-1.2), 1e-3));

    CPPUNIT_ASSERT(Verify<CloseTo>( log(att.GetTransmissionFactor().Y(345.*nanometer)),
                                  log(att.GetTransmissionFactor().Y(355.*nanometer)) * pow(345./355.,-1.2), 1e-3));
  }

  void
  testScattering()
  {
    Detector::GetInstance().Update(UTCDateTime(2005,1,1,0,0,0).GetTimeStamp());
    const Atmosphere& theAtm = Detector::GetInstance().GetAtmosphere();
    CoordinateSystemPtr
      CSll(Detector::GetInstance().GetFDetector().GetEye("Los Leones").GetLocalCoordinateSystem() );

    vector<double> wlength;
    wlength.push_back(300*nanometer);
    wlength.push_back(355*nanometer);
    wlength.push_back(400*nanometer);

    Point pointA(0, 0,     0 ,     CSll);
    Point pointB(0, 0,     1400*m, CSll);
    Point pointC(0, 600*m, 1600*m, CSll);

    // vary angle
    //
    ScatteringResult scatAB1 =
      theAtm.EvaluateMieScattering(pointA, pointB,
                                   20*degree, 1*km, wlength);
    CPPUNIT_ASSERT(Verify<CloseTo>(1e8*scatAB1.GetScatteringFactor().Y(400*nanometer), 0.856146, 1e-5));

    ScatteringResult scatAB2 =
      theAtm.EvaluateMieScattering(pointA, pointB,
                                   30*degree, 1*km, wlength);
    CPPUNIT_ASSERT(Verify<CloseTo>(1e8*scatAB2.GetScatteringFactor().Y(400*nanometer), 0.430129, 1e-5));

    ScatteringResult scatAB3 =
      theAtm.EvaluateMieScattering(pointA, pointB,
                                   40*degree, 1*km, wlength);
    CPPUNIT_ASSERT(Verify<CloseTo>(1e9*scatAB3.GetScatteringFactor().Y(400*nanometer), 2.33831, 1e-5));

    ScatteringResult scatAB4 =
      theAtm.EvaluateMieScattering(pointA, pointB,
                                   60*degree, 1*km, wlength);
    CPPUNIT_ASSERT(Verify<CloseTo>(1e9*scatAB4.GetScatteringFactor().Y(400*nanometer), 0.852745, 1e-5));

    // vary distance
    //
    ScatteringResult scatAB5 =
      theAtm.EvaluateMieScattering(pointA, pointB,
                                   30*degree, 10*km, wlength);
    CPPUNIT_ASSERT(Verify<CloseTo>(1e11*scatAB5.GetScatteringFactor().Y(400*nanometer),
                                 4.30129, 1e-5));

    ScatteringResult scatAB6 =
      theAtm.EvaluateMieScattering(pointA, pointB,
                                   60*degree, 10*km, wlength);
    CPPUNIT_ASSERT(Verify<CloseTo>(1e11*scatAB6.GetScatteringFactor().Y(400*nanometer),
                                 0.852745, 1e-5));


    // new pair of points (at same altitude)
    // These need to be checked - not sure if MeasuredDBMieModel is doing this
    // in the best way
    //
    ScatteringResult scatBC1 =
      theAtm.EvaluateMieScattering(pointB, pointC,
                                   40*degree, 1*km, wlength);
    CPPUNIT_ASSERT(Verify<CloseTo>(1e8*scatBC1.GetScatteringFactor().Y(400*nanometer),
                                 0.0441425, 1e-5));

    ScatteringResult scatBC2 =
      theAtm.EvaluateMieScattering(pointB, pointC,
                                   50*degree, 1*km, wlength);
    CPPUNIT_ASSERT(Verify<CloseTo>(1e8*scatBC2.GetScatteringFactor().Y(400*nanometer),
                                 0.0259341, 1e-5));

    ScatteringResult scatBC3 =
      theAtm.EvaluateMieScattering(pointB, pointC,
                                   60*degree, 1*km, wlength);
    CPPUNIT_ASSERT(Verify<CloseTo>(1e10*scatBC3.GetScatteringFactor().Y(400*nanometer),
                                 1.60981, 1e-5));

    ScatteringResult scatBC4 =
      theAtm.EvaluateMieScattering(pointB, pointC,
                                   90*degree, 1*km, wlength);
    CPPUNIT_ASSERT(Verify<CloseTo>(1e10*scatBC4.GetScatteringFactor().Y(400*nanometer),
                                 0.477644, 1e-5));

    // vary wavelengths
    //
    CPPUNIT_ASSERT(Verify<CloseTo>(1e8*scatAB1.GetScatteringFactor().Y(300*nanometer), 0.856146, 1e-5));
    CPPUNIT_ASSERT(Verify<CloseTo>(1e8*scatAB1.GetScatteringFactor().Y(355*nanometer), 0.856146, 1e-5));

    // check angle returned
    //
    CPPUNIT_ASSERT(Verify<CloseTo>(scatAB1.GetScatteringAngle(), 20*degree));
    CPPUNIT_ASSERT(Verify<CloseTo>(scatAB2.GetScatteringAngle(), 30*degree));
  }

  void
  testOutOfBound()
  {
    Detector::GetInstance().Update(UTCDateTime(2005,1,1,0,0,0).GetTimeStamp());
    const Atmosphere& theAtm = Detector::GetInstance().GetAtmosphere();
    CoordinateSystemPtr
      CSll(Detector::GetInstance().GetFDetector().GetEye("Los Leones").GetLocalCoordinateSystem() );

    vector<double> wlength;
    wlength.push_back(300*nanometer);
    wlength.push_back(355*nanometer);
    wlength.push_back(400*nanometer);

    Point A(0, 0,  0,      CSll);
    Point B(0, 0,  500*km, CSll);

    AttenuationResult AttABll = theAtm.EvaluateMieAttenuation(A, B, wlength);
    AttABll.GetTransmissionFactor().Y(355*nanometer);
  }

  void
  testNoData()
  {
    Detector::GetInstance().Update(TimeStamp(123));  // no data exists for this time
    const Atmosphere& theAtm = Detector::GetInstance().GetAtmosphere();
    CoordinateSystemPtr
      CSll(Detector::GetInstance().GetFDetector().GetEye("Los Leones").GetLocalCoordinateSystem() );

    vector<double> wlength;
    wlength.push_back(300*nanometer);
    wlength.push_back(355*nanometer);
    wlength.push_back(400*nanometer);

    Point A(0, 0, 0,    CSll);
    Point B(0, 0, 1*km, CSll);

    AttenuationResult AttABll = theAtm.EvaluateMieAttenuation(A, B, wlength);
  }

  void
  testHasData()
  {
    // Test HasData method at the level of the MieModel
    //

    Detector::GetInstance().Update(TimeStamp(123));
    MeasuredDBMieModel mieModelWithoutData;
    mieModelWithoutData.Init();
    CPPUNIT_ASSERT(Verify<Equal>(mieModelWithoutData.HasData(), false));

    Detector::GetInstance().Update(UTCDateTime(2005,1,1,0,0,0).GetTimeStamp());
    MeasuredDBMieModel mieModelWithData;
    mieModelWithData.Init();
    CPPUNIT_ASSERT(Verify<Equal>(mieModelWithData.HasData(), true));
  }

  void
  testNoDataMultipleCalls()
  {
    Detector::GetInstance().Update(TimeStamp(123));  // no data exists for this time
    const Atmosphere& theAtm = Detector::GetInstance().GetAtmosphere();
    CoordinateSystemPtr
      CSll(Detector::GetInstance().GetFDetector().GetEye("Los Leones").GetLocalCoordinateSystem() );

    vector<double> wlength;
    wlength.push_back(300*nanometer);
    wlength.push_back(355*nanometer);
    wlength.push_back(400*nanometer);

    Point A(0, 0, 0,    CSll);
    Point B(0, 0, 1*km, CSll);

    // check for bug 456
    //
    int numExceptions = 0;
    for (int i=0 ; i<5 ; i++){
      try {
        AttenuationResult AttABll = theAtm.EvaluateMieAttenuation(A, B, wlength);
      } catch (NoDataForModelException& ex) {
        numExceptions++;
        // should catch multiple NoDataForModelExceptions
      }
    }
    CPPUNIT_ASSERT(Verify<Equal>(numExceptions, 5));
  }

  void testUncertainty() {
    Detector::GetInstance().Update(UTCDateTime(2005,1,1,0,0,0).GetTimeStamp());

    CoordinateSystemPtr
      CSll(Detector::GetInstance().GetFDetector().GetEye("Los Leones").GetLocalCoordinateSystem() );

    ReferenceEllipsoid wgs84(ReferenceEllipsoid::Get(ReferenceEllipsoid::eWGS84));

    Point All(0, 0, 0,      CSll);
    Point Bll(0, 0, 2600*m, CSll);

    vector<double> wlength;
    wlength.push_back(300*nanometer);
    wlength.push_back(355*nanometer);
    wlength.push_back(400*nanometer);

    const Atmosphere& theAtm = Detector::GetInstance().GetAtmosphere();
    AttenuationResult attABll = theAtm.EvaluateMieAttenuation(All,
                                                              Bll,
                                                              wlength);
    CPPUNIT_ASSERT(Verify<CloseTo>(attABll.GetTransmissionFactor().Y(355*nanometer), 0.9814559317));
    theAtm.SetUncertaintyBound(Atmosphere::eMie, 1);
    attABll = theAtm.EvaluateMieAttenuation(All,
                                            Bll,
                                            wlength);
    CPPUNIT_ASSERT(Verify<CloseTo>(attABll.GetTransmissionFactor().Y(355*nanometer), 0.9796204908));
    theAtm.SetUncertaintyBound(Atmosphere::eMie, -1);
    attABll = theAtm.EvaluateMieAttenuation(All,
                                            Bll,
                                            wlength);
    CPPUNIT_ASSERT(Verify<CloseTo>(attABll.GetTransmissionFactor().Y(355*nanometer), 0.9832947194));

  }

};


CPPUNIT_TEST_SUITE_REGISTRATION(testMeasuredDBMieModel);


// Configure (x)emacs for this file ...
// Local Variables:
// mode:c++
// compile-command: "make -C .. -k"
// End: