+ subroutine ovtpvt
+ implicit double precision (a-h,o-z)
+c
+c
+c this routine generates the requested tabular listings of analysis
+c results. it calls plot to generate line-printer plots.
+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 /miscel/ atime,aprog(3),adate,atitle(10),defl,defw,defad,
+ 1 defas,rstats(50),iwidth,lwidth,nopage
+ common /dc/ tcstar(2),tcstop(2),tcincr(2),icvflg,itcelm(2),kssop,
+ 1 kinel,kidin,kovar,kidout
+ common /ac/ fstart,fstop,fincr,skw2,refprl,spw2,jacflg,idfreq,
+ 1 inoise,nosprt,nosout,nosin,idist,idprt
+ common /tran/ tstep,tstop,tstart,delmax,tdmax,forfre,jtrflg
+ common /outinf/ xincr,string(15),xstart,yvar(8),itab(8),itype(8),
+ 1 ilogy(8),npoint,numout,kntr,numdgt
+ common /blank/ value(1000)
+ integer nodplc(64)
+ complex*16 cvalue(32)
+ equivalence (value(1),nodplc(1),cvalue(1))
+c
+ complex*16 cval
+ dimension prform(3)
+ dimension subtit(4,3)
+ data subtit / 8hdc trans, 8hfer curv, 8hes , 8h ,
+ 1 8htransien, 8ht analys, 8his , 8h ,
+ 2 8hac analy, 8hsis , 8h , 8h /
+ data prform / 8h(1pe11.3, 8h,2x,8e00, 8h.00) /
+ data aper,rprn / 1h., 1h) /
+c
+ call second(t1)
+ if (icalc.le.0) go to 1000
+ call crunch
+ if (nogo.lt.0) go to 1000
+c
+c construct format statement to be used for printing the outputs
+c
+ ifract=max0(numdgt-1,0)
+ ifwdth=ifract+9
+ ipos=15
+ call alfnum(ifwdth,prform,ipos)
+ call move(prform,ipos,aper,1,1)
+ ipos=ipos+1
+ call alfnum(ifract,prform,ipos)
+ call move(prform,ipos,rprn,1,1)
+c
+ noprln=min0(8,(lwidth-12)/ifwdth)
+ if (mode-2) 50,60,300
+ 50 numout=jelcnt(41)+1
+ go to 70
+ 60 numout=jelcnt(42)+1
+c
+c dc and transient analysis printing
+c
+ 70 loc=locate(30+mode)
+ 80 if (loc.eq.0) go to 200
+ kntr=min0(noprln,nodplc(loc+3))
+ if (kntr.le.0) go to 120
+ call title(1,lwidth,1,subtit(1,mode))
+ call setprn(loc)
+c
+c get buffer space
+c
+ call getm8(locx,npoint)
+ call getm8(locy,kntr*npoint)
+c
+c interpolate outputs
+c
+ call ntrpl8(locx,locy,numpnt)
+c
+c print outputs
+c
+ do 100 i=1,numpnt
+ xvar=value(locx+i)
+ locyt=locy
+ do 90 k=1,kntr
+ yvar(k)=value(locyt+i)
+ locyt=locyt+npoint
+ 90 continue
+ write (6,prform) xvar,(yvar(k),k=1,kntr)
+ 100 continue
+ write (6,111)
+ 111 format(1hy)
+ call clrmem(locx)
+ call clrmem(locy)
+ 120 loc=nodplc(loc)
+ go to 80
+c
+c dc and transient analysis plotting
+c
+ 200 loc=locate(35+mode)
+ 210 if (loc.eq.0) go to 250
+ kntr=nodplc(loc+3)
+ if (kntr.le.0) go to 220
+ locv=nodplc(loc+1)
+ call title(1,lwidth,1,subtit(1,mode))
+ call setplt(loc)
+c
+c get buffer space
+c
+ call getm8(locx,npoint)
+ call getm8(locy,kntr*npoint)
+c
+c interpolate outputs and load plot buffers
+c
+ call ntrpl8(locx,locy,numpnt)
+ call plot(numpnt,locx,locy,locv)
+ call clrmem(locx)
+ call clrmem(locy)
+ 220 loc=nodplc(loc)
+ go to 210
+c
+c fourier analysis
+c
+ 250 if (mode.eq.1) go to 1000
+ if (nfour.eq.0) go to 1000
+ if (nogo.ne.0) go to 1000
+ call fouran
+ go to 1000
+c
+c ac analysis printing
+c
+ 300 numout=jelcnt(43)+jelcnt(44)+jelcnt(45)+1
+ do 599 id=33,35
+ loc=locate(id)
+ 320 if (loc.eq.0) go to 599
+ kntr=min0(noprln,nodplc(loc+3))
+ if (kntr.le.0) go to 595
+ call title(1,lwidth,1,subtit(1,mode))
+ call setprn(loc)
+c
+c print ac outputs
+c
+ lout=loutpt
+ do 590 i=1,icalc
+ xvar=dreal(cvalue(lout+1))
+ do 500 k=1,kntr
+ iseq=itab(k)
+ iseq=nodplc(iseq+4)
+ cval=cvalue(lout+iseq)
+ ktype=itype(k)
+ go to (450,450,430,440,450,450), ktype
+ 430 yvar(k)=dreal(cval)
+ go to 500
+ 440 yvar(k)=dimag(cval)
+ go to 500
+ 450 call magphs(cval,xmag,xphs)
+ go to (460,460,430,440,470,465), ktype
+ 460 yvar(k)=xmag
+ go to 500
+ 465 yvar(k)=20.0d0*dlog10(xmag)
+ go to 500
+ 470 yvar(k)=xphs
+ 500 continue
+ lout=lout+numout
+ 580 write (6,prform) xvar,(yvar(k),k=1,kntr)
+ 590 continue
+ write (6,111)
+ 595 loc=nodplc(loc)
+ go to 320
+ 599 continue
+c
+c ac analysis plotting
+c
+ do 760 id=38,40
+ loc=locate(id)
+ 610 if (loc.eq.0) go to 760
+ kntr=nodplc(loc+3)
+ if (kntr.le.0) go to 750
+ locv=nodplc(loc+1)
+ call title(1,lwidth,1,subtit(1,mode))
+ call setplt(loc)
+c
+ call getm8(locx,icalc)
+ call getm8(locy,kntr*icalc)
+c
+c load plot buffers
+c
+ lout=loutpt
+ do 710 i=1,icalc
+ xvar=dreal(cvalue(lout+1))
+ locyt=locy
+ do 700 k=1,kntr
+ iseq=itab(k)
+ iseq=nodplc(iseq+4)
+ cval=cvalue(lout+iseq)
+ ktype=itype(k)
+ go to (670,670,650,660,670,670), ktype
+ 650 yvr=dreal(cval)
+ go to 695
+ 660 yvr=dimag(cval)
+ go to 695
+ 670 call magphs(cval,xmag,xphs)
+ go to (680,680,650,660,690,685), ktype
+ 680 yvr=dlog10(xmag)
+ go to 695
+ 685 yvr=20.0d0*dlog10(xmag)
+ go to 695
+ 690 yvr=xphs
+ 695 value(locyt+i)=yvr
+ locyt=locyt+icalc
+ 700 continue
+ value(locx+i)=xvar
+ lout=lout+numout
+ 710 continue
+ call plot(icalc,locx,locy,locv)
+ call clrmem(locx)
+ call clrmem(locy)
+ 750 loc=nodplc(loc)
+ go to 610
+ 760 continue
+c
+c finished
+c
+ 1000 call clrmem(loutpt)
+ call second(t2)
+ rstats(11)=rstats(11)+t2-t1
+ return
+ end
+ subroutine ntrpl8(locx,locy,numpnt)
+ implicit double precision (a-h,o-z)
+c
+c this routine interpolates the analysis data to obtain the values
+c printed and/or plotted, using linear interpolation.
+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 /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 /outinf/ xincr,string(15),xstart,yvar(8),itab(8),itype(8),
+ 1 ilogy(8),npoint,numout,kntr,numdgt
+ common /blank/ value(1000)
+ integer nodplc(64)
+ complex*16 cvalue(32)
+ equivalence (value(1),nodplc(1),cvalue(1))
+ if(mode.ne.1) go to 4
+ numpnt=icalc
+ loco=loutpt
+ do 3 i=1,numpnt
+ locyt=locy
+ value(locx+i)=value(loco+1)
+ do 2 k=1,kntr
+ iseq=itab(k)
+ iseq=nodplc(iseq+4)
+ value(locyt+i)=value(loco+iseq)
+ locyt=locyt+npoint
+ 2 continue
+ loco=loco+numout
+ 3 continue
+ return
+ 4 continue
+ xvar=xstart
+ xvtol=xincr*1.0d-5
+ ippnt=0
+ icpnt=2
+ loco1=loutpt
+ loco2=loco1+numout
+ if (icalc.lt.2) go to 50
+ 10 x1=value(loco1+1)
+ x2=value(loco2+1)
+ dx1x2=x1-x2
+ 20 if (xincr.lt.0.0d0) go to 24
+ if (xvar.le.(x2+xvtol)) go to 30
+ go to 28
+ 24 if (xvar.ge.(x2+xvtol)) go to 30
+ 28 if (icpnt.ge.icalc) go to 100
+ icpnt=icpnt+1
+ loco1=loco2
+ loco2=loco1+numout
+ go to 10
+ 30 ippnt=ippnt+1
+ value(locx+ippnt)=xvar
+ dxx1=xvar-x1
+ locyt=locy
+ do 40 i=1,kntr
+ iseq=itab(i)
+ iseq=nodplc(iseq+4)
+ v1=value(loco1+iseq)
+ v2=value(loco2+iseq)
+ yvr=v1+(v1-v2)*dxx1/dx1x2
+ tol=dmin1(dabs(v1),dabs(v2))*1.0d-10
+ if (dabs(yvr).le.tol) yvr=0.0d0
+ value(locyt+ippnt)=yvr
+ locyt=locyt+npoint
+ 40 continue
+ if (ippnt.ge.npoint) go to 100
+ xvar=xstart+dfloat(ippnt)*xincr
+ if (dabs(xvar).ge.dabs(xvtol)) go to 20
+ xvar=0.0d0
+ go to 20
+c
+c special handling if icalc = 1
+c
+c... icalc=1; just copy over the single point and return
+ 50 ippnt=1
+ value(locx+ippnt)=xvar
+ locyt=locy
+ do 60 i=1,kntr
+ iseq=itab(i)
+ iseq=nodplc(iseq+4)
+ value(locyt+ippnt)=value(loco1+iseq)
+ locyt=locyt+npoint
+ 60 continue
+ go to 100
+c
+c return
+c
+ 100 numpnt=ippnt
+ return
+ end
+ subroutine setprn(loc)
+ implicit double precision (a-h,o-z)
+c
+c this routine formats the column headers for tabular listings of
+c output variables.
+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 /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 /dc/ tcstar(2),tcstop(2),tcincr(2),icvflg,itcelm(2),kssop,
+ 1 kinel,kidin,kovar,kidout
+ common /ac/ fstart,fstop,fincr,skw2,refprl,spw2,jacflg,idfreq,
+ 1 inoise,nosprt,nosout,nosin,idist,idprt
+ common /tran/ tstep,tstop,tstart,delmax,tdmax,forfre,jtrflg
+ common /outinf/ xincr,string(15),xstart,yvar(8),itab(8),itype(8),
+ 1 ilogy(8),npoint,numout,kntr,numdgt
+ common /blank/ value(1000)
+ integer nodplc(64)
+ complex*16 cvalue(32)
+ equivalence (value(1),nodplc(1),cvalue(1))
+c
+ data ablnk, atimex, afreq / 1h , 6h time, 6h freq /
+c
+c set limits depending upon the analysis mode
+c
+ if (mode-2) 10,20,30
+ 10 xstart=tcstar(1)
+ xincr=tcincr(1)
+ npoint=icvflg
+ itemp=itcelm(1)
+ loce=nodplc(itemp+1)
+ asweep=value(loce)
+ go to 40
+ 20 xstart=tstart
+ xincr=tstep
+ npoint=jtrflg
+ asweep=atimex
+ go to 40
+ 30 xstart=fstart
+ xincr=fincr
+ npoint=icalc
+ asweep=afreq
+c
+c construct and print the output variable names
+c
+ 40 loct=loc+2
+ ipos=1
+ npos=ipos+numdgt+8
+ do 90 i=1,kntr
+ loct=loct+2
+ itab(i)=nodplc(loct)
+ itype(i)=nodplc(loct+1)
+ call outnam(itab(i),itype(i),string,ipos)
+ if (ipos.ge.npos) go to 70
+ do 60 j=ipos,npos
+ call move(string,j,ablnk,1,1)
+ 60 continue
+ ipos=npos
+ go to 80
+ 70 call move(string,ipos,ablnk,1,1)
+ ipos=ipos+1
+ 80 npos=npos+numdgt+8
+ 90 continue
+ call move(string,ipos,ablnk,1,7)
+ jstop=(ipos+6)/8
+ write (6,91) asweep,(string(j),j=1,jstop)
+ 91 format(/3x,a8,5x,14a8,a4)
+ write (6,101)
+ 101 format(1hx/1h )
+ return
+ end
+ subroutine setplt(loc)
+ implicit double precision (a-h,o-z)
+c
+c this routine generates the 'legend' subheading used to identify
+c individual traces on multi-trace line-printer plots.
+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 /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 /dc/ tcstar(2),tcstop(2),tcincr(2),icvflg,itcelm(2),kssop,
+ 1 kinel,kidin,kovar,kidout
+ common /ac/ fstart,fstop,fincr,skw2,refprl,spw2,jacflg,idfreq,
+ 1 inoise,nosprt,nosout,nosin,idist,idprt
+ common /tran/ tstep,tstop,tstart,delmax,tdmax,forfre,jtrflg
+ common /outinf/ xincr,string(15),xstart,yvar(8),itab(8),itype(8),
+ 1 ilogy(8),npoint,numout,kntr,numdgt
+ common /blank/ value(1000)
+ integer nodplc(64)
+ complex*16 cvalue(32)
+ equivalence (value(1),nodplc(1),cvalue(1))
+c
+ dimension logopt(6)
+ data logopt / 2, 2, 1, 1, 1, 1 /
+ data ablnk, atimex, afreq / 1h , 6h time, 6h freq /
+ data pltsym / 8h*+=$0<>? /
+c
+c set limits depending upon the analysis mode
+c
+ if (mode-2) 10,20,30
+ 10 xstart=tcstar(1)
+ xincr=tcincr(1)
+ npoint=icvflg
+ itemp=itcelm(1)
+ loce=nodplc(itemp+1)
+ asweep=value(loce)
+ go to 40
+ 20 xstart=tstart
+ xincr=tstep
+ npoint=jtrflg
+ asweep=atimex
+ go to 40
+ 30 xstart=fstart
+ xincr=fincr
+ npoint=jacflg
+ asweep=afreq
+c
+c construct and print the output variables with corresponding plot
+c symbols
+c
+ 40 loct=loc+2
+ if (kntr.eq.1) go to 80
+ write (6,41)
+ 41 format('0legend:'/)
+ do 70 i=1,kntr
+ loct=loct+2
+ itab(i)=nodplc(loct)
+ ioutyp=nodplc(loct+1)
+ itype(i)=ioutyp
+ ilogy(i)=1
+ if (mode.le.2) go to 50
+ ilogy(i)=logopt(ioutyp)
+ 50 ipos=1
+ call outnam(itab(i),itype(i),string,ipos)
+ call move(string,ipos,ablnk,1,7)
+ jstop=(ipos+6)/8
+ call move(achar,1,pltsym,i,1)
+ write (6,61) achar,(string(j),j=1,jstop)
+ 61 format(1x,a1,2h: ,5a8)
+ 70 continue
+ 80 if (kntr.ge.2) go to 90
+ itab(1)=nodplc(loc+4)
+ ioutyp=nodplc(loc+5)
+ itype(1)=ioutyp
+ ilogy(1)=1
+ if (mode.le.2) go to 90
+ ilogy(1)=logopt(ioutyp)
+ 90 ipos=1
+ call outnam(itab(1),itype(1),string,ipos)
+ call move(string,ipos,ablnk,1,7)
+ jstop=(ipos+6)/8
+ write (6,101) asweep,(string(j),j=1,jstop)
+ 101 format(1hx/3x,a8,4x,5a8)
+ return
+ end
+ subroutine plot(numpnt,locx,locy,locv)
+ implicit double precision (a-h,o-z)
+c
+c this routine generates the line-printer plots.
+c
+ common /miscel/ atime,aprog(3),adate,atitle(10),defl,defw,defad,
+ 1 defas,rstats(50),iwidth,lwidth,nopage
+ 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 /outinf/ xincr,string(15),xstart,yvar(8),itab(8),itype(8),
+ 1 ilogy(8),npoint,numout,kntr,numdgt
+ common /blank/ value(1000)
+ integer nodplc(64)
+ complex*16 cvalue(32)
+ equivalence (value(1),nodplc(1),cvalue(1))
+c
+c
+ integer xxor
+ dimension ycoor(5,8),icoor(8),delplt(8)
+ dimension agraph(13),aplot(13)
+ dimension asym(2),pmin(8),jcoor(8)
+ data ablnk, aletx, aper / 1h , 1hx, 1h. /
+ data asym1, asym2, arprn / 8h(-------, 8h--------, 1h) /
+ data pltsym / 8h*+=$0<>? /
+c
+c
+ iwide=1
+ nwide=101
+ nwide4=25
+ if(lwidth.gt.80) go to 3
+ iwide=0
+ nwide=57
+ nwide4=14
+ 3 if (numpnt.le.0) go to 400
+ do 5 i=1,13
+ agraph(i)=ablnk
+ 5 continue
+ do 7 i=1,5
+ ispot=1+nwide4*(i-1)
+ call move(agraph,ispot,aper,1,1)
+ 7 continue
+ locyt=locy
+ lspot=locv-1
+ mltscl=0
+ if (value(locv).eq.0.0d0) mltscl=1
+ do 235 k=1,kntr
+ lspot=lspot+2
+ ymin=value(lspot)
+ ymax=value(lspot+1)
+ if (ymin.ne.0.0d0) go to 10
+ if (ymax.ne.0.0d0) go to 10
+ go to 100
+ 10 ymin1=dmin1(ymin,ymax)
+ ymax1=dmax1(ymin,ymax)
+ 30 if (ilogy(k).eq.1) go to 40
+ ymin1=dlog10(dmax1(ymin1,1.0d-20))
+ ymax1=dlog10(dmax1(ymax1,1.0d-20))
+ del=dmax1(ymax1-ymin1,0.0001d0)/4.0d0
+ go to 50
+ 40 del=dmax1(ymax1-ymin1,1.0d-20)/4.0d0
+ 50 ymin=ymin1
+ ymax=ymax1
+ go to 200
+c
+c determine max and min values
+c
+ 100 ymax1=value(locyt+1)
+ ymin1=ymax1
+ if (numpnt.eq.1) go to 150
+ do 110 i=2,numpnt
+ ymin1=dmin1(ymin1,value(locyt+i))
+ ymax1=dmax1(ymax1,value(locyt+i))
+ 110 continue
+c
+c scaling
+c
+ 150 call scale(ymin1,ymax1,4,ymin,ymax,del)
+c
+c determine coordinates
+c
+ 200 ycoor(1,k)=ymin
+ pmin(k)=ymin
+ small=del*1.0d-4
+ if (dabs(ycoor(1,k)).le.small) ycoor(1,k)=0.0d0
+ do 210 i=1,4
+ ycoor(i+1,k)=ycoor(i,k)+del
+ if (dabs(ycoor(i+1,k)).le.small) ycoor(i+1,k)=0.0d0
+ 210 continue
+ if (ilogy(k).eq.1) go to 230
+ do 220 i=1,5
+ 220 ycoor(i,k)=dexp(xlog10*ycoor(i,k))
+ 230 delplt(k)=del/dfloat(nwide4)
+ locyt=locyt+npoint
+ 235 continue
+c
+c count distinct coordinates
+c
+ icoor(1)=1
+ jcoor(1)=1
+ numcor=1
+ if (kntr.eq.1) go to 290
+ do 250 i=2,kntr
+ do 245 j=1,numcor
+ l=jcoor(j)
+c... coordinates are *equal* if the most significant 24 bits agree
+ do 240 k=1,5
+ y1=ycoor(k,i)
+ y2=ycoor(k,l)
+ if(y1.eq.0.0d0.and.y2.eq.0.0d0) go to 240
+ if(dabs((y1-y2)/dmax1(dabs(y1),dabs(y2))).ge.1.0d-7) go to 245
+ 240 continue
+ icoor(i)=l
+ go to 250
+ 245 continue
+ icoor(i)=i
+ numcor=numcor+1
+ jcoor(numcor)=i
+ 250 continue
+c
+c print coordinates
+c
+ 260 do 280 i=1,numcor
+ asym(1)=asym1
+ asym(2)=asym2
+ ipos=2
+ do 270 j=1,kntr
+ if (icoor(j).ne.jcoor(i)) go to 270
+ call move(asym,ipos,pltsym,j,1)
+ ipos=ipos+1
+ 270 continue
+ call move(asym,ipos,arprn,1,1)
+ k=jcoor(i)
+ if(iwide.ne.0) write(6,271) asym,(ycoor(j,k),j=1,5)
+ 271 format(/1hx,2a8,4h----,1pd12.3,4(15x,d10.3)/26x,51(2h -))
+ if(iwide.eq.0) write(6,273) asym,(ycoor(j,k),j=1,5)
+ 273 format(/1hx,2a8,1pd10.3,3(4x,d10.3),1x,d10.3/22x,29(2h -))
+ 280 continue
+ go to 300
+ 290 if(iwide.ne.0) write(6,291) (ycoor(j,1),j=1,5)
+ 291 format(/1hx,20x,1pd12.3,4(15x,d10.3)/26x,51(2h -))
+ if(iwide.eq.0) write(6,293) (ycoor(j,1),j=1,5)
+ 293 format(/1hx,14x,1pd12.3,3(4x,d10.3),1x,d10.3/22x,29(2h -))
+c
+c plotting
+c
+ 300 aspot=ablnk
+ do 320 i=1,numpnt
+ xvar=value(locx+i)
+ locyt=locy
+ call copy8(agraph,aplot,13)
+ do 310 k=1,kntr
+ yvr=value(locyt+i)
+ ktmp=icoor(k)
+ ymin1=pmin(ktmp)
+ jpoint=idint((yvr-ymin1)/delplt(k)+0.5d0)+1
+ if (jpoint.le.0) go to 306
+ if (jpoint.gt.nwide) go to 306
+ call move(aspot,1,aplot,jpoint,1)
+ if (aspot.eq.ablnk) go to 303
+ if (aspot.eq.aper) go to 303
+ call move(aplot,jpoint,aletx,1,1)
+ go to 306
+ 303 call move(aplot,jpoint,pltsym,k,1)
+ 306 locyt=locyt+npoint
+ 310 continue
+ yvr=value(locy+i)
+ if (ilogy(1).eq.1) go to 315
+ yvr=dexp(xlog10*yvr)
+ 315 if(iwide.ne.0) write(6,316) xvar,yvr,aplot
+ 316 format(1x,1pd10.3,2x,d10.3,2x,13a8)
+ if(iwide.eq.0) write(6,317) xvar,yvr,(aplot(k),k=1,8)
+ 317 format(1x,1pd10.3,1x,d10.3,7a8,a1)
+ 320 continue
+c
+c finished
+c
+ if(iwide.ne.0) write(6,331)
+ 331 format(26x,51(2h -)//)
+ if(iwide.eq.0) write(6,332)
+ 332 format(21x,29(2h -)//)
+ go to 500
+c
+c too few points
+c
+ 400 write (6,401)
+ 401 format('0warning: too few points for plotting'/)
+ 500 write (6,501)
+ 501 format(1hy)
+ return
+ end
+ subroutine scale(xmin,xmax,n,xminp,xmaxp,del)
+ implicit double precision (a-h,o-z)
+c
+c this routine determines the 'optimal' scale to use for the plot of
+c some output variable.
+c
+c
+c adapted from algorithm 463 of 'collected algorithms of the cacm'
+c
+ common /knstnt/ twopi,xlog2,xlog10,root2,rad,boltz,charge,ctok,
+ 1 gmin,reltol,abstol,vntol,trtol,chgtol,eps0,epssil,epsox
+ integer xxor
+ dimension vint(5)
+ data vint / 1.0d0,2.0d0,5.0d0,10.0d0,20.0d0 /
+ data eps / 1.0d-12 /
+c
+c
+c... trap too-small data spread
+ if(xmin.eq.0.0d0.and.xmax.eq.0.0d0) go to 4
+ if(dabs((xmax-xmin)/dmax1(dabs(xmin),dabs(xmax))).ge.1.0d-4)
+ 1 go to 10
+ 4 continue
+ if (xmin.ge.0.0d0) go to 5
+ xmax=0.5d0*xmin+eps
+ xmin=1.5d0*xmin-eps
+ go to 10
+ 5 xmax=1.5d0*xmin+eps
+ xmin=0.5d0*xmin-eps
+c... find approximate interval size, normalized to [1,10]
+ 10 a=(xmax-xmin)/dfloat(n)
+ nal=idint(dlog10(a))
+ if (a.lt.1.0d0) nal=nal-1
+ xfact=dexp(xlog10*dfloat(nal))
+ b=a/xfact
+c... find closest permissible interval size
+ do 20 i=1,3
+ if (b.lt.(vint(i)+eps)) go to 30
+ 20 continue
+ i=4
+c... compute interval size
+ 30 del=vint(i)*xfact
+ fm1=xmin/del
+ m1=fm1
+ if (fm1.lt.0.0d0) m1=m1-1
+ if (dabs(dfloat(m1)+1.0d0-fm1).lt.eps) m1=m1+1
+c... compute new maximum and minimum limits
+ xminp=del*dfloat(m1)
+ fm2=xmax/del
+ m2=fm2+1.0d0
+ if (fm2.lt.(-1.0d0)) m2=m2-1
+ if (dabs(fm2+1.0d0-dfloat(m2)).lt.eps) m2=m2-1
+ xmaxp=del*dfloat(m2)
+ np=m2-m1
+c... check whether another loop required
+ if (np.le.n) go to 40
+ i=i+1
+ go to 30
+c... do final adjustments and correct for roundoff error(s)
+ 40 nx=(n-np)/2
+ xminp=dmin1(xmin,xminp-dfloat(nx)*del)
+ xmaxp=dmax1(xmax,xminp+dfloat(n)*del)
+ return
+ end
+ subroutine fouran
+ implicit double precision (a-h,o-z)
+c
+c this routine determines the fourier coefficients of a transient
+c analysis waveform.
+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 /miscel/ atime,aprog(3),adate,atitle(10),defl,defw,defad,
+ 1 defas,rstats(50),iwidth,lwidth,nopage
+ 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 /tran/ tstep,tstop,tstart,delmax,tdmax,forfre,jtrflg
+ common /outinf/ xincr,string(15),xstart,yvar(8),itab(8),itype(8),
+ 1 ilogy(8),npoint,numout,kntr,numdgt
+ common /blank/ value(1000)
+ integer nodplc(64)
+ complex*16 cvalue(32)
+ equivalence (value(1),nodplc(1),cvalue(1))
+c
+c
+ dimension sinco(9),cosco(9)
+ dimension fortit(4)
+ data fortit / 8hfourier , 8hanalysis, 8h , 8h /
+ data ablnk / 1h /
+c
+c
+ forprd=1.0d0/forfre
+ xstart=tstop-forprd
+ kntr=1
+ xn=101.0d0
+ xincr=forprd/xn
+ npoint=xn
+ call getm8(locx,npoint)
+ call getm8(locy,npoint)
+ do 105 nknt=1,nfour
+ itab(1)=nodplc(ifour+nknt)
+ kfrout=itab(1)
+ call ntrpl8(locx,locy,numpnt)
+ dcco=0.0d0
+ call zero8(sinco,9)
+ call zero8(cosco,9)
+ loct=locy+1
+ ipnt=0
+ 10 yvr=value(loct+ipnt)
+ dcco=dcco+yvr
+ forfac=dfloat(ipnt)*twopi/xn
+ arg=0.0d0
+ do 20 k=1,9
+ arg=arg+forfac
+ sinco(k)=sinco(k)+yvr*dsin(arg)
+ cosco(k)=cosco(k)+yvr*dcos(arg)
+ 20 continue
+ ipnt=ipnt+1
+ if (ipnt.ne.npoint) go to 10
+ dcco=dcco/xn
+ forfac=2.0d0/xn
+ do 30 k=1,9
+ sinco(k)=sinco(k)*forfac
+ cosco(k)=cosco(k)*forfac
+ 30 continue
+ call title(0,72,1,fortit)
+ ipos=1
+ call outnam(kfrout,1,string,ipos)
+ call move(string,ipos,ablnk,1,7)
+ jstop=(ipos+6)/8
+ write (6,61) (string(j),j=1,jstop)
+ 61 format(' fourier components of transient response ',5a8///)
+ write (6,71) dcco
+ 71 format('0dc component =',1pd12.3/,
+ 1 '0harmonic frequency fourier normalized phase no
+ 2rmalized'/,
+ 3 ' no (hz) component component (deg) pha
+ 4se (deg)'//)
+ iknt=1
+ freq1=forfre
+ xnharm=1.0d0
+ call magphs(dcmplx(sinco(1),cosco(1)),xnorm,pnorm)
+ phasen=0.0d0
+ write (6,81) iknt,freq1,xnorm,xnharm,pnorm,phasen
+ 81 format(i6,1pd15.3,d12.3,0pf13.6,f10.3,f12.3/)
+ thd=0.0d0
+ do 90 iknt=2,9
+ freq1=dfloat(iknt)*forfre
+ call magphs(dcmplx(sinco(iknt),cosco(iknt)),
+ 1 harm,phase)
+ xnharm=harm/xnorm
+ phasen=phase-pnorm
+ thd=thd+xnharm*xnharm
+ write (6,81) iknt,freq1,harm,xnharm,phase,phasen
+ 90 continue
+ thd=100.0d0*dsqrt(thd)
+ write (6,101) thd
+ 101 format (//5x,'total harmonic distortion = ',f12.6,' percent')
+ 105 continue
+ call clrmem(locx)
+ call clrmem(locy)
+ 110 return
+ end
projects / unix-history / commitdiff