c       um-maad program from Gonzalez
        complex s11, s21, s12, s22, delta, c1, c2
        call enter( s11, s21, s12, s22, delta, dm, da, ak, c1, c2 )
1       continue
	call select( n )
	go to ( 10, 20, 30, 40, 50, 60, 70, 80, 90), n
10      call stabil( s11, s21, s12, s22, dm, da, ak, c1, c2 )
     	goto 1
20	call merit( s11, s21, s12, s22 )
	goto 1
30	call ungain( s11, s21, s22 )
	goto 1
40	call simat( s11, s22, dm, ak, c1, c2 )
	goto 1
50	call powerg( s11, s21, s12, s22, dm, ak, c1, c2 )
	goto 1
60	call noisef
	goto 1
70	call reflec( s11, s21, s12, s22 )
	goto 1
80	call conver
	goto 1
90	stop
	end


	subroutine enter( s11, s21, s12, s22, delta, dm, da, 
     +                    ak, c1, c2 )
        complex s11, s21, s12, s22, delta, c1, c2
9       print *,' Enter Mag. and Ang. of S11, S21, S12, S22'
	read *, s11m, s11a, s21m, s21a, s12m, s12a, s22m, s22a
	s11 = cmplx( s11m * cos( rad( s11a ) ), 
     +               s11m * sin( rad( s11a ) ) )
	s21 = cmplx( s21m * cos( rad( s21a ) ), 
     +               s21m * sin( rad( s21a ) ) )
	s12 = cmplx( s12m * cos( rad( s12a ) ), 
     +               s12m * sin( rad( s12a ) ) )
	s22 = cmplx( s22m * cos( rad( s22a ) ), 
     +               s22m * sin( rad( s22a ) ) )
        write( 6, 19 )
19      format( 6x, 's11=', 16x, 's21=', 16x,
     +          's12=', 16x, 's22=')
	write( 6, 29 ) s11m, s11a, s21m, s21a, s12m, s12a,
     +                 s22m, s22a
29      format( 8( 2x, f8.3 ) )
        delta = s11 * s22 - s12 * s21
	dm = cabs( delta )
	da = cang( delta )
	da = ang( da )
	c1 = s11 - delta * conjg( s22 )
	c2 = s22 - delta * conjg( s11 )
	h1 = cabs( s21 * s12 )
	if( h1 .eq. 0. ) goto 39
	ak = ( 1 - cabs( s11 ) **2 - cabs( s22 ) **2 
     +       + dm **2 ) / ( 2 * h1 )
        goto 49
39      ak = 1.e20
49      return
	end



	subroutine select( n )
	write( 6, 11 )
11      format( // '    Enter 1 for stability cirecles, k & delta '/
     +             '          2 for unilateral figure of merit' /
     +             '          3 for unilateral power gain circles'/
     +             '          4 for simultaneous conjugate match'/
     +             '          5 for power & avail. power gain circles'/
     +             '          6 for noise figure circles'/
     +             '          7 for i/o reflection coefficients'/
     +             '          8 for conversion between gamma & Z'/
     +             '          9 to Exit  ')
	read *, n
	return
	end



	subroutine stabil( s11, s21, s12, s22, dm, da, ak, c1, c2 )
	complex s11, s12, s21, s22, cl, cs, c1, c2
	write( 6, 21 ) dm, da, ak
21      format( //, 5x, 'Delta.Mag= ', f7.3, 5x, 'Delta.Ang= ',
     +          f8.3, 5x, '     k= ', f7.3,//)
	rld= cabs( s22 ) **2 - dm **2
	rl = cabs( s21 * s12 / rld )
	cl = conjg( c2 ) / rld
	clm = cabs( cl )
	cla = ang( cang( cl ) )
	rsd = cabs( s11 ) **2 - dm **2
	rs = cabs( s12 * s21 / rsd )
	cs = conjg( c1 ) / rsd
	csm = cabs( cs )
	csa = ang( cang( cs ) )
	if ( ( dm .lt. 1. ) .and. ( ak .ge. 1. ) ) goto 101
	print *,'  Potentially unstable transistor'
	print *,'  ==============================='
	goto 111
101     print *,'  Unconditionally stable transistor'
        print *,'  ================================'
111     print *,' **** Input stability circle **** '
        print *
	print *,'Radius = ', rs, 'Center.mag= ', csm, ' Center.ang= ', csa
	if( cabs( s22 ) - 1. ) 41, 51, 51
41      if( csm .gt. rs ) goto 121
        print *,' inside is stable '
	goto 131 
121     print *,' outside is stable'
131     goto 61
51      if( csm .gt. rs ) goto 141
	print *,' outside is stable'
	goto 61
141     print *,' inside is stable'
61      print *
	print *,' **** Output stability circle **** '
	print *
	print *,' radius= ', rl, 'center.mag= ', clm, 'center.ang= ', cla
	if( cabs( s11 ) - 1. ) 81, 91, 91
81      if( clm .gt. rl ) goto 211
	print *,' inside is stable '
	goto 221
211     print *,'  outside is stable'
221     goto 301
91      if( clm .gt. rl ) goto 231
	print *,' outside is stable'
	goto 301
231     print *,' inside is stable'
301     return
	end



       
       subroutine merit( s11, s21, s12, s22 )
       complex s11, s21, s12, s22
       un= cabs( s12 ) * cabs( s21 ) * cabs( s11 ) * cabs( s22 )
       ud= ( 1 - cabs( s11 ) **2 ) * ( 1 - cabs( s22 ) **2 )
       u = un / ud
       ul = 10 * alog10( ( 1 + u ) **(-2) )
       ur = 10 * alog10( ( 1 - u ) **(-2) )
       print *
       print *,' **** Unilateral figure of merit = ', u
       print *
       print *, ul, ' dB < GT/GTU < ', ur, ' dB'
       return
       end



       subroutine ungain( s11, s21, s22 )
       complex s11, s21, s22, s
       if( ( cabs( s11 ) .ge. 1. ) .or. ( cabs( s22 ) .ge. 1. ) ) 
     +    goto 82
       gsm = 10 * alog10( 1 / ( 1 - cabs( s11 ) **2 ) )
       glm = 10 * alog10( 1 / ( 1 - cabs( s22 ) **2 ) )
       go = 10 * alog10( cabs( s21 ) **2 )
       print *, ' s21, cabs(s21), log ', s21, cabs(s21), 
     +         log10( cabs(s21)**2 )
       gumax = gsm + glm + go
       print *
       print *,'   gsmax(dB)= ', gsm,'   glmax(dB)= ', glm
       print *,'   go   (dB)= ', go ,'   gtumax(dB)= ', gumax
       print *
       goto 2
82     print *
       print *,'   potentially unstable transistor'
       print *
2      print *,' enger 1 for Gs circle, 2 for Gl circle, 3 to exit'
       read *, m
       if( m .eq. 3 ) goto 92
       print *,' enter unilateral gain dBs -- start, stop, step '
       read *, start, stop, step
       goto ( 12, 22 ) m
       print *,' try again --- error!'
       goto 2
12     s = s11
       print *
       print *,'   **** Gs circles **** '
       goto 32
22     s = s22
       print *
       print *,'   **** Gl circles **** '
32     ainc = ang( cang( 1/ s ) )
       g= start
42     gi = 10 ** ( g / 10 )
       gi = gi * ( 1 - cabs( s ) **2 )
       risd = 1 - cabs( s ) **2 * ( 1 - gi )
       risn = sqrt( 1 - gi ) * ( 1 - cabs( s) **2 )
       ris = risn / risd
       dis = gi * cabs( s ) / ( 1 - cabs( s ) **2 * ( 1 - gi ) )
       print *
       print *, g, ' dB gain circle ', ' radius = ', ris
       print *,'   center.mag = ', dis,' center.ang = ', ainc
       print *
       g = g + step
       if( g .le. stop ) goto 42
       goto 2
92     return
       end



	function cang( w )
	complex w
	x = real( w )
	y = aimag( w )
	cang= atan2( y, x )
	return
	end



	function ang( rad )
	pi= 3.14159265
	ang= rad * 180 / pi
	return
	end



	function rad( ang )
	pi= 3.14159265
	rad = ang * pi / 180
	return
	end



	subroutine simat( s11, s22, dm, ak, c1, c2 )
	complex s11, s22, c1, c2, gammas, gammal
	if( ( dm .lt. 1. ) .and. ( ak .gt. 1. ) ) goto 37
	print *,' pot. unstable transistor, use operating or'
        print *,' available power gain design'
        goto 47
37      b1 = 1 + cabs( s11 ) **2 - cabs( s22 ) **2 - dm **2
        b2 = 1 + cabs( s22 ) **2 - cabs( s11 ) **2 - dm **2
	gammas = ( b1 - sqrt( b1 **2 - 4 * cabs( c1 ) **2 ) ) / 2 / c1
	gammal = ( b2 - sqrt( b2 **2 - 4 * cabs( c2 ) **2 ) ) / 2 / c2
	gsmag= cabs( gammas )
	gsang = ang( cang( gammas ) )
	glmag = cabs( gammal )
	glang = ang( cang( gammal ) )
	print *
	print *,' **** complex conj. match at the input port ****'
	print *,' **** GAMMAs.mag= ', gsmag, '  GAMMAS.ang= ', gsang
	print *
	print *,' **** Compl. conj. match at the output port ****'
	print *,' **** GAMMAl.mag= ', glmag, '  GAMMAl.ang= ', glang
47      return
	end




	subroutine powerg( s11, s21, s12, s22, dm, ak, c1, c2 )
	complex s, s11, s21, s12, s22, c1, c2, c, center
	print *
	if( ( ak .gt. 1. ) .and. ( dm .lt. 1. ) ) goto 93
	gmsg = cabs( s21 / s12 )
	gmsg = 10 * alog10( gmsg )
        print *,'  Potentially unst. transistor, MSG = ', gmsg, ' dB'
	goto 3
93      gpmax = ( ak - sqrt( ak **2 - 1 ) ) * cabs( s21 / s12 )
	gpmax = 10 * alog10( gpmax )
	print *,'  uncond. stable tran. gpmax= ', gpmax, ' dB'
3       print *
	print *,'  enter 1 for operating power gain circles'
	print *,'        2 for available gain circles'
	print *,'        3 for exit '
	read *, m
	if( m .eq. 3 ) goto 103
	print *,' enter gain in dBs --- start, stop, step '
	read *, start, stop, step
	goto ( 13, 23 ) m
	print *,' try again'
	goto 3
13      s = s22
	c = c2
	print *,' **** Operating power gain circles ****'
	goto 33
23      s = s11
	c = c1
	print *,' **** Available power gain circles ****'
33      tmp1 = cabs( s12 * s21 )
	tmp2 = cabs( s ) **2 - dm **2
	g = start
43      gi = 10 ** ( g / 10 )/ cabs( s21 ) **2
	center = gi * conjg( c ) / ( 1 + gi * tmp2 )
	distan = cabs( center )
	angle = ang( cang( center ) )
	rn = sqrt( 1 - 2 * ak * tmp1 * gi + ( tmp1 * gi ) **2 )
	print *,' ak, tmp1, tmp2, gi', ak, tmp1, tmp2, gi
	tt1= cabs(c)
	tt2= ang(cang(c))
	print *,' c, gi/(1+gi*tmp2)', tt1, tt2, gi/(1+gi*tmp2)
	radius = rn / ( 1 + gi * tmp2 )
	print *
	print *, g, ' dB gain Circle.  radius= ', radius
	print *,'   center.mag = ', distan,'   center.ang = ', angle
	g = g + step
	if ( g .le. stop ) goto 43
	goto 3
103     return
	end




	subroutine noisef
	complex go
	print *,' enter fimin(dB), go.mag, go.ang, rn '
	read *, fmin, gomag, goang, rn
	fmin = 10 ** ( fmin / 10 )
	r = rn / 50
4       print *
	print *,' enter noise fig. in dBs --- start, stop, step '
	read *, start, stop, step
	f = start
14      fi = 10 ** ( f / 10 )
	go = cmplx( gomag * cos( rad( goang ) ), 
     +              gomag * sin( rad( goang ) ) )
	tmp = cabs( 1 + go ) **2
	ani = ( fi - fmin ) * tmp / ( 4 * r )
	cfimag = gomag / ( 1 + ani )
	cfiang = goang
	rfi = sqrt( ( ani ) **2 + ani * ( 1 - gomag **2 ) ) / ( 1 + ani )
	print *, f, ' dB noise figure circle'
	print *
	print *,'  center.mag= ',cfimag,'  center.ang= ', cfiang
	print *,'    radius= ', rfi,' Ni = ', ani,' tmp = ', tmp
	print *
	f = f + step
	if( f .le. stop ) goto 14
	print *,'  enter 1 to continue, 0 to exit '
	read *, j
	if ( j .eq. 1 ) goto 4
	return
	end


	subroutine reflec( s11, s21, s12, s22 )
	complex s11, s21, s12, s22, sn, sd, gamma, g
	print *
5       print *,' enter  1 for input reflection coefficient '
	print *,'        2 for output reflection coefficient '
	print *,'        3 for exit '
	read *, l
	if( l .eq. 3 ) goto 55
	print *,'  enter gamma termination (mag,ang) '
	read *, gammag, gamang
	gamang= rad( gamang )
	gamma = cmplx( gammag * cos( gamang), gammag * sin( gamang ) )
	goto( 25, 35 ) l
	print *,'  try again '
	goto 5
25      sn = s11
        sd = s22
	print *,'   **** the input reflection coefficient **** '
	goto 45
35      sn= s22
	sd= s11
	print *,'   **** the output reflection coefficient **** '
45      g= sn + s12 * s21 * gamma / ( 1 - sd * gamma )
	gmag= cabs( g )
	gang= ang( cang( g ) )
	print *,'  gamma.mag = ', gmag,' gamma.ang = ', gang
	print *
	goto 5
55      return
	end



	subroutine conver
	complex z, zz, yy, gamma, gam
6	print *
	print *,'    enter 1 for Z to Gamma'
	print *,'          2 for Gamma to Z and Y'
	print *,'          3 for exit'
	read *, i
	if( i .eq. 3 ) goto 56
	goto ( 26, 36 ) i
	print *,'  try again'
	goto 6
26      print *,' enter normalized Z (real, imag)'
	read *, zre, zim
	z= cmplx( zre, zim )
	gamma= ( z - 1 ) / ( z + 1 )
	gmag= cabs( gamma )
	gang= ang( cang( gamma) )
	print *, '  gamma.mag = ', gmag, '  gamma.ang = ', gang
	print *
	goto 46
36      print *,'   enter gamm (mag, ang ) '
	read *, gamag, gaang
	gam = cmplx( gamag * cos( rad( gaang) ), 
     +     gamag * sin( rad( gaang ) ) )
	zz = ( 1 + gam ) / ( 1 - gam )
	yy = 1 / zz
	print *,'   impedance z = ', zz
	print *,'   admitance y = ', yy
	print *
46      goto 6
56      return
	end
