Vortex-Project/Plaqutte_Vortex_Counter.m at testing · mcawte/Vortex-Project · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
% Plaquette Vortex Counter
% This needs to calculate the phase around each lattice point. If the phase
% winds by 2pi, there is a vortex there.

 clear all
% close
 % clc

 Stirred_Vortices_84 = load('Data/Stirred_Vortices_84');

 PSI = Stirred_Vortices_84.PSI;

 PhaseofPSI = angle(PSI);

 [X,Y] = size(PSI);

 Vortex_Grid = zeros(X,Y);

tic

 for ii = 1:X-1
   for jj = 1:Y-1

       Alpha1 = PhaseofPSI(ii,jj);
       Beta1 = PhaseofPSI(ii,jj+1);

       m = 0;

       if Beta1 - Alpha1 > pi;
          m = m -1;
       elseif Beta1 - Alpha1 < -pi;
          m = m + 1;
       end

       Alpha2 = PhaseofPSI(ii,jj+1);
       Beta2 = PhaseofPSI(ii+1,jj+1);

         if Beta2 - Alpha2 > pi;
          m = m -1;
       elseif Beta2 - Alpha2 < -pi;
          m = m + 1;
       end


       Alpha3 = PhaseofPSI(ii+1,jj+1);
       Beta3 = PhaseofPSI(ii+1,jj);

         if Beta3 - Alpha3 > pi;
          m = m -1;
       elseif Beta3 - Alpha3 < -pi;
          m = m + 1;
       end


       Alpha4 = PhaseofPSI(ii+1,jj);
       Beta4 = PhaseofPSI(ii,jj);

         if Beta4 - Alpha4 > pi;
          m = m -1;
       elseif Beta4 - Alpha4 < -pi;
          m = m + 1;
       end


       Vortex_Grid(ii,jj) = m;

   end
 end

 Plaquette_Time = toc
 %Plaquette_cpu = cputime

 figure
 imagesc(X,Y,Vortex_Grid)
 set(gca,'ydir','normal')