home *** CD-ROM | disk | FTP | other *** search
/ Peanuts NeXT Software Archives / Peanuts-2.iso / Developer / hardware / dsp / drbub / analog / sspice.c < prev    next >
Encoding:
C/C++ Source or Header  |  1991-09-07  |  27.6 KB  |  935 lines
  1. /***************************************************************************
  2. *        STEVE SPICE            steve punte 10-24-87
  3. *
  4. *    This program is frequency domain only version of spice.  It
  5. *    was developed for the ece-149 class for the simulation of
  6. *    passive and active filters on standard IBM-ATs that normally
  7. *    have insufficient memory for even the reduced version of
  8. *    pspice.  This program is written very flexibbly to request
  9. *    memory from the operating system as it is needed.  Unfortionately
  10. *    the micro-soft C compiler requires that the maximum amount of
  11. *    memory ever needed be stated before hand.  Never the less
  12. *    this is a portable program to any system that compiles C.
  13. *
  14. *****************************************************************************
  15. *
  16. *    This is the header file.  All global variables are defined here.
  17. *    If this file is used in file "main", then the phrase "extern"
  18. *    is stripped of from all declarations.  Otherwise it remains for
  19. *    use in other files.
  20. *
  21. ****************************************************************************/
  22.  
  23. #include    <stdio.h>
  24. #include    <math.h>
  25. #define        TRUE        1
  26. #define        FALSE        0
  27. #define        PIE        3.141592653
  28. #define        STR_LNGTH    20  
  29.  
  30. #define        MAIN
  31. #ifdef        MAIN        /* Removes extern if used in main.c */
  32. #define        extern     
  33. #endif
  34. #define        ALLCODE        /* If all code is present, allows proper compliation*/ 
  35.                 /* Convenient pointers to dynamically allocated arrays */
  36. #define        addm(I, J)    (addm_mtrx_ptr + I + N * J)
  37. #define        curr(I )    (curr_vect_ptr + I )
  38.  
  39. extern struct component
  40.     {
  41.     char name[ STR_LNGTH ];        /* Name of component */
  42.     int node_i;            /* One connection of the component */
  43.     char name_i[ STR_LNGTH ];    /* Name of first node */
  44.     int node_j;            /* The second connection of the component */
  45.     char name_j[ STR_LNGTH ];    /* Name of second node */
  46.     int node_k;            /* A third connection, if existing */
  47.     char name_k[ STR_LNGTH ];    /* Name of third node */
  48.     int node_l;            /* A Fourth connection, if existing */
  49.     char name_l[ STR_LNGTH ];    /* Name of fourth node */
  50.     double value;            /* The value of the component */
  51.     struct component *link;        /* Link foward to next structure or NULL */
  52.     };
  53.  
  54. extern struct component *begin_caps ;    /* Start of capacitor link list */
  55. extern struct component *begin_inds ;    /* Start of inductor link list */
  56. extern struct component *begin_res  ;    /* Start of resistor link list */
  57. extern struct component *begin_vcv  ;    /* Start of VCV link list */
  58. extern struct component *begin_vci  ;    /* Start of VCI link list */
  59. extern struct component *end_caps   ;    /* Last entry of capacitor link list */
  60. extern struct component *end_inds   ;    /* Last entry of inductor link list */
  61. extern struct component *end_res    ;    /* Last entry of resistor link list */
  62. extern struct component *end_vcv    ;    /* Last entry of VCV link list */
  63. extern struct component *end_vci    ;    /* Last entry of VCI link list */
  64.  
  65. extern struct component VS_source;        /* Special entry for source */
  66. extern struct component OUT_node;        /* Special entry for output node marker */
  67.  
  68. extern int flag_FS    ;
  69. extern int flag_FE    ;
  70. extern int flag_FP    ;
  71. extern int flag_LOG    ;
  72. extern int flag_PHASE    ;
  73. extern int flag_RADS    ;
  74. extern int flag_PRTCMP  ;
  75. extern int flag_VS        ;
  76. extern int flag_OUT        ;
  77. extern int flag_MTRX       ;
  78.  
  79. extern double freq_start;        /* Starting Frequency */
  80. extern double freq_end;            /* Ending Frequency */
  81. extern int freq_part;            /* Number of partitions */
  82.  
  83. extern int N;                /* Number of distinct nodes */
  84. extern struct cmplx_num *addm_mtrx_ptr;    /* Addmittance Matrix Pointer */
  85. extern struct cmplx_num *curr_vect_ptr;    /* Impedance Matrix Pointer */
  86. /* -------------------------------------------------------------------*/
  87. /*************************************************************************
  88. *
  89. *        ENHANCED MATH PACKAGE
  90. *    An extra set of math functions, and in particular
  91. *    macros that perform complex number operations.  Note:
  92. *    all arguenments pass to the complex number macros MUST
  93. *    be pointers to a complex structure, and not an 
  94. *    actual structure.  Also the use of variable names "den"
  95. *    and "cmplx_tmp" will create problems.  These are macros,
  96. *    and not subroutines!
  97. *
  98. **************************************************************************/
  99.  
  100. struct cmplx_num    {        /* The Complex number structure */
  101.             double real;
  102.             double imag;
  103.             };
  104.  
  105. #define        abs(X)        ( X > 0 ? X : -X )
  106. #define        sqr(X)        ( X )*( X )
  107.  
  108. #define        cset(X, Y, Z)    X->real = Y; X->imag = Z
  109. #define        cmpl(X)        X->imag = - X->imag
  110. #define        cadd(X, Y, Z)    Z->real = X->real + Y->real;\
  111.                 Z->imag = X->imag + Y->imag
  112. #define        csub(X, Y, Z)    Z->real = X->real - Y->real;\
  113.                 Z->imag = X->imag - Y->imag
  114. #define        cmul(X, Y, Z)    { \
  115.                 struct cmplx_num cmplx_tmp; \
  116.                 cmplx_tmp.real = X->real * Y->real - X->imag * Y->imag ; \
  117.                 cmplx_tmp.imag = X->real * Y->imag + X->imag * Y->real ; \
  118.                 Z->real = cmplx_tmp.real; \
  119.                 Z->imag = cmplx_tmp.imag; \
  120.                 }
  121. #define        cdiv(X, Y, Z)    { \
  122.                 struct cmplx_num cmplx_tmp; \
  123.                 double den; \
  124.                 den = sqr( Y->real ) + sqr( Y->imag ); \
  125.                 cmplx_tmp.real = X->real * Y->real + X->imag * Y->imag ; \
  126.                 cmplx_tmp.imag = X->imag * Y->real - X->real * Y->imag ; \
  127.                 Z->imag = cmplx_tmp.imag / den; \
  128.                 Z->real = cmplx_tmp.real / den; \
  129.                 }
  130. #define        cmag(X)        sqrt( sqr( X->real) + sqr( X->imag))
  131. /* ----------------------------------------------------------------------------*/
  132. /************************************************************************
  133. *
  134. *    Main.c contains the primary loop of the spice simulation, namely
  135. *    the loop that increments frequency.  Main also first calls
  136. *    subroutines to read in all the net list, to check for critical
  137. *    information, and to possibly print the input data.  For each
  138. *    frequency, the addmitance matrix is recalculated, then 
  139. *    inverted.  The proper element is picked off, and printed.
  140. *
  141. *************************************************************************/
  142. #ifndef        ALLCODE
  143. #include    "header.h"
  144. #endif
  145.  
  146. main()
  147. {
  148. double freq_delta;    /* Frequency increment */
  149. double freq;        /* Current frequency of calculations */
  150. int i;
  151.  
  152. initialize();
  153. read_input();
  154. if( flag_PRTCMP == TRUE ) print_components();  /* Print all component values */
  155. check_input();
  156.  
  157.  
  158. /* Allocate Space for addmitance matrix */
  159. addm_mtrx_ptr = (struct cmplx_num *)malloc( N * N * sizeof( struct cmplx_num) );
  160. if( addm_mtrx_ptr <= NULL ) {
  161.     fprintf( stderr, "Insufficient Malloc memory allocations. Fatal error \n");
  162.     exit(1);
  163.     }
  164. curr_vect_ptr = (struct cmplx_num *)malloc( N * sizeof( struct cmplx_num) );
  165. if( curr_vect_ptr <= NULL ) {
  166.     fprintf( stderr, "Insufficient Malloc memory allocations. Fatal error \n");
  167.     exit(1);
  168.     }
  169.  
  170. #ifdef    IBM
  171. printf( "1 \n %d , SSPICE \n", freq_part );
  172. #endif
  173.  
  174. freq_delta = (freq_end - freq_start ) /freq_part;
  175. freq = freq_start ;
  176. for( i = 0; i <= freq_part ; i++ ) {
  177.     double omega;
  178.     fprintf( stderr, "*");
  179.  
  180.     if( flag_RADS == TRUE ) omega = freq ;
  181.     else omega = 2 * PIE * freq ;
  182.  
  183.     /* Fill addmittance matrix */
  184.     fill_matrix( omega);
  185.     if( flag_MTRX == TRUE ) print_mtrx( addm_mtrx_ptr );
  186.  
  187.     /* Take inverse of addmittance matrix */
  188.     inverse_matrix();
  189.     if( flag_MTRX == TRUE ) print_mtrx( addm_mtrx_ptr );
  190.  
  191.     print_output( freq );
  192.  
  193.     freq += freq_delta ;
  194.     }
  195.  
  196. } /* End of main */
  197. /* ----------------------------------------------------------------------- */
  198. /************************************************************************
  199. *
  200. *    This file contains two main subroutines: "fill_matrix" walks
  201. *    through all the component link lists, and makes numeric entries
  202. *    into the addmittance matrx.  This is only done for one frequency
  203. *    at a time.  "inverse_matrx" actually solves for the voltages
  204. *    of all nodes.  The solution ends up in the vector "curr", which
  205. *    originaly stands for the net current at each node.  The net current
  206. *    at each node is zero, with the exception of the source node.  It
  207. *    is assigned an external current of 1.0.
  208. *
  209. ***************************************************************************/
  210.  
  211. #ifndef        ALLCODE
  212. #include    "header.h"
  213. #endif
  214.  
  215. fill_matrix( omega)
  216. double omega;
  217. {
  218. struct component *x;
  219. struct cmplx_num addmitt;
  220. int i, j;
  221.  
  222. if( omega == 0.0 ) omega = 0.000001 ; /* To prevent division by zero */
  223.  
  224. /* Sets entire addmitance matrix to zero */
  225. for( i=0 ; i < N ; i++ ) {
  226.     for( j=0 ; j < N ; j++ ) {
  227.         addm( i, j)->real = 0.0;
  228.         addm( i, j)->imag = 0.0;
  229.         }
  230.     }
  231.  
  232. /* Enters all resistor addmittances to matrix */
  233. for( x = begin_res; x != NULL ; x = x->link ) {     /* Walks throught all res. */
  234.     cset( (&addmitt) , (1.0 / x->value), 0.0 );
  235.     add_element( &addmitt, x );
  236.     }
  237.  
  238. /* Enters all capacitor addmittances to matrix */
  239. for( x = begin_caps; x != NULL ; x = x->link ) {     /* Walks throught all cap. */
  240.     cset( (&addmitt) , 0.0 , (omega * x->value) );
  241.     add_element( &addmitt, x );
  242.     }
  243.  
  244. /* Enters all inductor addmittances to matrix */
  245. for( x = begin_inds; x != NULL ; x = x->link ) {     /* Walks throught all ind. */
  246.     cset( (&addmitt) , 0.0 , ( -1.0 / (omega * x->value)) );
  247.     add_element( &addmitt, x );
  248.     }
  249.  
  250. /* Enters all Voltage Controlled Current Sources addmittances to matrix */
  251. for( x = begin_vci; x != NULL ; x = x->link ) {
  252.     cset( (&addmitt) ,  x->value, 0.0 );
  253.  
  254.     if( x->node_k >= 0 ) {  /* Not a ground node */
  255.         if( x->node_i >= 0 ) {
  256.             addm( x->node_k, x->node_i)->real -= addmitt.real;
  257.             addm( x->node_k, x->node_i)->imag -= addmitt.imag;
  258.             }
  259.         if( x->node_j >= 0 ) {
  260.             addm( x->node_k, x->node_j)->real += addmitt.real;
  261.             addm( x->node_k, x->node_j)->imag += addmitt.imag;
  262.             }
  263.         }
  264.  
  265.     if( x->node_l >= 0 ) {  /* Not a ground node */
  266.         if( x->node_i >= 0 ) {
  267.             addm( x->node_l, x->node_i)->real += addmitt.real;
  268.             addm( x->node_l, x->node_i)->imag += addmitt.imag;
  269.             }
  270.         if( x->node_j >= 0 ) {
  271.             addm( x->node_l, x->node_j)->real -= addmitt.real;
  272.             addm( x->node_l, x->node_j)->imag -= addmitt.imag;
  273.             }
  274.         }
  275.     }
  276.  
  277. /* Enters all Voltage Controlled Voltage Sources addmittances to matrix */
  278. for( x = begin_vcv; x != NULL ; x = x->link ) {
  279.     cset( (&addmitt) ,  x->value, 0.0 );
  280.     for( i = 0 ; i < N ; i ++) {
  281.         addm( x->node_k, i)->real = 0.0 ;
  282.         addm( x->node_k, i)->imag = 0.0 ;
  283.         }
  284.     cset( addm( x->node_k, x->node_k), 1.0, 0.0 );
  285.  
  286.     if( x->node_j >= 0 ) {  /* Not a ground node */
  287.         addm( x->node_k, x->node_j)->real = x->value ;
  288.         }
  289.     if( x->node_i >= 0 ) {  /* Not a ground node */
  290.         addm( x->node_k, x->node_i)->real = -( x->value) ;
  291.         }
  292.     }
  293.  
  294. /* Enters source addmittance to matrix */
  295. for( i = 0 ; i < N ; i ++) {
  296.     addm( (&VS_source)->node_i, i)->real = 0.0 ;
  297.     addm( (&VS_source)->node_i, i)->imag = 0.0 ;
  298.     }
  299. cset( addm( (&VS_source)->node_i, (&VS_source)->node_i), 1.0, 0.0 );
  300.  
  301.  
  302. /* Initialize current vertor */
  303. for( i = 0; i < N; i++ ) { cset( curr( i ), 0.0 , 0.0 ); }
  304. cset( curr( (&VS_source)->node_i), 1.0, 0.0 );
  305. } /* End of enter_matrix */
  306.  
  307. add_element( a, x)
  308. struct cmplx_num *a;
  309. struct component *x;
  310. {
  311.  
  312. if( (x->node_i >= 0) && (x->node_j >=0) ) {
  313.     addm( x->node_i, x->node_i)->real += a->real;
  314.     addm( x->node_i, x->node_i)->imag += a->imag;
  315.  
  316.     addm( x->node_j, x->node_j)->real += a->real;
  317.     addm( x->node_j, x->node_j)->imag += a->imag;
  318.  
  319.     addm( x->node_i, x->node_j)->real -= a->real;
  320.     addm( x->node_i, x->node_j)->imag -= a->imag;
  321.  
  322.     addm( x->node_j, x->node_i)->real -= a->real;
  323.     addm( x->node_j, x->node_i)->imag -= a->imag;
  324.     }
  325.  
  326. else if( (x->node_i >= 0) && (x->node_j < 0) ) {
  327.     addm( x->node_i, x->node_i)->real += a->real;
  328.     addm( x->node_i, x->node_i)->imag += a->imag;
  329.     }
  330.  
  331. else if( (x->node_i < 0) && (x->node_j >= 0) ) {
  332.     addm( x->node_j, x->node_j)->real += a->real;
  333.     addm( x->node_j, x->node_j)->imag += a->imag;
  334.     }
  335.  
  336. else {
  337.     fprintf( stderr, "Component %s Seems to have no conections to any active nodes.  It is ignored \n", x->name );
  338.     }
  339.  
  340.  
  341. }
  342.  
  343.  
  344. #define        csubeq(X, Y)    { \
  345.                 X->real -= Y->real; \
  346.                 X->imag -= Y->imag; \
  347.                 }
  348.  
  349. /********************************************************************
  350. *  The following routing is an enhanced inversion gaussain elimination
  351. *  inversion routing. Is seeks the largest value in any column, and
  352. *  then performs pivoting and elimination using this entry.  In order to
  353. *  keep track of which columns have already been pivoted, the algorithm
  354. *  checks the sum of the magnitude of the elements along the row to
  355. *  just before the test pivot element.  If this sum is zero, or very
  356. *  close to it, it is assumed that this row has already been pivoited.
  357. **********************************************************************/
  358. inverse_matrix()
  359. {
  360. struct cmplx_num one_x, tmp_x, scale_x;
  361. struct cmplx_num *one, *tmp, *scale;
  362. int i, j, k, ii;
  363. double pre_max;        /* Previous Maximum */
  364. double cur_max;        /* Current Maximum */
  365. double par_row_sum ;     /* Partial Row Sum.  To detect if this row has already been cleared */
  366.  
  367. /* Need Temporary space that are pointers */
  368. one = &one_x;
  369. tmp = &tmp_x;
  370. scale = &scale_x;
  371.  
  372. cset( one, 1.0, 0.0);
  373.  
  374.  
  375. /* Scan across all columns */
  376. for( j = 0; j < N; j++) { 
  377.  
  378.     /* Find Maximun entry in this row. The magnitude of all previous entries in
  379.        this row is sumed in "par_row_sum".  Only if this sum is zero is it 
  380.        valid to pivot on this row. */
  381.     pre_max = 0.0 ;
  382.     for( i = 0; i < N; i++ ) {
  383.         par_row_sum = 0.0;
  384.         for( k = 0; k < j ; k++ ) par_row_sum += cmag( addm( i, k) ) ;
  385.         cur_max = cmag( addm( i, j) );
  386.  
  387.         /* test this pivot point */
  388.         if( ( cur_max > pre_max ) && ( cur_max > (1000000000 * par_row_sum ) )) {
  389.             ii = i;        /* Mark this N as ii */
  390.             pre_max = cur_max ;
  391.             }
  392.         }
  393.  
  394.     /* Normalize the row chosen above */
  395.     if( cmag( addm( ii, j )) == 0.0 ) {
  396.         fprintf( stderr,"##ERROR##\n, division by zero.\n");
  397.         exit(0);
  398.         }
  399.     cdiv( one , addm( ii, j), scale );
  400.     for( k = j; k < N; k++ ) {
  401.         cmul( addm( ii, k), scale, addm( ii, k) );
  402.         }
  403.     cmul( curr( ii), scale, curr( ii) );
  404.  
  405.  
  406.     /* Null out all other entries in this column, except
  407.        choses entry to zero. */
  408.     for( i = 0; i < N; i++ ) {
  409.         /* Don't do if pivot row, or if semi-pivot element is zero */
  410.         if( ( i != ii ) && ( cmag( addm( i , j ) ) != 0.0) ) {
  411.             cmul( addm( i, j), curr( ii ), tmp );
  412.             csubeq( curr( i ), tmp );
  413.             for( k = N-1; k >= j ; k-- ) {
  414.                 cmul( addm( i, j), addm( ii, k), tmp );
  415.                 csubeq( addm( i, k), tmp );
  416.                 }
  417.             }
  418.         }
  419.     } /* End of Column Scan */
  420. }  /* End of matrix inversion */
  421. /*****************************************************************************
  422. *
  423. *    This file contains subroutines for reading in the net list, processing
  424. *    the node names into node numbers, setting true appropriate flags, and
  425. *    recording misc. infomation, such as starting frequency, etc.
  426. *
  427. *****************************************************************************/
  428.  
  429. #ifndef        ALLCODE
  430. #include     "header.h"
  431. #endif
  432. #define        fscanf        my_fscanf
  433.  
  434. /* Primary input read loop.  Exits when EOF is detected */
  435. read_input()
  436. {
  437. char str[ STR_LNGTH ] ;    /* String buffer for input reading */
  438. FILE *fp_in    = stdin;
  439.  
  440. /* Loops over all input lines, and enters data into proper variables */
  441. while( fscanf( fp_in, "%s", str) != EOF ) {
  442.     double scan_modifier();
  443.     struct component *add_cmpt_link();
  444.  
  445.     if( !strcmp( str, "FS") ) {
  446.         fscanf( fp_in, "%f", &freq_start);
  447.         freq_start *= (int)scan_modifier( fp_in );
  448.         flag_FS = TRUE;
  449.         }
  450.  
  451.     else if( !strcmp( str, "FE") ) {
  452.         fscanf( fp_in, "%f", &freq_end);
  453.         freq_end *= (int)scan_modifier( fp_in );
  454.         flag_FE = TRUE;
  455.         }
  456.  
  457.     else if( !strcmp( str, "FP") ) {
  458.         fscanf( fp_in, "%d", &freq_part);
  459.         freq_part *= (int)scan_modifier( fp_in );
  460.         flag_FP = TRUE;
  461.         }
  462.  
  463.     else if( !strcmp( str, "PHASE") ) {
  464.         flag_PHASE = TRUE;
  465.         scan_modifier( fp_in );
  466.         }
  467.  
  468.     else if( !strcmp( str, "LOG") ) {
  469.         flag_LOG = TRUE;
  470.         scan_modifier( fp_in );
  471.         }
  472.  
  473.     else if( !strcmp( str, "MTRX") ) {
  474.         flag_MTRX = TRUE;
  475.         scan_modifier( fp_in );
  476.         }
  477.  
  478.     else if( !strcmp( str, "RADS") ) {
  479.         flag_RADS = TRUE;
  480.         scan_modifier( fp_in );
  481.         }
  482.  
  483.     else if( !strcmp( str, "PRTCMP") ) {
  484.         flag_PRTCMP = TRUE;
  485.         scan_modifier( fp_in );
  486.         }
  487.  
  488.     else if( !strcmp( str, "OUT") ) {
  489.         flag_OUT = TRUE;
  490.         fscanf( fp_in, "%s", OUT_node.name_i );
  491.         OUT_node.node_i = node_numb( OUT_node.name_i);
  492.         scan_modifier( fp_in );
  493.         }
  494.  
  495.     else if( *str == 'L' ) {
  496.         add_cmpt_link( &begin_inds, &end_inds );
  497.         strcpy( end_inds->name, str);
  498.         fscanf( fp_in, "%s", end_inds->name_i );
  499.         end_inds->node_i = node_numb( end_inds->name_i);
  500.         fscanf( fp_in, "%s", end_inds->name_j );
  501.         end_inds->node_j = node_numb( end_inds->name_j);
  502.         fscanf( fp_in, "%f", &end_inds->value );
  503.         end_inds->value *= scan_modifier( fp_in );
  504.         }
  505.  
  506.     else if( *str == 'C' ) {
  507.         add_cmpt_link( &begin_caps, &end_caps );
  508.         strcpy( end_caps->name, str);
  509.         fscanf( fp_in, "%s", end_caps->name_i );
  510.         end_caps->node_i = node_numb( end_caps->name_i);
  511.         fscanf( fp_in, "%s", end_caps->name_j );
  512.         end_caps->node_j = node_numb( end_caps->name_j);
  513.         fscanf( fp_in, "%f", &end_caps->value );
  514.         end_caps->value *= scan_modifier( fp_in );
  515.         }
  516.  
  517.     else if( *str == 'R' ) {
  518.         add_cmpt_link( &begin_res, &end_res );
  519.         strcpy( end_res->name, str);
  520.         fscanf( fp_in, "%s", end_res->name_i );
  521.         end_res->node_i = node_numb( end_res->name_i);
  522.         fscanf( fp_in, "%s", end_res->name_j );
  523.         end_res->node_j = node_numb( end_res->name_j);
  524.         fscanf( fp_in, "%f", &end_res->value );
  525.         end_res->value *= scan_modifier( fp_in );
  526.         }
  527.  
  528.     else if( !strcmp( str, "VCI") ) {
  529.         add_cmpt_link( &begin_vci, &end_vci );
  530.         strcpy( end_vci->name, str);
  531.         fscanf( fp_in, "%s", end_vci->name_i );
  532.         end_vci->node_i = node_numb( end_vci->name_i);
  533.         fscanf( fp_in, "%s", end_vci->name_j );
  534.         end_vci->node_j = node_numb( end_vci->name_j);
  535.         fscanf( fp_in, "%s", end_vci->name_k );
  536.         end_vci->node_k = node_numb( end_vci->name_k);
  537.         fscanf( fp_in, "%s", end_vci->name_l );
  538.         end_vci->node_l = node_numb( end_vci->name_l);
  539.         fscanf( fp_in, "%f", &end_vci->value );
  540.         end_vci->value *= scan_modifier( fp_in );
  541.         }
  542.  
  543.     else if( !strcmp( str, "VCV") ) {
  544.         add_cmpt_link( &begin_vcv, &end_vcv );
  545.         strcpy( end_vcv->name, str);
  546.         fscanf( fp_in, "%s", end_vcv->name_i );
  547.         end_vcv->node_i = node_numb( end_vcv->name_i);
  548.         fscanf( fp_in, "%s", end_vcv->name_j );
  549.         end_vcv->node_j = node_numb( end_vcv->name_j);
  550.         fscanf( fp_in, "%s", end_vcv->name_k );
  551.         end_vcv->node_k = node_numb( end_vcv->name_k);
  552.         fscanf( fp_in, "%f", &end_vcv->value );
  553.         end_vcv->value *= scan_modifier( fp_in );
  554.         }
  555.  
  556.     else if( !strcmp( str, "VS") ) {
  557.         flag_VS = TRUE;
  558.         fscanf( fp_in, "%s", VS_source.name_i );
  559.         VS_source.node_i = node_numb( VS_source.name_i);
  560.         scan_modifier( fp_in );
  561.         }
  562.  
  563.  
  564.     else {
  565.         fprintf( stderr, "Syntax Error: Unrecognized input string \"%s\". String ignored. \n", str );
  566.         }
  567.     } /* End of while */
  568. } /* End of read routine */
  569.  
  570.  
  571.  
  572.         
  573.  
  574. /*  Assigns a number to any ascii node name.  Number are chosen in
  575. *   sequential order beginning with zero.  If the name passed to
  576. *   this routing already exits, then the node number assigned
  577. *   to this name previously is returned. Gnd and several variations
  578. *   are special reserved names which always have the value -1
  579. *   associated with them.
  580. */
  581. node_numb( str )
  582. char *str;
  583. {
  584. struct name_list {         /* Link List structure to keep node names */
  585.     char name[ STR_LNGTH ];
  586.     int node_number;
  587.     struct name_list *link;
  588.     };
  589.  
  590. static struct name_list *bgn_nm_lst = NULL;
  591. static struct name_list *end_nm_lst = NULL;
  592. struct name_list *new_link, *x;
  593.  
  594. if( !strcmp( str, "GND") || !strcmp( str, "gnd") || !strcmp( str, "Gnd") ) return( -1 );
  595. for( x = bgn_nm_lst; x != NULL; x = x->link )
  596.     {
  597.     if( !strcmp( x->name, str)) return( x->node_number );
  598.     }
  599.  
  600. new_link = (struct name_list *)malloc( sizeof (struct name_list));
  601.  
  602. new_link->link = NULL;
  603. if( bgn_nm_lst == NULL ) bgn_nm_lst = new_link; /* IF first link, start chain */
  604. else end_nm_lst->link = new_link;         /* link head */
  605. end_nm_lst = new_link;                 /* Update end pointer */
  606.  
  607. strcpy( new_link->name, str);
  608. new_link->node_number = N++;
  609. return( new_link->node_number );
  610. }
  611.  
  612.  
  613. #undef        fscanf
  614. /* A lousy fucking patch for micro-soft C.  It doesn't
  615. *  seem to read floating point numbers properly.      
  616. */
  617. my_fscanf( fp, str, reg)
  618. FILE *fp;
  619. char *str;
  620. union {
  621.     int x;
  622.     double y;
  623.     char t;
  624.     } *reg;
  625.     
  626. {
  627. double atof();
  628. if( !strcmp( str, "%f") ) {
  629.     char temp[40];
  630.     fscanf( fp, "%s", temp);
  631.     reg->y = atof( temp);
  632.     }
  633. else fscanf( fp, str, reg);
  634. }
  635. #define        fscanf        my_fscanf
  636. /* ------------------------------------------------------------------------ */
  637. /********************************************************************
  638. *
  639. *    This file contains all routines which print anything to
  640. *    the user.
  641. *
  642. **********************************************************************/
  643.  
  644. #ifndef        ALLCODE
  645. #include     "header.h"
  646. #endif
  647.  
  648. #define     FLG_STR(X)    if( X == TRUE ) strcpy( str, "TRUE" ); \
  649.                         else strcpy( str, "FALSE")
  650.  
  651. /* Prints components and all other important information */
  652. print_components()
  653. {
  654. FILE *fp_out    = stdout;
  655. char str[ STR_LNGTH ] ;    /* String buffer for input reading */
  656. struct component *x;    /* a pointer to current component structure */
  657.  
  658. FLG_STR(flag_FS);
  659. fprintf( fp_out,"+Starting frequency flag is set %s \n", str);
  660. FLG_STR(flag_FE);
  661. fprintf( fp_out,"+Ending frequency flag is set %s \n", str);
  662. FLG_STR(flag_FP);
  663. fprintf( fp_out,"+Number of Points  flag is set %s \n", str);
  664. FLG_STR(flag_RADS);
  665. fprintf( fp_out,"+Radian flag is set %s \n", str);
  666. FLG_STR(flag_LOG);
  667. fprintf( fp_out,"+Logarithmitic flag is set %s \n", str);
  668. FLG_STR(flag_VS);
  669. fprintf( fp_out,"+Source designation flag is set %s \n", str);
  670. FLG_STR(flag_OUT);
  671. fprintf( fp_out,"+Output designation flag is set %s \n", str);
  672. fprintf( fp_out,"\n");
  673.  
  674. if( flag_FS == TRUE ) fprintf( fp_out,"+Starting Frequency = %f \n", freq_start);
  675. if( flag_FE == TRUE ) fprintf( fp_out,"+Ending Frequency = %f \n", freq_end);
  676. if( flag_FP == TRUE ) fprintf( fp_out,"+Number of partitions is = %d \n", freq_part);
  677.  
  678. fprintf( fp_out, "\n   Component   I Node   J Node   K Node   L Node     Value  \n");
  679.  
  680. for( x = begin_res; x != NULL ; x = x->link )
  681.     fprintf( fp_out, "+ %10s %8s %8s                       %10e Ohms \n", x->name, x->name_i, x->name_j, x->value);
  682.  
  683. for( x = begin_inds; x != NULL ; x = x->link ) 
  684.     fprintf( fp_out, "+ %10s %8s %8s                       %10e Heneries \n", x->name, x->name_i, x->name_j, x->value);
  685.  
  686. for( x = begin_caps; x != NULL ; x = x->link )
  687.     fprintf( fp_out, "+ %10s %8s %8s                       %10e Farads \n", x->name, x->name_i, x->name_j, x->value);
  688.  
  689. for( x = begin_vcv; x != NULL ; x = x->link )
  690.     fprintf( fp_out, "+ %10s %8s %8s %8s              %10e Av \n", x->name, x->name_i, x->name_j, x->name_k , x->value);
  691.  
  692. for( x = begin_vci; x != NULL ; x = x->link )
  693.     fprintf( fp_out, "+ %10s %8s %8s %8s %8s     %10e Mohs \n", x->name, x->name_i, x->name_j, x->name_k , x->name_l , x->value);
  694.  
  695. fprintf( fp_out, "+         VS %8s \n", VS_source.name_i);
  696. fprintf( fp_out, "+        OUT %8s \n", OUT_node.name_i);
  697.  
  698. } /* End of print routine */
  699.  
  700.  
  701.  
  702.  
  703. /*  Prints out all elements of a complex matrix 
  704. *   Primarily to be used for checking and debugging 
  705. */
  706. print_mtrx( mtrx )
  707. struct cmplx_num *mtrx;
  708. {
  709. FILE *fp_out    = stdout ;
  710. int i, j;
  711.  
  712. fprintf( fp_out,"------------------------------------------------------------ \n");
  713. for( i=0; i<N ; i++ ) {
  714.     for( j=0; j<N ; j++ ) fprintf( fp_out, "%6.2f  ", (mtrx + i + N*j)->real);
  715.     fprintf( fp_out, "   %6.2f\n", curr( i )->real);
  716.     for( j=0; j<N ; j++ ) fprintf( fp_out, "%6.2f  ", (mtrx + i + N*j)->imag);
  717.     fprintf( fp_out, "   %6.2f\n", curr( i )->imag);
  718.     fprintf( fp_out, "\n");
  719.     }
  720. fprintf( fp_out,"------------------------------------------------------------ \n");
  721. }
  722.         
  723.  
  724.  
  725. /* Prints the output transfer function */
  726. print_output( omega)
  727. double omega;
  728. {
  729. double magnitude, phase;
  730. int i, ii, out;
  731. double pre_max;
  732.  
  733. out = OUT_node.node_i ;
  734. pre_max = 0.0 ;
  735. for( i =0 ; i < N ; i++ ) {
  736.     if( cmag( addm( i, out )) > pre_max ) {
  737.         ii = i;
  738.         pre_max = cmag( addm( i, out ));
  739.         }
  740.     }
  741.     
  742.  
  743. phase = atan2( curr( ii )->imag, curr( ii )->real ); 
  744. magnitude = cmag( curr( ii ) );
  745. if( magnitude == 0.0 ) magnitude = 1000000000.0;
  746.  
  747. #ifdef    IBM
  748. if( flag_LOG == TRUE ) printf( "%f , %f \n", omega, ( -20 * log10( magnitude ) ) );
  749. else if( flag_PHASE == TRUE ) printf( "%f , %f \n", omega, ( 180 * phase / PIE) ); 
  750. else printf( "%f , %f \n", omega, magnitude );
  751.  
  752. #else
  753. if( flag_LOG == TRUE ) printf( "%f %f \n", omega, ( -20 * log10( magnitude ) ) );
  754. else if( flag_PHASE == TRUE ) printf( "%f %f \n", omega, ( 180 * phase / PIE ) ); 
  755. else printf( "%f %f \n", omega, magnitude );
  756. #endif
  757. }
  758. /* ----------------------------------------------------------------------------- */
  759. /****************************************************************************
  760. *    
  761. *    General miscellaneous subroutines that I didn't want to
  762. *    be cluttering up the more improtant files.
  763. *
  764. *****************************************************************************/
  765.  
  766. #ifndef        ALLCODE
  767. #include    "header.h"
  768. #endif
  769.  
  770.  
  771. /* Looks for some type of value modifie, such as kilo (x1000)
  772. *  or micro (x0.000001).  Then will be expecting a colon as end
  773. *  of line marker.  Returns a double persion number for
  774. *  scaling purposes.    
  775. */
  776.  
  777. double scan_modifier( input )
  778. FILE *input;
  779. {
  780. double scale = 1.0;
  781. char string[ STR_LNGTH ];
  782.  
  783. while( (fscanf( input, "%s", string) != EOF) && (*string != ';' ) ) {
  784.  
  785.     switch( *string) {
  786.  
  787.     case 'K':
  788.         scale = 1000;
  789.         break;
  790.  
  791.     case 'M':
  792.         scale = 1000000;
  793.         break;
  794.  
  795.     case 'G':
  796.         scale = 1000000000;
  797.         break;
  798.  
  799.     case 'm':
  800.         scale = 0.001;
  801.         break;
  802.  
  803.     case 'u':
  804.         scale = 0.000001;
  805.         break;
  806.  
  807.     case 'n':
  808.         scale = 0.000000001;
  809.         break;
  810.  
  811.     case 'p':
  812.         scale = 0.000000000001;
  813.         break;
  814.  
  815.     case 'E':
  816.         {
  817.         int K;
  818.         fscanf( input, "%d", &K);
  819.         if( K > 30 ) K = 30;
  820.         if( K < -30 ) K = -30;
  821.         while( K-- > 0 ) scale *= 10;
  822.         while( K++ < 0 ) scale *= 0.1;
  823.         }
  824.         break;
  825.  
  826.     default:
  827.         fprintf( stderr, "Unrecognized string \"%s\" \n", string);
  828.         fprintf( stderr, "Metric modifying unit expected. \n" );
  829.         fprintf( stderr, "String Ignored. \n");
  830.         
  831.  
  832.         }
  833.     }
  834. return( scale );
  835. }
  836.  
  837.  
  838. /* Add a Component structor to the Link.  A component structor is
  839. *  added to the link list.  The end_ptr_add arguement is the ADDRESS of a
  840. *  pointer which points to the last structure in the list.  In the
  841. *  case of a zero list, "pointer_address" is the ADDRESS of the
  842. *  variable "begin_list"
  843. */
  844.  
  845. struct component *add_cmpt_link( bgn_ptr_add, end_ptr_add )
  846. struct component **end_ptr_add;
  847. struct component **bgn_ptr_add;
  848. {
  849. struct component *new_link;
  850.  
  851. if( ( new_link = (struct component *)malloc( sizeof ( struct component ) ) ) == NULL )  {
  852.     fprintf( stderr, "Operating system fault Insufficient memory resources. \n"); 
  853.     fprintf( stderr, "Fatal Error. \n");
  854.     exit(1);
  855.     }
  856.  
  857. new_link->link = NULL;                    /* Link tail */
  858. if( *bgn_ptr_add == NULL ) *bgn_ptr_add = new_link;    /* IF first, Start Chain */
  859. else (*end_ptr_add)->link = new_link;            /* link head */
  860. *end_ptr_add = new_link;                /* Update end pointer */ 
  861. return( new_link );
  862. }
  863.  
  864.  
  865.  
  866. /*
  867. *  Initializes global variables only.
  868. *  This cannot be nicely done in header.h due to
  869. *  its dual nature use with a makefile
  870. */
  871.  
  872. initialize()
  873. {
  874. begin_caps = NULL;    /* Start of capacitor link list */
  875. begin_inds = NULL;    /* Start of inductor link list */
  876. begin_res  = NULL;    /* Start of resistor link list */
  877. begin_vcv  = NULL;    /* Start of VCV link list */
  878. begin_vci  = NULL;    /* Start of VCI link list */
  879. end_caps   = NULL;    /* Last entry of capacitor link list */
  880. end_inds   = NULL;    /* Last entry of inductor link list */
  881. end_res    = NULL;    /* Last entry of resistor link list */
  882. end_vcv    = NULL;    /* Last entry of VCV link list */
  883. end_vci    = NULL;    /* Last entry of VCI link list */
  884.  
  885. flag_FS        =    FALSE;
  886. flag_FE        =    FALSE;
  887. flag_FP        =    FALSE;
  888. flag_LOG    =    FALSE;
  889. flag_PHASE    =    FALSE;
  890. flag_RADS    =    FALSE;
  891. flag_PRTCMP     =    FALSE;
  892. flag_VS     =    FALSE;
  893. flag_OUT     =    FALSE;
  894. flag_MTRX     =    FALSE;
  895.  
  896. N        =     0;
  897. }
  898.  
  899.  
  900. /* Checks if essential input informaiton is present */
  901. check_input()
  902. {
  903.  
  904. if( flag_VS == FALSE ) {
  905.     fprintf( stderr, "No source node has been specified! \n" );
  906.     fprintf( stderr, "Fatal Error. \n");
  907.     exit(1);
  908.     }
  909.  
  910. if( flag_OUT == FALSE ) {
  911.     fprintf( stderr, "No output node has been specified! \n" );
  912.     fprintf( stderr, "Fatal Error. \n");
  913.     exit(1);
  914.     }
  915.  
  916. if( flag_FS == FALSE ) {
  917.     fprintf( stderr, "Starting frequency unspecified! \n" );
  918.     fprintf( stderr, "Fatal Error. \n");
  919.     exit(1);
  920.     }
  921.  
  922. if( flag_FE == FALSE ) {
  923.     fprintf( stderr, "Ending frequency unspecified! \n" );
  924.     fprintf( stderr, "Fatal Error. \n");
  925.     exit(1);
  926.     }
  927.  
  928. if( flag_FP == FALSE ) {
  929.     fprintf( stderr, "Number of frequency partitions unspecified! \n" );
  930.     fprintf( stderr, "Default of 20 will be used. \n");
  931.     freq_part = 20;
  932.     }
  933. }
  934. /* ------------------------------------------------------------------------- */
  935.