#!/usr/bin/python

import sys, math


group = None
g = 9.81	# gravitational acceleration, m/s2

# density, g/cm3

density = 11.34	# pure lead (Pb)
density =  9.31	# Pb50Sn50
density = 10.00	# Pb67Sn33
#density =  7.28	# pure tin (Sn)

#
# The z coordinate of the plane limiting the top of the counterweight. This is
# the altitude of the board surface minus the board clearance.
#
z_ceiling = 5.0	# mm

#
# The y coordinate of the axis around which our system rotates, i.e., the
# position of the center of the rear feet
#
y_axis = 16.0

off_x = -15+18
off_y = -46+5

channel_radius = 1	# mm
inlet_radius = 7
shaft = 3
ymax_piece = 50

total_mass = 0
total_torque = 0


#
# solve a quadratic equation of the form a*x^2+b*x+c = 0
#

def qeq(a, b, c):
    d = math.sqrt(b*b-4*a*c)
    return ((-b-d)/2/a, (-b+d)/2/a)


#
# find the x-coordinate of the center of mass of a trapezoid/trapezium with the
# four corners (0, 0), (x, 0), (0, y0), and (x, y0)
# we assume the mass distribution to be uniform
#

def cm_trap_a(x, y0, y1):
    if y0 == y1:
	return x/2.0
    f = float(y1-y0)/x/2
    return qeq(2*f, 2.0*y0, -x*(y0+y1)/2.0)[1]


#
# calculate a rectangle's contribution to mass and torque
#

def rect_calc(x0, y0, z0, x1, y1, z1):
    global total_mass, total_torque

    # mass, in g
    m = (x1-x0)*(y1-y0)*(z_ceiling-(z0+z1)/2.0)*density/1e3;

    # center of mass on y axis, in y coordinates (mm)
    y_center = y0+cm_trap_a(y1-y0, z_ceiling-z0, z_ceiling-z1)

    # weight, in N
    w = m*g/1000.0

    # torque, in Nm
    t = w*(y_center-y_axis)/1000.0

    total_mass += m
    total_torque += t

#
# gnuplot a rectangle
#

def rect_gnuplot(x0, y0, z0, x1, y1, z1):
    print x0, y0, z0
    print x1, y0, z0
    print x1, y1, z1
    print x0, y1, z1
    print x0, y0, z0
    print
    print


#
# add inlets and air escapes for gravitation casting
#

def channel(sk, x, y, r0, r1):
    if r0 == r1:
	cad.cylinder(x, y, 0, r0, ymax_piece-y)
        obj = cad.getlastobj()
    else:
	cad.cylinder(x, y, 0, r0, shaft)
	cyl = cad.getlastobj()
	cad.cone(x, y, shaft, r0, r1, ymax_piece-y-shaft)
	cone = cad.getlastobj()
	cad.fuse(cyl, cone)
	obj = cad.getlastobj()
    cad.rotate(obj, x, y, 0, 1, 0, 0, -math.pi/2)
    cad.cut(sk, obj)
    return cad.getlastobj()


def inlet(sk, x, y):
    return channel(sk, x+off_x, y+off_y, channel_radius, inlet_radius)
    pass


def escape(sk, x, y):
    return channel(sk, x+off_x, y+off_y, channel_radius, channel_radius)


#
# add a rectangle to the CAD model
#

def do_rect_cad(x0, y0, z0, x1, y1, z1):
    cad.sketch()
    sk = cad.getlastobj()

    cad.line3d(x0, y0, z0, x1, y0, z0)
    line = cad.getlastobj()
    cad.add(sk, line)

    cad.line3d(x1, y0, z0, x1, y1, z1)
    line = cad.getlastobj()
    cad.add(sk, line)

    cad.line3d(x1, y1, z1, x0, y1, z1)
    line = cad.getlastobj()
    cad.add(sk, line)

    cad.line3d(x0, y1, z1, x0, y0, z0)
    line = cad.getlastobj()
    cad.add(sk, line)

    cad.reorder(sk)

    return sk


def rect_cad(x0, y0, z0, x1, y1, z1):
    global group

    sk = do_rect_cad(x0, y0, z0, x1, y1, z1)

    cad.extrude(sk, 3)

    if group is None:
	group = cad.getlastobj()
    else:
	cad.fuse(group, cad.getlastobj())
	group = cad.getlastobj()


#
# add a rectangle with the following corners:
# (x0, y0, z0)
# (x1, y0, z0)
# (x0, y1, z1)
# (x1, y1, z1)
#

def rect(x0, y0, z0, x1, y1, z1):
    rect_calc(x0, y0, z0, x1, y1, z1)
    do(x0, y0, z0, x1, y1, z1)


#
# make the base
#

def make_base():
    rect(16,   46.0, 2.6,  22,   55,   2.6)	# left lateral, bottom
    rect(22,   46.0, 3.7,  24,   55,   3.7)	# on pedestal
    rect(89.5, 46.0, 2.6,  99.5, 55,   2.6)	# right lateral, bottom
    rect(88,   46.0, 3.7,  89.5, 55,   3.7)	# on pedestal

    rect(29.5, 50,   3.7,  36,   55,   3.7)	# left podium, to beam
    rect(38.5, 50,   3.7,  41,   55,   3.7)	# after beam
    rect(71,   50,   3.7,  82,   55,   3.7)	# right podium

    rect(16,   55,   2.6,  36,   60,   3.9)	# middle bar, to 1st beam
    rect(38.5, 55,   2.6,  46,   60,   3.9)	# between beams
    rect(48,   55,   2.6,  52.5, 60,   3.9)	# 2nd beam to battery lid
    rect(59.5, 55,   2.6,  99.5, 60,   3.9)	# right of battery lid "tongue"

    rect(36,   56,   3.9,  38.5, 60,   3.9)	# cover the beams
    rect(46,   56,   3.9,  48,   60,   3.9)

    rect(52.5, 56.5, 3.0,  59.5, 60,   3.9)	# cover the lid "tongue"

    rect(15,   60,   3.9,  60.5, 64,   3.9)	# end bar, left of beam
    rect(62.5, 60,   3.9, 100, 64,   3.9)	# right of beam

    rect(15,   64,   3.9,  34,   69.5, 3.9)	# left "ear"
    rect(89.5, 64,   3.9, 100,   69,   3.9)	# right "ear"

    rect(34,   64,   3.9,  60.5, 65,   3.9)	# extend inner area to sponge
    rect(62.5, 64,   3.9,  82.5, 65,   3.9)


if __name__ == "__main__":
    do = rect_gnuplot
else:
    import HeeksPython as cad
    do = rect_cad

make_base()

#
# for wax model
#
#if __name__ !=  "__main__":
#    sk = do_rect_cad(10, 40, z_ceiling, 110, 70, z_ceiling)
#    cad.extrude(sk, 3)
#    sk = cad.getlastobj()
#    cad.cut(group, sk)
#    group = cad.getlastobj()
#    cad.translate(group, -15, -69, -5)
#    cad.rotate(group, 0, 0, 0, 1, 0, 0, math.pi)

#
# add rectangular block for mold, then subtract the counterweight
#

if __name__ !=  "__main__":
    cad.translate(group, -15, -46, -5)
    cad.translate(group, 18, 5, 0)
    sk = do_rect_cad(0, 0, 0,  120, 45, 0)
    cad.extrude(sk, -10)
    sk = cad.getlastobj()
    cad.cut(sk, group)
    sk = cad.getlastobj()
    sk = escape(sk, 15+channel_radius, 69.5)
    sk = inlet(sk, 34-channel_radius, 69.5)
    sk = inlet(sk, 89.5+channel_radius, 69)
    sk = escape(sk, 100-channel_radius, 69)
    sk = escape(sk, 60.5-channel_radius, 65)
    sk = escape(sk, 62.5+channel_radius, 65)
    sk = inlet(sk, 50, 65)

print >>sys.stderr, "total mass =", total_mass, "g"
print >>sys.stderr, "total torque =", total_torque*1000.0, "mNm"