#declare r0=3       //Torus major radius r0
#declare r1=5       //Torus minor radius r1
#declare ri=9       //Inversion radius ri
#declare dx=3       //Inversion center displacement dx
#declare dy=0       //Inversion center displacement dy

#declare a0 =(pow(r1,2) - pow(dy,2) - pow(-r0 + dx,2))*
              (pow(r1,2) - pow(dy,2) - pow(r0 + dx,2))
#debug str(a0,0,-1)
#declare a1 = 0
#declare a2 = 0 
#declare a3 = (4*pow(r0,2)*pow(ri,2)*dx + 4*pow(r1,2)*
                pow(ri,2)*dx - 4*pow(ri,2)*pow(dy,2)*dx -
                4*pow(ri,2)*pow(dx,3)) 
#declare a4 = (2*(pow(r1,2) - pow(dy,2) - pow(-r0 + dx,2))*
               (pow(r1,2) - pow(dy,2) -    pow(r0 + dx,2))) 
#declare a5 = 0 
#declare a6 =(-4*pow(r0,2)*pow(ri,2)*dy + 4*pow(r1,2)*pow(ri,2)*
               dy - 4*pow(ri,2)*pow(dy,3) - 4*pow(ri,2)*dy*pow(dx,2)) 
#declare a7 =( 2*(pow(r1,2) - pow(dy,2) - pow(-r0 + dx,2))*
               (pow(r1,2) - pow(dy,2) -  pow(r0 + dx,2))) 
#declare a8 = 0 
#declare a9 =( -2*pow(r0,2)*pow(ri,4) -2*pow(r1,2)*pow(ri,4) +
               2*pow(ri,4)*pow(dy,2) +  6*pow(ri,4)*pow(dx,2)) 
#declare a10 =0 
#declare a11 = 0 
#declare a12 =( 4*pow(r0,2)*pow(ri,2)*dx +4*pow(r1,2)*pow(ri,2)*dx -
                4*pow(ri,2)*pow(dy,2)*dx - 4*pow(ri,2)*pow(dx,3)) 
#declare a13 = 0 
#declare a14 = 0 
#declare a15 = 0 
#declare a16 = 0 
#declare a17 = (4*pow(r0,2)*pow(ri,2)*dx +4*pow(r1,2)*pow(ri,2)*dx-
                4*pow(ri,2)*pow(dy,2)*dx -   4*pow(ri,2)*pow(dx,3)) 
#declare a18 = 0 
#declare a19 = (-4*pow(ri,6)*dx)   
#declare a20 = ((pow(r1,2) - pow(dy,2) - pow(-r0 + dx,2))*(pow(r1,2) -
                pow(dy,2) -  pow(r0 + dx,2))) 
#declare a21 = 0 
#declare a22 = (-4*pow(r0,2)*pow(ri,2)*dy + 4*pow(r1,2)*pow(ri,2)*dy -
                4*pow(ri,2)*pow(dy,3) -    4*pow(ri,2)*dy*pow(dx,2)) 
#declare a23 = (2*(pow(r1,2) - pow(dy,2) -pow(-r0 + dx,2))* (pow(r1,2) -
                pow(dy,2) -    pow(r0 + dx,2))) 
#declare a24 = 0 
#declare a25 = (2*pow(r0,2)*pow(ri,4) - 2*pow(r1,2)*pow(ri,4) +
                6*pow(ri,4)*pow(dy,2) + 2*pow(ri,4)*pow(dx,2)) 
#declare a26 = 0 
#declare a27 = (-4*pow(r0,2)*pow(ri,2)*dy +4*pow(r1,2)*pow(ri,2)*dy -
                 4*pow(ri,2)*pow(dy,3) - 4*pow(ri,2)*dy*pow(dx,2)) 
#declare a28 = 0 
#declare a29 = (-4*pow(ri,6)*dy)  
#declare a30 = (8*pow(ri,4)*dy*dx+ (pow(r1,2) - pow(dy,2) -
                  pow(-r0 + dx,2))*(pow(r1,2) - pow(dy,2) -
                  pow(r0 + dx,2))) 
#declare a31 = 0 
#declare a32 =(-2*pow(r0,2)*pow(ri,4) -2*pow(r1,2)*pow(ri,4) +
                2*pow(ri,4)*pow(dy,2) +  2*pow(ri,4)*pow(dx,2)) 
#declare a33 =0 
#declare a34 = (pow(ri,8))
#debug str(a34,0,-1)
#declare dup=poly{4 
               <a0, a1 ,a2 ,a3 ,a4 ,a5 ,a6 ,a7 ,a8 ,
               a9 ,a10 ,a11 ,a12 ,a13 ,a14 ,a15 ,a16 , 
               a17 ,a18 ,a19 ,a20 ,a21 ,a22 ,a23 ,a24 ,
               a25 ,a26 ,a27 ,a28 ,a29 ,a30 ,a31 ,a32 ,
               a33 ,a34 >}