+ subroutine acan
+ implicit double precision (a-h,o-z)
+c
+c
+c this routine drives the small-signal analyses.
+c
+ common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem,
+ 1 isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize,
+ 2 junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr,
+ 3 nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1,
+ 4 lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd,
+ 5 imynl,imvn,lcvn,loutpt,nsnod,nsmat,nsval,icnod,icmat,icval
+ common /miscel/ atime,aprog(3),adate,atitle(10),defl,defw,defad,
+ 1 defas,rstats(50),iwidth,lwidth,nopage
+ common /cirdat/ locate(50),jelcnt(50),nunods,ncnods,numnod,nstop,
+ 1 nut,nlt,nxtrm,ndist,ntlin,ibr,numvs
+ common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet,
+ 1 xmu,mode,modedc,icalc,initf,method,iord,maxord,noncon,iterno,
+ 2 itemno,nosolv,ipostp,iscrch
+ common /flags/ iprnta,iprntl,iprntm,iprntn,iprnto,limtim,limpts,
+ 1 lvlcod,lvltim,itl1,itl2,itl3,itl4,itl5,igoof,nogo,keof
+ common /knstnt/ twopi,xlog2,xlog10,root2,rad,boltz,charge,ctok,
+ 1 gmin,reltol,abstol,vntol,trtol,chgtol,eps0,epssil,epsox
+ common /ac/ fstart,fstop,fincr,skw2,refprl,spw2,jacflg,idfreq,
+ 1 inoise,nosprt,nosout,nosin,idist,idprt
+ common /cje/ maxtim,itime,icost
+ common /blank/ value(1000)
+ integer nodplc(64)
+ complex*16 cvalue(32)
+ equivalence (value(1),nodplc(1),cvalue(1))
+ data aendor /9.87654321d0/
+ call second(t1)
+c.. post-processor initialization
+ if(ipostp.eq.0) go to 1
+ numcur=jelcnt(9)
+ numpos=nunods+numcur
+ call getm16(ibuff,numpos)
+ numpos=numpos*4
+ if(numcur.eq.0) go to 1
+ loc=locate(9)
+ loccur=nodplc(loc+6)-1
+c
+c allocate storage
+c
+ 1 call getm8(ndiag,2*nstop)
+ call getm8(lvn,nstop)
+ call getm8(imvn,nstop)
+ call getm16(lcvn,nstop)
+ call getm8(lynl,nstop+nut+nlt)
+ call getm8(imynl,nstop+nut+nlt)
+ if (idist.ne.0) call dinit
+ nandd=0
+ if (inoise.eq.0) go to 10
+ if (idist.eq.0) go to 10
+ nandd=1
+ call getm16(lvnctc,nstop)
+ 10 call getm16(loutpt,0)
+ call crunch
+ lyu=lynl+nstop
+ lyl=lyu+nut
+ numout=jelcnt(43)+jelcnt(44)+jelcnt(45)+1
+c if(ipostp.ne.0) call pheadr(atitle)
+ icalc=0
+ freq=fstart
+c
+c load y matrix and c vector, solve for v vector
+c
+ 100 call getcje
+ if ((maxtim-itime).le.limtim) go to 900
+ omega=twopi*freq
+ call acload
+ call acdcmp
+ call acsol
+ if (igoof.eq.0) go to 200
+ write (6,121) igoof,freq
+ 121 format('0warning: underflow ',i4,' time(s) in ac analysis at freq
+ 1 = ',1pd9.3,' hz')
+ igoof=0
+c
+c store outputs
+c
+ 200 call extmem(loutpt,numout)
+ loco=loutpt+icalc*numout
+ icalc=icalc+1
+ cvalue(loco+1)=dcmplx(freq,omega)
+ loc=locate(43)
+ 310 if (loc.eq.0) go to 350
+ if (nodplc(loc+5).ne.0) go to 320
+ node1=nodplc(loc+2)
+ node2=nodplc(loc+3)
+ iseq=nodplc(loc+4)
+ cvalue(loco+iseq)=cvalue(lcvn+node1)-cvalue(lcvn+node2)
+ loc=nodplc(loc)
+ go to 310
+ 320 iptr=nodplc(loc+2)
+ iptr=nodplc(iptr+6)
+ iseq=nodplc(loc+4)
+ cvalue(loco+iseq)=cvalue(lcvn+iptr)
+ loc=nodplc(loc)
+ go to 310
+ 350 if(ipostp.eq.0) go to 400
+ cvalue(ibuff+1)=dcmplx(freq,0.0d0)
+ call copy16(cvalue(lcvn+2),cvalue(ibuff+2),nunods-1)
+ if(numcur.ne.0) call copy16(cvalue(lcvn+loccur+1),
+ 1 cvalue(ibuff+nunods+1),numcur)
+ call dblsgl(cvalue(ibuff+1),numpos)
+c call fwrite(cvalue(ibuff+1),numpos)
+c
+c noise and distortion analyses
+c
+ 400 if (nandd.eq.0) go to 410
+ call copy16(cvalue(lcvn+1),cvalue(lvnctc+1),nstop)
+ 410 if (inoise.ne.0) call noise(loco)
+ if (nandd.eq.0) go to 420
+ call copy16(cvalue(lvnctc+1),cvalue(lcvn+1),nstop)
+ 420 if (idist.ne.0) call disto(loco)
+c
+c increment frequency
+c
+ if (icalc.ge.jacflg) go to 1000
+ if (idfreq.ge.3) go to 510
+ freq=freq*fincr
+ go to 100
+ 510 freq=freq+fincr
+ go to 100
+c
+c finished
+c
+ 900 write (6,901)
+ 901 format('0*error*: cpu time limit exceeded ... analysis stopped'/)
+ nogo=1
+ 1000 if(ipostp.eq.0) go to 1010
+ cvalue(ibuff+1)=aendor
+c call fwrite(cvalue(ibuff+1),numpos)
+ if(ipostp.ne.0) call clrmem(ibuff)
+ 1010 call clrmem(lvnim1)
+ call clrmem(lx0)
+ call clrmem(lvn)
+ call clrmem(imvn)
+ call clrmem(lcvn)
+ call clrmem(imynl)
+ call clrmem(lynl)
+ call clrmem(ndiag)
+ if (idist.eq.0) go to 1020
+ call clrmem(ld0)
+ call clrmem(ld1)
+ 1020 if (nandd.eq.0) go to 1040
+ call clrmem(lvnctc)
+ 1040 call second(t2)
+ rstats(7)=rstats(7)+t2-t1
+ rstats(8)=rstats(8)+icalc
+ return
+ end
+ subroutine cdiv(xr,xi,yr,yi,cr,ci)
+c.. ok if cr and ci are really xr and xi or yr and yi
+ implicit double precision (a-h,o-z)
+ xrtemp=xr
+ xitemp=xi
+ yrtemp=yr
+ yitemp=yi
+ amag2=yrtemp*yrtemp+yitemp*yitemp
+ cr=(xrtemp*yrtemp+xitemp*yitemp)/amag2
+ ci=(xitemp*yrtemp-xrtemp*yitemp)/amag2
+ return
+ end
+ subroutine cmult(xr,xi,yr,yi,cr,ci)
+c.. ok if cr and ci are really xr and xi or yr and yi
+ implicit double precision (a-h,o-z)
+ xrtemp=xr
+ xitemp=xi
+ yrtemp=yr
+ yitemp=yi
+ cr=xrtemp*yrtemp-xitemp*yitemp
+ ci=xitemp*yrtemp+xrtemp*yitemp
+ return
+ end
+ subroutine dblsgl(carray,nwds)
+c.. note that carray really contains complex*16
+ complex*8 carray(1)
+ complex*8 ctemp8(2),ctemp3
+ complex*16 ctemp
+ equivalence (ctemp,ctemp8(1))
+c
+c.. gather up numpos 16-bit words from dbl-complex to
+c.. make up sgl-complex
+c
+ num=nwds/4
+ do 10 i=1,num
+ ndex=2*i-1
+ ctemp8(1)=carray(ndex)
+ ctemp8(2)=carray(ndex+1)
+ ctemp3=ctemp
+ carray(i)=ctemp3
+ 10 continue
+ return
+ end
+ subroutine acdcmp
+ implicit double precision (a-h,o-z)
+c
+c this routine performs an lu factorization of the circuit equation
+c coefficient matrix.
+c
+ common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem,
+ 1 isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize,
+ 2 junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr,
+ 3 nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1,
+ 4 lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd,
+ 5 imynl,imvn,lcvn,loutpt,nsnod,nsmat,nsval,icnod,icmat,icval
+ common /cirdat/ locate(50),jelcnt(50),nunods,ncnods,numnod,nstop,
+ 1 nut,nlt,nxtrm,ndist,ntlin,ibr,numvs
+ common /flags/ iprnta,iprntl,iprntm,iprntn,iprnto,limtim,limpts,
+ 1 lvlcod,lvltim,itl1,itl2,itl3,itl4,itl5,igoof,nogo,keof
+ common /knstnt/ twopi,xlog2,xlog10,root2,rad,boltz,charge,ctok,
+ 1 gmin,reltol,abstol,vntol,trtol,chgtol,eps0,epssil,epsox
+ common /blank/ value(1000)
+ integer nodplc(64)
+ complex*16 cvalue(32)
+ equivalence (value(1),nodplc(1),cvalue(1))
+c
+c
+ do 100 i=2,nstop
+ io=nodplc(iorder+i)
+ gdiag=dabs(value(lynl+io))+dabs(value(imynl+io))
+ if (gdiag.ge.gmin) go to 10
+ value(lynl+io)=gmin
+ value(imynl+io)=0.0d0
+ igoof=igoof+1
+ 10 jstart=nodplc(ilc+i)
+ jstop=nodplc(ilc+i+1)-1
+ if (jstart.gt.jstop) go to 100
+ do 90 j=jstart,jstop
+ call cdiv(value(lyl+j),value(imynl+nstop+nut+j),value(lynl+io),
+ 1 value(imynl+io),value(lyl+j),value(imynl+nstop+nut+j))
+ icol=nodplc(ilr+j)
+ kstart=nodplc(iur+i)
+ kstop=nodplc(iur+i+1)-1
+ if (kstart.gt.kstop) go to 90
+ do 80 k=kstart,kstop
+ irow=nodplc(iuc+k)
+ if (icol-irow) 20,60,40
+c
+c find (icol,irow) matrix term (upper triangle)
+c
+ 20 l=nodplc(iur+icol+1)
+ 30 l=l-1
+ if (nodplc(iuc+l).ne.irow) go to 30
+ ispot=lyu+l
+ ispot2=imynl+nstop+l
+ go to 70
+c
+c find (icol,irow) matrix term (lower triangle)
+c
+ 40 l=nodplc(ilc+irow+1)
+ 50 l=l-1
+ if (nodplc(ilr+l).ne.icol) go to 50
+ ispot=lyl+l
+ ispot2=imynl+nstop+nut+l
+ go to 70
+c
+c find (icol,irow) matrix term (diagonal)
+c
+ 60 ispot=lynl+nodplc(iorder+irow)
+ ispot2=imynl+nodplc(iorder+irow)
+c
+ 70 call cmult(value(lyl+j),value(imynl+nstop+nut+j),
+ 1 value(lyu+k),value(imynl+nstop+k),xreal,ximag)
+ value(ispot)=value(ispot)-xreal
+ value(ispot2)=value(ispot2)-ximag
+ 80 continue
+ 90 continue
+ 100 continue
+ return
+ end
+ subroutine acsol
+ implicit double precision (a-h,o-z)
+c
+c this routine solves the circuit equations by performing a forward
+c and backward substitution using the previously-computed lu factors.
+c
+ common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem,
+ 1 isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize,
+ 2 junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr,
+ 3 nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1,
+ 4 lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd,
+ 5 imynl,imvn,lcvn,loutpt,nsnod,nsmat,nsval,icnod,icmat,icval
+ common /cirdat/ locate(50),jelcnt(50),nunods,ncnods,numnod,nstop,
+ 1 nut,nlt,nxtrm,ndist,ntlin,ibr,numvs
+ common /blank/ value(1000)
+ integer nodplc(64)
+ complex*16 cvalue(32)
+ equivalence (value(1),nodplc(1),cvalue(1))
+c
+c forward substitution
+c
+ do 20 i=2,nstop
+ jstart=nodplc(ilc+i)
+ jstop=nodplc(ilc+i+1)-1
+ if (jstart.gt.jstop) go to 20
+ io=nodplc(iorder+i)
+ if (value(lvn+io).ne.0.0d0) go to 5
+ if (value(imvn+io).eq.0.0d0) go to 20
+ 5 do 10 j=jstart,jstop
+ jo=nodplc(ilr+j)
+ jo=nodplc(iorder+jo)
+ call cmult(value(lyl+j),value(imynl+nstop+nut+j),value(lvn+io),
+ 1 value(imvn+io),xreal,ximag)
+ value(lvn+jo)=value(lvn+jo)-xreal
+ value(imvn+jo)=value(imvn+jo)-ximag
+ 10 continue
+ 20 continue
+c
+c back substitution
+c
+ k=nstop+1
+ do 50 i=2,nstop
+ k=k-1
+ io=nodplc(iorder+k)
+ jstart=nodplc(iur+k)
+ jstop=nodplc(iur+k+1)-1
+ if (jstart.gt.jstop) go to 40
+ do 30 j=jstart,jstop
+ jo=nodplc(iuc+j)
+ jo=nodplc(iorder+jo)
+ call cmult(value(lyu+j),value(imynl+nstop+j),value(lvn+jo),
+ 1 value(imvn+jo),xreal,ximag)
+ value(lvn+io)=value(lvn+io)-xreal
+ value(imvn+io)=value(imvn+io)-ximag
+ 30 continue
+ 40 call cdiv(value(lvn+io),value(imvn+io),value(lynl+io),
+ 1 value(imynl+io),value(lvn+io),value(imvn+io))
+ 50 continue
+ do 100 i=2,nstop
+ cvalue(lcvn+i)=dcmplx(value(lvn+i),value(imvn+i))
+ 100 continue
+ cvalue(lcvn+1)=dcmplx(0.0d0,0.0d0)
+ return
+ end
+ subroutine acload
+ implicit double precision (a-h,o-z)
+c
+c this routine zeroes-out and then loads the complex coefficient
+c matrix.
+c
+ common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem,
+ 1 isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize,
+ 2 junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr,
+ 3 nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1,
+ 4 lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd,
+ 5 imynl,imvn,lcvn,loutpt,nsnod,nsmat,nsval,icnod,icmat,icval
+ common /cirdat/ locate(50),jelcnt(50),nunods,ncnods,numnod,nstop,
+ 1 nut,nlt,nxtrm,ndist,ntlin,ibr,numvs
+ common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet,
+ 1 xmu,mode,modedc,icalc,initf,method,iord,maxord,noncon,iterno,
+ 2 itemno,nosolv,ipostp,iscrch
+ common /blank/ value(1000)
+ integer nodplc(64)
+ complex*16 cvalue(32)
+ equivalence (value(1),nodplc(1),cvalue(1))
+c
+c
+ complex*16 cval
+c
+c zero y matrix and current vector
+c
+ call zero8(value(lvn+1),nstop)
+ call zero8(value(imvn+1),nstop)
+ call zero8(value(lynl+1),nstop+nut+nlt)
+ call zero8(value(imynl+1),nstop+nut+nlt)
+c
+c resistors
+c
+ loc=locate(1)
+ 20 if (loc.eq.0) go to 30
+ locv=nodplc(loc+1)
+ val=value(locv+1)
+ locy=lynl+nodplc(loc+6)
+ value(locy)=value(locy)+val
+ locy=lynl+nodplc(loc+7)
+ value(locy)=value(locy)+val
+ locy=lynl+nodplc(loc+4)
+ value(locy)=value(locy)-val
+ locy=lynl+nodplc(loc+5)
+ value(locy)=value(locy)-val
+ loc=nodplc(loc)
+ go to 20
+c
+c capacitors
+c
+ 30 loc=locate(2)
+ 40 if (loc.eq.0) go to 50
+ locv=nodplc(loc+1)
+ val=omega*value(locv+1)
+ locyi=imynl+nodplc(loc+10)
+ value(locyi)=value(locyi)+val
+ locyi=imynl+nodplc(loc+11)
+ value(locyi)=value(locyi)+val
+ locyi=imynl+nodplc(loc+5)
+ value(locyi)=value(locyi)-val
+ locyi=imynl+nodplc(loc+6)
+ value(locyi)=value(locyi)-val
+ loc=nodplc(loc)
+ go to 40
+c
+c inductors
+c
+ 50 loc=locate(3)
+ 60 if (loc.eq.0) go to 70
+ locv=nodplc(loc+1)
+ val=omega*value(locv+1)
+ locyi=imynl+nodplc(loc+13)
+ locy=lynl+nodplc(loc+13)
+ value(locy)=0.0d0
+ value(locyi)=-val
+ locy=lynl+nodplc(loc+6)
+ locyi=imynl+nodplc(loc+6)
+ value(locy)=1.0d0
+ value(locyi)=0.0d0
+ locy=lynl+nodplc(loc+7)
+ locyi=imynl+nodplc(loc+7)
+ value(locy)=-1.0d0
+ value(locyi)=0.0d0
+ locy=lynl+nodplc(loc+8)
+ locyi=imynl+nodplc(loc+8)
+ value(locy)=1.0d0
+ value(locyi)=0.0d0
+ locy=lynl+nodplc(loc+9)
+ locyi=imynl+nodplc(loc+9)
+ value(locy)=-1.0d0
+ value(locyi)=0.0d0
+ loc=nodplc(loc)
+ go to 60
+c
+c mutual inductors
+c
+ 70 loc=locate(4)
+ 80 if (loc.eq.0) go to 90
+ locv=nodplc(loc+1)
+ val=omega*value(locv+1)
+ locy=lynl+nodplc(loc+4)
+ locyi=imynl+nodplc(loc+4)
+ value(locy)=0.0d0
+ value(locyi)=-val
+ locy=lynl+nodplc(loc+5)
+ locyi=imynl+nodplc(loc+5)
+ value(locy)=0.0d0
+ value(locyi)=-val
+ loc=nodplc(loc)
+ go to 80
+c
+c nonlinear voltage controlled current sources
+c
+ 90 loc=locate(5)
+ 95 if (loc.eq.0) go to 100
+ ndim=nodplc(loc+4)
+ lmat=nodplc(loc+7)
+ loct=lx0+nodplc(loc+12)+2
+ do 97 i=1,ndim
+ val=value(loct)
+ loct=loct+2
+ locy=lynl+nodplc(lmat+1)
+ value(locy)=value(locy)+val
+ locy=lynl+nodplc(lmat+2)
+ value(locy)=value(locy)-val
+ locy=lynl+nodplc(lmat+3)
+ value(locy)=value(locy)-val
+ locy=lynl+nodplc(lmat+4)
+ value(locy)=value(locy)+val
+ lmat=lmat+4
+ 97 continue
+ loc=nodplc(loc)
+ go to 95
+c
+c nonlinear voltage controlled voltage sources
+c
+ 100 loc=locate(6)
+ 105 if (loc.eq.0) go to 110
+ ndim=nodplc(loc+4)
+ lmat=nodplc(loc+8)
+ loct=lx0+nodplc(loc+13)+3
+ locy=lynl+nodplc(lmat+1)
+ locyi=imynl+nodplc(lmat+1)
+ value(locy)=+1.0d0
+ value(locyi)=0.0d0
+ locy=lynl+nodplc(lmat+2)
+ locyi=imynl+nodplc(lmat+2)
+ value(locy)=-1.0d0
+ value(locyi)=0.0d0
+ locy=lynl+nodplc(lmat+3)
+ locyi=imynl+nodplc(lmat+3)
+ value(locy)=+1.0d0
+ value(locyi)=0.0d0
+ locy=lynl+nodplc(lmat+4)
+ locyi=imynl+nodplc(lmat+4)
+ value(locy)=-1.0d0
+ value(locyi)=0.0d0
+ lmat=lmat+4
+ do 107 i=1,ndim
+ val=value(loct)
+ loct=loct+2
+ locy=lynl+nodplc(lmat+1)
+ value(locy)=value(locy)-val
+ locy=lynl+nodplc(lmat+2)
+ value(locy)=value(locy)+val
+ lmat=lmat+2
+ 107 continue
+ loc=nodplc(loc)
+ go to 105
+c
+c nonlinear current controlled current sources
+c
+ 110 loc=locate(7)
+ 115 if (loc.eq.0) go to 120
+ ndim=nodplc(loc+4)
+ lmat=nodplc(loc+7)
+ loct=lx0+nodplc(loc+12)+2
+ do 117 i=1,ndim
+ val=value(loct)
+ loct=loct+2
+ locy=lynl+nodplc(lmat+1)
+ locyi=imynl+nodplc(lmat+1)
+ value(locy)=+val
+ value(locyi)=0.0d0
+ locy=lynl+nodplc(lmat+2)
+ locyi=imynl+nodplc(lmat+2)
+ value(locy)=-val
+ value(locyi)=0.0d0
+ lmat=lmat+2
+ 117 continue
+ loc=nodplc(loc)
+ go to 115
+c
+c nonlinear current controlled voltage sources
+c
+ 120 loc=locate(8)
+ 125 if (loc.eq.0) go to 140
+ ndim=nodplc(loc+4)
+ lmat=nodplc(loc+8)
+ loct=lx0+nodplc(loc+13)+3
+ locy=lynl+nodplc(lmat+1)
+ locyi=imynl+nodplc(lmat+1)
+ value(locy)=+1.0d0
+ value(locyi)=0.0d0
+ locy=lynl+nodplc(lmat+2)
+ locyi=imynl+nodplc(lmat+2)
+ value(locy)=-1.0d0
+ value(locyi)=0.0d0
+ locy=lynl+nodplc(lmat+3)
+ locyi=imynl+nodplc(lmat+3)
+ value(locy)=+1.0d0
+ value(locyi)=0.0d0
+ locy=lynl+nodplc(lmat+4)
+ locyi=imynl+nodplc(lmat+4)
+ value(locy)=-1.0d0
+ value(locyi)=0.0d0
+ lmat=lmat+4
+ do 127 i=1,ndim
+ val=value(loct)
+ loct=loct+2
+ locy=lynl+nodplc(lmat+i)
+ value(locy)=value(locy)-val
+ 127 continue
+ loc=nodplc(loc)
+ go to 125
+c
+c voltage sources
+c
+ 140 loc=locate(9)
+ 150 if (loc.eq.0) go to 160
+ locv=nodplc(loc+1)
+ iptr=nodplc(loc+6)
+ value(lvn+iptr)=value(locv+2)
+ value(imvn+iptr)=value(locv+3)
+ locy=lynl+nodplc(loc+7)
+ value(locy)=value(locy)+1.0d0
+ locy=lynl+nodplc(loc+8)
+ value(locy)=value(locy)-1.0d0
+ locy=lynl+nodplc(loc+9)
+ value(locy)=value(locy)+1.0d0
+ locy=lynl+nodplc(loc+10)
+ value(locy)=value(locy)-1.0d0
+ loc=nodplc(loc)
+ go to 150
+c
+c current sources
+c
+ 160 loc=locate(10)
+ 170 if (loc.eq.0) go to 200
+ locv=nodplc(loc+1)
+ node1=nodplc(loc+2)
+ node2=nodplc(loc+3)
+ value(lvn+node1)=value(lvn+node1)-value(locv+2)
+ value(imvn+node1)=value(imvn+node1)-value(locv+3)
+ value(lvn+node2)=value(lvn+node2)+value(locv+2)
+ value(imvn+node2)=value(imvn+node2)+value(locv+3)
+ loc=nodplc(loc)
+ go to 170
+c
+c diodes
+c
+ 200 loc=locate(11)
+ 210 if (loc.eq.0) go to 250
+ locv=nodplc(loc+1)
+ area=value(locv+1)
+ locm=nodplc(loc+5)
+ locm=nodplc(locm+1)
+ loct=lx0+nodplc(loc+11)
+ gspr=value(locm+2)*area
+ geq=value(loct+2)
+ xceq=value(loct+4)*omega
+ locy=lynl+nodplc(loc+13)
+ value(locy)=value(locy)+gspr
+ locy=lynl+nodplc(loc+14)
+ locyi=imynl+nodplc(loc+14)
+ value(locy)=value(locy)+geq
+ value(locyi)=value(locyi)+xceq
+ locy=lynl+nodplc(loc+15)
+ locyi=imynl+nodplc(loc+15)
+ value(locy)=value(locy)+geq+gspr
+ value(locyi)=value(locyi)+xceq
+ locy=lynl+nodplc(loc+7)
+ value(locy)=value(locy)-gspr
+ locy=lynl+nodplc(loc+8)
+ locyi=imynl+nodplc(loc+8)
+ value(locy)=value(locy)-geq
+ value(locyi)=value(locyi)-xceq
+ locy=lynl+nodplc(loc+9)
+ value(locy)=value(locy)-gspr
+ locy=lynl+nodplc(loc+10)
+ locyi=imynl+nodplc(loc+10)
+ value(locy)=value(locy)-geq
+ value(locyi)=value(locyi)-xceq
+ loc=nodplc(loc)
+ go to 210
+c
+c bjts
+c
+ 250 loc=locate(12)
+ 260 if (loc.eq.0) go to 300
+ locv=nodplc(loc+1)
+ area=value(locv+1)
+ locm=nodplc(loc+8)
+ locm=nodplc(locm+1)
+ loct=lx0+nodplc(loc+22)
+ gcpr=value(locm+20)*area
+ gepr=value(locm+19)*area
+ gpi=value(loct+4)
+ gmu=value(loct+5)
+ gm=value(loct+6)
+ go=value(loct+7)
+ xgm=0.0d0
+ td=value(locm+28)
+ if(td.eq.0.0d0) go to 270
+ arg=td*omega
+ gm=gm+go
+ xgm=-gm*dsin(arg)
+ gm=gm*dcos(arg)-go
+ 270 gx=value(loct+16)
+ xcpi=value(loct+9)*omega
+ xcmu=value(loct+11)*omega
+ xcbx=value(loct+15)*omega
+ xccs=value(loct+13)*omega
+ xcmcb=value(loct+17)*omega
+ locy=lynl+nodplc(loc+24)
+ value(locy)=value(locy)+gcpr
+ locy=lynl+nodplc(loc+25)
+ locyi=imynl+nodplc(loc+25)
+ value(locy)=value(locy)+gx
+ value(locyi)=value(locyi)+xcbx
+ locy=lynl+nodplc(loc+26)
+ value(locy)=value(locy)+gepr
+ locy=lynl+nodplc(loc+27)
+ locyi=imynl+nodplc(loc+27)
+ value(locy)=value(locy)+gmu+go+gcpr
+ value(locyi)=value(locyi)+xcmu+xccs+xcbx
+ locy=lynl+nodplc(loc+28)
+ locyi=imynl+nodplc(loc+28)
+ value(locy)=value(locy)+gx+gpi+gmu
+ value(locyi)=value(locyi)+xcpi+xcmu+xcmcb
+ locy=lynl+nodplc(loc+29)
+ locyi=imynl+nodplc(loc+29)
+ value(locy)=value(locy)+gpi+gepr+gm+go
+ value(locyi)=value(locyi)+xcpi+xgm
+ locy=lynl+nodplc(loc+10)
+ value(locy)=value(locy)-gcpr
+ locy=lynl+nodplc(loc+11)
+ value(locy)=value(locy)-gx
+ locy=lynl+nodplc(loc+12)
+ value(locy)=value(locy)-gepr
+ locy=lynl+nodplc(loc+13)
+ value(locy)=value(locy)-gcpr
+ locy=lynl+nodplc(loc+14)
+ locyi=imynl+nodplc(loc+14)
+ value(locy)=value(locy)-gmu+gm
+ value(locyi)=value(locyi)-xcmu+xgm
+ locy=lynl+nodplc(loc+15)
+ locyi=imynl+nodplc(loc+15)
+ value(locy)=value(locy)-gm-go
+ value(locyi)=value(locyi)-xgm
+ locy=lynl+nodplc(loc+16)
+ value(locy)=value(locy)-gx
+ locy=lynl+nodplc(loc+17)
+ locyi=imynl+nodplc(loc+17)
+ value(locy)=value(locy)-gmu
+ value(locyi)=value(locyi)-xcmu-xcmcb
+ locy=lynl+nodplc(loc+18)
+ locyi=imynl+nodplc(loc+18)
+ value(locy)=value(locy)-gpi
+ value(locyi)=value(locyi)-xcpi
+ locy=lynl+nodplc(loc+19)
+ value(locy)=value(locy)-gepr
+ locy=lynl+nodplc(loc+20)
+ locyi=imynl+nodplc(loc+20)
+ value(locy)=value(locy)-go
+ value(locyi)=value(locyi)+xcmcb
+ locy=lynl+nodplc(loc+21)
+ locyi=imynl+nodplc(loc+21)
+ value(locy)=value(locy)-gpi-gm
+ value(locyi)=value(locyi)-xcpi-xgm-xcmcb
+ locyi=imynl+nodplc(loc+31)
+ value(locyi)=value(locyi)+xccs
+ locyi=imynl+nodplc(loc+32)
+ value(locyi)=value(locyi)-xccs
+ locyi=imynl+nodplc(loc+33)
+ value(locyi)=value(locyi)-xccs
+ locyi=imynl+nodplc(loc+34)
+ value(locyi)=value(locyi)-xcbx
+ locyi=imynl+nodplc(loc+35)
+ value(locyi)=value(locyi)-xcbx
+ loc=nodplc(loc)
+ go to 260
+c
+c jfets
+c
+ 300 loc=locate(13)
+ 310 if (loc.eq.0) go to 350
+ locv=nodplc(loc+1)
+ area=value(locv+1)
+ locm=nodplc(loc+7)
+ locm=nodplc(locm+1)
+ loct=lx0+nodplc(loc+19)
+ gdpr=value(locm+4)*area
+ gspr=value(locm+5)*area
+ gm=value(loct+5)
+ gds=value(loct+6)
+ ggs=value(loct+7)
+ xgs=value(loct+9)*omega
+ ggd=value(loct+8)
+ xgd=value(loct+11)*omega
+ locy=lynl+nodplc(loc+20)
+ value(locy)=value(locy)+gdpr
+ locy=lynl+nodplc(loc+21)
+ locyi=imynl+nodplc(loc+21)
+ value(locy)=value(locy)+ggd+ggs
+ value(locyi)=value(locyi)+xgd+xgs
+ locy=lynl+nodplc(loc+22)
+ value(locy)=value(locy)+gspr
+ locy=lynl+nodplc(loc+23)
+ locyi=imynl+nodplc(loc+23)
+ value(locy)=value(locy)+gdpr+gds+ggd
+ value(locyi)=value(locyi)+xgd
+ locy=lynl+nodplc(loc+24)
+ locyi=imynl+nodplc(loc+24)
+ value(locy)=value(locy)+gspr+gds+gm+ggs
+ value(locyi)=value(locyi)+xgs
+ locy=lynl+nodplc(loc+9)
+ value(locy)=value(locy)-gdpr
+ locy=lynl+nodplc(loc+10)
+ locyi=imynl+nodplc(loc+10)
+ value(locy)=value(locy)-ggd
+ value(locyi)=value(locyi)-xgd
+ locy=lynl+nodplc(loc+11)
+ locyi=imynl+nodplc(loc+11)
+ value(locy)=value(locy)-ggs
+ value(locyi)=value(locyi)-xgs
+ locy=lynl+nodplc(loc+12)
+ value(locy)=value(locy)-gspr
+ locy=lynl+nodplc(loc+13)
+ value(locy)=value(locy)-gdpr
+ locy=lynl+nodplc(loc+14)
+ locyi=imynl+nodplc(loc+14)
+ value(locy)=value(locy)-ggd+gm
+ value(locyi)=value(locyi)-xgd
+ locy=lynl+nodplc(loc+15)
+ value(locy)=value(locy)-gds-gm
+ locy=lynl+nodplc(loc+16)
+ locyi=imynl+nodplc(loc+16)
+ value(locy)=value(locy)-ggs-gm
+ value(locyi)=value(locyi)-xgs
+ locy=lynl+nodplc(loc+17)
+ value(locy)=value(locy)-gspr
+ locy=lynl+nodplc(loc+18)
+ value(locy)=value(locy)-gds
+ loc=nodplc(loc)
+ go to 310
+c
+c mosfets
+c
+ 350 loc=locate(14)
+ 360 if (loc.eq.0) go to 400
+ locv=nodplc(loc+1)
+ locm=nodplc(loc+8)
+ itype=nodplc(locm+2)
+ locm=nodplc(locm+1)
+ loct=lx0+nodplc(loc+26)
+ gdpr=value(locv+11)
+ gspr=value(locv+12)
+ if(itype.eq.0) go to 380
+ xl=value(locv+1)-2.0d0*value(locm+20)
+ xw=value(locv+2)-2.0d0*value(locm+36)
+ covlgs=value(locm+8)*xw
+ covlgd=value(locm+9)*xw
+ covlgb=value(locm+10)*xl
+ didvg=value(loct)
+ didvd=value(loct+1)
+ didvs=value(loct+2)
+ geqbd=value(loct+3)
+ geqbs=value(loct+5)
+ ccgg=value(loct+8)
+ ccgd=value(loct+9)
+ ccgs=value(loct+10)
+ ccbg=value(loct+11)
+ ccbd=value(loct+12)
+ ccbs=value(loct+13)
+ capbd=value(loct+14)
+ capbs=value(loct+15)
+ xccg2=-0.5d0*(ccgg+ccbg)*omega
+ xccd2=-0.5d0*(ccgd+ccbd)*omega
+ xccs2=-0.5d0*(ccgs+ccbs)*omega
+ xccdg=xccg2-covlgd*omega
+ xccdd=xccd2+(capbd+covlgd)*omega
+ xccds=xccs2
+ xccsg=xccg2-covlgs*omega
+ xccsd=xccd2
+ xccss=xccs2+(capbs+covlgs)*omega
+ xccgg=(ccgg+covlgd+covlgs+covlgb)*omega
+ xccgd=(ccgd-covlgd)*omega
+ xccgs=(ccgs-covlgs)*omega
+ xccbg=(ccbg-covlgb)*omega
+ xccbd=xccbd+(ccbd-capbd)*omega
+ xccbs=xccbs+(ccbs-capbs)*omega
+ gccdg=didvg
+ gccdd=geqbd+didvd
+ gccds=didvs
+ gccsg=-didvg
+ gccsd=-didvd
+ gccss=geqbs-didvs
+ gccbd=-geqbd
+ gccbs=-geqbs
+ locyi=imynl+nodplc(loc+28)
+ value(locyi)=value(locyi)+xccgg
+ locyi=imynl+nodplc(loc+30)
+ value(locyi)=value(locyi)-xccbg-xccbd-xccbs
+ locyi=imynl+nodplc(loc+31)
+ value(locyi)=value(locyi)+xccdd
+ locyi=imynl+nodplc(loc+32)
+ value(locyi)=value(locyi)+xccss
+ locyi=imynl+nodplc(loc+11)
+ value(locyi)=value(locyi)-xccgg-xccgd-xccgs
+ locyi=imynl+nodplc(loc+12)
+ value(locyi)=value(locyi)+xccgd
+ locyi=imynl+nodplc(loc+13)
+ value(locyi)=value(locyi)+xccgs
+ locyi=imynl+nodplc(loc+15)
+ value(locyi)=value(locyi)+xccbg
+ locyi=imynl+nodplc(loc+16)
+ value(locyi)=value(locyi)+xccbd
+ locyi=imynl+nodplc(loc+17)
+ value(locyi)=value(locyi)+xccbs
+ locyi=imynl+nodplc(loc+19)
+ value(locyi)=value(locyi)+xccdg
+ locyi=imynl+nodplc(loc+20)
+ value(locyi)=value(locyi)-xccdg-xccdd-xccds
+ locyi=imynl+nodplc(loc+21)
+ value(locyi)=value(locyi)+xccds
+ locyi=imynl+nodplc(loc+22)
+ value(locyi)=value(locyi)+xccsg
+ locyi=imynl+nodplc(loc+24)
+ value(locyi)=value(locyi)-xccsg-xccsd-xccss
+ locyi=imynl+nodplc(loc+25)
+ value(locyi)=value(locyi)+xccsd
+ locy=lynl+nodplc(loc+27)
+ value(locy)=value(locy)+gdpr
+ locy=lynl+nodplc(loc+29)
+ value(locy)=value(locy)+gspr
+ locy=lynl+nodplc(loc+30)
+ value(locy)=value(locy)-gccbd-gccbs
+ locy=lynl+nodplc(loc+31)
+ value(locy)=value(locy)+gdpr+gccdd
+ locy=lynl+nodplc(loc+32)
+ value(locy)=value(locy)+gspr+gccss
+ locy=lynl+nodplc(loc+10)
+ value(locy)=value(locy)-gdpr
+ locy=lynl+nodplc(loc+14)
+ value(locy)=value(locy)-gspr
+ locy=lynl+nodplc(loc+16)
+ value(locy)=value(locy)+gccbd
+ locy=lynl+nodplc(loc+17)
+ value(locy)=value(locy)+gccbs
+ locy=lynl+nodplc(loc+18)
+ value(locy)=value(locy)-gdpr
+ locy=lynl+nodplc(loc+19)
+ value(locy)=value(locy)+gccdg
+ locy=lynl+nodplc(loc+20)
+ value(locy)=value(locy)-gccdg-gccdd-gccds
+ locy=lynl+nodplc(loc+21)
+ value(locy)=value(locy)+gccds
+ locy=lynl+nodplc(loc+22)
+ value(locy)=value(locy)+gccsg
+ locy=lynl+nodplc(loc+23)
+ value(locy)=value(locy)-gspr
+ locy=lynl+nodplc(loc+24)
+ value(locy)=value(locy)-gccsg-gccsd-gccss
+ locy=lynl+nodplc(loc+25)
+ value(locy)=value(locy)+gccsd
+ loc=nodplc(loc)
+ go to 360
+c... ga-as fets
+ 380 continue
+ devmod=value(locv+8)
+ ggd=value(loct+6)
+ gm=0.0d0
+ gds=0.0d0
+ gmj1=value(loct+7)
+ gdb=value(loct+8)
+ ggs=value(loct+9)
+ xcds=value(loct+10)*omega
+ gsb=value(loct+11)
+ xcgs=value(loct+12)*omega
+ gmj2=value(loct+13)
+ xcgd=value(loct+14)*omega
+ xcgb=value(loct+16)*omega
+c write(6,1001) ggd,gm,gds,gmj1,gdb,ggs,gsb,gmj2
+ 1001 format(' ggd gm gds gmj1 gdb ggs gsb gmj2'/,1x,1p8e9.1)
+ if(devmod.gt.0.0d0) go to 385
+ gmrev=gm
+ gmnrm=0.0d0
+ gm=0.0d0
+ gmj1r=gmj1
+ gmj1=0.0d0
+ gmj1n=0.0d0
+ gmj2n=0.0d0
+ gmj2r=0.0d0
+ go to 390
+ 385 gmnrm=gm
+ gmrev=0.0d0
+ gm=0.0d0
+ gmj2n=gmj2
+ gmj2r=0.0d0
+ gmj2=0.0d0
+ gmj1r=0.0d0
+ gmj1n=0.0d0
+ 390 locy=lynl+nodplc(loc+27)
+ value(locy)=value(locy)+gdpr
+ locy=lynl+nodplc(loc+28)
+ locyi=imynl+nodplc(loc+28)
+ value(locy)=value(locy)+ggd+ggs
+ value(locyi)=value(locyi)+xcgd+xcgs+xcgb
+ locy=lynl+nodplc(loc+29)
+ value(locy)=value(locy)+gspr
+ locy=lynl+nodplc(loc+30)
+ locyi=imynl+nodplc(loc+30)
+ value(locy)=value(locy)+gdb+gsb+gmj1+gmj2
+ value(locyi)=value(locyi)+xcgb
+ locy=lynl+nodplc(loc+31)
+ locyi=imynl+nodplc(loc+31)
+ value(locy)=value(locy)+gdpr+gds+gdb+ggd-gmrev-gmj1r
+ value(locyi)=value(locyi)+xcds+xcgd
+ locy=lynl+nodplc(loc+32)
+ locyi=imynl+nodplc(loc+32)
+ value(locy)=value(locy)+gspr+gds+gsb+ggs+gmnrm-gmj2n
+ value(locyi)=value(locyi)+xcds+xcgs
+ locy=lynl+nodplc(loc+10)
+ value(locy)=value(locy)-gdpr
+ locyi=imynl+nodplc(loc+11)
+ value(locyi)=value(locyi)-xcgb
+ locy=lynl+nodplc(loc+12)
+ locyi=imynl+nodplc(loc+12)
+ value(locy)=value(locy)-ggd
+ value(locyi)=value(locyi)-xcgd
+ locy=lynl+nodplc(loc+13)
+ locyi=imynl+nodplc(loc+13)
+ value(locy)=value(locy)-ggs
+ value(locyi)=value(locyi)-xcgs
+ locy=lynl+nodplc(loc+14)
+ value(locy)=value(locy)-gspr
+ locy=lynl+nodplc(loc+15)
+ locyi=imynl+nodplc(loc+15)
+ value(locy)=value(locy)-gmj2n-gmj1n-gmj1r-gmj2r-gmj1-gmj2
+ value(locyi)=value(locyi)-xcgb
+ locy=lynl+nodplc(loc+16)
+ value(locy)=value(locy)-gdb+gmj2r+gmj1r
+ locy=lynl+nodplc(loc+17)
+ value(locy)=value(locy)-gsb+gmj1n+gmj2n
+ locy=lynl+nodplc(loc+18)
+ value(locy)=value(locy)-gdpr
+ locy=lynl+nodplc(loc+19)
+ locyi=imynl+nodplc(loc+19)
+ value(locy)=value(locy)-ggd+gmnrm+gmrev+gmj1r+gmj1n+gmj1+gm
+ value(locyi)=value(locyi)-xcgd
+ locy=lynl+nodplc(loc+20)
+ value(locy)=value(locy)-gdb-gmj1-gm
+ locy=lynl+nodplc(loc+21)
+ locyi=imynl+nodplc(loc+21)
+ value(locy)=value(locy)-gds-gmnrm-gmj1n
+ value(locyi)=value(locyi)-xcds
+ locy=lynl+nodplc(loc+22)
+ locyi=imynl+nodplc(loc+22)
+ value(locy)=value(locy)-ggs-gmnrm-gmrev+gmj2r+gmj2n+gmj2-gm
+ value(locyi)=value(locyi)-xcgs
+ locy=lynl+nodplc(loc+23)
+ value(locy)=value(locy)-gspr
+ locy=lynl+nodplc(loc+24)
+ value(locy)=value(locy)-gsb-gmj2+gm
+ locy=lynl+nodplc(loc+25)
+ locyi=imynl+nodplc(loc+25)
+ value(locy)=value(locy)-gds+gmrev-gmj2r
+ value(locyi)=value(locyi)-xcds
+ loc=nodplc(loc)
+ go to 360
+c
+c transmission lines
+c
+ 400 loc=locate(17)
+ 410 if (loc.eq.0) go to 1000
+ locv=nodplc(loc+1)
+ z0=value(locv+1)
+ y0=1.0d0/z0
+ td=value(locv+2)
+ arg=-omega*td
+ rval=dcos(arg)
+ xval=dsin(arg)
+ locy=lynl+nodplc(loc+10)
+ value(locy)=value(locy)+y0
+ locy=lynl+nodplc(loc+11)
+ locyi=imynl+nodplc(loc+11)
+ value(locy)=-y0
+ value(locyi)=0.0d0
+ locy=lynl+nodplc(loc+12)
+ locyi=imynl+nodplc(loc+12)
+ value(locy)=-1.0d0
+ value(locyi)=0.0d0
+ locy=lynl+nodplc(loc+13)
+ value(locy)=value(locy)+y0
+ locy=lynl+nodplc(loc+14)
+ locyi=imynl+nodplc(loc+14)
+ value(locy)=-1.0d0
+ value(locyi)=0.0d0
+ locy=lynl+nodplc(loc+15)
+ locyi=imynl+nodplc(loc+15)
+ value(locy)=-y0
+ value(locyi)=0.0d0
+ locy=lynl+nodplc(loc+16)
+ locyi=imynl+nodplc(loc+16)
+ value(locy)=+y0
+ value(locyi)=0.0d0
+ locy=lynl+nodplc(loc+17)
+ locyi=imynl+nodplc(loc+17)
+ value(locy)=+1.0d0
+ value(locyi)=0.0d0
+ locy=lynl+nodplc(loc+18)
+ locyi=imynl+nodplc(loc+18)
+ value(locy)=+y0
+ value(locyi)=0.0d0
+ locy=lynl+nodplc(loc+19)
+ locyi=imynl+nodplc(loc+19)
+ value(locy)=+1.0d0
+ value(locyi)=0.0d0
+ locy=lynl+nodplc(loc+20)
+ locyi=imynl+nodplc(loc+20)
+ value(locy)=-1.0d0
+ value(locyi)=0.0d0
+ locy=lynl+nodplc(loc+21)
+ locyi=imynl+nodplc(loc+21)
+ value(locy)=-rval
+ value(locyi)=-xval
+ locy=lynl+nodplc(loc+22)
+ locyi=imynl+nodplc(loc+22)
+ value(locy)=+rval
+ value(locyi)=+xval
+ locy=lynl+nodplc(loc+23)
+ locyi=imynl+nodplc(loc+23)
+ value(locy)=+1.0d0
+ value(locyi)=0.0d0
+ locy=lynl+nodplc(loc+24)
+ locyi=imynl+nodplc(loc+24)
+ value(locy)=-rval*z0
+ value(locyi)=-xval*z0
+ locy=lynl+nodplc(loc+25)
+ locyi=imynl+nodplc(loc+25)
+ value(locy)=-rval
+ value(locyi)=-xval
+ locy=lynl+nodplc(loc+26)
+ locyi=imynl+nodplc(loc+26)
+ value(locy)=+rval
+ value(locyi)=+xval
+ locy=lynl+nodplc(loc+27)
+ locyi=imynl+nodplc(loc+27)
+ value(locy)=-1.0d0
+ value(locyi)=0.0d0
+ locy=lynl+nodplc(loc+28)
+ locyi=imynl+nodplc(loc+28)
+ value(locy)=+1.0d0
+ value(locyi)=0.0d0
+ locy=lynl+nodplc(loc+29)
+ locyi=imynl+nodplc(loc+29)
+ value(locy)=-rval*z0
+ value(locyi)=-xval*z0
+ locy=lynl+nodplc(loc+31)
+ locyi=imynl+nodplc(loc+31)
+ value(locy)=-y0
+ value(locyi)=0.0d0
+ locy=lynl+nodplc(loc+32)
+ locyi=imynl+nodplc(loc+32)
+ value(locy)=-y0
+ value(locyi)=0.0d0
+ loc=nodplc(loc)
+ go to 410
+c
+c reorder right-hand side
+c
+ 1000 do 1010 i=2,nstop
+ j=nodplc(iswap+i)
+ value(ndiag+i)=value(lvn+j)
+ value(ndiag+i+nstop)=value(imvn+j)
+ 1010 continue
+ call copy8(value(ndiag+1),value(lvn+1),nstop)
+ call copy8(value(ndiag+nstop+1),value(imvn+1),nstop)
+c
+c finished
+c
+ return
+ end
+ subroutine noise(loco)
+ implicit double precision (a-h,o-z)
+c
+c this routine computes the noise due to various circuit elements.
+c
+ common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem,
+ 1 isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize,
+ 2 junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr,
+ 3 nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1,
+ 4 lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd,
+ 5 imynl,imvn,lcvn,loutpt,nsnod,nsmat,nsval,icnod,icmat,icval
+ common /cirdat/ locate(50),jelcnt(50),nunods,ncnods,numnod,nstop,
+ 1 nut,nlt,nxtrm,ndist,ntlin,ibr,numvs
+ common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet,
+ 1 xmu,mode,modedc,icalc,initf,method,iord,maxord,noncon,iterno,
+ 2 itemno,nosolv,ipostp,iscrch
+ common /miscel/ atime,aprog(3),adate,atitle(10),defl,defw,defad,
+ 1 defas,rstats(50),iwidth,lwidth,nopage
+ common /knstnt/ twopi,xlog2,xlog10,root2,rad,boltz,charge,ctok,
+ 1 gmin,reltol,abstol,vntol,trtol,chgtol,eps0,epssil,epsox
+ common /ac/ fstart,fstop,fincr,skw2,refprl,spw2,jacflg,idfreq,
+ 1 inoise,nosprt,nosout,nosin,idist,idprt
+ common /blank/ value(1000)
+ integer nodplc(64)
+ complex*16 cvalue(32)
+ equivalence (value(1),nodplc(1),cvalue(1))
+c
+c
+ dimension vno1(12),vno2(12),vno3(12),vno4(12),vno5(12),vno6(12)
+ dimension vntot(12),anam(12),string(5)
+ dimension titln(4),v(2)
+ dimension afmt1(3),afmt2(3)
+ complex*16 cval,c(1)
+ equivalence (c(1),v(1),cval)
+ equivalence (v(1),vreal),(v(2),vimag)
+ data titln / 8hnoise an, 8halysis , 8h , 8h /
+ data alsrb,alsrc,alsre,alsrs,alsrd / 2hrb,2hrc,2hre,2hrs,2hrd /
+ data alsib,alsic,alsid,alsfn / 2hib,2hic,2hid,2hfn /
+ data alstot / 5htotal /
+ data aslash,ablnk / 1h/, 1h /
+ data afmt1 /8h(////,11,8hx, (2x,,8ha8)) /
+ data afmt2 /8h(1h0,a8,,8h1p d10.,8h3) /
+ kntr=12
+ if(lwidth.le.80) kntr=7
+ ipos=11
+ call move(afmt1,ipos,ablnk,1,2)
+ call alfnum(kntr,afmt1,ipos)
+ ipos=11
+ call move(afmt2,ipos,ablnk,1,2)
+ call alfnum(kntr,afmt2,ipos)
+ nprnt=0
+ freq=omega/twopi
+ if (icalc.ge.2) go to 10
+ fourkt=4.0d0*charge*vt
+ twoq=2.0d0*charge
+ noposo=nodplc(nosout+2)
+ nonego=nodplc(nosout+3)
+ kntlim=lwidth/11
+ nkntr=1
+ 10 if (nosprt.eq.0) go to 30
+ if (nkntr.gt.icalc) go to 30
+ nprnt=1
+ nkntr=nkntr+nosprt
+ call title(0,lwidth,1,titln)
+ write (6,16) freq
+ 16 format('0 frequency = ',1pd10.3,' hz'/)
+c
+c obtain adjoint circuit solution
+c
+ 30 vnrms=0.0d0
+ cval=cvalue(lcvn+noposo)-cvalue(lcvn+nonego)
+ vout=dsqrt(vreal*vreal+vimag*vimag)
+ vout=dmax1(vout,1.0d-20)
+ call zero8(value(lvn+1),nstop)
+ call zero8(value(imvn+1),nstop)
+ value(lvn+noposo)=-1.0d0
+ value(lvn+nonego)=+1.0d0
+ call acasol
+c
+c resistors
+c
+ if (jelcnt(1).eq.0) go to 200
+ ititle=0
+ 91 format(//'0**** resistor squared noise voltages (sq v/hz)')
+ 100 loc=locate(1)
+ kntr=0
+ 110 if (loc.eq.0) go to 130
+ kntr=kntr+1
+ locv=nodplc(loc+1)
+ anam(kntr)=value(locv)
+ node1=nodplc(loc+2)
+ node2=nodplc(loc+3)
+ cval=cvalue(lcvn+node1)-cvalue(lcvn+node2)
+ vntot(kntr)=(vreal*vreal+vimag*vimag)*fourkt*value(locv+1)
+ vnrms=vnrms+vntot(kntr)
+ if (kntr.ge.kntlim) go to 140
+ 120 loc=nodplc(loc)
+ go to 110
+ 130 if (kntr.eq.0) go to 200
+ 140 if (nprnt.eq.0) go to 160
+ if (ititle.eq.0) write (6,91)
+ ititle=1
+ write (6,afmt1) (anam(i),i=1,kntr)
+ write (6,afmt2) alstot,(vntot(i),i=1,kntr)
+ 160 kntr=0
+ if (loc.ne.0) go to 120
+c
+c diodes
+c
+ 200 if (jelcnt(11).eq.0) go to 300
+ ititle=0
+ 201 format(//'0**** diode squared noise voltages (sq v/hz)')
+ 210 loc=locate(11)
+ kntr=0
+ 220 if (loc.eq.0) go to 240
+ kntr=kntr+1
+ locv=nodplc(loc+1)
+ anam(kntr)=value(locv)
+ node1=nodplc(loc+2)
+ node2=nodplc(loc+3)
+ node3=nodplc(loc+4)
+ locm=nodplc(loc+5)
+ locm=nodplc(locm+1)
+ loct=nodplc(loc+11)
+ area=value(locv+1)
+ fnk=value(locm+10)
+ fna=value(locm+11)
+c
+c ohmic resistance
+c
+ cval=cvalue(lcvn+node1)-cvalue(lcvn+node3)
+ vno1(kntr)=(vreal*vreal+vimag*vimag)*fourkt*value(locm+2)*area
+c
+c junction shot noise and flicker noise
+c
+ cval=cvalue(lcvn+node3)-cvalue(lcvn+node2)
+ vtemp=vreal*vreal+vimag*vimag
+ arg=dmax1(dabs(value(lx0+loct+1)),1.0d-20)
+ vno2(kntr)=vtemp*twoq*arg
+ vno3(kntr)=vtemp*fnk*dexp(fna*dlog(arg))/freq
+ vntot(kntr)=vno1(kntr)+vno2(kntr)+vno3(kntr)
+ vnrms=vnrms+vntot(kntr)
+ if (kntr.ge.kntlim) go to 250
+ 230 loc=nodplc(loc)
+ go to 220
+ 240 if (kntr.eq.0) go to 300
+ 250 if (nprnt.eq.0) go to 260
+ if (ititle.eq.0) write (6,201)
+ ititle=1
+ write (6,afmt1) (anam(i),i=1,kntr)
+ write (6,afmt2) alsrs,(vno1(i),i=1,kntr)
+ write (6,afmt2) alsid,(vno2(i),i=1,kntr)
+ write (6,afmt2) alsfn,(vno3(i),i=1,kntr)
+ write (6,afmt2) alstot,(vntot(i),i=1,kntr)
+ 260 kntr=0
+ if (loc.ne.0) go to 230
+c
+c bipolar junction transistors
+c
+ 300 if (jelcnt(12).eq.0) go to 400
+ ititle=0
+ 301 format(//'0**** transistor squared noise voltages (sq v/hz)')
+ 310 loc=locate(12)
+ kntr=0
+ 320 if (loc.eq.0) go to 340
+ kntr=kntr+1
+ locv=nodplc(loc+1)
+ anam(kntr)=value(locv)
+ node1=nodplc(loc+2)
+ node2=nodplc(loc+3)
+ node3=nodplc(loc+4)
+ node4=nodplc(loc+5)
+ node5=nodplc(loc+6)
+ node6=nodplc(loc+7)
+ locm=nodplc(loc+8)
+ locm=nodplc(locm+1)
+ loct=nodplc(loc+22)
+ area=value(locv+1)
+ fnk=value(locm+44)
+ fna=value(locm+45)
+c
+c extrinsic resistances
+c
+c... base resistance
+ cval=cvalue(lcvn+node2)-cvalue(lcvn+node5)
+ vno1(kntr)=(vreal*vreal+vimag*vimag)*fourkt*value(lx0+loct+16)
+ 1 *area
+c... collector resistance
+ cval=cvalue(lcvn+node1)-cvalue(lcvn+node4)
+ vno2(kntr)=(vreal*vreal+vimag*vimag)*fourkt*value(locm+20)*area
+c... emitter resistance
+ cval=cvalue(lcvn+node3)-cvalue(lcvn+node6)
+ vno3(kntr)=(vreal*vreal+vimag*vimag)*fourkt*value(locm+19)*area
+c
+c base current shot noise and flicker noise
+c
+ cval=cvalue(lcvn+node5)-cvalue(lcvn+node6)
+ vtemp=vreal*vreal+vimag*vimag
+ arg=dmax1(dabs(value(lx0+loct+3)),1.0d-20)
+ vno4(kntr)=vtemp*twoq*arg
+ vno5(kntr)=vtemp*fnk*dexp(fna*dlog(arg))/freq
+c
+c collector current shot noise
+c
+ cval=cvalue(lcvn+node4)-cvalue(lcvn+node6)
+ vno6(kntr)=(vreal*vreal+vimag*vimag)*twoq*dabs(value(lx0+loct+2))
+ vntot(kntr)=vno1(kntr)+vno2(kntr)+vno3(kntr)+vno4(kntr)+vno5(kntr)
+ 1 +vno6(kntr)
+ vnrms=vnrms+vntot(kntr)
+ if (kntr.ge.kntlim) go to 350
+ 330 loc=nodplc(loc)
+ go to 320
+ 340 if (kntr.eq.0) go to 400
+ 350 if (nprnt.eq.0) go to 360
+ if (ititle.eq.0) write (6,301)
+ ititle=1
+ write (6,afmt1) (anam(i),i=1,kntr)
+ write (6,afmt2) alsrb,(vno1(i),i=1,kntr)
+ write (6,afmt2) alsrc,(vno2(i),i=1,kntr)
+ write (6,afmt2) alsre,(vno3(i),i=1,kntr)
+ write (6,afmt2) alsib,(vno4(i),i=1,kntr)
+ write (6,afmt2) alsic,(vno6(i),i=1,kntr)
+ write (6,afmt2) alsfn,(vno5(i),i=1,kntr)
+ write (6,afmt2) alstot,(vntot(i),i=1,kntr)
+ 360 kntr=0
+ if (loc.ne.0) go to 330
+c
+c jfets
+c
+ 400 if (jelcnt(13).eq.0) go to 500
+ ititle=0
+ 401 format(//'0**** jfet squared noise voltages (sq v/hz)')
+ 410 loc=locate(13)
+ kntr=0
+ 420 if (loc.eq.0) go to 440
+ kntr=kntr+1
+ locv=nodplc(loc+1)
+ anam(kntr)=value(locv)
+ node1=nodplc(loc+2)
+ node2=nodplc(loc+3)
+ node3=nodplc(loc+4)
+ node4=nodplc(loc+5)
+ node5=nodplc(loc+6)
+ locm=nodplc(loc+7)
+ locm=nodplc(locm+1)
+ loct=nodplc(loc+19)
+ area=value(locv+1)
+ fnk=value(locm+10)
+ fna=value(locm+11)
+c
+c extrinsic resistances
+c
+c... drain resistance
+ cval=cvalue(lcvn+node1)-cvalue(lcvn+node4)
+ vno1(kntr)=(vreal*vreal+vimag*vimag)*fourkt*value(locm+4)*area
+c... source resistance
+ cval=cvalue(lcvn+node3)-cvalue(lcvn+node5)
+ vno2(kntr)=(vreal*vreal+vimag*vimag)*fourkt*value(locm+5)*area
+c
+c drain current shot noise and flicker noise
+c
+ cval=cvalue(lcvn+node4)-cvalue(lcvn+node5)
+ vtemp=vreal*vreal+vimag*vimag
+ vno3(kntr)=vtemp*fourkt*2.0d0*dabs(value(lx0+loct+5))/3.0d0
+ arg=dmax1(dabs(value(lx0+loct+3)),1.0d-20)
+ vno4(kntr)=vtemp*fnk*dexp(fna*dlog(arg))/freq
+ vntot(kntr)=vno1(kntr)+vno2(kntr)+vno3(kntr)+vno4(kntr)
+ vnrms=vnrms+vntot(kntr)
+ if (kntr.ge.kntlim) go to 450
+ 430 loc=nodplc(loc)
+ go to 420
+ 440 if (kntr.eq.0) go to 500
+ 450 if (nprnt.eq.0) go to 460
+ if (ititle.eq.0) write (6,401)
+ ititle=1
+ write (6,afmt1) (anam(i),i=1,kntr)
+ write (6,afmt2) alsrd,(vno1(i),i=1,kntr)
+ write (6,afmt2) alsrs,(vno2(i),i=1,kntr)
+ write (6,afmt2) alsid,(vno3(i),i=1,kntr)
+ write (6,afmt2) alsfn,(vno4(i),i=1,kntr)
+ write (6,afmt2) alstot,(vntot(i),i=1,kntr)
+ 460 kntr=0
+ if (loc.ne.0) go to 430
+c
+c mosfets
+c
+ 500 if (jelcnt(14).eq.0) go to 600
+ ititle=0
+ 501 format(//'0**** mosfet squared noise voltages (sq v/hz)')
+ 510 loc=locate(14)
+ kntr=0
+ 520 if (loc.eq.0) go to 540
+ kntr=kntr+1
+ locv=nodplc(loc+1)
+ anam(kntr)=value(locv)
+ node1=nodplc(loc+2)
+ node2=nodplc(loc+3)
+ node3=nodplc(loc+4)
+ node4=nodplc(loc+5)
+ node5=nodplc(loc+6)
+ node6=nodplc(loc+7)
+ locm=nodplc(loc+8)
+ itype=nodplc(locm+2)
+ loct=nodplc(loc+26)
+ locm=nodplc(locm+1)
+ if(itype.eq.0) go to 522
+ xl=value(locv+1)-2.0d0*value(locm+20)
+ xw=value(locm+2)-2.0d0*value(locm+36)
+ cox=value(locm+13)*xl*xw
+ fnk=value(locm+27)
+ fna=value(locm+28)
+c
+c extrinsic resistances
+c
+c... drain resistance
+ 522 cval=cvalue(lcvn+node1)-cvalue(lcvn+node5)
+ vno1(kntr)=(vreal*vreal+vimag*vimag)*fourkt*value(locv+11)
+c... source resistance
+ cval=cvalue(lcvn+node3)-cvalue(lcvn+node6)
+ vno2(kntr)=(vreal*vreal+vimag*vimag)*fourkt*value(locv+12)
+c
+c drain current shot noise and flicker noise
+c
+ cval=cvalue(lcvn+node5)-cvalue(lcvn+node6)
+ vtemp=vreal*vreal+vimag*vimag
+ gm=value(lx0+loct+7)
+ arg=dmax1(dabs(value(lx0+loct+4)),1.0d-20)
+ if(itype.ne.0) go to 524
+ modeop=value(locv+8)
+ if(modeop.le.0) gm=value(lx0+loct+13)
+ if(value(locm+10).ne.0.0d0) gm=value(locm+10)
+ xnexp=value(locm+11)
+ fnk=value(locm+12)
+ fna=value(locm+13)
+ vno4(kntr)=vtemp*fnk*dexp(fna*dlog(arg))/(freq**xnexp)
+ 524 vno3(kntr)=vtemp*fourkt*dabs(gm)/1.5d0
+ if(itype.eq.0) go to 525
+ vno4(kntr)=vtemp*fnk*dexp(fna*dlog(arg))/(freq*cox)
+ 525 vntot(kntr)=vno1(kntr)+vno2(kntr)+vno3(kntr)+vno4(kntr)
+ vnrms=vnrms+vntot(kntr)
+ if (kntr.ge.kntlim) go to 550
+ 530 loc=nodplc(loc)
+ go to 520
+ 540 if (kntr.eq.0) go to 600
+ 550 if (nprnt.eq.0) go to 560
+ if (ititle.eq.0) write (6,501)
+ ititle=1
+ write (6,afmt1) (anam(i),i=1,kntr)
+ write (6,afmt2) alsrd,(vno1(i),i=1,kntr)
+ write (6,afmt2) alsrs,(vno2(i),i=1,kntr)
+ write (6,afmt2) alsid,(vno3(i),i=1,kntr)
+ write (6,afmt2) alsfn,(vno4(i),i=1,kntr)
+ write (6,afmt2) alstot,(vntot(i),i=1,kntr)
+ 560 kntr=0
+ if (loc.ne.0) go to 530
+c
+c compute equivalent input noise voltage
+c
+ 600 vnout=dsqrt(vnrms)
+ vnin=vnout/vout
+ if (nprnt.eq.0) go to 620
+ do 610 i=1,5
+ string(i)=ablnk
+ 610 continue
+ ioutyp=1
+ ipos=1
+ call outnam(nosout,ioutyp,string,ipos)
+ call move(string,ipos,aslash,1,1)
+ ipos=ipos+1
+ locv=nodplc(nosin+1)
+ anam1=value(locv)
+ call move(string,ipos,anam1,1,8)
+ write (6,611) vnrms,vnout,string,vout,anam1,vnin
+ 611 format(////,
+ 1 '0**** total output noise voltage',9x,'= ',1pd10.3,' sq v/hz'/,
+ 2 1h0,40x,'= ',d10.3,' v/rt hz'/,
+ 3 '0 transfer function value:',
+ 4 1h0,7x,4a8,a1,'= ',d10.3,/,
+ 5 '0 equivalent input noise at ',a8,' = ',d10.3,' /rt hz')
+c
+c save noise outputs
+c
+ 620 loc=locate(44)
+ 630 if (loc.eq.0) go to 1000
+ iseq=nodplc(loc+4)
+ if (nodplc(loc+5).ne.2) go to 640
+ cvalue(loco+iseq)=vnout
+ go to 650
+ 640 cvalue(loco+iseq)=vnin
+ 650 loc=nodplc(loc)
+ go to 630
+c
+c finished
+c
+ 1000 return
+ end
+ subroutine acasol
+ implicit double precision (a-h,o-z)
+c
+c this routine evaluates the response of the adjoint circuit by
+c doing a forward/backward substitution step using the transpose of the
+c circuit equation coefficient matrix.
+c
+ common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem,
+ 1 isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize,
+ 2 junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr,
+ 3 nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1,
+ 4 lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd,
+ 5 imynl,imvn,lcvn,loutpt,nsnod,nsmat,nsval,icnod,icmat,icval
+ common /cirdat/ locate(50),jelcnt(50),nunods,ncnods,numnod,nstop,
+ 1 nut,nlt,nxtrm,ndist,ntlin,ibr,numvs
+ common /blank/ value(1000)
+ integer nodplc(64)
+ complex*16 cvalue(32)
+ equivalence (value(1),nodplc(1),cvalue(1))
+c
+c evaluates adjoint response by doing forward/backward substitution on
+c the transpose of the y matrix
+c
+c
+c forward substitution
+c
+ do 20 i=2,nstop
+ io=nodplc(iorder+i)
+ call cdiv(value(lvn+io),value(imvn+io),value(lynl+io),
+ 1 value(imynl+io),value(lvn+io),value(imvn+io))
+ jstart=nodplc(iur+i)
+ jstop=nodplc(iur+i+1)-1
+ if (jstart.gt.jstop) go to 20
+ if (value(lvn+io).ne.0.0d0) go to 5
+ if (value(imvn+io).eq.0.0d0) go to 20
+ 5 do 10 j=jstart,jstop
+ jo=nodplc(iuc+j)
+ jo=nodplc(iorder+jo)
+ call cmult(value(lyu+j),value(imynl+nstop+j),value(lvn+io),
+ 1 value(imvn+io),xreal,ximag)
+ value(lvn+jo)=value(lvn+jo)-xreal
+ value(imvn+jo)=value(imvn+jo)-ximag
+ 10 continue
+ 20 continue
+c
+c backward substitution
+c
+ k=nstop+1
+ do 40 i=2,nstop
+ k=k-1
+ io=nodplc(iorder+k)
+ jstart=nodplc(ilc+k)
+ jstop=nodplc(ilc+k+1)-1
+ if (jstart.gt.jstop) go to 40
+ do 30 j=jstart,jstop
+ jo=nodplc(ilr+j)
+ jo=nodplc(iorder+jo)
+ call cmult(value(lyl+j),value(imynl+nstop+nut+j),value(lvn+jo),
+ 1 value(imvn+jo),xreal,ximag)
+ value(lvn+io)=value(lvn+io)-xreal
+ value(imvn+io)=value(imvn+io)-ximag
+ 30 continue
+ 40 continue
+c
+c reorder right-hand side
+c
+ do 50 i=2,nstop
+ j=nodplc(iswap+i)
+ value(ndiag+i)=value(lvn+j)
+ value(ndiag+i+nstop)=value(imvn+j)
+ 50 continue
+ call copy8(value(ndiag+1),value(lvn+1),nstop)
+ call copy8(value(ndiag+nstop+1),value(imvn+1),nstop)
+ do 120 i=2,nstop
+ cvalue(lcvn+i)=dcmplx(value(lvn+i),value(imvn+i))
+ 120 continue
+ cvalue(lcvn+1)=dcmplx(0.0d0,0.0d0)
+c
+c finished
+c
+ return
+ end
+ subroutine dinit
+ implicit double precision (a-h,o-z)
+c
+c this routine performs storage-allocation and one-time computation
+c needed to do the small-signal distortion analysis.
+c
+ common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem,
+ 1 isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize,
+ 2 junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr,
+ 3 nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1,
+ 4 lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd,
+ 5 imynl,imvn,lcvn,loutpt,nsnod,nsmat,nsval,icnod,icmat,icval
+ common /cirdat/ locate(50),jelcnt(50),nunods,ncnods,numnod,nstop,
+ 1 nut,nlt,nxtrm,ndist,ntlin,ibr,numvs
+ common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet,
+ 1 xmu,mode,modedc,icalc,initf,method,iord,maxord,noncon,iterno,
+ 2 itemno,nosolv,ipostp,iscrch
+ common /flags/ iprnta,iprntl,iprntm,iprntn,iprnto,limtim,limpts,
+ 1 lvlcod,lvltim,itl1,itl2,itl3,itl4,itl5,igoof,nogo,keof
+ common /blank/ value(1000)
+ integer nodplc(64)
+ complex*16 cvalue(32)
+ equivalence (value(1),nodplc(1),cvalue(1))
+c
+c
+ call getm8(ld0,ndist)
+ call getm16(ld1,5*nstop)
+c
+c bipolar junction transistors
+c
+ loc=locate(12)
+ 100 if (loc.eq.0) go to 200
+ locv=nodplc(loc+1)
+ area=value(locv+1)
+ locm=nodplc(loc+8)
+ locm=nodplc(locm+1)
+ loct=lx0+nodplc(loc+22)
+ locd=ld0+nodplc(loc+23)
+ csat=value(locm+1)*area
+ ova=value(locm+4)
+ tf=value(locm+24)
+ tr=value(locm+33)
+ czbe=value(locm+21)*area
+ czbc=value(locm+29)*area
+ pe=value(locm+22)
+ xme=value(locm+23)
+ pc=value(locm+30)
+ xmc=value(locm+31)
+ fcpe=value(locm+46)
+ fcpc=value(locm+50)
+ vbe=value(loct)
+ vbc=value(loct+1)
+ gpi=value(loct+4)
+ go=value(loct+7)
+ gm=value(loct+6)
+ gmu=value(loct+5)
+ if (vbe.gt.0.0d0) go to 110
+ evbe=1.0d0
+ cbe=csat*vbe/vt
+ go to 120
+ 110 evbe=dexp(vbe/vt)
+ cbe=csat*(evbe-1.0d0)
+ 120 if (vbc.gt.0.0d0) go to 130
+ evbc=1.0d0
+ cbc=csat*vbc/vt
+ arg=1.0d0-vbc/pc
+ go to 140
+ 130 evbc=dexp(vbc/vt)
+ cbc=csat*(evbc-1.0d0)
+ 140 if (vbe.ge.fcpe) go to 150
+ arg=1.0d0-vbe/pe
+ sarg=dexp(xme*dlog(arg))
+ cjeo=czbe/sarg
+ argbe=pe-vbe
+ cje1=xme*cjeo/argbe
+ cje2=xme*(1.0d0+xme)*cje1/argbe
+ go to 160
+ 150 denom=dexp((1.0d0+xme)*dlog(1.0d0-fcpe))
+ cjeo=czbe*(1.0d0-fcpe*(1.0d0+xme)+xme*vbe/pe)/denom
+ cje1=czbe*xme/(denom*pe)
+ cje2=0.0d0
+ 160 if (vbc.ge.fcpc) go to 170
+ arg=1.0d0-vbc/pc
+ sarg=dexp(xmc*dlog(arg))
+ cjco=czbc/sarg
+ argbc=pc-vbc
+ cjc1=xmc*cjco/argbc
+ cjc2=xmc*(1.0d0+xmc)*cjc1/argbc
+ go to 180
+ 170 denom=dexp((1.0d0+xmc)*dlog(1.0d0-fcpc))
+ cjco=czbc*(1.0d0-fcpc*(1.0d0+xmc)+xmc*vbc/pc)/denom
+ cjc1=czbc*xmc/(denom*pc)
+ cjc2=0.0d0
+ 180 twovt=vt+vt
+ go2=(-go+csat*(evbe+evbc)*ova)/twovt
+ gmo2=(cbe+csat)*ova/vt-2.0d0*go2
+ gm2=(gm+go)/twovt-gmo2-go2
+ gmu2=gmu/twovt
+ if (vbc.le.0.0d0) gmu2=0.0d0
+ gpi2=gpi/twovt
+ if (vbe.le.0.0d0) gpi2=0.0d0
+ cbo=tf*csat*evbe/vt
+ cbor=tr*csat*evbc/vt
+ cb1=cbo/vt
+ cb1r=cbor/vt
+ trivt=3.0d0*vt
+ go3=-(go2+(cbc+csat)*ova/twovt)/trivt
+ gmo23=-3.0d0*go3
+ gm2o3=-gmo23+(cbe+csat)*ova/(vt*twovt)
+ gm3=(gm2-(cbe-cbc)*ova/twovt)/trivt
+ gmu3=gmu2/trivt
+ gpi3=gpi2/trivt
+ cb2=cb1/twovt
+ cb2r=cb1r/twovt
+ value(locd)=cje1
+ value(locd+1)=cje2
+ value(locd+2)=cjc1
+ value(locd+3)=cjc2
+ value(locd+4)=go2
+ value(locd+5)=gmo2
+ value(locd+6)=gm2
+ value(locd+7)=gmu2
+ value(locd+8)=gpi2
+ value(locd+9)=cbo
+ value(locd+10)=cbor
+ value(locd+11)=cb1
+ value(locd+12)=cb1r
+ value(locd+13)=go3
+ value(locd+14)=gmo23
+ value(locd+15)=gm2o3
+ value(locd+16)=gm3
+ value(locd+17)=gmu3
+ value(locd+18)=gpi3
+ value(locd+19)=cb2
+ value(locd+20)=cb2r
+ loc=nodplc(loc)
+ go to 100
+c
+c diodes
+c
+ 200 loc=locate(11)
+ 210 if (loc.eq.0) go to 300
+ locv=nodplc(loc+1)
+ area=value(locv+1)
+ locm=nodplc(loc+5)
+ locm=nodplc(locm+1)
+ loct=lx0+nodplc(loc+11)
+ locd=ld0+nodplc(loc+12)
+ csat=value(locm+1)*area
+ vte=value(locm+3)*vt
+ tau=value(locm+4)
+ czero=value(locm+5)*area
+ phib=value(locm+6)
+ xm=value(locm+7)
+ fcpb=value(locm+12)
+ vd=value(loct)
+ geq=value(loct+2)
+ evd=1.0d0
+ if (vd.ge.0.0d0) evd=dexp(vd/vte)
+ if (vd.ge.fcpb) go to 220
+ arg=1.0d0-vd/phib
+ sarg=dexp(xm*dlog(arg))
+ cdjo=czero/sarg
+ argd=phib-vd
+ cdj1=xm*czero/argd
+ cdj2=xm*(1.0d0+xm)*cdj1/argd
+ go to 230
+ 220 denom=dexp((1.0d0+xm)*dlog(1.0d0-fcpb))
+ cdjo=czero*(1.0d0-fcpb*(1.0d0+xm)+xm*vd/phib)/denom
+ cdj1=czero*xm/(denom*phib)
+ cdj2=0.0d0
+ cdj2=0.0d0
+ 230 cdbo=tau*csat*evd/vte
+ cdb1=cdbo/vte
+ twovte=2.0d0*vte
+ geq2=geq/twovte
+ if (vd.le.0.0d0) geq2=0.0d0
+ trivte=3.0d0*vte
+ geq3=geq2/trivte
+ cdb2=cdb1/twovte
+ value(locd)=cdj1
+ value(locd+1)=cdj2
+ value(locd+2)=cdbo
+ value(locd+3)=cdb1
+ value(locd+4)=geq2
+ value(locd+5)=geq3
+ value(locd+6)=cdb2
+ loc=nodplc(loc)
+ go to 210
+c
+c finished
+c
+ 300 return
+ end
+ subroutine disto(loco)
+ implicit double precision (a-h,o-z)
+c
+c this routine performs the small-signal distortion analysis.
+c
+ common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem,
+ 1 isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize,
+ 2 junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr,
+ 3 nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1,
+ 4 lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd,
+ 5 imynl,imvn,lcvn,loutpt,nsnod,nsmat,nsval,icnod,icmat,icval
+ common /miscel/ atime,aprog(3),adate,atitle(10),defl,defw,defad,
+ 1 defas,rstats(50),iwidth,lwidth,nopage
+ common /cirdat/ locate(50),jelcnt(50),nunods,ncnods,numnod,nstop,
+ 1 nut,nlt,nxtrm,ndist,ntlin,ibr,numvs
+ common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet,
+ 1 xmu,mode,modedc,icalc,initf,method,iord,maxord,noncon,iterno,
+ 2 itemno,nosolv,ipostp,iscrch
+ common /knstnt/ twopi,xlog2,xlog10,root2,rad,boltz,charge,ctok,
+ 1 gmin,reltol,abstol,vntol,trtol,chgtol,eps0,epssil,epsox
+ common /flags/ iprnta,iprntl,iprntm,iprntn,iprnto,limtim,limpts,
+ 1 lvlcod,lvltim,itl1,itl2,itl3,itl4,itl5,igoof,nogo,keof
+ common /ac/ fstart,fstop,fincr,skw2,refprl,spw2,jacflg,idfreq,
+ 1 inoise,nosprt,nosout,nosin,idist,idprt
+ common /blank/ value(1000)
+ integer nodplc(64)
+ complex*16 cvalue(32)
+ equivalence (value(1),nodplc(1),cvalue(1))
+c
+c
+ complex*16 difvn1,difvn2,difvn3,difvi1,difvi2,difvi3,dsgo2,dsgm2,
+ 1 dsgmu2,dsgpi2,dscb1,dscb1r,dscje1,dscjc1,disto1,disto2,disto3,
+ 2 dsgmo2,dgm2o3,dgmo23,bew,cew,bcw,be2w,ce2w,bc2w,bew2,cew2,
+ 3 bcw2,bew12,cew12,bcw12,dscdb1,dscdj1,dsg2,cvabe,cvabc,cvace,
+ 4 cvout,cvdist
+ dimension distit(4)
+ dimension vdo(2,12)
+ complex*16 cvdo(12)
+ equivalence (cvdo(1),vdo(1,1))
+ data distit / 8hdistorti, 8hon analy, 8hsis , 8h /
+ icvw1=ld1
+ icv2w1=icvw1+nstop
+ icvw2=icv2w1+nstop
+ icvw12=icvw2+nstop
+ icvadj=icvw12+nstop
+ iprnt=0
+ if (icalc.ge.2) go to 10
+ idnp=nodplc(idist+2)
+ idnn=nodplc(idist+3)
+ locv=nodplc(idist+1)
+ rload=1.0d0/value(locv+1)
+ kntr=1
+ 10 if (idprt.eq.0) go to 30
+ if (kntr.gt.icalc) go to 30
+ iprnt=1
+ kntr=kntr+idprt
+ call title(0,lwidth,1,distit)
+ 30 freq1=dreal(cvalue(loco+1))
+ freq2=skw2*freq1
+ call copy16(cvalue(lcvn+1),cvalue(icvw1+1),nstop)
+ cvout=cvalue(icvw1+idnp)-cvalue(icvw1+idnn)
+ call magphs(cvout,omag,ophase)
+c
+c begin the distortion analysis
+c
+ do 1000 kdisto=1,7
+ cvdist=dcmplx(0.0d0,0.0d0)
+ go to (1000,110,120,130,140,160,170),kdisto
+ 110 freqd=2.0d0*freq1
+ arg=dsqrt(2.0d0*rload*refprl)/(omag*omag)
+ if (iprnt.eq.0) go to 200
+ write (6,111) freq1,freqd,omag,ophase
+ 111 format (///5x,'2nd harmonic distortion',30x,'freq1 = ',1pd9.2,
+ 1 ' hz'//5x,'distortion frequency ',d9.2,' hz',16x,
+ 2 'mag ',d9.3,3x,'phs ',0pf7.2)
+ go to 200
+ 120 freqd=3.0d0*freq1
+ arg=2.0d0*rload*refprl/(omag*omag*omag)
+ if (iprnt.eq.0) go to 200
+ write (6,121) freq1,freqd,omag,ophase
+ 121 format (1h1,4x,'3rd harmonic distortion',30x,'freq1 = ',1pd9.2,
+ 1 ' hz'//5x,'distortion frequency ',d9.2,' hz',16x,
+ 2 'mag ',d9.3,3x,'phs ',0pf7.2)
+ go to 200
+ 130 freqd=freq2
+ go to 200
+ 140 freqd=freq1-freq2
+ arg=dsqrt(2.0d0*rload*refprl)*spw2/(omag*omag)
+ if (iprnt.eq.0) go to 200
+ write (6,151) freq1,freq2,freqd,omag,ophase,ow2mag,ow2phs
+ 151 format (1h1,4x,'2nd order intermodulation difference component',
+ 1 7x,'freq1 = ',1pd9.2,' hz',15x,'freq2 = ',d9.2,' hz'//
+ 2 5x,'distortion frequency ',d9.2,' hz',16x,'mag ',
+ 3 d9.3,3x,'phs ',0pf7.2,9x,'mag ',1pd9.3,3x,'phs ',0pf7.2)
+ go to 200
+ 160 freqd=freq1+freq2
+ arg=dsqrt(2.0d0*rload*refprl)*spw2/(omag*omag)
+ if (iprnt.eq.0) go to 200
+ write (6,161) freq1,freq2,freqd,omag,ophase,ow2mag,ow2phs
+ 161 format (1h1,4x,'2nd order intermodulation sum component',
+ 1 14x,'freq1 = ',1pd9.2,' hz',15x,'freq2 = ',d9.2,' hz'//
+ 2 5x,'distortion frequency ',d9.2,' hz',16x,'mag ',
+ 3 d9.3,3x,'phs ',0pf7.2,9x,'mag ',1pd9.3,3x,'phs ',0pf7.2)
+ go to 200
+ 170 freqd=2.0d0*freq1-freq2
+ arg=2.0d0*rload*refprl*spw2/(omag*omag*omag)
+ if (iprnt.eq.0) go to 200
+ write (6,171) freq1,freq2,freqd,omag,ophase,ow2mag,ow2phs
+ 171 format (1h1,4x,'3rd order intermodulation difference component',
+ 1 7x,'freq1 = ',1pd9.2,' hz',15x,'freq2 = ',d9.2,' hz'//
+ 2 5x,'distortion frequency ',d9.2,' hz',16x,'mag ',
+ 3 d9.3,3x,'phs ',0pf7.2,9x,'mag ',1pd9.3,3x,'phs ',0pf7.2)
+c
+c load and decompose y matrix
+c
+ 200 omega=twopi*freqd
+ igoof=0
+ call acload
+ call acdcmp
+ if (igoof.eq.0) go to 220
+ write (6,211) igoof,freqd
+ 211 format('0warning: underflow ',i4,' time(s) in distortion analysis
+ 1 at freq = ',1pd9.3,' hz')
+ igoof=0
+ 220 if (kdisto.eq.4) go to 710
+c
+c obtain adjoint solution
+c
+ call zero8(value(lvn+1),nstop)
+ call zero8(value(imvn+1),nstop)
+ value(lvn+idnp)=-1.0d0
+ value(lvn+idnn)=+1.0d0
+ call acasol
+ call copy16(cvalue(lcvn+1),cvalue(icvadj+1),nstop)
+ call zero8(value(lvn+1),nstop)
+ call zero8(value(imvn+1),nstop)
+c
+c bjts
+c
+ if (jelcnt(12).eq.0) go to 500
+ ititle=0
+ 301 format (////1x,'bjt distortion components'//1x,'name',11x,'gm',
+ 1 8x,'gpi',7x,'go',8x,'gmu',6x,'gmo2',7x,'cb',8x,'cbr',7x,'cje',
+ 2 7x,'cjc',6x,'total')
+ 311 format (////1x,'bjt distortion components'//1x,'name',11x,'gm',
+ 1 8x,'gpi',7x,'go',8x,'gmu',6x,'gmo2',7x,'cb',8x,'cbr',7x,'cje',
+ 2 7x,'cjc',6x,'gm203',5x,'gmo23',5x,'total')
+ 320 loc=locate(12)
+ 330 if (loc.eq.0) go to 500
+ locv=nodplc(loc+1)
+ loct=lx0+nodplc(loc+22)
+ locd=ld0+nodplc(loc+23)
+ node1=nodplc(loc+5)
+ node2=nodplc(loc+6)
+ node3=nodplc(loc+7)
+ cje1=value(locd)
+ cje2=value(locd+1)
+ cjc1=value(locd+2)
+ cjc2=value(locd+3)
+ go2=value(locd+4)
+ gmo2=value(locd+5)
+ gm2=value(locd+6)
+ gmu2=value(locd+7)
+ gpi2=value(locd+8)
+ cb1=value(locd+11)
+ cb1r=value(locd+12)
+ go3=value(locd+13)
+ gmo23=value(locd+14)
+ gm2o3=value(locd+15)
+ gm3=value(locd+16)
+ gmu3=value(locd+17)
+ gpi3=value(locd+18)
+ cb2=value(locd+19)
+ cb2r=value(locd+20)
+ bew=cvalue(icvw1+node2)-cvalue(icvw1+node3)
+ cew=cvalue(icvw1+node1)-cvalue(icvw1+node3)
+ bcw=cvalue(icvw1+node2)-cvalue(icvw1+node1)
+ if (kdisto.eq.2) go to 370
+ be2w=cvalue(icv2w1+node2)-cvalue(icv2w1+node3)
+ ce2w=cvalue(icv2w1+node1)-cvalue(icv2w1+node3)
+ bc2w=cvalue(icv2w1+node2)-cvalue(icv2w1+node1)
+ if (kdisto.eq.3) go to 380
+ bew2=cvalue(icvw2+node2)-cvalue(icvw2+node3)
+ cew2=cvalue(icvw2+node1)-cvalue(icvw2+node3)
+ bcw2=cvalue(icvw2+node2)-cvalue(icvw2+node1)
+ if (kdisto.eq.5) go to 390
+ if (kdisto.eq.6) go to 400
+ bew12=cvalue(icvw12+node2)-cvalue(icvw12+node3)
+ cew12=cvalue(icvw12+node1)-cvalue(icvw12+node3)
+ bcw12=cvalue(icvw12+node2)-cvalue(icvw12+node1)
+ go to 410
+c
+c calculate hd2 current generators
+c
+ 370 difvn1=0.5d0*cew*cew
+ difvn2=0.5d0*bew*bew
+ difvn3=0.5d0*bcw*bcw
+ dsgmo2=gmo2*0.5d0*bew*cew
+ go to 420
+c
+c calculate hd3 current generators
+c
+ 380 difvi1=0.50d0*cew*ce2w
+ difvn1=0.25d0*cew*cew*cew
+ difvi2=0.50d0*bew*be2w
+ difvn2=0.25d0*bew*bew*bew
+ difvi3=0.50d0*bcw*bc2w
+ difvn3=0.25d0*bcw*bcw*bcw
+ dsgmo2=gmo2*(bew*ce2w+be2w*cew)*0.5d0
+ go to 430
+c
+c calculate im2d current generators
+c
+ 390 difvn1=cew*dconjg(cew2)
+ difvn2=bew*dconjg(bew2)
+ difvn3=bcw*dconjg(bcw2)
+ dsgmo2=gmo2*0.5d0*(bew*dconjg(cew2)+cew*dconjg(bew2))
+ go to 420
+c
+c calculate im2s current generators
+c
+ 400 difvn1=cew*cew2
+ difvn2=bew*bew2
+ difvn3=bcw*bcw2
+ dsgmo2=gmo2*0.5d0*(bew*cew2+bew2*cew)
+ go to 420
+c
+c calculate im3 current generators
+c
+ 410 difvi1=0.5d0*(ce2w*dconjg(cew2)+cew*cew12)
+ difvi2=0.5d0*(be2w*dconjg(bew2)+bew*bew12)
+ difvi3=0.5d0*(bc2w*dconjg(bcw2)+bcw*bcw12)
+ difvn1=cew*cew*dconjg(cew2)*0.75d0
+ difvn2=bew*bew*dconjg(bew2)*0.75d0
+ difvn3=bcw*bcw*dconjg(bcw2)*0.75d0
+ dsgmo2=gmo2*0.5d0*(dconjg(bew2)*ce2w+bew*cew12+dconjg(cew2)*be2w+
+ 1 cew*bew12)
+ go to 430
+c
+ 420 dsgo2=go2*difvn1
+ dsgm2=gm2*difvn2
+ dsgmu2=gmu2*difvn3
+ dsgpi2=gpi2*difvn2
+ dscb1=0.5d0*cb1*omega*dcmplx(-dimag(difvn2),dreal(difvn2))
+ dscb1r=0.5d0*cb1r*omega*dcmplx(-dimag(difvn3),dreal(difvn3))
+ dscje1=0.5d0*cje1*omega*dcmplx(-dimag(difvn2),dreal(difvn2))
+ dscjc1=0.5d0*cjc1*omega*dcmplx(-dimag(difvn3),dreal(difvn3))
+ go to 440
+c
+ 430 dsgo2=2.0d0*go2*difvi1+go3*difvn1
+ dsgm2=2.0d0*gm2*difvi2+gm3*difvn2
+ dsgmu2=2.0d0*gmu2*difvi3+gmu3*difvn3
+ dsgpi2=2.0d0*gpi2*difvi2+gpi3*difvn2
+ dscb1=omega*(cb1*difvi2+cb2*difvn2/3.0d0)
+ dscb1=dcmplx(-dimag(dscb1),dreal(dscb1))
+ dscb1r=omega*(cb1r*difvi3+cb2r*difvn3/3.0d0)
+ dscb1r=dcmplx(-dimag(dscb1r),dreal(dscb1r))
+ dscje1=omega*(cje1*difvi2+cje2*difvn2/3.0d0)
+ dscje1=dcmplx(-dimag(dscje1),dreal(dscje1))
+ dscjc1=omega*(cjc1*difvi3+cjc2*difvn3/3.0d0)
+ dscjc1=dcmplx(-dimag(dscjc1),dreal(dscjc1))
+c
+c determine contribution of each distortion source
+c
+ 440 cvabe=cvalue(icvadj+node2)-cvalue(icvadj+node3)
+ cvabc=cvalue(icvadj+node2)-cvalue(icvadj+node1)
+ cvace=cvalue(icvadj+node1)-cvalue(icvadj+node3)
+ disto1=dsgm2+dsgo2+dsgmo2
+ disto2=dsgpi2+dscb1+dscje1
+ disto3=dsgmu2+dscb1r+dscjc1
+ cvdo(1)=dsgm2*cvace*arg
+ cvdo(2)=dsgpi2*cvabe*arg
+ cvdo(3)=dsgo2*cvace*arg
+ cvdo(4)=dsgmu2*cvabc*arg
+ cvdo(5)=dsgmo2*cvace*arg
+ cvdo(6)=dscb1*cvabe*arg
+ cvdo(7)=dscb1r*cvabc*arg
+ cvdo(8)=dscje1*cvabe*arg
+ cvdo(9)=dscjc1*cvabc*arg
+ if (kdisto.eq.3) go to 450
+ if (kdisto.eq.7) go to 460
+ cvdo(10)=cvdo(1)+cvdo(2)+cvdo(3)+cvdo(4)+cvdo(5)+cvdo(6)+cvdo(7)+
+ 1 cvdo(8)+cvdo(9)
+ cvdist=cvdist+cvdo(10)
+ if (iprnt.eq.0) go to 480
+ do 445 j=1,10
+ call magphs(cvdo(j),xmag,xphs)
+ cvdo(j)=dcmplx(xmag,xphs)
+ 445 continue
+ if (ititle.eq.0) write (6,301)
+ ititle=1
+ write (6,446) value(locv),(vdo(1,j),j=1,10)
+ 446 format(1h0,a8,'mag',1p12d10.3)
+ write (6,447) (vdo(2,j),j=1,10)
+ 447 format(9x,'phs',12(1x,f7.2,2x))
+ go to 480
+ 450 dgm2o3=gm2o3*cew*bew*bew*0.25d0
+ dgmo23=gmo23*bew*cew*cew*0.25d0
+ go to 470
+ 460 dgm2o3=gm2o3*(0.5d0*bew*dconjg(bew2)*cew+0.25d0*bew*bew*
+ 1 dconjg(cew2))
+ dgmo23=gmo23*(0.5d0*cew*dconjg(cew2)*bew+0.25d0*cew*cew*
+ 1 dconjg(bew2))
+ 470 disto1=disto1+dgm2o3+dgmo23
+ cvdo(10)=dgm2o3*cvace*arg
+ cvdo(11)=dgmo23*cvace*arg
+ cvdo(12)=cvdo(1)+cvdo(2)+cvdo(3)+cvdo(4)+cvdo(5)+cvdo(6)+cvdo(7)+
+ 1 cvdo(8)+cvdo(9)+cvdo(10)+cvdo(11)
+ cvdist=cvdist+cvdo(12)
+ if (iprnt.eq.0) go to 480
+ do 475 j=1,12
+ call magphs(cvdo(j),xmag,xphs)
+ cvdo(j)=dcmplx(xmag,xphs)
+ 475 continue
+ if (ititle.eq.0) write (6,311)
+ ititle=1
+ write (6,446) value(locv),(vdo(1,j),j=1,12)
+ write (6,447) (vdo(2,j),j=1,12)
+ 480 value(lvn+node1)=value(lvn+node1)
+ 1 -dreal(disto1-disto3)
+ value(lvn+node2)=value(lvn+node2)
+ 1 -dreal(disto2+disto3)
+ value(lvn+node3)=value(lvn+node3)
+ 1 +dreal(disto1+disto2)
+ value(imvn+node1)=value(imvn+node1)
+ 1 -dimag(disto1-disto3)
+ value(imvn+node2)=value(imvn+node2)
+ 1 -dimag(disto2+disto3)
+ value(imvn+node3)=value(imvn+node3)
+ 1 +dimag(disto1+disto2)
+ loc=nodplc(loc)
+ go to 330
+c
+c junction diodes
+c
+ 500 if (jelcnt(11).eq.0) go to 700
+ ititle=0
+ 501 format (////1x,'diode distortion components'//1x,'name',
+ 1 11x,'geq',7x,'cb',8x,'cj',7x,'total')
+ 510 loc=locate(11)
+ 520 if (loc.eq.0) go to 700
+ locv=nodplc(loc+1)
+ node1=nodplc(loc+2)
+ node2=nodplc(loc+3)
+ node3=nodplc(loc+4)
+ locm=nodplc(loc+5)
+ locm=nodplc(locm+1)
+ loct=lx0+nodplc(loc+11)
+ locd=ld0+nodplc(loc+12)
+ cdj1=value(locd)
+ cdj2=value(locd+1)
+ cdb1=value(locd+3)
+ geq2=value(locd+4)
+ geq3=value(locd+5)
+ cdb2=value(locd+6)
+ bew=cvalue(icvw1+node3)-cvalue(icvw1+node2)
+ if (kdisto.eq.2) go to 540
+ be2w=cvalue(icv2w1+node3)-cvalue(icv2w1+node2)
+ if (kdisto.eq.3) go to 550
+ bew2=cvalue(icvw2+node3)-cvalue(icvw2+node2)
+ if (kdisto.eq.5) go to 560
+ if (kdisto.eq.6) go to 570
+ bew12=cvalue(icvw12+node3)-cvalue(icvw12+node2)
+ go to 580
+c
+c calculate hd2 current generators
+c
+ 540 difvn1=0.5d0*bew*bew
+ go to 590
+c
+c calculate hd3 current generators
+c
+ 550 difvi1=0.5d0*bew*be2w
+ difvn1=0.25d0*bew*bew*bew
+ go to 600
+c
+c calculate im2d current generators
+c
+ 560 difvn1=bew*dconjg(bew2)
+ go to 590
+c
+c calculate im2s current generators
+c
+ 570 difvn1=bew*bew2
+ go to 590
+c
+c calculate im3 current generators
+c
+ 580 difvi1=0.5d0*(be2w*dconjg(bew2)+bew*bew12)
+ difvn1=bew*bew*dconjg(bew2)*0.75d0
+ go to 600
+ 590 dsg2=geq2*difvn1
+ dscdb1=0.5d0*cdb1*omega*dcmplx(-dimag(difvn1),dreal(difvn1))
+ dscdj1=0.5d0*cdj1*omega*dcmplx(-dimag(difvn1),dreal(difvn1))
+ go to 610
+c
+ 600 dsg2=2.0d0*geq2*difvi1+geq3*difvn1
+ dscdb1=omega*(cdb1*difvi1+cdb2*difvn1/3.0d0)
+ dscdb1=dcmplx(-dimag(dscdb1),dreal(dscdb1))
+ dscdj1=omega*(cdj1*difvi1+cdj2*difvn1/3.0d0)
+ dscdj1=dcmplx(-dimag(dscdj1),dreal(dscdj1))
+c
+c determine contribution of each distortion source
+c
+ 610 cvabe=cvalue(icvadj+node3)-cvalue(icvadj+node2)
+ disto1=dsg2+dscdb1+dscdj1
+ cvdo(1)=dsg2*cvabe*arg
+ cvdo(2)=dscdb1*cvabe*arg
+ cvdo(3)=dscdj1*cvabe*arg
+ cvdo(4)=cvdo(1)+cvdo(2)+cvdo(3)
+ cvdist=cvdist+cvdo(4)
+ if (iprnt.eq.0) go to 680
+ do 670 j=1,4
+ call magphs(cvdo(j),xmag,xphs)
+ cvdo(j)=dcmplx(xmag,xphs)
+ 670 continue
+ if (ititle.eq.0) write (6,501)
+ ititle=1
+ write (6,446) value(locv),(vdo(1,j),j=1,4)
+ write (6,447) (vdo(2,j),j=1,4)
+ 680 value(lvn+node2)=value(lvn+node2)+dreal(disto1)
+ value(lvn+node3)=value(lvn+node3)-dreal(disto1)
+ value(imvn+node2)=value(imvn+node2)+dimag(disto1)
+ value(imvn+node3)=value(imvn+node3)-dimag(disto1)
+ loc=nodplc(loc)
+ go to 520
+c
+c obtain total distortion solution if necessary
+c
+ 700 go to (1000,710,790,710,710,840,860),kdisto
+ 710 call acsol
+c
+c store solution, print and store answers
+c
+ 760 go to (1000,770,790,800,820,840,860),kdisto
+ 770 call copy16(cvalue(lcvn+1),cvalue(icv2w1+1),nstop)
+ call magphs(cvdist,o2mag,o2phs)
+ if (iprnt.eq.0) go to 900
+ o2log=20.0d0*dlog10(o2mag)
+ write (6,781) o2mag,o2phs,o2log
+ 781 format (///5x,'hd2 magnitude ',1pd10.3,5x,'phase ',0pf7.2,
+ 1 5x,'= ',f7.2,' db')
+ go to 900
+ 790 call magphs(cvdist,o3mag,o3phs)
+ if (iprnt.eq.0) go to 900
+ o3log=20.0d0*dlog10(o3mag)
+ write (6,791) o3mag,o3phs,o3log
+ 791 format (///5x,'hd3 magnitude ',1pd10.3,5x,'phase ',0pf7.2,
+ 1 5x,'= ',f7.2,' db')
+ go to 900
+ 800 call copy16(cvalue(lcvn+1),cvalue(icvw2+1),nstop)
+ cvout=cvalue(icvw2+idnp)-cvalue(icvw2+idnn)
+ call magphs(cvout,ow2mag,ow2phs)
+ go to 1000
+ 820 call copy16(cvalue(lcvn+1),cvalue(icvw12+1),nstop)
+ 840 call magphs(cvdist,o12mag,o12phs)
+ if (iprnt.eq.0) go to 900
+ o12log=20.0d0*dlog10(o12mag)
+ if (kdisto.eq.6) go to 850
+ write (6,841) o12mag,o12phs,o12log
+ 841 format (///5x,'im2d magnitude ',1pd10.3,5x,'phase ',0pf7.2,
+ 1 5x,'= ',f7.2,' db')
+ go to 900
+ 850 write (6,851) o12mag,o12phs,o12log
+ 851 format (///5x,'im2s magnitude ',1pd10.3,5x,'phase ',0pf7.2,
+ 1 5x,'= ',f7.2,' db')
+ go to 900
+ 860 call magphs(cvdist,o21mag,o21phs)
+ if (iprnt.eq.0) go to 900
+ o21log=20.0d0*dlog10(o21mag)
+ write (6,861) o21mag,o21phs,o21log
+ 861 format (///5x,'im3 magnitude ',1pd10.3,5x,'phase ',0pf7.2,
+ 1 5x,'= ',f7.2,' db')
+ cma=dabs(4.0d0*o21mag*dcos((o21phs-ophase)/rad))
+ cma=dmax1(cma,1.0d-20)
+ cmp=dabs(4.0d0*o21mag*dsin((o21phs-ophase)/rad))
+ cmp=dmax1(cmp,1.0d-20)
+ cmalog=20.0d0*dlog10(cma)
+ cmplog=20.0d0*dlog10(cmp)
+ write (6,866)
+ 866 format (////5x,'approximate cross modulation components')
+ write (6,871) cma,cmalog
+ 871 format (/5x,'cma magnitude ',1pd10.3,24x,'= ',0pf7.2,' db')
+ write (6,881) cmp,cmplog
+ 881 format (/5x,'cmp magnitude ',1pd10.3,24x,'= ',0pf7.2,' db')
+c
+c save distortion outputs
+c
+ 900 iflag=kdisto+2
+ if (iflag.ge.7) iflag=iflag-1
+ loc=locate(45)
+ 910 if (loc.eq.0) go to 1000
+ if (nodplc(loc+5).ne.iflag) go to 920
+ iseq=nodplc(loc+4)
+ cvalue(loco+iseq)=cvdist
+ 920 loc=nodplc(loc)
+ go to 910
+ 1000 continue
+c
+c finished
+c
+ 2000 return
+ end
projects / unix-history / commitdiff