Corrected code for polylog for use in custom_wf. There are 2 variables N & zpow (correspond to a & b in custom_wf) ---------------------------------------------------- import org.jwildfire.create.tina.base.XForm; import org.jwildfire.create.tina.variation.FlameTransformationContext; import org.jwildfire.create.tina.base.XYZPoint; import static org.jwildfire.base.mathlib.MathLib.*; public void init(FlameTransformationContext pContext, XForm pXForm) { N = a; // ********* USER PARAM ********* Riem = new RiemannInt(); HSTerm = HarmonicS( (int)N -1 ); } RiemannInt Riem; double HSTerm; double N; public double HarmonicS(int N) { // Defined in Crandalls paper if (N < 1) return 0.; double Hs=0.0 ; int i=1; for (; i 0) return 1.; return -1.; } public double Sig2() { // avoids returning 0 if (re>=0) return 1.; return -1.; } public double Mag2() { return re*re+im*im; } public double MagInv() { double M2 = this.Mag2(); return (M2 < 1e-8?1:1./M2); } public double Radius() { // slower than Mag2 return Math.sqrt(this.Mag2()); } public double Arg() { return Math.atan2(im, re); } public Complex ToP() { return new Complex( this.Radius() , this.Arg() ); } public Complex UnP() { return new Complex( re*Math.cos(im) , re*Math.sin(im) ); } // Save / Restore utils public void Save() { // saves the current state save_re = re; save_im = im; } public void Restore() { // revert to prev state, keeping into save re & im the current re = save_re; im = save_im; } public void Switch() { // revert to prev state, keeping into save re & im the current double r2 = save_re; double i2 = save_im; save_re = re; save_im = im; re = r2; im = i2; } public void Keep( Complex zz ) { // saves zz save_re = zz.re; save_im = zz.im; } public Complex Recall() { // gives you what was saved return new Complex( save_re, save_im); } public void NextPow() { // super useful in series!!! Do save before calling this or it wont work... this.Mul( this.Recall() ); } // Arith public void Sqr() { double r2 = re*re - im*im; double i2 = 2*re*im; re = r2; im = i2; } public void Scale( double mul ) { re = re*mul; im = im*mul; } public void Mul( Complex zz) { double r2 = re*zz.re - im*zz.im; double i2 = re*zz.im + im*zz.re; re = r2; im = i2; } public void Div( Complex zz) { // divides by zz double r2 = im*zz.im + re*zz.re; double i2 = im*zz.re - re*zz.im; double M2 = zz.MagInv(); re = r2*M2; im = i2*M2; } public void DivR( Complex zz) { // changes to zz / ( current value) double r2 = zz.im*im + zz.re*re; double i2 = zz.im*re - zz.re*im; double M2 = this.MagInv(); re = r2*M2; im = i2*M2; } public void Add( Complex zz) { re += zz.re; im += zz.im; } public void AMean( Complex zz) { this.Add(zz); this.Scale(0.5); } public void GMean( Complex zz) { Complex PF = new Complex(zz); PF.Sqr(); this.Sqr(); this.Add(PF); this.Pow(0.5); this.Scale(0.5); } public void Sub( Complex zz) { re -= zz.re; im -= zz.im; } public void SubR( Complex zz) { re = zz.re-re; im = zz.im-im; } public void Inc() { re += 1.0; } public void Dec() { re -= 1.0; } public void Pow( double ex) { Complex PF = this.ToP(); PF.re = Math.pow(PF.re,ex); PF.im = PF.im * ex; this.Copy( PF.UnP() ); } public void Exp() { re = Math.exp(re); this.Copy( this.UnP() ); } public void SinH() { double rr = 0.0; double ri = 0.0; double er = 1.0; re = Math.exp(re); er /= re; rr = 0.5 * (re - er); ri = rr + er; re = cos(im) * rr; im = sin(im) * ri; } public void CosH() { double rr = 0.0; double ri = 0.0; double er = 1.0; re = Math.exp(re); er /= re; rr = 0.5 * (re - er); ri = rr + er; re = cos(im) * ri; im = sin(im) * rr; } public void Log() { Complex PF = new Complex(this.Arg()); PF.im = 0.5 * Math.log( this.Mag2() ); // faster than log(mag()) PF.Flip(); // assign arg to cmplx part, log to real this.Copy( PF ); } public void Pow( Complex ex ) { this.Log(); this.Mul(ex); this.Exp(); } } public void transform(FlameTransformationContext pContext, XForm pXForm, XYZPoint pAffineTP, XYZPoint pVarTP, double pAmount) { // approx (very good) of Li[n](z) for n > 1 double vv = pAmount; double zpow = b; // ********* USER PARAM ********* it is a z exponent really Complex z = new Complex(pAffineTP.x, pAffineTP.y); z.Pow(zpow); z.Save(); if (z.Mag2() > 250000.0 || N >= 20 ) { // no convergence, or N too big... When N is big then Li tends to z pVarTP.x += vv * z.re; pVarTP.y += vv *z.im; return; } Complex LiN = new Complex(); int i; Complex T = new Complex(); Complex zPl1 = new Complex(z); if (z.Mag2() < 0.07 ) { // normal series. Li = sum((z^k)/(k^n)) for (i=1; i < 20; i++) { T.Copy( new Complex(pow(i,N) ) ); T.DivR(z); LiN.Add(T); z.NextPow(); } pVarTP.x += vv * LiN.re; pVarTP.y += vv * LiN.im; return; } // Crandall method (very simple and fast!) that uses Erdelyi series // from now on we will use ln(z) only so switch to it z.Log(); z.Save(); z.One(); for (i=0; i < 20; i++) { double zee = Riem.Z( (int)N-i ); if (zee != 0.0) { T.Copy(z); T.Scale( zee / (cern.jet.math.Arithmetic.longFactorial(i)) ); LiN.Add(T); } if (i == N-1) { zPl1.Copy(z); } z.NextPow(); } z.Restore(); // back to log(z) again... z.Neg(); z.Log(); z.Neg(); // -log(-log(z)) must be added now... z.re += HSTerm; T.Copy(z); z.Copy(zPl1); z.Scale(1.0 / cern.jet.math.Arithmetic.longFactorial((int)N-1) ); z.Mul(T); LiN.Add(z); pVarTP.x += vv * LiN.re; pVarTP.y += vv * LiN.im; return; } NOTE: this seems to only work with certain other variations. It works well with synth & julian. Others are hit & miss. This test is with hypertile3d2 & synth mode=3
Posted on: Sun, 17 Aug 2014 01:30:07 +0000
Recently Viewed Topics
© 2015