% Multi-Quadrant Dynamic PV model
% Setup file for SimPowerSystems in Matlab Simulink
% Created by Katherine A. Kim 
% University of Illinois at Urbana-Champaign
% First created 12-5-2012
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
close all
clc

%Load Test Data
load IVdat_Cell8    % Vdat, Idat
load Cdat_Cell8     % C, VdatC, IdatC, Vbd
load LRdat_Cell8    % Rs, Ls

Vbd = -20

%PV panel parameters 
%test parameters
Gin = 0;
Tin = 25;

%Datasheet parameters
Imp = 7.89/12;   %Max. Power Current Ipm (A) 7.89 - 8.16
Vmp = 0.510;   %Max Power Voltage Vpm (V) 0.510 - 0.521
Iscn = 8.49/12;    %Short Circuit Current Isc (A) 8.49 - 8.73
Vocn = 0.6230;  %Open Circuit Voltage Voc (V) 0.623 - 0.625
Ki = 0.00065*Iscn;     %Current Temperature Coefficient [I/K]
Kv = -0.0036*Vocn;   %Voltage Temperature Coefficient [V/K]
Ns = 1;          %Number of series-strung cells

%plot(Vdat,Vdat.*(Idat+Imp))
%[Pmax,imax] = max(Vdat.*(Idat+Imp))
%Vdat(imax)

% Constants
k = 1.3806503e-23;   %Boltzmann [J/K]
q = 1.60217646e-19;  %Electron charge [C]
%a = 1.4;             %Diode constant
Tn = Tin + 273.15;    % Nominal operating temperature [K]
Vtn = k*Tn*Ns/q

figure(1)
plot(Vdat,Idat,'*k')
hold on

%separate the function into three sections: breakdown, middle, and forward
iVbd_srt = find(Vdat<=Vbd,1,'first')
iVbd_lst = find(Vdat<=Vbd,1,'last')
iVmid_srt = find(Vdat>Vbd & Vdat<Vmp,1,'first')
iVmid_lst = find(Vdat>Vbd & Vdat<Vmp,1,'last')
iVf_srt = find(Vdat>=Vmp,1,'first')
iVf_lst = find(Vdat>=Vmp,1,'last')
iRshmeas_srt = find(Vdat>Vbd*0.6 & Vdat<Vmp*0.6,1,'first')
iRshmeas_lst = find(Vdat>Vbd*0.6 & Vdat<Vmp*0.6,1,'last')

Rsh = (Vdat(iRshmeas_lst)-Vdat(iRshmeas_srt))./-(Idat(iRshmeas_lst)-Idat(iRshmeas_srt))
plot([Vdat(iRshmeas_srt) Vdat(iRshmeas_lst)],[Idat(iRshmeas_srt) Idat(iRshmeas_lst)],'rp')
plot(Vdat,-Vdat/Rsh)

Ifit = Idat+Vdat/Rsh;
%plot(Vdat,Ifit,'b*')
Vfit = Vdat + Idat*Rs;
plot(Vfit,Ifit,'gs')

s0 = fitoptions('Method','NonlinearLeastSquares',...
               'Lower',[1.0],...
               'Upper',[2.0],...
               'Startpoint',[1.2]);
 
f0 = fittype(['-' num2str(Iscn) '/(exp(' num2str(Vocn) '/(a*' num2str(Vtn) '))-1)*(exp(x/(a*' num2str(Vtn) '))-1)'],'options',s0); %a=Isr, b=Kr, -Isr*(exp(Krb*(Vbd-Vd)/aVt)-exp(Krb*Vbd/aVt))

srt = min([iVf_srt iRshmeas_srt]);
lst = max([iVf_lst iRshmeas_lst]);
[cfit0,stats] = fit(transpose(Vfit(srt:lst)),transpose(Ifit(srt:lst)),f0)
plot(Vfit(srt:lst),cfit0(Vfit(srt:lst)),'<m')
%axis([0 Vocn -inf 0])
 
a = cfit0.a
aVtn = a*k*Tn*Ns/q;

s0 = fitoptions('Method','NonlinearLeastSquares',...
               'Lower',[0,-Inf],...
               'Upper',[Inf,Inf],...
               'Startpoint',[1,1e-1]);
 
% f0 = fittype('a*(exp(b*(-18.0220-x)/0.0334)-exp(b*(-18.0220)/0.0334))','options',s0); %a=Isr, b=Kr, -Isr*(exp(Krb*(Vbd-Vd)/aVt)-exp(Krb*Vbd/aVt))
f0 = fittype(['a*(exp(b*(' num2str(Vbd) '-x)/' num2str(aVtn) ')-exp(b*(' num2str(Vbd) ')/' num2str(aVtn) '))'],'options',s0); %a=Isr, b=Kr, -Isr*(exp(Krb*(Vbd-Vd)/aVt)-exp(Krb*Vbd/aVt))

srt = min([iVbd_srt iVmid_srt]);
lst = max([iVbd_lst iVmid_lst]);
[cfit0,stats] = fit(transpose(Vfit(srt:lst)),transpose(Ifit(srt:lst)),f0)
plot(Vdat(srt:lst),Ifit(srt:lst),'ob')
plot(Vdat(srt:lst),cfit0(Vdat(srt:lst)),'pr')

Isr = cfit0.a
Kr = cfit0.b
%%
%Capactior Data
VdatC = VdatC(2:end);
C = C(2:end);
IdatC = -IdatC(2:end);

Vmp = 0.412 %hack due to lack of data
iVbd_srt = find(VdatC<=Vbd,1,'first')
iVbd_lst = find(VdatC<=Vbd,1,'last')
iVmid_srt = find(VdatC>Vbd & VdatC<Vmp,1,'first')
iVmid_lst = find(VdatC>Vbd & VdatC<Vmp,1,'last')
iVf_srt = find(VdatC>=Vmp,1,'first')
iVf_lst = find(VdatC>=Vmp,1,'last')

%Vd = Vdat+(-Idat.*Rs);

%Junction Cap
% s1 = fitoptions('Method','NonlinearLeastSquares',...
%                'Lower',[-Inf],...
%                'Upper',[Inf],...
%                'Startpoint',[1e-4]);
% %Determine appropriate phio value (either Vocn or largest meas C voltage)
% phio = ceil(max([Vocn/Ns max(VdatC)])*100)/100+0.2
% 
% % f1 = fittype('a/(1-x/0.638)^(1/2)','options',s1); %a=Cjo, x=Vd
% f1 = fittype(['a/(1-x/' num2str(phio) ')^(1/2)'],'options',s1); %a=Cjo, x=Vd

s1 = fitoptions('Method','NonlinearLeastSquares',...
               'Lower',[-Inf,0],...
               'Upper',[Inf,10],...
               'Startpoint',[1e-5,Vocn]);
f1 = fittype(['a/(1-x/b)^(1/2)'],'options',s1); %a=Cjo, x=Vd, b = phio

[cfit1,stats] = fit(transpose(VdatC(iVmid_srt:iVmid_lst)),transpose(C(iVmid_srt:iVmid_lst)),f1)

figure(2)
srt = min([iVmid_srt iVf_srt]);
lst = max([iVmid_lst iVf_lst]);
semilogy(VdatC,C,'k*')
hold on
semilogy(VdatC(srt:lst),cfit1(VdatC(srt:lst)),'rp--')
Cjo = cfit1.a
phio = cfit1.b

Cdiff = C - transpose(cfit1(VdatC));
plot(VdatC,Cdiff,'m^')

s2 = fitoptions('Method','NonlinearLeastSquares',...
               'Lower',[-Inf],...
               'Upper',[Inf],...
               'Startpoint',[1e-5]);
% f2 = fittype('a*x/0.0334','options',s2); %a=tau, x=Id
f2 = fittype(['a*x/' num2str(aVtn)],'options',s2); %a=tau, x=Id
[cfit2,stats] = fit(transpose(IdatC(iVf_srt:iVf_lst)),transpose(Cdiff(iVf_srt:iVf_lst)),f2)
tau = cfit2.a
plot(VdatC(iVf_srt:iVf_lst),cfit2(IdatC(iVf_srt:iVf_lst)))

% phir = floor(min([min(Vdat) min(VdatC)])*10)/10
% phir = Vbd-aVtn/Kr*log(Iscn*1.5/Isr+1)
%phir = floor((Vbd-aVtn/Kr*log(Iscn*1.5/Isr+1))*10)/10
%phir = -aVtn/Kr*log(exp(-Vbd*Kr/aVtn)*1.5*Iscn/Isr+1)
% phir = floor((-aVtn/Kr*log(exp(-Vbd*Kr/aVtn)*1.5*Iscn/Isr+1))*10)/10
% Xr = 6;

s3 = fitoptions('Method','NonlinearLeastSquares',...
               'Lower',[-Inf],...
               'Upper',[Inf],...
               'Startpoint',[1e-6]);
% f3 = fittype('a/(1-x/-26.5)^(6)','options',s3); %a=Cjr, x=Vd
%f3 = fittype(['a/(1-x/' num2str(phir) ')^(' num2str(Xr) ')'],'options',s3); %a=Cjr, x=Vd
f3 = fittype(['a*(-x)/' num2str(aVtn)],'options',s2); %a=taur, x=Id
% [cfit3,stats] = fit(transpose(VdatC(iVbd_srt:iVbd_lst)),transpose(Cdiff(iVbd_srt:iVbd_lst)),f3)
% Cjr = cfit3.a
% 
% semilogy(VdatC(iVbd_srt:iVbd_lst),Cdiff(iVbd_srt:iVbd_lst),'b*')
% hold on
% plot(VdatC(iVbd_srt:iVbd_lst),cfit3(VdatC(iVbd_srt:iVbd_lst)))
try
    %[cfit3,stats] = fit(transpose(VdatC(iVbd_srt:iVbd_lst)),transpose(Cdiff(iVbd_srt:iVbd_lst)),f3)
    [cfit3,stats] = fit(transpose(IdatC(iVbd_srt:iVbd_lst)),transpose(Cdiff(iVbd_srt:iVbd_lst)),f3)
%     Cjr = cfit3.a
    taubd = cfit3.a
    semilogy(VdatC(iVbd_srt:iVbd_lst),Cdiff(iVbd_srt:iVbd_lst),'b*')
    plot(VdatC(iVbd_srt:iVbd_lst),cfit3(IdatC(iVbd_srt:iVbd_lst)))
catch
%     Cjr = 0
    taubd = 0
end
%% Print out parameters
% clc
display('--------------------')
display('Input Parameters:')
display(['Iscn: ' num2str(Iscn) ' A'])
display(['Vocn: ' num2str(Vocn) ' V'])
display(['Imp: ' num2str(Imp) ' A'])
display(['Vmp: ' num2str(Vmp) ' V'])
display(['Rs: ' num2str(Rs) ' Ohms'])
display(['Rsh: ' num2str(Rsh) ' Ohms'])
display(['Ki: ' num2str(Ki) ' I/K'])
display(['Kv: ' num2str(Kv) ' V/K'])
display(['Ns: ' num2str(Ns) ' '])
display(['Vbd: ' num2str(Vbd) ' V'])
display(' ')
display('Inherent Parameters:')
display(['seires inductance (Ls): ' num2str(Ls) ' H'])
display(['diode ideality factor (a): ' num2str(a) ''])
display(['Isr: ' num2str(Isr) ' A'])
display(['Kr: ' num2str(Kr) ' '])
display(['Cjo: ' num2str(Cjo) ' F'])
display(['phio: ' num2str(phio) ' V'])
display(['tau: ' num2str(tau) ' s'])
display(['taubd: ' num2str(taubd) ' s'])
% display(['Cjr: ' num2str(Cjr) ' F'])
% display(['phir: ' num2str(phir) ' V'])
% display(['Xr: ' num2str(Xr) ''])
display(' ')
display('Copy and Paste Below code into Mask Edit Initialization:')
display('% Fit parameters from measured data')
display(['Ls = ' num2str(Ls) '; %H'])
display(['a = ' num2str(a) ';'])
display('% reverse breakdown parameters')
display(['Isr = ' num2str(Isr) '; %A'])
display(['Kr = ' num2str(Kr) ';'])
display('% capacitance parameters')
display(['Cjo = ' num2str(Cjo) '; %F'])
display(['phio = ' num2str(phio) '; %V'])
display(['tau = ' num2str(tau) '; %s'])
display(['taubd = ' num2str(taubd) '; %s'])
% display(['Cjr = ' num2str(Cjr) '; %F'])
% display(['phir = ' num2str(phir) '; %V'])
% display(['Xr = ' num2str(Xr) ';'])

%% Save all model parameters
 save('simulink_PV_dyn_param.mat','Iscn','Vocn','Rs','Rsh','Ki','Kv','Ns','Vbd','Isr','Ls','Kr','Cjo','tau','taubd');
