def z (x:var[]..[], y:var[]..[])
r = Math.Sqrt(Math.Pow(x,2)+Math.Pow(y,2));
eta = (1-(x/b)) * (1+(x/b)) * (1-(y/c)) * (1+(y/d))/
(1-((a*x)/(r*b))) * (1+((a*x)/(r*b))) * (1-((a*y)/(r*c))) * (1+((a*y)/(r*d)));
psi = (1-(x/b)) * (1+(x/b)) * (1-(y/c)) * (1+(y/d));
(((Math.Sqrt(Math.Pow(b-x,2) + Math.Pow(c-y,2))) / ((b-x)*(c-y))) +
((Math.Sqrt(Math.Pow(b-x,2) + Math.Pow(d+y,2))) / ((b-x) * (d+y))) +
((Math.Sqrt(Math.Pow(b+x,2) + Math.Pow(c-y,2))) / ((b+x) * (c-y))) +
((Math.Sqrt(Math.Pow(b+x,2) + Math.Pow(d+y,2))) / ((b+x)*(d+y))));
z1 = ((hcen - hedg)*eta) + hedg;
z2 = alp * (((1-lm)*(((35.0+(10.0*psi))*0.5*(1+Math.Cos(2*the)))+(12.0*(0.5*(1-Math.Cos(2*the))+Math.Sin(the)))+((7.5+(12.0*psi))*(0.5*(1-Math.Cos(2*the))-Math.Sin(the)))-1.6))
- ((5.0*(1+Math.Cos(2*the)))+(10.0*(Math.Pow((0.5*(0.5*(1-Math.Cos(2*the))+Math.Sin(the))),2))*(1.0-(3.0*alp))))
+ (2.5*(Math.Pow((0.5*(0.5*(1-Math.Cos(2*the)-Math.Sin(the)))),2))*(Math.Pow(((r/a)-1),2))));
z3 = bet * ((lm*((1.75*(1+Math.Cos(2*the)))+(1.5*(1-Math.Cos(2*the)))+(0.3*Math.Sin(the))))
+ (1.05*(Math.Pow(Math.E,-mu*(1-(x/b)))+Math.Pow(Math.E,-mu*(1+(x/b))))*(Math.Pow(Math.E,-mu*(1-(y/c)))+Math.Pow(Math.E,-mu*(1+(y/d))))));
y = -51.125..46.025..#30;
p = NurbsSurface.ByPoints(Point.ByCoordinates(x<1>,y<2>).Translate(Vector.ByCoordinates(0,0,z(x<1>,y<2>)))).Trim(Plane.ByOriginNormal(Point.ByCoordinates(0,0,19.71),Vector.ZAxis()),Point.Origin());