ThI = -2*gamma*Sign[y]*ArcTan[Rmx/Rpx] + Sign[x*y]*(-Rma+2*gamma*ArcTan[Rmalp/Rpalp] )

Rpx = (1/Sqrt[2])*Sqrt[Sqrt[x^2+y^2]+x]

Rmx = (1/Sqrt[2])*Sqrt[Sqrt[x^2+y^2]-x]

Rpa = (1/Sqrt[2])*Sqrt[Sqrt[a^2+b^2]+a]

Rma = (1/Sqrt[2])*Sqrt[Sqrt[a^2+b^2]-a]

Rpalp = (1/Sqrt[2])*Sqrt[Sqrt[alp^2+bet^2]+alp]

Rmalp = (1/Sqrt[2])*Sqrt[Sqrt[alp^2+bet^2]-alp]

a = x^2-y^2+gamma^2

b = 2*x*y

alp = Rpa+gamma

bet = Sign[b]*Rma

(*  We'll look for Stokes lines in the lower half z plane.   I'll define a function f = ThI - Pi*gamma/2 and then look for its zeros. *)

ThI/.{x->0,y->-Pi/2}

f = Expand[ThI - Pi*gamma/2];

(* Set gamma = 1 *)
gamma = 1

(* Note that the equation satisfied by the Stokes lines is f = 0.  So we find the values of y for various values of x for which f = 0.  Here x will be positive *)

St1 = Table[{N[i/100],Part[FindRoot[f/.{x->i/100},{y,-1}],1,2]},{i,1,1000}];

St1 = Prepend[St1,{0,-1}];

ThR = -Rpa-(gamma/2)*Log[x^2+y^2] + (gamma/2)*Log[alp^2+bet^2]

FindRoot[ThR/.{x->0/10},{y,-1}]

FindRoot[ThR/.{x->1/100},{y,-1}]

St2 = Table[{N[i/1000],Part[FindRoot[ThR/.{x->i/1000},{y,-1+i/1000}],1,2]},{i,0,662}];

St2append = Part[FindRoot[ThR/.{y->0},{x,0.66}],1,2]

St2 = Append[St2,{St2append,0}];

St3 = Table[{0,N[-1+i/10]},{i,0,10}]

p1 = ListPlot[{St1,St2,St3},PlotRange->{{-1,5},{-2,0}},Axes->False,Frame->True,Joined->True]

Export["C:\\Users\\anderson\\Research\\de Sitter States\\Dillon\\Stokes_dS_1.dat",Grid[St1]]

Export["C:\\Users\\anderson\\Research\\de Sitter States\\Dillon\\Stokes_dS_2.dat",Grid[St2]]

Export["C:\\Users\\anderson\\Research\\de Sitter States\\Dillon\\Stokes_dS_3.dat",Grid[St3]]

(* OK, let's try the other side of the axis.  *)

g = Expand[ThI + Pi*gamma/2];

St4 = Table[{N[i/100],Part[FindRoot[g/.{x->i/100},{y,1}],1,2]},{i,1,1000}];

St4 = Prepend[St4,{0,1}];

St5 = Table[{N[i/1000],Part[FindRoot[ThR/.{x->i/1000},{y,1-i/1000}],1,2]},{i,1,660}];

St5append = Part[FindRoot[ThR/.{y->0},{x,0.66}],1,2]

St5 = Append[St5,{St5append,0}];

St6 = Table[{0,N[1-i/10]},{i,0,10}]

p2 = ListPlot[{St4,St5,St6},PlotRange->{{-1,5},{-2,2}},Axes->False,Frame->True,Joined->True]

p3 = Show[p1,p2,PlotRange->{{-1,5},{-2,2}},Axes->True]

Export["C:\\Users\\anderson\\Research\\de Sitter States\\Dillon\\Stokes_z_plot.pdf",p3]

Export["C:\\Users\\anderson\\Research\\de Sitter States\\Dillon\\Stokes_dS_4.dat",Grid[St4]]

Export["C:\\Users\\anderson\\Research\\de Sitter States\\Dillon\\Stokes_dS_5.dat",Grid[St5]]

Export["C:\\Users\\anderson\\Research\\de Sitter States\\Dillon\\Stokes_dS_6.dat",Grid[St6]]

(* OK now try the u coordinate *)

uR=-(1/2) *Log[ x^2+y^2]
uI=-ArcSin[y/Sqrt[x^2+y^2]]

Su1= Table[{N[-(1/2)*Log[St1[[i,1]]^2+St1[[i,2]]^2]],N[-ArcSin[St1[[i,2]]/Sqrt[St1[[i,1]]^2+St1[[i,2]]^2]]]},{i,1,1001}];

Su2= Table[{-(1/2)*Log[St2[[i,1]]^2+St2[[i,2]]^2],-ArcSin[St2[[i,2]]/Sqrt[St2[[i,1]]^2+St2[[i,2]]^2]]},{i,1,662}];

(* For the next one x = 0 so *)

uR/.{x->0}

uI/.{x->0}

(* So uI = Pi/2  if y <0 and uI = -Pi/2 if y > 0.  Also since |y| ≤1, uR >0 *)

Su3= Table[{N[i/10],N[Pi/2]},{i,0,21}]

p4 = ListPlot[{Su1,Su2,Su3},Joined->True]

Su4 = Table[{Su1[[i,1]],N[Su1[[i,2]]+2*Pi]},{i,1001}];

Su5 = Table[{Su2[[i,1]],N[Su2[[i,2]]+2*Pi]},{i,662}];

Su6 = Table[{Su3[[i,1]],N[Su3[[i,2]]+2*Pi]},{i,21}];

p5 = ListPlot[{Su4,Su5,Su6},Joined->True]

Su7 = Table[{Su4[[i,1]],N[Su4[[i,2]]+2*Pi]},{i,1001}];

Su8 = Table[{Su5[[i,1]],N[Su5[[i,2]]+2*Pi]},{i,662}];

Su9 = Table[{Su6[[i,1]],N[Su6[[i,2]]+2*Pi]},{i,21}];

p6 = ListPlot[{Su7,Su8,Su9},Joined->True]

Show[p4,p5,p6,PlotRange->All]

Export["C:\\Users\\anderson\\Research\\de Sitter States\\Dillon\\Stokes_dS_u_1.dat",Grid[Su1]]

Export["C:\\Users\\anderson\\Research\\de Sitter States\\Dillon\\Stokes_dS_u_2.dat",Grid[Su2]]

Export["C:\\Users\\anderson\\Research\\de Sitter States\\Dillon\\Stokes_dS_u_3.dat",Grid[Su3]]

Export["C:\\Users\\anderson\\Research\\de Sitter States\\Dillon\\Stokes_dS_u_4.dat",Grid[Su4]]

Export["C:\\Users\\anderson\\Research\\de Sitter States\\Dillon\\Stokes_dS_u_5.dat",Grid[Su5]]

Export["C:\\Users\\anderson\\Research\\de Sitter States\\Dillon\\Stokes_dS_u_6.dat",Grid[Su6]]

Export["C:\\Users\\anderson\\Research\\de Sitter States\\Dillon\\Stokes_dS_u_7.dat",Grid[Su7]]

Export["C:\\Users\\anderson\\Research\\de Sitter States\\Dillon\\Stokes_dS_u_8.dat",Grid[Su8]]

Export["C:\\Users\\anderson\\Research\\de Sitter States\\Dillon\\Stokes_dS_u_9.dat",Grid[Su9]]

(* Now we'll include the ones from the other side of the z-plot. *)

Su10= Table[{N[-(1/2)*Log[St4[[i,1]]^2+St4[[i,2]]^2]],N[-ArcSin[St4[[i,2]]/Sqrt[St4[[i,1]]^2+St4[[i,2]]^2]]]},{i,1,1001}];

Su11= Table[{-(1/2)*Log[St5[[i,1]]^2+St5[[i,2]]^2],-ArcSin[St5[[i,2]]/Sqrt[St5[[i,1]]^2+St5[[i,2]]^2]]},{i,1,661}];

Su12= Table[{N[i/10],N[-Pi/2]},{i,0,21}]

p7 = ListPlot[{Su10,Su11,Su12},Joined->True]

Show[p4,p7,PlotRange->All]

(* Now do one at uI = -Pi/2 + 2Pi = 3 Pi/2 *)

Su13 = Table[{Su10[[i,1]],N[Su10[[i,2]]+2*Pi]},{i,1001}];

Su14 = Table[{Su11[[i,1]],N[Su11[[i,2]]+2*Pi]},{i,661}];

Su15 = Table[{Su12[[i,1]],N[Su12[[i,2]]+2*Pi]},{i,21}];

p8= ListPlot[{Su13,Su14,Su15},Joined->True]

Show[p4,p5,p7,p8,PlotRange->All]

(* OK, let's find the anti-Stokes and Stokes lines for x < 0 in the z-coordinate. *)

St7 = Table[{N[-i/100],Part[FindRoot[f/.{x->-i/100},{y,-1}],1,2]},{i,1,1000}];

St7= Prepend[St7,{0,-1}];

St8 = Table[{N[-i/1000],Part[FindRoot[ThR/.{x->-i/1000},{y,-1+i/1000}],1,2]},{i,0,662}];

St8append = Part[FindRoot[ThR/.{y->0},{x,-0.66}],1,2]

St8 = Append[St8,{-0.6627434193491816,0}];

p9 = ListPlot[{St7,St8},PlotRange->All]

Show[p1,p2,p9,PlotRange->{{-1,5},{-2,2}},Axes->True]

St9 = Table[{N[-i/100],Part[FindRoot[g/.{x->-i/100},{y,1}],1,2]},{i,1,1000}];

St9= Prepend[St9,{0,1}];

St10 = Table[{N[-i/1000],Part[FindRoot[ThR/.{x->-i/1000},{y,1-i/1000}],1,2]},{i,1,662}];

St8append

St10 = Append[St10,{St8append ,0.}];

p10 = ListPlot[{St9,St10},PlotRange->{{-1,0},{0,1}},Joined->True]

Show[p1,p2,p9,p10,PlotRange->{{-1,5},{-2,2}},Axes->True]

(* Now translate to the u coordinates *)

(* NOTE THAT I SHOULD HAVE USED A DIFFERENT TRANSFORMATION FOR THE IMAGINARY PART OF U.  IT SHOULD HAVE BEEN + ArcSin[y/Sqrt[x^2+y^2] + Pi.
  HOWEVER IT DOESN'T MAKE ANY DIFFERENCE BECAUSE I CONSIDER BOTH POSITIVE VALUES OF y AND NEGATIVE VALUES OF y WHEN x < 0.  So I will leave the 
program as it is.  *)

Su16= Table[{N[-(1/2)*Log[St7[[i,1]]^2+St7[[i,2]]^2]],N[-ArcSin[St7[[i,2]]/Sqrt[St7[[i,1]]^2+St7[[i,2]]^2]]]+Pi},{i,1,1001}];

Su17= Table[{-(1/2)*Log[St8[[i,1]]^2+St8[[i,2]]^2],-ArcSin[St8[[i,2]]/Sqrt[St8[[i,1]]^2+St8[[i,2]]^2]]+Pi},{i,1,663}];

p11= ListPlot[{Su16,Su17},Joined->True]

Su18= Table[{N[-(1/2)*Log[St9[[i,1]]^2+St9[[i,2]]^2]],N[-ArcSin[St9[[i,2]]/Sqrt[St9[[i,1]]^2+St9[[i,2]]^2]]]+Pi},{i,1,1001}];

Su19= Table[{-(1/2)*Log[St10[[i,1]]^2+St10[[i,2]]^2],-ArcSin[St10[[i,2]]/Sqrt[St10[[i,1]]^2+St10[[i,2]]^2]]+Pi},{i,1,663}];

p12 = ListPlot[{Su18,Su19},Joined->True,PlotRange->All]

Show[p4,p5,p7,p8,p11,p12,PlotRange->All]

(* So for the u plot above I'm not getting any extra lines for some reason - maybe the way the translation is occurring.  *)

Show[p1,p2,p9,p10,PlotRange->{{-1,5},{-2,2}},Axes->True]

(* For the z plot above I do get the extra lines *)

Export["C:\\Users\\anderson\\Research\\de Sitter States\\Dillon\\Stokes_dS_7.dat",Grid[St7]]

Export["C:\\Users\\anderson\\Research\\de Sitter States\\Dillon\\Stokes_dS_8.dat",Grid[St8]]

Export["C:\\Users\\anderson\\Research\\de Sitter States\\Dillon\\Stokes_dS_9.dat",Grid[St9]]

Export["C:\\Users\\anderson\\Research\\de Sitter States\\Dillon\\Stokes_dS_10.dat",Grid[St10]]

Export["C:\\Users\\anderson\\Research\\de Sitter States\\Dillon\\Stokes_dS_u_10.dat",Grid[Su10]]

Export["C:\\Users\\anderson\\Research\\de Sitter States\\Dillon\\Stokes_dS_u_11.dat",Grid[Su11]]

Export["C:\\Users\\anderson\\Research\\de Sitter States\\Dillon\\Stokes_dS_u_12.dat",Grid[Su12]]

Export["C:\\Users\\anderson\\Research\\de Sitter States\\Dillon\\Stokes_dS_u_13.dat",Grid[Su13]]

Export["C:\\Users\\anderson\\Research\\de Sitter States\\Dillon\\Stokes_dS_u_14.dat",Grid[Su14]]

Export["C:\\Users\\anderson\\Research\\de Sitter States\\Dillon\\Stokes_dS_u_15.dat",Grid[Su15]]

Export["C:\\Users\\anderson\\Research\\de Sitter States\\Dillon\\Stokes_dS_u_16.dat",Grid[Su16]]

Export["C:\\Users\\anderson\\Research\\de Sitter States\\Dillon\\Stokes_dS_u_17.dat",Grid[Su17]]

Export["C:\\Users\\anderson\\Research\\de Sitter States\\Dillon\\Stokes_dS_u_18.dat",Grid[Su18]]

Export["C:\\Users\\anderson\\Research\\de Sitter States\\Dillon\\Stokes_dS_u_19.dat",Grid[Su19]]

N[-Pi/2,16]

N[3*Pi/2,16]

N[5*Pi/2,16]