Dim Nx, pi

Public Sub Struc(k, N, jw, w, jC, C, So)
'        common /MANY/Nx
'        complex So(2, 2), S0(2, 2), ph
Dim jW0(100), jW1(100), W0(100), W1(100)
Dim jCW(100), CW(100)
Dim S0(7), Bt(1)
pi = 3.14159265358979
pi2 = 1.5707963267949
Nx = 20
'        E=(1.,0.)
'        j=(0.,1.)

        For i = 0 To N
                For ii = 0 To 2
                jW0(ii) = jw(i, ii)
                W0(ii) = w(i, ii)
                CW(ii) = C(i, ii)
                jCW(ii) = jC(i, ii)
                jW1(ii) = jw(i + 1, ii)
                W1(ii) = w(i + 1, ii)
                Next
'        if(i.eq.0)then
         If i = 0 Then
        Call EL(k, jW0, W0, jCW, CW, jW1, W1, So)
        Else
        Call EL(k, jW0, W0, jCW, CW, jW1, W1, S0)
Bt(0) = k: Bt(1) = 0: x = w(i, 0)
'        ph = k * w(i, 0)
        Call Cas(Bt, x, So, S0, So)
        End If
        Next
'1       format(10f8.3)

End Sub

Sub EL(k, jW0, W0, jC, C, jW1, W1, So)
'        subroutine EL(k, jw0, w0, Jc, C, jw1, w1, So)
'        real A0(0:30),X0(0:30)
'        real A1(0:30),X1(0:30),k
'        integer jW0(0:2),JW1(0:2),JC(0:2)
'        real W0(0:2),W1(0:2),C(0:2)
'        complex So(2, 2)
Dim A0(30), X0(30), A1(30), X1(30)

Call StripJunc(jW0(1), jW0(2), W0(1), W0(2), jC(1), jC(2), C(1), C(2), X0, A0)
Call StripJunc(jW1(1), jW1(2), W1(1), W1(2), jC(1), jC(2), C(1), C(2), X1, A1)
        Call Js(k, C(0), X0, A0, X1, A1, So)
'        Return
'        End

End Sub

Sub Js(k, x, X0(), A0, X1(), A1, So)
'        subroutine Js()
'        common /MANY/Nx
'        real A0(0:30),X0(0:30)
'        real A1(0:30),X1(0:30)
'        Real mu, k
'        complex R0, R, R1, Y0, Y1, E, So(2, 2)
Dim S0(1), S1(1), S2(1), S3(1)
pi = 3.14159265358979
pi2 = 1.5707963267949
        R = 0: R0 = 0: R1 = 0
        For N = 0 To Nx
        g2 = k * k - X0(N) * X1(N)
        If g2 < 0 Then
        w = Sqr(-g2)
        O0 = k / w / tanh(w * x)
           If w * x > 15 Then
           O = 0
           Else
           O = k / w / sinh(w * x)
           End If
        Else
        w = Sqr(g2)
        O0 = -k / w / Tan(w * x)
        O = -k / w / Sin(w * x)
        End If
        R0 = R0 + O0 * A0(N) * A0(N)
        R1 = R1 + O0 * A1(N) * A1(N)
        R = R + O * A0(N) * A1(N)
        Next
        
'        R0 = R0 + cmplx(0, O0 * A0(n) * A0(n))
'        R1 = R1 + cmplx(0, O0 * A1(n) * A1(n))
'        R = R + cmplx(0, O * A0(n) * A1(n))
'        enddo
Call Spath(R0, R, R1, S0, S1)
Call Spath(R1, R, R0, S3, S2)
Call Smatrix(S0, S1, S2, S3, So)
'        Y0 = E / (R0 - R * R / (E + R1))
'        Y1 = E / (R1 - R * R / (E + R0))
'        So(1, 1) = (E - Y0) / (E + Y0)
'        So(2, 2) = (E - Y1) / (E + Y1)
'        So(1, 2) = (E - So(1, 1)) * R / (E + R1)
'        So(2, 1) = (E - So(2, 2)) * R / (E + R0)
'        Return
'        End

End Sub
Sub Spath(R0, R, R1, S0, S1)
Dim y(1)
Call Yin(R0, R, R1, y)
M = (1 + y(0)) ^ 2 + (y(1)) ^ 2
S0(0) = (1 - y(0) * y(0) - y(1) * y(1)) / M
S0(1) = -2 * y(1) / M
MR = 1 + R1 * R1
S1(0) = R * (S0(1) + (1 - S0(0)) * R1) / MR
S1(1) = R * (1 - S0(0) - R1 * S0(1)) / MR

End Sub
Sub Yin(R0, R, R1, y)
Rs = R * R - R0 * R1: Ms = Rs * Rs + R0 * R0
y(0) = R * R / Ms: y(1) = (R1 * Rs - R0) / Ms
End Sub

Sub StripJunc(is0, is1, xs0, xs1, i0, i1, X0, X1, Vx, ax)
'        subroutine StripJunc()
'        common /MANY/Nx
'        real mu,Ax(0:30),Vx(0:30)
        pi = 3.14159265358979
        pi2 = 1.5707963267949
        ep = 1
        If is0 + is1 = 0 Then ep = 2
        Aso = xs1 - xs0
        Ao = Sqr(2 / ep / Aso)
        v = pi2 * (is0 + is1) / (xs1 - xs0)
        a = X1 - X0
        TS1 = Tn(Ao, is0, v, Aso)
        TS0 = Tn(Ao, is0, v, 0)
        RS1 = Rn(Ao, is0, v, Aso)
        RS0 = Rn(Ao, is0, v, 0)
        An0 = Sqr(2 / a)
        AN1 = Sqr(1 / a)
        For N = 0 To Nx
        An = An0
        If N + i0 + i1 = 0 Then An = AN1
        mu = pi2 * (2 * N + i0 + i1) / a
        Vx(N) = mu
        T1 = Tn(An, i0, mu, xs1 - X0)
        T0 = Tn(An, i0, mu, xs0 - X0)
        R1 = Rn(An, i0, mu, xs1 - X0)
        R0 = Rn(An, i0, mu, xs0 - X0)
        If mu = v Then
                If v = 0 Then
                ax(N) = T0 * TS0 * Aso
                Else
        ax(N) = 0.5 * (Aso * (T1 * TS1 + R1 * RS1) + (T1 * RS1 - T0 * RS0) / v)
                End If
        Else
        D1 = v * T1 * RS1 - mu * R1 * TS1
        D0 = v * T0 * RS0 - mu * R0 * TS0
        ax(N) = -(D1 - D0) / (mu * mu - v * v)
        End If
        Next

End Sub

        Function Tn(An, i, mu, x)
        pi2 = 1.5707963267949
        Tn = An * Cos(mu * x - pi2 * i)
        End Function

        Function Rn(An, i, mu, x)
        pi2 = 1.5707963267949
        Rn = An * Sin(mu * x - pi2 * i)
        End Function



Sub Cas(Bt, x, C0, C1, Co)
Dim z(2), E(1), E2(1), C(3, 3), G(3, 2), DZ(1)
pi = 3.14159265358979
EE = Exp(-Bt(1) * x): E(0) = Cos(Bt(0) * x) * EE: E(1) = Sin(Bt(0) * x) * EE
E2(0) = E(0) * E(0) - E(1) * E(1): E2(1) = 2 * E(0) * E(1)
For i0 = 0 To 1: z(i0) = z(i0) + E2(i0): G(0, i0) = C0(i0): G(3, i0) = C1(6 + i0)
For i1 = 0 To 1: z(i0 + i1) = z(i0 + i1) - C0(6 + i0) * C1(i1):
For i2 = 0 To 1: j = i0 + i1 + i2
C(0, j) = C(0, j) + C0(2 + i0) * C1(i1) * C0(4 + i2)
C(1, j) = C(1, j) + C0(2 + i0) * C1(2 + i1) * E(i2)
C(2, j) = C(2, j) + C1(4 + i0) * C0(4 + i1) * E(i2)
C(3, j) = C(3, j) + C1(4 + i0) * C1(2 + i1) * C0(6 + i2)
Next: Next
Next
For j = 0 To 3: For i = 0 To 1: C(j, i) = C(j, i) - C(j, i + 2): Next: Next
DZ(0) = z(0) - z(2): DZ(1) = -z(1): ZM = DZ(0) * DZ(0) + DZ(1) * DZ(1)
For j = 0 To 3
For i0 = 0 To 1: For i1 = 0 To 1:
G(j, i0 + i1) = G(j, i0 + i1) + DZ(i0) * C(j, i1) / ZM
Next: Next
Co(2 * j) = G(j, 0) - G(j, 2): Co(2 * j + 1) = G(j, 1): Next
End Sub


Sub Sinvert(S, Sinv)
For i = 0 To 3: For j = 0 To 1: Sinv(2 * i + j) = S(2 * (3 - i) + j): Next: Next
End Sub

Sub Smatrix(s11, s12, s21, s22, S)
S(0) = s11(0): S(1) = s11(1)
S(2) = s12(0): S(3) = s12(1)
S(4) = s21(0): S(5) = s21(1)
S(6) = s22(0): S(7) = s22(1)
End Sub

Sub Scomp(S, s11, s12, s21, s22)
s11(0) = S(0): s11(1) = S(1)
s12(0) = S(2): s12(1) = S(3)
s21(0) = S(4): s21(1) = S(5)
s22(0) = S(6): s22(1) = S(7)
End Sub


Sub xoy(x, y, xy)
xy(0) = x(0) * y(0) - x(1) * y(1):   xy(1) = x(0) * y(1) + x(1) * y(0)
End Sub

Sub xdy(x, y, xy)
Dim M
M = y(0) * y(0) + y(1) * y(1)
xy(0) = (x(0) * y(0) + x(1) * y(1)) / M:  xy(1) = (-x(0) * y(1) + x(1) * y(0)) / M
End Sub

Sub cexp(x, ex)
Dim a
a = Exp(x(0)): ex(0) = a * Cos(x(1)): ex(1) = a * Sin(x(1))
End Sub

Function cabs(x)
cabs = Sqr(x(0) * x(0) + x(1) * x(1))
End Function

Function arg(x)
Dim up, um
pi = 3.14159265358979
pi2 = 1.5707963267949
up = x(1) + x(0): um = x(1) - x(0)
If up > 0 And um > 0 Then
arg = pi2 - Atn(x(0) / x(1))
End If
If up >= 0 And um <= 0 Then
arg = Atn(x(1) / x(0))
End If
If up < 0 And um < 0 Then
arg = -pi2 - Atn(x(0) / x(1))
End If
If up <= 0 And um >= 0 Then
If x(1) >= 0 Then arg = pi + Atn(x(1) / x(0)) Else arg = -pi + Atn(x(1) / x(0))
End If
End Function

Sub xpy(x, y, xy)
xy(0) = x(0) + y(0):  xy(1) = x(1) + y(1)
End Sub

Sub xmy(x, y, xy)
xy(0) = x(0) - y(0):   xy(1) = x(1) - y(1)
End Sub


Sub conj(x, conjx)
conjx(0) = x(0): conjx(1) = -x(1)
End Sub

Sub clog(x, lnx)
Dim absx, argx
absx = cabs(x)
argx = arg(x)
lnx(0) = Log(absx): lnx(1) = argx
End Sub



Function ch(x)
Dim y
 y = Exp(x): ch = 0.5 * (y + 1 / y)
End Function

Function sh(x)
Dim y
 y = Exp(x): sh = 0.5 * (y - 1 / y)
End Function

Function csqr(x, sqrx)
Dim z, res, ims
z = Sqr(x(0) * x(0) + x(1) * x(1))
res = Sqr((z + x(0)) / 2)
ims = Sqr((z - x(0)) / 2)
sqx(0) = res: sqrx(1) = ims
End Function


Function csin(x, sinx)
sinx(0) = Sin(xre) * ch(xim): sinx(1) = Cos(x(0)) * sh(x(1))
End Function

Function ccos(x, cosx)
cosx(0) = Cos(x(0)) * ch(x(1)): cosx(1) = -Sin(x(0)) * sh(x(1))
End Function

Function x_2(x, x2)
x2(0) = x(0) * x(0) - x(1) * x(1): x2(1) = 2 * x(0) * x(1)
End Function

Function aox(a, x, ax)
aox(0) = a * x(0): aox(1) = a * x(1)
End Function

Function sinh(x)
Dim y
y = Exp(x): sinh = 0.5 * (y - 1 / y)
End Function

Function cosh(x)
Dim y
y = Exp(x): cosh = 0.5 * (y + 1 / y)
End Function

Function tanh(x)
Dim y
If x < 0 Then
y = Exp(2 * x): tanh = (y - 1) / (y + 1)
Else
y = Exp(-2 * x): tanh = (1 - y) / (1 + y)
End If
End Function

Public Sub SdB(So, RdB, TdB)
Dim R(1), T(1)
R(0) = So(0, 0, 0): R(1) = So(0, 0, 1)
T(0) = So(0, 1, 0): T(1) = So(0, 1, 1)
RdB = dB(R): TdB = dB(T)
End Sub


Public Function dB(x)
R = cabs(x)
If R > 1E-20 Then
dB = 20 * Log(R) / Log(10)
Else
dB = -400
End If
End Function


        Function Wave(IndUnit, f)
        pi = 3.14159265358979
        If IndUnit = 0 Then L = 11.80315 / f
        If IndUnit = 1 Then L = 299.792858 / f
        If IndUnit = 2 Then L = 299792.858 / f
        Wave = 2 * pi / L
        End Function
        
        Function dB_Line(f, So)
	spac="       "
        dBL = FormatNumber(f, 5)
        For i = 0 To 3
        R = Sqr((So(2 * i)) ^ 2 + (So(2 * i + 1)) ^ 2)
        If R < 1E-20 Then R = 1E-20
        dBs = 20 * Log(R) / Log(10)
         dBL = dBL & spac&FormatNumber(dBs, 5)
        Next
        dB_Line = dBL
        End Function