% gives the distribution of times spent in each min.
function thermal_dist(e,g, v0, deltaE,len)

% square this
%e = (4*kT/gamma) *delta_t;
% g = 1/4kT

x = zeros(len,1);
y = zeros(len,1);
x(1) = sqrt(e/2)*randn();
y(1) = sqrt(e/2)*randn();
xmin = -0.2;
xmax = 0.7;
deltax = xmax-xmin;
ymin = xmin;
ymax = xmax;
deltay=ymax-ymin;

   for i=2:len
   x(i) =x(i-1) -e*g*dVdx(x(i-1),y(i-1),v0,deltaE) + sqrt(e/2)*randn() ;
   if(x(i)>xmax) 
       %s = floor(x(i)/deltax);
       x(i) = 2*xmax - x(i);
   elseif (x(i)<xmin)
       x(i) = 2*xmin - x(i); 
   end
       
       y(i) = y(i-1) -e*g*dVdy(x(i-1),y(i-1),v0,deltaE) + sqrt(e/2)*randn();
   
   if(y(i)>ymax) 
       %s = floor(x(i)/deltax);
       y(i) = 2*ymax - y(i);
   elseif(y(i)<ymin)
       y(i) = 2*ymin - y(i); 
   end
   
   
   %x(i) = x(i-1) + e*g*x(i-1)*(1-x(i-1)^2) + sqrt(e/2)*randn();
   end
   %o = occupancy(x,100);
   figure(1);
   %imagesc(x,y,,100)
   hist2d(x,y,100);
   figure(2);
   plot(x)
   figure(3);
   plot(y)
   Es = V(x,y,v0,deltaE);
   figure(4);
   hist(Es,1000)
   a = linspace(xmin,xmax,100);
   b = linspace(ymin,ymax,100);
   f = zeros(100);
   for i=1:100
      f(:,i) = V(a,b(i),v0,deltaE);
   end
   figure(5);
   imagesc(a,b,f);
end

function E = V(x,y, v0, de)

%E = x.^2.*( e1*(x-1).^2 + e2*(x+1).^2)/4;
%E = v0*(1 - x.^2-y.^2).^2 +0.75*de*(x - (x.^3)./3);
E = v0*cos(2*pi*x).*cos(2*pi*y) + de*3*(x-0.25).*(1 - 16*(x-0.25).^2/3);
end

function dEdx = dVdx(x,y,v0,de)

%dEdx = 0.5*x.*(e1*(x-1).^2 + e2*(x+1).^2) + 0.5*x.^2.*( e1*(x-1) + e2*(x+1));
%dEdx = -4*v0*x.*(1-x.^2) + 0.75*de*(1 - x.^2);
dEdx = -v0*2*pi*sin(2*pi*x).*cos(2*pi*y) + de*3*(1 - 16*(x-0.25).^2);
end

function dEdy = dVdy(x,y,v0,de)

%dEdx = 0.5*x.*(e1*(x-1).^2 + e2*(x+1).^2) + 0.5*x.^2.*( e1*(x-1) + e2*(x+1));
%dEdx = -4*v0*x.*(1-x.^2) + 0.75*de*(1 - x.^2);
dEdy = -v0*2*pi*sin(2*pi*y).*cos(2*pi*x);
end





%plot(traj)
