Sorry, just noticed something missing. This uses atanf,
which is a function that does what I think atan2 should
(given two coordinates y and x it computes the angle of
the polar coordinates of the point (x,y)). Here it is:


#declare atan_aux1=function(a,b) {atan2(min(a,b),max(a,b))}
#declare atan_aux2=function(a,b,aa,ba) {
  (2*(aa>ba)-1)*
  (2*(a>0)-1)*
  (2*(b>0)-1)*
  atan_aux1(aa,ba)
}
#declare atan_aux3=function(a,b,aa,ba) {
  (1+(aa>ba)*(2*(a<=0)-1))*(2*(b>0)-1)
}
#declare atan_aux4=function(a,b,aa,ba) {
  atan_aux3(a,b,aa,ba)*pi/2+atan_aux2(a,b,aa,ba)
}
#declare atanf=function(b,a) {
  atan_aux4(a,b,abs(a),abs(b))
}


And while I am looking at the code: This may be faster:


#macro pre_coil (l,r2,r1,steepness)
  #local rmax=(r1)+(r2);
  #local dmax=rmax*sqrt(2);
  #local omeg=(steepness)/(r2);
  #local pre_plane_stretch=1/(1+omeg*omeg);
  #local plane_stretch=sqrt(pre_plane_stretch);
  #local sqrsum=function(a,b) {a*a+b*b}
  #local sqrsum_=function(a,b,c) {a*a+b*b+c*c}
  #local coil_aux1=function(ca,sa,b,c) {atanf(c*ca+b*sa,b*ca-c*sa)}
  #local coil_aux2=function(a_,b,c) {coil_aux1(cos(a_),sin(a_),b,c)}
  #local coil_aux3=function(r,phi) {sqrsum(r-(r2),plane_stretch*r*phi)}
  isosurface {
    function {coil_aux3(sqrt(sqrsum(y,z)),coil_aux2(omeg*x,y,z))}
    contained_by {box {<0,-rmax,-rmax> <l,rmax,rmax>}}
    threshold (r1)*(r1)
    max_gradient 2*max(r2,sqrt(sqrsum_(
      (dmax-(r2)+pi*pi*pre_plane_stretch*dmax),
      (pi*pre_plane_stretch*dmax)
      (pi*omeg*pre_plane_stretch*dmax*dmax))))
    all_intersections
  }
#end

Post a reply to this message