%Calculate the Forward current through the Vactrol LED based on the input
%voltage(0-10.7) and offset(0-10.7)
%Added a "smoothed" model which provides a more realistic response
function If = getCurrent(Vc,offset,smoothed)
A = 3.4645912;
B = 1136.2129956;
Ifmin = 10.1e-6;
Ifmax = 40e-3;
R2max = 10e3;
R6max = 20e3;
R3 = 150e3;
R5 = 100e3;
R4 = 470e3;
R7 = 33e3;
R8 = 4.7e3;
R9 = 470;
VB = 3.9; % Zener Diode Reverse Breakdown Voltage
VF = 0.7; %Zener Diode Forward Voltage
VT = 26e-3; 
n = 3.9696;
kl = 6.3862;
G = 2e5; %Op-amp model open loop gain
Vs = 15; 
k0 = 1.468e2;
k1 = 4.9202e-1;
k2 = 4.1667e-4;
k3 = 7.3915e-9;
gamma = 0.0001; %correction factor
Is =10.1E-5;
beta = 87;
%max voltage in the model is ~ 10.7V. Clip here
if (offset+Vc) > 10.7
    Vc = 10.7-offset;
end
if smoothed == "false"
%Use Parker-D'angelo model
    R6 = R6max;%scale * R6max;
    R1 = (1-offset)*R2max;
    R2 = offset*R2max;

    alpha = 1 + (R6+R7) * (1/R3 + 1/R5);
    beta = ((1/alpha) - 1)/(R6 + R7) - 1/R8;

    bound1 = 600* alpha *n*VT/(G*(R6+R7-1/(alpha*beta))); 

    Ia = Vc/R5 + Vs/(R3*(1+R1/R2));

    if (Ia <= -bound1)

        V3 = -Ia/(alpha*beta);

    elseif ( Ia < bound1)
        x = G*Ia*(R6+R7-1/(alpha*beta))/(alpha*n*VT);
        w = k0 + k1*x + k2*x^2 + k3*x^3;
        V3 = -(alpha/G)*n*VT*w - Ia/(alpha*beta);
     else 
       V3 = kl*alpha/G*n*VT-Ia*(R6+R7);
     end

    Ifbound1 = alpha*(Ifmin- beta*V3);
    Ifbound2 = VB/(R6+R7);
    Ifbound3 = (gamma*G*VB + alpha*R9*(VB*beta+Ifmax))/(gamma*G*(R6+R7) + R9);

    if (Ia <= Ifbound1) 
        If = Ifmin;

    elseif (Ia <= Ifbound2) 
        If = beta * V3 + Ia/alpha;

    elseif (Ia <= Ifbound3) 
        If = gamma * G *(Ia*(R6+R7) - VB)/(alpha*R9) - beta*VB + Ia/alpha;
    else 
        If = Ifmax;

    end
else
    %Use smoothed model
    If = Is*(exp((Vc+offset)*beta/n*VT)-1);
    if If > 0.04
        If = 0.04;
    end
end

