{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "769e731f",
   "metadata": {
    "hide_input": false
   },
   "outputs": [],
   "source": [
    "from numpy import *\n",
    "import matplotlib.pyplot as plt\n",
    "import time\n",
    "from scipy.spatial import cKDTree\n",
    "import pandas as pd\n",
    "\n",
    "\n",
    "def arc(radius=0,start_angle=0,end_angle=0,cp=[0,0],s=20):\n",
    "    cp=array(cp)\n",
    "    r=linspace(start_angle,end_angle,s+1)\n",
    "    x=radius*cos(pi/180*r)\n",
    "    y=radius*sin(pi/180*r)\n",
    "    c=(cp+array([x,y]).swapaxes(0,1))\n",
    "    return c.tolist()        \n",
    "\n",
    "def pts(p):\n",
    "    return array(p)[:,0:2].cumsum(axis=0).tolist()\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "def pts1(p):\n",
    "    p=[[a[0],a[1],0] if len(a)==2 else a for a in p]\n",
    "    b=array(p)[:,0:2].cumsum(axis=0)\n",
    "    c=array([array(p)[:,2].tolist()])\n",
    "    return concatenate((b,c.T),1).tolist()\n",
    "\n",
    "def cw(p):\n",
    "    p=array(p)[:,0:2]\n",
    "    q=p[1:].tolist()+[p[0].tolist()]\n",
    "    r=[p[len(p)-1].tolist()]+p[0:len(p)-1].tolist()\n",
    "    a=array(q)-p\n",
    "    b=p-array(r)\n",
    "    c=where(cross(b,a)>0,1,0).sum()\n",
    "    d=where(cross(b,a)<0,1,0).sum()\n",
    "    e=1 if c<d else -1\n",
    "    return e    \n",
    "\n",
    "def cwv(p):\n",
    "    p=sec\n",
    "    p0=[p[len(p)-1]]+p[:-1]\n",
    "    p1=p\n",
    "    p2=p[1:]+[p[0]]\n",
    "    p0,p1,p2=array([p0,p1,p2])\n",
    "    p=array([p0,p1,p2]).transpose(1,0,2).tolist()\n",
    "    return [cw(p1) for p1 in p]\n",
    "\n",
    "\n",
    "def ang(x,y):\n",
    "    if x>=0 and y>=0:\n",
    "        return arctan(y/(0.000001 if x==0 else x))*180/pi\n",
    "    elif x<0 and y>=0:\n",
    "        return 180-abs(arctan(y/x))*180/pi\n",
    "    elif  x<0 and y<0:\n",
    "        return 180+abs(arctan(y/x))*180/pi\n",
    "    else:\n",
    "        return 360-abs(arctan(y/(0.000001 if x==0 else x)))*180/pi\n",
    "\n",
    "\n",
    "\n",
    "def q(vector=[1,0,0],point=[0,5,0],theta=0):\n",
    "\n",
    "    t=theta\n",
    "    v=vector/(linalg.norm(vector)+.00001)\n",
    "    a=t/2*pi/180\n",
    "    p=[cos(a),multiply(v,sin(a))]\n",
    "    p1=[p[0],-p[1]]\n",
    "    q=[0,[point[0],point[1],0] if len(point)==2 else point]\n",
    "    pq=[p[0]*q[0]-p[1]@q[1],multiply(p[0],q[1])+p[1]*q[0]+cross(p[1],q[1])]\n",
    "    pqp1=[pq[0]*p1[0]-pq[1]@p1[1],pq[0]*p1[1]+pq[1]*p1[0]+cross(pq[1],p1[1])]\n",
    "    transformation=pqp1[1].tolist()\n",
    "    return transformation\n",
    "\n",
    "def uv(v):\n",
    "    v=array(v)\n",
    "    return (v/(linalg.norm(v))).tolist()\n",
    "\n",
    "def norm(v):\n",
    "    return linalg.norm(v)\n",
    "\n",
    "def fillet2d(pl,rl,s):\n",
    "    p0=array(array(pl)[len(pl)-2:len(pl)].tolist()+array(pl)[0:len(pl)-2].tolist())\n",
    "    p1=array([array(pl)[len(pl)-1].tolist()]+array(pl)[0:len(pl)-1].tolist())\n",
    "    p2=array(pl)\n",
    "    p3=array(array(pl)[1:len(pl)].tolist()+[array(pl)[0].tolist()])\n",
    "    p4=array(array(pl)[2:len(pl)].tolist()+array(pl)[0:2].tolist())\n",
    "    r0=array([array(rl)[len(rl)-1].tolist()]+array(rl)[0:len(rl)-1].tolist())\n",
    "    r1=array(rl)\n",
    "    r2=array(array(rl)[1:len(rl)].tolist()+[array(rl)[0].tolist()])\n",
    "    u0=(p0-p1)/linalg.norm(p0-p1,axis=1).reshape(-1,1)\n",
    "    u1=(p2-p1)/linalg.norm(p2-p1,axis=1).reshape(-1,1)\n",
    "    u2=(p1-p2)/linalg.norm(p1-p2,axis=1).reshape(-1,1)\n",
    "    u3=(p3-p2)/linalg.norm(p3-p2,axis=1).reshape(-1,1)\n",
    "    u4=(p2-p3)/linalg.norm(p2-p3,axis=1).reshape(-1,1)\n",
    "    u5=(p4-p3)/linalg.norm(p4-p3,axis=1).reshape(-1,1)\n",
    "    theta0= (180-arccos(einsum('ij,ij->i',u0,u1))*180/pi)/2\n",
    "    theta1= (180-arccos(einsum('ij,ij->i',u2,u3))*180/pi)/2\n",
    "    theta2= (180-arccos(einsum('ij,ij->i',u4,u5))*180/pi)/2\n",
    "    return f2d(p1,p2,p3,r0,r1,r2,theta0,theta1,theta2,u2,u3,s)\n",
    "    \n",
    "\n",
    "\n",
    "def each(a):\n",
    "    c=[]\n",
    "    for p in a:\n",
    "        for p1 in p:\n",
    "            c.append(p1)\n",
    "    return c\n",
    "\n",
    "def cr1(pl,s=20):\n",
    "    pl1=array(pl)[:,0:2].tolist()\n",
    "    rl=[0 if len(p)==2 else p[2] for p in pl]\n",
    "    return fillet2d(pl1,rl,s)\n",
    "\n",
    "def cr(pl,s=20):\n",
    "    sec=array(cr1(pl,s)).round(8)\n",
    "    s1=sec[sort(unique(sec,axis=0,return_index=True)[1])].tolist()\n",
    "    return s1\n",
    "\n",
    "def cr_c(pl,s=20):\n",
    "    sec=array(cr1(pl,s)).round(8)\n",
    "    s1=sec[sort(unique(sec,axis=0,return_index=True)[1])].tolist()\n",
    "    p0,p1=array([s1[len(s1)-1],s1[0]])\n",
    "    v=p1-p0\n",
    "    p=(p0+v*.999).tolist()\n",
    "    \n",
    "    return s1+[p]\n",
    "\n",
    "def f2d(p1,p2,p3,r0,r1,r2,theta0,theta1,theta2,u2,u3,s):\n",
    "    l1=linalg.norm(p1-p2,axis=1)\n",
    "    l2=r0*tan(theta0*pi/180)+r1*tan(theta1*pi/180)\n",
    "    l3=linalg.norm(p3-p2,axis=1)\n",
    "    l4=r1*tan(theta1*pi/180)+r2*tan(theta2*pi/180)\n",
    "    rf1=[r1[i] if l1[i]>l2[i] else 0 if l2[i]==0 else l1[i]/l2[i]*r1[i] for i in range(len(l1))]\n",
    "    rf2=[r1[i] if l3[i]>l4[i] else 0 if l4[i]==0 else l3[i]/l4[i]*r1[i] for i in range(len(l3))]\n",
    "    rf=swapaxes([rf1,rf2],0,1).min(axis=1)\n",
    "    p=p2+u2*(rf*tan(theta1*pi/180)).reshape(-1,1)\n",
    "    q=swapaxes([p1,p2,p3],0,1)\n",
    "    r=[]\n",
    "    for i in range(len(q)):\n",
    "        r.append(cw(q[i]))\n",
    "    r=array(r)\n",
    "    n=r==-1\n",
    "    n1=p-u2@array(rm(90))*rf.reshape(-1,1)\n",
    "    n2=p-u2@array(rm(-90))*rf.reshape(-1,1)\n",
    "    cp=[]\n",
    "    for i in range(len(n)):\n",
    "        if n[i]==True:\n",
    "            cp.append(n1[i])\n",
    "        else:\n",
    "            cp.append(n2[i])\n",
    "\n",
    "    cp=array(cp)\n",
    "    a1=[]\n",
    "#     alpha=(p-cp)/linalg.norm(p-cp,axis=1).reshape(-1,1)\n",
    "    alpha=[ [0,0] if linalg.norm(p[i]-cp[i])==0 else (p[i]-cp[i])/linalg.norm(p[i]-cp[i]) for i in range(len(p))]\n",
    "    for i in range(len(alpha)):\n",
    "        a1.append(ang(alpha[i][0],alpha[i][1]))\n",
    "    a1=array(a1)\n",
    "    boo=[]\n",
    "    for i in range(len(p1)):\n",
    "        boo.append(cw([p1[i],p2[i],p3[i]]))\n",
    "    boo=array(boo)   \n",
    "    a2=where(boo==-1,a1+2*theta1,a1-2*theta1)\n",
    "    ar=[]\n",
    "    for i in range(len(rf)):\n",
    "        ar.append(arc(rf[i],a1[i],a2[i],cp[i],s))\n",
    "    ar=array(ar)\n",
    "    c1=r1==0\n",
    "    c2=linalg.norm(u2-u3,axis=1)<.2\n",
    "    d=[]\n",
    "    for i in range(len(c1)):\n",
    "        if c1[i] or c2[i]:\n",
    "            d.append([p2[i].tolist()])\n",
    "        else:\n",
    "            d.append(ar[i].tolist())\n",
    "    return concatenate(d).tolist()\n",
    "\n",
    "\n",
    "def flip(sec): \n",
    "    return sec[::-1]\n",
    "    \n",
    "\n",
    "def r_3p(p1,p2,p3):\n",
    "    p4=add(p1,divide(subtract(p2,p1),2)).tolist()\n",
    "    p5=add(p2,divide(subtract(p3,p2),2)).tolist()\n",
    "    u1=uv(subtract(p2,p4))\n",
    "    u2=uv(subtract(p3,p5))\n",
    "    p6=add(p4,dot(u1,rm(90))).tolist()\n",
    "    p7=add(p5,dot(u2,rm(90))).tolist()\n",
    "    cp=i_p2d([p4,p6],[p5,p7])\n",
    "    r=norm(subtract(p1,cp))\n",
    "    return r\n",
    "\n",
    "\n",
    "def max_r(sec):\n",
    "    c=[]\n",
    "    for i in range(len(sec)):\n",
    "        i_2minus=len(sec)-2 if i==0 else len(sec)-1 if i==1 else i-2\n",
    "        i_minus=len(sec)-1 if i==0 else i-1\n",
    "        i_plus=i+1 if i<len(sec)-1 else 0\n",
    "        i_2plus=i+2 if i<len(sec)-2 else 0 if i<len(sec)-1 else 1\n",
    "        pi_2minus=sec[i_2minus]\n",
    "        pi_minus=sec[i_minus]\n",
    "        pi=sec[i]\n",
    "        pi_plus=sec[i_plus]\n",
    "        pi_2plus=sec[i_2plus]\n",
    "        v1=subtract(pi_minus,pi_2minus)\n",
    "        v2=subtract(pi,pi_minus)\n",
    "        v3=subtract(pi_plus,pi)\n",
    "        v4=subtract(pi_2plus,pi_plus)\n",
    "        l1=norm(v1).round(3)\n",
    "        l2=norm(v2).round(3)\n",
    "        l3=norm(v3).round(3)\n",
    "        l4=norm(v4).round(3)\n",
    "        r1=r_3p([pi_2minus,pi_minus,pi]).round(3)\n",
    "        r2=r_3p([pi_minus,pi,pi_plus]).round(3)\n",
    "        r3=r_3p([pi,pi_plus,pi_2plus]).round(3)\n",
    "        c.append(0 if l2!=l3 and (r1!=r2 or r2!=r3) else r2)\n",
    "    return max(c)\n",
    "        \n",
    "\n",
    "def offset_l(l,d):\n",
    "    u=uv(subtract(l[1],l[0]))\n",
    "    p0=add(l[0],dot(u,multiply(d,rm(-90)))).tolist()\n",
    "    p1=add(l[1],dot(u,multiply(d,rm(-90)))).tolist()\n",
    "    return [p0,p1]\n",
    "\n",
    "def seg(sec):\n",
    "    c=[]\n",
    "    for i in range(len(sec)):\n",
    "        i_plus=i+1 if i<len(sec)-1 else 0\n",
    "        p0=sec[i]\n",
    "        p1=sec[i_plus]\n",
    "        l=[p0,p1]\n",
    "        c.append(l)\n",
    "    return c\n",
    "\n",
    "def offset_seg(sec,r):\n",
    "    s=seg(sec)\n",
    "    c=[]\n",
    "    for p in s:\n",
    "        c.append(offset_l(p,r))\n",
    "    return c\n",
    "\n",
    "def offset_segv(sec,d):\n",
    "    s=sec\n",
    "    s1=s[1:]+[s[0]]\n",
    "    x=(array(s1)-array(s))\n",
    "    y=linalg.norm(x,axis=1)\n",
    "    u=x/y.reshape(-1,1)\n",
    "    p0=array(s)+u@array(rm(-90))*d\n",
    "    p1=array(s1)+u@array(rm(-90))*d\n",
    "    return swapaxes([p0,p1],0,1).tolist()\n",
    "\n",
    "def offset_points(sec,r):\n",
    "    s=seg(sec)\n",
    "    c=[]\n",
    "    for p in s:\n",
    "        c.append(offset_l(p,r)[0])\n",
    "    return array(c).tolist()\n",
    "\n",
    "def offset_pointsv(sec,r):\n",
    "    return array(offset_segv(sec,r))[:,0].tolist()\n",
    "\n",
    "def offset_seg_cw(sec,r):\n",
    "    c=[]\n",
    "    for i in range(len(sec)):\n",
    "        i_minus=len(sec)-1 if i==0 else i-1\n",
    "        i_plus=i+1 if i<len(sec)-1 else 0\n",
    "        p0=sec[i_minus]\n",
    "        p1=sec[i]\n",
    "        p2=sec[i_plus]\n",
    "        clock=cw([p0,p1,p2])\n",
    "        if clock==1:\n",
    "            c.append(offset_l([p1,p2],r))\n",
    "    d=[]\n",
    "    for a in c:\n",
    "        for b in a:\n",
    "            d.append(b)\n",
    "    return d\n",
    "\n",
    "def lim(t,s,e):\n",
    "    return t>=s and t<=e\n",
    "\n",
    "def remove_extra_points(points_list):\n",
    "    return array(points_list)[sort(unique(points_list,axis=0,return_index=True)[1])].tolist()\n",
    "\n",
    "def convert_secv(sec,d):\n",
    "    pi_2minus=sec[-2:]+sec[:-2]\n",
    "    pi_minus=[sec[-1]]+sec[:-1]\n",
    "    p_i=sec\n",
    "    pi_plus=sec[1:]+[sec[0]]\n",
    "    pi_2plus=sec[2:]+sec[:2]\n",
    "\n",
    "    v1=array(pi_minus)-array(pi_2minus)\n",
    "    v2=array(p_i)-array(pi_minus)\n",
    "    v3=array(pi_plus)-array(p_i)\n",
    "    v4=array(pi_2plus)-array(pi_plus)\n",
    "\n",
    "    l1=linalg.norm(v1,axis=1).round(3)\n",
    "    l2=linalg.norm(v2,axis=1).round(3)\n",
    "    l3=linalg.norm(v3,axis=1).round(3)\n",
    "    l4=linalg.norm(v4,axis=1).round(3)\n",
    "\n",
    "    p4=array(pi_2minus)+(array(pi_minus)-array(pi_2minus))/2\n",
    "    p5=array(pi_minus)+(array(p_i)-array(pi_minus))/2\n",
    "\n",
    "    u1=(array(pi_minus)-p4)/linalg.norm(array(pi_minus)-p4,axis=1).reshape(-1,1)\n",
    "    u2=(array(p_i)-p5)/linalg.norm(array(p_i)-p5).reshape(-1,1)\n",
    "\n",
    "    v5=array(pi_minus)-p4\n",
    "    v6=(v5/linalg.norm(v5,axis=1).reshape(-1,1))\n",
    "    r1=r_3pv(array(pi_2minus),array(pi_minus),array(p_i)).round(3)\n",
    "    r2=r_3pv(array(pi_minus),array(p_i),array(pi_plus)).round(3)\n",
    "    r3=r_3pv(array(p_i),array(pi_plus),array(pi_2plus)).round(3)\n",
    "    r=where((l2!=l3) & ((r1!=r2) | (r2!=r3)),0,r2)\n",
    "    arr=swapaxes([pi_minus,p_i,pi_plus],0,1)\n",
    "    clock=array(list(map(cw,arr)))\n",
    "    c1=where(r==0,True,False)\n",
    "    c2=where(r>=d,True,False)\n",
    "    c3=where(clock==1,True,False)\n",
    "    p=array(sec)[c1 | c2 | c3].round(6)\n",
    "    p=p[sort(unique(p,axis=0,return_index=True)[1])]\n",
    "    p1=cKDTree(array(sec)).query(p)[1].tolist()\n",
    "    p2=[p1[len(p1)-1]]+p1[0:len(p1)-1]\n",
    "    p3=p1[1:len(p1)]+[p1[0]]\n",
    "    p4=p1[2:len(p1)]+p1[0:2]\n",
    "    a=i_p2dv(array(sec)[p2],array(sec)[p1],array(sec)[p3],array(sec)[p4])\n",
    "    b=array(sec)[p1]\n",
    "    c=array(p3)-array(p1)>1\n",
    "    d=[]\n",
    "    for i in range(len(c)):\n",
    "        if c[i]==True:\n",
    "            d.append(a[i].tolist())\n",
    "        else:\n",
    "            d.append(b[i].tolist())\n",
    "    d_minus=[d[len(d)-1]]+d[0:len(d)-1]\n",
    "    d_plus=d[1:len(d)]+[d[0]]\n",
    "    va=array(d)-array(d_minus)\n",
    "    vb=array(d_plus)-array(d_minus)\n",
    "    normva=1/linalg.norm(va,axis=1)\n",
    "    normvb=1/linalg.norm(vb,axis=1)\n",
    "    ua=einsum('ij,i->ij',va,normva)\n",
    "    ub=einsum('ij,i->ij',vb,normvb)\n",
    "    return array(d)[(ua!=ub).all(axis=1)].tolist()           \n",
    "\n",
    "\n",
    "def convert_secv1(sec,d):\n",
    "    pi_2minus=sec[-2:]+sec[:-2]\n",
    "    pi_minus=[sec[-1]]+sec[:-1]\n",
    "    p_i=sec\n",
    "    pi_plus=sec[1:]+[sec[0]]\n",
    "    pi_2plus=sec[2:]+sec[:2]\n",
    "\n",
    "    v1=array(pi_minus)-array(pi_2minus)\n",
    "    v2=array(p_i)-array(pi_minus)\n",
    "    v3=array(pi_plus)-array(p_i)\n",
    "    v4=array(pi_2plus)-array(pi_plus)\n",
    "\n",
    "    l1=linalg.norm(v1,axis=1).round(3)\n",
    "    l2=linalg.norm(v2,axis=1).round(3)\n",
    "    l3=linalg.norm(v3,axis=1).round(3)\n",
    "    l4=linalg.norm(v4,axis=1).round(3)\n",
    "\n",
    "    p4=array(pi_2minus)+(array(pi_minus)-array(pi_2minus))/2\n",
    "    p5=array(pi_minus)+(array(p_i)-array(pi_minus))/2\n",
    "\n",
    "    u1=(array(pi_minus)-p4)/linalg.norm(array(pi_minus)-p4,axis=1).reshape(-1,1)\n",
    "    u2=(array(p_i)-p5)/linalg.norm(array(p_i)-p5).reshape(-1,1)\n",
    "\n",
    "    v5=array(pi_minus)-p4\n",
    "    v6=(v5/linalg.norm(v5,axis=1).reshape(-1,1))\n",
    "    r1=r_3pv(array(pi_2minus),array(pi_minus),array(p_i)).round(3)\n",
    "    r2=r_3pv(array(pi_minus),array(p_i),array(pi_plus)).round(3)\n",
    "    r3=r_3pv(array(p_i),array(pi_plus),array(pi_2plus)).round(3)\n",
    "    r=where((l2!=l3) & ((r1!=r2) | (r2!=r3)),0,r2)\n",
    "    arr=swapaxes([pi_minus,p_i,pi_plus],0,1)\n",
    "    clock=array(list(map(cw,arr)))\n",
    "    c1=where(r==0,True,False)\n",
    "    c2=where(r>=d,True,False)\n",
    "    c3=where(clock==-1,True,False)\n",
    "    p=array(sec)[c1 | c2 | c3]\n",
    "    p1=cKDTree(array(sec)).query(p)[1].tolist()\n",
    "    p2=[p1[len(p1)-1]]+p1[0:len(p1)-1]\n",
    "    p3=p1[1:len(p1)]+[p1[0]]\n",
    "    p4=p1[2:len(p1)]+p1[0:2]\n",
    "    a=i_p2dv(array(sec)[p2],array(sec)[p1],array(sec)[p3],array(sec)[p4])\n",
    "    b=array(sec)[p1]\n",
    "    c=array(p3)-array(p1)>1\n",
    "    d=[]\n",
    "    for i in range(len(c)):\n",
    "        if c[i]==True:\n",
    "            d.append(a[i].tolist())\n",
    "        else:\n",
    "            d.append(b[i].tolist())\n",
    "    d_minus=[d[len(d)-1]]+d[0:len(d)-1]\n",
    "    d_plus=d[1:len(d)]+[d[0]]\n",
    "    va=array(d)-array(d_minus)\n",
    "    vb=array(d_plus)-array(d_minus)\n",
    "    normva=1/linalg.norm(va,axis=1)\n",
    "    normvb=1/linalg.norm(vb,axis=1)\n",
    "    ua=einsum('ij,i->ij',va,normva)\n",
    "    ub=einsum('ij,i->ij',vb,normvb)\n",
    "    return array(d)[(ua!=ub).all(axis=1)].tolist()           \n",
    "\n",
    "\n",
    "def list_r(sec):\n",
    "    pi_2minus=sec[-2:]+sec[:-2]\n",
    "    pi_minus=[sec[-1]]+sec[:-1]\n",
    "    p_i=sec\n",
    "    pi_plus=sec[1:]+[sec[0]]\n",
    "    pi_2plus=sec[2:]+sec[:2]\n",
    "\n",
    "    v1=array(pi_minus)-array(pi_2minus)\n",
    "    v2=array(p_i)-array(pi_minus)\n",
    "    v3=array(pi_plus)-array(p_i)\n",
    "    v4=array(pi_2plus)-array(pi_plus)\n",
    "\n",
    "    l1=linalg.norm(v1,axis=1).round(3)\n",
    "    l2=linalg.norm(v2,axis=1).round(3)\n",
    "    l3=linalg.norm(v3,axis=1).round(3)\n",
    "    l4=linalg.norm(v4,axis=1).round(3)\n",
    "\n",
    "    p4=array(pi_2minus)+(array(pi_minus)-array(pi_2minus))/2\n",
    "    p5=array(pi_minus)+(array(p_i)-array(pi_minus))/2\n",
    "\n",
    "    u1=(array(pi_minus)-p4)/linalg.norm(array(pi_minus)-p4,axis=1).reshape(-1,1)\n",
    "    u2=(array(p_i)-p5)/linalg.norm(array(p_i)-p5).reshape(-1,1)\n",
    "\n",
    "    v5=array(pi_minus)-p4\n",
    "    v6=(v5/linalg.norm(v5,axis=1).reshape(-1,1))\n",
    "    r1=r_3pv(array(pi_2minus),array(pi_minus),array(p_i)).round(3)\n",
    "    r2=r_3pv(array(pi_minus),array(p_i),array(pi_plus)).round(3)\n",
    "    r3=r_3pv(array(p_i),array(pi_plus),array(pi_2plus)).round(3)\n",
    "    r=where((l2!=l3) & ((r1!=r2) | (r2!=r3)),0,r2)\n",
    "    return r\n",
    "\n",
    "def list_ra(sec):\n",
    "    pi_2minus=sec[-2:]+sec[:-2]\n",
    "    pi_minus=[sec[-1]]+sec[:-1]\n",
    "    p_i=sec\n",
    "    pi_plus=sec[1:]+[sec[0]]\n",
    "    pi_2plus=sec[2:]+sec[:2]\n",
    "\n",
    "    v1=array(pi_minus)-array(pi_2minus)\n",
    "    v2=array(p_i)-array(pi_minus)\n",
    "    v3=array(pi_plus)-array(p_i)\n",
    "    v4=array(pi_2plus)-array(pi_plus)\n",
    "\n",
    "    l1=linalg.norm(v1,axis=1).round(3)\n",
    "    l2=linalg.norm(v2,axis=1).round(3)\n",
    "    l3=linalg.norm(v3,axis=1).round(3)\n",
    "    l4=linalg.norm(v4,axis=1).round(3)\n",
    "\n",
    "    p4=array(pi_2minus)+(array(pi_minus)-array(pi_2minus))/2\n",
    "    p5=array(pi_minus)+(array(p_i)-array(pi_minus))/2\n",
    "\n",
    "    u1=(array(pi_minus)-p4)/linalg.norm(array(pi_minus)-p4,axis=1).reshape(-1,1)\n",
    "    u2=(array(p_i)-p5)/linalg.norm(array(p_i)-p5).reshape(-1,1)\n",
    "\n",
    "    v5=array(pi_minus)-p4\n",
    "    v6=(v5/linalg.norm(v5,axis=1).reshape(-1,1))\n",
    "    r1=r_3pv(array(pi_2minus),array(pi_minus),array(p_i)).round(3)\n",
    "    r2=r_3pv(array(pi_minus),array(p_i),array(pi_plus)).round(3)\n",
    "    r3=r_3pv(array(p_i),array(pi_plus),array(pi_2plus)).round(3)\n",
    "    r=where((l2!=l3) & ((r1!=r2) | (r2!=r3)),0,r2)\n",
    "    return r2\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "def rnd_v(v,n):\n",
    "    b=[]\n",
    "    for i in v:\n",
    "        b.append(round(i,n))\n",
    "    return b\n",
    "\n",
    "def i_m2d(m):\n",
    "    return linalg.pinv(transpose(m)).tolist()\n",
    "\n",
    "def rm(theta):\n",
    "    pi=3.141592653589793\n",
    "    return [[cos(theta * pi/180),sin(theta * pi/180)],[-sin(theta * pi/180),cos(theta * pi/180)]]\n",
    "\n",
    "def max_rv(sec):\n",
    "    pi_2minus=sec[-2:]+sec[:-2]\n",
    "    pi_minus=[sec[-1]]+sec[:-1]\n",
    "    p_i=sec\n",
    "    pi_plus=sec[1:]+[sec[0]]\n",
    "    pi_2plus=sec[2:]+sec[:2]\n",
    "\n",
    "    v1=array(pi_minus)-array(pi_2minus)\n",
    "    v2=array(p_i)-array(pi_minus)\n",
    "    v3=array(pi_plus)-array(p_i)\n",
    "    v4=array(pi_2plus)-array(pi_plus)\n",
    "\n",
    "    l1=linalg.norm(v1,axis=1).round(3)\n",
    "    l2=linalg.norm(v2,axis=1).round(3)\n",
    "    l3=linalg.norm(v3,axis=1).round(3)\n",
    "    l4=linalg.norm(v4,axis=1).round(3)\n",
    "\n",
    "    p4=array(pi_2minus)+(array(pi_minus)-array(pi_2minus))/2\n",
    "    p5=array(pi_minus)+(array(p_i)-array(pi_minus))/2\n",
    "\n",
    "    u1=(array(pi_minus)-p4)/linalg.norm(array(pi_minus)-p4,axis=1).reshape(-1,1)\n",
    "    u2=(array(p_i)-p5)/linalg.norm(array(p_i)-p5).reshape(-1,1)\n",
    "\n",
    "    v5=array(pi_minus)-p4\n",
    "    v6=(v5/linalg.norm(v5,axis=1).reshape(-1,1))\n",
    "    r1=r_3pv(array(pi_2minus),array(pi_minus),array(p_i)).round(3)\n",
    "    r2=r_3pv(array(pi_minus),array(p_i),array(pi_plus)).round(3)\n",
    "    r3=r_3pv(array(p_i),array(pi_plus),array(pi_2plus)).round(3)\n",
    "    return max(where((l2!=l3) & ((r1!=r2) | (r2!=r3)),0,r2))\n",
    "\n",
    "def r_3p(p):\n",
    "    p4=add(p[0],divide(subtract(p[1],p[0]),2)).tolist()\n",
    "    p5=add(p[1],divide(subtract(p[2],p[1]),2)).tolist()\n",
    "    u1=uv(subtract(p[1],p4))\n",
    "    u2=uv(subtract(p[2],p5))\n",
    "    p6=add(p4,dot(u1,rm(90))).tolist()\n",
    "    p7=add(p5,dot(u2,rm(90))).tolist()\n",
    "    cp=i_p2d([p4,p6],[p5,p7])\n",
    "    r=norm(subtract(p[0],cp))\n",
    "    return r\n",
    "\n",
    "\n",
    "def i_p2d(l1,l2):\n",
    "    p0,p1,p2,p3=l1[0],l1[1],l2[0],l2[1]\n",
    "    p0,p1,p2,p3=array([p0,p1,p2,p3])\n",
    "    v1=p1-p0\n",
    "    v2=p3-p2\n",
    "    im=linalg.pinv(array([v1,-v2]).T)\n",
    "    t1=(im@(p2-p0))[0]\n",
    "    ip=p0+v1*t1\n",
    "    \n",
    "    return ip.tolist()\n",
    "\n",
    "\n",
    "def offset_seg_cwv(sec,r):\n",
    "    pi_minus=[sec[-1]]+sec[:-1]\n",
    "    p_i=sec\n",
    "    pi_plus=sec[1:]+[sec[0]]\n",
    "    c=array(list(map(cw,swapaxes([pi_minus,p_i,pi_plus],0,1))))\n",
    "    return array(offset_segv(sec,r))[c==1].reshape(-1,2)                    \n",
    "            \n",
    "\n",
    "\n",
    "def s_intv(s):\n",
    "    c=[]\n",
    "    for i in range(len(s)):\n",
    "        p0=array([s[i]]*len(s))[:,0]\n",
    "        p1=array([s[i]]*len(s))[:,1]\n",
    "        v1=p1-p0\n",
    "        p2=array(s)[:,0]\n",
    "        p3=array(s)[:,1]\n",
    "        v2=p3-p2\n",
    "        m=swapaxes([swapaxes([v1.T[0],-v2.T[0]],0,1),swapaxes([v1.T[1],-v2.T[1]],0,1)],0,1)\n",
    "        n=m[where(linalg.det(m)!=0)]\n",
    "        pa=p0[where(linalg.det(m)!=0)]\n",
    "        pb=p2[where(linalg.det(m)!=0)]\n",
    "        v=v1[where(linalg.det(m)!=0)]\n",
    "        A=linalg.inv(n)\n",
    "        B=pb-pa\n",
    "        def mul(a,b):\n",
    "            return a@b\n",
    "        t=einsum('ijk,ik->ij',A,B)[:,0].round(4)\n",
    "        u=einsum('ijk,ik->ij',A,B)[:,1].round(4)\n",
    "        t1=where(t>=0,where(t<=1,True,False),False)\n",
    "        u1=where(u>=0,where(u<=1,True,False),False)\n",
    "        d=(pa+v*t.reshape(-1,1))[where(t1&u1==True)].tolist()\n",
    "        if d!=[]:\n",
    "            c=c+d\n",
    "    return c\n",
    "\n",
    "def s_intv1(s):\n",
    "    c=[]\n",
    "    for i in range(len(s)):\n",
    "        p0=array([s[i]]*len(s))[:,0]\n",
    "        p1=array([s[i]]*len(s))[:,1]\n",
    "        v1=p1-p0\n",
    "        p2=array(s)[:,0]\n",
    "        p3=array(s)[:,1]\n",
    "        v2=p3-p2\n",
    "        m=swapaxes([swapaxes([v1.T[0],-v2.T[0]],0,1),swapaxes([v1.T[1],-v2.T[1]],0,1)],0,1)\n",
    "        n=m[where(linalg.det(m)!=0)]\n",
    "        pa=p0[where(linalg.det(m)!=0)]\n",
    "        pb=p2[where(linalg.det(m)!=0)]\n",
    "        v=v1[where(linalg.det(m)!=0)]\n",
    "        A=linalg.inv(n)\n",
    "        B=pb-pa\n",
    "        def mul(a,b):\n",
    "            return a@b\n",
    "        t=einsum('ijk,ik->ij',A,B)[:,0].round(4)\n",
    "        u=einsum('ijk,ik->ij',A,B)[:,1].round(4)\n",
    "        t1=where(t>0,where(t<1,True,False),False)\n",
    "        u1=where(u>0,where(u<1,True,False),False)\n",
    "        d=(pa+v*t.reshape(-1,1))[where(t1&u1==True)].tolist()\n",
    "        if d!=[]:\n",
    "            c=c+d\n",
    "    return c\n",
    "\n",
    "\n",
    "def r_3pv(p1,p2,p3):\n",
    "    p4=p1+(p2-p1)/2\n",
    "    p5=p2+(p3-p2)/2\n",
    "    u1=(p2-p4)/linalg.norm(p2-p4,axis=1).reshape(-1,1)\n",
    "    u2=(p3-p5)/linalg.norm(p3-p5,axis=1).reshape(-1,1)\n",
    "    p6=p4+u1@array([[0,1],[-1,0]])\n",
    "    p7=p5+u2@array([[0,1],[-1,0]])\n",
    "    cp=i_p2dv(p4,p6,p5,p7)\n",
    "    r=linalg.norm(p1-cp,axis=1)\n",
    "    return r\n",
    "\n",
    "def i_p2dv(p0,p1,p2,p3):\n",
    "    v1=p1-p0\n",
    "    v2=p3-p2\n",
    "    a=linalg.pinv(swapaxes(transpose(array([v1,-v2])),0,1))\n",
    "    b=p2-p0\n",
    "    t=einsum('ijk,ik->ij',a,b)[:,0]\n",
    "    return p0+einsum('ij,i->ij',v1,t)\n",
    "\n",
    "def sort_points(sec,list):\n",
    "    if list!=[]:\n",
    "        b=[]\n",
    "        for p in sec:\n",
    "            a=[]\n",
    "            for i in range(len(list)):\n",
    "                a.append(norm(subtract(list[i],p)))\n",
    "            for i,x in enumerate(a):\n",
    "                if x==min(a):\n",
    "                    b.append(list[i])\n",
    "            \n",
    "        return b\n",
    "            \n",
    "def sort_pointsv(sec,sec1):\n",
    "    a=array(sec)\n",
    "    b=array(sec1)\n",
    "    c=[]\n",
    "    for p in a:\n",
    "        d=linalg.norm(b-p,axis=1)\n",
    "        c.append(b[where(d==min(d))][0])\n",
    "    return array(c).tolist()\n",
    "\n",
    "\n",
    "\n",
    "def m_points(sec,sl=20):# multiple points within straight lines of a closed section 'sec' with equal segment length 'sl' in the straight line segments\n",
    "    p0=array(sec)\n",
    "    p1=array(sec)[1:].tolist()+[sec[0]]\n",
    "    lnth=linalg.norm(array(p1)-array(p0),axis=1)\n",
    "    sec1=concatenate([array(l([p0[i],p1[i]],lnth[i]/sl)) if lnth[i]>=sl*2 else [p0[i]] for i in range(len(p0))])\n",
    "    return sec1.tolist()\n",
    "\n",
    "def m_points_o(sec,sl=20):# multiple points within straight lines of an open section 'sec' with equal segment length 'sl' in the straight line segments\n",
    "    p0=array(sec)\n",
    "    p1=array(sec)[1:].tolist()+[sec[0]]\n",
    "    lnth=linalg.norm(array(p1)-array(p0),axis=1)\n",
    "    sec1=concatenate([array(l([p0[i],p1[i]],lnth[i]/sl)) if lnth[i]>=sl*2 else [p0[i]] for i in range(len(p0)-1)])\n",
    "    return sec1.tolist()\n",
    "\n",
    "\n",
    "def l(l,s=20):# line 'l' with number of segments 's'\n",
    "    p0,p1=array(l[0]),array(l[1])\n",
    "    v=p1-p0\n",
    "    u=[v/linalg.norm(v)]\n",
    "    length=linalg.norm(v)\n",
    "    r=arange(0,length,length/s)\n",
    "    return (p0+einsum('ij,k->kj',u,r)).tolist()\n",
    "\n",
    "def l_len(l):# length of a line 'l'\n",
    "    p0,p1=array(l[0]),array(l[1])\n",
    "    v=p1-p0\n",
    "    u=[v/(linalg.norm(v)+.00001)]\n",
    "    length=linalg.norm(v)\n",
    "    return length.tolist()\n",
    "\n",
    "def arc_2p(p1,p2,r,cw=1,s=20):#arc with 2 points 'p1,p2' with radius 'r' and with orientation clockwise (1) or counterclock wise(-1)\n",
    "    p1,p2=array([p1,p2])\n",
    "    p3=p1+(p2-p1)/2\n",
    "    d=linalg.norm(p3-p1)\n",
    "    l=sqrt(abs(r**2-d**2))\n",
    "    v=p1-p3\n",
    "    u=v/linalg.norm(v)\n",
    "    cp=p3+(u*l)@rm(-90 if cw==-1 else 90)\n",
    "    v1,v2=p1-cp,p2-cp\n",
    "    a1,a2=ang(v1[0],v1[1]),ang(v2[0],v2[1])\n",
    "    a3= (a2+360 if a2<a1 else a2) if cw==-1 else (a2 if a2<a1 else a2-360)\n",
    "    return arc(r,a1,a3,cp,s)\n",
    "\n",
    "def arc_long_2p(p1,p2,r,cw=1,s=20):#long arc with 2 points 'p1,p2' with radius 'r' and with orientation clockwise (1) or counterclock wise(-1)\n",
    "    p1,p2=array([p1,p2])\n",
    "    p3=p1+(p2-p1)/2\n",
    "    d=linalg.norm(p3-p1)\n",
    "    l=sqrt(abs(r**2-d**2))\n",
    "    v=p1-p3\n",
    "    u=v/linalg.norm(v)\n",
    "    cp=p3+(u*l)@rm(90 if cw==-1 else -90)\n",
    "    v1,v2=p1-cp,p2-cp\n",
    "    a1,a2=ang(v1[0],v1[1]),ang(v2[0],v2[1])\n",
    "    a3=(a2+360 if a2<a1 else a2) if cw==-1 else (a2 if a2<a1 else a2-360)\n",
    "    return arc(r,a1,a3,cp,s)\n",
    "\n",
    "def arc_2p_cp(p1,p2,r,cw=-1):# center point of an arc with 2 points 'p1,p2' with radius 'r' and with orientation clockwise (1) or counterclock wise(-1)\n",
    "    p1,p2=array([p1,p2])\n",
    "    p3=p1+(p2-p1)/2\n",
    "    d=linalg.norm(p3-p1)\n",
    "    l=sqrt(abs(r**2-d**2))\n",
    "    v=p1-p3\n",
    "    u=v/linalg.norm(v)\n",
    "    cp=p3+(u*l)@rm(-90 if cw==-1 else 90)\n",
    "    return cp\n",
    "\n",
    "def offset(sec,r):# offset for a section 'sec' by amount 'r'\n",
    "#     return io(sec,r) if r<0 else sec if r==0 else oo_convex(sec,r) if convex(sec)==True else outer_offset(sec,r)\n",
    "    return inner_offset(sec,r) if r<0 else sec if r==0 else oo_convex(sec,r) if convex(sec)==True else out_offset(sec,r)\n",
    "\n",
    "\n",
    "def prism(sec,path):# prism with a section 'sec' and sculpting 'path'\n",
    "    return [array(trns([0,0,y],offset(sec,round(x,3)))).tolist() for (x,y) in path]\n",
    "\n",
    "def trns(p,sec):#translates a prism or section by [x,y,z] distance\n",
    "    return [ (array([p1[0],p1[1],0])+array(p)).tolist() if len(p1)==2 else (array(p1)+array(p)).tolist() for p1 in sec]\n",
    "\n",
    "def prism1(sec,path,n):\n",
    "        a=m_points(sec,n)\n",
    "        return [ trns([0,0,y], array(m_points(offset(sec,x),n))[cKDTree(m_points(offset(sec,x),n)).query(a)[1]]) for (x,y) in path ]\n",
    "\n",
    "def offset_points_cw(sec,r):\n",
    "    s=seg(sec)\n",
    "    c=[]\n",
    "    for i in range(len(sec)):\n",
    "        i_minus=len(sec)-1 if i==0 else i-1\n",
    "        i_plus=i+1 if i<len(sec)-1 else 0\n",
    "        p0=sec[i_minus]\n",
    "        p1=sec[i]\n",
    "        p2=sec[i_plus]\n",
    "        if cw([p0,p1,p2])==1:\n",
    "            c.append(offset_l([p1,p2],r)[0])\n",
    "    return c\n",
    "\n",
    "\n",
    "def cytz(path):# converts 'y' points to 'z' points in a 2d list of points\n",
    "    return [[p[0],0,p[1]] for p in path]\n",
    "\n",
    "def surf_extrude(sec,path):# extrudes an open section 'sec' to a 'path' to create surface\n",
    "    p0=path\n",
    "    p1=p0[1:]+[p0[0]]\n",
    "    p0,p1=array(p0),array(p1)\n",
    "    v=p1-p0\n",
    "    a1=vectorize(ang)(v[:,0],v[:,1])\n",
    "    b=sqrt(v[:,0]**2+v[:,1]**2)\n",
    "    a2=vectorize(ang)(b,v[:,2])\n",
    "    c=[]\n",
    "    for i in range(len(path)-1):\n",
    "        sec1=trns(p0[i],q_rot(['x90','z-90',f'y{-a2[i]}',f'z{a1[i]}'],sec))\n",
    "        sec2=trns(p1[i],q_rot(['x90','z-90',f'y{-a2[i]}',f'z{a1[i]}'],sec))\n",
    "        if i<len(path)-2:\n",
    "            c.append([sec1])\n",
    "        else:\n",
    "            c.append([sec1,sec2])\n",
    "    return concatenate(c).tolist()\n",
    "\n",
    "def cpo(prism): # changes the orientation of points of a prism\n",
    "    return swapaxes(array(prism),0,1).tolist()\n",
    "\n",
    "def c2t3(p):# converts 2d list to 3d\n",
    "    if len(array(p).shape)>2:\n",
    "        return [trns([0,0,0],p) for p in p]\n",
    "    else:\n",
    "        return trns([0,0,0],p)\n",
    "\n",
    "def c3t2(a): # converts 3d list to 2d list \n",
    "    if len(array(a).shape)==3:\n",
    "        return array([ swapaxes([p[:,0],p[:,1]],0,1) for p in array(a)]).tolist()\n",
    "    else:\n",
    "        p=array(a)\n",
    "        return swapaxes([p[:,0],p[:,1]],0,1).tolist()\n",
    "\n",
    "def nv(p):# normal vector to the plane 'p' with atleast 3 known points\n",
    "    p0,p1,p2=array(trns([0,0,0],[p[0],p[1],p[2]]))\n",
    "    nv=cross(p0-p1,p2-p1)\n",
    "    m=1/linalg.norm(nv) if linalg.norm(nv)>0 else 1e5\n",
    "    return (nv*m).tolist()\n",
    "\n",
    "def fillet_3p_3d(p0,p1,p2,r,s):# fillet with 3 known points 'p0,p1,p2' in 3d space. 'r' is the radius of fillet and 's' is the number of segments in the fillet\n",
    "    p0,p1,p2=array(trns([0,0,0],[p0,p1,p2]))\n",
    "    n=array(nv([p0,p1,p2]))\n",
    "    u1=(p0-p1)/(linalg.norm(p0-p1)+.00001)\n",
    "    u2=(p2-p1)/(linalg.norm(p2-p1)+.00001)\n",
    "    theta=(180-arccos(u1@u2)*180/pi)/2\n",
    "    alpha=arccos(u1@u2)*180/pi\n",
    "    l=r*tan(theta*pi/180)\n",
    "    cp=p1+q(n,u1*r/cos(theta*pi/180),alpha/2)\n",
    "    pa=p1+u1*l\n",
    "    arc=[ cp+q(n,pa-cp,-i) for i in linspace(0,theta*2,s)]\n",
    "    a,b,c=arc[0],arc[1:s-1],arc[s-1]\n",
    "    return concatenate([[p1],arc]).tolist()\n",
    "\n",
    "def fillet_3p_3d_cp(p0,p1,p2,r):# center point 'cp' of the fillet with 3 known points 'p0,p1,p2' in 3d space. 'r' is the radius of fillet\n",
    "    p0,p1,p2=array(trns([0,0,0],[p0,p1,p2]))\n",
    "    n=array(nv([p0,p1,p2]))\n",
    "    u1=(p0-p1)/(linalg.norm(p0-p1)+.00001)\n",
    "    u2=(p2-p1)/(linalg.norm(p2-p1)+.00001)\n",
    "    theta=(180-arccos(u1@u2)*180/pi)/2\n",
    "    alpha=arccos(u1@u2)*180/pi\n",
    "    l=r*tan(theta*pi/180)\n",
    "    cp=p1+q(n,u1*r/cos(theta*pi/180),alpha/2)\n",
    "    return cp.tolist()\n",
    "\n",
    "def i_p3d(l1,l2): # intersection point between 2 lines 'l1' and 'l2' in 3d space where both the lines are in the same plane\n",
    "    l1,l2=array(l1),array(l2)\n",
    "    v1=l1[1]-l1[0]\n",
    "    v2=l2[1]-l2[0]\n",
    "    u1=v1/(linalg.norm(v1)+.00001)\n",
    "    u2=v2/(linalg.norm(v2)+.00001)\n",
    "    v3=l2[0]-l1[0]\n",
    "    t1= (linalg.pinv(array([v1,-v2,[1,1,1]]).T)@array(v3))[0]\n",
    "    ip=l1[0]+v1*t1\n",
    "    return ip.tolist()\n",
    "\n",
    "def arc_3p_3d(points,s): # arc with 3 known list of 'points' in 3d space where 's' is the number of segments in the arc\n",
    "    points=array(points)\n",
    "    v1=points[0]-points[1]\n",
    "    v2=points[2]-points[1]\n",
    "    u1=v1/linalg.norm(v1)\n",
    "    u2=v2/linalg.norm(v2)\n",
    "    n=cross(u1,u2)\n",
    "    alpha=arccos(u1@u2)*180/pi\n",
    "    pa=v1/2\n",
    "    pb=v2/2\n",
    "    pap=pa+q(n,u1,90)\n",
    "    pbp=pb+q(n,u2,-90)\n",
    "    l1=[pa,pap]\n",
    "    l2=[pb,pbp]\n",
    "    cp=i_p3d(l1,l2)\n",
    "    v3=points[0]-(points[1]+cp)\n",
    "    u3=v3/linalg.norm(v3)\n",
    "    v4=points[2]-(points[1]+cp)\n",
    "    u4=v4/linalg.norm(v4)\n",
    "    theta= 360-arccos(u3@u4)*180/pi if alpha<90 else arccos(u3@u4)*180/pi\n",
    "    radius=linalg.norm(pa-cp)\n",
    "    arc=trns(points[1]+cp,[ q(n,points[0]-(points[1]+cp),-i)  for i in linspace(0,theta,s) ])\n",
    "    return array(arc).tolist()\n",
    "\n",
    "def r_3p_3d(points):# radius of the circle with 3 known list of 'points' in 3d space\n",
    "    points=array(points)\n",
    "    v1=points[0]-points[1]\n",
    "    v2=points[2]-points[1]\n",
    "    u1=v1/(linalg.norm(v1)+.00001)\n",
    "    u2=v2/(linalg.norm(v2)+.00001)\n",
    "    n=cross(u1,u2)\n",
    "    alpha=arccos(u1@u2)*180/pi\n",
    "    pa=v1/2\n",
    "    pb=v2/2\n",
    "    pap=pa+q(n,u1,90)\n",
    "    pbp=pb+q(n,u2,-90)\n",
    "    l1=[pa,pap]\n",
    "    l2=[pb,pbp]\n",
    "    cp=i_p3d(l1,l2)\n",
    "    v3=points[0]-(points[1]+cp)\n",
    "    u3=v3/(linalg.norm(v3)+.00001)\n",
    "    v4=points[2]-(points[1]+cp)\n",
    "    u4=v4/(linalg.norm(v4)+.00001)\n",
    "    theta= 360-arccos(u3@u4)*180/pi if alpha<90 else arccos(u3@u4)*180/pi\n",
    "    radius=linalg.norm(pa-cp)\n",
    "    return radius\n",
    "\n",
    "def cir_3p_3d(points,s):#circle with 3 known list of 'points' in 3d space where 's' is the number of segments in the circle \n",
    "    points=array(points)\n",
    "    v1=points[0]-points[1]\n",
    "    v2=points[2]-points[1]\n",
    "    u1=v1/linalg.norm(v1)\n",
    "    u2=v2/linalg.norm(v2)\n",
    "    n=cross(u1,u2)\n",
    "    alpha=arccos(u1@u2)*180/pi\n",
    "    pa=v1/2\n",
    "    pb=v2/2\n",
    "    pap=pa+q(n,u1,90)\n",
    "    pbp=pb+q(n,u2,-90)\n",
    "    l1=[pa,pap]\n",
    "    l2=[pb,pbp]\n",
    "    cp=i_p3d(l1,l2)\n",
    "    v3=points[0]-(points[1]+cp)\n",
    "    u3=v3/linalg.norm(v3)\n",
    "    v4=points[2]-(points[1]+cp)\n",
    "    u4=v4/linalg.norm(v4)\n",
    "    theta= 360-arccos(u3@u4)*180/pi if alpha<90 else arccos(u3@u4)*180/pi\n",
    "    radius=linalg.norm(pa-cp)\n",
    "    arc=trns(points[1]+cp,[ q(n,points[0]-(points[1]+cp),-i)  for i in linspace(0,360,s) ])\n",
    "    return array(arc).tolist()\n",
    "\n",
    "def scl2d(sec,sl):# scale the 2d section 'sec' by a scaling factor 'sl'. this places the scaled section in the bottom center of the original section\n",
    "    s1=array(trns([0,0,0],sec))\n",
    "    cp=array(s1).mean(axis=0)\n",
    "    rev=array(s1).mean(axis=0)+(array(s1)-array(s1).mean(axis=0))*sl\n",
    "    y1=cp-array([0,array(s1)[:,1].min(),0])\n",
    "    y2=cp-array([0,rev[:,1].min(),0])\n",
    "    d=y2-y1\n",
    "    return c3t2(trns(d,rev))\n",
    "\n",
    "def scl2d_c(sec,sl):# scale the 2d section 'sec' with scaling factor 'sl'. this places the scaled section in the center of original section or the center of both original and scaled section remains the same.\n",
    "    s1=array(trns([0,0,0],sec))\n",
    "    cp=array(s1).mean(axis=0)\n",
    "    rev=array(s1).mean(axis=0)+(array(s1)-array(s1).mean(axis=0))*sl\n",
    "    return c3t2(rev)\n",
    "\n",
    "def scl3d(p,s):# scale 3d prism 'p' with scaling factor 's'. This places the scaled prism at the same bottom of the original prism\n",
    "    p=array(p)\n",
    "    cp=p.reshape(-1,3).mean(axis=0)\n",
    "    rev=cp+(p-cp)*s\n",
    "    z1=p.reshape(-1,3)[:,2].min()\n",
    "    z2=rev.reshape(-1,3)[:,2].min()\n",
    "    d=z1-z2\n",
    "    return trns([0,0,d],rev)\n",
    "\n",
    "def scl3dc(p,s):# scale a 3d prism 'p' with scaling factor 's'. This places the scaled prism in the center of the original prism or the center of both the prism is same\n",
    "    p=array(p)\n",
    "    cp=p.reshape(-1,3).mean(axis=0)\n",
    "    rev=cp+(p-cp)*s\n",
    "    return rev.tolist()\n",
    "\n",
    "\n",
    "def io(sec,r):# used for inner offset in offset function\n",
    "    if r<0:\n",
    "        s=flip(sec) if cw(sec)==1 else sec\n",
    "        s1=s\n",
    "#         s1=convert_secv(s,max_r(s)+1 if abs(r)>=max_r(s) else abs(r))\n",
    "        s2=offset_seg(s1,r)\n",
    "        s3=offset_seg_cw(s1,r)\n",
    "        s4=s_int(s2)\n",
    "        s5=sec_clean(s1,s4+s3,abs(r))\n",
    "        s6=array(s5)[cKDTree(s5).query(s)[1]]\n",
    "        return s6.tolist()\n",
    "\n",
    "\n",
    "    \n",
    "def outer_offset(sec,r):# used for offset function\n",
    "    s1=flip(sec) if cw(sec)==1 else sec\n",
    "    p0=[sec[len(sec)-1]]+sec[:len(sec)-1]\n",
    "    p1=sec\n",
    "    p2=sec[1:]+[sec[0]]\n",
    "    p0,p1,p2=array([p0,p1,p2])\n",
    "    v1=p0-p1\n",
    "    u1=v1/linalg.norm(v1,axis=1).reshape(-1,1)\n",
    "    v2=p2-p1\n",
    "    u2=v2/linalg.norm(v2,axis=1).reshape(-1,1)\n",
    "    theta=arccos(einsum('ij,ij->i',u1,u2))*180/pi\n",
    "    alpha=180-theta\n",
    "    pa=p1+einsum('ij,i->ij',(u1*r),tan(alpha/2*pi/180))\n",
    "    pb=p1+einsum('ij,i->ij',(u2*r),tan(alpha/2*pi/180))\n",
    "    cp=array([ arc_2p_cp(pa[i],pb[i],r,1) for i in range(len(p1))])\n",
    "    pc=p1+(u1@rm(90))*r\n",
    "    pd=p1+(u2@rm(-90))*r\n",
    "    op=[ array(arc_2p(pc[i],pd[i],r,-1,0 if linalg.norm(pc[i]-pd[i])<1 else 5)) if cw([p0[i],p1[i],p2[i]])==-1 else [cp[i]] for i in range(len(p1))]\n",
    "    radius=r_3pv(p0,p1,p2)\n",
    "    op01=concatenate([op[i] for i in range(len(sec)) if (cw([p0[i],p1[i],p2[i]])==-1) | (radius[i]>=r)]).tolist()\n",
    "    p0=op01\n",
    "    p1=op01[1:]+[op01[0]]\n",
    "    p2=op01[len(op01)-2:len(op01)]+op01[:len(op01)-2]\n",
    "    p3=[op01[len(op01)-1]]+op01[:len(op01)-1]\n",
    "    p4=op01[2:]+op01[0:2]\n",
    "    p5=op01[3:]+op01[0:3]\n",
    "\n",
    "    p0,p1,p2,p3,p4,p5=array([p0,p1,p2,p3,p4,p5])\n",
    "    v1=p1-p0\n",
    "    u1=v1/linalg.norm(v1,axis=1).reshape(-1,1)\n",
    "    ip=swapaxes(array([i_p2dv(p0,p1,p2,p3),i_p2dv(p0,p1,p4,p5)]),0,1)\n",
    "    l1=linalg.norm(p1-p0,axis=1)\n",
    "    a=ip-p0[:,None]\n",
    "    b=1/sqrt(einsum('ijk,ijk->ij',a,a))\n",
    "    c=sqrt(einsum('ijk,ijk->ij',a,a))\n",
    "    u2=einsum('ijk,ij->ijk',a,b)\n",
    "    c1=c<l1[:,None]\n",
    "    c2=((u2<0)==(u1<0)[:,None]).all(axis=1)\n",
    "    op02=[ ip[i][0].tolist() if (c1&c2)[i][0]==True else (ip[i][1].tolist() if (c1&c2)[i][1]==True else op01[i]) for i in range(len(ip))]\n",
    "    p=array(op02)\n",
    "    p1=array(m_points1(s1,10))\n",
    "    p.shape,p1.shape\n",
    "    p2=p[:,None]-p1\n",
    "    p3=sqrt(einsum('ijk,ijk->ij',p2,p2)).min(axis=1)\n",
    "    p4=p[(p3>=r-.02)]\n",
    "\n",
    "    return p4[cKDTree(p4).query(s1)[1]].tolist()\n",
    "\n",
    "def m_points1(sec,s):# multiple points with in the straight lines in the closed section 'sec'. 's' is the number of segments between each straight line\n",
    "    s1=sec\n",
    "    s2=sec[1:]+[sec[0]]\n",
    "    s1,s2=array([s1,s2])\n",
    "    u=(s2-s1)/linalg.norm(s2-s1,axis=1).reshape(-1,1)\n",
    "    l=linalg.norm(s2-s1,axis=1)\n",
    "    n=(l/s).round(0)+1\n",
    "    p=linspace(zeros(len(l)),l,s,axis=1)\n",
    "    q=einsum('ij,ik->ikj',u,p)\n",
    "    s1.shape,q.shape\n",
    "    return (s1[:,None]+q).reshape(-1,2).tolist()\n",
    "\n",
    "def ibsap(sec,pnt):# intersection between section and a point. used to find whether the poin is inside the section or outside the section\n",
    "    p0=array(pnt)\n",
    "    p2=sec\n",
    "    p3=sec[1:]+[sec[0]]\n",
    "    p2,p3=array([p2,p3])\n",
    "    v1=[1,0]\n",
    "    v2=(p3-p2)+[0,.00001]\n",
    "    im=linalg.pinv(array([[v1]*len(v2),-v2]).transpose(1,0,2).transpose(0,2,1))\n",
    "    p=p2-p0\n",
    "    t1=einsum('ijk,ik->ij',im,p)[:,0]\n",
    "    t2=einsum('ijk,ik->ij',im,p)[:,1]\n",
    "    c1=(t2>=0)&(t2<=1)\n",
    "    c2=t1>=0\n",
    "    t=t1[c1&c2]\n",
    "    p4=p0[None,:]+array(v1)*t.reshape(-1,1)\n",
    "    return p4.tolist()\n",
    "\n",
    "def sec_clean(sec,sec1,r):\n",
    "    sec1=array([p for p in sec1 if len(ibsap(sec,p))%2==1])\n",
    "    p0=sec\n",
    "    p1=sec[1:]+[sec[0]]\n",
    "    sec6=swapaxes(array([p0,p1]),0,1)\n",
    "    p0,p1=array([p0,p1])\n",
    "    v1=p1-p0\n",
    "    v2=sec1[:,None]-p0\n",
    "    v3=sec1[:,None]-p1\n",
    "    u1=v1/linalg.norm(v1,axis=1).reshape(-1,1)\n",
    "    n=1/linalg.norm(v1,axis=1)\n",
    "    u1.shape,v2.shape\n",
    "    d=einsum('jk,ijk->ij',u1,v2)\n",
    "    t=einsum('ij,j->ij',d,n).round(3)\n",
    "    u1.shape,d.shape\n",
    "    n1=einsum('jk,ij->ijk',u1,d)\n",
    "    p1=p0+n1\n",
    "    sec1.shape,p1.shape\n",
    "    n2=sec1[:,None]-p1\n",
    "    n3=sqrt(einsum('ijk,ijk->ij',n2,n2)).round(3)\n",
    "    n4=where((t>=0)&(t<=1),n3,1e5).min(axis=1)\n",
    "    m=sec1[(n4>=abs(r)-.02)&(n4<=abs(r)+.02)].tolist()\n",
    "    return array(m)[cKDTree(m).query(sec)[1]].tolist()\n",
    "\n",
    "\n",
    "def sec_clean1(sec,sec1,r):\n",
    "#     sec1=array([p for p in sec1 if len(ibsap(sec,p))%2==1])\n",
    "    p0=sec\n",
    "    p1=sec[1:]+[sec[0]]\n",
    "    sec6=swapaxes(array([p0,p1]),0,1)\n",
    "    p0,p1=array([p0,p1])\n",
    "    v1=p1-p0\n",
    "    v2=sec1[:,None]-p0\n",
    "    v3=sec1[:,None]-p1\n",
    "    u1=v1/linalg.norm(v1,axis=1).reshape(-1,1)\n",
    "    n=1/linalg.norm(v1,axis=1)\n",
    "    u1.shape,v2.shape\n",
    "    d=einsum('jk,ijk->ij',u1,v2)\n",
    "    t=einsum('ij,j->ij',d,n).round(3)\n",
    "    u1.shape,d.shape\n",
    "    n1=einsum('jk,ij->ijk',u1,d)\n",
    "    p1=p0+n1\n",
    "    sec1.shape,p1.shape\n",
    "    n2=sec1[:,None]-p1\n",
    "    n3=sqrt(einsum('ijk,ijk->ij',n2,n2)).round(3)\n",
    "    n4=where((t>=0)&(t<=1),n3,1e5).min(axis=1)\n",
    "    m=sec1[(n4>=abs(r)-.02)&(n4<=abs(r)+.02)].tolist()\n",
    "    return array(m)[cKDTree(m).query(sec)[1]].tolist()\n",
    "\n",
    "\n",
    "\n",
    "def fillet_2cir(r1,r2,c1,c2,r): # fillet between 2 circles with radius 'r1' and 'r2' and center points 'c1' and 'c2' and 'r' is the radius of the fillet\n",
    "    c1,c2=array([c1,c2])\n",
    "    l1=linalg.norm(c2-c1)\n",
    "    l2=r1+r\n",
    "    l3=r2+r\n",
    "    t=(l1**2+l2**2-l3**2)/(2*l1)\n",
    "    h=sqrt(l2**2-t**2)\n",
    "    v=c2-c1\n",
    "    u=v/linalg.norm(v)\n",
    "    p1=c1+u*t+(u@rm(90))*h\n",
    "    a1=ang((c1-p1)[0],(c1-p1)[1])\n",
    "    a2=ang((c2-p1)[0],(c2-p1)[1])\n",
    "    p2=c1+u*t+u@rm(-90)*h\n",
    "    a3=ang((c2-p2)[0],(c2-p2)[1])\n",
    "    a4=ang((c1-p2)[0],(c1-p2)[1])\n",
    "    a5=ang((p1-c1)[0],(p1-c1)[1])\n",
    "    a6=ang((p2-c1)[0] ,(p2-c1)[1])\n",
    "    a7=ang((p1-c2)[0] ,(p1-c2)[1])\n",
    "    a8=ang((p2-c2)[0] ,(p2-c2)[1])\n",
    "\n",
    "    arc1=arc(r,360+a2 if a2<a1 else a2,a1,p1)\n",
    "    arc2=arc(r,360+a4 if a4<a3 else a4,a3,p2)\n",
    "    arc3=arc(r2,360+a7 if a7<a8 else a7,a8,c2)\n",
    "    arc4=arc(r1,a5,360+a6 if a6<a5 else a6,c1)\n",
    "\n",
    "    return arc2+arc1\n",
    "\n",
    "def filleto_2cir(r1,r2,c1,c2,r): # fillet between 2 circles with radius 'r1' and 'r2' and center points 'c1' and 'c2' and 'r' is the radius of the fillet. This is an open fillet where first or the second fillet can be called based on requirement\n",
    "    c1,c2=array([c1,c2])\n",
    "    l1=linalg.norm(c2-c1)\n",
    "    l2=r1+r\n",
    "    l3=r2+r\n",
    "    t=(l1**2+l2**2-l3**2)/(2*l1)\n",
    "    h=sqrt(l2**2-t**2)\n",
    "    v=c2-c1\n",
    "    u=v/linalg.norm(v)\n",
    "    p1=c1+u*t+(u@rm(90))*h\n",
    "    a1=ang((c1-p1)[0],(c1-p1)[1])\n",
    "    a2=ang((c2-p1)[0],(c2-p1)[1])\n",
    "    p2=c1+u*t+u@rm(-90)*h\n",
    "    a3=ang((c2-p2)[0],(c2-p2)[1])\n",
    "    a4=ang((c1-p2)[0],(c1-p2)[1])\n",
    "    a5=ang((p1-c1)[0],(p1-c1)[1])\n",
    "    a6=ang((p2-c1)[0] ,(p2-c1)[1])\n",
    "    a7=ang((p1-c2)[0] ,(p1-c2)[1])\n",
    "    a8=ang((p2-c2)[0] ,(p2-c2)[1])\n",
    "\n",
    "    arc1=arc(r,360+a2 if a2<a1 else a2,a1,p1)\n",
    "    arc2=arc(r,360+a4 if a4<a3 else a4,a3,p2)\n",
    "    arc3=arc(r2,360+a7 if a7<a8 else a7,a8,c2)\n",
    "    arc4=arc(r1,a5,360+a6 if a6<a5 else a6,c1)\n",
    "\n",
    "    return [arc2,arc1]\n",
    "\n",
    "def tctp(r1,r2,cp1,cp2): # 2 circle tangent points (one side) r1 and r2 are the radius of 2 circles and cp1 and cp2 are the center points\n",
    "    cp1,cp2=array([cp1,cp2])\n",
    "    v1=cp2-cp1,\n",
    "    u1=v1/linalg.norm(v1)\n",
    "    ang1=arcsin((r2-r1)/linalg.norm(cp2-cp1))*180/pi\n",
    "\n",
    "    t1=cp1+u1@rm(90+ang1)*r1\n",
    "    t2=cp2+u1@rm(90+ang1)*r2\n",
    "\n",
    "    t3=cp1+u1@rm(-90-ang1)*r1\n",
    "    t4=cp2+u1@rm(-90-ang1)*r2\n",
    "    return [t1[0].tolist(),t2[0].tolist()]\n",
    "\n",
    "def tctpf(r1,r2,cp1,cp2): #2 circle tangent point full (both the sides)\n",
    "    cp1,cp2=array([cp1,cp2])\n",
    "    v1=cp2-cp1,\n",
    "    u1=v1/linalg.norm(v1)\n",
    "    ang1=arcsin((r2-r1)/linalg.norm(cp2-cp1))*180/pi\n",
    "\n",
    "    t1=cp1+u1@rm(90+ang1)*r1\n",
    "    t2=cp2+u1@rm(90+ang1)*r2\n",
    "\n",
    "    t3=cp1+u1@rm(-90-ang1)*r1\n",
    "    t4=cp2+u1@rm(-90-ang1)*r2\n",
    "    return [t1[0].tolist(),t2[0].tolist(),t4[0].tolist(),t3[0].tolist()]\n",
    "\n",
    "def circle(r,cp=[0,0],s=50): # circle with radius r and center point cp, s is the number of segments in the circle\n",
    "    return array([ [cp[0]+r*cos(i*pi/180),cp[1]+r*sin(i*pi/180)] for i in linspace(0,360,s)][0:-1]).tolist()\n",
    "\n",
    "def circle_c(r,cp=[0,0],s=50):\n",
    "    c=array([ [cp[0]+r*cos(i*pi/180),cp[1]+r*sin(i*pi/180)] for i in linspace(0,360,s)][0:-1]).tolist()\n",
    "    p0,p1=array([c[len(c)-1],c[0]])\n",
    "    v=p1-p0\n",
    "    p=(p0+v*.999).tolist()\n",
    "    return c+[p]\n",
    "\n",
    "def qmr1(s,r,pl):\n",
    "    for i in range(len(s)):\n",
    "        a=[1,0,0] if s[i]=='x' else [0,1,0] if s[i]=='y' else [0,0,1]\n",
    "        b=r[i]\n",
    "        pl=[q(a,p,b) for p in pl]\n",
    "    return pl\n",
    "\n",
    "def qmr2(s,r,pl):\n",
    "    for i in range(len(s)):\n",
    "        a=[1,0,0] if s[i]=='x' else [0,1,0] if s[i]=='y' else [0,0,1]\n",
    "        b=r[i]\n",
    "        pl=[[q(a,p1,b) for p1 in p]for p in pl]\n",
    "    return pl\n",
    "\n",
    "def q_rot(s,pl):\n",
    "    if len(array(pl).shape)==2:\n",
    "        return qmr1([p[0] for p in s],[0 if len(p)==1 else float(p[1:]) for p in s],pl)\n",
    "    else:\n",
    "        return qmr2([p[0] for p in s],[0 if len(p)==1 else float(p[1:]) for p in s],pl)\n",
    "    \n",
    "def l_extrude(sec,h=1,a=0,steps=1):\n",
    "    s=2 if a==0 else steps\n",
    "    return [trns([0,0,h*i if a==0 else h/a*i],q_rot([f\"z{0 if a==0 else i}\"],sec)) for i in linspace(0,1 if a==0 else a,s)]\n",
    "\n",
    "def cylinder(r1=1,r2=1,h=1,cp=[0,0],s=50,r=0,d=0,d1=0,d2=0,center=False):\n",
    "    ra=r if r>0 else d/2 if d>0 else d1/2 if d1>0 else r1\n",
    "    rb=r if r>0 else d/2 if d>0 else d2/2 if d2>0 else r2\n",
    "    sec=circle(ra,cp,s)\n",
    "    \n",
    "    path=pts([[-ra+.1,0],[ra-.1,0],[rb-ra,h],[-rb+.1,0]])\n",
    "    p= trns([0,0,-h/2],prism(sec,path)) if center==True else prism(sec,path)\n",
    "    return p\n",
    "\n",
    "def square(s=0,center=False):\n",
    "    m= s if type(s)==int or type(s)==float else s[0]\n",
    "    n= s if type(s)==int or type(s)==float else s[1]\n",
    "    sec=cr(pts1([[0,0,.001],[m,0,.001],[0,n,.001],[-m,0,.001]]),10)\n",
    "    sec1= [[p[0]-m/2,p[1]-n/2] for p in sec] if center==True else sec\n",
    "    return sec1\n",
    "\n",
    "def rsz3d(prism,rsz):\n",
    "    prism1=array(prism).reshape(-1,3)\n",
    "    max_x=prism1[:,0].max()\n",
    "    max_y=prism1[:,1].max()\n",
    "    max_z=prism1[:,2].max()\n",
    "    min_x=prism1[:,0].min()\n",
    "    min_y=prism1[:,1].min()\n",
    "    min_z=prism1[:,2].min()\n",
    "    avg=prism1.mean(axis=0)\n",
    "    \n",
    "    r_x=rsz[0]/(max_x-min_x)\n",
    "    r_y=rsz[1]/(max_y-min_y)\n",
    "    r_z=rsz[2]/(max_z-min_z)\n",
    "    \n",
    "    rev_prism=[[[avg[0]+r_x*(p[0]-avg[0]),avg[1]+r_y*(p[1]-avg[1]),avg[2]+r_z*(p[2]-avg[2])] for p in prism[i]] \n",
    "               for i in range(len(prism))]\n",
    "    t=((array(bb(rev_prism))-array(bb(prism)))/2).tolist()\n",
    "    return trns(t,rev_prism)\n",
    "\n",
    "def rsz3dc(prism,rsz):\n",
    "    prism1=array(prism).reshape(-1,3)\n",
    "    max_x=prism1[:,0].max()\n",
    "    max_y=prism1[:,1].max()\n",
    "    max_z=prism1[:,2].max()\n",
    "    min_x=prism1[:,0].min()\n",
    "    min_y=prism1[:,1].min()\n",
    "    min_z=prism1[:,2].min()\n",
    "    avg=prism1.mean(axis=0)\n",
    "    \n",
    "    r_x=rsz[0]/(max_x-min_x)\n",
    "    r_y=rsz[1]/(max_y-min_y)\n",
    "    r_z=rsz[2]/(max_z-min_z)\n",
    "    \n",
    "    rev_prism=[[[avg[0]+r_x*(p[0]-avg[0]),avg[1]+r_y*(p[1]-avg[1]),avg[2]+r_z*(p[2]-avg[2])] for p in prism[i]] \n",
    "               for i in range(len(prism))]\n",
    "    return rev_prism\n",
    "\n",
    "\n",
    "def bb(prism):\n",
    "    prism1=array(prism).reshape(-1,3)\n",
    "    max_x=prism1[:,0].max()\n",
    "    max_y=prism1[:,1].max()\n",
    "    max_z=prism1[:,2].max()\n",
    "    min_x=prism1[:,0].min()\n",
    "    min_y=prism1[:,1].min()\n",
    "    min_z=prism1[:,2].min()\n",
    "    return [max_x-min_x,max_y-min_y,max_z-min_z]\n",
    "\n",
    "# def cube(s,center=False):\n",
    "#     m=s if type(s)==int or type(s)==float else s[0]\n",
    "#     n=s if type(s)==int or type(s)==float else s[1]\n",
    "#     o=s if type(s)==int or type(s)==float else s[2]\n",
    "#     path=cr(pts1([[-m/2,0],[m/2,0],[0,o],[-m/2,0]]),1)\n",
    "#     p=trns([-m/2,-n/2,-o/2],rsz3d(prism(square(m),path),[m,n,o])) if center==True else rsz3d(prism(square(m),path),[m,n,o])\n",
    "#     return array(p).tolist()\n",
    "\n",
    "def cube(s,center=False):\n",
    "    if center==False:\n",
    "        return l_extrude(square([s[0],s[1]]),s[2])\n",
    "    elif center==True:\n",
    "        return trns([0,0,-s[2]/2],l_extrude(square([s[0],s[1]],True),s[2]))\n",
    "\n",
    "\n",
    "def sphere(r=0,cp=[0,0,0],s=50):\n",
    "    path=arc(r,-90,90,s=s)\n",
    "    p=[ trns([cp[0],cp[1],p[1]+cp[2]],circle(p[0],s=s)) for p in path]\n",
    "    return array(p).tolist()\n",
    "\n",
    "def rsz2d(sec,rsz):\n",
    "    avg=array(sec).mean(axis=0)\n",
    "    max_x=array(sec)[:,0].max()\n",
    "    min_x=array(sec)[:,0].min()\n",
    "    max_y=array(sec)[:,1].max()\n",
    "    min_y=array(sec)[:,1].min()\n",
    "    r_x=rsz[0]/(max_x-min_x)\n",
    "    r_y=rsz[1]/(max_y-min_y)\n",
    "    s=array([ avg+array([r_x*(sec[i][0]-avg[0]),r_y*(sec[i][1]-avg[1])-((min_y-avg[1])*r_y-(min_y-avg[1]))]) for i in range(len(sec))]).round(4)\n",
    "    return s[sort(unique(s,axis=0,return_index=True)[1])].tolist()\n",
    "    \n",
    "def rsz2dc(sec,rsz):\n",
    "    avg=array(sec).mean(axis=0)\n",
    "    max_x=array(sec)[:,0].max()\n",
    "    min_x=array(sec)[:,0].min()\n",
    "    max_y=array(sec)[:,1].max()\n",
    "    min_y=array(sec)[:,1].min()\n",
    "    r_x=rsz[0]/(max_x-min_x)\n",
    "    r_y=rsz[1]/(max_y-min_y)\n",
    "    s=array([ avg+array([r_x*(sec[i][0]-avg[0]),r_y*(sec[i][1]-avg[1])]) for i in range(len(sec))]).round(4)\n",
    "    return s[sort(unique(s,axis=0,return_index=True)[1])].tolist()\n",
    "\n",
    "def ip(prism,prism1):\n",
    "    pa=prism\n",
    "    pb=prism1\n",
    "    p1=array([[ [[pa[i][j],pa[i][j+1],pa[i+1][j]],[pa[i+1][j+1],pa[i+1][j],pa[i][j+1]]] if j<len(pa[i])-1 \n",
    "     else [[pa[i][j],pa[i][0],pa[i+1][j]],[pa[i+1][0],pa[i+1][j],pa[i][0]]] \n",
    "     for j in range(len(pa[i]))] \n",
    "              for i in range(len(pa)-1)]).reshape(-1,3,3)\n",
    "    p2=array([[[pb[i][j],pb[i+1][j]] for j in range(len(pb[i]))] for i in range(len(pb)-1)]).reshape(-1,2,3)\n",
    "    pm=p1[:,0]\n",
    "    pn=p1[:,1]\n",
    "    po=p1[:,2]\n",
    "    px=p2[:,0]\n",
    "    py=p2[:,1]\n",
    "    v1,v2,v3=py-px,pn-pm,po-pm\n",
    "    t1=einsum('ijk,jk->ij',px[:,None]-pm,cross(v2,v3))/einsum('ik,jk->ij',-v1,cross(v2,v3))\n",
    "    t2=einsum('ijk,ijk->ij',px[:,None]-pm,cross(v3,-v1[:,None]))/einsum('ik,jk->ij',-v1,cross(v2,v3))\n",
    "    t3=einsum('ijk,ijk->ij',px[:,None]-pm,cross(-v1[:,None],v2))/einsum('ik,jk->ij',-v1,cross(v2,v3))\n",
    "    p=px[:,None]+einsum('ik,ij->ijk',v1,t1)\n",
    "    condition=(t1>=0)&(t1<=1)&(t2>=0)&(t2<=1)&(t3>=0)&(t3<=1)&((t2+t3)>=0)&((t2+t3)<=1)\n",
    "    p=p[condition]\n",
    "#     p=p[unique(p,return_index=True)[1]]\n",
    "    return p.tolist()\n",
    "\n",
    "def ipf(prism,prism1,r,s,o=0):\n",
    "    pa=prism\n",
    "    pb=prism1\n",
    "    p1=array([[ [[pa[i][j],pa[i][j+1],pa[i+1][j]],[pa[i+1][j+1],pa[i+1][j],pa[i][j+1]]] if j<len(pa[i])-1 \n",
    "     else [[pa[i][j],pa[i][0],pa[i+1][j]],[pa[i+1][0],pa[i+1][j],pa[i][0]]] \n",
    "     for j in range(len(pa[i])-1)] \n",
    "              for i in range(len(pa)-1)]).reshape(-1,3,3)\n",
    "    p2=array([[[pb[i][j],pb[i+1][j]] for j in range(len(pb[i]))] for i in range(len(pb)-1)]).reshape(-1,2,3)\n",
    "    pm=p1[:,0]\n",
    "    pn=p1[:,1]\n",
    "    po=p1[:,2]\n",
    "    px=p2[:,0]\n",
    "    py=p2[:,1]\n",
    "    v1,v2,v3=py-px,pn-pm,po-pm\n",
    "#     px+v1*t1=pm+v2*t2+v3*t3\n",
    "#     v1*t1-v2*t2-v3*t3=pm-px\n",
    "    u1=v1/(linalg.norm(v1,axis=1).reshape(-1,1)+.0001)\n",
    "    t1=einsum('ijk,jk->ij',px[:,None]-pm,cross(v2,v3))/(einsum('ik,jk->ij',-v1,cross(v2,v3))+.0001)\n",
    "    t2=einsum('ijk,ijk->ij',px[:,None]-pm,cross(v3,-v1[:,None]))/(einsum('ik,jk->ij',-v1,cross(v2,v3))+.0001)\n",
    "    t3=einsum('ijk,ijk->ij',px[:,None]-pm,cross(-v1[:,None],v2))/(einsum('ik,jk->ij',-v1,cross(v2,v3))+.0001)\n",
    "    p=px[:,None]+einsum('ik,ij->ijk',v1,t1)\n",
    "    pq=p+(u1*r)[:,None]\n",
    "    p.shape,pq.shape\n",
    "    condition=(t1>=0)&(t1<=1)&(t2>=0)&(t2<=1)&(t3>=0)&(t3<=1)&((t2+t3)>=0)&((t2+t3)<=1)\n",
    "    p=p[condition].tolist()\n",
    "    pp=p[1:]+[p[0]]\n",
    "    pq=pq[condition].tolist()\n",
    "    v4=array(pp)-array(p)\n",
    "    pnt=array(pq)-array(p)\n",
    "    n=cross(v4,pnt)\n",
    "    n=n/(linalg.norm(n,axis=1).reshape(-1,1)+.0001)*r\n",
    "    pnt=n\n",
    "    cir=[[(p[i]+array(q(v4[i],pnt[i],t))).tolist() for t in linspace(-90,90,10)]for i in range(len(v4))] if o==0 else \\\n",
    "    [[(p[i]+array(q(v4[i],pnt[i],-t))).tolist() for t in linspace(90,270,10)]for i in range(len(v4))] \n",
    "    p2=[[ [cir[i][j],cir[i][j+1]] for j in range(len(cir[i])-1)] for i in range(len(cir))]\n",
    "    p2=array(p2).reshape(-1,2,3)\n",
    "    px=p2[:,0]\n",
    "    py=p2[:,1]\n",
    "    v1,v2,v3=py-px,pn-pm,po-pm\n",
    "    u1=v1/(linalg.norm(v1,axis=1).reshape(-1,1)+.0001)\n",
    "    t1=einsum('ijk,jk->ij',px[:,None]-pm,cross(v2,v3))/(einsum('ik,jk->ij',-v1,cross(v2,v3))+.0001)\n",
    "    t2=einsum('ijk,ijk->ij',px[:,None]-pm,cross(v3,-v1[:,None]))/(einsum('ik,jk->ij',-v1,cross(v2,v3))+.0001)\n",
    "    t3=einsum('ijk,ijk->ij',px[:,None]-pm,cross(-v1[:,None],v2))/(einsum('ik,jk->ij',-v1,cross(v2,v3))+.0001)\n",
    "    m=px[:,None]+einsum('ik,ij->ijk',v1,t1)\n",
    "    condition=(t1>=0)&(t1<=1)&(t2>=0)&(t2<=1)&(t3>=0)&(t3<=1)&((t2+t3)>=0)&((t2+t3)<=1)\n",
    "    m=m[condition]\n",
    "    m=unique(m,axis=0)[:-1]\n",
    "    m=m[cKDTree(m).query(p)[1]].tolist()\n",
    "    p=swapaxes(array([m,p,pq]),0,1)\n",
    "    p=[[fillet_3p_3d(p2,p1,p0,r_3p_3d([p0,p1,p2]),s)]for (p0,p1,p2) in p]\n",
    "    p=array(p).reshape(-1,s+1,3).tolist()\n",
    "    return p+[p[0]]\n",
    "\n",
    "def ipe(prism,prism1,r,s,o):\n",
    "    a=array(cpo(ipf(prism,prism1,r,s,o)))[1:].tolist()\n",
    "    return a\n",
    "\n",
    "\n",
    "def s_int(s): #creates intersection between all the segments of a section\n",
    "    p0=array([array(s)[:,0]]*len(s)).transpose(1,0,2)\n",
    "    p1=array([array(s)[:,1]]*len(s)).transpose(1,0,2)\n",
    "    v1=p1-p0\n",
    "    p2=array([array(s)[:,0]]*len(s))\n",
    "    p3=array([array(s)[:,1]]*len(s))\n",
    "    v2=p3-p2\n",
    "    v1.shape,v2.shape\n",
    "    A=linalg.pinv(array([v1,-v2]).transpose(1,0,2,3).transpose(0,2,1,3).transpose(0,1,3,2))\n",
    "    B=p2-p0\n",
    "    t=einsum('ijkl,ijl->ijk',A,B)[:,:,0].round(4)\n",
    "    u=einsum('ijkl,ijl->ijk',A,B)[:,:,1].round(4)\n",
    "    condition=(t>=0)&(t<=1)&(u>=0)&(u<=1)\n",
    "    d=(p0+einsum('ijk,ij->ijk',v1,t))[condition].tolist()\n",
    "    return d\n",
    "\n",
    "def comb(n,i): # calculates number of possible combinations for \"n\" items with \"i\" selected items\n",
    "    return int(math.factorial(n)/(math.factorial(i)*math.factorial(n-i)))\n",
    "\n",
    "def bezier(p,s=10): # bezier curve defined by points 'p' and number of segments 's'\n",
    "    return array([ array([ comb((len(p)-1),i)*(1-t)**((len(p)-1)-i)*t**i*array(p[i])  for i in range(len(p))]).sum(0) for t in linspace(0,1,s)]).tolist()\n",
    "\n",
    "def arc_3d(v=[0,0,1],r=1,theta1=0,theta2=360,cw=-1,s=50): # 3d arc defined by normal vector 'v', radius 'r1', start angle 'theta1', end angle 'theta2' , clockwise(1) or counter clockwise(-1) and number of segments 's'\n",
    "    theta=0 if v[:2]==[0,0] else ang(v[0],v[1])\n",
    "    v=q([0,0,1],v,-theta)\n",
    "    alpha=ang(v[0],v[2])\n",
    "    arc1=arc(r,theta1,theta2,[0,0],s=s) if cw==-1 else flip(arc(r,theta1,theta2,[0,0],s=s))\n",
    "    arc2=q_rot(['x90','z90'],arc1)\n",
    "    return array(q_rot([f'z{theta}',f'y{-alpha}'],arc2)).tolist()\n",
    "\n",
    "def plane(nv,radius): # plane defined by normal 'nv' and 'radius'\n",
    "    sec1=arc_3d(nv,.01,0,360,-1)\n",
    "    sec2=arc_3d(nv,radius,0,360,-1)\n",
    "    plane=[sec1,sec2]\n",
    "    return plane\n",
    "\n",
    "def l_cir_ip(line,cir): # line circle intersection point\n",
    "    p0,p1=array(line)\n",
    "    p2=array(cir)\n",
    "    p3=array(cir[1:]+[cir[0]])\n",
    "    v1=p1-p0\n",
    "    v2=p3-p2\n",
    "    im=linalg.pinv(array([[v1]*len(v2),-v2]).transpose(1,0,2).transpose(0,2,1))\n",
    "    pnt=p2-p0\n",
    "    t=einsum('ijk,ik->ij',im,pnt)\n",
    "    condition=((t>=0)&(t<=1)).all(1)\n",
    "    ip=p2+v2*t[:,1].reshape(-1,1)\n",
    "    return ip[condition].tolist()\n",
    "\n",
    "def s_pnt(pnt): # starting point for calculating convex hull (bottom left point)\n",
    "    pnt=array(pnt)\n",
    "    c1=pnt[:,1]==pnt[:,1].min()\n",
    "    s1=pnt[c1]\n",
    "    c2=s1[:,0]==s1[:,0].min()\n",
    "    return s1[c2][0].tolist()\n",
    "\n",
    "def n_pnt(pnt,sp,an):\n",
    "    pnt,sp=array(pnt),array(sp)\n",
    "    pnt=pnt[(pnt!=sp).all(1)]\n",
    "    a=pnt-sp\n",
    "    a1=vectorize(ang)(a[:,0],a[:,1])\n",
    "    n_pnt=pnt[a1==a1[a1>=an].min()][0].tolist()\n",
    "    return [n_pnt,a1[a1>=an].min().tolist()]\n",
    "\n",
    "def c_hull(pnt): # convex hull for an array of points\n",
    "    c=[]\n",
    "    np=n_pnt(pnt,s_pnt(pnt),0)\n",
    "    for i in range(len(pnt)):\n",
    "        c.append(np[0])\n",
    "        np=n_pnt(pnt,np[0],np[1])\n",
    "        if np[0]==s_pnt(pnt):\n",
    "            break\n",
    "    return [s_pnt(pnt)]+c\n",
    "\n",
    "def convex(sec): # to check whether a closed section is convex\n",
    "    s=flip(sec) if cw(sec)==1 else sec\n",
    "    return True if offset_points_cw(s,-1)==[] else False\n",
    "\n",
    "def oo_convex(sec,r): #outer offset of a convex section\n",
    "    s=flip(sec) if cw(sec)==1 else sec\n",
    "    return offset_points(sec,r)\n",
    "\n",
    "def cir_p_t(cir,pnt): # circle to point tangent line (point should be outside the circle)\n",
    "    p0=cir\n",
    "    p1=cir[1:]+[cir[0]]\n",
    "    p0,p1=array([p0,p1])\n",
    "    v=p1-p0\n",
    "    a1=vectorize(ang)(v[:,0],v[:,1])\n",
    "    v1=array(pnt)-p0\n",
    "    a2=vectorize(ang)(v1[:,0],v1[:,1])\n",
    "    an=abs(a1-a2).round(4)\n",
    "    a=360/len(cir)/2\n",
    "    cond=abs(a1-a2)<a\n",
    "    an1=abs(a1-a2)[cond].round(4)\n",
    "    return array(cir)[an==an1][0].tolist()\n",
    "\n",
    "def p_cir_t(pnt,cir): # point to circle tangent line (point should be outside the circle)\n",
    "    p0=cir\n",
    "    p1=cir[1:]+[cir[0]]\n",
    "    p0,p1=array([p0,p1])\n",
    "    v=p1-p0\n",
    "    a1=vectorize(ang)(v[:,0],v[:,1])\n",
    "    v1=p0-array(pnt)\n",
    "    a2=vectorize(ang)(v1[:,0],v1[:,1])\n",
    "    an=abs(a1-a2).round(4)\n",
    "    a=360/len(cir)/2\n",
    "    cond=abs(a1-a2)<a\n",
    "    an1=abs(a1-a2)[cond].round(4)\n",
    "    return array(cir)[an==an1][0].tolist()\n",
    "\n",
    "def p_extrude(sec,path): # section extrude through a path\n",
    "    p0=path\n",
    "    p1=p0[1:]+[p0[0]]\n",
    "    p0,p1=array(p0),array(p1)\n",
    "    v=p1-p0\n",
    "    a1=vectorize(ang)(v[:,0],v[:,1])\n",
    "    b=sqrt(v[:,0]**2+v[:,1]**2)\n",
    "    a2=vectorize(ang)(b,v[:,2])\n",
    "    c=[]\n",
    "    for i in range(len(path)-1):\n",
    "        sec1=trns(p0[i],q_rot(['x90','z-90',f'y{-a2[i]}',f'z{a1[i]}'],sec))\n",
    "        sec2=trns(p1[i],q_rot(['x90','z-90',f'y{-a2[i]}',f'z{a1[i]}'],sec))\n",
    "        if i<len(path)-2:\n",
    "            c.append([sec1])\n",
    "        else:\n",
    "            c.append([sec1,sec2])\n",
    "    return flip(concatenate(c).tolist())\n",
    "\n",
    "def p_extrudec(sec,path): # section extrude through a path (closed path)\n",
    "    p0=path\n",
    "    p1=p0[1:]+[p0[0]]\n",
    "    p0,p1=array(p0),array(p1)\n",
    "    v=p1-p0\n",
    "    a1=vectorize(ang)(v[:,0],v[:,1])\n",
    "    b=sqrt(v[:,0]**2+v[:,1]**2)\n",
    "    a2=vectorize(ang)(b,v[:,2])\n",
    "    c=[]\n",
    "    for i in range(len(path)-1):\n",
    "        sec1=trns(p0[i],q_rot(['x90','z-90',f'y{-a2[i]}',f'z{a1[i]}'],sec))\n",
    "        sec2=trns(p1[i],q_rot(['x90','z-90',f'y{-a2[i]}',f'z{a1[i]}'],sec))\n",
    "        if i<len(path)-2:\n",
    "            c.append([sec1])\n",
    "        else:\n",
    "            c.append([sec1,sec2])\n",
    "    a=concatenate(c).tolist()\n",
    "    return flip(a+[a[0]])\n",
    "\n",
    "def v_sec_extrude(sec,path,o): #variable section extrude through a given path\n",
    "    sec=[offset(sec,i) for i in linspace(0,o,len(path))]\n",
    "    p0=path\n",
    "    p1=p0[1:]+[p0[0]]\n",
    "    p0,p1=array(p0),array(p1)\n",
    "    v=p1-p0\n",
    "    a1=vectorize(ang)(v[:,0],v[:,1])\n",
    "    b=sqrt(v[:,0]**2+v[:,1]**2)\n",
    "    a2=vectorize(ang)(b,v[:,2])\n",
    "    c=[]\n",
    "    for i in range(len(path)-1):\n",
    "        sec1=trns(p0[i],q_rot(['x90','z-90',f'y{-a2[i]}',f'z{a1[i]}'],sec[i]))\n",
    "        sec2=trns(p1[i],q_rot(['x90','z-90',f'y{-a2[i]}',f'z{a1[i]}'],sec[i]))\n",
    "        if i<len(path)-2:\n",
    "            c.append([sec1])\n",
    "        else:\n",
    "            c.append([sec1,sec2])\n",
    "    return concatenate(c).tolist()\n",
    "\n",
    "def t_cir_tarc(r1,r2,cp1,cp2,r,s=50): #two circle tangent arc\n",
    "    cp1,cp2=array([cp1,cp2])\n",
    "    l1=linalg.norm(cp2-cp1)\n",
    "    l2=r-r1\n",
    "    l3=r-r2\n",
    "    x=(l2**2-l3**2+l1**2)/(2*l1)\n",
    "    h=sqrt(l2**2-x**2)\n",
    "    v1=cp2-cp1\n",
    "    u1=v1/linalg.norm(v1)\n",
    "    p0=cp1+u1*x\n",
    "    cp3=p0-(u1@rm(90))*h\n",
    "    v2=cp2-cp3\n",
    "    u2=v2/linalg.norm(v2)\n",
    "    v3=cp1-cp3\n",
    "    u3=v3/linalg.norm(v3)\n",
    "    ang1=ang(u2[0],u2[1])\n",
    "    ang2=ang(u3[0],u3[1])\n",
    "    return array(arc(r,ang1,ang2,cp3,s)).tolist()\n",
    "\n",
    "def tcct(r1,r2,cp1,cp2,cw=-1): # two circle cross tangent\n",
    "    v1=[1,1]\n",
    "    v2=[-r2,r1]\n",
    "    cp1,cp2=array([cp1,cp2])\n",
    "    d=linalg.norm(cp2-cp1)\n",
    "    d1=(linalg.inv(array([v1,v2]).T)@array([d,0]))[0]\n",
    "    d2=(linalg.inv(array([v1,v2]).T)@array([d,0]))[1]\n",
    "    a=arcsin(r1/d1)*180/pi\n",
    "    v3=cp2-cp1\n",
    "    u3=v3/linalg.norm(v3)\n",
    "    b=arccos(u3@array([1,0]))*180/pi\n",
    "    if cw==-1:\n",
    "        if v3[0]>0 and v3[1]<=0:\n",
    "            theta1=270+a-b\n",
    "            theta2=90+a-b\n",
    "        elif v3[0]>=0 and v3[1]>0:\n",
    "            theta1=270+a+b\n",
    "            theta2=90+a+b\n",
    "        elif v3[0]<0 and v3[1]>=0:\n",
    "            theta1=270+a+b\n",
    "            theta2=90+a+b\n",
    "        else:\n",
    "            theta1=270+a-b\n",
    "            theta2=90+a-b\n",
    "    else:\n",
    "        if v3[0]>0 and v3[1]<=0:\n",
    "            theta2=270-a-b\n",
    "            theta1=90-a-b\n",
    "        elif v3[0]>=0 and v3[1]>0:\n",
    "            theta2=270-a+b\n",
    "            theta1=90-a+b\n",
    "        elif v3[0]<0 and v3[1]>=0:\n",
    "            theta2=270-a+b\n",
    "            theta1=90-a+b\n",
    "        else:\n",
    "            theta2=270-a-b\n",
    "            theta1=90-a-b\n",
    "        \n",
    "    p0=(cp1+array([r1*cos(theta1*pi/180),r1*sin(theta1*pi/180)])).tolist()\n",
    "    p1=(cp2+array([r2*cos(theta2*pi/180),r2*sin(theta2*pi/180)])).tolist()\n",
    "    return [p0,p1]\n",
    "\n",
    "def arc_3p(p1,p2,p3,s=30):\n",
    "    p1,p2,p3=array([p1,p2,p3])\n",
    "    p4=p1+(p2-p1)/2\n",
    "    p5=p2+(p3-p2)/2\n",
    "    v1=p2-p4\n",
    "    u1=v1/linalg.norm(v1)\n",
    "    v2=p3-p5\n",
    "    u2=v2/linalg.norm(v2)\n",
    "    p6=p4+u1@rm(90)\n",
    "    p7=p5+u2@rm(90)\n",
    "    cp=i_p2d([p4,p6],[p5,p7])\n",
    "    r=linalg.norm(p1-cp)\n",
    "    v3=p1-cp\n",
    "    v4=p2-cp\n",
    "    v5=p3-cp\n",
    "    a1=ang(v3[0],v3[1])\n",
    "    a2=ang(v4[0],v4[1])\n",
    "    a3=ang(v5[0],v5[1])\n",
    "    a4=(a3+360 if a3<a1 else a3) if cw([p1,p2,p3])==-1 else (a3 if a3<a1 else a3-360)\n",
    "    return arc(r,a1,a4,cp,s)\n",
    "\n",
    "def cir_3p(p1,p2,p3,s=30):\n",
    "    p1,p2,p3=array([p1,p2,p3])\n",
    "    p4=p1+(p2-p1)/2\n",
    "    p5=p2+(p3-p2)/2\n",
    "    v1=p2-p4\n",
    "    u1=v1/linalg.norm(v1)\n",
    "    v2=p3-p5\n",
    "    u2=v2/linalg.norm(v2)\n",
    "    p6=p4+u1@rm(90)\n",
    "    p7=p5+u2@rm(90)\n",
    "    cp=i_p2d([p4,p6],[p5,p7])\n",
    "    r=linalg.norm(p1-cp)\n",
    "#     v3=p1-cp\n",
    "#     v4=p2-cp\n",
    "#     v5=p3-cp\n",
    "#     a1=ang(v3[0],v3[1])\n",
    "#     a2=ang(v4[0],v4[1])\n",
    "#     a3=ang(v5[0],v5[1])\n",
    "#     a4=(a3+360 if a3<a1 else a3) if cw([p1,p2,p3])==-1 else (a3 if a3<a1 else a3-360)\n",
    "    return circle(r,cp,s)\n",
    "\n",
    "def cp_3p(p1,p2,p3):\n",
    "    p1,p2,p3=array([p1,p2,p3])\n",
    "    p4=p1+(p2-p1)/2\n",
    "    p5=p2+(p3-p2)/2\n",
    "    v1=p2-p4\n",
    "    u1=v1/linalg.norm(v1)\n",
    "    v2=p3-p5\n",
    "    u2=v2/linalg.norm(v2)\n",
    "    p6=p4+u1@rm(90)\n",
    "    p7=p5+u2@rm(90)\n",
    "    cp=i_p2d([p4,p6],[p5,p7])\n",
    "    return array(cp).tolist()\n",
    "\n",
    "\n",
    "\n",
    "def ip_surf(surf2,surf1):\n",
    "    i,j,_=array(surf2).shape\n",
    "    a=surf2\n",
    "    b=surf1\n",
    "    p1=array([[[[a[i][j],a[i+1][j],a[i][j+1]],[a[i+1][j+1],a[i][j+1],a[i+1][j]]] \n",
    "            for j in range(j-1)]  for i in range(i-1)]).reshape(-1,3,3)\n",
    "    p2=array([[[b[i][j],b[i+1][j]] for j in range(len(b[i]))] for i in range(len(b)-1)]).reshape(-1,2,3)\n",
    "    pm=p1[:,0]\n",
    "    pn=p1[:,1]\n",
    "    po=p1[:,2]\n",
    "    px=p2[:,0]\n",
    "    py=p2[:,1]\n",
    "    v1,v2,v3=py-px,pn-pm,po-pm\n",
    "    t1=einsum('ijk,jk->ij',px[:,None]-pm,cross(v2,v3))/(einsum('ik,jk->ij',-v1,cross(v2,v3))+.00001)\n",
    "    t2=einsum('ijk,ijk->ij',px[:,None]-pm,cross(v3,-v1[:,None]))/(einsum('ik,jk->ij',-v1,cross(v2,v3))+.00001)\n",
    "    t3=einsum('ijk,ijk->ij',px[:,None]-pm,cross(-v1[:,None],v2))/(einsum('ik,jk->ij',-v1,cross(v2,v3))+.00001)\n",
    "    p=px[:,None]+einsum('ik,ij->ijk',v1,t1)\n",
    "    condition=(t1>=0)&(t1<=1)&(t2>=0)&(t2<=1)&(t3>=0)&(t3<=1)&((t2+t3)>=0)&((t2+t3)<=1)\n",
    "    return p[condition].tolist()\n",
    "\n",
    "def perp(sec,point,radius):\n",
    "    sec=array(seg(sec))\n",
    "    p0=sec[:,0]\n",
    "    p1=sec[:,1]\n",
    "    v1=p1-p0\n",
    "    u1=v1/(linalg.norm(v1,axis=1).reshape(-1,1)+.00001)\n",
    "    v2=array(point)-p0\n",
    "    v1norm=linalg.norm(v1,axis=1)\n",
    "    v2norm=linalg.norm(v2,axis=1)\n",
    "    v2cost=einsum('ij,ij->i',u1,v2)\n",
    "    cond1=v2cost>=0\n",
    "    cond2=v2cost<=v1norm\n",
    "    d=sqrt(v2norm**2-v2cost**2)\n",
    "    d=min(d[(cond1)&(cond2)]).round(4)\n",
    "    cond3=d==round(abs(radius),3)\n",
    "    return point if cond3 else []\n",
    "\n",
    "def perp_point(line,point,distance):\n",
    "    p0=line[0]\n",
    "    p1=line[1]\n",
    "    p0,p1=array([p0,p1])\n",
    "    v1=p1-p0\n",
    "    u1=v1/(linalg.norm(v1)+.00001)\n",
    "    v2=array(point)-p0\n",
    "    v1norm=linalg.norm(v1)\n",
    "    v2norm=linalg.norm(v2)\n",
    "    v2cost=u1@v2\n",
    "    cond1=v2cost>=0\n",
    "    cond2=v2cost<=v1norm\n",
    "    d=sqrt(v2norm**2-v2cost**2)\n",
    "    cond3=d<=distance\n",
    "    return point if cond1 & cond2 & cond3  else []\n",
    "\n",
    "def perp_dist(line,point):\n",
    "    p0=line[0]\n",
    "    p1=line[1]\n",
    "    p0,p1=array([p0,p1])\n",
    "    v1=p1-p0\n",
    "    u1=v1/(linalg.norm(v1)+.00001)\n",
    "    v2=array(point)-p0\n",
    "    v1norm=linalg.norm(v1)\n",
    "    v2norm=linalg.norm(v2)\n",
    "    v2cost=u1@v2\n",
    "    d=sqrt(v2norm**2-v2cost**2)\n",
    "    return d\n",
    "\n",
    "\n",
    "def pies(sec,pnt):\n",
    "    sec1=array([p for p in pnt if len(ibsap(sec,p))%2==1])\n",
    "    return sec1.tolist()\n",
    "\n",
    "def sq(d,cp=[0,0]):\n",
    "    cp=array(cp)-d/2\n",
    "    cp=[cp[0],cp[1],0]\n",
    "    return c3t2(trns(cp,[[0,0],[d,0],[d,d],[0,d]]))\n",
    "\n",
    "def near_points(points,s_p,n):\n",
    "    l=array([ linalg.norm(array(p)-array(s_p)) for p in points])\n",
    "    l1=sort(l)[0:n+1]\n",
    "    index=array([[i for i in range(len(l)) if p==l[i]]for p in l1]).reshape(-1)\n",
    "    p1=array(points)[index].tolist()\n",
    "    return p1[1:]\n",
    "\n",
    "def next_point(points,s_p):\n",
    "    a1=[270+(360-ang((array(p)-array(s_p))[0],(array(p)-array(s_p))[1]))\n",
    "        if ang((array(p)-array(s_p))[0],(array(p)-array(s_p))[1])>270 else\n",
    "        270-ang((array(p)-array(s_p))[0],(array(p)-array(s_p))[1])\n",
    "        for p in points]\n",
    "    n_p=array(points)[a1==max(a1)][0].tolist()\n",
    "    return n_p\n",
    "\n",
    "def exclude_points(points,pnts):\n",
    "    return [p for p in points if p not in pnts]\n",
    "\n",
    "def i_p2dw(l1,l):\n",
    "    p0,p1=array(l1)\n",
    "    p2,p3=array(l)\n",
    "    v1=p1-p0\n",
    "    v2=p3-p2\n",
    "#                     p0+v1*t1=p2+v2*t2\n",
    "#                     v1*t1-v2*t2=p2-p0\n",
    "    im=linalg.inv(array([v1,-v2]).transpose(1,0)\n",
    "                  +array([[.000001,.000002],[.000002,.000003]]))\n",
    "    t=(im@(p2-p0))[0]\n",
    "    u=(im@(p2-p0))[1]\n",
    "    return  (p0+v1*t).tolist() if (0<t<1)& (0<u<1) else []\n",
    "\n",
    "\n",
    "def pies1(s8,s4):\n",
    "    p0=array(s4)\n",
    "    p2=s8\n",
    "    p3=s8[1:]+[s8[0]]\n",
    "    p2,p3=array([p2,p3])\n",
    "    v1=array([[[1,0]]*len(p2)]*len(p0))\n",
    "    v2=array([((p3-p2)+[0,.00001]).tolist()]*len(p0))\n",
    "    # im=linalg.pinv(array([[v1]*len(v2),-v2]).transpose(1,0,2).transpose(0,2,1))\n",
    "    # im=array([im.tolist()]*len(p0))\n",
    "    p=p2-p0[:,None]\n",
    "    # t=einsum('ijkl,ijl->ijl',im,p)\n",
    "    # s10=[p0[i].tolist() for i in range(len(p0)) if \\\n",
    "    #     t[i][(t[i][:,0]>=0)&(t[i][:,1]>=0)&(t[i][:,1]<=1)].shape[0]%2 \\\n",
    "    #  ==1]\n",
    "\n",
    "    im=linalg.pinv(array([v1,-v2]).transpose(1,0,2,3).transpose(0,2,1,3))\n",
    "    im.shape,p.shape\n",
    "    t=einsum('ijkl,ijk->ijl',im,p)\n",
    "    s10=[p0[i].tolist() for i in range(len(p0)) if \\\n",
    "        t[i][(t[i][:,0]>=0)&(t[i][:,1]>=0)&(t[i][:,1]<=1)].shape[0]%2 \\\n",
    "     ==1]\n",
    "    return s10\n",
    "\n",
    "def rsec(line,radius):\n",
    "    p0=line[0]\n",
    "    p1=line[1]\n",
    "    p0,p1=array([p0,p1])\n",
    "    v=p1-p0\n",
    "    a=ang(v[0],v[1])\n",
    "    return arc(radius,a+90,a+270,p0,int(round(10+log10(radius+1)**6,0)))+arc(radius,a-90,a+90,p1,int(round(10+log10(radius+1)**6,0)))\n",
    "\n",
    "\n",
    "\n",
    "def cleaning_seg(sec):\n",
    "    r=-max_r(sec)-1\n",
    "    s=seg(sec)\n",
    "    s1=offset_points(sec,r)\n",
    "    s2=seg(s1)\n",
    "    u=array([(array(p[1])-array(p[0]))/linalg.norm(array(p[1])-array(p[0])) for p in s])\n",
    "    u1=array([(array(p[1])-array(p[0]))/linalg.norm(array(p[1])-array(p[0])) for p in s2])\n",
    "    s3=array(s)[linalg.norm(u-u1,axis=1)<1].tolist()\n",
    "    return s3\n",
    "\n",
    "def cleaning_sec_inner(sec,r):\n",
    "    s=cleaning_seg(sec)\n",
    "    s1=[rsec(p,abs(r)-.01) for p in s]\n",
    "    return s1\n",
    "\n",
    "def cleaning_sec_outer(sec,r):\n",
    "    s=cleaning_seg(sec)\n",
    "    s1=[rsec(p,abs(r)-.1) for p in s]\n",
    "    return s1\n",
    "\n",
    "def inner_offset(sec,r):\n",
    "    sec=flip(sec) if cw(sec)==1 else sec\n",
    "    s=offset_points(sec,r)\n",
    "    if s_intv1(seg(s))!=[]:\n",
    "        s1=unique(s_intv(seg(s)),axis=0).tolist()\n",
    "        for p in cleaning_sec_inner(sec,r):\n",
    "            s2=pies1(p,s1)\n",
    "            s1=exclude_points(s1,s2)\n",
    "        s1=array(s1)[cKDTree(s1).query(sec)[1]].tolist()\n",
    "        return s1\n",
    "    else:\n",
    "        return s\n",
    "\n",
    "def out_offset(sec,r):\n",
    "    sec=flip(sec) if cw(sec)==1 else sec\n",
    "    s=offset_points(sec,r)\n",
    "    if s_intv1(seg(s))!=[]:\n",
    "        s1=unique(s_intv(seg(s)),axis=0).tolist()\n",
    "        for p in cleaning_sec_outer(sec,r):\n",
    "            s2=pies1(p,s1)\n",
    "            s1=exclude_points(s1,s2)\n",
    "        s1=array(s1)[cKDTree(s1).query(sec)[1]].tolist()\n",
    "        return s1\n",
    "    else:\n",
    "        return s\n",
    "\n",
    "\n",
    "def swp(bead2):\n",
    "    n1=arange(len(bead2[0])).tolist()\n",
    "    n2=array([[[[(j+1)+i*len(bead2[0]),j+i*len(bead2[0]),j+(i+1)*len(bead2[0])],[(j+1)+i*len(bead2[0]),j+(i+1)*len(bead2[0]),(j+1)+(i+1)*len(bead2[0])]] \\\n",
    "             if j<len(bead2[0])-1 else \\\n",
    "             [[0+i*len(bead2[0]),j+i*len(bead2[0]),j+(i+1)*len(bead2[0])],[0+i*len(bead2[0]),j+(i+1)*len(bead2[0]),0+(i+1)*len(bead2[0])]] \\\n",
    "                 for j in range(len(bead2[0]))] for i in range(len(bead2)-1)]).reshape(-1,3).tolist()\n",
    "    n3=(array(flip(arange(len(bead2[0]))))+(len(bead2)-1)*len(bead2[0])).tolist()\n",
    "    n=[n1]+n2+[n3]\n",
    "    pnt=array(bead2).reshape(-1,3).round(4).tolist()\n",
    "    return f'polyhedron({pnt},{n},convexity=10);'\n",
    "\n",
    "def swp_c(bead2):\n",
    "    n1=arange(len(bead2[0])).tolist()\n",
    "    n2=array([[[[(j+1)+i*len(bead2[0]),j+i*len(bead2[0]),j+(i+1)*len(bead2[0])],[(j+1)+i*len(bead2[0]),j+(i+1)*len(bead2[0]),(j+1)+(i+1)*len(bead2[0])]] \\\n",
    "             if j<len(bead2[0])-1 else \\\n",
    "             [[0+i*len(bead2[0]),j+i*len(bead2[0]),j+(i+1)*len(bead2[0])],[0+i*len(bead2[0]),j+(i+1)*len(bead2[0]),0+(i+1)*len(bead2[0])]] \\\n",
    "                 for j in range(len(bead2[0]))] for i in range(len(bead2)-1)]).reshape(-1,3).tolist()\n",
    "    n3=(array(flip(arange(len(bead2[0]))))+(len(bead2)-1)*len(bead2[0])).tolist()\n",
    "    n=[n1]+n2+[n3]\n",
    "    pnt=array(bead2).reshape(-1,3).round(4).tolist()\n",
    "    return f'polyhedron({pnt},{n2},convexity=10);'\n",
    "\n",
    "def resurf(surf,f):\n",
    "    base=array(c3t2(surf)).reshape(-1,2).tolist()\n",
    "    c=[]\n",
    "    for i in range(len(surf)):\n",
    "        if len(base)<=2:\n",
    "            break\n",
    "        else:\n",
    "            base1=concave_hull(base,f)\n",
    "            base=exclude_points(base,base1)\n",
    "            c.append(base1)\n",
    "    base=concave_hull(array(c3t2(surf)).reshape(-1,2).tolist(),f)\n",
    "    c=[array(p)[cKDTree(p).query(base)[1]].tolist() for p in c]\n",
    "    base=array(c3t2(surf)).reshape(-1,2).tolist()\n",
    "    surf=array(surf).reshape(-1,3)\n",
    "    c= [surf[cKDTree(base).query(p)[1]].tolist() for p in c]\n",
    "    return c\n",
    "\n",
    "def surf_extrudef(surf,t=-.05):\n",
    "    s=cpo(surf)\n",
    "    s1=trns([0,0,t],[flip(p) for p in s])\n",
    "    s2=array([s,s1]).transpose(1,0,2,3)\n",
    "    \n",
    "    i,j,k,l=s2.shape\n",
    "    s2=s2.reshape(i,j*k,l).tolist()\n",
    "    return s2 if t<0 else flip(s2)\n",
    "\n",
    "\n",
    "\n",
    "def swp_prism_h(prism_big,prism_small):\n",
    "    p1=prism_big\n",
    "    p2=flip(prism_small)\n",
    "    p3=p1+p2+[p1[0]]\n",
    "    return p3\n",
    "    \n",
    "def pmdp(line,pnts): #perpendicular minimum distance point\n",
    "    if pnts==[]:\n",
    "        return line\n",
    "    else:\n",
    "        a=[perp_dist(line,p) for p in pnts]\n",
    "        b=array(pnts)[min(a)==array(a)][0].tolist()\n",
    "        return [line[0],b,line[1]]\n",
    "    \n",
    "\n",
    "\n",
    "def surf_base(surf,h=0):\n",
    "    s=cpo(surf)\n",
    "    s1=trns([0,0,h],c2t3(c3t2([flip(p) for p in s])))\n",
    "    s2=array([s,s1]).transpose(1,0,2,3)\n",
    "    \n",
    "    i,j,k,l=s2.shape\n",
    "    s2=s2.reshape(i,j*k,l).tolist()\n",
    "    return s2 if t<0 else flip(s2)\n",
    "\n",
    "def cr_3d(p,s=5): # Corner radius 3d where 'p' are the list of points (turtle movement) and 's' is number of segments for each arc\n",
    "    pnts=array(p)[:,0:3]\n",
    "    pnts=pnts.cumsum(0)\n",
    "\n",
    "    rds=array(p)[:,3]\n",
    "    c=[]\n",
    "    for i in range(len(pnts)):\n",
    "        if i==0:\n",
    "            p0=pnts[len(pnts)-1]\n",
    "            p1=pnts[i]\n",
    "            p2=pnts[i+1]\n",
    "        elif i<len(pnts)-1:\n",
    "            p0=pnts[i-1]\n",
    "            p1=pnts[i]\n",
    "            p2=pnts[i+1]\n",
    "        else:\n",
    "            p0=pnts[i-1]\n",
    "            p1=pnts[i]\n",
    "            p2=pnts[0]\n",
    "        c.append(fillet_3p_3d(p0,p1,p2,rds[i],s)[1:])\n",
    "    c=array(c).reshape(-1,3).tolist()\n",
    "    return c\n",
    "\n",
    "def p_exc(sec,path,option=0):\n",
    "    p=array(path)\n",
    "    c,d,e=[],[],[]\n",
    "    a2,j=0,0\n",
    "    for i in range(len(p)):\n",
    "        i_plus=i+1 if i<len(p)-1 else 0\n",
    "        p0=p[i]\n",
    "        p1=p[i_plus]\n",
    "        v1=uv(p1-p0)\n",
    "        vz=[0,0,1]\n",
    "        v2=cross(vz,v1).tolist()\n",
    "        theta=0 if v2==[0,0,0] else arccos(array(vz)@array(v1))*180/pi\n",
    "        a=0 if theta==0 else ang(v2[0],v2[1])\n",
    "        a1=0 if v1[2]==0 else ang(v1[0],v1[1]) if v1[2]<0 else -ang(v1[0],v1[1])\n",
    "        d.append(a1)\n",
    "        if i>0:\n",
    "            j= j+1 if abs(a1-d[i-1])>179 else j\n",
    "            e.append(j)\n",
    "            a2=a1+j*180\n",
    "        if option==0:\n",
    "            sec1=trns(path[i],[q(v2,q([0,0,1],p,a+a2),theta) for p in sec])\n",
    "        else:\n",
    "            sec1=trns(path[i],[q(v2,q([0,0,1],p,a+a1),theta) for p in sec])\n",
    "            \n",
    "        c.append(sec1)\n",
    "    c=c+[c[0]]\n",
    "    return c\n",
    "\n",
    "def p_ex(sec,path,option=0):\n",
    "    p=array(path)\n",
    "    c,d,e=[],[],[]\n",
    "    a2,j=0,0\n",
    "    for i in range(len(p)-1):\n",
    "        i_plus=i+1 if i<len(p)-1 else 0\n",
    "        p0=p[i]\n",
    "        p1=p[i_plus]\n",
    "        v1=uv(p1-p0)\n",
    "        vz=[0,0,1]\n",
    "        v2=cross(vz,v1).tolist()\n",
    "        theta=0 if v2==[0,0,0] else arccos(array(vz)@array(v1))*180/pi\n",
    "        a=0 if theta==0 else ang(v2[0],v2[1])\n",
    "        a1=0 if v1[2]==0 else ang(v1[0],v1[1]) if v1[2]<0 else -ang(v1[0],v1[1])\n",
    "        d.append(a1)\n",
    "        if i>0:\n",
    "            j= j+1 if abs(abs(a1)-d[i-1])>100 else j\n",
    "            e.append(j)\n",
    "            a2=a1+j*180\n",
    "        if option==0:\n",
    "            sec1=trns(path[i],[q(v2,q([0,0,1],p,a+a2),theta) for p in sec])\n",
    "        else:\n",
    "            sec1=trns(path[i],[q(v2,q([0,0,1],p,a+a1),theta) for p in sec])\n",
    "            \n",
    "        c.append(sec1)\n",
    "    c=c+[trns(array(path[-1])-array(path[-2]),c[-1])]\n",
    "    return c\n",
    "\n",
    "def helix(radius=10,pitch=10, number_of_coils=1, step_angle=1):\n",
    "    return [[radius*cos(i*pi/180),radius*sin(i*pi/180),i/360*pitch] for i in arange(0,360*number_of_coils,step_angle)]\n",
    "\n",
    "def surf_offset(surf,o):\n",
    "    c=[]\n",
    "    for i in range(len(surf)):\n",
    "        for j in range(len(surf[0])):\n",
    "            j_plus=j+1 if j<len(surf[0])-1 else 0\n",
    "            p0=surf[i][j]\n",
    "            p1=surf[i][j_plus] if i<len(surf)-1 else surf[i-1][j]\n",
    "            p2=surf[i+1][j] if i<len(surf)-1 else surf[i][j_plus]\n",
    "            p0,p1,p2=array([p0,p1,p2])\n",
    "            v1=p1-p0\n",
    "            v2=p2-p0\n",
    "            p=p0+array(uv(cross(v1,v2)))*o\n",
    "            c.append(p.tolist())\n",
    "    l,m,n=array(surf).shape\n",
    "    return array(c).reshape(l,m,n).tolist()\n",
    "\n",
    "def path_to_vectors(path):\n",
    "    c=[]\n",
    "    for i in range(len(path)):\n",
    "        i_plus=i+1 if i<len(path)-1 else 0\n",
    "        p0=path[i]\n",
    "        p1=path[i_plus]\n",
    "        p0,p1=array([p0,p1])\n",
    "        v1=p1-p0\n",
    "        c.append(v1.round(4).tolist())\n",
    "    return vector_correct(c)\n",
    "\n",
    "def vector_correct(c):\n",
    "    for i in range(len(c)):\n",
    "        if i>0 and i<len(c)-1:\n",
    "            if c[i][0]==0 and abs(c[i-1][0])>0 and abs(c[i+1][0])>0:\n",
    "                c[i][0]=.001\n",
    "            if c[i][1]==0 and abs(c[i-1][1])>0 and abs(c[i+1][1])>0:\n",
    "                c[i][1]=.001\n",
    "            if c[i][2]==0 and abs(c[i-1][2])>0 and abs(c[i+1][2])>0:\n",
    "                c[i][2]=.001\n",
    "    path=c\n",
    "    return path\n",
    "\n",
    "def concave_hull(pnts,x=1,loops=10):#x is sensitivity 1 is max and 100 is almost convex hull, loops can be any number less than the number of points\n",
    "    c=c_hull(pnts)\n",
    "    for j in range(loops):\n",
    "        c1=seg(c)\n",
    "        pnts1=exclude_points(pnts,c)\n",
    "        c2=[]\n",
    "        for i in range(len(c1)):\n",
    "            p0,p1=array(c1[i])\n",
    "            v1=p1-p0\n",
    "            u1=array(uv(v1))\n",
    "            v1norm=linalg.norm(v1)\n",
    "            pnts2=[p for p in array(pnts1) if ((u1@(p-p0))>=0)&((u1@(p-p0))<=v1norm) & (abs(cross(v1,p-p0))/v1norm <= v1norm/x) ]\n",
    "            if pnts2!=[]:\n",
    "                lengths=[cross(v1,(p-p0))/v1norm for p in array(pnts2)]\n",
    "                pnt=array(pnts2)[lengths==min(lengths)][0]\n",
    "                pnts1=exclude_points(pnts1,pnt)\n",
    "                c2.append([p0.tolist(),pnt.tolist(),p1.tolist()])\n",
    "            else:\n",
    "                c2.append(c1[i])\n",
    "\n",
    "\n",
    "        c3=remove_extra_points(concatenate(c2).tolist())\n",
    "        n=s_intv1(seg(c3))\n",
    "        if n!=[]:\n",
    "            d=[[p[1] for p1 in array(n) if (array(uv(p1-p[0])).round(4)==array(uv(p[1]-p[0])).round(4)).all() ] for p in array(seg(c3))]\n",
    "            d=concatenate([p for p in d if p!=[]]).tolist()\n",
    "            c3=exclude_points(c3,d)\n",
    "        c=c3\n",
    "    while s_intv1(seg(c))!=[]:\n",
    "        n=s_intv1(seg(c3))\n",
    "        if n==[]:\n",
    "            break\n",
    "        else:\n",
    "            d=[[p[1] for p1 in array(n) if (array(uv(p1-p[0])).round(4)==array(uv(p[1]-p[0])).round(4)).all() ] for p in array(seg(c3))]\n",
    "            d=concatenate([p for p in d if p!=[]]).tolist()\n",
    "            c3=exclude_points(c3,d)\n",
    "        \n",
    "        \n",
    "    return c3\n",
    "\n",
    "# def ipfillet(p,p1,r=1,s=5,o=0):\n",
    "#     pa=[[[[p[i][j],p[i][j+1],p[i+1][j]],[p[i+1][j+1],p[i+1][j],p[i][j+1]]] if j<len(p[0])-1 else \\\n",
    "#          [[p[i][j],p[i][0],p[i+1][j]],[p[i+1][0],p[i+1][j],p[i][0]]] \\\n",
    "#          for j in range(len(p[0]))] for i in range(len(p)-1)]\n",
    "#     pa=array(pa).reshape(-1,3,3)\n",
    "    \n",
    "#     pb=[[[p1[i][j],p1[i+1][j]] for j in range(len(p1[0]))] for i in range(len(p1)-1)]\n",
    "#     pb=array(pb).reshape(-1,2,3)\n",
    "#     ip=[]\n",
    "#     for a in pb:\n",
    "#         for b in pa:\n",
    "#             v1=a[1]-a[0]\n",
    "#             v2=b[1]-b[0]\n",
    "#             v3=b[2]-b[0]\n",
    "#             c=array([v1,-v2,-v3]).transpose(1,0)\n",
    "#             d=b[0]-a[0]\n",
    "#             t=linalg.pinv(c)@d\n",
    "#             if (0<=t[0]<=1)&(0<=t[1]<=1)&(0<=t[2]<=1)&(0<=(t[1]+t[2])<=1):\n",
    "#                 pnt=a[0]+v1*t[0]\n",
    "#                 u1=array(uv(v1))\n",
    "#                 pnt1=pnt+u1*r\n",
    "#                 ip.append([pnt,pnt1])\n",
    "    \n",
    "#     if o==0:\n",
    "#         e=[ip[i][0]+array(uv(cross(ip[i+1][0]-ip[i][0],ip[i][1]-ip[i][0])))*r if i<len(ip)-1 else \\\n",
    "#            ip[i][0]+array(uv(cross(ip[0][0]-ip[i][0],ip[i][1]-ip[i][0])))*r \\\n",
    "#            for i in range(len(ip))]\n",
    "#     elif o==1:\n",
    "#         e=[ip[i][0]+array(uv(cross(ip[i][1]-ip[i][0],ip[i+1][0]-ip[i][0])))*r if i<len(ip)-1 else \\\n",
    "#            ip[i][0]+array(uv(cross(ip[i][1]-ip[i][0],ip[0][0]-ip[i][0])))*r \\\n",
    "#            for i in range(len(ip))]\n",
    "        \n",
    "#     f=[[ip[i][0]+q(ip[i+1][0]-ip[i][0],e[i]-ip[i][0], theta) if i<len(ip)-1 else \\\n",
    "#         ip[i][0]+q(ip[0][0]-ip[i][0],e[i]-ip[i][0] , theta) \\\n",
    "#         for theta in linspace(0,90,3)] for i in range(len(ip))]\n",
    "#     f=array([seg(array(p).tolist())[0:-1] for p in f]).reshape(-1,2,3)\n",
    "#     ip1=[]\n",
    "    \n",
    "#     for a in f:\n",
    "#         for b in pa:\n",
    "#             v1=a[1]-a[0]\n",
    "#             v2=b[1]-b[0]\n",
    "#             v3=b[2]-b[0]\n",
    "#             c=array([v1,-v2,-v3]).transpose(1,0)\n",
    "#             d=b[0]-a[0]\n",
    "#             t=linalg.pinv(c)@d\n",
    "#             if (0<=t[0]<=1)&(0<=t[1]<=1)&(0<=t[2]<=1)&(0<=(t[1]+t[2])<=1):\n",
    "#                 pnt=a[0]+v1*t[0]\n",
    "                \n",
    "#                 ip1.append(pnt)\n",
    "    \n",
    "#     g=array([ip1,array(ip)[:,0],array(ip)[:,1]]).transpose(1,0,2)\n",
    "#     h=[fillet_3p_3d(p0,p1,p2,r_3p_3d([p0,p1,p2])*1.5,s) for (p0,p1,p2) in g]\n",
    "            \n",
    "#     return h+[h[0]]\n",
    "\n",
    "def ipfillet(p,p1,r=1,s=5,o=0):\n",
    "  \n",
    "    pa=[[[[p[i][j],p[i][j+1],p[i+1][j]],[p[i+1][j+1],p[i+1][j],p[i][j+1]]] if j<len(p[0])-1 else \\\n",
    "         [[p[i][j],p[i][0],p[i+1][j]],[p[i+1][0],p[i+1][j],p[i][0]]] \\\n",
    "         for j in range(len(p[0]))] for i in range(len(p)-1)]\n",
    "    pa=array(pa).reshape(-1,3,3)\n",
    "\n",
    "    pb=[[[p1[i][j],p1[i+1][j]] for j in range(len(p1[0]))] for i in range(len(p1)-1)]\n",
    "    pb=array(pb).reshape(-1,2,3)\n",
    "    a1,a2,a3=pa[:,0],pa[:,1],pa[:,2]\n",
    "    b1,b2=pb[:,0],pb[:,1]\n",
    "    v1,v2,v3=b2-b1,a2-a1,a3-a1\n",
    "    i,j=len(v1),len(v2)\n",
    "    v1=v1.repeat(j,0)\n",
    "    v2=array((v2).tolist()*i)\n",
    "    v3=array((v3).tolist()*i)\n",
    "    c=linalg.pinv(array([v1,-v2,-v3]).transpose(1,0,2).transpose(0,2,1))\n",
    "    d=array(a1.tolist()*i)-b1.repeat(j,0)\n",
    "    t=einsum('ijk,ik->ij',c,d)\n",
    "    p0=b1.repeat(j,0)\n",
    "    pnt=p0+einsum('ij,i->ij',v1,t[:,0])\n",
    "    v1norm=1/sqrt(einsum('ij,ij->i',v1,v1))\n",
    "    u1=einsum('ij,i->ij',v1,v1norm)\n",
    "    pnt1=pnt+u1*r\n",
    "    cond=(t[:,0]>=0)&(t[:,0]<=1)&(t[:,1]>=0)&(t[:,1]<=1)&(t[:,2]>=0)&(t[:,2]<=1)&((t[:,1]+t[:,2])>=0)&((t[:,1]+t[:,2])<=1)\n",
    "    pnt=pnt[cond]\n",
    "    pnt1=pnt1[cond]\n",
    "    ip=array([pnt,pnt1]).transpose(1,0,2)\n",
    "    if o==0:\n",
    "        e=[ip[i][0]+array(uv(cross(ip[i+1][0]-ip[i][0],ip[i][1]-ip[i][0])))*r if i<len(ip)-1 else \\\n",
    "           ip[i][0]+array(uv(cross(ip[0][0]-ip[i][0],ip[i][1]-ip[i][0])))*r \\\n",
    "           for i in range(len(ip))]\n",
    "    elif o==1:\n",
    "        e=[ip[i][0]+array(uv(cross(ip[i][1]-ip[i][0],ip[i+1][0]-ip[i][0])))*r if i<len(ip)-1 else \\\n",
    "           ip[i][0]+array(uv(cross(ip[i][1]-ip[i][0],ip[0][0]-ip[i][0])))*r \\\n",
    "           for i in range(len(ip))]\n",
    "\n",
    "    f=[[ip[i][0]+q(ip[i+1][0]-ip[i][0],e[i]-ip[i][0], theta) if i<len(ip)-1 else \\\n",
    "        ip[i][0]+q(ip[0][0]-ip[i][0],e[i]-ip[i][0] , theta) \\\n",
    "        for theta in linspace(-90,90,3)] for i in range(len(ip))]\n",
    "    f=array([seg(array(p).tolist())[0:-1] for p in f]).reshape(-1,2,3)\n",
    "\n",
    "    a1,a2,a3=pa[:,0],pa[:,1],pa[:,2]\n",
    "    b1,b2=f[:,0],f[:,1]\n",
    "    v1,v2,v3=b2-b1,a2-a1,a3-a1\n",
    "    i,j=len(v1),len(v2)\n",
    "    v1=v1.repeat(j,0)\n",
    "    v2=array((v2).tolist()*i)\n",
    "    v3=array((v3).tolist()*i)\n",
    "    c=linalg.pinv(array([v1,-v2,-v3]).transpose(1,0,2).transpose(0,2,1))\n",
    "    d=array(a1.tolist()*i)-b1.repeat(j,0)\n",
    "    t=einsum('ijk,ik->ij',c,d)\n",
    "    p0=b1.repeat(j,0)\n",
    "    pnt2=p0+einsum('ij,i->ij',v1,t[:,0])\n",
    "    cond=(t[:,0]>=0)&(t[:,0]<=1)&(t[:,1]>=0)&(t[:,1]<=1)&(t[:,2]>=0)&(t[:,2]<=1)&((t[:,1]+t[:,2])>=0)&((t[:,1]+t[:,2])<=1)\n",
    "    pnt2=pnt2[cond]\n",
    "\n",
    "    g=array([pnt2,pnt,pnt1]).transpose(1,0,2)\n",
    "    h=[fillet_3p_3d(p0,p1,p2,r_3p_3d([p0,p1,p2]),s) for (p0,p1,p2) in g]\n",
    "\n",
    "    return h+[h[0]]\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 260,
   "id": "bd441777",
   "metadata": {},
   "outputs": [],
   "source": [
    "# sec=cr_c(pts1([[0,0,.2],[8,3,3],[5,7,1],[-8,0,2],[-5,20,1]]),20)\n",
    "# sec=cr_c(pts1([[0,0,.1],[7,5,2],[5,7,3],[-5,7,5],[-7,5,5]]),20)\n",
    "# sec=cr_c(pts1([[-15,0,2.5],[0,15,3],[30,0,3],[0,-15,2.5],[5,0,2.5],[0,20,7],[-40,0,7],[0,-20,2.5]]),20)\n",
    "# sec=cr(pts1([[0,0,0.01],[5,0,0.01],[0,3,0.01],[-5,0,0.01]]),10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3694,
   "id": "fbb74ea5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# path=cr(pts1([[2,0],[-2,0,2],[0,7,4],[-4.2,0]]),40)\n",
    "# path=cr(pts1([[2,0],[-2,0,2],[0,7,5.2],[-5.2,0]]),20)\n",
    "# path=cr(pts1([[2,0],[-2,0,2],[0,7,2],[-2,0]]),20)\n",
    "# path=cr(pts1([[-5/2,0],[5/2,0],[0,5],[-5/2,0]]),1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "id": "33f9760d",
   "metadata": {},
   "outputs": [],
   "source": [
    "set_printoptions(suppress=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "1a735837",
   "metadata": {},
   "outputs": [],
   "source": [
    "cir=circle(10)\n",
    "sol1=l_extrude(cir,20)\n",
    "sketch1=[[-20,0],[20,0]]\n",
    "sketch2=[[-20,0,5],[20,0,15]]\n",
    "surf1=surf_extrude(sketch1,sketch2)\n",
    "ip=ip_surf(surf1,sol1)\n",
    "ip1=trns([0,0,3],ip)\n",
    "sol2=[ip,ip1]\n",
    "with open('/Users/sanjeevprabhakar/iCloud Drive (Archive)/Documents/Download/openscad/trial.scad','w+')as f:\n",
    "    f.write(f'''\n",
    "    include<dependencies.scad>\n",
    "    %{swp(surf1)}\n",
    "    %{swp(sol1)}\n",
    "    color(\"blue\")points({ip},.5);\n",
    "    color(\"green\")points({ip1},.5);\n",
    "    {swp(sol2)}\n",
    "    ''')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}