mos.jl

using ModelingToolkit
using StructuralTransformations
using OrdinaryDiffEq
using IfElse: ifelse
using Plots


function NMOS_L1(Vg, Vd, Vs, Vth, W, L, uCox)
    Vgs = Vg-Vs
    Vds = Vd-Vs
    Vo = Vgs-Vth # overdrive
    # ModelingToolkit doesn't handle Boolean expressions, this sometimes works
    ifelse(signbit(Vo),
           zero(Vds), # below threshold is 0
           ifelse(signbit(Vds-Vo),
                  uCox*W/L*(Vo*Vds-Vds^2/2),
                  uCox/2*W/L*Vo^2))
end

mymax(x, y) = ifelse(signbit(x-y), y, x)

# https://www2.eecs.berkeley.edu/Pubs/TechRpts/1990/ERL-90-19.pdf
function NMOS_L6(Vg, Vd, Vs, Vb,
                 B, n, K, m, λ₀, λ₁, Vt0, γ, ϕf, W, Leff)
    Vbs = Vb-Vs
    Vgs = Vg-Vs
    Vds = Vd-Vs
    Vth=Vt0+γ*(√(2ϕf-Vbs)-√(2ϕf))
    Vdsat=K*mymax(1e-9, Vgs-Vth)^m
    Idsat=W/Leff*B*mymax(1e-9, Vgs-Vth)^n
    λ=λ₀-λ₁*Vbs
    Id5=Idsat*(1+λ*Vds) # Vds > Vdsat (saturated region)
    Id3=Id5*(2-(Vds/Vdsat))*Vds/Vdsat
    sat = Vds-Vdsat
    ifelse(signbit(sat), Id3, Id5)
end
# x = LinRange(0, 2, 100)
# res = NMOS_L6.(1, x, 0, 0,
#     4.9721e-5, 1.0484, 0.83496, 0.6193, 0.066265,
#     0.0038573, 0.85502, 0.29648, 0.20556/2, 100e-6, 1e-6)
# plot(x, hcat(res...)')

@parameters t W L C R Vcc f
@variables Vg(t) Vc(t) Id(t) Ic(t) Ir(t)
@derivatives D'~t

# It does not like it AT ALL if you put functions on LHS or derivatives on RHS
eqs = [
    0 ~ -Vg+1+sin(2*pi*f*t),
    0 ~ -Id+NMOS_L6(Vg, Vc, 0, 0,
                    4.9721e-5, 1.0484, 0.83496, 0.6193, 0.066265,
                    0.0038573, 0.85502, 0.29648, 0.20556/2, W, L),
    D(Vc) ~ Ic/C,
    0 ~ -Ir+(Vcc-Vc)/R,
    0 ~ Ir - Id - Ic,
]

sys = ODESystem(eqs, t, [Vg, Vc, Id, Ic , Ir], [W, L, C, R, Vcc, f])
u0 = zeros(5)
tspan = (0.0, 10e-3)
params = [
    W => 100e-6,
    L => 1e-6,
    C => 1e-6,
    R => 10e3,
    Vcc => 5,
    f => 1e3
]
# Turn into a first order differential equation system
# sys = ModelingToolkit.ode_order_lowering(sys)

# Perform index reduction to get an index 1 DAE
# sys = StructuralTransformations.dae_index_lowering(sys)
prob = ODEProblem(sys, u0, tspan, params, jac=true)
# soln = solve(prob,abstol=1e-8, reltol=1e-8)
soln = solve(prob, Trapezoid(), abstol=1e-8, reltol=1e-8)
plot(soln)

using AutoOptimize

_prob, _ = auto_optimize(prob)

_soln = solve(prob, Trapezoid(), abstol=1e-8, reltol=1e-8)

plot(_soln) #differs from soln

# Turn into a first order differential equation system
sys2 = ModelingToolkit.ode_order_lowering(sys)

       # Perform index reduction to get an index 1 DAE
sys3 = StructuralTransformations.dae_index_lowering(sys2) ## fails

prob2 = ODEProblem(sys2, u0, tspan, params, jac=true)

sol2 = solve(prob2, Trapezoid(), abstol=1e-8, reltol=1e-8) ##fails