Re: [boost] [units] Crazy idea/question

From: Matthias Schabel (boost_at_[hidden])
Date: 2007-02-16 23:29:10

Hi Michael,

> Same here, I was only wondering aloud.  In fact I hope Matthias will
> *only* provide the compile-time version in this first iteration.

So do I ;^) Thanks for the support.

> Well, I have a question to Matthias concerning similar things.  How
> does one use a unitless number?  Do you use dimensionless?

There is a dimensionless tag, but quantities of dimensionless units are
specialized so that they behave just like regular value_types - by
definition if something is dimensionless, there is no difference in its
value, irrespective of unit system. It is mainly there to allow checking
for results of dimensional calculations - if you get a dimensionless
unit, it came from a dimensional analysis calculation. But there is no
real benefit to explicitly constructing dimensionless quantities in
general; for unit computations (without associated values) it can be
nice to be able to output explicitly that the result of some calculation
is dimensionless.

> I posted two formulas that are from a radar handbook.  The first is:
>
> h = R^2 / 2ka
>
> where h is the beam height, R is radar range, k is the refractivity
> (usually 4/3), and a is the Earth's radius. R and a must be the same
> unit, and h is returned in that unit.  How would you code this using
> mcs::units?

For lengths in meters, this would be :

template<typename T>
quantity<SI::length,T>
radar_beam_height(const quantity<SI::length,T>& radar_range,
                                     const quantity<SI::length,T>& earth_radius,
                                     T k = 4.0/3.0)
{
        return quantity<SI::length,T>(pow<2>(radar_range)/ 
(2.0*k*earth_radius));
}

To be more general, we can specialize for arbitrary unit system:

template<class System,typename T>
quantity<unit<System,length_type>,T>
radar_beam_height(const quantity<unit<System,length_type>,T>&  
radar_range,
                                     const quantity<unit<System,length_type>,T>& earth_radius,
                                     T k = 4.0/3.0)
{
        return quantity<unit<System,length_type>,T>(pow<2>(radar_range)/ 
(2.0*k*earth_radius));
}

This function works as long as all quantities are of the same system.  
You could further
generalize it to work with three different types of lengths : one for  
the radar range,
one for the earth radius, and one for the beam height:

template<class return_type,class System1,class System2,typename T>
return_type
radar_beam_height(const quantity<unit<System1,length_type>,T>&  
radar_range,
                                     const quantity<unit<System2,length_type>,T>& earth_radius,
                                     T k = 4.0/3.0)
{
        // need to decide which system to use for calculation
        const return_type	rr(radar_range),
                                        er(earth_radius);

        return return_type (pow<2>(rr)/(2.0*k*er));
}

This way, if your code was like this:

quantity<nautical::length>		radar_range(24.5*miles);
quantity<SI::length>			earth_radius(6371e3*meters);

you could do this:

quantity<CGS::length>		beam_height = radar_beam_height<CGS::system> 
(radar_range,earth_radius);

and get the value computed in centimeters.

> template <typename T>
> T radar_beam_height(T range, T earth_radius,

I would return a quantity of length and pass arguments of length to  
retain type safety

> quantity<SomeSystem::dimensionless> refractivity =
> quantity<SomeSystem::dimensionless>(4/3))

No need, dimensionless is equivalent to T. Watch out for integer  
division in 4/3...

>    const quantity<SomeSystem::dimensionless> two(2);

Again, just a scalar T is fine here.

>    return quantity<T>( pow(range, 2) / ( two * refractivity *  
> range ) );

You need to use the static power function since powers/roots change  
the dimensions
of the units, and a return quantity in the units you like:

return quantity<SI::length>(pow<2>(range)/(two*refractivity*range));

> With calling code like:
> // What's the radar beam height at 300 nautical miles using the mean
> Earth radius?
> quantity<SI::length> height =
> radar_beam_height(quantity<SI::length>(1028971200),
> quantity<SI::length>(6371008.7714));

You can make this simpler:

quantity<SI::length>	height =
        radar_beam_height(quantity<SI::length>(quantity<nautical::length> 
(300.0)),
                                             quantity<SI::length>(6371.0087714*kilo*meters));

> The other question I have is very open-ended.  An approximation for
> the above formula where R is given in nautical miles, and h is
> returned as feet is:
>
> h = (R / 1.23)^2
>
> How would you go about implementing that with mcs::units?  I have no
> idea what unit 1.23 is, or how they even derived that number, other
> than it has the Earth radius and conversion from nautical miles to
> feet built-in.

Engineering approximations can be problematic in the sense that they are
often completely empirical and, therefore, don't particularly respect  
the rules
of dimensional analysis. In these cases it's best to just wrap the  
approximation
in a function that takes the correct arguments and returns the  
appropriate value:

quantity<imperial::length> radar_beam_height(const  
quantity<nautical::length>& range)
{
        return quantity<imperial::length>(std::pow(range.value()/1.23,2.0) 
*imperial::feet);
}

If you really wanted to think about the units of your approximation  
in great detail,
you have :

feet = (nautical miles/(something))^2

so

something = nautical miles/feet^(1/2)

and your constant is

C = 1.23 nm/ft^1/2

This is actually a good example of where heterogeneous units would be  
useful...
Anyway, check out the attached example code for functional versions  
of the
above. You will also need to patch conversion.hpp - there is a minor  
bug in it
that I needed to fix.

> quantity<imperial::feet> radar_beam_height 
> (quantity<nautical::miles> range)

It's important to maintain a semantic distinction in this library :  
by convention, I
have the template arguments refer to the type of unit (length, mass,  
etc...) rather
than the specific unit (meters, kilograms, etc...)...

An important proviso, though : because you cannot use floating point  
numbers as
template arguments, there will be some engineering approximations  
that cannot
be expressed in dimensionally correct form (those with floating point  
powers and/or
roots). However, those formulas are best wrapped in a dimensionally  
correct
wrapper function in any case, to guarantee that the intended unit is  
correct.

Matthias