### ----------------------------------------------------
### -- Geo::Ellipsoid
### -- Licenses: Artistic-2.0
### -- Authors: Tom Browder <tbrowder@acm.org> (translator Perl 5 to Raku, current maintainer), Jim Gibson (CPAN Perl 5 author, v1.12)
### -- File: lib/Geo/Ellipsoid.rakumod
### ----------------------------------------------------
### -- Chunk 2 of 2
 $cz * $d/6.0) - $x) *
    $d /4.0) + $cz) * $sy * $d) + $y) * $c * $a * $r;

    printf "s=%.8f\n", $s if $DEBUG;

    # REPLACE THIS WITH FUNCTION CALL!!
    # adjust azimuth to (0,360) or (-180,180) as specified
    # units MUST be radians at this point
    if self.bearing_sym {
        $faz += $twopi if $faz < -(pi);
        $faz -= $twopi if $faz >= pi;
    } 
    else {
        $faz += $twopi if $faz < 0;
        $faz -= $twopi if $faz >= $twopi;
    }

    if self.units eq 'degrees' {
        $faz = rad2deg $faz;
    } 

    # return result
    #my @disp = (($s/self.conversion), $faz);
    my @disp = ($s/self.conversion, $faz);
    print "disp = (@disp)\n" if $DEBUG;
    return (|@disp);
} # _inverse

#	_forward
#
#	Calculate the location (latitude, longitude) of a point
#	given a starting point and a displacement from that
#	point as (range, bearing) where range is in the class's
#       current units and bearing is in degrees from true north.
#
# pseudo "private" method
method _forward($lat1 is copy, $lon1 is copy, $range, $bearing is copy)
{
  #die "FATAL: need to use radians for trig funcs, units are: {self.units}" if self.units eq 'degrees';
  if self.units eq 'degrees' {
      $lat1    = deg2rad $lat1;
      $lon1    = deg2rad $lon1;
      $bearing = deg2rad $bearing;
  }

  if $DEBUG {
    printf "_forward(lat1=%.8f,lon1=%.8f,range=%.8f,bearing=%.8f)\n",
      $lat1, $lon1, $range, $bearing;
  }

  my $eps = 0.5e-13;

  my $a = self.equatorial;
  my $f = self.flattening;
  my $r = 1.0 - $f;
  my $tu = $r * sin($lat1) / cos($lat1);
  my $faz = $bearing;
  my $s = self.conversion * $range;

  if $DEBUG {
      say "DEBUG: \$faz = $faz";
      say $faz.WHAT;
  }

  my $sf = sin($faz);
  my $cf = cos($faz);

  my $baz = 0.0;
  $baz = 2.0 * atan2($tu,$cf) if $cf != 0.0;

  my $cu = 1.0 / sqrt(1.0 + $tu*$tu);
  my $su = $tu * $cu;
  my $sa = $cu * $sf;
  my $c2a = 1.0 - ($sa*$sa);
  my $x = 1.0 + sqrt((((1.0/($r*$r)) - 1.0)*$c2a) +1.0);
  $x = ($x-2.0)/$x;
  my $c = 1.0 - $x;
  $c = ((($x*$x)/4.0) + 1.0)/$c;
  my $d = $x * ((0.375*$x*$x)-1.0);
  $tu = (($s/$r)/$a)/$c;
  my $y = $tu;

  if $DEBUG {
    printf "r=%.8f, tu=%.8f, faz=%.8f\n", $r, $tu, $faz;
    printf "baz=%.8f, sf=%.8f, cf=%.8f\n", $baz, $sf, $cf;
    printf "cu=%.8f, su=%.8f, sa=%.8f\n", $cu, $su, $sa;
    printf "x=%.8f, c=%.8f, y=%.8f\n", $x, $c, $y;
  }

  my ($cy, $cz, $e, $sy);
  repeat {
    $sy = sin($y);
    $cy = cos($y);
    $cz = cos($baz+$y);
    $e = (2.0*$cz*$cz)-1.0;
    $c = $y;
    $x = $e * $cy;
    $y = (2.0 * $e) - 1.0;
    $y = ((((((((($sy*$sy*4.0)-3.0)*$y*$cz*$d)/6.0)+$x)*$d)/4.0)-$cz)*$sy*$d) +
      $tu;
    } while (abs($y-$c) > $eps);

  $baz = ($cu*$cy*$cf) - ($su*$sy);
  $c = $r*sqrt(($sa*$sa) + ($baz*$baz));
  $d = $su*$cy + $cu*$sy*$cf;
  my $lat2 = atan2($d,$c);
  $c = $cu*$cy - $su*$sy*$cf;
  $x = atan2($sy*$sf,$c);
  $c = (((((-3.0*$c2a)+4.0)*$f)+4.0)*$c2a*$f)/16.0;
  $d = (((($e*$cy*$c) + $cz)*$sy*$c)+$y)*$sa;
  my $lon2 = $lon1 + $x - (1.0-$c)*$d*$f;
  #$baz = atan2($sa,$baz) + pi;

  if self.units eq 'degrees' {
      $lat2 = rad2deg $lat2;
      $lon2 = rad2deg $lon2;
  }

  # return result
  return ($lat2, $lon2);

} # _forward

=begin pod

=head1 DEFINED ELLIPSOIDS

The following ellipsoids are defined in Geo::Ellipsoid, with the
semi-major axis in meters and the reciprocal flattening as shown.
The default ellipsoid is WGS84.

    Ellipsoid        Semi-Major Axis (m.)     1/Flattening
    ---------        -------------------     ---------------
    AIRY                 6377563.396         299.3249646
    AIRY-MODIFIED        6377340.189         299.3249646
    AUSTRALIAN           6378160.0           298.25
    BESSEL-1841          6377397.155         299.1528128
    CLARKE-1880          6378249.145         293.465
    EVEREST-1830         6377276.345         290.8017
    EVEREST-MODIFIED     6377304.063         290.8017
    FISHER-1960          6378166.0           298.3
    FISHER-1968          6378150.0           298.3
    GRS80                6378137.0           298.25722210088
    HOUGH-1956           6378270.0           297.0
    HAYFORD              6378388.0           297.0
    IAU76                6378140.0           298.257
    KRASSOVSKY-1938      6378245.0           298.3
    NAD27                6378206.4           294.9786982138
    NWL-9D               6378145.0           298.25
    SOUTHAMERICAN-1969   6378160.0           298.25
    SOVIET-1985          6378136.0           298.257
    WGS72                6378135.0           298.26
    WGS84                6378137.0           298.257223563

=head1 LIMITATIONS

The methods should not be used on points which are too near the poles
(above or below 89 degrees), and should not be used on points which
are antipodal, i.e., exactly on opposite sides of the ellipsoid. The
methods will not return valid results in these cases.

=head1 ACKNOWLEDGEMENTS

The conversion algorithms used here are Perl translations of Fortran
routines written by LCDR S<L. Pfeifer> NGS Rockville MD that implement
S<T. Vincenty's> Modified Rainsford's method with Helmert's elliptical
terms as published in "Direct and Inverse Solutions of Ellipsoid on
the Ellipsoid with Application of Nested Equations", S<T. Vincenty,>
Survey Review, April 1975.

The Fortran source code files inverse.for and forward.for
may be obtained from

    ftp://ftp.ngs.noaa.gov/pub/pcsoft/for_inv.3d/source/

=head1 AUTHOR

Jim Gibson, C<< <Jim@Gibson.org> >>
Tom Browder, C<< <tom.browder@gmail.com> >>

=head1 BUGS

See LIMITATIONS, above.

Please enter any bugs or feature requests to the issues 
 through the web interface at
L<https://github.com/tbrowder/Geo-Ellipsoid-Perl6>.

=head1 COPYRIGHT & LICENSE

Copyright (c) 2005-2008 Jim Gibson, all rights reserved.
Copyright (c) 2015-2017 Tom Browder, all rights reserved.

This program is free software; you can redistribute it and/or modify it
under the same terms as Perl itself.

=head1 SEE ALSO

Geo::Distance, Geo::Ellipsoids

=end pod


