projects / unix-history / blame
| Commit | Line | Data |
|---|---|---|
| e0fe7396 D |
1 | subroutine ovtpvt |
| 2 | implicit double precision (a-h,o-z) | |
| 3 | c | |
| 4 | c | |
| 5 | c this routine generates the requested tabular listings of analysis | |
| 6 | c results. it calls plot to generate line-printer plots. | |
| 7 | c | |
| 8 | common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem, | |
| 9 | 1 isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize, | |
| 10 | 2 junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr, | |
| 11 | 3 nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1, | |
| 12 | 4 lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd, | |
| 13 | 5 imynl,imvn,lcvn,loutpt,nsnod,nsmat,nsval,icnod,icmat,icval | |
| 14 | common /cirdat/ locate(50),jelcnt(50),nunods,ncnods,numnod,nstop, | |
| 15 | 1 nut,nlt,nxtrm,ndist,ntlin,ibr,numvs | |
| 16 | common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet, | |
| 17 | 1 xmu,mode,modedc,icalc,initf,method,iord,maxord,noncon,iterno, | |
| 18 | 2 itemno,nosolv,ipostp,iscrch | |
| 19 | common /flags/ iprnta,iprntl,iprntm,iprntn,iprnto,limtim,limpts, | |
| 20 | 1 lvlcod,lvltim,itl1,itl2,itl3,itl4,itl5,igoof,nogo,keof | |
| 21 | common /miscel/ atime,aprog(3),adate,atitle(10),defl,defw,defad, | |
| 22 | 1 defas,rstats(50),iwidth,lwidth,nopage | |
| 23 | common /dc/ tcstar(2),tcstop(2),tcincr(2),icvflg,itcelm(2),kssop, | |
| 24 | 1 kinel,kidin,kovar,kidout | |
| 25 | common /ac/ fstart,fstop,fincr,skw2,refprl,spw2,jacflg,idfreq, | |
| 26 | 1 inoise,nosprt,nosout,nosin,idist,idprt | |
| 27 | common /tran/ tstep,tstop,tstart,delmax,tdmax,forfre,jtrflg | |
| 28 | common /outinf/ xincr,string(15),xstart,yvar(8),itab(8),itype(8), | |
| 29 | 1 ilogy(8),npoint,numout,kntr,numdgt | |
| 30 | common /blank/ value(1000) | |
| 31 | integer nodplc(64) | |
| 32 | complex*16 cvalue(32) | |
| 33 | equivalence (value(1),nodplc(1),cvalue(1)) | |
| 34 | c | |
| 35 | complex*16 cval | |
| 36 | dimension prform(3) | |
| 37 | dimension subtit(4,3) | |
| 38 | data subtit / 8hdc trans, 8hfer curv, 8hes , 8h , | |
| 39 | 1 8htransien, 8ht analys, 8his , 8h , | |
| 40 | 2 8hac analy, 8hsis , 8h , 8h / | |
| 41 | data prform / 8h(1pe11.3, 8h,2x,8e00, 8h.00) / | |
| 42 | data aper,rprn / 1h., 1h) / | |
| 43 | c | |
| 44 | call second(t1) | |
| 45 | if (icalc.le.0) go to 1000 | |
| 46 | call crunch | |
| 47 | if (nogo.lt.0) go to 1000 | |
| 48 | c | |
| 49 | c construct format statement to be used for printing the outputs | |
| 50 | c | |
| 51 | ifract=max0(numdgt-1,0) | |
| 52 | ifwdth=ifract+9 | |
| 53 | ipos=15 | |
| 54 | call alfnum(ifwdth,prform,ipos) | |
| 55 | call move(prform,ipos,aper,1,1) | |
| 56 | ipos=ipos+1 | |
| 57 | call alfnum(ifract,prform,ipos) | |
| 58 | call move(prform,ipos,rprn,1,1) | |
| 59 | c | |
| 60 | noprln=min0(8,(lwidth-12)/ifwdth) | |
| 61 | if (mode-2) 50,60,300 | |
| 62 | 50 numout=jelcnt(41)+1 | |
| 63 | go to 70 | |
| 64 | 60 numout=jelcnt(42)+1 | |
| 65 | c | |
| 66 | c dc and transient analysis printing | |
| 67 | c | |
| 68 | 70 loc=locate(30+mode) | |
| 69 | 80 if (loc.eq.0) go to 200 | |
| 70 | kntr=min0(noprln,nodplc(loc+3)) | |
| 71 | if (kntr.le.0) go to 120 | |
| 72 | call title(1,lwidth,1,subtit(1,mode)) | |
| 73 | call setprn(loc) | |
| 74 | c | |
| 75 | c get buffer space | |
| 76 | c | |
| 77 | call getm8(locx,npoint) | |
| 78 | call getm8(locy,kntr*npoint) | |
| 79 | c | |
| 80 | c interpolate outputs | |
| 81 | c | |
| 82 | call ntrpl8(locx,locy,numpnt) | |
| 83 | c | |
| 84 | c print outputs | |
| 85 | c | |
| 86 | do 100 i=1,numpnt | |
| 87 | xvar=value(locx+i) | |
| 88 | locyt=locy | |
| 89 | do 90 k=1,kntr | |
| 90 | yvar(k)=value(locyt+i) | |
| 91 | locyt=locyt+npoint | |
| 92 | 90 continue | |
| 93 | write (6,prform) xvar,(yvar(k),k=1,kntr) | |
| 94 | 100 continue | |
| 95 | write (6,111) | |
| 96 | 111 format(1hy) | |
| 97 | call clrmem(locx) | |
| 98 | call clrmem(locy) | |
| 99 | 120 loc=nodplc(loc) | |
| 100 | go to 80 | |
| 101 | c | |
| 102 | c dc and transient analysis plotting | |
| 103 | c | |
| 104 | 200 loc=locate(35+mode) | |
| 105 | 210 if (loc.eq.0) go to 250 | |
| 106 | kntr=nodplc(loc+3) | |
| 107 | if (kntr.le.0) go to 220 | |
| 108 | locv=nodplc(loc+1) | |
| 109 | call title(1,lwidth,1,subtit(1,mode)) | |
| 110 | call setplt(loc) | |
| 111 | c | |
| 112 | c get buffer space | |
| 113 | c | |
| 114 | call getm8(locx,npoint) | |
| 115 | call getm8(locy,kntr*npoint) | |
| 116 | c | |
| 117 | c interpolate outputs and load plot buffers | |
| 118 | c | |
| 119 | call ntrpl8(locx,locy,numpnt) | |
| 120 | call plot(numpnt,locx,locy,locv) | |
| 121 | call clrmem(locx) | |
| 122 | call clrmem(locy) | |
| 123 | 220 loc=nodplc(loc) | |
| 124 | go to 210 | |
| 125 | c | |
| 126 | c fourier analysis | |
| 127 | c | |
| 128 | 250 if (mode.eq.1) go to 1000 | |
| 129 | if (nfour.eq.0) go to 1000 | |
| 130 | if (nogo.ne.0) go to 1000 | |
| 131 | call fouran | |
| 132 | go to 1000 | |
| 133 | c | |
| 134 | c ac analysis printing | |
| 135 | c | |
| 136 | 300 numout=jelcnt(43)+jelcnt(44)+jelcnt(45)+1 | |
| 137 | do 599 id=33,35 | |
| 138 | loc=locate(id) | |
| 139 | 320 if (loc.eq.0) go to 599 | |
| 140 | kntr=min0(noprln,nodplc(loc+3)) | |
| 141 | if (kntr.le.0) go to 595 | |
| 142 | call title(1,lwidth,1,subtit(1,mode)) | |
| 143 | call setprn(loc) | |
| 144 | c | |
| 145 | c print ac outputs | |
| 146 | c | |
| 147 | lout=loutpt | |
| 148 | do 590 i=1,icalc | |
| 149 | xvar=dreal(cvalue(lout+1)) | |
| 150 | do 500 k=1,kntr | |
| 151 | iseq=itab(k) | |
| 152 | iseq=nodplc(iseq+4) | |
| 153 | cval=cvalue(lout+iseq) | |
| 154 | ktype=itype(k) | |
| 155 | go to (450,450,430,440,450,450), ktype | |
| 156 | 430 yvar(k)=dreal(cval) | |
| 157 | go to 500 | |
| 158 | 440 yvar(k)=dimag(cval) | |
| 159 | go to 500 | |
| 160 | 450 call magphs(cval,xmag,xphs) | |
| 161 | go to (460,460,430,440,470,465), ktype | |
| 162 | 460 yvar(k)=xmag | |
| 163 | go to 500 | |
| 164 | 465 yvar(k)=20.0d0*dlog10(xmag) | |
| 165 | go to 500 | |
| 166 | 470 yvar(k)=xphs | |
| 167 | 500 continue | |
| 168 | lout=lout+numout | |
| 169 | 580 write (6,prform) xvar,(yvar(k),k=1,kntr) | |
| 170 | 590 continue | |
| 171 | write (6,111) | |
| 172 | 595 loc=nodplc(loc) | |
| 173 | go to 320 | |
| 174 | 599 continue | |
| 175 | c | |
| 176 | c ac analysis plotting | |
| 177 | c | |
| 178 | do 760 id=38,40 | |
| 179 | loc=locate(id) | |
| 180 | 610 if (loc.eq.0) go to 760 | |
| 181 | kntr=nodplc(loc+3) | |
| 182 | if (kntr.le.0) go to 750 | |
| 183 | locv=nodplc(loc+1) | |
| 184 | call title(1,lwidth,1,subtit(1,mode)) | |
| 185 | call setplt(loc) | |
| 186 | c | |
| 187 | call getm8(locx,icalc) | |
| 188 | call getm8(locy,kntr*icalc) | |
| 189 | c | |
| 190 | c load plot buffers | |
| 191 | c | |
| 192 | lout=loutpt | |
| 193 | do 710 i=1,icalc | |
| 194 | xvar=dreal(cvalue(lout+1)) | |
| 195 | locyt=locy | |
| 196 | do 700 k=1,kntr | |
| 197 | iseq=itab(k) | |
| 198 | iseq=nodplc(iseq+4) | |
| 199 | cval=cvalue(lout+iseq) | |
| 200 | ktype=itype(k) | |
| 201 | go to (670,670,650,660,670,670), ktype | |
| 202 | 650 yvr=dreal(cval) | |
| 203 | go to 695 | |
| 204 | 660 yvr=dimag(cval) | |
| 205 | go to 695 | |
| 206 | 670 call magphs(cval,xmag,xphs) | |
| 207 | go to (680,680,650,660,690,685), ktype | |
| 208 | 680 yvr=dlog10(xmag) | |
| 209 | go to 695 | |
| 210 | 685 yvr=20.0d0*dlog10(xmag) | |
| 211 | go to 695 | |
| 212 | 690 yvr=xphs | |
| 213 | 695 value(locyt+i)=yvr | |
| 214 | locyt=locyt+icalc | |
| 215 | 700 continue | |
| 216 | value(locx+i)=xvar | |
| 217 | lout=lout+numout | |
| 218 | 710 continue | |
| 219 | call plot(icalc,locx,locy,locv) | |
| 220 | call clrmem(locx) | |
| 221 | call clrmem(locy) | |
| 222 | 750 loc=nodplc(loc) | |
| 223 | go to 610 | |
| 224 | 760 continue | |
| 225 | c | |
| 226 | c finished | |
| 227 | c | |
| 228 | 1000 call clrmem(loutpt) | |
| 229 | call second(t2) | |
| 230 | rstats(11)=rstats(11)+t2-t1 | |
| 231 | return | |
| 232 | end | |
| 233 | subroutine ntrpl8(locx,locy,numpnt) | |
| 234 | implicit double precision (a-h,o-z) | |
| 235 | c | |
| 236 | c this routine interpolates the analysis data to obtain the values | |
| 237 | c printed and/or plotted, using linear interpolation. | |
| 238 | c | |
| 239 | common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem, | |
| 240 | 1 isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize, | |
| 241 | 2 junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr, | |
| 242 | 3 nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1, | |
| 243 | 4 lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd, | |
| 244 | 5 imynl,imvn,lcvn,loutpt,nsnod,nsmat,nsval,icnod,icmat,icval | |
| 245 | common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet, | |
| 246 | 1 xmu,mode,modedc,icalc,initf,method,iord,maxord,noncon,iterno, | |
| 247 | 2 itemno,nosolv,ipostp,iscrch | |
| 248 | common /outinf/ xincr,string(15),xstart,yvar(8),itab(8),itype(8), | |
| 249 | 1 ilogy(8),npoint,numout,kntr,numdgt | |
| 250 | common /blank/ value(1000) | |
| 251 | integer nodplc(64) | |
| 252 | complex*16 cvalue(32) | |
| 253 | equivalence (value(1),nodplc(1),cvalue(1)) | |
| 254 | if(mode.ne.1) go to 4 | |
| 255 | numpnt=icalc | |
| 256 | loco=loutpt | |
| 257 | do 3 i=1,numpnt | |
| 258 | locyt=locy | |
| 259 | value(locx+i)=value(loco+1) | |
| 260 | do 2 k=1,kntr | |
| 261 | iseq=itab(k) | |
| 262 | iseq=nodplc(iseq+4) | |
| 263 | value(locyt+i)=value(loco+iseq) | |
| 264 | locyt=locyt+npoint | |
| 265 | 2 continue | |
| 266 | loco=loco+numout | |
| 267 | 3 continue | |
| 268 | return | |
| 269 | 4 continue | |
| 270 | xvar=xstart | |
| 271 | xvtol=xincr*1.0d-5 | |
| 272 | ippnt=0 | |
| 273 | icpnt=2 | |
| 274 | loco1=loutpt | |
| 275 | loco2=loco1+numout | |
| 276 | if (icalc.lt.2) go to 50 | |
| 277 | 10 x1=value(loco1+1) | |
| 278 | x2=value(loco2+1) | |
| 279 | dx1x2=x1-x2 | |
| 280 | 20 if (xincr.lt.0.0d0) go to 24 | |
| 281 | if (xvar.le.(x2+xvtol)) go to 30 | |
| 282 | go to 28 | |
| 283 | 24 if (xvar.ge.(x2+xvtol)) go to 30 | |
| 284 | 28 if (icpnt.ge.icalc) go to 100 | |
| 285 | icpnt=icpnt+1 | |
| 286 | loco1=loco2 | |
| 287 | loco2=loco1+numout | |
| 288 | go to 10 | |
| 289 | 30 ippnt=ippnt+1 | |
| 290 | value(locx+ippnt)=xvar | |
| 291 | dxx1=xvar-x1 | |
| 292 | locyt=locy | |
| 293 | do 40 i=1,kntr | |
| 294 | iseq=itab(i) | |
| 295 | iseq=nodplc(iseq+4) | |
| 296 | v1=value(loco1+iseq) | |
| 297 | v2=value(loco2+iseq) | |
| 298 | yvr=v1+(v1-v2)*dxx1/dx1x2 | |
| 299 | tol=dmin1(dabs(v1),dabs(v2))*1.0d-10 | |
| 300 | if (dabs(yvr).le.tol) yvr=0.0d0 | |
| 301 | value(locyt+ippnt)=yvr | |
| 302 | locyt=locyt+npoint | |
| 303 | 40 continue | |
| 304 | if (ippnt.ge.npoint) go to 100 | |
| 305 | xvar=xstart+dfloat(ippnt)*xincr | |
| 306 | if (dabs(xvar).ge.dabs(xvtol)) go to 20 | |
| 307 | xvar=0.0d0 | |
| 308 | go to 20 | |
| 309 | c | |
| 310 | c special handling if icalc = 1 | |
| 311 | c | |
| 312 | c... icalc=1; just copy over the single point and return | |
| 313 | 50 ippnt=1 | |
| 314 | value(locx+ippnt)=xvar | |
| 315 | locyt=locy | |
| 316 | do 60 i=1,kntr | |
| 317 | iseq=itab(i) | |
| 318 | iseq=nodplc(iseq+4) | |
| 319 | value(locyt+ippnt)=value(loco1+iseq) | |
| 320 | locyt=locyt+npoint | |
| 321 | 60 continue | |
| 322 | go to 100 | |
| 323 | c | |
| 324 | c return | |
| 325 | c | |
| 326 | 100 numpnt=ippnt | |
| 327 | return | |
| 328 | end | |
| 329 | subroutine setprn(loc) | |
| 330 | implicit double precision (a-h,o-z) | |
| 331 | c | |
| 332 | c this routine formats the column headers for tabular listings of | |
| 333 | c output variables. | |
| 334 | c | |
| 335 | common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem, | |
| 336 | 1 isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize, | |
| 337 | 2 junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr, | |
| 338 | 3 nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1, | |
| 339 | 4 lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd, | |
| 340 | 5 imynl,imvn,lcvn,loutpt,nsnod,nsmat,nsval,icnod,icmat,icval | |
| 341 | common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet, | |
| 342 | 1 xmu,mode,modedc,icalc,initf,method,iord,maxord,noncon,iterno, | |
| 343 | 2 itemno,nosolv,ipostp,iscrch | |
| 344 | common /miscel/ atime,aprog(3),adate,atitle(10),defl,defw,defad, | |
| 345 | 1 defas,rstats(50),iwidth,lwidth,nopage | |
| 346 | common /dc/ tcstar(2),tcstop(2),tcincr(2),icvflg,itcelm(2),kssop, | |
| 347 | 1 kinel,kidin,kovar,kidout | |
| 348 | common /ac/ fstart,fstop,fincr,skw2,refprl,spw2,jacflg,idfreq, | |
| 349 | 1 inoise,nosprt,nosout,nosin,idist,idprt | |
| 350 | common /tran/ tstep,tstop,tstart,delmax,tdmax,forfre,jtrflg | |
| 351 | common /outinf/ xincr,string(15),xstart,yvar(8),itab(8),itype(8), | |
| 352 | 1 ilogy(8),npoint,numout,kntr,numdgt | |
| 353 | common /blank/ value(1000) | |
| 354 | integer nodplc(64) | |
| 355 | complex*16 cvalue(32) | |
| 356 | equivalence (value(1),nodplc(1),cvalue(1)) | |
| 357 | c | |
| 358 | data ablnk, atimex, afreq / 1h , 6h time, 6h freq / | |
| 359 | c | |
| 360 | c set limits depending upon the analysis mode | |
| 361 | c | |
| 362 | if (mode-2) 10,20,30 | |
| 363 | 10 xstart=tcstar(1) | |
| 364 | xincr=tcincr(1) | |
| 365 | npoint=icvflg | |
| 366 | itemp=itcelm(1) | |
| 367 | loce=nodplc(itemp+1) | |
| 368 | asweep=value(loce) | |
| 369 | go to 40 | |
| 370 | 20 xstart=tstart | |
| 371 | xincr=tstep | |
| 372 | npoint=jtrflg | |
| 373 | asweep=atimex | |
| 374 | go to 40 | |
| 375 | 30 xstart=fstart | |
| 376 | xincr=fincr | |
| 377 | npoint=icalc | |
| 378 | asweep=afreq | |
| 379 | c | |
| 380 | c construct and print the output variable names | |
| 381 | c | |
| 382 | 40 loct=loc+2 | |
| 383 | ipos=1 | |
| 384 | npos=ipos+numdgt+8 | |
| 385 | do 90 i=1,kntr | |
| 386 | loct=loct+2 | |
| 387 | itab(i)=nodplc(loct) | |
| 388 | itype(i)=nodplc(loct+1) | |
| 389 | call outnam(itab(i),itype(i),string,ipos) | |
| 390 | if (ipos.ge.npos) go to 70 | |
| 391 | do 60 j=ipos,npos | |
| 392 | call move(string,j,ablnk,1,1) | |
| 393 | 60 continue | |
| 394 | ipos=npos | |
| 395 | go to 80 | |
| 396 | 70 call move(string,ipos,ablnk,1,1) | |
| 397 | ipos=ipos+1 | |
| 398 | 80 npos=npos+numdgt+8 | |
| 399 | 90 continue | |
| 400 | call move(string,ipos,ablnk,1,7) | |
| 401 | jstop=(ipos+6)/8 | |
| 402 | write (6,91) asweep,(string(j),j=1,jstop) | |
| 403 | 91 format(/3x,a8,5x,14a8,a4) | |
| 404 | write (6,101) | |
| 405 | 101 format(1hx/1h ) | |
| 406 | return | |
| 407 | end | |
| 408 | subroutine setplt(loc) | |
| 409 | implicit double precision (a-h,o-z) | |
| 410 | c | |
| 411 | c this routine generates the 'legend' subheading used to identify | |
| 412 | c individual traces on multi-trace line-printer plots. | |
| 413 | c | |
| 414 | common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem, | |
| 415 | 1 isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize, | |
| 416 | 2 junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr, | |
| 417 | 3 nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1, | |
| 418 | 4 lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd, | |
| 419 | 5 imynl,imvn,lcvn,loutpt,nsnod,nsmat,nsval,icnod,icmat,icval | |
| 420 | common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet, | |
| 421 | 1 xmu,mode,modedc,icalc,initf,method,iord,maxord,noncon,iterno, | |
| 422 | 2 itemno,nosolv,ipostp,iscrch | |
| 423 | common /miscel/ atime,aprog(3),adate,atitle(10),defl,defw,defad, | |
| 424 | 1 defas,rstats(50),iwidth,lwidth,nopage | |
| 425 | common /dc/ tcstar(2),tcstop(2),tcincr(2),icvflg,itcelm(2),kssop, | |
| 426 | 1 kinel,kidin,kovar,kidout | |
| 427 | common /ac/ fstart,fstop,fincr,skw2,refprl,spw2,jacflg,idfreq, | |
| 428 | 1 inoise,nosprt,nosout,nosin,idist,idprt | |
| 429 | common /tran/ tstep,tstop,tstart,delmax,tdmax,forfre,jtrflg | |
| 430 | common /outinf/ xincr,string(15),xstart,yvar(8),itab(8),itype(8), | |
| 431 | 1 ilogy(8),npoint,numout,kntr,numdgt | |
| 432 | common /blank/ value(1000) | |
| 433 | integer nodplc(64) | |
| 434 | complex*16 cvalue(32) | |
| 435 | equivalence (value(1),nodplc(1),cvalue(1)) | |
| 436 | c | |
| 437 | dimension logopt(6) | |
| 438 | data logopt / 2, 2, 1, 1, 1, 1 / | |
| 439 | data ablnk, atimex, afreq / 1h , 6h time, 6h freq / | |
| 440 | data pltsym / 8h*+=$0<>? / | |
| 441 | c | |
| 442 | c set limits depending upon the analysis mode | |
| 443 | c | |
| 444 | if (mode-2) 10,20,30 | |
| 445 | 10 xstart=tcstar(1) | |
| 446 | xincr=tcincr(1) | |
| 447 | npoint=icvflg | |
| 448 | itemp=itcelm(1) | |
| 449 | loce=nodplc(itemp+1) | |
| 450 | asweep=value(loce) | |
| 451 | go to 40 | |
| 452 | 20 xstart=tstart | |
| 453 | xincr=tstep | |
| 454 | npoint=jtrflg | |
| 455 | asweep=atimex | |
| 456 | go to 40 | |
| 457 | 30 xstart=fstart | |
| 458 | xincr=fincr | |
| 459 | npoint=jacflg | |
| 460 | asweep=afreq | |
| 461 | c | |
| 462 | c construct and print the output variables with corresponding plot | |
| 463 | c symbols | |
| 464 | c | |
| 465 | 40 loct=loc+2 | |
| 466 | if (kntr.eq.1) go to 80 | |
| 467 | write (6,41) | |
| 468 | 41 format('0legend:'/) | |
| 469 | do 70 i=1,kntr | |
| 470 | loct=loct+2 | |
| 471 | itab(i)=nodplc(loct) | |
| 472 | ioutyp=nodplc(loct+1) | |
| 473 | itype(i)=ioutyp | |
| 474 | ilogy(i)=1 | |
| 475 | if (mode.le.2) go to 50 | |
| 476 | ilogy(i)=logopt(ioutyp) | |
| 477 | 50 ipos=1 | |
| 478 | call outnam(itab(i),itype(i),string,ipos) | |
| 479 | call move(string,ipos,ablnk,1,7) | |
| 480 | jstop=(ipos+6)/8 | |
| 481 | call move(achar,1,pltsym,i,1) | |
| 482 | write (6,61) achar,(string(j),j=1,jstop) | |
| 483 | 61 format(1x,a1,2h: ,5a8) | |
| 484 | 70 continue | |
| 485 | 80 if (kntr.ge.2) go to 90 | |
| 486 | itab(1)=nodplc(loc+4) | |
| 487 | ioutyp=nodplc(loc+5) | |
| 488 | itype(1)=ioutyp | |
| 489 | ilogy(1)=1 | |
| 490 | if (mode.le.2) go to 90 | |
| 491 | ilogy(1)=logopt(ioutyp) | |
| 492 | 90 ipos=1 | |
| 493 | call outnam(itab(1),itype(1),string,ipos) | |
| 494 | call move(string,ipos,ablnk,1,7) | |
| 495 | jstop=(ipos+6)/8 | |
| 496 | write (6,101) asweep,(string(j),j=1,jstop) | |
| 497 | 101 format(1hx/3x,a8,4x,5a8) | |
| 498 | return | |
| 499 | end | |
| 500 | subroutine plot(numpnt,locx,locy,locv) | |
| 501 | implicit double precision (a-h,o-z) | |
| 502 | c | |
| 503 | c this routine generates the line-printer plots. | |
| 504 | c | |
| 505 | common /miscel/ atime,aprog(3),adate,atitle(10),defl,defw,defad, | |
| 506 | 1 defas,rstats(50),iwidth,lwidth,nopage | |
| 507 | common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet, | |
| 508 | 1 xmu,mode,modedc,icalc,initf,method,iord,maxord,noncon,iterno, | |
| 509 | 2 itemno,nosolv,ipostp,iscrch | |
| 510 | common /knstnt/ twopi,xlog2,xlog10,root2,rad,boltz,charge,ctok, | |
| 511 | 1 gmin,reltol,abstol,vntol,trtol,chgtol,eps0,epssil,epsox | |
| 512 | common /outinf/ xincr,string(15),xstart,yvar(8),itab(8),itype(8), | |
| 513 | 1 ilogy(8),npoint,numout,kntr,numdgt | |
| 514 | common /blank/ value(1000) | |
| 515 | integer nodplc(64) | |
| 516 | complex*16 cvalue(32) | |
| 517 | equivalence (value(1),nodplc(1),cvalue(1)) | |
| 518 | c | |
| 519 | c | |
| 520 | integer xxor | |
| 521 | dimension ycoor(5,8),icoor(8),delplt(8) | |
| 522 | dimension agraph(13),aplot(13) | |
| 523 | dimension asym(2),pmin(8),jcoor(8) | |
| 524 | data ablnk, aletx, aper / 1h , 1hx, 1h. / | |
| 525 | data asym1, asym2, arprn / 8h(-------, 8h--------, 1h) / | |
| 526 | data pltsym / 8h*+=$0<>? / | |
| 527 | c | |
| 528 | c | |
| 529 | iwide=1 | |
| 530 | nwide=101 | |
| 531 | nwide4=25 | |
| 532 | if(lwidth.gt.80) go to 3 | |
| 533 | iwide=0 | |
| 534 | nwide=57 | |
| 535 | nwide4=14 | |
| 536 | 3 if (numpnt.le.0) go to 400 | |
| 537 | do 5 i=1,13 | |
| 538 | agraph(i)=ablnk | |
| 539 | 5 continue | |
| 540 | do 7 i=1,5 | |
| 541 | ispot=1+nwide4*(i-1) | |
| 542 | call move(agraph,ispot,aper,1,1) | |
| 543 | 7 continue | |
| 544 | locyt=locy | |
| 545 | lspot=locv-1 | |
| 546 | mltscl=0 | |
| 547 | if (value(locv).eq.0.0d0) mltscl=1 | |
| 548 | do 235 k=1,kntr | |
| 549 | lspot=lspot+2 | |
| 550 | ymin=value(lspot) | |
| 551 | ymax=value(lspot+1) | |
| 552 | if (ymin.ne.0.0d0) go to 10 | |
| 553 | if (ymax.ne.0.0d0) go to 10 | |
| 554 | go to 100 | |
| 555 | 10 ymin1=dmin1(ymin,ymax) | |
| 556 | ymax1=dmax1(ymin,ymax) | |
| 557 | 30 if (ilogy(k).eq.1) go to 40 | |
| 558 | ymin1=dlog10(dmax1(ymin1,1.0d-20)) | |
| 559 | ymax1=dlog10(dmax1(ymax1,1.0d-20)) | |
| 560 | del=dmax1(ymax1-ymin1,0.0001d0)/4.0d0 | |
| 561 | go to 50 | |
| 562 | 40 del=dmax1(ymax1-ymin1,1.0d-20)/4.0d0 | |
| 563 | 50 ymin=ymin1 | |
| 564 | ymax=ymax1 | |
| 565 | go to 200 | |
| 566 | c | |
| 567 | c determine max and min values | |
| 568 | c | |
| 569 | 100 ymax1=value(locyt+1) | |
| 570 | ymin1=ymax1 | |
| 571 | if (numpnt.eq.1) go to 150 | |
| 572 | do 110 i=2,numpnt | |
| 573 | ymin1=dmin1(ymin1,value(locyt+i)) | |
| 574 | ymax1=dmax1(ymax1,value(locyt+i)) | |
| 575 | 110 continue | |
| 576 | c | |
| 577 | c scaling | |
| 578 | c | |
| 579 | 150 call scale(ymin1,ymax1,4,ymin,ymax,del) | |
| 580 | c | |
| 581 | c determine coordinates | |
| 582 | c | |
| 583 | 200 ycoor(1,k)=ymin | |
| 584 | pmin(k)=ymin | |
| 585 | small=del*1.0d-4 | |
| 586 | if (dabs(ycoor(1,k)).le.small) ycoor(1,k)=0.0d0 | |
| 587 | do 210 i=1,4 | |
| 588 | ycoor(i+1,k)=ycoor(i,k)+del | |
| 589 | if (dabs(ycoor(i+1,k)).le.small) ycoor(i+1,k)=0.0d0 | |
| 590 | 210 continue | |
| 591 | if (ilogy(k).eq.1) go to 230 | |
| 592 | do 220 i=1,5 | |
| 593 | 220 ycoor(i,k)=dexp(xlog10*ycoor(i,k)) | |
| 594 | 230 delplt(k)=del/dfloat(nwide4) | |
| 595 | locyt=locyt+npoint | |
| 596 | 235 continue | |
| 597 | c | |
| 598 | c count distinct coordinates | |
| 599 | c | |
| 600 | icoor(1)=1 | |
| 601 | jcoor(1)=1 | |
| 602 | numcor=1 | |
| 603 | if (kntr.eq.1) go to 290 | |
| 604 | do 250 i=2,kntr | |
| 605 | do 245 j=1,numcor | |
| 606 | l=jcoor(j) | |
| 607 | c... coordinates are *equal* if the most significant 24 bits agree | |
| 608 | do 240 k=1,5 | |
| 609 | y1=ycoor(k,i) | |
| 610 | y2=ycoor(k,l) | |
| 611 | if(y1.eq.0.0d0.and.y2.eq.0.0d0) go to 240 | |
| 612 | if(dabs((y1-y2)/dmax1(dabs(y1),dabs(y2))).ge.1.0d-7) go to 245 | |
| 613 | 240 continue | |
| 614 | icoor(i)=l | |
| 615 | go to 250 | |
| 616 | 245 continue | |
| 617 | icoor(i)=i | |
| 618 | numcor=numcor+1 | |
| 619 | jcoor(numcor)=i | |
| 620 | 250 continue | |
| 621 | c | |
| 622 | c print coordinates | |
| 623 | c | |
| 624 | 260 do 280 i=1,numcor | |
| 625 | asym(1)=asym1 | |
| 626 | asym(2)=asym2 | |
| 627 | ipos=2 | |
| 628 | do 270 j=1,kntr | |
| 629 | if (icoor(j).ne.jcoor(i)) go to 270 | |
| 630 | call move(asym,ipos,pltsym,j,1) | |
| 631 | ipos=ipos+1 | |
| 632 | 270 continue | |
| 633 | call move(asym,ipos,arprn,1,1) | |
| 634 | k=jcoor(i) | |
| 635 | if(iwide.ne.0) write(6,271) asym,(ycoor(j,k),j=1,5) | |
| 636 | 271 format(/1hx,2a8,4h----,1pd12.3,4(15x,d10.3)/26x,51(2h -)) | |
| 637 | if(iwide.eq.0) write(6,273) asym,(ycoor(j,k),j=1,5) | |
| 638 | 273 format(/1hx,2a8,1pd10.3,3(4x,d10.3),1x,d10.3/22x,29(2h -)) | |
| 639 | 280 continue | |
| 640 | go to 300 | |
| 641 | 290 if(iwide.ne.0) write(6,291) (ycoor(j,1),j=1,5) | |
| 642 | 291 format(/1hx,20x,1pd12.3,4(15x,d10.3)/26x,51(2h -)) | |
| 643 | if(iwide.eq.0) write(6,293) (ycoor(j,1),j=1,5) | |
| 644 | 293 format(/1hx,14x,1pd12.3,3(4x,d10.3),1x,d10.3/22x,29(2h -)) | |
| 645 | c | |
| 646 | c plotting | |
| 647 | c | |
| 648 | 300 aspot=ablnk | |
| 649 | do 320 i=1,numpnt | |
| 650 | xvar=value(locx+i) | |
| 651 | locyt=locy | |
| 652 | call copy8(agraph,aplot,13) | |
| 653 | do 310 k=1,kntr | |
| 654 | yvr=value(locyt+i) | |
| 655 | ktmp=icoor(k) | |
| 656 | ymin1=pmin(ktmp) | |
| 657 | jpoint=idint((yvr-ymin1)/delplt(k)+0.5d0)+1 | |
| 658 | if (jpoint.le.0) go to 306 | |
| 659 | if (jpoint.gt.nwide) go to 306 | |
| 660 | call move(aspot,1,aplot,jpoint,1) | |
| 661 | if (aspot.eq.ablnk) go to 303 | |
| 662 | if (aspot.eq.aper) go to 303 | |
| 663 | call move(aplot,jpoint,aletx,1,1) | |
| 664 | go to 306 | |
| 665 | 303 call move(aplot,jpoint,pltsym,k,1) | |
| 666 | 306 locyt=locyt+npoint | |
| 667 | 310 continue | |
| 668 | yvr=value(locy+i) | |
| 669 | if (ilogy(1).eq.1) go to 315 | |
| 670 | yvr=dexp(xlog10*yvr) | |
| 671 | 315 if(iwide.ne.0) write(6,316) xvar,yvr,aplot | |
| 672 | 316 format(1x,1pd10.3,2x,d10.3,2x,13a8) | |
| 673 | if(iwide.eq.0) write(6,317) xvar,yvr,(aplot(k),k=1,8) | |
| 674 | 317 format(1x,1pd10.3,1x,d10.3,7a8,a1) | |
| 675 | 320 continue | |
| 676 | c | |
| 677 | c finished | |
| 678 | c | |
| 679 | if(iwide.ne.0) write(6,331) | |
| 680 | 331 format(26x,51(2h -)//) | |
| 681 | if(iwide.eq.0) write(6,332) | |
| 682 | 332 format(21x,29(2h -)//) | |
| 683 | go to 500 | |
| 684 | c | |
| 685 | c too few points | |
| 686 | c | |
| 687 | 400 write (6,401) | |
| 688 | 401 format('0warning: too few points for plotting'/) | |
| 689 | 500 write (6,501) | |
| 690 | 501 format(1hy) | |
| 691 | return | |
| 692 | end | |
| 693 | subroutine scale(xmin,xmax,n,xminp,xmaxp,del) | |
| 694 | implicit double precision (a-h,o-z) | |
| 695 | c | |
| 696 | c this routine determines the 'optimal' scale to use for the plot of | |
| 697 | c some output variable. | |
| 698 | c | |
| 699 | c | |
| 700 | c adapted from algorithm 463 of 'collected algorithms of the cacm' | |
| 701 | c | |
| 702 | common /knstnt/ twopi,xlog2,xlog10,root2,rad,boltz,charge,ctok, | |
| 703 | 1 gmin,reltol,abstol,vntol,trtol,chgtol,eps0,epssil,epsox | |
| 704 | integer xxor | |
| 705 | dimension vint(5) | |
| 706 | data vint / 1.0d0,2.0d0,5.0d0,10.0d0,20.0d0 / | |
| 707 | data eps / 1.0d-12 / | |
| 708 | c | |
| 709 | c | |
| 710 | c... trap too-small data spread | |
| 711 | if(xmin.eq.0.0d0.and.xmax.eq.0.0d0) go to 4 | |
| 712 | if(dabs((xmax-xmin)/dmax1(dabs(xmin),dabs(xmax))).ge.1.0d-4) | |
| 713 | 1 go to 10 | |
| 714 | 4 continue | |
| 715 | if (xmin.ge.0.0d0) go to 5 | |
| 716 | xmax=0.5d0*xmin+eps | |
| 717 | xmin=1.5d0*xmin-eps | |
| 718 | go to 10 | |
| 719 | 5 xmax=1.5d0*xmin+eps | |
| 720 | xmin=0.5d0*xmin-eps | |
| 721 | c... find approximate interval size, normalized to [1,10] | |
| 722 | 10 a=(xmax-xmin)/dfloat(n) | |
| 723 | nal=idint(dlog10(a)) | |
| 724 | if (a.lt.1.0d0) nal=nal-1 | |
| 725 | xfact=dexp(xlog10*dfloat(nal)) | |
| 726 | b=a/xfact | |
| 727 | c... find closest permissible interval size | |
| 728 | do 20 i=1,3 | |
| 729 | if (b.lt.(vint(i)+eps)) go to 30 | |
| 730 | 20 continue | |
| 731 | i=4 | |
| 732 | c... compute interval size | |
| 733 | 30 del=vint(i)*xfact | |
| 734 | fm1=xmin/del | |
| 735 | m1=fm1 | |
| 736 | if (fm1.lt.0.0d0) m1=m1-1 | |
| 737 | if (dabs(dfloat(m1)+1.0d0-fm1).lt.eps) m1=m1+1 | |
| 738 | c... compute new maximum and minimum limits | |
| 739 | xminp=del*dfloat(m1) | |
| 740 | fm2=xmax/del | |
| 741 | m2=fm2+1.0d0 | |
| 742 | if (fm2.lt.(-1.0d0)) m2=m2-1 | |
| 743 | if (dabs(fm2+1.0d0-dfloat(m2)).lt.eps) m2=m2-1 | |
| 744 | xmaxp=del*dfloat(m2) | |
| 745 | np=m2-m1 | |
| 746 | c... check whether another loop required | |
| 747 | if (np.le.n) go to 40 | |
| 748 | i=i+1 | |
| 749 | go to 30 | |
| 750 | c... do final adjustments and correct for roundoff error(s) | |
| 751 | 40 nx=(n-np)/2 | |
| 752 | xminp=dmin1(xmin,xminp-dfloat(nx)*del) | |
| 753 | xmaxp=dmax1(xmax,xminp+dfloat(n)*del) | |
| 754 | return | |
| 755 | end | |
| 756 | subroutine fouran | |
| 757 | implicit double precision (a-h,o-z) | |
| 758 | c | |
| 759 | c this routine determines the fourier coefficients of a transient | |
| 760 | c analysis waveform. | |
| 761 | c | |
| 762 | common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem, | |
| 763 | 1 isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize, | |
| 764 | 2 junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr, | |
| 765 | 3 nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1, | |
| 766 | 4 lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd, | |
| 767 | 5 imynl,imvn,lcvn,loutpt,nsnod,nsmat,nsval,icnod,icmat,icval | |
| 768 | common /cirdat/ locate(50),jelcnt(50),nunods,ncnods,numnod,nstop, | |
| 769 | 1 nut,nlt,nxtrm,ndist,ntlin,ibr,numvs | |
| 770 | common /flags/ iprnta,iprntl,iprntm,iprntn,iprnto,limtim,limpts, | |
| 771 | 1 lvlcod,lvltim,itl1,itl2,itl3,itl4,itl5,igoof,nogo,keof | |
| 772 | common /miscel/ atime,aprog(3),adate,atitle(10),defl,defw,defad, | |
| 773 | 1 defas,rstats(50),iwidth,lwidth,nopage | |
| 774 | common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet, | |
| 775 | 1 xmu,mode,modedc,icalc,initf,method,iord,maxord,noncon,iterno, | |
| 776 | 2 itemno,nosolv,ipostp,iscrch | |
| 777 | common /knstnt/ twopi,xlog2,xlog10,root2,rad,boltz,charge,ctok, | |
| 778 | 1 gmin,reltol,abstol,vntol,trtol,chgtol,eps0,epssil,epsox | |
| 779 | common /tran/ tstep,tstop,tstart,delmax,tdmax,forfre,jtrflg | |
| 780 | common /outinf/ xincr,string(15),xstart,yvar(8),itab(8),itype(8), | |
| 781 | 1 ilogy(8),npoint,numout,kntr,numdgt | |
| 782 | common /blank/ value(1000) | |
| 783 | integer nodplc(64) | |
| 784 | complex*16 cvalue(32) | |
| 785 | equivalence (value(1),nodplc(1),cvalue(1)) | |
| 786 | c | |
| 787 | c | |
| 788 | dimension sinco(9),cosco(9) | |
| 789 | dimension fortit(4) | |
| 790 | data fortit / 8hfourier , 8hanalysis, 8h , 8h / | |
| 791 | data ablnk / 1h / | |
| 792 | c | |
| 793 | c | |
| 794 | forprd=1.0d0/forfre | |
| 795 | xstart=tstop-forprd | |
| 796 | kntr=1 | |
| 797 | xn=101.0d0 | |
| 798 | xincr=forprd/xn | |
| 799 | npoint=xn | |
| 800 | call getm8(locx,npoint) | |
| 801 | call getm8(locy,npoint) | |
| 802 | do 105 nknt=1,nfour | |
| 803 | itab(1)=nodplc(ifour+nknt) | |
| 804 | kfrout=itab(1) | |
| 805 | call ntrpl8(locx,locy,numpnt) | |
| 806 | dcco=0.0d0 | |
| 807 | call zero8(sinco,9) | |
| 808 | call zero8(cosco,9) | |
| 809 | loct=locy+1 | |
| 810 | ipnt=0 | |
| 811 | 10 yvr=value(loct+ipnt) | |
| 812 | dcco=dcco+yvr | |
| 813 | forfac=dfloat(ipnt)*twopi/xn | |
| 814 | arg=0.0d0 | |
| 815 | do 20 k=1,9 | |
| 816 | arg=arg+forfac | |
| 817 | sinco(k)=sinco(k)+yvr*dsin(arg) | |
| 818 | cosco(k)=cosco(k)+yvr*dcos(arg) | |
| 819 | 20 continue | |
| 820 | ipnt=ipnt+1 | |
| 821 | if (ipnt.ne.npoint) go to 10 | |
| 822 | dcco=dcco/xn | |
| 823 | forfac=2.0d0/xn | |
| 824 | do 30 k=1,9 | |
| 825 | sinco(k)=sinco(k)*forfac | |
| 826 | cosco(k)=cosco(k)*forfac | |
| 827 | 30 continue | |
| 828 | call title(0,72,1,fortit) | |
| 829 | ipos=1 | |
| 830 | call outnam(kfrout,1,string,ipos) | |
| 831 | call move(string,ipos,ablnk,1,7) | |
| 832 | jstop=(ipos+6)/8 | |
| 833 | write (6,61) (string(j),j=1,jstop) | |
| 834 | 61 format(' fourier components of transient response ',5a8///) | |
| 835 | write (6,71) dcco | |
| 836 | 71 format('0dc component =',1pd12.3/, | |
| 837 | 1 '0harmonic frequency fourier normalized phase no | |
| 838 | 2rmalized'/, | |
| 839 | 3 ' no (hz) component component (deg) pha | |
| 840 | 4se (deg)'//) | |
| 841 | iknt=1 | |
| 842 | freq1=forfre | |
| 843 | xnharm=1.0d0 | |
| 844 | call magphs(dcmplx(sinco(1),cosco(1)),xnorm,pnorm) | |
| 845 | phasen=0.0d0 | |
| 846 | write (6,81) iknt,freq1,xnorm,xnharm,pnorm,phasen | |
| 847 | 81 format(i6,1pd15.3,d12.3,0pf13.6,f10.3,f12.3/) | |
| 848 | thd=0.0d0 | |
| 849 | do 90 iknt=2,9 | |
| 850 | freq1=dfloat(iknt)*forfre | |
| 851 | call magphs(dcmplx(sinco(iknt),cosco(iknt)), | |
| 852 | 1 harm,phase) | |
| 853 | xnharm=harm/xnorm | |
| 854 | phasen=phase-pnorm | |
| 855 | thd=thd+xnharm*xnharm | |
| 856 | write (6,81) iknt,freq1,harm,xnharm,phase,phasen | |
| 857 | 90 continue | |
| 858 | thd=100.0d0*dsqrt(thd) | |
| 859 | write (6,101) thd | |
| 860 | 101 format (//5x,'total harmonic distortion = ',f12.6,' percent') | |
| 861 | 105 continue | |
| 862 | call clrmem(locx) | |
| 863 | call clrmem(locy) | |
| 864 | 110 return | |
| 865 | end |