How to use the andes.consts.Gy function in andes

def gycall(self, dae):
        k = zeros(self.n, 1)
        if self.system.config.forcepq:
            return
        elif self.system.config.forcez:
            if self.v0:
                k = div(2 * dae.y[self.v], self.v0**2)
        else:
            k += mul(self.below, div(2 * dae.y[self.v], self.vmin**2))
            k += mul(self.above, div(2 * dae.y[self.v], self.vmax**2))
        k = mul(self.u, k)

        dae.add_jac(Gy, mul(self.p, k), self.a, self.v)
        dae.add_jac(Gy, mul(self.q, k), self.v, self.v)

def gyisland(self, dae):
        """Reset gy(x) for islanded buses and areas"""
        if self.system.Bus.islanded_buses:
            a = self.system.Bus.islanded_buses
            v = [self.system.Bus.n + item for item in a]
            dae.set_jac(Gy, 1e-6, a, a)
            dae.set_jac(Gy, 1e-6, v, v)

def gycall(self, dae):
        dV2 = mul(self.u, 2 * dae.y[self.v])
        dPdV = mul(dV2, matrix(self.g, (self.n, 1), 'd'))
        dQdV = -mul(dV2, matrix(self.b, (self.n, 1), 'd'))

        dae.add_jac(Gy, dPdV, self.a, self.v)
        dae.add_jac(Gy, dQdV, self.v, self.v)

def gycall(self, dae):
        dae.add_jac(
            Gy,
            mul(0.5, self.ngen, pi, self.rho, (self.R)**2, (self.Vwn)**3,
                div(1, self.mva_mega), (dae.x[self.vw])**3), self.pw, self.cp)
        dae.add_jac(
            Gy,
            mul(25.52, exp(mul(-12.5, dae.y[self.ilamb]))) + mul(
                -12.5, -1.1 + mul(25.52, dae.y[self.ilamb]) + mul(
                    -0.08800000000000001, dae.x[self.theta_p]),
                exp(mul(-12.5, dae.y[self.ilamb]))), self.cp, self.ilamb)
        dae.add_jac(Gy,
                    -(dae.y[self.lamb] + mul(0.08, dae.x[self.theta_p]))**-2,
                    self.ilamb, self.lamb)

def build_gy(self, dae):
        """Build line Jacobian matrix"""
        if not self.n:
            idx = range(dae.m)
            dae.set_jac(Gy, 1e-6, idx, idx)
            return

        Vn = polar(1.0, dae.y[self.a])
        Vc = mul(dae.y[self.v], Vn)
        Ic = self.Y * Vc

        diagVn = spdiag(Vn)
        diagVc = spdiag(Vc)
        diagIc = spdiag(Ic)

        dS = self.Y * diagVn
        dS = diagVc * conj(dS)
        dS += conj(diagIc) * diagVn

        dR = diagIc
        dR -= self.Y * diagVc

def gycall(self, dae):
        dV2 = mul(self.u, 2 * dae.y[self.v])
        dPdV = mul(dV2, matrix(self.g, (self.n, 1), 'd'))
        dQdV = -mul(dV2, matrix(self.b, (self.n, 1), 'd'))

        dae.add_jac(Gy, dPdV, self.a, self.v)
        dae.add_jac(Gy, dQdV, self.v, self.v)

dae.add_jac(Gy, 2 * mul(bsh, V) + mul(gsh, Vsh, sin(theta - thetash)) -
                    mul(bsh, Vsh, cos(theta - thetash)), self.vsh, self.v)
        dae.add_jac(Gy, mul(gsh, V, sin(theta - thetash)) - mul(bsh, V, cos(theta - thetash)), self.vsh, self.vsh)
        dae.add_jac(Gy, mul(gsh, V, Vsh, cos(theta - thetash)) + mul(bsh, V, Vsh, sin(theta - thetash)),
                    self.vsh, self.a)
        dae.add_jac(Gy, -mul(gsh, V, Vsh, cos(theta - thetash)) - mul(bsh, V, Vsh, sin(theta - thetash)),
                    self.vsh, self.ash)

        dae.add_jac(Gy, 0.5 * mul(self.u, 2 * V - 2 * mul(Vsh, cos(theta - thetash)), abs(iIsh), abs(iZsh) ** 2),
                    self.Ish, self.v)
        dae.add_jac(Gy, 0.5 * mul(self.u, 2 * Vsh - 2 * mul(V, cos(theta - thetash)), abs(iIsh), abs(iZsh) ** 2),
                    self.Ish, self.vsh)
        dae.add_jac(Gy, 0.5 * mul(self.u, 2 * V, Vsh, sin(theta - thetash), abs(iIsh), abs(iZsh) ** 2),
                    self.Ish, self.a)
        dae.add_jac(Gy, 0.5 * mul(self.u, 2 * V, Vsh, -sin(theta - thetash), abs(iIsh), abs(iZsh) ** 2),
                    self.Ish, self.ash)

        dae.add_jac(Gy, -2 * mul(self.u, self.k2, dae.y[self.Ish]), self.pdc, self.Ish)

        dae.add_jac(Gy, mul(2 * gsh, Vsh) - mul(gsh, V, cos(theta - thetash)) + mul(bsh, V, sin(theta - thetash)),
                    self.pdc, self.vsh)
        dae.add_jac(Gy, -mul(gsh, Vsh, cos(theta - thetash)) + mul(bsh, Vsh, sin(theta - thetash)),
                    self.pdc, self.v)
        dae.add_jac(Gy, mul(gsh, V, Vsh, sin(theta - thetash)) + mul(bsh, V, Vsh, cos(theta - thetash)),
                    self.pdc, self.a)
        dae.add_jac(Gy, -mul(gsh, V, Vsh, sin(theta - thetash)) - mul(bsh, V, Vsh, cos(theta - thetash)),
                    self.pdc, self.ash)

        for gidx, yidx in zip(self.glim, self.ylim):
            dae.set_jac(Gy, 0.0, [gidx] * dae.m, range(dae.m))
            dae.set_jac(Gy, 1.0, [gidx], [yidx])

def gycall(self, dae):
        Turbine.gycall(self, dae)
        MPPT.gycall(self, dae)
        dae.add_jac(Gy, -dae.x[self.isq], self.isd, self.vsd)
        dae.add_jac(Gy, dae.y[self.vsq], self.isd, self.isd)
        dae.add_jac(Gy, dae.y[self.isd], self.isd, self.vsq)
        dae.add_jac(Gy, -mul(dae.x[self.omega_m], self.xd), self.vsq, self.isd)
        dae.add_jac(Gy, dae.y[self.isd], self.ps, self.vsd)
        dae.add_jac(Gy, dae.y[self.vsd], self.ps, self.isd)
        dae.add_jac(Gy, dae.x[self.isq], self.ps, self.vsq)
        dae.add_jac(Gy, mul(dae.x[self.isq], self.xq - self.xd), self.te,
                    self.isd)
        dae.add_jac(Gy, -div(1, dae.y[self.v1] - dae.y[self.v2]), self.v1,
                    self.ps)
        dae.add_jac(
            Gy, -mul(dae.y[self.ps], (dae.y[self.v1] - dae.y[self.v2])**-2),
            self.v1, self.v2)
        dae.add_jac(Gy,
                    mul(dae.y[self.ps], (dae.y[self.v1] - dae.y[self.v2])**-2),
                    self.v1, self.v1)

def gycall(self, dae):
        if not self.active:
            return
        V = mul(2, self.u, dae.y[self.v])
        dae.add_jac(Gy, mul(self.gf, V), self.a, self.v)
        dae.add_jac(Gy, -mul(self.bf, V), self.v, self.v)

iZsh = div(self.u, abs(Zsh))
        Vh = polar(dae.y[self.v], dae.y[self.a] * 1j)
        Vsh = polar(dae.y[self.vsh], dae.y[self.ash] * 1j)
        Ish = div(Vsh - Vh + 1e-6, Zsh)
        iIsh = div(self.u, Ish)

        gsh = div(self.u, Zsh).real()
        bsh = div(self.u, Zsh).imag()
        V = dae.y[self.v]
        theta = dae.y[self.a]
        Vsh = dae.y[self.vsh]
        thetash = dae.y[self.ash]
        Vdc = dae.y[self.v1] - dae.y[self.v2]
        iVdc2 = div(self.u, Vdc**2)

        dae.add_jac(Gy, div(self.u, Vdc), self.v1, self.pdc)
        dae.add_jac(Gy, -mul(self.u, dae.y[self.pdc], iVdc2), self.v1, self.v1)
        dae.add_jac(Gy, mul(self.u, dae.y[self.pdc], iVdc2), self.v1, self.v2)

        dae.add_jac(Gy, -div(self.u, Vdc), self.v2, self.pdc)
        dae.add_jac(Gy, mul(self.u, dae.y[self.pdc], iVdc2), self.v2, self.v1)
        dae.add_jac(Gy, -mul(self.u, dae.y[self.pdc], iVdc2), self.v2, self.v2)

        dae.add_jac(Gy, -2 * mul(gsh, V) + mul(gsh, Vsh, cos(theta - thetash)) +
                    mul(bsh, Vsh, sin(theta - thetash)), self.ash, self.v)
        dae.add_jac(Gy, mul(gsh, V, cos(theta - thetash)) + mul(bsh, V, sin(theta - thetash)), self.ash, self.vsh)
        dae.add_jac(Gy, -mul(gsh, V, Vsh, sin(theta - thetash)) + mul(bsh, V, Vsh, cos(theta - thetash)),
                    self.ash, self.a)
        dae.add_jac(Gy, mul(gsh, V, Vsh, sin(theta - thetash)) - mul(bsh, V, Vsh, cos(theta - thetash)),
                    self.ash, self.ash)

        dae.add_jac(Gy, 2 * mul(bsh, V) + mul(gsh, Vsh, sin(theta - thetash)) -