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 *--------------------------------------------------------------------* subroutine 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