Page MenuHome
Files F7872

implicit_implicit.py
Public
Actions

implicit_implicit.py
View Options

import math
from math import pi,sin,cos,acos,atan,tan,fabs,sqrt,log,tanh
from cmath import atanh
import implicit_vefm
from implicit_vefm import *
''' All classes of implicit have an init function to set them up and a calculate function to sample the field at a given point'''
class cauchypoint:
def __init__(self,centre,strength):
self.centre=centre
self.strength=strength
self.neg=0
if self.strength<0:
self.neg=1
def calculate(self,point):
length=point-self.centre
length.findlength()
l2=length.length*length.length
s2=self.strength*self.strength
calc=(1.0+s2*l2)**2
answer=1.0/calc
if self.neg:
answer=-answer
# return (answer*self.strength)-self.level#-.5#01
return answer
class distanceline:
def __init__(self,start,end,strength):
self.base=start
self.end=end
self.axis=end-start
self.axis.normalize()
self.axis.findlength()
line=end-start
line.findlength()
self.length=line.length
self.strength=strength
self.l2=self.length*self.length
self.s2=self.strength*self.strength
if self.strength<0:
self.neg=1
else:
self.neg=0
def calculate(self,point):
test=self.findclosepoint(point)
newcauch=cauchypoint(test,self.strength)
answer=newcauch.calculate(point) #one+two+three
return answer
def findclosepoint(self,point):
u=(point.x-self.base.x)*(self.end.x-self.base.x)+(point.y-self.base.y)*(self.end.y-self.base.y)+(point.z-self.base.z)*(self.end.z-self.base.z)
u=u/(self.length**2)
if u<0.0:
point=self.base
elif u>1.0:
point=self.end
else:
newx=self.base.x+(u*(self.end.x-self.base.x))
newy=self.base.y+(u*(self.end.y-self.base.y))
newz=self.base.z+(u*(self.end.z-self.base.z))
point=vertex(newx,newy,newz)
return point
class cauchyline:
def __init__(self,start,end,strength):
self.base=start
self.end=end
self.axis=end-start
self.axis.normalize()
self.axis.findlength()
line=end-start
line.findlength()
self.length=line.length
self.strength=strength
self.l2=self.length*self.length
self.s2=self.strength*self.strength
if self.strength<0:
self.neg=1
else:
self.neg=0
def calculate(self,r):
d=r-self.base
d.findlength()
x=self.dotproduct(d,self.axis)
lmx=self.length-x
d2=d.length*d.length
x2=x**2
p,p2,p3=self.generatep(d2,x2)
q2=self.generateq(d2,x)
a=x/(2*p2*(p2+(self.s2*x2)))
b=lmx/(2*(p2*q2))
term=atan((self.strength/p)*x)+atan((self.strength/p)*lmx)
c=1/(2*self.strength*p3)*term
answer=a+b+c
if self.neg:
answer=-answer
# return (answer*self.strength)-self.level#-.5#01
return answer
def dotproduct(self,a,b):
dot=(a.x*b.x)+(a.y*b.y)+(a.z*b.z)
return dot
def generatep(self,d2,x2):
p2=1+(self.s2*(d2-x2))
p=sqrt(p2)
p3=p**3
return p,p2,p3
def generateq(self,d2,x):
q2=1+(self.s2*(d2+self.l2-(2*self.length*x)))
return q2
class cauchylinevariable:
def __init__(self,start,end,strength,controlpoints):
self.base=start
self.end=end
self.axis=end-start
self.axis.normalize()
self.axis.findlength()
line=end-start
line.findlength()
self.length=line.length
self.strength=strength
self.l2=self.length*self.length
self.s2=self.strength*self.strength
if self.strength<0:
self.neg=1
else:
self.neg=0
self.q0=controlpoints[0]
self.q1=controlpoints[1]
self.q2=controlpoints[2]
self.q3=controlpoints[3]
def calculate(self,r):
d=edge(self.base,r)
d.findlength()
d2=d.vect.x*d.vect.x+d.vect.y*d.vect.y+d.vect.z*d.vect.z
h=d.vect.x*self.axis.x+d.vect.y*self.axis.y+d.vect.z*self.axis.z
s2=self.strength*self.strength
h2=h*h
lmh=self.length-h
pp= 1.0+s2*(d2-h2)
ww=pp+s2*h2
qq=pp+s2*lmh*lmh
p=sqrt(fabs(pp))
sdp=self.strength/p
tpp=2.0*pp
sp=self.strength*p
term=atan(sdp*h)+atan(sdp*lmh)
fline1=(h/ww+lmh/qq)/tpp+term/(tpp*sp)
ts2=2.0*s2
flinet=h*fline1+(1.0/ww-1.0/qq)/ts2
flinet2=2.0*h*flinet-(h2+pp/s2)*fline1+term/(s2*sp)
ts4=ts2*s2
flinet3=3.0*h*flinet2-3.0*h2*flinet+h2*h*fline1+log(qq/ww)/ts4-pp*(1.0/ww-1.0/qq)/ts4
p0=self.q0
p1=(-3.0*self.q0+3.0*self.q1)/self.length
p2=(3.0*self.q0-6.0*self.q1+3.0*self.q2)/(self.length*self.length)
p3=(-self.q0+3.0*self.q1-3.0*self.q2+self.q3)/(self.length*self.length*self.length)
result=p0*fline1+p1*flinet+p2*flinet2+p3*flinet3
return result#-0.01
class cauchyarc:
## TOTALLY BROKEN, DON'T WORK, NUFFIN
def __init__(self,origin,angle,radius,strength):
self.origin=origin
self.angle=angle
self.radius=radius
self.strength=strength
self.neg=0
self.r2=self.radius*self.radius
self.r4=self.radius*self.radius*self.radius*self.radius
self.s2=self.strength*self.strength
self.s4=self.strength*self.strength*self.strength*self.strength
self.k=2.0*self.radius*self.s2
# self.k=(2.0*self.radius*self.strength)**2
self.s2r2=self.r2*self.s2
self.sinan=sin(angle)
self.cosan=cos(angle)
self.tanan=tan(angle/2.0)
if self.strength<0:
self.neg=1
def calculate(self,point):
k=self.k
r2=self.r2
s2=self.s2
r4=self.r4
s4=self.s4
s2r2=self.s2r2
z2=point.z*point.z
y2=point.y*point.y
x2=point.x*point.x
sinan=self.sinan
cosan=self.cosan
tanan=self.tanan
d2=(point.x*point.x)+(point.y*point.y)+(point.z*point.z)
s2d2=s2*d2
b=1.0+s2r2+s2d2
p2=(-r4*s4)+(2.0*r2*s2)*(s2*(d2-2.0*point.z*point.z)-1.0)-((1.0+s2d2)**2)
w=(1.0+s2d2)**2
x=s2*(d2-(2.0*z2))#-1.0
y=2.0*s2r2
z=-(r4*s4)
# p2=z+(y*x-w)
# print "p2",p2,"w ",w,"x ",x,"z ",z
xp2=point.x*p2
p=sqrt(abs(p2))
# p=sqrt(p2)
p3=p2*p2
A=(b*point.y)/(xp2*(k*point.x-b))
m=k*(x2*y2)*sinan-(b*point.y)
n=xp2*(k*(point.x*cosan+point.y*sinan)-b)
B=m/n
i=atanh((k*point.y)/p)
j=((k*point.x+b)*tanan)-(k*point.y)
o=atanh(j/p)
C=(2*b/p3)*(i+j)
return abs(A+B+C)
class cauchytriangle:
pass
def __init__(self,triangle,strength,radius):
self.strength=strength
self.tri=triangle
self.radius=1.0/radius
self.e0=edge(self.tri.vertices[0],self.tri.vertices[1])
self.e1=edge(self.tri.vertices[1],self.tri.vertices[2])
self.e2=edge(self.tri.vertices[2],self.tri.vertices[0])
self.e0.findlength()
self.e1.findlength()
self.e2.findlength()
# print "lengths",self.e0.vect.length,self.e1.vect.length,self.e2.vect.length
self.longest,self.apex,self.base,self.hypot=self.findconfig()
self.b,self.a1,self.a2,self.h,self.U,self.V=self.findmeasures()
# print "self.longest.vect.length",self.longest.vect.length
# print "self.apex",self.apex
# print "self.base",self.base
# print "self.hypot.vect.length",self.hypot.vect.length
# print "self.b",self.b.x,self.b.y,self.b.z
# print "self.a1",self.a1
# print "self.a2",self.a2
# print "self.a1+self.a2",self.a1+self.a2
# print "self.h",self.h
# print "self.U",self.U.length
# print "self.V",self.V.length
#
#rightangle=1.5707963267948966
def findconfig(self):
if self.e0.vect.length>self.e1.vect.length:
longest=self.e0
else:
longest=self.e1
# print "first length test",longest.vect.length
if longest.vect.length<self.e2.vect.length:
longest=self.e2
# print "second length test",longest.vect.length
if longest is self.e0:
apex=self.tri.vertices[2]
base=self.tri.vertices[0]
hypot=self.e2
# print "longest is e0"
if longest is self.e1:
apex=self.tri.vertices[0]
base=self.tri.vertices[1]
hypot=self.e0
# print "longest is e1"
if longest is self.e2:
apex=self.tri.vertices[1]
base=self.tri.vertices[2]
hypot=self.e1
# print "longest is e2"
return longest,apex,base,hypot
def findmeasures(self):
# baseangle=self.tri.dotproduct(self.longest.vect,self.hypot.vectb)
dot=self.dotproduct(self.longest.vect,self.hypot.vectb)
costheta=dot/(self.longest.vect.length*self.hypot.vect.length)
baseangle=acos(costheta)
# baseangle=acos(self.dotproduct(self.longest.vect,self.hypot.vectb))
# print "baseangle ",baseangle
apexangle=pi*0.5-baseangle
# print "apexangle",apexangle
h=sin(baseangle)*self.hypot.vect.length
# print "h",h
a1=sin(apexangle)*self.hypot.vect.length
# print "a1",a1
a2=self.longest.vect.length-a1
# print "a2",a2
# print "longest",self.longest,self.longest.vect
bnorm=vertex(self.longest.vect.x,self.longest.vect.y,self.longest.vect.z)
bnorm.normalize()
b=self.base+(bnorm*a1)
# print "b",b.x,b.y,b.z
U=self.base-b
V=b-self.apex
U.normalize()
V.normalize()
return b,a1,a2,h,U,V
def calculate(self,point):
point=point*self.radius
d=self.b-point
d.findlength()
u=self.dotproduct(d,self.U)
v=self.dotproduct(d,self.V)
s2=self.strength*self.strength
d2=d.length*d.length
u2=u*u
v2=v*v
h2=self.h*self.h
a12=self.a1*self.a1
a22=self.a2*self.a2
C2=(1.0/s2)+(d2-u2)
g=v-self.h
q=C2-v2
w=C2-(2.0*v*self.h)+h2
m=(self.a2*g)+(u*self.h)
n=(u*self.h)-(self.a1*g)
A2=(a12*w)+(h2*(q+u2))-2.0*self.a1*self.h*u*g
B2=(a22*w)+(h2*(q+u2))+2.0*self.a2*self.h*u*g
A=sqrt(A2)
B=sqrt(B2)
C=sqrt(C2)
atanA1=atan(((v*self.h)+(self.a1*(self.a1+u)))/A)
atanA2=atan(((g*self.h)+(self.a1*u))/-A)
atanB1=atan(((v*self.h)+(self.a2*(self.a2-u)))/-B)
atanB2=atan(((g*self.h)-(self.a2*u))/B)
atanC1=atan((self.a1+u)/C)
atanC2=atan((self.a2-u)/C)
one=(n/A)*(atanA1+atanA2)
two=(m/B)*(atanB1+atanB2)
three=(v/C)*(atanC1+atanC2)
div=1.0/(2.0*q*self.strength)
whole=div*(one+two+three)
return whole
def dotproduct(self,a,b):
dot=(a.x*b.x)+(a.y*b.y)+(a.z*b.z)
return dot

File Metadata

Mime Type
text/x-python
Storage Engine
local-disk
Storage Format
Raw Data
Storage Handle
24/9a/1bb710118b72090bb6cc7d52ff1f

Event Timeline