using MultivariatePolynomials
using SemialgebraicSets
using TypedPolynomials
function milnorpov(p,q,theta,eps=1, filename="/path/to/somefile.pov")
	#p has to be greater than q. If p and q are not coprime, i.e. if this produces a link rather than a knot, the code as it is only plot one of the components.
  @polyvar x y z
    P=(2*x*eps+im*2*y*eps)^p-(x^2+y^2+z^2+eps^2)^(p-q)*((2*z*eps+im *(x^2+y^2+z^2-eps^2)))^q
    RP=0
    IP=0
    for i in terms(P)
        RP+= real(i.coefficient)*i.monomial
        IP+=imag(i.coefficient)*i.monomial
    end
    d=max(2*p,2*q+2)
    strI=""
    strR=""
    IP=IP-tan(theta)*RP
    for i in d:-1:0
        for j in d:-1:0
            for k in d:-1:0
                if i+j+k<=d
                    strI*=string(coefficient(IP,x^i*y^j*z^k))
                    strI*=","
                    if coefficient(RP,x^i*y^j*z^k) >0
                        strR*="+"*string(coefficient(RP,x^i*y^j*z^k))*"*"*"pow(x,$i)"*"*"*"pow(y,$j)"*"*"*"pow(z,$k)"
                    end
                    if coefficient(RP,x^i*y^j*z^k) <0
                        strR*=string(coefficient(RP,x^i*y^j*z^k))*"*"*"pow(x,$i)"*"*"*"pow(y,$j)"*"*"*"pow(z,$k)"
                    end
                end
            end
        end
    end
    strI=chop(strI)
    a=@set x^p-y^q==0 && x^2+y^2==eps^2
    m=a[1][1]
    n=a[1][2]
    strK="$eps *$m*cos($q*ctr)/($eps- $n*sin($p*ctr)) , $eps*$m*sin($q*ctr)/($eps- $n*sin($p*ctr)) , $eps*$n*cos($p*ctr)/($eps- $n*sin($p*ctr))"
    out="
    #version version;
#include \"colors.inc\"
#include \"rad_def.inc\"
background { color White}
    #declare maxG=100000000;
#declare maxC=10;

global_settings { assumed_gamma 1.0
radiosity {
//	      Rad_Settings(Radiosity_IndoorHQ ,off,off) //disable this for faster rendering
	         }
	 }

camera
  {
  location <0, 0, -6.5>
  right x*image_width/image_height up <0,1,0>
  look_at  <0.0, 0.0, 0.0>
  angle 20
  }


light_source
  {
  2*<5, 10, -40>
  color White
  area_light <5, 0, 0>, <0, 5, 0>, 5, 5
  adaptive 1
  }

plane
  {
  <0,0,1>, 1.5
  pigment { color White }
  finish
    {
    ambient 0.0
    diffuse 0.5
    }
}
        #declare plastic= texture{
	    finish {
    ambient 0.0
   diffuse 0.3
    specular 0.2
    roughness 0.1
    reflection {0.0 metallic}
    }
    normal {
                        bumps 0.15 scale 0.05
                }
}

        #declare surf= texture{
	    finish {
    ambient 0.0
   diffuse 0.3
    specular 0.2
    roughness 0.1
    reflection {0.0 metallic}
    }
}
    	      union {
  #declare r_tube = 0.033366382419178446;  // thickness (radius) of tube

  #declare ctr=0;
  #while(ctr <= 6.28)
	  sphere{-0.4*<$strK>,r_tube

	  texture {plastic

pigment {color rgb <30/255,144/255,255/255>}
}

	      }
	      #declare ctr=ctr+0.001;
      #end

poly{$d,
	<$strI>
sturm on
clipped_by {isosurface{
function{$strR}
	threshold 0
	contained_by {box{maxC*<-1,-1,-1>, maxC*<1,1,1>}}
max_gradient maxG
	      }}
 texture {surf}
 pigment {color rgbf <255/255,143/255,31/255,0.3>}
scale -0.4
	      }
rotate <0,20,0> //change this to view the result from a different angle
translate <0,0,-0.5>
scale -0.9}
    "
    io=open(filename, "w")
    write(io, out)
    close(io)
end