testLowLevelAtmInterface.cc

Go to the documentation of this file.

#include <det/Detector.h>

#include <fdet/FDetector.h>
#include <fdet/Eye.h>
#include <fdet/Telescope.h>
#include <fdet/Pixel.h>

#include <atm/AerosolDB.h>
#include <atm/OverallQualityDB.h>
#include <atm/AerosolZone.h>
#include <atm/AttSlice.h>
#include <atm/PFSlice.h>
#include <atm/MolecularZone.h>
#include <atm/LidarDB.h>
#include <atm/LidarZone.h>
#include <atm/VAODSlice.h>
#include <atm/MolecularIds.h>
#include <atm/MolecularDB.h>
#include <atm/GOESDB.h>

#include <fwk/CentralConfig.h>

#include <utl/ErrorLogger.h>
#include <utl/AugerUnits.h>
#include <utl/AugerException.h>
#include <utl/Point.h>
#include <fwk/CoordinateSystemRegistry.h>

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

#include <iostream>

// For some reason, CLHEP SystemOfUnits.h
// #defines pascal to hep_pascal, putting it in
// the global namespace and screwing up our own
// definition of pascal from AugerUnits.h.  It appears
// the CLHEP units are being pulled in via the geometry
// package.  The following line undoes the CLHEP sabotage.
#undef pascal

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

namespace {
  const double gNCloudProbabilityLevels = 5; // DB returns from [0, 4]
}

class testAtmLowLevelAtmInterface : public CppUnit::TestFixture {


  CPPUNIT_TEST_SUITE(testAtmLowLevelAtmInterface);

  CPPUNIT_TEST(testGOESDB);

  //
  // All the Checks are performed against quantities stored in the
  // validate_Atm_xyz databases.
  //

  // Test AerosolDB interface (zone iterators and whatnot).
  CPPUNIT_TEST(testAerosolDB);

  // Test iterating over and retrieving data from attenuation slices
  CPPUNIT_TEST(testAttSlice);

  // Test iterating over and retrieving data from phase function slices
  CPPUNIT_TEST(testPFSlice);

  // Test MolecularDB interface
  // THIS TEST HAS BEEN COMMENTED OUT AS IT WAS PREPARED FOR VERSION 0_A OF
  // THE MOLECULAR DATABASE. THIS VERSION USED SEASONAL PROFILES, WHICH HAVE
  // SINCE BEEN REPLACED BY GDAS AND WEATHER STATION MEASUREMENTS, WHICH ARE
  // WRITTEN IN VERSION 1_A OF THE MOLECULAR DB. THIS TEST SHOULD BE REWRITTEN
  // BY SOMEONE TO CHECK VERSION 1_A. THIS WOULD INVOLVE ALSO CREATING A SIMPLE
  // TEST DATABASE FOLLOWING THE SAME CONVENTIONS AS VERSION 1_A OF THE DATABASE.
  //  CPPUNIT_TEST(testMolecularDB);

  // Test iterating over and retrieveing data from molecular layers
  CPPUNIT_TEST(testMolecularLayer);

  // Test failure to find data in molecular db
  CPPUNIT_TEST(testMolecularDBNotFound);

  // Test case of no data in aerosol db
  CPPUNIT_TEST(testAerosolDBNotFound);

  // Test OverallQualityDB interface
  CPPUNIT_TEST(testOverallQualityDB);
  CPPUNIT_TEST_EXCEPTION(testOverallQualityDBNotFound, DataNotFoundInDBException);

  // Test LidarDB interface
  CPPUNIT_TEST(testLidarDB);

  // Test iterating over and retrieveing data from VAOD slices in Lidar db
  CPPUNIT_TEST(testVAODSlice);

  // Test the GOES cloud db
  CPPUNIT_TEST_EXCEPTION(testGOESDBNotFound1, DataNotFoundInDBException);
  CPPUNIT_TEST_EXCEPTION(testGOESDBNotFound2, DataNotFoundInDBException);
  CPPUNIT_TEST_EXCEPTION(testGOESDBNotFound3, NonExistentComponentException);

  // Test the cloud camera db.
  // Somebody should write a cloud db test.
  //  CPPUNIT_TEST(testCloudDB);

  CPPUNIT_TEST_SUITE_END();

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

  void tearDown() { }

  // Test interface to Atm_Aerosol database
  void
  testAerosolDB()
  {
    Detector& det = Detector::GetInstance();
    det.Update(1111112);

    const AerosolDB& aerosol1 = det.GetAtmosphere().GetAerosolDB();

    int numZonesAerosol1 = 0;
    for (AerosolDB::ZoneIterator zIt = aerosol1.ZonesBegin();
         zIt != aerosol1.ZonesEnd(); ++zIt) {

      if (zIt->GetName() == "Coihueco") {
        CPPUNIT_ASSERT(Verify<CloseTo>(zIt->GetNorthing(), 6114141.8*meter));
        CPPUNIT_ASSERT(Verify<CloseTo>(zIt->GetEasting(),  445346.*meter));
        ++numZonesAerosol1;
      }

      if (zIt->GetName() == "twilightZone") {
        CPPUNIT_ASSERT(Verify<CloseTo>(zIt->GetNorthing(), 144924.*meter));
        CPPUNIT_ASSERT(Verify<CloseTo>(zIt->GetEasting(),  1233122.*meter));
        ++numZonesAerosol1;
      }

      if (zIt->GetName() == "Los Morados") {
        CPPUNIT_ASSERT(Verify<CloseTo>(zIt->GetNorthing(), 6094571.*meter));
        CPPUNIT_ASSERT(Verify<CloseTo>(zIt->GetEasting(),  498885.*meter));
        ++numZonesAerosol1;
      }

    }
    CPPUNIT_ASSERT(Verify<Equal>(numZonesAerosol1, 3));  // should be three entries for time 1111112

    det.Update(2222223);

    const AerosolDB& aerosol2 = det.GetAtmosphere().GetAerosolDB();

    int numZonesAerosol2 = 0;
    for (AerosolDB::ZoneIterator zIt = aerosol2.ZonesBegin();
         zIt != aerosol2.ZonesEnd(); ++zIt)
      ++numZonesAerosol2;
    CPPUNIT_ASSERT(Verify<Equal>(numZonesAerosol2, 2));
  }

  void
  testAttSlice()
  {
    Detector& det = Detector::GetInstance();
    det.Update(1111112);            // the second aerosol object !
    const AerosolDB& aerosolDB = det.GetAtmosphere().GetAerosolDB();

    // iterate over zones
    for (AerosolDB::ZoneIterator zIt = aerosolDB.ZonesBegin();
         zIt != aerosolDB.ZonesEnd(); ++zIt) {

      double bogusValue = 0;
      int iSliceNumber = 1;

      // iterate over attenuation slices for this zone
      for (AerosolZone::AttSliceIterator attIt = zIt->AttSlicesBegin();
           attIt != zIt->AttSlicesEnd(); ++attIt) {

        if (zIt->GetName() == "Coihueco") {
          CPPUNIT_ASSERT(Verify<CloseTo>(attIt->GetMinHeight(), bogusValue*meter));
          CPPUNIT_ASSERT(Verify<CloseTo>(attIt->GetMaxHeight(), (bogusValue+10.)*meter));

          CPPUNIT_ASSERT(Verify<CloseTo>(attIt->GetVAOD(), (bogusValue+11.)));
          CPPUNIT_ASSERT(Verify<CloseTo>(attIt->GetMinVAOD(), (bogusValue+12.)));
          CPPUNIT_ASSERT(Verify<CloseTo>(attIt->GetMaxVAOD(), (bogusValue+13.)));

          CPPUNIT_ASSERT(Verify<CloseTo>(attIt->GetAttAlpha(), (77.0 + bogusValue)/meter));
          CPPUNIT_ASSERT(Verify<CloseTo>(attIt->GetMinAttAlpha(), (0.77 - bogusValue/100.)/meter));
          CPPUNIT_ASSERT(Verify<CloseTo>(attIt->GetMaxAttAlpha(), (0.77 + bogusValue/100.)/meter));

          CPPUNIT_ASSERT(Verify<CloseTo>(attIt->GetAttLambda(), 78.0 + bogusValue));
          CPPUNIT_ASSERT(Verify<CloseTo>(attIt->GetMinAttLambda(), 0.78 - bogusValue/100.));
          CPPUNIT_ASSERT(Verify<CloseTo>(attIt->GetMaxAttLambda(), 0.77 + bogusValue/100.));
        }

        if (zIt->GetName() == "twilightZone") {
          CPPUNIT_ASSERT(Verify<CloseTo>(attIt->GetMinHeight(), 0.*m));
          CPPUNIT_ASSERT(Verify<CloseTo>(attIt->GetMaxHeight(), 10000.*m));
          CPPUNIT_ASSERT(Verify<CloseTo>(attIt->GetVAOD(), 10000.));
          CPPUNIT_ASSERT(Verify<CloseTo>(attIt->GetMinVAOD(), 10.));
          CPPUNIT_ASSERT(Verify<CloseTo>(attIt->GetMaxVAOD(), 10.));


          CPPUNIT_ASSERT(Verify<CloseTo>(attIt->GetAttAlpha(), 10./m));
          CPPUNIT_ASSERT(Verify<CloseTo>(attIt->GetMinAttAlpha(), 1./m));
          CPPUNIT_ASSERT(Verify<CloseTo>(attIt->GetMaxAttAlpha(), 2./m));

          CPPUNIT_ASSERT(Verify<CloseTo>(attIt->GetAttLambda(), 20.0));
          CPPUNIT_ASSERT(Verify<CloseTo>(attIt->GetMinAttLambda(), 2.));
          CPPUNIT_ASSERT(Verify<CloseTo>(attIt->GetMaxAttLambda(), 3.));
        }

        if (zIt->GetName() == "Los Morados") {
          CPPUNIT_ASSERT(Verify<CloseTo>(attIt->GetMinHeight(), 0.*meter));
          CPPUNIT_ASSERT(Verify<CloseTo>(attIt->GetMaxHeight(), 1000.*meter));
          CPPUNIT_ASSERT(Verify<CloseTo>(attIt->GetVAOD(), 1000. ));
          CPPUNIT_ASSERT(Verify<CloseTo>(attIt->GetMinVAOD(), 10.));
          CPPUNIT_ASSERT(Verify<CloseTo>(attIt->GetMaxVAOD(), 10.));

          CPPUNIT_ASSERT(Verify<CloseTo>(attIt->GetAttAlpha(), 100./m));
          CPPUNIT_ASSERT(Verify<CloseTo>(attIt->GetMinAttAlpha(), 10./m));
          CPPUNIT_ASSERT(Verify<CloseTo>(attIt->GetMaxAttAlpha(), 11./m));

          CPPUNIT_ASSERT(Verify<CloseTo>(attIt->GetAttLambda(), 200.));
          CPPUNIT_ASSERT(Verify<CloseTo>(attIt->GetMinAttLambda(), 20.));
          CPPUNIT_ASSERT(Verify<CloseTo>(attIt->GetMaxAttLambda(), 21.));
        }
        ++iSliceNumber;
        bogusValue += 10;
      }

    }  // iterating on zones
  }

  void
  testPFSlice()
  {
    Detector& theDet = Detector::GetInstance();

    theDet.Update(1111112); // the second aerosol object !

    const Atmosphere& theAtm = theDet.GetAtmosphere();
    const AerosolDB& aerosolDB = theAtm.GetAerosolDB();

    for (AerosolDB::ZoneIterator zIt = aerosolDB.ZonesBegin();
         zIt != aerosolDB.ZonesEnd(); ++zIt) {

      for (AerosolZone::PFSliceIterator pfIt = zIt->PFSlicesBegin();
           pfIt != zIt->PFSlicesEnd(); ++pfIt) {

        if (zIt->GetName() == "Coihueco") {
          CPPUNIT_ASSERT(Verify<CloseTo>(pfIt->GetMinHeight(), 0.*meter));
          CPPUNIT_ASSERT(Verify<CloseTo>(pfIt->GetMaxHeight(), 10000.*meter));

          CPPUNIT_ASSERT(Verify<CloseTo>(pfIt->GetF(), 1.0));
          CPPUNIT_ASSERT(Verify<CloseTo>(pfIt->GetFError(), 0.2));
          CPPUNIT_ASSERT(Verify<CloseTo>(pfIt->GetG(), 2.0));
          CPPUNIT_ASSERT(Verify<CloseTo>(pfIt->GetGError(), 0.3));

          CPPUNIT_ASSERT(Verify<CloseTo>(pfIt->GetFLambda(), 4.));
          CPPUNIT_ASSERT(Verify<CloseTo>(pfIt->GetFLambdaError(), 0.4));
          CPPUNIT_ASSERT(Verify<CloseTo>(pfIt->GetGLambda(), 5.));
          CPPUNIT_ASSERT(Verify<CloseTo>(pfIt->GetGLambdaError(), 0.5));
        }

        if (zIt->GetName() == "twilightZone") {
          CPPUNIT_ASSERT(Verify<CloseTo>(pfIt->GetMinHeight(), 0.*m));
          CPPUNIT_ASSERT(Verify<CloseTo>(pfIt->GetMaxHeight(), 11000.*m));

          CPPUNIT_ASSERT(Verify<CloseTo>(pfIt->GetF(), 66.));
          CPPUNIT_ASSERT(Verify<CloseTo>(pfIt->GetFError(), 0.66));
          CPPUNIT_ASSERT(Verify<CloseTo>(pfIt->GetG(), 77.));
          CPPUNIT_ASSERT(Verify<CloseTo>(pfIt->GetGError(), 0.77));

          CPPUNIT_ASSERT(Verify<CloseTo>(pfIt->GetFLambda(), 88.));
          CPPUNIT_ASSERT(Verify<CloseTo>(pfIt->GetFLambdaError(), 0.88));
          CPPUNIT_ASSERT(Verify<CloseTo>(pfIt->GetGLambda(), 99.));
          CPPUNIT_ASSERT(Verify<CloseTo>(pfIt->GetGLambdaError(), 0.99));
        }

        if (zIt->GetName() == "Los Morados") {
          CPPUNIT_ASSERT(Verify<CloseTo>(pfIt->GetMinHeight(), 0.*m));
          CPPUNIT_ASSERT(Verify<CloseTo>(pfIt->GetMaxHeight(), 1100.*m));

          CPPUNIT_ASSERT(Verify<CloseTo>(pfIt->GetF(), 66.));
          CPPUNIT_ASSERT(Verify<CloseTo>(pfIt->GetFError(), 0.66));
          CPPUNIT_ASSERT(Verify<CloseTo>(pfIt->GetG(), 77.));
          CPPUNIT_ASSERT(Verify<CloseTo>(pfIt->GetGError(), 0.77));

          CPPUNIT_ASSERT(Verify<CloseTo>(pfIt->GetFLambda(), 88.));
          CPPUNIT_ASSERT(Verify<CloseTo>(pfIt->GetFLambdaError(), 0.88));
          CPPUNIT_ASSERT(Verify<CloseTo>(pfIt->GetGLambda(), 99.));
          CPPUNIT_ASSERT(Verify<CloseTo>(pfIt->GetGLambdaError(), 0.99));
        }
      }
    }
  }

  void
  testMolecularDB()
  {
    Detector& det = Detector::GetInstance();
    det.Update(1111112);

    const MolecularDB& molecular1 = det.GetAtmosphere().GetMolecularDB(MolecularIds::eWeatherStation);

    CPPUNIT_ASSERT(molecular1.GetProfileName() == "weather station");
    int numZonesMolecular1 = 0;
    for (MolecularDB::ZoneIterator zIt = molecular1.ZonesBegin();
         zIt != molecular1.ZonesEnd(); ++zIt) {

      if (zIt->GetName() == "CLF") {
        CPPUNIT_ASSERT(Verify<CloseTo>(zIt->GetNorthing(), 6095769.0*meter));
        CPPUNIT_ASSERT(Verify<CloseTo>(zIt->GetEasting(),  469378.0*meter));
      }
      ++numZonesMolecular1;
    }

    CPPUNIT_ASSERT(Verify<Equal>(numZonesMolecular1, 1));

    det.Update(2222223);
    const MolecularDB& molecular2 = det.GetAtmosphere().GetMolecularDB(MolecularIds::eUSStdA);

    CPPUNIT_ASSERT(molecular2.GetProfileName() == "US-StdA");
    int numZonesMolecular2 = 0;
    for (MolecularDB::ZoneIterator zIt = molecular2.ZonesBegin();
         zIt != molecular2.ZonesEnd(); ++zIt)
      ++numZonesMolecular2;

    CPPUNIT_ASSERT(Verify<Equal>(numZonesMolecular2, 2));

    // test lookup for different profile types
    //
    det.Update(3333334);

    const MolecularDB& m1 = det.GetAtmosphere().GetMolecularDB(MolecularIds::eMalargueJanuaryModel);
    CPPUNIT_ASSERT(m1.GetProfileName() == string("Malargue January model"));
    const MolecularDB& m2 = det.GetAtmosphere().GetMolecularDB(MolecularIds::eMalargueFebruaryModel);
    CPPUNIT_ASSERT(m2.GetProfileName() == string("Malargue February model"));
    const MolecularDB& m3 = det.GetAtmosphere().GetMolecularDB(MolecularIds::eMalargueMarchModel);
    CPPUNIT_ASSERT(m3.GetProfileName() == string("Malargue March model"));
    const MolecularDB& m4 = det.GetAtmosphere().GetMolecularDB(MolecularIds::eMalargueAprilModel);
    CPPUNIT_ASSERT(m4.GetProfileName() == string("Malargue April model"));
    const MolecularDB& m5 = det.GetAtmosphere().GetMolecularDB(MolecularIds::eMalargueMayModel);
    CPPUNIT_ASSERT(m5.GetProfileName() == string("Malargue May model"));
    const MolecularDB& m6 = det.GetAtmosphere().GetMolecularDB(MolecularIds::eMalargueJuneModel);
    CPPUNIT_ASSERT(m6.GetProfileName() == string("Malargue June model"));
    const MolecularDB& m7 = det.GetAtmosphere().GetMolecularDB(MolecularIds::eMalargueJulyModel);
    CPPUNIT_ASSERT(m7.GetProfileName() == string("Malargue July model"));
    const MolecularDB& m8 = det.GetAtmosphere().GetMolecularDB(MolecularIds::eMalargueAugustModel);
    CPPUNIT_ASSERT(m8.GetProfileName() == string("Malargue August model"));
    const MolecularDB& m9 = det.GetAtmosphere().GetMolecularDB(MolecularIds::eMalargueSeptemberModel);
    CPPUNIT_ASSERT(m9.GetProfileName() == string("Malargue September model"));
    const MolecularDB& m10 = det.GetAtmosphere().GetMolecularDB(MolecularIds::eMalargueOctoberModel);
    CPPUNIT_ASSERT(m10.GetProfileName() == string("Malargue October model"));
    const MolecularDB& m11 = det.GetAtmosphere().GetMolecularDB(MolecularIds::eMalargueNovemberModel);
    CPPUNIT_ASSERT(m11.GetProfileName() == string("Malargue November model"));
    const MolecularDB& m12 = det.GetAtmosphere().GetMolecularDB(MolecularIds::eMalargueDecemberModel);
    CPPUNIT_ASSERT(m12.GetProfileName() == string("Malargue December model"));
  }

  void
  testMolecularLayer()
  {
    Detector& det = Detector::GetInstance();
    det.Update(1111112);
    const MolecularDB& molecularDB = det.GetAtmosphere().GetMolecularDB(MolecularIds::eWeatherStation);

    // iterate over zones
    for (MolecularDB::ZoneIterator zIt = molecularDB.ZonesBegin();
         zIt != molecularDB.ZonesEnd(); ++zIt) {

      [[maybe_unused]] double bogusValue = 0;
      int iLayerNumber  = 1;

      // iterate over molecular layers for this zone
      for (MolecularZone::LayerIterator sIt = zIt->LayersBegin();
           sIt != zIt->LayersEnd(); ++sIt) {

        if (zIt->GetName() == "Coihueco") {

          const double dummy = iLayerNumber;
          const double dummy_err = dummy/10.;

          CPPUNIT_ASSERT(Verify<CloseTo>(sIt->GetHeight(), dummy*meter));
          CPPUNIT_ASSERT(Verify<CloseTo>(sIt->GetHeightError(), dummy_err*meter));
          CPPUNIT_ASSERT(Verify<CloseTo>(sIt->GetDepth(), dummy*g/cm2));
          CPPUNIT_ASSERT(Verify<CloseTo>(sIt->GetDepthError(), dummy_err*g/cm2));
          CPPUNIT_ASSERT(Verify<CloseTo>(sIt->GetPressure(), dummy*100.*pascal));
          CPPUNIT_ASSERT(Verify<CloseTo>(sIt->GetPressureError(), dummy_err*100.*pascal));
          CPPUNIT_ASSERT(Verify<CloseTo>(sIt->GetTemperature(), dummy*kelvin));
          CPPUNIT_ASSERT(Verify<CloseTo>(sIt->GetTemperatureError(), dummy_err*kelvin));
          CPPUNIT_ASSERT(Verify<CloseTo>(sIt->GetHumidity(), dummy*percent));
          CPPUNIT_ASSERT(Verify<CloseTo>(sIt->GetHumidityError(), dummy_err*percent));
          CPPUNIT_ASSERT(Verify<CloseTo>(sIt->GetWindSpeed(), dummy*m/s));
          CPPUNIT_ASSERT(Verify<CloseTo>(sIt->GetWindSpeedError(), dummy_err*m/s));
          CPPUNIT_ASSERT(Verify<CloseTo>(sIt->GetAirDensity(), dummy*kg/m3));
          CPPUNIT_ASSERT(Verify<CloseTo>(sIt->GetAirDensityError(), dummy_err*kg/m3));
        }

        if (zIt->GetName() == "twilightZone") {
          CPPUNIT_ASSERT(Verify<CloseTo>(sIt->GetHeight(), 1.*m));
          CPPUNIT_ASSERT(Verify<CloseTo>(sIt->GetHeightError(), .1*m));
          CPPUNIT_ASSERT(Verify<CloseTo>(sIt->GetDepth(), 2.*g/cm2));
          CPPUNIT_ASSERT(Verify<CloseTo>(sIt->GetDepthError(), .2*g/cm2));
          CPPUNIT_ASSERT(Verify<CloseTo>(sIt->GetPressure(), 3.*100.*pascal));
          CPPUNIT_ASSERT(Verify<CloseTo>(sIt->GetPressureError(), .3*100.*pascal));
          CPPUNIT_ASSERT(Verify<CloseTo>(sIt->GetTemperature(), 4.*kelvin));
          CPPUNIT_ASSERT(Verify<CloseTo>(sIt->GetTemperatureError(), .4*kelvin));
          CPPUNIT_ASSERT(Verify<CloseTo>(sIt->GetHumidity(), 5.*percent));
          CPPUNIT_ASSERT(Verify<CloseTo>(sIt->GetHumidityError(), .5*percent));
          CPPUNIT_ASSERT(Verify<CloseTo>(sIt->GetWindSpeed(), 6.*m/s));
          CPPUNIT_ASSERT(Verify<CloseTo>(sIt->GetWindSpeedError(), .6*m/s));
          CPPUNIT_ASSERT(Verify<CloseTo>(sIt->GetAirDensity(), 7.*kg/m3));
          CPPUNIT_ASSERT(Verify<CloseTo>(sIt->GetAirDensityError(), .7*kg/m3));
        }

        if (zIt->GetName() == "Los Morados") {
          CPPUNIT_ASSERT(Verify<CloseTo>(sIt->GetHeight(), 11.*m));
          CPPUNIT_ASSERT(Verify<CloseTo>(sIt->GetHeightError(), 1.*m));
          CPPUNIT_ASSERT(Verify<CloseTo>(sIt->GetDepth(), 12.*g/cm2));
          CPPUNIT_ASSERT(Verify<CloseTo>(sIt->GetDepthError(), 2.*g/cm2));
          CPPUNIT_ASSERT(Verify<CloseTo>(sIt->GetPressure(), 0.*100.*pascal));
          CPPUNIT_ASSERT(Verify<CloseTo>(sIt->GetPressureError(), 0.*100.*pascal));
          CPPUNIT_ASSERT(Verify<CloseTo>(sIt->GetTemperature(), 0.*kelvin));
          CPPUNIT_ASSERT(Verify<CloseTo>(sIt->GetTemperatureError(), 0.*kelvin));
          CPPUNIT_ASSERT(Verify<CloseTo>(sIt->GetHumidity(), 0.));
          CPPUNIT_ASSERT(Verify<CloseTo>(sIt->GetHumidityError(), 0.));
          CPPUNIT_ASSERT(Verify<CloseTo>(sIt->GetWindSpeed(), 0.*m/s));
          CPPUNIT_ASSERT(Verify<CloseTo>(sIt->GetWindSpeedError(), 0.*m/s));
          CPPUNIT_ASSERT(Verify<CloseTo>(sIt->GetAirDensity(), 0.*kg/m3));
          CPPUNIT_ASSERT(Verify<CloseTo>(sIt->GetAirDensityError(), 0.*kg/m3));
        }
        ++iLayerNumber;
        bogusValue += 10;
      }
    }
  }

  void
  testMolecularDBNotFound()
  {
    Detector& det = Detector::GetInstance();
    det.Update(0);
    const MolecularDB& molecular1 =
      det.GetAtmosphere().GetMolecularDB(MolecularIds::eWeatherStation);

    int nzones = 0;
    for (MolecularDB::ZoneIterator zIt = molecular1.ZonesBegin();
         zIt != molecular1.ZonesEnd(); ++zIt)
      ++nzones;

    CPPUNIT_ASSERT(Verify<Equal>(nzones, 0));
  }

  void
  testAerosolDBNotFound()
  {
    Detector& det = Detector::GetInstance();
    det.Update(0);

    const AerosolDB& aerosolDB = det.GetAtmosphere().GetAerosolDB();

    // confirm that no zones exist

    int nzones = 0;
    for (AerosolDB::ZoneIterator zIt = aerosolDB.ZonesBegin();
         zIt != aerosolDB.ZonesEnd(); ++zIt)
      ++nzones;

    CPPUNIT_ASSERT(Verify<Equal>(nzones, 0));
  }

  void
  testOverallQualityDB()
  {
    Detector& det = Detector::GetInstance();
    det.Update(691901819);

    const OverallQualityDB& overallDB =
      det.GetAtmosphere().GetOverallQualityDB();

    CPPUNIT_ASSERT(Verify<CloseTo>(overallDB.GetHorizontalUniformity(), 0.9));
    CPPUNIT_ASSERT(Verify<CloseTo>(overallDB.GetCloudCoverage(), 0.1));
    CPPUNIT_ASSERT(Verify<CloseTo>(overallDB.GetMinCloudBase(), 1000*meter));

    det.Update(1383803637);
    const OverallQualityDB& overallDB2 = det.GetAtmosphere().GetOverallQualityDB();

    CPPUNIT_ASSERT(Verify<CloseTo>(overallDB2.GetHorizontalUniformity(), 0.8));
    CPPUNIT_ASSERT(Verify<CloseTo>(overallDB2.GetCloudCoverage(),        0.2));
    CPPUNIT_ASSERT(Verify<CloseTo>(overallDB2.GetMinCloudBase(),         1000*meter));
  }

  void
  testOverallQualityDBNotFound()
  {
    Detector& det = Detector::GetInstance();
    det.Update(0);

    const OverallQualityDB& overallDB =
      det.GetAtmosphere().GetOverallQualityDB();
    double cc = overallDB.GetCloudCoverage();
    ++cc; // avoid gcc warning;
  }

  void
  testLidarDB()
  {
    Detector& det = Detector::GetInstance();
    det.Update(2222223);

    const LidarDB& lidar = det.GetAtmosphere().GetLidarDB();

    int numLidarZones = 0;

    for (LidarDB::ZoneIterator zIt = lidar.ZonesBegin() ;
         zIt != lidar.ZonesEnd(); ++zIt) {

      if (zIt->GetName() == "Coihueco") {
        CPPUNIT_ASSERT(Verify<CloseTo>(zIt->GetNorthing(), 6114141.8*meter));
        CPPUNIT_ASSERT(Verify<CloseTo>(zIt->GetEasting(),  445346.*meter));

        CPPUNIT_ASSERT(Verify<CloseTo>(zIt->GetCloudCoverage(), .2));
        CPPUNIT_ASSERT(Verify<CloseTo>(zIt->GetLowestCloudHeight(), 2482*meter));
        CPPUNIT_ASSERT(Verify<CloseTo>(zIt->GetLowestCloudThickness(), 449*meter));
        CPPUNIT_ASSERT(Verify<CloseTo>(zIt->GetLowestCloudVAOD(), 0.168));
        CPPUNIT_ASSERT(Verify<CloseTo>(zIt->GetLidarMaxHeight(), 8000*meter));

        ++numLidarZones;
      }

      if (zIt->GetName() == "myZone") {
        CPPUNIT_ASSERT(Verify<CloseTo>(zIt->GetNorthing(), 6083819*meter));
        CPPUNIT_ASSERT(Verify<CloseTo>(zIt->GetEasting(),  464292.5*meter));

        CPPUNIT_ASSERT(Verify<CloseTo>(zIt->GetCloudCoverage(), .3));
        CPPUNIT_ASSERT(Verify<CloseTo>(zIt->GetLowestCloudHeight(), 1000*meter));
        CPPUNIT_ASSERT(Verify<CloseTo>(zIt->GetLowestCloudThickness(), 200*meter));
        CPPUNIT_ASSERT(Verify<CloseTo>(zIt->GetLowestCloudVAOD(), 0.2));

        ++numLidarZones;
      }
    }
    CPPUNIT_ASSERT(Verify<Equal>(numLidarZones, 2));
  }

  void
  testVAODSlice()
  {
    Detector& det = Detector::GetInstance();
    det.Update(2222223);
    const LidarDB& lidarDB = det.GetAtmosphere().GetLidarDB();

    // iterate over zones
    for (LidarDB::ZoneIterator zIt = lidarDB.ZonesBegin() ;
         zIt != lidarDB.ZonesEnd(); ++zIt) {

      // ref the numbers written by the test program testAtmLidarHandler
      double bogusValue = 0.;
      int iSliceNumber  = 1;

      // iterate over VAOD slices for this zone
      for (LidarZone::VAODSliceIterator vaodIt = zIt->VAODSlicesBegin();
           vaodIt != zIt->VAODSlicesEnd(); ++vaodIt) {

        if (zIt->GetName() == "Coihueco") {
          CPPUNIT_ASSERT(Verify<CloseTo>(vaodIt->GetMinHeight(), bogusValue*meter));
          CPPUNIT_ASSERT(Verify<CloseTo>(vaodIt->GetMaxHeight(), (bogusValue+10.)*meter));
          CPPUNIT_ASSERT(Verify<CloseTo>(vaodIt->GetVAOD()     , bogusValue+11.));
          CPPUNIT_ASSERT(Verify<CloseTo>(vaodIt->GetMinVAOD()  , bogusValue+12.));
          CPPUNIT_ASSERT(Verify<CloseTo>(vaodIt->GetMaxVAOD()  , bogusValue+13.));
        }

        if (zIt->GetName() == "myZone") {
          CPPUNIT_ASSERT(Verify<CloseTo>(vaodIt->GetMinHeight(), 1.*meter));
          CPPUNIT_ASSERT(Verify<CloseTo>(vaodIt->GetMaxHeight(), 1000.*meter));
          CPPUNIT_ASSERT(Verify<CloseTo>(vaodIt->GetVAOD(), 0.5));
          CPPUNIT_ASSERT(Verify<CloseTo>(vaodIt->GetMinVAOD(), 0.4));
          CPPUNIT_ASSERT(Verify<CloseTo>(vaodIt->GetMaxVAOD(), 0.6));
        }
        ++iSliceNumber;
        bogusValue += 10;
      }
    }
  }


  void
  testCloudDB()
  {
    Detector& det = Detector::GetInstance();

    det.Update(691901819);
    {
      const fdet::Pixel& pix1 =
        Detector::GetInstance().GetFDetector().GetEye("Los Leones").GetTelescope(1).GetPixel(1);
      if (pix1.HasCloudFraction())
        CPPUNIT_ASSERT(Verify<CloseTo>(double(pix1.GetCloudFraction()), 0.4));
    }

    det.Update(691904871);
    {
      const fdet::Pixel& pix2 =
        Detector::GetInstance().GetFDetector().GetEye("Coihueco").GetTelescope(3).GetPixel(3);
      if (pix2.HasCloudFraction())
        CPPUNIT_ASSERT(Verify<CloseTo>(double(pix2.GetCloudFraction()), 1.));
    }
  }

  void
  testGOESDB()
  {
    Detector& det = Detector::GetInstance();

    det.Update(1);
    {
      const GOESDB& goesdb = det.GetAtmosphere().GetGOESDB();
      CPPUNIT_ASSERT(Verify<Equal>(goesdb.HasData(0), false));
    }

    int t = 692904871;
    det.Update(t);
    {
      const GOESDB& goesDB = det.GetAtmosphere().GetGOESDB();

      // get a pixel by id
      for (unsigned int idx = 0; idx <= 100; ++idx) {
        const double cloudP = goesDB.GetCloudProbability(idx);
        CPPUNIT_ASSERT(Verify<Equal>(cloudP, int(idx*4/100)/gNCloudProbabilityLevels));
      }

      // do it again (now pixels should be cached)
      for (unsigned int idx = 0; idx <= 100; ++idx) {
        const double cloudP = goesDB.GetCloudProbability(idx);
        CPPUNIT_ASSERT(Verify<Equal>(cloudP, int(idx*4/100)/gNCloudProbabilityLevels));
      }
    }

    // reset detector and get a pixel by id (at new time)
    t += 3700;
    det.Update(t);
    {
      const GOESDB& goesDB = det.GetAtmosphere().GetGOESDB();

      for (unsigned int idx = 0; idx <= 100; ++idx) {
        double cloudP = goesDB.GetCloudProbability(idx);
        CPPUNIT_ASSERT(Verify<Equal>(cloudP, 3./gNCloudProbabilityLevels));
      }
    }

    // get a pixel by location.  This uses the same test
    // data used to fill the GOES validation database.
    // The code to fill the db is in
    // OfflineTools/DBTools/Documentation/Examples/WriteToAtmGOESDB

    CoordinateSystemPtr csCLF = CoordinateSystemRegistry::Get("CLF");
    const Point pCLF = Point(0, 0, 0, csCLF);

    vector<int> prob;
    vector<unsigned int> pixelId;

    det.Update(700);
    {
      const GOESDB& goesDB = det.GetAtmosphere().GetGOESDB();

      Reader r(TESTGOES, Reader::eNONE);
      Branch tb = r.GetTopBranch();

      tb.GetChild("pixelId").GetData(pixelId);
      tb.GetChild("prob").GetData(prob);

      for (vector<unsigned int>::iterator i = pixelId.begin(); i != pixelId.end(); ++i) {
        CPPUNIT_ASSERT(Verify<Equal>(prob[*i]/gNCloudProbabilityLevels, goesDB.GetCloudProbability(*i)));
      }

      // get a pixel over the CLF (to confirm I have the basic code correct in
      // GetCloudProbability(Point) - that is, there should be some cloud data
      // over the CLF

      CPPUNIT_ASSERT(Verify<Equal>(goesDB.GetCloudProbability(pCLF), 2./gNCloudProbabilityLevels));
    }

    // starting from Malargue point and the magic numbers in GetCloudProbability(Point)
    // step through all the pixels and confirm I read back the numbers in the test file

    // new db so we can check both the cached and the not-yet cached
    // versions of GetCloudProbability
    det.Update(701);
    {
      const GOESDB& goesDB = det.GetAtmosphere().GetGOESDB();

      // magic numbers to compute pixel id
      const int fgNumEasting = 30; // pixels in easting
      const double fgStartEasting = 553.2; // start of easting coverage
      const double fgStartNorthing = 9736.34; // start of northing coverage
      const double fgDeltaEasting = 2394.9; // size of pixel in easting
      const double fgDeltaNorthing = 5493.03; // size of pixel in northing

      // get site coordinate sytem (usually Malargue)
      CoordinateSystemPtr csDetectorSite = det.GetSiteCoordinateSystem();

      int northPosition = 0;
      for (vector<unsigned int>::iterator i = pixelId.begin(); i != pixelId.end(); ++i) {
        const double eastPosition = *i % fgNumEasting;

        // compute a position x,y somewhere in the current pixel.
        // the fgDeltaEasting/(something) at the end is just a little margin to be sure
        // we are inside the pixel. Note from Figure 2.1 in GAP 2011_091, the east pixels
        // shift by 1 per 6 north pixels.  So for x (below), fgDeltaEasting/6.1 should
        // get us into the pixel without wandering off into the wrong pixel as we move
        // to the north.
        const double x = eastPosition*fgDeltaEasting + fgStartEasting + fgDeltaEasting/6.1;
        const double y = northPosition*fgDeltaNorthing + fgStartNorthing + fgDeltaNorthing/2;
        const Point p(x, y, 0, csDetectorSite);

        const double cloudP4 = goesDB.GetCloudProbability(p);
        const double cloudP3 = goesDB.GetCloudProbability(p);

        CPPUNIT_ASSERT(Verify<Equal>(cloudP4, prob[*i]/gNCloudProbabilityLevels));
        CPPUNIT_ASSERT(Verify<Equal>(cloudP3, prob[*i]/gNCloudProbabilityLevels));

        if (eastPosition == fgNumEasting - 1) {
          northPosition += 1;
        }
      }

      // test maximum cloud probability along line of sight between CLF and XLF
      // first calculate it by stepping between the two positions, then use
      // GOESDB::GetMaximumCloudProbability(Point, Point)
      const Point pXLF(0, 0, 0, CoordinateSystemRegistry::Get("XLF"));
      const Vector clfToXLF = Normalized(pXLF - pCLF);
      const double distance = (pXLF - pCLF).GetMag();
      double currDistance = 0;
      double maxPropBetweenCLFandXLF = 0;
      while (currDistance < distance) {
        const double currProb = goesDB.GetCloudProbability(pCLF + currDistance * clfToXLF);
        if (currProb > maxPropBetweenCLFandXLF)
          maxPropBetweenCLFandXLF = currProb;
        currDistance += 10*m;
      }

      const double cloudProbCLFXLF = goesDB.GetMaximumCloudProbability(pXLF, pCLF);
      CPPUNIT_ASSERT(Verify<Equal>(cloudProbCLFXLF, maxPropBetweenCLFandXLF));

      // order of points shound not matter
      const double cloudProbXLFCLF = goesDB.GetMaximumCloudProbability(pCLF, pXLF);
      CPPUNIT_ASSERT(Verify<Equal>(cloudProbCLFXLF, cloudProbXLFCLF));
    }
  }

  void
  testGOESDBNotFound1()
  {
    Detector& det = Detector::GetInstance();

    det.Update(1);
    const GOESDB& goesdb = det.GetAtmosphere().GetGOESDB();
    goesdb.GetCloudProbability(0);
  }

  void
  testGOESDBNotFound2()
  {
    Detector& det = Detector::GetInstance();

    det.Update(692904871);
    const GOESDB& goesdb = det.GetAtmosphere().GetGOESDB();
    goesdb.GetCloudProbability(12345);
  }

  void
  testGOESDBNotFound3()
  {
    // confirm there is no cloud data over the malargue point
    // (there is no GOES pixel in the database for that point)

    Detector& det = Detector::GetInstance();
    det.Update(700);
    const GOESDB& goesdb = det.GetAtmosphere().GetGOESDB();

    // get site coordinate sytem (usually Malargue)
    CoordinateSystemPtr csDetectorSite = det.GetSiteCoordinateSystem();
    Point pCLF = Point(0, 0, 0, csDetectorSite);

    goesdb.GetCloudProbability(pCLF);
  }

};


CPPUNIT_TEST_SUITE_REGISTRATION(testAtmLowLevelAtmInterface);


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