whittemore1
psi=If[x^2+y^2< 1, 1/2(x^2+y^2), -1/Sqrt[x^2+y^2]];
f[x_, y_] = D[psi, y]; g[x_, y_]=-D[psi, x];
NumOfPoints=100; (*Number of points on the boundary of parcels*)
radius = 0.31; (*Radius of each initial parcel*)
(* Initial Centers of individual parcels*)
centers={{0.3, 0.3},{0.9, 0.9},{1.8, 1.8}};
tfin=10; (*Total Time*)
interval = 0.2; (*Time between snapshots*)
diffeqn[tfinal_,aa_,bb_]:=NDSolve[{x'[t]==f[x[t],y[t]],
y'[t]==g[x[t],y[t]],
x[0]==aa, y[0]==bb},
{x,y}, {t,0,tfinal}, MaxSteps->10000];
Do[initdata[i]=N[Table[{centers[[i,1]]+radius*Cos[w],
centers[[i,2]]+radius*Sin[w]},
{w, 0, 2*Pi, 2*Pi/NumOfPoints}]];
oldsolution[i]=Table[diffeqn[tfin,initdata[i][[j,1]],
initdata[i][[j,2]]], {j, Length[initdata[i]]}];
data[i]=Table[{x[t], y[t]} /. oldsolution[i], {t, 0, tfin, interval}];
ColoredSnapshots[i]=Table[Graphics[{RGBColor[1-0.2*i, 0.2*i, 0],
Map[Polygon, Flatten[data[i][[j]],1], {0}]}], {j, Length[data[i]]}],
{i,Length[centers]}];
graph2=ParametricPlot[{Cos[t], Sin[t]}, {t, 0, 2Pi},
PlotStyle->{Thickness[0.01]}];
graph1=ContourPlot[psi, {x, -3, 3}, {y, -3, 3}, PlotPoints->250,
Contours->20];
Do[
giffile=Show[graph1,graph2,Table[ColoredSnapshots[k][[i]],
{k,Length[centers]}],
PlotRange->{{-3,3}, {-3,3}},
AspectRatio->Automatic];
llabel=StringJoin[ToString[i],".gif"];
Display[llabel,giffile,"GIF"],
{i, Length[ColoredSnapshots[1]]}];