# -*- coding: utf-8 -*-
#
# e-PRTR XML-3 Export
# (c) 2009, 2010 Matthias Lüttgert, ENDA GmbH & Co. KG
# (c) 2009 Torsten Lüdtke, RISA GmbH
#
# This program is free software; you can redistribute it and/or modify it under
# the terms of the GNU General Public License as published by the Free Software
# Foundation; either version 3 of the License, or (at your option) any later
# version.
#
# This program is distributed in the hope that it will be useful, but WITHOUT ANY
# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
# PARTICULAR PURPOSE. See the GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License along with
# this program; if not, see <http://www.gnu.org/licenses/>.
#
# DEUTSCH:
#
# Dieses Programm ist freie Software. Sie können es unter den Bedingungen der GNU
# General Public License, wie von der Free Software Foundation veröffentlicht,
# weitergeben und/oder modifizieren, entweder gemäß Version 3 der Lizenz oder
# (nach Ihrer Option) jeder späteren Version.
#
# Die Veröffentlichung dieses Programms erfolgt in der Hoffnung, daß es Ihnen von
# Nutzen sein wird, aber OHNE IRGENDEINE GARANTIE, sogar ohne die implizite
# Garantie der MARKTREIFE oder der VERWENDBARKEIT FÜR EINEN BESTIMMTEN
# ZWECK. Details finden Sie in der GNU General Public License.
#
# Sie sollten ein Exemplar der GNU General Public License zusammen mit diesem
# Programm erhalten haben. Falls nicht, siehe <http://www.gnu.org/licenses/>.

  class Gausskrueger2Georef
  include Math
  
  def initialize (rechtswert, hochwert)
    rho  = 180 / PI
    e2   = 0.0067192188
    c    = 6398786.849
    b2   = (hochwert / 10000855.7646) * (hochwert / 10000855.7646)
    bf   = 325632.08677 * (hochwert / 10000855.7646) * ((((((0.00000562025 * b2 + 0.00022976983) * b2 - 0.00113566119) * b2 + 0.00424914906) * b2 - 0.00831729565) * b2 + 1))
    bf   /= 3600 * rho
    co   = Math.cos(bf)
    g2   = e2 * (co * co)
    g1   = c / Math.sqrt(1 + g2)
    t    = Math.tan(bf)
    fa   = (rechtswert - (rechtswert / 1000000).floor * 1000000 - 500000) / g1
    geor = ((bf - fa * fa * t * (1 + g2) / 2 + fa * fa * fa * fa * t * (5 + 3 * t * t + 6 * g2 - 6 * g2 * t * t) / 24) * rho)
    geoh = fa - fa * fa * fa * (1 + 2 * t * t + g2) / 6 + fa * fa * fa * fa * fa * (1 + 28 * t * t + 24 * t * t * t * t) / 120
    geoh = geoh * rho / co + (rechtswert / 1000000).floor * 3

    abessel = 6377397.155
    eebessel = 0.0066743722296294277832
    scalefactor = 0.00000982
    rotxrad = -7.16069806998785E-06
    rotyrad = 3.56822869296619E-07
    rotzrad = 7.06858347057704E-06
    #shiftxmeters = 591.28
    #shiftymeters = 81.35
    #shiftzmeters = 396.39
    shiftxmeters = 598.1 # nach WIKIPEDIA Transfortmationskorrektur für DHDN/Potsdam 2001
    shiftymeters = 73.7  # nach WIKIPEDIA Transfortmationskorrektur für DHDN/Potsdam 2001
    shiftzmeters = 418.2 # nach WIKIPEDIA Transfortmationskorrektur für DHDN/Potsdam 2001
    awgs84 = 6378137
    eewgs84 = 0.0066943799

    geor = (geor / 180) * PI
	  geoh = (geoh / 180) * PI

  	n = eebessel * Math.sin(geor) * Math.sin(geor)
    n = 1 - n
    n = Math.sqrt(n)
    n = abessel / n

    cartesianxmeters = n * Math.cos(geor) * Math.cos(geoh)
    cartesianymeters = n * Math.cos(geor) * Math.sin(geoh)
	  cartesianzmeters = n * (1 - eebessel) * Math.sin(geor)

    cartoutputxmeters = (1 + scalefactor) * cartesianxmeters + rotzrad * cartesianymeters - rotyrad * cartesianzmeters + shiftxmeters
    cartoutputymeters = -rotzrad * cartesianxmeters + (1 + scalefactor) * cartesianymeters + rotxrad * cartesianzmeters + shiftymeters
    cartoutputzmeters = rotyrad * cartesianxmeters - rotxrad * cartesianymeters + (1 + scalefactor) * cartesianzmeters + shiftzmeters

    geoh = Math.atan((cartoutputymeters / cartoutputxmeters))

    latitude   = (cartoutputxmeters * cartoutputxmeters) + (cartoutputymeters * cartoutputymeters)
    latitude   = Math.sqrt(latitude)
    latitude   = cartoutputzmeters / latitude
    latitude   = Math.atan(latitude)
    latitudeit = 99999999

    loop do
    	latitudeit = latitude

      n        = 1 - eewgs84 * Math.sin(latitude) * Math.sin(latitude)
      n        = Math.sqrt(n)
      n        = awgs84 / n
      latitude = cartoutputxmeters * cartoutputxmeters + cartoutputymeters * cartoutputymeters
      latitude = Math.sqrt(latitude)
      latitude = (cartoutputzmeters + eewgs84 * n * Math.sin(latitudeit)) / latitude
      latitude = Math.atan(latitude)

      break unless (latitude - latitudeit).abs >= 0.000000000000001
    end

    @breitengrad = (latitude / PI) * 180
    @laengengrad = (geoh / PI) * 180
  end

  def breitengrad
    return @breitengrad
  end

  def laengengrad
    return @laengengrad
  end
end
