![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to acans.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [CSRG BSD Unix] / 3BSD / cmd / spice|
Return to acans.f CVS log | Up to [CSRG BSD Unix] / 3BSD / cmd / spice |
BSD 3.0
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
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.