function [Id]=mvs_2_0_0_etsoi(coeff,bias_data)
% Symmetrical Short-Channel MOSFET model for ETSOI devices.
% Returns the drain current, Id [A]. This model is only valid for Vg >~
% Vg(psis=phif) where psis is the surface potential.  I.e range of validity
% is from onset of weak inversion through strong inversion. This model uses
% a new inversion charge model; Implemented by Shaloo Rakheja, MIT,

version = 2.00; % Model version
%% input parameters known and not fitted.
global input_parms;

type=input_parms(1);    % Type of transistor. nFET type=1; pFET type=-1
W=input_parms(2);       % Transistor width [m]
Lgdr=input_parms(3);    % Physical gate length [m]. This is the designed gate length for litho printing.
dLg = input_parms(4);   % Overlap length including both source and drain sides [m].
Cox = input_parms(5);   % Insulator capacitance [F/m^2] {calibrated from C-V}
Tjun=input_parms(6);    % Junction temperature [K].

% constants
phit = 8.617e-5*Tjun;   % Thermal voltage [V]
qe = 1.602e-19;         % Electron charge [Col/]
kT = phit*qe;           % Thermal energy [Joules]
hbar = 6.62e-34/(2*pi); % Reduced Planck's constant [J-s]
eps0 = 8.85e-12;        % Permittivity of free space [F/m]
nq = 1/3;               % QM corr. exponent. Theoretical = 1/3

%fitted coefficients
Rs0=coeff(1);                           % Source access resistance [Ohms-meter]
Rd0 = Rs0;                              % Assumed symmetric source and drain. If asymmetric, use coeff(2) for Rd0.
delta=coeff(3);                         % Drain induced barrier lowering (DIBL) [V/V]
n0=coeff(4);                            % Subthreshold swing factor [unit-less] {typically between 1.0 and 2.0}
nd = coeff(5);                          % Punch-through [1/V]
theta = coeff(6);                       % Saturation voltage for critical length
beta = coeff(7);                        % Parameter to control slope of saturation function, Fsat, in trasport
beta_c = coeff(8);                      % Parameter to control slope of critical length saturation function


energy_diff_volt = coeff(9);            % Threshold voltage [V]

mt = coeff(10);                         % Transverse electron mass in Silicon [Kg]
ml = coeff(11);                         % Longitudinal effective mass in Silicon [Kg]

mu_eff=coeff(12);                       % Long-channel effective mobility [m^2/Vs]
ksee = coeff(13);                       % Coefficient for calculating VS injection velocity
B = coeff(14);                          % Stern QM correction numerator
dqm0 = coeff(15);                       % Distance of charge centroid from the interface when channel charge is negligible [m] 
eps = coeff(16);                        % Relative permittivity of semiconductor channel material
nu = coeff(17);                         % Relative occupancy of delta2 subbands
relax0 = coeff(18);                     % Successive iteration relaxation coefficient

Leff = Lgdr-dLg;                        % Effective channel length [m]
QB = (B/dqm0)^(1/nq);                   % QM corr. based on fitted distance of centroid

% Conductivity and DOS electron masses in delta2 and delta4 sub-bands
mC_delta2 = 4*mt;
mD_delta2 = 2*mt;
mC_delta4 = 4*(sqrt(mt)+sqrt(ml)).^2;
mD_delta4 = 4*sqrt(mt*ml);

% Calculate non-degenerate values of thermal velocity and mean free path
vT_delta2_int = sqrt(2*kT/(pi)*mC_delta2/mD_delta2.^2);
vT_delta4_int = sqrt(2*kT/(pi)*mC_delta4/mD_delta4.^2);
vT_avg_int = nu*vT_delta2_int+(1-nu)*vT_delta4_int;
lambda_int = 2*phit*mu_eff./vT_avg_int;

%%%%% ####### model file begins
% Direction of current flow: dir=+1 when "x" terminal is the source dir=-1
% when "y" terminal is the source |

%% bias values
Vd_pre= bias_data(:,1);
Vg_pre= bias_data(:,2);
Vs_pre=bias_data(:,3);

for len_bias=1:length(Vd_pre)
    Vd=Vd_pre(len_bias);
    Vg=Vg_pre(len_bias);
    Vs=Vs_pre(len_bias);
    dir=type*sign(Vd-Vs);
    
    Vds=abs(Vd-Vs);
    Vgs=max(type*(Vg-Vs),type*(Vg-Vd));
    
    %%% Initialization of variables-------------------------------
    psi_solx=0; psi_solxx = 1; Qn=0;  Vdsat = phit;
    Idx = 0;
    Rs = Rs0/W;
    Rd = Rd0/W;
    
    dvg=Idx.*Rs;
    dvd=Idx.*(Rd);
    count=1;
    
    
    %%% Current calculation loop %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    while max(abs((psi_solx-psi_solxx)./abs(psi_solxx+psi_solx)))>1e-12;
        count=count+1;
        if count>1000, break, end
        Idxx=Idx;
        psi_solxx = psi_solx;
        
        
        % The degenerate MFP and thermal velocity are the same as their
        % respective NDG values; since the effects of degeneracy are very
        % low. See Fig.7 in Analysis of Carrier Transport in Short-Channel
        % MOSFETs DOI: 10.1109/TED.2013.2294380
        
        lambda_dg = lambda_int; 
        vT_dg=vT_avg_int; 
        
        dvg=(Idx.*Rs+dvg*0)/1;
        dvd=(Idx.*(Rd)+dvd*0)/1;
        dvds=dvg+dvd;             % total voltage drop from Vds
        Vdsi=Vds-dvds;
        Vgsi = Vgs-dvg;
 
        Es=(energy_diff_volt+psi_solx)/phit;
        expEs=exp(Es);
        FDhalfs=FD_half_integral(Es, expEs);  % After Blakemore
        extr_coef=FDhalfs./log(1+expEs); % Will be used for evaluating charge.
        Ed=(energy_diff_volt+psi_solx-Vdsi)/phit; 
        expEd=exp(Ed);
        
        % calculation of source and drain flux for both delta2 and delta4
        % sub-bands
        flux_source_delta2 = log(1+expEs);
        flux_source_delta4 = flux_source_delta2*(1-nu)./nu;
        flux_drain_delta2 = log(1+expEd);
        flux_drain_delta4 = flux_drain_delta2*(1-nu)./nu;
        
        flux_source = flux_source_delta2+flux_source_delta4; % net source flux
        flux_drain = flux_drain_delta2+flux_drain_delta4; % net drain flux
        
        %%% Calculate Critical length and channel transmission
        Lcrit_lin = Leff;
        Lcrit_sat = ksee*Leff;
        Vdsat_c = theta*phit; % Saturation voltage only for the critical length
        Fsat_c = (abs(Vdsi)./Vdsat_c)./((1+(abs(Vdsi)./Vdsat_c).^beta_c).^(1./beta_c));
        Lcrit = Lcrit_lin.*(1-Fsat_c)+Lcrit_sat.*Fsat_c;
        Tx = lambda_dg./(lambda_dg+Lcrit);
        
        %%% calculate saturation voltage
        coef_sat = (2*flux_source./((2-Tx).*flux_source+Tx.*flux_drain)).^(-1);
        Vdsat = 2*phit*(lambda_dg+Leff)./(lambda_dg+2*ksee*Leff)*coef_sat; % Saturation voltage 
        Fsat=(abs(Vdsi)./Vdsat)./((1+(abs(Vdsi)./Vdsat).^beta).^(1./beta));
        
        %%% Calculate VS injection velocity--------------
        vx0 = vT_dg.*(lambda_dg./(lambda_dg+2*ksee*Leff)); 
        v=vx0.*Fsat;
        
        %%% Calculate charge at VS point --------------------------------
        Qn = -qe*kT/(2*pi*hbar.^2)*mD_delta2.*((2-Tx)*flux_source+Tx*flux_drain);
        n = n0+abs(nd*Vdsi);
        dqm = B/(QB+11/32*abs(Qn)).^nq;
        Cstern = eps*eps0/dqm;
        Cgc=Cstern.*Cox./(Cstern+Cox);    % Effective capacitance [F/m^2]
        
        %%% Calculate currrent
        %%% -------------------------------------------------
        xx=1+((abs(Vdsi)+1e-15)./(abs(Vds)+1e-15))/2;  
        relax=relax0*xx.*(extr_coef).^2;% Increase relax for high Rs,d
        Idx = (W.*abs(Qn).*v + 5*relax.*Idxx)./(5*relax+1);
        psi_solx = ((Vgsi+delta*Vdsi+Qn/Cgc)/n+4*relax.*psi_solxx)./(4*relax+1);
    end
    %%% Current, positive into terminal y  %%%%%%%%%%%%%%%%%%%%%%%%
    Id(len_bias,1)=type*dir.*Idx; % in (Amperes)
end