#!/usr/bin/perl -w

use strict;
use POSIX;

# The vertical tube stands upright. and is 35x35 with rounded edges of a diameter of:
# diagonal measurment gives me 46.8mm.
my $edge_35 = 35;
# this would be 3.256mm radius
my $edge_radius_35 = (sqrt(2) * $edge_35 - 46.8) / (2 *(sqrt(2) - 1));

# for first trike (snow-safety) 30x30x2, diagonal is 40
my $edge_30 = 30 * (116/115);
my $edge_radius_30 = (sqrt(2) * $edge_30 - 40*(116/115)) / (2 *(sqrt(2) - 1));

# The parameters of the 10x30
my $h_10x30 = 30;
my $w_10x30 = 10;
my $edge_radius_10x30 = 1.47; # estimation  + calibration

# we need a trace along the xy plane for the vertical tube.
# we use an angle rotating in the xy plane with stepsize radius
my $stepsize = 0.005;

my $PI = 3.1415926535;

my $a45 = $PI / 4;
my $a90 = $PI / 2;
my $a135 = 3 * $PI / 4;
my $a180 = $PI;
my $a225 = 5 * $PI / 4;
my $a270 = 6 * $PI / 4;
my $a315 = 7 * $PI / 4;
my $a360 = 2 * $PI;

print "%!\n";

print "0.5 setlinewidth\n";
draw_section_kantelknik_voorbouw();
#draw_section_tube_tube(51.6*$PI/180, 1.5, $edge_35, $edge_35, $edge_radius_35);
#draw_section_tube_tube(9.9*$PI/180, 1.5, $w_10x30, $h_10x30, $edge_radius_10x30);
#draw_section_tretlager();

#draw_section_cilinder_cilinder(0, 28.6/2, 28.6/2);

#draw_section_tube_tube(90*$PI/180, 2, $edge_30, $edge_30, $edge_radius_30);

print "showpage\n";

# the horizontal tube is 45mm diameter. What we need is a Z position
# depending on the XY positions. The tube is horizontally put with its
# axis along X.
sub Z_kantelknik_tube {
    my ($x, $y) = @_;

    # radius:
    my $r = 45/2;
    # only the y value is interesting:
    # y2 + z2 = r2
    my $z2 = $r*$r - $y * $y;
    die "no solution for $y\n" if($z2 < 0);
    return $r - sqrt($z2) + 30;
}

# give the point (x,y) to calculate z for, given ax+by+cz = d
# so z = d - ax - by / c
sub Z_plane {
    my ($x, $y, $a, $b, $c, $d) = @_;

    if($c < 0.00001) {
        # vertical plane: always give same answer....
        return 0;
    } else {
        return ($d - $a * $x - $b * $y) / $c;
    }
}

# The cilinder is at the $angle to the xy-plane, axis of the cilinder
# through the zero and in the zy-plane
sub Z_cilinder {
    my ($x, $y, $radius, $angle) = @_;

    # we cut of for to large:
    if($x >= $radius) {
        return $radius;
    } else {
        # x decides what angle in the crosssection we are going to hit:
        my $t = cos($angle) * sqrt(sqr($radius) - sqr($x));
        return 30 + $radius - ($t + $y * tan($angle));
    }
}

sub sqr {
    return $_[0] * $_[0];
}

sub min {
    return $_[0] < $_[1]? $_[0] : $_[1];
}

# make mm into 1/72 inch (what postscripts uses)
sub corps {
    return $_[0] * 72 / 25.4;
    #return $_[0];
}

sub moveto {
    my ($x, $y) = @_;
    print corps($x), " ", corps($y), " moveto\n";
}

sub lineto {
    my ($x, $y) = @_;
    print corps($x), " ", corps($y), " lineto\n";
}

sub line {
    my ($x1, $y1, $x2, $y2) = @_;
    moveto($x1, $y1);
    lineto($x2, $y2);
    print "stroke\n";
}

# calculates the trail around a square tube. parameter is the angle.
#sub calc_xyztrail {
#    my ($angle) = @_;
#    my ($a, $kxx, $kxy, $kyx, $kyy, $trail_parts, $ktrail);
#    if(0 <= $angle && $angle < $a45) {
#        ($a, $kxx, $kxy, $kyx, $kyy, $trail_parts, $ktrail) = ($angle, 1, 0, 0, 1, 0, 1);
#    } elsif($angle < $a90) {
#        ($a, $kxx, $kxy, $kyx, $kyy, $trail_parts, $ktrail) = ($a90 - $angle, 0, 1, 1, 0, 2, -1);
#    } elsif($angle < $a135) {
#        ($a, $kxx, $kxy, $kyx, $kyy, $trail_parts, $ktrail) = ($angle - $a90, 0, -1, 1, 0, 2, 1);
#    } elsif($angle < $a180) {
#        ($a, $kxx, $kxy, $kyx, $kyy, $trail_parts, $ktrail) = ($a180 - $angle, -1, 0, 0, 1, 4, -1);
#    } elsif($angle < $a225) {
#        ($a, $kxx, $kxy, $kyx, $kyy, $trail_parts, $ktrail) = ($angle - $a180, -1, 0, 0, -1, 4, 1);
#    } elsif($angle < $a270) {
#        ($a, $kxx, $kxy, $kyx, $kyy, $trail_parts, $ktrail) = ($a270 - $angle, 0, -1, -1, 0, 6, -1);
#    } elsif($angle < $a315) {
#        ($a, $kxx, $kxy, $kyx, $kyy, $trail_parts, $ktrail) = ($angle - $a270, 0, 1, -1, 0, 6, 1);
#    } else {
#        ($a, $kxx, $kxy, $kyx, $kyy, $trail_parts, $ktrail) = ($a360 - $angle, 1, 0, 0, -1, 8, -1);
#    }
#    # suppose first we are on the quadrant: 0 -> 45
#    my $x = $edge_35 / 2;
#    my $y = $edge_35 * tan($a) / 2;
#    # the trail is the number of complete quadrants
#    # plus a strait piece
#    my $trail = ($trail_parts * ($edge_35/2 - $edge_radius_35 + $edge_radius_35*$PI/4) +
#                 $ktrail * min($y, $edge_35/2 - $edge_radius_35));
#    if($y > ($edge_35/2 - $edge_radius_35)) {
#        # we are in the edge circle, correct x and y.
#        # its a bit tricky here:
#        my $R = sqrt(2) * ($edge_35/2 - $edge_radius_35); # the inner circle
#        # rules of the sines
#        my $c = asin(($R * sin($a45 - $a)) / $edge_radius_35);
#        # the angle under which the small edge circle is:
#        my $b = $a - $c;
#        $x = $edge_35/2 - $edge_radius_35 + $edge_radius_35 * cos($b);
#        $y = $edge_35/2 - $edge_radius_35 + $edge_radius_35 * sin($b);
#        $trail += $ktrail * $edge_radius_35 * $b;
#    }
#
#    # correcto for all quadrants
#    my $rx = $kxx * $x + $kxy * $y;
#    my $ry = $kyx * $x + $kyy * $y;
#
#    return ($rx, $ry, $trail);
#}

# calculates the trail around a square tube.
# parameters are the angle plus tube parameters
# tube_parameters:
# - width (x direction)
# - height (y direction)
# - edge radius

# note the angle used is not really from the centerpoint. We part the
# tube in 8 quadrants. Each quadrant gets its own centerpoint. The
# centerpoint of a quadrant is the point of the square from out of the
# edge that lies in the quadrant.

sub calc_xyztrail_tube {
    my ($angle, $width, $height, $edge_radius) = @_;

    my ($a, $kxx, $kxy, $kyx, $kyy, $trail_v_parts, $trail_h_parts, $ktrail);
    if(0 <= $angle && $angle < $a45) {
        ($a, $kxx, $kxy, $kyx, $kyy, $trail_v_parts, $trail_h_parts, $ktrail) = ($angle, 1, 0, 0, 1, 0, 0, 1);
    } elsif($angle < $a90) {
        ($a, $kxx, $kxy, $kyx, $kyy, $trail_v_parts, $trail_h_parts, $ktrail) = ($a90 - $angle, 0, 1, 1, 0, 1, 1, -1);
    } elsif($angle < $a135) {
        ($a, $kxx, $kxy, $kyx, $kyy, $trail_v_parts, $trail_h_parts, $ktrail) = ($angle - $a90, 0, -1, 1, 0, 1, 1, 1);
    } elsif($angle < $a180) {
        ($a, $kxx, $kxy, $kyx, $kyy, $trail_v_parts, $trail_h_parts, $ktrail) = ($a180 - $angle, -1, 0, 0, 1, 2, 2, -1);
    } elsif($angle < $a225) {
        ($a, $kxx, $kxy, $kyx, $kyy, $trail_v_parts, $trail_h_parts, $ktrail) = ($angle - $a180, -1, 0, 0, -1, 2, 2, 1);
    } elsif($angle < $a270) {
        ($a, $kxx, $kxy, $kyx, $kyy, $trail_v_parts, $trail_h_parts, $ktrail) = ($a270 - $angle, 0, -1, -1, 0, 3, 3, -1);
    } elsif($angle < $a315) {
        ($a, $kxx, $kxy, $kyx, $kyy, $trail_v_parts, $trail_h_parts, $ktrail) = ($angle - $a270, 0, 1, -1, 0, 3, 3, 1);
    } else {
        ($a, $kxx, $kxy, $kyx, $kyy, $trail_v_parts, $trail_h_parts, $ktrail) = ($a360 - $angle, 1, 0, 0, -1, 4, 4, -1);
    }
    # lenght of the strait part
    my $transformed_width = sqr($kxx)*$width + sqr($kxy)*$height;
    my $transformed_height = sqr($kxx)*$height + sqr($kxy)*$width;
    my $strait_part = $transformed_height/2 - $edge_radius;
    # suppose first we are on the quadrant: 0 -> 45
    my $x = $transformed_width / 2;
    my $y = $transformed_height * tan($a) / 2;
    # the trail is the number of complete quadrants
    # plus a strait piece
    my $trail = ($trail_v_parts * ($height/2 - $edge_radius + $edge_radius*$PI/4) +
                 $trail_h_parts * ($width/2 - $edge_radius + $edge_radius*$PI/4) +
                 $ktrail * min($y, $strait_part));
    if($y > $strait_part) {
        # we are in the edge circle, correct x and y.
        # its a bit tricky here:
        my $R = sqrt(2) * $strait_part; # the inner circle
        # rules of the sines
        my $c = asin(($R * sin($a45 - $a)) / $edge_radius);
        # the angle under which the small edge circle is:
        my $b = $a - $c;
        $x = $transformed_width/2 - $edge_radius + $edge_radius * cos($b);
        $y = $strait_part + $edge_radius * sin($b);
        $trail += $ktrail * $edge_radius * $b;
    }

    # correcto for all quadrants
    my $rx = $kxx * $x + $kxy * $y;
    my $ry = $kyx * $x + $kyy * $y;

    return ($rx, $ry, $trail);
}

sub calc_xyztrail_cilinder {
    my ($angle, $radius) = @_;

    my $trail = $radius * $angle;
    my $rx = $radius * cos($angle);
    my $ry = $radius * sin($angle);
    return ($rx, $ry, $trail);
}

sub draw_section_kantelknik_voorbouw {
    # the angle rotates between 0 and 360 degrees
    my $angle;
    
    # move the pic in a visable area
    print corps(50), " ", corps(50), " translate\n";
    # print some vertical alignment lines:
    for($angle = 0; $angle < ($a360+0.1); $angle += $a45) {
        my ($rx, $ry, $trail) = calc_xyztrail_tube($angle, $edge_35, $edge_35, $edge_radius_35);
        my $rz = Z_kantelknik_tube($rx,$ry);
        line($trail, 0, $trail, $rz);
    }
    # now comes the real graph
    moveto(0,0);
    my ($rx, $ry, $rz, $trail);
    
    for($angle = 0; $angle < $a360; $angle += $stepsize) {
        ($rx, $ry, $trail) = calc_xyztrail_tube($angle, $edge_35, $edge_35, $edge_radius_35);
        my $rz = Z_kantelknik_tube($rx,$ry);
        # the uitgerolde plot
        lineto($trail, $rz);
    }
    # output the thing
    print "stroke\n";
    line(0, 0, $trail, 0);
}

sub draw_section_cilinder_cilinder {
    # the angle rotates between 0 and 360 degrees
    my ($tube_angle, $radius1, $radius2) = @_;
    
    my $angle;
    my ($rx, $ry, $rz, $trail);
    
    # move the pic in a visable area
    print corps(50), " ", corps(50), " translate\n";
    # print some vertical alignment lines:
    for($angle = 0; $angle < ($a360+0.1); $angle += $a45) {
        my ($rx, $ry, $trail) = calc_xyztrail_cilinder($angle, $radius1);
        $rz = Z_cilinder($rx,$ry, $radius2, $tube_angle);
#        print "coord: $rx, $ry, $rz, angle $angle\n";
        line($trail, 0, $trail, $rz);
    }
    
    moveto(0,0);
    for($angle = 0; $angle < $a360; $angle += $stepsize) {
        ($rx, $ry, $trail) = calc_xyztrail_cilinder($angle, $radius1);
        $rz = Z_cilinder($rx,$ry, $radius2, $tube_angle);
        # the uitgerolde plot
        lineto($trail, $rz);
    }
    # output the thing
    print "stroke\n";
    line(0, 0, $trail, 0);
}

# parameters:
# - the angel between the tubes. ( 0 is strait through, a90 is loodrecht)
# - the thickness of the tube
# plus the tube parameters:
# - width,
# - height,
# - edge radius

# I suppose at the moment that you do not cut through the tube, but
# rather bend on face of the tube, I suppose that the inside of the wall
# will stay constant in lenght.

sub draw_section_tube_tube {
    my ($tube_angle, $thickness, $width, $height, $edge_radius) = @_;

    # move the pic in a visable area
#    
    # calculate a,b,c for ax+by+cz=1 (we rotate around Y axis)
    my @plane = (tan($tube_angle/2), 0, 1, 50);
    print corps(50), " ", corps(50), " translate\n";
    draw_section_tube_tube_half($width, $height, $edge_radius, @plane);
    my ($rx, $ry, $trail) = calc_xyztrail_tube($a180, $width, $height, $edge_radius);
    my $rz = Z_plane($rx,$ry, @plane);
    print corps(0), " ", corps(2*$rz), " translate\n";
    print "180 rotate\n";
    ($rx, $ry, $trail) = calc_xyztrail_tube($a360, $width, $height, $edge_radius);
    print -corps($trail), " 0 translate\n";
    # translate to effect the thickness of the wall
    print "0 ", corps(2*tan($tube_angle/2)*$thickness), " translate\n";
    draw_section_tube_tube_half($width, $height, $edge_radius, @plane);
}

sub draw_section_tube_tube_half {
    my ($width, $height, $edge_radius, @plane) = @_;
    my ($rx, $ry, $rz, $trail);
    my $angle;

    for($angle = 0; $angle < ($a360+0.1); $angle += $a45) {
        my ($rx, $ry, $trail) = calc_xyztrail_tube($angle, $width, $height, $edge_radius);
        my $rz = Z_plane($rx,$ry, @plane);
#        print "coord: $rx, $ry, $rz\n";
        line($trail, 0, $trail, $rz);
    }
    
    moveto(0,0);
    for($angle = 0; $angle < $a360; $angle += $stepsize) {
        ($rx, $ry, $trail) = calc_xyztrail_tube($angle, $width, $height, $edge_radius);
        $rz = Z_plane($rx, $ry, @plane);
        # the uitgerolde plot
        lineto($trail, $rz);
    }
    # output the thing
    print "stroke\n";
    line(0, 0, $trail, 0);
}