from math import sin, cos, acos, tan, atan, sqrt
from math import degrees, radians, pi
import Materials;

materialManager = Materials.MaterialManager()
wood = materialManager.getMaterial('ba2474ee-f62c-45f5-b388-823ea105847f')
wood.Name
wood.Description

ToolDiameter = FreeCAD.Units.Quantity('3 mm')
ToolNumberOfFlutes = 2
ToolRakeAngle = FreeCAD.Units.Quantity('30°')
ToolHelixAngle = FreeCAD.Units.Quantity('15°')

ToolMaxChipLoad = FreeCAD.Units.Quantity('0.030 mm') # not a tool setting; differs per material! (ToolMaxTorque would be nice but no vendor specifies this. And for soft materials large chips jam the bit before max torque is reached)

ae = FreeCAD.Units.Quantity('3 mm') # width of cut (radial)
ap = FreeCAD.Units.Quantity('5 mm') # depth of cut (axial)

kc11 = FreeCAD.Units.Quantity(wood.PhysicalProperties['UnitCuttingForce'])
h0 = FreeCAD.Units.Quantity('1 mm') # unit chip thickness, per definition 1mm for k_c1.1
mc = float(wood.PhysicalProperties['ChipThicknessExponent'])

# project angle: https://math.stackexchange.com/questions/2207665/projecting-an-angle-from-one-plane-to-another-plane
# not really worth taking the helix into account here; below 40° the effect is neglectable
gamma_eff = degrees(atan(tan(ToolRakeAngle.getValueAs("rad")/cos(ToolHelixAngle.getValueAs("rad")))))
gamma = ToolRakeAngle.getValueAs("deg")
Kg = 1 - 0.01 * gamma_eff # correction factor for rake angle

kapr = radians(90) # straight milling cutter, i.e. chamfer=90° aka no chamfer

D = ToolDiameter

phie = acos(1 - (2 * ae / D)) # engangement angle

# honor chip-thinning: calculate fz from h_max when phi_e < 90°
hmax = ToolMaxChipLoad # maximum undeformed chip thickness
if ae >  D:
    raise ValueError("Width of cut must not exceed the tool diameter")
elif ae >= D/2:
    fz = hmax # feed per tooth
else:
    # due to the curvature of the cutter's cross section (i.e. a circle with diameter D)
    # we must increase the feed per tooth to still reach h_max for shallow width of cut (radial)
    # reference: https://www.machiningdoctor.com/calculators/chip-thinning-calculator/
    fz = hmax  / sqrt(1 - (1 - 2*ae/D)**2)

Sb = D * pi * (phie / (2*pi)) # chip arc length
hm = fz * (ae/Sb) * sin(kapr) # mean undeformed chip thickness using Cavalieri's principle

Kw = 1.2 # correction factor for tool wear: 1 for new sharp tools, 1.2 for used tools, 1.5 for dull tools that need to be replaced

kc = kc11 * (hm/h0)**-mc * Kg * Kw # specific cutting force

Fcz = ap * hm * kc # cutting force per flute

z = ToolNumberOfFlutes

ze = phie * z / (2*pi) # engaged flutes

Fc = Fcz * ze # cutting force

vc_set = FreeCAD.Units.Quantity(wood.PhysicalProperties['SurfaceSpeedCarbide'])
vc_set.getValueAs("m/min")

n_set = vc_set / (pi * D) # spindle speed
n_set.getValueAs("1/min")

n_max = FreeCAD.Units.Quantity("30000/min")
n = min(n_set, n_max)
n.getValueAs("1/min")

vc = n * (pi * D)
vc.getValueAs("m/min")

Pc = Fc * vc # mechanical cutting power

eff = 0.85 # machine efficiency: 
P = Pc / eff # electrical spindle power
P.getValueAs("kW")

Mc = Fc * D / 2 # cutting torque (maybe better base this on h_max instead of h_mean?)
Mc.getValueAs("Nm")

vf = n * z * fz # feed rate
vf.getValueAs("mm/min")