package Geo::Coordinates::WGS84toCH1903;

use strict;
use warnings;

our $VERSION = '0.01';

sub new {
	my $class = shift;	
	
	my $self = {	
		_WGSLatDec   => undef,
		_WGSLongDec  => undef,
		_CH1903X     => undef,        
		_CH1903Y     => undef	
	};

	bless($self, $class);
	return $self;
}

sub attach_WGS84_dec {

	# Read params
	my @params = @_;
	my ($self, $lat, $long) = @params;
	
	# Sanity check
	if (defined($lat) && defined($long) && (scalar(@params) == 3)){
		
		# Process if defined
		$self->{_WGSLatDec}  = $lat;
		$self->{_WGSLongDec} = $long;
	}
	else{
		# Give warning if not defined
		print "Did not attach: attach_WGS84_Dec needs a decimal latitude and longitude as input (2 values)!\n"
	}
}	

sub attach_WGS84_sex {

	# Read params
	my @params = @_;
	my ($self, $latD, $latM, $latS, $longD, $longM, $longS) = @params;
	
	# Sanity check
	if (defined($latD) && defined($latM) && defined($latS) && defined($longD) && defined($longM) && defined($longS) && (scalar(@params) == 7)){
		
		# Process if defined
		$self->{_WGSLatDec}  = SEXtoDEC($latD,$latM,$latS);
		$self->{_WGSLongDec} = SEXtoDEC($longD,$longM,$longS);
	}
	else{
		# Give warning if not defined
		print "Did not attach: attach_WGS84_Sex needs latitude and longitude DMS as input (6 values)!\n"
	}
}

sub attach_CH1903 {

	# Read params
	my @params = @_;
	my ($self, $y, $x) = @params;
	
	# Sanity check
	if (defined($x) && defined($y) && (scalar(@params) == 3)){
		
		# Process if defined
		$self->{_CH1903X}  = $x;
		$self->{_CH1903Y}  = $y;
	}
	else{
		# Give warning if not defined
		print "Did not attach: attach_CH1903 needs a Y and X coordinate!\n"
	}
}

sub fetch_CH1903 {
	
	# Read params
	my $self = shift;
	
	# Perform checks
	if (defined($self->{_CH1903X}) && defined($self->{_CH1903X})){
		
		# Report if already defined
		return($self->{_CH1903Y},$self->{_CH1903X});
	}
	elsif(defined($self->{_WGSLatDec}) && defined($self->{_WGSLongDec})){
	
		# Calculate if angle is defined
		$self->{_CH1903X} = int(WGStoCH1903x($self->{_WGSLatDec},$self->{_WGSLongDec}) + 0.5);
		$self->{_CH1903Y} = int(WGStoCH1903y($self->{_WGSLatDec},$self->{_WGSLongDec}) + 0.5);
		
		# Return
		return($self->{_CH1903Y},$self->{_CH1903X});
	}
	else{
	
		# Give warning if nothing attached
		print "Did not fetch CH1903 data: no data has been attached!\n";
	}
}

sub fetch_WGS84_dec {

	# Read params
	my $self = shift;
	
	# Perform checks
	if (defined($self->{_WGSLatDec}) && defined($self->{_WGSLongDec})){
	
		# Report if already defined
		return($self->{_WGSLatDec},$self->{_WGSLongDec});
	}
	elsif(defined($self->{_CH1903X}) && defined($self->{_CH1903X})){
		
		# Calculate if CH1903 is defined
		$self->{_WGSLatDec} = int((CH1903toWGSlat($self->{_CH1903Y},$self->{_CH1903X}))*10000 + 0.5)/10000;
		$self->{_WGSLongDec}= int((CH1903toWGSlong($self->{_CH1903Y},$self->{_CH1903X}))*10000 + 0.5)/10000;
		
		# Return
		return($self->{_WGSLatDec},$self->{_WGSLongDec});
	}
	else{
	
		# Give warning if nothing attached
		print "Did not fetch WGS84 data: no data has been attached!\n";
	}
}

sub fetch_WGS84_sex {

	# Read params
	my $self = shift;
	
	# Fetch decimal
	my @dec = $self->fetch_WGS84_dec();
	
	# Proceed if fetch ok
	if (scalar(@dec) > 0){
	
		# Convert to sexagesimal
		my @sex1 = DECtoSEX($dec[0]);
		my @sex2 = DECtoSEX($dec[1]);
		
		# Join array
		my @sex;
		foreach(@sex1){
			push(@sex,$_)
		}
		foreach(@sex2){
			push(@sex,$_)
		}

		# Return
		return(@sex);
	}
}

sub unattach {

	# Read params
	my $self = shift;
	
	# Empty object
	$self = {	
		_WGSLatDec   => undef,
		_WGSLongDec  => undef,
		_CH1903X     => undef,        
		_CH1903Y     => undef,	
	}
}

#############
# FUNCTIONS #
#############

###################################
## Converts WGS coordinates to   ##
## CH1903 Y coordinate           ##
##                               ##
## Input:  dd.mmss     [decimal] ##
## Output: CH1903 Y    [integer] ##
##                               ##
###################################

sub WGStoCH1903y {

	# Read params
	my @params = @_;
	my ($lat, $long) = @params;

	# Convert sexagesimal to seconds
	my @lat  = DECtoSEX($lat);
	my @long = DECtoSEX($long);
	$lat  = DEGtoSEC(@lat);
	$long = DEGtoSEC(@long);

	# Set auxiliaries
	my $lat_aux = ($lat - 169028.66)/10000;
	my $long_aux = ($long - 26782.5)/10000;

	# Process Y
	my $y = 600072.37 
	        + (211455.93 * $long_aux)
	        -  (10938.51 * $long_aux * $lat_aux)
	        -  (    0.36 * $long_aux * $lat_aux * $lat_aux)
	        -  (   44.54 * $long_aux * $long_aux*  $long_aux);

	# Return Y coordinate
	return $y;
}

###################################
## Converts WGS coordinates to   ##
## CH1903 X coordinate           ##
##                               ##
## Input:  dd.mmss     [decimal] ##
## Output: CH1903 X    [integer] ##
##                               ##
###################################

sub WGStoCH1903x {

	# Read params
	my @params = @_;
	my ($lat, $long) = @params;

	# Convert sexagesimal to seconds
	my @lat  = DECtoSEX($lat);
	my @long = DECtoSEX($long);
	$lat  = DEGtoSEC(@lat);
	$long = DEGtoSEC(@long);

	# Set auxiliaries
	my $lat_aux = ($lat - 169028.66)/10000;
	my $long_aux = ($long - 26782.5)/10000;

	# Process X
	my $x = 200147.07
	        + (308807.95 * $lat_aux) 
	        + (  3745.25 * $long_aux * $long_aux)
	        + (    76.63 * $lat_aux  * $lat_aux) 
	        - (   194.56 * $long_aux * $long_aux * $lat_aux)
	        + (   119.79 * $lat_aux  * $lat_aux  * $lat_aux);

	# Return X coordinate
	return $x; 
}

###################################
## Converts CH1903 coordinates to##
## WGS decimal latitude          ##
##                               ##
## Input:  CH1903 coor [integer] ##
## Output: dd.mmss     [decimal] ##
##                               ##
###################################

sub CH1903toWGSlat {

	# Read params
	my @params = @_;
	my ($y,$x) = @params;	

	# Set auxiliaries
	my $y_aux = ($y - 600000)/1000000;
	my $x_aux = ($x - 200000)/1000000;
  
	# Process latitude
	my $lat = 16.9023892
		+ (3.238272 * $x_aux)
		- (0.270978 * $y_aux * $y_aux)
		- (0.002528 * $x_aux * $x_aux)
		- (  0.0447 * $y_aux * $y_aux * $x_aux)
		- (  0.0140 * $x_aux * $x_aux * $x_aux);
    
	# Unit 10000" to 1 " and converts seconds to degrees
	$lat = $lat * 100/36;
  
	# Return dd.mmss
	return $lat;
}

###################################
## Converts CH1903 coordinates to##
## WGS decimal longitude         ##
##                               ##
## Input:  CH1903 coor [integer] ##
## Output: dd.mmss     [decimal] ##
##                               ##
###################################

sub CH1903toWGSlong {

	# Read params
	my @params = @_;
	my ($y,$x) = @params;	

	# Set auxiliaries
	my $y_aux = ($y - 600000)/1000000;
	my $x_aux = ($x - 200000)/1000000;
  
	# Process longitude
	my $long = 2.6779094
		+ (4.728982 * $y_aux)
		+ (0.791484 * $y_aux * $x_aux)
		+ (  0.1306 * $y_aux * $x_aux * $x_aux)
		- (  0.0436 * $y_aux * $y_aux * $y_aux);
     
	# Unit 10000" to 1 " and converts seconds to degrees
	$long = $long * 100/36;
   
	# Return dd.mmss
	return $long;
}

###################################
## Converts a sexagesimal (base 60#
## angle to a decimal angle      ##
##                               ##
## Input:  dd mm ss      [array] ##
## Output: dd.mmss     [decimal] ##
##                               ##
###################################

sub SEXtoDEC {

	# Read params
	my @params = @_;
	my ($deg,$min,$sec) = @params;	
	
	# Extract decimal
	my $dec = $deg + $min/60 + $sec/3600;

	# return dd.mmss
	return $dec;
}

###################################
## Converts a decimal angle to a ##
## sexagesimal (base 60) angle   ##
##                               ##
## Input:  dd.mmss     [decimal] ##
## Output: dd mm ss      [array] ##
##                               ##
###################################

sub DECtoSEX {

	# Read params
	my @params = @_;
	my ($angle) = @params;	

	# Extract DMS
	my $deg = int($angle);
	my $min = int(($angle-$deg)*60);
	my $sec = ((($angle-$deg)*60)-$min)*60; 

	# Clean up seconds
	$sec = int($sec + 0.5);			# Round
	if (length($sec)<2){ $sec = "0".$sec};	# Make 2 digit number

	# Return dd mm ss array
	return ($deg,$min,$sec);
}

###################################
## Converts a sexagesimal (base 60#
## angle to seconds              ##
##                               ##
## Input:  dd mm ss      [array] ##
## Output: seconds     [integer] ##
##                               ##
###################################

sub DEGtoSEC {

	# Read params
	my @params = @_;
	my ($deg,$min,$sec) = @params;	

	# Calculate seconds
	my $seconds = $sec + $min*60 + $deg*3600;
	
	# Return seconds
	return $seconds;
}

1;

__END__

=head1 NAME

Geo::Coordinates::WGS84toCH1903 - converts WGS84 coordinate system to/from the Swiss CH1903 coordinate system 

=head1 SYNOPSIS

  use Geo::Coordinates::WGS84toCH1903;

  # Constructor
  my $coordinates = new Geo::Coordinates::WGS84toCH1903;
  
  # Attach decimal WGS84 coordinates to the object
  $coordinates->attach_WGS84_dec(45.966667,7.65);
  
  # Attach sexagesimal WGS84 coordinates to the object
  $coordinates->attach_WGS84_sex(45,58,0,7,39,0);
  
  # Attach CH1903 coordinates to the object
  $coordinates->attach_CH1903(4617049,91670);

  # Return decimal WGS84 coordinates from the object (regardless of the input type)
  my ($lat,$long) = $coordinates->fetch_WGS84_dec();
  
  # Return sexagesimal WGS84 coordinates from the object (regardless of the input type)
  my @DMS_coordinates = $coordinates->fetch_WGS84_sex(); 
  
  # Return CH1903 coordinates from the object (regardless of the input type)
  my ($x,$y) = $coordinates->fetch_CH1903();
  
  # Clear the object
  $coordinates->unattach();
  
=head1 DESCRIPTION

This module converts WGS84 coordinate system to/from the Swiss CH1903 coordinate system. The WGS84 coordinate 
system is what many people know from maps and considered 'the default'. The CH1903 coordinate system covers 
Switzerland and surroundings. The coordinate system and the transformation functions are defined by the 
Official Federal Office of Topography of Switzerland (swisstopo).

The 0 / 0 coordinate is located near Bordeaux, France. This definition invokes that CH1903 coordinates are only
usefull for places, such as Switzerland, located in the 1st quadrant (i.e. NE of this point) of the coordinate
system.

=head1 METHODS

=over 8

=item attach_WGS84_dec(latitude, longitude)

    Attaches a decimal WGS84 coordinate to the object.
    
    Needs a two-dimensional array as input where both latitude and longitude are decimal. 
    As the coordinate system is designed to be used in Switzerland, input only N and E locations.
    
    This method has no return values.
    
=item attach_WGS84_sex(degrees, minutes, seconds, degrees, minutes, seconds)

    Attaches a sexagesimal (base 60) WGS84 coordinate to the object.
    
    Needs a six-dimensional array as input where both latitude and longitude are described by degree, minutes 
    and seconds. E.g. 45°58'35''N / 7°39'31''E would be (45,58,35,7,39,31) as input.
    As the coordinate system is designed to be used in Switzerland, input only N and E locations
    
    This method has no return values.

=item attach_CH1903(y, x)

    Attaches a CH1903 coordinate to the object.
    
    Needs a two-dimensional array input where both coordinates are integers.
    
    This method has no return values.

=item fetch_WGS84_dec()

    This method has no input values.

    Returns a decimal WGS84 coordinate if a decimal/sexagesimal WGS84 coordinate was attached to the object.
    If a CH1903 coordinate was attached to the object, the coordinate is converted to WGS84, stored in the object 
    , and returned.
    
    The return value is a two-dimensional array with latitude and longitude.

=item fetch_WGS84_sex()

    This method has no input values.

    Returns a sexagesimal WGS84 coordinate if a decimal/sexagesimal WGS84 coordinate was attached to the object.
    If a CH1903 coordinate was attached to the object, the coordinate is converted to WGS84, stored in the object 
    , and returned.
    
    The return value is a six-dimensional array with both the latitude and longitude in sexagesimal DMS coordinates.

=item fetch_CH1903()

   This method has no input values.
   
   Return a CH1903 coordinate if a CH1903 coordinate was attached to the object.
   If a decimal/sexagesimal WGS84 coordinate was attached to the object, the coordinate is converted to CH1903,
   stored in the object, and returned.
   
   The return value is a two-dimensional array with both the y and x CH1903 coordinate.


=item unattach()

   This method has no input and no return values.
   
   Empties the object.

=back

=head1 BUGS

Please report to author.

=head1 SEE ALSO

Official Federal Office of Topography swisstopo website
http://www.swisstopo.admin.ch/internet/swisstopo/en/home/products/software/products/skripts.html

Wikipedia page on the CH1903 coordinates
http://en.wikipedia.org/wiki/Swiss_coordinate_system

Online conversion tool for WGS84<>CH1903 coordinats
http://www.swisstopo.admin.ch/internet/swisstopo/en/home/apps/calc/navref.html

=head1 AUTHOR

CPAN package: Joachim De Schrijver, <joachim.deschrijver@gmail.com>

=head1 ACKNOWLEDGEMENTS

The Perl realization is based on the implementations provided by Federal Office of Topography swisstopo. 

The implementations are available on the abovementioned swisstopo website:

=head1 LICENSE AND COPYRIGHT

This library is free software; you can redistribute it and/or modify
it under the same terms as Perl itself, either Perl version 5.10.1 or,
at your option, any later version of Perl 5 you may have available.

Copyright (C) 2012 Joachim De Schrijver

=cut