# The Rayleigh-Benard Flow

eps=0.1;omega=6;
psi[x_, z_,t_] = Sin[Pi*x]*Sin[Pi*z]+eps*Cos[omega*t]*Cos[Pi*x]*Sin[Pi*z];
f[x_,z_,t_]=D[psi[x,z,t],z]; g[x_,z_,t_]=-D[psi[x,z,t],x];

sol1[a_, b_,tfinal_]:=NDSolve[{x'[t]==f[x[t],z[t],t],z'[t]==g[x[t],z[t],t],
x[0]==a, z[0]==b}, {x,z}, {t,0,tfinal}, MaxSteps->5000];
rightexit={};leftexit={};
Do[a=0.01+i/20; b=0.01;solution=sol1[a,b, 40];
aa=ToString[a]; bb=ToString[b];
llabel=StringJoin["(",aa,",",bb,")"];
Clear[c,average,x1, domain,range,tpositive,tnegative];
x1[t_]=First[x[t]-1.000001/. solution];
domain=Table[t,{t,0,40, 0.1}];
range=x1[domain];
testrange=Table[range[[i]]*range[[i+1]], {i, Length[range]-1}];
Catch[Do[index=i;
If[testrange[[i]]<0, Throw[i]]; , {i, Length[testrange]}]];
If[range[[index]]<0,tnegative=domain[[index]],tnegative=domain[[index+1]]];
If[range[[index]]<0,tpositive=domain[[index+1]],tpositive=domain[[index]]];
Do[average=(tpositive+tnegative)/2;c=x1[average];
If[c < 0, tnegative=average, tpositive=average],{i, 50}];
ExitToRight=average;
rightexit=Append[rightexit, {a, b, ExitToRight}];
x2[t_]=First[x[t]-0.00001/. solution];
Plot[x2[t], {t, 0, 40}, PlotPoints->500,
PlotLabel->StringJoin["x(t) with x(0) = ", aa]];
domain=Table[t,{t,0,40, 0.1}];
range=x2[domain];
testrange=Table[range[[i]]*range[[i+1]], {i, Length[range]-1}];Catch[Do[index=i;
If[testrange[[i]]<0, Throw[i]]; , {i, Length[testrange]}]];
If[range[[index]]<0,tnegative=domain[[index]],tnegative=domain[[index+1]]];
If[range[[index]]<0,tpositive=domain[[index+1]],tpositive=domain[[index]]];
Do[average=(tpositive+tnegative)/2;c=x2[average];
If[c < 0, tnegative=average, tpositive=average],{i, 50}];
ExitToLeft=average;
Print[llabel," exits to the cell on the right at t = ", ExitToRight];
Print[llabel," exits to the cell on the left at t = ", ExitToLeft];
leftexit=Append[leftexit, {a, b, ExitToLeft}],
{i, 1, 19}];
ExitTime=Table[Min[rightexit[[i,3]],leftexit[[i,3]]], {i,Length[leftexit]}];
graph1=ListPlot[ExitTime, PlotJoined->True, PlotRange->All, AxesOrigin->{0,0},
PlotLabel-> StringJoin["Exit time of particles located at z = ", bb]]