%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Laboratory of Renewable Energy Science and Engineering, EPFL  %
%                             2017                               %
%               Yannick Gaudy, Sophia Haussener                  %
%                        Renewable energy                        %
%                           Solution 2                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;
close all;
clc;

% Exercise 2 %

% Constant declaration
h=6.626070040E-34; % J*s Planck constant
c0=299792458; % m/s Speed of light
kB=1.3806488E-23; % J/K
T=6000; %Sun temperature in K
eV=1.602176565E-19; % eV to J

% Efficiency for different bandgap

Egap=0:0.1:5; %eV
Efficiency=zeros(1,length(Egap));
lambda=200:1:2500; % in nm
interval=1e-9; % wavelength interval in m
lambda=lambda/1e9; % in m
E_lambda=h*c0./lambda/eV;
el=zeros(1,length(lambda));
eb=zeros(1,length(lambda));

for k=1:length(Egap)
lambda_limit=h*c0/(Egap(k)*eV);


for i=1:length(lambda)
    eb(i)=2*h*c0^2/(lambda(i)^5*(exp(h*c0/(kB*lambda(i)*T))-1)); % Black body emission power
end
Sum_eb=sum(eb*interval);

for i=1:length(lambda)
    if lambda(i)<lambda_limit
        el(i)=Egap(k)/E_lambda(i)*2*h*c0^2/(lambda(i)^5*(exp(h*c0/(kB*lambda(i)*T))-1));
    end
    if lambda(i)>lambda_limit
        el(i)=0;
    end
end
Sum_el=sum(el*interval);
Efficiency(k)=Sum_el/Sum_eb*100;

end

figure(1);
plot(Egap,Efficiency,'color','k','LineWidth',2.5);hold on; grid on;
set(gca,'Fontsize',14);
xlabel('Bandgap (eV)','Fontsize',14);
ylabel('Solar cell efficiency (%)','Fontsize',14);

fig1=figure(1);
print(fig1,'Solar_cell_efficiency','-dpng');

%%%%%%  Exercise 4   %%%%%%

Tamb=25; % Temperature in Celsius
T1=Tamb+273.15; % Temperature in K
Isc=33; % mA/cm2
Voc=0.55; % V
FF=0.7; % Fill factor
Pmax=FF*Isc*Voc; %mW/cm2
I0=1E-8; % dark current in mA/cm2
IL=Isc+I0; % Light current in mA/cm2
Rs=2; % Initial guess of serie resistance in Ohm*cm2
V=0:0.005:0.6; % Fixed potetential potential
I=zeros(1,length(V));

err_limit=1e-5; % error limit
cu_limit=1E2; % number of loop limit
I(1)=30; %mA/cm2 initial guess
err=1; % error
cu=0; % for loop counting
while(err>err_limit && cu<cu_limit)
I1=IL-I0*(exp(eV*(V(1)+I(1)*Rs/1000)/(kB*T1))-1);
err=abs(I(1)-I1);
I(1)=I1;
cu=cu+1;
end

for i=2:length(V)
    i;
    I(i)=I(i-1); %mA/cm2 initial guess
    err=1; % error
    cu=0; % for loop counting
        while(err>err_limit && cu<cu_limit)
        I1=IL-I0*(exp(eV*(V(i)+I(i)*Rs/1000)/(kB*T1))-1);
        I(i)=I1;
        err=abs(I(i)-I1);
        cu=cu+1;
        end
end

P=I.*V;
Pmax1=max(P);
k=0;

% Run loop

while(Pmax1 > Pmax)
    k=k+1;
    Rs=Rs+0.01;
    err=1; % error
    cu=0; % for loop counting
    while(err>err_limit && cu<cu_limit) % Initial value at V=0
    I1=IL-I0*(exp(eV*(V(1)+I(1)*Rs/1000)/(kB*T1))-1);
    err=abs(I(1)-I1);
    I(1)=I1;
    cu=cu+1;
    end

    for i=2:length(V)
        i;
        I(i)=I(i-1); %mA/cm2 initial guess
        err=1; % error
        cu=0; % for loop counting
            while(err>err_limit && cu<cu_limit)
            I1=IL-I0*(exp(eV*(V(i)+I(i)*Rs/1000)/(kB*T1))-1);
            I(i)=I1;
            err=abs(I(i)-I1);
            cu=cu+1;
            end
    end

P=I.*V;
Pmax1=max(P);
end;

figure(2);
plot(V,I,'color','k','LineWidth',2.5);hold on; grid on;
set(gca,'Fontsize',14);
ylim([0 35]);
xlim([0 V(end)]);
xlabel('Potential (V)','Fontsize',14);
ylabel('Current density (mA/cm^2)','Fontsize',14);


figure(3);
plot(V,P,'color','k','LineWidth',2.5);hold on; grid on;
set(gca,'Fontsize',14);
ylim([0 15]);
xlim([0 V(end)]);
xlabel('Potential (V)','Fontsize',14);
ylabel('Solar cell power (mW/cm^2)','Fontsize',14);

fig2=figure(2);
print(fig2,'IV_Rs','-dpng');

fig3=figure(3);
print(fig3,'P_Rs','-dpng');