CUMCM2024ProblemA-T3

Python代码

import numpy as np
import math
import sympy as sy
from sqlalchemy.sql.operators import truediv
from sympy import symbols,diff
import pandas as pd

def Proj1(di):
    return (2 * np.pi)*(16-450/di)
    C = []  # 各个时间龙头的角度
    S = []  # 各个时间龙头的路程（单位：厘米）
    c = 0.0
    pre_s=0.0
    while 1:
        cita = c
        R=di*(16 - c/(2 * np.pi))
        if R - 450 < 0:
            return c
        l=0.001
        r=1.0
        while l<=r:
            mid=(l+r)/2
            c=cita+mid
            now_s = di * ((np.log(((c - 32 * np.pi) ** 2 + 1) ** 0.5 + c - np.pi * 32) - np.log(((32 * np.pi) ** 2 + 1) ** 0.5 - 32 * np.pi)) + ((c - 32 * np.pi) * ((c - 32 * np.pi) ** 2 + 1) ** 0.5) + 32 * np.pi * (((32 * np.pi) ** 2 + 1) ** 0.5)) / (4 * np.pi)
            if now_s-pre_s>100:
                r=mid-0.00001
            else:
                l=mid+0.00001

        c=cita+l
        Final_s = di * ((np.log(((c - 32 * np.pi) ** 2 + 1) ** 0.5 + c - np.pi * 32) - np.log(
            ((32 * np.pi) ** 2 + 1) ** 0.5 - 32 * np.pi)) + (
                                  (c - 32 * np.pi) * ((c - 32 * np.pi) ** 2 + 1) ** 0.5) + 32 * np.pi * (
                                  ((32 * np.pi) ** 2 + 1) ** 0.5)) / (4 * np.pi)
        #print("NOWCITA ", c)
        C.append(c)
        S.append(Final_s)
        #print("NOWS ", Final_s)
        pre_s=Final_s


traversal = 3




def Proj2(cita,di):
    def d(x):
        nonlocal y
        Z = (di ** 2 * ((16 - x / 2 / np.pi) ** 2 + (16 - y / 2 / np.pi) ** 2 - 2 * (16 - x / 2 / np.pi) * (
                    16 - y / 2 / np.pi) * sy.cos(x - y)) - 165 ** 2) / 1000000
        return Z

    def d2(x):
        nonlocal y
        Z = (di ** 2 * ((16 - x / 2 / np.pi) ** 2 + (16 - y / 2 / np.pi) ** 2 - 2 * (16 - x / 2 / np.pi) * (
                    16 - y / 2 / np.pi) * sy.cos(x - y)) - 286 ** 2) / 1000000
        return Z

    a = symbols('a', real=True)
    y=cita
    y0 = y
    X = []
    for i in range(223):
        if i == 0:
            x1 = y - 0.1
            for t in range(traversal):
                x = x1
                Delt = d2(x) / diff(d2(a), a).subs({a: x})
                if abs(Delt) < np.pi/32:
                    x1=x-Delt
                else:
                    x1=x-np.pi/32
            y = x1
            if x1 < 0:
                break
            X.append(x1)
        else:
            x1 = y - 0.1
            for t in range(traversal):
                x = x1
                Delt = d(x) / diff(d(a), a).subs({a: x})
                if abs(Delt) < np.pi / 32:
                    x1 = x - Delt
                else:
                    x1 = x - np.pi / 32
            y = x1
            if x1 < 0:
                break
            X.append(x1)
    print(f'各位置角度{[y0] + X}')
    return [y0] + X




def position(c,di):
    x = di * math.cos(c) * (16 - c / (2 * math.pi))
    y = -di * math.sin(c) * (16 - c / (2 * math.pi))
    return [x, y]


def k_of_v(c):
    x, y = position(c)
    return -x / y


def k_of_stick(c1, c2):
    x1, y1 = position(c1)
    x2, y2 = position(c2)
    return (y1 - y2) / (x1 - x2)


def tan_c(k1, k2):
    return (k1 - k2) / (1 + k1 * k2)


def cos_c(tan_c):
    return (1 / (1 + tan_c ** 2)) ** 0.5


def C_to_V(same_second_c):
    ssc = same_second_c
    V = [1]
    for i in range(len(ssc) - 1):
        k1 = k_of_v(ssc[i])
        k2 = k_of_v(ssc[i + 1])
        k3 = k_of_stick(ssc[i], ssc[i + 1])
        cos1 = cos_c(tan_c(k1, k3))
        cos2 = cos_c(tan_c(k2, k3))
        v = V[i] * cos1 / cos2
        V.append(v)
    return V


def C_to_xy(same_second_c):
    ssc = same_second_c
    XY = []
    for c in ssc:
        XY = XY + position(c)
    return XY


class ExtendedAndParallelLines:
    def __init__(self, point1, point2, extension=27.5, distance=15):
        self.point1 = point1
        self.point2 = point2
        self.extension = extension
        self.distance = distance

        # 计算延长后的线段
        self.extended_point1, self.extended_point2 = self.extend_segment()

        # 计算延长线段的斜率和y轴截距
        self.slope, self.intercept = self.calculate_line_equation(self.extended_point1, self.extended_point2)

        # 计算垂直线段的斜率和y轴截距
        self.perpendicular_slope = -1 / self.slope
        self.intercept1 = self.calculate_new_intercept(self.extended_point1, self.distance)
        self.intercept2 = self.calculate_new_intercept(self.extended_point1, -self.distance)

        # 存储平移后两条线段的端点坐标

        self.moved_perpendicular_point1 ,self.moved_perpendicular_point2 = self.calculate_moved_point12()


        self.moved_perpendicular_point3, self.moved_perpendicular_point4 = self.calculate_moved_point34()

    def extend_segment(self):
        """延长线段"""
        x1, y1 = self.point1
        x2, y2 = self.point2
        length = math.sqrt((x2 - x1)**2 + (y2 - y1)**2)
        direction = ((x2 - x1) / length, (y2 - y1) / length)
        extended_point1 = (x1 - self.extension * direction[0], y1 - self.extension * direction[1])
        extended_point2 = (x2 + self.extension * direction[0], y2 + self.extension * direction[1])
        return extended_point1, extended_point2

    def calculate_line_equation(self, point1, point2):
        """计算通过两点直线的斜率和y轴截距"""
        x1, y1 = point1
        x2, y2 = point2
        slope = (y2 - y1) / (x2 - x1)
        intercept = y1 - slope * x1
        return slope, intercept

    def calculate_new_intercept(self, point, distance):
        """计算垂直线段的y轴截距"""
        x, y = point
        # 使用点到直线的距离公式计算新的截距
        intercept = y - self.perpendicular_slope * x + distance / math.sqrt(1 + self.perpendicular_slope**2)
        return intercept

    def calculate_moved_point12(self):
        """计算平移后的点12"""
        x1, y1 = self.extended_point1
        x2, y2 = self.extended_point2
        length = math.sqrt((x2 - x1) ** 2 + (y2 - y1) ** 2)
        direction = ((x2 - x1) / length, (y2 - y1) / length)
        moved_extended_point1 = (x1 - self.distance * direction[1], y1 + self.distance * direction[0])
        moved_extended_point2 = (x2 - self.distance * direction[1], y2 + self.distance * direction[0])
        return moved_extended_point1, moved_extended_point2

    def calculate_moved_point34(self):
        """计算平移后的点34"""
        x1, y1 = self.extended_point1
        x2, y2 = self.extended_point2
        length = math.sqrt((x2 - x1) ** 2 + (y2 - y1) ** 2)
        direction = ((x2 - x1) / length, (y2 - y1) / length)
        moved_extended_point1 = (x1 + self.distance * direction[1], y1 - self.distance * direction[0])
        moved_extended_point2 = (x2 + self.distance * direction[1], y2 - self.distance * direction[0])
        return moved_extended_point1, moved_extended_point2

    def __str__(self):
        return (f"原始线段: 通过点{self.point1}和{self.point2}\n"
                f"延长后的线段: 通过点{self.extended_point1}和{self.extended_point2}\n"
                f"线段: y = {self.slope}x + {self.intercept1}\n"
                f"平移后垂直线段的端点1_1: {self.moved_perpendicular_point1}\n"
                f"平移后垂直线段的端点1_2: {self.moved_perpendicular_point2}\n"
                f"平移后垂直线段的端点2_1: {self.moved_perpendicular_point3}\n"
                f"平移后垂直线段的端点2_2: {self.moved_perpendicular_point4}\n")


def Intersect(l1 , l2):
    v1 = (l1[0] - l2[0], l1[1] - l2[1])
    v2 = (l1[0] - l2[2], l1[1] - l2[3])
    v0 = (l1[0] - l1[2], l1[1] - l1[3])
    a = v0[0] * v1[1] - v0[1] * v1[0]
    b = v0[0] * v2[1] - v0[1] * v2[0]

    temp = l1
    l1 = l2
    l2 = temp
    v1 = (l1[0] - l2[0], l1[1] - l2[1])
    v2 = (l1[0] - l2[2], l1[1] - l2[3])
    v0 = (l1[0] - l1[2], l1[1] - l1[3])
    c = v0[0] * v1[1] - v0[1] * v1[0]
    d = v0[0] * v2[1] - v0[1] * v2[0]

    if a * b < 0 and c * d < 0:
        return True
    else:
        return False


def C_to_xy2(same_second_c,di):
    ssc = same_second_c
    XY = []
    for c in ssc:
        XY.append(position(c,di))
    return XY

def ifcross(Line1, Line2):
    Flag1=Intersect((Line1.moved_perpendicular_point1[0],Line1.moved_perpendicular_point1[1],
                     Line1.moved_perpendicular_point2[0],Line1.moved_perpendicular_point2[1]),
                    (Line2.moved_perpendicular_point1[0],Line2.moved_perpendicular_point1[1],
                     Line2.moved_perpendicular_point2[0],Line2.moved_perpendicular_point2[1]))
    Flag2=Intersect((Line1.moved_perpendicular_point1[0],Line1.moved_perpendicular_point1[1],
                     Line1.moved_perpendicular_point2[0],Line1.moved_perpendicular_point2[1]),
                    (Line2.moved_perpendicular_point3[0],Line2.moved_perpendicular_point3[1],
                     Line2.moved_perpendicular_point4[0],Line2.moved_perpendicular_point4[1]))
    Flag3=Intersect((Line1.moved_perpendicular_point3[0],Line1.moved_perpendicular_point3[1],
                     Line1.moved_perpendicular_point4[0],Line1.moved_perpendicular_point4[1]),
                    (Line2.moved_perpendicular_point1[0],Line2.moved_perpendicular_point1[1],
                     Line2.moved_perpendicular_point2[0],Line2.moved_perpendicular_point2[1]))
    Flag4=Intersect((Line1.moved_perpendicular_point3[0],Line1.moved_perpendicular_point3[1],
                     Line1.moved_perpendicular_point4[0],Line1.moved_perpendicular_point4[1]),
                    (Line2.moved_perpendicular_point3[0],Line2.moved_perpendicular_point3[1],
                     Line2.moved_perpendicular_point4[0],Line2.moved_perpendicular_point4[1]))
    return Flag1+Flag2+Flag3+Flag4

def judge(Lines):
    print("NUMlines", len(Lines))
    for Line1n in range(len(Lines)):
        for Line2n in range(Line1n+2,len(Lines)):
            if ifcross(Lines[Line1n],Lines[Line2n])>0:
                print("cross:", Lines[Line1n],Lines[Line2n])
                print("cross:", Line1n, Line2n)
                return True
    return False

l=42.3
r=43.3#42.8
while l<=r:
    #l=l+1
    #mid=l
    mid=(l+r)/2
    ansflag = 0
    for i in range(0,60):
        cita=Proj1(mid)-float(i)/20
        all2_xy = C_to_xy2(Proj2(cita,mid),mid)
        All_lines=[]
        print(all2_xy)
        flag=0
        pre=[0,0]
        for xy in all2_xy:
            if flag==0:
                pre=xy
                flag=1
            else:
                All_lines.append(ExtendedAndParallelLines(pre,xy))
                pre=xy
        judgecross=judge(All_lines)
        print("mid=",mid,"相交：",judgecross)
        if judgecross==1:
            #l=mid+0.01
            break
        else:
            #r=mid-0.01
            print("ANS =",mid)
            ansflag=ansflag+1
            #print("ANSflag =",ansflag)
            pass
    if ansflag==60:
        r=mid-0.01
    else:
        l=mid+0.01
#print("最终结果=",l)