GammaLib

/* __ Includes ___________________________________________________________ */

#include <iostream>

#include <cmath>

#include "GSkyDir.hpp"

/***********************************************************************//**

 * @brief GSkyDir::dist

 ***************************************************************************/

double GSkyDir::dist(const GSkyDir& dir) const

{

    double cosdis = std::sin(m_dec)*std::sin(dir.m_dec) +

                    std::cos(m_dec)*std::cos(dir.m_dec) *

                    std::cos(dir.m_ra - m_ra);

    double dist   = std::acos(cosdis);

    return dist;

}

/***********************************************************************//**

 * @brief GSkyDir2::dist

 ***************************************************************************/

double GSkyDir2::dist(const GSkyDir2& dir) const

{

    double cosdis = sin_dec() * dir.sin_dec() +

                    cos_dec() * dir.cos_dec() *

                    std::cos(dir.m_ra - m_ra);

    double dist   = std::acos(cosdis);

    return dist;

}

/***********************************************************************//**

 * @brief GSkyDir3::dist

 ***************************************************************************/

double GSkyDir3::dist(const GSkyDir3& dir) const

{

    if (!m_has_cache) {

        m_has_cache = true;

        m_sin_dec   = std::sin(m_dec);

        m_cos_dec   = std::cos(m_dec);

    }

    if (!dir.m_has_cache) {

        dir.m_has_cache = true;

        dir.m_sin_dec   = std::sin(dir.m_dec);

        dir.m_cos_dec   = std::cos(dir.m_dec);

    }

    double cosdis = m_sin_dec * dir.m_sin_dec +

                    m_cos_dec * dir.m_cos_dec *

                    std::cos(dir.m_ra - m_ra);

    double dist   = std::acos(cosdis);

    return dist;

}

/***********************************************************************//**

 * @brief GSkyDir4::dist

 ***************************************************************************/

double GSkyDir4::dist(const GSkyDir4& dir) const

{

    if (m_cache == NULL) {

        m_cache          = new cache;

        m_cache->sin_dec = std::sin(m_dec);

        m_cache->cos_dec = std::cos(m_dec);

    }

    if (dir.m_cache == NULL) {

        dir.m_cache          = new cache;

        dir.m_cache->sin_dec = std::sin(dir.m_dec);

        dir.m_cache->cos_dec = std::cos(dir.m_dec);

    }

    double cosdis = m_cache->sin_dec * dir.m_cache->sin_dec +

                    m_cache->cos_dec * dir.m_cache->cos_dec *

                    std::cos(dir.m_ra - m_ra);

    double dist   = std::acos(cosdis);

    return dist;

}

/***********************************************************************//**

 * @brief test_GSkyDir

 ***************************************************************************/

void test_GSkyDir(int nouter, int ninner)

{

    // Start timer

    clock_t t_start = clock();

    // Loop

    double sum = 0.0;

    for (int outer = 0; outer < nouter; ++outer) {

        GSkyDir crab( 83.6331,  22.0145);

        GSkyDir vela(128.8361, -45.1764);

        for (int inner = 0; inner < ninner; ++inner) {

            double dist = crab.dist(vela);

            sum += dist;

        }

    }

    // Stop timer

    double t_elapse = (double)(clock() - t_start) / (double)CLOCKS_PER_SEC;

    std::cout << "  test_GSkyDir1: CPU usage = " << t_elapse << " sec";

    std::cout << " sum=" << sum << std::endl;

    // Return

    return;

}

/***********************************************************************//**

 * @brief test_GSkyDir2

 ***************************************************************************/

void test_GSkyDir2(int nouter, int ninner)

{

    // Start timer

    clock_t t_start = clock();

    // Loop

    double sum = 0.0;

    for (int outer = 0; outer < nouter; ++outer) {

        GSkyDir2 crab( 83.6331,  22.0145);

        GSkyDir2 vela(128.8361, -45.1764);

        for (int inner = 0; inner < ninner; ++inner) {

            double dist = crab.dist(vela);

            sum += dist;

        }

    }

    // Stop timer

    double t_elapse = (double)(clock() - t_start) / (double)CLOCKS_PER_SEC;

    std::cout << "  test_GSkyDir2: CPU usage = " << t_elapse << " sec";

    std::cout << " sum=" << sum << std::endl;

    // Return

    return;

}

/***********************************************************************//**

 * @brief test_GSkyDir3

 ***************************************************************************/

void test_GSkyDir3(int nouter, int ninner)

{

    // Start timer

    clock_t t_start = clock();

    // Loop

    double sum = 0.0;

    for (int outer = 0; outer < nouter; ++outer) {

        GSkyDir3 crab( 83.6331,  22.0145);

        GSkyDir3 vela(128.8361, -45.1764);

        for (int inner = 0; inner < ninner; ++inner) {

            double dist = crab.dist(vela);

            sum += dist;

        }

    }

    // Stop timer

    double t_elapse = (double)(clock() - t_start) / (double)CLOCKS_PER_SEC;

    std::cout << "  test_GSkyDir3: CPU usage = " << t_elapse << " sec";

    std::cout << " sum=" << sum << std::endl;

    // Return

    return;

}

/***********************************************************************//**

 * @brief test_GSkyDir4

 ***************************************************************************/

void test_GSkyDir4(int nouter, int ninner)

{

    // Start timer

    clock_t t_start = clock();

    // Loop

    double sum = 0.0;

    for (int outer = 0; outer < nouter; ++outer) {

        GSkyDir4 crab( 83.6331,  22.0145);

        GSkyDir4 vela(128.8361, -45.1764);

        for (int inner = 0; inner < ninner; ++inner) {

            double dist = crab.dist(vela);

            sum += dist;

        }

    }

    // Stop timer

    double t_elapse = (double)(clock() - t_start) / (double)CLOCKS_PER_SEC;

    std::cout << "  test_GSkyDir4: CPU usage = " << t_elapse << " sec";

    std::cout << " sum=" << sum << std::endl;

    // Return

    return;

}

/***********************************************************************//**

 * @brief Main

 ***************************************************************************/

int main(void)

{

    // Test 10000000 outer, 1 inner

    std::cout << "Test 10000000 outer, 1 inner" << std::endl;

    test_GSkyDir(10000000, 1);

    test_GSkyDir2(10000000, 1);

    test_GSkyDir3(10000000, 1);

    test_GSkyDir4(10000000, 1);

    // Test 5000000 outer, 2 inner

    std::cout << "Test 5000000 outer, 2 inner" << std::endl;

    test_GSkyDir(5000000, 2);

    test_GSkyDir2(5000000, 2);

    test_GSkyDir3(5000000, 2);

    test_GSkyDir4(5000000, 2);

    // Test 3333333 outer, 3 inner

    std::cout << "Test 3333333 outer, 3 inner" << std::endl;

    test_GSkyDir(3333333, 3);

    test_GSkyDir2(3333333, 3);

    test_GSkyDir3(3333333, 3);

    test_GSkyDir4(3333333, 3);

    // Test 1000000 outer, 10 inner

    std::cout << "Test 1000000 outer, 10 inner" << std::endl;

    test_GSkyDir(1000000, 10);

    test_GSkyDir2(1000000, 10);

    test_GSkyDir3(1000000, 10);

    test_GSkyDir4(1000000, 10);

    // Test 100000 outer, 100 inner

    std::cout << "Test 100000 outer, 100 inner" << std::endl;

    test_GSkyDir(100000, 100);

    test_GSkyDir2(100000, 100);

    test_GSkyDir3(100000, 100);

    test_GSkyDir4(100000, 100);

    // Test 10000 outer, 1000 inner

    std::cout << "Test 10000 outer, 1000 inner" << std::endl;

    test_GSkyDir(10000, 1000);

    test_GSkyDir2(10000, 1000);

    test_GSkyDir3(10000, 1000);

    test_GSkyDir4(10000, 1000);

    // Exit

    return 0;

}