/* ac_fix.c  94.03.13
 * Copyright 1983-1992   Albert Davis
 * Takes care of nonlinearities, behavioral modeling, etc. in ac analysis.
 */
#include "ecah.h"
#include "ac.h"
#include "branch.h"
#include "error.h"
#include "expr.h"
#include "declare.h"
/*--------------------------------------------------------------------------*/
	complex_t	acfix(branch_t*);
static	void	acf_ac(const branch_t*,double**,complex_t*);
static	void	acf_dc(const branch_t*,double**,complex_t*);
static	void	acf_dctran(const branch_t*,double**,complex_t*);
static	void	acf_frequency(const branch_t*,double**,complex_t*);
static	void	acf_period(const branch_t*,double**,complex_t*);
static	void	acf_ramp(const branch_t*,double**,complex_t*);	
static	void	acf_time(const branch_t*,double**,complex_t*);	
static	void	acf_tran(const branch_t*,double**,complex_t*);	
static	void	acf_bandwidth(const branch_t*,double**,complex_t*);
static	void	acf_complex(const branch_t*,double**,complex_t*);
static	void	acf_cornerdown(const branch_t*,double**,complex_t*);
static	void	acf_cornerup(const branch_t*,double**,complex_t*);
static	void	acf_delay(const branch_t*,double**,complex_t*);	
static	void	acf_exp(const branch_t*,double**,complex_t*);	
static	void	acf_expterm(const branch_t*,double**,complex_t*);
static	void	acf_generator(const branch_t*,double**,complex_t*);
static	void	acf_max(const branch_t*,double**,complex_t*);	
static	void	acf_netfunc(const branch_t*,double**,complex_t*);	
static	void	acf_notch(const branch_t*,double**,complex_t*);	
static	void	acf_numeric(const branch_t*,double**,complex_t*);	
static	void	acf_offset(const branch_t*,double**,complex_t*);
static	void	acf_polar(const branch_t*,double**,complex_t*);	
static	void	acf_polyterm(const branch_t*,double**,complex_t*);
static	void	acf_pulse(const branch_t*,double**,complex_t*);	
static	void	acf_pwl(const branch_t*,double**,complex_t*);	
static	void	acf_sffm(const branch_t*,double**,complex_t*);
static	void	acf_sin(const branch_t*,double**,complex_t*);	
static	void	acf_tanh(const branch_t*,double**,complex_t*);	
static	void	acf_ic(const branch_t*,double**,complex_t*);	
static	void	acf_ii(const branch_t*,double**,complex_t*);	
static	void	acf_iv(const branch_t*,double**,complex_t*);	
static	void	acf_tempco(const branch_t*,double**,complex_t*);
/*--------------------------------------------------------------------------*/
#define arg0 (arg[0][0])
#define arg1 (arg[0][1])
#define arg2 (arg[0][2])
#define arg3 (arg[0][3])
#define arg4 (arg[0][4])
#define arg5 (arg[0][5])
#define arg6 (arg[0][6])
#define arg7 (arg[0][7])
#define xx0  (brh->acbias)

extern const ac_t ac;
extern const char e_int[];	/* error message	*/
static int skip;
/*--------------------------------------------------------------------------*/
complex_t acfix(branch_t *brh)
{
 complex_t y;			/* return value				    */

 if (brh->f->devclass == rnSOURCE  ||  brh->f->devclass == rnNISOURCE){
    skip = YES;
    y.x = 0.;
    y.y = 0.;
 }else{
    skip = NO;
    y.x = brh->val;
    y.y = 0.;
 }

 if (brh->x){
    struct expr *x;
    double *arg;
    int *key;

    x = (struct expr*)brh->x;
    arg = x->args->args;
    for (key = x->keys->args;  *key;  key++){
       switch (*key){
	  case eAC:		acf_ac(brh,&arg,&y);		 break;
	  case eDC:		acf_dc(brh,&arg,&y);	 	break;
	  case eDCTRAN:		acf_dctran(brh,&arg,&y);	 break;
	  case eFREQUENCY:	acf_frequency(brh,&arg,&y);	break;
	  case ePERIOD:		acf_period(brh,&arg,&y);	 break;
	  case eRAMP:		acf_ramp(brh,&arg,&y);		 break;
	  case eTIME:		acf_time(brh,&arg,&y);		 break;
	  case eTRAN:		acf_tran(brh,&arg,&y);	 	break;

	  case eBANDWIDTH:	acf_bandwidth(brh,&arg,&y);	 break;
	  case eCOMPLEX:	acf_complex(brh,&arg,&y);	 break;
	  case eCORNERDOWN:	acf_cornerdown(brh,&arg,&y);	 break;
	  case eCORNERUP:	acf_cornerup(brh,&arg,&y);	 break;
	  case eDELAY:		acf_delay(brh,&arg,&y);		 break;
	  case eEXP:		acf_exp(brh,&arg,&y);		 break;
	  case eEXPTERM:	acf_expterm(brh,&arg,&y);	 break;
	  case eGENERATOR:	acf_generator(brh,&arg,&y);	 break;
	  case eMAX:		acf_max(brh,&arg,&y);		 break;
	  case eNETFUNC:	acf_netfunc(brh,&arg,&y);	 break;
	  case eNOTCH:		acf_notch(brh,&arg,&y);		 break;
	  case eNUMERIC:	acf_numeric(brh,&arg,&y);	 break;
	  case eOFFSET:		acf_offset(brh,&arg,&y);	 break;
	  case ePOLAR:		acf_polar(brh,&arg,&y);		 break;
	  case ePOLYTERM:	acf_polyterm(brh,&arg,&y);	 break;
	  case ePULSE:		acf_pulse(brh,&arg,&y);		 break;
	  case ePWL:		acf_pwl(brh,&arg,&y);		 break;
	  case eSFFM:		acf_sffm(brh,&arg,&y);		 break;
	  case eSIN:		acf_sin(brh,&arg,&y);		 break;
	  case eTANH:		acf_tanh(brh,&arg,&y);		 break;

	  case eIC:		acf_ic(brh,&arg,&y);		 break;
	  case eII:		acf_ii(brh,&arg,&y);		 break;
	  case eIV:		acf_iv(brh,&arg,&y);		 break;
	  case eTEMPCO:		acf_tempco(brh,&arg,&y);	 break;

	  default:  error(bWARNING, "%s: undefined function: %d\n", printlabel(brh,NO), *key); break;
       }
    }
 }
 return y;
}
/*--------------------------------------------------------------------------*/
/* acf_ac:  keyword: ac
 * following args for ac analysis only
 * here : do it
 * no args
 */
/*ARGSUSED*/
static void acf_ac(const branch_t *brh, double **arg, complex_t* y)
{
 y->x = y->y = 0.;
 skip = NO;
 *arg += aAC;
}
/*--------------------------------------------------------------------------*/
/* acf_dc:  keyword: dc
 * following args for dc analysis only
 * here: skip them
 * no args
 */
/*ARGSUSED*/
static void acf_dc(const branch_t *brh, double **arg, complex_t* y)
{
 skip = YES;
 *arg += aDC;
}
/*--------------------------------------------------------------------------*/
/* acf_dctran:  keyword: dctran
 * following args for dc and transient analysis
 * here: skip them
 * no args
 */
/*ARGSUSED*/
static void acf_dctran(const branch_t *brh, double **arg, complex_t* y)
{
 skip = YES;
 *arg += aDCTRAN;
}
/*--------------------------------------------------------------------------*/
/* acf_frequency:  keyword: frequency
 * the value is frequency dependent (ac only)
 * works only in ac analysis, otherwise could be non-causal.
 * arg0 = test freq
 */
/*ARGSUSED*/
static void acf_frequency(const branch_t *brh, double **arg, complex_t* y)
{
 error(bWARNING,"frequency not implemented: %s\n", printlabel(brh,NO));
 *arg += aFREQUENCY;
}
/*--------------------------------------------------------------------------*/
/* acf_period:  keyword: period
 * periodic in time
 * implies that following args for transient analysis only.
 * here: skip them
 * arg0 = period
 */
/*ARGSUSED*/
static void acf_period(const branch_t *brh, double **arg, complex_t* y)
{
 skip = YES;
 *arg += aPERIOD;
}
/*--------------------------------------------------------------------------*/
/* acf_ramp: keyword: ramp
 * ramp value between plotted points
 * transient only: skip here
 */
/*ARGSUSED*/
static void acf_ramp(const branch_t *brh, double **arg, complex_t* y)
{
 skip = YES;
 *arg += aRAMP;
}
/*--------------------------------------------------------------------------*/
/* acf_time:  keyword: time
 * following args take effect at the given time (switch)
 * obviously, transient only.
 * arg0 = switch time
 */
/*ARGSUSED*/
static void acf_time(const branch_t *brh, double **arg, complex_t* y)
{
 skip = YES;
 *arg += aTIME;
}
/*--------------------------------------------------------------------------*/
/* acf_tran:  keyword: transient
 * following args for transient analysis only.
 * no args
 */
/*ARGSUSED*/
static void acf_tran(const branch_t *brh, double **arg, complex_t* y)
{
 skip = YES;
 *arg += aTRAN;
}
/*--------------------------------------------------------------------------*/
/* acf_bandwidth:  function: bandwidth
 * gain block, with a bandwidth. (ac only)
 * bad design: should be keyword instead.
 * arg0 = dc gain
 * arg1 = 3 db freq.
 */
/*ARGSUSED*/
static void acf_bandwidth(const branch_t *brh, double **arg, complex_t* y)
{
 if (!skip  &&  arg1 != 0.){
    double ratio;
    double coeff;
    ratio = ac.freq / arg1;
    coeff = arg0 / (1.+(ratio*ratio));
    y->x += coeff;
    y->y -= coeff * ratio;
 }
 *arg += aBANDWIDTH;
}
/*--------------------------------------------------------------------------*/
/* acf_complex:  function: complex
 * complex value, cartesian coordinates, in frequency domain (ac only)
 * arg0 = real part
 * arg1 = imaginary part
 */
/*ARGSUSED*/
static void acf_complex(const branch_t *brh, double **arg, complex_t* y)
{
 if (!skip){
    y->x += arg0;
    y->y += arg1;
 }
 *arg += aCOMPLEX;
}
/*--------------------------------------------------------------------------*/
/* acf_cornerdown:  function: cornerdown
 * piecewise linear function:
 * output is 0 for bias >= arg1
 * derivative is arg0 for bias <= arg1
 * arg0 = active slope
 * arg1 = break point
 */
static void acf_cornerdown(const branch_t *brh, double **arg, complex_t* y)
{
 if (!skip  &&  xx0 < arg1)
    y->x += arg0;
 *arg += aCORNERDOWN;
}
/*--------------------------------------------------------------------------*/
/* acf_cornerup:  function: cornerup
 * piecewise linear function:
 * output is 0 for bias <= arg1
 * derivative is arg0 for bias >= arg1
 * arg0 = active slope
 * arg1 = break point
 */
static void acf_cornerup(const branch_t *brh, double **arg, complex_t* y)
{
 if (!skip  &&  xx0 > arg1)
    y->x += arg0;
 *arg += aCORNERUP;
}
/*--------------------------------------------------------------------------*/
/* acf_delay:  function: delay
 * time delay (ac only)
 * arg0 = magnitude
 * arg1 = delay
 */
static void acf_delay(const branch_t *brh, double **arg, complex_t* y)
{
 if (!skip){
    double ratio;
    double dphase;
    ratio = ac.freq * arg1;
    if (ratio > 100000.){
       error(bPICKY, "delay too long: %s\n", printlabel(brh,NO));
       ratio = 0.;
    }
    dphase = -360.0 * ratio;
    y->x += real(arg0, dphase);
    y->y += imag(arg0, dphase);
 }
 *arg += aDELAY;
}
/*--------------------------------------------------------------------------*/
/* acf_exp:  spice compatible function: exp
 * spice source exponential function: exponential function of time
 * no-op in ac analysis, with a complaint
 * arg0 = initial value
 * arg1 = pulsed value
 * arg2 = rise delay
 * arg3 = rise time const
 * arg4 = fall delay
 * arg5 = fall time const
 */
/*ARGSUSED*/
static void acf_exp(const branch_t *brh, double **arg, complex_t* y)
{
 if (!skip)
    error(bDEBUG,"%s: exp not supported in ac analysis\n",printlabel(brh,NO));
 *arg += aEXP;
}
/*--------------------------------------------------------------------------*/
/* acf_expterm:  function: expterm
 * exponent term, non-integer power term (not exponential)
 * like polyterm, but for non-integers
 * only works for positive input, because negative number raised to
 * non-integer power is usually complex.
 * arg0 = coefficient
 * arg1 = exponent
 */
static void acf_expterm(const branch_t *brh, double **arg, complex_t* y)
{
 if (!skip  &&  xx0 > 0.){
    double coeff;
    coeff = arg0 * pow(xx0,arg1-1);
    y->x += coeff * arg1;
 }
 *arg += aEXPTERM;
}
/*--------------------------------------------------------------------------*/
/* acf_generator:  function: generator
 * value is derived from the "signal generator" (generator command)
 * intended for fixed sources, as circuit input.
 * arg0 = scale factor
 */
/*ARGSUSED*/
static void acf_generator(const branch_t *brh, double **arg, complex_t* y)
{
 if (!skip){
    y->x += arg0;
    /* y->y += 0.; */
 }
 *arg += aGENERATOR;
}
/*--------------------------------------------------------------------------*/
/* acf_max:  function: max
 * piecewise linear function: block that clips
 * bad design: should be keyword, so it can apply to the total part
 *	       should be able to specify upper and lower limits
 * arg0 = normal value
 * arg1 = output clip pt
 */
static void acf_max(const branch_t *brh, double **arg, complex_t* y)
{
 if (!skip  &&  arg0 != 0.){
    double clip_input;
    clip_input = arg1/arg0;
    if (xx0 >= -clip_input   &&   xx0 <= clip_input)
       y->x += arg0;
 }
 *arg += aMAX;
}
/*--------------------------------------------------------------------------*/
/* acf_netfunc:  function: netfunction
 * gain block, nth order response. (ac only)
 * BUG: only partial checking for number of args
 * arg0 = arg count
 * arg1 = dc gain
 * arg2 = ref frequency
 * arg3 = order of denom
 * arg* = coefs of denom
 * arg* = coefs of numer
 */
/*ARGSUSED*/
static void acf_netfunc(const branch_t *brh, double **arg, complex_t* y)
{
 int argcount;
 argcount = (int)arg0;

 if (!skip  &&  arg1 != 0.){
    double *denomcoeff;	/* coefficients of denominator			*/
    double *numercoeff;	/* coefficients of numerator			*/
    int denomorder;	/* order of denominator (#coefs = order + 1)	*/
    int numerorder;	/* order of numerator				*/
    double wacc;	/* stash for powers of f			*/
    complex_t denom;	/* denominator					*/
    complex_t numer;	/* numerator					*/
    complex_t result;	/* result of evaluating the function		*/
    int ii;		/* generic loop index				*/
    double normfreq;	/* normalized frequency				*/
    double nf2;		/*	"	"	squared			*/

    normfreq = ac.freq / ((arg2 != 0.) ? arg2 : 1./kPIx2);
    nf2 = normfreq * normfreq;

    denomorder = (int)arg3;
    numerorder = argcount - denomorder - 6;
    if (numerorder > denomorder){
       numerorder = denomorder;
       error(bWARNING, "%s: too many args\n", printlabel(brh,NO));
    }else if (numerorder < 0){
       denomorder += numerorder;
       numerorder = 0;
       error(bWARNING, "%s: too few args\n", printlabel(brh,NO));
    }

    denomcoeff = &(arg[0][4]);
    numercoeff = &(arg[0][5+denomorder]);

    wacc = 1.;
    denom.x = denomcoeff[0]; /* * wacc */
    for (ii = 2;   ii <= denomorder;   ii += 2){
       wacc *= -nf2;
       denom.x += denomcoeff[ii] * wacc;
    }

    wacc = normfreq;
    denom.y = denomcoeff[1] * wacc;
    for (ii = 3;   ii <= denomorder;   ii += 2){
       wacc *= -nf2;
       denom.y += denomcoeff[ii] * wacc;
    }

    wacc = 1.;
    numer.x = numercoeff[0]; /* * wacc */
    for (ii = 2;   ii <= numerorder;   ii += 2){
       wacc *= -nf2;
       numer.x += numercoeff[ii] * wacc;
    }

    wacc = normfreq;
    numer.y = numercoeff[1] * wacc;
    for (ii = 3;   ii <= numerorder;   ii += 2){
       wacc *= -nf2;
       numer.y += numercoeff[ii] * wacc;
    }

    result = cdiv(numer, denom);
    y->x += result.x * arg1;
    y->y += result.y * arg1;
 }
 *arg += argcount;
}
/*--------------------------------------------------------------------------*/
/* acf_notch:  function: notch
 * piecewise linear function: crossover notch, dead zone
 * symmetric around zero
 * arg0 = normal value
 * arg1 = dead zone size
 */
static void acf_notch(const branch_t *brh, double **arg, complex_t* y)
{
 if (!skip){
    if (xx0 > arg1   ||   xx0 < -arg1)
       y->x += arg0;
 }
 *arg += aNOTCH;
}
/*--------------------------------------------------------------------------*/
/* acf_numeric:  simple numeric argument
 * arg0 = value
 */
/*ARGSUSED*/
static void acf_numeric(const branch_t *brh, double **arg, complex_t* y)
{
 if (!skip)
    y->x += arg0;
 *arg += aNUMERIC;
}
/*--------------------------------------------------------------------------*/
/* acf_offset:  function: offset
 * fixed dc offset
 * bad design: should be keyword
 * arg0 = gain
 * arg1 = output offset
 */
/*ARGSUSED*/
static void acf_offset(const branch_t *brh, double **arg, complex_t* y)
{
 if (!skip)
    y->x += arg0;
 *arg += aOFFSET;
}
/*--------------------------------------------------------------------------*/
/* acf_polar:  function: polar
 * complex value in polar coordinates, in frequency domain (ac only)
 * arg0 = magnitude
 * arg1 = phase
 */
/*ARGSUSED*/
static void acf_polar(const branch_t *brh, double **arg, complex_t* y)
{
 if (!skip){
    y->x += real(arg0, arg1);
    y->y += imag(arg0, arg1);
 }
 *arg += aPOLAR;
}
/*--------------------------------------------------------------------------*/
/* acf_polyterm:  function: polyterm
 * polynomial term, one term in a polynomial
 * power must be integer, is rounded to nearest integer
 * caution: will divide by zero with zero input and negative exponent
 * arg0 = coefficient
 * arg1 = exponent
 */
static void acf_polyterm(const branch_t *brh, double **arg, complex_t* y)
{
 if (!skip){
    int expo;
    expo = (int)floor(arg1+.5);
    if (expo != 0){
       double coeff;
       coeff = arg0 * ipow(xx0,expo-1);
       y->x += coeff * expo;
    }
 }
 *arg += aPOLYTERM;
}
/*--------------------------------------------------------------------------*/
/* acf_pulse: spice compatible function: pulse
 * spice pulse function (for sources)
 * no-op in ac analysis, with a complaint
 * arg0 = initial value
 * arg1 = pulsed value
 * arg2 = delay time
 * arg3 = rise time
 * arg4 = fall time
 * arg5 = pulse width
 * arg6 = period
 */
/*ARGSUSED*/
static void acf_pulse(const branch_t *brh, double **arg, complex_t* y)
{
 if (!skip)
   error(bDEBUG,"%s: pulse not supported in ac analysis\n",printlabel(brh,NO));
 *arg += aPULSE;
}
/*--------------------------------------------------------------------------*/
/* acf_pwl:  spice compatible function: pwl
 * "piece-wise linear" point-plotting function of time.
 * no-op in ac analysis, with a complaint
 * arg0 = arg count
 * no other args matter
 */
/*ARGSUSED*/
static void acf_pwl(const branch_t *brh, double **arg, complex_t* y)
{
 if (!skip)
    error(bDEBUG,"%s: pwl not supported in ac analysis\n",printlabel(brh,NO));
 *arg += (int)arg0;
}
/*--------------------------------------------------------------------------*/
/* acf_sffm:  spice compatible function: sffm
 * single frequency frequency modulation
 * no-op in ac analysis, with a complaint
 * arg0 = dc offset
 * arg1 = amplitude
 * arg2 = carrier freq
 * arg3 = mod index
 * arg4 = signal freq
 */
/*ARGSUSED*/
static void acf_sffm(const branch_t *brh, double **arg, complex_t* y)
{
 if (!skip)
    error(bDEBUG,"%s: sffm not supported in ac analysis\n",printlabel(brh,NO));
 *arg += aSFFM;
}
/*--------------------------------------------------------------------------*/
/* acf_sin:  spice compatible function: sin
 * sinusoidal function, mainly for sources
 * no-op in ac analysis, with a complaint
 * arg0 = offset
 * arg1 = amplitude
 * arg2 = frequency
 * arg3 = delay
 * arg4 = damping
 */
/*ARGSUSED*/
static void acf_sin(const branch_t *brh, double **arg, complex_t* y)
{
 if (!skip)
    error(bDEBUG,"%s: sin not supported in ac analysis\n",printlabel(brh,NO));
 *arg += aSIN;
}
/*--------------------------------------------------------------------------*/
/* acf_tanh:  function: tanh
 * for now, a copy of max
 * piecewise linear function: block that clips
 * bad design: should be keyword, so it can apply to the total part
 *	       should be able to specify upper and lower limits
 * arg0 = normal value
 * arg1 = output clip pt
 */
static void acf_tanh(const branch_t *brh, double **arg, complex_t* y)
{
 if (!skip){
    if (arg1 != 0.){
       double cosine;
       cosine = cosh(xx0 * arg0/arg1);
       y->x += arg0 / (cosine*cosine);
    }
    /* else 0 */
 }
 *arg += aTANH;
}
/*--------------------------------------------------------------------------*/
/*ARGSUSED*/
static void acf_ic(const branch_t *brh, double **arg, complex_t* y)
{
 *arg += 1;
}
/*--------------------------------------------------------------------------*/
/*ARGSUSED*/
static void acf_ii(const branch_t *brh, double **arg, complex_t* y)
{
 *arg += 1;
}
/*--------------------------------------------------------------------------*/
/*ARGSUSED*/
static void acf_iv(const branch_t *brh, double **arg, complex_t* y)
{
 *arg += 1;
}
/*--------------------------------------------------------------------------*/
/*ARGSUSED*/
static void acf_tempco(const branch_t *brh, double **arg, complex_t* y)
{
#ifdef NEVER
 tempco = arg0;
#endif
 *arg += 1;
}
/*--------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------*/
