2444.sagews

# 定义给定高度height和半径radius的圆缺的面积
def HemiCircle(radius,height):
    var('r y h')
    assume(r>0)
    assume(0<h<=2*r)
    f=2*( r^2 - (y-r)^2 )^(1/2)
    area=f.integrate(y,0,h)
    if height>2*radius:
        height=2*radius
    if height<0:
        height=0
    return area(r=radius,h=height)

# 定义给定高度height和长半径longradius、短半径shortradius的椭圆面积
def HemiEllipse(longradius,shortradius,height):
    var('a b y h')
    assume(a>0)
    assume(b>0)
    assume(0<h<2*b)
    f=2*( ( 1 - (y-b)^2/(b^2) )*(a^2) )^(1/2)
    area=f.integrate(y,0,h)
    if height>2*shortradius:
        height=2*shortradius
    if height<0:
        height=0
    return area(a=longradius,b=shortradius,h=height)

# 定义给定油罐侧切面椭圆长半径longradius、短半径shortradius、显示油高height、油浮左侧长度leftlength、油浮右侧长度rightlength、和倾斜角度alpha时的体积
def TankVolume(longradius,shortradius,leftlength,rightlength,height,alpha):
    var('x h')
    assume(x>0)
    assume(h>0)
    f=HemiEllipse(longradius,shortradius,height-x*tan(alpha))
    if height<=(leftlength+rightlength)*tan(alpha):
        area=f.integrate(x,0,height/tan(alpha))
    else:
        area=f.integrate(x,0,leftlength+rightlength)
    return area(h=height+leftlength*tan(alpha))

# 定义第一问的体积V和高度h、旋转角度alpha的函数关系
var('h alpha')
assume(tan(alpha/180*pi)>0)
assume(h<0.6)#V和V0的表达式中惠出现arcsin，SAGE要求明确，
V=TankVolume(1.78/2,1.2/2,0.4,2.05,h,alpha/180*pi)
V0=TankVolume(1.78/2,1.2/2,0.4,2.05,h,0)#零度时单独定义，原因为积分式子中包含tan的倒数

# 水平入油
hh=[0.15902,0.17614,0.19259,0.20850,0.22393,0.23897,0.25366,0.26804,0.28216,0.29603,0.30969,0.32315,0.33644,0.34957,0.36256,0.37542,0.38816,0.40079,0.41332,0.42576,0.43812,0.45040,0.46262,0.47478,0.48689,0.49895,0.51097,0.52295,0.53490,0.54682,0.55872,0.57061,0.58248,0.59435,0.60622,0.61809,0.62996,0.64185,0.65375,0.66567,0.67763,0.67854,0.69053,0.69082,0.70285,0.71491,0.72703,0.73919,0.75142,0.76370,0.76416,0.77653,0.78899,0.80154,0.81419,0.82695,0.83983,0.85284,0.86600,0.87932,0.89282,0.89284,0.90653,0.92045,0.93461,0.94905,0.96380,0.97891,0.99443,1.01043,1.02699,1.04425,1.06237,1.08159,1.10233,1.12532,1.15236,1.19349]
vv=[0.312,0.362,0.412,0.462,0.512,0.562,0.612,0.662,0.712,0.762,0.812,0.862,0.912,0.962,1.012,1.062,1.112,1.162,1.212,1.262,1.312,1.362,1.412,1.462,1.512,1.562,1.612,1.662,1.712,1.762,1.812,1.862,1.912,1.962,2.012,2.062,2.112,2.162,2.212,2.262,2.312,2.31583,2.36583,2.36706,2.41706,2.46706,2.51706,2.56706,2.61706,2.66698,2.66883,2.71883,2.76883,2.81883,2.86883,2.91883,2.96883,3.01883,3.06883,3.11883,3.16883,3.16891,3.21891,3.26891,3.31891,3.36891,3.41891,3.46891,3.51891,3.56891,3.61891,3.66891,3.71891,3.76891,3.81891,3.86891,3.91891,3.96891]
p=list_plot(zip(hh,vv))
p+=plot(V0(h=x),(x,0,1.2),color='red')
p.show(frame=true)

#绝对误差
verror=[(V0(h=hh[each])-vv[each]) for each in range(len(hh))]
p=list_plot(zip(hh,verror))
p.show(frame=true)