#########################################################################
#####  This script can be used to plot J-V and P-V characteristics  #####
######   Run the script after running a voltage sweep in DEVICE   #######
## The script also calculates Jsc, Voc, FF and photovoltaic efficiency ##
#########################################################################
#############  Copyright 2015 Lumerical Solutions, Inc.    ##############
#########################################################################

clear;
closeall;

Psolar = 100;  # mW/cm^2;

# Get results from voltage sweep

base = getresult("CHARGE","base");

I = pinch(base.I);
V = base.V_base;
Lx = getnamed("CHARGE simulation region","x span")*100;
Ly = getnamed("CHARGE","norm length")*100;

J = I/Lx/Ly*1000;  # mA/cm^2

kern = find(J<0);

kern = length(kern);

slope = (J(kern+1) - J(kern))/(V(kern+1)-V(kern));

# Get Voc and Jsc

Voc = V(kern) + abs(J(kern))/slope;
Jsc = -J(1);

# truncate the J, V arrays to sweep from 0 V to Voc

J = [J(1:kern); 0];
V = [V(1:kern); Voc];

plot(V,J,"Voltage (V)","Current density (mA/cm^2)");
legend("");

# Calculate power voltage relationship

P = abs(J*V);
plot(V,P,"Voltage (V)","Power (mW/cm^2)");
legend("");


?"Voc = " + num2str(Voc) + "V";  #  Open-circuit voltage
?"Jsc = " + num2str(Jsc) + " mA/cm^2";  #  Short-circuit current
?"Pmax = " + num2str(max(P)) + " mW/cm^2";  # Maximum output power
?"FF = " + num2str(max(P)/(Voc*Jsc));  #  Fill-factor
?"eta = " + num2str(max(P)/Psolar*100) + "%";  #  Photovoltaic efficiency