clearvars
clc
Import = importdata('tideHeight1s.xlsx');

Times=Import.data(1:end,1);
TideH=Import.data(1:end,2);

[iterT, redundant1] = size(Times); %Iteration limit, predefined array size

Turbineshaft=7.2; %meters
TurbineArea=pi()*(Turbineshaft/2)^2;
LagoonArea=11.5*(10^6); %m^2
DischargeCoef=0.8;
TurbineCoef=0.5;
Gravconst=9.807; %m/s^2
seawaterdensity=1029; %kg/m^3
maxHdiff=7; %Max difference of lagoon and sea water level to open sluice, max makes sense is 7
%(difference between max sea level and turbine height) and lowering through
%optimisation to limit variable size
maxturbine=50; %Number max turbine shafts for iteration, array size (based on largest sensible value)
%lowered through optimisation to limit variable size
heightdif=[0.1:0.1:maxHdiff]; %Set of sluice height for Discharge Calc and Optimisation iteration
[redundant2, iterH]=size(heightdif); %Iteration limit, array size

%Criteria for start max water level
I=find(TideH == max(TideH(:)),1);

H=zeros(iterT,iterH,maxturbine);
Q=zeros(iterT,iterH,maxturbine);
V=zeros(iterT,iterH,maxturbine);
P=zeros(iterT,iterH,maxturbine);
optimalmodeP=zeros(iterH,maxturbine);
optimalsumP=zeros(iterH,maxturbine);
optimalsumdevP=zeros(iterH,maxturbine);
optimalmodedevP=zeros(iterH,maxturbine);

%Loop for turbine number
%Loop for a range of height differences for sluice opening
%Discharge calculation loop

for x=1:maxturbine   %number of turbines
    y=1;
    for y=1:iterH %
        %Loop calculation of Discharge and Lagoon water level (H)
        z=1;
        for z=1:iterT
            if z<=I %Set maximum Lagoon water level
            H(z,y,x)=TideH(z);
            Q(z,y,x)=0;
            elseif z>I && (TideH(I)-TideH(z)) < heightdif(y)   %Until open sluice    
            H(z,y,x)=H(z-1,y,x);
            Q(z,y,x)=0;
            elseif z>I && H(z-1,y,x) > Turbineshaft %Calculate Discharge and Lagoon water level, per x turbines
            Q(z,y,x)=DischargeCoef*x*TurbineArea*sqrt(2*Gravconst*(H(z-1,y,x)-TideH(z-1)));    
            H(z,y,x)=H(z-1,y,x)-(Q(z,y,x)/LagoonArea);
            elseif z>I && H(z,y,x)<=TideH(z) && H(z,y,x)>Turbineshaft %If lagoon water falls too quickly, match sea level
            Q(z,y,x)=DischargeCoef*x*TurbineArea*sqrt(2*Gravconst*(H(z-1,y,x)-TideH(z-1)));    
            H(z,y,x)=TideH(z);    
            elseif z>I && H(z,y,x)<=Turbineshaft %Point where sea level falls to height of shaft
            Q(z,y,x)=0;
            H(z,y,x)=H(z-1,y,x);  
            end
        end
    end
end

%Loop for calculating velocity PER TURBINE (because x(V^3) does not equal
%(xV)^3
for d=1:maxturbine   %number of turbines
    
V(:,:,d)=Q(:,:,d)/(TurbineArea*d);
P(:,:,d)=(TurbineCoef*seawaterdensity*TurbineArea*d*(V(:,:,d).^3))/2;
end

%Optimising solutions for power, but what is optimal?


for n=1:maxturbine
   for m=1:iterH
   S=P(:,m,n);
   time=nnz(S); %Total time of non-zero power
   S=nonzeros(S);
   sumP=sum(S); %Sum over all elements, since 1 step is 1s, it's an integral
   optimalsumP(m,n)=sumP*time; %Weighed for duration time
   end
end

optimalsumP=rmmissing(optimalsumP); %remove NaNs
O = max(max(optimalsumP));
[X, Y] = find(optimalsumP==O,1);
%Cheching for the largest weighed integral
%Because the above loop matches values of main loop, X and Y are indices of
%desired data column

for n=1:maxturbine
   for m=1:iterH
   S=P(:,m,n);
   time=nnz(S); %Total time of non-zero power
   S=nonzeros(S);
   modeP=mode(S); %Mode is which element is most common
   optimalmodeP(m,n)=modeP*time; 
   %Looking for highest intersection of mode and time
   end
end

optimalmodeP=rmmissing(optimalmodeP);
O = max(max(optimalmodeP));
[M, N] = find(optimalmodeP==O,1);

for n=1:maxturbine
   for m=1:iterH
   S=P(:,m,n);
   time=nnz(S); %Total time of non-zero power
   S=nonzeros(S);
   sumP=sum(S); %Sum over all elements
   devP=mad(S); % Max deviation from average
   optimalsumdevP(m,n)=sumP/devP*time; %Integral wighed for time and deviation
   end
end

optimalsumdevP=rmmissing(optimalsumdevP);
O = max(max(optimalsumdevP));
[R, T] = find(optimalsumdevP==O,1);


for n=1:maxturbine
   for m=1:iterH
   S=P(:,m,n);
   time=nnz(S); %Total time of non-zero power
   S=nonzeros(S);
   devP=mad(S);
   modeP=mode(S); %Mode is which element is most common
   optimalmodedevP(m,n)=modeP/devP*time; 
   %Intersection of mode, deviation and time
   end
end

optimalmodedevP=rmmissing(optimalmodedevP);
O = max(max(optimalmodedevP));
[W, E] = find(optimalmodedevP==O);

%Just to see the values
optimalsumP(X,Y)
optimalmodeP(M,N)
optimalsumdevP(R,T)
optimalmodedevP(W,E)

Timeh=Times/3600;
        
close all
figure (1)
xlabel('Time (hours)')
ylabel('Water level (m)')
grid on
hold on
plot(Timeh, TideH)
plot(Timeh, H(:,X,Y))
plot(Timeh, H(:,M,N))
plot(Timeh, H(:,W,E))
plot(Timeh, H(:,R,T))
legend('Sea', 'Sum P', 'Mode P', 'Sum/Dev P', 'Mode/Dev P' )

Zero=zeros(1,21601);

figure (2)
xlabel('Time (hours)')
ylabel('Power (MW)')
grid on
hold on
plot(Timeh, Zero)
plot(Timeh, P(:,X,Y)/10^9)
plot(Timeh, P(:,M,N)/10^9)
plot(Timeh, P(:,W,E)/10^9)
plot(Timeh, P(:,R,T)/10^9)
legend('Zero', 'Sum P', 'Mode P', 'Sum/Dev P', 'Mode/Dev P')

M10=M/10;
X10=X/10;
W10=W/10;
R10=R/10;
fprintf('Optimal height to open sluice for Sum is %1.1f', X10)
fprintf('meters')
fprintf('\n')
fprintf('Optimal number of turbines for Sum is %d', Y)
fprintf('\n')
fprintf('Optimal height to open sluice for Mode is %1.1f', M10)
fprintf('meters')
fprintf('\n')
fprintf('Optimal number of turbines for Mode is %d', N)
fprintf('\n')
fprintf('Optimal height to open sluice for Sum/Dev is %1.1f', W10)
fprintf('meters')
fprintf('\n')
fprintf('Optimal number of turbines for Sum/Dev is %d', E)
fprintf('\n')
fprintf('Optimal height to open sluice for Mode/Dev is %1.1f', R10)
fprintf('meters')
fprintf('\n')
fprintf('Optimal number of turbines for Mode/Dev is %d', T)
fprintf('\n')