#!/usr/bin/perl
# Las Animas Co. CO, but only between PLSS correction lines.
# Here we recognize that Range lines are parallel to meridians of longitude,
# and Township lines are parallel to parallels of latitude.
# Therefore we will no longer do any map projections!
#
# Author: Dan Jacobson https://www.jidanni.org/
# Copyright: https://www.gnu.org/licenses/gpl.html
# Created: 2023-11-11T03:22:37+0000
# Last-Updated: 2024-05-13T20:18:10+0000
#     Update #: 99
#

use strict;
use warnings "all";

# OK, let's start from the corner of the
# 7th Guide Meridian West, (R 56/57 W) and
# 6th Standard Parallel South (T 30/31 S).
my %P7GW6PS = ( x => -103.7353009, y => 37.3798999 );

# We must stay between the standard parallels (correction lines).
# Townships: T26-30S. And county extent: R51-67W.

#Populate lookup table
my %adr = ( 0 => { x => ( address("R57W") )[0], y => ( address("T31S") )[0] } );

# It is five townships (30 miles) to the 5th Standard Parallel South.
my %center = %P7GW6PS;
my %range  = ( x => [ -11 * 6 .. 6 * 6 ], y => [ 0 .. 5 * 6 ] );

#use Data::Dumper;print Dumper \%range;exit;
my $increment        = 6;
my $numbers_per_mile = 1000;
my %radius           = ( a => 6378206.4, b => 6356583.8 );    # Clarke 1866

# From analyzing the maps:
sub address {    #uses the "beginning edge" of a township or range.
    my %origin = ( T => -42, R => -70 );
    my $n;
    for (@_) {
        /^(T)(\d+)(S)$/
          || /^(R)(\d+)(W)$/
          || die "\"$_[0]\" should be like T57S, R37W.";
        $n = $2;
        if ( $3 eq "S" || $3 eq "W" ) { $n = -$n; }
        my $adr = -6000 * ( $origin{$1} - $n );
        my $des;
        if    ( $1 eq "T" ) { $des = "N" }
        elsif ( $1 eq "R" ) { $des = "E" }
        else                { die "$1 not T and not R?" }
        return $adr, $des;
    }
}
use Math::Trig;
use Math::Trig ':pi';
sub dlon { 360 * 1609.344 / ( cos deg2rad(pop) ) / $radius{a} / pi2; }
sub dlat { 360 * 1609.344 / $radius{b} / pi2; }

print "WKT,Name\n";
for my $y ( @{ $range{y} } ) {
    for my $x ( @{ $range{x} } ) {
        if ( $y % $increment || $x % $increment ) { next; }
        my %nos;
        my %v = ( x => $x, y => $y );
        for (qw /x y/) {
            $nos{adr}{$_} = $adr{0}{$_} + $numbers_per_mile * $v{$_};
            $nos{cr}{$_}  = ( $_ eq "y" ? -1 : 0 ) + $nos{adr}{$_} / 500;
            die "leftover in CR $_" unless $nos{cr}{$_} = int $nos{cr}{$_};
        }
        printf qq{"POINT (%f %f)",%d/%d\n},
          $center{x} - dlon($y) * $x,
          $center{y} + dlat() * $y,
          $nos{cr}{x}, $nos{cr}{y};
        next;
        printf qq{"POINT (%f %f)",%dE/%d %dN/%d\n},
          $center{x} - dlon($y) * $x,
          $center{y} + dlat() * $y,
          $nos{adr}{x}, $nos{cr}{x}, $nos{adr}{y}, $nos{cr}{y};
    }
}