testParametricXMLRayleighModel.cc

Go to the documentation of this file.

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

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

#include <fwk/CentralConfig.h>
#include <fwk/CoordinateSystemRegistry.h>

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

#include <cppunit/extensions/HelperMacros.h>
#include <tst/Verify.h>

#include <iostream>

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


class testParametricXMLRayleighModel : public CppUnit::TestFixture{

  CPPUNIT_TEST_SUITE(testParametricXMLRayleighModel);

  INFO("PLEASE IMPLEMENT ME");

//   CPPUNIT_TEST(testEvaluateAttenuation);
//   CPPUNIT_TEST(testEvaluateScattering);

  CPPUNIT_TEST_SUITE_END();

private:

public:

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

  void tearDown(){
  }

  void testEvaluateAttenuation() {

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

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

    //
    // 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.,   10000.*m,  CSll);
    Point Dll(0., 2000.*m,   2600.*m,   CSll);
    Point Ell(0.,20000.*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(350.*nanometer);
    wlength.push_back(400.*nanometer);


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


    AttenuationResult attABll = theAtm.EvaluateRayleighAttenuation(All,
                                                                   Bll,
                                                                   wlength);
    CPPUNIT_ASSERT(Verify<CloseTo>(attABll.GetTransmissionFactor().Y(300.*nanometer),0.771771));


    AttenuationResult attACll = theAtm.EvaluateRayleighAttenuation(All,
                                                                   Cll,
                                                                   wlength);
    CPPUNIT_ASSERT(Verify<CloseTo>(attACll.GetTransmissionFactor().Y(300.*nanometer), 0.498456));

    AttenuationResult attADll = theAtm.EvaluateRayleighAttenuation(All,
                                                                   Dll,
                                                                   wlength);
    CPPUNIT_ASSERT(Verify<CloseTo>(attADll.GetTransmissionFactor().Y(300.*nanometer),0.721197));

    AttenuationResult attBDll = theAtm.EvaluateRayleighAttenuation(Bll,
                                                                   Dll,
                                                                   wlength);
    CPPUNIT_ASSERT(Verify<CloseTo>(attBDll.GetTransmissionFactor().Y(300.*nanometer), 0.838598));

    AttenuationResult attBEll = theAtm.EvaluateRayleighAttenuation(Bll,
                                                                   Ell,
                                                                   wlength);
      CPPUNIT_ASSERT(Verify<CloseTo>(attBEll.GetTransmissionFactor().Y(300.*nanometer), 0.172627));

    //
    //
    // Calculate distance between sea level and a fixed height above sea level using two different
    // coordinate systems.  The results should be identical (modulo difference in
    // the vertical due to earth curvature).  This test turns up error if, for eg.,
    // you use GetZ(some CS) rather than the height above reference ellipsoid
    // in ParametricXMLRayleighModel::GDistance
    //
    double heightSLll = ( (ReferenceEllipsoid::Get(ReferenceEllipsoid::eWGS84)).
                           PointToLatitudeLongitudeHeight( Point(0., 0., 0., CSll)) ).get<2>();
    double heightSLst = ( (ReferenceEllipsoid::Get(ReferenceEllipsoid::eWGS84)).
                           PointToLatitudeLongitudeHeight( Point(0., 0., 0., CSst)) ).get<2>();

    AttenuationResult attCGroundLl
      = theAtm.EvaluateRayleighAttenuation(Point(0., 0., 0., CSll),
                                           Point(0., 0., 10000*m, CSll),
                                           wlength);
    AttenuationResult attCGroundSt
      = theAtm.EvaluateRayleighAttenuation(Point(0., 0., heightSLll - heightSLst, CSst),
                                           Point(0., 0., heightSLll - heightSLst + 10000.*m, CSst),
                                           wlength);

    CPPUNIT_ASSERT(Verify<CloseTo>( attCGroundLl.GetTransmissionFactor().Y(300.*nanometer),
                                  attCGroundSt.GetTransmissionFactor().Y(300.*nanometer)));

  }

  void testEvaluateScattering() {

    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(350*nanometer);
    wlength.push_back(400*nanometer);

    Point pointA(0., 0.    ,  0.,     CSll);
    Point pointB(0., 0.    ,  500.*m, CSll);
    Point pointC(0., 500.*m,  500.*m, CSll);

    // vary angle
    //
    ScatteringResult scatAB1 =
      theAtm.EvaluateRayleighScattering(pointA, pointB,
                                        0., 1.*km, wlength);
    CPPUNIT_ASSERT(Verify<CloseTo>(1.e9*scatAB1.GetScatteringFactor().Y(400.*nanometer), 2.06171, 1e-5));

    ScatteringResult scatAB2 =
      theAtm.EvaluateRayleighScattering(pointA, pointB,
                                        10.*degree, 1.*km, wlength);
    CPPUNIT_ASSERT(Verify<CloseTo>(1e9*scatAB2.GetScatteringFactor().Y(400.*nanometer), 2.03062, 1e-5));

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

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


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

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


    // new pair of points (at same altitude)
    //
    ScatteringResult scatBC1 =
      theAtm.EvaluateRayleighScattering(pointB, pointC,
                                   0., 1.*km, wlength);
    CPPUNIT_ASSERT(Verify<CloseTo>(1e9*scatBC1.GetScatteringFactor().Y(400*nanometer),
                                 1.99968, 1e-5));

    ScatteringResult scatBC2 =
      theAtm.EvaluateRayleighScattering(pointB, pointC,
                                   10.*degree, 1.*km, wlength);
     CPPUNIT_ASSERT(Verify<CloseTo>(1e9*scatBC2.GetScatteringFactor().Y(400*nanometer),
                                 1.96953, 1e-5));

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

    ScatteringResult scatBC4 =
      theAtm.EvaluateRayleighScattering(pointB, pointC,
                                   120.*degree, 1.*km, wlength);
    CPPUNIT_ASSERT(Verify<CloseTo>(1e9*scatBC4.GetScatteringFactor().Y(400*nanometer),
                                 1.2498, 1e-5));

    // TO DO : vary wavelengths
    //

    // check angle returned by ScatteringResult
    //
    CPPUNIT_ASSERT(Verify<CloseTo>(scatAB1.GetScatteringAngle(), 0.));
    CPPUNIT_ASSERT(Verify<CloseTo>(scatAB2.GetScatteringAngle(), 10.*degree));


  }

};

CPPUNIT_TEST_SUITE_REGISTRATION(testParametricXMLRayleighModel);



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