%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Berakningsprogram för solcellsanlaggningar
% 
% Syftet är att beräkna effekt (per sträng/tracker/växelriktare) och
% ström (per sträng och tracker) för godtycklig anläggning (norra halvklotet) vid godtycklig
% tidpunkt.
% 
% Varje sträng är beräknad som en punkt, (samt molnfritt väder förutsätts, steg #1), har även inplementerat
% globalstrålningshistorik från SMHI för Göteborg för att kunna prediktera elproduktion.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
clc;
% Horisonthöjden för respektive sträng. N=0, E=90, S=180, W=270
% Befintlig
Horisonthojd=[
 0	10	20	30	40	50	60	70	80	90 100 110 120 130	140	150	160	170	180	190	200	210	220	230	240	250	260	270	280	290	300	310	320	330	340	350	360
30	30	30	30	30	30	22	22	22	22	22	22	 7	 6	  5	  5	 35	 35	 12	 12	 12	 12	 20	 20	 38	 38	 38	 35	 35	 35	 35	 35	 35	 35	 35	 35	 30
30	30	30	30	30	30	20	20	 7	 7	 7	 7	 7	 6	  5	  5	 35	 35	 12	 12	 12	 12	 20	 20	 36	 36  36	 35	 35	 35	 35	 35	 35	 35	 35	 35	 30
30	30	30	30	30	30	20	20	 7	 7	 7	 7 	 7	 6	  5	  5	 35	 35	 12	 12	 12	 12	 20	 20	 35	 35	 38	 35	 35	 35	 35	 35	 35	 35	 35	 35	 30
30	30	30	30	30	30	30	20	20	 7	 7	 7	 7	 7	  7	  7	  7	  6	  5	  5	 35	 35	 12	 12	 12	 12	 20	 20	 30	 10	 10	 15	 35	 35	 35	 35	 30
30	30	30	60	60	60	60	60	60	60	60	60	60	 5	  5 	5	  5  	5  25	 25	 25	  8	  8	 20	 20	 20	 20	 20	 20	 15	 15	 35	 35	 35	 35	 35	 30
30	30	30	30	30	30	30	20	20	 7	 7	 7	 7	 6	  5	  5	 35	 35	 12	 12	 12	 12	 20	 20	 35	 35	 38	 35	 35	 35	 35	 35	 35	 35	 35	 35	 30];

% Stora eken och träd mot grannar i syd är borta
%Horisonthojd=[
% 0	10	20	30	40	50	60	70	80	90 100 110 120 130	140	150	160	170	180	190	200	210	220	230	240	250	260	270	280	290	300	310	320	330	340	350	360
%30	30	30	30	30	30	22	22	22	22	22	22	 7	 6	  5	  5	 12	 12	 12	 12	 12	 12	 20	 20	 38	 38	 38	 35	 35	 35	 35	 35	 35	 35	 35	 35	 30
%30	30	30	30	30	30	20	20	 7	 7	 7	 7	 7	 6	  5	  5	 12	 12	 12	 12	 12	 12	 20	 20	 36	 36  36	 35	 35	 35	 35	 35	 35	 35	 35	 35	 30
%30	30	30	30	30	30	20	20	 7	 7	 7	 7 	 7	 6	  5	  5	 12	 12	 12	 12	 12	 12	 20	 20	 35	 35	 38	 35	 35	 35	 35	 35	 35	 35	 35	 35	 30
%30	30	30	30	30	30	30	20	20	 7	 7	 7	 7	 7	  7	  7	  7	  6	  5	  5	 12	 12	 12	 12	 12	 12	 20	 20	 30	 10	 10	 15	 35	 35	 35	 35	 30
%30	30	30	60	60	60	60	60	60	60	60	60	60	 5	  5 	5	  5  	5  25	 25	 25	  8	  8	 20	 20	 20	 20	 20	 20	 15	 15	 35	 35	 35	 35	 35	 30
%30	30	30	30	30	30	30	20	20	 7	 7	 7	 7	 6	  5	  5	 12	 12	 12	 12	 12	 12	 20	 20	 35	 35	 38	 35	 35	 35	 35	 35	 35	 35	 35	 35	 30];

% Ange for respektive sträng: sträng#, tracker#, azimuth (orientering N=0, E=90, S=180, W=270), lutning (horisontalplan=0), antal paneler
Strang=[
  1   2  3   4   5 6    % Sträng
  1   1  2   2   2 3    % Tracker#, får inte vara 0
140 140 73 253 210 0    % Azimuth. N=0, E=90, S=180, W=270
 27  27 27  27  60 0    % Lutning, Horisontalplan=0
 22  22 18  18   0 0];  % Antal paneler, om inga paneler finnes skall det vara 0


% Vilken ort: longitud, latitud
% Här är det Kungälv
Ort=[11.83 57.85];
 
% Data solpaneler
Panel_b=1.649;
Panel_l=0.991;
Panel_p=290;
Pkoeff=-0.0039;
Mppt_i=9.2;
Panel_t=45;
Ref_t=25;

% Solinstrålningsparametrar for modell
Rjord=40000/2/pi;           % Radie jorden, km
Solin_dir=960;              % Direkt solinstralning, W/m2
Hatm=200;                   % Atmosfärshöjd i modellen, km
Atm_abs=0.25;               % Atmosfärsabsorbtion W/km
Solin_diff_rel=0.15;        % Solinstrålning diffus relativt

% SMHI-data för göteborg
% Predikterad globalstrålning i (kWh/m2 per månad) linjärt extrapolerad till 2018
Soldata=[
 11.5 
 24.7 
 83.3 
 123.9 
 156.3 
 172.2 
 161.5 
 125.7 
 88.5 
 39.5 
 13.6 
 7.4 
]';

% Inverterdata
Pmax_inv_DC=26300;
Pmax_inv_AC=17500;
Imax_track=[33 27 0];
Imax_track_1_2=51;

% Tidpunkt 
Year=2018;
%Month=6;
%Day=29;

% Loop for dagnummer, hela året, 365 dagar
%Dagnummer=datenum(Year,Month,Day)-datenum(Year,1,0);
Dagar_i_ar=365;
n=1;
o=1:Dagar_i_ar;
for n=1:Dagar_i_ar;
Dagnummer=o(n);
J_prim=360*Dagnummer/Dagar_i_ar;

% Loop for timmar
k=1;
l=0:0.5:23.5;
%l=0:0.5:24;
m=size(l);
for k=1:m(2);


%Local_time=12;
Local_time=l(k);
Rel_UTC=2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%        Ingen mer indata nedan               %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Var solen befinner sig
%Dagnummer=datenum(Year,Month,Day)-datenum(Year,1,0);
%Dagar_i_ar=datenum(Year,12,31)-datenum(Year,0,0);
%J_prim=360*Dagnummer/Dagar_i_ar;
Soldeklinationen=(0.3948-23.2559*cos((J_prim+9.1)*pi/180)-0.3915*cos((2*J_prim+5.4)*pi/180)-0.1764*cos((3*J_prim+26)*pi/180));
Tidsekvation=(0.0066+7.3525*cos((J_prim+85.9)*pi/180)+9.9359*cos((2*J_prim+108.9)*pi/180)+0.3387*cos((3*J_prim+105.2)*pi/180))/60;
MOZ=Local_time-Rel_UTC+Ort(1)/15;
WOZ=MOZ+Tidsekvation;
Timvinkeln=(12-WOZ)*15;
Solhojden=asin(cos(Timvinkeln*pi/180)*cos(Ort(2)*pi/180)*cos(Soldeklinationen*pi/180)+sin(Ort(2)*pi/180)*sin(Soldeklinationen*pi/180))*180/pi;
Solens_azimuth=ifelse (WOZ<12,
  180-acos((sin(Solhojden*pi/180)*sin(Ort(2)*pi/180)-sin(Soldeklinationen*pi/180))/(cos(Solhojden*pi/180)*cos(Ort(2)*pi/180)))*180/pi,
  180+acos((sin(Solhojden*pi/180)*sin(Ort(2)*pi/180)-sin(Soldeklinationen*pi/180))/(cos(Solhojden*pi/180)*cos(Ort(2)*pi/180)))*180/pi);
Solhojden_inv=90-Solhojden;

% Aktuell horisonthöjd, dvs skugga eller ej
% Solinstrålning
Satm=(Rjord+Hatm)*sin((180-(Solhojden+90)-asin(sin((Solhojden+90)*pi/180)*Rjord/(Rjord+Hatm))*180/pi)*pi/180)/sin((Solhojden+90)*pi/180);
Solin_dir_atm=Solin_dir-Satm*Atm_abs;
Panelata=Panel_p/(Panel_b*Panel_l*1000);
Solin_hor=ifelse (Solhojden<0,0,Solin_dir_atm*sin(Solhojden*pi/180));
Solin_diff_hor=Solin_hor*Solin_diff_rel;

% Loop for alla strängar
i=1;
j=size(Horisonthojd(:,1));
for i=1:(j(1)-1)
  Azimuth_diff(i)=Solens_azimuth-Strang(3,i);
  Hjalpvariabel_cos(i)=acos(cos(Strang(4,i)*pi/180)*cos(Solhojden_inv*pi/180)+sin(Strang(4,i)*pi/180)*sin(Solhojden_inv*pi/180)*cos(Azimuth_diff(i)*pi/180))*180/pi;
  Solin_glob(i)=Solin_dir_atm*cos(Hjalpvariabel_cos(i)*pi/180);
  Hor_interpol(i)=interp1(Horisonthojd(1,:),Horisonthojd(i+1,:),Solens_azimuth);
  Solin_glob_hor(i)=ifelse (Solhojden>Hor_interpol(i),Solin_glob(i),0);
  Solin_diff_hor(i)=ifelse (Solhojden<0,0,Solin_dir_atm*sin(Solhojden*pi/180)*Solin_diff_rel);
  Pdir(i)=Solin_glob_hor(i)*Panelata*Strang(5,i)*Panel_b*Panel_l*(1+Pkoeff*(Panel_t-Ref_t));
  Pdiff(i)=Solin_diff_hor(i)*Panelata*Strang(5,i)*Panel_b*Panel_l*(1+Pkoeff*(Panel_t-Ref_t));
  Ptot(i)=ifelse (Pdir(i)>Pdiff(i),Pdir(i),Pdiff(i));
  Idir(i)=Solin_glob_hor(i)*Mppt_i/1000;
  Idiff(i)=Solin_diff_hor(i)*Mppt_i/1000;
  Itot(i)=Idir(i)+Idiff(i);
  Icorr(i)=ifelse (Strang(5,i)>0,1,0);
  Trackertemp(i,Strang(2,i),1)=Ptot(i);
  Trackertemp(i,Strang(2,i),2)=Itot(i)*Icorr(i);
    i=i+1;
endfor

Tracker=Trackertemp(1,:,:)+Trackertemp(2,:,:)+Trackertemp(3,:,:)+Trackertemp(4,:,:)+Trackertemp(5,:,:)+Trackertemp(6,:,:);
Pmax_Inverter=Tracker(1,1)+Tracker(1,2);
Plim_Inverter=ifelse (Pmax_Inverter>Pmax_inv_AC,Pmax_inv_AC,Pmax_Inverter);

% från W till kW
Inverter(k,1)=Plim_Inverter/1000;
Inverter(k,2)=Pmax_Inverter/1000;

% från W/m2 till kW/m2
Solstralning(k)=Solin_hor/1000;
k=k+1;
endfor

% momentant, kW
Inverter_dag(n,:,:)=Inverter;

% ackumulerat över dagen, kWh/m2
Solstralning_dag(n,:)=sum(Solstralning)/2;
n
n=n+1;
endfor

% ackumulerad globalstrålning respektive månad
Solstralning_manad(1)=sum(Solstralning_dag(1:31, :));
Solstralning_manad(2)=sum(Solstralning_dag(32:61, :));
Solstralning_manad(3)=sum(Solstralning_dag(62:90, :));
Solstralning_manad(4)=sum(Solstralning_dag(91:120, :));
Solstralning_manad(5)=sum(Solstralning_dag(121:151, :));
Solstralning_manad(6)=sum(Solstralning_dag(152:181, :));
Solstralning_manad(7)=sum(Solstralning_dag(182:213, :));
Solstralning_manad(8)=sum(Solstralning_dag(214:243, :));
Solstralning_manad(9)=sum(Solstralning_dag(244:274, :));
Solstralning_manad(10)=sum(Solstralning_dag(275:304, :));
Solstralning_manad(11)=sum(Solstralning_dag(305:335, :));
Solstralning_manad(12)=sum(Solstralning_dag(336:365, :));

% solfaktor=1 är molnfritt
Solfaktor=Soldata./Solstralning_manad(1,:);
Solstralning_manad(2,:)=Soldata;

% ackumulerad eleffekt i kWh
% _unlim=invertern inte effektbegränsad
% _lim=invertern är effektbegränsad, se gränser i indatan
% _moln=statistisk molnighet pålagd
Inverter_manad(1,:)=sum(sum(Inverter_dag(1:31,:,:)))/2;
Inverter_manad(2,:)=sum(sum(Inverter_dag(32:61,:,:)))/2;
Inverter_manad(3,:)=sum(sum(Inverter_dag(62:90,:,:)))/2;
Inverter_manad(4,:)=sum(sum(Inverter_dag(91:120,:,:)))/2;
Inverter_manad(5,:)=sum(sum(Inverter_dag(121:151,:,:)))/2;
Inverter_manad(6,:)=sum(sum(Inverter_dag(152:181,:,:)))/2;
Inverter_manad(7,:)=sum(sum(Inverter_dag(182:213,:,:)))/2;
Inverter_manad(8,:)=sum(sum(Inverter_dag(214:243,:,:)))/2;
Inverter_manad(9,:)=sum(sum(Inverter_dag(244:274,:,:)))/2;
Inverter_manad(10,:)=sum(sum(Inverter_dag(275:304,:,:)))/2;
Inverter_manad(11,:)=sum(sum(Inverter_dag(305:335,:,:)))/2;
Inverter_manad(12,:)=sum(sum(Inverter_dag(336:365,:,:)))/2;
Inverter_manad_moln=Inverter_manad.*Solfaktor';
Inverter_ar_unlim=sum(Inverter_manad(:,2))
Inverter_ar_moln_unlim=sum(Inverter_manad_moln(:,2))
Inverter_ar_lim=sum(Inverter_manad(:,1))
Inverter_ar_moln_lim=sum(Inverter_manad_moln(:,1))
Inverter_manad_moln(:,1)

% i kW
Rekord=sum(sum(Inverter_dag(180,:,2)))/2
%Rekord_zzz=sum(sum(Inverter_dag(322,:,2)))/2

% "Plot-Loop" för veckonummer, endast effekt i W, ej ackumulerad
p=1;
q=1:52;
for p=1:52;
  Inverter_vecka(p,:,:)=Inverter_dag(7*p-6,:,:,:);
endfor

Manad=[1 2 3 4 5 6 7 8 9 10 11 12];

figure(10);
plot(l,Inverter_dag(180,:,2));
title('Rekorddagen 2018-06-29');
xlabel('Klockslag');
ylabel('Effekt, kW');
grid on;
xticks([0:24]);
yticks([0:1:18]);
axis([0 24 0 18]);

figure(11);
plot(l,Inverter_dag(322,:,2));
title('Hostdagen 2018-11-18');
xlabel('Klockslag');
ylabel('Effekt, kW');
grid on;
xticks([0:24]);
yticks([0:1:18]);
axis([0 24 0 18]);

figure(20);
plot(Manad,Solfaktor);
title('Globalstralningsfaktor (1=helt molnfritt) for Goteborg 2018');
xlabel('Månad');
ylabel('Solfaktor');
grid on;
xticks([1:12]);
yticks([0.4:0.05:1]);
axis([1 12 0.4 1]);


figure(30);
%[ll,oo]=meshgrid(l,o);
[ll,qq]=meshgrid(l,q);
%mesh(ll,oo,Inverter_dag(:,:,1));
mesh(ll,qq,Inverter_vecka(:,:,1));
%surf(ll,oo,Inverter_dag(:,:,1));
%axis([0 24 1 365 0 17.5]);
axis([0 24 1 52 0 17.5]);
xticks([0:24]);
yticks([0:2:52]);
%yticks([1:365]);
zticks([0:2.5:17.5]);
xlabel('Tid, h');
ylabel('Veckonummer');
%ylabel('Dagnummer');
zlabel('Effekt AC, kW');
title('Solelanlaggning Kungalv vid molnfritt');
