Reza Malek-Madani

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]]}];
go to Top