#! perl -sw
use strict;
use PDL;
use PDL::Constants qw(PI DEGRAD);
use PDL::IO::GD;
use constant FNAME => 'CGScrewPDL';
use constant {
    COS30   => 0.86602540378443864676372317075294,
    TAN30   => 0.57735026918962576450914878050196,
};
use enum qw[ X Y A ];
$|=1;

## Construct thread profile segment
sub constructThreadProfile {
    my ($dia, $pitch, $yInc) = @_;
    my $H =  $pitch * COS30;
    # draw flat crest
    my $x = $dia / 2;
    my @yLim = (0, $pitch / 16);
    my $yQuant = ($yLim[-1]-$yLim[-2]) / $yInc;
    my $ySeq = zeroes($yQuant)->xlinvals(0, $yLim[-1]);
    my $topcrest = cat(pdl($x), $ySeq, pdl(1))->mv(-1,0);

    ## draw upper 30deg flank.
    my $xd = $yInc / TAN30;
    push @yLim, $yLim[-1] + 5 / 16 * $pitch;
    $yQuant = ($yLim[-1]-$yLim[-2]) / $yInc;
    my $xSeq = zeroes($yQuant)->xlinvals($x, $x - $xd*$yQuant);
    $ySeq = zeroes($yQuant)->xlinvals($yLim[-2], $yLim[-1]);
    my $topflank = cat($xSeq, $ySeq, pdl(0.5))->mv(-1,0);

    ## draw root radius
    push @yLim, $yLim[-1] + 4 / 16 * $pitch;
    my ($cx, $cy, $r) = ( ( $dia/2 - 7/8*$H + $H/3 ), $pitch / 2, $H / 6 );
    $yQuant = ($yLim[-1]-$yLim[-2]) / $yInc;
    $ySeq = zeroes($yQuant)->xlinvals($yLim[-2], $yLim[-1]);
    my $dy = $cy - $ySeq;
    my $dx = sqrt( $r**2 - $dy**2 );
    my $rootradius = cat($xSeq = $cx - $dx, $ySeq, $dx / $r)->mv(-1,0);
    $x = $xSeq->at(-1);

    ## draw lower 30deg flank
    push @yLim, $yLim[-1] + 5 / 16 * $pitch;
    $yQuant = ($yLim[-1]-$yLim[-2]) / $yInc;
    $xSeq = zeroes($yQuant)->xlinvals($x, $x + $xd*$yQuant);
    $ySeq = zeroes($yQuant)->xlinvals($yLim[-2], $yLim[-1]);
    my $botflank = cat($xSeq, $ySeq, pdl(-0.5))->mv(-1,0);
    $x = $xSeq->at(-1);

    ## the bottom crest
    push @yLim, $pitch;
    $yQuant = ($yLim[-1]-$yLim[-2]) / $yInc;
    $ySeq = zeroes($yQuant)->xlinvals($yLim[-2], $yLim[-1]);
    my $botcrest = cat(pdl($x), $ySeq, pdl(1))->mv(-1,0);
    $topcrest->glue(1, $topflank, $rootradius, $botflank, $botcrest);
}

our $M //= 10;
our $P //= 1.5;
our $L //= 2;
our $S //= 100;
my $profilePDL = constructThreadProfile( $M, $P, 1 / ( 10* $S ) );
my ($w, $h) = ( $M * $S + 200, int( $L / $P + 1 ) * $P * $S + 200 );
my $imPDL = zeroes(byte, $w, $h) + 128; # no RGB as all same, dummy at end

my $xc = $w / 2;
my $yTrans = $P / 360;
my $maxRad = $M * $S / 2;
my $yOffs = 100 + sequence($L / $P) * $P * $S;
for my $yOff ( $yOffs->list ) {
print "*";
    my $p = zeroes(4*$maxRad+1)->xlinvals(-$maxRad, $maxRad); # 4mR
    my $rot = atan2(sqrt($maxRad**2 - $p**2), $p)->dummy(0);  # 1 4mR
    my $this_profile = $profilePDL                            # 3 pE
      + cat(pdl(0), $yTrans*DEGRAD*$rot, pdl(0))->mv(-1,0)    # 3 1  4mR
      ;                                                       # 3 pE 4mR
    $this_profile *= cat(cos($rot), pdl(1), cos( PI/2 - $rot))->mv(-1,0);
    my $newX = $xc + $this_profile->slice([X,0,0]) * $S;
    my $newY = $yOff + $this_profile->slice([Y,0,0]) * $S;
    my $color = $this_profile->slice([A,0,0])->abs * 256 + (100 - $newY);
    $imPDL->indexND( $newX->cat($newY)->mv(-1,0) ) .= $color;
}
print "\n";

my $fname = sprintf "%sM%.2fxP%.2fxL%.2fxS%d.png", FNAME, $M, $P, $L, $S;
my $im = PDL::IO::GD->new({pdl=> $imPDL->dummy(2,3)});
$im->String( gdFontGetSmall, 0,0, $fname, 0 );
$im->write_Png($fname);