Public Sub hu000(k,Ao,Bo,At,Bt,Xt,Nc,Ac,Bc,Sc,Dc,Hc,Shu000)
dim wK(1), ph(1),S0(1,1,1),S1(1,1,1),S2(1,1,1)
call hu006(k,pi/Ac,wK)     
ph(0)=wK(0)*Dc(0): ph(1)=wK(1)*Dc(0)
call hu002(k,Ao,Bo,At,Bt,Xt,Ac,Bc,S2)
call Sinvert(S2,S0)
call hu001(wK,Nc,Ac,Bc,Sc,Dc,Hc,S1)
call hu004(ph, S1, S0,S1)
call hu004(ph, S2, S1,Shu000)
End Sub

Sub hu001(wK,Nc,Ac,Bc,Sc,Dc,Hc, Shu001)
dim SS(1,1,1),ph(1)
for i=Nc to 0 step -1
if i=Nc then
call hu007(wK,Bc,Sc,Hc(i),Hc(i), Shu001)
else
ph(0)=wK(0)*Dc(i+1): ph(1)=wK(1)*Dc(i+1)
call hu007(wK,Bc,Sc,Hc(i),Hc(i),SS)
call hu004(ph, SS, Shu001,Shu001)
end if
next
End Sub

Sub hu002(k,Ao,Bo,At,Bt,Xt,Ac,Bc,Shu002)
dim SS(1,1,1), So(1,1,1),wK(1),ph(1)
call hu006(k, pi/At,wK)     
ph(0)=wK(0)*Xt: ph(1)=wK(1)*Xt
call EH(k, Ao, Bo, At, Bt,SS)
call EH(k, At, Bt, Ac, Bc,So)
Call hu004(ph, SS, So,Shu002)
End Sub



Sub hu003(y0, jy, y1,S)
Dim yy0(1), yy1(1), u(1), s11(1), s22(1), s12(1), s21(1), ymy0(1),ymy1(1), ypy0(1),ypy1(1)
u(0) = 1. : u(1)=0.
call xpy(y1,jy,yy0): call xpy(y0, jy,yy1)
call xmy(y0,yy0,ymy0): call xpy(y0,yy0,ypy0)
call xmy(y1,yy1,ymy1): call xpy(y1,yy1,ypy1)
call xdy(ymy0, ypy0,S11)
call xdy(ymy1, ypy1,S22)
call xpy(u, S11,S12)
call xpy(u, S22,S21)
call Smatrix(s11,s12,s21,s22,S)
End Sub

Sub Sinvert(S,Sinv)
for i=0 to 1: for j=0 to 1
Sinv(1-i,1-j,0)=S(i,j,0): Sinv(1-i,1-j,1)=S(i,j,1)
next: next
End Sub

Sub Smatrix(s11,s12,s21,s22,S)
S(0,0,0)=S11(0):S(0,0,1)=S11(1)
S(0,1,0)=S12(0):S(0,1,1)=S12(1)
S(1,0,0)=S21(0):S(1,0,1)=S21(1)
S(1,1,0)=S22(0):S(1,1,1)=S22(1)
End Sub

Sub Scomp(S,s11,s12,s21,s22)
S11(0)=S(0,0,0): S11(1)=S(0,0,1)
S12(0)=S(0,1,0): S12(1)=S(0,1,1)
S21(0)=S(1,0,0): S21(1)=S(1,0,1)
S22(0)=S(1,1,0): S22(1)=S(1,1,1)
End Sub

Function hu004(phase, S0, S1,So)
Dim ex(1), ex2(1),s011(1),s012(1),s021(1),s022(1) 
Dim zn(1),s111(1),s112(1),s121(1),s122(1), ph(1) 
Dim ss11(1),ss12(1),ss21(1),ss22(1),j(1)  
Dim gg11(1),gg12(1),gg21(1),gg22(1)  
j(0)=0.: j(1)=1.
call Scomp(S0,s011,s012,s021,s022)
call Scomp(S1,s111,s112,s121,s122)
call xoy(j,phase,ph)
call cexp(ph,ex): call x_2(ex,ex2)
call xoy(s022,s111,zn): call xmy(ex2,zn,zn)
call xoy(s012,s111,gg11): call xoy(gg11,s021,ss11):call xdy(ss11,zn,gg11): call xpy(s011,gg11,ss11)
call xoy(s012,s112,ss12): call xoy(ss12,ex,gg12): call xdy(gg12,zn,ss12)
call xoy(s121,s021,ss21): call xoy(ss21,ex,gg21): call xdy(gg21,zn,ss21)
call xoy(s121,s022,gg22): call xoy(gg22,s112,ss22): call xdy(ss22,zn,gg22): call xpy(s122,gg22,ss22)
call Smatrix(ss11,ss12,ss21,ss22,So)
End Function

Sub hu005(k, kc, x, tau)
Dim g2
g2 = k * k - kc * kc
If g2 >= 0 Then
tau(0)=sqr(g2)*x:   tau(1)=0. 
Else
tau(0)=0. :  tau(1) =  -Sqr(-g2) * X
End If
End Sub 

Sub hu006(k, kc, betta)
Dim g2
g2 = k * k - kc * kc
If g2 >= 0 Then
betta(0)=sqr(g2):   betta(1)=0. 
Else
betta(0)=0. :  betta(1) =  -Sqr(-g2) 
End If
End Sub 

        Sub hu007(wK,B,S,H0,H1,So)
        dim s0(1),s1(1),j(1),E(1),Yin(1),R(1),jopa(1), T(1), EpY(1),EmY(1), EpR(1)
         call term_sums(10,wK,B,S,H0,H1,S0,S1)
         jopa2=(1.+S0(0))^2+S0(1)*S0(1)
         jopa(0)=(s1(0)*(1.+s0(0))+s1(1)*s0(1))/jopa2: jopa(1)=(-s1(0)*s0(1)+s1(1)*(1.+s0(0)))/jopa2
         Yin(0)=s0(0)-s1(0)*jopa(0)+s1(1)*jopa(1): Yin(1)=s0(1)-s1(0)*jopa(1)-s1(1)*jopa(0):
        EpY(0)=1.+Yin(0): EpY(1)=Yin(1)
        EmY(0)=1.-Yin(0): EmY(1)=-Yin(1)
        call xdy(EmY,EpY,R)
        EpR(0)=1.+R(0):  EpR(1)=R(1)
        call xoy(jopa,EpR,T)
        call Smatrix(R,T,T,R,So)
        End Sub

        sub term_sums(N,wk,B,S,H0,H1,S0,S1)
        dim j(1),ro_bet(1),tu_bet(1), ss0(1),ss1(1)  
        dim mu 
        j(0)=0.: j(1)=1.
        AM0=1./sqr(B*(H0+B+H1))
        AM1=sqr(2.)*AM0
        ss0(0)=0.: ss0(1)=0.   
        ss1(0)=0.: ss1(1)=0. 
        for m=0 to N step 2
            mu=pi*m/(h0+b+h1)
            if m = 0 then
            Am=AM0*b
            else
            Am=AM1*(sin(mu*(b+h0))-sin(mu*h0))/mu
            end if

        G2=wK(0)*wK(0)-wK(1)*wK(1)-mu*mu
            if G2 > .0 then
            Gm=sqr(G2)
            ro_bet(0)=0.: ro_bet(1)=-0.5/gm/sin(gm*S)
            tu_bet(0)=0.: tu_bet(1)=-1./gm/tan(gm*S)
            else
            Gm=sqr(-G2)
            ro_bet(0)=0.: ro_bet(1)=0.5/gm/sinh(gm*S)
            tu_bet(0)=0.: tu_bet(1)=1./gm/tanh(gm*S)
            end if
        ss0(0)=ss0(0)+Am*Am*tu_bet(0):  ss0(1)=ss0(1)+Am*Am*tu_bet(1):       
        ss1(0)=ss1(0)+Am*Am*ro_bet(0):   ss1(1)=ss1(1)+Am*Am*ro_bet(1):      
        next 
        s0(0)=wK(0)*ss0(0)-wK(1)*ss0(1):  s0(1)=wK(0)*ss0(1)+wK(1)*ss0(0): 
        s1(0)=2.*(wK(0)*ss1(0)-wK(1)*ss1(1)):  s1(1)=2.*(wK(0)*ss1(1)+wK(1)*ss1(0)): 
      end sub

        Sub EH(k, a0, b0, a1, b1,So)
        Dim y(1), jB(1)
        Dim Sum0(1), Sum1(1)
        Dim Y0(1), Sm(1)
        Dim u(1)
        Dim A, b, ai, bi
        u(0)=1. : u(1)=0. 
        If a0 <= a1 And b0 <= b1 Then
        ai = a0
        bi = b0
        A = a1
        b = b1
        call EH0(k,A,B,Ai,Bi,Sum0)
        call EH1(k,A,B,Ai,Bi,Sum1)
        call Yi(k,A,B,Ai,Bi,Y)
        jB(0) = Sum0(0)+Sum1(0)    
        jB(1) = Sum0(1)+Sum1(1)   
  
        Else
        ai = a1
        bi = b1
        A = a0
        b = b0
        call EH0(k,A,B,Ai,Bi,Sum0)
        call EH1(k,A,B,Ai,Bi,Sum1)
        call Yi(k,A,B,Ai,Bi,Y0)
        Sm(0) = Sum0(0)+Sum1(0)    :         Sm(1) = Sum0(1)+Sum1(1)     
        call xdy(u,Y0,Y)
        call xoy(Sm,Y,jB)
        End If

        call hu003(u, jB, Y, So)
        End Sub

        Sub Yi(k, A, b, ai, bi,Y)
        Dim Ko(1), kno(1), Sm(1)  
        Sm(0)=0. : Sm(1)=0.
        xa = pi / A
        xsi = pi / ai
        call hu006(k,xsi,ko) 
        call hu006(k,xa,kno)
        ga2= (ga(1, A, ai)) ^ 2     
        Sm(0)= kno(0)*ga2   
        Sm(1)= kno(1)*ga2   
        Rm = 4. * bi / (ai * A * b)
        ko2=ko(0)*ko(0)+ko(1)*ko(1)
        Y(0) = Rm*(Sm(0)*ko(0)+Sm(1)*ko(1))/ko2
        Y(1)=  Rm*(Sm(1)*ko(0)-Sm(0)*ko(1))/ko2
        End Sub

        Sub EH0(k, A, b, ai, bi,Sum0)
        Dim Ko(1), kno(1), Sm(1)  
        Dim ko2, kno2, ga2, xn,xa,xsi
         If A = ai Then
        Sum0(0) = 0. : Sum0(1)=0. 
        Else
        Sm(0) = 0.: Sm(1)=0. 
        xa = pi / A
        xsi = pi / ai
        call hu006(k,xsi,ko) 
         ko2=ko(0)*ko(0)+ko(1)*ko(1)
       For n = 3 To 21 Step 2
        xn = xa * n
        call hu006(k,xn,kno)  
        kno2=kno(0)*kno(0)+kno(1)*kno(1)
        ga2=(ga(n,A,ai))^2
        Sm(0)=Sm(0)+ga2*kno(0)
        Sm(1)=Sm(1)+ga2*kno(1)
        Next
        Rm = 4. * bi / (ai * A * b)
        Sum0(0)=Rm*(Sm(0)*ko(0)+Sm(1)*ko(1))/ko2
        Sum0(1)=Rm*(Sm(1)*ko(0)-Sm(0)*ko(1))/ko2
        End If

        End Sub

        Sub EH1(k, A, b, ai, bi,Sum1)
        Dim Ko(1), knm(1), Sm(1)   
        Dim knm2
         Sm(0) = 0.: Sm(1)=0.
        xa = pi / A
        xb = pi / b
        xsi = pi / ai
        call hu006(k,xsi,ko)
        Nmax = 11
        Mmax = 5
        If A = ai Then
        Nmax = 1
        Mmax = 10
        End If
        If b = bi Then
        Sum1(0) = 0.: Sum1(1)=0.
        Else
        For n = 1 To Nmax Step 2
        For M = 2 To Mmax Step 2
        xn = xa * n
        xm = xb * M
        xnm=sqr(xn*xn+xm*xm)
        call hu006(k,xnm,knm)

        Anm = (k * k - xn * xn) * (ga(n, A, ai) * gb(M, b, bi)) ^ 2
        knm2=knm(0)*knm(0)+knm(1)*knm(1)
        Sm(0)=Sm(0)+Anm*knm(0)/knm2
        Sm(1)=Sm(1)-Anm*knm(1)/knm2
        Next
        Next
        Rm = 8. / (ai * bi * A * b)
        ko2=ko(0)*ko(0)+ko(1)*ko(1)
        Sum1(0)=Rm*(Sm(0)*ko(0)+Sm(1)*ko(1))/ko2
        Sum1(1)=Rm*(Sm(1)*ko(0)-Sm(0)*ko(1))/ko2
        End If
        End Sub

        Function gb(M, b, bi)
        pi2 = 1.570796326795
        xn = pi2 * M * bi / b
        If Abs(xn) < 0.005 Then
        gb = bi * (1. - xn * xn / 6.)
        Else
        gb = bi * Sin(xn) / xn
        End If
        End Function
        

        Function ga(n, A, ai)
        pi2 = 1.570796326795
        afs = 0.4112335167121
        xn = ai * n / A
        dxn = 1. - xn
        If Abs(dxn) < 0.005 Then
        ga = ai * (1. - afs * dxn * dxn) / (1. + xn)
        Else
        ga = ai * Cos(pi2 * xn) / (1. - xn * xn) / pi2
        End If
        End Function

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
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 y <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 > 1.e-20 then
dB=20.*log(R)/log(10.)
else
dB=-400
end if
End Function