projects / unix-history / blame
? search:
summary | tags | clone url | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Start development on BSD-SCCS
[unix-history] / .ref-BSD-3 / usr / src / cmd / spice / dctrans.f
CommitLineData
7316c58a
D		 1 subroutine dctran
2 implicit double precision (a-h,o-z)
3c
4c
5c this routine controls the dc transfer curve, dc operating point,
6c and transient analyses. the variables mode and modedc (defined below)
7c determine exactly which analysis is performed.
8c
9 common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem,
10 1 isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize,
11 2 junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr,
12 3 nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1,
13 4 lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd,
14 5 imynl,imvn,lcvn,loutpt,nsnod,nsmat,nsval,icnod,icmat,icval
15 common /miscel/ atime,aprog(3),adate,atitle(10),defl,defw,defad,
16 1 defas,rstats(50),iwidth,lwidth,nopage
17 common /cirdat/ locate(50),jelcnt(50),nunods,ncnods,numnod,nstop,
18 1 nut,nlt,nxtrm,ndist,ntlin,ibr,numvs
19 common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet,
20 1 xmu,mode,modedc,icalc,initf,method,iord,maxord,noncon,iterno,
21 2 itemno,nosolv,ipostp,iscrch
22 common /flags/ iprnta,iprntl,iprntm,iprntn,iprnto,limtim,limpts,
23 1 lvlcod,lvltim,itl1,itl2,itl3,itl4,itl5,igoof,nogo,keof
24 common /dc/ tcstar(2),tcstop(2),tcincr(2),icvflg,itcelm(2),kssop,
25 1 kinel,kidin,kovar,kidout
26 common /tran/ tstep,tstop,tstart,delmax,tdmax,forfre,jtrflg
27 common /cje/ maxtim,itime,icost
28 common /blank/ value(1000)
29 integer nodplc(64)
30 complex*16 cvalue(32)
31 equivalence (value(1),nodplc(1),cvalue(1))
32c
33c
34 dimension subtit(4,2)
35 dimension avhdr(3),avfrm(4)
36 data aendor /9.87654321d+27/
37 data avhdr / 8h( (2x,a4, 8h,3x,a7,3, 5hx)//) /
38 data avfrm / 8h( (1h ,a, 8h1,i3,1h), 8h,f10.4,3, 4hx)/) /
39 data anode, avltg / 4hnode, 7hvoltage /
40 data subtit / 8hsmall si, 8hgnal bia, 8hs soluti, 8hon ,
41 1 8hinitial , 8htransien, 8ht soluti, 8hon /
42 data lprn /1h(/
43c
44c the variables *mode*, *modedc*, and *initf* are used by spice to
45c keep track of the state of the analysis. the values of these flags
46c (and the corresponding meanings) are as follows:
47c
48c flag value meaning
49c ---- ----- -------
50c
51c mode 1 dc analysis (subtype defined by *modedc*)
52c 2 transient analysis
53c 3 ac analysis (small signal)
54c
55c modedc 1 dc operating point
56c 2 initial operating point for transient analysis
57c 3 dc transfer curve computation
58c
59c initf 1 converge with 'off' devices allowed to float
60c 2 initialize junction voltages
61c 3 converge with 'off' devices held 'off'
62c 4 store small-signal parameters away
63c 5 first timepoint in transient analysis
64c 6 prediction step
65c
66c note: *modedc* is only significant if *mode* = 1.
67c
68c
69c initialize
70c
71 call second(t1)
72c.. don't take any chances with lx3, set to large number
73 lx3=20000000
74 lx2=20000000
75 nolx2=0
76 nolx3=0
77 loctim=5
78c.. post-processing initialization
79 if(ipostp.eq.0) go to 1
80 numcur=jelcnt(9)
81 numpos=nunods+numcur
82 call getm8(ibuff,numpos)
83 numpos=numpos*4
84 if(numcur.eq.0) go to 1
85 loc=locate(9)
86 loccur=nodplc(loc+6)-1
87c... set up format
88 1 nvprln=4+(lwidth-72)/19
89 nvprln=min0(nvprln,ncnods-1)
90 ipos=2
91 call alfnum(nvprln,avfrm,ipos)
92 ipos=2
93 call alfnum(nvprln,avhdr,ipos)
94c... allocate storage
95 if (mode.eq.2) go to 5
96 need=4*nstop+nut+nlt+nxtrm
97 call avlm8(navl)
98 if(need.le.navl) go to 4
99c... not enough memory for dc operating point analysis
100 write(6,3) need,navl
101 3 format('0insufficient memory available for dc analysis.',/
102 1' memory required ',i6,', memory available ',i6,'.')
103 nogo=1
104 go to 1100
105 4 call getm8(lvnim1,nstop)
106 call getm8(ndiag,nstop)
107 call getm8(lvn,nstop+nstop+nut+nlt)
108 call getm8(lx0,nxtrm)
109 if (modedc.ne.3) go to 15
110 5 call getm8(lx1,nxtrm)
111 if(nolx2.eq.0) call getm8(lx2,nxtrm)
112 if (mode.ne.2) go to 12
113 if(nolx3.eq.0) call getm8(lx3,nxtrm)
114 call getm8(ltd,0)
115 12 call getm8(loutpt,0)
116 15 call crunch
117 lynl=lvn+nstop
118 lyu=lynl+nstop
119 lyl=lyu+nut
120 20 if (mode.eq.2) go to 500
121 time=0.0d0
122 ag(1)=0.0d0
123 call sorupd
124 if (modedc.eq.3) go to 300
125c
126c .... single point dc analysis
127c
128c
129c compute dc operating point
130c
131 100 initf=2
132 call iter8(itl1)
133 rstats(6)=rstats(6)+iterno
134 if (igoof.ne.0) go to 150
135 if (modedc.ne.1) go to 120
136 initf=4
137 call diode
138 call bjt
139 call jfet
140 call mosfet
141c
142c print operating point
143c
144 120 if ((mode.eq.1).and.(modedc.eq.2).and.(nosolv.ne.0)) go to 1000
145 call title(-1,lwidth,1,subtit(1,modedc))
146 write (6,avhdr) (anode,avltg,i=1,nvprln)
147 write (6,avfrm) (lprn,nodplc(junode+i),value(lvnim1+i),i=2,ncnods)
148 go to 1000
149c
150c no convergence
151c
152 150 nogo=1
153 write (6,151)
154 151 format('1*error*: no convergence in dc analysis'/'0last node vol'
155 1 ,'tages:'/)
156 write (6,avhdr) (anode,avltg,i=1,nvprln)
157 write (6,avfrm) (lprn,nodplc(junode+i),value(lvnim1+i),i=2,ncnods)
158 go to 1000
159c
160c .... dc transfer curves
161c
162 300 numout=jelcnt(41)+1
163 if(ipostp.ne.0) call pheadr(atitle)
164 itemp=itcelm(1)
165 locs=nodplc(itemp+1)
166 temval=value(locs+1)
167 icvfl2=1
168 if(itcelm(2).eq.0) go to 310
169 itemp=itcelm(2)
170 locs2=nodplc(itemp+1)
171 temv2=value(locs2+1)
172 value(locs2+1)=tcstar(2)
173 temp=dabs((tcstop(2)-tcstar(2))/tcincr(2))+0.5d0
174 icvfl2=idint(temp)+1
175 icvfl2=max0(icvfl2,1)
176 310 delta=tcincr(1)
177 do 320 i=1,7
178 delold(i)=delta
179 320 continue
180 icvfl1=icvflg/icvfl2
181 value(locs+1)=tcstar(1)
182 icalc=0
183 ical2=0
184 loctim=3
185 340 initf=2
186 call iter8(itl1)
187 rstats(4)=rstats(4)+iterno
188 call copy8(value(lx0+1),value(lx1+1),nxtrm)
189 if(nolx2.eq.0) call copy8(value(lx0+1),value(lx2+1),nxtrm)
190 if (igoof.ne.0) go to 450
191 go to 360
192 350 call getcje
193 if ((maxtim-itime).le.limtim) go to 460
194 initf=6
195 call iter8(itl2)
196 rstats(4)=rstats(4)+iterno
197 if (igoof.ne.0) go to 340
198c
199c store outputs
200c
201 360 call extmem(loutpt,numout)
202 loco=loutpt+icalc*numout
203 icalc=icalc+1
204 ical2=ical2+1
205 value(loco+1)=value(locs+1)
206 loc=locate(41)
207 370 if (loc.eq.0) go to 400
208 if (nodplc(loc+5).ne.0) go to 380
209 node1=nodplc(loc+2)
210 node2=nodplc(loc+3)
211 iseq=nodplc(loc+4)
212 value(loco+iseq)=value(lvnim1+node1)-value(lvnim1+node2)
213 loc=nodplc(loc)
214 go to 370
215 380 iptr=nodplc(loc+2)
216 iptr=nodplc(iptr+6)
217 iseq=nodplc(loc+4)
218 value(loco+iseq)=value(lvnim1+iptr)
219 loc=nodplc(loc)
220 go to 370
221c
222c increment source value
223c
224 400 if(ipostp.eq.0) go to 410
225 value(ibuff+1)=value(locs+1)
226 call copy8(value(lvnim1+2),value(ibuff+2),nunods-1)
227 if(numcur.ne.0) call copy8(value(lvnim1+loccur+1),
228 1 value(ibuff+nunods+1),numcur)
229 call fwrite(value(ibuff+1),numpos)
230 410 if (icalc.ge.icvflg) go to 490
231 if(ical2.ge.icvfl1) go to 480
232 if(nolx2.ne.0) go to 420
233 call ptrmem(lx2,itemp)
234 call ptrmem(lx1,lx2)
235 go to 430
236 420 call ptrmem(lx1,itemp)
237 430 call ptrmem(lx0,lx1)
238 call ptrmem(itemp,lx0)
239 value(locs+1)=tcstar(1)+dfloat(ical2)*delta
240 go to 350
241c
242c no convergence
243c
244 450 itemp=itcelm(1)
245 loce=nodplc(itemp+1)
246 write (6,451) value(loce),value(locs+1)
247 451 format('1*error*: no convergence in dc transfer curves at ',a8,
248 1 ' = ',1pd10.3/'0last node voltages:'/)
249 write (6,avhdr) (anode,avltg,i=1,nvprln)
250 write (6,avfrm) (lprn,nodplc(junode+i),value(lvnim1+i),i=2,ncnods)
251 go to 470
252 460 write (6,461)
253 461 format('0*error*: cpu time limit exceeded ... analysis stopped'/)
254 470 nogo=1
255 go to 490
256c... reset first sweep variable ... step second
257 480 ical2=0
258 value(locs+1)=tcstar(1)
259 value(locs2+1)=value(locs2+1)+tcincr(2)
260 go to 340
261c
262c finished with dc transfer curves
263c
264 490 value(locs+1)=temval
265 if(itcelm(2).ne.0) value(locs2+1)=temv2
266 if(ipostp.eq.0) go to 1000
267 value(ibuff+1)=aendor
268 call fwrite(value(ibuff+1),numpos)
269 go to 1000
270c
271c .... transient analysis
272c
273 500 numout=jelcnt(42)+1
274 if(ipostp.ne.0) call pheadr(atitle)
275 numese=jelcnt(2)+jelcnt(3)+jelcnt(11)+jelcnt(12)+jelcnt(13)
276 1 +jelcnt(14)
277 if (numese.eq.0) delmax=dmin1(delmax,tstep)
278 initf=5
279 iord=1
280 loctim=9
281 icalc=0
282 numtp=0
283 numrtp=0
284 numnit=0
285 time=0.0d0
286 ibkflg=1
287 delbkp=delmax
288 nbkpt=1
289 delta=delmax
290 do 510 i=1,7
291 delold(i)=delta
292 510 continue
293 delmin=1.0d-9*delmax
294 go to 650
295c
296c increment time, update sources, and solve next timepoint
297c
298 600 time=time+delta
299 call sorupd
300 if (nogo.ne.0) go to 950
301 call getcje
302 if ((maxtim-itime).le.limtim) go to 920
303 if ((itl5.ne.0).and.(numnit.ge.itl5)) go to 905
304 call comcof
305 if (initf.ne.5) initf=6
306 itrlim=itl4
307 if ((numtp.eq.0).and.(nosolv.ne.0)) itrlim=itl1
308 call iter8(itrlim)
309 if(itl5.ne.0) numnit=numnit+iterno
310 numtp=numtp+1
311 if (numtp.ne.1) go to 605
312 if(nolx2.eq.0) call copy8(value(lx1+1),value(lx2+1),nxtrm)
313 if(nolx3.eq.0) call copy8(value(lx1+1),value(lx3+1),nxtrm)
314c.. note that time-point is cut when itrlim exceeded regardless
315c.. of which time-step contol is specified thru 'lvltim'.
316 605 if (igoof.eq.0) go to 610
317 jord=iord
318 iord=1
319 if (jord.ge.5) call clrmem(lx7)
320 if (jord.ge.4) call clrmem(lx6)
321 if (jord.ge.3) call clrmem(lx5)
322 if ((jord.ge.2).and.(method.ne.1)) call clrmem(lx4)
323 igoof=0
324 time=time-delta
325 delta=delta/8.0d0
326 go to 620
327 610 delnew=delta
328 if (numtp.eq.1) go to 630
329 call trunc(delnew)
330 if (delnew.ge.(0.9d0*delta)) go to 630
331 time=time-delta
332 delta=delnew
333 620 numrtp=numrtp+1
334 ibkflg=0
335 delold(1)=delta
336 if (delta.ge.delmin) go to 600
337 time=time+delta
338 go to 900
339c.. time-point accepted
340 630 call copy8(delold(1),delold(2),6)
341 delta=delnew
342 delold(1)=delta
343c
344c determine order of integration method
345c
346c... skip if trapezoidal algorithm used
347 if ((method.eq.1).and.(iord.eq.2)) go to 650
348 if (numtp.eq.1) go to 650
349 ordrat=1.05d0
350 if (iord.gt.1) go to 635
351 iord=2
352 call trunc(delnew)
353 iord=1
354 if ((delnew/delta).le.ordrat) go to 650
355 if (maxord.le.1) go to 650
356 iord=2
357 if (method.eq.1) go to 650
358 call getm8(lx4,nxtrm)
359 go to 650
360 635 if (iord.lt.maxord) go to 640
361 iord=iord-1
362 call trunc(delnew)
363 iord=iord+1
364 if ((delnew/delta).le.ordrat) go to 650
365 go to 642
366 640 iord=iord-1
367 call trunc(delnew)
368 iord=iord+1
369 if ((delnew/delta).le.ordrat) go to 645
370 642 iord=iord-1
371 if (iord.eq.1) call clrmem(lx4)
372 if (iord.eq.2) call clrmem(lx5)
373 if (iord.eq.3) call clrmem(lx6)
374 if (iord.eq.4) call clrmem(lx7)
375 go to 650
376 645 iord=iord+1
377 call trunc(delnew)
378 iord=iord-1
379 if ((delnew/delta).le.ordrat) go to 650
380 iord=iord+1
381 if (iord.eq.2) call getm8(lx4,nxtrm)
382 if (iord.eq.3) call getm8(lx5,nxtrm)
383 if (iord.eq.4) call getm8(lx6,nxtrm)
384 if (iord.eq.5) call getm8(lx7,nxtrm)
385c
386c store outputs
387c
388 650 if ((time+delta).le.tstart) go to 685
389 if ((numtp.eq.0).and.(nosolv.ne.0)) go to 685
390 call extmem(loutpt,numout)
391 loco=loutpt+icalc*numout
392 icalc=icalc+1
393 value(loco+1)=time
394 loc=locate(42)
395 670 if (loc.eq.0) go to 682
396 if (nodplc(loc+5).ne.0) go to 680
397 node1=nodplc(loc+2)
398 node2=nodplc(loc+3)
399 iseq=nodplc(loc+4)
400 value(loco+iseq)=value(lvnim1+node1)-value(lvnim1+node2)
401 loc=nodplc(loc)
402 go to 670
403 680 iptr=nodplc(loc+2)
404 iptr=nodplc(iptr+6)
405 iseq=nodplc(loc+4)
406 value(loco+iseq)=value(lvnim1+iptr)
407 loc=nodplc(loc)
408 go to 670
409 682 if(ipostp.eq.0) go to 684
410 value(ibuff+1)=time
411 call copy8(value(lvnim1+2),value(ibuff+2),nunods-1)
412 if(numcur.ne.0) call copy8(value(lvnim1+loccur+1),
413 1 value(ibuff+nunods+1),numcur)
414 call fwrite(value(ibuff+1),numpos)
415 684 continue
416 685 if (jelcnt(17).eq.0) go to 694
417 call sizmem(ltd,ltdsiz)
418 numtd=ltdsiz/ntlin
419 if (numtd.le.3) go to 689
420 baktim=time-tdmax
421 if (baktim.lt.0.0d0) go to 689
422 lcntr=0
423 ltemp=ltd
424 do 686 i=1,numtd
425 if (value(ltemp+1).ge.baktim) go to 687
426 ltemp=ltemp+ntlin
427 lcntr=lcntr+1
428 686 continue
429 go to 689
430 687 if (lcntr.le.2) go to 689
431 lcntr=lcntr-2
432 nwords=lcntr*ntlin
433 ltemp=ltemp-ntlin-ntlin
434 call copy8(value(ltemp+1),value(ltd+1),ltdsiz-nwords)
435 call relmem(ltd,nwords)
436 call sizmem(ltd,ltdsiz)
437 689 call extmem(ltd,ntlin)
438 ltdptr=ltd+ltdsiz
439 value(ltdptr+1)=time
440 loc=locate(17)
441 690 if (loc.eq.0) go to 693
442 locv=nodplc(loc+1)
443 z0=value(locv+1)
444 node1=nodplc(loc+2)
445 node2=nodplc(loc+3)
446 node3=nodplc(loc+4)
447 node4=nodplc(loc+5)
448 ibr1=nodplc(loc+8)
449 ibr2=nodplc(loc+9)
450 lspot=nodplc(loc+30)+ltdptr
451 if ((initf.eq.5).and.(nosolv.ne.0)) go to 691
452 value(lspot)=value(lvnim1+node3)-value(lvnim1+node4)
453 1 +value(lvnim1+ibr2)*z0
454 value(lspot+1)=value(lvnim1+node1)-value(lvnim1+node2)
455 1 +value(lvnim1+ibr1)*z0
456 go to 692
457 691 value(lspot)=value(locv+7)+value(locv+8)*z0
458 value(lspot+1)=value(locv+5)+value(locv+6)*z0
459 692 loc=nodplc(loc)
460 go to 690
461c
462c add two *fake* backpoints to ltd for interpolation near time=0.0d0
463c
464 693 if (numtd.ne.0) go to 694
465 call extmem(ltd,ntlin+ntlin)
466 call copy8(value(ltd+1),value(ltd+ntlin+1),ntlin)
467 call copy8(value(ltd+1),value(ltd+2*ntlin+1),ntlin)
468 value(ltd+2*ntlin+1)=time
469 value(ltd+ntlin+1)=time-delta
470 value(ltd+1)=time-delta-delta
471c
472c rotate state vector storage
473c
474 694 go to (710,706,702,698,696,696), iord
475 696 call ptrmem(lx7,itemp)
476 call ptrmem(lx6,lx7)
477 go to 700
478 698 call ptrmem(lx6,itemp)
479 700 call ptrmem(lx5,lx6)
480 go to 704
481 702 call ptrmem(lx5,itemp)
482 704 call ptrmem(lx4,lx5)
483 go to 708
484 706 if (method.eq.1) go to 710
485 call ptrmem(lx4,itemp)
486 708 call ptrmem(lx3,lx4)
487 go to 713
488 710 if(nolx3.eq.0) go to 712
489 if(nolx2.eq.0) go to 711
490 call ptrmem(lx1,itemp)
491 go to 714
492 711 call ptrmem(lx2,itemp)
493 call ptrmem(lx1,lx2)
494 go to 714
495 712 call ptrmem(lx3,itemp)
496 713 call ptrmem(lx2,lx3)
497 call ptrmem(lx1,lx2)
498 714 call ptrmem(lx0,lx1)
499 call ptrmem(itemp,lx0)
500c
501c check breakpoints
502c
503 750 if (ibkflg.eq.0) go to 760
504c.. just accepted analysis at breakpoint
505 jord=iord
506 iord=1
507 if (jord.ge.5) call clrmem(lx7)
508 if (jord.ge.4) call clrmem(lx6)
509 if (jord.ge.3) call clrmem(lx5)
510 if ((jord.ge.2).and.(method.ne.1)) call clrmem(lx4)
511 ibkflg=0
512 nbkpt=nbkpt+1
513 if (nbkpt.gt.numbkp) go to 950
514 temp=dmin1(delbkp,value(lsbkpt+nbkpt)-time)
515 delta=dmin1(delta,0.1d0*temp,delmax)
516 if (numtp.eq.0) delta=delta/10.0d0
517 delold(1)=delta
518 go to 600
519 760 del1=value(lsbkpt+nbkpt)-time
520 if ((1.01d0*delta).le.del1) go to 600
521 ibkflg=1
522 delbkp=delta
523 delta=del1
524 delold(1)=delta
525 go to 600
526c
527c transient analysis failed
528c
529 900 write (6,901)
530 901 format('1*error*: internal timestep too small in transient analys
531 1is'/)
532 go to 910
533 905 write (6,906) itl5
534 906 format('1*error*: transient analysis iterations exceed limit of '
535 1,i5,/'0this limit may be overridden using the itl5 parameter on th
536 2e .option card')
537 910 write (6,911) time,delta,numnit
538 911 format(1h0,10x,'time = ',1pd12.5,'; delta = ',d12.5,'; numnit =
539 1',i6/)
540 write (6,916)
541 916 format(1h0/'0last node voltages:'/)
542 write (6,avhdr) (anode,avltg,i=1,nvprln)
543 write (6,avfrm) (lprn,nodplc(junode+i),value(lvnim1+i),i=2,ncnods)
544 go to 930
545 920 write (6,921) time
546 921 format('0*error*: cpu time limit exceeded in transient analysis '
547 1 ,'at time = ',1pd13.6/)
548 930 nogo=1
549c
550c finished with transient analysis
551c
552 950 rstats(10)=rstats(10)+numnit
553 rstats(30)=rstats(30)+numtp
554 rstats(31)=rstats(31)+numrtp
555 rstats(32)=rstats(32)+numnit
556 if(ipostp.eq.0) go to 1000
557 value(ibuff+1)=aendor
558 call fwrite(value(ibuff+1),numpos)
559c
560c return unneeded memory
561c
562 1000 if (mode.eq.2) go to 1010
563 if (modedc.ne.3) go to 1100
564 1010 call clrmem(lvnim1)
565 call clrmem(lx0)
566 call clrmem(ndiag)
567 call clrmem(lvn)
568 call clrmem(lx1)
569 if(nolx2.eq.0) call clrmem(lx2)
570 if ((mode.eq.1).and.(modedc.eq.3)) go to 1020
571 if(nolx3.eq.0) call clrmem(lx3)
572 if (mode.eq.1) go to 1020
573 call clrmem(ltd)
574 if (iord.eq.1) go to 1020
575 if (method.eq.1) go to 1020
576 call clrmem(lx4)
577 if (iord.eq.2) go to 1020
578 call clrmem(lx5)
579 if (iord.eq.3) go to 1020
580 call clrmem(lx6)
581 if (iord.eq.4) go to 1020
582 call clrmem(lx7)
583 1020 call extmem(loutpt,2*numout)
584 1100 if(ipostp.ne.0) call clrmem(ibuff)
585 call second(t2)
586 rstats(loctim)=rstats(loctim)+t2-t1
587 return
588 end
589 subroutine fwrite(iarray,numwds)
590c
591c write onto 'punch' file numwds 16-bit words starting
592c with location iarray(/1/)
593c
594 integer iarray(1)
595c
596c... data must be written out in 40 word (80 byte) chunks
597c
598 integer idata(20)
599 numwd4=(numwds+1)/2
600 numblk=(numwd4-1)/20+1
601 kntr=1
602 numlft=numwd4
603 do 10 i=1,numblk
604 kstop=min0(numlft,20)
605 call copy4(iarray(kntr),idata(1),kstop)
606 write(7) idata
607 kntr=kntr+20
608 numlft=numlft-20
609 10 continue
610 return
611 end
612 subroutine pheadr(aheadr)
613 implicit double precision (a-h,o-z)
614 common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem,
615 1 isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize,
616 2 junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr,
617 3 nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1,
618 4 lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd,
619 5 imynl,imvn,lcvn,loutpt,nsnod,nsmat,nsval,icnod,icmat,icval
620 common /cirdat/ locate(50),jelcnt(50),nunods,ncnods,numnod,nstop,
621 1 nut,nlt,nxtrm,ndist,ntlin,ibr,numvs
622 common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet,
623 1 xmu,mode,modedc,icalc,initf,method,iord,maxord,noncon,iterno,
624 2 itemno,nosolv,ipostp,iscrch
625 common /dc/ tcstar(2),tcstop(2),tcincr(2),icvflg,itcelm(2),kssop,
626 1 kinel,kidin,kovar,kidout
627 common /miscel/ atime,aprog(3),adate,atitle(10),defl,defw,defad,
628 1 defas,rstats(50),iwidth,lwidth,nopage
629 common /blank/ value(1000)
630 integer nodplc(64)
631 complex*16 cvalue(32)
632 integer*2 int2,nodpl2(128)
633 equivalence (value(1),nodpl2(1))
634 equivalence (value(1),nodplc(1),cvalue(1))
635 dimension aheadr(10)
636c
637c put out the header records onto the post-processing file
638c routine is used for all analysis modes (mode=1,2,3)
639c
640 dimension xtype(2)
641 data xtype /4htime,4hfreq/
642 data ablnk,aletv,aleti /1h ,1hv,1hi/
643 call getm8(ibuff,12)
644 call copy8(aheadr(1),value(ibuff+1),10)
645 value(ibuff+11)=adate
646 value(ibuff+12)=atime
647 call fwrite(value(ibuff+1),48)
648 numout=nunods+jelcnt(9)
649 info=3
650 call getm8(inames,numout)
651 call getm4(itypes,numout)
652 call getm4(iseqs,numout)
653 itype2=itypes*2
654 iseq2=iseqs*2
655 iknt=1
656 nodpl2(iseq2+1)=1
657 if(mode.ne.1) go to 10
658 loc=itcelm(1)
659 locv=nodplc(loc+1)
660 value(inames+1)=value(locv)
661 anam=ablnk
662 call move(anam,1,value(locv),1,1)
663 ityp=0
664 if(anam.eq.aletv) ityp=3
665 if(anam.eq.aleti) ityp=4
666 nodpl2(itype2+1)=ityp
667 go to 20
668 10 value(inames+1)=xtype(mode-1)
669 nodpl2(itype2+1)=mode-1
670 20 do 30 i=2,nunods
671 nodpl2(itype2+i)=3
672 nodpl2(iseq2+i)=i
673 value(inames+i)=ablnk
674 ipos=1
675 call alfnum(nodplc(junode+i),value(inames+i),ipos)
676 30 continue
677 loc=locate(9)
678 iknt=nunods
679 40 if(loc.eq.0) go to 50
680 iknt=iknt+1
681 nodpl2(itype2+iknt)=4
682 nodpl2(iseq2+iknt)=iknt
683 locv=nodplc(loc+1)
684 value(inames+iknt)=value(locv)
685 loc=nodplc(loc)
686 go to 40
687 50 int2=numout
688 call fwrite(int2,1)
689 int2=info
690 call fwrite(int2,1)
691 nwds=numout*4
692 call fwrite(value(inames+1),nwds)
693 call fwrite(nodpl2(itype2+1),numout)
694 call fwrite(nodpl2(iseq2+1),numout)
695 call clrmem(ibuff)
696 call clrmem(inames)
697 call clrmem(itypes)
698 call clrmem(iseqs)
699 return
700 end
701 subroutine comcof
702 implicit double precision (a-h,o-z)
703c
704c this routine calculates the timestep-dependent terms used in the
705c numerical integration.
706c
707 common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem,
708 1 isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize,
709 2 junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr,
710 3 nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1,
711 4 lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd,
712 5 imynl,imvn,lcvn,loutpt,nsnod,nsmat,nsval,icnod,icmat,icval
713 common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet,
714 1 xmu,mode,modedc,icalc,initf,method,iord,maxord,noncon,iterno,
715 2 itemno,nosolv,ipostp,iscrch
716 common /flags/ iprnta,iprntl,iprntm,iprntn,iprnto,limtim,limpts,
717 1 lvlcod,lvltim,itl1,itl2,itl3,itl4,itl5,igoof,nogo,keof
718 common /blank/ value(1000)
719 integer nodplc(64)
720 complex*16 cvalue(32)
721 equivalence (value(1),nodplc(1),cvalue(1))
722 dimension gmat(7,7)
723c
724c compute coefficients for particular integration method
725c
726 if (method.ne.1) go to 5
727 if (iord.eq.1) go to 5
728c... trapezoidal method
729 ag(1)=1.0d0/delta/(1.0d0-xmu)
730 ag(2)=xmu/(1.0d0-xmu)
731 go to 200
732c
733c construct gear coefficient matrix
734c
735 5 istop=iord+1
736 call zero8(ag,istop)
737 ag(2)=-1.0d0
738 do 10 i=1,istop
739 gmat(1,i)=1.0d0
740 10 continue
741 do 20 i=2,istop
742 gmat(i,1)=0.0d0
743 20 continue
744 arg=0.0d0
745 do 40 i=2,istop
746 arg=arg+delold(i-1)
747 arg1=1.0d0
748 do 30 j=2,istop
749 arg1=arg1*arg
750 gmat(j,i)=arg1
751 30 continue
752 40 continue
753c
754c solve for gear coefficients ag(*)
755c
756c
757c lu decomposition
758c
759 do 70 i=2,istop
760 jstart=i+1
761 if (jstart.gt.istop) go to 70
762 do 60 j=jstart,istop
763 gmat(j,i)=gmat(j,i)/gmat(i,i)
764 do 50 k=jstart,istop
765 gmat(j,k)=gmat(j,k)-gmat(j,i)*gmat(i,k)
766 50 continue
767 60 continue
768 70 continue
769c
770c forward substitution
771c
772 do 90 i=2,istop
773 jstart=i+1
774 if (jstart.gt.istop) go to 90
775 do 80 j=jstart,istop
776 ag(j)=ag(j)-gmat(j,i)*ag(i)
777 80 continue
778 90 continue
779c
780c backward substitution
781c
782 ag(istop)=ag(istop)/gmat(istop,istop)
783 ir=istop
784 do 110 i=2,istop
785 jstart=ir
786 ir=ir-1
787 do 100 j=jstart,istop
788 ag(ir)=ag(ir)-gmat(ir,j)*ag(j)
789 100 continue
790 ag(ir)=ag(ir)/gmat(ir,ir)
791 110 continue
792c
793c finished
794c
795 200 return
796 end
797 subroutine trunc(delnew)
798 implicit double precision (a-h,o-z)
799c
800c this routine determines the new transient stepsize by either
801c calling terr to estimate the local truncation error, or by checking
802c on the number of iterations needed to converge at the last timepoint.
803c
804 common /cirdat/ locate(50),jelcnt(50),nunods,ncnods,numnod,nstop,
805 1 nut,nlt,nxtrm,ndist,ntlin,ibr,numvs
806 common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet,
807 1 xmu,mode,modedc,icalc,initf,method,iord,maxord,noncon,iterno,
808 2 itemno,nosolv,ipostp,iscrch
809 common /flags/ iprnta,iprntl,iprntm,iprntn,iprnto,limtim,limpts,
810 1 lvlcod,lvltim,itl1,itl2,itl3,itl4,itl5,igoof,nogo,keof
811 common /tran/ tstep,tstop,tstart,delmax,tdmax,forfre,jtrflg
812 common /blank/ value(1000)
813 integer nodplc(64)
814 complex*16 cvalue(32)
815 equivalence (value(1),nodplc(1),cvalue(1))
816c
817c
818 if (lvltim.ne.0) go to 5
819 delnew=dmin1(tstep,delmax)
820 return
821 5 if (lvltim.ne.1) go to 10
822 delnew=delta
823 if (iterno.gt.itl3) return
824 delnew=dmin1(2.0d0*delta,tstep,delmax)
825 return
826c
827c capacitors
828c
829 10 delnew=1.0d20
830 loc=locate(2)
831 20 if (loc.eq.0) go to 30
832 loct=nodplc(loc+8)
833 call terr(loct,delnew)
834 loc=nodplc(loc)
835 go to 20
836c
837c inductors
838c
839 30 loc=locate(3)
840 40 if (loc.eq.0) go to 50
841 loct=nodplc(loc+11)
842 call terr(loct,delnew)
843 loc=nodplc(loc)
844 go to 40
845c
846c diodes
847c
848 50 loc=locate(11)
849 60 if (loc.eq.0) go to 70
850 loct=nodplc(loc+11)
851 call terr(loct+3,delnew)
852 loc=nodplc(loc)
853 go to 60
854c
855c bjts
856c
857 70 loc=locate(12)
858 80 if (loc.eq.0) go to 90
859 loct=nodplc(loc+22)
860 call terr(loct+8,delnew)
861 call terr(loct+10,delnew)
862 call terr(loct+12,delnew)
863 loc=nodplc(loc)
864 go to 80
865c
866c jfets
867c
868 90 loc=locate(13)
869 100 if (loc.eq.0) go to 110
870 loct=nodplc(loc+19)
871 call terr(loct+9,delnew)
872 call terr(loct+11,delnew)
873 loc=nodplc(loc)
874 go to 100
875c
876c mosfets
877c
878 110 loc=locate(14)
879 120 if (loc.eq.0) go to 200
880 loct=nodplc(loc+26)
881 call terr(loct+12,delnew)
882 call terr(loct+14,delnew)
883 call terr(loct+16,delnew)
884 call terr(loct+18,delnew)
885 call terr(loct+20,delnew)
886 loc=nodplc(loc)
887 go to 120
888c
889c delta is allowed only to double at each timepoint
890c
891 200 delnew=dmin1(2.0d0*delta,delnew,delmax)
892 return
893 end
894 subroutine terr(loct,delnew)
895 implicit double precision (a-h,o-z)
896c
897c this routine estimates the local truncation error for a particular
898c circuit element. it then computes the appropriate stepsize which
899c should be used.
900c
901 common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem,
902 1 isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize,
903 2 junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr,
904 3 nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1,
905 4 lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd,
906 5 imynl,imvn,lcvn,loutpt,nsnod,nsmat,nsval,icnod,icmat,icval
907 common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet,
908 1 xmu,mode,modedc,icalc,initf,method,iord,maxord,noncon,iterno,
909 2 itemno,nosolv,ipostp,iscrch
910 common /knstnt/ twopi,xlog2,xlog10,root2,rad,boltz,charge,ctok,
911 1 gmin,reltol,abstol,vntol,trtol,chgtol,eps0,epssil,epsox
912 common /blank/ value(1000)
913 integer nodplc(64)
914 complex*16 cvalue(32)
915 equivalence (value(1),nodplc(1),cvalue(1))
916c
917c
918 dimension qcap(1),ccap(1),diff(8),deltmp(7),coef(6)
919 equivalence (qcap(1),value(1)),(ccap(1),value(2))
920 data coef / 5.000000000d-1, 2.222222222d-1, 1.363636364d-1,
921 1 9.600000000d-2, 7.299270073d-2, 5.830903790d-2 /
922 data xtwelv / 8.333333333d-2 /
923c
924c
925 tol=reltol*dmax1(dabs(ccap(lx0+loct)),dabs(ccap(lx1+loct)))+abstol
926 ctol=reltol*dmax1(dabs(qcap(lx0+loct)),dabs(qcap(lx1+loct)),
927 1 chgtol)/delta
928 tol=dmax1(tol,ctol)
929c
930c determine divided differences
931c
932 go to (6,5,4,3,2,1), iord
933 1 diff(8)=qcap(lx7+loct)
934 2 diff(7)=qcap(lx6+loct)
935 3 diff(6)=qcap(lx5+loct)
936 4 diff(5)=qcap(lx4+loct)
937 5 diff(4)=qcap(lx3+loct)
938 6 diff(3)=qcap(lx2+loct)
939 diff(2)=qcap(lx1+loct)
940 diff(1)=qcap(lx0+loct)
941 istop=iord+1
942 do 10 i=1,istop
943 deltmp(i)=delold(i)
944 10 continue
945 20 do 30 i=1,istop
946 diff(i)=(diff(i)-diff(i+1))/deltmp(i)
947 30 continue
948 istop=istop-1
949 if (istop.eq.0) go to 100
950 do 40 i=1,istop
951 deltmp(i)=deltmp(i+1)+delold(i)
952 40 continue
953 go to 20
954c
955c diff(1) contains divided difference
956c
957 100 const=coef(iord)
958 if ((method.eq.1).and.(iord.eq.2)) const=xtwelv
959 del=trtol*tol/dmax1(abstol,const*dabs(diff(1)))
960 if (iord.eq.1) go to 200
961 if (iord.ge.3) go to 150
962 del=dsqrt(del)
963 go to 200
964 150 del=dexp(dlog(del)/dfloat(iord))
965 200 delnew=dmin1(delnew,del)
966 return
967 end
968 subroutine sorupd
969 implicit double precision (a-h,o-z)
970c
971c this routine updates the independent voltage and current sources
972c used in the circuit. it also updates the ltd table (which contains
973c previous (delayed) values of the sources used to model transmission
974c lines).
975c
976 common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem,
977 1 isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize,
978 2 junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr,
979 3 nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1,
980 4 lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd,
981 5 imynl,imvn,lcvn,loutpt,nsnod,nsmat,nsval,icnod,icmat,icval
982 common /cirdat/ locate(50),jelcnt(50),nunods,ncnods,numnod,nstop,
983 1 nut,nlt,nxtrm,ndist,ntlin,ibr,numvs
984 common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet,
985 1 xmu,mode,modedc,icalc,initf,method,iord,maxord,noncon,iterno,
986 2 itemno,nosolv,ipostp,iscrch
987 common /flags/ iprnta,iprntl,iprntm,iprntn,iprnto,limtim,limpts,
988 1 lvlcod,lvltim,itl1,itl2,itl3,itl4,itl5,igoof,nogo,keof
989 common /blank/ value(1000)
990 integer nodplc(64)
991 complex*16 cvalue(32)
992 equivalence (value(1),nodplc(1),cvalue(1))
993c
994c
995 do 500 id=9,10
996 loc=locate(id)
997 10 if (loc.eq.0) go to 500
998 locv=nodplc(loc+1)
999 locp=nodplc(loc+5)
1000 itype=nodplc(loc+4)+1
1001 go to (490,100,200,300,400,450), itype
1002c
1003c pulse source
1004c
1005 100 v1=value(locp+1)
1006 v2=value(locp+2)
1007 t1=value(locp+3)
1008 t2=value(locp+4)
1009 t3=value(locp+5)
1010 t4=value(locp+6)
1011 period=value(locp+7)
1012 time1=time
1013 if (time1.le.0.0d0) go to 160
1014 110 if (time1.lt.t1+period) go to 120
1015 time1=time1-period
1016 go to 110
1017 120 if (time1.lt.t4) go to 130
1018 value(locv+1)=v1
1019 go to 490
1020 130 if (time1.lt.t3) go to 140
1021 value(locv+1)=v2+(time1-t3)*(v1-v2)/(t4-t3)
1022 go to 490
1023 140 if (time1.lt.t2) go to 150
1024 value(locv+1)=v2
1025 go to 490
1026 150 if (time1.lt.t1) go to 160
1027 value(locv+1)=v1+(time1-t1)*(v2-v1)/(t2-t1)
1028 go to 490
1029 160 value(locv+1)=v1
1030 go to 490
1031c
1032c sinusoidal source
1033c
1034 200 v1=value(locp+1)
1035 v2=value(locp+2)
1036 omeg=value(locp+3)
1037 t1=value(locp+4)
1038 theta=value(locp+5)
1039 time1=time-t1
1040 if (time1.gt.0.0d0) go to 210
1041 value(locv+1)=v1
1042 go to 490
1043 210 if (theta.ne.0.0d0) go to 220
1044 value(locv+1)=v1+v2*dsin(omeg*time1)
1045 go to 490
1046 220 value(locv+1)=v1+v2*dsin(omeg*time1)*dexp(-time1*theta)
1047 go to 490
1048c
1049c exponential source
1050c
1051 300 v1=value(locp+1)
1052 v2=value(locp+2)
1053 t1=value(locp+3)
1054 tau1=value(locp+4)
1055 t2=value(locp+5)
1056 tau2=value(locp+6)
1057 time1=time
1058 if (time1.gt.t1) go to 310
1059 value(locv+1)=v1
1060 go to 490
1061 310 if (time1.gt.t2) go to 320
1062 value(locv+1)=v1+(v2-v1)*(1.0d0-dexp((t1-time1)/tau1))
1063 go to 490
1064 320 value(locv+1)=v1+(v2-v1)*(1.0d0-dexp((t1-time1)/tau1))
1065 1 +(v1-v2)*(1.0d0-dexp((t2-time1)/tau2))
1066 go to 490
1067c
1068c piecewise-linear source
1069c
1070 400 t1=value(locp+1)
1071 v1=value(locp+2)
1072 t2=value(locp+3)
1073 v2=value(locp+4)
1074 iknt=4
1075 410 if (time.le.t2) go to 420
1076 t1=t2
1077 v1=v2
1078 t2=value(locp+iknt+1)
1079 v2=value(locp+iknt+2)
1080 iknt=iknt+2
1081 go to 410
1082 420 value(locv+1)=v1+((time-t1)/(t2-t1))*(v2-v1)
1083 go to 490
1084c
1085c single-frequency fm
1086c
1087 450 v1=value(locp+1)
1088 v2=value(locp+2)
1089 omegc=value(locp+3)
1090 xmod=value(locp+4)
1091 omegs=value(locp+5)
1092 value(locv+1)=v1+v2*dsin(omegc*time+xmod*dsin(omegs*time))
1093 490 loc=nodplc(loc)
1094 go to 10
1095 500 continue
1096c
1097c update transmission line sources
1098c
1099 if (jelcnt(17).eq.0) go to 1000
1100 if (mode.ne.2) go to 1000
1101 call sizmem(ltd,ltdsiz)
1102 numtd=ltdsiz/ntlin
1103 if (numtd.lt.3) go to 900
1104 loc=locate(17)
1105 610 if (loc.eq.0) go to 1000
1106 locv=nodplc(loc+1)
1107 td=value(locv+2)
1108 baktim=time-td
1109 if (baktim.lt.0.0d0) go to 640
1110 ltdptr=nodplc(loc+30)
1111 icntr=2
1112 l1=ltd
1113 l2=l1+ntlin
1114 l3=l2+ntlin
1115 t1=value(l1+1)
1116 t2=value(l2+1)
1117 620 t3=value(l3+1)
1118 icntr=icntr+1
1119 if (baktim.le.t3) go to 630
1120 if (icntr.eq.numtd) go to 900
1121 l1=l2
1122 l2=l3
1123 l3=l2+ntlin
1124 t1=t2
1125 t2=t3
1126 go to 620
1127 630 dt1t2=t1-t2
1128 dt1t3=t1-t3
1129 dt2t3=t2-t3
1130 tdnom1=1.0d0/(dt1t2*dt1t3)
1131 tdnom2=-1.0d0/(dt1t2*dt2t3)
1132 tdnom3=1.0d0/(dt2t3*dt1t3)
1133 dtt1=baktim-t1
1134 dtt2=baktim-t2
1135 dtt3=baktim-t3
1136 tfact1=dtt2*dtt3*tdnom1
1137 tfact2=dtt1*dtt3*tdnom2
1138 tfact3=dtt1*dtt2*tdnom3
1139 value(locv+3)=value(l1+ltdptr+0)*tfact1+value(l2+ltdptr+0)*tfact2
1140 1 +value(l3+ltdptr+0)*tfact3
1141 value(locv+4)=value(l1+ltdptr+1)*tfact1+value(l2+ltdptr+1)*tfact2
1142 1 +value(l3+ltdptr+1)*tfact3
1143 640 loc=nodplc(loc)
1144 go to 610
1145c
1146c internal logic error: less than 3 entries in ltd
1147c
1148 900 nogo=1
1149 write (6,901) numtd,icntr
1150 901 format('0*abort*: internal spice error: sorupd: ',2i5/)
1151c
1152c finished
1153c
1154 1000 return
1155 end
1156 subroutine iter8(itlim)
1157 implicit double precision (a-h,o-z)
1158c
1159c this routine drives the newton-raphson iteration technique used to
1160c solve the set of nonlinear circuit equations.
1161c
1162 common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem,
1163 1 isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize,
1164 2 junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr,
1165 3 nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1,
1166 4 lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd,
1167 5 imynl,imvn,lcvn,loutpt,nsnod,nsmat,nsval,icnod,icmat,icval
1168 common /knstnt/ twopi,xlog2,xlog10,root2,rad,boltz,charge,ctok,
1169 1 gmin,reltol,abstol,vntol,trtol,chgtol,eps0,epssil,epsox
1170 common /cirdat/ locate(50),jelcnt(50),nunods,ncnods,numnod,nstop,
1171 1 nut,nlt,nxtrm,ndist,ntlin,ibr,numvs
1172 common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet,
1173 1 xmu,mode,modedc,icalc,initf,method,iord,maxord,noncon,iterno,
1174 2 itemno,nosolv,ipostp,iscrch
1175 common /flags/ iprnta,iprntl,iprntm,iprntn,iprnto,limtim,limpts,
1176 1 lvlcod,lvltim,itl1,itl2,itl3,itl4,itl5,igoof,nogo,keof
1177 common /blank/ value(1000)
1178 integer nodplc(64)
1179 complex*16 cvalue(32)
1180 equivalence (value(1),nodplc(1),cvalue(1))
1181c
1182c
1183 igoof=0
1184 iterno=0
1185 ndrflo=0
1186 noncon=0
1187 ipass=0
1188c
1189c construct linear equations and check convergence
1190c
1191 10 call load
1192 if ((mode.eq.1).and.(modedc.eq.2).and.(nosolv.ne.0)) go to 300
1193 iterno=iterno+1
1194 go to (20,30,40,50,50,50),initf
1195 20 if(mode.ne.1) go to 22
1196 call sizmem(nsnod,nic)
1197 if(nic.eq.0) go to 22
1198 if(ipass.ne.0) noncon=ipass
1199 ipass=0
1200 22 if(noncon.eq.0) go to 300
1201 go to 100
1202 30 initf=3
1203 40 if (noncon.eq.0) initf=1
1204 ipass=1
1205 go to 100
1206 50 initf=1
1207c
1208c solve equations for next iteration
1209c
1210 100 if (iterno.ge.itlim) go to 200
1211 110 call dcdcmp
1212 call dcsol
1213 120 if (igoof.eq.0) go to 130
1214 ndrflo=ndrflo+1
1215 igoof=0
1216 130 value(lvn+1)=0.0d0
1217 ntemp=noncon
1218 noncon=0
1219 if (ntemp.gt.0) go to 150
1220 if (iterno.eq.1) go to 150
1221 do 140 i=2,numnod
1222 vold=value(lvnim1+i)
1223 vnew=value(lvn+i)
1224 tol=reltol*dmax1(dabs(vold),dabs(vnew))+vntol
1225 if (dabs(vold-vnew).le.tol) go to 140
1226 noncon=noncon+1
1227 140 continue
1228 150 call copy8(value(lvn+1),value(lvnim1+1),nstop)
1229 go to 10
1230c
1231c no convergence
1232c
1233 200 igoof=1
1234 300 if (ndrflo.eq.0) go to 400
1235 write (6,301) ndrflo
1236 301 format('0warning: underflow occurred ',i4,' time(s)')
1237c
1238c finished
1239c
1240 400 return
1241 end
1242 subroutine load
1243 implicit double precision (a-h,o-z)
1244c
1245c this routine zeroes-out and then loads the coefficient matrix.
1246c the active devices and the controlled sources are loaded by separate
1247c subroutines.
1248c
1249 common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem,
1250 1 isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize,
1251 2 junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr,
1252 3 nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1,
1253 4 lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd,
1254 5 imynl,imvn,lcvn,loutpt,nsnod,nsmat,nsval,icnod,icmat,icval
1255 common /cirdat/ locate(50),jelcnt(50),nunods,ncnods,numnod,nstop,
1256 1 nut,nlt,nxtrm,ndist,ntlin,ibr,numvs
1257 common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet,
1258 1 xmu,mode,modedc,icalc,initf,method,iord,maxord,noncon,iterno,
1259 2 itemno,nosolv,ipostp,iscrch
1260 common /flags/ iprnta,iprntl,iprntm,iprntn,iprnto,limtim,limpts,
1261 1 lvlcod,lvltim,itl1,itl2,itl3,itl4,itl5,igoof,nogo,keof
1262 common /knstnt/ twopi,xlog2,xlog10,root2,rad,boltz,charge,ctok,
1263 1 gmin,reltol,abstol,vntol,trtol,chgtol,eps0,epssil,epsox
1264 common /blank/ value(1000)
1265 integer nodplc(64)
1266 complex*16 cvalue(32)
1267 equivalence (value(1),nodplc(1),cvalue(1))
1268c
1269c
1270 dimension qcap(1),ccap(1)
1271 equivalence (qcap(1),value(1)),(ccap(1),value(2))
1272 dimension find(1),vind(1)
1273 equivalence (find(1),value(1)),(vind(1),value(2))
1274c
1275c zero y matrix and current vector
1276c
1277 call zero8(value(lvn+1),nstop+nstop+nut+nlt)
1278c
1279c resistors
1280c
1281 loc=locate(1)
1282 20 if (loc.eq.0) go to 30
1283 locv=nodplc(loc+1)
1284 val=value(locv+1)
1285 locy=lynl+nodplc(loc+6)
1286 value(locy)=value(locy)+val
1287 locy=lynl+nodplc(loc+7)
1288 value(locy)=value(locy)+val
1289 locy=lynl+nodplc(loc+4)
1290 value(locy)=value(locy)-val
1291 locy=lynl+nodplc(loc+5)
1292 value(locy)=value(locy)-val
1293 loc=nodplc(loc)
1294 go to 20
1295c
1296c capacitors
1297c
1298 30 loc=locate(2)
1299 40 if (loc.eq.0) go to 100
1300 locv=nodplc(loc+1)
1301 node1=nodplc(loc+2)
1302 node2=nodplc(loc+3)
1303 lcoef=nodplc(loc+7)
1304 call sizmem(nodplc(loc+7),ncoef)
1305 loct=nodplc(loc+8)
1306 vcap=value(locv+2)
1307 if ((mode.eq.1).and.(initf.eq.2)) go to 45
1308 if ((nosolv.ne.0).and.(initf.eq.5)) go to 45
1309 vcap=value(lvnim1+node1)-value(lvnim1+node2)
1310 45 value(locv+3)=vcap
1311 if (mode.eq.1) go to 60
1312 if (initf.ne.6) go to 50
1313 qcap(lx0+loct)=qcap(lx1+loct)
1314 go to 60
1315 50 call evpoly(qcap(lx0+loct),-1,lcoef,ncoef,locv+2,1,loc+8)
1316 if (initf.ne.5) go to 60
1317 if (nosolv.eq.0) go to 55
1318 vcap=value(locv+2)
1319 value(locv+3)=vcap
1320 call evpoly(qcap(lx0+loct),-1,lcoef,ncoef,locv+2,1,loc+8)
1321 55 qcap(lx1+loct)=qcap(lx0+loct)
1322 60 call evpoly(value(locv+1),0,lcoef,ncoef,locv+2,1,loc+8)
1323 if (mode.eq.1) go to 90
1324 call intgr8(geq,ceq,value(locv+1),loct)
1325 ceq=ceq-geq*vcap+ag(1)*qcap(lx0+loct)
1326 if(initf.ne.5) go to 70
1327 ccap(lx1+loct)=ccap(lx0+loct)
1328 70 locy=lynl+nodplc(loc+10)
1329 value(locy)=value(locy)+geq
1330 locy=lynl+nodplc(loc+11)
1331 value(locy)=value(locy)+geq
1332 locy=lynl+nodplc(loc+5)
1333 value(locy)=value(locy)-geq
1334 locy=lynl+nodplc(loc+6)
1335 value(locy)=value(locy)-geq
1336 value(lvn+node1)=value(lvn+node1)-ceq
1337 value(lvn+node2)=value(lvn+node2)+ceq
1338 90 loc=nodplc(loc)
1339 go to 40
1340c
1341c inductors
1342c
1343 100 if (jelcnt(3).eq.0) go to 400
1344 if (mode.eq.1) go to 150
1345 if (initf.eq.6) go to 150
1346 loc=locate(3)
1347 110 if (loc.eq.0) go to 120
1348 locv=nodplc(loc+1)
1349 iptr=nodplc(loc+5)
1350 lcoef=nodplc(loc+10)
1351 call sizmem(nodplc(loc+10),ncoef)
1352 loct=nodplc(loc+11)
1353 cind=value(lvnim1+iptr)
1354 if ((mode.eq.1).and.(initf.eq.2)) cind=value(locv+2)
1355 if ((initf.eq.5).and.(nosolv.ne.0)) cind=value(locv+2)
1356 value(locv+3)=cind
1357 call evpoly(find(lx0+loct),-1,lcoef,ncoef,locv+2,1,loc+11)
1358 loc=nodplc(loc)
1359 go to 110
1360 120 loc=locate(4)
1361 130 if (loc.eq.0) go to 150
1362 locv=nodplc(loc+1)
1363 nl1=nodplc(loc+2)
1364 nl2=nodplc(loc+3)
1365 iptr1=nodplc(nl1+5)
1366 iptr2=nodplc(nl2+5)
1367 loct1=nodplc(nl1+11)
1368 loct2=nodplc(nl2+11)
1369 find(lx0+loct1)=find(lx0+loct1)+value(locv+1)*value(lvnim1+iptr2)
1370 find(lx0+loct2)=find(lx0+loct2)+value(locv+1)*value(lvnim1+iptr1)
1371 loc=nodplc(loc)
1372 go to 130
1373 150 loc=locate(3)
1374 160 if (loc.eq.0) go to 300
1375 locv=nodplc(loc+1)
1376 iptr=nodplc(loc+5)
1377 lcoef=nodplc(loc+10)
1378 call sizmem(nodplc(loc+10),ncoef)
1379 loct=nodplc(loc+11)
1380 cind=value(lvnim1+iptr)
1381 if ((mode.eq.1).and.(initf.eq.2)) cind=value(locv+2)
1382 if ((nosolv.ne.0).and.(initf.eq.5)) cind=value(locv+2)
1383 value(locv+3)=cind
1384 if (mode.ne.1) go to 200
1385 veq=0.0d0
1386 req=0.0d0
1387 go to 210
1388 200 if (initf.ne.6) go to 205
1389 find(lx0+loct)=find(lx1+loct)
1390 go to 210
1391 205 if (initf.ne.5) go to 210
1392 find(lx1+loct)=find(lx0+loct)
1393 210 call evpoly(value(locv+1),0,lcoef,ncoef,locv+2,1,loc+11)
1394 if (mode.eq.1) go to 250
1395 call intgr8(req,veq,value(locv+1),loct)
1396 if (ncoef.eq.1) go to 250
1397 veq=veq-req*cind+ag(1)*find(lx0+loct)
1398 250 value(lvn+iptr)=veq
1399 if(initf.ne.5) go to 260
1400 vind(lx1+loct)=vind(lx0+loct)
1401 260 locy=lynl+nodplc(loc+13)
1402 value(locy)=-req
1403 locy=lynl+nodplc(loc+6)
1404 value(locy)=1.0d0
1405 locy=lynl+nodplc(loc+7)
1406 value(locy)=-1.0d0
1407 locy=lynl+nodplc(loc+8)
1408 value(locy)=1.0d0
1409 locy=lynl+nodplc(loc+9)
1410 value(locy)=-1.0d0
1411 loc=nodplc(loc)
1412 go to 160
1413c
1414c mutual inductances
1415c
1416 300 loc=locate(4)
1417 310 if (loc.eq.0) go to 400
1418 locv=nodplc(loc+1)
1419 req=ag(1)*value(locv+1)
1420 locy=lynl+nodplc(loc+4)
1421 value(locy)=-req
1422 locy=lynl+nodplc(loc+5)
1423 value(locy)=-req
1424 loc=nodplc(loc)
1425 go to 310
1426c
1427c nonlinear controlled sources
1428c
1429 400 call nlcsrc
1430c
1431c voltage sources
1432c
1433 loc=locate(9)
1434 610 if (loc.eq.0) go to 700
1435 locv=nodplc(loc+1)
1436 iptr=nodplc(loc+6)
1437 value(lvn+iptr)=value(locv+1)
1438 locy=lynl+nodplc(loc+7)
1439 value(locy)=value(locy)+1.0d0
1440 locy=lynl+nodplc(loc+8)
1441 value(locy)=value(locy)-1.0d0
1442 locy=lynl+nodplc(loc+9)
1443 value(locy)=value(locy)+1.0d0
1444 locy=lynl+nodplc(loc+10)
1445 value(locy)=value(locy)-1.0d0
1446 loc=nodplc(loc)
1447 go to 610
1448c
1449c current sources
1450c
1451 700 loc=locate(10)
1452 710 if (loc.eq.0) go to 800
1453 locv=nodplc(loc+1)
1454 node1=nodplc(loc+2)
1455 node2=nodplc(loc+3)
1456 val=value(locv+1)
1457 value(lvn+node1)=value(lvn+node1)-val
1458 value(lvn+node2)=value(lvn+node2)+val
1459 loc=nodplc(loc)
1460 go to 710
1461c
1462c call device model routines
1463c
1464 800 call diode
1465 call bjt
1466 call jfet
1467 call mosfet
1468c
1469c transmission lines
1470c
1471 loc=locate(17)
1472 910 if (loc.eq.0) go to 980
1473 locv=nodplc(loc+1)
1474 z0=value(locv+1)
1475 y0=1.0d0/z0
1476 node1=nodplc(loc+2)
1477 node2=nodplc(loc+3)
1478 node3=nodplc(loc+4)
1479 node4=nodplc(loc+5)
1480 ibr1=nodplc(loc+8)
1481 ibr2=nodplc(loc+9)
1482 locy=lynl+nodplc(loc+10)
1483 value(locy)=value(locy)+y0
1484 locy=lynl+nodplc(loc+11)
1485 value(locy)=-y0
1486 locy=lynl+nodplc(loc+12)
1487 value(locy)=-1.0d0
1488 locy=lynl+nodplc(loc+13)
1489 value(locy)=value(locy)+y0
1490 locy=lynl+nodplc(loc+14)
1491 value(locy)=-1.0d0
1492 locy=lynl+nodplc(loc+15)
1493 value(locy)=-y0
1494 locy=lynl+nodplc(loc+16)
1495 value(locy)=+y0
1496 locy=lynl+nodplc(loc+17)
1497 value(locy)=+1.0d0
1498 locy=lynl+nodplc(loc+18)
1499 value(locy)=+y0
1500 locy=lynl+nodplc(loc+19)
1501 value(locy)=+1.0d0
1502 locy=lynl+nodplc(loc+20)
1503 value(locy)=-1.0d0
1504 locy=lynl+nodplc(loc+23)
1505 value(locy)=+1.0d0
1506 locy=lynl+nodplc(loc+27)
1507 value(locy)=-1.0d0
1508 locy=lynl+nodplc(loc+28)
1509 value(locy)=+1.0d0
1510 locy=lynl+nodplc(loc+31)
1511 value(locy)=-y0
1512 locy=lynl+nodplc(loc+32)
1513 value(locy)=-y0
1514 if (mode.ne.1) go to 920
1515 locy=lynl+nodplc(loc+21)
1516 value(locy)=-1.0d0
1517 locy=lynl+nodplc(loc+22)
1518 value(locy)=+1.0d0
1519 locy=lynl+nodplc(loc+24)
1520 value(locy)=-(1.0d0-gmin)*z0
1521 locy=lynl+nodplc(loc+25)
1522 value(locy)=-1.0d0
1523 locy=lynl+nodplc(loc+26)
1524 value(locy)=+1.0d0
1525 locy=lynl+nodplc(loc+29)
1526 value(locy)=-(1.0d0-gmin)*z0
1527 go to 950
1528 920 if (initf.ne.5) go to 930
1529 if (nosolv.ne.0) go to 925
1530 value(locv+3)=value(lvnim1+node3)-value(lvnim1+node4)
1531 1 +value(lvnim1+ibr2)*z0
1532 value(locv+4)=value(lvnim1+node1)-value(lvnim1+node2)
1533 1 +value(lvnim1+ibr1)*z0
1534 go to 930
1535 925 value(locv+3)=value(locv+7)+value(locv+8)*z0
1536 value(locv+4)=value(locv+5)+value(locv+6)*z0
1537 930 value(lvn+ibr1)=value(locv+3)
1538 value(lvn+ibr2)=value(locv+4)
1539 950 loc=nodplc(loc)
1540 go to 910
1541c
1542c initialize nodes
1543c
1544 980 if(mode.ne.1) go to 995
1545 if(initf.ne.3.and.initf.ne.2) go to 1000
1546 call sizmem(nsnod,nic)
1547 if(nic.eq.0) go to 995
1548 call sizmem(icnod,ntest)
1549 if(modedc.eq.2.and.ntest.ne.0) go to 995
1550 g=1.0d0
1551 do 990 i=1,nic
1552 locy=lynl+nodplc(nsmat+i)
1553 value(locy)=value(locy)+g
1554 node=nodplc(nsnod+i)
1555 value(lvn+node)=value(lvn+node)+value(nsval+i)*g
1556 990 continue
1557c
1558c transient initial conditions (uic not specified)
1559c
1560 995 if(mode.ne.1) go to 1000
1561 if(modedc.ne.2) go to 1000
1562 if(nosolv.ne.0) go to 1000
1563 call sizmem(icnod,nic)
1564 if(nic.eq.0) go to 1000
1565 g=1.0d0
1566 do 996 i=1,nic
1567 locy=lynl+nodplc(icmat+i)
1568 value(locy)=value(locy)+g
1569 node=nodplc(icnod+i)
1570 value(lvn+node)=value(lvn+node)+value(icval+i)*g
1571 996 continue
1572c
1573c reorder right-hand side
1574c
1575 1000 do 1010 i=2,nstop
1576 j=nodplc(iswap+i)
1577 value(ndiag+i)=value(lvn+j)
1578 1010 continue
1579 call copy8(value(ndiag+1),value(lvn+1),nstop)
1580c
1581c finished
1582c
1583 return
1584 end
1585 subroutine nlcsrc
1586 implicit double precision (a-h,o-z)
1587c
1588c this routine loads the nonlinear controlled sources into the
1589c coefficient matrix.
1590c
1591 common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem,
1592 1 isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize,
1593 2 junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr,
1594 3 nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1,
1595 4 lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd,
1596 5 imynl,imvn,lcvn,loutpt,nsnod,nsmat,nsval,icnod,icmat,icval
1597 common /cirdat/ locate(50),jelcnt(50),nunods,ncnods,numnod,nstop,
1598 1 nut,nlt,nxtrm,ndist,ntlin,ibr,numvs
1599 common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet,
1600 1 xmu,mode,modedc,icalc,initf,method,iord,maxord,noncon,iterno,
1601 2 itemno,nosolv,ipostp,iscrch
1602 common /flags/ iprnta,iprntl,iprntm,iprntn,iprnto,limtim,limpts,
1603 1 lvlcod,lvltim,itl1,itl2,itl3,itl4,itl5,igoof,nogo,keof
1604 common /knstnt/ twopi,xlog2,xlog10,root2,rad,boltz,charge,ctok,
1605 1 gmin,reltol,abstol,vntol,trtol,chgtol,eps0,epssil,epsox
1606 common /blank/ value(1000)
1607 integer nodplc(64)
1608 complex*16 cvalue(32)
1609 equivalence (value(1),nodplc(1),cvalue(1))
1610c
1611c nonlinear voltage-controlled current sources
1612c
1613 loc=locate(5)
1614 10 if (loc.eq.0) go to 100
1615 node1=nodplc(loc+2)
1616 node2=nodplc(loc+3)
1617 ndim=nodplc(loc+4)
1618 lnod=nodplc(loc+6)
1619 lmat=nodplc(loc+7)
1620 lcoef=nodplc(loc+8)
1621 call sizmem(nodplc(loc+8),ncoef)
1622 larg=nodplc(loc+9)
1623 lexp=nodplc(loc+10)
1624 lic=nodplc(loc+11)
1625 loct=nodplc(loc+12)+1
1626 icheck=0
1627 do 20 i=1,ndim
1628 call update(value(lic+i),loct,nodplc(lnod+1),nodplc(lnod+2),2,
1629 1 icheck)
1630 value(larg+i)=value(lx0+loct)
1631 loct=loct+2
1632 lnod=lnod+2
1633 20 continue
1634 call evpoly(cold,0,lcoef,ncoef,larg,ndim,lexp)
1635 loct=nodplc(loc+12)
1636 if (icheck.eq.1) go to 30
1637 if (initf.eq.6) go to 30
1638 tol=reltol*dmax1(dabs(cold),dabs(value(lx0+loct)))+abstol
1639 if (dabs(cold-value(lx0+loct)).lt.tol) go to 40
1640 30 noncon=noncon+1
1641 40 value(lx0+loct)=cold
1642 ceq=cold
1643 do 50 i=1,ndim
1644 call evpoly(geq,i,lcoef,ncoef,larg,ndim,lexp)
1645 loct=loct+2
1646 value(lx0+loct)=geq
1647 ceq=ceq-geq*value(larg+i)
1648 locy=lynl+nodplc(lmat+1)
1649 value(locy)=value(locy)+geq
1650 locy=lynl+nodplc(lmat+2)
1651 value(locy)=value(locy)-geq
1652 locy=lynl+nodplc(lmat+3)
1653 value(locy)=value(locy)-geq
1654 locy=lynl+nodplc(lmat+4)
1655 value(locy)=value(locy)+geq
1656 lmat=lmat+4
1657 50 continue
1658 value(lvn+node1)=value(lvn+node1)-ceq
1659 value(lvn+node2)=value(lvn+node2)+ceq
1660 loc=nodplc(loc)
1661 go to 10
1662c
1663c nonlinear voltage controlled voltage sources
1664c
1665 100 loc=locate(6)
1666 110 if (loc.eq.0) go to 200
1667 node1=nodplc(loc+2)
1668 node2=nodplc(loc+3)
1669 ndim=nodplc(loc+4)
1670 iptr=nodplc(loc+6)
1671 lnod=nodplc(loc+7)
1672 lmat=nodplc(loc+8)
1673 lcoef=nodplc(loc+9)
1674 call sizmem(nodplc(loc+9),ncoef)
1675 larg=nodplc(loc+10)
1676 lexp=nodplc(loc+11)
1677 lic=nodplc(loc+12)
1678 loct=nodplc(loc+13)+2
1679 icheck=0
1680 do 120 i=1,ndim
1681 call update(value(lic+i),loct,nodplc(lnod+1),nodplc(lnod+2),2,
1682 1 icheck)
1683 value(larg+i)=value(lx0+loct)
1684 loct=loct+2
1685 lnod=lnod+2
1686 120 continue
1687 call evpoly(volt,0,lcoef,ncoef,larg,ndim,lexp)
1688 loct=nodplc(loc+13)
1689 if (icheck.eq.1) go to 130
1690 if (initf.eq.6) go to 130
1691 tol=reltol*dmax1(dabs(volt),dabs(value(lx0+loct)))+vntol
1692 if (dabs(volt-value(lx0+loct)).lt.tol) go to 140
1693 130 noncon=noncon+1
1694 140 value(lx0+loct)=volt
1695 value(lx0+loct+1)=value(lvnim1+iptr)
1696 veq=volt
1697 locy=lynl+nodplc(lmat+1)
1698 value(locy)=+1.0d0
1699 locy=lynl+nodplc(lmat+2)
1700 value(locy)=-1.0d0
1701 locy=lynl+nodplc(lmat+3)
1702 value(locy)=+1.0d0
1703 locy=lynl+nodplc(lmat+4)
1704 value(locy)=-1.0d0
1705 lmat=lmat+4
1706 loct=loct+1
1707 do 150 i=1,ndim
1708 call evpoly(vgain,i,lcoef,ncoef,larg,ndim,lexp)
1709 loct=loct+2
1710 value(lx0+loct)=vgain
1711 veq=veq-vgain*value(larg+i)
1712 locy=lynl+nodplc(lmat+1)
1713 value(locy)=value(locy)-vgain
1714 locy=lynl+nodplc(lmat+2)
1715 value(locy)=value(locy)+vgain
1716 lmat=lmat+2
1717 150 continue
1718 value(lvn+iptr)=veq
1719 loc=nodplc(loc)
1720 go to 110
1721c
1722c nonlinear current-controlled current sources
1723c
1724 200 loc=locate(7)
1725 210 if (loc.eq.0) go to 300
1726 node1=nodplc(loc+2)
1727 node2=nodplc(loc+3)
1728 ndim=nodplc(loc+4)
1729 lvs=nodplc(loc+6)
1730 lmat=nodplc(loc+7)
1731 lcoef=nodplc(loc+8)
1732 call sizmem(nodplc(loc+8),ncoef)
1733 larg=nodplc(loc+9)
1734 lexp=nodplc(loc+10)
1735 lic=nodplc(loc+11)
1736 loct=nodplc(loc+12)+1
1737 icheck=0
1738 do 220 i=1,ndim
1739 iptr=nodplc(lvs+i)
1740 iptr=nodplc(iptr+6)
1741 call update(value(lic+i),loct,iptr,1,2,icheck)
1742 value(larg+i)=value(lx0+loct)
1743 loct=loct+2
1744 220 continue
1745 call evpoly(csrc,0,lcoef,ncoef,larg,ndim,lexp)
1746 loct=nodplc(loc+12)
1747 if (icheck.eq.1) go to 230
1748 if (initf.eq.6) go to 230
1749 tol=reltol*dmax1(dabs(csrc),dabs(value(lx0+loct)))+abstol
1750 if (dabs(csrc-value(lx0+loct)).lt.tol) go to 240
1751 230 noncon=noncon+1
1752 240 value(lx0+loct)=csrc
1753 ceq=csrc
1754 do 250 i=1,ndim
1755 call evpoly(cgain,i,lcoef,ncoef,larg,ndim,lexp)
1756 loct=loct+2
1757 value(lx0+loct)=cgain
1758 ceq=ceq-cgain*value(larg+i)
1759 locy=lynl+nodplc(lmat+1)
1760 value(locy)=value(locy)+cgain
1761 locy=lynl+nodplc(lmat+2)
1762 value(locy)=value(locy)-cgain
1763 lmat=lmat+2
1764 250 continue
1765 value(lvn+node1)=value(lvn+node1)-ceq
1766 value(lvn+node2)=value(lvn+node2)+ceq
1767 loc=nodplc(loc)
1768 go to 210
1769c
1770c nonlinear current controlled voltage sources
1771c
1772 300 loc=locate(8)
1773 310 if (loc.eq.0) go to 1000
1774 node1=nodplc(loc+2)
1775 node2=nodplc(loc+3)
1776 ndim=nodplc(loc+4)
1777 ibr=nodplc(loc+6)
1778 lvs=nodplc(loc+7)
1779 lmat=nodplc(loc+8)
1780 lcoef=nodplc(loc+9)
1781 call sizmem(nodplc(loc+9),ncoef)
1782 larg=nodplc(loc+10)
1783 lexp=nodplc(loc+11)
1784 lic=nodplc(loc+12)
1785 loct=nodplc(loc+13)+2
1786 icheck=0
1787 do 320 i=1,ndim
1788 iptr=nodplc(lvs+i)
1789 iptr=nodplc(iptr+6)
1790 call update(value(lic+i),loct,iptr,1,2,icheck)
1791 value(larg+i)=value(lx0+loct)
1792 loct=loct+2
1793 320 continue
1794 call evpoly(volt,0,lcoef,ncoef,larg,ndim,lexp)
1795 loct=nodplc(loc+13)
1796 if (icheck.eq.1) go to 330
1797 if (initf.eq.6) go to 330
1798 tol=reltol*dmax1(dabs(volt),dabs(value(lx0+loct)))+vntol
1799 if (dabs(volt-value(lx0+loct)).lt.tol) go to 340
1800 330 noncon=noncon+1
1801 340 value(lx0+loct)=volt
1802 value(lx0+loct+1)=value(lvnim1+ibr)
1803 veq=volt
1804 locy=lynl+nodplc(lmat+1)
1805 value(locy)=+1.0d0
1806 locy=lynl+nodplc(lmat+2)
1807 value(locy)=-1.0d0
1808 locy=lynl+nodplc(lmat+3)
1809 value(locy)=+1.0d0
1810 locy=lynl+nodplc(lmat+4)
1811 value(locy)=-1.0d0
1812 lmat=lmat+4
1813 loct=loct+1
1814 do 350 i=1,ndim
1815 call evpoly(transr,i,lcoef,ncoef,larg,ndim,lexp)
1816 loct=loct+2
1817 value(lx0+loct)=transr
1818 veq=veq-transr*value(larg+i)
1819 locy=lynl+nodplc(lmat+i)
1820 value(locy)=value(locy)-transr
1821 350 continue
1822 value(lvn+ibr)=veq
1823 loc=nodplc(loc)
1824 go to 310
1825c
1826c finished
1827c
1828 1000 return
1829 end
1830 subroutine update(vinit,loct,node1,node2,nupda,icheck)
1831 implicit double precision (a-h,o-z)
1832c
1833c this routine updates and limits the controlling variables for the
1834c nonlinear controlled sources.
1835c
1836 common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem,
1837 1 isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize,
1838 2 junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr,
1839 3 nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1,
1840 4 lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd,
1841 5 imynl,imvn,lcvn,loutpt,nsnod,nsmat,nsval,icnod,icmat,icval
1842 common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet,
1843 1 xmu,mode,modedc,icalc,initf,method,iord,maxord,noncon,iterno,
1844 2 itemno,nosolv,ipostp,iscrch
1845 common /blank/ value(1000)
1846 integer nodplc(64)
1847 complex*16 cvalue(32)
1848 equivalence (value(1),nodplc(1),cvalue(1))
1849c
1850c
1851 go to (40,10,40,20,30,50), initf
1852 10 vnew=vinit
1853 go to 70
1854 20 vnew=value(lx0+loct)
1855 go to 70
1856 30 vnew=value(lx1+loct)
1857 go to 70
1858 40 vnew=value(lvnim1+node1)-value(lvnim1+node2)
1859 go to 60
1860 50 call copy8(value(lx1+loct),value(lx0+loct),nupda)
1861 xfact=delta/delold(2)
1862 vnew=(1.0d0+xfact)*value(lx1+loct)-xfact*value(lx2+loct)
1863 60 if (dabs(vnew).le.1.0d0) go to 80
1864 delv=vnew-value(lx0+loct)
1865 if (dabs(delv).le.0.1d0) go to 80
1866 vlim=dmax1(dabs(0.1d0*value(lx0+loct)),0.1d0)
1867 vnew=value(lx0+loct)+dsign(dmin1(dabs(delv),vlim),delv)
1868 go to 70
1869 70 icheck=1
1870 80 value(lx0+loct)=vnew
1871 return
1872 end
1873 subroutine evpoly(result,itype,lcoef,ncoef,larg,
1874 1 narg,lexp)
1875 implicit double precision (a-h,o-z)
1876c
1877c this routine evaluates a polynomial. lcoef points to the coef-
1878c ficients, and larg points to the values of the polynomial argument(s).
1879c
1880 common /blank/ value(1000)
1881 integer nodplc(64)
1882 complex*16 cvalue(32)
1883 equivalence (value(1),nodplc(1),cvalue(1))
1884c
1885c
1886 if (itype) 100,200,300
1887c
1888c integration (polynomial *must* be one-dimensional)
1889c
1890 100 result=0.0d0
1891 arg=1.0d0
1892 arg1=value(larg+1)
1893 do 110 i=1,ncoef
1894 arg=arg*arg1
1895 result=result+value(lcoef+i)*arg/dfloat(i)
1896 110 continue
1897 go to 1000
1898c
1899c evaluation of the polynomial
1900c
1901 200 result=value(lcoef+1)
1902 if (ncoef.eq.1) go to 1000
1903 call zero4(nodplc(lexp+1),narg)
1904 do 220 i=2,ncoef
1905 call nxtpwr(nodplc(lexp+1),narg)
1906 if (value(lcoef+i).eq.0.0d0) go to 220
1907 arg=1.0d0
1908 do 210 j=1,narg
1909 call evterm(val,value(larg+j),nodplc(lexp+j))
1910 arg=arg*val
1911 210 continue
1912 result=result+value(lcoef+i)*arg
1913 220 continue
1914 go to 1000
1915c
1916c partial derivative with respect to the itype*th variable
1917c
1918 300 result=0.0d0
1919 if (ncoef.eq.1) go to 1000
1920 call zero4(nodplc(lexp+1),narg)
1921 do 330 i=2,ncoef
1922 call nxtpwr(nodplc(lexp+1),narg)
1923 if (nodplc(lexp+itype).eq.0) go to 330
1924 if (value(lcoef+i).eq.0.0d0) go to 330
1925 arg=1.0d0
1926 do 320 j=1,narg
1927 if (j.eq.itype) go to 310
1928 call evterm(val,value(larg+j),nodplc(lexp+j))
1929 arg=arg*val
1930 go to 320
1931 310 call evterm(val,value(larg+j),nodplc(lexp+j)-1)
1932 arg=arg*dfloat(nodplc(lexp+j))*val
1933 320 continue
1934 result=result+value(lcoef+i)*arg
1935 330 continue
1936c
1937c finished
1938c
1939 1000 return
1940 end
1941 subroutine evterm(val,arg,iexp)
1942 implicit double precision (a-h,o-z)
1943c
1944c this routine evaluates one term of a polynomial.
1945c
1946 jexp=iexp+1
1947 if (jexp.ge.6) go to 60
1948 go to (10,20,30,40,50), jexp
1949 10 val=1.0d0
1950 go to 100
1951 20 val=arg
1952 go to 100
1953 30 val=arg*arg
1954 go to 100
1955 40 val=arg*arg*arg
1956 go to 100
1957 50 val=arg*arg
1958 val=val*val
1959 go to 100
1960 60 if (arg.eq.0.0d0) go to 70
1961 argexp=dfloat(iexp)*dlog(dabs(arg))
1962 if (argexp.lt.-200.0d0) go to 70
1963 val=dexp(argexp)
1964 if((iexp/2)*2.eq.iexp) go to 100
1965 val=dsign(val,arg)
1966 go to 100
1967 70 val=0.0d0
1968c
1969c finished
1970c
1971 100 return
1972 end
1973 subroutine nxtpwr(pwrseq,pdim)
1974 implicit double precision (a-h,o-z)
1975c
1976c this routine determines the 'next' set of exponents for the
1977c different dimensions of a polynomial.
1978c
1979 integer pwrseq(1),pdim,psum
1980c
1981c
1982 if (pdim.eq.1) go to 80
1983 k=pdim
1984 10 if (pwrseq(k).ne.0) go to 20
1985 k=k-1
1986 if (k.ne.0) go to 10
1987 go to 80
1988 20 if (k.eq.pdim) go to 30
1989 pwrseq(k)=pwrseq(k)-1
1990 pwrseq(k+1)=pwrseq(k+1)+1
1991 go to 100
1992 30 km1=k-1
1993 do 40 i=1,km1
1994 if (pwrseq(i).ne.0) go to 50
1995 40 continue
1996 pwrseq(1)=pwrseq(pdim)+1
1997 pwrseq(pdim)=0
1998 go to 100
1999 50 psum=1
2000 k=pdim
2001 60 if (pwrseq(k-1).ge.1) go to 70
2002 psum=psum+pwrseq(k)
2003 pwrseq(k)=0
2004 k=k-1
2005 go to 60
2006 70 pwrseq(k)=pwrseq(k)+psum
2007 pwrseq(k-1)=pwrseq(k-1)-1
2008 go to 100
2009 80 pwrseq(1)=pwrseq(1)+1
2010c
2011c finished
2012c
2013 100 return
2014 end
2015 subroutine intgr8(geq,ceq,capval,loct)
2016 implicit double precision (a-h,o-z)
2017c
2018c this routine performs the actual numerical integration for each
2019c circuit element.
2020c
2021 common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem,
2022 1 isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize,
2023 2 junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr,
2024 3 nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1,
2025 4 lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd,
2026 5 imynl,imvn,lcvn,loutpt,nsnod,nsmat,nsval,icnod,icmat,icval
2027 common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet,
2028 1 xmu,mode,modedc,icalc,initf,method,iord,maxord,noncon,iterno,
2029 2 itemno,nosolv,ipostp,iscrch
2030 common /blank/ value(1000)
2031 integer nodplc(64)
2032 complex*16 cvalue(32)
2033 equivalence (value(1),nodplc(1),cvalue(1))
2034c
2035c
2036 dimension qcap(1),ccap(1)
2037 equivalence (qcap(1),value(1)),(ccap(1),value(2))
2038c
2039c
2040 if (method.eq.2) go to 100
2041c
2042c trapezoidal algorithm
2043c
2044 if (iord.eq.1) go to 100
2045 ccap(lx0+loct)=-ccap(lx1+loct)*ag(2)
2046 1 +ag(1)*(qcap(lx0+loct)-qcap(lx1+loct))
2047 go to 190
2048c
2049c gears algorithm
2050c
2051 100 go to (110,120,130,140,150,160), iord
2052 110 ccap(lx0+loct)=ag(1)*qcap(lx0+loct)+ag(2)*qcap(lx1+loct)
2053 go to 190
2054 120 ccap(lx0+loct)=ag(1)*qcap(lx0+loct)+ag(2)*qcap(lx1+loct)
2055 1 +ag(3)*qcap(lx2+loct)
2056 go to 190
2057 130 ccap(lx0+loct)=ag(1)*qcap(lx0+loct)+ag(2)*qcap(lx1+loct)
2058 1 +ag(3)*qcap(lx2+loct)+ag(4)*qcap(lx3+loct)
2059 go to 190
2060 140 ccap(lx0+loct)=ag(1)*qcap(lx0+loct)+ag(2)*qcap(lx1+loct)
2061 1 +ag(3)*qcap(lx2+loct)+ag(4)*qcap(lx3+loct)
2062 2 +ag(5)*qcap(lx4+loct)
2063 go to 190
2064 150 ccap(lx0+loct)=ag(1)*qcap(lx0+loct)+ag(2)*qcap(lx1+loct)
2065 1 +ag(3)*qcap(lx2+loct)+ag(4)*qcap(lx3+loct)
2066 2 +ag(5)*qcap(lx4+loct)+ag(6)*qcap(lx5+loct)
2067 go to 190
2068 160 ccap(lx0+loct)=ag(1)*qcap(lx0+loct)+ag(2)*qcap(lx1+loct)
2069 1 +ag(3)*qcap(lx2+loct)+ag(4)*qcap(lx3+loct)
2070 2 +ag(5)*qcap(lx4+loct)+ag(6)*qcap(lx5+loct)
2071 3 +ag(7)*qcap(lx6+loct)
2072c... ceq is the equivalent current applicable to linear capacitance
2073c (inductance) only, i.e. q=c*v
2074 190 ceq=ccap(lx0+loct)-ag(1)*qcap(lx0+loct)
2075 geq=ag(1)*capval
2076 return
2077 end
2078 subroutine pnjlim(vnew,vold,vt,vcrit,icheck)
2079 implicit double precision (a-h,o-z)
2080c
2081c this routine limits the change-per-iteration of device pn-junction
2082c voltages.
2083c
2084 if (vnew.le.vcrit) go to 30
2085 vlim=vt+vt
2086 delv=vnew-vold
2087 if (dabs(delv).le.vlim) go to 30
2088 if (vold.le.0.0d0) go to 20
2089 arg=1.0d0+delv/vt
2090 if (arg.le.0.0d0) go to 10
2091 vnew=vold+vt*dlog(arg)
2092 icheck=1
2093 go to 100
2094 10 vnew=vcrit
2095 icheck=1
2096 go to 100
2097 20 vnew=vt*dlog(vnew/vt)
2098 icheck=1
2099 go to 100
2100 30 continue
2101c
2102c finished
2103c
2104 100 return
2105 end
2106 subroutine diode
2107 implicit double precision (a-h,o-z)
2108c
2109c this routine processes diodes for dc and transient analyses.
2110c
2111 common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem,
2112 1 isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize,
2113 2 junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr,
2114 3 nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1,
2115 4 lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd,
2116 5 imynl,imvn,lcvn,loutpt,nsnod,nsmat,nsval,icnod,icmat,icval
2117 common /cirdat/ locate(50),jelcnt(50),nunods,ncnods,numnod,nstop,
2118 1 nut,nlt,nxtrm,ndist,ntlin,ibr,numvs
2119 common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet,
2120 1 xmu,mode,modedc,icalc,initf,method,iord,maxord,noncon,iterno,
2121 2 itemno,nosolv,ipostp,iscrch
2122 common /knstnt/ twopi,xlog2,xlog10,root2,rad,boltz,charge,ctok,
2123 1 gmin,reltol,abstol,vntol,trtol,chgtol,eps0,epssil,epsox
2124 common /blank/ value(1000)
2125 integer nodplc(64)
2126 complex*16 cvalue(32)
2127 equivalence (value(1),nodplc(1),cvalue(1))
2128c
2129c
2130 dimension vdo(1),cdo(1),gdo(1),qd(1),cqd(1)
2131 equivalence (vdo(1),value(1)),(cdo(1),value(2)),
2132 1 (gdo(1),value(3)),(qd(1),value(4)),(cqd(1),value(5))
2133c
2134c
2135 loc=locate(11)
2136 10 if (loc.eq.0) return
2137 locv=nodplc(loc+1)
2138 node1=nodplc(loc+2)
2139 node2=nodplc(loc+3)
2140 node3=nodplc(loc+4)
2141 locm=nodplc(loc+5)
2142 ioff=nodplc(loc+6)
2143 locm=nodplc(locm+1)
2144 loct=nodplc(loc+11)
2145c
2146c dc model parameters
2147c
2148 area=value(locv+1)
2149 csat=value(locm+1)*area
2150 gspr=value(locm+2)*area
2151 vte=value(locm+3)*vt
2152 bv=value(locm+13)
2153 vcrit=value(locm+18)
2154c
2155c initialization
2156c
2157 icheck=1
2158 go to (100,20,30,50,60,70),initf
2159 20 if(mode.ne.1.or.modedc.ne.2.or.nosolv.eq.0) go to 25
2160 vd=value(locv+2)
2161 go to 300
2162 25 if(ioff.ne.0) go to 40
2163 vd=vcrit
2164 go to 300
2165 30 if (ioff.eq.0) go to 100
2166 40 vd=0.0d0
2167 go to 300
2168 50 vd=vdo(lx0+loct)
2169 go to 300
2170 60 vd=vdo(lx1+loct)
2171 go to 300
2172 70 xfact=delta/delold(2)
2173 vdo(lx0+loct)=vdo(lx1+loct)
2174 vd=(1.0d0+xfact)*vdo(lx1+loct)-xfact*vdo(lx2+loct)
2175 cdo(lx0+loct)=cdo(lx1+loct)
2176 gdo(lx0+loct)=gdo(lx1+loct)
2177 go to 110
2178c
2179c compute new nonlinear branch voltage
2180c
2181 100 vd=value(lvnim1+node3)-value(lvnim1+node2)
2182 110 delvd=vd-vdo(lx0+loct)
2183 cdhat=cdo(lx0+loct)+gdo(lx0+loct)*delvd
2184c
2185c bypass if solution has not changed
2186c
2187 if (6 .eq.6) go to 200
2188 tol=reltol*dmax1(dabs(vd),dabs(vdo(lx0+loct)))+vntol
2189 if (dabs(delvd).ge.tol) go to 200
2190 tol=reltol*dmax1(dabs(cdhat),dabs(cdo(lx0+loct)))+abstol
2191 if (dabs(cdhat-cdo(lx0+loct)).ge.tol) go to 200
2192 vd=vdo(lx0+loct)
2193 cd=cdo(lx0+loct)
2194 gd=gdo(lx0+loct)
2195 go to 800
2196c
2197c limit new junction voltage
2198c
2199 200 vlim=vte+vte
2200 if(bv.eq.0.0d0) go to 205
2201 if (vd.lt.dmin1(0.0d0,-bv+10.0d0*vte)) go to 210
2202 205 icheck=0
2203 call pnjlim(vd,vdo(lx0+loct),vte,vcrit,icheck)
2204 go to 300
2205 210 vdtemp=-(vd+bv)
2206 call pnjlim(vdtemp,-(vdo(lx0+loct)+bv),vte,vcrit,icheck)
2207 vd=-(vdtemp+bv)
2208c
2209c compute dc current and derivitives
2210c
2211 300 if (vd.lt.-5.0d0*vte) go to 310
2212 evd=dexp(vd/vte)
2213 cd=csat*(evd-1.0d0)+gmin*vd
2214 gd=csat*evd/vte+gmin
2215 go to 330
2216 310 if(bv.eq.0.0d0) go to 315
2217 if(vd.lt.-bv) go to 320
2218 315 gd=-csat/vd+gmin
2219 cd=gd*vd
2220 go to 330
2221 320 evrev=dexp(-(bv+vd)/vt)
2222 cd=-csat*(evrev-1.0d0+bv/vt)
2223 gd=csat*evrev/vt
2224 330 if (mode.ne.1) go to 500
2225 if ((modedc.eq.2).and.(nosolv.ne.0)) go to 500
2226 if (initf.eq.4) go to 500
2227 go to 700
2228c
2229c charge storage elements
2230c
2231 500 tau=value(locm+4)
2232 czero=value(locm+5)*area
2233 pb=value(locm+6)
2234 xm=value(locm+7)
2235 fcpb=value(locm+12)
2236 if (vd.ge.fcpb) go to 510
2237 arg=1.0d0-vd/pb
2238 sarg=dexp(-xm*dlog(arg))
2239 qd(lx0+loct)=tau*cd+pb*czero*(1.0d0-arg*sarg)/(1.0d0-xm)
2240 capd=tau*gd+czero*sarg
2241 go to 520
2242 510 f1=value(locm+15)
2243 f2=value(locm+16)
2244 f3=value(locm+17)
2245 czof2=czero/f2
2246 qd(lx0+loct)=tau*cd+czero*f1+czof2*(f3*(vd-fcpb)
2247 1 +(xm/(pb+pb))*(vd*vd-fcpb*fcpb))
2248 capd=tau*gd+czof2*(f3+xm*vd/pb)
2249c
2250c store small-signal parameters
2251c
2252 520 if ((mode.eq.1).and.(modedc.eq.2).and.(nosolv.ne.0)) go to 700
2253 if (initf.ne.4) go to 600
2254 value(lx0+loct+4)=capd
2255 go to 1000
2256c
2257c transient analysis
2258c
2259 600 if (initf.ne.5) go to 610
2260 qd(lx1+loct)=qd(lx0+loct)
2261 610 call intgr8(geq,ceq,capd,loct+3)
2262 gd=gd+geq
2263 cd=cd+cqd(lx0+loct)
2264 if (initf.ne.5) go to 700
2265 cqd(lx1+loct)=cqd(lx0+loct)
2266c
2267c check convergence
2268c
2269 700 if (initf.ne.3) go to 710
2270 if (ioff.eq.0) go to 710
2271 go to 750
2272 710 if (icheck.eq.1) go to 720
2273 tol=reltol*dmax1(dabs(cdhat),dabs(cd))+abstol
2274 if (dabs(cdhat-cd).le.tol) go to 750
2275 720 noncon=noncon+1
2276 750 vdo(lx0+loct)=vd
2277 cdo(lx0+loct)=cd
2278 gdo(lx0+loct)=gd
2279c
2280c load current vector
2281c
2282 800 cdeq=cd-gd*vd
2283 value(lvn+node2)=value(lvn+node2)+cdeq
2284 value(lvn+node3)=value(lvn+node3)-cdeq
2285c
2286c load matrix
2287c
2288 locy=lynl+nodplc(loc+13)
2289 value(locy)=value(locy)+gspr
2290 locy=lynl+nodplc(loc+14)
2291 value(locy)=value(locy)+gd
2292 locy=lynl+nodplc(loc+15)
2293 value(locy)=value(locy)+gd+gspr
2294 locy=lynl+nodplc(loc+7)
2295 value(locy)=value(locy)-gspr
2296 locy=lynl+nodplc(loc+8)
2297 value(locy)=value(locy)-gd
2298 locy=lynl+nodplc(loc+9)
2299 value(locy)=value(locy)-gspr
2300 locy=lynl+nodplc(loc+10)
2301 value(locy)=value(locy)-gd
2302 1000 loc=nodplc(loc)
2303 go to 10
2304 end
2305 subroutine bjt
2306 implicit double precision (a-h,o-z)
2307c
2308c this routine processes bjts for dc and transient analyses.
2309c
2310 common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem,
2311 1 isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize,
2312 2 junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr,
2313 3 nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1,
2314 4 lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd,
2315 5 imynl,imvn,lcvn,loutpt,nsnod,nsmat,nsval,icnod,icmat,icval
2316 common /cirdat/ locate(50),jelcnt(50),nunods,ncnods,numnod,nstop,
2317 1 nut,nlt,nxtrm,ndist,ntlin,ibr,numvs
2318 common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet,
2319 1 xmu,mode,modedc,icalc,initf,method,iord,maxord,noncon,iterno,
2320 2 itemno,nosolv,ipostp,iscrch
2321 common /knstnt/ twopi,xlog2,xlog10,root2,rad,boltz,charge,ctok,
2322 1 gmin,reltol,abstol,vntol,trtol,chgtol,eps0,epssil,epsox
2323 common /blank/ value(1000)
2324 integer nodplc(64)
2325 complex*16 cvalue(32)
2326 equivalence (value(1),nodplc(1),cvalue(1))
2327c
2328c
2329 dimension vbeo(1),vbco(1),cco(1),cbo(1),gpio(1),gmuo(1),gmo(1),
2330 1 goo(1),qbe(1),cqbe(1),qbc(1),cqbc(1),qcs(1),cqcs(1),qbx(1),
2331 2 cqbx(1),gxo(1),cexbc(1)
2332 equivalence (vbeo(1),value(1)),(vbco(1),value(2)),
2333 1 (cco(1),value(3)),(cbo(1),value(4)),(gpio(1),value(5)),
2334 2 (gmuo(1),value(6)),(gmo(1),value(7)),(goo(1),value(8)),
2335 3 (qbe(1),value(9)),(cqbe(1),value(10)),(qbc(1),value(11)),
2336 4 (cqbc(1),value(12)),(qcs(1),value(13)),(cqcs(1),value(14)),
2337 5 (qbx(1),value(15)),(cqbx(1),value(16)),(gxo(1),value(17)),
2338 6 (cexbc(1),value(18))
2339c
2340c
2341 loc=locate(12)
2342 10 if (loc.eq.0) return
2343 locv=nodplc(loc+1)
2344 node1=nodplc(loc+2)
2345 node2=nodplc(loc+3)
2346 node3=nodplc(loc+4)
2347 node4=nodplc(loc+5)
2348 node5=nodplc(loc+6)
2349 node6=nodplc(loc+7)
2350 node7=nodplc(loc+30)
2351 locm=nodplc(loc+8)
2352 ioff=nodplc(loc+9)
2353 type=nodplc(locm+2)
2354 locm=nodplc(locm+1)
2355 loct=nodplc(loc+22)
2356 gccs=0.0d0
2357 ceqcs=0.0d0
2358 geqbx=0.0d0
2359 ceqbx=0.0d0
2360 geqcb=0.0d0
2361c
2362c dc model paramters
2363c
2364 area=value(locv+1)
2365 bfm=value(locm+2)
2366 brm=value(locm+8)
2367 csat=value(locm+1)*area
2368 rbpr=value(locm+18)/area
2369 rbpi=value(locm+16)/area-rbpr
2370 gcpr=value(locm+20)*area
2371 gepr=value(locm+19)*area
2372 ova=value(locm+4)
2373 ovb=value(locm+10)
2374 oik=value(locm+5)/area
2375 c2=value(locm+6)*area
2376 vte=value(locm+7)*vt
2377 oikr=value(locm+11)/area
2378 c4=value(locm+12)*area
2379 vtc=value(locm+13)*vt
2380 vcrit=value(locm+54)
2381 td=value(locm+28)
2382 xjrb=value(locm+17)*area
2383c
2384c initialization
2385c
2386 icheck=1
2387 go to (100,20,30,50,60,70),initf
2388 20 if(mode.ne.1.or.modedc.ne.2.or.nosolv.eq.0) go to 25
2389 vbe=type*value(locv+2)
2390 vce=type*value(locv+3)
2391 vbc=vbe-vce
2392 vbx=vbc
2393 vcs=0.0d0
2394 go to 300
2395 25 if(ioff.ne.0) go to 40
2396 vbe=vcrit
2397 vbc=0.0d0
2398 go to 300
2399 30 if (ioff.eq.0) go to 100
2400 40 vbe=0.0d0
2401 vbc=0.0d0
2402 go to 300
2403 50 vbe=vbeo(lx0+loct)
2404 vbc=vbco(lx0+loct)
2405 vbx=type*(value(lvnim1+node2)-value(lvnim1+node4))
2406 vcs=type*(value(lvnim1+node7)-value(lvnim1+node4))
2407 go to 300
2408 60 vbe=vbeo(lx1+loct)
2409 vbc=vbco(lx1+loct)
2410 vbx=type*(value(lvnim1+node2)-value(lvnim1+node4))
2411 vcs=type*(value(lvnim1+node7)-value(lvnim1+node4))
2412 if(mode.ne.2.or.nosolv.eq.0) go to 300
2413 vbx=type*(value(locv+2)-value(locv+3))
2414 vcs=0.0d0
2415 go to 300
2416 70 xfact=delta/delold(2)
2417 vbeo(lx0+loct)=vbeo(lx1+loct)
2418 vbe=(1.0d0+xfact)*vbeo(lx1+loct)-xfact*vbeo(lx2+loct)
2419 vbco(lx0+loct)=vbco(lx1+loct)
2420 vbc=(1.0d0+xfact)*vbco(lx1+loct)-xfact*vbco(lx2+loct)
2421 cco(lx0+loct)=cco(lx1+loct)
2422 cbo(lx0+loct)=cbo(lx1+loct)
2423 gpio(lx0+loct)=gpio(lx1+loct)
2424 gmuo(lx0+loct)=gmuo(lx1+loct)
2425 gmo(lx0+loct)=gmo(lx1+loct)
2426 goo(lx0+loct)=goo(lx1+loct)
2427 go to 110
2428c
2429c compute new nonlinear branch voltages
2430c
2431 100 vbe=type*(value(lvnim1+node5)-value(lvnim1+node6))
2432 vbc=type*(value(lvnim1+node5)-value(lvnim1+node4))
2433 110 delvbe=vbe-vbeo(lx0+loct)
2434 delvbc=vbc-vbco(lx0+loct)
2435 vbx=type*(value(lvnim1+node2)-value(lvnim1+node4))
2436 vcs=type*(value(lvnim1+node7)-value(lvnim1+node4))
2437 cchat=cco(lx0+loct)+(gmo(lx0+loct)+goo(lx0+loct))*delvbe
2438 1 -(goo(lx0+loct)+gmuo(lx0+loct))*delvbc
2439 cbhat=cbo(lx0+loct)+gpio(lx0+loct)*delvbe+gmuo(lx0+loct)*delvbc
2440c
2441c limit nonlinear branch voltages
2442c
2443 200 icheck=0
2444 call pnjlim(vbe,vbeo(lx0+loct),vt,vcrit,icheck)
2445 call pnjlim(vbc,vbco(lx0+loct),vt,vcrit,icheck)
2446c
2447c determine dc current and derivitives
2448c
2449 300 vtn=vt*value(locm+3)
2450 if(vbe.le.-5.0d0*vtn) go to 320
2451 evbe=dexp(vbe/vtn)
2452 cbe=csat*(evbe-1.0d0)+gmin*vbe
2453 gbe=csat*evbe/vtn+gmin
2454 if (c2.ne.0.0d0) go to 310
2455 cben=0.0d0
2456 gben=0.0d0
2457 go to 350
2458 310 evben=dexp(vbe/vte)
2459 cben=c2*(evben-1.0d0)
2460 gben=c2*evben/vte
2461 go to 350
2462 320 gbe=-csat/vbe+gmin
2463 cbe=gbe*vbe
2464 gben=-c2/vbe
2465 cben=gben*vbe
2466 350 vtn=vt*value(locm+9)
2467 if(vbc.le.-5.0d0*vtn) go to 370
2468 evbc=dexp(vbc/vtn)
2469 cbc=csat*(evbc-1.0d0)+gmin*vbc
2470 gbc=csat*evbc/vtn+gmin
2471 if (c4.ne.0.0d0) go to 360
2472 cbcn=0.0d0
2473 gbcn=0.0d0
2474 go to 400
2475 360 evbcn=dexp(vbc/vtc)
2476 cbcn=c4*(evbcn-1.0d0)
2477 gbcn=c4*evbcn/vtc
2478 go to 400
2479 370 gbc=-csat/vbc+gmin
2480 cbc=gbc*vbc
2481 gbcn=-c4/vbc
2482 cbcn=gbcn*vbc
2483c
2484c determine base charge terms
2485c
2486 400 q1=1.0d0/(1.0d0-ova*vbc-ovb*vbe)
2487 if (oik.ne.0.0d0) go to 405
2488 if (oikr.ne.0.0d0) go to 405
2489 qb=q1
2490 dqbdve=q1*qb*ovb
2491 dqbdvc=q1*qb*ova
2492 go to 410
2493 405 q2=oik*cbe+oikr*cbc
2494 arg=dmax1(0.0d0,1.0d0+4.0d0*q2)
2495 sqarg=1.0d0
2496 if(arg.ne.0.0d0) sqarg=dsqrt(arg)
2497 qb=q1*(1.0d0+sqarg)/2.0d0
2498 dqbdve=q1*(qb*ovb+oik*gbe/sqarg)
2499 dqbdvc=q1*(qb*ova+oikr*gbc/sqarg)
2500c
2501c weil's approx. for excess phase applied with backward-
2502c euler integration
2503c
2504 410 cc=0.0d0
2505 cex=cbe
2506 gex=gbe
2507 if(mode.eq.1) go to 420
2508 if(td.eq.0.0d0) go to 420
2509 arg1=delta/td
2510 arg2=3.0d0*arg1
2511 arg1=arg2*arg1
2512 denom=1.0d0+arg1+arg2
2513 arg3=arg1/denom
2514 if(initf.ne.5) go to 411
2515 cexbc(lx1+loct)=cbe/qb
2516 cexbc(lx2+loct)=cexbc(lx1+loct)
2517 411 cc=(cexbc(lx1+loct)*(1.0d0+delta/delold(2)+arg2)
2518 1 -cexbc(lx2+loct)*delta/delold(2))/denom
2519 cex=cbe*arg3
2520 gex=gbe*arg3
2521 cexbc(lx0+loct)=cc+cex/qb
2522c
2523c determine dc incremental conductances
2524c
2525 420 cc=cc+(cex-cbc)/qb-cbc/brm-cbcn
2526 cb=cbe/bfm+cben+cbc/brm+cbcn
2527 gx=rbpr+rbpi/qb
2528 if(xjrb.eq.0.0d0) go to 430
2529 arg1=dmax1(cb/xjrb,1.0d-9)
2530 arg2=(-1.0d0+dsqrt(1.0d0+14.59025d0*arg1))/2.4317d0/dsqrt(arg1)
2531 arg1=dtan(arg2)
2532 gx=rbpr+3.0d0*rbpi*(arg1-arg2)/arg2/arg1/arg1
2533 430 if(gx.ne.0.0d0) gx=1.0d0/gx
2534 gpi=gbe/bfm+gben
2535 gmu=gbc/brm+gbcn
2536 go=(gbc+(cex-cbc)*dqbdvc/qb)/qb
2537 gm=(gex-(cex-cbc)*dqbdve/qb)/qb-go
2538 if (mode.ne.1) go to 500
2539 if ((modedc.eq.2).and.(nosolv.ne.0)) go to 500
2540 if (initf.eq.4) go to 500
2541 go to 700
2542c
2543c charge storage elements
2544c
2545 500 tf=value(locm+24)
2546 tr=value(locm+33)
2547 czbe=value(locm+21)*area
2548 pe=value(locm+22)
2549 xme=value(locm+23)
2550 cdis=value(locm+32)
2551 ctot=value(locm+29)*area
2552 czbc=ctot*cdis
2553 czbx=ctot-czbc
2554 pc=value(locm+30)
2555 xmc=value(locm+31)
2556 fcpe=value(locm+46)
2557 czcs=value(locm+38)*area
2558 ps=value(locm+39)
2559 xms=value(locm+40)
2560 xtf=value(locm+25)
2561 ovtf=value(locm+26)
2562 xjtf=value(locm+27)*area
2563 if(tf.eq.0.0d0) go to 505
2564 if(vbe.le.0.0d0) go to 505
2565 argtf=0.0d0
2566 arg2=0.0d0
2567 arg3=0.0d0
2568 if(xtf.eq.0.0d0) go to 504
2569 argtf=xtf
2570 if(ovtf.ne.0.0d0) argtf=argtf*dexp(vbc*ovtf)
2571 arg2=argtf
2572 if(xjtf.eq.0.0d0) go to 503
2573 temp=cbe/(cbe+xjtf)
2574 argtf=argtf*temp*temp
2575 arg2=argtf*(3.0d0-temp-temp)
2576 503 arg3=cbe*argtf*ovtf
2577 504 cbe=cbe*(1.0d0+argtf)/qb
2578 gbe=(gbe*(1.0d0+arg2)-cbe*dqbdve)/qb
2579 geqcb=tf*(arg3-cbe*dqbdvc)/qb
2580 505 if (vbe.ge.fcpe) go to 510
2581 arg=1.0d0-vbe/pe
2582 sarg=dexp(-xme*dlog(arg))
2583 qbe(lx0+loct)=tf*cbe+pe*czbe*(1.0d0-arg*sarg)/(1.0d0-xme)
2584 capbe=tf*gbe+czbe*sarg
2585 go to 520
2586 510 f1=value(locm+47)
2587 f2=value(locm+48)
2588 f3=value(locm+49)
2589 czbef2=czbe/f2
2590 qbe(lx0+loct)=tf*cbe+czbe*f1+czbef2*(f3*(vbe-fcpe)
2591 1 +(xme/(pe+pe))*(vbe*vbe-fcpe*fcpe))
2592 capbe=tf*gbe+czbef2*(f3+xme*vbe/pe)
2593 520 fcpc=value(locm+50)
2594 f1=value(locm+51)
2595 f2=value(locm+52)
2596 f3=value(locm+53)
2597 if (vbc.ge.fcpc) go to 530
2598 arg=1.0d0-vbc/pc
2599 sarg=dexp(-xmc*dlog(arg))
2600 qbc(lx0+loct)=tr*cbc+pc*czbc*(1.0d0-arg*sarg)/(1.0d0-xmc)
2601 capbc=tr*gbc+czbc*sarg
2602 go to 540
2603 530 czbcf2=czbc/f2
2604 qbc(lx0+loct)=tr*cbc+czbc*f1+czbcf2*(f3*(vbc-fcpc)
2605 1 +(xmc/(pc+pc))*(vbc*vbc-fcpc*fcpc))
2606 capbc=tr*gbc+czbcf2*(f3+xmc*vbc/pc)
2607 540 if(vbx.ge.fcpc) go to 550
2608 arg=1.0d0-vbx/pc
2609 sarg=dexp(-xmc*dlog(arg))
2610 qbx(lx0+loct)=pc*czbx*(1.0d0-arg*sarg)/(1.0d0-xmc)
2611 capbx=czbx*sarg
2612 go to 560
2613 550 czbxf2=czbx/f2
2614 qbx(lx0+loct)=czbx*f1+czbxf2*(f3*(vbx-fcpc)+(xmc/(pc+pc))*
2615 1 (vbx*vbx-fcpc*fcpc))
2616 capbx=czbxf2*(f3+xmc*vbx/pc)
2617 560 if(vcs.ge.0.0d0) go to 570
2618 arg=1.0d0-vcs/ps
2619 sarg=dexp(-xms*dlog(arg))
2620 qcs(lx0+loct)=ps*czcs*(1.0d0-arg*sarg)/(1.0d0-xms)
2621 capcs=czcs*sarg
2622 go to 580
2623 570 qcs(lx0+loct)=vcs*czcs*(1.0d0+xms*vcs/(2.0d0*ps))
2624 capcs=czcs*(1.0d0+xms*vcs/ps)
2625c
2626c store small-signal parameters
2627c
2628 580 if ((mode.eq.1).and.(modedc.eq.2).and.(nosolv.ne.0)) go to 700
2629 if (initf.ne.4) go to 600
2630 value(lx0+loct+9)=capbe
2631 value(lx0+loct+11)=capbc
2632 value(lx0+loct+13)=capcs
2633 value(lx0+loct+15)=capbx
2634 value(lx0+loct+17)=geqcb
2635 go to 1000
2636c
2637c transient analysis
2638c
2639 600 if (initf.ne.5) go to 610
2640 qbe(lx1+loct)=qbe(lx0+loct)
2641 qbc(lx1+loct)=qbc(lx0+loct)
2642 qbx(lx1+loct)=qbx(lx0+loct)
2643 qcs(lx1+loct)=qcs(lx0+loct)
2644 610 call intgr8(geq,ceq,capbe,loct+8)
2645 geqcb=geqcb*ag(1)
2646 gpi=gpi+geq
2647 cb=cb+cqbe(lx0+loct)
2648 call intgr8(geq,ceq,capbc,loct+10)
2649 gmu=gmu+geq
2650 cb=cb+cqbc(lx0+loct)
2651 cc=cc-cqbc(lx0+loct)
2652 call intgr8(gccs,ceq,capcs,loct+12)
2653 ceqcs=type*(cqcs(lx0+loct)-vcs*gccs)
2654 call intgr8(geqbx,ceq,capbx,loct+14)
2655 ceqbx=type*(cqbx(lx0+loct)-vbx*geqbx)
2656 if (initf.ne.5) go to 700
2657 cqbe(lx1+loct)=cqbe(lx0+loct)
2658 cqbc(lx1+loct)=cqbc(lx0+loct)
2659 cqbx(lx1+loct)=cqbx(lx0+loct)
2660 cqcs(lx1+loct)=cqcs(lx0+loct)
2661c
2662c check convergence
2663c
2664 700 if (initf.ne.3) go to 710
2665 if (ioff.eq.0) go to 710
2666 go to 750
2667 710 if (icheck.eq.1) go to 720
2668 tol=reltol*dmax1(dabs(cchat),dabs(cc))+abstol
2669 if (dabs(cchat-cc).gt.tol) go to 720
2670 tol=reltol*dmax1(dabs(cbhat),dabs(cb))+abstol
2671 if (dabs(cbhat-cb).le.tol) go to 750
2672 720 noncon=noncon+1
2673 750 vbeo(lx0+loct)=vbe
2674 vbco(lx0+loct)=vbc
2675 cco(lx0+loct)=cc
2676 cbo(lx0+loct)=cb
2677 gpio(lx0+loct)=gpi
2678 gmuo(lx0+loct)=gmu
2679 gmo(lx0+loct)=gm
2680 goo(lx0+loct)=go
2681 gxo(lx0+loct)=gx
2682c
2683c load current excitation vector
2684c
2685 900 ceqbe=type*(cc+cb-vbe*(gm+go+gpi)+vbc*(go-geqcb))
2686 ceqbc=type*(-cc+vbe*(gm+go)-vbc*(gmu+go))
2687 value(lvn+node2)=value(lvn+node2)-ceqbx
2688 value(lvn+node4)=value(lvn+node4)+ceqcs+ceqbx+ceqbc
2689 value(lvn+node5)=value(lvn+node5)-ceqbe-ceqbc
2690 value(lvn+node6)=value(lvn+node6)+ceqbe
2691 value(lvn+node7)=value(lvn+node7)-ceqcs
2692c
2693c load y matrix
2694c
2695 locy=lynl+nodplc(loc+24)
2696 value(locy)=value(locy)+gcpr
2697 locy=lynl+nodplc(loc+25)
2698 value(locy)=value(locy)+gx+geqbx
2699 locy=lynl+nodplc(loc+26)
2700 value(locy)=value(locy)+gepr
2701 locy=lynl+nodplc(loc+27)
2702 value(locy)=value(locy)+gmu+go+gcpr+gccs+geqbx
2703 locy=lynl+nodplc(loc+28)
2704 value(locy)=value(locy)+gx +gpi+gmu+geqcb
2705 locy=lynl+nodplc(loc+29)
2706 value(locy)=value(locy)+gpi+gepr+gm+go
2707 locy=lynl+nodplc(loc+10)
2708 value(locy)=value(locy)-gcpr
2709 locy=lynl+nodplc(loc+11)
2710 value(locy)=value(locy)-gx
2711 locy=lynl+nodplc(loc+12)
2712 value(locy)=value(locy)-gepr
2713 locy=lynl+nodplc(loc+13)
2714 value(locy)=value(locy)-gcpr
2715 locy=lynl+nodplc(loc+14)
2716 value(locy)=value(locy)-gmu+gm
2717 locy=lynl+nodplc(loc+15)
2718 value(locy)=value(locy)-gm-go
2719 locy=lynl+nodplc(loc+16)
2720 value(locy)=value(locy)-gx
2721 locy=lynl+nodplc(loc+17)
2722 value(locy)=value(locy)-gmu-geqcb
2723 locy=lynl+nodplc(loc+18)
2724 value(locy)=value(locy)-gpi
2725 locy=lynl+nodplc(loc+19)
2726 value(locy)=value(locy)-gepr
2727 locy=lynl+nodplc(loc+20)
2728 value(locy)=value(locy)-go+geqcb
2729 locy=lynl+nodplc(loc+21)
2730 value(locy)=value(locy)-gpi-gm-geqcb
2731 locy=lynl+nodplc(loc+31)
2732 value(locy)=value(locy)+gccs
2733 locy=lynl+nodplc(loc+32)
2734 value(locy)=value(locy)-gccs
2735 locy=lynl+nodplc(loc+33)
2736 value(locy)=value(locy)-gccs
2737 locy=lynl+nodplc(loc+34)
2738 value(locy)=value(locy)-geqbx
2739 locy=lynl+nodplc(loc+35)
2740 value(locy)=value(locy)-geqbx
2741 1000 loc=nodplc(loc)
2742 go to 10
2743 end
2744 subroutine fetlim(vnew,vold,vto,icheck)
2745c
2746c *** fetlim is not used in this version ***
2747c * if problems arrise with the conver- *
2748c * gence of MOSFET circuit it should be *
2749c * re-installed. R.Newton. *
2750c *** ***
2751c
2752 return
2753 end
2754 subroutine jfet
2755 implicit double precision (a-h,o-z)
2756c
2757c this routine processes jfets for dc and transient analyses.
2758c
2759 common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem,
2760 1 isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize,
2761 2 junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr,
2762 3 nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1,
2763 4 lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd,
2764 5 imynl,imvn,lcvn,loutpt,nsnod,nsmat,nsval,icnod,icmat,icval
2765 common /cirdat/ locate(50),jelcnt(50),nunods,ncnods,numnod,nstop,
2766 1 nut,nlt,nxtrm,ndist,ntlin,ibr,numvs
2767 common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet,
2768 1 xmu,mode,modedc,icalc,initf,method,iord,maxord,noncon,iterno,
2769 2 itemno,nosolv,ipostp,iscrch
2770 common /knstnt/ twopi,xlog2,xlog10,root2,rad,boltz,charge,ctok,
2771 1 gmin,reltol,abstol,vntol,trtol,chgtol,eps0,epssil,epsox
2772 common /blank/ value(1000)
2773 integer nodplc(64)
2774 complex*16 cvalue(32)
2775 equivalence (value(1),nodplc(1),cvalue(1))
2776c
2777c
2778 dimension vgso(1),vgdo(1),cgo(1),cdo(1),cgdo(1),gmo(1),gdso(1),
2779 1 ggso(1),ggdo(1),qgs(1),cqgs(1),qgd(1),cqgd(1)
2780 equivalence (vgso(1),value( 1)),(vgdo(1),value( 2)),
2781 1 (cgo (1),value( 3)),(cdo (1),value( 4)),
2782 2 (cgdo(1),value( 5)),(gmo (1),value( 6)),
2783 3 (gdso(1),value( 7)),(ggso(1),value( 8)),
2784 4 (ggdo(1),value( 9)),(qgs (1),value(10)),
2785 5 (cqgs(1),value(11)),(qgd (1),value(12)),
2786 6 (cqgd(1),value(13))
2787c
2788c
2789 loc=locate(13)
2790 10 if (loc.eq.0) return
2791 locv=nodplc(loc+1)
2792 node1=nodplc(loc+2)
2793 node2=nodplc(loc+3)
2794 node3=nodplc(loc+4)
2795 node4=nodplc(loc+5)
2796 node5=nodplc(loc+6)
2797 locm=nodplc(loc+7)
2798 ioff=nodplc(loc+8)
2799 type=nodplc(locm+2)
2800 locm=nodplc(locm+1)
2801 loct=nodplc(loc+19)
2802c
2803c dc model parameters
2804c
2805 area=value(locv+1)
2806 vto=value(locm+1)
2807 beta=value(locm+2)*area
2808 xlamb=value(locm+3)
2809 gdpr=value(locm+4)*area
2810 gspr=value(locm+5)*area
2811 csat=value(locm+9)*area
2812 vcrit=value(locm+16)
2813c
2814c initialization
2815c
2816 icheck=1
2817 go to (100,20,30,50,60,70), initf
2818 20 if(mode.ne.1.or.modedc.ne.2.or.nosolv.eq.0) go to 25
2819 vds=type*value(locv+2)
2820 vgs=type*value(locv+3)
2821 vgd=vgs-vds
2822 go to 300
2823 25 if(ioff.ne.0) go to 40
2824 vgs=-1.0d0
2825 vgd=-1.0d0
2826 go to 300
2827 30 if (ioff.eq.0) go to 100
2828 40 vgs=0.0d0
2829 vgd=0.0d0
2830 go to 300
2831 50 vgs=vgso(lx0+loct)
2832 vgd=vgdo(lx0+loct)
2833 go to 300
2834 60 vgs=vgso(lx1+loct)
2835 vgd=vgdo(lx1+loct)
2836 go to 300
2837 70 xfact=delta/delold(2)
2838 vgso(lx0+loct)=vgso(lx1+loct)
2839 vgs=(1.0d0+xfact)*vgso(lx1+loct)-xfact*vgso(lx2+loct)
2840 vgdo(lx0+loct)=vgdo(lx1+loct)
2841 vgd=(1.0d0+xfact)*vgdo(lx1+loct)-xfact*vgdo(lx2+loct)
2842 cgo(lx0+loct)=cgo(lx1+loct)
2843 cdo(lx0+loct)=cdo(lx1+loct)
2844 cgdo(lx0+loct)=cgdo(lx1+loct)
2845 gmo(lx0+loct)=gmo(lx1+loct)
2846 gdso(lx0+loct)=gdso(lx1+loct)
2847 ggso(lx0+loct)=ggso(lx1+loct)
2848 ggdo(lx0+loct)=ggdo(lx1+loct)
2849 go to 110
2850c
2851c compute new nonlinear branch voltages
2852c
2853 100 vgs=type*(value(lvnim1+node2)-value(lvnim1+node5))
2854 vgd=type*(value(lvnim1+node2)-value(lvnim1+node4))
2855 110 delvgs=vgs-vgso(lx0+loct)
2856 delvgd=vgd-vgdo(lx0+loct)
2857 delvds=delvgs-delvgd
2858 cghat=cgo(lx0+loct)+ggdo(lx0+loct)*delvgd+ggso(lx0+loct)*delvgs
2859 cdhat=cdo(lx0+loct)+gmo(lx0+loct)*delvgs+gdso(lx0+loct)*delvds
2860 1 -ggdo(lx0+loct)*delvgd
2861c
2862c bypass if solution has not changed
2863c
2864 if (initf.eq.6) go to 200
2865 tol=reltol*dmax1(dabs(vgs),dabs(vgso(lx0+loct)))+vntol
2866 if (dabs(delvgs).ge.tol) go to 200
2867 tol=reltol*dmax1(dabs(vgd),dabs(vgdo(lx0+loct)))+vntol
2868 if (dabs(delvgd).ge.tol) go to 200
2869 tol=reltol*dmax1(dabs(cghat),dabs(cgo(lx0+loct)))+abstol
2870 if (dabs(cghat-cgo(lx0+loct)).ge.tol) go to 200
2871 tol=reltol*dmax1(dabs(cdhat),dabs(cdo(lx0+loct)))+abstol
2872 if (dabs(cdhat-cdo(lx0+loct)).ge.tol) go to 200
2873 vgs=vgso(lx0+loct)
2874 vgd=vgdo(lx0+loct)
2875 vds=vgs-vgd
2876 cg=cgo(lx0+loct)
2877 cd=cdo(lx0+loct)
2878 cgd=cgdo(lx0+loct)
2879 gm=gmo(lx0+loct)
2880 gds=gdso(lx0+loct)
2881 ggs=ggso(lx0+loct)
2882 ggd=ggdo(lx0+loct)
2883 go to 900
2884c
2885c limit nonlinear branch voltages
2886c
2887 200 icheck=0
2888 call pnjlim(vgs,vgso(lx0+loct),vt,vcrit,icheck)
2889 call pnjlim(vgd,vgdo(lx0+loct),vt,vcrit,icheck)
2890 call fetlim(vgs,vgso(lx0+loct),vto,icheck)
2891 call fetlim(vgd,vgdo(lx0+loct),vto,icheck)
2892c
2893c determine dc current and derivatives
2894c
2895 300 vds=vgs-vgd
2896 if (vgs.gt.-5.0d0*vt) go to 310
2897 ggs=-csat/vgs+gmin
2898 cg=ggs*vgs
2899 go to 320
2900 310 evgs=dexp(vgs/vt)
2901 ggs=csat*evgs/vt+gmin
2902 cg=csat*(evgs-1.0d0)+gmin*vgs
2903 320 if (vgd.gt.-5.0d0*vt) go to 330
2904 ggd=-csat/vgd+gmin
2905 cgd=ggd*vgd
2906 go to 340
2907 330 evgd=dexp(vgd/vt)
2908 ggd=csat*evgd/vt+gmin
2909 cgd=csat*(evgd-1.0d0)+gmin*vgd
2910 340 cg=cg+cgd
2911c
2912c compute drain current and derivitives for normal mode
2913c
2914 400 if (vds.lt.0.0d0) go to 450
2915 vgst=vgs-vto
2916c
2917c normal mode, cutoff region
2918c
2919 if (vgst.gt.0.0d0) go to 410
2920 cdrain=0.0d0
2921 gm=0.0d0
2922 gds=0.0d0
2923 go to 490
2924c
2925c normal mode, saturation region
2926c
2927 410 betap=beta*(1.0d0+xlamb*vds)
2928 twob=betap+betap
2929 if (vgst.gt.vds) go to 420
2930 cdrain=betap*vgst*vgst
2931 gm=twob*vgst
2932 gds=xlamb*beta*vgst*vgst
2933 go to 490
2934c
2935c normal mode, linear region
2936c
2937 420 cdrain=betap*vds*(vgst+vgst-vds)
2938 gm=twob*vds
2939 gds=twob*(vgst-vds)+xlamb*beta*vds*(vgst+vgst-vds)
2940 go to 490
2941c
2942c compute drain current and derivitives for inverse mode
2943c
2944 450 vgdt=vgd-vto
2945c
2946c inverse mode, cutoff region
2947c
2948 if (vgdt.gt.0.0d0) go to 460
2949 cdrain=0.0d0
2950 gm=0.0d0
2951 gds=0.0d0
2952 go to 490
2953c
2954c inverse mode, saturation region
2955c
2956 460 betap=beta*(1.0d0-xlamb*vds)
2957 twob=betap+betap
2958 if (vgdt.gt.-vds) go to 470
2959 cdrain=-betap*vgdt*vgdt
2960 gm=-twob*vgdt
2961 gds=xlamb*beta*vgdt*vgdt-gm
2962 go to 490
2963c
2964c inverse mode, linear region
2965c
2966 470 cdrain=betap*vds*(vgdt+vgdt+vds)
2967 gm=twob*vds
2968 gds=twob*vgdt-xlamb*beta*vds*(vgdt+vgdt+vds)
2969c
2970c compute equivalent drain current source
2971c
2972 490 cd=cdrain-cgd
2973 if (mode.ne.1) go to 500
2974 if ((modedc.eq.2).and.(nosolv.ne.0)) go to 500
2975 if (initf.eq.4) go to 500
2976 go to 700
2977c
2978c charge storage elements
2979c
2980 500 czgs=value(locm+6)*area
2981 czgd=value(locm+7)*area
2982 phib=value(locm+8)
2983 twop=phib+phib
2984 fcpb=value(locm+12)
2985 fcpb2=fcpb*fcpb
2986 f1=value(locm+13)
2987 f2=value(locm+14)
2988 f3=value(locm+15)
2989 czgsf2=czgs/f2
2990 czgdf2=czgd/f2
2991 if (vgs.ge.fcpb) go to 510
2992 sarg=dsqrt(1.0d0-vgs/phib)
2993 qgs(lx0+loct)=twop*czgs*(1.0d0-sarg)
2994 capgs=czgs/sarg
2995 go to 520
2996 510 qgs(lx0+loct)=czgs*f1+czgsf2*(f3*(vgs-fcpb)
2997 1 +(vgs*vgs-fcpb2)/(twop+twop))
2998 capgs=czgsf2*(f3+vgs/twop)
2999 520 if (vgd.ge.fcpb) go to 530
3000 sarg=dsqrt(1.0d0-vgd/phib)
3001 qgd(lx0+loct)=twop*czgd*(1.0d0-sarg)
3002 capgd=czgd/sarg
3003 go to 560
3004 530 qgd(lx0+loct)=czgd*f1+czgdf2*(f3*(vgd-fcpb)
3005 1 +(vgd*vgd-fcpb2)/(twop+twop))
3006 capgd=czgdf2*(f3+vgd/twop)
3007c
3008c store small-signal parameters
3009c
3010 560 if ((mode.eq.1).and.(modedc.eq.2).and.(nosolv.ne.0)) go to 700
3011 if (initf.ne.4) go to 600
3012 value(lx0+loct+9)=capgs
3013 value(lx0+loct+11)=capgd
3014 go to 1000
3015c
3016c transient analysis
3017c
3018 600 if (initf.ne.5) go to 610
3019 qgs(lx1+loct)=qgs(lx0+loct)
3020 qgd(lx1+loct)=qgd(lx0+loct)
3021 610 call intgr8(geq,ceq,capgs,loct+9)
3022 ggs=ggs+geq
3023 cg=cg+cqgs(lx0+loct)
3024 call intgr8(geq,ceq,capgd,loct+11)
3025 ggd=ggd+geq
3026 cg=cg+cqgd(lx0+loct)
3027 cd=cd-cqgd(lx0+loct)
3028 cgd=cgd+cqgd(lx0+loct)
3029 if (initf.ne.5) go to 700
3030 cqgs(lx1+loct)=cqgs(lx0+loct)
3031 cqgd(lx1+loct)=cqgd(lx0+loct)
3032c
3033c check convergence
3034c
3035 700 if (initf.ne.3) go to 710
3036 if (ioff.eq.0) go to 710
3037 go to 750
3038 710 if (icheck.eq.1) go to 720
3039 tol=reltol*dmax1(dabs(cghat),dabs(cg))+abstol
3040 if (dabs(cghat-cg).ge.tol) go to 720
3041 tol=reltol*dmax1(dabs(cdhat),dabs(cd))+abstol
3042 if (dabs(cdhat-cd).le.tol) go to 750
3043 720 noncon=noncon+1
3044 750 vgso(lx0+loct)=vgs
3045 vgdo(lx0+loct)=vgd
3046 cgo(lx0+loct)=cg
3047 cdo(lx0+loct)=cd
3048 cgdo(lx0+loct)=cgd
3049 gmo(lx0+loct)=gm
3050 gdso(lx0+loct)=gds
3051 ggso(lx0+loct)=ggs
3052 ggdo(lx0+loct)=ggd
3053c
3054c load current vector
3055c
3056 900 ceqgd=type*(cgd-ggd*vgd)
3057 ceqgs=type*((cg-cgd)-ggs*vgs)
3058 cdreq=type*((cd+cgd)-gds*vds-gm*vgs)
3059 value(lvn+node2)=value(lvn+node2)-ceqgs-ceqgd
3060 value(lvn+node4)=value(lvn+node4)-cdreq+ceqgd
3061 value(lvn+node5)=value(lvn+node5)+cdreq+ceqgs
3062c
3063c load y matrix
3064c
3065 locy=lynl+nodplc(loc+20)
3066 value(locy)=value(locy)+gdpr
3067 locy=lynl+nodplc(loc+21)
3068 value(locy)=value(locy)+ggd+ggs
3069 locy=lynl+nodplc(loc+22)
3070 value(locy)=value(locy)+gspr
3071 locy=lynl+nodplc(loc+23)
3072 value(locy)=value(locy)+gdpr+gds+ggd
3073 locy=lynl+nodplc(loc+24)
3074 value(locy)=value(locy)+gspr+gds+gm+ggs
3075 locy=lynl+nodplc(loc+9)
3076 value(locy)=value(locy)-gdpr
3077 locy=lynl+nodplc(loc+10)
3078 value(locy)=value(locy)-ggd
3079 locy=lynl+nodplc(loc+11)
3080 value(locy)=value(locy)-ggs
3081 locy=lynl+nodplc(loc+12)
3082 value(locy)=value(locy)-gspr
3083 locy=lynl+nodplc(loc+13)
3084 value(locy)=value(locy)-gdpr
3085 locy=lynl+nodplc(loc+14)
3086 value(locy)=value(locy)+gm-ggd
3087 locy=lynl+nodplc(loc+15)
3088 value(locy)=value(locy)-gds-gm
3089 locy=lynl+nodplc(loc+16)
3090 value(locy)=value(locy)-ggs-gm
3091 locy=lynl+nodplc(loc+17)
3092 value(locy)=value(locy)-gspr
3093 locy=lynl+nodplc(loc+18)
3094 value(locy)=value(locy)-gds
3095 1000 loc=nodplc(loc)
3096 go to 10
3097 end
3098 subroutine mosfet
3099 implicit double precision (a-h,o-z)
3100c
3101c this routine processes mosfets for dc and transient analyses.
3102c
3103 common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem,
3104 1 isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize,
3105 2 junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr,
3106 3 nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1,
3107 4 lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd,
3108 5 imynl,imvn,lcvn,loutpt,nsnod,nsmat,nsval,icnod,icmat,icval
3109 common /cirdat/ locate(50),jelcnt(50),nunods,ncnods,numnod,nstop,
3110 1 nut,nlt,nxtrm,ndist,ntlin,ibr,numvs
3111 common /mosarg/ gamma,beta,vto,phi,cox,vbi,xnfs,xnsub,xd,xj,xl,
3112 1 xlamda,utra,uexp,vbp,von,vdsat,theta,vcrit,vtra,gleff,cdrain
3113 common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet,
3114 1 xmu,mode,modedc,icalc,initf,method,iord,maxord,noncon,iterno,
3115 2 itemno,nosolv,ipostp,iscrch
3116 common /knstnt/ twopi,xlog2,xlog10,root2,rad,boltz,charge,ctok,
3117 1 gmin,reltol,abstol,vntol,trtol,chgtol,eps0,epssil,epsox
3118 common /blank/ value(1000)
3119 integer nodplc(64)
3120 complex*16 cvalue(32)
3121 equivalence (value(1),nodplc(1),cvalue(1))
3122c
3123c
3124 dimension vbdo(1),vbso(1),vgbo(1),gdgo(1),cdo(1),cbo(1),
3125 1 gddo(1),gdso(1),gbgo(1),gbdo(1),gbso(1),qb(1),cqb(1),qg(1),
3126 2 cqg(1),qd(1),cqd(1)
3127c.. note: for direct indexing with 'value', use, e.g. loct+2 to find vgbo
3128 equivalence (vbdo (1),value( 1)),(vbso (1),value( 2)),
3129 1 (vgbo (1),value( 3)),
3130 2 (cdo (1),value( 5)),(cbo (1),value( 6)),
3131 3 (gddo (1),value( 7)),(gdgo (1),value(8)),
3132 4 (gdso (1),value( 9)),(gbgo (1),value(10)),
3133 5 (gbdo (1),value(11)),(gbso (1),value(12)),
3134 6 (qb (1),value(13)),(cqb (1),value(14)),
3135 7 (qg (1),value(15)),(cqg (1),value(16)),
3136 8 (qd (1),value(17)),(cqd (1),value(18))
3137c
3138c
3139 loc=locate(14)
3140 10 if (loc.eq.0) return
3141 locm=nodplc(loc+8)
3142 if(nodplc(locm+2).ne.0) go to 15
3143 call gasfet(loc)
3144 go to 1000
3145 15 locv=nodplc(loc+1)
3146 node1=nodplc(loc+2)
3147 node2=nodplc(loc+3)
3148 node3=nodplc(loc+4)
3149 node4=nodplc(loc+5)
3150 node5=nodplc(loc+6)
3151 node6=nodplc(loc+7)
3152 ioff=nodplc(loc+9)
3153 type=nodplc(locm+2)
3154 locm=nodplc(locm+1)
3155 loct=nodplc(loc+26)
3156c
3157c dc model parameters
3158c
3159 xj=value(locm+19)
3160 xl=value(locv+1)-2.0d0*value(locm+20)
3161 xw=value(locv+2)-2.0d0*value(locm+36)
3162 devmod=value(locv+8)
3163 vto=type*value(locm+1)
3164 beta=value(locm+2)*xw/xl
3165 gamma=value(locm+3)
3166 phi=value(locm+4)
3167 xlamda=value(locm+5)
3168 csat=value(locm+15)
3169 ad=value(locv+3)
3170 as=value(locv+4)
3171 cdsat=csat*ad
3172 cssat=csat*as
3173 gdpr=value(locv+11)
3174 gspr=value(locv+12)
3175 covlgs=value(locm+8)*xw
3176 covlgd=value(locm+9)*xw
3177 covlgb=value(locm+10)*xl
3178 cox=value(locm+13)
3179 xnsub=value(locm+16)
3180 xnfs=value(locm+18)
3181 vbp=value(locm+24)
3182 uexp=value(locm+25)
3183 utra=value(locm+26)
3184 vbi=type*value(locm+34)
3185 xd=value(locm+35)
3186 vcrit=value(locm+37)*xl
3187 vtra=value(locm+38)*xl
3188 theta=value(locm+39)
3189 gleff=value(locm+40)
3190c
3191c initialization
3192c
3193 icheck=1
3194 go to (100,20,30,50,60,70), initf
3195 20 if(mode.ne.1.or.modedc.ne.2.or.nosolv.eq.0) go to 25
3196 vbs=type*value(locv+7)
3197 vbd=type*value(locv+5)-vbs
3198 vgb=type*value(locv+6)-vbs
3199 go to 300
3200 25 if(ioff.ne.0) go to 40
3201 vbs=0.0d0
3202 vbd=0.0d0
3203 vgb=vto
3204 go to 300
3205 30 if (ioff.eq.0) go to 100
3206 40 vbs=0.0d0
3207 vbd=0.0d0
3208 vgb=0.0d0
3209 go to 300
3210 50 vbs=vbso(lx0+loct)
3211 vbd=vbdo(lx0+loct)
3212 vgb=vgbo(lx0+loct)
3213 go to 300
3214 60 vbs=vbso(lx1+loct)
3215 vbd=vbdo(lx1+loct)
3216 vgb=vgbo(lx1+loct)
3217 go to 300
3218 70 xfact=delta/delold(2)
3219 vbso(lx0+loct)=vbso(lx1+loct)
3220 vbs=(1.0d0+xfact)*vbso(lx1+loct)-xfact*vbso(lx2+loct)
3221 vbdo(lx0+loct)=vbdo(lx1+loct)
3222 vbd=(1.0d0+xfact)*vbdo(lx1+loct)-xfact*vbdo(lx2+loct)
3223 vgbo(lx0+loct)=vgbo(lx1+loct)
3224 vgb=(1.0d0+xfact)*vgbo(lx1+loct)-xfact*vgbo(lx2+loct)
3225 cdo(lx0+loct)=cdo(lx1+loct)
3226 cbo(lx0+loct)=cbo(lx1+loct)
3227 gdgo(lx0+loct)=gdgo(lx1+loct)
3228 gddo(lx0+loct)=gddo(lx1+loct)
3229 gdso(lx0+loct)=gdso(lx1+loct)
3230 gbgo(lx0+loct)=gbgo(lx1+loct)
3231 gbdo(lx0+loct)=gbdo(lx1+loct)
3232 gbso(lx0+loct)=gbso(lx1+loct)
3233 go to 110
3234c
3235c compute new nonlinear branch voltages
3236c
3237 100 vbs=type*(value(lvnim1+node4)-value(lvnim1+node6))
3238 vbd=type*(value(lvnim1+node4)-value(lvnim1+node5))
3239 vgb=type*(value(lvnim1+node2)-value(lvnim1+node4))
3240 110 delvbs=vbs-vbso(lx0+loct)
3241 delvbd=vbd-vbdo(lx0+loct)
3242 delvgb=vgb-vgbo(lx0+loct)
3243 cdhat=cdo(lx0+loct)+gdgo(lx0+loct)*delvgb-gddo(lx0+loct)*delvbd
3244 1 -gdso(lx0+loct)*delvbs
3245 cbhat=cbo(lx0+loct)+gbgo(lx0+loct)*delvgb-gbdo(lx0+loct)*delvbd
3246 1 -gbso(lx0+loct)*delvbs
3247c
3248c bypass if solution has not changed
3249c
3250c********** kill bypass for now!!!!!
3251 if (6 .eq.6) go to 200
3252 tol=reltol*dmax1(dabs(vbs),dabs(vbso(lx0+loct)))+vntol
3253 if (dabs(delvbs).ge.tol) go to 200
3254 tol=reltol*dmax1(dabs(vbd),dabs(vbdo(lx0+loct)))+vntol
3255 if (dabs(delvbd).ge.tol) go to 200
3256 tol=reltol*dmax1(dabs(vgb),dabs(vgbo(lx0+loct)))+vntol
3257 if (dabs(delvgb).ge.tol) go to 200
3258 tol=reltol*dmax1(dabs(cdhat),dabs(cdo(lx0+loct)))+abstol
3259 if (dabs(cdhat-cdo(lx0+loct)).ge.tol) go to 200
3260 tol=reltol*dmax1(dabs(cbhat),dabs(cbo(lx0+loct)))+abstol
3261 if (dabs(cbhat-(cbo(lx0+loct))).ge.tol) go to 200
3262 vbd=vbdo(lx0+loct)
3263 vbs=vbso(lx0+loct)
3264 vgb=vgbo(lx0+loct)
3265 cdrain=cdo(lx0+loct)
3266 cbulk=cbo(lx0+loct)
3267 gccdg=gdgo(lx0+loct)
3268 gccdd=gddo(lx0+loct)
3269 gccds=gdso(lx0+loct)
3270 gccbg=gbgo(lx0+loct)
3271 gccbd=gbdo(lx0+loct)
3272 gccbs=gbso(lx0+loct)
3273 go to 900
3274c
3275c limit nonlinear branch voltages
3276c
3277 200 von=type*value(locv+9)
3278 icheck=0
3279 call fetlim(vgb,vgbo(lx0+loct),von,icheck)
3280 vcornr=0.0d0
3281c if(vbs.gt.0.0d0) vcornr=vt*dlog(vt/(root2*cssat))
3282 call pnjlim(vbs,vbso(lx0+loct),vt,vcornr,icheck)
3283c vbs=dmax1(vbso(lx0+loct)-10.0d0,vbs)
3284 vcornr=0.0d0
3285c if(vbd.gt.0.0d0) vcornr=vt*dlog(vt/(root2*cdsat))
3286 call pnjlim(vbd,vbdo(lx0+loct),vt,vcornr,icheck)
3287c vbd=dmax1(vbdo(lx0+loct)-10.0d0,vbd)
3288c
3289c determine bulk-drain and bulk-source diode terms
3290c
3291 300 fivevt=-5.0d0*vt
3292 if (vbs.gt.fivevt) go to 310
3293 geqbs=-cssat/vbs+gmin
3294 cbulk=geqbs*vbs
3295 go to 320
3296 310 evbs=dexp(vbs/vt)
3297 geqbs=cssat*evbs/vt+gmin
3298 cbulk=cssat*(evbs-1.0d0)+gmin*vbs
3299 320 if (vbd.gt.fivevt) go to 330
3300 geqbd=-cdsat/vbd+gmin
3301 cbd=geqbd*vbd
3302 cbulk=cbulk+cbd
3303 go to 340
3304 330 evbd=dexp(vbd/vt)
3305 geqbd=cdsat*evbd/vt+gmin
3306 cbd=cdsat*(evbd-1.0d0)+gmin*vbd
3307 cbulk=cbulk+cbd
3308c.. cbd must also be subtracted from drain current
3309 340 continue
3310 gccdd=geqbd
3311 gccbd=-geqbd
3312 gccss=geqbs
3313 gccbs=-geqbs
3314 if (mode.ne.1) go to 350
3315c.. zero out some conductances and cgate
3316 cgate=0.0d0
3317 gccgg=0.0d0
3318 gccgd=0.0d0
3319 gccgs=0.0d0
3320 gccbg=0.0d0
3321c
3322c compute drain current and derivatives
3323c
3324 350 cox=cox*xl*xw
3325 if(vbd.gt.vbs) go to 360
3326c.. normal operation
3327 devmod=1.0d0
3328 call calcq(vgb,vbd,vbs,qgate,qchan,qbulk,
3329 1 ccgg,ccgd,ccgs,ccbg,ccbd,ccbs,didvg,didvd,didvs)
3330 didvg=beta*didvg
3331 didvd=beta*didvd
3332 didvs=beta*didvs
3333 go to 370
3334c.. inverted operation
3335 360 devmod=-1.0d0
3336 call calcq(vgb,vbs,vbd,qgate,qchan,qbulk,
3337 1 ccgg,ccgs,ccgd,ccbg,ccbs,ccbd,didvg,didvs,didvd)
3338 didvg=-beta*didvg
3339 didvd=-beta*didvd
3340 didvs=-beta*didvs
3341 cdrain=-cdrain
3342 370 cdrain=beta*cdrain-cbd
3343c if(mode.ne.1) write(6,6666) qgate,qchan,qbulk,time,delta
3344 6666 format(' qg qc qb',1p3e11.3,' time delta ',2e11.3)
3345 value(locv+8)=devmod
3346 value(locv+9)=type*von
3347 value(locv+10)=type*vdsat
3348 if(mode.ne.1) go to 500
3349 if ((modedc.eq.2).and.(nosolv.ne.0)) go to 500
3350 if (initf.eq.4) go to 500
3351 gccdg=didvg
3352 gccdd=gccdd+didvd
3353 gccds=didvs
3354 gccsg=-didvg
3355 gccsd=-didvd
3356 gccss=gccss-didvs
3357 go to 700
3358c
3359c charge storage elements
3360c
3361c.. bulk-drain and bulk-source depletion capacitances
3362 500 czbd=value(locm+11)*ad
3363 czbs=value(locm+12)*as
3364 phib=value(locm+14)
3365 twop=phib+phib
3366 fcpb=value(locm+29)
3367 if(vbs.lt.fcpb.and.vbd.lt.fcpb) go to 505
3368 fcpb2=fcpb*fcpb
3369 f1=value(locm+30)
3370 f2=value(locm+31)
3371 f3=value(locm+32)
3372 czbsf2=czbs/f2
3373 czbdf2=czbd/f2
3374 if (vbs.ge.fcpb) go to 510
3375 505 sarg=dsqrt(1.0d0-vbs/phib)
3376 qbs=twop*czbs*(1.0d0-sarg)
3377 capbs=czbs/sarg
3378 go to 520
3379 510 qbs=czbs*f1+czbsf2*(f3*(vbs-fcpb)
3380 1 +(vbs*vbs-fcpb2)/(twop+twop))
3381 capbs=czbsf2*(f3+vbs/twop)
3382 520 if (vbd.ge.fcpb) go to 530
3383 sarg=dsqrt(1.0d0-vbd/phib)
3384 qbd=twop*czbd*(1.0d0-sarg)
3385 capbd=czbd/sarg
3386 go to 540
3387 530 qbd=czbd*f1+czbdf2*(f3*(vbd-fcpb)
3388 1 +(vbd*vbd-fcpb2)/(twop+twop))
3389 capbd=czbdf2*(f3+vbd/twop)
3390c.. bulk and channel charge (plus overlaps)
3391 540 qgd=covlgd*(vgb+vbd)
3392 qgs=covlgs*(vgb+vbs)
3393 qgb=covlgb*vgb
3394 qg(lx0+loct)=qgate+qgb+qgd+qgs
3395 qd(lx0+loct)=qchan*0.5d0-qgd-qbd
3396 qb(lx0+loct)=qbulk+qbd+qbs-qgb
3397c
3398c store small-signal parameters
3399c
3400 590 if ((mode.eq.1).and.(modedc.eq.2).and.(nosolv.ne.0)) go to 900
3401 if (initf.ne.4) go to 600
3402 value(lx0+loct)=didvg
3403 value(lx0+loct+1)=didvd
3404 value(lx0+loct+2)=didvs
3405 value(lx0+loct+3)=geqbd
3406c.. cdrain is used in printing as well as noise calculation
3407 value(lx0+loct+4)=cdrain
3408 value(lx0+loct+5)=geqbs
3409c.. (loct+6 not used)
3410c.. this is the 'gm' term used in the noise calculation
3411 value(lx0+loct+7)=didvg
3412 value(lx0+loct+8)=ccgg
3413 value(lx0+loct+9)=ccgd
3414 value(lx0+loct+10)=ccgs
3415 value(lx0+loct+11)=ccbg
3416 value(lx0+loct+12)=ccbd
3417 value(lx0+loct+13)=ccbs
3418 value(lx0+loct+14)=capbd
3419 value(lx0+loct+15)=capbs
3420 go to 1000
3421c
3422c transient analysis
3423c
3424 600 if (initf.ne.5) go to 610
3425 qb(lx1+loct)=qb(lx0+loct)
3426 qg(lx1+loct)=qg(lx0+loct)
3427 qd(lx1+loct)=qd(lx0+loct)
3428c.. integrate qb
3429 610 call intgr8(geq,ceq,0.0d0,loct+12)
3430c.. integrate qg
3431 call intgr8(geq,ceq,0.0d0,loct+14)
3432c.. integrate qd
3433 call intgr8(geq,ceq,0.0d0,loct+16)
3434c.. divvey up the channel charge 50/50 to source and drain
3435c.. note that symmetry also precludes need for 'devmod' decisions
3436 gccg2=-0.5d0*(ccgg+ccbg)*ag(1)
3437 gccd2=-0.5d0*(ccgd+ccbd)*ag(1)
3438 gccs2=-0.5d0*(ccgs+ccbs)*ag(1)
3439 gccdg=gccg2+didvg-covlgd*ag(1)
3440 gccdd=gccdd+gccd2+didvd+(capbd+covlgd)*ag(1)
3441 gccds=gccs2+didvs
3442 gccsg=gccg2-didvg-covlgs*ag(1)
3443 gccsd=gccd2-didvd
3444 gccss=gccss+gccs2-didvs+(capbs+covlgs)*ag(1)
3445 gccgg=(ccgg+covlgd+covlgs+covlgb)*ag(1)
3446 gccgd=(ccgd-covlgd)*ag(1)
3447 gccgs=(ccgs-covlgs)*ag(1)
3448 gccbg=(ccbg-covlgb)*ag(1)
3449 gccbd=gccbd+(ccbd-capbd)*ag(1)
3450 gccbs=gccbs+(ccbs-capbs)*ag(1)
3451 cgate=cqg(lx0+loct)
3452 cbulk=cbulk+cqb(lx0+loct)
3453 cdrain=cdrain+cqd(lx0+loct)
3454 if (initf.ne.5) go to 700
3455 cqb(lx1+loct)=cqb(lx0+loct)
3456 cqg(lx1+loct)=cqg(lx0+loct)
3457 cqd(lx1+loct)=cqd(lx0+loct)
3458c
3459c check convergence
3460c
3461 700 if (initf.ne.3) go to 710
3462 if (ioff.ne.0) go to 750
3463 710 if (icheck.eq.1) go to 720
3464c tol=reltol*dmax1(dabs(cdhat),dabs(cdrain))+abstol
3465c if (dabs(cdhat-cdrain).ge.tol) go to 720
3466c tol=reltol*dmax1(dabs(cbhat),dabs(cbulk))+abstol
3467c if (dabs(cbhat-cbulk).le.tol) go to 750
3468 tol=reltol*dabs(vgb)+vntol
3469 if(dabs(delvgb).ge.tol) go to 720
3470 tol=reltol*dabs(vbd)+vntol
3471 if(dabs(delvbd).ge.tol) go to 720
3472 tol=reltol*dabs(vbs)+vntol
3473 if(dabs(delvbs).lt.tol) go to 750
3474 720 noncon=noncon+1
3475 750 vbdo(lx0+loct)=vbd
3476 vbso(lx0+loct)=vbs
3477 vgbo(lx0+loct)=vgb
3478 cdo(lx0+loct)=cdrain
3479 cbo(lx0+loct)=cbulk
3480 gdgo(lx0+loct)=gccdg
3481 gddo(lx0+loct)=gccdd
3482 gdso(lx0+loct)=gccds
3483 gbgo(lx0+loct)=gccbg
3484 gbdo(lx0+loct)=gccbd
3485 gbso(lx0+loct)=gccbs
3486c
3487c load current vector
3488c
3489 900 ceqg=type*(cgate-gccgg*vgb+gccgd*vbd+gccgs*vbs)
3490 ceqb=type*(cbulk-gccbg*vgb+gccbd*vbd+gccbs*vbs)
3491 ceqd=type*(cdrain-gccdg*vgb+gccdd*vbd+gccds*vbs)
3492 value(lvn+node2)=value(lvn+node2)-ceqg
3493 value(lvn+node4)=value(lvn+node4)-ceqb
3494 value(lvn+node5)=value(lvn+node5)-ceqd
3495 value(lvn+node6)=value(lvn+node6)+ceqd+ceqg+ceqb
3496c
3497c load y matrix
3498c
3499 locy=lynl+nodplc(loc+27)
3500 value(locy)=value(locy)+gdpr
3501 locy=lynl+nodplc(loc+28)
3502 value(locy)=value(locy)+gccgg
3503 locy=lynl+nodplc(loc+29)
3504 value(locy)=value(locy)+gspr
3505 locy=lynl+nodplc(loc+30)
3506 value(locy)=value(locy)-gccbg-gccbd-gccbs
3507 locy=lynl+nodplc(loc+31)
3508 value(locy)=value(locy)+gdpr+gccdd
3509 locy=lynl+nodplc(loc+32)
3510 value(locy)=value(locy)+gspr+gccss
3511 locy=lynl+nodplc(loc+10)
3512 value(locy)=value(locy)-gdpr
3513 locy=lynl+nodplc(loc+11)
3514 value(locy)=value(locy)-gccgg-gccgd-gccgs
3515 locy=lynl+nodplc(loc+12)
3516 value(locy)=value(locy)+gccgd
3517 locy=lynl+nodplc(loc+13)
3518 value(locy)=value(locy)+gccgs
3519 locy=lynl+nodplc(loc+14)
3520 value(locy)=value(locy)-gspr
3521 locy=lynl+nodplc(loc+15)
3522 value(locy)=value(locy)+gccbg
3523 locy=lynl+nodplc(loc+16)
3524 value(locy)=value(locy)+gccbd
3525 locy=lynl+nodplc(loc+17)
3526 value(locy)=value(locy)+gccbs
3527 locy=lynl+nodplc(loc+18)
3528 value(locy)=value(locy)-gdpr
3529 locy=lynl+nodplc(loc+19)
3530 value(locy)=value(locy)+gccdg
3531 locy=lynl+nodplc(loc+20)
3532 value(locy)=value(locy)-gccdg-gccdd-gccds
3533 locy=lynl+nodplc(loc+21)
3534 value(locy)=value(locy)+gccds
3535 locy=lynl+nodplc(loc+22)
3536 value(locy)=value(locy)+gccsg
3537 locy=lynl+nodplc(loc+23)
3538 value(locy)=value(locy)-gspr
3539 locy=lynl+nodplc(loc+24)
3540 value(locy)=value(locy)-gccsg-gccsd-gccss
3541 locy=lynl+nodplc(loc+25)
3542 value(locy)=value(locy)+gccsd
3543 1000 loc=nodplc(loc)
3544 go to 10
3545 end
3546 subroutine calcq(vgb,vbd,vbs,qg,qc,qb,
3547 1 ccgg,ccgd,ccgs,ccbg,ccbd,ccbs,
3548 2 didvg,didvd,didvs)
3549 implicit double precision (a-h,o-z)
3550 common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet,
3551 1 xmu,mode,modedc,icalc,initf,method,iord,maxord,noncon,iterno,
3552 2 itemno,nosolv,ipostp,iscrch
3553 common /mosarg/ gamma,beta,vto,phi,cox,vbi,xnfs,xnsub,xd,xj,xl,
3554 1 xlamda,utra,uexp,vbp,von,vdsat,theta,vcrit,vtra,gleff,cdrain
3555 common /knstnt/ twopi,xlog2,xlog10,root2,rad,boltz,charge,ctok,
3556 1 gmin,reltol,abstol,vntol,trtol,chgtol,eps0,epssil,epsox
3557 iflag=1
3558 if(mode.ne.1) go to 5
3559 iflag=0
3560 if(modedc.eq.2.and.nosolv.ne.0) iflag=1
3561 if(initf.eq.4) iflag=1
3562 5 vd=dmax1(phi-vbd,1.0d-8)
3563 vg=vgb-vbi+phi
3564 vs=dmax1(phi-vbs,1.0d-8)
3565 vsp5=dsqrt(vs)
3566 gammad=gamma
3567 if(gamma.eq.0.0d0) go to 7
3568 if(xj.eq.0.0d0) go to 7
3569 arg=dsqrt(1.0d0+xd*2.0d0*vsp5/xj)
3570 gfact=1.0d0-xj/xl*(arg-1.0d0)
3571 gfact=dmax1(0.5d0,gfact)
3572 gammad=gamma*gfact
3573 7 vth=gammad*vsp5+vs
3574c.. von is referenced to vgb for 'fetlim'
3575c.. change mosfet to reference to vgs (von=vth+vbi-vs) for
3576c.. printing
3577 von=vth+vbi-phi
3578 vdsat=0.0d0
3579 if(vg.lt.vth) go to 100
3580c
3581c 'on' region (linear and saturated)
3582c
3583 gamma2=gammad*0.5d0
3584 sqarg=dsqrt(gamma2*gamma2+vg)
3585 vsat=(sqarg-gamma2)**2
3586 vsatcl=vsat
3587 vs2=vs*vs
3588 vs3=vs2*vs
3589 vs5=vs3*vs2
3590 vs1p5=vs*vsp5
3591 vs2p5=vs1p5*vs
3592 if(vcrit.eq.0.0d0) go to 9
3593c...... iterate to new vsat
3594 iter=1
3595 8 ve2=vsat*vsat
3596c write(6,8878) iter,vsat
3597 8878 format(' iter vsat ',i4,1pd11.2)
3598 vep5=dsqrt(vsat)
3599 ve1p5=vsat*vep5
3600 arg2=0.5d0*(ve2-vs2)
3601 arg1p5=gammad*(ve1p5-vs1p5)/1.5d0
3602 cdrain=vg*(vsat-vs)-arg1p5-arg2
3603 didve=vg-gammad*vep5-vsat
3604 d2idve=-0.5d0*gammad/dmax1(vep5,1.0d-5)-1.0d0
3605 if(vtra.eq.0.0d0) go to 88
3606 trafac=1.0d0/(1.0d0+(vsat-vs)/vtra)
3607 dtrdve=-trafac*trafac/vtra
3608 d2idve=d2idve*trafac+(dtrdve+dtrdve)*(didve-cdrain*trafac/vtra)
3609 didve=didve*trafac+dtrdve*cdrain
3610 cdrain=cdrain*trafac
3611 88 delv=(didve*vcrit-cdrain)/dabs(didve-vcrit*d2idve)
3612c.. limit voltage excursion to 1/2 old vsat
3613 if(dabs(delv).gt.0.5d0*vsat) delv=vsat*dsign(0.5d0,delv)
3614 vsat=vsat+delv
3615c vsat=dmax1(vsat,1.0d-5)
3616c vsat=dmin1(vsat,vsatcl)
3617 if(dabs(delv).lt.1.0d-6) go to 9
3618 iter=iter+1
3619 if(iter.gt.20) write(6,7777) vg,vs,vsat,delv
3620 7777 format(' iteration count for vsat hit limit of 20'/,
3621 1 ' vg vs vsat delv ',1p4d11.3)
3622 if(iter.gt.20) go to 9
3623 go to 8
3624c.. end of iteration loop
3625c.. vdsat is referenced to vds for printing only
3626 9 vdsat=vsat-vs
3627 if(vsat.gt.vsatcl) write(6,9989) vsat,vsatcl
3628 9989 format(' ********error****** vsat is larger than classical vsat',/
3629 1' vsat ',1pd11.3,' classical vsat ',d11.3)
3630 9999 format(' vsat ',1pd11.3,' classical vsat ',d11.3)
3631 if(vd.ge.vsat) go to 10
3632c.. linear region
3633 ve=vd
3634 dvedvd=1.0d0
3635 dvedvg=0.0d0
3636 go to 15
3637c.. saturated region
3638 10 ve=vsat
3639 dvedvd=0.0d0
3640c**************** zero dvedvg!!!
3641c dvedvg=1.0d0-gamma2/sqarg
3642 dvedvg=0.0d0
3643c
3644 15 ve2=ve*ve
3645 ve3=ve2*ve
3646 ve5=ve3*ve2
3647 vep5=dsqrt(ve)
3648 ve1p5=ve*vep5
3649 ve2p5=ve1p5*ve
3650 arg2=0.5d0*(ve2-vs2)
3651 arg1p5=gammad*(ve1p5-vs1p5)/1.5d0
3652 cdrain=vg*(ve-vs)-arg1p5-arg2
3653 didve=vg-gammad*vep5-ve
3654 didvg=ve-vs+didve*dvedvg
3655 didvs=-vg+gammad*vsp5+vs
3656 if(iflag.eq.0) go to 30
3657 if(dabs(cdrain).gt.1.0d-5) go to 20
3658c .. special case when ve almost equals vs and regular formulas don't work
3659c write(6,5475) time,cdrain
3660 5475 format(' time = ',1pd15.6,' cdrain = ',e11.3)
3661 qg=cox*(vg-vs)
3662 ccgg=cox
3663 ccgd=-0.5d0*cox
3664 ccgs=ccgd
3665 qb=-cox*gammad*vsp5
3666 ccbg=0.0d0
3667 ccbd=-cox*0.25d0*gammad/dmax1(vsp5,1.0d-2)
3668 ccbs=ccbd
3669 go to 30
3670c
3671 20 arg2p5=gammad*0.4d0*(ve2p5-vs2p5)
3672 varg=(vg*arg2-arg2p5-(ve3-vs3)/3.0d0)/cdrain
3673 qg=cox*(vg-varg)
3674 dqgdve=cox/cdrain*(varg-ve)*didve
3675 ccgg=cox*(1.0d0-(arg2-varg*(ve-vs))/cdrain)
3676 1 +dqgdve*dvedvg
3677 ccgd=dqgdve*dvedvd
3678 ccgs=cox/cdrain*(varg-vs)*didvs
3679 qb=-cox/cdrain*(vg*arg1p5-gammad*gammad*arg2-arg2p5)
3680 dqbdve=-cox/cdrain*(gammad*vep5+qb/cox)*didve
3681 ccbd=dqbdve*dvedvd
3682 ccbs=-cox/cdrain*(gammad*vsp5+qb/cox)*didvs
3683 ccbg=-cox/cdrain*(arg1p5+qb/cox*(ve-vs))
3684 1 +dqbdve*dvedvg
3685c.. mobility factor (a-la bdm)
3686 30 if(uexp.eq.0.0d0) go to 35
3687 vdenom=vg-vth-utra*(ve-vs)
3688 if(vdenom.le.vbp) go to 45
3689 arg=vbp/vdenom
3690 ufact=dexp(uexp*dlog(arg))
3691 dcoef=-uexp*ufact*arg/vbp
3692 didvg=ufact*didvg+cdrain*dcoef
3693 didvs=ufact*didvs-cdrain*dcoef*(0.5d0*gammad/vsp5+1.0d0-utra)
3694 didve=ufact*didve-cdrain*dcoef*utra
3695 cdrain=cdrain*ufact
3696 go to 45
3697c.. lateral field effects
3698 35 if(vtra.eq.0.0d0) go to 40
3699 trafac=1.0d0/(1.0d0+(ve-vs)/vtra)
3700 dtrdve=-trafac*trafac/vtra
3701 didve=didve*trafac+dtrdve*cdrain
3702c.. note that dtrdvs=-dtrdve
3703 didvs=didvs*trafac-dtrdve*cdrain
3704 didvg=didvg*trafac
3705 cdrain=cdrain*trafac
3706c.. mobility variation a-la sun-daseking
3707 40 if(theta.eq.0.0d0) go to 45
3708 ufact=1.0d0/(1.0d0+theta*(vg-vth))
3709 dufact=-theta*ufact*ufact
3710 didve=ufact*didve
3711 didvg=ufact*didvg+cdrain*dufact
3712 didvs=ufact*didvs-cdrain*dufact*(0.5d0*gammad/vsp5+1.0d0)
3713 cdrain=cdrain*ufact
3714c .. done with 've', use it
3715 45 didvd=didve*dvedvd
3716c.. channel length modulation
3717 if(vcrit.ne.0.0d0) go to 80
3718 if (xlamda.gt.0.0d0) go to 50
3719 if (xnsub.eq.0.0d0) go to 50
3720c.. frohman-grove (lousy) formulation modified a-la newton
3721 arg1=(vd-vsat)/4.0d0
3722 arg2=dsqrt(1.0d0+arg1*arg1)
3723 arg3=dsqrt(arg1+arg2)
3724 clfact=1.0d0/(1.0d0-xd/xl*arg3)
3725 if(clfact.le.0.0d0) go to 60
3726 dclfct=0.125d0*clfact*clfact*xd/xl*(1.0d0+arg1/arg2)/arg3
3727 didvd=clfact*didvd+cdrain*dclfct
3728 didvg=clfact*didvg-cdrain*dclfct*(1.0d0-gamma2/sqarg)
3729 didvs=clfact*didvs
3730 cdrain=cdrain*clfact
3731 go to 200
3732c.. simple (1+vds*lambda/l) formulation
3733 50 xlfact=xlamda/xl
3734 clfact=1.0d0+xlfact*(vd-vs)
3735 didvd=clfact*didvd+cdrain*xlfact
3736 didvs=clfact*didvs-cdrain*xlfact
3737 didvg=clfact*didvg
3738 cdrain=cdrain*clfact
3739 go to 200
3740c.. device punched thru
3741 60 clfact=1000.0d0
3742 if(ipunch.gt.50) go to 200
3743 ipunch=ipunch+1
3744 write(6,61)
3745 61 format('0warning: channel length reduced to zero in mosfet')
3746 go to 200
3747c.. into saturation with vcrit ne 0
3748 80 if(vd.le.vsat) go to 200
3749 xk1=vcrit/2.0d0/xl/gleff
3750 temp=dsqrt(xk1*xk1+(vd-vsat)/gleff)
3751 clfact=1.0d0+(temp-xk1)/xl
3752 dclfct=0.5d0/xl/temp/gleff
3753 didvd=didvd*clfact+cdrain*dclfct
3754 didvs=didvs*clfact
3755 didvg=didvg*clfact
3756 cdrain=cdrain*clfact
3757c
3758 go to 200
3759c
3760c.. cut-off region (vg<vth)
3761c
3762 100 continue
3763 cdrain=0.0d0
3764 didvg=0.0d0
3765 didvd=0.0d0
3766 didvs=0.0d0
3767 if(iflag.eq.0) return
3768 if(vg.gt.0.0d0) go to 120
3769 qg=cox*vg
3770 ccgg=cox
3771 go to 130
3772 120 gamma2=gammad*0.5d0
3773 sqarg=dsqrt(gamma2*gamma2+vg)
3774 qg=gammad*cox*(sqarg-gamma2)
3775 ccgg=0.5d0*cox*gammad/sqarg
3776 130 qb=-qg
3777 ccbg=-ccgg
3778 ccgd=0.0d0
3779 ccgs=0.0d0
3780 ccbd=0.0d0
3781 ccbs=0.0d0
3782 200 qc=-(qg+qb)
3783c write(6,7657) time,vg,vs,ve,qg,ccgg,ccgd,ccgs,qb,ccbg,ccbd,ccbs
3784 7657 format(' time = ',1pd12.5,' vg vs ve ',3d11.3,
3785 1 /,' qg ',4d11.3,' qb ',4d11.3)
3786 return
3787 end
3788 subroutine gasfet(loc)
3789c *** a gasfet (or any other) model may ***
3790c * be inserted here... if you happen *
3791c * to have a good one! *
3792c *** ***
3793c
3794 return
3795 end
Continuous source code history from original PDP-7 UNIX to FreeBSD 2, the first fully unencumbered FreeBSD release.