Wikel Tripel Map Projection in Ruby

MapProjection_WinkelTripel.rb

# Usage
# projection = MapProjection::WinkelTripel.new(w, h)
# projections.project( latitude, longitude)
#   => [ x, y ] # cartesian coordinates
# Where x and y will be between 0, 0 (upper-left) and w, h (lower-right)
#   suitable for compositing images on a bitmap

module Math
  DEG2RAD = PI / 180
  
  # Provide sinc function to Math: http://en.wikipedia.org/wiki/Sinc_function
  def self.sinc x
    return 1 if 0 == x
    sin(x) / x
  end
  
  def self.deg_to_rad degrees
    degrees * DEG2RAD
  end
end

# Modified implementation of Winkel Tripel projection, lat, lng -> cartesian cooridinates
# Returns cartesian coordinates scaled to the dimensions you provide
module MapProjection
  NORTH_POLE   = [ 90.0, 0.0 ]
  SOUTH_POLE   = [ -90.0, 0.0 ]
  COTE_DIVOIRE = [ 0.0, 0.0 ]
  EAST_EDGE    = [ 0, 180 ]
  PORTLAND     = [ 45.5184, -122.6554 ]
  
  class WinkelTripel
    attr_reader :width, :height
    def initialize width, height
      @width, @height = width, height
      @x_scale = width / WIDTH
      @y_scale = height / HEIGHT
    end
    
    def project lat, lng
      x, y = self.class.unscaled_project(lat, lng)
      [ (x + X_OFFSET) * @x_scale, height - (y + Y_OFFSET) * @y_scale ]
    end
    
    private
      def self.unscaled_project lat, lng
        lat = Math.deg_to_rad(lat)
        lng = Math.deg_to_rad(lng)
        
        alpha = Math.acos(Math.cos(lat) * Math.cos(lng / 2))
        
        x = (lng * COS_PHI_SUB_ONE + (2 * Math.cos(lat) *
            Math.sin(lng / 2)) / Math.sinc(alpha)) / 2
        y = (lat + (Math.sin(lat) / Math.sinc(alpha))) / 2
        
        [ x, y ]
      end
    
      PHI_SUB_ONE     = Math.acos(2 / Math::PI)
      COS_PHI_SUB_ONE = 2 / Math::PI # Math.cos(PHI_SUB_ONE)
      
      X_OFFSET        = unscaled_project(*EAST_EDGE)[0]
      Y_OFFSET        = unscaled_project(*NORTH_POLE)[1]
      
      WIDTH           = X_OFFSET * 2
      HEIGHT          = Y_OFFSET * 2
  end
end