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 = -10;
Gamma2 = 8;

U = 4;
a = 2;

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

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

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

%time step
ts = 0.000001;

%Initial vortex center location
yc1 = -1;
xc1= -1;
yc2 = 1;
xc2 = 1;

%Initialize zero matrices
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);
u2 = zeros(ysize,xsize);
v2 = zeros(ysize,xsize);
u3 = zeros(ysize,xsize);
v3 = zeros(ysize,xsize);

dPsidy1 = zeros(ysize,xsize);
dPsidx1 = zeros(ysize,xsize);
dPsidy2 = zeros(ysize,xsize);
dPsidx2 = zeros(ysize,xsize);
dPsidy3 = zeros(ysize,xsize);
dPsidx3 = zeros(ysize,xsize);


%%%%%%%%%%%%Calculate velocity caused by cross flow%%%%%%%%%%%%%%
Psi3 = U*((-sind(a).*x) + (cosd(a).*y));
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
        dPsidy3(i,j) = (Psi3(i+1,j) - Psi3(i-1,j) )/ (2*spacing);
        dPsidx3(i,j) = (Psi3(i,j+1) - Psi3(i,j-1)) / (2*spacing);
    end
end

u3 = dPsidy3;   %matrix of x velocity
v3 = -1* dPsidx3; %matrix of y velocity

u3(u1>1000000) = 0;
u3(u1<-1000000) = 0;
v3(v1>1000000) = 0;
v3(v1<-1000000) = 0;
%[%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

        %Psi defined for vortices
    Psi1 =(Gamma1*log(r1)./(2*pi));
    Psi2 =(Gamma2*log(r2)./(2*pi));

        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

        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
        u2 = dPsidy2;   %matrix of x velocity
        v2 = -1* dPsidx2; %matrix of y velocity

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

        u2(u2>1000000) = 0;
        u2(u2<-1000000) = 0;
        v2(v2>1000000) = 0;
        v2(v2<-1000000) = 0;


        xc1 = xc1 + (u2(n1, m1).*ts) + (u3(n1, m1).*ts);
        yc1 = yc1 + (v2(n1, m1).*ts) + (v3(n1, m1).*ts);
        xc2 = xc2 + (u1(n1, m1).*ts) + (u3(n1, m1).*ts);
        yc2 = yc2 + (v1(n1, m1).*ts) + (v3(n1, m1).*ts);


        xloc1(1, m1) = xc1;
        yloc1(1, m1) = yc1;
        xloc2(1, m1) = xc2;
        yloc2(1, m1) = yc2;



%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

figure;

plot(xloc1, yloc1);
hold on;
plot(xloc2, yloc2);
ylim([-5 5]);
xlim([-5 5]);