/* POV-Ray 3.1 macros for using complex numbers
	contents include:
		Re(w), Im(w), 
		cAdd(w1,w2), cSub(w1,w2), 
		cMult(w1,w2), cDiv(w1,w2), cOneOver(w),
		cSquare(w), cCube(w), cSqrt(w),
		cCos(w), cSin(w), cTan(w),
		cACos(w), cASin(w), cAtan(w),
		cPow(w1,w2), cExp(w), cLog(w),
		cToThe(w,R)

	Mostly, these just return the result of complex
	arithmetic on complex arguments.

	There are sure to be major errors in the
	following code. Please send bug reports to:
		dsharp@interport.net
		( mailto://dsharp@interport.net )
*/

// some definitions for convenience

#declare E = 2.718281828;

#macro Cosh(b)
   		(exp(b)+exp(-b))/2
#end

#macro Sinh(b)
   		(exp(b)-exp(-b))/2
#end

#macro cString(w)
	concat("<",str(w.u,0,-1),", ",str(w.v,0,-1),">")
#end

#declare cOne=<1,0>;
#declare pOne=<1,0>;

#macro Polar2Cart(r,theta)
	#local uu=r*cos(theta);
	#local vv=r*sin(theta);
	<uu,vv>
#end
#macro cPolar2Cart(cpolar)
	#local uu=cpolar.u*cos(cpolar.v);
	#local vv=cpolar.u*sin(cpolar.v);
	<uu,vv>
#end
#macro cCart2Polar(w)
	<sqrt(w.u*w.u+w.v*w.v),atan2(w.v,w.u) >
#end


// w assumed "complex"  and represented as
// a two-d vector:
//		 w=<w.u,w.v> = w.u + i* w.v
// represented in 'cartesian' form)

#macro Re(w) // return the real part of w
	w.u
#end

#macro Im(w) // return the imaginary part of w
	w.v
#end

#macro cTimesI(w) // i*w
	<-w.v,w.u>
#end

#macro cConj(w)   // complex conjugate of w
	<w.u,-w.v>
#end

// magnitude of w
#macro cMag(w)
	sqrt(w.u*w.u+w.v*w.v)
#end

// 'argument' of w
#macro cArg(w)
	atan2(w.v,w.u)
#end

// cAdd and cSub are pretty superfluous since
// POV-Ray already adds and subtracts
// <uu,vv> vectors in the right way for
// our purposes
#macro cAdd(w1,w2) // returns w1+w2
	<w1.u+w2.u ,w1.v+w2.v>
#end
#macro cSub(w1,w2) // returns w1-w2
	<w1.u-w2.u ,w1.v-w2.v>
#end

#macro cMult(w1,w2) // returns w1*w2
	<w1.u*w2.u - w1.v*w2.v,w1.v*w2.u + w1.u*w2.v>
#end

#macro cDiv(w1,w2) // returns w1/w2
	#local absz2=w2.u*w2.u+w2.v*w2.v;
	<(w1.u*w2.u+w1.v*w2.v)/absz2, (w1.v*w2.u-w1.u*w2.v)/absz2>
#end

#macro cOneOver(w)  // 1/w
	#local den=w.u*w.u+w.v*w.v;
	<w.u/den,-w.v/den>
#end

#macro cSquare(w) // w^2
   #local ru=w.u*w.u-w.v*w.v;
   #local rv=2*w.u*w.v;
   <ru,rv>
#end

#macro cCube(w) // w^3
	cMult(cSquare(w),w)
#end

// square of <r,theta> form
#macro pSquare(p) // 
	<p.u*p.u,2*p.v>
#end

#macro cSqrt(w) // square root of w
	#if(w.v=0)
                #if(w.u>0=0)
                    #local retc=<sqrt(w.u),0>;
                #else
                     #local retc=<0,sqrt(-w.u)>;
                #end
	#else
		#local sqrtr=sqrt(sqrt(w.u*w.u+w.v*w.v));
		#local theta=atan2(w.v,w.u);
	
		#local retc=<sqrtr*cos(theta/2),sqrtr*sin(theta/2)>;
	#end
	retc
#end
// cosine of complex w
#macro cCos(w)
	#local ru=cos(w.u)*Cosh(w.v);
	#local rv=-sin(w.u)*Sinh(w.v);
	<ru,rv>
#end

// sine of complex w
#macro cSin(w)
	#local ru=sin(w.u)*Cosh(w.v);
	#local rv=cos(w.u)*Sinh(w.v);
	<ru,rv>
#end

// tangent of complex w
#macro cTan(w)
	cDiv(cSin(w),cCos(w))	
#end

// natural logarithm of w
#macro cLog(w)
     #if(w.u=0 & w.v=0)	//
		#local ru=0;	// these values for cLog(0) should be changed to taste
		#local rv=0;	// or just eliminate the "if... else...end"
					// so POV will generate an error if cLog is passed <0,0>.
	#else
		#local ru=log(sqrt(w.u*w.u+w.v*w.v));
		#local rv=atan2(w.v,w.u);
	#end
	<ru,rv>
#end



#macro cACos(w) //Arccos(z)=-i*Log{z+sqrt(z*z-1)}
       #local zz=w+cSqrt(cSquare(w)-cOne);
       #local ww=cLog(zz);
       <ww.v,-ww.u>
#end


#macro cASin(w) //Arcsin(z)=-i*Log{i*z+sqrt(1-z*z)}
	  #local zz=<-w.v,w.u>+cSqrt(cOne-cSquare(w));
       #local ww=cLog(zz);
       <ww.v,-ww.u>
#end

#macro cAtan(w) // Arctan(z)=i/2*Log{(1-i*z)/(1+i*z)}
       #local zz=cLog(cDiv(<1+w.v,-w.u>,<1-w.v,w.u>));
       #local zz=<-zz.v,zz.u>/2;
       zz
#end

#macro cExp(w)	
	<exp(w.u)*cos(w.v),exp(w.u)*sin(w.v)>
#end

#macro cPow(w1,w2) // w1 to the w2 power
	#local A=cLog(w1);
	#local B=cMult(w2,A);
	cExp(B)
#end



#macro cToThe(w,R) // complex number to real power
	#local powr=pow(sqrt(w.u*w.u+w.v*w.v),R);
	#local theta=atan2(w.v,w.u);
	<powr*cos(theta*R),powr*sin(theta*R)>
#end