British Museum Great Court

def z (x:var[]..[], y:var[]..[])
{
a = 22.245;
b = 36.625;
c = 46.025;
d = 51.125;
​
//λ
lm = 0.5;
​
//μ
mu = 14.0;
​
hcen = 20.955;
hedg = 19.71;
r = Math.Sqrt(Math.Pow(x,2)+Math.Pow(y,2));
​
//θ
the = Math.Asin(y/r);
​
//η
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));
​
//α
alp = ((r/a)-1) * psi;
​
//β
bet = (1-(a/r))/
(((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))))));
return = z1 + z2 +z3;
};
​
x =-36.6..36.6..#30;
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());

cen1 = Circle.ByCenterPointRadius(Point.Origin(),15);
edg1 = Rectangle.ByWidthLength(73,97).Translate(Vector.YAxis(),3);
pnt1 = (edg1.Explode())<1>.PointAtParameter((0..1..#10)<2>);
pnt2 = cen1.ClosestPointTo(pnt1);
pnt3 = Line.ByStartPointEndPoint(pnt1,pnt2).PointAtParameter(0.5).Translate(Vector.ZAxis(),6);
pnt4 = List.Transpose(Arc.ByThreePoints(pnt1,pnt3,pnt2))<1>.PointAtParameter((0..1..#10)<2>);
pnt5 = Point.PruneDuplicates(List.Flatten(pnt4<1>,-1),0.001);
srf1 = Surface.ByLoft(List.AddItemToFront(edg1,List.RestOfItems(NurbsCurve.ByControlPoints(pnt5,3,true))));