|||
import math
import random
import sys
import numpy as np
if (len(sys.argv)<3):
print 'too few input parameters, format:'
print 'python distparabola.py p x1 y1 '
print 'example:'
print 'python distparabola.py 1 1 0'
#sys.exit()
p = 1
x1 = 1
y1 = 0
print 'default values used.'
print ' '
else:
p = float(sys.argv[1])
x1 = float(sys.argv[2])
y1 = float(sys.argv[3])
TwoPi=math.pi*2
def sign(x):
if abs(x)<1e-6:
return 0
elif x>0:
return 1
else:
return -1
# Function from https://www.douban.com/note/523859088/
def getCubicRoot(a,b,c,d):
a=float(a)
b=b/(3.0*a)
c=c/(6.0*a)
d=d/(2.0*a)
solve=[None,None,None]
Alph=-b**3+3.0*b*c-d
Beta=b**2-2.0*c
Delt=Alph**2-Beta**3
if (Delt)>0.0:
tht=Alph+math.sqrt(Delt)
R1=abs(abs(tht)**(1.0/3.0))*sign(tht)
tht=Alph-math.sqrt(Delt)
R2=abs(abs(tht)**(1.0/3.0))*sign(tht)
solve[0]=-b+R1+R2
elif abs(Delt)<1e-6:
R1=abs(abs(Alph)**(1.0/3.0))*sign(Alph)
if abs(R1)<1e-6:
solve[0]=-b
else:
solve[0]=-b+2.0*R1
solve[1]=-b-R1
elif Delt<0.0:
tht=math.acos(Alph/(math.sqrt(Beta)*Beta))
solve[0]=-b+2.0*math.sqrt(Beta)*math.cos(tht/3.0)
solve[1]=-b+2.0*math.sqrt(Beta)*math.cos((tht+TwoPi)/3.0)
solve[2]=-b+2.0*math.sqrt(Beta)*math.cos((tht-TwoPi)/3.0)
return solve
c=-2.0*p*(y1-p)
d=-2.0*p**2*x1
a=1
b=0
ans=getCubicRoot(a,b,c,d)
#print ans
x0=ans[0]
y0=2.0*p*x0**2
dist=np.sqrt((x1-x0)**2+(y1-y0)**2)
print dist
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-9-26 13:35
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社