testParametricXMLMieModel.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/ParametricXMLMieModel.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 atm;
using namespace fwk;
using namespace utl;
using namespace tst;


class testParametricXMLMieModel : public CppUnit::TestFixture {

  CPPUNIT_TEST_SUITE(testParametricXMLMieModel);
  CPPUNIT_TEST(testAttenuation);
  CPPUNIT_TEST(testScattering);
  CPPUNIT_TEST(testExpertMode);
  CPPUNIT_TEST_SUITE_END();

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

  void tearDown() { }

  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(350.*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(300.*nanometer), 0.915745));

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

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

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

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

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

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

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

    // vary angle
    //
    AttenuationResult attAB1=
      theAtm.EvaluateMieAttenuation(pointA, pointB, wlength);

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

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

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

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

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

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


    // new pair of points (at same altitude)
    //
    ScatteringResult scatBC1 =
      theAtm.EvaluateMieScattering(pointB, pointC,
                                   40.*degree, 1.*km, wlength);
    CPPUNIT_ASSERT(Verify<CloseTo>(1.e8*scatBC1.GetScatteringFactor().Y(400*nanometer),
                                 0.214726, 1.e-5));

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

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

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

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

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

  void
  testExpertMode()
  {
    ParametricXMLMieModel myModel;
    myModel.Init();

    CPPUNIT_ASSERT(Verify<Equal>(myModel.GetHorizAttLength(),24000.*meter));
    CPPUNIT_ASSERT(Verify<Equal>(myModel.GetMixHeight(),1000.*meter));
    CPPUNIT_ASSERT(Verify<Equal>(myModel.GetScaleHeight(),1000.*meter));
    CPPUNIT_ASSERT(Verify<Equal>(myModel.GetAngstromCoeff(),1.));

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

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

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

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

    myModel.SetHorizAttLength(100.*meter);
    myModel.SetScaleHeight(1000.*meter);

    CPPUNIT_ASSERT(Verify<Equal>(myModel.GetHorizAttLength(),100.*meter));
    CPPUNIT_ASSERT(Verify<Equal>(myModel.GetScaleHeight(),1000.*meter));
  }

};


CPPUNIT_TEST_SUITE_REGISTRATION(testParametricXMLMieModel);


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