Делим четырехугольник на дачные участки. 1

Еще в 50-годы моя мама при дачном правлении состояла в Земельной Комиссии. В те времена могли выделять участки для учителей Москвы. Делалось так: колышками обозначались вершины четырехугольника с таким расчетом, чтобы площадь оказывалась 40 соток. В точности прямоугольник часто не получался, так как геология не позволяла (овраги, дороги, ручьи и т.д.). Известны координаты вершин (x1,y1),....,(x4,y4). Нужно было геометрически или алгебраически наложить на чертеж крестовину (она голубого цвета) так, чтобы все четыре участка оказались площадью 10 соток. Найти нужно координаты пяти точек, которые обозначены буквами в кружочках. Мама попросила меня, уже школьника, такую задачу решить. Решил ее или нет - хоть убейте, не помню. Но в архиве мамы этот чертеж нашел, и попытался два дня назад с задачей справиться. Весь вчерашний день осмысливал ходы рассуждений, понял, что теоретически мне не одолеть. Поэтому стал писать программу для компа, чтобы получить хоть какой-нибудь результат методом Монте Карло. Параллельно обратился к математическому форуму. Скажу честно - мне было сильно трудно, хотел даже бросить всё. Но ради мамы, ради ее памяти... В конце-концов уже поздно вечером буквально проломил сложности и получил великолепный инструмент для расчета координат крестовины внутри четырехугольника. Пример дан на рисунке в центре. А спустя всего час мой давний помощник по форуму под ником michel дал свое решение! При помощи  Mathcad сделал и всё-всё совпало с моими потугами! О его результатах расскажу во второй части.

Текст моей проги:

rem делим четырехугольник на 4 равновеликие части
x1=1:y1=1:x2=3:y2=11:x3=13:y3=15:x4=15:y4=5
x5=(x1+x2)/2:y5=(y1+y2)/2
x6=(x2+x3)/2:y6=(y2+y3)/2
x7=(x3+x4)/2:y7=(y3+y4)/2
x8=(x1+x4)/2:y8=(y1+y4)/2
x9=(x5+x7)/2:y9=(y5+y7)/2
xo=x5:yo=y5:xt=x7:yt=y7
p()
a57=a:b57=b
b68=y9+x9/a57
z=0.1:smax=10^10
a570=a57:b570=b57:b680=b68
for i=1 to 300000
a57=a570*(1+z*(ran()-.5))
b57=b570*(1+z*(ran()-.5))
b68=b680*(1+z*(ran()-.5))
a68=-1/a57
xo=x1:yo=y1:xt=x2:yt=y2
p()
a12=a:b12=b
aa1=a57:bb1=b57
aa2=a12:bb2=b12
x()
x5=x0:y5=y0
xo=x3:yo=y3:xt=x4:yt=y4
p()
a34=a:b34=b
aa1=a57:bb1=b57
aa2=a34:bb2=b34
x()
x7=x0:y7=y0
aa1=a57:bb1=b57
aa2=a68:bb2=b68
x()
x9=x0:y9=y0
xo=x2:yo=y2:xt=x3:yt=y3
p()
a23=a:b23=b
aa1=a23:bb1=b23:aa2=a68:bb2=b68
x()
x6=x0:y6=y0
xo=x1:yo=y1:xt=x4:yt=y4
p()
a14=a:b14=b
aa1=a14:bb1=b14:aa2=a68:bb2=b68
x()
x8=x0:y8=y0
xa=x1:ya=y1:xb=x8:yb=y8:xc=x9:yc=y9:xd=x5:yd=y5
s()
s1=s
xa=x5:ya=y5:xb=x9:yb=y9:xc=x6:yc=y6:xd=x2:yd=y2
s()
s2=s
xa=x9:ya=y9:xb=x7:yb=y7:xc=x3:yc=y3:xd=x6:yd=y6
s()
s3=s
xa=x8:ya=y8:xb=x4:yb=y4:xc=x7:yc=y7:xd=x9:yd=y9
s()
s4=s
d1=(s1-s2)^2+(s1-s3)^2+(s1-s4)^2
if d1<smax then
print x5,y5,x6,y6,x7,y7,x8,y8,x9,y9,d1
smax=d1
a570=a57:b570=b57:b680=b68
s10=s1:s20=s2:s30=s3:s40=s4
if d1<0.02 then z=0.00001:fi
fi
next i
print s10,s20,s30,s40
sub p()
a=(yt-yo)/(xt-xo):b=yo-a*xo
end sub
sub x()
x0=(bb1-bb2)/(aa2-aa1):y0=aa1*x0+bb1
end sub
sub s()
sP=xa*yb+xb*yc+xc*yd+xd*ya
sM=ya*xb+yb*xc+yc*xd+yd*xa
s=1/2*(sP-sM)
end sub

4 сентября 2021 г.