clear all
close all
clc
% -------------------------------------------------------------

% -------------------------------------------------------------
% ucitavanje fileova
v=load('bavaniste.txt'); % merni podaci
S93=load('GE 3.6-137.txt'); % kriva snage
% proracun koeficijenta snage agregata, u tabeli je koef. pritiska dat
for Ht= 85:1:132
% podaci o agregatu
PN=3.6*10^6; %Nazivna snaga vetroturbine
A=14741; %porvsina koju zahvata vetroturbina
 % visina glavcine
nv=150+Ht; %nadmorska visina glavcine
% -------------------------------------------------------------

% -------------------------------------------------------------
% ucitavanje podataka o brzini vetra na 20 i 40m , temperaturi i smeru
v20=v(:,6);  %srednje desetominutne brzine vetra na 20m
v40=v(:,3);  %srednje desetominutne brzine vetra na 40m
T=v(:,8);    %temperatura
S=v(:,7);    % smer vetra
% usrednjavanje podataka o brzini vetra i prikaz rezultata
vsr20=mean(v20);
vsr40=mean(v40);
disp('Srednja godisnja brzina vetra na visini od 20 metara: ');
disp(vsr20);
disp('Srednja godisnja brzina vetra na visini od 40 metara: ');
disp(vsr40);
%proracun srednje godisnje vrednosti koeficijenta smicanja vetra
alfasr=log(vsr40/vsr20)/log(60/10);
disp('Srednja godisnja vrednost koeficijenta smicanja vetra:');
disp((alfasr));
% -------------------------------------------------------------



% -------------------------------------------------------------
% relativna cestina javljanja brzina vetra po smerovima
sm=zeros(1,360);
for i=0:1:359
    for k=1:1:length(S)
        if S(k)==i
            sm(i+1)=sm(i+1)+1;
        end
    end
end
sm=sm*100/length(S);
% -------------------------------------------------------------

% -------------------------------------------------------------
%proracun brzine vetra na visini turbine
ln20=log(10);
ln40=log(60);
lnHt=log(Ht);
% ekstrapolacija brzina vetra na visinu osovine vetroturbine i proracun
% gustine vazduha
for i=1:length(v40)
    if v40(i)>0
        if v40(i)-v20(i)<=0
            vHt(i)=v40(i);
            ka=exp(-0.000118*nv);
            kt(i)=288.15/(273.15+T(i));
            ro(i)=1.225*ka*kt(i);
            pHt(i)=0.5*vHt(i)^3*ro(i);
        else
             lnz0(i)=(v40(i)*ln20-v20(i)*ln40)/(v40(i)-v20(i));
             vHt(i)=(lnHt-lnz0(i))*v40(i)/(ln40-lnz0(i));
             ka=exp(-0.000118*nv);
             kt(i)=288.15/(273.15+T(i));
             ro(i)=1.225*ka*kt(i);
             if  vHt(i)>25
                 vHt(i)=0;
                 ka=exp(-0.000118*nv);
                 kt(i)=288.15/(273.15+T(i));
                 ro(i)=1.225*ka*kt(i);
             end
             pHt(i)=0.5*vHt(i)^3*ro(i);
        end
    else
        ka=exp(-0.000118*nv);
        kt(i)=288.15/(273.15+T(i));
        ro(i)=1.225*ka*kt(i);
        vHt(i)=0;
        pHt(i)=0;
        lnz0(i)=0.00001;
    end
end

%srednja godisnja brzina vetra na visini turbine
vHtsr=mean(vHt);
disp('Srednja godisnja brzina vetra na visini turbine:');
disp((vHtsr));
%srednja godisnja gustina vazduha na visini osovine vetroturbine
gustinavazduha=mean(ro);
disp('Srednja godisnja gustina vazduha na visini turbine:');
disp((gustinavazduha));
%srednja godisnja gustina snaga vetra na Htm
pHtsr=mean(pHt);
disp('Srednja godisnja gustina snage vetra na visini turbine:');
disp((pHtsr));
% -------------------------------------------------------------

% -------------------------------------------------------------
% graficki prikaz diskretnog histograma brzina
vHtr=(vHt);
for i=1:max(vHtr)
    vHts(i)=length(find(vHtr==i))/(6*8760);
end


%Proracun parametara Weibullove statistike
Epf=mean(vHt.^3)/(mean(vHt))^3;
k=1+(3.69/Epf^2)
c=mean(vHt)/gamma(1+(1/k))
%Graficki prikaz Weibullove funkcije gustine raspodele brzina vetra
hold on
fplot(@(v) (k/c)*((v/c)^(k-1))*exp(-(v/c)^k),[0 25],'r')

%Proracun gustine snage na osnovu Weibullove statistike
pHtsrweibull=0.5*gustinavazduha*c^3*gamma(1+(3/k));
disp('Srednja godisnja gustina snage vetra na visini turbine, dobijena pomocu Weibullove raspodele:');
disp((pHtsrweibull));
% -------------------------------------------------------------

% -------------------------------------------------------------
%proracun elektricne snage generisanja za odabranu vetroturbinu
vHtx=ceil(vHt)-1;
for i=1:length(vHtx)-1*0
    if vHtx(i)>1
    cpHt(i)=S93(vHtx(i),3)+(vHt(i)-vHtx(i))*(S93(vHtx(i)+1,3)-S93(vHtx(i),3)); %linearna interpolacija faktora iskoriscenja vetroturbine
    else
      cpHt(i)=0;
    end
    PHtv(i)=0.5*ro(i)*vHt(i)^3*A; %proracun snage vjetra koju zahvata vetroturbina
    PHt(i)=cpHt(i)*PHtv(i); %proracun elektricne snage vetrogeneratora sa odabranom vetroturbinom
    PN=3.6*10^6; %Nazivna snaga vetroturbine
    if PHt(i)>PN
    PHt(i)=PN;
    end
end

%proracun ukupne godisnje generisane elektricne energije [MWh]
W=sum(PHt)/6000000;
disp('Ukupna godisnja proizvodnja elektricne energije:');
disp((W)); % u MWh
%proracun faktora iskoriscenja kapaciteta
CPF=W/(3.6*8760);
disp('CF:');
disp((CPF*100));
% -------------------------------------------------------------
Wge(Ht)=W;
CFge(Ht)=CPF;

%ekonomska analiza
ICC85=1.7*PN/(10^6); % pocetni investicioni troskovi
Itot=ICC85*(1+0.095*((Ht-85)/85));%promena sa visinom
ICC(Ht)=Itot;

n=20; % period amortizacije
m=0.01; % operativni troskovi u Euro/kWh
A=0.97;  % raspolozivost vetroelektrane
Fid_in_tariff=92; %fid in tarifa po Mwh
j=0;
for i=0.0:0.001:0.2
    j=j+1;
    x(j)=i;
    a(j)=(i*(1+i)^n)/((1+i)^n-1);
    c(j)=m+(a(j)*Itot*1000000/(A*W*1000));
end
Ctr(Ht)=c(61);
figure1 = figure('Color',[1 1 1]);
axes('Parent',figure1,'YGrid','on','XGrid','on','FontWeight','bold',...
    'FontSize',12);
box('on');
hold('all');
plot(x*100,c*1000,'LineWidth',2,'Color',[0.4784 0.06275 0.8941]);
title('Ekonomska Analiza');
xlabel(' Interesna stopa [%]','FontWeight','demi','FontSize',14);
ylabel('Troskovi proizvodnje [Euro/MWh]','FontWeight','bold','FontSize',14);
hold on
fid=(Fid_in_tariff+0.0001)*ones(1,length(x));
plot(x*100,fid,'LineWidth',2,'LineStyle','--','Color',[1 0 0]);
% -------------------------------------------------------------

end