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