      program smith
*
*  Written by: Louis Wilson
*              louis@wizard.lasc.lockheed.com
*  Program to plot Impdeance as a function of frequency
*  on a smith chart.
*  Input file is a title line followed by
*  data pairs of frequency and complex impedance
*
*  Warning - calls sind, cosd, atan2d degree argument
*            trigonometric functions. I your fortran
*            does not have these you will have to write
*            them.
*
      character*45 title1,cdate
      complex er
      real r(200),phi(200),freq(200)
      common /stuff/ er,r,np,phi,freq,n,cdate,title1
*
      call read_file
*
      call setup
      call smith(r,phi,freq,n,1)
      call finish
*
      end
*--------------------------------------------------------------------*
      subroutine read_file
      character*45 title1,cdate
      character*80 fname
      complex er,z(16),gam
      real r(200),phi(200),freq(200)
      common /stuff/ er,r,np,phi,freq,n,cdate,title1
*
*  Sun fortran date/time call - returns a 24 character string
*  with date and time.
*     call fdate(cdate)
      cdate =' '
      print*,'enter the data file name'
      read(*,'(a)')fname
      open(8,file=fname,status='old')
      i = 1
      read(8,'(a)')title1
1     read(8,*,err = 2,end = 2)freq(i),z(i)
      i = i + 1
      goto 1
2     continue
      n = i - 1
      do i = 1, n
         gam    =  (z(i)-1)/(z(i)+1)
         r(i)   =  cabs(gam)
         xr = real(gam)
         if(xr.eq.0) xr = 1e-6
         xi = aimag(gam)
         phi(i) =  atan2d( xi,xr )
      enddo
      return
      end
*--------------------------------------------------------------------*
      swbroutine smith(r,phi,freq,n,imark)
*
*  Plot a Smith Chart
*
*  Inputs:
*     r(n)    =   magnitude of the reflection coefficient (unitless)
*     phi(n)  =   phase angle of the reflection coefficient (degrees)
*     freq(n) =   frequency of each element of r and phi (ghz)
*     n       =   number of data points
*
      real  phi(n),r(n),freq(n)
*
      xdist = 4.25
      ydist = 4.0
      radius = 3.0
*
*   Plot grid
*
      call smgrid(xdist,ydist,radius)
*
*   Plot and label curve
*
      call reg_titles
      call polar(xdist,ydist)
      write(10,*)'1.0 setgray'
      call pcurve(r, phi, n, radius,3.0)
      write(10,*)'0.0 setgray'
      call pcurve(r, phi, n, radius,1.0)
*
      return
      end
*-----------------------------------------------------------------*
      subroutine smgrid(xdist,ydist,r)
*
*  Plot the Smith Chart Grid Lines
*
*  Inputs:
*     xdist =     x offset (in inches) of smith chart origin
*     ydist =     y offset (in inches) of smith chart origin
*     r =         radius (in inches) of smith chart
*
      real xdist,ydist, r,z(6)
      data z/0.,.2,.5,1.,2.,5./
*
*  Zero Reactance line
*
      call line (xdist+r,ydist,xdist-r,ydist)
*
*  Resistance circles
*
      do i=6,1,-1
         r1  = i*r/6 
         call circle(xdist - r1 + r,ydist,r1)
      enddo
      a1 = 90.
      a2 = 270.
*
*  Reactance circles
*
      do i= 2,6
         rc = r/z(i)
         call get_angle(a3,cmplx(0.,z(i)),r,rc)
         call arc(xdist + r,ydist - rc,rc,a1,a1 + a3)
         call arc(xdist + r,ydist + rc,rc,a2-a3,a2)
      enddo
*
*  plot resistance labels
*
      call setfont(6)
      dr = r/3.
      ycen = ydist
      do i = -3,2,1
         j = 4 + i
         if(z(j).eq.0)then
            xcen = xdist + i*dr +.05
           call label(int(z(j)),xcen,ycen + .05,0.0)
         else
            xcen = xdist + i*dr -.12
           call label(z(j),xcen,ycen,0.0)
         endif
      enddo
*
*  plot reactance labels
*
      rr = r - .12
      do i=2,6
*
*  Upper labels
*
         call get_rot(a,cmplx(0.,z(i)))
         xcen = xdist + rr*cosd(a)
         ycen = ydist + rr*sind(a)
         call label(z(i),xcen,ycen,a)
*
*  Lower labels
*
         xcen = xdist + r*cosd(a)
         ycen = ydist - r*sind(a)
         call label(z(i),xcen,ycen,-a + 180)
      enddo
      return
      end
*--------------------------------------------------------------------*
      subroutine reg_titles
      complex er
      character*45 title1,title2,title3,title4,title5, cdate
      real r(200),phi(200),freq(200)
      character*45 getlog
      common /stuff/ er,r,np,phi,freq,n,cdate,title1
*
      title2 = 'Smith Chart - normalized coordinates'
      title3 = 'Impedance Plot'
      title4 = cdate
*
* Sun fortran function to get login name
*
*     title5 = getlog()
      title5 = ' '
*
*   Label polarization
*
*
*   Label angle of incidence
*
      call setfont(20)
      call text (title1, 1.00, 9.75)
      call text (title2, 1.00, 9.50)
      call text (title3, 1.00, 9.25)
      call setfont(10)
      call text (title4, 1.00, 0.50)
      call text (title5, 1.00, 0.35)
      return
      end
*-----------------------------------------------------------------*
      subroutine get_rot(a,z)
      complex z,rho
      rho = (z-1)/(z+1)
      a = atan2d(aimag(rho),real(rho))
      return
      end
*-----------------------------------------------------------------*
      subroutine get_angle(ang,z,r,rc)
*
*  gets angle for plotting an arc
*
      complex z,rho
      rho = (z-1)/(z+1)
      rl = r*sind(atan2d(aimag(rho),real(rho))/2)
      ang = 2*asind(rl/rc)
      return
      end
*-----------------------------------------------------------------*
*  PostScript Routines
*-----------------------------------------------------------------*
      subroutine setup
*
*  This subroutine sets up the Post-Script header
*  and startup definitions
*
      open(10,file='smith.ps',status='unknown')
      write(10,*)'%!PS-Adobe-2.0'
      write(10,*)'%%Title: Smith Chart'
      write(10,*)'%%Creator: Louis Wilson'
      write(10,*)'%%BoundingBox: 0 0 612 792'
      write(10,*)'%%EndComments'
      write(10,*)'%%BeginProlog'
      write(10,*)'/Circle{'
      write(10,*)'    newpath 1 361 arc'
      write(10,*)'    stroke}def'
      write(10,*)'/Line{'
      write(10,*)'     newpath moveto'
      write(10,*)'     lineto'
      write(10,*)'     stroke }def'
      write(10,*)'%%EndProlog'
      write(10,*)'%%BeginScript'
      write(10,*)'.1 setlinewidth'
      write(10,*)'0 0 moveto'
      return
      end
*-----------------------------------------------------------------*
      subroutine circle(x,y,r)
      real x,y,r
      xx = x*72
      yy = y*72
      rr = r*72
      write(10,10)xx,yy,rr
 10   format(f8.2,2x,f8.2,2x,f8.2,2x,'Circle')
      return
      end
*-----------------------------------------------------------------*
      subroutine arc(x,y,r,a,b)
      xx = x*72
      yy = y*72
      rr = r*72
      write(10,*)'newpath'
      write(10,10)xx,yy,rr,a,b
      write(10,*)'stroke'
 10   format(f8.2,2x,f8.2,2x,f8.2,2x,f7.2,2x,f7.2,2x,'arc')
      return
      end
*-----------------------------------------------------------------*
      subroutine finish
      write(10,*)'%%EndScript'
      write(10,*)'showpage'
      close(10)
      return
      end
*-----------------------------------------------------------------*
      subroutine line(x1,y1,x2,y2)
      xx1 = x1*72
      xx2 = x2*72
      yy1 = y1*72
      yy2 = y2*72
      write(10,10)xx1,yy1,xx2,yy2 
10    format(f7.2,2x,f7.2,2x,f7.2,2x,f7.2,2x,'Line')
      return
      end
*-----------------------------------------------------------------*
      subroutine label(a,x,y,r)
      xx = x*72
      yy = y*72
      write(10,10)xx,yy
      write(10,13)r
      if(a.eq.0)then
         write(10,11)int(a)
      else
         write(10,12)a
      endif
      write(10,13)-r
 10   format(f7.2,2x,f7.2,2x,'moveto')
 11   format('(',i1,')',2x,'show')
 12   format('(',f3.1,')',2x,'show')
 13   format(f7.2,2x,'rotate')
      return
      end
*-----------------------------------------------------------------*
      subroutine setfont(iscale)
      write(10,10)iscale
 10   format('/Helvetica findfont',1x,,i2,1x,'scalefont setfont')
      return
      end
*-----------------------------------------------------------------*
      subroutine text(a,x,y)
*
*  Prints character string at location specified
*  by x and y.
*
      character*(*) a
      xx = x*72
      yy = y*72
      call count(a,n)
      write(10,10)xx,yy
      write(10,*)'(',a(1:n),') show'
 10   format(f7.2,2x,f7.2,2x,'moveto')
      return
      end
*-----------------------------------------------------------------*
      subroutine count(a,n)
*
*  Counts characters in a title
*
      character*(*) a
      character*1 b
      b = ' '
      nl = len(a)
 1    if(a(nl:nl).eq.b)then
         nl = nl - 1
         goto 1
      else
         goto 2
      endif
 2    n = nl
      return
      end
*-----------------------------------------------------------------*
      subroutine polar(x,y)
      xx = x*72
      yy = y*72
      write(10,10)xx,yy
 10   format(f7.2,2x,f7.2,2x,'translate')
      return
      end
*-----------------------------------------------------------------*
      subroutine pcurve(r,p,n,radius,xw)
      parameter(dtr = 3.14159265/180.)
      real r(*),p(*),radius
      write(10,*)'newpath'
      write(10,12)xw
      rad = radius*72
      ang = p(1)*dtr
      rr  = r(1)*rad
      x   = rr*cos(ang)
      y   = rr*sin(ang)
      write(10,10)x,y
      do i = 2, n
         ang = p(i)*dtr
         rr  = r(i)*rad
         x   = rr*cos(ang)
         y   = rr*sin(ang)
         write(10,11)x,y
      enddo
      write(10,*)'stroke'
 10   format(f7.2,2x,f7.2,2x,'moveto')
 11   format(f7.2,2x,f7.2,2x,'lineto')
 12   format(f5.2,2x,'setlinewidth')
      return
      end
