### ----------------------------------------------------
### -- Math::Trig
### -- Licenses: Artistic-2.0
### -- Authors: Tom Browder <tbrowder@acm.org> [since 2020-02-14], Jonathan Scott Duff <duff@pobox.com> [deceased]
### -- File: lib/Math/Trig.rakumod
### ----------------------------------------------------
unit module Math::Trig;

sub rad2rad ($rad) is export
{
    $rad % tau;
}

sub deg2deg ($deg) is export
{
    $deg % 360;
}

sub grad2grad ($grad) is export
{
    $grad % 400;
}

sub rad2deg ($rad) is export
{
    $rad * 180 / pi;
}

sub deg2rad ($deg) is export
{
    $deg * pi / 180;
}

sub grad2deg ($grad) is export
{
    360 / 400 * $grad;
}

sub deg2grad ($deg) is export
{
    400 / 360 * $deg;
}

sub rad2grad ($rad) is export
{
    400 / tau * $rad;
}

sub grad2rad ($grad) is export
{
    tau / 400 * $grad;
}

sub cartesian-to-spherical($x,$y,$z) is export(:radial)
{
    my $rho = sqrt($x*$x + $y*$y + $z*$z);
    $rho, atan2($y, $x), acos($z / $rho);
}

sub spherical-to-cartesian($rho, $theta, $phi) is export(:radial)
{
    ( $rho * cos( $theta ) * sin( $phi ),
      $rho * sin( $theta ) * sin( $phi ),
      $rho * cos( $phi   ) );
}

sub spherical-to-cylindrical($rho, $theta, $phi) is export(:radial)
{
    my ($x, $y, $z) = spherical-to-cartesian($rho,$theta,$phi);
    ( sqrt( $x * $x + $y * $y ), $theta, $z );
}

sub cartesian-to-cylindrical($x,$y,$z) is export(:radial)
{
    ( sqrt( $x * $x + $y * $y ), atan2( $y, $x ), $z );
}

sub cylindrical-to-cartesian($rho, $theta, $z) is export(:radial)
{
    ( $rho * cos( $theta ), $rho * sin( $theta ), $z );
}

sub cylindrical-to-spherical($rho, $theta, $phi)  is export(:radial)
{
    cartesian-to-spherical( |cylindrical-to-cartesian( $rho, $theta, $phi ) );
}

sub great-circle-distance($theta0, $phi0, $theta1, $phi1, $rho = 1) is export(:great-circle)
{
    my $lat0 = pi/2 - $phi0;
    my $lat1 = pi/2 - $phi1;
    $rho * acos( cos( $lat0 ) * cos( $lat1 ) * cos( $theta0 - $theta1 ) +
                 sin( $lat0 ) * sin( $lat1 ) );
}

sub great-circle-direction($theta0, $phi0, $theta1, $phi1) is export(:great-circle)
{
    my $lat0 = pi/2 - $phi0;
    my $lat1 = pi/2 - $phi1;

    rad2rad(2 * pi -
        atan2(sin($theta0-$theta1) * cos($lat1),
                cos($lat0) * sin($lat1) -
                    sin($lat0) * cos($lat1) * cos($theta0-$theta1)));
}

our &great-circle-bearing is export(:great-circle) = &great-circle-direction;

sub great-circle-waypoint($theta0, $phi0, $theta1, $phi1, $point = 0.5) is export(:great-circle)
{
    my $d = great-circle-distance( $theta0, $phi0, $theta1, $phi1 );

    return if $d == pi;

    my $sd = sin($d);

    return ($theta0, $phi0) if $sd == 0;

    my $A = sin((1 - $point) * $d) / $sd;
    my $B = sin(     $point  * $d) / $sd;

    my $lat0 = pi/2 - $phi0;
    my $lat1 = pi/2 - $phi1;

    my $x = $A * cos($lat0) * cos($theta0) + $B * cos($lat1) * cos($theta1);
    my $y = $A * cos($lat0) * sin($theta0) + $B * cos($lat1) * sin($theta1);
    my $z = $A * sin($lat0)                + $B * sin($lat1);

    my $theta = atan2($y, $x);
    my $phi   = acos($z);

    ($theta, $phi);
}

our &great-circle-midpoint is export(:great-circle) = &great-circle-waypoint.assuming(*,*,*,*,0.5);

sub great-circle-destination( $theta0, $phi0, $dir0, $dst ) is export(:great-circle)
{
    my $lat0 = pi/2 - $phi0;

    my $phi1   = asin(sin($lat0)*cos($dst) +
                      cos($lat0)*sin($dst)*cos($dir0));

    my $theta1 = $theta0 + atan2(sin($dir0)*sin($dst)*cos($lat0),
                                 cos($dst)-sin($lat0)*sin($phi1));

    my $dir1 = great-circle-bearing($theta1, $phi1, $theta0, $phi0) + pi;

    $dir1 -= 2*pi if $dir1 > 2*pi;

    ($theta1, $phi1, $dir1);
}

# Trigonemetric functions with arguments in degrees:
#   rad = deg × π / 180°
#   deg = rad × 180° / π
constant $deg2rad is export(:deg2rad :TRIG :ALL) = pi / 180.0;
constant $rad2deg is export(:rad2deg :TRIG :ALL) = 180.0 / pi;
constant deg2rad is export(:deg2rad :TRIG :ALL) = pi / 180.0;
constant rad2deg is export(:rad2deg :TRIG :ALL) = 180.0 / pi;

constant pi2 is export(:pi2 :TRIG :ALL) = pi/2;
constant pi4 is export(:pi2 :TRIG :ALL) = pi/4;
constant _2pi is export(:_2pi :TRIG :ALL) = pi*2;
constant _4pi is export(:_4pi :TRIG :ALL) = pi*4;

# arguments in degrees
sub sind($degrees) is export(:sind :TRIG :ALL) {
    return sin $degrees * $deg2rad;
}
sub cosd($degrees) is export(:cosd :TRIG :ALL) {
    return cos $degrees * $deg2rad;
}
sub tand($degrees) is export(:tand :TRIG :ALL) {
    return tan $degrees * $deg2rad;
}

# return results in degrees
sub asind($x) is export(:asind :TRIG :ALL) {
    return $rad2deg * asin $x;
}
sub acosd($x) is export(:acosd :TRIG :ALL) {
    return $rad2deg * acos $x;
}
sub atand($x) is export(:atand :TRIG :ALL) {
    return $rad2deg * atan $x;
}
sub atan2d($y, $x) is export(:atan2d :TRIG :ALL) {
    return $rad2deg * atan2 $y, $x;
}


