/* This file is part of the source code for 3D-XplorMath-J, Version 1.0 (January 2008). * Copyright (c) 2008 The 3D-XplorMath Consortium (http://3d-xplormath.org). * This source code is released under a BSD License, which allows redistribution * in source and binary form, with or without modification, provided copyright * and license information are included, and with no warranty or guarantee of * any kind. For details, see http://3d-xplormath.org/j/source/BSDLicense.txt */ package vmm.surface.parametric; import vmm.core.Complex; import vmm.core.IntegerParam; import vmm.core.RealParamAnimateable; import vmm.core3D.ComplexVector3D; import vmm.core3D.GridTransformMatrix; import vmm.core3D.Vector3D; public class TwistedScherk extends WeierstrassMinimalSurface { private IntegerParam exponent = new IntegerParam("vmm.surface.parametric.TwistedScherk.MainEx",2); private RealParamAnimateable aa = new RealParamAnimateable("vmm.surface.parametric.TwistedScherk.aa",0.1,0,0.44); //private RealParamAnimateable bb = new RealParamAnimateable("vmm.surface.parametric.TwistedScherk.bb",0.9318,0,0.9); //private IntegerParam pieces = new IntegerParam("vmm.surface.parametric.Scherk_Weierstrass.pieces",1); //ComplexVector3D[][] integrationGrid; //boolean needsIntegrationGrid = true; public TwistedScherk() { super(); setFramesForMorphing(18); afp.reset(0,0,0); //addParameter(bb); addParameter(aa); aa.setMinimumValueForInput(0.0); aa.setMaximumValueForInput(0.4999); addParameter(exponent); //setDefaultViewpoint(new Vector3D(-44,9,25)) setDefaultViewpoint(new Vector3D(-3.5, 45, 22)); setDefaultViewUp(new Vector3D(-0.12, -0.45, 0.89)); setDefaultWindow(-5.2,5.2,-4,4); uPatchCount.setValueAndDefault(32); vPatchCount.setValueAndDefault(12); //SurfaceView.setGridSpacing(6); umin.reset(-4.8); umax.reset(4.8); vmin.reset(0); vmax.reset(1); removeParameter(umax); removeParameter(vmax); removeParameter(vmin); umin.setMaximumValueForInput(-0.1); //iFirstInHelper = false; multipleCopyOptions = new int[] {2}; // indicates that we can show twice the usual number of copies of fundamntal piece canShowConjugateSurface = true; wantsToSeeDomain = false; // true;// wantsToSeeGaussImage = false;//true;// if (wantsToSeeGaussImage) wantsToSeeDomain = true; if (wantsToSeeDomain) {setDefaultViewUp(new Vector3D(0,0,1)); setDefaultViewpoint(new Vector3D(0,0,40));} } private int ex,i0; private double twist, rs, rd, R, RR, rl, u0, r1, udiff; private Complex rotTwist, gs, cSym; private Vector3D transZ, normal, surf; private GridTransformMatrix Pxy, Rtw, Rtw_; // ============================= // ============================== // protected void createData() { //System.out.println("Before super.create"); super.createData(); data.discardGridTransforms(); if ((flag0)&&(!wantsToSeeDomain)) computeHalfPeriod(); if ((!inAssociateMorph)&&(flag0||flag05)&&(!wantsToSeeDomain) ) { GridTransformMatrix[] trList = new GridTransformMatrix[ex*8]; trList[0] = new GridTransformMatrix(); if (flag05){ trList[1] = new GridTransformMatrix().scale(1,-1,1); // reflect in x-z-plane data.addGridTransform(trList[1]); } else if (flag0) { trList[1] = new GridTransformMatrix().scale(-1,1,-1).reverse(); // rotate around y-axis data.addGridTransform(trList[1]); } for (int e = 1; e < ex; e++){ trList[2*e] = new GridTransformMatrix(trList[2*e-2]).rotateZ(360.0/ex); data.addGridTransform(trList[2*e]); trList[2*e+1] = new GridTransformMatrix(trList[2*e-1]).rotateZ(360.0/ex); data.addGridTransform(trList[2*e+1]); } if (flag0){ for (int e = 0; e < 2*ex; e++){ trList[2*ex + e] = new GridTransformMatrix(trList[e]).leftMultiplyBy(Pxy); data.addGridTransform(trList[2*ex + e]); } if (getNumberOfPieces()==2) { for (int e = 0; e < 2*ex; e++){ trList[4*ex + e] = new GridTransformMatrix(trList[e]).leftMultiplyBy(Rtw); data.addGridTransform(trList[4*ex + e]); } } // if pieces } // if flag0 } // if } public static double paramRescaleEnd(double x){ double y; y = Math.sin(Math.PI/2*x); return y; } public static double paramRescaleStart(double x){ double y; y = 1-Math.cos(Math.PI/2*x); return y; } public static double paramRescaleBoth(double x){ double y; y = (1.0-Math.cos(Math.PI/2*x))/2.0; return y; } /** * This grid is precisely adapted to the symmetry lines and ends of the surface, * overrides default Cartesian grid */ protected Complex domainGrid(double u, double v){ // return ( ); z = exp(u+i*v) Complex z,w; /*double a; double ui = umin.getValue(); double ua = umax.getValue(); if (u <= u0) a = ui + (u0-ui)*paramRescaleEnd((u-ui)/(u0-ui)); else a = u0 + (ua-u0)*paramRescaleStart((u-u0)/(ua-u0)); double r = Math.exp(a); */ double r = Math.exp(u); double ph = Math.PI/2*( 1 + paramRescaleEnd(v)); z = new Complex(r*Math.cos(ph), r*Math.sin(ph)); w = domainGridz(z); if (wantsToSeeGaussImage) w = gauss(w); return w; } protected Complex domainGridz(Complex z){ //double pSym = Math.PI/(2.0*ex); //cSym = new Complex(Math.cos(pSym),Math.sin(pSym)); //gs = gauss(cSym); Complex aux,w; w = new Complex(z.re + 1, z.im); aux = new Complex(1-z.re, -z.im); w.assignDivide(aux); // Polar centers at +1,-1 from w = (z+1)/(1-z), i --> i w.re = rs*w.re + rd; w.im = rs*w.im; // w --> rs*w + rd, Polar centers at RR, -1/RR, i --> rs*i+rd aux = w.log(); //if (aux.im < 0) aux.im = aux.im +2*Math.PI; aux.assignTimes(1.0/ex); w = aux.exponential(); // w --> w^{1/ex}, Polar centers at R, exp(i*pi/ex)/R //double test = (w.minus(cSym)).r(); //if (test < 0.0005) w.assignTimes(0.9995); // This is to show cSym in the grid return w; } // Index of grid point closest to zero and on real line protected void zeroIndex(){ double aux; double test = 10; Complex z; for (int i = 0; i < ucount ; i++){ z = domainGrid(umin.getValue() +i*du,vmax.getValue()); aux = z.r(); if ((aux <= test)&&(Math.abs(z.im)< 0.00001)) { test = aux; i0=i; } } } /** * The following two functions are the Weierstrass data that * define this surface. It is best shown on the above domainGrid. */ protected Complex gauss(Complex z) { // k = ex. gauss = z^{k-1}[(z^k - R^k)/i / (R^k*z^k + 1)]^twist Complex w = new Complex(z.integerPower(ex-1)); Complex aux = w.times(z); // z^k Complex cl = new Complex(aux.im, -aux.re + RR); // (z^k - R^k)/i aux.re = RR*aux.re + 1; aux.im = RR*aux.im; // R^k*z^k + 1 (deviates from notes?) cl.assignDivide(aux); //= cl.dividedBy(aux); cl = cl.log(); cl.assignTimes(twist); // = cl.times(twist); cl = cl.exponential(); w.assignTimes(cl); w.assignTimes(rotTwist); return w; } protected Complex hPrime(Complex z) { // return i/z/(z^ex - z^{-ex} -R^ex + R^{-ex}) // The domain Grid has polar centers at the poles! Complex aux = z.integerPower(ex); Complex w = aux.inverse(); aux.re = aux.re - w.re - 2*rd; aux.im = aux.im - w.im; aux.assignInvert(); w.re = -r1*aux.im; w.im = r1*aux.re; w.assignDivide(z); return w; } protected void redoConstants(){ // global variables independent of the points, hence thread safe super.redoConstants(); ex = exponent.getValue(); twist = (aa.getValue()); if (twist > 0.999/ex) { twist = 0.999/ex; aa.reset(twist,0,twist); } R = 0.93; if (!wantsToSeeDomain){ // The period computation gives no sign change for trivial gauss if (twist==0) R = 1.0; else { double Rtop = Math.max(0.03, Math.min(1.0, 1.0 - (2.0*ex/Math.PI )*twist )); Rtop = Math.exp(Math.log(Rtop)/ex); double Rbot = 0.8*Rtop; R = search(Rbot,Rtop, 1e-8);} System.out.println("Periods Above. Period Closing Parameter Below:"); System.out.println(R); } redoConstantsHere(); udiff = (umax.getValue() - umin.getValue())/2.0; // This is constant under repeated calls /* if (rl0*(ucount-1) > udiff/2.0){ int ir0 = (int)Math.round( 2.0*rl0/udiff*(ucount-1) ); udiff = 2.0*rl0*(ucount-1)/(double)ir0; }*/ umin.reset(rl - udiff); // How can one set umin dependent on ex? umax.reset(rl + udiff); zeroIndex(); // needs du //System.out.println(i0); } protected void redoConstantsHere(){ double tau = 0; rotTwist = new Complex(Math.cos(tau),Math.sin(tau)); r1 = Math.sqrt((ex - 0.5)*Math.sqrt(ex)); RR = Math.pow(R,ex); rs = 0.5*(RR + 1.0/RR); rd = 0.5*(RR - 1.0/RR); rl = Math.log(RR); // The middle line passes through the symmetry point u0 = Math.log(RR/(rs-rd)); // u0: The middle line passes the vertex at 0 double rl0 = Math.abs(rl-u0); // Make the y-axis a symmetry line: Complex gy = gauss(domainGridz(new Complex(-0.05,0))); rotTwist = rotTwist.times(new Complex(gy.re/gy.r(), -gy.im/gy.r() ) ); double pSym = Math.PI/(2.0*ex); cSym = new Complex(Math.cos(pSym),Math.sin(pSym)); } /** * We want to center the surface already at the helper Level. We cannot use the symmetry * lines of the surface and therefore need to integrate towards the image of ZERO_C. * getCenter is called from inside createHelperArray and not in Associate Morph. */ protected ComplexVector3D getCenter(){ double phSym = Math.PI/(2.0*ex); double eps_ = 0.00001; Complex eps_C = new Complex(eps_*Math.cos(phSym),eps_*Math.sin(phSym)); // detour ZERO_C Complex zi0 = domainGrid(umin.getValue()+i0*du,vmax.getValue()); //System.out.println(zi0.re); ComplexVector3D gC = new ComplexVector3D(helperArray[i0][vcount-1]); gC = gC.plus(ComplexVectorOneStepIntegrator(zi0, eps_C)); // place the origin on the symmetry axis //gC = new ComplexVector3D(ZERO_C,ZERO_C,ZERO_C); return gC; } public void computeHalfPeriod(){ if (!wantsToSeeDomain) { halfPeriod = new ComplexVector3D(helperArray[ucount-1][0]); gs = gauss(cSym); //= normal = surfaceNormal( rl,0); Complex gss = gs.times(gs); double cg = gss.re; double sg = gss.im; //Complex gss= gauss(domainGrid((int)Math.floor(ucount/2.0), (int)Math.floor(vcount/2.0))); //System.out.println("Before surf"); surf = surfacePoint( rl,0 ); Pxy = GridTransformMatrix.SetGridTransformMatrix(cg,sg,0,0, sg,-cg,0,0, 0,0,-1,2*surf.z, false); Vector3D yrot = Pxy.applyToNormal(new Vector3D(0,1,0)); //System.out.println(yrot.x);System.out.println(yrot.y);System.out.println(yrot.z); Rtw = GridTransformMatrix.SetGridTransformMatrix(yrot.y,-yrot.x,0,0, yrot.x,yrot.y,0,0, 0,0,1,-2*surf.z, true); /* System.out.println(' '); // This checks the Gauss map at the symmetry point System.out.println(cSym.re);System.out.println(cSym.im);System.out.println(cSym.r()); System.out.println(gs.re);System.out.println(gs.im);System.out.println(gs.r()); System.out.println(normal.x);System.out.println(normal.y);System.out.println(normal.z); System.out.println(surf.x);System.out.println(surf.y);System.out.println(surf.z); System.out.println(normal.y/normal.x); System.out.println(surf.y/surf.x); // period condition */ /* System.out.println(' '); // This checks the Gauss map on the straight lines Complex gy = gauss(domainGridz(new Complex(-0.05,0))); Complex gx = gauss(domainGridz(new Complex(-20.0,0))); System.out.println(gy.re);System.out.println(gy.im);//System.out.println(gs.r()); System.out.println(gx.re);System.out.println(gx.im); */ } //System.out.println("After CompHalfPeriod"); } /** * Override to close the hole around the image of ZERO_C. */ public Vector3D surfacePoint(double u, double v) { int i = (int) Math.floor(0.25 + (u-umin.getValue())/du ); int j = (int) Math.floor(0.25 + (v-vmin.getValue())/dv ); //Complex zInitial = domainGrid(umin.getValue() + i*du, vmin.getValue() + j*dv); //Complex z = domainGrid(u,v); ComplexVector3D auxW = new ComplexVector3D(helperArray[i][j].plus( ComplexVectorOneStepIntegrator(domainGrid(umin.getValue() + i*du, vmin.getValue() + j*dv), domainGrid(u,v)) )); //ComplexVector3D eW = new ComplexVector3D(auxW.y.minus(auxW.x), (auxW.y.plus(auxW.x)).times(I_C), auxW.z ) ; ComplexVector3D eW = new ComplexVector3D(helperToMinimal(auxW)); if (((i==i0)||(i==(i0+1)))&&(j==vcount-1)) eW = new ComplexVector3D(ZERO_C,ZERO_C,eW.z); if (AFP==0) if (wantsToSeeDomain) return new Vector3D(eW.z.re, eW.z.im, 0); else return eW.re(); else return (eW.re().times(Math.cos(AFP))).plus((eW.im().times(Math.sin(AFP)))); } /** * This function computes the period for given exponent and twist as function * of the position R etc of the poles of hPrime * @param s */ protected double comperiod(double s){ R = s; redoConstantsHere(); Complex[] periodpath = new Complex[2*ex+8]; ComplexVector3D[] vW = new ComplexVector3D[2*ex+8]; Vector3D[] sM = new Vector3D[2*ex+8]; double ph = Math.PI/ex; Complex eph = new Complex(Math.cos(ph), Math.sin(ph)); Complex gss = gauss(cSym); double period = 0.0; double sx = 0; double sy = 0; periodpath[0] = cSym.times(0.01); for (int i = 1; i < 2*ex+1; i++){ periodpath[i] = periodpath[i-1].times(eph); } periodpath[2*ex+1] = cSym.times(0.25); periodpath[2*ex+2] = cSym.times(0.5); periodpath[2*ex+3] = cSym.times(0.75); periodpath[2*ex+4] = cSym.times(1.0); vW[0] = new ComplexVector3D(ZERO_C, ZERO_C, ZERO_C); for (int i = 1; i < 2*ex+1; i++){ vW[i] = vW[i-1].plus( ComplexVectorOneStepIntegrator(periodpath[i-1], periodpath[i]) ); } for (int i = 2*ex+1; i < 2*ex+5; i++){ vW[i] = vW[i-1].plus( ComplexVectorIntegrator(periodpath[i-1], periodpath[i], 8) ); } for (int i = 0; i < 2*ex+5; i++){ sM[i] = helperToMinimal(vW[i]).re(); } for (int i = 0; i < 2*ex+1; i++){ sx = sx + sM[i].x; sy = sy + sM[i].y; } period = (sM[2*ex+4].y - sy/(2.0*ex))*gss.re - (sM[2*ex+4].x - sx/(2.0*ex))*gss.im; //System.out.println(sM[2*ex+4].x - sx/(2.0*ex));System.out.println(sM[2*ex+4].y - sy/(2.0*ex)); System.out.println(period); return period; } // =========== // The regula falsi starts here // ============ // //final static int ITMAX = 10000; // not used //final static double EPS = 1e-10; // not used //public static int stepCount; //static double sign( double v ) { return v<0?-1:v>0?1:0; } //static double abs( double v ) { return Math.abs(v); } //public static double search(RealFunctionOfOneVariable f, double left, /** * Searchs the root of function f in brackets [left,right] with * prescribed precision. * left and right must really bracket a root, e.g. f(left)*f(right)<0. * @param f real function * @param left bracket * @param right bracket * @param precision of root * @return a root of x */ double search(double left,double right, double precision) { // a and b are the endponts of the bracketing interval // and stay not orderd during the algorithm double a = left; double b = right; //double f_a = f.eval(a); double f_a = comperiod(a); //double f_b = f.eval(b); double f_b = comperiod(b); double u = 0.0; double v = 0.0; if (Math.abs(f_b) >= 1e-13){ if (sign(f_a) * sign(f_b) > 0) //TO DO: discuss this { b = a; f_b = f_a; a = 0.7*a; f_a = comperiod(a); if (sign(f_a) * sign(f_b) > 0){ b = a; f_b = f_a; a = 0.5*a; f_a = comperiod(a); if (sign(f_a) * sign(f_b) > 0){ b = a; f_b = f_a; a = 0.3*a; f_a = comperiod(a); if (sign(f_a) * sign(f_b) > 0) throw new IllegalArgumentException("root is not bracketed"); } } } } int count = 1; int iteration_steps = 0; if (f_b==0.0) { u = b; } else{ u = (a * f_b - b * f_a) / (f_b - f_a); if (Math.abs(f_b)< 1e-13){ // Deals only with very small f_a double d = (b-a)/1024.0; a = b - d; b = b + d; f_a = comperiod(a); //f_a = f.eval(a); f_b = comperiod(b); //f_b = f.eval(b); v = comperiod(u); //v = f.eval(u); if (sign(f_b) * sign(v) > 0) { if (sign(f_a) * sign(f_b) > 0) //to do: discuss this throw new IllegalArgumentException("root is not bracketed"); else u = (a * f_b - b * f_a) / (f_b - f_a); } else{ u = (b * v - u * f_b) / (v - f_b); } } // small f_a is treated. // .................. MAIN Iteration while (abs(a - b) > 2 * precision ) { //v = f.eval(u); v = comperiod(u); if (v == 0.0){ a = u; b = u; } else { if (sign(v) == sign(f_a)) { // BAD case, same endpoint a is improved count = count + 1; } else { // GOOD case, count was not increased, the other endpoint b improved double tmp1 = f_a; f_a = f_b; f_b = tmp1; double tmp2 = a; a = b; b = tmp2; count = 1; } final double weight = count * Math.min(0.25, (v / f_a) * (v / f_a)); // The weight (##) a = u; f_a = v; u = (a * f_b - b * f_a) / (f_b - f_a); // zero of secant //.......... The improved guess with the help of weight: .......... (##) u = (u + weight * b) / (1 + weight); iteration_steps++; } // else } // while END of iteration, error bound achieved. } // if (!(f_a==0)) stepCount = iteration_steps; return u; } }