clear all
close all


global defpara drtdef actpara drtact drtrec 

% Jsc=21;                     %mA/cm^2
qq=1.602e-19;               %C
%% define all parameters relevant for your cell
for i=1  % variable to loop a parameter
    

defpara.WFback=4.2 ;        %eV
defpara.WFfront=5.5;        %eV
defpara.Sn_back=1.000e+5;          %m/s
defpara.Sp_back=1.000e+8;          %m/s
defpara.Sn_front=1.000e+8;          %m/s
defpara.Sp_front=1.000e-1;          %m/s


defpara.PVK_d=500.*10^(-9) ; %m
defpara.PVK_vth_n=1.000e+05 ; %m/s
defpara.PVK_vth_p=1.000e+05 ; %m/s
defpara.PVK_eps=68 ; %F/m
defpara.PVK_chi=4.15 ; %eV
defpara.PVK_Eg=2.4 ; %eV
defpara.PVK_Nc=2.000000e+24 ; %/m^3
defpara.PVK_Nv=2.000000e+24 ; %/m^3
defpara.PVK_mu_n=2.000000e-06; %m^2/Vs
defpara.PVK_mu_p=2.000000e-06; %m^2/Vs
defpara.PVK_kdir=6.000000e-17 ;%m^3/s
defpara.PVK_Na=0; %m^-3
defpara.PVK_Nd=5e23; %m^-3




%% define working point and uncomment measurements you want to execute
actpara.T=300;              %K
actpara.WPvolt=0;           %V
actpara.WPfreq=1e3;         %Hz
actpara.dark=1;             %0 -> illuminated measurements, 1 -> no illumination
% actpara.Grate=Jsc/qq/1000;  %1/s/cm^2

% IV
% actpara.IV='Calculate!';
% actpara.IVstart=0;          %V
% actpara.IVstop=1.4;         %V
% actpara.IVinc=0.02;         %V

% %CV
% actpara.CV='Calculate!';
% actpara.CVstart=-1;          %V
% actpara.CVstop=1.5;         %V
% actpara.CVinc=0.02;         %V

% % %Cf
actpara.Cf='Calculate!';
actpara.Cfstart=1e-1;       %Hz
actpara.Cfstop=1e2;         %Hz
actpara.Cfppd=10;           %points per decade
% % 
% %QE
% actpara.QE='Calculate!';
% actpara.QEstart=300;        %nm
% actpara.QEstop=900;         %nm
% actpara.QEinc=10;           %nm



drtdef='C:\Users\s.ravi.shankar\sciebo2\Sandheep\Mott-Schottky\code\deffile.def';            %put your definition filename here
drtact='C:\Users\s.ravi.shankar\sciebo2\Sandheep\Mott-Schottky\code\actfile.act';            %put your action filename here
drtrec='C:\Users\s.ravi.shankar\sciebo2\Sandheep\Mott-Schottky\code\recfile.srf';            %put your recorder filename here
scriptfile='C:\Users\s.ravi.shankar\sciebo2\Sandheep\Mott-Schottky\code\scriptfile.script'; %put your Scriptfilename here

%% write the your parameters and measurements in the files that SCAPS reads

writedeffile_doped_no_contacts

writeactfile

writerecfile


%% write your script file

fileID=fopen(scriptfile,'w');                                %open your script file
fprintf(fileID,'set quitscript.quitscaps\n ')                 %closes SCAPS when the simulations are done
fprintf(fileID,'load definitionfile %s\n',drtdef)               %makes SCAPS load your other files
fprintf(fileID,'load actionlistfile %s\n', drtact)
fprintf(fileID,'load recordersettingsfile  %s\n', drtrec)
% fprintf(fileID,'load batchsettingsfile %s\n', batfile)
fprintf(fileID,'calculate singleshot\n')                      %calculates all measurements once
fprintf(fileID,'calculate recorder\n')                          %excecutes measurement and records the quantities specified in recorder file
% fprintf(fileID,'save results.iv results%d.iv\n',i)                 %saves IV results to results.iv
% fprintf(fileID,'save results.cf results.cf\n,i')                 %saves Cf results to results.cf
% fprintf(fileID,'save results.cv results%d.cv\n',i)                 %saves CV results to results.cv
%fprintf(fileID,'save results.qe results.qe\n')                 %saves QE results to results.qe
fprintf(fileID,'save results.recorder results%d.eb\n',i)           %saves recorder results, file type of recorder results depends on recorder input
fclose(fileID)                                                  %closes your scriptfile


orderToSCAPS=['! C:\Scaps3307\scaps3307.exe ' scriptfile];      %creates the commandline: 'SCAPS, run scriptfile'
eval(orderToSCAPS)                                              %executes the commandline for SCAPS
%% j-V curve
% filename1=sprintf('C:\\Scaps3307\\results\\results%d.iv',i);
% [header1, data1] = mhdrload(filename1);
% voltage=data1(:,1,1);
% j=data1(:,2,1);
% figure(i+1);
% plot(voltage,-j,'LineWidth',2);
% hold on
% xlabel('V (V)') 
% ylabel('j (mA.cm^{-2})')
% set(gca,'YLim',[-5 20])
% set(gca,'XMinorTick','on','YMinorTick','on')
% set(gca,'fontsize',12)
% xL = get(gca, 'XLim');
% plot(xL, [0 0],'HandleVisibility','off') 
% hold off
%% data for band diagram and carrier concentrations
filename2=sprintf('C:\\Scaps3307\\results\\results%d.eb',i);
[header2, data2] = mhdrload(filename2);
xpos=data2(:,2,1).*1000;   % position in nm
e_conc=data2(:,2,2) ;        % electron concentration in cm^-3
h_conc=data2(:,2,3);         % hole concentration in cm^-3
E_cb=data2(:,2,7);          % conduction band energy in eV
E_vb=data2(:,2,8);          % valence band energy in eV
E_fn=data2(:,2,9);          % electron Fermi energy in eV
E_fp=data2(:,2,10);          % electron Fermi energy in eV
%% concentrations figure
figure(i+2);
c = {'b' ; 'r' ; 'k'} ;
hold on
semilogy(xpos, e_conc,'LineWidth',2,'color',c{i+1});
semilogy(xpos, h_conc,'LineWidth',2,'color',c{i+2});
legend({'n','p'},'fontsize',12)
set(gca, 'YScale', 'log');
% set(gca,'XLim',[0 392])
% set(gca,'YLim',[10^(0) 10^(20)])
xlabel('x(nm)')
ylabel('Concentration (cm^{-3})')
hold off
%% Band diagram
figure(i+3);
plot(xpos, E_cb,'LineWidth',2);
hold on
plot(xpos, E_vb,'LineWidth',2) ;
plot(xpos, E_fn,'--','LineWidth',2) ;
plot(xpos, E_fp, '--','LineWidth',2) ;
legend({'E_{cb}','E_{vb}','E_{fn}','E_{fp}'},'fontsize',12)
% set(gca,'XLim',[0 392])
xlabel('x(nm)')
ylabel('Energy (eV)')   
hold off
%% Mott Schottky
% filename3=sprintf('C:\\Scaps3307\\results\\results%d.cv', i);
% [header3, data3] = mhdrload(filename3);
% V=data3(:,1);
% Cap=data3(:,2).*10.^(-9);                               %F/cm^2
% export_cap(:,i)=Cap;
% pp=13;                                                  % scaling power, change line 186 when changing g
% MS=1./(10.^(pp).*Cap.^2);                                % 1/C^2 (F^-2.cm^4)
% MS_export(:,i)=1./Cap.^(2);  
% % 
% % % % figure(i+4);
% figure(1);
% semilogy(V, Cap);
% hold on
% % set(gca,'XLim',[-0.2 1.4])                                    % x-axis range
% % set(gca,'YLim',[1e-8 1e-7])                                    % y-axis range
% % 
% set(gca,'XMinorTick','on','YMinorTick','on')
% xlabel('V (V)') 
% % % ylabel('C^{-2} x 10^{15} (F^{-2}.cm^4)')                  % 10 raised to the power pp
% ylabel('Capacitance (F.cm^{-2})')
% set(gca,'fontsize',13)
% % outlier = excludedata(V,MS,'domain',[0.9 1.06]);            % specify fitting domain
% % V_fit=V(~outlier);                                           % x data for fiting
% % MS_fit=MS(~outlier);                                         % y data for fitting
% % p=polyfit(V_fit, MS_fit, 1)     ;                             %fit
% % scaled_slope=abs(p(1));
% % fit=polyval(p, V);                                           %plot fit
% % Vfb=abs(p(2)./p(1))                                         %flatband potential
% % plot(V, fit);     
% figure(2);
% plot(V, MS);
% hold on
% set(gca,'XMinorTick','on','YMinorTick','on')
% xlabel('V (V)') 
% ylabel('C^{-2} x 10^{13} (F^{-2}.cm^4)')                  % 10 raised to the power pp
% set(gca,'fontsize',13)

%% Doping profile
% filename4=sprintf('C:\\Scaps3307\\results\\results%d.cv', i);
% [header4, data4] = mhdrload(filename4);
% V=data4(:,1);
% dep_width=data4(:,4); %micrometer
% dop_den=data4(:,5) ;  % cm^-3
% M(:,i)=dop_den ;
% N(:,i)=dep_width;
% semilogy(dep_width, dop_den,'LineWidth',2);
% hold on
% set(gca,'XLim',[0 40])                                    % x-axis range
% set(gca,'YLim',[10^(11) 10^(18)])                                    % y-axis range
% set(gca,'fontsize',13)
% set(gca,'XMinorTick','on','YMinorTick','on')
% set(gca,'YTick',[10^(11) 10^(12) 10^(13) 10^(14) 10^(15) 10^(16) 10^(17) 10^(18)]);
% xlabel('apparent depletion width (µm)') 
% ylabel('N_d (cm^{-3})') 
%% Dopant density calculation
% q=1.602.*10.^(-19); % C
% eps_r=75;
% eps_vac=8.854.*10.^(-14); %% F.cm^-1
% slope=scaled_slope.*(10.^(g)); %% slope of mott schottky plot rescaled 
% Nd=2./(q.*slope.*eps_r.*eps_vac) % cm^-3

end
%  xL = get(gca, 'XLim');                                     
% legend({'1 kHz','5 kHz','10 kHz','50 kHz','100 kHz','500 kHz'},'fontsize',13)
% set(gca,'fontsize',12) 
% plot(xL, [0 0],'HandleVisibility','off')                 % reference line at y=0
% export_data=[M N];
export_data=[xpos E_cb E_vb E_fn E_fp e_conc h_conc];
% save result.dat export_data -ascii
