code:
/*
 [auto_generated]
 boost/numeric/odeint/stepper/implicit_euler.hpp

 [begin_description]
 Impementation of the implicit Euler method. Works with ublas::vector as state type.
 [end_description]

 Copyright 2010-2012 Mario Mulansky
 Copyright 2010-2012 Karsten Ahnert
 Copyright 2012 Christoph Koke

 Distributed under the Boost Software License, Version 1.0.
 (See accompanying file LICENSE_1_0.txt or
 copy at http://www.boost.org/LICENSE_1_0.txt)
 */

#ifndef BOOST_NUMERIC_ODEINT_STEPPER_IMPLICIT_EULER_HPP_INCLUDED
#define BOOST_NUMERIC_ODEINT_STEPPER_IMPLICIT_EULER_HPP_INCLUDED

#include <utility>

#include <boost/numeric/odeint/util/bind.hpp>
#include <boost/numeric/odeint/util/unwrap_reference.hpp>
#include <boost/numeric/odeint/stepper/stepper_categories.hpp>

#include <boost/numeric/odeint/util/ublas_wrapper.hpp>
#include <boost/numeric/odeint/util/is_resizeable.hpp>
#include <boost/numeric/odeint/util/resizer.hpp>

#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/lu.hpp>

namespace boost {
namespace numeric {
namespace odeint {

template< class ValueType , class Resizer = initially_resizer >
class implicit_euler
{

public:

    typedef ValueType value_type;
    typedef value_type time_type;
    typedef boost::numeric::ublas::vector< value_type > state_type;
    typedef state_wrapper< state_type > wrapped_state_type;
    typedef state_type deriv_type;
    typedef state_wrapper< deriv_type > wrapped_deriv_type;
    typedef boost::numeric::ublas::matrix< value_type > matrix_type;
    typedef state_wrapper< matrix_type > wrapped_matrix_type;
    typedef boost::numeric::ublas::permutation_matrix< size_t > pmatrix_type;
    typedef state_wrapper< pmatrix_type > wrapped_pmatrix_type;
    typedef Resizer resizer_type;
    typedef stepper_tag stepper_category;
    typedef implicit_euler< ValueType , Resizer > stepper_type;

    implicit_euler( value_type epsilon = 1E-6 )
    : m_epsilon( epsilon )
    { }

    template< class System >
    void do_step( System system , state_type &x , time_type t , time_type dt )
    {
        typedef typename odeint::unwrap_reference< System >::type system_type;
        typedef typename odeint::unwrap_reference< typename system_type::first_type >::type deriv_func_type;
        typedef typename odeint::unwrap_reference< typename system_type::second_type >::type jacobi_func_type;
        system_type &sys = system;
        deriv_func_type &deriv_func = sys.first;
        jacobi_func_type &jacobi_func = sys.second;

        m_resizer.adjust_size( x , detail::bind( &stepper_type::template resize_impl<state_type> , detail::ref( *this ) , detail::_1 ) );

        for( size_t i=0 ; i<x.size() ; ++i )
            m_pm.m_v[i] = i;

        t += dt;

        // apply first Newton step
        deriv_func( x , m_dxdt.m_v , t );

        m_b.m_v = dt * m_dxdt.m_v;

        jacobi_func( x , m_jacobi.m_v  , t );
        m_jacobi.m_v *= dt;
        m_jacobi.m_v -= boost::numeric::ublas::identity_matrix< value_type >( x.size() );

        solve( m_b.m_v , m_jacobi.m_v );

        m_x.m_v = x - m_b.m_v;

        // iterate Newton until some precision is reached
        // ToDo: maybe we should apply only one Newton step -> linear implicit one-step scheme
        while( boost::numeric::ublas::norm_2( m_b.m_v ) > m_epsilon )
        {
            deriv_func( m_x.m_v , m_dxdt.m_v , t );
            m_b.m_v = x - m_x.m_v + dt*m_dxdt.m_v;

            // simplified version, only the first Jacobian is used
            //            jacobi( m_x , m_jacobi , t );
            //            m_jacobi *= dt;
            //            m_jacobi -= boost::numeric::ublas::identity_matrix< value_type >( x.size() );

            solve( m_b.m_v , m_jacobi.m_v );

            m_x.m_v -= m_b.m_v;
        }
        x = m_x.m_v;
    }

    template< class StateType >
    void adjust_size( const StateType &x )
    {
        resize_impl( x );
    }

private:

    template< class StateIn >
    bool resize_impl( const StateIn &x )
    {
        bool resized = false;
        resized |= adjust_size_by_resizeability( m_dxdt , x , typename is_resizeable<deriv_type>::type() );
        resized |= adjust_size_by_resizeability( m_x , x , typename is_resizeable<state_type>::type() );
        resized |= adjust_size_by_resizeability( m_b , x , typename is_resizeable<deriv_type>::type() );
        resized |= adjust_size_by_resizeability( m_jacobi , x , typename is_resizeable<matrix_type>::type() );
        resized |= adjust_size_by_resizeability( m_pm , x , typename is_resizeable<pmatrix_type>::type() );
        return resized;
    }

    void solve( state_type &x , matrix_type &m )
    {
        int res = boost::numeric::ublas::lu_factorize( m , m_pm.m_v );
        if( res != 0 ) std::exit(0);
        boost::numeric::ublas::lu_substitute( m , m_pm.m_v , x );
    }

private:

    value_type m_epsilon;
    resizer_type m_resizer;
    wrapped_deriv_type m_dxdt;
    wrapped_state_type m_x;
    wrapped_deriv_type m_b;
    wrapped_matrix_type m_jacobi;
    wrapped_pmatrix_type m_pm;

};

} // odeint
} // numeric
} // boost

#endif // BOOST_NUMERIC_ODEINT_STEPPER_IMPLICIT_EULER_HPP_INCLUDED


code:
// (C) Copyright Tobias Schwinger
//
// Use modification and distribution are subject to the boost Software License,
// Version 1.0. (See http://www.boost.org/LICENSE_1_0.txt).

//------------------------------------------------------------------------------

// no include guards, this file is intended for multiple inclusions

// this file has been generated from the master.hpp file in the same directory
# define BOOST_FT_cc_id 1
# define BOOST_FT_cc_name implicit_cc
# define BOOST_FT_cc BOOST_PP_EMPTY
# define BOOST_FT_cond BOOST_FT_CC_IMPLICIT
# if BOOST_FT_cond
# define BOOST_FT_config_valid 1
# include BOOST_FT_cc_file
# endif
# undef BOOST_FT_cond
# undef BOOST_FT_cc_name
# undef BOOST_FT_cc
# undef BOOST_FT_cc_id
# define BOOST_FT_cc_id 2
# define BOOST_FT_cc_name cdecl_cc
# define BOOST_FT_cc BOOST_PP_IDENTITY(__cdecl )
# define BOOST_FT_cond BOOST_FT_CC_CDECL
# if BOOST_FT_cond
# define BOOST_FT_config_valid 1
# include BOOST_FT_cc_file
# endif
# undef BOOST_FT_cond
# undef BOOST_FT_cc_name
# undef BOOST_FT_cc
# undef BOOST_FT_cc_id
# define BOOST_FT_cc_id 3
# define BOOST_FT_cc_name stdcall_cc
# define BOOST_FT_cc BOOST_PP_IDENTITY(__stdcall )
# define BOOST_FT_cond BOOST_FT_CC_STDCALL
# if BOOST_FT_cond
# define BOOST_FT_config_valid 1
# include BOOST_FT_cc_file
# endif
# undef BOOST_FT_cond
# undef BOOST_FT_cc_name
# undef BOOST_FT_cc
# undef BOOST_FT_cc_id
# define BOOST_FT_cc_id 4
# define BOOST_FT_cc_name pascal_cc
# define BOOST_FT_cc BOOST_PP_IDENTITY(pascal )
# define BOOST_FT_cond BOOST_FT_CC_PASCAL
# if BOOST_FT_cond
# define BOOST_FT_config_valid 1
# include BOOST_FT_cc_file
# endif
# undef BOOST_FT_cond
# undef BOOST_FT_cc_name
# undef BOOST_FT_cc
# undef BOOST_FT_cc_id
# define BOOST_FT_cc_id 5
# define BOOST_FT_cc_name fastcall_cc
# define BOOST_FT_cc BOOST_PP_IDENTITY(__fastcall)
# define BOOST_FT_cond BOOST_FT_CC_FASTCALL
# if BOOST_FT_cond
# define BOOST_FT_config_valid 1
# include BOOST_FT_cc_file
# endif
# undef BOOST_FT_cond
# undef BOOST_FT_cc_name
# undef BOOST_FT_cc
# undef BOOST_FT_cc_id
# define BOOST_FT_cc_id 6
# define BOOST_FT_cc_name clrcall_cc
# define BOOST_FT_cc BOOST_PP_IDENTITY(__clrcall )
# define BOOST_FT_cond BOOST_FT_CC_CLRCALL
# if BOOST_FT_cond
# define BOOST_FT_config_valid 1
# include BOOST_FT_cc_file
# endif
# undef BOOST_FT_cond
# undef BOOST_FT_cc_name
# undef BOOST_FT_cc
# undef BOOST_FT_cc_id
# define BOOST_FT_cc_id 7
# define BOOST_FT_cc_name thiscall_cc
# define BOOST_FT_cc BOOST_PP_IDENTITY(__thiscall)
# define BOOST_FT_cond BOOST_FT_CC_THISCALL
# if BOOST_FT_cond
# define BOOST_FT_config_valid 1
# include BOOST_FT_cc_file
# endif
# undef BOOST_FT_cond
# undef BOOST_FT_cc_name
# undef BOOST_FT_cc
# undef BOOST_FT_cc_id
# define BOOST_FT_cc_id 8
# define BOOST_FT_cc_name thiscall_cc
# define BOOST_FT_cc BOOST_PP_EMPTY
# define BOOST_FT_cond BOOST_FT_CC_IMPLICIT_THISCALL
# if BOOST_FT_cond
# define BOOST_FT_config_valid 1
# include BOOST_FT_cc_file
# endif
# undef BOOST_FT_cond
# undef BOOST_FT_cc_name
# undef BOOST_FT_cc
# undef BOOST_FT_cc_id
# ifndef BOOST_FT_config_valid
# define BOOST_FT_cc_id 1
# define BOOST_FT_cc_name implicit_cc
# define BOOST_FT_cc BOOST_PP_EMPTY
# define BOOST_FT_cond 0x00000001
# include BOOST_FT_cc_file
# undef BOOST_FT_cond
# undef BOOST_FT_cc_name
# undef BOOST_FT_cc
# undef BOOST_FT_cc_id
# else
# undef BOOST_FT_config_valid
# endif

license: mit


code:
// (C) Copyright Tobias Schwinger
//
// Use modification and distribution are subject to the boost Software License,
// Version 1.0. (See http://www.boost.org/LICENSE_1_0.txt).

//------------------------------------------------------------------------------

#ifndef BOOST_FT_DETAIL_TO_SEQUENCE_HPP_INCLUDED
#define BOOST_FT_DETAIL_TO_SEQUENCE_HPP_INCLUDED

#include <boost/mpl/eval_if.hpp>
#include <boost/mpl/identity.hpp>
#include <boost/mpl/is_sequence.hpp>
#include <boost/mpl/placeholders.hpp>
#include <boost/type_traits/add_reference.hpp>

#include <boost/function_types/is_callable_builtin.hpp>

namespace boost { namespace function_types { namespace detail {

// wrap first arguments in components, if callable builtin type
template<typename T>
struct to_sequence
{
  typedef typename
   mpl::eval_if
   < is_callable_builtin<T>
   , to_sequence< components<T> >
   , mpl::identity< T >
   >::type
  type;
};

// reduce template instantiations, if possible
template<typename T, typename U>
struct to_sequence< components<T,U> >
{
  typedef typename components<T,U>::types type;
};

} } } // namespace ::boost::function_types::detail

#endif


code:
#include "FileNode.hpp"

#include <QVariant>
#include <QUrl>
#include <QDir>
#include <QFileInfo>

#include "Directory.hpp"
#include "MessageBox.hpp"

FileNode::FileNode (QTreeWidgetItem * parent
                    , QNetworkAccessManager * network_manager
                    , QString const& local_file_path
                    , QUrl const& url
                    , bool http_only)
  : QTreeWidgetItem {parent, Type}
  , remote_file_ {this, network_manager, local_file_path, http_only}
  , block_sync_ {false}
{
  sync_blocker b {this};
  setFlags (flags () | Qt::ItemIsUserCheckable);
  setText (0, QFileInfo {local_file_path}.fileName ()); // display
  setData (0, Qt::UserRole, url);
  setData (0, Qt::UserRole + 1, local_file_path); // local absolute path
  setCheckState (0, Qt::Unchecked);
}

void FileNode::error (QString const& title, QString const& message)
{
  MessageBox::warning_message (treeWidget (), title, message);
}

bool FileNode::sync (bool local)
{
  if (block_sync_)
    {
      return true;
    }
  return remote_file_.sync (data (0, Qt::UserRole).toUrl (), local);
}

void FileNode::download_progress (qint64 bytes_received, qint64 total_bytes)
{
  sync_blocker b {this};
  setData (1, Qt::UserRole, total_bytes);
  if (bytes_received < 0)
    {
      setData (1, Qt::DisplayRole, 0ll);
      setCheckState (0, Qt::Unchecked);
    }
  else
    {
      setData (1, Qt::DisplayRole, bytes_received);
    }
  static_cast<Directory *> (treeWidget ())->update (parent ());
}

void FileNode::download_finished (bool success)
{
  sync_blocker b {this};
  if (!success)
    {
      setData (1, Qt::UserRole, 0ll);
      setData (1, Qt::DisplayRole, 0ll);
    }
  setCheckState (0, success ? Qt::Checked : Qt::Unchecked);
  static_cast<Directory *> (treeWidget ())->update (parent ());
}

void FileNode::abort ()
{
  sync_blocker b {this};
  remote_file_.abort ();
}


code:
/*==============================================================================
    Copyright (c) 2010-2011 Bryce Lelbach

    Distributed under the Boost Software License, Version 1.0. (See accompanying
    file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
==============================================================================*/

#ifndef BOOST_DETAIL_SORTED_HPP
#define BOOST_DETAIL_SORTED_HPP

#include <boost/detail/iterator.hpp>

#include <functional>

namespace boost {
namespace detail {

template<class Iterator, class Comp>
inline Iterator is_sorted_until (Iterator first, Iterator last, Comp c) {
  if (first == last)
    return last;

  Iterator it = first; ++it;

  for (; it != last; first = it, ++it)
    if (c(*it, *first))
      return it;

  return it;
}

template<class Iterator>
inline Iterator is_sorted_until (Iterator first, Iterator last) {
  typedef typename boost::detail::iterator_traits<Iterator>::value_type
    value_type;

  typedef std::less<value_type> c;

  return ::boost::detail::is_sorted_until(first, last, c());
}

template<class Iterator, class Comp>
inline bool is_sorted (Iterator first, Iterator last, Comp c) {
  return ::boost::detail::is_sorted_until(first, last, c) == last;
}

template<class Iterator>
inline bool is_sorted (Iterator first, Iterator last) {
  return ::boost::detail::is_sorted_until(first, last) == last;
}

} // detail
} // boost

#endif // BOOST_DETAIL_SORTED_HPP


code:
// Boost.Units - A C++ library for zero-overhead dimensional analysis and
// unit/quantity manipulation and conversion
//
// Copyright (C) 2003-2008 Matthias Christian Schabel
// Copyright (C) 2008 Steven Watanabe
//
// Distributed under the Boost Software License, Version 1.0. (See
// accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)

#ifndef BOOST_UNITS_ACCELERATION_DERIVED_DIMENSION_HPP
#define BOOST_UNITS_ACCELERATION_DERIVED_DIMENSION_HPP

#include <boost/units/derived_dimension.hpp>
#include <boost/units/physical_dimensions/length.hpp>
#include <boost/units/physical_dimensions/time.hpp>

namespace boost {

namespace units {

/// derived dimension for acceleration : L T^-2
typedef derived_dimension<length_base_dimension,1,
                          time_base_dimension,-2>::type acceleration_dimension;

} // namespace units

} // namespace boost

#endif // BOOST_UNITS_ACCELERATION_DERIVED_DIMENSION_HPP


code:
// (C) Copyright Steve Cleary, Beman Dawes, Howard Hinnant & John Maddock 2000.
//  Use, modification and distribution are subject to the Boost Software License,
//  Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
//  http://www.boost.org/LICENSE_1_0.txt).
//
//  See http://www.boost.org/libs/type_traits for most recent version including documentation.

#ifndef BOOST_TT_ADD_POINTER_HPP_INCLUDED
#define BOOST_TT_ADD_POINTER_HPP_INCLUDED

#include <boost/type_traits/remove_reference.hpp>

namespace boost {

#if defined(__BORLANDC__) && (__BORLANDC__ < 0x5A0)
//
// For some reason this implementation stops Borlands compiler
// from dropping cv-qualifiers, it still fails with references
// to arrays for some reason though (shrug...) (JM 20021104)
//
template <typename T>
struct add_pointer
{
    typedef T* type;
};
template <typename T>
struct add_pointer<T&>
{
    typedef T* type;
};
template <typename T>
struct add_pointer<T&const>
{
    typedef T* type;
};
template <typename T>
struct add_pointer<T&volatile>
{
    typedef T* type;
};
template <typename T>
struct add_pointer<T&const volatile>
{
    typedef T* type;
};

#else

template <typename T>
struct add_pointer
{
    typedef typename remove_reference<T>::type no_ref_type;
    typedef no_ref_type* type;
};

#endif

} // namespace boost

#endif // BOOST_TT_ADD_POINTER_HPP_INCLUDED


code:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define RADS 0.0174532925199433
#define DEGS 57.2957795130823
#define TPI 6.28318530717959
#define PI 3.1415927

/* ratio of     earth radius to astronomical unit */
#define ER_OVER_AU 0.0000426352325194252

/* all prototypes here */

double getcoord(int coord);
void getargs(int argc, char *argv[], int *y, int *m, double *tz, double *glong, double *glat);
double range(double y);
double rangerad(double y);
double days(int y, int m, int dn, double hour);
double days_(int *y, int *m, int *dn, double *hour);
void moonpos(double, double *, double *, double *);
void sunpos(double , double *, double *, double *);
double moontransit(int y, int m, int d, double timezone, double glat, double glong, int *nt);
double atan22(double y, double x);
double epsilon(double d);
void equatorial(double d, double *lon, double *lat, double *r);
void ecliptic(double d, double *lon, double *lat, double *r);
double gst(double d);
void topo(double lst, double glat, double *alp, double *dec, double *r);
double alt(double glat, double ha, double dec);
void libration(double day, double lambda, double beta, double alpha, double *l, double *b, double *p);
void illumination(double day, double lra, double ldec, double dr, double sra, double sdec, double *pabl, double *ill);
int daysinmonth(int y, int m);
int isleap(int y);
void tmoonsub_(double *day, double *glat, double *glong, double *moonalt,
   double *mrv, double *l, double *b, double *paxis);

static const char
usage[] = "  Usage: tmoon date[yyyymm] timz[+/-h.hh] long[+/-dddmm] lat[+/-ddmm]\n"
            "example: tmoon 200009 0 -00155 5230\n";

/*
  getargs() gets the arguments from the command line, does some basic error
  checking, and converts arguments into numerical form. Arguments are passed
  back in pointers. Error messages print to stderr so re-direction of output
  to file won't leave users blind. Error checking prints list of all errors
  in a command line before quitting.
*/
void getargs(int argc, char *argv[], int *y, int *m, double *tz,
             double *glong, double *glat) {

  int date, latitude, longitude;
  int mflag = 0, yflag = 0, longflag = 0, latflag = 0, tzflag = 0;
  int longminflag = 0, latminflag = 0, dflag = 0;

  /* if not right number of arguments, then print example command line */

  if (argc !=5) {
    fprintf(stderr, usage);
    exit(EXIT_FAILURE);
  }

  date = atoi(argv[1]);
  *y = date / 100;
  *m = date - *y * 100;
  *tz = (double) atof(argv[2]);
  longitude = atoi(argv[3]);
  latitude = atoi(argv[4]);
  *glong = RADS * getcoord(longitude);
  *glat = RADS * getcoord(latitude);

  /* set a flag for each error found */

  if (*m > 12 || *m < 1) mflag = 1;
  if (*y > 2500) yflag = 1;
  if (date < 150001) dflag = 1;
  if (fabs((float) *glong) > 180 * RADS) longflag = 1;
  if (abs(longitude) % 100 > 59) longminflag = 1;
  if (fabs((float) *glat) > 90 * RADS) latflag = 1;
  if (abs(latitude) % 100 > 59) latminflag = 1;
  if (fabs((float) *tz) > 12) tzflag = 1;

  /* print all the errors found */

  if (dflag == 1) {
    fprintf(stderr, "date: dates must be in form yyyymm, gregorian, and later than 1500 AD\n");
  }
  if (yflag == 1) {
    fprintf(stderr, "date: too far in future - accurate from 1500 to 2500\n");
  }
  if (mflag == 1) {
    fprintf(stderr, "date: month must be in range 0 to 12, eg - August 2000 is entered as 200008\n");
  }
  if (tzflag == 1) {
    fprintf(stderr, "timz: must be in range +/- 12 hours, eg -6 for Chicago\n");
  }
  if (longflag == 1) {
    fprintf(stderr, "long: must be in range +/- 180 degrees\n");
  }
  if (longminflag == 1) {
    fprintf(stderr, "long: last two digits are arcmin - max 59\n");
  }
  if (latflag == 1) {
    fprintf(stderr, " lat: must be in range +/- 90 degrees\n");
  }
  if (latminflag == 1) {
    fprintf(stderr, " lat: last two digits are arcmin - max 59\n");
  }

  /* quits if one or more flags set */

  if (dflag + mflag + yflag + longflag + latflag + tzflag + longminflag + latminflag > 0) {
    exit(EXIT_FAILURE);
  }

}

/*
   returns coordinates in decimal degrees given the
   coord as a ddmm value stored in an integer.
*/
double getcoord(int coord) {
  int west = 1;
  double glg, deg;
  if (coord < 0) west = -1;
  glg = fabs((double) coord/100);
  deg = floor(glg);
  glg = west* (deg + (glg - deg)*100 / 60);
  return(glg);
}

/*
  days() takes the year, month, day in the month and decimal hours
  in the day and returns the number of days since J2000.0.
  Assumes Gregorian calendar.
*/
double days(int y, int m, int d, double h) {
  int a, b;
  double day;

  /*
    The lines below work from 1900 march to feb 2100
    a = 367 * y - 7 * (y + (m + 9) / 12) / 4 + 275 * m / 9 + d;
    day = (double)a - 730531.5 + hour / 24;
  */

  /*  These lines work for any Gregorian date since 0 AD */
  if (m ==1 || m==2) {
    m +=12;
    y -= 1;
  }
  a = y / 100;
  b = 2 - a + a/4;
  day = floor(365.25*(y + 4716)) + floor(30.6001*(m + 1))
    + d + b - 1524.5 - 2451545 + h/24;
  return(day);
}
double days_(int *y0, int *m0, int *d0, double *h0)
{
  return days(*y0,*m0,*d0,*h0);
}

/*
Returns 1 if y a leap year, and 0 otherwise, according
to the Gregorian calendar
*/
int isleap(int y) {
  int a = 0;
  if(y % 4 == 0) a = 1;
  if(y % 100 == 0) a = 0;
  if(y % 400 == 0) a = 1;
  return(a);
}

/*
Given the year and the month, function returns the
number of days in the month. Valid for Gregorian
calendar.
*/
int daysinmonth(int y, int m) {
  int b = 31;
  if(m == 2) {
    if(isleap(y) == 1) b= 29;
    else b = 28;
  }
  if(m == 4 || m == 6 || m == 9 || m == 11) b = 30;
  return(b);
}

/*
moonpos() takes days from J2000.0 and returns ecliptic coordinates
of moon in the pointers. Note call by reference.
This function is within a couple of arcminutes most of the time,
and is truncated from the Meeus Ch45 series, themselves truncations of
ELP-2000. Returns moon distance in earth radii.
Terms have been written out explicitly rather than using the
table based method as only a small number of terms is
retained.
*/
void moonpos(double d, double *lambda, double *beta, double *rvec) {
  double dl, dB, dR, L, D, M, M1, F, e, lm, bm, rm, t;

  t = d / 36525;

  L = range(218.3164591  + 481267.88134236  * t) * RADS;
  D = range(297.8502042  + 445267.1115168  * t) * RADS;
  M = range(357.5291092  + 35999.0502909  * t) * RADS;
  M1 = range(134.9634114  + 477198.8676313  * t - .008997 * t * t) * RADS;
  F = range(93.27209929999999  + 483202.0175273  * t - .0034029*t*t)*RADS;
  e = 1 - .002516 * t;

  dl =      6288774 * sin(M1);
  dl +=     1274027 * sin(2 * D - M1);
  dl +=      658314 * sin(2 * D);
  dl +=      213618 * sin(2 * M1);
  dl -=  e * 185116 * sin(M);
  dl -=      114332 * sin(2 * F) ;
  dl +=       58793 * sin(2 * D - 2 * M1);
  dl +=   e * 57066 * sin(2 * D - M - M1) ;
  dl +=       53322 * sin(2 * D + M1);
  dl +=   e * 45758 * sin(2 * D - M);
  dl -=   e * 40923 * sin(M - M1);
  dl -=       34720 * sin(D) ;
  dl -=   e * 30383 * sin(M + M1) ;
  dl +=       15327 * sin(2 * D - 2 * F) ;
  dl -=       12528 * sin(M1 + 2 * F);
  dl +=       10980 * sin(M1 - 2 * F);
  lm = rangerad(L + dl / 1000000 * RADS);

  dB =   5128122 * sin(F);
  dB +=   280602 * sin(M1 + F);
  dB +=   277693 * sin(M1 - F);
  dB +=   173237 * sin(2 * D - F);
  dB +=    55413 * sin(2 * D - M1 + F);
  dB +=    46271 * sin(2 * D - M1 - F);
  dB +=    32573 * sin(2 * D + F);
  dB +=    17198 * sin(2 * M1 + F);
  dB +=     9266 * sin(2 * D + M1 - F);
  dB +=     8822 * sin(2 * M1 - F);
  dB += e * 8216 * sin(2 * D - M - F);
  dB +=     4324 * sin(2 * D - 2 * M1 - F);
  bm = dB / 1000000 * RADS;

  dR =    -20905355 * cos(M1);
  dR -=     3699111 * cos(2 * D - M1);
  dR -=     2955968 * cos(2 * D);
  dR -=      569925 * cos(2 * M1);
  dR +=   e * 48888 * cos(M);
  dR -=        3149 * cos(2 * F);
  dR +=      246158 * cos(2 * D - 2 * M1);
  dR -=  e * 152138 * cos(2 * D - M - M1) ;
  dR -=      170733 * cos(2 * D + M1);
  dR -=  e * 204586 * cos(2 * D - M);
  dR -=  e * 129620 * cos(M - M1);
  dR +=      108743 * cos(D);
  dR +=  e * 104755 * cos(M + M1);
  dR +=       79661 * cos(M1 - 2 * F);
  rm = 385000.56  + dR / 1000;

  *lambda = lm;
  *beta = bm;
  /* distance to Moon must be in Earth radii */
  *rvec = rm / 6378.14;
}

/*
topomoon() takes the local siderial time, the geographical
latitude of the observer, and pointers to the geocentric
equatorial coordinates. The function overwrites the geocentric
coordinates with topocentric coordinates on a simple spherical
earth model (no polar flattening). Expects Moon-Earth distance in
Earth radii.    Formulas scavenged from Astronomical Almanac 'low
precision formulae for Moon position' page D46.
*/

void topo(double lst, double glat, double *alp, double *dec, double *r) {
  double x, y, z, r1;
  x = *r * cos(*dec) * cos(*alp) - cos(glat) * cos(lst);
  y = *r * cos(*dec) * sin(*alp) - cos(glat) * sin(lst);
  z = *r * sin(*dec)  - sin(glat);
  r1 = sqrt(x*x + y*y + z*z);
  *alp = atan22(y, x);
  *dec = asin(z / r1);
  *r = r1;
}

/*
moontransit() takes date, the time zone and geographic longitude
of observer and returns the time (decimal hours) of lunar transit
on that day if there is one, and sets the notransit flag if there
isn't. See Explanatory Supplement to Astronomical Almanac
section 9.32 and 9.31 for the method.
*/

double moontransit(int y, int m, int d, double tz, double glat, double glong, int *notransit) {
  double hm, ht, ht1, lon, lat, rv, dnew, lst;
  int itcount;

  ht1 = 180 * RADS;
  ht = 0;
  itcount = 0;
  *notransit = 0;
  do {
    ht = ht1;
    itcount++;
    dnew = days(y, m, d, ht * DEGS/15) - tz/24;
    lst = gst(dnew) + glong;
    /* find the topocentric Moon ra (hence hour angle) and dec */
    moonpos(dnew, &lon, &lat, &rv);
    equatorial(dnew, &lon, &lat, &rv);
    topo(lst, glat, &lon, &lat, &rv);
    hm = rangerad(lst -  lon);
    ht1 = rangerad(ht - hm);
    /* if no convergence, then no transit on that day */
    if (itcount > 30) {
      *notransit = 1;
      break;
    }
  }
  while (fabs(ht - ht1) > 0.04 * RADS);
  return(ht1);
}

/*
  Calculates the selenographic coordinates of either the sub Earth point
  (optical libration) or the sub-solar point (selen. coords of centre of
  bright hemisphere).  Based on Meeus chapter 51 but neglects physical
  libration and nutation, with some simplification of the formulas.
*/
void libration(double day, double lambda, double beta, double alpha, double *l, double *b, double *p) {
  double i, f, omega, w, y, x, a, t, eps;
  t = day / 36525;
  i = 1.54242 * RADS;
  eps = epsilon(day);
  f = range(93.2720993 + 483202.0175273 * t - .0034029 * t * t) * RADS;
  omega = range(125.044555 - 1934.1361849 * t + .0020762 * t * t) * RADS;
  w = lambda - omega;
  y = sin(w) * cos(beta) * cos(i) - sin(beta) * sin(i);
  x = cos(w) * cos(beta);
  a = atan22(y, x);
  *l = a - f;

  /*  kludge to catch cases of 'round the back' angles  */
  if (*l < -90 * RADS) *l += TPI;
  if (*l > 90 * RADS)  *l -= TPI;
  *b = asin(-sin(w) * cos(beta) * sin(i) - sin(beta) * cos(i));

  /*  pa pole axis - not used for Sun stuff */
  x = sin(i) * sin(omega);
  y = sin(i) * cos(omega) * cos(eps) - cos(i) * sin(eps);
  w = atan22(x, y);
  *p = rangerad(asin(sqrt(x*x + y*y) * cos(alpha - w) / cos(*b)));
}

/*
  Takes: days since J2000.0, eq coords Moon, ratio of moon to sun distance,
  eq coords Sun
  Returns: position angle of bright limb wrt NCP, percentage illumination
  of Sun
*/
void illumination(double day , double lra, double ldec, double dr, double sra, double sdec, double *pabl, double *ill) {
  double x, y, phi, i;
  (void)day;
  y = cos(sdec) * sin(sra - lra);
  x = sin(sdec) * cos(ldec) - cos(sdec) * sin(ldec) * cos (sra - lra);
  *pabl = atan22(y, x);
  phi = acos(sin(sdec) * sin(ldec) + cos(sdec) * cos(ldec) * cos(sra-lra));
  i = atan22(sin(phi) , (dr - cos(phi)));
  *ill = 0.5*(1 + cos(i));
}

/*
sunpos() takes days from J2000.0 and returns ecliptic longitude
of Sun in the pointers. Latitude is zero at this level of precision,
but pointer left in for consistency in number of arguments.
This function is within 0.01 degree (1 arcmin) almost all the time
for a century either side of J2000.0. This is from the 'low precision
fomulas for the Sun' from C24 of Astronomical Alamanac
*/
void sunpos(double d, double *lambda, double *beta, double *rvec) {
  double L, g, ls, bs, rs;

  L = range(280.461 + .9856474 * d) * RADS;
  g = range(357.528 + .9856003 * d) * RADS;
  ls = L + (1.915 * sin(g) + .02 * sin(2 * g)) * RADS;
  bs = 0;
  rs = 1.00014 - .01671 * cos(g) - .00014 * cos(2 * g);
  *lambda = ls;
  *beta = bs;
  *rvec = rs;
}

/*
this routine returns the altitude given the days since J2000.0
the hour angle and declination of the object and the latitude
of the observer. Used to find the Sun's altitude to put a letter
code on the transit time, and to find the Moon's altitude at
transit just to make sure that the Moon is visible.
*/
double alt(double glat, double ha, double dec) {
  return(asin(sin(dec) * sin(glat) + cos(dec) * cos(glat) * cos(ha)));
}

/* returns an angle in degrees in the range 0 to 360 */
double range(double x) {
  double a, b;
  b = x / 360;
  a = 360 * (b - floor(b));
  if (a < 0)
    a = 360 + a;
  return(a);
}

/* returns an angle in rads in the range 0 to two pi */
double rangerad(double x) {
  double a, b;
  b = x / TPI;
  a = TPI * (b - floor(b));
  if (a < 0)
    a = TPI + a;
  return(a);
}

/*
gets the atan2 function returning angles in the right
order and  range
*/
double atan22(double y, double x) {
  double a;

  a = atan2(y, x);
  if (a < 0) a += TPI;
  return(a);
}

/*
returns mean obliquity of ecliptic in radians given days since
J2000.0.
*/
double epsilon(double d) {
  double t = d/ 36525;
  return((23.4392911111111 - (t* (46.8150 + 0.00059*t)/3600)) *RADS);
}

/*
replaces ecliptic coordinates with equatorial coordinates
note: call by reference destroys original values
R is unchanged.
*/
void equatorial(double d, double *lon, double *lat, double * r) {
  double  eps, ceps, seps, l, b;
  (void)r;

  l = *lon;
  b = * lat;
  eps = epsilon(d);
  ceps = cos(eps);
  seps = sin(eps);
  *lon = atan22(sin(l)*ceps - tan(b)*seps, cos(l));
  *lat = asin(sin(b)*ceps + cos(b)*seps*sin(l));
}

/*
replaces equatorial coordinates with ecliptic ones. Inverse
of above, but used to find topocentric ecliptic coords.
*/
void ecliptic(double d, double *lon, double *lat, double * r) {
  double  eps, ceps, seps, alp, dec;
  (void)r;

  alp = *lon;
  dec = *lat;
  eps = epsilon(d);
  ceps = cos(eps);
  seps = sin(eps);
  *lon = atan22(sin(alp)*ceps + tan(dec)*seps, cos(alp));
  *lat = asin(sin(dec)*ceps - cos(dec)*seps*sin(alp));
}

/*
returns the siderial time at greenwich meridian as
an angle in radians given the days since J2000.0
*/
double gst( double d) {
  double t = d / 36525;
  double theta;
  theta = range(280.46061837 + 360.98564736629 * d + 0.000387933 * t * t);
  return(theta * RADS);
}

void tmoonsub_(double *day, double *glat, double *glong, double *moonalt,
   double *mrv, double *l, double *b, double *paxis)
{
  double mlambda, mbeta;
  double malpha, mdelta;
  double lst, mhr;
  double tlambda, tbeta, trv;

  lst = gst(*day) + *glong;

  /* find Moon topocentric coordinates for libration calculations */

  moonpos(*day, &mlambda, &mbeta, mrv);
  malpha = mlambda;
  mdelta = mbeta;
  equatorial(*day, &malpha, &mdelta, mrv);
  topo(lst, *glat, &malpha, &mdelta, mrv);
  mhr = rangerad(lst - malpha);
  *moonalt = alt(*glat, mhr, mdelta);

  /* Optical libration and Position angle of the Pole */

  tlambda = malpha;
  tbeta = mdelta;
  trv = *mrv;
  ecliptic(*day, &tlambda, &tbeta, &trv);
  libration(*day, tlambda, tbeta, malpha,  l, b, paxis);
}


code:
/*=============================================================================
    Copyright (c) 1999-2003 Jaakko Jarvi
    Copyright (c) 1999-2003 Jeremiah Willcock
    Copyright (c) 2001-2011 Joel de Guzman

    Distributed under the Boost Software License, Version 1.0. (See accompanying
    file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
==============================================================================*/
#if !defined(FUSION_OUT_05052005_0121)
#define FUSION_OUT_05052005_0121

#include <boost/fusion/support/config.hpp>
#include <ostream>
#include <boost/fusion/sequence/io/detail/manip.hpp>

#include <boost/mpl/bool.hpp>
#include <boost/fusion/sequence/intrinsic/begin.hpp>
#include <boost/fusion/sequence/intrinsic/end.hpp>
#include <boost/fusion/iterator/deref.hpp>
#include <boost/fusion/iterator/next.hpp>
#include <boost/fusion/iterator/equal_to.hpp>

namespace boost { namespace fusion { namespace detail
{
    template <typename Tag>
    struct delimiter_out
    {
        // print a delimiter
        template <typename OS>
        static void
        print(OS& os, char const* delim, mpl::false_ = mpl::false_())
        {
            detail::string_ios_manip<Tag, OS> manip(os);
            manip.print(delim);
        }

        template <typename OS>
        static void
        print(OS&, char const*, mpl::true_)
        {
        }
    };

    struct print_sequence_loop
    {
        template <typename OS, typename First, typename Last>
        static void
        call(OS&, First const&, Last const&, mpl::true_)
        {
        }

        template <typename OS, typename First, typename Last>
        static void
        call(OS& os, First const& first, Last const& last, mpl::false_)
        {
            result_of::equal_to<
                typename result_of::next<First>::type
              , Last
            >
            is_last;

            os << *first;
            delimiter_out<tuple_delimiter_tag>::print(os, " ", is_last);
            call(os, fusion::next(first), last, is_last);
        }

        template <typename OS, typename First, typename Last>
        static void
        call(OS& os, First const& first, Last const& last)
        {
            result_of::equal_to<First, Last> eq;
            call(os, first, last, eq);
        }
    };

    template <typename OS, typename Sequence>
    inline void
    print_sequence(OS& os, Sequence const& seq)
    {
        delimiter_out<tuple_open_tag>::print(os, "(");
        print_sequence_loop::call(os, fusion::begin(seq), fusion::end(seq));
        delimiter_out<tuple_close_tag>::print(os, ")");
    }
}}}

#endif


code:
//  Copyright 2000 John Maddock (john@johnmaddock.co.uk)
//  Copyright 2002 Aleksey Gurtovoy (agurtovoy@meta-comm.com)
//
//  Use, modification and distribution are subject to the Boost Software License,
//  Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
//  http://www.boost.org/LICENSE_1_0.txt).
//
//  See http://www.boost.org/libs/type_traits for most recent version including documentation.

#ifndef BOOST_TT_IS_FUNCTION_HPP_INCLUDED
#define BOOST_TT_IS_FUNCTION_HPP_INCLUDED

#include <boost/type_traits/is_reference.hpp>
#include <boost/type_traits/detail/config.hpp>

#if !defined(BOOST_TT_TEST_MS_FUNC_SIGS)
#   include <boost/type_traits/detail/is_function_ptr_helper.hpp>
#else
#   include <boost/type_traits/detail/is_function_ptr_tester.hpp>
#   include <boost/type_traits/detail/yes_no_type.hpp>
#endif

// is a type a function?
// Please note that this implementation is unnecessarily complex:
// we could just use !is_convertible<T*, const volatile void*>::value,
// except that some compilers erroneously allow conversions from
// function pointers to void*.

namespace boost {

#if !defined( __CODEGEARC__ )

namespace detail {

#if !defined(BOOST_TT_TEST_MS_FUNC_SIGS)
template<bool is_ref = true>
struct is_function_chooser
{
   template< typename T > struct result_
      : public false_type {};
};

template <>
struct is_function_chooser<false>
{
    template< typename T > struct result_
        : public ::boost::type_traits::is_function_ptr_helper<T*> {};
};

template <typename T>
struct is_function_impl
    : public is_function_chooser< ::boost::is_reference<T>::value >
        ::BOOST_NESTED_TEMPLATE result_<T>
{
};

#else

template <typename T>
struct is_function_impl
{
#if BOOST_WORKAROUND(BOOST_MSVC_FULL_VER, >= 140050000)
#pragma warning(push)
#pragma warning(disable:6334)
#endif
    static T* t;
    BOOST_STATIC_CONSTANT(
        bool, value = sizeof(::boost::type_traits::is_function_ptr_tester(t))
        == sizeof(::boost::type_traits::yes_type)
        );
#if BOOST_WORKAROUND(BOOST_MSVC_FULL_VER, >= 140050000)
#pragma warning(pop)
#endif
};

template <typename T>
struct is_function_impl<T&> : public false_type
{};
#ifndef BOOST_NO_CXX11_RVALUE_REFERENCES
template <typename T>
struct is_function_impl<T&&> : public false_type
{};
#endif

#endif

} // namespace detail

#endif // !defined( __CODEGEARC__ )

#if defined( __CODEGEARC__ )
template <class T> struct is_function : integral_constant<bool, __is_function(T)> {};
#else
template <class T> struct is_function : integral_constant<bool, ::boost::detail::is_function_impl<T>::value> {};
#ifndef BOOST_NO_CXX11_RVALUE_REFERENCES
template <class T> struct is_function<T&&> : public false_type {};
#endif
#endif
} // namespace boost

#endif // BOOST_TT_IS_FUNCTION_HPP_INCLUDED


code:
// Copyright Kevlin Henney, 2000-2005.
// Copyright Alexander Nasonov, 2006-2010.
// Copyright Antony Polukhin, 2011-2014.
//
// Distributed under the Boost Software License, Version 1.0. (See
// accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)
//
// what:  lexical_cast custom keyword cast
// who:   contributed by Kevlin Henney,
//        enhanced with contributions from Terje Slettebo,
//        with additional fixes and suggestions from Gennaro Prota,
//        Beman Dawes, Dave Abrahams, Daryle Walker, Peter Dimov,
//        Alexander Nasonov, Antony Polukhin, Justin Viiret, Michael Hofmann,
//        Cheng Yang, Matthew Bradbury, David W. Birdsall, Pavel Korzh and other Boosters
// when:  November 2000, March 2003, June 2005, June 2006, March 2011 - 2014

#ifndef BOOST_LEXICAL_CAST_DETAIL_CONVERTER_LEXICAL_HPP
#define BOOST_LEXICAL_CAST_DETAIL_CONVERTER_LEXICAL_HPP

#include <boost/config.hpp>
#ifdef BOOST_HAS_PRAGMA_ONCE
#   pragma once
#endif

#if defined(BOOST_NO_STRINGSTREAM) || defined(BOOST_NO_STD_WSTRING)
#define BOOST_LCAST_NO_WCHAR_T
#endif

#include <cstddef>
#include <string>
#include <boost/limits.hpp>
#include <boost/mpl/bool.hpp>
#include <boost/mpl/identity.hpp>
#include <boost/mpl/if.hpp>
#include <boost/type_traits/is_integral.hpp>
#include <boost/type_traits/is_float.hpp>
#include <boost/type_traits/has_left_shift.hpp>
#include <boost/type_traits/has_right_shift.hpp>
#include <boost/static_assert.hpp>
#include <boost/detail/lcast_precision.hpp>

#include <boost/lexical_cast/detail/widest_char.hpp>
#include <boost/lexical_cast/detail/is_character.hpp>

#ifndef BOOST_NO_CXX11_HDR_ARRAY
#include <array>
#endif

#include <boost/array.hpp>
#include <boost/range/iterator_range_core.hpp>
#include <boost/container/container_fwd.hpp>

#include <boost/lexical_cast/detail/converter_lexical_streams.hpp>

namespace boost {

    namespace detail // normalize_single_byte_char<Char>
    {
        // Converts signed/unsigned char to char
        template < class Char >
        struct normalize_single_byte_char
        {
            typedef Char type;
        };

        template <>
        struct normalize_single_byte_char< signed char >
        {
            typedef char type;
        };

        template <>
        struct normalize_single_byte_char< unsigned char >
        {
            typedef char type;
        };
    }

    namespace detail // deduce_character_type_later<T>
    {
        // Helper type, meaning that stram character for T must be deduced
        // at Stage 2 (See deduce_source_char<T> and deduce_target_char<T>)
        template < class T > struct deduce_character_type_later {};
    }

    namespace detail // stream_char_common<T>
    {
        // Selectors to choose stream character type (common for Source and Target)
        // Returns one of char, wchar_t, char16_t, char32_t or deduce_character_type_later<T> types
        // Executed on Stage 1 (See deduce_source_char<T> and deduce_target_char<T>)
        template < typename Type >
        struct stream_char_common: public boost::mpl::if_c<
            boost::detail::is_character< Type >::value,
            Type,
            boost::detail::deduce_character_type_later< Type >
        > {};

        template < typename Char >
        struct stream_char_common< Char* >: public boost::mpl::if_c<
            boost::detail::is_character< Char >::value,
            Char,
            boost::detail::deduce_character_type_later< Char* >
        > {};

        template < typename Char >
        struct stream_char_common< const Char* >: public boost::mpl::if_c<
            boost::detail::is_character< Char >::value,
            Char,
            boost::detail::deduce_character_type_later< const Char* >
        > {};

        template < typename Char >
        struct stream_char_common< boost::iterator_range< Char* > >: public boost::mpl::if_c<
            boost::detail::is_character< Char >::value,
            Char,
            boost::detail::deduce_character_type_later< boost::iterator_range< Char* > >
        > {};

        template < typename Char >
        struct stream_char_common< boost::iterator_range< const Char* > >: public boost::mpl::if_c<
            boost::detail::is_character< Char >::value,
            Char,
            boost::detail::deduce_character_type_later< boost::iterator_range< const Char* > >
        > {};

        template < class Char, class Traits, class Alloc >
        struct stream_char_common< std::basic_string< Char, Traits, Alloc > >
        {
            typedef Char type;
        };

        template < class Char, class Traits, class Alloc >
        struct stream_char_common< boost::container::basic_string< Char, Traits, Alloc > >
        {
            typedef Char type;
        };

        template < typename Char, std::size_t N >
        struct stream_char_common< boost::array< Char, N > >: public boost::mpl::if_c<
            boost::detail::is_character< Char >::value,
            Char,
            boost::detail::deduce_character_type_later< boost::array< Char, N > >
        > {};

        template < typename Char, std::size_t N >
        struct stream_char_common< boost::array< const Char, N > >: public boost::mpl::if_c<
            boost::detail::is_character< Char >::value,
            Char,
            boost::detail::deduce_character_type_later< boost::array< const Char, N > >
        > {};

#ifndef BOOST_NO_CXX11_HDR_ARRAY
        template < typename Char, std::size_t N >
        struct stream_char_common< std::array<Char, N > >: public boost::mpl::if_c<
            boost::detail::is_character< Char >::value,
            Char,
            boost::detail::deduce_character_type_later< std::array< Char, N > >
        > {};

        template < typename Char, std::size_t N >
        struct stream_char_common< std::array< const Char, N > >: public boost::mpl::if_c<
            boost::detail::is_character< Char >::value,
            Char,
            boost::detail::deduce_character_type_later< std::array< const Char, N > >
        > {};
#endif

#ifdef BOOST_HAS_INT128
        template <> struct stream_char_common< boost::int128_type >: public boost::mpl::identity< char > {};
        template <> struct stream_char_common< boost::uint128_type >: public boost::mpl::identity< char > {};
#endif

#if !defined(BOOST_LCAST_NO_WCHAR_T) && defined(BOOST_NO_INTRINSIC_WCHAR_T)
        template <>
        struct stream_char_common< wchar_t >
        {
            typedef char type;
        };
#endif
    }

    namespace detail // deduce_source_char_impl<T>
    {
        // If type T is `deduce_character_type_later` type, then tries to deduce
        // character type using boost::has_left_shift<T> metafunction.
        // Otherwise supplied type T is a character type, that must be normalized
        // using normalize_single_byte_char<Char>.
        // Executed at Stage 2  (See deduce_source_char<T> and deduce_target_char<T>)
        template < class Char >
        struct deduce_source_char_impl
        {
            typedef BOOST_DEDUCED_TYPENAME boost::detail::normalize_single_byte_char< Char >::type type;
        };

        template < class T >
        struct deduce_source_char_impl< deduce_character_type_later< T > >
        {
            typedef boost::has_left_shift< std::basic_ostream< char >, T > result_t;

#if defined(BOOST_LCAST_NO_WCHAR_T)
            BOOST_STATIC_ASSERT_MSG((result_t::value),
                "Source type is not std::ostream`able and std::wostream`s are not supported by your STL implementation");
            typedef char type;
#else
            typedef BOOST_DEDUCED_TYPENAME boost::mpl::if_c<
                result_t::value, char, wchar_t
            >::type type;

            BOOST_STATIC_ASSERT_MSG((result_t::value || boost::has_left_shift< std::basic_ostream< type >, T >::value),
                "Source type is neither std::ostream`able nor std::wostream`able");
#endif
        };
    }

    namespace detail  // deduce_target_char_impl<T>
    {
        // If type T is `deduce_character_type_later` type, then tries to deduce
        // character type using boost::has_right_shift<T> metafunction.
        // Otherwise supplied type T is a character type, that must be normalized
        // using normalize_single_byte_char<Char>.
        // Executed at Stage 2  (See deduce_source_char<T> and deduce_target_char<T>)
        template < class Char >
        struct deduce_target_char_impl
        {
            typedef BOOST_DEDUCED_TYPENAME normalize_single_byte_char< Char >::type type;
        };

        template < class T >
        struct deduce_target_char_impl< deduce_character_type_later<T> >
        {
            typedef boost::has_right_shift<std::basic_istream<char>, T > result_t;

#if defined(BOOST_LCAST_NO_WCHAR_T)
            BOOST_STATIC_ASSERT_MSG((result_t::value),
                "Target type is not std::istream`able and std::wistream`s are not supported by your STL implementation");
            typedef char type;
#else
            typedef BOOST_DEDUCED_TYPENAME boost::mpl::if_c<
                result_t::value, char, wchar_t
            >::type type;

            BOOST_STATIC_ASSERT_MSG((result_t::value || boost::has_right_shift<std::basic_istream<wchar_t>, T >::value),
                "Target type is neither std::istream`able nor std::wistream`able");
#endif
        };
    }

    namespace detail  // deduce_target_char<T> and deduce_source_char<T>
    {
        // We deduce stream character types in two stages.
        //
        // Stage 1 is common for Target and Source. At Stage 1 we get
        // non normalized character type (may contain unsigned/signed char)
        // or deduce_character_type_later<T> where T is the original type.
        // Stage 1 is executed by stream_char_common<T>
        //
        // At Stage 2 we normalize character types or try to deduce character
        // type using metafunctions.
        // Stage 2 is executed by deduce_target_char_impl<T> and
        // deduce_source_char_impl<T>
        //
        // deduce_target_char<T> and deduce_source_char<T> functions combine
        // both stages

        template < class T >
        struct deduce_target_char
        {
            typedef BOOST_DEDUCED_TYPENAME stream_char_common< T >::type stage1_type;
            typedef BOOST_DEDUCED_TYPENAME deduce_target_char_impl< stage1_type >::type stage2_type;

            typedef stage2_type type;
        };

        template < class T >
        struct deduce_source_char
        {
            typedef BOOST_DEDUCED_TYPENAME stream_char_common< T >::type stage1_type;
            typedef BOOST_DEDUCED_TYPENAME deduce_source_char_impl< stage1_type >::type stage2_type;

            typedef stage2_type type;
        };
    }

    namespace detail // extract_char_traits template
    {
        // We are attempting to get char_traits<> from T
        // template parameter. Otherwise we'll be using std::char_traits<Char>
        template < class Char, class T >
        struct extract_char_traits
                : boost::false_type
        {
            typedef std::char_traits< Char > trait_t;
        };

        template < class Char, class Traits, class Alloc >
        struct extract_char_traits< Char, std::basic_string< Char, Traits, Alloc > >
            : boost::true_type
        {
            typedef Traits trait_t;
        };

        template < class Char, class Traits, class Alloc>
        struct extract_char_traits< Char, boost::container::basic_string< Char, Traits, Alloc > >
            : boost::true_type
        {
            typedef Traits trait_t;
        };
    }

    namespace detail // array_to_pointer_decay<T>
    {
        template<class T>
        struct array_to_pointer_decay
        {
            typedef T type;
        };

        template<class T, std::size_t N>
        struct array_to_pointer_decay<T[N]>
        {
            typedef const T * type;
        };
    }

    namespace detail // lcast_src_length
    {
        // Return max. length of string representation of Source;
        template< class Source,         // Source type of lexical_cast.
                  class Enable = void   // helper type
                >
        struct lcast_src_length
        {
            BOOST_STATIC_CONSTANT(std::size_t, value = 1);
        };

        // Helper for integral types.
        // Notes on length calculation:
        // Max length for 32bit int with grouping "\1" and thousands_sep ',':
        // "-2,1,4,7,4,8,3,6,4,7"
        //  ^                    - is_signed
        //   ^                   - 1 digit not counted by digits10
        //    ^^^^^^^^^^^^^^^^^^ - digits10 * 2
        //
        // Constant is_specialized is used instead of constant 1
        // to prevent buffer overflow in a rare case when
        // <boost/limits.hpp> doesn't add missing specialization for
        // numeric_limits<T> for some integral type T.
        // When is_specialized is false, the whole expression is 0.
        template <class Source>
        struct lcast_src_length<
                    Source, BOOST_DEDUCED_TYPENAME boost::enable_if<boost::is_integral<Source> >::type
                >
        {
#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
            BOOST_STATIC_CONSTANT(std::size_t, value =
                  std::numeric_limits<Source>::is_signed +
                  std::numeric_limits<Source>::is_specialized + /* == 1 */
                  std::numeric_limits<Source>::digits10 * 2
              );
#else
            BOOST_STATIC_CONSTANT(std::size_t, value = 156);
            BOOST_STATIC_ASSERT(sizeof(Source) * CHAR_BIT <= 256);
#endif
        };

        // Helper for floating point types.
        // -1.23456789e-123456
        // ^                   sign
        //  ^                  leading digit
        //   ^                 decimal point
        //    ^^^^^^^^         lcast_precision<Source>::value
        //            ^        "e"
        //             ^       exponent sign
        //              ^^^^^^ exponent (assumed 6 or less digits)
        // sign + leading digit + decimal point + "e" + exponent sign == 5
        template<class Source>
        struct lcast_src_length<
                Source, BOOST_DEDUCED_TYPENAME boost::enable_if<boost::is_float<Source> >::type
            >
        {

#ifndef BOOST_LCAST_NO_COMPILE_TIME_PRECISION
            BOOST_STATIC_ASSERT(
                    std::numeric_limits<Source>::max_exponent10 <=  999999L &&
                    std::numeric_limits<Source>::min_exponent10 >= -999999L
                );

            BOOST_STATIC_CONSTANT(std::size_t, value =
                    5 + lcast_precision<Source>::value + 6
                );
#else // #ifndef BOOST_LCAST_NO_COMPILE_TIME_PRECISION
            BOOST_STATIC_CONSTANT(std::size_t, value = 156);
#endif // #ifndef BOOST_LCAST_NO_COMPILE_TIME_PRECISION
        };
    }

    namespace detail // lexical_cast_stream_traits<Source, Target>
    {
        template <class Source, class Target>
        struct lexical_cast_stream_traits {
            typedef BOOST_DEDUCED_TYPENAME boost::detail::array_to_pointer_decay<Source>::type src;
            typedef BOOST_DEDUCED_TYPENAME boost::remove_cv<src>::type            no_cv_src;

            typedef boost::detail::deduce_source_char<no_cv_src>                           deduce_src_char_metafunc;
            typedef BOOST_DEDUCED_TYPENAME deduce_src_char_metafunc::type           src_char_t;
            typedef BOOST_DEDUCED_TYPENAME boost::detail::deduce_target_char<Target>::type target_char_t;

            typedef BOOST_DEDUCED_TYPENAME boost::detail::widest_char<
                target_char_t, src_char_t
            >::type char_type;

#if !defined(BOOST_NO_CXX11_CHAR16_T) && defined(BOOST_NO_CXX11_UNICODE_LITERALS)
            BOOST_STATIC_ASSERT_MSG(( !boost::is_same<char16_t, src_char_t>::value
                                        && !boost::is_same<char16_t, target_char_t>::value),
                "Your compiler does not have full support for char16_t" );
#endif
#if !defined(BOOST_NO_CXX11_CHAR32_T) && defined(BOOST_NO_CXX11_UNICODE_LITERALS)
            BOOST_STATIC_ASSERT_MSG(( !boost::is_same<char32_t, src_char_t>::value
                                        && !boost::is_same<char32_t, target_char_t>::value),
                "Your compiler does not have full support for char32_t" );
#endif

            typedef BOOST_DEDUCED_TYPENAME boost::mpl::if_c<
                boost::detail::extract_char_traits<char_type, Target>::value,
                BOOST_DEDUCED_TYPENAME boost::detail::extract_char_traits<char_type, Target>,
                BOOST_DEDUCED_TYPENAME boost::detail::extract_char_traits<char_type, no_cv_src>
            >::type::trait_t traits;

            typedef boost::mpl::bool_
            	<
                boost::is_same<char, src_char_t>::value &&                                 // source is not a wide character based type
                (sizeof(char) != sizeof(target_char_t)) &&  // target type is based on wide character
                (!(boost::detail::is_character<no_cv_src>::value))
            	> is_string_widening_required_t;

            typedef boost::mpl::bool_
            	<
            	!(boost::is_integral<no_cv_src>::value ||
                  boost::detail::is_character<
                    BOOST_DEDUCED_TYPENAME deduce_src_char_metafunc::stage1_type          // if we did not get character type at stage1
                  >::value                                                           // then we have no optimization for that type
            	 )
            	> is_source_input_not_optimized_t;

            // If we have an optimized conversion for
            // Source, we do not need to construct stringbuf.
            BOOST_STATIC_CONSTANT(bool, requires_stringbuf =
            	(is_string_widening_required_t::value || is_source_input_not_optimized_t::value)
            );

            typedef boost::detail::lcast_src_length<no_cv_src> len_t;
        };
    }

    namespace detail
    {
        template<typename Target, typename Source>
        struct lexical_converter_impl
        {
            typedef lexical_cast_stream_traits<Source, Target>  stream_trait;

            typedef detail::lexical_istream_limited_src<
                BOOST_DEDUCED_TYPENAME stream_trait::char_type,
                BOOST_DEDUCED_TYPENAME stream_trait::traits,
                stream_trait::requires_stringbuf,
                stream_trait::len_t::value + 1
            > i_interpreter_type;

            typedef detail::lexical_ostream_limited_src<
                BOOST_DEDUCED_TYPENAME stream_trait::char_type,
                BOOST_DEDUCED_TYPENAME stream_trait::traits
            > o_interpreter_type;

            static inline bool try_convert(const Source& arg, Target& result) {
                i_interpreter_type i_interpreter;

                // Disabling ADL, by directly specifying operators.
                if (!(i_interpreter.operator <<(arg)))
                    return false;

                o_interpreter_type out(i_interpreter.cbegin(), i_interpreter.cend());

                // Disabling ADL, by directly specifying operators.
                if(!(out.operator >>(result)))
                    return false;

                return true;
            }
        };
    }

} // namespace boost

#undef BOOST_LCAST_NO_WCHAR_T

#endif // BOOST_LEXICAL_CAST_DETAIL_CONVERTER_LEXICAL_HPP


code:
//---------------------------------------------------------------------------//
// Copyright (c) 2013 Kyle Lutz <kyle.r.lutz@gmail.com>
//
// Distributed under the Boost Software License, Version 1.0
// See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt
//
// See http://boostorg.github.com/compute for more information.
//---------------------------------------------------------------------------//

#ifndef BOOST_COMPUTE_LAMBDA_CONTEXT_HPP
#define BOOST_COMPUTE_LAMBDA_CONTEXT_HPP

#include <boost/proto/core.hpp>
#include <boost/proto/context.hpp>
#include <boost/type_traits.hpp>
#include <boost/preprocessor/repetition.hpp>

#include <boost/compute/config.hpp>
#include <boost/compute/function.hpp>
#include <boost/compute/lambda/result_of.hpp>
#include <boost/compute/lambda/functional.hpp>
#include <boost/compute/type_traits/result_of.hpp>
#include <boost/compute/type_traits/type_name.hpp>
#include <boost/compute/detail/meta_kernel.hpp>

namespace boost {
namespace compute {
namespace lambda {

namespace mpl = boost::mpl;
namespace proto = boost::proto;

#define BOOST_COMPUTE_LAMBDA_CONTEXT_DEFINE_BINARY_OPERATOR(tag, op) \
    template<class LHS, class RHS> \
    void operator()(tag, const LHS &lhs, const RHS &rhs) \
    { \
        if(proto::arity_of<LHS>::value > 0){ \
            stream << '('; \
            proto::eval(lhs, *this); \
            stream << ')'; \
        } \
        else { \
            proto::eval(lhs, *this); \
        } \
        \
        stream << op; \
        \
        if(proto::arity_of<RHS>::value > 0){ \
            stream << '('; \
            proto::eval(rhs, *this); \
            stream << ')'; \
        } \
        else { \
            proto::eval(rhs, *this); \
        } \
    }

// lambda expression context
template<class Args>
struct context : proto::callable_context<context<Args> >
{
    typedef void result_type;
    typedef Args args_tuple;

    // create a lambda context for kernel with args
    context(boost::compute::detail::meta_kernel &kernel, const Args &args_)
        : stream(kernel),
          args(args_)
    {
    }

    // handle terminals
    template<class T>
    void operator()(proto::tag::terminal, const T &x)
    {
        // terminal values in lambda expressions are always literals
        stream << stream.lit(x);
    }

    // handle placeholders
    template<int I>
    void operator()(proto::tag::terminal, placeholder<I>)
    {
        stream << boost::get<I>(args);
    }

    // handle functions
    #define BOOST_COMPUTE_LAMBDA_CONTEXT_FUNCTION_ARG(z, n, unused) \
        BOOST_PP_COMMA_IF(n) BOOST_PP_CAT(const Arg, n) BOOST_PP_CAT(&arg, n)

    #define BOOST_COMPUTE_LAMBDA_CONTEXT_FUNCTION(z, n, unused) \
    template<class F, BOOST_PP_ENUM_PARAMS(n, class Arg)> \
    void operator()( \
        proto::tag::function, \
        const F &function, \
        BOOST_PP_REPEAT(n, BOOST_COMPUTE_LAMBDA_CONTEXT_FUNCTION_ARG, ~) \
    ) \
    { \
        proto::value(function).apply(*this, BOOST_PP_ENUM_PARAMS(n, arg)); \
    }

    BOOST_PP_REPEAT_FROM_TO(1, BOOST_COMPUTE_MAX_ARITY, BOOST_COMPUTE_LAMBDA_CONTEXT_FUNCTION, ~)

    #undef BOOST_COMPUTE_LAMBDA_CONTEXT_FUNCTION

    // operators
    BOOST_COMPUTE_LAMBDA_CONTEXT_DEFINE_BINARY_OPERATOR(proto::tag::plus, '+')
    BOOST_COMPUTE_LAMBDA_CONTEXT_DEFINE_BINARY_OPERATOR(proto::tag::minus, '-')
    BOOST_COMPUTE_LAMBDA_CONTEXT_DEFINE_BINARY_OPERATOR(proto::tag::multiplies, '*')
    BOOST_COMPUTE_LAMBDA_CONTEXT_DEFINE_BINARY_OPERATOR(proto::tag::divides, '/')
    BOOST_COMPUTE_LAMBDA_CONTEXT_DEFINE_BINARY_OPERATOR(proto::tag::modulus, '%')
    BOOST_COMPUTE_LAMBDA_CONTEXT_DEFINE_BINARY_OPERATOR(proto::tag::less, '<')
    BOOST_COMPUTE_LAMBDA_CONTEXT_DEFINE_BINARY_OPERATOR(proto::tag::greater, '>')
    BOOST_COMPUTE_LAMBDA_CONTEXT_DEFINE_BINARY_OPERATOR(proto::tag::less_equal, "<=")
    BOOST_COMPUTE_LAMBDA_CONTEXT_DEFINE_BINARY_OPERATOR(proto::tag::greater_equal, ">=")
    BOOST_COMPUTE_LAMBDA_CONTEXT_DEFINE_BINARY_OPERATOR(proto::tag::equal_to, "==")
    BOOST_COMPUTE_LAMBDA_CONTEXT_DEFINE_BINARY_OPERATOR(proto::tag::not_equal_to, "!=")
    BOOST_COMPUTE_LAMBDA_CONTEXT_DEFINE_BINARY_OPERATOR(proto::tag::logical_and, "&&")
    BOOST_COMPUTE_LAMBDA_CONTEXT_DEFINE_BINARY_OPERATOR(proto::tag::logical_or, "||")
    BOOST_COMPUTE_LAMBDA_CONTEXT_DEFINE_BINARY_OPERATOR(proto::tag::bitwise_and, '&')
    BOOST_COMPUTE_LAMBDA_CONTEXT_DEFINE_BINARY_OPERATOR(proto::tag::bitwise_or, '|')
    BOOST_COMPUTE_LAMBDA_CONTEXT_DEFINE_BINARY_OPERATOR(proto::tag::bitwise_xor, '^')
    BOOST_COMPUTE_LAMBDA_CONTEXT_DEFINE_BINARY_OPERATOR(proto::tag::assign, '=')

    // subscript operator
    template<class LHS, class RHS>
    void operator()(proto::tag::subscript, const LHS &lhs, const RHS &rhs)
    {
        proto::eval(lhs, *this);
        stream << '[';
        proto::eval(rhs, *this);
        stream << ']';
    }

    // ternary conditional operator
    template<class Pred, class Arg1, class Arg2>
    void operator()(proto::tag::if_else_, const Pred &p, const Arg1 &x, const Arg2 &y)
    {
        proto::eval(p, *this);
        stream << '?';
        proto::eval(x, *this);
        stream << ':';
        proto::eval(y, *this);
    }

    boost::compute::detail::meta_kernel &stream;
    Args args;
};

namespace detail {

template<class Expr, class Arg>
struct invoked_unary_expression
{
    typedef typename ::boost::compute::result_of<Expr(Arg)>::type result_type;

    invoked_unary_expression(const Expr &expr, const Arg &arg)
        : m_expr(expr),
          m_arg(arg)
    {
    }

    Expr m_expr;
    Arg m_arg;
};

template<class Expr, class Arg>
boost::compute::detail::meta_kernel&
operator<<(boost::compute::detail::meta_kernel &kernel,
           const invoked_unary_expression<Expr, Arg> &expr)
{
    context<boost::tuple<Arg> > ctx(kernel, boost::make_tuple(expr.m_arg));
    proto::eval(expr.m_expr, ctx);

    return kernel;
}

template<class Expr, class Arg1, class Arg2>
struct invoked_binary_expression
{
    typedef typename ::boost::compute::result_of<Expr(Arg1, Arg2)>::type result_type;

    invoked_binary_expression(const Expr &expr,
                              const Arg1 &arg1,
                              const Arg2 &arg2)
        : m_expr(expr),
          m_arg1(arg1),
          m_arg2(arg2)
    {
    }

    Expr m_expr;
    Arg1 m_arg1;
    Arg2 m_arg2;
};

template<class Expr, class Arg1, class Arg2>
boost::compute::detail::meta_kernel&
operator<<(boost::compute::detail::meta_kernel &kernel,
           const invoked_binary_expression<Expr, Arg1, Arg2> &expr)
{
    context<boost::tuple<Arg1, Arg2> > ctx(
        kernel,
        boost::make_tuple(expr.m_arg1, expr.m_arg2)
    );
    proto::eval(expr.m_expr, ctx);

    return kernel;
}

} // end detail namespace

// forward declare domain
struct domain;

// lambda expression wrapper
template<class Expr>
struct expression : proto::extends<Expr, expression<Expr>, domain>
{
    typedef proto::extends<Expr, expression<Expr>, domain> base_type;

    BOOST_PROTO_EXTENDS_USING_ASSIGN(expression)

    expression(const Expr &expr = Expr())
        : base_type(expr)
    {
    }

    // result_of protocol
    template<class Signature>
    struct result
    {
    };

    template<class This>
    struct result<This()>
    {
        typedef
            typename ::boost::compute::lambda::result_of<Expr>::type type;
    };

    template<class This, class Arg>
    struct result<This(Arg)>
    {
        typedef
            typename ::boost::compute::lambda::result_of<
                Expr,
                typename boost::tuple<Arg>
            >::type type;
    };

    template<class This, class Arg1, class Arg2>
    struct result<This(Arg1, Arg2)>
    {
        typedef typename
            ::boost::compute::lambda::result_of<
                Expr,
                typename boost::tuple<Arg1, Arg2>
            >::type type;
    };

    template<class Arg>
    detail::invoked_unary_expression<expression<Expr>, Arg>
    operator()(const Arg &x) const
    {
        return detail::invoked_unary_expression<expression<Expr>, Arg>(*this, x);
    }

    template<class Arg1, class Arg2>
    detail::invoked_binary_expression<expression<Expr>, Arg1, Arg2>
    operator()(const Arg1 &x, const Arg2 &y) const
    {
        return detail::invoked_binary_expression<
                   expression<Expr>,
                   Arg1,
                   Arg2
                >(*this, x, y);
    }

    // function<> conversion operator
    template<class R, class A1>
    operator function<R(A1)>() const
    {
        using ::boost::compute::detail::meta_kernel;

        std::stringstream source;

        ::boost::compute::detail::meta_kernel_variable<A1> arg1("x");

        source << "inline " << type_name<R>() << " lambda"
               << ::boost::compute::detail::generate_argument_list<R(A1)>('x')
               << "{\n"
               << "    return " << meta_kernel::expr_to_string((*this)(arg1)) << ";\n"
               << "}\n";

        return make_function_from_source<R(A1)>("lambda", source.str());
    }

    template<class R, class A1, class A2>
    operator function<R(A1, A2)>() const
    {
        using ::boost::compute::detail::meta_kernel;

        std::stringstream source;

        ::boost::compute::detail::meta_kernel_variable<A1> arg1("x");
        ::boost::compute::detail::meta_kernel_variable<A1> arg2("y");

        source << "inline " << type_name<R>() << " lambda"
               << ::boost::compute::detail::generate_argument_list<R(A1, A2)>('x')
               << "{\n"
               << "    return " << meta_kernel::expr_to_string((*this)(arg1, arg2)) << ";\n"
               << "}\n";

        return make_function_from_source<R(A1, A2)>("lambda", source.str());
    }
};

// lambda expression domain
struct domain : proto::domain<proto::generator<expression> >
{
};

} // end lambda namespace
} // end compute namespace
} // end boost namespace

#endif // BOOST_COMPUTE_LAMBDA_CONTEXT_HPP


code:
program QRA64code

! Provides examples of message packing, bit and symbol ordering,
! QRA (63,12) encoding, and other necessary details of the QRA64
! protocol.

  use packjt
  character*22 msg,msg0,msg1,decoded,cok*3,msgtype*10,arg*12
  character*6 mycall
  logical ltext
  integer dgen(12),sent(63),dec(12)
  integer icos7(0:6)
  data icos7/2,5,6,0,4,1,3/     !Defines a 7x7 Costas array

  include 'testmsg.f90'

  nargs=iargc()
  if(nargs.lt.1) then
     print*,'Usage: qra64code "message"'
     print*,'       qra64code -t'
     go to 999
  endif

  call getarg(1,msg)                     !Get message from command line
  nmsg=1
  if(msg(1:2).eq."-t") nmsg=NTEST

  write(*,1010)
1010 format("     Message                 Decoded                Err? Type"/74("-"))

  do imsg=1,nmsg
     if(nmsg.gt.1) msg=testmsg(imsg)
    call fmtmsg(msg,iz)                    !To upper, collapse mult blanks
     msg0=msg                               !Input message
     call chkmsg(msg,cok,nspecial,flip)     !See if it includes "OOO" report
     msg1=msg                               !Message without "OOO"
     call packmsg(msg1,dgen,itype,.false.) !Pack message into 12 six-bit bytes
     msgtype=""
     if(itype.eq.1) msgtype="Std Msg"
     if(itype.eq.2) msgtype="Type 1 pfx"
     if(itype.eq.3) msgtype="Type 1 sfx"
     if(itype.eq.4) msgtype="Type 2 pfx"
     if(itype.eq.5) msgtype="Type 2 sfx"
     if(itype.eq.6) msgtype="Free text"

     call qra64_enc(dgen,sent)              !Encode using QRA64

     call unpackmsg(dgen,decoded,.false.,'      ') !Unpack the user message
     call fmtmsg(decoded,iz)
     ii=imsg
     write(*,1020) ii,msg0,decoded,itype,msgtype
1020 format(i4,1x,a22,2x,a22,4x,i3,": ",a13)
  enddo

  if(nmsg.eq.1) then
     write(*,1030) dgen
1030 format(/'Packed message, 6-bit symbols ',12i3) !Display packed symbols

     write(*,1040) sent
1040 format(/'Information-carrying channel symbols'/(i5,29i3))

     write(*,1050) 10*icos7,sent(1:32),10*icos7,sent(33:63),10*icos7
1050 format(/'Channel symbols including sync'/(i5,29i3))
  endif

999 end program QRA64code


code:
#ifndef BOOST_WEAK_PTR_HPP_INCLUDED
#define BOOST_WEAK_PTR_HPP_INCLUDED

//
//  weak_ptr.hpp
//
//  Copyright (c) 2001, 2002, 2003 Peter Dimov
//
//  Distributed under the Boost Software License, Version 1.0.
//  See accompanying file LICENSE_1_0.txt or copy at
//  http://www.boost.org/LICENSE_1_0.txt
//
//  See http://www.boost.org/libs/smart_ptr/weak_ptr.htm for documentation.
//

#include <boost/smart_ptr/weak_ptr.hpp>

#endif  // #ifndef BOOST_WEAK_PTR_HPP_INCLUDED


code:
subroutine watterson(c,npts,fs,delay,fspread)

  complex c(0:npts-1)
  complex c2(0:npts-1)
  complex cs1(0:npts-1)
  complex cs2(0:npts-1)

  nonzero=0
  df=fs/npts
  if(fspread.gt.0.0) then
     do i=0,npts-1
        xx=gran()
        yy=gran()
        cs1(i)=0.707*cmplx(xx,yy)
        xx=gran()
        yy=gran()
        cs2(i)=0.707*cmplx(xx,yy)
     enddo
     call four2a(cs1,npts,1,-1,1)     !To freq domain
     call four2a(cs2,npts,1,-1,1)
     do i=0,npts-1
        f=i*df
        if(i.gt.npts/2) f=(i-npts)*df
        x=(f/(0.5*fspread))**2
        a=0.
        if(x.le.50.0) then
           a=exp(-x)
        endif
        cs1(i)=a*cs1(i)
        cs2(i)=a*cs2(i)
        if(abs(f).lt.10.0) then
           p1=real(cs1(i))**2 + aimag(cs1(i))**2
           p2=real(cs2(i))**2 + aimag(cs2(i))**2
           if(p1.gt.0.0) nonzero=nonzero+1
!           write(62,3101) f,p1,p2,db(p1+1.e-12)-60,db(p2+1.e-12)-60
!3101       format(f10.3,2f12.3,2f10.3)
        endif
     enddo
     call four2a(cs1,npts,1,1,1)     !Back to time domain
     call four2a(cs2,npts,1,1,1)
     cs1(0:npts-1)=cs1(0:npts-1)/npts
     cs2(0:npts-1)=cs2(0:npts-1)/npts
  endif

  nshift=nint(0.001*delay*fs)
  c2(0:npts-1)=cshift(c(0:npts-1),nshift)
  sq=0.
  do i=0,npts-1
     if(nonzero.gt.1) then
        c(i)=0.5*(cs1(i)*c(i) + cs2(i)*c2(i))
     else
        c(i)=0.5*(c(i) + c2(i))
     endif
     sq=sq + real(c(i))**2 + aimag(c(i))**2
!     write(61,3001) i/12000.0,c(i)
!3001 format(3f12.6)
  enddo
  rms=sqrt(sq/npts)
  c=c/rms

  return
end subroutine watterson