# -*- python -*-
# The first line tells emacs to use python mode for editing


#
# mold.cad
#
# molding and casting example
#
# Neil Gershenfeld
# CBA MIT 7/28/07
#
# (c) Massachusetts Institute of Technology 2007
# Permission granted for experimental and personal use;
# license for commercial sale available from MIT.
#

#
# define shapes and transformation
#
# circle(x0, y0, r)
# cylinder(x0, y0, z0, z1, r)
# sphere(x0, y0, z0, r)
# torus(x0, y0, z0, r0, r1)
# rectangle(x0, x1, y0, y1)
# cube(x0, x1, y0, y1, z0, z1)
# function(Z_of_XY)
# functions(upper_Z_of_XY,lower_Z_of_XY)
# add(part1, part2)
# subtract(part1, part2)
# intersect(part1, part2)
# move(part,dx,dy)
# translate(part,dx,dy,dz)
# rotate(part, angle)
# rotate_x(part, angle)
# rotate_y(part, angle)
# rotate_z(part, angle)
# rotate_z_90(part)
# rotate_z_180(part)
# rotate_z_270(part)
# reflect_x(part)
# reflect_y(part)
# reflect_z(part)
# reflect_xy(part)
# reflect_xz(part)
# reflect_yz(part)
# scale_x(part, x0, sx)
# scale_y(part, y0, sy)
# scale_z(part, z0, sz)
# scale_xy(part, x0, y0, sxy)
# scale_xyz(part, x0, y0, z0, sxyz)
# coscale_x_y(part, x0, y0, y1, angle0, angle1, amplitude, offset)
# coscale_x_z(part, x0, z0, z1, angle0, angle1, amplitude, offset)
# coscale_xy_z(part, x0, y0, z0, z1, angle0, angle1, amplitude, offset)
# taper_x_y(part, x0, y0, y1, s0, s1)
# taper_x_z(part, x0, z0, z1, s0, s1)
# taper_xy_z(part, x0, y0, z0, z1, s0, s1)
# shear_x_y(part, y0, y1, dx0, dx1)
# shear_x_z(part, z0, z1, dx0, dx1)
# (more to come)

# coshear

def circle(x0, y0, r):
   part = "(((X-x0)**2 + (Y-y0)**2) <= r**2)"
   part = replace(part,'x0',str(x0))
   part = replace(part,'y0',str(y0))
   part = replace(part,'r',str(r))
   return part

def cylinder(x0, y0, z0, z1, r):
   part = "(((X-x0)**2 + (Y-y0)**2 <= r**2) & (Z >= z0) & (Z <= z1))"
   part = replace(part,'x0',str(x0))
   part = replace(part,'y0',str(y0))
   part = replace(part,'z0',str(z0))
   part = replace(part,'z1',str(z1))
   part = replace(part,'r',str(r))
   return part

def sphere(x0, y0, z0, r):
   part = "(((X-x0)**2 + (Y-y0)**2 + (Z-z0)**2) <= r**2)"
   part = replace(part,'x0',str(x0))
   part = replace(part,'y0',str(y0))
   part = replace(part,'z0',str(z0))
   part = replace(part,'r',str(r))
   return part

def torus(x0, y0, z0, r0, r1):
   part = "(((r0 - sqrt((X-x0)**2 + (Y-y0)**2))**2 + (Z-z0)**2) <= r1**2)"
   part = replace(part,'x0',str(x0))
   part = replace(part,'y0',str(y0))
   part = replace(part,'z0',str(z0))
   part = replace(part,'r0',str(r0))
   part = replace(part,'r1',str(r1))
   return part

def rectangle(x0, x1, y0, y1):
   part = "((X >= x0) & (X <= x1) & (Y >= y0) & (Y <= y1))"
   part = replace(part,'x0',str(x0))
   part = replace(part,'x1',str(x1))
   part = replace(part,'y0',str(y0))
   part = replace(part,'y1',str(y1))
   return part

def cube(x0, x1, y0, y1, z0, z1):
   part = "((X >= x0) & (X <= x1) & (Y >= y0) & (Y <= y1) & (Z >= z0) & (Z <= z1))"
   part = replace(part,'x0',str(x0))
   part = replace(part,'x1',str(x1))
   part = replace(part,'y0',str(y0))
   part = replace(part,'y1',str(y1))
   part = replace(part,'z0',str(z0))
   part = replace(part,'z1',str(z1))
   return part

def function(Z_of_XY):
   part = '(Z <= '+Z_of_XY+')'
   return part

def functions(upper_Z_of_XY,lower_Z_of_XY):
   part = '(Z <= '+upper_Z_of_XY+') & (Z >= '+lower_Z_of_XY+')'
   return part

def add(part1, part2):
   part = "(part1) | (part2)"
   part = replace(part,'part1',part1)
   part = replace(part,'part2',part2)
   return part

def subtract(part1, part2):
   part = "(part1) & ~(part2)"
   part = replace(part,'part1',part1)
   part = replace(part,'part2',part2)
   return part

def intersect(part1, part2):
   part = "(part1) & (part2)"
   part = replace(part,'part1',part1)
   part = replace(part,'part2',part2)
   return part

def move(part,dx,dy):
   part = replace(part,'X','(X-'+str(dx)+')')
   part = replace(part,'Y','(Y-'+str(dy)+')')
   return part

def translate(part,dx,dy,dz):
   part = replace(part,'X','(X-'+str(dx)+')')
   part = replace(part,'Y','(Y-'+str(dy)+')')
   part = replace(part,'Z','(Z-'+str(dz)+')')
   return part

def rotate(part, angle):
   angle = angle*pi/180
   part = replace(part,'X','(cos(angle)*X+sin(angle)*y)')
   part = replace(part,'Y','(-sin(angle)*X+cos(angle)*y)')
   part = replace(part,'y','Y')
   part = replace(part,'angle',str(angle))
   return part

def rotate_x(part, angle):
   angle = angle*pi/180
   part = replace(part,'Y','(cos(angle)*Y+sin(angle)*z)')
   part = replace(part,'Z','(-sin(angle)*Y+cos(angle)*z)')
   part = replace(part,'z','Z')
   part = replace(part,'angle',str(angle))
   return part

def rotate_y(part, angle):
   angle = angle*pi/180
   part = replace(part,'X','(cos(angle)*X+sin(angle)*z)')
   part = replace(part,'Z','(-sin(angle)*X+cos(angle)*z)')
   part = replace(part,'z','Z')
   part = replace(part,'angle',str(angle))
   return part

def rotate_z(part, angle):
   angle = angle*pi/180
   part = replace(part,'X','(cos(angle)*X+sin(angle)*y)')
   part = replace(part,'Y','(-sin(angle)*X+cos(angle)*y)')
   part = replace(part,'y','Y')
   part = replace(part,'angle',str(angle))
   return part

def rotate_z_90(part):
   part = reflect_xy(part)
   part = reflect_y(part)
   return part

def rotate_z_180(part):
   part = reflect_xy(part)
   part = reflect_y(part)
   part = reflect_xy(part)
   part = reflect_y(part)
   return part

def rotate_z_270(part):
   part = reflect_xy(part)
   part = reflect_y(part)
   part = reflect_xy(part)
   part = reflect_y(part)
   part = reflect_xy(part)
   part = reflect_y(part)
   return part

def reflect_x(part):
   part = replace(part,'X','-X')
   return part

def reflect_y(part):
   part = replace(part,'Y','-Y')
   return part

def reflect_z(part):
   part = replace(part,'Z','-Z')
   return part

def reflect_xy(part):
   part = replace(part,'X','temp')
   part = replace(part,'Y','X')
   part = replace(part,'temp','Y')
   return part

def reflect_xz(part):
   part = replace(part,'X','temp')
   part = replace(part,'Z','X')
   part = replace(part,'temp','Z')
   return part

def reflect_yz(part):
   part = replace(part,'Y','temp')
   part = replace(part,'Z','Y')
   part = replace(part,'temp','Z')
   return part

def scale_x(part, x0, sx):
   part = replace(part,'X','(x0 + (X-x0)/sx)')
   part = replace(part,'x0',str(x0))
   part = replace(part,'sx',str(sx))
   return part

def scale_y(part, y0, sy):
   part = replace(part,'Y','(y0 + (Y-y0)/sy)')
   part = replace(part,'y0',str(y0))
   part = replace(part,'sy',str(sy))
   return part

def scale_z(part, z0, sz):
   part = replace(part,'Z','(z0 + (Z-z0)/sz)')
   part = replace(part,'z0',str(z0))
   part = replace(part,'sz',str(sz))
   return part

def scale_xy(part, x0, y0, sxy):
   part = replace(part,'X','(x0 + (X-x0)/sx)')
   part = replace(part,'Y','(y0 + (Y-y0)/sy)')
   part = replace(part,'x0',str(x0))
   part = replace(part,'y0',str(y0))
   part = replace(part,'sxy',str(sxy))
   return part

def scale_xyz(part, x0, y0, z0, sxy):
   part = replace(part,'X','(x0 + (X-x0)/sx)')
   part = replace(part,'Y','(y0 + (Y-y0)/sy)')
   part = replace(part,'Z','(z0 + (Z-z0)/sz)')
   part = replace(part,'x0',str(x0))
   part = replace(part,'y0',str(y0))
   part = replace(part,'z0',str(z0))
   part = replace(part,'sxyz',str(sxyz))
   return part

def coscale_x_y(part, x0, y0, y1, angle0, angle1, amplitude, offset):
   phase0 = pi*angle0/180.
   phase1 = pi*angle1/180.
   part = replace(part,'X','(x0 + (X-x0)/(offset + amplitude*cos(phase0 + (phase1-phase0)*(Y-y0)/(y1-y0))))')
   part = replace(part,'x0',str(x0))
   part = replace(part,'y0',str(y0))
   part = replace(part,'y1',str(y1))
   part = replace(part,'phase0',str(phase0))
   part = replace(part,'phase1',str(phase1))
   part = replace(part,'amplitude',str(amplitude))
   part = replace(part,'offset',str(offset))
   return part

def coscale_x_z(part, x0, z0, z1, angle0, angle1, amplitude, offset):
   phase0 = pi*angle0/180.
   phase1 = pi*angle1/180.
   part = replace(part,'X','(x0 + (X-x0)/(offset + amplitude*cos(phase0 + (phase1-phase0)*(Z-z0)/(z1-z0))))')
   part = replace(part,'x0',str(x0))
   part = replace(part,'z0',str(z0))
   part = replace(part,'z1',str(z1))
   part = replace(part,'phase0',str(phase0))
   part = replace(part,'phase1',str(phase1))
   part = replace(part,'amplitude',str(amplitude))
   part = replace(part,'offset',str(offset))
   return part

# starting with: 'X<20'
# end up with: 'X-A*(sin(Y*B/pi))<20';


# curvify_z_x modifies the Z component of a part to be moved by some
# amount relative to sin(x)
def curvify_z_x(part, coeff, phi, offset):
   part = replace(part,'Z','Z-coeff*(sin((X-offset)*phi/pi))')
   part = replace(part,'coeff',str(coeff))
   part = replace(part,'phi',str(phi))
   part = replace(part,'offset',str(offset))
   return part


# curvify_z_xy modifies the Z component of a part to be moved by some
# amount relative to sin(d), where d is the distance from the origin
# in the x-y plane. In other words, if you start with a flat surface
# (constant z value), this will distort it such that there are
# radiating waves forming concentric circles around the origin.
def curvify_z_xy(part, coeff, phi, offset):
   part = replace(part,'Z','Z-coeff*(sin((sqrt(X**2+Y**2)-offset)*phi/pi))')
   part = replace(part,'coeff',str(coeff))
   part = replace(part,'phi',str(phi))
   part = replace(part,'offset',str(offset))
   return part


def coscale_xy_z(part, x0, y0, z0, z1, angle0, angle1, amplitude, offset):
   phase0 = pi*angle0/180.
   phase1 = pi*angle1/180.
   part = replace(part,'X','(x0 + (X-x0)/(offset + amplitude*cos(phase0 + (phase1-phase0)*(Z-z0)/(z1-z0))))')
   part = replace(part,'Y','(y0 + (Y-y0)/(offset + amplitude*cos(phase0 + (phase1-phase0)*(Z-z0)/(z1-z0))))')
   part = replace(part,'x0',str(x0))
   part = replace(part,'y0',str(y0))
   part = replace(part,'z0',str(z0))
   part = replace(part,'z1',str(z1))
   part = replace(part,'phase0',str(phase0))
   part = replace(part,'phase1',str(phase1))
   part = replace(part,'amplitude',str(amplitude))
   part = replace(part,'offset',str(offset))
   return part

def taper_x_y(part, x0, y0, y1, s0, s1):
   part = replace(part,'X','(x0 + (X-x0)*(y1-y0)/(s1*(Y-y0) + s0*(y1-Y)))')
   part = replace(part,'x0',str(x0))
   part = replace(part,'y0',str(y0))
   part = replace(part,'y1',str(y1))
   part = replace(part,'s0',str(s0))
   part = replace(part,'s1',str(s1))
   return part

def taper_x_z(part, x0, z0, z1, s0, s1):
   part = replace(part,'X','(x0 + (X-x0)*(z1-z0)/(s1*(Z-z0) + s0*(z1-Z)))')
   part = replace(part,'x0',str(x0))
   part = replace(part,'z0',str(z0))
   part = replace(part,'z1',str(z1))
   part = replace(part,'s0',str(s0))
   part = replace(part,'s1',str(s1))
   return part

def taper_xy_z(part, x0, y0, z0, z1, s0, s1):
   part = replace(part,'X','(x0 + (X-x0)*(z1-z0)/(s1*(Z-z0) + s0*(z1-Z)))')
   part = replace(part,'Y','(y0 + (Y-y0)*(z1-z0)/(s1*(Z-z0) + s0*(z1-Z)))')
   part = replace(part,'x0',str(x0))
   part = replace(part,'y0',str(y0))
   part = replace(part,'z0',str(z0))
   part = replace(part,'z1',str(z1))
   part = replace(part,'s0',str(s0))
   part = replace(part,'s1',str(s1))
   return part

def shear_x_y(part, y0, y1, dx0, dx1):
   part = replace(part,'X','(X - dx0 - (dx1-dx0)*(Y-y0)/(y1-y0))')
   part = replace(part,'y0',str(y0))
   part = replace(part,'y1',str(y1))
   part = replace(part,'dx0',str(dx0))
   part = replace(part,'dx1',str(dx1))
   return part

def shear_x_z(part, z0, z1, dx0, dx1):
   part = replace(part,'X','(X - dx0 - (dx1-dx0)*(Z-z0)/(z1-z0))')
   part = replace(part,'z0',str(z0))
   part = replace(part,'z1',str(z1))
   part = replace(part,'dx0',str(dx0))
   part = replace(part,'dx1',str(dx1))
   return part

def coshear_x_z(part, z0, z1, angle0, angle1, amplitude, offset):
   phase0 = pi*angle0/180.
   phase1 = pi*angle1/180.
   part = replace(part,'X','(X - offset - amplitude*cos(phase0 + (phase1-phase0)*(Z-z0)/(z1-z0)))')
   part = replace(part,'z0',str(z0))
   part = replace(part,'z1',str(z1))
   part = replace(part,'phase0',str(phase0))
   part = replace(part,'phase1',str(phase1))
   part = replace(part,'amplitude',str(amplitude))
   part = replace(part,'offset',str(offset))
   return part

#
# define part
#

bodyRadius = 0.4
tubeRadius=0.125
depth = 1.5
toolD=0.125
sweepRadius = bodyRadius*2+toolD
ZOffset = -.21
#borderSlop = tubeRadius*2
borderSlop = toolD*2.5
bigNum = 10


########
# Define body
########
tophalf = sphere(0,0,0,bodyRadius)
tophalf = scale_z(tophalf, 0, 1.5)
tophalf = subtract(tophalf,cube(-bigNum,bigNum,-bigNum,bigNum,-bigNum,0))
bottomhalf = sphere(0,0,0,bodyRadius)
bottomhalf = scale_z(bottomhalf, 0, 0.5)
bottomhalf = subtract(bottomhalf,cube(-bigNum,bigNum,-bigNum,bigNum,0,bigNum))
body = add(tophalf, bottomhalf)



upper_mold = 'True'

eyeVent = cylinder(0,0,0,1.0+ZOffset,tubeRadius)
eyeVentR = translate(eyeVent,tubeRadius*2,0,0)
eyeVentL = rotate_z(eyeVentR,70)

leg1 = cylinder(0,0,0,10,tubeRadius)
leg1 = rotate_y(leg1,90)
#leg1 = curvify_z_x(leg1,0.2,30,0)
# Don't rotate the leg downward for now... this introduces overhang, and would require a 3-part mold
#leg1 = rotate_y(leg1,20)


topmask = cube(-bigNum, bigNum, -bigNum, bigNum, 0, bigNum)

legs = add(leg1, rotate_z(leg1, 1*360/8))
legs = add(legs, rotate_z(leg1, 2*360/8))
legs = add(legs, rotate_z(leg1, 3*360/8))
legs = add(legs, rotate_z(leg1, 4*360/8))
legs = add(legs, rotate_z(leg1, 5*360/8))
legs = add(legs, rotate_z(leg1, 6*360/8))
legs = add(legs, rotate_z(leg1, 7*360/8))

topmask = curvify_z_xy(topmask,0.2,30,bodyRadius)
legs = curvify_z_xy(legs,0.2,30,bodyRadius)


lower_mold = 'True'


octopus = add(body,eyeVentR)
octopus = add(octopus, eyeVentL)

octopus = add(octopus, legs)

wrapper = 'True'
#xx = sweepRadius+2*toolD
xx = sweepRadius
wrapper = subtract(wrapper,cube(-xx,xx,-xx,xx,-xx,xx))
topwrapper = subtract(wrapper,cube(-bigNum,bigNum,-bigNum,bigNum,-bigNum,-tubeRadius/2))
bottomwrapper = subtract(wrapper,cube(-bigNum,bigNum,-bigNum,bigNum,-tubeRadius/2,bigNum))
ditch = cube(-(xx+toolD),(xx+toolD),-(xx+toolD),(xx+toolD),-tubeRadius/2-toolD/2,0)
ditch = subtract(ditch,cube(-xx,xx,-xx,xx,-tubeRadius/2-toolD/2,0))

upper_mold = subtract(topmask,octopus)
upper_mold = add(upper_mold,topwrapper)
upper_mold = subtract(upper_mold,bottomwrapper)

lower_mold = 'True'
lower_mold = subtract(lower_mold, topmask)
lower_mold = subtract(lower_mold, octopus)
lower_mold = add(lower_mold, bottomwrapper)
lower_mold = subtract(lower_mold, topwrapper)
lower_mold = subtract(lower_mold,ditch)

universe=cube(-bigNum,bigNum,-bigNum,bigNum,-bodyRadius,bodyRadius)
universe = subtract(universe,cube(-(xx+1.5*toolD),(xx+1.5*toolD),-(xx+1.5*toolD),(xx+1.5*toolD),-bodyRadius*1.5,bodyRadius*1.5))

upper_mold = subtract(upper_mold,universe)
lower_mold = subtract(lower_mold,universe)


upper_mold = reflect_z(upper_mold)
topmask = reflect_z(topmask)


#part = brick
part = upper_mold
#part = lower_mold

# The following line is just a cut-out to help visualization.
# Comment it out when actually cutting.
#part = subtract(part,cube(0,10,-10,0,-10,10))

#part = translate(part,1+sweepRadius,1+sweepRadius,-depth/2)
part = translate(part,1+sweepRadius,1+sweepRadius,ZOffset)

# torus experiment:
#part = torus(0, 0, 0, 0.4, 0.1)
#part = translate(part,2,2,-1)

#
# define limits and parameters
#

cad.xmin = 1-borderSlop
cad.xmax = 1+2*sweepRadius+borderSlop
cad.ymin = 1-borderSlop
cad.ymax = 1+2*sweepRadius+borderSlop
cad.zmin = -depth
cad.zmax = 0
cad.rx = 20 # x view rotation (degrees)
cad.rz = 20 # z view rotation (degrees)
#dpi = 20 # low resolution for previewing
dpi = 70 # high resolution for machining
nxy = int(dpi*(cad.xmax-cad.xmin))
cad.nx = nxy
cad.ny = nxy
dz = 0.02
cad.nz = int((cad.zmax-cad.zmin)/dz)
cad.inches_per_unit = 1.0 # use inch units

#
# assign part to cad.function
#

cad.function = part