(*** Input the expressions by Goldman, Cox, and Brenner ***)

(* Shifted the lubrication curves to match the numerical data by O'Neill for ζ = δ/a to order O(δ/a) *) 

Ftx[ζ_]:=(8/15)*Log[ζ]-0.9588;
Tty[ζ_]:=-(1/10)*Log[ζ]-0.1895;
Frx[ζ_]:=-(2/15)*Log[ζ]-0.2526;
Try[ζ_]:=(2/5)*Log[ζ]-0.3817;

(* Plots for the translating sphere case, these match the lines in the Goldman paper *)

LogLinearPlot[-Ftx[d],{d,0.001,10},PlotStyle->{Thick,Red},TicksStyle->Directive[Black,Thick,30],AxesLabel->{Style["δ/a",Black,Thick,30],Style["-(Subscript[F, x])/(6πμaU)",Black,Thick,30]},ImageSize->750]
LogLinearPlot[Tty[d],{d,0.001,10},PlotStyle->{Thick,Red},TicksStyle->Directive[Black,Thick,30],AxesLabel->{Style["δ/a",Black,Thick,30],Style["(Subscript[T, y])/(8((πμa)^2)U)",Black,Thick,30]},ImageSize->750]

(* Plots for the rotating sphere case, these too match the lines in the paper *)

LogLinearPlot[Frx[d],{d,0.001,10},PlotStyle->{Thick,Red},TicksStyle->Directive[Black,Thick,30],AxesLabel->{Style["δ/a",Black,Thick,30],Style["(Subscript[F, x])/(6((πμa)^2)Ω)",Black,Thick,30]},ImageSize->750]
LogLinearPlot[-Try[d],{d,0.00001,10},PlotStyle->{Thick,Red},TicksStyle->Directive[Black,Thick,30],AxesLabel->{Style["δ/a",Black,Thick,30],Style["-(Subscript[T, y])/(8((πμa)^3)Ω)",Black,Thick,30]},ImageSize->750]

(*** Extract the velocity and angular velocity for the applied torque ***)

(* Set of equations for a sphere both translating and rotating *)

fx=6*π*μ*a*(U*ftx+a*Ω*frx)
ty=8*π*μ*a^2*(U*tty+a*Ω*try)

(* Solve speed/rotation ratio for no external force (fx==0) and extract the velocity *)

solU=Solve[fx==0,Ω][[1]]
rhs=ty/.solU//Simplify;
vel=Solve[Ty==rhs,U][[1,1,2]]

(* Now plug the Goldman expressions in and obtain the velocity and angular velocity in terms of the gap size *)

subs={frx->Frx[δ/a],ftx->Ftx[δ/a],try->Try[δ/a],tty->Tty[δ/a]};
VV=vel/.subs
ΩΩ=solU[[1,2]]/.{U->vel}/.subs

(* Determine the relevant expression for the reduced velocity and study the limiting behavior *)

frac=VV/(ΩΩ*a);
frac//Simplify

Limit[frac,δ->0]

(* Lastly plot the reduced velocity as a function of swimmer-wall separation *)

pfrac=frac/.{δ->x*a}//Simplify;
LogLinearPlot[pfrac,{x,10^(-10),0.1},PlotStyle->{Thick,Red},TicksStyle->Directive[Black,Thick,30],AxesLabel->{Style["δ/a",Black,Thick,30],Style["V/(Ωa)",Black,Thick,30]},ImageSize->750]