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: 1985 Wayne A. Christopher
5 1992 Stephen R. Whiteley
6 ****************************************************************************/
7
8 /*
9 *
10 * Routines to do complex mathematical functions. These routines require
11 * the -lm libraries. We sacrifice a lot of space to be able
12 * to avoid having to do a seperate call for every vector element,
13 * but it pays off in time savings. These routines should never
14 * allow FPE's to happen.
15 *
16 * Complex functions are called as follows:
17 * cx_something(data, type, length, &newlength, &newtype),
18 * and return a char * that is cast to complex or double.
19 *
20 */
21
22 #include "spice.h"
23 #include "ftedefs.h"
24 #include "ftecmath.h"
25
26 /* This flag determines whether degrees or radians are used. The radtodeg
27 * and degtorad macros are no-ops if this is false.
28 */
29
30 bool cx_degrees = false;
31
32 char *
cx_mag(data,type,length,newlength,newtype)33 cx_mag(data, type, length, newlength, newtype)
34
35 char *data;
36 short type;
37 int length;
38 int *newlength;
39 short *newtype;
40 {
41 double *d = alloc_d(length);
42 double *dd = (double *) data;
43 complex *cc = (complex *) data;
44 int i;
45
46 *newlength = length;
47 *newtype = VF_REAL;
48 if (type == VF_REAL)
49 for (i = 0; i < length; i++)
50 d[i] = FTEcabs(dd[i]);
51 else
52 for (i = 0; i < length; i++)
53 d[i] = cmag(&cc[i]);
54 return ((char *) d);
55 }
56
57
58 char *
cx_ph(data,type,length,newlength,newtype)59 cx_ph(data, type, length, newlength, newtype)
60
61 char *data;
62 short type;
63 int length;
64 int *newlength;
65 short *newtype;
66 {
67 double *d = alloc_d(length);
68 complex *cc = (complex *) data;
69 int i;
70
71 *newlength = length;
72 *newtype = VF_REAL;
73 if (type == VF_COMPLEX)
74 for (i = 0; i < length; i++) {
75 d[i] = radtodeg(cph(&cc[i]));
76 }
77 /* Otherwise it is 0, but tmalloc zeros the stuff already. */
78 return ((char *) d);
79 }
80
81 /* If this is pure imaginary we might get real, but never mind... */
82
83 char *
cx_j(data,type,length,newlength,newtype)84 cx_j(data, type, length, newlength, newtype)
85
86 char *data;
87 short type;
88 int length;
89 int *newlength;
90 short *newtype;
91 {
92 complex *c = alloc_c(length);
93 complex *cc = (complex *) data;
94 double *dd = (double *) data;
95 int i;
96
97 *newlength = length;
98 *newtype = VF_COMPLEX;
99 if (type == VF_COMPLEX)
100 for (i = 0; i < length; i++) {
101 realpart(&c[i]) = - imagpart(&cc[i]);
102 imagpart(&c[i]) = realpart(&cc[i]);
103 }
104 else
105 for (i = 0; i < length; i++) {
106 imagpart(&c[i]) = dd[i];
107 /* Real part is already 0. */
108 }
109 return ((char *) c);
110 }
111
112
113 char *
cx_real(data,type,length,newlength,newtype)114 cx_real(data, type, length, newlength, newtype)
115
116 char *data;
117 short type;
118 int length;
119 int *newlength;
120 short *newtype;
121 {
122 double *d = alloc_d(length);
123 double *dd = (double *) data;
124 complex *cc = (complex *) data;
125 int i;
126
127 *newlength = length;
128 *newtype = VF_REAL;
129 if (type == VF_COMPLEX)
130 for (i = 0; i < length; i++)
131 d[i] = realpart(&cc[i]);
132 else
133 for (i = 0; i < length; i++)
134 d[i] = dd[i];
135 return ((char *) d);
136 }
137
138
139 char *
cx_imag(data,type,length,newlength,newtype)140 cx_imag(data, type, length, newlength, newtype)
141
142 char *data;
143 short type;
144 int length;
145 int *newlength;
146 short *newtype;
147 {
148 double *d = alloc_d(length);
149 double *dd = (double *) data;
150 complex *cc = (complex *) data;
151 int i;
152
153 *newlength = length;
154 *newtype = VF_REAL;
155 if (type == VF_COMPLEX)
156 for (i = 0; i < length; i++)
157 d[i] = imagpart(&cc[i]);
158 else
159 for (i = 0; i < length; i++)
160 d[i] = dd[i];
161 return ((char *) d);
162 }
163
164 /* This is obsolete... */
165
166 char *
cx_pos(data,type,length,newlength,newtype)167 cx_pos(data, type, length, newlength, newtype)
168
169 char *data;
170 short type;
171 int length;
172 int *newlength;
173 short *newtype;
174 {
175 double *d = alloc_d(length);
176 double *dd = (double *) data;
177 complex *cc = (complex *) data;
178 int i;
179
180 *newlength = length;
181 *newtype = VF_REAL;
182 if (type == VF_COMPLEX)
183 for (i = 0; i < length; i++)
184 d[i] = ((realpart(&cc[i]) > 0.0) ? 1.0 : 0.0);
185 else
186 for (i = 0; i < length; i++)
187 d[i] = ((dd[i] > 0.0) ? 1.0 : 0.0);
188 return ((char *) d);
189 }
190
191
192 char *
cx_db(data,type,length,newlength,newtype)193 cx_db(data, type, length, newlength, newtype)
194
195 char *data;
196 short type;
197 int length;
198 int *newlength;
199 short *newtype;
200 {
201 double *d = alloc_d(length);
202 double *dd = (double *) data;
203 complex *cc = (complex *) data;
204 double tt;
205 int i;
206
207 *newlength = length;
208 *newtype = VF_REAL;
209 if (type == VF_COMPLEX)
210 for (i = 0; i < length; i++) {
211 tt = cmag(&cc[i]);
212 rcheck(tt >= 0, "db");
213 if (tt == 0.0)
214 d[i] = 20.0 * - log(HUGE);
215 else
216 d[i] = 20.0 * log10(tt);
217 }
218 else
219 for (i = 0; i < length; i++) {
220 rcheck(dd[i] >= 0, "db");
221 if (dd[i] == 0.0)
222 d[i] = 20.0 * - log(HUGE);
223 else
224 d[i] = 20.0 * log10(dd[i]);
225 }
226 return ((char *) d);
227 }
228
229
230 char *
cx_log(data,type,length,newlength,newtype)231 cx_log(data, type, length, newlength, newtype)
232
233 char *data;
234 short type;
235 int length;
236 int *newlength;
237 short *newtype;
238 {
239 double *d, td;
240 complex *c;
241 double *dd = (double *) data;
242 complex *cc = (complex *) data;
243 int i;
244
245 if (type == VF_REAL) {
246 d = alloc_d(length);
247 *newtype = VF_REAL;
248 }
249 else {
250 c = alloc_c(length);
251 *newtype = VF_COMPLEX;
252 }
253 *newlength = length;
254 if (type == VF_COMPLEX) {
255 for (i = 0; i < length; i++) {
256 td = cmag(&cc[i]);
257 /* Perhaps we should trap when td = 0.0, but Ken wants
258 * this to be possible...
259 */
260 rcheck(td >= 0, "log");
261 if (td == 0.0) {
262 realpart(&c[i]) = - log10(HUGE);
263 imagpart(&c[i]) = 0.0;
264 }
265 else {
266 realpart(&c[i]) = log10(td);
267 imagpart(&c[i]) = atan2(imagpart(&cc[i]),
268 realpart(&cc[i]));
269 }
270 }
271 return ((char *) c);
272 }
273 else {
274 for (i = 0; i < length; i++) {
275 rcheck(dd[i] >= 0, "log");
276 if (dd[i] == 0.0)
277 d[i] = - log10(HUGE);
278 else
279 d[i] = log10(dd[i]);
280 }
281 return ((char *) d);
282 }
283 }
284
285
286 char *
cx_ln(data,type,length,newlength,newtype)287 cx_ln(data, type, length, newlength, newtype)
288
289 char *data;
290 short type;
291 int length;
292 int *newlength;
293 short *newtype;
294 {
295 double *d, td;
296 complex *c;
297 double *dd = (double *) data;
298 complex *cc = (complex *) data;
299 int i;
300
301 if (type == VF_REAL) {
302 d = alloc_d(length);
303 *newtype = VF_REAL;
304 }
305 else {
306 c = alloc_c(length);
307 *newtype = VF_COMPLEX;
308 }
309 *newlength = length;
310 if (type == VF_COMPLEX) {
311 for (i = 0; i < length; i++) {
312 td = cmag(&cc[i]);
313 rcheck(td >= 0, "ln");
314 if (td == 0.0) {
315 realpart(&c[i]) = - log(HUGE);
316 imagpart(&c[i]) = 0.0;
317 }
318 else {
319 realpart(&c[i]) = log(td);
320 imagpart(&c[i]) = atan2(imagpart(&cc[i]),
321 realpart(&cc[i]));
322 }
323 }
324 return ((char *) c);
325 }
326 else {
327 for (i = 0; i < length; i++) {
328 rcheck(dd[i] >= 0, "ln");
329 if (dd[i] == 0.0)
330 d[i] = - log(HUGE);
331 else
332 d[i] = log(dd[i]);
333 }
334 return ((char *) d);
335 }
336 }
337
338
339 char *
cx_exp(data,type,length,newlength,newtype)340 cx_exp(data, type, length, newlength, newtype)
341
342 char *data;
343 short type;
344 int length;
345 int *newlength;
346 short *newtype;
347 {
348 double *d, td;
349 complex *c;
350 double *dd = (double *) data;
351 complex *cc = (complex *) data;
352 int i;
353
354 if (type == VF_REAL) {
355 d = alloc_d(length);
356 *newtype = VF_REAL;
357 }
358 else {
359 c = alloc_c(length);
360 *newtype = VF_COMPLEX;
361 }
362 *newlength = length;
363 if (type == VF_COMPLEX) {
364 for (i = 0; i < length; i++) {
365 td = exp(realpart(&cc[i]));
366 realpart(&c[i]) = td * cos(realpart(&cc[i]));
367 imagpart(&c[i]) = td * sin(imagpart(&cc[i]));
368 }
369 return ((char *) c);
370 }
371 else {
372 for (i = 0; i < length; i++)
373 d[i] = exp(dd[i]);
374 return ((char *) d);
375 }
376 }
377
378
379 char *
cx_sqrt(data,type,length,newlength,newtype)380 cx_sqrt(data, type, length, newlength, newtype)
381
382 char *data;
383 short type;
384 int length;
385 int *newlength;
386 short *newtype;
387 {
388 double *d;
389 complex *c;
390 double *dd = (double *) data;
391 complex *cc = (complex *) data;
392 int i, cres = (type == VF_REAL) ? 0 : 1;
393
394 if (type == VF_REAL)
395 for (i = 0; i < length; i++)
396 if (dd[i] < 0.0)
397 cres = 1;
398 if (cres) {
399 c = alloc_c(length);
400 *newtype = VF_COMPLEX;
401 }
402 else {
403 d = alloc_d(length);
404 *newtype = VF_REAL;
405 }
406 *newlength = length;
407 if (type == VF_COMPLEX) {
408 for (i = 0; i < length; i++) {
409 if (realpart(&cc[i]) == 0.0) {
410 if (imagpart(&cc[i]) == 0.0) {
411 realpart(&c[i]) = 0.0;
412 imagpart(&c[i]) = 0.0;
413 }
414 else if (imagpart(&cc[i]) > 0.0) {
415 realpart(&c[i]) = sqrt (0.5 *
416 imagpart(&cc[i]));
417 imagpart(&c[i]) = realpart(&c[i]);
418 }
419 else {
420 imagpart(&c[i]) = sqrt( -0.5 *
421 imagpart(&cc[i]));
422 realpart(&c[i]) = - imagpart(&c[i]);
423 }
424 }
425 else if (realpart(&cc[i]) > 0.0) {
426 if (imagpart(&cc[i]) == 0.0) {
427 realpart(&c[i]) =
428 sqrt(realpart(&cc[i]));
429 imagpart(&c[i]) = 0.0;
430 }
431 else if (imagpart(&cc[i]) < 0.0) {
432 realpart(&c[i]) = -sqrt(0.5 *
433 (cmag(&cc[i]) + realpart(&cc[i])));
434 }
435 else {
436 realpart(&c[i]) = sqrt(0.5 *
437 (cmag(&cc[i]) + realpart(&cc[i])));
438 }
439 imagpart(&c[i]) = imagpart(&cc[i]) / (2.0 *
440 realpart(&c[i]));
441 }
442 else { /* realpart(&cc[i]) < 0.0) */
443 if (imagpart(&cc[i]) == 0.0) {
444 realpart(&c[i]) = 0.0;
445 imagpart(&c[i]) =
446 sqrt(- realpart(&cc[i]));
447 }
448 else {
449 if (imagpart(&cc[i]) < 0.0)
450 imagpart(&c[i]) = - sqrt(0.5 *
451 (cmag(&cc[i]) -
452 realpart(&cc[i])));
453 else
454 imagpart(&c[i]) = sqrt(0.5 *
455 (cmag(&cc[i]) -
456 realpart(&cc[i])));
457 realpart(&c[i]) = imagpart(&cc[i]) /
458 (2.0 * imagpart(&c[i]));
459 }
460 }
461 }
462 return ((char *) c);
463 }
464 else if (cres) {
465 for (i = 0; i < length; i++)
466 if (dd[i] < 0.0)
467 imagpart(&c[i]) = sqrt(- dd[i]);
468 else
469 realpart(&c[i]) = sqrt(dd[i]);
470 return ((char *) c);
471 }
472 else {
473 for (i = 0; i < length; i++)
474 d[i] = sqrt(dd[i]);
475 return ((char *) d);
476 }
477 }
478
479
480 char *
cx_sin(data,type,length,newlength,newtype)481 cx_sin(data, type, length, newlength, newtype)
482
483 char *data;
484 short type;
485 int length;
486 int *newlength;
487 short *newtype;
488 {
489 double *d;
490 complex *c;
491 double *dd = (double *) data;
492 complex *cc = (complex *) data;
493 int i;
494
495 if (type == VF_REAL) {
496 d = alloc_d(length);
497 *newtype = VF_REAL;
498 }
499 else {
500 c = alloc_c(length);
501 *newtype = VF_COMPLEX;
502 }
503 *newlength = length;
504 if (type == VF_COMPLEX) {
505 for (i = 0; i < length; i++) {
506 realpart(&c[i]) = sin(degtorad(realpart(&cc[i]))) *
507 cosh(degtorad(imagpart(&cc[i])));
508 imagpart(&c[i]) = cos(degtorad(realpart(&cc[i]))) *
509 sinh(degtorad(imagpart(&cc[i])));
510 }
511 return ((char *) c);
512 }
513 else {
514 for (i = 0; i < length; i++)
515 d[i] = sin(degtorad(dd[i]));
516 return ((char *) d);
517 }
518 }
519
520
521 char *
cx_cos(data,type,length,newlength,newtype)522 cx_cos(data, type, length, newlength, newtype)
523
524 char *data;
525 short type;
526 int length;
527 int *newlength;
528 short *newtype;
529 {
530 double *d;
531 complex *c;
532 double *dd = (double *) data;
533 complex *cc = (complex *) data;
534 int i;
535
536 if (type == VF_REAL) {
537 d = alloc_d(length);
538 *newtype = VF_REAL;
539 }
540 else {
541 c = alloc_c(length);
542 *newtype = VF_COMPLEX;
543 }
544 *newlength = length;
545 if (type == VF_COMPLEX) {
546 for (i = 0; i < length; i++) {
547 realpart(&c[i]) = cos(degtorad(realpart(&cc[i]))) *
548 cosh(degtorad(imagpart(&cc[i])));
549 imagpart(&c[i]) = - sin(degtorad(realpart(&cc[i]))) *
550 sinh(degtorad(imagpart(&cc[i])));
551 }
552 return ((char *) c);
553 }
554 else {
555 for (i = 0; i < length; i++)
556 d[i] = cos(degtorad(dd[i]));
557 return ((char *) d);
558 }
559 }
560
561