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