home *** CD-ROM | disk | FTP | other *** search
/ NetNews Usenet Archive 1993 #3 / NN_1993_3.iso / spool / alt / sources / 3112 < prev    next >
Encoding:
Text File  |  1993-01-28  |  60.6 KB  |  1,983 lines
  1. Path: sparky!uunet!ferkel.ucsb.edu!taco!rock!stanford.edu!ames!haven.umd.edu!darwin.sura.net!bogus.sura.net!ukma!cs.widener.edu!dsinc!ub!galileo.cc.rochester.edu!ee.rochester.edu!rbc!al
  2. From: al@rbc.uucp (Al Davis)
  3. Newsgroups: alt.sources
  4. Subject: ACS  circuit simulator  part 10/20
  5. Message-ID: <1993Jan27.040639.11348@rbc.uucp>
  6. Date: 27 Jan 93 04:06:39 GMT
  7. Sender: al@rbc.uucp (Al Davis)
  8. Organization: Huh?
  9. Lines: 1972
  10.  
  11. #! /bin/sh
  12. # This is a shell archive, meaning:
  13. # 1. Remove everything above the #! /bin/sh line.
  14. # 2. Save the resulting text in a file.
  15. # 3. Execute the file with /bin/sh (not csh) to create the files:
  16. #    src/ac.c
  17. #    src/ac_clear.c
  18. #    src/ac_fill.c
  19. #    src/ac_fix.c
  20. #    src/ac_load.c
  21. #    src/ac_out.c
  22. #    src/ac_probe.c
  23. #    src/ac_setup.c
  24. #    src/ac_sweep.c
  25. #    src/ac_z.c
  26. #    src/allocate.c
  27. #    src/argparse.c
  28. #    src/array.c
  29. # This archive created: Tue Jan 26 22:50:59 1993
  30. export PATH; PATH=/bin:$PATH
  31. if test -f 'src/ac.c'
  32. then
  33.     echo shar: will not over-write existing file "'src/ac.c'"
  34. else
  35. cat << \SHAR_EOF > 'src/ac.c'
  36. /* ac  12/29/92
  37.  * Copyright 1983-1992   Albert Davis
  38.  * ac analysis top
  39.  */
  40. #include "ecah.h"
  41. #include "ac.h"
  42. #include "mode.h"
  43. #include "status.h"
  44. #include "declare.h"
  45. /*--------------------------------------------------------------------------*/
  46.     void    cmd_ac(const char*,int*);
  47. /*--------------------------------------------------------------------------*/
  48. extern struct status stats;
  49. extern ac_t ac;
  50. extern int run_mode;    /* variations on handling of dot cmds            */
  51. extern int sim_mode;    /* simulation type (AC, DC, ...)            */
  52. /*--------------------------------------------------------------------------*/
  53. void cmd_ac(cmd,cnt)
  54. const char *cmd;
  55. int  *cnt;
  56. {
  57.  sim_mode = sAC;
  58.  time_zstart(&(stats.total));
  59.  time_zstart(&(stats.ac));
  60.  
  61.  allocate(sim_mode);
  62.  ac_setup(cmd,cnt,&ac);
  63.  if (run_mode != rEXECUTE)
  64.     return;
  65.  ac_sweep(&ac);
  66.  
  67.  time_stop(&(stats.ac));
  68.  time_stop(&(stats.total));
  69. }
  70. /*--------------------------------------------------------------------------*/
  71. /*--------------------------------------------------------------------------*/
  72. SHAR_EOF
  73. fi # end of overwriting check
  74. if test -f 'src/ac_clear.c'
  75. then
  76.     echo shar: will not over-write existing file "'src/ac_clear.c'"
  77. else
  78. cat << \SHAR_EOF > 'src/ac_clear.c'
  79. /* acclear  04/02/92
  80.  * Copyright 1983-1992   Albert Davis
  81.     Clears working arrays allocated by acalloc.     It is
  82. necessary to clear them separately because they must be cleared for each new
  83. frequency, but not reallocated.     We cannot use free/calloc because it would
  84. require setting new pointer tables.
  85.     */
  86. #include "ecah.h"
  87. #include "status.h"
  88. #include "declare.h"
  89. /*--------------------------------------------------------------------------*/
  90.     void    ac_clear(void);
  91. /*--------------------------------------------------------------------------*/
  92. extern double *recon, *imcon;        /* constant terms                */
  93. extern double *reals, *imags;        /* big array allocation bases        */
  94. extern int decomp;
  95.  
  96. extern const unsigned aspace;        /* count of non-zero array elements        */
  97. extern const struct status stats;
  98. /*--------------------------------------------------------------------------*/
  99. void ac_clear()
  100. {
  101.  (void)memset((void*)reals, 0, (size_t)aspace      * sizeof(double));
  102.  (void)memset((void*)imags, 0, (size_t)aspace      * sizeof(double));
  103.  (void)memset((void*)recon, 0, (size_t)(stats.total_nodes+2)* sizeof(double));
  104.  (void)memset((void*)imcon, 0, (size_t)(stats.total_nodes+2)* sizeof(double));
  105.  decomp = 0;
  106. }
  107. /*--------------------------------------------------------------------------*/
  108. /*--------------------------------------------------------------------------*/
  109. SHAR_EOF
  110. fi # end of overwriting check
  111. if test -f 'src/ac_fill.c'
  112. then
  113.     echo shar: will not over-write existing file "'src/ac_fill.c'"
  114. else
  115. cat << \SHAR_EOF > 'src/ac_fill.c'
  116. /* ac_fill  12/29/92
  117.  * Copyright 1983-1992   Albert Davis
  118.  * Fill working matrix for AC analysis
  119.  */
  120. #include "ecah.h"
  121. #include "branch.h"
  122. #include "dev.h"
  123. #include "status.h"
  124. #include "declare.h"
  125. /*--------------------------------------------------------------------------*/
  126.     void    ac_solve(void);
  127.     void    ac_fill_rl(branch_t*);
  128. /*--------------------------------------------------------------------------*/
  129. extern struct status stats;
  130. /*--------------------------------------------------------------------------*/
  131. void ac_solve()
  132. {
  133.  ac_clear();
  134.  time_start(&(stats.load));
  135.  ac_fill_rl(firstbranch_dev());
  136.  time_stop(&(stats.load));
  137.  xsolve();
  138. }
  139. /*--------------------------------------------------------------------------*/
  140. void ac_fill_rl(brh)
  141. branch_t *brh;
  142. {
  143.  branch_t *stop;
  144.  if (exists(brh)){
  145.     stop = brh;
  146.     do {
  147.        doac_branch(brh);
  148.     } while (brh=nextbranch_dev(brh),  brh != stop);
  149.  }
  150. }
  151. /*--------------------------------------------------------------------------*/
  152. /*--------------------------------------------------------------------------*/
  153. SHAR_EOF
  154. fi # end of overwriting check
  155. if test -f 'src/ac_fix.c'
  156. then
  157.     echo shar: will not over-write existing file "'src/ac_fix.c'"
  158. else
  159. cat << \SHAR_EOF > 'src/ac_fix.c'
  160. /* acfix  01/18/93
  161.  * Copyright 1983-1992   Albert Davis
  162.  * Takes care of nonlinearities, behavioral modeling, etc. in ac analysis.
  163.  */
  164. #include "ecah.h"
  165. #include "ac.h"
  166. #include "branch.h"
  167. #include "error.h"
  168. #include "expr.h"
  169. #include "declare.h"
  170. /*--------------------------------------------------------------------------*/
  171.     complex_t    acfix(branch_t*);
  172. static    void    acf_ac(const branch_t*,double**,int*,complex_t*);
  173. static    void    acf_dc(const branch_t*,double**,int*,complex_t*);
  174. static    void    acf_dctran(const branch_t*,double**,int*,complex_t*);
  175. static    void    acf_frequency(const branch_t*,double**,int*,complex_t*);
  176. static    void    acf_period(const branch_t*,double**,int*,complex_t*);
  177. static    void    acf_ramp(const branch_t*,double**,int*,complex_t*);    
  178. static    void    acf_time(const branch_t*,double**,int*,complex_t*);    
  179. static    void    acf_tran(const branch_t*,double**,int*,complex_t*);    
  180. static    void    acf_bandwidth(const branch_t*,double**,complex_t*);
  181. static    void    acf_complex(const branch_t*,double**,complex_t*);
  182. static    void    acf_cornerdown(const branch_t*,double**,complex_t*);
  183. static    void    acf_cornerup(const branch_t*,double**,complex_t*);
  184. static    void    acf_delay(const branch_t*,double**,complex_t*);    
  185. static    void    acf_exp(const branch_t*,double**,complex_t*);    
  186. static    void    acf_expterm(const branch_t*,double**,complex_t*);
  187. static    void    acf_generator(const branch_t*,double**,complex_t*);
  188. static    void    acf_max(const branch_t*,double**,complex_t*);    
  189. static    void    acf_netfunc(const branch_t*,double**,complex_t*);    
  190. static    void    acf_notch(const branch_t*,double**,complex_t*);    
  191. static    void    acf_numeric(const branch_t*,double**,complex_t*);    
  192. static    void    acf_offset(const branch_t*,double**,complex_t*);
  193. static    void    acf_polar(const branch_t*,double**,complex_t*);    
  194. static    void    acf_polyterm(const branch_t*,double**,complex_t*);
  195. static    void    acf_pulse(const branch_t*,double**,complex_t*);    
  196. static    void    acf_pwl(const branch_t*,double**,complex_t*);    
  197. static    void    acf_sffm(const branch_t*,double**,complex_t*);
  198. static    void    acf_sin(const branch_t*,double**,complex_t*);    
  199. static    void    acf_tanh(const branch_t*,double**,complex_t*);    
  200. static    void    acf_ic(const branch_t*,double**,complex_t*);    
  201. static    void    acf_ii(const branch_t*,double**,complex_t*);    
  202. static    void    acf_iv(const branch_t*,double**,complex_t*);    
  203. static    void    acf_tempco(const branch_t*,double**,complex_t*);
  204. /*--------------------------------------------------------------------------*/
  205. #define arg0 (arg[0][0])
  206. #define arg1 (arg[0][1])
  207. #define arg2 (arg[0][2])
  208. #define arg3 (arg[0][3])
  209. #define arg4 (arg[0][4])
  210. #define arg5 (arg[0][5])
  211. #define arg6 (arg[0][6])
  212. #define arg7 (arg[0][7])
  213. #define xx0  (brh->acbias)
  214.  
  215. extern const ac_t ac;
  216. extern const char e_int[];    /* error message    */
  217. static int skip;
  218. /*--------------------------------------------------------------------------*/
  219. complex_t acfix(brh)
  220. branch_t *brh;
  221. {
  222.  complex_t y;            /* return value                    */
  223.  
  224.  skip = NO;
  225.  y.x = brh->val;
  226.  y.y = 0.;
  227.  
  228.  if (brh->x){
  229.     struct expr *x;
  230.     double *arg;
  231.     int *key;
  232.  
  233.     x = (struct expr*)brh->x;
  234.     arg = x->args->args;
  235.     for (key = x->keys->args;  *key;  key++){
  236.        switch (*key){
  237.       case eAC:        acf_ac(brh,&arg,&skip,&y);     break;
  238.       case eDC:        acf_dc(brh,&arg,&skip,&y);     break;
  239.       case eDCTRAN:        acf_dctran(brh,&arg,&skip,&y);     break;
  240.       case eFREQUENCY:    acf_frequency(brh,&arg,&skip,&y);break;
  241.       case ePERIOD:        acf_period(brh,&arg,&skip,&y);     break;
  242.       case eRAMP:        acf_ramp(brh,&arg,&skip,&y);     break;
  243.       case eTIME:        acf_time(brh,&arg,&skip,&y);     break;
  244.       case eTRAN:        acf_tran(brh,&arg,&skip,&y);     break;
  245.  
  246.       case eBANDWIDTH:    acf_bandwidth(brh,&arg,&y);     break;
  247.       case eCOMPLEX:    acf_complex(brh,&arg,&y);     break;
  248.       case eCORNERDOWN:    acf_cornerdown(brh,&arg,&y);     break;
  249.       case eCORNERUP:    acf_cornerup(brh,&arg,&y);     break;
  250.       case eDELAY:        acf_delay(brh,&arg,&y);         break;
  251.       case eEXP:        acf_exp(brh,&arg,&y);         break;
  252.       case eEXPTERM:    acf_expterm(brh,&arg,&y);     break;
  253.       case eGENERATOR:    acf_generator(brh,&arg,&y);     break;
  254.       case eMAX:        acf_max(brh,&arg,&y);         break;
  255.       case eNETFUNC:    acf_netfunc(brh,&arg,&y);     break;
  256.       case eNOTCH:        acf_notch(brh,&arg,&y);         break;
  257.       case eNUMERIC:    acf_numeric(brh,&arg,&y);     break;
  258.       case eOFFSET:        acf_offset(brh,&arg,&y);     break;
  259.       case ePOLAR:        acf_polar(brh,&arg,&y);         break;
  260.       case ePOLYTERM:    acf_polyterm(brh,&arg,&y);     break;
  261.       case ePULSE:        acf_pulse(brh,&arg,&y);         break;
  262.       case ePWL:        acf_pwl(brh,&arg,&y);         break;
  263.       case eSFFM:        acf_sffm(brh,&arg,&y);         break;
  264.       case eSIN:        acf_sin(brh,&arg,&y);         break;
  265.       case eTANH:        acf_tanh(brh,&arg,&y);         break;
  266.  
  267.       case eIC:        acf_ic(brh,&arg,&y);         break;
  268.       case eII:        acf_ii(brh,&arg,&y);         break;
  269.       case eIV:        acf_iv(brh,&arg,&y);         break;
  270.       case eTEMPCO:        acf_tempco(brh,&arg,&y);     break;
  271.  
  272.       default:  error(bWARNING, "%s: undefined function: %d\n", printlabel(brh,NO), *key); break;
  273.        }
  274.     }
  275.  }
  276.  return y;
  277. }
  278. /*--------------------------------------------------------------------------*/
  279. /* acf_ac:  keyword: ac
  280.  * following args for ac analysis only
  281.  * here : do it
  282.  */
  283. /*ARGSUSED*/
  284. static void acf_ac(brh,arg,skip,y)        /* no args        */
  285. const branch_t *brh;
  286. double **arg;
  287. int *skip;
  288. complex_t* y;
  289. {
  290.  y->x = y->y = 0.;
  291.  *skip = NO;
  292.  *arg += aAC;
  293. }
  294. /*--------------------------------------------------------------------------*/
  295. /* acf_dc:  keyword: dc
  296.  * following args for dc analysis only
  297.  * here: skip them
  298.  */
  299. /*ARGSUSED*/
  300. static void acf_dc(brh,arg,skip,y)        /* no args        */
  301. const branch_t *brh;
  302. double **arg;
  303. int *skip;
  304. complex_t *y;
  305. {
  306.  *skip = YES;
  307.  *arg += aDC;
  308. }
  309. /*--------------------------------------------------------------------------*/
  310. /* acf_dctran:  keyword: dctran
  311.  * following args for dc and transient analysis
  312.  * here: skip them
  313.  */
  314. /*ARGSUSED*/
  315. static void acf_dctran(brh,arg,skip,y)        /* no args        */
  316. const branch_t *brh;
  317. double **arg;
  318. int *skip;
  319. complex_t* y;
  320. {
  321.  *skip = YES;
  322.  *arg += aDCTRAN;
  323. }
  324. /*--------------------------------------------------------------------------*/
  325. /* acf_frequency:  keyword: frequency
  326.  * the value is frequency dependent (ac only)
  327.  * works only in ac analysis, otherwise could be non-causal.
  328.  */
  329. /*ARGSUSED*/
  330. static void acf_frequency(brh,arg,skip,y)    /* arg0 = test freq    */
  331. const branch_t *brh;
  332. double **arg;
  333. int *skip;
  334. complex_t *y;
  335. {
  336.  error(bWARNING,"frequency not implemented: %s\n", printlabel(brh,NO));
  337.  *arg += aFREQUENCY;
  338. }
  339. /*--------------------------------------------------------------------------*/
  340. /* acf_period:  keyword: period
  341.  * periodic in time
  342.  * implies that following args for transient analysis only.
  343.  * here: skip them
  344.  */
  345. /*ARGSUSED*/
  346. static void acf_period(brh,arg,skip,y)        /* arg0 = period    */
  347. const branch_t *brh;
  348. double **arg;
  349. int *skip;
  350. complex_t *y;
  351. {
  352.  *skip = YES;
  353.  *arg += aPERIOD;
  354. }
  355. /*--------------------------------------------------------------------------*/
  356. /* acf_ramp: keyword: ramp
  357.  * ramp value between plotted points
  358.  * transient only: skip here
  359.  */
  360. /*ARGSUSED*/
  361. static void acf_ramp(brh,arg,skip,y)
  362. const branch_t *brh;
  363. double **arg;
  364. int *skip;
  365. complex_t *y;
  366. {
  367.  *skip = YES;
  368.  *arg += aRAMP;
  369. }
  370. /*--------------------------------------------------------------------------*/
  371. /* acf_time:  keyword: time
  372.  * following args take effect at the given time (switch)
  373.  * obviously, transient only.
  374.  */
  375. /*ARGSUSED*/
  376. static void acf_time(brh,arg,skip,y)        /* arg0 = switch time    */
  377. const branch_t *brh;
  378. double **arg;
  379. int *skip;
  380. complex_t *y;
  381. {
  382.  *skip = YES;
  383.  *arg += aTIME;
  384. }
  385. /*--------------------------------------------------------------------------*/
  386. /* acf_tran:  keyword: transient
  387.  * following args for transient analysis only.
  388.  */
  389. /*ARGSUSED*/
  390. static void acf_tran(brh,arg,skip,y)        /* no args        */
  391. const branch_t *brh;
  392. double **arg;
  393. int *skip;
  394. complex_t *y;
  395. {
  396.  *skip = YES;
  397.  *arg += aTRAN;
  398. }
  399. /*--------------------------------------------------------------------------*/
  400. /* acf_bandwidth:  function: bandwidth
  401.  * gain block, with a bandwidth. (ac only)
  402.  * bad design: should be keyword instead.
  403.  */
  404. /*ARGSUSED*/
  405. static void acf_bandwidth(brh,arg,y)        /* arg0 = dc gain    */
  406. const branch_t *brh;                /* arg1 = 3 db freq.    */
  407. double **arg;
  408. complex_t *y;
  409. {
  410.  if (!skip  &&  arg1 != 0.){
  411.     double ratio;
  412.     double coeff;
  413.     ratio = ac.freq / arg1;
  414.     coeff = arg0 / (1.+(ratio*ratio));
  415.     y->x += coeff;
  416.     y->y -= coeff * ratio;
  417.  }
  418.  *arg += aBANDWIDTH;
  419. }
  420. /*--------------------------------------------------------------------------*/
  421. /* acf_complex:  function: complex
  422.  * complex value, cartesian coordinates, in frequency domain (ac only)
  423.  */
  424. /*ARGSUSED*/
  425. static void acf_complex(brh,arg,y)        /* arg0 = real part     */
  426. const branch_t *brh;                /* arg1 = imaginary part */
  427. double **arg;
  428. complex_t *y;
  429. {
  430.  if (!skip){
  431.     y->x += arg0;
  432.     y->y += arg1;
  433.  }
  434.  *arg += aCOMPLEX;
  435. }
  436. /*--------------------------------------------------------------------------*/
  437. /* acf_cornerdown:  function: cornerdown
  438.  * piecewise linear function:
  439.  * output is 0 for bias >= arg1
  440.  * derivative is arg0 for bias <= arg1
  441.  */
  442. static void acf_cornerdown(brh,arg,y)        /* arg0 = active slope    */
  443. const branch_t *brh;                /* arg1 = break point    */
  444. double **arg;
  445. complex_t *y;
  446. {
  447.  if (!skip  &&  xx0 < arg1)
  448.     y->x += arg0;
  449.  *arg += aCORNERDOWN;
  450. }
  451. /*--------------------------------------------------------------------------*/
  452. /* acf_cornerup:  function: cornerup
  453.  * piecewise linear function:
  454.  * output is 0 for bias <= arg1
  455.  * derivative is arg0 for bias >= arg1
  456.  */
  457. static void acf_cornerup(brh,arg,y)        /* arg0 = active slope    */
  458. const branch_t *brh;                /* arg1 = break point    */
  459. double **arg;
  460. complex_t *y;
  461. {
  462.  if (!skip  &&  xx0 > arg1)
  463.     y->x += arg0;
  464.  *arg += aCORNERUP;
  465. }
  466. /*--------------------------------------------------------------------------*/
  467. /* acf_delay:  function: delay
  468.  * time delay (ac only)
  469.  */
  470. static void acf_delay(brh,arg,y)        /* arg0 = magnitude    */
  471. const branch_t *brh;                /* arg1 = delay        */
  472. double **arg;
  473. complex_t *y;
  474. {
  475.  if (!skip){
  476.     double ratio;
  477.     double dphase;
  478.     ratio = ac.freq * arg1;
  479.     if (ratio > 100000.){
  480.        error(bPICKY, "delay too long: %s\n", printlabel(brh,NO));
  481.        ratio = 0.;
  482.     }
  483.     dphase = -360.0 * ratio;
  484.     y->x += real(arg0, dphase);
  485.     y->y += imag(arg0, dphase);
  486.  }
  487.  *arg += aDELAY;
  488. }
  489. /*--------------------------------------------------------------------------*/
  490. /* acf_exp:  spice compatible function: exp
  491.  * spice source exponential function: exponential function of time
  492.  * no-op in ac analysis, with a complaint
  493.  */
  494. /*ARGSUSED*/
  495. static void acf_exp(brh,arg,y)            /* arg0 = initial value      */
  496. const branch_t *brh;                /* arg1 = pulsed value      */
  497. double **arg;                    /* arg2 = rise delay      */
  498. complex_t *y;                    /* arg3 = rise time const */
  499. {                        /* arg4 = fall delay      */
  500.                         /* arg5 = fall time const */
  501.  if (!skip)
  502.     error(bDEBUG,"%s: exp not supported in ac analysis\n",printlabel(brh,NO));
  503.  *arg += aEXP;
  504. }
  505. /*--------------------------------------------------------------------------*/
  506. /* acf_expterm:  function: expterm
  507.  * exponent term, non-integer power term (not exponential)
  508.  * like polyterm, but for non-integers
  509.  * only works for positive input, because negative number raised to
  510.  * non-integer power is usually complex.
  511.  */
  512. static void acf_expterm(brh,arg,y)        /* arg0 = coefficient    */
  513. const branch_t *brh;                /* arg1 = exponent    */
  514. double **arg;
  515. complex_t *y;
  516. {
  517.  if (!skip  &&  xx0 > 0.){
  518.     double coeff;
  519.     coeff = arg0 * pow(xx0,arg1-1);
  520.     y->x += coeff * arg1;
  521.  }
  522.  *arg += aEXPTERM;
  523. }
  524. /*--------------------------------------------------------------------------*/
  525. /* acf_generator:  function: generator
  526.  * value is derived from the "signal generator" (generator command)
  527.  * intended for fixed sources, as circuit input.
  528.  */
  529. /*ARGSUSED*/
  530. static void acf_generator(brh,arg,y)        /* arg0 = scale factor    */
  531. const branch_t *brh;
  532. double **arg;
  533. complex_t *y;
  534. {
  535.  if (!skip){
  536.     y->x += arg0;
  537.     /* y->y += 0.; */
  538.  }
  539.  *arg += aGENERATOR;
  540. }
  541. /*--------------------------------------------------------------------------*/
  542. /* acf_max:  function: max
  543.  * piecewise linear function: block that clips
  544.  * bad design: should be keyword, so it can apply to the total part
  545.  *           should be able to specify upper and lower limits
  546.  */
  547. static void acf_max(brh,arg,y)            /* arg0 = normal value     */
  548. const branch_t *brh;                /* arg1 = output clip pt */
  549. double **arg;
  550. complex_t *y;
  551. {
  552.  if (!skip  &&  arg0 != 0.){
  553.     double clip_input;
  554.     clip_input = arg1/arg0;
  555.     if (xx0 >= -clip_input   &&   xx0 <= clip_input)
  556.        y->x += arg0;
  557.  }
  558.  *arg += aMAX;
  559. }
  560. /*--------------------------------------------------------------------------*/
  561. /* acf_netfunc:  function: netfunction
  562.  * gain block, nth order response. (ac only)
  563.  * BUG: only partial checking for number of args
  564.  */
  565. /*ARGSUSED*/                    /* arg0 = arg count    */
  566. static void acf_netfunc(brh,arg,y)        /* arg1 = dc gain    */
  567. const branch_t *brh;                /* arg2 = ref frequency    */
  568. double **arg;                    /* arg3 = order of denom*/
  569. complex_t *y;                    /* arg* = coefs of denom*/
  570. {                        /* arg* = coefs of numer*/
  571.  int argcount;
  572.  argcount = (int)arg0;
  573.  
  574.  if (!skip  &&  arg1 != 0.){
  575.     double *denomcoeff;    /* coefficients of denominator            */
  576.     double *numercoeff;    /* coefficients of numerator            */
  577.     int denomorder;    /* order of denominator (#coefs = order + 1)    */
  578.     int numerorder;    /* order of numerator                */
  579.     double wacc;    /* stash for powers of f            */
  580.     complex_t denom;    /* denominator                    */
  581.     complex_t numer;    /* numerator                    */
  582.     complex_t result;    /* result of evaluating the function        */
  583.     int ii;        /* generic loop index                */
  584.     double normfreq;    /* normalized frequency                */
  585.     double nf2;        /*    "    "    squared            */
  586.  
  587.     normfreq = ac.freq / ((arg2 != 0.) ? arg2 : 1./PI2);
  588.     nf2 = normfreq * normfreq;
  589.  
  590.     denomorder = (int)arg3;
  591.     numerorder = argcount - denomorder - 6;
  592.     if (numerorder > denomorder){
  593.        numerorder = denomorder;
  594.        error(bWARNING, "%s: too many args\n", printlabel(brh,NO));
  595.     }else if (numerorder < 0){
  596.        denomorder += numerorder;
  597.        numerorder = 0;
  598.        error(bWARNING, "%s: too few args\n", printlabel(brh,NO));
  599.     }
  600.  
  601.     denomcoeff = &(arg[0][4]);
  602.     numercoeff = &(arg[0][5+denomorder]);
  603.  
  604.     wacc = 1.;
  605.     denom.x = denomcoeff[0]; /* * wacc */
  606.     for (ii = 2;   ii <= denomorder;   ii += 2){
  607.        wacc *= -nf2;
  608.        denom.x += denomcoeff[ii] * wacc;
  609.     }
  610.  
  611.     wacc = normfreq;
  612.     denom.y = denomcoeff[1] * wacc;
  613.     for (ii = 3;   ii <= denomorder;   ii += 2){
  614.        wacc *= -nf2;
  615.        denom.y += denomcoeff[ii] * wacc;
  616.     }
  617.  
  618.     wacc = 1.;
  619.     numer.x = numercoeff[0]; /* * wacc */
  620.     for (ii = 2;   ii <= numerorder;   ii += 2){
  621.        wacc *= -nf2;
  622.        numer.x += numercoeff[ii] * wacc;
  623.     }
  624.  
  625.     wacc = normfreq;
  626.     numer.y = numercoeff[1] * wacc;
  627.     for (ii = 3;   ii <= numerorder;   ii += 2){
  628.        wacc *= -nf2;
  629.        numer.y += numercoeff[ii] * wacc;
  630.     }
  631.  
  632.     result = cdiv(numer, denom);
  633.     y->x += result.x * arg1;
  634.     y->y += result.y * arg1;
  635.  }
  636.  *arg += argcount;
  637. }
  638. /*--------------------------------------------------------------------------*/
  639. /* acf_notch:  function: notch
  640.  * piecewise linear function: crossover notch, dead zone
  641.  * symmetric around zero
  642.  */
  643. static void acf_notch(brh,arg,y)        /* arg0 = normal value     */
  644. const branch_t *brh;                /* arg1 = dead zone size */
  645. double **arg;
  646. complex_t *y;
  647. {
  648.  if (!skip){
  649.     if (xx0 > arg1   ||   xx0 < -arg1)
  650.        y->x += arg0;
  651.  }
  652.  *arg += aNOTCH;
  653. }
  654. /*--------------------------------------------------------------------------*/
  655. /* acf_numeric:  simple numeric argument
  656.  */
  657. /*ARGSUSED*/
  658. static void acf_numeric(brh,arg,y)        /* arg0 = value        */
  659. const branch_t *brh;
  660. double **arg;
  661. complex_t *y;
  662. {
  663.  if (!skip)
  664.     y->x += arg0;
  665.  *arg += aNUMERIC;
  666. }
  667. /*--------------------------------------------------------------------------*/
  668. /* acf_offset:  function: offset
  669.  * fixed dc offset
  670.  * bad design: should be keyword
  671.  */
  672. /*ARGSUSED*/
  673. static void acf_offset(brh,arg,y)        /* arg0 = gain        */
  674. const branch_t *brh;                /* arg1 = output offset    */
  675. double **arg;
  676. complex_t *y;
  677. {
  678.  if (!skip)
  679.     y->x += arg0;
  680.  *arg += aOFFSET;
  681. }
  682. /*--------------------------------------------------------------------------*/
  683. /* acf_polar:  function: polar
  684.  * complex value in polar coordinates, in frequency domain (ac only)
  685.  */
  686. /*ARGSUSED*/
  687. static void acf_polar(brh,arg,y)        /* arg0 = magnitude    */
  688. const branch_t *brh;                /* arg1 = phase        */
  689. double **arg;
  690. complex_t *y;
  691. {
  692.  if (!skip){
  693.     y->x += real(arg0, arg1);
  694.     y->y += imag(arg0, arg1);
  695.  }
  696.  *arg += aPOLAR;
  697. }
  698. /*--------------------------------------------------------------------------*/
  699. /* acf_polyterm:  function: polyterm
  700.  * polynomial term, one term in a polynomial
  701.  * power must be integer, is rounded to nearest integer
  702.  * caution: will divide by zero with zero input and negative exponent
  703.  */
  704. static void acf_polyterm(brh,arg,y)        /* arg0 = coefficient    */
  705. const branch_t *brh;                /* arg1 = exponent    */
  706. double **arg;
  707. complex_t *y;
  708. {
  709.  if (!skip){
  710.     int expo;
  711.     expo = (int)floor(arg1+.5);
  712.     if (expo != 0){
  713.        double coeff;
  714.        coeff = arg0 * ipow(xx0,expo-1);
  715.        y->x += coeff * expo;
  716.     }
  717.  }
  718.  *arg += aPOLYTERM;
  719. }
  720. /*--------------------------------------------------------------------------*/
  721. /* acf_pulse: spice compatible function: pulse
  722.  * spice pulse function (for sources)
  723.  * no-op in ac analysis, with a complaint
  724.  */
  725. /*ARGSUSED*/
  726. static void acf_pulse(brh,arg,y)        /* arg0 = initial value    */
  727. const branch_t *brh;                /* arg1 = pulsed value    */
  728. double **arg;                    /* arg2 = delay time    */
  729. complex_t *y;                    /* arg3 = rise time    */
  730. {                        /* arg4 = fall time    */
  731.                         /* arg5 = pulse width    */
  732.  if (!skip)                    /* arg6 = period    */
  733.    error(bDEBUG,"%s: pulse not supported in ac analysis\n",printlabel(brh,NO));
  734.  *arg += aPULSE;
  735. }
  736. /*--------------------------------------------------------------------------*/
  737. /* acf_pwl:  spice compatible function: pwl
  738.  * "piece-wise linear" point-plotting function of time.
  739.  * no-op in ac analysis, with a complaint
  740.  */
  741. /*ARGSUSED*/
  742. static void acf_pwl(brh,arg,y)            /* arg0 = time        */
  743. const branch_t *brh;                /* arg1 = value        */
  744. double **arg;                    /* etc.            */
  745. complex_t *y;
  746. {
  747.  if (!skip)
  748.     error(bDEBUG,"%s: pwl not supported in ac analysis\n",printlabel(brh,NO));
  749.  
  750.  if (aPWL == aVARIABLE){
  751.     *arg += (int)(**arg);
  752.  }else{
  753.     *arg += aPWL;
  754.     error(bWARNING, e_int, "acf_pwl");
  755.  }
  756. }
  757. /*--------------------------------------------------------------------------*/
  758. /* acf_sffm:  spice compatible function: sffm
  759.  * single frequency frequency modulation
  760.  * no-op in ac analysis, with a complaint
  761.  */
  762. /*ARGSUSED*/
  763. static void acf_sffm(brh,arg,y)            /* arg0 = dc offset    */
  764. const branch_t *brh;                /* arg1 = amplitude    */
  765. double **arg;                    /* arg2 = carrier freq    */
  766. complex_t *y;                    /* arg3 = mod index    */
  767. {                        /* arg4 = signal freq    */
  768.  if (!skip)
  769.     error(bDEBUG,"%s: sffm not supported in ac analysis\n",printlabel(brh,NO));
  770.  *arg += aSFFM;
  771. }
  772. /*--------------------------------------------------------------------------*/
  773. /* acf_sin:  spice compatible function: sin
  774.  * sinusoidal function, mainly for sources
  775.  * no-op in ac analysis, with a complaint
  776.  */
  777. /*ARGSUSED*/
  778. static void acf_sin(brh,arg,y)            /* arg0 = offset    */
  779. const branch_t *brh;                /* arg1 = amplitude    */
  780. double **arg;                    /* arg2 = frequency    */
  781. complex_t *y;                    /* arg3 = delay        */
  782. {                        /* arg4 = damping    */
  783.  if (!skip)
  784.     error(bDEBUG,"%s: sin not supported in ac analysis\n",printlabel(brh,NO));
  785.  *arg += aSIN;
  786. }
  787. /*--------------------------------------------------------------------------*/
  788. /* acf_tanh:  function: tanh
  789.  * for now, a copy of max
  790.  * piecewise linear function: block that clips
  791.  * bad design: should be keyword, so it can apply to the total part
  792.  *           should be able to specify upper and lower limits
  793.  */
  794. static void acf_tanh(brh,arg,y)            /* arg0 = normal value     */
  795. const branch_t *brh;                /* arg1 = output clip pt */
  796. double **arg;
  797. complex_t *y;
  798. {
  799.  if (!skip){
  800.     if (arg1 != 0.){
  801.        double cosine;
  802.        cosine = cosh(xx0 * arg0/arg1);
  803.        y->x += arg0 / (cosine*cosine);
  804.     }
  805.     /* else 0 */
  806.  }
  807.  *arg += aTANH;
  808. }
  809. /*--------------------------------------------------------------------------*/
  810. /*ARGSUSED*/
  811. static void acf_ic(brh,arg,y)
  812. const branch_t *brh;
  813. double **arg;
  814. complex_t *y;
  815. {
  816.  *arg += 1;
  817. }
  818. /*--------------------------------------------------------------------------*/
  819. /*ARGSUSED*/
  820. static void acf_ii(brh,arg,y)
  821. const branch_t *brh;
  822. double **arg;
  823. complex_t *y;
  824. {
  825.  *arg += 1;
  826. }
  827. /*--------------------------------------------------------------------------*/
  828. /*ARGSUSED*/
  829. static void acf_iv(brh,arg,y)
  830. const branch_t *brh;
  831. double **arg;
  832. complex_t *y;
  833. {
  834.  *arg += 1;
  835. }
  836. /*--------------------------------------------------------------------------*/
  837. /*ARGSUSED*/
  838. static void acf_tempco(brh,arg,y)
  839. const branch_t *brh;
  840. double **arg;
  841. complex_t *y;
  842. {
  843. #ifdef NEVER
  844.  tempco = arg0;
  845. #endif
  846.  *arg += 1;
  847. }
  848. /*--------------------------------------------------------------------------*/
  849. /*--------------------------------------------------------------------------*/
  850. SHAR_EOF
  851. fi # end of overwriting check
  852. if test -f 'src/ac_load.c'
  853. then
  854.     echo shar: will not over-write existing file "'src/ac_load.c'"
  855. else
  856. cat << \SHAR_EOF > 'src/ac_load.c'
  857. /* ac_load  04/02/92
  858.  * Copyright 1983-1992   Albert Davis
  859.  * Load AC matrix from pre-computed values
  860.  */
  861. #include "ecah.h"
  862. #include "array.h"
  863. #include "branch.h"
  864. #include "declare.h"
  865. /*--------------------------------------------------------------------------*/
  866.     void    acloadsource(const branch_t*);
  867.     void    acloadpassive(const branch_t*);
  868.     void    acloadpassivereal(const branch_t*);
  869.     void    acloadpassiveimaginary(const branch_t*);
  870.     void    acloadactive(const branch_t*);
  871. /*--------------------------------------------------------------------------*/
  872. extern double *recon, *imcon;    /* right side vector, 2 = complex        */
  873. /*--------------------------------------------------------------------------*/
  874. void acloadsource(brh)
  875. const branch_t *brh;
  876. {
  877.  if (brh->n[OUT2].m != 0){
  878.     recon[brh->n[OUT2].m] += brh->acg.x;
  879.     imcon[brh->n[OUT2].m] += brh->acg.y;
  880.  }
  881.  if (brh->n[OUT1].m != 0){
  882.     recon[brh->n[OUT1].m] -= brh->acg.x;
  883.     imcon[brh->n[OUT1].m] -= brh->acg.y;
  884.  }
  885. }
  886. /*--------------------------------------------------------------------------*/
  887. void acloadpassive(brh)
  888. const branch_t *brh;
  889. {
  890.  if (brh->acg.x != 0.)
  891.     acloadpassivereal(brh);
  892.  if (brh->acg.y != 0.)
  893.     acloadpassiveimaginary(brh);
  894. }
  895. /*--------------------------------------------------------------------------*/
  896. void acloadpassivereal(brh)
  897. const branch_t *brh;
  898. {
  899.  if (brh->n[OUT2].m != 0){
  900.     *Red(brh->n[OUT2].m,brh->n[OUT2].m) += brh->acg.x;
  901.     if (brh->n[OUT1].m != 0){
  902.        *Red(brh->n[OUT1].m,brh->n[OUT1].m) += brh->acg.x;
  903.        *Re(brh->n[OUT1].m,brh->n[OUT2].m)  -= brh->acg.x;
  904.        *Re(brh->n[OUT2].m,brh->n[OUT1].m)  -= brh->acg.x;
  905.     }
  906.  }else if (brh->n[OUT1].m != 0){
  907.     *Red(brh->n[OUT1].m,brh->n[OUT1].m) += brh->acg.x;
  908.  }
  909. }
  910. /*--------------------------------------------------------------------------*/
  911. void acloadpassiveimaginary(brh)
  912. const branch_t *brh;
  913. {
  914.  if (brh->n[OUT2].m != 0){
  915.     *Imd(brh->n[OUT2].m,brh->n[OUT2].m) += brh->acg.y;
  916.     if (brh->n[OUT1].m != 0){
  917.        *Imd(brh->n[OUT1].m,brh->n[OUT1].m) += brh->acg.y;
  918.        *Im(brh->n[OUT1].m,brh->n[OUT2].m)  -= brh->acg.y;
  919.        *Im(brh->n[OUT2].m,brh->n[OUT1].m)  -= brh->acg.y;
  920.     }
  921.  }else if (brh->n[OUT1].m != 0){
  922.     *Imd(brh->n[OUT1].m,brh->n[OUT1].m) += brh->acg.y;
  923.  }
  924. }
  925. /*--------------------------------------------------------------------------*/
  926. void acloadactive(brh)
  927. const branch_t *brh;
  928. {
  929.  if (brh->acg.x != 0.){
  930.     *re(brh->n[OUT1].m,brh->n[IN1].m) += brh->acg.x;
  931.     *re(brh->n[OUT2].m,brh->n[IN2].m) += brh->acg.x;
  932.     *re(brh->n[OUT1].m,brh->n[IN2].m) -= brh->acg.x;
  933.     *re(brh->n[OUT2].m,brh->n[IN1].m) -= brh->acg.x;
  934.  }
  935.  if (brh->acg.y != 0.){
  936.     *im(brh->n[OUT1].m,brh->n[IN1].m) += brh->acg.y;
  937.     *im(brh->n[OUT2].m,brh->n[IN2].m) += brh->acg.y;
  938.     *im(brh->n[OUT1].m,brh->n[IN2].m) -= brh->acg.y;
  939.     *im(brh->n[OUT2].m,brh->n[IN1].m) -= brh->acg.y;
  940.  }
  941. }
  942. /*--------------------------------------------------------------------------*/
  943. /*--------------------------------------------------------------------------*/
  944. SHAR_EOF
  945. fi # end of overwriting check
  946. if test -f 'src/ac_out.c'
  947. then
  948.     echo shar: will not over-write existing file "'src/ac_out.c'"
  949. else
  950. cat << \SHAR_EOF > 'src/ac_out.c'
  951. /* ac_out  01/26/93
  952.  * Copyright 1983-1992   Albert Davis
  953.  * ac analysis output functions
  954.  */
  955. #include "ecah.h"
  956. #include "ac.h"
  957. #include "branch.h"
  958. #include "dc.h"
  959. #include "io.h"
  960. #include "mode.h"
  961. #include "probh.h"
  962. #include "status.h"
  963. #include "tr.h"
  964. #include "worst.h"
  965. #include "declare.h"
  966. /*--------------------------------------------------------------------------*/
  967.     void    ac_out(const ac_t*);
  968. static    void    ac_print(double);
  969. /*--------------------------------------------------------------------------*/
  970. extern const struct ioctrl io;
  971. extern       struct status stats;
  972.  
  973. extern const probe_t *probelist[];
  974. extern const int sim_mode;
  975. extern const int worstcase;
  976. extern const complex_t *fodata;            /* used as a flag */
  977. /*--------------------------------------------------------------------------*/
  978. void ac_out(ac)
  979. const ac_t *ac;
  980. {
  981.  time_start(&(stats.output));
  982.  plotac(ac->freq);
  983.  ac_print(ac->freq);
  984.  time_stop(&(stats.output));
  985. }
  986. /*--------------------------------------------------------------------------*/
  987. /* ac_print: print the list of results (text form) to io.where
  988.  * The argument is the first column (independent variable, aka "x")
  989.  */
  990. static void ac_print(x)
  991. double x;
  992. {
  993.  if (!io.ploton){
  994.     if (x != NOT_VALID)
  995.        mprintf(io.where,ftos(x,"           ",5,io.formaat));
  996.     if (probelist[sim_mode]){
  997.        int ii;
  998.        for (ii = 0;  *probelist[sim_mode][ii].what;  ii++){
  999.       double value;
  1000.       probe_t prb;    /* to hide MSC BUG */
  1001.       prb = probelist[sim_mode][ii];
  1002.       value = acprobe(prb);
  1003.       mprintf(io.where,ftos(value,"           ",5,io.formaat));
  1004.        }
  1005.     }
  1006.     mprintf(io.where,"\n");
  1007.  }
  1008. }
  1009. /*--------------------------------------------------------------------------*/
  1010. /*--------------------------------------------------------------------------*/
  1011. SHAR_EOF
  1012. fi # end of overwriting check
  1013. if test -f 'src/ac_probe.c'
  1014. then
  1015.     echo shar: will not over-write existing file "'src/ac_probe.c'"
  1016. else
  1017. cat << \SHAR_EOF > 'src/ac_probe.c'
  1018. /* acprobe  01/11/93
  1019.  * Copyright 1983-1992   Albert Davis
  1020.  * collection of probe functions for ac analysis
  1021.  */
  1022. #include "ecah.h"
  1023. #include "branch.h"
  1024. #include "error.h"
  1025. #include "probh.h"
  1026. #include "status.h"
  1027. #include "declare.h"
  1028. /*--------------------------------------------------------------------------*/
  1029.     double    acprobe(probe_t);
  1030. static    double    acprobe_node(probe_t);
  1031. /*--------------------------------------------------------------------------*/
  1032. extern const struct status stats;
  1033. extern const double temp;
  1034. extern const double trtime;
  1035. extern const double *recon;    /* constant vector, real part            */
  1036. extern const double *imcon;    /* constant vector, imaginary part        */
  1037. extern const int *nm;        /* node map (map internal to user nodes)    */
  1038. /*--------------------------------------------------------------------------*/
  1039. double acprobe(prb)
  1040. probe_t prb;
  1041. {
  1042.  if (prb.p.node <= stats.total_nodes){
  1043.     return acprobe_node(prb);
  1044.  }else{
  1045.     return acprobe_branch(prb.p.brh,prb.what);
  1046.  }
  1047. }
  1048. /*--------------------------------------------------------------------------*/
  1049. static double acprobe_node(prb)
  1050. probe_t prb;
  1051. {
  1052.  int dummy = 0;
  1053.  
  1054.  setmatch(prb.what,&dummy);
  1055.  if (rematch("Temperature")){
  1056.     return temp + ABS_ZERO;
  1057.  }else if (rematch("TIme")){
  1058.     return trtime;
  1059.  }else if (prb.p.node > stats.total_nodes){
  1060.     return NOT_VALID;
  1061.  }else if (rematch("V") || rematch("VM")){
  1062.     return hypot(recon[nm[prb.p.node]],imcon[nm[prb.p.node]]);
  1063.  }else if (rematch("VDB")){
  1064.     double v = hypot(recon[nm[prb.p.node]],imcon[nm[prb.p.node]]);
  1065.     return (v > VMIN)  ?  20. * log10(v)  :  DBVMIN;
  1066.  }else if (rematch("VP")){
  1067.     return phase(recon[nm[prb.p.node]],imcon[nm[prb.p.node]]);
  1068.  }else{ /* bad parameter */
  1069.     return NOT_VALID;
  1070.  }
  1071. }
  1072. /*--------------------------------------------------------------------------*/
  1073. /*--------------------------------------------------------------------------*/
  1074. SHAR_EOF
  1075. fi # end of overwriting check
  1076. if test -f 'src/ac_setup.c'
  1077. then
  1078.     echo shar: will not over-write existing file "'src/ac_setup.c'"
  1079. else
  1080. cat << \SHAR_EOF > 'src/ac_setup.c'
  1081. /* ac_setup  12/31/92
  1082.  * Copyright 1983-1992   Albert Davis
  1083.  * ac analysis setup
  1084.  */
  1085. #include "ecah.h"
  1086. #include "ac.h"
  1087. #include "argparse.h"
  1088. #include "error.h"
  1089. #include "io.h"
  1090. #include "mode.h"
  1091. #include "options.h"
  1092. #include "worst.h"
  1093. #include "declare.h"
  1094. /*--------------------------------------------------------------------------*/
  1095.     void    ac_setup(const char*,int*,ac_t*);
  1096. static    void    ac_optby(const char*,int*);
  1097. static    void    ac_optdecade(const char*,int*);
  1098. static    void    ac_optlin(const char*,int*);
  1099. static    void    ac_optoctave(const char*,int*);
  1100. static    void    ac_opttimes(const char*,int*);
  1101. /*--------------------------------------------------------------------------*/
  1102. extern       struct ioctrl io;
  1103. extern const struct options opt;
  1104. extern int worstcase;    /* worst case, monte carlo mode    (enum type, worst.h)*/
  1105. extern double temp;    /* actual ambient temperature, kelvin            */
  1106. static int needslinfix;    /* flag: lin option needs patch later (spice compat)*/
  1107. static ac_t *ac;
  1108. /*--------------------------------------------------------------------------*/
  1109. void ac_setup(cmd,cnt,params)
  1110. const char *cmd;
  1111. int  *cnt;
  1112. ac_t *params;
  1113. {
  1114.  ac = params;
  1115.  io.where |= io.mstdout;
  1116.  temp = opt.tempamb;
  1117.  io.ploton = io.plotset;
  1118.  ac->echo = worstcase = NO;
  1119.  for (;;){
  1120.     if (argparse(cmd,cnt,REPEAT,
  1121.     "*",        aFUNCTION,  ac_opttimes,
  1122.         "ACMAx",    aENUM,        &worstcase,    wMAXAC,
  1123.     "ACMIn",    aENUM,        &worstcase,    wMINAC,
  1124.     "Ambient",    aODOUBLE,   &temp,    opt.tempamb,
  1125.     "By",        aFUNCTION,  ac_optby,
  1126.         "DCMAx",    aENUM,        &worstcase,    wMAXDC,
  1127.     "DCMIn",    aENUM,        &worstcase,    wMINDC,
  1128.     "Decade",    aFUNCTION,  ac_optdecade,
  1129.     "Echo",        aENUM,        &ac->echo,    YES,
  1130.     "LAg",        aENUM,        &worstcase, wLAG,
  1131.     "LEad",        aENUM,        &worstcase, wLEAD,
  1132.         "LIn",        aFUNCTION,  ac_optlin,
  1133.     "MAx",        aENUM,        &worstcase, wMAXAC,
  1134.     "MIn",        aENUM,        &worstcase, wMINAC,
  1135.     "NOPlot",    aENUM,        &io.ploton, NO,
  1136.     "Octave",    aFUNCTION,  ac_optoctave,
  1137.     "PLot",        aENUM,        &io.ploton, YES,
  1138.     "Reftemp",    aODOUBLE,   &temp,    opt.tnom,
  1139.     "SEnsitivity",    aENUM,        &worstcase, wSENS,
  1140.     "Temperature",    aODOUBLE,   &temp,    -ABS_ZERO,
  1141.     "TImes",    aFUNCTION,  ac_opttimes,
  1142.     ""))
  1143.     ;    
  1144.     else if (isfloat(cmd[*cnt])){
  1145.        ac->start = ctof(cmd,cnt);
  1146.        ac->stop  = ctof(cmd,cnt);
  1147.        if (ac->stop==0.)
  1148.           ac->stop = ac->start;
  1149.        if (isfloat(cmd[*cnt]))
  1150.           ac_optby(cmd,cnt);
  1151.     }else if (outset(cmd,cnt,(char*)NULL,"ac ")){
  1152.     ;
  1153.     }else{
  1154.     syntax(cmd,cnt,bWARNING);
  1155.     break;
  1156.     }
  1157.  }
  1158.  initio(io.where,io.whence);
  1159.  if (worstcase==wRAND  ||  worstcase==wWORST)
  1160.     io.ploton = NO;
  1161.  if (needslinfix){                /* LIN option is # of steps.   */
  1162.      ac->step=(ac->stop-ac->start)/ac->step;/* Must compute ac->step after */
  1163.      needslinfix = NO;                /* reading start and stop,       */
  1164.  }                        /* but step must be read first */
  1165.  if (ac->step==0.){                /* for Spice compatibility       */
  1166.     ac->step = ac->stop - ac->start;
  1167.     ac->linswp = YES;
  1168.  }
  1169. }
  1170. /*--------------------------------------------------------------------------*/
  1171. static void ac_optby(cmd,cnt)
  1172. const char *cmd;
  1173. int *cnt;
  1174. {
  1175.  ac->step = ctof(cmd,cnt);
  1176.  needslinfix = NO;
  1177.  ac->linswp = YES;
  1178. }
  1179. /*--------------------------------------------------------------------------*/
  1180. static void ac_optdecade(cmd,cnt)
  1181. const char *cmd;
  1182. int *cnt;
  1183. {
  1184.  ac->step = fabs(ctof(cmd,cnt));
  1185.  if (ac->step == 0.)
  1186.     ac->step = 1.;
  1187.  ac->step = pow(10., 1./ac->step);
  1188.  needslinfix = NO;
  1189.  ac->linswp = NO;
  1190. }
  1191. /*--------------------------------------------------------------------------*/
  1192. static void ac_optlin(cmd,cnt)
  1193. const char *cmd;
  1194. int *cnt;
  1195. {
  1196.  ac->step = fabs(ctof(cmd,cnt));    /* need to fix ac->step, later    */
  1197.  if (ac->step == 0.)            /* do it at the end of ac_setup    */
  1198.     ac->step = 1.;            /* a kluge, but this is a patch    */
  1199.  needslinfix = YES;            /* and I am too lazy to do it    */
  1200.  ac->linswp = YES;            /* right.            */
  1201. }
  1202. /*--------------------------------------------------------------------------*/
  1203. static void ac_optoctave(cmd,cnt)
  1204. const char *cmd;
  1205. int *cnt;
  1206. {
  1207.  ac->step = fabs(ctof(cmd,cnt));
  1208.  if (ac->step == 0.)
  1209.     ac->step = 1.;
  1210.  ac->step = pow(2.00000001, 1./ac->step);
  1211.  needslinfix = NO;
  1212.  ac->linswp = NO;
  1213. }
  1214. /*--------------------------------------------------------------------------*/
  1215. static void ac_opttimes(cmd,cnt)
  1216. const char *cmd;
  1217. int *cnt;
  1218. {
  1219.  ac->step = fabs(ctof(cmd,cnt));
  1220.  if (ac->step == 0.   &&   ac->start != 0.)
  1221.     ac->step = ac->stop / ac->start;
  1222.  needslinfix = NO;
  1223.  ac->linswp = NO;
  1224. }
  1225. /*--------------------------------------------------------------------------*/
  1226. /*--------------------------------------------------------------------------*/
  1227. SHAR_EOF
  1228. fi # end of overwriting check
  1229. if test -f 'src/ac_sweep.c'
  1230. then
  1231.     echo shar: will not over-write existing file "'src/ac_sweep.c'"
  1232. else
  1233. cat << \SHAR_EOF > 'src/ac_sweep.c'
  1234. /* ac_sweep  01/26/93
  1235.  * Copyright 1983-1992   Albert Davis
  1236.  * ac analysis sweep
  1237.  */
  1238. #include "ecah.h"
  1239. #include "ac.h"
  1240. #include "io.h"
  1241. #include "declare.h"
  1242. /*--------------------------------------------------------------------------*/
  1243.     void    ac_sweep(ac_t*);
  1244. static    void    ac_first(ac_t*);
  1245. static    int    ac_next(ac_t*);
  1246. /*--------------------------------------------------------------------------*/
  1247. extern const struct ioctrl io;
  1248. /*--------------------------------------------------------------------------*/
  1249. void ac_sweep(ac)
  1250. ac_t *ac;
  1251. {
  1252.  tr_head(ac->start, ac->stop, ac->linswp, "Freq");
  1253.  ac_first(ac);
  1254.  do {
  1255.     ac_solve();
  1256.     ac_out(ac);
  1257.  } while (ac_next(ac));
  1258. }
  1259. /*--------------------------------------------------------------------------*/
  1260. static void ac_first(ac)
  1261. ac_t *ac;
  1262. {
  1263.  ac->freq = ac->start;
  1264.  ac->omega = ac->freq * PI2;
  1265. }
  1266. /*--------------------------------------------------------------------------*/
  1267. static int ac_next(ac)
  1268. ac_t *ac;
  1269. {
  1270. double stop;
  1271.  
  1272.  if (io.whence)
  1273.     return YES;
  1274.  stop = (ac->linswp)
  1275.     ? ac->stop - ac->step/100.
  1276.     : ac->stop / pow(ac->step,.01);
  1277.  if (!inorder(ac->start, ac->freq, stop))
  1278.     return NO;
  1279.  ac->freq = (ac->linswp)
  1280.     ? ac->freq + ac->step
  1281.     : ac->freq * ac->step;
  1282.  ac->omega = ac->freq * PI2;
  1283.  if (inorder(ac->freq, ac->start, ac->stop))
  1284.     return NO;
  1285.  return YES;
  1286. }
  1287. /*--------------------------------------------------------------------------*/
  1288. /*--------------------------------------------------------------------------*/
  1289. SHAR_EOF
  1290. fi # end of overwriting check
  1291. if test -f 'src/ac_z.c'
  1292. then
  1293.     echo shar: will not over-write existing file "'src/ac_z.c'"
  1294. else
  1295. cat << \SHAR_EOF > 'src/ac_z.c'
  1296. /* acz    01/11/93
  1297.  * Copyright 1983-1992   Albert Davis
  1298.  * impedance at a port
  1299.  */
  1300. #include "ecah.h"
  1301. #include "error.h"
  1302. #include "status.h"
  1303. #include "declare.h"
  1304. /*--------------------------------------------------------------------------*/
  1305.     complex_t    acz(int,int,complex_t);
  1306. /*--------------------------------------------------------------------------*/
  1307. extern const struct status stats;
  1308. extern const char e_om[];
  1309.  
  1310. extern double *recon;
  1311. extern double *imcon;
  1312. /*--------------------------------------------------------------------------*/
  1313. complex_t acz(n1,n2,parallel)
  1314. int n1,n2;                /* node pair                */
  1315. complex_t parallel;            /* parallel admittance to remove        */
  1316. {
  1317.  static double *rezap, *imzap;
  1318.  static unsigned gotsize;
  1319.  unsigned newsize;
  1320.  double *resave, *imsave;            /* temporary storage: normal right side */
  1321.  complex_t raw_z;            /* raw impedance (includes par branch)  */
  1322.  
  1323.  newsize = ((unsigned)stats.total_nodes+2) * sizeof(double);
  1324.  if (!rezap || !imzap){
  1325.     rezap = (double*)calloc(newsize,1);
  1326.     imzap = (double*)calloc(newsize,1);
  1327.     gotsize = newsize;
  1328.  }else if (newsize > gotsize){
  1329.     rezap = (double*)realloc((void*)rezap,newsize);
  1330.     imzap = (double*)realloc((void*)imzap,newsize);
  1331.     gotsize = newsize;
  1332.  }
  1333.  if (!rezap || !imzap)
  1334.     error(bERROR, e_om, "acz");
  1335.  
  1336.  resave = recon;
  1337.  imsave = imcon;
  1338.  
  1339.  (void)memset((void*)recon,0,newsize);
  1340.  (void)memset((void*)imcon,0,newsize);
  1341.  if (n1)
  1342.     recon[n1] =     1.;
  1343.  if (n2)
  1344.     recon[n2] = -1.;
  1345.  
  1346.  xsolve();
  1347.  raw_z.x = recon[n1] - recon[n2];
  1348.  raw_z.y = imcon[n1] - imcon[n2];
  1349.  
  1350.  recon = resave;
  1351.  imcon = imsave;
  1352.  
  1353.  /* return  ckt_z */
  1354.  /* return  1 / ((1/raw_z) - parallel) */
  1355.  return  cflip(csub(cflip(raw_z),parallel));
  1356. }
  1357. /*--------------------------------------------------------------------------*/
  1358. /*--------------------------------------------------------------------------*/
  1359. SHAR_EOF
  1360. fi # end of overwriting check
  1361. if test -f 'src/allocate.c'
  1362. then
  1363.     echo shar: will not over-write existing file "'src/allocate.c'"
  1364. else
  1365. cat << \SHAR_EOF > 'src/allocate.c'
  1366. /* allocate  01/11/93
  1367.  * Copyright 1983-1992   Albert Davis
  1368.  * Allocates space for the admittance matrix.  Allocation below the
  1369.  * diagonal is by row, above the diagonal is by column, and stored backwards.
  1370.  * This broken vector is stored.  The length of it increases with increasing
  1371.  * position.  The maximum length of the nth vector is 2n-1.  For a band matrix
  1372.  * only those elements that are non-zero or are nearer to the diagonal than a
  1373.  * non-zero element are stored.
  1374.  */
  1375. #include "ecah.h"
  1376. #include "branch.h"
  1377. #include "dev.h"
  1378. #include "error.h"
  1379. #include "mode.h"
  1380. #include "nodestat.h"
  1381. #include "options.h"
  1382. #include "status.h"
  1383. #include "declare.h"
  1384. /*--------------------------------------------------------------------------*/
  1385.     void     allocate(int);
  1386. static    void     count_nodes(void);
  1387. static    void     expand_list(branch_t*);
  1388. static    void     map_nodes(void);
  1389. static    void     order_reverse(void);
  1390. static    void     order_forward(void);
  1391. static    void     order_auto(void);
  1392. static    void     map_list(branch_t*);
  1393. static    void     build_basenode_table(void);
  1394. static    void     addto_basenode_table(const branch_t*);
  1395. static    void     lownode(int,int);
  1396. static    void     alloc_hold_vectors(void);
  1397. static    unsigned size_arrays(void);
  1398. static    void     alloc_matrix(int,unsigned);
  1399. static    void     alloc_constants(int);
  1400. static    void     alloc_pointer_table(int);
  1401. static    void     fill_pointer_table(int);
  1402. /*--------------------------------------------------------------------------*/
  1403. extern const struct options opt;
  1404. extern       struct status stats;
  1405.  
  1406. extern const char e_om[];
  1407.  
  1408. extern int decomp;            /* how far the array been decomposed    */
  1409. extern int *basnode;            /* lowest node connected to this node   */
  1410. extern int *nm;                /* node map table                */
  1411. extern unsigned aspace;            /* count of non-zero array elements        */
  1412. extern double **rerow, **imrow;        /* array of row pointers            */
  1413. extern double **recol, **imcol;        /* array of column pointers            */
  1414. extern double **redia, **imdia;        /* array of diagonal pointers        */
  1415. extern double *recon,*imcon,*oldcon,*fwcon;/* constant terms            */
  1416. extern double *reals, *imags;        /* big array allocation bases        */
  1417. extern struct nodestat *nstat;        /* node status flags            */
  1418. extern double *volts;
  1419. /*--------------------------------------------------------------------------*/
  1420. /* allocate: memory allocation for simulation
  1421.  */
  1422. void allocate(mode)
  1423. int mode;
  1424. {
  1425.  decomp = 0;                    /* start decomposition over        */
  1426.  
  1427.  if (mode != sSTATUS){
  1428.     time_reset(&(stats.load));
  1429.     time_reset(&(stats.lu));
  1430.     time_reset(&(stats.back));
  1431.     time_reset(&(stats.review));
  1432.     time_reset(&(stats.output));
  1433.  
  1434.     time_zstart(&(stats.setup));
  1435.  }
  1436.  if (!volts){                        /* if "volts" not allocated,        */
  1437.     dealloc(YES);                /* must allocate everything        */
  1438.     count_nodes();             /* and start solution over.        */
  1439.     expand_list(firstbranch_dev());
  1440.     map_nodes();
  1441.     build_basenode_table();
  1442.     alloc_hold_vectors();
  1443.     aspace = size_arrays();
  1444.     if (mode==sSTATUS)
  1445.        dealloc(YES);
  1446.  }
  1447.  if (mode==sSTATUS)
  1448.     return;                        /* done if status   */
  1449.  if (!exists(firstbranch_all()))
  1450.     error(bERROR, "no circuit\n");
  1451.  
  1452.                 /* lost between commands.  re-- every time  */
  1453.  alloc_matrix(mode,aspace);
  1454.  alloc_constants(mode);
  1455.  alloc_pointer_table(mode);
  1456.  fill_pointer_table(mode);
  1457.  if (mode != sSTATUS)
  1458.     time_stop(&(stats.setup));
  1459. }
  1460. /*--------------------------------------------------------------------------*/
  1461. /* count_nodes: count nodes in main ckt (not subckts)
  1462.  * update the variables "stats.total_nodes" and "stats.user_nodes"
  1463.  * zeros "stats.subckt_nodes" and "stats.model_nodes"
  1464.  */
  1465. static void count_nodes()
  1466. {
  1467.  const branch_t *stop;
  1468.  
  1469.  stats.user_nodes = 0;
  1470.  stop = firstbranch_dev();
  1471.  if (isdevice(stop)){
  1472.     const branch_t *brh;
  1473.     brh = stop;
  1474.     do {
  1475.        int ii;
  1476.        for (ii = 0;  brh->n[ii].e != INVALIDNODE;  ii++){
  1477.           if (brh->n[ii].e > stats.user_nodes)
  1478.          stats.user_nodes = brh->n[ii].e;
  1479.        }
  1480.     } while (brh = nextbranch_dev(brh),  brh != stop);
  1481.  }
  1482.  stats.total_nodes = stats.user_nodes;
  1483.  stats.subckt_nodes = stats.model_nodes = 0;
  1484. }
  1485. /*--------------------------------------------------------------------------*/
  1486. /* expand_list: expand (flatten) a list of components (subckts)
  1487.  * Scan component list.  Expand each subckt: create actual elements
  1488.  * for flat representation to use for simulation.
  1489.  * Recursive to allow for nested subckts.
  1490.  */
  1491. static void expand_list(stop)
  1492. branch_t *stop;
  1493. {
  1494.  if (isdevice(stop)){
  1495.     branch_t *brh;
  1496.     brh = stop;
  1497.     do {
  1498.        expand_branch(brh);
  1499.        if (brh->subckt){
  1500.           expand_list(brh->subckt);
  1501.        }
  1502.     } while (brh=nextbranch_dev(brh),  brh != stop);
  1503.  }
  1504. }
  1505. /*--------------------------------------------------------------------------*/
  1506. /* map_nodes: map intermediate node number to internal node number.
  1507.  * Ideally, this function would find some near-optimal order
  1508.  * and squash out gaps.
  1509.  */
  1510. static void map_nodes()
  1511. {
  1512.  nm = (int*)calloc((unsigned)stats.total_nodes+1, sizeof(int));
  1513.  if (!nm)                /* node map table            */
  1514.     error(bERROR, e_om, "nm");
  1515.  time_zstart(&(stats.order));
  1516.  switch (opt.order){
  1517.     default:
  1518.        error(bWARNING, "invalid order spec: %d\n", opt.order);
  1519.        /* fall through to order_reverse. */
  1520.     case oAUTO:
  1521.        order_auto();
  1522.        break;
  1523.     case oREVERSE:
  1524.        order_reverse();
  1525.        break;
  1526.     case oFORWARD:
  1527.        order_forward();
  1528.        break;
  1529.  }
  1530.  time_stop(&(stats.order));
  1531.  map_list(firstbranch_dev());
  1532. }
  1533. /*--------------------------------------------------------------------------*/
  1534. /* order_reverse: force ordering to reverse of user ordering
  1535.  *  subcircuits at beginning, results on border at the bottom
  1536.  */
  1537. static void order_reverse()
  1538. {
  1539.  int node;
  1540.  nm[0] = 0;
  1541.  for (node=1;  node<=stats.total_nodes;  ++node)
  1542.     nm[node] = stats.total_nodes - node + 1;
  1543. }
  1544. /*--------------------------------------------------------------------------*/
  1545. /* order_forward: use user ordering, with subcircuits added to end
  1546.  * results in border at the top (worst possible if lots of subcircuits)
  1547.  */
  1548. static void order_forward()
  1549. {
  1550.  int node;
  1551.  nm[0] = 0;
  1552.  for (node=1;  node<=stats.total_nodes;  ++node)
  1553.     nm[node] = node;
  1554. }
  1555. /*--------------------------------------------------------------------------*/
  1556. /* order_auto: full automatic ordering
  1557.  * reverse, for now
  1558.  */
  1559. static void order_auto()
  1560. {
  1561.  int node;
  1562.  nm[0] = 0;
  1563.  for (node=1;  node<=stats.total_nodes;  ++node)
  1564.     nm[node] = stats.total_nodes - node + 1;
  1565. }
  1566. /*--------------------------------------------------------------------------*/
  1567. /* map_list: map intermediate node to internal (working) node number.
  1568.  * recursive: actual mapping done here
  1569.  */
  1570. static void map_list(stop)
  1571. branch_t *stop;
  1572. {
  1573.  if (isdevice(stop)){
  1574.     branch_t *brh;
  1575.     brh = stop;
  1576.     do {
  1577.        int ii;
  1578.        for (ii = 0;  brh->n[ii].t != INVALIDNODE;  ii++){
  1579.       brh->n[ii].m = nm[brh->n[ii].t];
  1580.        }
  1581.        brh->n[ii].m = INVALIDNODE;
  1582.        if (brh->subckt)
  1583.           map_list(brh->subckt);
  1584.     } while (brh = nextbranch_dev(brh),  brh != stop);
  1585.  }
  1586. }
  1587. /*--------------------------------------------------------------------------*/
  1588. /* build_basenode_table: build a table containing the lowest node number
  1589.  * connected to each node.  This determines the bandwidth of the matrix,
  1590.  * or the length of the vectors (in the matrix) corresponding to each node.
  1591.  * Actual work is done by addto_basenode_table.
  1592.  */
  1593. static void build_basenode_table()
  1594. {
  1595.  int node;
  1596.  
  1597.  basnode = (int*)calloc((unsigned)stats.total_nodes+1, sizeof(int));
  1598.  if (!basnode)                /* base-node, (vector size) table        */
  1599.     error(bERROR, e_om, "basnode");
  1600.  
  1601.  for (node = 0;  node <= stats.total_nodes;  ++node)
  1602.     basnode[node] = node;        /* initialize each to itself        */
  1603.  
  1604.  addto_basenode_table(firstbranch_dev());
  1605. }
  1606. /*--------------------------------------------------------------------------*/
  1607. /* addto_basenode_table: actual work of build_basenode_table.
  1608.  * Find and store "base-node" for each node.  This is the lowest numbered
  1609.  * node that it connects to.  We don't do anything with higher numbered nodes.
  1610.  * (Unnecessary because of symmetry.)   Recursive, for subckts.
  1611.  */
  1612. static void addto_basenode_table(stop)
  1613. const branch_t *stop;
  1614. {
  1615.  if (isdevice(stop)){
  1616.     const branch_t *brh;
  1617.     brh = stop;
  1618.     do {
  1619.        int ii,jj;
  1620.        for (ii = 0;  brh->n[ii].m != INVALIDNODE;  ii++){
  1621.           if (brh->n[ii].m != 0)
  1622.          for (jj = 0;  jj < ii ;  jj++)
  1623.             lownode(brh->n[ii].m,brh->n[jj].m);
  1624.        }
  1625.        if (brh->subckt){
  1626.           addto_basenode_table(brh->subckt);   /* do subckts */
  1627.        }
  1628.     } while (brh = nextbranch_dev(brh), brh != stop);
  1629.  }
  1630. }
  1631. /*--------------------------------------------------------------------------*/
  1632. /* lownode: find the lowest node that each node connects to.
  1633.  * (connects to == has distance 1 from)
  1634.  */
  1635. static void lownode(node1,node2)
  1636. int node1,node2;
  1637. {
  1638.  if (node1 == 0  ||  node2 == 0)    /* node 0 is ground, and doesn't    */
  1639.     ;                    /* count as a connection        */
  1640.  else if ( node1 < basnode[node2] ) 
  1641.     basnode[node2]=node1;
  1642.  else if ( node2 < basnode[node1] ) 
  1643.     basnode[node1]=node2;
  1644. }
  1645. /*--------------------------------------------------------------------------*/
  1646. /* alloc_hold_vectors: allocate space to hold data between commands.
  1647.  * store voltage between commands, for restart and convergence assistance
  1648.  */
  1649. static void alloc_hold_vectors()
  1650. {
  1651.  volts = (double*)calloc((unsigned)stats.total_nodes+2, sizeof(double));
  1652.  if (!volts)
  1653.     error(bERROR, e_om, "volts");
  1654. }
  1655. /*--------------------------------------------------------------------------*/
  1656. /* size_arrays: determine memory needed for main matrix
  1657.  * (count non-zero array elememts)
  1658.  */
  1659. static unsigned size_arrays()
  1660. {
  1661.  int node;
  1662.  unsigned space;
  1663.  
  1664.  space = 0;
  1665.  for (node = 0;   node <= stats.total_nodes;   ++node)
  1666.     space += 2 * ( node-basnode[node] ) + 1;
  1667.  return space;
  1668. }
  1669. /*--------------------------------------------------------------------------*/
  1670. /* alloc_matrix: allocate main matrix
  1671.  * (the thing we do LU decomposition to, to solve)
  1672.  */
  1673. static void alloc_matrix(mode,space)
  1674. int mode;
  1675. unsigned space;
  1676. {
  1677.  if (space > SEGSIZ/sizeof(double))    /* attempt to trap allocation of    */
  1678.     error(bERROR, e_om, "seg size");    /* single matrix large enough to    */
  1679.                         /* cause segmentation problems        */
  1680.                     /* also traps space*sizeof(double) */
  1681.                     /* overflow inside calloc.        */
  1682.  reals = (double*)calloc(space, sizeof(double));
  1683.  if (!reals)
  1684.     error(bERROR, e_om, "reals");
  1685.  imags = (double*)calloc(space, sizeof(double));
  1686.  if (!imags)
  1687.     error(bERROR, e_om, "imags");
  1688. }
  1689. /*--------------------------------------------------------------------------*/
  1690. /* alloc_constants: allocate constant vectors
  1691.  * (constant is not the right word...)  allocate space for the right-side
  1692.  * vector (initially current sources, on solution becomes voltages) and copies
  1693.  * used for one-time-ago and convergence checking
  1694.  */
  1695. static void alloc_constants(mode)
  1696. int mode;
  1697. {
  1698.  recon = (double*)calloc((unsigned)stats.total_nodes+2, sizeof(double));
  1699.  if (!recon)                        /* constants terms  */
  1700.     error(bERROR, e_om, "recon");
  1701.  imcon = (double*)calloc((unsigned)stats.total_nodes+2, sizeof(double));
  1702.  if (!imcon)                        /* imaginary (AC)  */
  1703.     error(bERROR, e_om, "imcon");            /* or next (DC, tr) */
  1704.  fwcon = (double*)calloc((unsigned)stats.total_nodes+2, sizeof(double));
  1705.  if (!fwcon)                        /* fwd sub, (dctr)  */
  1706.     error(bERROR, e_om, "fwcon");
  1707.  oldcon = (double*)calloc((unsigned)stats.total_nodes+2, sizeof(double));
  1708.  if (!oldcon)                        /* last time (tr)   */
  1709.     error(bERROR, e_om, "oldcon");
  1710.  nstat = (struct nodestat*)calloc((unsigned)stats.total_nodes+2,
  1711.                 sizeof(struct nodestat));
  1712.  if (!nstat)
  1713.     error(bERROR, e_om, "nstat");
  1714. }
  1715. /*--------------------------------------------------------------------------*/
  1716. /* alloc_pointer_table: allocate tables for sparse matrix indexing
  1717.  */
  1718. static void alloc_pointer_table(mode)
  1719. int mode;
  1720. {
  1721.  redia = (double**)calloc((unsigned)stats.total_nodes+1, sizeof(double*));
  1722.  if (!redia)        /* pointers to diagonal */
  1723.     error(bERROR, e_om, "redia");
  1724.  rerow = (double**)calloc((unsigned)stats.total_nodes+1, sizeof(double*));
  1725.  if (!rerow)        /* row pointers: points to col 0 of each row */
  1726.     error(bERROR, e_om, "rerow");
  1727.  recol = (double**)calloc((unsigned)stats.total_nodes+1, sizeof(double*));
  1728.  if (!recol)        /* column pointers: points to row 0 */
  1729.     error(bERROR, e_om, "recol");
  1730.  
  1731.  imdia = (double**)calloc((unsigned)stats.total_nodes+1, sizeof(double*));
  1732.  if (!imdia)
  1733.     error(bERROR, e_om, "imdia");
  1734.  imrow = (double**)calloc((unsigned)stats.total_nodes+1, sizeof(double*));
  1735.  if (!imrow)
  1736.     error(bERROR, e_om, "imrow");
  1737.  imcol = (double**)calloc((unsigned)stats.total_nodes+1, sizeof(double*));
  1738.  if (!imcol)
  1739.     error(bERROR, e_om, "imcol");
  1740. }
  1741. /*--------------------------------------------------------------------------*/
  1742. /* fill_pointer_table: build tables for sparse matrix access
  1743.  * there are 3 tables: row, column and diagonal
  1744.  * row pointer points to each row, where column 0 would be if it existed,
  1745.  * even though col 0 doesn't actually exist.
  1746.  * Matrix is stored by row below the diagonal, and by column, backwards, above,
  1747.  * in such a way that the diagonal is in both the row and column.
  1748.  */
  1749. static void fill_pointer_table(mode)
  1750. int mode;
  1751. {
  1752.  double *repnt, *impnt;            /* running location of place in array   */
  1753.  int node;
  1754.  
  1755.  repnt = reals;                    /* fill in pointer tables   */
  1756.  impnt = imags;
  1757.  for (node = 0;   node <= stats.total_nodes;   ++node){
  1758.     recol[node] = repnt - basnode[node];    /* pointers to the bogus    */
  1759.     rerow[node] = recol[node] + 2*node;        /* row or column 0        */
  1760.     redia[node]    = recol[node] + node;
  1761.     repnt = &repnt[2*(node-basnode[node])+1];    /* Next free slot.        */
  1762.  
  1763.     imcol[node] = impnt - basnode[node];
  1764.     imrow[node] = imcol[node] + 2*node;
  1765.     imdia[node] = imcol[node] + node;
  1766.     impnt = &impnt[2*(node-basnode[node])+1];
  1767.  }
  1768. }
  1769. /*--------------------------------------------------------------------------*/
  1770. /*--------------------------------------------------------------------------*/
  1771. SHAR_EOF
  1772. fi # end of overwriting check
  1773. if test -f 'src/argparse.c'
  1774. then
  1775.     echo shar: will not over-write existing file "'src/argparse.c'"
  1776. else
  1777. cat << \SHAR_EOF > 'src/argparse.c'
  1778. /* argparse  04/02/92
  1779.  * Copyright 1983-1992   Albert Davis
  1780.  * parse name followed by numeric argument from command line
  1781.  * to call...  list the arguments:
  1782.  *   each is:
  1783.  *    string to match (upper case part must match, lower may match
  1784.  *        user string is not case sensitive
  1785.  *    variable type (see argparse.h)
  1786.  *    pointer to variable to fill
  1787.  *    if enumerated type: value
  1788.  *    if "2double": second argument
  1789.  *   end with either null pointer or null string
  1790.  *
  1791.  * cnt is incrmented to the next argument in the input string
  1792.  * returns YES if it did something, NO if not.
  1793.  *
  1794.  * will try to handle as many arguments as possible.
  1795.  */
  1796. #include "ecah.h"
  1797. #include "argparse.h"
  1798. #include "error.h"
  1799. #include "declare.h"
  1800. /*--------------------------------------------------------------------------*/
  1801.     int    argparse();
  1802. /*--------------------------------------------------------------------------*/
  1803. extern    const char e_int[];
  1804. /*--------------------------------------------------------------------------*/
  1805. /*VARARGS3*/
  1806. int argparse(cmd,cnt,mode,va_alist)
  1807. const char *cmd;
  1808. int *cnt;
  1809. int mode;
  1810. va_dcl /* ; */
  1811. {
  1812.  char *key;
  1813.  va_list marker;
  1814.  int type;
  1815.  int did = NO;
  1816.  
  1817.  setmatch(cmd,cnt);
  1818.  do {
  1819.     va_start(marker);
  1820.     while (key = va_arg(marker,char*),  key && *key){
  1821.        type = va_arg(marker,int);
  1822.        if (type == aDOUBLE || type == aUDOUBLE || type == aODOUBLE
  1823.                || type == a2DOUBLE || type == aSDOUBLE){
  1824.       double *arg1;
  1825.       double *arg2;
  1826.       double offset;
  1827.       double scale;
  1828.       arg1 = va_arg(marker,double*);
  1829.       arg2 = (type==a2DOUBLE) ? va_arg(marker,double*) : (double*)NULL;
  1830.       offset = (type==aODOUBLE  ||  type==aOIDOUBLE)
  1831.                       ? va_arg(marker,double)  : 0.;
  1832.       scale  = (type==aSDOUBLE  ||  type==aSIDOUBLE)
  1833.                       ? va_arg(marker,double)  : 1.;
  1834.       if (rematch(key)){
  1835.          *arg1 = (type==aUDOUBLE) ? fabs(ctof(cmd,cnt)) : ctof(cmd,cnt);
  1836.          *arg1 += offset;
  1837.          *arg1 *= scale;
  1838.          if (type==aIDOUBLE || type==aOIDOUBLE || type==aSIDOUBLE)
  1839.         if (*arg1 != 0.)
  1840.            *arg1 = 1. / *arg1;
  1841.          if (arg2)
  1842.         *arg2 = ctof(cmd,cnt);
  1843.          did = YES;
  1844.          break;
  1845.       }
  1846.        }else if (type == aINT || type == aUINT  ||  type == aFINT){
  1847.       int *arg;
  1848.       arg = va_arg(marker,int*);
  1849.       if (rematch(key)){
  1850.          switch (type){
  1851.         case aUINT:    *arg = abs(ctoi(cmd,cnt));    break;
  1852.         case aINT:    *arg = ctoi(cmd,cnt);        break;
  1853.         case aFINT:    *arg = (int)ctof(cmd,cnt);    break;
  1854.          }
  1855.          did = YES;
  1856.          break;
  1857.       }
  1858.        }else if (type == aENUM || type == aORENUM || type == aANDENUM){
  1859.       int *arg;
  1860.       int value;
  1861.       arg = va_arg(marker,int*);
  1862.       value = va_arg(marker,int);
  1863.       if (rematch(key)){
  1864.          switch (type){
  1865.         case aENUM:    *arg = value;    break;
  1866.         case aORENUM:   *arg |= value;    break;
  1867.         case aANDENUM:  *arg &= value;    break;
  1868.          }
  1869.          did = YES;
  1870.          break;
  1871.       }
  1872.        }else if (type == aFUNCTION || type == a2FUNCTION){
  1873.       void (*arg1)(const char*,int*);
  1874.       void (*arg2)(const char*,int*);
  1875.       arg1 = va_arg(marker,void*);
  1876.       arg2 = (type==a2FUNCTION) ? va_arg(marker,void*) : NULL;
  1877.       if (rematch(key)){
  1878.          (*arg1)(cmd,cnt);
  1879.          if (arg2)
  1880.         (*arg2)(cmd,cnt);
  1881.          did = YES;
  1882.          break;
  1883.       }
  1884.        }else{
  1885.       error(bERROR, e_int, "argparse");
  1886.        }
  1887.     }
  1888.     va_end(marker);
  1889.  } while (key && *key && mode==REPEAT);
  1890.  return did;
  1891. }
  1892. /*--------------------------------------------------------------------------*/
  1893. /*--------------------------------------------------------------------------*/
  1894. SHAR_EOF
  1895. fi # end of overwriting check
  1896. if test -f 'src/array.c'
  1897. then
  1898.     echo shar: will not over-write existing file "'src/array.c'"
  1899. else
  1900. cat << \SHAR_EOF > 'src/array.c'
  1901. /* array  04/02/92
  1902.  * Copyright 1983-1992   Albert Davis
  1903.  * Access main arrays.
  1904.  * Returns pointer to the requested array element.
  1905.  * If the place is not allocated, returns pointer to a zero.
  1906.  * There are several of these, to access different arrays.
  1907.  * Robust versions, with bounds checking, of sorts.
  1908.  * For fast versions (macros) see array.h.
  1909.  */
  1910. #include "ecah.h"
  1911. #include "error.h"
  1912. #include "status.h"
  1913. #include "declare.h"
  1914. /*--------------------------------------------------------------------------*/
  1915.     double    *re(int,int);
  1916.     double    *im(int,int);
  1917. /*--------------------------------------------------------------------------*/
  1918. extern const struct status stats;
  1919. extern const int *basnode;        /* lowest node connected to this    */
  1920. extern double **rerow, **imrow;        /* array of row pointers        */
  1921. extern double **recol, **imcol;        /* array of column pointers        */
  1922. extern double **redia, **imdia;        /* array of diagonal pointers        */
  1923. extern double *recon, *imcon;        /* constant terms            */
  1924. extern double zero;            /* really a constant, i hope        */
  1925. static double trash;            /* place to deposit node 0        */
  1926. extern const char e_int[];
  1927. /*--------------------------------------------------------------------------*/
  1928. double *re(row,col)
  1929. int row,col;
  1930. {
  1931.  if (col==row)
  1932.     return (redia[row]);
  1933.  if (col>row){        /* above the diagonal */
  1934.     if (row==0)
  1935.        return &trash;
  1936.     if (col>stats.total_nodes)
  1937.        return &recon[row];
  1938.     if (row<basnode[col])
  1939.        return &zero;
  1940.     return (recol[col]+row);
  1941.  }
  1942.  /* if (col<row) */{    /* below the diagonal */
  1943.     if (col==0)
  1944.        return &trash;
  1945.     if (row>stats.total_nodes)
  1946.        error(bERROR, e_int, "array");
  1947.     if (col<basnode[row])
  1948.        return &zero;
  1949.     return (rerow[row]-col);
  1950.  }
  1951. }
  1952. /*--------------------------------------------------------------------------*/
  1953. double *im(row,col)
  1954. int row,col;
  1955. {
  1956.  if (col==row)
  1957.     return (imdia[row]);
  1958.  if (col>row){        /* above the diagonal */
  1959.     if (row==0)
  1960.        return &trash;
  1961.     if (col>stats.total_nodes)
  1962.        return &imcon[row];
  1963.     if (row<basnode[col])
  1964.        return &zero;
  1965.     return (imcol[col]+row);
  1966.  }
  1967.  /* if (col<row) */{    /* below the diagonal */
  1968.     if (col==0)
  1969.        return &trash;
  1970.     if (row>stats.total_nodes)
  1971.        error(bERROR, e_int, "array");
  1972.     if (col<basnode[row])
  1973.        return &zero;
  1974.     return (imrow[row]-col);
  1975.  }
  1976. }
  1977. /*--------------------------------------------------------------------------*/
  1978. /*--------------------------------------------------------------------------*/
  1979. SHAR_EOF
  1980. fi # end of overwriting check
  1981. #    End of shell archive
  1982. exit 0
  1983.