(* This program draws general CMC 1 surfaces in H3, using the Weierstrass data and 
              solving the Bryant equation 
              F^(-1) dF=2by2matrix(g,-g^2,1,-g) dz numerically.  Then the immersion 
              is F times its conjugate transpose.   The variable "dual" being either 
              0 or 1 determines whether the original surface (F times its conjugate 
              transpose) or the dual surface (F^(-1) times its conjugate transpose) 
              is made.   *)

r=9;  (* r=1 for a catenoid cousin 
         r=2 for an Enneper cousin, dual (dual=1) or original (dual=0)
         r=3 for a horosphere 
         r=4,5,6,7 for a genus 1 Enneper cousin dual 
               (then patch the four graphics figures together) 
         r=8 for a warped catenoid cousin
         r=9 for a trinoid cousin dual
         r=11 for example 7.4 of Umehara-Yamada paper in Annals of Math.
         r=13 for example 7.3 of Umehara-Yamada paper in Annals of Math.
         r=15 for a genus 0 4-noid with low total absolute curvature
         r=17 most symmetric 4-noid cousin
       *)

       (* BEGIN SURFACE DEPENDENT PART *)

(* first specify planar grid we use for triangulation *)

(* catenoid cousin: *)

If[r==1, 

h = 0.5; (* 1 must be integer multiple of h(=step size in algorithm) *)

xsteps = 15;
ysteps = 16;
xstepsize = -0.066;
ystepsize = 1.57075/4;
xbegin = 1.0; (* x is r, 0 < x < 1 *)
xend = 0.01; 
ybegin = 0.0; (* y goes from 0 to 2pi *)
yend = 6.283;
x0 = 1.0;
y0 = 0.0; 

mu = 0.8; 
kkk = -((mu^2.0)-1.0)/(4.0);

init11real = (1.0+mu)/(2*Sqrt[mu]);
init11imag = 0.0;
init12real = -(-1.0+mu)/(2*Sqrt[mu]);
init12imag = 0.0;
init21real = -(-1.0+mu)/(2*Sqrt[mu]);
init21imag = 0.0;
init22real = (1.0+mu)/(2*Sqrt[mu]);
init22imag = 0.0;

dual = 1;  (* take dual surface *)
];

(* Enneper cousin dual: *)

If[r==2,

h = 0.5;   (*  1 must be an integer multiple of h,
          it specifies step size in the algorithm  *)

xsteps = 40;
ysteps = 16;
xstepsize = 0.125;
ystepsize = 0.3925/4;
xbegin = 0.0; (* x is radial length *)
xend = 5.0;
ybegin = 0.0; (* y is theta *)
yend = 1.57;
x0 = 0.0;
y0 = 0.0; 

kkk = 1.0; (* coef in front of omega, can vary it *)

init11real = 1.0;
init11imag = 0.0;
init12real = 0.0;
init12imag = 0.0;
init21real = 0.0;
init21imag = 0.0;
init22real = 1.0;
init22imag = 0.0;

dual = 0;  (* 1 for dual surface, 0 for original surface *)

(*       alternate faster option, for use below:
v[i,j,1] = Cosh[Sqrt[kkk] z1];
v[i,j,2] = (1/Sqrt[kkk]) Sinh[Sqrt[kkk] z1]-z1 Cosh[Sqrt[kkk] z1];
v[i,j,3] = Sqrt[kkk] Sinh[Sqrt[kkk] z1];
v[i,j,4] = Cosh[Sqrt[kkk] z1] - Sqrt[kkk] z1 Sinh[Sqrt[kkk] z1];   *)

];

(* horosphere: *)

If[r==3,

h = 0.01;   (*  1 must be an integer multiple of h,
          it specifies step size in the algorithm  *)

xsteps = 64;
ysteps = 64;
xstepsize = 0.125/2;
ystepsize = 3.14159/32;
xbegin = 0.0;
xend = 4.0;
ybegin = 0.0;
yend = 6.28318;
x0 = 0.0;
y0 = 0.0;

kkk = 1.0; (* coef in front of omega, can vary it *)

init11real = 1.0;
init11imag = 0.0;
init12real = 0.0;
init12imag = 0.0;
init21real = 0.0;
init21imag = 0.0;
init22real = 1.0;
init22imag = 0.0;

dual = 1;  (* take dual surface *)
];

(* genus 1 Enneper cousin dual: *)

If[7.1>r>3.9,

h = 0.1;   (*  1 must be an integer multiple of h,
          it specifies step size in the algorithm  *)

avoiddividebyzero = 0.000000001; (* this should 
                   be a small, positive real *)
avoidnegativeimaginaryparts = avoiddividebyzero;
ybegin = 1.0+avoiddividebyzero; (* y=radial length *)
xbegin = 1.57+avoidnegativeimaginaryparts; 
                                (* x is theta *)
y0 = 1.0+avoiddividebyzero;
x0 = 1.57+avoidnegativeimaginaryparts; 

If[r==4,
ysteps = 20;
xsteps = 12;
ystepsize = -0.05;
xstepsize = -0.3925/3;
yend = 0.0+avoiddividebyzero;
xend = 0.0+avoidnegativeimaginaryparts;
];

If[r==5,
ysteps = 20;
xsteps = 12;
ystepsize = -0.05;
xstepsize = 0.3925/3;
yend = 0.0+avoiddividebyzero;
xend = 3.14+avoidnegativeimaginaryparts;
];

If[r==6,
ysteps = 60;
xsteps = 12;
ystepsize = 0.2;
xstepsize = -0.3925/3;
yend = 13.0+avoiddividebyzero;
xend = 0.0+avoidnegativeimaginaryparts;
];

If[r==7,
ysteps = 60;
xsteps = 12;
ystepsize = 0.2;
xstepsize = 0.3925/3;
yend = 13.0+avoiddividebyzero;
xend = 3.14+avoidnegativeimaginaryparts;
];

kkk = 0.25; (* coef in front of omega, can vary it *)
bbb = 0+Sqrt[NIntegrate[Sqrt[x]/Sqrt[1-x*x],{x,0,1}]/
NIntegrate[Sqrt[1-x*x]/Sqrt[x],{x,0,1}]];
           (* coef in front of g, can vary it *)
aaa = 1.0; (* variable in omega and g, can vary it *)

init11real = 1.00215;
init11imag = 0.0827759;
init12real = 0.0982006;
init12imag = -0.102215;
init21real = 0.0407553;
init21imag = 0.0397232;
init22real = 0.999062;
init22imag = -0.0827852;

dual = 1;  (* take dual surface *)

r=4;
];

(* warped catenoid cousins: *)

If[r==8, 

h = 0.05;   (*  1 must be an integer multiple of h,
          it specifies step size in the algorithm  *)

ell = 4; 
m = 1;
b = 0.5;

xx00 = (16 ell^2 (1+b^2)/((m^2-ell^2)^2))^0.25;
Print["xx00 is ",xx00];

xsteps = 36;
ysteps = 32;
xstepsize = -0.025;
ystepsize = 1.57075/16;
xbegin = xx00; (* x is r, 0 < x < 1 *)
xend = xx00-0.9; 
ybegin = 0.0; (* y goes from 0 to 2pi *)
yend = 6.283/2;
x0 = xbegin;
y0 = 0.0; 

rrr = (((ell-m)^2/(ell+m)^2) ((16 ell^2 xx00^(-2 m-2 ell)+
           16 ell^2 b^2 xx00^(-2 m-2 ell)+(ell+m)^4 xx00^(-2 m+2 ell)+
           8 ell b (ell+m)^2 xx00^(-2 m))/
           (16 ell^2 xx00^(2 m-2 ell)+16 ell^2 b^2 xx00^(2 m-2 ell)+
           (ell-m)^4 xx00^(2 m+2 ell)+8 ell b (ell-m)^2 xx00^(2 m) )))^0.25;

(* alternate faster method, for use below:   
v[i,j,1] = rrr Sqrt[(m^2-ell^2)/m] (1/(ell-m)) z1^((m-ell)/2);
v[i,j,2] = rrr Sqrt[(m^2-ell^2)/m] ((-b) (1/(ell-m)) z1^((m-ell)/2) + 
                        ((m-ell)/(4 ell)) z1^((ell+m)/2));
v[i,j,3] = Sqrt[(m^2-ell^2)/m] (1/(ell+m)) z1^((-m-ell)/2)/rrr;
v[i,j,4] = Sqrt[(m^2-ell^2)/m] ((-b) (1/(ell+m)) z1^((-m-ell)/2) + 
                        ((-m-ell)/(4 ell)) z1^((ell-m)/2))/rrr;      *)
         
init11real = Re[rrr Sqrt[(m^2-ell^2)/m] (1/(ell-m)) (xx00^2)^((m-ell)/2)];
init11imag = Im[rrr Sqrt[(m^2-ell^2)/m] (1/(ell-m)) (xx00^2)^((m-ell)/2)];
init12real = Re[rrr Sqrt[(m^2-ell^2)/m] ((-b) (1/(ell-m)) (xx00^2)^((m-ell)/2) + 
                        ((m-ell)/(4 ell)) (xx00^2)^((ell+m)/2))];
init12imag = Im[rrr Sqrt[(m^2-ell^2)/m] ((-b) (1/(ell-m)) (xx00^2)^((m-ell)/2) + 
                        ((m-ell)/(4 ell)) (xx00^2)^((ell+m)/2))];
init21real = Re[Sqrt[(m^2-ell^2)/m] (1/(ell+m)) (xx00^2)^((-m-ell)/2)/rrr];
init21imag = Im[Sqrt[(m^2-ell^2)/m] (1/(ell+m)) (xx00^2)^((-m-ell)/2)/rrr];
init22real = Re[Sqrt[(m^2-ell^2)/m] ((-b) (1/(ell+m)) (xx00^2)^((-m-ell)/2) + 
                        ((-m-ell)/(4 ell)) (xx00^2)^((ell-m)/2))/rrr];
init22imag = Im[Sqrt[(m^2-ell^2)/m] ((-b) (1/(ell+m)) (xx00^2)^((-m-ell)/2) + 
                        ((-m-ell)/(4 ell)) (xx00^2)^((ell-m)/2))/rrr];

dual = 0;  (* take original (not dual) surface *)
];

(* trinoid cousin: *)

If[r==9,

h = 0.1;

xsteps = 36;
ysteps = 32;
xstepsize = 0.1111111/4;
ystepsize = -0.1309/4;
xbegin = 0.0; (* x is r, 0 < x < 1 *)
xend = 0.9999999; 
ybegin = 1.0472; (* y is theta, 0 < y < pi/3= 1.0472 *)
yend = 0.0;
x0 = 0.0;
y0 = 0.0;

(* kkk = -0.3;   (* coefficient of omega *)
verticalmove = 0.94; *)

kkk = 0.8;   (* coefficient of omega *)
verticalmove = 1.05;

init11real = verticalmove;
init11imag = 0.0;
init12real = 0.0;
init12imag = 0.0;
init21real = 0.0;
init21imag = 0.0;
init22real = 1/verticalmove;
init22imag = 0.0;

dual = 1;
];

(* UY1, Example 7.4: *)

If[r==11, 

h = 0.05;   (*  1 must be an integer multiple of h,
          it specifies step size in the algorithm  *)

xsteps = 88;
ysteps = 32;
xstepsize = 0.0125;
ystepsize = 1.57075/32;
xbegin = 1.0; (* x is r, 0 < x < 1 *)
xend = 2.1; 
ybegin = 0.0; (* y goes from 0 to 2pi *)
yend = 6.283/4;
x0 = 1.0;
y0 = 0.0; 

m = 3; 

init11real = (-0.988/0.204813)/2.5;
init11imag = 0.0;
init12real = (Sqrt[1-0.988 0.988]/0.204813)/2.5;
init12imag = 0.0;
init21real = 2.5 (0.988 - Sqrt[1-0.988 0.988] 0.204813);
init21imag = 0.0;
init22real = 2.5 (-Sqrt[1-0.988 0.988] - 0.988 0.204813);
init22imag = 0.0;

dual = 0;  (* don't take dual surface *)
];

(* Umehara-Yamada paper in Annals of Math., example 7.3 *)

If[r==13,

h = 0.1;

xsteps = 20;
ysteps = 32;
xstepsize = -0.025;
ystepsize = -0.1309/4;
xbegin = 0.5; (* x is r, 0 < x < 1 *)
xend = 0.0; 
ybegin = 2*1.0472; (* y is theta, 0 < y < pi/3= 1.0472 *)
yend = 1.0472;
x0 = xbegin;
y0 = 2*1.0472;

m = 3;   
mu = 2.5;
a=(mu+m)/(mu-m);

rrr = 1.32;
theta = -1.0472;
init11real = (rrr/Sqrt[Abs[rrr^2-1]]) N[Cos[theta]];
init11imag = (rrr/Sqrt[Abs[rrr^2-1]]) N[Sin[theta]];
init12real = (1/Sqrt[Abs[rrr^2-1]]) N[Cos[theta]];
init12imag = -(1/Sqrt[Abs[rrr^2-1]]) N[Sin[theta]];
init21real = (1/Sqrt[Abs[rrr^2-1]]) N[Cos[theta]];
init21imag = (1/Sqrt[Abs[rrr^2-1]]) N[Sin[theta]];
init22real = (rrr/Sqrt[Abs[rrr^2-1]]) N[Cos[theta]];
init22imag = -(rrr/Sqrt[Abs[rrr^2-1]]) N[Sin[theta]];

init11real = 0.594378034539528599;
init11imag = -1.02949684247990047;
init12real = 0.294885854414141093;
init12imag = 0.510776142173220293;
init21real = 0.566165152442706176;
init21imag = -0.98062197480503972;
init22real = 0.701490421782231798;
init22imag = 1.21504200321622257;

dual = 0;
];

(* genus 0 4-noid cousin with low total absolute curvature: *)

If[r==15, 

h = 0.1; (* 1 must be integer multiple of h(=step size in algorithm) *)

ysteps = 30;
xsteps = 36;
ystepsize = 1.0/ysteps;
xstepsize = -1.57074/xsteps;
ybegin = 0.0; (* y is r, 0 < x < 1 *)
yend = 1.0; 
xbegin = 1.57075; (* x goes from pi/2 to 0 *)
xend = 0.00001;
y0 = 0.0;
x0 = 1.57075; 

aaaa=0.8;
bbbb=0.5;
beta=0.5;
(* good example is a=0.8 b=0.564 beta=-0.5, i11r=1.08737, i22r=0.919651 *)
(*  good example is a=0.9 b=0.763 beta=-0.75, i11r=1.12822, i22r=0.88635 *)
(*  good example is a=0.9 b=0.804 beta=-0.75, i11r=1.19561, i22r=0.83639 *)
(*  good example is a=0.97  beta=-6/7   b approx 0.875 I say, b in range 0.92, 0.93 Ume says,
            with my b, i11r=1.04942, i22r=0.952909 *)
(* Ume says good example is a=0.99  beta=-7/8,  b in range 0.945,0.9455  *)
(* Yam says good example is a=0.99  beta=-0.9,  b approx 0.952  *)
(* a=0.5, b=0.6, beta=0.8095, i11r=2.72403, i22r=0.367103 *)

pppp=Sqrt[(1/(6 bbbb bbbb)) (1-bbbb^4+Sqrt[1-14 bbbb^4+bbbb^8])];
lambda=-aaaa aaaa beta (beta+2) (aaaa^2+aaaa^(-2))^2/((aaaa aaaa+bbbb bbbb) (-1+aaaa aaaa bbbb bbbb));

init11real = 0.78888;
init11imag = 0.0;
init12real = 0.0;
init12imag = 0.0;
init21real = 0.0;
init21imag = 0.0;
init22real = 1.26762;
init22imag = 0.0;

dual = 1;  (* take dual surface *)
];

(* most symmetric 4-noid cousin: *)

If[r==17,

h = 0.1;

xsteps = 36;
ysteps = 32;
xstepsize = 0.9999999/xsteps;
ystepsize = -0.78539/ysteps;
xbegin = 0.0; (* x is r, 0 < x < 1 *)
xend = 0.9999999; 
ybegin = 0.78539; (* y is theta, 0 < y < pi/4= 0.78539 *)
yend = 0.0;
x0 = 0.0;
y0 = 0.0;

beta=-0.48;
kkk = 4 beta (beta+2);   (* coefficient of omega *)
verticalmove = 0.4;

init11real = verticalmove;
init11imag = 0.0;
init12real = 0.0;
init12imag = 0.0;
init21real = 0.0;
init21imag = 0.0;
init22real = 1/verticalmove;
init22imag = 0.0;

dual = 1;
];

         (* END SURFACE DEPENDENT PART *)

       (* BEGIN SURFACE INDEPENDENT PART *)

(* The array v is of the form v[i,j,k].  i and j 
specify grid pts. k specifies type of data in the entry:

k = 1  ------------- F11 (or F22 for dual)
k = 2  ------------- F12 (or -F12 for dual)
k = 3  ------------- F21 (or -F21 for dual)
k = 4  ------------- F22 (or F11 for dual)
k = 13 ------------- x1, actual projected R3 coordinate
k = 14 ------------- x2, actual projected R3 coordinate
k = 15 ------------- x3, actual projected R3 coordinate
*)

Array[v, {xsteps+1, ysteps+1, 15}];

(* Now we do first loop thru grid 1 <= i,j <= n .  
We define v[i,j,k] for k <= 8.  *)

oldx = x0;
oldy = y0;
i11r = init11real;
i11i = init11imag;
i12r = init12real;
i12i = init12imag;
i21r = init21real;
i21i = init21imag;
i22r = init22real;
i22i = init22imag;
 
j=0;
Do[ 
{i=0;
j=j+1;
Print["next y = ",y];
Do[ 
{i=i+1;

(* Now we specify the initial pt, end pt, and 
straight line path used in "Solve" fct, and 
do the numerical estimate algorithm *) 

s = .;

       (* END SURFACE INDEPENDENT PART *)

       (* BEGIN SURFACE DEPENDENT PART *)

(* Catenoid W data is: *)

If[r==1, 
z0 = (oldx) * (oldx) * Exp[I oldy];
z1 = (x) * (x) * Exp[I (y)];
path[s_] = (1-s+0. I)*z0 + (s + 0. I)*z1;
g[s_] = (1.0/path[s]); 
omega[s_] = kkk; 
];

(* Enneper W data is: *)

If[r==2, 
z0 = (oldx) * Exp[I oldy];
z1 = (x) * Exp[I (y)];
path[s_] = (1-s+0. I)*z0 + (s + 0. I)*z1;
g[s_] = path[s];
omega[s_] = kkk;
];

(* Horosphere W data is: *)

If[r==3, 
z0 = (oldx) * Exp[I oldy];
z1 = (x) * Exp[I (y)];
path[s_] = (1-s+0. I)*z0 + (s + 0. I)*z1;
g[s_] = 0.0;
omega[s_] = kkk;
];

(* genus 1 Enneper W data is: *)

If[r==4, 
z0 = (oldy) * Exp[I oldx];
z1 = (y) * Exp[I (x)];
path[s_] = (1-s+0. I)*z0 + (s + 0. I)*z1;
g[s_] = 
 bbb*Sqrt[path[s]-1]*Sqrt[path[s]+aaa]/Sqrt[path[s]];
omega[s_] = 
 kkk*Sqrt[path[s]]/(Sqrt[path[s]-1]*Sqrt[path[s]+aaa]);
];

(* warped catenoid cousin W data is: *)

If[r==8, 
z0 = (oldx) * (oldx) * Exp[I oldy];
z1 = (x) * (x) * Exp[I (y)];
path[s_] = (1-s+0. I)*z0 + (s + 0. I)*z1;
g[s_] = ((m m-ell ell)/(4 ell)) path[s]^ell + b; 
omega[s_] = 1/(path[s]^(ell+1)); 
];

(* Trinoid W data is: *)

If[r==9, 
z0 = (oldx)^0.15 * Exp[I (1/(1.0472^4)) oldy^5];
z1 = (x)^0.15 * Exp[I (1/(1.0472^4)) y^5];
path[s_] = (1-s+0. I)*z0 + (s + 0. I)*z1;
g[s_] = (path[s]^(2.0)); 
omega[s_] = (kkk)/(((path[s]^3.0)-1.0)^(2.0));
];

(* UY1 example 7.4 W data is: *)

If[r==11, 
z0 = (oldx) * (oldx) * Exp[I oldy];
z1 = (x) * (x) * Exp[I (y)];
path[s_] = (1-s+0. I)*z0 + (s + 0. I)*z1;
g[s_] = path[s]^m; 
omega[s_] = 1/(path[s]^2); 
];

(* UY1 example 7.3 W data is: *)

If[r==13, 
z0 = (oldx) * Exp[I oldy];
z1 = (x) * Exp[I y];
path[s_] = (1-s+0. I)*z0 + (s + 0. I)*z1;
g[s_] = (path[s]^mu) ((path[s]^m+a)/(a path[s]^m+1)); 
omega[s_] = ((1-mu^2)/(4 a mu)) path[s]^(-mu-1) ((a path[s]^m+1)/(path[s]^m+1))^2;
];

(* 4-noid with low TA's W data is: *)

If[r==15, 
z0 = Sqrt[oldy] * Exp[I oldx];
z1 = Sqrt[y] * Exp[I x];
path[s_] = (1-s+0. I)*z0 + (s + 0. I)*z1;
g[s_] = -(path[s] (pppp pppp path[s]^2-1))/(pppp^2-path[s]^2); 
omega[s_] = -lambda (((path[s]^2-
                  pppp pppp)^2)/((path[s]^2-
                    aaaa aaaa)^2 pppp pppp (path[s]^2-(aaaa^(-2)))^2 )); 
];

(* most symmetric 4-noid cousin W data is: *)

If[r==17, 
z0 = (oldx)^0.5 * Exp[I (1/(0.78539^4)) oldy^5];
z1 = (x)^0.5 * Exp[I (1/(0.78539^4)) y^5];
path[s_] = (1-s+0. I)*z0 + (s + 0. I)*z1;
g[s_] = (path[s]^(3.0)); 
omega[s_] = (kkk/3)/(((path[s]^4.0)-1.0)^(2.0));
];

       (* END SURFACE DEPENDENT PART *)

     (* BEGIN SURFACE INDEPENDENT PART *)

s = .;
oldar = i11r;
oldai = i11i;
oldbr = i12r;
oldbi = i12i;
Do[{
newar = oldar + h * 
      Re[(z1-z0)*(((oldar+oldai I)*g[s])+(oldbr+oldbi I))*omega[s]];
newai = oldai + h * 
      Im[(z1-z0)*(((oldar+oldai I)*g[s])+(oldbr+oldbi I))*omega[s]];
newbr = oldbr + h * 
      -Re[(z1-z0)*g[s]*(((oldar+oldai I)*g[s])+(oldbr+oldbi I))*omega[s]];
newbi = oldbi + h * 
      -Im[(z1-z0)*g[s]*(((oldar+oldai I)*g[s])+(oldbr+oldbi I))*omega[s]];
oldar = newar;
oldai = newai;
oldbr = newbr;
oldbi = newbi;
s = .;
},
{s, 0.0, 1.0 - h, h} ];
Ar = oldar;
Ai = oldai;
Br = oldbr;
Bi = oldbi;
v[i,j,1] = Ar + I Ai;
v[i,j,2] = Br + I Bi;

s = .;
oldar = i21r;
oldai = i21i;
oldbr = i22r;
oldbi = i22i;
Do[{
newar = oldar + h * 
      Re[(z1-z0)*(((oldar+oldai I)*g[s])+(oldbr+oldbi I))*omega[s]];
newai = oldai + h * 
      Im[(z1-z0)*(((oldar+oldai I)*g[s])+(oldbr+oldbi I))*omega[s]];
newbr = oldbr + h * 
      -Re[(z1-z0)*g[s]*(((oldar+oldai I)*g[s])+(oldbr+oldbi I))*omega[s]];
newbi = oldbi + h * 
      -Im[(z1-z0)*g[s]*(((oldar+oldai I)*g[s])+(oldbr+oldbi I))*omega[s]];
oldar = newar;
oldai = newai;
oldbr = newbr;
oldbi = newbi;
s = .;
},
{s, 0.0, 1.0 - h, h} ];
Cr = oldar;
Ci = oldai;
Dr = oldbr;
Di = oldbi;
v[i,j,3] = Cr + I Ci;
v[i,j,4] = Dr + I Di; 

s = .; 

(* end of numerical estimate algorithm *)

i11r = Ar;
i11i = Ai;
i12r = Br;
i12i = Bi;
i21r = Cr;
i21i = Ci;
i22r = Dr;
i22i = Di;
oldx = x;
oldy = y},

{x, xbegin, xend, xstepsize} ];

i11r = Re[v[1,j,1]];
i11i = Im[v[1,j,1]];
i12r = Re[v[1,j,2]];
i12i = Im[v[1,j,2]];
i21r = Re[v[1,j,3]];
i21i = Im[v[1,j,3]];
i22r = Re[v[1,j,4]];
i22i = Im[v[1,j,4]];
oldx = x0;
oldy = y},

{y, ybegin, yend, ystepsize} ];

(* Now do second loop thru grid 1 <= i,j <= n .  We 
define remaining v[i,j,k].  *)

j=0;Do[{i=0;j=j+1;Do[ {i=i+1;

If[dual==1,{temp = v[i,j,1];
v[i,j,1] = v[i,j,4];
v[i,j,2] = -v[i,j,2];
v[i,j,3] = -v[i,j,3];
v[i,j,4] = temp;}];

(* Now convert to SL(2,C) model for H3.  *)

f11 = (v[i,j,1])*Conjugate[v[i,j,1]] + 
        (v[i,j,2])*Conjugate[v[i,j,2]];
f12 = (v[i,j,1])*Conjugate[v[i,j,3]] + 
        (v[i,j,2])*Conjugate[v[i,j,4]];
f21 = (v[i,j,3])*Conjugate[v[i,j,1]] + 
        (v[i,j,4])*Conjugate[v[i,j,2]];
f22 = (v[i,j,3])*Conjugate[v[i,j,3]] + 
        (v[i,j,4])*Conjugate[v[i,j,4]];

(* Now convert to Minkowski model for H3.  *)

t = (f11 + f22)*(0.5);
x1 = (f12+f21)*(0.5);
x2 = (f12-f21)*(-0.5 I);
x3 = (f11-f22)*(0.5);

(* Now convert to Poincare model for H3.  *)

v[i,j,13] = Re[x1/(1+t)];
v[i,j,14] = Re[x2/(1+t)];
v[i,j,15] = Re[x3/(1+t)]; 


},
{x, xbegin, xend, xstepsize} ]},
{y, ybegin, yend, ystepsize} ];

(* Now we make the graphics.*)

forjavaview=Show[Table[Graphics3D[{EdgeForm[ ],Polygon[{Table[v[i,j,k],{k,13,15}],Table[v[i+1,j,k],
  {k,13,15}],Table[v[i,j+1,k],{k,13,15}]}]}],{i,1,xsteps},{j,1,ysteps}],
  Table[Graphics3D[{EdgeForm[ ],Polygon[{Table[v[i+1,j+1,k],{k,13,15}],Table[v[i+1,j,k],
  {k,13,15}],Table[v[i,j+1,k],{k,13,15}]}]}],{i,1,xsteps},{j,1,ysteps}],
  Graphics3D[{PointSize[0],Point[{v[1,1,13],v[1,1,14],v[1,1,15]}]}],
  ParametricPlot3D[{Cos[t],Sin[t],0},{t,0,2 Pi}],
  ParametricPlot3D[{Sqrt[3] Cos[t]/2,Sqrt[3] Sin[t]/2,1/2},{t,0,2 Pi}],
  ParametricPlot3D[{Sqrt[3] Cos[t]/2,Sqrt[3] Sin[t]/2,-1/2},{t,0,2 Pi}],
  ParametricPlot3D[{Cos[t],0,Sin[t]},{t,0,2 Pi}],
  ParametricPlot3D[{0,Cos[t],Sin[t]},{t,0,2 Pi}],
  PlotRange -> {{-1,1},{-1,1},{-1,1}}, Boxed -> False];

<< JavaView`RunThrough`

RunJavaView[forjavaview]