use strict;
use warnings;

use Math::Trig;


my $my_debug  = 0;

#we need set the following parameters

#metric meter
my $ellipse_a = 6378137;
my $ellipse_inverse_f = 298.2572221010002;

my $easting = 629980.000;
my $northing = 228527.250;
my $origin_easting = 629980.000;
my $origin_northing = 228527.250;
my $size_easting = 526;
my $size_northing = 474;
#my $size_easting = 5;
#my $size_northing = 4;
my $step_easting = 28.500000000000000;
my $step_northing = -28.500000000000000;

my $latitude_false_orign = deg2rad(33.75);
my $longitude_false_orign = deg2rad(-79);

my $latitude_first_parallel = deg2rad(34.33333333333334);
my $latitude_second_parallel = deg2rad(36.16666666666666);

#metric foot
my $easting_false = 609601.22;
my $northing_false = 0;

my @latitude_array;
my @longitude_array;



#END

#$my_debug and print cos(pi/2) . "\n";

my $ellipse_f = 1 / $ellipse_inverse_f;
my $ellipse_e = sqrt(2 * $ellipse_f - $ellipse_f * $ellipse_f);

$my_debug and print "ellipse_e ". $ellipse_e . "\n";

my $latitude_first_parallel_sin = sin($latitude_first_parallel);
my $latitude_second_parallel_sin = sin($latitude_first_parallel);
my $latitude_false_orign_sin = sin($latitude_false_orign);

my $m1 = cos($latitude_first_parallel) / sqrt( 1 - ($ellipse_e **2) * ($latitude_first_parallel_sin **2));
my $m2 = cos($latitude_second_parallel) / sqrt( 1 - ($ellipse_e **2) * ($latitude_second_parallel_sin **2));

$my_debug and print "m1 : " . $m1 . "\n";
$my_debug and print "m2 : " . $m2 . "\n";

my $t1 = tan(pi/4-$latitude_first_parallel/2) / 
    (
      (
        ( 1 - $ellipse_e * $latitude_first_parallel_sin ) / 
        ( 1 + $ellipse_e * $latitude_first_parallel_sin )  
      ) **     
      ($ellipse_e/2) 
    );
    
my $t2 = tan(pi/4-$latitude_second_parallel/2) / 
    (
      (
        ( 1 - $ellipse_e * $latitude_second_parallel_sin ) / 
        ( 1 + $ellipse_e * $latitude_second_parallel_sin )  
      ) **     
      ($ellipse_e/2) 
    );
    
my $tf = tan(pi/4-$latitude_false_orign/2) / 
    (
      (
        ( 1 - $ellipse_e * $latitude_false_orign_sin ) / 
        ( 1 + $ellipse_e * $latitude_false_orign_sin )  
      ) **     
      ($ellipse_e/2) 
    );
    
$my_debug and print "t1 : ". $t1 . "\n";
$my_debug and print "t2 : ". $t2 . "\n";
$my_debug and print "tf : ". $tf . "\n";

my $n  = ( log($m1) - log($m2) )/ (log($t1) - log($t2));
$my_debug and print "n: ". $n . "\n";

my $F = $m1 / ($n * ($t1**$n)) ;
$my_debug and print "F : ". $F . "\n";


# declare the cooridnates to converted to longitude and latitude


my $r_false= $ellipse_a * $F * ( $tf ** $n);

$my_debug and print "r_false : ". $r_false . "\n";



#declare the subroutine

sub latitude_cal
{
   my ($r_prime,$t_prime,$theta_prime) = @_;
   my $latitude_rad = pi/2 - 2 * atan($t_prime);
   
   # iterata 4 times
   for(my $count =4 ; $count>0 ; $count--){
   $latitude_rad = pi/2 - 2 * atan($t_prime * (
                 (
                  (1-$ellipse_e * sin($latitude_rad)) / (1+$ellipse_e * sin($latitude_rad))
                 )**($ellipse_e/2)));
   }
   my $latitude_deg=rad2deg($latitude_rad);
  
}

sub longitude_cal
{
  my ($r_prime,$t_prime,$theta_prime) = @_;
  my $longitude_rad = $theta_prime / $n + $longitude_false_orign;
  my $longitude_deg=rad2deg($longitude_rad);
}



#$my_debug and print "latitude : ". latitude_cal . "\n";
#$my_debug and print "longitude : ". longitude_cal . "\n";



$easting = $origin_easting;

for(my $count_easting=0;$count_easting<$size_easting+1;$count_easting++)
{
  $northing =  $origin_northing;
  
  for(my $count_northing=0;$count_northing<$size_northing+1;$count_northing++)
  { 
      
      
      my $r_prime = sqrt(($easting-$easting_false)**2+($r_false-$northing+$northing_false)**2);

      if( $n < 0 )
      {
         $r_prime = -$r_prime;
      }

      $my_debug and print "r_prime : ". $r_prime . "\n";

      my $t_prime = ( $r_prime / ($ellipse_a * $F)) ** (1 / $n);
      $my_debug and print "t_prime : ". $t_prime . "\n";
      
      my $theta_prime = atan(($easting-$easting_false)/($r_false-$northing+$northing_false));
      $my_debug and print "theta_prime ". $theta_prime . "\n";
      
      $latitude_array[$count_northing][$count_easting]=latitude_cal($r_prime,$t_prime,$theta_prime);
      $longitude_array[$count_northing][$count_easting]=longitude_cal($r_prime,$t_prime,$theta_prime);
      
      $my_debug and print "$count_easting,$count_northing \n";
      
      $northing += $step_northing;
   }
   
   $easting += $step_easting;
 }
 
 my $data_file="geo.info";
 open(my $IN,">$data_file") or die("cannot open file $data_file \n");
 

 for(my $count_northing=0;$count_northing<$size_northing+1;$count_northing++)
{ 
   for(my $count_easting=0;$count_easting<$size_easting+1;$count_easting++)
  { 
      print $IN "( $count_northing, $count_easting  ) : ";
      print $IN $longitude_array[$count_northing][$count_easting] . " " ;
      print $IN $latitude_array[$count_northing][$count_easting] . " \n" ;
   }
 }
 
 close($IN);