#! perl -slw
no warnings 'pack';
use strict;
use Data::Dump qw[ pp ];
use GD;
use constant FNAME => 'CGScrew';
use constant {
    COS30   => 0.86602540378443864676372317075294,
    TAN30   => 0.57735026918962576450914878050196,
    DEG2RAD => 0.017453292519943295769236907684886,
    RAD2DEG => 57.295779513082320876798154814105,
};
use enum qw[ X Y Z A ];
sub rgb2n{ unpack 'N', pack 'CCCC', 0, @_ }

## Construct thread profile segment
sub constructThreadProfile {
    my @points;
    my( $dia, $pitch, $yInc ) = @_;
    my $H =  $pitch * COS30;
    # draw flat crest
    my $x = $dia / 2;
    my $y = -$yInc;
    push @points, [ $x, $y += $yInc, 0, 1 ] while $y < ( $pitch / 16 );

    ## draw upper 30deg flank.
    my $xd = $yInc / TAN30;
    my $yLim = $y + 5 / 16 * $pitch;
    push @points, [ $x -= $xd, $y += $yInc, 0, 0.5 ] while $y < $yLim;

    ## draw root radius
    $yLim = $y + $pitch / 4;
    ## cx = $dia /2 - 7/8*$H +$H/3
    my( $cx, $cy, $r ) = ( ( $dia/2 - 7/8*$H + $H/3 ), $pitch / 2, $H / 6 );
    while( $y < $yLim ) {
        my $dy = $cy - $y;
        my $dx = sqrt( $r**2 - $dy**2 );
        push @points, [ $cx - $dx, $cy - $dy, 0,  $dx / $r ];
        $y += $yInc;
    }
    $y -= $yInc;
    ## draw lower 30deg flank
    $yLim = $y + 5 / 16 * $pitch;
    push @points, [ $x += $xd, $y +=  $yInc, 0, - 0.5 ] while $y < $yLim;

    push @points, [ $x, $y += $yInc, 0, 1 ] while $y < $pitch;
    return \@points;
}

our $M //= 10;
our $P //= 1.5;
our $L //= 2;
our $S //= 100;
my $fname = sprintf "%sM%.2fxP%.2fxL%.2fxS%d.png", FNAME, $M, $P, $L, $S;

my $profile = constructThreadProfile( $M, $P, 1 / ( 10* $S ) );
#pp $profile;

my( $w, $h ) = ( $M * $S + 200, int( $L / $P + 1 ) * $P * $S + 200 );
my $xc = $w / 2;

my $im =  GD::Image->new( $w, $h, 1 );
$im->fill( 0,0, rgb2n( 128, 128, 128 ) );

sub xformPoint {
    my( $point, $rot, $yTrans ) = @_;
    $rot *= DEG2RAD;
    my $x = $point->[X] * cos( $rot );
    my $y = $point->[Y] + $yTrans;
    my $z = $point->[X] * sin( $rot );
    my $a = $point->[A] * cos( 90 * DEG2RAD - $rot );
    return [ $x, $y, $z, $a ];
}

my $yTrans = $P / 360;
my $maxRad = $M * $S / 2;
my $yOff = 100;
for my $turn ( 1 .. $L / $P ) {
    for my $p ( map $_/2, -$maxRad*2 .. $maxRad*2 ) {
        my $rot = RAD2DEG * atan2( sqrt( $maxRad**2 - $p**2 ), $p );
        for my $point ( @$profile ) {
            my $newPoint = xformPoint( $point, $rot, $yTrans * $rot );
            my( $newX, $newY ) = ( $xc + $newPoint->[X] * $S, $yOff + $newPoint->[Y] * $S );
            my $color = ( abs( $newPoint->[A] ) * 256 + (100 - $newY) );
            $im->setPixel( $newX, $newY, rgb2n( ( $color ) x 3 ) );
        }
    }
    $yOff += $P * $S;
}

$im->string( gdSmallFont, 0,0, $fname, 0 );
open O, '>:raw',  $fname or die $!;
print O $im->png;
close O;

system $fname;