1. Laser Cut at the 'Critical-Section'

(********************************************************)
(* Clear variables so that they no longer have old, junk values *)
(********************************************************)
Quiet[
Clear[
x,y,z,θ,α,
ρ,(Subscript[ρ, n]),ψ,
(Subscript[d, o]),
(OverBar[t]),(Subscript[(OverBar[d]), h]),(Subscript[(OverBar[d]), c]),(OverBar[k]),
(Subscript[r, o]),(Subscript[r, i]),(Subscript[r, h]),(Subscript[r, c]),t,k,(Subscript[k, s]),ϕ,(Subscript[ϕ, deg])
]
];
(********************************************************)


(*******************************)
(* Assign values to the parameters *)
(*******************************)

(OverBar[k])  = 0.0029972;
(Subscript[ϕ, deg]) = 0;
(Subscript[d, o])= 5;
(OverBar[t])  = 0.1;
(Subscript[(OverBar[d]), c])= 1.0;
(Subscript[(OverBar[d]), h])= 0.5;

(Subscript[r, o])=((Subscript[d, o])/2);
(*******************************)


(**********************************************************************)
(* Make rules for converting from normalized parameters to absolute quantities *)
(**********************************************************************)
substitutionRules = {

t      ->    (OverBar[t]) (Subscript[d, o]),

(Subscript[r, i])    ->    (Subscript[r, o])-(OverBar[t]) (Subscript[d, o]),

(Subscript[r, h])    ->    (1/2) (Subscript[(OverBar[d]), h]) (Subscript[d, o]),

(Subscript[r, c])    ->    (1/2) (Subscript[(OverBar[d]), c]) (Subscript[d, o]),

k      ->    (OverBar[k]) (Subscript[d, o]),

(Subscript[k, s])    ->    (((OverBar[k]) (Subscript[d, o]))/2),

ϕ      ->    (((Subscript[ϕ, deg]) π)/180.0)

};
 (**********************************************************************)

(*********************************)
(* Define equations and inequalities *)
(*********************************)

(Subscript[pathOfCutAtCriticalSection, implicit])=
(x^2)+(z^2)==(Subscript[r, h]^2)   && (x^2)+(y^2)==(Subscript[r, c]^2);

(Subscript[pathOfCutAtCriticalSection, parametric]) = (({{(Subscript[r, h]) Cos[ψ]
}, {(Sqrt[(Subscript[r, c]^2)-((((Subscript[r, h]) Cos[ψ]))^2)])
}, {(Subscript[r, h]) Sin[ψ]}}));

(*********************************)

Print["(Subscript[pathOfCutAtCriticalSection, implicit]) : ", (Subscript[pathOfCutAtCriticalSection, implicit])];
Print["(Subscript[pathOfCutAtCriticalSection, parametric]) : ", MatrixForm[(Subscript[pathOfCutAtCriticalSection, parametric])]];

(*******)
(* Plot *)
(*******)
AbsoluteTiming[
(Subscript[pathOfCutAtCriticalSection, parametric-plot]) = 
ParametricPlot3D[
ReplaceAll[
(Subscript[pathOfCutAtCriticalSection, parametric]),
substitutionRules
],
{ψ, 0, 2*Pi},
Axes->Automatic,
AxesLabel->{"X", "Y", "Z"}
]
]

2. Laser Cut at a distance 'r' from the tube axis

(********************************************************)
(* Clear variables so that they no longer have old, junk values *)
(********************************************************)
Quiet[
Clear[
x,y,z,θ,α,
ρ,(Subscript[ρ, n]),ψ,
(Subscript[d, o]),
(OverBar[t]),(Subscript[(OverBar[d]), h]),(Subscript[(OverBar[d]), c]),(OverBar[k]),
(Subscript[r, o]),(Subscript[r, i]),(Subscript[r, h]),(Subscript[r, c]),t,k,(Subscript[k, s]),ϕ,(Subscript[ϕ, deg])
]
];
(********************************************************)

(*******************************)
(* Assign values to the parameters *)
(*******************************)

(OverBar[k])  = 0.0029972;
(Subscript[ϕ, deg]) = 0;
(Subscript[d, o])= 5;
(OverBar[t])  = 0.1;
(Subscript[(OverBar[d]), c])= 1.0;
(Subscript[(OverBar[d]), h])= 0.5;

(Subscript[r, o])=((Subscript[d, o])/2);
(*******************************)

(**********************************************************************)
(* Make rules for converting from normalized parameters to absolute quantities *)
(**********************************************************************)
substitutionRules = {

t      ->    (OverBar[t]) (Subscript[d, o]),

(Subscript[r, i])    ->    (Subscript[r, o])-(OverBar[t]) (Subscript[d, o]),

(Subscript[r, h])    ->    (1/2) (Subscript[(OverBar[d]), h]) (Subscript[d, o]),

(Subscript[r, c])    ->    (1/2) (Subscript[(OverBar[d]), c]) (Subscript[d, o]),

k      ->    (OverBar[k]) (Subscript[d, o]),

(Subscript[k, s])    ->    (((OverBar[k]) (Subscript[d, o]))/2),

ϕ      ->    (((Subscript[ϕ, deg]) π)/180.0)

};
 (**********************************************************************)

(*********************************)
(* Define equations and inequalities *)
(*********************************)

(Subscript[pathOfCutAtAnySection, implicit])=
((x^2)/(((((Subscript[r, c])/ρ)))^2))+ ((z^2)/(1^2))==(Subscript[r, h]^2) && (x^2)+(y^2)==(ρ^2);

(Subscript[pathOfCutAtAnySection, parametric]) = (({{(((Subscript[r, h]) Cos[ψ])/((Subscript[r, c])/ρ))
}, {(Sqrt[(ρ^2)-((((((Subscript[r, h]) Cos[ψ])/((Subscript[r, c])/ρ))))^2)])
}, {(Subscript[r, h]) Sin[ψ]}}));

(*********************************)

Print["(Subscript[pathOfCutAtAnySection, implicit]) : ", (Subscript[pathOfCutAtAnySection, implicit])];
Print["(Subscript[pathOfCutAtAnySection, parametric]) : ", MatrixForm[(Subscript[pathOfCutAtAnySection, parametric])]];

(*******)
(* Plot *)
(*******)
AbsoluteTiming[
ParametricPlot3D[
ReplaceAll[
(Subscript[pathOfCutAtAnySection, parametric]),
substitutionRules
],
{ψ, 0, 2*Pi},
{ρ, ReplaceAll[(Subscript[r, i]),substitutionRules], (Subscript[r, o])},
Axes->Automatic,
AxesLabel->{"X", "Y", "Z"},
PlotPoints -> 50
]
]

3. Equation of the Laser-Beam Axis

(********************************************************)
(* Clear variables so that they no longer have old, junk values *)
(********************************************************)
Quiet[
Clear[
x,y,z,θ,α,
ρ,(Subscript[ρ, n]),ψ,
(Subscript[d, o]),
(OverBar[t]),(Subscript[(OverBar[d]), h]),(Subscript[(OverBar[d]), c]),(OverBar[k]),
(Subscript[r, o]),(Subscript[r, i]),(Subscript[r, h]),(Subscript[r, c]),t,k,(Subscript[k, s]),ϕ,(Subscript[ϕ, deg])
]
];
(********************************************************)

(*******************************)
(* Assign values to the parameters *)
(*******************************)

(OverBar[k])  = 0.0029972;
(Subscript[ϕ, deg]) = 0;
(Subscript[d, o])= 5;
(OverBar[t])  = 0.1;
(Subscript[(OverBar[d]), c])= 1.0;
(Subscript[(OverBar[d]), h])= 0.5;

(Subscript[r, o])=((Subscript[d, o])/2);
(*******************************)

(**********************************************************************)
(* Make rules for converting from normalized parameters to absolute quantities *)
(**********************************************************************)
substitutionRules = {

t      ->    (OverBar[t]) (Subscript[d, o]),

(Subscript[r, i])    ->    (Subscript[r, o])-(OverBar[t]) (Subscript[d, o]),

(Subscript[r, h])    ->    (1/2) (Subscript[(OverBar[d]), h]) (Subscript[d, o]),

(Subscript[r, c])    ->    (1/2) (Subscript[(OverBar[d]), c]) (Subscript[d, o]),

k      ->    (OverBar[k]) (Subscript[d, o]),

(Subscript[k, s])    ->    (((OverBar[k]) (Subscript[d, o]))/2),

ϕ      ->    (((Subscript[ϕ, deg]) π)/180.0),

ψ      ->    (π/6)

};
 (**********************************************************************)

(*********************************)
(* Define equations and inequalities *)
(*********************************)

(Subscript[pathOfCutAtAnySection, parametric]) = (({{(((Subscript[r, h]) Cos[ψ])/((Subscript[r, c])/ρ))
}, {(Sqrt[(ρ^2)-((((((Subscript[r, h]) Cos[ψ])/((Subscript[r, c])/ρ))))^2)])
}, {(Subscript[r, h]) Sin[ψ]}}));

(Subscript[equationOfLaserBeamAxis, implicit])=
((x^2)/(((((Subscript[r, c])/ρ)))^2))+ ((z^2)/(1^2))==(Subscript[r, h]^2) && (x^2)+(y^2)==(ρ^2);
(Subscript[equationOfLaserBeamAxis, parametric]) = (({{(((Subscript[r, h]) Cos[ψ])/((((Subscript[r, c])/ρ))))}, {(Sqrt[(ρ^2)-((((((Subscript[r, h]) Cos[ψ])/((((Subscript[r, c])/ρ)))) ))^2)])}, {(Subscript[r, h]) Sin[ψ]}}));
(*********************************)

Print["(Subscript[equationOfLaserBeamAxis, implicit]) : ", (Subscript[equationOfLaserBeamAxis, implicit])];
Print["(Subscript[equationOfLaserBeamAxis, parametric]) : ", MatrixForm[(Subscript[equationOfLaserBeamAxis, parametric])]];

(*******)
(* Plot *)
(*******)
AbsoluteTiming[

(Subscript[equationOfLaserBeamAxis, parametric-plot]) = ParametricPlot3D[
ReplaceAll[(Subscript[equationOfLaserBeamAxis, parametric]), substitutionRules],
{ρ, 0, (Subscript[r, o])},
Axes->Automatic,
AxesLabel->{"X", "Y", "Z"}
];

Quiet[
Clear[
x,y,z,θ,α,
ρ,(Subscript[ρ, n]),ψ,
(Subscript[d, o]),
(OverBar[t]),(Subscript[(OverBar[d]), h]),(Subscript[(OverBar[d]), c]),(OverBar[k]),
(Subscript[r, o]),(Subscript[r, i]),(Subscript[r, h]),(Subscript[r, c]),t,k,(Subscript[k, s]),ϕ,(Subscript[ϕ, deg])
]
];

substitutionRules = {
t->(OverBar[t]) (Subscript[d, o]),
(Subscript[r, i])->(Subscript[r, o])-(OverBar[t]) (Subscript[d, o]),
(Subscript[r, h])->(1/2) (Subscript[(OverBar[d]), h]) (Subscript[d, o]),
(Subscript[r, c])->(1/2) (Subscript[(OverBar[d]), c]) (Subscript[d, o]),
k->(OverBar[k]) (Subscript[d, o]),
(Subscript[k, s])->(((OverBar[k]) (Subscript[d, o]))/2),
ϕ -> (((Subscript[ϕ, deg]) π)/180.0)
};

(Subscript[pathOfCutAtAnySection, parametric-plot]) = ParametricPlot3D[
ReplaceAll[(Subscript[pathOfCutAtAnySection, parametric]), substitutionRules],
{ψ, 0, 2*Pi},
{ρ, ReplaceAll[(Subscript[r, i]),substitutionRules], (Subscript[r, o])},
Axes->Automatic,
AxesLabel->{"X", "Y", "Z"},
PlotPoints -> 50
];

bounds = ParametricPlot3D[
ReplaceAll[(Subscript[pathOfCutAtAnySection, parametric]), substitutionRules],
{ψ, 0, 2*Pi},
{ρ, 0, (Subscript[r, o])},
Axes->Automatic,
AxesLabel->{"X", "Y", "Z"},
Mesh->None,
PlotStyle->Directive[Orange,Opacity[0.5]]
];

Show[bounds, (Subscript[pathOfCutAtAnySection, parametric-plot]), (Subscript[equationOfLaserBeamAxis, parametric-plot])]

]

4. Locus of the normal lines

(********************************************************)
(* Clear variables so that they no longer have old, junk values *)
(********************************************************)
Quiet[
Clear[
x,y,z,θ,α,
ρ,(Subscript[ρ, n]),ψ,
(Subscript[d, o]),
(OverBar[t]),(Subscript[(OverBar[d]), h]),(Subscript[(OverBar[d]), c]),(OverBar[k]),
(Subscript[r, o]),(Subscript[r, i]),(Subscript[r, h]),(Subscript[r, c]),t,k,(Subscript[k, s]),ϕ,(Subscript[ϕ, deg])
]
];
(********************************************************)

(*******************************)
(* Assign values to the parameters *)
(*******************************)

(OverBar[k])  = 0.0029972;
(Subscript[ϕ, deg]) = 0;
(Subscript[d, o])= 5;
(OverBar[t])  = 0.1;
(Subscript[(OverBar[d]), c])= 1.0;
(Subscript[(OverBar[d]), h])= 0.5;

(Subscript[r, o])=((Subscript[d, o])/2);
(*******************************)

(**********************************************************************)
(* Make rules for converting from normalized parameters to absolute quantities *)
(**********************************************************************)
substitutionRules = {

t      ->    (OverBar[t]) (Subscript[d, o]),

(Subscript[r, i])    ->    (Subscript[r, o])-(OverBar[t]) (Subscript[d, o]),

(Subscript[r, h])    ->    (1/2) (Subscript[(OverBar[d]), h]) (Subscript[d, o]),

(Subscript[r, c])    ->    (1/2) (Subscript[(OverBar[d]), c]) (Subscript[d, o]),

k      ->    (OverBar[k]) (Subscript[d, o]),

(Subscript[k, s])    ->    (((OverBar[k]) (Subscript[d, o]))/2),

ϕ      ->    (((Subscript[ϕ, deg]) π)/180.0),

ψ      ->    (π/6)

};
 (**********************************************************************)

(*********************************)
(* Define equations and inequalities *)
(*********************************)

(Subscript[pathOfCutAtCriticalSection, parametric]) = (({{(Subscript[r, h]) Cos[ψ]
}, {(Sqrt[(Subscript[r, c]^2)-((((Subscript[r, h]) Cos[ψ]))^2)])
}, {(Subscript[r, h]) Sin[ψ]}}));

(Subscript[equationOfLaserBeamAxis, parametric]) = (({{(((Subscript[r, h]) Cos[ψ])/((((Subscript[r, c])/ρ))))}, {(Sqrt[(ρ^2)-((((((Subscript[r, h]) Cos[ψ])/((((Subscript[r, c])/ρ)))) ))^2)])}, {(Subscript[r, h]) Sin[ψ]}}));

(Subscript[equationOfNormalLines, parametric]) = (({{(((Subscript[ρ, n]) (ρ ν))/(Subscript[r, h]))
}, {((ρ (1-(((Subscript[ρ, n]) (ν^2))/(Subscript[r, h]))))/(Sqrt[1-(ν^2)]))
}, {(((Subscript[ρ, n]) (Subscript[r, h]) Sin[ψ])/(Subscript[r, h]))}}));

(Subscript[equationOfNormalLines, parametric]) = ReplaceAll[
(Subscript[equationOfNormalLines, parametric]),
 ν -> (((Subscript[r, h])/(Subscript[r, c])))Cos[ψ]
];

(*********************************)

Print["(Subscript[equationOfNormalLines, parametric]) : ", MatrixForm[(Subscript[equationOfNormalLines, parametric])]];

(*******)
(* Plot *)
(*******)
AbsoluteTiming[

(Subscript[equationOfNormalLines, parametric-plot]) = ParametricPlot3D[
ReplaceAll[(Subscript[equationOfNormalLines, parametric]), substitutionRules],
{ρ, 0, (Subscript[r, o])},
{(Subscript[ρ, n]), 0, ReplaceAll[(Subscript[r, h ]),substitutionRules]},
Axes->Automatic,
AxesLabel->{"X", "Y", "Z"},
PlotStyle->Directive[Cyan,Opacity[1]]
];

(Subscript[equationOfLaserBeamAxis, parametric-plot]) = ParametricPlot3D[
ReplaceAll[(Subscript[equationOfLaserBeamAxis, parametric]), substitutionRules],
{ρ, 0, (Subscript[r, o])},
Axes->Automatic,
AxesLabel->{"X", "Y", "Z"}
];

Quiet[
Clear[
x,y,z,θ,α,
ρ,(Subscript[ρ, n]),ψ,
(Subscript[d, o]),
(OverBar[t]),(Subscript[(OverBar[d]), h]),(Subscript[(OverBar[d]), c]),(OverBar[k]),
(Subscript[r, o]),(Subscript[r, i]),(Subscript[r, h]),(Subscript[r, c]),t,k,(Subscript[k, s]),ϕ,(Subscript[ϕ, deg])
]
];

substitutionRules = {
t->(OverBar[t]) (Subscript[d, o]),
(Subscript[r, i])->(Subscript[r, o])-(OverBar[t]) (Subscript[d, o]),
(Subscript[r, h])->(1/2) (Subscript[(OverBar[d]), h]) (Subscript[d, o]),
(Subscript[r, c])->(1/2) (Subscript[(OverBar[d]), c]) (Subscript[d, o]),
k->(OverBar[k]) (Subscript[d, o]),
(Subscript[k, s])->(((OverBar[k]) (Subscript[d, o]))/2),
ϕ -> (((Subscript[ϕ, deg]) π)/180.0)
};

(Subscript[pathOfCutAtAnySection, parametric-plot]) = ParametricPlot3D[
ReplaceAll[(Subscript[pathOfCutAtAnySection, parametric]), substitutionRules],
{ψ, 0, 2*Pi},
{ρ, ReplaceAll[(Subscript[r, i]),substitutionRules], (Subscript[r, o])},
Axes->Automatic,
AxesLabel->{"X", "Y", "Z"},
PlotPoints -> 50
];

bounds = ParametricPlot3D[
ReplaceAll[(Subscript[pathOfCutAtAnySection, parametric]), substitutionRules],
{ψ, 0, 2*Pi},
{ρ, 0, (Subscript[r, o])},
Axes->Automatic,
AxesLabel->{"X", "Y", "Z"},
Mesh->None,
PlotStyle->Directive[Orange,Opacity[1]],
PlotRange->{{ReplaceAll[-1.25 (Subscript[r, h]), substitutionRules], ReplaceAll[1.25 (Subscript[r, h]), substitutionRules]}, {0, 1.25 (Subscript[r, o])}, {ReplaceAll[-1.25 (Subscript[r, h]), substitutionRules], ReplaceAll[1.25 (Subscript[r, h]), substitutionRules]}},
PlotPoints -> 100
];

Show[bounds, (Subscript[pathOfCutAtAnySection, parametric-plot]), (Subscript[equationOfNormalLines, parametric-plot]), (Subscript[equationOfLaserBeamAxis, parametric-plot])]

]

5. Equation of the Laser Kerf Cone

(********************************************************)
(* Clear variables so that they no longer have old, junk values *)
(********************************************************)
Quiet[
Clear[
x,y,z,θ,α,
ρ,(Subscript[ρ, n]),ψ,
(Subscript[d, o]),
(OverBar[t]),(Subscript[(OverBar[d]), h]),(Subscript[(OverBar[d]), c]),(OverBar[k]),
(Subscript[r, o]),(Subscript[r, i]),(Subscript[r, h]),(Subscript[r, c]),t,k,(Subscript[k, s]),ϕ,(Subscript[ϕ, deg])
]
];
(********************************************************)

(*******************************)
(* Assign values to the parameters *)
(*******************************)

(OverBar[k])  = 0.0029972;
(Subscript[ϕ, deg]) = 0;
(Subscript[d, o])= 5;
(OverBar[t])  = 0.1;
(Subscript[(OverBar[d]), c])= 1.0;
(Subscript[(OverBar[d]), h])= 0.5;

(Subscript[r, o])=((Subscript[d, o])/2);
(*******************************)

(**********************************************************************)
(* Make rules for converting from normalized parameters to absolute quantities *)
(**********************************************************************)
substitutionRules = {

t      ->    (OverBar[t]) (Subscript[d, o]),

(Subscript[r, i])    ->    (Subscript[r, o])-(OverBar[t]) (Subscript[d, o]),

(Subscript[r, h])    ->    (1/2) (Subscript[(OverBar[d]), h]) (Subscript[d, o]),

(Subscript[r, c])    ->    (1/2) (Subscript[(OverBar[d]), c]) (Subscript[d, o]),

k      ->    (OverBar[k]) (Subscript[d, o]),

(Subscript[k, s])    ->    (((OverBar[k]) (Subscript[d, o]))/2),

ϕ      ->    (((Subscript[ϕ, deg]) π)/180.0),

ψ      ->    (π/6)

};
 (**********************************************************************)

(*********************************)
(* Define equations and inequalities *)
(*********************************)

(Subscript[pathOfCutAtAnySection, parametric]) = (({{(((Subscript[r, h]) Cos[ψ])/((Subscript[r, c])/ρ))
}, {(Sqrt[(ρ^2)-((((((Subscript[r, h]) Cos[ψ])/((Subscript[r, c])/ρ))))^2)])
}, {(Subscript[r, h]) Sin[ψ]}}));

(*The cone radius (ρ*Tan[ϕ] - (Subscript[r, c])*Tan[ϕ] + (Subscript[k, s])) is a linear interpolation from (Subscript[k, s]) at ρ=(Subscript[r, c]) to 0 at ρ=(Subscript[r, c])-((Subscript[k, s])/(Tan[ϕ]))*)
(* Parametric form *)
(Subscript[equationOfLaserCone, parametric]) = (({{(ρ Tan[ϕ]-(Subscript[r, c]) Tan[ϕ]+(Subscript[k, s])) Cos[α]}, {ρ}, {(ρ Tan[ϕ]-(Subscript[r, c]) Tan[ϕ]+(Subscript[k, s])) Sin[α]}}));

(Subscript[equationOfLaserCone, parametric]) = (Subscript[equationOfLaserCone, parametric]) + (({{0}, {0}, {(Subscript[r, h]) Sin[ψ]}}));

(Subscript[equationOfLaserCone, parametric]) = 
RotationMatrix[
-ArcTan[(Sqrt[((Subscript[r, c])^2)-((((Subscript[r, h]) Cos[ψ]))^2)]),(Subscript[r, h]) Cos[ψ]],
{0,0,1}
].(Subscript[equationOfLaserCone, parametric]);

(*********************************)

Print["(Subscript[equationOfLaserCone, parametric]) : ", MatrixForm[(Subscript[equationOfLaserCone, parametric])]];

(*******)
(* Plot *)
(*******)
AbsoluteTiming[

(Subscript[equationOfLaserCone, parametric-plot]) = ParametricPlot3D[
ReplaceAll[(Subscript[equationOfLaserCone, parametric]), substitutionRules],
{α, -Pi, Pi},
{ρ, 0, ReplaceAll[(Subscript[r, c]), substitutionRules]},
Axes->Automatic,
Mesh -> None,
PlotStyle->Directive[Magenta,Opacity[0.5]],
AxesLabel->{"X", "Y", "Z"}
];

Quiet[
Clear[
x,y,z,θ,α,
ρ,(Subscript[ρ, n]),ψ,
(Subscript[d, o]),
(OverBar[t]),(Subscript[(OverBar[d]), h]),(Subscript[(OverBar[d]), c]),(OverBar[k]),
(Subscript[r, o]),(Subscript[r, i]),(Subscript[r, h]),(Subscript[r, c]),t,k,(Subscript[k, s]),ϕ,(Subscript[ϕ, deg])
]
];

substitutionRules = {
t->(OverBar[t]) (Subscript[d, o]),
(Subscript[r, i])->(Subscript[r, o])-(OverBar[t]) (Subscript[d, o]),
(Subscript[r, h])->(1/2) (Subscript[(OverBar[d]), h]) (Subscript[d, o]),
(Subscript[r, c])->(1/2) (Subscript[(OverBar[d]), c]) (Subscript[d, o]),
k->(OverBar[k]) (Subscript[d, o]),
(Subscript[k, s])->(((OverBar[k]) (Subscript[d, o]))/2),
ϕ -> (((Subscript[ϕ, deg]) π)/180.0)
};

(Subscript[pathOfCutAtAnySection, parametric-plot]) = ParametricPlot3D[
ReplaceAll[(Subscript[pathOfCutAtAnySection, parametric]), substitutionRules],
{ψ, 0, 2*Pi},
{ρ, ReplaceAll[(Subscript[r, i]),substitutionRules], (Subscript[r, o])},
Axes->Automatic,
AxesLabel->{"X", "Y", "Z"},
PlotPoints -> 50
];

bounds = ParametricPlot3D[
ReplaceAll[(Subscript[pathOfCutAtAnySection, parametric]), substitutionRules],
{ψ, 0, 2*Pi},
{ρ, 0, (Subscript[r, o])},
Axes->Automatic,
AxesLabel->{"X", "Y", "Z"},
Mesh->None,
PlotStyle->Directive[Orange,Opacity[1]],
PlotRange->{{ReplaceAll[-1.25 (Subscript[r, h]), substitutionRules], ReplaceAll[1.25 (Subscript[r, h]), substitutionRules]}, {0, 1.25 (Subscript[r, o])}, {ReplaceAll[-1.25 (Subscript[r, h]), substitutionRules], ReplaceAll[1.25 (Subscript[r, h]), substitutionRules]}},
PlotPoints -> 100
];

Show[bounds, (Subscript[pathOfCutAtAnySection, parametric-plot]), (Subscript[equationOfNormalLines, parametric-plot]), (Subscript[equationOfLaserCone, parametric-plot])]

]

6. Pin and Socket Surfaces

(********************************************************)
(* Clear variables so that they no longer have old, junk values *)
(********************************************************)
Quiet[
Clear[
x,y,z,θ,α,
ρ,(Subscript[ρ, n]),ψ,
(Subscript[d, o]),
(OverBar[t]),(Subscript[(OverBar[d]), h]),(Subscript[(OverBar[d]), c]),(OverBar[k]),
(Subscript[r, o]),(Subscript[r, i]),(Subscript[r, h]),(Subscript[r, c]),t,k,(Subscript[k, s]),ϕ,(Subscript[ϕ, deg])
]
];
(********************************************************)

(*******************************)
(* Assign values to the parameters *)
(*******************************)

(OverBar[k])  = 0.0029972;
(Subscript[ϕ, deg]) = 0;
(Subscript[d, o])= 5;
(OverBar[t])  = 0.1;
(Subscript[(OverBar[d]), c])= 1.0;
(Subscript[(OverBar[d]), h])= 0.5;

(Subscript[r, o])=((Subscript[d, o])/2);
(*******************************)

(**********************************************************************)
(* Make rules for converting from normalized parameters to absolute quantities *)
(**********************************************************************)
substitutionRules = {

t      ->    (OverBar[t]) (Subscript[d, o]),

(Subscript[r, i])    ->    (Subscript[r, o])-(OverBar[t]) (Subscript[d, o]),

(Subscript[r, h])    ->    (1/2) (Subscript[(OverBar[d]), h]) (Subscript[d, o]),

(Subscript[r, c])    ->    (1/2) (Subscript[(OverBar[d]), c]) (Subscript[d, o]),

k      ->    (OverBar[k]) (Subscript[d, o]),

(Subscript[k, s])    ->    (((OverBar[k]) (Subscript[d, o]))/2),

ϕ      ->    (((Subscript[ϕ, deg]) π)/180.0)

};
 (**********************************************************************)

(*********************************)
(* Define equations and inequalities *)
(*********************************)

(Subscript[equationOfNormalLines, parametric]) = (({{(((Subscript[ρ, n]) (ρ ν))/(Subscript[r, h]))
}, {((ρ (1-(((Subscript[ρ, n]) (ν^2))/(Subscript[r, h]))))/(Sqrt[1-(ν^2)]))
}, {(((Subscript[ρ, n]) (Subscript[r, h]) Sin[ψ])/(Subscript[r, h]))}}));

(Subscript[equationOfNormalLines, parametric]) = ReplaceAll[
(Subscript[equationOfNormalLines, parametric]),
 ν -> (((Subscript[r, h])/(Subscript[r, c])))Cos[ψ]
];

(Subscript[equationOfLaserCone, parametric]) = (({{((ρ Cos[ψ] (Subscript[r, h])+Cos[α] (Sqrt[(Subscript[r, c]^2)-((Cos[ψ])^2) (Subscript[r, h]^2)]) ((Subscript[k, s])+(ρ-(Subscript[r, c])) Tan[ϕ]))/(Subscript[r, c]))}, {((ρ (Sqrt[(Subscript[r, c]^2)-((Cos[ψ])^2) (Subscript[r, h]^2)])-Cos[α] Cos[ψ] (Subscript[r, h]) ((Subscript[k, s])+(ρ-(Subscript[r, c])) Tan[ϕ]))/(Subscript[r, c]))}, {Sin[ψ] (Subscript[r, h])+Sin[α] ((Subscript[k, s])+(ρ-(Subscript[r, c])) Tan[ϕ])}}));

AbsoluteTiming[
(Subscript[coneNormalSolutions, parametric]) = Solve[
(Subscript[equationOfNormalLines, parametric])[[1]] == (Subscript[equationOfLaserCone, parametric])[[1]]
 && (Subscript[equationOfNormalLines, parametric])[[3]] == (Subscript[equationOfLaserCone, parametric])[[3]]
 ,{(Subscript[ρ, n]),α}
];
]
Print["(Subscript[coneNormalSolutions, parametric]) : ", MatrixForm[(Subscript[coneNormalSolutions, parametric])]];

(Subscript[equationOfSocketSurface, parametric]) = ReplaceAll[
(Subscript[equationOfNormalLines, parametric]),
 (Subscript[coneNormalSolutions, parametric])[[1]]
];
Print["(Subscript[equationOfSocketSurface, parametric]) : ", MatrixForm[(Subscript[equationOfSocketSurface, parametric])]];

(Subscript[equationOfPinSurface, parametric]) = ReplaceAll[
(Subscript[equationOfNormalLines, parametric]),
 (Subscript[coneNormalSolutions, parametric])[[3]]
];
Print["(Subscript[equationOfPinSurface, parametric]) : ", MatrixForm[(Subscript[equationOfPinSurface, parametric])]];

eq3 = ReplaceAll[
(Subscript[equationOfNormalLines, parametric]),
 (Subscript[coneNormalSolutions, parametric])[[2]]
];
Print["eq3 : ", MatrixForm[eq3]];

eq4 = ReplaceAll[
(Subscript[equationOfNormalLines, parametric]),
 (Subscript[coneNormalSolutions, parametric])[[4]]
];
Print["eq4 : ", MatrixForm[eq4]];

(*********************************)

AbsoluteTiming[

(Subscript[equationOfPinSurface, parametric-plot]) = ParametricPlot3D[
ReplaceAll[(Subscript[equationOfPinSurface, parametric]), substitutionRules],
{ρ,ReplaceAll[(Subscript[r, i]), substitutionRules], (Subscript[r, o])},
{ψ, 0, 2 π},
Axes->Automatic,
AxesLabel->{"X", "Y", "Z"},
PlotStyle->Directive[Blue,Opacity[1]],
PlotPoints->25
];

(Subscript[equationOfSocketSurface, parametric-plot]) = ParametricPlot3D[
ReplaceAll[(Subscript[equationOfSocketSurface, parametric]), substitutionRules],
{ρ,ReplaceAll[(Subscript[r, i]), substitutionRules], (Subscript[r, o])},
{ψ, 0, 2 π},
Axes->Automatic,
AxesLabel->{"X", "Y", "Z"},
PlotStyle->Directive[Orange,Opacity[1]],
PlotPoints->25
];

(*
(Subscript[eq3, parametric-plot]) = ParametricPlot3D[
ReplaceAll[eq3, substitutionRules],
{ρ,ReplaceAll[(Subscript[r, i]), substitutionRules], (Subscript[r, o])},
{ψ, Pi, 2*Pi},
Axes->Automatic,
AxesLabel->{"X", "Y", "Z"},
PlotStyle->Directive[Magenta,Opacity[1]],
PlotPoints->25
];

(Subscript[eq4, parametric-plot]) = ParametricPlot3D[
ReplaceAll[eq4, substitutionRules],
{ρ,ReplaceAll[(Subscript[r, i]), substitutionRules], (Subscript[r, o])},
{ψ, Pi, 2*Pi},
Axes->Automatic,
AxesLabel->{"X", "Y", "Z"},
PlotStyle->Directive[Green,Opacity[1]],
PlotPoints->25
];
*)

Show[(Subscript[equationOfPinSurface, parametric-plot]),(Subscript[equationOfSocketSurface, parametric-plot])
(*, (Subscript[eq3, parametric-plot]), (Subscript[eq4, parametric-plot])*)]

]

7. Rotated Pin and Socket Surfaces

(********************************************************)
(* Clear variables so that they no longer have old, junk values *)
(********************************************************)
Quiet[
Clear[
x,y,z,θ,α,
ρ,(Subscript[ρ, n]),ψ,
(Subscript[d, o]),
(OverBar[t]),(Subscript[(OverBar[d]), h]),(Subscript[(OverBar[d]), c]),(OverBar[k]),
(Subscript[r, o]),(Subscript[r, i]),(Subscript[r, h]),(Subscript[r, c]),t,k,(Subscript[k, s]),ϕ,(Subscript[ϕ, deg])
]
];
(********************************************************)

(*******************************)
(* Assign values to the parameters *)
(*******************************)

(OverBar[k])  = 0.0029972;
(Subscript[ϕ, deg]) = 0;
(Subscript[θ, deg]) = 60;
(Subscript[d, o])= 5;
(OverBar[t])  = 0.025;
(Subscript[(OverBar[d]), c])= 1.0;
(Subscript[(OverBar[d]), h])= 0.4;

(Subscript[r, o])=((Subscript[d, o])/2);
(*******************************)

(**********************************************************************)
substitutionRules = {
t      ->    (OverBar[t]) (Subscript[d, o]),
(Subscript[r, i])    ->    (Subscript[r, o])-(OverBar[t]) (Subscript[d, o]),
(Subscript[r, h])    ->    (1/2) (Subscript[(OverBar[d]), h]) (Subscript[d, o]),
(Subscript[r, c])    ->    (1/2) (Subscript[(OverBar[d]), c]) (Subscript[d, o]),
k      ->    (OverBar[k]) (Subscript[d, o]),
(Subscript[k, s])    ->    (((OverBar[k]) (Subscript[d, o]))/2),
ϕ      ->    (((Subscript[ϕ, deg]) π)/180.0),
θ -> (((Subscript[θ, deg]) π)/180.0)
};
 (**********************************************************************)

(*********************************)
(* Define equations and inequalities *)
(*********************************)

(Subscript[equationOfPinSurface, parametric]) = (({{((ρ Cos[ψ] ((Subscript[r, h])-(((Subscript[k, s]) (Sqrt[((Subscript[r, c])-Cos[ψ] (Subscript[r, h])) ((Subscript[r, c])+Cos[ψ] (Subscript[r, h]))]) (Sqrt[(Subscript[k, s])+ρ Tan[ϕ]-(Subscript[r, c]) Tan[ϕ]]))/(Sqrt[(ρ^2) ((Cos[ψ])^2) (Subscript[k, s])+((Sin[ψ])^2) (Subscript[k, s]) (Subscript[r, c]^2)-((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[k, s]) (Subscript[r, h]^2)+(ρ^3) ((Cos[ψ])^2) Tan[ϕ]-(ρ^2) ((Cos[ψ])^2) (Subscript[r, c]) Tan[ϕ]+ρ ((Sin[ψ])^2) (Subscript[r, c]^2) Tan[ϕ]-((Sin[ψ])^2) (Subscript[r, c]^3) Tan[ϕ]-ρ ((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[r, h]^2) Tan[ϕ]+((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[r, c]) (Subscript[r, h]^2) Tan[ϕ]]))-((ρ (Sqrt[((Subscript[r, c])-Cos[ψ] (Subscript[r, h])) ((Subscript[r, c])+Cos[ψ] (Subscript[r, h]))]) Tan[ϕ] (Sqrt[(Subscript[k, s])+ρ Tan[ϕ]-(Subscript[r, c]) Tan[ϕ]]))/(Sqrt[(ρ^2) ((Cos[ψ])^2) (Subscript[k, s])+((Sin[ψ])^2) (Subscript[k, s]) (Subscript[r, c]^2)-((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[k, s]) (Subscript[r, h]^2)+(ρ^3) ((Cos[ψ])^2) Tan[ϕ]-(ρ^2) ((Cos[ψ])^2) (Subscript[r, c]) Tan[ϕ]+ρ ((Sin[ψ])^2) (Subscript[r, c]^2) Tan[ϕ]-((Sin[ψ])^2) (Subscript[r, c]^3) Tan[ϕ]-ρ ((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[r, h]^2) Tan[ϕ]+((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[r, c]) (Subscript[r, h]^2) Tan[ϕ]]))+(((Subscript[r, c]) (Sqrt[((Subscript[r, c])-Cos[ψ] (Subscript[r, h])) ((Subscript[r, c])+Cos[ψ] (Subscript[r, h]))]) Tan[ϕ] (Sqrt[(Subscript[k, s])+ρ Tan[ϕ]-(Subscript[r, c]) Tan[ϕ]]))/(Sqrt[(ρ^2) ((Cos[ψ])^2) (Subscript[k, s])+((Sin[ψ])^2) (Subscript[k, s]) (Subscript[r, c]^2)-((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[k, s]) (Subscript[r, h]^2)+(ρ^3) ((Cos[ψ])^2) Tan[ϕ]-(ρ^2) ((Cos[ψ])^2) (Subscript[r, c]) Tan[ϕ]+ρ ((Sin[ψ])^2) (Subscript[r, c]^2) Tan[ϕ]-((Sin[ψ])^2) (Subscript[r, c]^3) Tan[ϕ]-ρ ((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[r, h]^2) Tan[ϕ]+((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[r, c]) (Subscript[r, h]^2) Tan[ϕ]]))))/(Subscript[r, c]))
}, {((ρ (1-((((Cos[ψ])^2) (Subscript[r, h]) ((Subscript[r, h])-(((Subscript[k, s]) (Sqrt[((Subscript[r, c])-Cos[ψ] (Subscript[r, h])) ((Subscript[r, c])+Cos[ψ] (Subscript[r, h]))]) (Sqrt[(Subscript[k, s])+ρ Tan[ϕ]-(Subscript[r, c]) Tan[ϕ]]))/(Sqrt[(ρ^2) ((Cos[ψ])^2) (Subscript[k, s])+((Sin[ψ])^2) (Subscript[k, s]) (Subscript[r, c]^2)-((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[k, s]) (Subscript[r, h]^2)+(ρ^3) ((Cos[ψ])^2) Tan[ϕ]-(ρ^2) ((Cos[ψ])^2) (Subscript[r, c]) Tan[ϕ]+ρ ((Sin[ψ])^2) (Subscript[r, c]^2) Tan[ϕ]-((Sin[ψ])^2) (Subscript[r, c]^3) Tan[ϕ]-ρ ((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[r, h]^2) Tan[ϕ]+((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[r, c]) (Subscript[r, h]^2) Tan[ϕ]]))-((ρ (Sqrt[((Subscript[r, c])-Cos[ψ] (Subscript[r, h])) ((Subscript[r, c])+Cos[ψ] (Subscript[r, h]))]) Tan[ϕ] (Sqrt[(Subscript[k, s])+ρ Tan[ϕ]-(Subscript[r, c]) Tan[ϕ]]))/(Sqrt[(ρ^2) ((Cos[ψ])^2) (Subscript[k, s])+((Sin[ψ])^2) (Subscript[k, s]) (Subscript[r, c]^2)-((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[k, s]) (Subscript[r, h]^2)+(ρ^3) ((Cos[ψ])^2) Tan[ϕ]-(ρ^2) ((Cos[ψ])^2) (Subscript[r, c]) Tan[ϕ]+ρ ((Sin[ψ])^2) (Subscript[r, c]^2) Tan[ϕ]-((Sin[ψ])^2) (Subscript[r, c]^3) Tan[ϕ]-ρ ((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[r, h]^2) Tan[ϕ]+((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[r, c]) (Subscript[r, h]^2) Tan[ϕ]]))+(((Subscript[r, c]) (Sqrt[((Subscript[r, c])-Cos[ψ] (Subscript[r, h])) ((Subscript[r, c])+Cos[ψ] (Subscript[r, h]))]) Tan[ϕ] (Sqrt[(Subscript[k, s])+ρ Tan[ϕ]-(Subscript[r, c]) Tan[ϕ]]))/(Sqrt[(ρ^2) ((Cos[ψ])^2) (Subscript[k, s])+((Sin[ψ])^2) (Subscript[k, s]) (Subscript[r, c]^2)-((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[k, s]) (Subscript[r, h]^2)+(ρ^3) ((Cos[ψ])^2) Tan[ϕ]-(ρ^2) ((Cos[ψ])^2) (Subscript[r, c]) Tan[ϕ]+ρ ((Sin[ψ])^2) (Subscript[r, c]^2) Tan[ϕ]-((Sin[ψ])^2) (Subscript[r, c]^3) Tan[ϕ]-ρ ((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[r, h]^2) Tan[ϕ]+((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[r, c]) (Subscript[r, h]^2) Tan[ϕ]]))))/(Subscript[r, c]^2))))/(Sqrt[1-((((Cos[ψ])^2) (Subscript[r, h]^2))/(Subscript[r, c]^2))]))
}, {Sin[ψ] ((Subscript[r, h])-(((Subscript[k, s]) (Sqrt[((Subscript[r, c])-Cos[ψ] (Subscript[r, h])) ((Subscript[r, c])+Cos[ψ] (Subscript[r, h]))]) (Sqrt[(Subscript[k, s])+ρ Tan[ϕ]-(Subscript[r, c]) Tan[ϕ]]))/(Sqrt[(ρ^2) ((Cos[ψ])^2) (Subscript[k, s])+((Sin[ψ])^2) (Subscript[k, s]) (Subscript[r, c]^2)-((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[k, s]) (Subscript[r, h]^2)+(ρ^3) ((Cos[ψ])^2) Tan[ϕ]-(ρ^2) ((Cos[ψ])^2) (Subscript[r, c]) Tan[ϕ]+ρ ((Sin[ψ])^2) (Subscript[r, c]^2) Tan[ϕ]-((Sin[ψ])^2) (Subscript[r, c]^3) Tan[ϕ]-ρ ((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[r, h]^2) Tan[ϕ]+((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[r, c]) (Subscript[r, h]^2) Tan[ϕ]]))-((ρ (Sqrt[((Subscript[r, c])-Cos[ψ] (Subscript[r, h])) ((Subscript[r, c])+Cos[ψ] (Subscript[r, h]))]) Tan[ϕ] (Sqrt[(Subscript[k, s])+ρ Tan[ϕ]-(Subscript[r, c]) Tan[ϕ]]))/(Sqrt[(ρ^2) ((Cos[ψ])^2) (Subscript[k, s])+((Sin[ψ])^2) (Subscript[k, s]) (Subscript[r, c]^2)-((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[k, s]) (Subscript[r, h]^2)+(ρ^3) ((Cos[ψ])^2) Tan[ϕ]-(ρ^2) ((Cos[ψ])^2) (Subscript[r, c]) Tan[ϕ]+ρ ((Sin[ψ])^2) (Subscript[r, c]^2) Tan[ϕ]-((Sin[ψ])^2) (Subscript[r, c]^3) Tan[ϕ]-ρ ((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[r, h]^2) Tan[ϕ]+((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[r, c]) (Subscript[r, h]^2) Tan[ϕ]]))+(((Subscript[r, c]) (Sqrt[((Subscript[r, c])-Cos[ψ] (Subscript[r, h])) ((Subscript[r, c])+Cos[ψ] (Subscript[r, h]))]) Tan[ϕ] (Sqrt[(Subscript[k, s])+ρ Tan[ϕ]-(Subscript[r, c]) Tan[ϕ]]))/(Sqrt[(ρ^2) ((Cos[ψ])^2) (Subscript[k, s])+((Sin[ψ])^2) (Subscript[k, s]) (Subscript[r, c]^2)-((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[k, s]) (Subscript[r, h]^2)+(ρ^3) ((Cos[ψ])^2) Tan[ϕ]-(ρ^2) ((Cos[ψ])^2) (Subscript[r, c]) Tan[ϕ]+ρ ((Sin[ψ])^2) (Subscript[r, c]^2) Tan[ϕ]-((Sin[ψ])^2) (Subscript[r, c]^3) Tan[ϕ]-ρ ((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[r, h]^2) Tan[ϕ]+((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[r, c]) (Subscript[r, h]^2) Tan[ϕ]])))}}));
(Subscript[equationOfSocketSurface, parametric]) = (({{((ρ Cos[ψ] ((Subscript[r, h])+(((Subscript[k, s]) (Sqrt[((Subscript[r, c])-Cos[ψ] (Subscript[r, h])) ((Subscript[r, c])+Cos[ψ] (Subscript[r, h]))]) (Sqrt[(Subscript[k, s])+ρ Tan[ϕ]-(Subscript[r, c]) Tan[ϕ]]))/(Sqrt[(ρ^2) ((Cos[ψ])^2) (Subscript[k, s])+((Sin[ψ])^2) (Subscript[k, s]) (Subscript[r, c]^2)-((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[k, s]) (Subscript[r, h]^2)+(ρ^3) ((Cos[ψ])^2) Tan[ϕ]-(ρ^2) ((Cos[ψ])^2) (Subscript[r, c]) Tan[ϕ]+ρ ((Sin[ψ])^2) (Subscript[r, c]^2) Tan[ϕ]-((Sin[ψ])^2) (Subscript[r, c]^3) Tan[ϕ]-ρ ((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[r, h]^2) Tan[ϕ]+((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[r, c]) (Subscript[r, h]^2) Tan[ϕ]]))+((ρ (Sqrt[((Subscript[r, c])-Cos[ψ] (Subscript[r, h])) ((Subscript[r, c])+Cos[ψ] (Subscript[r, h]))]) Tan[ϕ] (Sqrt[(Subscript[k, s])+ρ Tan[ϕ]-(Subscript[r, c]) Tan[ϕ]]))/(Sqrt[(ρ^2) ((Cos[ψ])^2) (Subscript[k, s])+((Sin[ψ])^2) (Subscript[k, s]) (Subscript[r, c]^2)-((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[k, s]) (Subscript[r, h]^2)+(ρ^3) ((Cos[ψ])^2) Tan[ϕ]-(ρ^2) ((Cos[ψ])^2) (Subscript[r, c]) Tan[ϕ]+ρ ((Sin[ψ])^2) (Subscript[r, c]^2) Tan[ϕ]-((Sin[ψ])^2) (Subscript[r, c]^3) Tan[ϕ]-ρ ((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[r, h]^2) Tan[ϕ]+((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[r, c]) (Subscript[r, h]^2) Tan[ϕ]]))-(((Subscript[r, c]) (Sqrt[((Subscript[r, c])-Cos[ψ] (Subscript[r, h])) ((Subscript[r, c])+Cos[ψ] (Subscript[r, h]))]) Tan[ϕ] (Sqrt[(Subscript[k, s])+ρ Tan[ϕ]-(Subscript[r, c]) Tan[ϕ]]))/(Sqrt[(ρ^2) ((Cos[ψ])^2) (Subscript[k, s])+((Sin[ψ])^2) (Subscript[k, s]) (Subscript[r, c]^2)-((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[k, s]) (Subscript[r, h]^2)+(ρ^3) ((Cos[ψ])^2) Tan[ϕ]-(ρ^2) ((Cos[ψ])^2) (Subscript[r, c]) Tan[ϕ]+ρ ((Sin[ψ])^2) (Subscript[r, c]^2) Tan[ϕ]-((Sin[ψ])^2) (Subscript[r, c]^3) Tan[ϕ]-ρ ((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[r, h]^2) Tan[ϕ]+((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[r, c]) (Subscript[r, h]^2) Tan[ϕ]]))))/(Subscript[r, c]))}, {((ρ (1-((((Cos[ψ])^2) (Subscript[r, h]) ((Subscript[r, h])+(((Subscript[k, s]) (Sqrt[((Subscript[r, c])-Cos[ψ] (Subscript[r, h])) ((Subscript[r, c])+Cos[ψ] (Subscript[r, h]))]) (Sqrt[(Subscript[k, s])+ρ Tan[ϕ]-(Subscript[r, c]) Tan[ϕ]]))/(Sqrt[(ρ^2) ((Cos[ψ])^2) (Subscript[k, s])+((Sin[ψ])^2) (Subscript[k, s]) (Subscript[r, c]^2)-((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[k, s]) (Subscript[r, h]^2)+(ρ^3) ((Cos[ψ])^2) Tan[ϕ]-(ρ^2) ((Cos[ψ])^2) (Subscript[r, c]) Tan[ϕ]+ρ ((Sin[ψ])^2) (Subscript[r, c]^2) Tan[ϕ]-((Sin[ψ])^2) (Subscript[r, c]^3) Tan[ϕ]-ρ ((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[r, h]^2) Tan[ϕ]+((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[r, c]) (Subscript[r, h]^2) Tan[ϕ]]))+((ρ (Sqrt[((Subscript[r, c])-Cos[ψ] (Subscript[r, h])) ((Subscript[r, c])+Cos[ψ] (Subscript[r, h]))]) Tan[ϕ] (Sqrt[(Subscript[k, s])+ρ Tan[ϕ]-(Subscript[r, c]) Tan[ϕ]]))/(Sqrt[(ρ^2) ((Cos[ψ])^2) (Subscript[k, s])+((Sin[ψ])^2) (Subscript[k, s]) (Subscript[r, c]^2)-((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[k, s]) (Subscript[r, h]^2)+(ρ^3) ((Cos[ψ])^2) Tan[ϕ]-(ρ^2) ((Cos[ψ])^2) (Subscript[r, c]) Tan[ϕ]+ρ ((Sin[ψ])^2) (Subscript[r, c]^2) Tan[ϕ]-((Sin[ψ])^2) (Subscript[r, c]^3) Tan[ϕ]-ρ ((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[r, h]^2) Tan[ϕ]+((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[r, c]) (Subscript[r, h]^2) Tan[ϕ]]))-(((Subscript[r, c]) (Sqrt[((Subscript[r, c])-Cos[ψ] (Subscript[r, h])) ((Subscript[r, c])+Cos[ψ] (Subscript[r, h]))]) Tan[ϕ] (Sqrt[(Subscript[k, s])+ρ Tan[ϕ]-(Subscript[r, c]) Tan[ϕ]]))/(Sqrt[(ρ^2) ((Cos[ψ])^2) (Subscript[k, s])+((Sin[ψ])^2) (Subscript[k, s]) (Subscript[r, c]^2)-((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[k, s]) (Subscript[r, h]^2)+(ρ^3) ((Cos[ψ])^2) Tan[ϕ]-(ρ^2) ((Cos[ψ])^2) (Subscript[r, c]) Tan[ϕ]+ρ ((Sin[ψ])^2) (Subscript[r, c]^2) Tan[ϕ]-((Sin[ψ])^2) (Subscript[r, c]^3) Tan[ϕ]-ρ ((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[r, h]^2) Tan[ϕ]+((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[r, c]) (Subscript[r, h]^2) Tan[ϕ]]))))/(Subscript[r, c]^2))))/(Sqrt[1-((((Cos[ψ])^2) (Subscript[r, h]^2))/(Subscript[r, c]^2))]))}, {Sin[ψ] ((Subscript[r, h])+(((Subscript[k, s]) (Sqrt[((Subscript[r, c])-Cos[ψ] (Subscript[r, h])) ((Subscript[r, c])+Cos[ψ] (Subscript[r, h]))]) (Sqrt[(Subscript[k, s])+ρ Tan[ϕ]-(Subscript[r, c]) Tan[ϕ]]))/(Sqrt[(ρ^2) ((Cos[ψ])^2) (Subscript[k, s])+((Sin[ψ])^2) (Subscript[k, s]) (Subscript[r, c]^2)-((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[k, s]) (Subscript[r, h]^2)+(ρ^3) ((Cos[ψ])^2) Tan[ϕ]-(ρ^2) ((Cos[ψ])^2) (Subscript[r, c]) Tan[ϕ]+ρ ((Sin[ψ])^2) (Subscript[r, c]^2) Tan[ϕ]-((Sin[ψ])^2) (Subscript[r, c]^3) Tan[ϕ]-ρ ((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[r, h]^2) Tan[ϕ]+((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[r, c]) (Subscript[r, h]^2) Tan[ϕ]]))+((ρ (Sqrt[((Subscript[r, c])-Cos[ψ] (Subscript[r, h])) ((Subscript[r, c])+Cos[ψ] (Subscript[r, h]))]) Tan[ϕ] (Sqrt[(Subscript[k, s])+ρ Tan[ϕ]-(Subscript[r, c]) Tan[ϕ]]))/(Sqrt[(ρ^2) ((Cos[ψ])^2) (Subscript[k, s])+((Sin[ψ])^2) (Subscript[k, s]) (Subscript[r, c]^2)-((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[k, s]) (Subscript[r, h]^2)+(ρ^3) ((Cos[ψ])^2) Tan[ϕ]-(ρ^2) ((Cos[ψ])^2) (Subscript[r, c]) Tan[ϕ]+ρ ((Sin[ψ])^2) (Subscript[r, c]^2) Tan[ϕ]-((Sin[ψ])^2) (Subscript[r, c]^3) Tan[ϕ]-ρ ((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[r, h]^2) Tan[ϕ]+((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[r, c]) (Subscript[r, h]^2) Tan[ϕ]]))-(((Subscript[r, c]) (Sqrt[((Subscript[r, c])-Cos[ψ] (Subscript[r, h])) ((Subscript[r, c])+Cos[ψ] (Subscript[r, h]))]) Tan[ϕ] (Sqrt[(Subscript[k, s])+ρ Tan[ϕ]-(Subscript[r, c]) Tan[ϕ]]))/(Sqrt[(ρ^2) ((Cos[ψ])^2) (Subscript[k, s])+((Sin[ψ])^2) (Subscript[k, s]) (Subscript[r, c]^2)-((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[k, s]) (Subscript[r, h]^2)+(ρ^3) ((Cos[ψ])^2) Tan[ϕ]-(ρ^2) ((Cos[ψ])^2) (Subscript[r, c]) Tan[ϕ]+ρ ((Sin[ψ])^2) (Subscript[r, c]^2) Tan[ϕ]-((Sin[ψ])^2) (Subscript[r, c]^3) Tan[ϕ]-ρ ((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[r, h]^2) Tan[ϕ]+((Cos[ψ])^2) ((Sin[ψ])^2) (Subscript[r, c]) (Subscript[r, h]^2) Tan[ϕ]])))}}));

(Subscript[rotatedEquationOfPinSurface, parametric]) = RotationMatrix[θ, {0, 1, 0}].(Subscript[equationOfPinSurface, parametric]);
Print[MatrixForm[ReplaceAll[(Subscript[rotatedEquationOfPinSurface, parametric]), substitutionRules]]];
(*********************************)

(********)
(* Plot *)
(********)
AbsoluteTiming[
(* Test Plot *)
(Subscript[rotatedEquationOfPinSurface, parametric-plot]) = ParametricPlot3D[
ReplaceAll[(Subscript[rotatedEquationOfPinSurface, parametric]), substitutionRules],
{ρ,ReplaceAll[(Subscript[r, i]), substitutionRules], (Subscript[r, o])},
{ψ, 0, 2*Pi},
Axes->Automatic,
AxesLabel->{"X", "Y", "Z"},
PlotStyle->Directive[Blue,Opacity[1]],
PlotPoints->50,
Mesh->None
];
(Subscript[equationOfSocketSurface, parametric-plot]) = ParametricPlot3D[
ReplaceAll[(Subscript[equationOfSocketSurface, parametric]), substitutionRules],
{ρ,ReplaceAll[(Subscript[r, i]), substitutionRules], (Subscript[r, o])},
{ψ, 0, 2*Pi},
Axes->Automatic,
AxesLabel->{"X", "Y", "Z"},
PlotStyle->Directive[Orange,Opacity[1]],
PlotPoints->50,
Mesh->None
];
Show[(Subscript[rotatedEquationOfPinSurface, parametric-plot]),(Subscript[equationOfSocketSurface, parametric-plot])]
]

8. Trials

SetOptions[ParametricPlot3D,ImageSize->{1000,1000}];
(********************************************************)
(* Clear variables so that they no longer have old, junk values *)
(********************************************************)
Quiet[
Clear[
x,y,z,θ,α,
ρ,(Subscript[ρ, n]),ψ,
(Subscript[d, o]),
(OverBar[t]),(Subscript[(OverBar[d]), h]),(Subscript[(OverBar[d]), c]),(OverBar[k]),
(Subscript[r, o]),(Subscript[r, i]),(Subscript[r, h]),(Subscript[r, c]),t,k,(Subscript[k, s]),ϕ,(Subscript[ϕ, deg])
]
];
(********************************************************)

(*******************************)

(* Assign values to the parameters *)
(*******************************)

(Subscript[d, o])= 5;
(Subscript[(OverBar[d]), c])= 1.0;
(Subscript[ϕ, deg]) = 0;
(OverBar[k])  = 0.00381;
(Subscript[(OverBar[d]), h])= 0.40;
(OverBar[t])  = 0.1;
(Subscript[θ, deg]) = 510;

(Subscript[r, o])=((Subscript[d, o])/2);
(*******************************)

(**********************************************************************)
substitutionRules = {
t      ->    (OverBar[t]) (Subscript[d, o]),
(Subscript[r, i])    ->    (Subscript[r, o])-(OverBar[t]) (Subscript[d, o]),
(Subscript[r, h])    ->    (1/2) (Subscript[(OverBar[d]), h]) (Subscript[d, o]),
(Subscript[r, c])    ->    (1/2) (Subscript[(OverBar[d]), c]) (Subscript[d, o]),
k      ->    (OverBar[k]) (Subscript[d, o]),
(Subscript[k, s])    ->    (((OverBar[k]) (Subscript[d, o]))/2),
ϕ      ->    (((Subscript[ϕ, deg]) π)/180.0),
θ -> (((Subscript[θ, deg]) π)/180.0)
};
 (**********************************************************************)

(*********************************)
(* Define equations and inequalities *)
(*********************************)
 (**********************************************************************)

(*
parameterEliminationRules = {
ρ->(Sqrt[(x^2)+(y^2)]),
ψ->ArcTan[(Subscript[r, c])x,z (Sqrt[(x^2)+(y^2)])]
};

(Subscript[equationOfPin, implicit]) = z == ReplaceAll[(Subscript[equationOfPinSurface, parametric])[[3]], parameterEliminationRules];
AbsoluteTiming[
(Subscript[equationOfPin, implicit]) = FullSimplify[(Subscript[equationOfPin, implicit])];
]

Print["(Subscript[equationOfPin, implicit]): ", (Subscript[equationOfPin, implicit])];
Print["(Subscript[equationOfPin, implicit]): ", ReplaceAll[(Subscript[equationOfPin, implicit]), substitutionRules]];

(Subscript[equationOfSocket, implicit]) = z == ReplaceAll[(Subscript[equationOfSocketSurface, parametric])[[3]], parameterEliminationRules];
AbsoluteTiming[
(Subscript[equationOfSocket, implicit]) = FullSimplify[(Subscript[equationOfSocket, implicit])];
]

Print["(Subscript[equationOfSocket, implicit]): ", (Subscript[equationOfSocket, implicit])];
Print["(Subscript[equationOfSocket, implicit]): ", ReplaceAll[(Subscript[equationOfSocket, implicit]), substitutionRules]];
*)

Needs["NDSolve`FEM`"]
socketInequality=1>(((Subscript[r, h]) (Sqrt[(x^2)+(y^2)]))/(Sqrt[((x^2)+(y^2)) (z^2)+(x^2) (Subscript[r, c]^2)]))+(((Sqrt[((x^2)+(y^2)) (z^2)+(x^2) (Subscript[r, c]^2)-(x^2) (Subscript[r, h]^2)]) ((Subscript[k, s])+((Sqrt[(x^2)+(y^2)])-(Subscript[r, c])) Tan[ϕ]))/(Sqrt[((x^2)+(z^2)) (((x^2)+(y^2)) (z^2)+(x^2) (Subscript[r, c]^2))-(x^2) (z^2) (Subscript[r, h]^2)]))&&(Subscript[r, i]^2)≤(x^2)+(y^2)≤(Subscript[r, o]^2);
pinInequality=1<(((Subscript[r, h]) (Sqrt[(x^2)+(y^2)]))/(Sqrt[((x^2)+(y^2)) (z^2)+(x^2) (Subscript[r, c]^2)]))-(((Sqrt[((x^2)+(y^2)) (z^2)+(x^2) (Subscript[r, c]^2)-(x^2) (Subscript[r, h]^2)]) ((Subscript[k, s])+((Sqrt[(x^2)+(y^2)])-(Subscript[r, c])) Tan[ϕ]))/(Sqrt[((x^2)+(z^2)) (((x^2)+(y^2)) (z^2)+(x^2) (Subscript[r, c]^2))-(x^2) (z^2) (Subscript[r, h]^2)]))&&(Subscript[r, i]^2)≤(x^2)+(y^2)≤(Subscript[r, o]^2);

                                                                                                          
socketInequalitySub = ReplaceAll[socketInequality, substitutionRules];
pinInequalitySub = ReplaceAll[pinInequality, substitutionRules];
socketRegion = ImplicitRegion[socketInequalitySub,{{x,-2.5,2.5},{y,0,2.5},{z,-2.75,2.75}}];
pinRegion = ImplicitRegion[pinInequalitySub,{{x,-2.5,2.5},{y,0,2.5},{z,-2.75,2.75}}];
rotatedPinRegion = ImplicitRegion[pinInequalitySub/.Thread[{x,y,z}->RotationTransform[-(Subscript[θ, deg]) * Pi/180.0,{0,1,0}][{x,y,z}]]//Evaluate,{{x,-2.5,2.5},{y,0,2.5},{z,-2.75,2.75}}];
commonRegion=RegionIntersection[socketRegion,rotatedPinRegion];
discretizedCommonRegion1 = DiscretizeRegion[commonRegion, MaxCellMeasure->{"Area" -> 0.001}]
discretizedCommonRegion2 = MeshRegion[ToElementMesh[commonRegion, MaxCellMeasure->{"Area" -> 0.001}]]
Volume[discretizedCommonRegion1]
Volume[discretizedCommonRegion2]
RegionMeasure[discretizedCommonRegion2]

(*********************************)

9. Pin and Socket Spaces

SetOptions[ParametricPlot3D,ImageSize->{1000,1000}];
(********************************************************)
(* Clear variables so that they no longer have old, junk values *)
(********************************************************)
Quiet[
Clear[
x,y,z,θ,α,
ρ,(Subscript[ρ, n]),ψ,
(Subscript[d, o]),
(OverBar[t]),(Subscript[(OverBar[d]), h]),(Subscript[(OverBar[d]), c]),(OverBar[k]),
(Subscript[r, o]),(Subscript[r, i]),(Subscript[r, h]),(Subscript[r, c]),t,k,(Subscript[k, s]),ϕ,(Subscript[ϕ, deg])
]
];
(********************************************************)

(*******************************)
(* Assign values to the parameters *)
(*******************************)

(OverBar[k])  = 0.0029972;
(Subscript[ϕ, deg]) = 0;
(Subscript[θ, deg]) = 45;
(Subscript[d, o])= 5;
(OverBar[t])  = 0.1;
(Subscript[(OverBar[d]), c])= 1.0;
(Subscript[(OverBar[d]), h])= 0.5;

(Subscript[r, o])=((Subscript[d, o])/2);
(*******************************)

(**********************************************************************)
substitutionRules = {

t      ->    (OverBar[t]) (Subscript[d, o]),

(Subscript[r, i])    ->    (Subscript[r, o])-(OverBar[t]) (Subscript[d, o]),

(Subscript[r, h])    ->    (1/2) (Subscript[(OverBar[d]), h]) (Subscript[d, o]),

(Subscript[r, c])    ->    (1/2) (Subscript[(OverBar[d]), c]) (Subscript[d, o]),

k      ->    (OverBar[k]) (Subscript[d, o]),

(Subscript[k, s])    ->    (((OverBar[k]) (Subscript[d, o]))/2),

ϕ      ->    (((Subscript[ϕ, deg]) π)/180.0),


θ -> (((Subscript[θ, deg]) π)/180.0)

};
 (**********************************************************************)

(*
parameterEliminationRules = {
ρ->(Sqrt[(x^2)+(y^2)]),
ψ->ArcTan[(Subscript[r, c])x,z (Sqrt[(x^2)+(y^2)])]
};

(Subscript[equationOfPin, implicit]) = z == ReplaceAll[(Subscript[equationOfPinSurface, parametric])[[3]], parameterEliminationRules];
AbsoluteTiming[
(Subscript[equationOfPin, implicit]) = FullSimplify[(Subscript[equationOfPin, implicit])];
]

Print["(Subscript[equationOfPin, implicit]): ", (Subscript[equationOfPin, implicit])];
Print["(Subscript[equationOfPin, implicit]): ", ReplaceAll[(Subscript[equationOfPin, implicit]), substitutionRules]];

(Subscript[equationOfSocket, implicit]) = z == ReplaceAll[(Subscript[equationOfSocketSurface, parametric])[[3]], parameterEliminationRules];
AbsoluteTiming[
(Subscript[equationOfSocket, implicit]) = FullSimplify[(Subscript[equationOfSocket, implicit])];
]

Print["(Subscript[equationOfSocket, implicit]): ", (Subscript[equationOfSocket, implicit])];
Print["(Subscript[equationOfSocket, implicit]): ", ReplaceAll[(Subscript[equationOfSocket, implicit]), substitutionRules]];
*)
socketInequality=1>(((Subscript[r, h]) (Sqrt[(x^2)+(y^2)]))/(Sqrt[((x^2)+(y^2)) (z^2)+(x^2) (Subscript[r, c]^2)]))+(((Sqrt[((x^2)+(y^2)) (z^2)+(x^2) (Subscript[r, c]^2)-(x^2) (Subscript[r, h]^2)]) ((Subscript[k, s])+((Sqrt[(x^2)+(y^2)])-(Subscript[r, c])) Tan[ϕ]))/(Sqrt[((x^2)+(z^2)) (((x^2)+(y^2)) (z^2)+(x^2) (Subscript[r, c]^2))-(x^2) (z^2) (Subscript[r, h]^2)]))&&(Subscript[r, i]^2)≤(x^2)+(y^2)≤(Subscript[r, o]^2);
pinInequality=1<(((Subscript[r, h]) (Sqrt[(x^2)+(y^2)]))/(Sqrt[((x^2)+(y^2)) (z^2)+(x^2) (Subscript[r, c]^2)]))-(((Sqrt[((x^2)+(y^2)) (z^2)+(x^2) (Subscript[r, c]^2)-(x^2) (Subscript[r, h]^2)]) ((Subscript[k, s])+((Sqrt[(x^2)+(y^2)])-(Subscript[r, c])) Tan[ϕ]))/(Sqrt[((x^2)+(z^2)) (((x^2)+(y^2)) (z^2)+(x^2) (Subscript[r, c]^2))-(x^2) (z^2) (Subscript[r, h]^2)]))&&(Subscript[r, i]^2)≤(x^2)+(y^2)≤(Subscript[r, o]^2);
(*
socketInequalitySub = ReplaceAll[socketInequality, substitutionRules];
pinInequalitySub = ReplaceAll[pinInequality, substitutionRules];
socketRegion = ImplicitRegion[socketInequalitySub,{{x,-2.5,2.5},{y,0,2.5},{z,-2.75,2.75}}];
pinRegion = ImplicitRegion[pinInequalitySub,{{x,-2.5,2.5},{y,0,2.5},{z,-2.75,2.75}}];
rotatedPinRegion = ImplicitRegion[pinInequalitySub/.Thread[{x,y,z}->RotationTransform[-(Subscript[θ, deg]),{0,1,0}][{x,y,z}]]//Evaluate,{{x,-2.5,2.5},{y,0,2.5},{z,-2.75,2.75}}];
commonRegion=RegionIntersection[socketRegion,rotatedPinRegion];
(*AbsoluteTiming[discretizedCommonRegion = DiscretizeRegion[commonRegion, MaxCellMeasure->{"Area" -> 0.001}]]*)
AbsoluteTiming[discretizedCommonRegion = DiscretizeRegion[rotatedPinRegion]]
(*volume = Volume[discretizedCommonRegion];*)
*)

m1 = Table[
NaN (* Default value*),
 4 * 7  (* Number or rows = Total number of iterations *),
{3} (* 1 column each for storing (Subscript[(OverBar[d]), h]), (OverBar[k]), θ, and volume*)
];
executionTimes = Table[
0,
4 * 7 * 180,
{1}
];

iterationCount = 1;
execCount = 1;
lastExecTime = 0;
averageExecTime = 0;
lastRow = {0,0,0,0};
(Subscript[(OverBar[d]), h-set]) = {0.3, 0.4, 0.5, 0.75};
(Subscript[(OverBar[d]), h-counter]) = 0;
(Subscript[θ, deg]) = 0;
progressPanel = Panel[
Column[{
Row[{"\!\(\*SubscriptBox[OverscriptBox[\(d\), \(_\)], \(\(h\)\(\\ \\ \\ \\ \\ \\ \)\)]\)",Spacer[20],ProgressIndicator[Dynamic[(Subscript[(OverBar[d]), h-counter])],{1, 4}], Spacer[20],Dynamic[(Subscript[(OverBar[d]), h-counter]) / 4.0* 100.0], "%"}],
Row[{"(OverBar[k])     ",Spacer[20],ProgressIndicator[Dynamic[(OverBar[k])],{0.001, 0.020}], Spacer[20],Dynamic[(OverBar[k]) / 0.07 * 100.0], "%"}],
Row[{"(Subscript[θ, deg])",Spacer[20],ProgressIndicator[Dynamic[(Subscript[θ, deg])],{0, 180}], Spacer[20],Dynamic[(Subscript[θ, deg]) / 90.0 * 100.0], "%"}],
Row[{"Last Iteration Time (s): ",Spacer[20],Dynamic[lastExecTime]}],
Row[{"Average Iteration Time (s): ",Spacer[20],Dynamic[averageExecTime]}],
Row[{"Worst-case Estimated Time (s): ",Spacer[20],Dynamic[averageExecTime * (4 * 7 * 90 - execCount)]}],
Row[{"Last Row: ",Spacer[20],Dynamic[lastRow]}],
Row[{MatrixForm[Dynamic[m1]]}]
}],
"Progress"
]
execCount = 1;
AbsoluteTiming[
For[(Subscript[(OverBar[d]), h-counter]) = 1, (Subscript[(OverBar[d]), h-counter]) ≤ Length[(Subscript[(OverBar[d]), h-set])], (Subscript[(OverBar[d]), h-counter]) =(Subscript[(OverBar[d]), h-counter]) + 1,
(Subscript[(OverBar[d]), h]) = (Subscript[(OverBar[d]), h-set])[[(Subscript[(OverBar[d]), h-counter])]];
For[ (OverBar[k])  =  0.001, (OverBar[k]) ≤ 0.007, (OverBar[k]) = (OverBar[k]) + 0.001,
substitutionRules = {
t      ->    (OverBar[t]) (Subscript[d, o]),
(Subscript[r, i])    ->    (Subscript[r, o])-(OverBar[t]) (Subscript[d, o]),
(Subscript[r, h])    ->    (1/2) (Subscript[(OverBar[d]), h]) (Subscript[d, o]),
(Subscript[r, c])    ->    (1/2) (Subscript[(OverBar[d]), c]) (Subscript[d, o]),
k      ->    (OverBar[k]) (Subscript[d, o]),
(Subscript[k, s])    ->    (((OverBar[k]) (Subscript[d, o]))/2),
ϕ      ->    (((Subscript[ϕ, deg]) π)/180.0),
θ      ->    (((Subscript[θ, deg]) π)/180.0)
};
socketInequalitySub = ReplaceAll[socketInequality, substitutionRules];
pinInequalitySub = ReplaceAll[pinInequality, substitutionRules];
socketRegion = ImplicitRegion[socketInequalitySub,{{x,-2.5,2.5},{y,0,2.5},{z,-2.75,2.75}}];
(*pinRegion = ImplicitRegion[pinInequalitySub,{{x,-2.5,2.5},{y,0,2.5},{z,-2.75,2.75}}];*)
For[(Subscript[θ, deg]) =  0.1, (Subscript[θ, deg]) <= 90, (Subscript[θ, deg]) = (Subscript[θ, deg]) + 0.1, 
executionTimes[[execCount, 1]] = First[AbsoluteTiming[
volume = 0;
rotatedPinRegion = ImplicitRegion[pinInequalitySub/.Thread[{x,y,z}->RotationTransform[-(Subscript[θ, deg]) * Pi/180.0,{0,1,0}][{x,y,z}]]//Evaluate,{{x,-2.5,2.5},{y,0,2.5},{z,-2.75,2.75}}];
commonRegion=RegionIntersection[socketRegion,rotatedPinRegion];
(*discretizedCommonRegion = DiscretizeRegion[commonRegion, MaxCellMeasure->{"Area" -> 0.001}];*)
discretizedCommonRegion = MeshRegion[ToElementMesh[commonRegion, MaxCellMeasure->{"Area" -> 0.001}]];
volume = Volume[discretizedCommonRegion];
]];
lastExecTime = executionTimes[[execCount, 1]];
averageExecTime = Mean[executionTimes[[1;;execCount,1]]];
execCount = execCount + 1;
If[ !NumericQ[volume],
lastRow = { (Subscript[(OverBar[d]), h]), (OverBar[k]), (Subscript[θ, deg]), 0};
Continue[];
];
If[volume > 0.0001,
m1[[iterationCount, 1]] = (Subscript[(OverBar[d]), h]);
m1[[iterationCount, 2]] = (OverBar[k]);
m1[[iterationCount, 3]] = (Subscript[θ, deg]);
lastRow = m1[[iterationCount,  All]];
iterationCount = iterationCount + 1;
(*Print[MatrixForm[m1]];*)
Break[];
];                         
];
If[m1[[iterationCount - 1, 3]] ≥ 90, 
Break[];
];
];
];
]

AbsoluteTiming[
Plot1 = ListPlot[
Table[m1[[(i-1)*7+1;;i*7, 2;;3]], {i, 1, 4, 1}],
AxesLabel->{"Normalized Kerf Width","θ (in degrees)"},
Joined->False,
InterpolationOrder->1,
 PlotLegends-> Placed[LineLegend[Table[m1[[(j-1)*7+1,1]],{j, 1, 4, 1}], LegendLabel-> "Normalized Hinge Size"], Right],
 BaseStyle->{FontWeight->"Bold",FontSize->22}
];
]
Show[Plot1]