' pntfillet : 2d fillet of 2 lines through a point

function ptw(a)
ptw=" w"&cstr(a(0))&","&cstr(a(1))&",0 "
end function
 
function v2prodi(a,b)
v2prodi=a(0)*b(0)+a(1)*b(1)
end function

function v2prode(a,b)
v2prode=a(0)*b(1)-a(1)*b(0)
end function

function v2ver(v)
m=sqr(v(0)^2+v(1)^2)
if m=0 then 
  v2ver=null
else
  v2ver=array(v(0)/m,v(1)/m)
end if
end function

function v2sub(a,b)
v2sub=array(a(0)-b(0),a(1)-b(1))
end function

function v2sum(a,b)
v2sum=array(a(0)+b(0),a(1)+b(1))
end function

function v2mul(v,a)
v2mul=array(v(0)*a,v(1)*a)
end function

function v2psuml(p,u,l)
v2psuml=array(p(0)+u(0)*l,p(1)+u(1)*l)
end function

function v2verpp(a,b)
v=v2sub(b,a)
v2verpp=v2ver(v)
end function

function v2mid(pa,pb)
v2mid=array((pa(0)+pb(0))/2,(pa(1)+pb(1))/2)
end function

function v2dist(a,b)
v2dist=sqr((a(0)-b(0))^2+(a(1)-b(1))^2)
end function

function prjptlin(pt,pl,ul)
l=v2prodi(v2sub(pt,pl),ul)
prjptlin=v2psuml(pl,ul,l)
end function

function intlinlin(pa,ua,pb,ub)
if abs(ua(0))>abs(ua(1)) then
  aa=(pb(0)-pa(0))/ua(0)
  bb=ub(0)/ua(0)
  h=(pb(1)-pa(1)-aa*ua(1))/(bb*ua(1)-ub(1))
  intlinlin=v2psuml(pb,ub,h)
else
  cc=(pa(0)-pb(0))/ub(0)
  dd=ua(0)/ub(0)
  k=(pb(1)-pa(1)+cc*ub(1))/(ua(1)-dd*ub(1))
  intlinlin=v2psuml(pa,ua,k)
end if
end function

sub main
la=rhino.getobject("First line ?")
lb=rhino.getobject("Second line ?")
pt=rhino.getpoint("Point ?")
rad=rhino.getstring("Radius Min or Max ?","Max",array("Min","Max"))
pa0=rhino.curvestartpoint(la)
pa1=rhino.curveendpoint(la)
pb0=rhino.curvestartpoint(lb)
pb1=rhino.curveendpoint(lb)
ua=v2verpp(pa0,pa1)
ub=v2verpp(pb0,pb1)
pint=intlinlin(pa0,ua,pb0,ub)
' rhino.addpoint(pint)
if v2dist(pint,pa0) > v2dist(pint,pa1) then
  tm=pa0
  pa0=pa1
  pa1=tm
  ua=v2verpp(pa0,pa1)
end if
if v2dist(pint,pb0) > v2dist(pint,pb1) then
  tm=pb0
  pb0=pb1
  pb1=tm
  ub=v2verpp(pb0,pb1)
end if
pa2=v2psuml(pint,ua,10)
pb2=v2psuml(pint,ub,10)
pm=v2mid(pa2,pb2)
' rhino.addpoint(pm)
um=v2verpp(pint,pm)
cosal=v2prodi(ua,um)
sinal=sqr(1-cosal*cosal)
kb=v2sub(pint,pt)
aa=um(0)^2+um(1)^2-sinal^2
bb=kb(0)*um(0)+kb(1)*um(1)
cc=kb(0)^2+kb(1)^2
k0=(-bb+sqr(bb*bb-aa*cc))/aa
k1=(-bb-sqr(bb*bb-aa*cc))/aa
c0=v2psuml(pint,um,k0)
c1=v2psuml(pint,um,k1)
' rhino.addpoint(c0)
' rhino.addpoint(c1)
r0=v2dist(c0,pt)
r1=v2dist(c1,pt)
if r0>r1 then
  tm=r0
  r0=r1
  r1=tm
  tm=c0
  c0=c1
  c1=tm
end if
if rad="Min" then
  r=r0
  c=c0
else
  r=r1
  c=c1
end if
pra=prjptlin(c,pint,ua)
prb=prjptlin(c,pint,ub)
' rhino.addpoint(pra)
' rhino.addpoint(prb)
rhino.command"arc startpoint"&ptw(pra)&ptw(prb)&ptw(pt)
' rhino.command"fillet "&cstr(r)&" selid "&la&" selid "&lb
end sub

pi=3.14159265358979323846

main