Lambert.cxx

Go to the documentation of this file.

#include "Lambert.h"

#include "UnaryTransform.h"

#include "axes/AxisModelBase.h"
#include "axes/AxisTick.h"

#include "graphics/Rectangle.h"

#include <cmath>
#include <cstdio>

#include <cassert>

using namespace hippodraw;

using std::max;
using std::abs;
using std::vector;

Lambert::Lambert ( UnaryTransform * z )
  : PeriodicBinaryTransform ( z, true, true, false, true,
                              -180.0, +180.0,
                              - 90.0, + 90.0 )
{
  m_name = "Lambert";
}

Lambert::Lambert ( const Lambert & t )
  :  PeriodicBinaryTransform ( t )
{
}

Lambert::~Lambert ()
{
}

#ifdef CLONE_DEFECT
TransformBase * Lambert::clone () const
#else
Lambert    * Lambert::clone () const
#endif
{
  return new Lambert ( *this );
}

bool
Lambert::
isLinearInXY () const
{
  return false;
}

/* virtual */
void
Lambert::
transform ( double & lon, double & lat ) const
{
  double lon_r = M_PI * lon / 180.0; // convert to radians from degrees
  double lat_r = M_PI * lat / 180.0; // convert to radians from degrees

  double q = 2 * sin ( ( M_PI / 2 - lat_r ) / 2 );
  
  lon = ( 180.0 / M_PI ) * q * sin ( lon_r );// Tranform and convert to degrees
  lat = ( 180.0 / M_PI ) * q * cos ( lon_r );// Tranform and convert to degrees
}

void
Lambert::
inverseTransform ( double & lon, double & lat ) const
{
  double x  = lon * M_PI / 180.0; // Convert to radians
  double y  = lat * M_PI / 180.0; // Convert to radians
  double lon_r = atan2(x ,  y);   // Gets longitude in radians
                                  // (atan2 ensure value is in -M_PI to M_PI)  

  double lat_r =  M_PI / 2 - 2 * asin( 0.5 * sqrt ( x * x + y * y ) );

  lon = ( 180.0 / M_PI ) * lon_r;
  lat = ( 180.0 / M_PI ) * lat_r;
}

/* virtual */
void
Lambert::
transform ( std::vector< double > & lon,
            std::vector< double > & lat ) const
{
  assert ( lat.size() == lon.size() );

  std::vector < double >:: iterator it_lat = lat.begin ();
  std::vector < double >:: iterator it_lon = lon.begin ();

  for ( ; it_lat != lat.end(); ++it_lat, ++it_lon ) {
    transform ( *it_lon, *it_lat );
  }
}

/* virtual */
double Lambert::aspectRatio () const
{
  return 1.0;
}

HippoRectangle Lambert::calcRectangle ( const Range & lat,
                                        const Range & lon )
{
  double x_lo = lat.low ();
  double x_hi = lat.high ();
  
  double y_lo = lon.low ();
  double y_hi = lon.high ();

  double x, y;
  
  double x_min =  1000;
  double x_max = -1000;
  double y_min =  1000;
  double y_max = -1000;

  int n = 50;
  double dx = ( x_hi - x_lo ) / n;
  double dy = ( y_hi - y_lo ) / n;

  // Finding out xmin, xmax, ymin, ymax along line y =  y_lo
  for ( int i = 0; i <= n; i++)
    {
      x = x_lo + i * dx; 
      y = y_lo;
      transform ( x, y );
      x_min = ( x_min < x ) ? x_min : x;
      x_max = ( x_max > x ) ? x_max : x;
      y_min = ( y_min < y ) ? y_min : y;
      y_max = ( y_max > y ) ? y_max : y;
    }

  // Finding out xmin, xmax, ymin, ymax along line y =  y_hi
  for ( int i = 0; i <= n; i++)
    {
      x = x_lo + i * dx; 
      y = y_hi;
      transform ( x, y );
      x_min = ( x_min < x ) ? x_min : x;
      x_max = ( x_max > x ) ? x_max : x;
      y_min = ( y_min < y ) ? y_min : y;
      y_max = ( y_max > y ) ? y_max : y;
    }

  // Finding out xmin, xmax, ymin, ymax along line x =  x_lo
  for ( int i = 0; i <= n; i++)
    {
      x = x_lo; 
      y = y_lo + i * dy;
      transform ( x, y );
      x_min = ( x_min < x ) ? x_min : x;
      x_max = ( x_max > x ) ? x_max : x;
      y_min = ( y_min < y ) ? y_min : y;
      y_max = ( y_max > y ) ? y_max : y;
    }

  // Finding out xmin, xmax, ymin, ymax along line x =  x_hi
  for ( int i = 0; i <= n; i++)
    {
      x = x_hi; 
      y = y_lo + i * dy;
      transform ( x, y );
      x_min = ( x_min < x ) ? x_min : x;
      x_max = ( x_max > x ) ? x_max : x;
      y_min = ( y_min < y ) ? y_min : y;
      y_max = ( y_max > y ) ? y_max : y;
    }
    
  return HippoRectangle ( x_min, y_min, x_max - x_min, y_max - y_min );
}

/* virtual */
void Lambert::validate ( Range & lat, Range & lon ) const
{
  if ( lat.low () < -180.0 ) lat.setLow ( -180.0 );
  if ( lat.high () > 180.0 ) lat.setHigh ( 180.0 );

  if ( lon.low () < -90.0 ) lon.setLow ( -90.0 );
  if ( lon.high () > 90.0 ) lon.setHigh ( 90.0 );
}


const vector < AxisTick > &
Lambert::
setTicks ( AxisModelBase & model, hippodraw::Axes::Type axis )
{ 
  if ( axis != Axes:: Z ) {
    setTickStep ( model );
    setFirstTick ( model );
    return genTicks ( model, axis );
  }
  else {
    return m_z -> setTicks ( model );
  }
}

void
Lambert::
adjustValues ( AxisModelBase & model,
               hippodraw::Axes::Type,
               const Range & limit )
{
  // Does not do anything as of now 
  return;
}

inline double FLT_EQUAL( double x, double y )
{
  return ( (double)abs( x - y ) <= 2.0 * ( y * FLT_EPSILON + FLT_MIN ) );
}


void Lambert::setTickStep( AxisModelBase & axis )
{
  const Range & range = axis.getRange(false);
  double rangeLength = range.length();

  double scale_factor = axis.getScaleFactor();
  rangeLength *= scale_factor;
  
  // The following algorithm determines the magnitude of the range...
  double rmag = floor( log10( rangeLength ) );
  axis.setRMag( rmag );
  
  double scalelow = range.low() * scale_factor;
  double scalehigh = range.high() * scale_factor;

  // We will also need the largest magnitude for labels.
  double pmag = max( floor( log10( abs ( scalehigh ) ) ), 
                     floor( log10( abs ( scalelow ) ) ) );
  axis.setPMag( pmag );
  
  axis.setTickStep( rangeLength / 4.0 );
}

void Lambert::setFirstTick( AxisModelBase & axis )
{
  const Range & range = axis.getRange(false);
  
  axis.setFirstTick ( range.low() );
}


const vector < AxisTick > &
Lambert::
genTicks( AxisModelBase & axis, hippodraw::Axes::Type axistype )
{ 
  double y = 0.0, ylabel;
  
  int num_ticks = 0;
  m_ticks.clear();
  double pmag = axis.getPMag();
  double rmag = axis.getRMag();
  double first_tick = axis.getFirstTick();
  double tick_step  = axis.getTickStep();
  double scale_factor = axis.getScaleFactor();
  
  // pmag will get set to 0 if it is less than or equal to 3.  This
  // is used later to determine scientific notation.  However, m_rmag
  // is still needed as the original magnitude for calculations such
  // as decimal place notation, and rounding to nice numbers.
  
  bool use_pmag = abs ( pmag ) > 3.0;

  axis.setUsePMag ( use_pmag );
  
  char pstr[10];
  char labl[10];

  int decimals = 0;

  // The following if-else block decides the pstr string, which holds
  // the number of decimal places in the label.

  //   if( fabs( m_pmag ) > 3.0 ) {
  if ( use_pmag ) {  
    // If we are greater than mag 3, we are showing scientific
    // notation.  How many decimals we show is determined by the
    // difference between the range magnitude and the power magnitude.
  
    decimals = static_cast<int>( pmag - rmag );
    // minumum 1 decimal in scientific notation
    
    if( !decimals ) decimals++;
  
  } else {
    
    if( rmag > 0.0 ){
    
      // If we are less than mag 3 and positive, then no decimal
      // accuracy is needed.
      
      decimals = 0;
      
    } else {
    
      // If we are less than mag 3 and negative, then we are suddenly
      // looking at decimal numbers not in scientific notation.
      // Therefore we hold as many decimal places as the magnitude.
      
      decimals = static_cast<int>( abs( rmag ) );
      
    }
  
  }
  // @todo decimals should never be negative here, but it does end up
  //    negative in some cases. See the "dirty fix" in Range.cxx, that
  //    dirty-fixed this problem too. But a better fix is needed. 
  if (decimals < 0) 
    decimals = 0;
   
  sprintf( pstr, "%%1.%df", decimals );

  y = first_tick;
  const Range & range = axis.getRange(false);
  double range_high = range.high();
  range_high *= scale_factor;
  
  while( y <= range_high || FLT_EQUAL( range_high, y ) )
    {
      double value = 0.0;
      
      if ( axistype == Axes::X ) {
        value = moduloAddX( y, xOffset() );
      }
      else if ( axistype == Axes::Y ) {
        value = moduloAddY( y, yOffset() );
      }
      
      // Now we either keep the original magnitude
      // or reduce it in order to express it in scientific notation.
      
      if ( use_pmag ) ylabel = value / pow( 10.0, pmag );
      else ylabel = value;
      
      value /= scale_factor;
      sprintf( labl, pstr, ylabel );
      m_ticks.push_back( AxisTick ( y, labl ) );
      
      num_ticks++;
      y += tick_step;
      
    }
  
  return m_ticks;
}

---