clearvars
close all

OP = [.25 .25 7 7,];

xmin = -5;
xmax = 5;
ymin = -5;
ymax = 5;
spacing = 0.01;

%Define value for vortices and cross flow
Gamma1 = -20;
Gamma2 = 8;

U = 5;
a = 0;

x = (xmin:spacing:xmax);
y = (ymin:spacing:ymax);

xsize = numel(x);
ysize = numel(y);

[x,y] = meshgrid(x,y);

%time step
ts = 0.00001;
dur = 0;

%Vortex center location
yc1 = -1;
xc1= -1;
yc2 = 1;
xc2 = 1;


r1 = sqrt(((x-xc1).^2)+((y-yc1).^2));
r2 = sqrt(((x-xc2).^2)+((y-yc2).^2));

%Psi defined for vortices and cross flow
Psi1= (Gamma1*log(r1)./(2*pi));

Psi2 = Gamma2*log(r2)./(2*pi);

Psi3 = U*((-sind(a).*x) + (cosd(a).*y));

xloc1 = zeros(1,xsize); %x location of vortex 1
yloc1 = zeros(1,ysize); %y location of vortex 1

xloc2 = zeros(1,xsize); %x location of vortex 2
yloc2 = zeros(1,ysize); %y location of vortex 2



u1 = zeros(ysize,xsize);
v1 = zeros(ysize,xsize);

dPsidy1 = zeros(ysize,xsize);
dPsidx1 = zeros(ysize,xsize);

for i = 2:(ysize-1)
    for j = 2:xsize-1
        %The y value changes as you increase the row (so i)
        %Essentially this is delta y over delta x
        dPsidy1(i,j) = (Psi1(i+1,j) - Psi1(i-1,j) )/ (2*spacing);
        dPsidx1(i,j) = (Psi1(i,j+1) - Psi1(i,j-1)) / (2*spacing);
    end
end

dPsidy2 = zeros(ysize,xsize);
dPsidx2 = zeros(ysize,xsize);

for i = 2:(ysize-1)
    for j = 2:xsize-1
        %The y value changes as you increase the row (so i)
        %Essentially this is delta y over delta x
        dPsidy2(i,j) = (Psi2(i+1,j) - Psi2(i-1,j) )/ (2*spacing);
        dPsidx2(i,j) = (Psi2(i,j+1) - Psi2(i,j-1)) / (2*spacing);
    end
end

u1 = dPsidy1;   %matrix of x velocity
v1 = -1* dPsidx1; %matrix of y velocity

u1(u1>1000000) = 0;
u1(u1<-1000000) = 0;
v1(v1>1000000) = 0;
v1(v1<-1000000) = 0;

for n1 = 1:(ysize)
    for m1 = 1:(xsize)
        xc2 = xc2 + (u1(n1, m1).*ts);
        yc2 = yc2 + (v1(n1, m1).*ts);
    end
end
%totalv = (u.^2 + v.^2).^.5;

%hold on
%contourf(x,y,Psi, 1000, 'LineStyle', 'none');
%colormap
%colorbar
%h = colorbar;
%title(h, '\psi')
%hold on;


%ph = streamslice(x,y,u1,v1);

%set(ph,'color','k');

%axis([xmin xmax ymin ymax]);
%set(gcf,'Units','inches','OuterPosition',OP);

%title('Streamlines for vortices');
%xlabel('x location');
%ylabel('y location');
%hold off