1 /***************************************************************************
2 JSPICE3 adaptation of Spice3e2 - Copyright (c) Stephen R. Whiteley 1992
3 Copyright 1990 Regents of the University of California.  All rights reserved.
4 Authors: 1988 Thomas L. Quarles
5          1993 Stephen R. Whiteley
6 ****************************************************************************/
7 
8 #define ANAL_EXT
9 
10 /* subroutine to do Transfer Function analysis */
11 
12 #include "spice.h"
13 #include <stdio.h>
14 #include <math.h>
15 #include "srcdefs.h"
16 #include "util.h"
17 #include "sperror.h"
18 #include "outdata.h"
19 #include "tfdefs.h"
20 #include "cktext.h"
21 
22 #ifdef __STDC__
23 static int tf_acoperation(CKTcircuit*,struct sOUTdata*,int);
24 static int tf_dcoperation(CKTcircuit*,struct sOUTdata*,int);
25 #else
26 static int tf_acoperation();
27 static int tf_dcoperation();
28 #endif
29 
30 
31 /* ARGSUSED */
32 int
TFan(cktp,restart)33 TFan(cktp,restart)
34 
35 GENERIC *cktp;
36 int restart;    /* forced restart flag */
37 {
38     CKTcircuit *ckt = (CKTcircuit *)cktp;
39     static struct sOUTdata outd;
40     int error;
41     int i;
42     IFuid uids[3];
43     int code;
44     TFAN *job = (TFAN*)ckt->CKTcurJob;
45 
46     if (ckt->CKTjjPresent) {
47         (*(SPfrontEnd->IFerror))(ERR_FATAL,
48             "Transfer function not possible with Josephson junctions", NULL);
49             return (OK);
50     }
51 
52     if (restart) {
53 
54         /* make sure sources are present */
55         code = CKTtypelook("Source");
56         error = CKTfndDev((GENERIC*)ckt,&code,&job->TFinSrcDev,
57             job->TFinSrc, (GENERIC *)NULL,(IFuid)NULL);
58         if (error != 0) {
59             (*(SPfrontEnd->IFerror))(ERR_FATAL,
60                 "Transfer function source %s not in circuit",
61                 &job->TFinSrc);
62             return (E_NOTFOUND);
63         }
64         if (((SRCinstance*)job->TFinSrcDev)->SRCtype == SRC_V) {
65             job->TFinIsV = 1;
66             job->TFinIsI = 0;
67         }
68         else {
69             job->TFinIsI = 1;
70             job->TFinIsV = 0;
71         }
72 
73         if (job->TFoutIsI) {
74             error = CKTfndDev((GENERIC*)ckt,&code,&job->TFoutSrcDev,
75                 job->TFoutSrc,(GENERIC *)NULL,(IFuid)NULL);
76             if (error != 0) {
77                 (*(SPfrontEnd->IFerror))(ERR_FATAL,
78                     "Transfer function source %s not in circuit",
79                     &job->TFoutSrc);
80                 return (E_NOTFOUND);
81             }
82             if (((SRCinstance*)job->TFoutSrcDev)->SRCtype != SRC_V) {
83                 (*(SPfrontEnd->IFerror))(ERR_FATAL,
84                     "Transfer function source %s not a voltage source",
85                     &job->TFoutSrc);
86                 return (E_NOTFOUND);
87             }
88         }
89 
90         if (job->DC.eltName[0] != NULL) {
91             /* DC source was given */
92             SRCinstance *here;
93 
94             here = (SRCinstance*)NULL;
95             for (i = 0; i <= job->DC.nestLevel; i++) {
96                 error = CKTfndDev((GENERIC*)ckt,&code,(GENERIC**)&here,
97                     job->DC.eltName[i], (GENERIC *)NULL,(IFuid)NULL);
98 
99                 if (error) {
100                     (*(SPfrontEnd->IFerror))(ERR_FATAL,
101                         "DCtrCurv: source %s not in circuit",
102                         &(job->DC.eltName[i]));
103                     return (E_NOTFOUND);
104                 }
105                 job->DC.elt[i] = (GENinstance*)here;
106                 job->DC.vsave[i] = here->SRCdcValue;
107             }
108         }
109 
110         /* make a UID for the transfer function output */
111         (*(SPfrontEnd->IFnewUid))(ckt,&uids[0],(IFuid)NULL,"tranfunc",
112                 UID_OTHER,(GENERIC **)NULL);
113 
114         /* make a UID for the input impedance */
115         (*(SPfrontEnd->IFnewUid))(ckt,&uids[1],job->TFinSrc,
116                 "Zi", UID_OTHER,(GENERIC **)NULL);
117 
118         /* make a UID for the output impedance */
119         if (job->TFoutIsI) {
120             (*(SPfrontEnd->IFnewUid))
121                 (ckt,&uids[2],job->TFoutSrc,
122                 "Zo", UID_OTHER,(GENERIC **)NULL);
123         }
124         else {
125             char *name = (char *)
126                 MALLOC(sizeof(char)*(strlen(job->TFoutName)+22));
127             char *s = strchr(job->TFoutName,'(');
128             if (s)
129                 (void)sprintf(name,"Zo%s",s);
130             else
131                 (void)sprintf(name,"Zo(%s)",job->TFoutName);
132             (*(SPfrontEnd->IFnewUid))(ckt,&uids[2],(IFuid)NULL,
133                     name, UID_OTHER,(GENERIC **)NULL);
134         }
135 
136         if (job->AC.stepType == DCSTEP) {
137             outd.dataType = IF_REAL;
138             outd.refName = (IFuid)NULL;
139             outd.refType = 0;
140         }
141         else {
142             outd.dataType = IF_COMPLEX;
143             (*(SPfrontEnd->IFnewUid))((GENERIC *)ckt,&outd.refName,
144                 (IFuid)NULL, "frequency", UID_OTHER,(GENERIC **)NULL);
145             outd.refType = IF_REAL;
146         }
147 
148         outd.circuitPtr = (GENERIC *)ckt;
149         outd.analysisPtr = (GENERIC*)ckt->CKTcurJob;
150         outd.analName = ckt->CKTcurJob->JOBname;
151         outd.numNames = 3;
152         outd.dataNames = uids;
153         outd.plotPtr = &job->TFplot;
154         outd.numPts = 1;
155         outd.initValue = 0;
156         outd.finalValue = 0;
157         outd.step = 0;
158         outd.count = 0;
159         (*(SPfrontEnd->OUTbeginPlot))((GENERIC*)&outd);
160     }
161     error = DCTloop(tf_dcoperation,ckt,restart,&job->DC,&outd);
162 
163     if (error == E_PAUSE)
164         return (error);
165 
166     (*(SPfrontEnd->OUTendPlot))(*outd.plotPtr);
167     return (OK);
168 }
169 
170 
171 static int
tf_dcoperation(ckt,outd,restart)172 tf_dcoperation(ckt,outd,restart)
173 
174 CKTcircuit *ckt;
175 struct sOUTdata *outd;
176 int restart;
177 {
178     TFAN *job = (TFAN*)ckt->CKTcurJob;
179     int error;
180     int size;
181     int i;
182     double outputs[3];
183     double A;
184     IFvalue outdata;
185     IFvalue refval;
186     SRCinstance *ptr;
187 
188     /* find the operating point */
189 
190     error = CKTic(ckt);
191     if (error)
192         return (error);
193 
194     error = CKTop(ckt,
195             MODEDCOP | MODEINITJCT,
196             MODEDCOP | MODEINITFLOAT,
197             ckt->CKTdcMaxIter);
198     if (error)
199         return (error);
200 
201     if (job->AC.stepType != DCSTEP) {
202         error = ACloop(tf_acoperation,ckt,restart,&job->AC,outd);
203         return (error);
204     }
205 
206     size = spGetSize(ckt->CKTmatrix,1);
207 
208     for (i = 0; i <= size; i++) {
209         ckt->CKTrhs[i] = 0;
210     }
211 
212     ptr = (SRCinstance*)job->TFinSrcDev;
213     if (job->TFinIsI) {
214         ckt->CKTrhs[ptr->SRCposNode] -= 1;
215         ckt->CKTrhs[ptr->SRCnegNode] += 1;
216     }
217     else {
218         ckt->CKTrhs[ptr->SRCbranch] += 1;
219     }
220 
221     spSolve(ckt->CKTmatrix,ckt->CKTrhs,ckt->CKTrhs,NULL,NULL);
222     ckt->CKTrhs[0] = 0;
223 
224     /* find transfer function */
225     if (job->TFoutIsV) {
226         outputs[0] =
227             ckt->CKTrhs[job->TFoutPos->number] -
228                 ckt->CKTrhs[job->TFoutNeg->number] ;
229     }
230     else {
231         ptr = (SRCinstance*)job->TFoutSrcDev;
232         outputs[0] = ckt->CKTrhs[ptr->SRCbranch];
233     }
234 
235     /* now for input impedance */
236     ptr = (SRCinstance*)job->TFinSrcDev;
237     if (job->TFinIsI) {
238         outputs[1] =
239             ckt->CKTrhs[ptr->SRCnegNode] - ckt->CKTrhs[ptr->SRCposNode];
240     }
241     else {
242         A = ckt->CKTrhs[ptr->SRCbranch];
243         if (FABS(A) < 1e-20)
244             outputs[1] = 1e20;
245         else
246             outputs[1] = -1/A;
247     }
248 
249     if (job->TFoutIsI && (job->TFoutSrc == job->TFinSrc)) {
250         outputs[2] = outputs[1];
251         /* no need to compute output impedance when it is the same as
252          * the input
253          */
254     }
255     else {
256         /* now for output resistance */
257         for (i = 0; i <= size; i++) {
258             ckt->CKTrhs[i] = 0;
259         }
260         if (job->TFoutIsV) {
261             ckt->CKTrhs[job->TFoutPos->number] -= 1;
262             ckt->CKTrhs[job->TFoutNeg->number] += 1;
263         }
264         else {
265             ptr = (SRCinstance*)job->TFoutSrcDev;
266             ckt->CKTrhs[ptr->SRCbranch] += 1;
267         }
268         spSolve(ckt->CKTmatrix,ckt->CKTrhs,ckt->CKTrhs,
269             ckt->CKTirhs,ckt->CKTirhs);
270         ckt->CKTrhs[0] = 0;
271 
272         if (job->TFoutIsV) {
273             outputs[2] =
274                 ckt->CKTrhs[job->TFoutNeg->number] -
275                 ckt->CKTrhs[job->TFoutPos->number];
276         }
277         else {
278             A = ckt->CKTrhs[ptr->SRCbranch];
279             if (FABS(A) < 1e-20)
280                 outputs[2] = 1e20;
281             else
282                 outputs[2] = 1/A;
283         }
284     }
285 
286     outdata.v.numValue = 3;
287     outdata.v.vec.rVec = outputs;
288     refval.rValue = ckt->CKTomega/(2*M_PI);
289     (*(SPfrontEnd->OUTdata))(*outd->plotPtr,&refval,&outdata);
290     outd->count++;
291     return (OK);
292 }
293 
294 
295 static int
tf_acoperation(ckt,outd,restart)296 tf_acoperation(ckt,outd,restart)
297 
298 CKTcircuit *ckt;
299 struct sOUTdata *outd;
300 int restart;
301 {
302     int size;
303     double A, B, MAG;
304     IFcomplex coutputs[3];
305     IFvalue outdata;
306     IFvalue refval;
307     int error;
308     int i;
309     SRCinstance *ptr;
310     TFAN *job = (TFAN*)ckt->CKTcurJob;
311 
312     ckt->CKTmode = MODEDCOP | MODEINITSMSIG;
313     error = CKTload(ckt);
314     if (error)
315         return (error);
316     ckt->CKTmode = MODEAC;
317     error = NIacIter(ckt);
318     if (error)
319         return (error);
320 
321     size = spGetSize(ckt->CKTmatrix,1);
322 
323     for (i = 0; i <= size; i++) {
324         ckt->CKTrhs[i] = 0;
325         ckt->CKTirhs[i] = 0;
326     }
327 
328     ptr = (SRCinstance*)job->TFinSrcDev;
329     if (job->TFinIsI) {
330         ckt->CKTrhs[ptr->SRCposNode] -= 1;
331         ckt->CKTrhs[ptr->SRCnegNode] += 1;
332     }
333     else {
334         ckt->CKTrhs[ptr->SRCbranch] += 1;
335     }
336 
337     spSolve(ckt->CKTmatrix,ckt->CKTrhs,ckt->CKTrhs,
338         ckt->CKTirhs,ckt->CKTirhs);
339     ckt->CKTrhs[0] = 0;
340     ckt->CKTirhs[0] = 0;
341 
342     /* find transfer function */
343     if (job->TFoutIsV) {
344         coutputs[0].real =
345             ckt->CKTrhs[job->TFoutPos->number] -
346                 ckt->CKTrhs[job->TFoutNeg->number] ;
347         coutputs[0].imag =
348             ckt->CKTirhs[job->TFoutPos->number] -
349                 ckt->CKTirhs[job->TFoutNeg->number] ;
350     }
351     else {
352         ptr = (SRCinstance*)job->TFoutSrcDev;
353         coutputs[0].real = ckt->CKTrhs[ptr->SRCbranch];
354         coutputs[0].imag = ckt->CKTirhs[ptr->SRCbranch];
355     }
356 
357     /* now for input impedance */
358     ptr = (SRCinstance*)job->TFinSrcDev;
359     if (job->TFinIsI) {
360         coutputs[1].real =
361             ckt->CKTrhs[ptr->SRCnegNode] - ckt->CKTrhs[ptr->SRCposNode];
362         coutputs[1].imag =
363             ckt->CKTirhs[ptr->SRCnegNode] - ckt->CKTirhs[ptr->SRCposNode];
364     }
365     else {
366         A = ckt->CKTrhs[ptr->SRCbranch];
367         B = ckt->CKTirhs[ptr->SRCbranch];
368         MAG = A*A + B*B;
369         if (MAG < 1e-40)
370             MAG = 1e-40;
371         coutputs[1].real = -A/MAG;
372         coutputs[1].imag = B/MAG;
373     }
374 
375     if (job->TFoutIsI && (job->TFoutSrc == job->TFinSrc)) {
376         coutputs[2].real = coutputs[1].real;
377         coutputs[2].imag = coutputs[1].imag;
378         /* no need to compute output impedance when it is the same as
379          * the input
380          */
381     }
382     else {
383         /* now for output resistance */
384         for (i = 0; i <= size; i++) {
385             ckt->CKTrhs[i] = 0;
386             ckt->CKTirhs[i] = 0;
387         }
388         if (job->TFoutIsV) {
389             ckt->CKTrhs[job->TFoutPos->number] -= 1;
390             ckt->CKTrhs[job->TFoutNeg->number] += 1;
391         }
392         else {
393             ptr = (SRCinstance*)job->TFoutSrcDev;
394             ckt->CKTrhs[ptr->SRCbranch] += 1;
395         }
396         spSolve(ckt->CKTmatrix,ckt->CKTrhs,ckt->CKTrhs,
397             ckt->CKTirhs,ckt->CKTirhs);
398         ckt->CKTrhs[0] = 0;
399         ckt->CKTirhs[0] = 0;
400 
401         if (job->TFoutIsV) {
402             coutputs[2].real =
403                 ckt->CKTrhs[job->TFoutNeg->number] -
404                 ckt->CKTrhs[job->TFoutPos->number];
405             coutputs[2].imag =
406                 ckt->CKTirhs[job->TFoutNeg->number] -
407                 ckt->CKTirhs[job->TFoutPos->number];
408         }
409         else {
410             A = ckt->CKTrhs[ptr->SRCbranch];
411             B = ckt->CKTirhs[ptr->SRCbranch];
412             MAG = A*A + B*B;
413             if (MAG < 1e-40)
414                 MAG = 1e-40;
415             coutputs[2].real = A/MAG;
416             coutputs[2].imag = -B/MAG;
417         }
418     }
419 
420     outdata.v.numValue = 3;
421     outdata.v.vec.cVec = coutputs;
422     refval.rValue = ckt->CKTomega/(2*M_PI);
423     (*(SPfrontEnd->OUTdata))(*outd->plotPtr,&refval,&outdata);
424     outd->count++;
425     return (OK);
426 }
427