# using Pkg
# Pkg.add("Modia")
# Pkg.add("ModiaMath")

using Modia
using Modia.Synchronous
using Modia.Electric
using Modia.Rotational
using Modia.Blocks
using ModiaMath


# Modification of https://www.mathworks.com/help/simulink/slref/pulsegenerator.html
# Uses a mod to get the periodicity of the pulse signal.
@model PulseSignal begin
    amplitude = Parameter(start= 1.0, info = "Amplitude of pulse wave")
    period = Parameter(start = 0.0, info = "Period of pulse wave")
    width = Parameter(start = 0.0, info = "width% of period for pulse wave")

    offset = Parameter(start = 0.0, info = "Offset of output signal y")
    startTime = Parameter(start = 0.0, info = "Output is y=offset for time < startTime")

    y = Variable(start = 0.0, info = "Output signal")
    t = Float(start=0.0) # Time

    @equations begin
        der(t) = 1
        y = offset + (((t-startTime) >= 0 && ((mod(t - startTime, period)/period) <= width/100)) ? amplitude : 0)
    end
end

# https://www.maplesoft.com/documentation_center/online_manuals/modelica/Modelica_Electrical_Analog_Sensors.html#Modelica.Electrical.Analog.Sensors.VoltageSensor
@model VoltageSensor begin
    p = Pin()
    n = Pin()
    v = Voltage()

 @equations begin
     p.i = 0
     n.i = 0
     v = p.v - n.v
 end
end

# https://www.maplesoft.com/documentation_center/online_manuals/modelica/Modelica_Electrical_Analog_Sources.html#Modelica.Electrical.Analog.Sources.SignalCurrent
@model SignalCurrent begin
    p = Pin()
    n = Pin()
    v = Float()
    i = Float()
    @equations begin
        v = p.v - n.v
        0 = p.i + n.i
        i = p.i
    end
end

# https://www.maplesoft.com/documentation_center/online_manuals/modelica/Modelica_Electrical_Analog_Ideal.html#Modelica.Electrical.Analog.Ideal.IdealDiode
@model IdealDiode3 begin # Ideal diode
    v = Voltage(info = "Voltage between `p` and `n`")
    i = Current(info = "Current from `p` to `n`")
    p = Pin(info = "Positive pin")
    n = Pin(info = "Negative pin")
    Ron = Parameter() # 1.0E-5 # Forward state-on differential resistance (closed diode resistance)
    Goff = Parameter() # 1.0E-5 # Backward state-off conductance (opened diode conductance)
    Vknee = Parameter() # Forward threshold voltage
    # off = Variable(start=true) # Switching state
    s = Float(start = 0.0, info = "Auxiliary variable") # For position on the diode characteristic


    #=
    s = 0: knee point
    s < 0: below knee point, diode conducting
    s > 0: above knee point, diode locking
    =#
    @equations begin
       # i = if (v > Vknee); (v - Vknee*(1-Ron*Goff))/Ron else v*Goff end
       i = if Synchronous.positive(s); (1 + Vknee) else (Ron + Vknee) end
       v = if Synchronous.positive(s); (Goff + Goff*Vknee) else (1 + Goff*Vknee) end
       v = p.v - n.v
       0 = p.i + n.i
       i = p.i
   end
end

@model model1 begin
    # Conversions
    mmHgToPa = 133.32
    PaTommHg = 1/133.32
    mLTom3 = 1e-6
    m3TomL = 1e6


    MmPulseCurrentSource = SignalCurrent()
    MmPulse = PulseSignal(amplitude = 0.000171,
            width= 60,
            period = 0.857,
            offset=0,
            startTime=0.341628959) # See matlab
    MmGroundSource = Ground()
    LlCapacitor = Capacitor(C = 2*mLTom3/mmHgToPa)

    AaCurrent = CurrentSensor()


    # AaDiode = IdealDiode2(Vknee=0.1, Ron=0.001, Goff=1e-18)
    AaDiode = IdealDiode()


    AaResistor = Resistor(R=mmHgToPa*m3TomL*.05*3)
    AaInductor = Inductor(L=289210)

    BbResistor = Resistor(R=0.1*mmHgToPa/mLTom3)
    SsResistor = Resistor(R=0.85*mmHgToPa/mLTom3)
    SsCapacitor = Capacitor(C=2*mLTom3/mmHgToPa)
    CcPressure = ConstantVoltage(V=4*mmHgToPa)
    CcGroundSource = Ground()

    OoResistor = Resistor(R=0.05*3*mmHgToPa/mLTom3)
    OoCapacitor = Capacitor(C=0.5*mLTom3/mmHgToPa)

    HhSource = Ground()

    LlVoltageSensorForPressure = VoltageSensor()
    AaVoltageSensorForPressure = VoltageSensor()
    OoVoltageSensorForPressure = VoltageSensor()


  @equations begin

      connect(MmPulse.y, MmPulseCurrentSource.i) # convert physical to current signal
      connect(MmPulseCurrentSource.p, MmGroundSource.p)
      connect(MmPulseCurrentSource.n, AaCurrent.p)
      connect(MmPulseCurrentSource.n, LlCapacitor.p)


      # Current sensor connects directly to resistor is there is no diode
      # Else, current sensor -> diode -> resistor
      # connect(AaCurrent.n, AaResistor.p) # For no diode
      connect(AaCurrent.n, AaDiode.p) # With diode
      connect(AaDiode.n, AaResistor.p) # With diode


      connect(AaResistor.n, AaInductor.p)

      connect(AaInductor.n, BbResistor.p)
      connect(BbResistor.n, SsCapacitor.p)

      connect(BbResistor.n, SsResistor.p)
      connect(SsResistor.n, CcPressure.n)
      connect(CcPressure.p, CcGroundSource.p) # ground

      connect(AaInductor.n, OoResistor.p)
      connect(OoResistor.n, OoCapacitor.p)


      ## Ground in "middle" of circuit
      connect(LlCapacitor.n, HhSource.p) # elastance
      connect(OoCapacitor.n, HhSource.p) # Oo
      connect(SsCapacitor.n, HhSource.p) # ground

      # Voltage Sensors
      connect(OoResistor.p, OoVoltageSensorForPressure.p)
      connect(OoVoltageSensorForPressure.n, OoCapacitor.n)

      connect(AaResistor.p, AaVoltageSensorForPressure.p)
      connect(AaVoltageSensorForPressure.n, AaInductor.n)

      connect(LlCapacitor.p, LlVoltageSensorForPressure.p)
      connect(LlVoltageSensorForPressure.n, LlCapacitor.n)

  end
end;

nSteps = 10000; startTime = 0; stopTime = 150; t = range(startTime, stop=stopTime, length=nSteps)

result1 = simulateModel(model1, t, logTiming = true, logTranslation=true,
                    removeSingularities = true, relTol = 1e-4)

plot(result1, ["OoVoltageSensorForPressure.v", "AaVoltageSensorForPressure.v", "LlVoltageSensorForPressure.v",  "AaCurrent.i"])