Clear[ inverseRadialSmoothedCap ] 

inverseRadialSmoothedCap[rr_, ee_, a_, b_] = 
Module[{vPlusCubic, rrn},
vPlusCubic[r_, e_, alpha_, beta_] =
Module[{vv, ww, s, v, aa, bb},
vv[r_] = aa ( r - e)^3/e^4 + bb(r - e)^2/ e^3 + beta(r - e)/e^2 - beta/e ;
ww[r_] = vv'[r] ;
s= Solve[ {-(beta + alpha)/e == vv[0], 0 == ww[0]}, {aa, bb}]  // Flatten ;

v[r_, e_, alpha_, beta_] = vv[r] /. s  ;

(*w[r_, e_, alpha_, beta_] = D[v[r, e, alpha, beta], r] ;
{v[e, e, α, β], w[e,e, α, β], v[0, e, α, β], w[0, e, α, β]} // Simplify*)

v[r, e, alpha, beta]
] ;

rrn = Norm[ rr ] ;
Piecewise[{{vPlusCubic[rrn,ee,a,b],rrn^2<ee^2}},-b/rrn]
] ;

inverseRadialSmoothedCap[{x, y}, ϵ, α, β] // TraditionalForm

Timing[Module[{e, a, b, range},
e = 0.5 ;
a = 1 ;
b = 1 ;
range = 5 ;
Plot3D[ 
inverseRadialSmoothedCap[{x, y}, e, a, b],
{x, -range, range}, 
{y, -range, range}, 
PlotRange -> {Full, Full, {-(a+b)/e-0.5, e}},
(*Ticks -> None,*)
(* http://mathematica.stackexchange.com/questions/13657/plot-showing-discontinuity-where-it-shouldnt *)
Exclusions->None
(*, Epilog->{
Point[ { 0,e, -b/e} ],
Point[ { 0, 0, -(b+a)/e} ],
Text["(0, ϵ, -(β/(Abs[r])))", {0, 2.5 e, -b/e} ],
Text["(0, 0, -((α + β)/(Abs[r])))", { 0,0.9, -(a + b)/e } ]
}*)
]
]]

(*Clear[ inverseRadialSmoothedPeriodicExtension ]
inverseRadialSmoothedPeriodicExtension[r_, a_, b_, e_, alpha_, beta_] =
Sum[ inverseRadialSmoothedCap[r + i a  + j b, e, alpha, beta], {i, -2, 2},{j, -2, 2}] ;
*)

(*Clear[ periodicProjectionToBaseUnitCell ] 
periodicProjectionToBaseUnitCell[ v_, am_, gm_ ] :=
Module[{a, b, g, h},
{a, b} = { am[[1]], am[[2]] } ;
{g, h} = Floor[{gm[[All, 1]], gm[[All, 2]]} . v / (2 Pi)] ;
v - (a g + b h )
]*)

(*Module[{a, b, v, am, gm},
a = { 1, 13 } ;
b = { 15, 3 } ;
v = { 75, 99 } ;

am = { a, b } ;
gm =2 Pi Inverse[ am ] ;

periodicProjectionToBaseUnitCell[ v, am, gm ]
]*)

ff[{x_, y_}] =Module[{a, b, am, gm, g1, g2, e, alpha, beta, fourierTerms, sampleV, sampleVunitParams,range, terms, p, maxFourierTerms, uOfx, vOfx, f},
a = { 5, 13 } ;
b = { 15, 3 } ;
(*a = { 20, 0 } ;
b = { 0, 20 } ;*)
am = { a, b } ;
gm =2 Pi Inverse[ am ] ;
g1 = gm[[All, 1]] ;
g2 = gm[[All, 2]] ;

uOfx[r_] = r . g1/(2 Pi) ;
vOfx[r_] = r . g2 /(2 Pi) ;

e = 0.5 ; 
alpha = 1 ;
beta = 1 ;
range = 25 ;
terms = 1 ;
maxFourierTerms = 2 ;

sampleV[{s_, t_}] = Sum[inverseRadialSmoothedCap[{s, t} -( i a  + j b), e, alpha, beta], {i, -terms, terms}, {j, -terms, terms}] ;
sampleVunitParams[{u_, v_}] = sampleV[ u a + v b ] ;

f[{xx_, yy_}] =FourierSeries[ sampleVunitParams[{xx, yy}], {xx, yy}, {maxFourierTerms, maxFourierTerms}, FourierParameters->{1,2Pi}] ;

{
sampleV[x, y], 
f[x, y], 
f[r, s] /. { r -> uOfx[{{x},{y}}], s -> vOfx[{{x},{y}}] } 
}
]

(*p = Plot3D[ /
{sampleV[{x, y}]},
{x, -range, range},
{y, -range, range}, 
PlotRange -> Full, 
PlotPoints -> 50,
Exclusions -> None ] ;*)

(*fourierTerms[  nFourierTerms_ ] := 
Table[
 NIntegrate[ 
sampleV[ uu a + vv b]  Exp[- 2 Pi I( i uu  + j vv)],
{uu, -1/2, 1/2},
{vv, -1/2, 1/2}
,Method->"ClenshawCurtisRule"
,WorkingPrecision->50
,AccuracyGoal->5 
],
{i, 0, nFourierTerms - 1},
{j, 0, nFourierTerms - 1 }
] ;

fv = fourierTerms[ maxFourierTerms ]  ;

fSeries[{r_, s_}, n_] := 
fv[[1]][[1]] 
+ 2 Re[ Sum[fv[[1]][[j+1]]Exp[ I  j g2. {0, s} ], {j, 1, n-1}]]  
+ 2 Re[ Sum[fv[[i+1]][[j+1]]Exp[ I (i g1 + j g2). {r, s} ], {i, 1, n-1}, {j, 1, n-1}]]  ;
(*fSeries[r, 4]*)*)

(* Switched from Row to GraphicsRow: 
http://mathematica.stackexchange.com/questions/2058/sizing-cells-in-a-graphicsgrid-graphicsrow *)
(*Manipulate[
,
{{n, 1}, 1, maxFourierTerms, 1},
{{m, 1}, 1, maxFourierTerms, 1},
SaveDefinitions->True
]*)


Clear[sampleV, fSeries]
sampleV[{x_, y_}] = ff[{x, y}][[1]] ;
fSeries[{x_, y_}] = ff[{x, y}][[3]] ;

(*Graphics*)(*Row[{*)

fSeries[{x, y}] // TraditionalForm
,Plot3D[ 
sampleV[{r, s}],
{r, -20, 20},
{s, -20, 20},
PlotRange -> Full,
Exclusions->None,
PlotPoints -> 50
]

Plot3D[ 
fSeries[{r, s}],
{r, -20, 20},
{s, -20, 20},
PlotRange -> Full,
Exclusions->None,
PlotPoints -> 50
]

(*,Plot3D[ 
2 Re[ fv[[n+1]] [[m+1]]Exp[ I (n g1 + m g2). {r, s} ]],
{r, -20, 20},
{s, -20, 20},
PlotRange -> {Full, Full,{-2, 2}},
PlotPoints -> 50,
Exclusions->None
]*)

(*Module[{n, t},
n = 3 ;
t = Table[
 f[i] g[j]
,
{i, 0, n},
{j, 0, n}
] ;
{t, t // MatrixForm}
] *)

scratch notes for text

(*HoldForm[(Subscript[a, rs]) = Integrate[ V[u,v] E^(- 2 Pi I(r u + s v)), {u, 0, 1}, {v, 0, 1}]] // TraditionalForm
HoldForm[ s(x, y) = (∑_(r,s)) (Subscript[a, rs])((E^(I r x+I s y))) = 2 Re((∑_(r,s > 0)) (Subscript[a, rs])(E^(I r x+I s y)) + (∑_(s > 0)) (Subscript[a, 0s])((E^(I s y))) )+ (Subscript[a, 00])] *)
(*HoldForm[ (∑_(r,s)) (Subscript[a, rs])E^(Ir x + I s y)  *)(*// TraditionalForm*)

(Subscript[a, rs])=(∫_0^1)(∫_0^1)V(u,v) (E^(-2 π I(r u+s v)))\[DifferentialD]v\[DifferentialD]u
s(x, y) = (∑_(r,s)) (Subscript[a, rs])((E^(I r x+I s y)))

Note on Fourier sums.  Our sum is of the form

s(x, y) = (∑_(r,s)) (Subscript[a, rs])((e^(i r x+i s y))),

where for real V(u,v)

(Subscript[a, rs])=(∫_0^1)(∫_0^1)V(u,v) (e^(-2 π i(r u+s v)))dvdu.

We can avoid evaluating the terms in the left hand plane, by writing:

s(x, y) = 2 Re((∑_(r,s > 0)) (Subscript[a, rs])(e^(i r x+i s y)) + (∑_(s > 0)) (Subscript[a, 0s])((e^(i s y))) )+ (Subscript[a, 00]),

Clear[ir]
ir = Compile[{{uParam, _Real}, {vParam, _Real}, {eParam, _Real}, {alphaParam, _Real}, {betaParam, _Real}},inverseRadialSmoothedCap[{uParam ,vParam } , eParam, alphaParam, betaParam]] 

Timing[Module[{e, a, b, range},
e = 0.5 ;
a = 1 ;
b = 1 ;
range = 5 ;
Plot3D[ 
ir[x, y, 0.5, 1, 1],
{x, -range, range}, 
{y, -range, range}, 
PlotRange -> {Full, Full, {-(a+b)/e-0.5, e}},
(*Ticks -> None,*)
(* http://mathematica.stackexchange.com/questions/13657/plot-showing-discontinuity-where-it-shouldnt *)
Exclusions->None
(*, Epilog->{
Point[ { 0,e, -b/e} ],
Point[ { 0, 0, -(b+a)/e} ],
Text["(0, ϵ, -(β/(Abs[r])))", {0, 2.5 e, -b/e} ],
Text["(0, 0, -((α + β)/(Abs[r])))", { 0,0.9, -(a + b)/e } ]
}*)
]
]]