C  ---------------------------------------------------------------------
C  REFLECTION AND TRANSMISSION COEFFICIENTS FOR A STACK OF PLANAR LAYERS
C                  WITH OPTIONAL INTERFACE IMPEDANCES
C  ---------------------------------------------------------------------
C
C                      v 1.7   16 December 2000
C
C           (layers may now have permeability and permittivity)
C
C           copyright 1997,8,9,2000     all rights reserved
C
C           B. Chambers
C           Communications and Radar Group
C           Dept of Electronic & Electrical Eng
C           University of Sheffield, S1 3JD, UK
C
C           e-mail:  b.chambers@sheffield.ac.uk
C           WWW:     http://www.shef.ac.uk/~el1bc
C
C           NOTE: Some input and graphical commands are specific
C                 to Salford Fortran FTN77
C
C           Ref:  "Radar Cross-Section", Knott, Shaeffer and Tuley
C                 Artech House, 2nd edition, 1993, pp 303 - 308
C
C           NOTE: non-standard time convention for jw
C     ----------------------------------------------------------------
C
C     n      = number of dielectric layers in stack
C     np2    = number of layers in calculation (n+2) (1 and np2 are air)
C     np1    = number of interfaces between layers
C     ityp() = type of interface, 1..parallel RC, 2..series RC
C     nmax   = max allowed number of layers in stack (already set to 10)
C     yn()   = admittance of a layer
C     yj()   = admittance of an interface
C     xj()   = x coordinate of an interface (xj(1)=0 for back of stack)
C     d()    = thickness of a layer in mm
C     th     = angle of incidence in degs
C     p      = polarisation, 1...parallel, 2...perpendicular
C
C     Input data stored in 'fname.dat' file
C     Computed transmission and reflection data stored in 'fname.plt'
C
C     -----------------------------------------------------------------

      parameter (nmax=10,np1max=nmax+1,np2max=nmax+2)
      real refa(20000),tlab(20000),gj(np1max),cj(np1max),xj(np1max),
     $ d(np2max),er(np2max),tde(np2max),lj(np1max),
     $ mu(np2max),tdm(np2max)

      complex gamma(np2max),a(np2max),b(np2max),yn(np2max),yj(np1max),
     $ erc(np2max),rhoa,tab,ce1,ce2,ce3,ce4,ci,cst(np2max),cct(np2max),
     $ cmj,muc(np2max)

      integer*2 iq,fora,morc
      integer p,ityp(np1max)
      character*8 fname
      character*12 fnplt,fndat

      pi=4.*atan(1.)
      rad=pi/180.
      iendflag=0
      scale=-40.  ! initial y axis = 0 to -40 dB on plots
      cmj=cmplx(0.,1.)

      call gk2ecl   !Salford FTN77 command - clears screen
      write(6,*)
      write(6,*)'Analysis of a dielectric stack with optional interface
     $impedances'
      write(6,*)'_______________________________________________________
     $__________'
      write(6,*)
      write(6,*)
      call coua@(' Enter system characteristic impedance in ohms ')
      read(5,*)z00
      write(6,*)


      call coua@(' Enter output data file name ')
      read(5,300)fname
300   format(a8)


      fndat=fname//'.dat'
      fnplt=fname//'.plt'
      open(10,file=fndat,status='unknown')
      open(11,file=fnplt,status='unknown')

      write(10,597)
597   format('Dielectric stack with optional interface impedances')
      write(10,598)
598   format('---------------------------------------------------')
      write(10,*)
      write(10,700)z00
700   format('System characteristic impedance = ',f6.1,' ohms')

      write(6,*)
      call coua@(' Enter freq start, stop and step in GHz  ')
      read(5,*)freq1,freq2,freq3
      write(10,600)freq1,freq3,freq2
600   format('Start, step and stop freq = ',3f8.3,' GHz')
      ifreq1=nint(1.e3*freq1)
      ifreq2=nint(1.e3*freq2)
      ifreq3=nint(1.e3*freq3)

      write(6,*)
      call coua@(' Enter angle of incidence in degs ')
      read(5,*)th
      write(10,550)th
550   format('Angle of incidence = ',f5.1,' degs')

2005  write(6,*)
      call coua@(' Specify polarisation, 1...parallel, 2...perpendicular
     $  ')
      read(5,*)p
      if(p.lt.1.or.p.gt.2)goto2005
      if(p.eq.1)write(10,551)
      if(p.eq.2)write(10,552)
551   format('Polarisation is parallel')
552   format('Polarisation is perpendicular')

2001  write(6,*)
      call coua@(' Enter number of dielectric layers in stack  ')
      read(5,*)n
      if(n.gt.nmax)goto 2001

      write(10,601)n
601   format('No of dielectric layers in stack = ',i3)

      np1=n+1
      np2=n+2
      xj(1)=0.
      y0=1./z00
      erc(1)=cmplx(1.,0.)
      muc(1)=cmplx(1.,0.)
      erc(np2)=erc(1)
      muc(np2)=muc(1)
      ci=cmplx(0.,1.)

      do 10 i=1,np1
      gj(i)=0.
      cj(i)=0.
      yj(i)=cmplx(0.,0.)
10    continue

      write(10,599)
599   format('Layer data (from back) d(mm), Er, tanDe, Ur, tanDm:')
      write(10,*)
      write(6,*)
      write(6,*)'Enter data for each layer in the stack starting from th
     $e back (layer no 1)'
      write(6,*)
      do 1 i=2,np1
      im=i-1
      write(6,1000)im
1000  format(' For layer number ',i2,' enter d(mm),Er,tanDe,Ur,tanDm')
      read(5,*)d(im),er(i),tde(i),mu(i),tdm(i)
      write(10,602)d(im),er(i),tde(i)
602   format(3f8.4)
      dmm=1.e-3*d(im)
      xj(i)=xj(im)+dmm
      erc(i)=er(i)*cmplx(1.,tde(i))
      muc(i)=mu(i)*cmplx(1.,tdm(i))
      yn(i)=y0*csqrt(erc(i)/muc(i))

1     continue

      cst(np2)=cmplx(sin(rad*th),0.)
      cct(np2)=csqrt(1.-cst(np2)*cst(np2))

C     calculate incidence angles in each layer using Snell's Law

      do 333 i=np1,1,-1
      cst(i)=cst(i+1)*csqrt((erc(i+1)*muc(i+1))/(erc(i)*muc(i)))
      cct(i)=csqrt(1.-cst(i)*cst(i))
333   continue

      write(6,*)
      write(6,*)'Is there an interface impedance between the air and the
     $ back layer ?  '
      call yesno(iq)
      if(iq.eq.1)then

2002  write(6,*)
      write(6,*)' Choose type'
      write(6,*)' 1....parallel RC'
      write(6,*)' 2....series RC'
      write(6,*)' 3....parallel RL'
      write(6,*)' 4....series RL'
      read(5,*)ityp(1)
      if(ityp(1).lt.1.or.ityp(1).gt.4)goto 2002

      write(6,*)
      call coua@(' Enter interface sheet resistance in ohms/sq ')
      read(5,*)sr
      if(sr.lt.0.1)then
      sr=0.001
      iendflag=1
      endif
      gj(1)=1./sr
      write(10,*)
      write(10,603)sr
603   format('Air-last layer interface,  R = ',f10.1,' ohms/sq')

      if(ityp(1).lt.3)then
      write(6,*)
      call coua@(' Enter interface capacitance in fF/sq ')
      read(5,*)sc
      cj(1)=1.e-15*sc
      write(10,613)sc
613   format('                           C = ',f10.1,' fF/sq')
      endif

      if(ityp(1).gt.2)then
      write(6,*)
      call coua@(' Enter interface inductance in nH/sq ')
      read(5,*)sc
      lj(1)=1.e-9*sc
      write(10,6613)sc
6613  format('                           L = ',f10.1,' nH/sq')
      endif

      endif

      if(n.gt.1)then
      do 20 i=2,n
      im=i-1
      write(6,*)
      write(6,2010)im,i
2010  format(' Is there an interface impedance between layers ',i2,
     $' and ',i2,' ?')
      call yesno(iq)
      if(iq.eq.1)then

2003  write(6,*)
      write(6,*)' Choose type'
      write(6,*)' 1....parallel RC'
      write(6,*)' 2....series RC'
      write(6,*)' 3....parallel RL'
      write(6,*)' 4....series RL'
      read(5,*)ityp(i)
      if(ityp(i).lt.1.or.ityp(i).gt.4)goto 2003

      write(6,*)
      call coua@(' Enter interface sheet resistance in ohms/sq ')
      read(5,*)sr
      write(10,*)
      write(10,604)im,i,sr
604   format('Interface between layers',i3,' and',i3,', R = ',f10.1,' oh
     $ms/sq')
      gj(i)=1./sr

      if(ityp(i).lt.3)then
      write(6,*)
      call coua@(' Enter interface capacitance in fF/sq ')
      read(5,*)sc
      cj(i)=1.e-15*sc
      write(10,6130)sc
6130  format('                                    C = ',f10.1,' fF/sq')
      endif

      if(ityp(i).gt.2)then
      write(6,*)
      call coua@(' Enter interface inductance in nH/sq ')
      read(5,*)sc
      lj(i)=1.e-9*sc
      write(10,66130)sc
66130 format('                                    L = ',f10.1,' nH/sq')
      endif

      endif
20    continue

      endif
      write(6,*)
      write(6,*)'Is there an interface impedance between the front layer
     $ and air ? '
      call yesno(iq)
      if(iq.eq.1)then

2004  write(6,*)
      write(6,*)' Choose type'
      write(6,*)' 1....parallel RC'
      write(6,*)' 2....series RC'
      write(6,*)' 3....parallel RL'
      write(6,*)' 4....series RL'
      read(5,*)ityp(np1)
      if(ityp(np1).lt.1.or.ityp(np1).gt.4)goto 2004

      write(6,*)
      call coua@(' Enter interface sheet resistance in ohms/sq ')
      read(5,*)sr
      write(10,*)
      write(10,605)sr
605   format('First layer-air interface, R = ',f10.1,' ohms/sq')
      gj(np1)=1./sr

      if(ityp(np1).lt.3)then
      write(6,*)
      call coua@(' Enter interface capacitance in fF/sq ')
      read(5,*)sc
      cj(np1)=1.e-15*sc
      write(10,6133)sc
6133  format('                           C = ',f10.1,' fF/sq')

      endif

      if(ityp(np1).gt.2)then
      write(6,*)
      call coua@(' Enter interface inductance in nH/sq ')
      read(5,*)sc
      lj(np1)=1.e-9*sc
      write(10,66134)sc
66134 format('                           L = ',f10.1,' nH/sq')
      endif

      endif

      if(iendflag.eq.0)then
      a(1)=cmplx(1.,0.)
      b(1)=cmplx(0.,0.)
      else
      a(2)=cmplx(1.,0.)
      b(2)=cmplx(-1.,0.)
      endif

      yn(1)=cmplx(y0,0.)
      yn(np2)=yn(1)

C     step frequency

950   lll=0
      do 2 k=ifreq1,ifreq2,ifreq3
      f=1.e6*float(k)
      fghz=1.e-3*float(k)
      lll=lll+1			!count number of points to be plotted
      do 3 i=1,np2
      gamma(i)=ci*2.e-8*pi*f*csqrt(erc(i)*muc(i))/3.
3     continue

      do 33 i=1,np1
      if(ityp(i).eq.1) yj(i)=cmplx(gj(i),-2.*pi*f*cj(i))
      if(ityp(i).eq.2) yj(i)=2.*pi*f*cj(i)*gj(i)/cmplx(2.*pi*f*cj(i),
     $+gj(i))
      if(ityp(i).eq.3) yj(i)=cmplx(gj(i),1./(2.*pi*f*lj(i)))
      if(ityp(i).eq.4) yj(i)=gj(i)/(1.-cmj*2.*pi*f*lj(i)*gj(i))
33    continue

      if(iendflag.eq.0)then
      jss=1
      else
      jss=2
      endif

      do 4 j=jss,np1
      ce1=cexp(gamma(j+1)*xj(j)*cct(j+1))
      ce2=cexp(-gamma(j)*xj(j)*cct(j))
      ce3=1./ce2
      ce4=1./ce1

C     parallel polarisation case

      if(p.eq.1)then
      a(j+1)=0.5*ce1*(a(j)*(yn(j+1)*cct(j)+yn(j)*cct(j+1)
     $       +yj(j)*cct(j)*cct(j+1))*ce2
     $       +b(j)*(yn(j+1)*cct(j)-yn(j)*cct(j+1)
     $       +yj(j)*cct(j)*cct(j+1))*ce3)/(yn(j+1)*cct(j+1))

      b(j+1)=0.5*ce4*(a(j)*(yn(j+1)*cct(j)-yn(j)*cct(j+1)
     $       -yj(j)*cct(j)*cct(j+1))*ce2
     $       +b(j)*(yn(j+1)*cct(j)+yn(j)*cct(j+1)
     $       -yj(j)*cct(j)*cct(j+1))*ce3)/(yn(j+1)*cct(j+1))
      endif

C     perpendicular polarisation case

      if(p.eq.2)then
      a(j+1)=0.5*ce1*(a(j)*(yn(j+1)*cct(j+1)+yn(j)*cct(j)+yj(j))*ce2
     $       +b(j)*(yn(j+1)*cct(j+1)-yn(j)*cct(j)+yj(j))*ce3)
     $       /(yn(j+1)*cct(j+1))

      b(j+1)=0.5*ce4*(a(j)*(yn(j+1)*cct(j+1)-yn(j)*cct(j)-yj(j))*ce2
     $       +b(j)*(yn(j+1)*cct(j+1)+yn(j)*cct(j)-yj(j))*ce3)
     $       /(yn(j+1)*cct(j+1))
      endif

4     continue

C     correct phase angles for stack thickness in air

      angcorr=360.*xj(np1)*float(k)/300.

      rhoa=b(np2)/a(np2)
      refa(lll)=20.*alog10(cabs(rhoa))
      refang=-atan2(aimag(rhoa),real(rhoa))/rad-2.*angcorr
2000  if(refang.lt.-360.)then
      refang=refang+360.
      goto 2000
      end if
      if(iendflag.eq.0)then
      tab=1./a(np2)
      tabang=-atan2(aimag(tab),real(tab))/rad-angcorr
3000  if(tabang.lt.-360.)then
      tabang=tabang+360.
      goto 3000
      end if
      tlab(lll)=20.*alog10(cabs(tab))
      else
      tlab(lll)=-1000.
      endif
      write(11,100)fghz,refa(lll),refang,tlab(lll),tabang
2     continue
100   format(f8.4,4f11.4)

      morc=2      ! 1 for monochrome, 2 for colour plot
      fora=1      ! 1 for x axis = frequency

      call plot(lll,refa,tlab,freq1,freq2,fname,fora,morc,scale)
      write(6,*)
555   call coua@('Plot again with different scale ?  (Y/N) ')
      call yesno(iq)
      if (iq.eq.1)then
      write(6,*)
      call coua@('Enter new scale in -dB ')
      read(5,*)scale
      call plot(lll,refa,tlab,freq1,freq2,fname,fora,morc,scale)
      goto 555
      endif
      write(6,*)
      call coua@('Change any interface impedances ? (Y/N) ')
      call yesno(iq)
      if (iq.eq.1)then
998   write(6,*)
      call coua@('Enter interface number (1 is back, 0 to exit) ')
      read(5,*)in
      if (in.eq.0)goto 999
2006  write(6,*)
      write(6,*)' Choose type'
      write(6,*)' 1....parallel RC'
      write(6,*)' 2....series RC'
      write(6,*)' 3....parallel RL'
      write(6,*)' 4....series RL'
      read(5,*)ityp(in)
      if(ityp(in).lt.1.or.ityp(in).gt.4)goto 2006

      write(6,*)
      call coua@('Enter resistance value in ohms/sq ')
      read(5,*)sr
      gj(in)=1./sr

      if(ityp(in).lt.3)then
      write(6,*)
      call coua@(' Enter interface capacitance in fF/sq ')
      read(5,*)sc
      cj(in)=1.e-15*sc
      endif

      if(ityp(in).gt.2)then
      write(6,*)
      call coua@(' Enter interface inductance in nH/sq ')
      read(5,*)sc
      lj(in)=1.e-9*sc
      endif

      goto 998

999   goto 950
      endif
      close(10)
      close(11)

      stop
      end

      subroutine yesno(ia)
      integer*2 ik,ia
      ia=-1
1     call get_key@(ik)
      if(ik.eq.89.or.ik.eq.121)ia=1
      if(ik.eq.78.or.ik.eq.110)ia=0
      if(ia.eq.-1)goto 1
      write(6,101)char(ik)
101   format(a1)
      return
      end


      SUBROUTINE PLOT(lll,ref,tl,fstart,fstop,name,fora,morc,scale)
      dimension icol(8),ref(20000),tl(20000),refp(20000),tlp(20000)
      integer c1,c2,c3,c4,c5,iscale,bcol,bcol2,bcol3,fcol
      integer*2 icont,fora,morc
      character*3 xsf,xff,cnum,chscale
      character*8 edate,date
      character *8 name

      LLLM=LLL-1
      if(lllm.lt.1)then
      write(6,*)
      write(6,*)'NO DATA TO BE PLOTTED'
      call get_key@(icont)
      return
      endif

22    iscale=nint(scale)
      write(chscale,7654)iscale
7654  format(i3)
      date=edate()
      do2 i=1,lll
      if(ref(i).lt.scale)then
      refp(i)=scale
      else
      refp(i)=ref(i)
      endif

      if(tl(i).lt.scale)then
      tlp(i)=scale
      else
      tlp(i)=tl(i)
      endif

2     continue

C     set up start and stop frequencies on x axis

      if(fstart.lt.10.0.and.fstop.lt.10.0)then
      write(xsf,10)fstart
      write(xff,10)fstop
10    format(f3.1)
      else
      xsf=cnum(nint(fstart))
      xff=cnum(nint(fstop))
      endif

C     mono graphics
      if(morc.eq.1)then
      bcol=0          ! black
      bcol2=0         ! black
      bcol3=0         ! black
      fcol=15
      c1=15           ! white
      c2=15
      c3=8            ! dark grey
      c4=8            ! dark grey
      c5=15
      icol(1)=15
      icol(2)=7
      icol(3)=15
      icol(4)=15
      icol(5)=15
      icol(6)=7
      icol(7)=15      ! white
      icol(8)=15
      endif

C     16 colour graphics
      if(morc.eq.2)then
      bcol=8          ! dark grey
      bcol2=0         ! black
      bcol3=9         ! light blue
      fcol=7          ! light grey
      c1=10           ! light turquoise
      c2=10           ! light green
      c3=0            ! black
      c4=8            ! dark grey
      c5=12           ! light red
      icol(1)=12      ! red
      icol(2)=9       ! blue
      icol(3)=10      ! green
      icol(4)=6       ! brown
      icol(5)=13      ! purple
      icol(6)=11      ! turquoise
      icol(7)=14      ! yellow
      icol(8)=15      ! white
      endif

      call VGA@
      call clear_screen_area@(0,0,639,479,bcol)
      call clear_screen_area@(60,40,559,440,bcol2)
      call clear_screen_area@(0,0,639,37,bcol3)
      call clear_screen_area@(0,36,639,38,bcol2)

      call set_text_attribute@(0,0.0,0.0,0.0)
      call draw_text@(name,60,4,c2)
      call draw_text@(date,495,22,c2)
      call draw_text@('0',35,38,c2)
      call draw_text@(chscale,20,428,c2)

      if(fora.eq.1)then
      CALL draw_text@(xsf,61,455,c2)
      call draw_text@('Frequency (GHz)',250,455,c2)
      CALL draw_text@(xff,545,455,c2)
      else
      call draw_text@('0',61,455,c2)
      call draw_text@('Angle of incidence (deg)',215,455,c2)
      call draw_text@('90',545,455,c2)
      endif

      call set_text_attribute@(0,0.0,90.0,0.0)
      call draw_text@('Reflectivity',40,375,c5)
      call draw_text@('/Transmission Loss (dB)',40,275,c2)

      call set_text_attribute@(0,0.0,0.0,0.0)

      DO9876 IR=1,10
      ifir=440-nint(40.*FLOAT(IR))
      call draw_line@(60,ifir,559,ifir,c4)
9876  CONTINUE

      if(scale.eq.-30.0)then
      do98766 ir=1,5
c      if(ir.eq.3)goto 98766
      ifir=440-nint(400.*float(ir)/6.)
      icol2=c4
      do 98767 jj=1,100
      jst=5*(jj-1)+60
      jsf=5*jj+60
      call draw_line@(jst,ifir,jsf,ifir,icol2)
      if(icol2.eq.c4)then
      icol2=bcol2
      else
      icol2=c4
      endif
98767 continue
98766 continue
      endif

c     sgt up number of frequency axis graticule lines on plot

      if(fora.eq.1)then
      xinc=(fstop-fstart)
      if(xinc.le.1.or.xinc.gt.20.)xinc=10.
      endif

c     set up number of angle axis graticule lines on plot

      if(fora.eq.2)xinc=18.

      nxinc=int(xinc+0.00001)
      DO9877 IR=1,NXINC
      IIX=60+nint(499.*FLOAT(IR)/XINC)
      call draw_line@(IIX,440,IIX,40,c4)
9877  CONTINUE

      call draw_line@(560,440,60,440,c3)
      call draw_line@(60,440,60,40,c3)
      call draw_line@(60,40,560,40,c3)
      call draw_line@(560,40,560,440,c3)

      dx=499./float(LLLM)
      dy=abs(400./scale)

      DO7 I=1,LLLM
      ix1=60+nint(dx*float(i-1))
      ix2=60+nint(dx*float(i))
      iy1=40-nint(refp(i)*dy)
      iy2=40-nint(refp(i+1)*dy)
      call draw_line@(ix1,iy1,ix2,iy2,icol(1))

7     CONTINUE

      DO77 I=1,LLLM
      ix1=60+nint(dx*float(i-1))
      ix2=60+nint(dx*float(i))
      iy1=40-nint(tlp(i)*dy)
      iy2=40-nint(tlp(i+1)*dy)
      call draw_line@(ix1,iy1,ix2,iy2,icol(3))

77     CONTINUE


      call get_key@(icont)

      if(icont.eq.67.or.icont.eq.99)then
      if(morc.eq.1)then
      morc=2
      else
      morc=1
      endif
      goto22
      endif

      call text_mode@
      return
      END
