/* d_cap.c  94.12.09
 * Copyright 1983-1992   Albert Davis
 * capacitor models
 * two levels: linear (lin) and nonlinear (nl)
 * y.x = volts, y.f0 = coulombs, ev = y.f1 = farads
 */
#include "ecah.h"
#include "ac.h"
#include "branch.h"
#include "error.h"
#include "mode.h"
#include "options.h"
#include "tr.h"
#include "declare.h"
/*--------------------------------------------------------------------------*/
static	void	expand_cap(branch_t*);
static	double	probe_cap(const branch_t*,const char*);
static 	int	tr_cap_lin(branch_t*);
static 	int	tr_cap_nl(branch_t*);
static	void	integrate(branch_t*);
static 	void	ac_cap_lin(branch_t*);
static 	void	ac_cap_nl(branch_t*);
static	double	tr_review_cap(branch_t*);
/*--------------------------------------------------------------------------*/
extern const double init_curr;	/* initial current			    */
extern const double init_volt;	/* initial voltage			    */
extern const double trtime0;	/* transient analysis time		    */
extern const int sim_mode;	/* what type simulation is this		    */
extern const int sim_phase;	/* how far it got			    */
extern const ac_t ac;
extern const transient_t *tr;
extern const struct options opt;
extern const char e_int[];
extern int in_curr;		/* flag: initial current		    */
extern int in_volt;		/* flag: initial voltage		    */
/*--------------------------------------------------------------------------*/
functions_t dev_cap = {
   (generic_t*)NULL,	/* x */
   sizeof(functions_t),	/* ssize */
   sizeof(branch_t),	/* elementsize */
   (functions_t*)NULL,	/* super */
   2, 			/* numnodes */
   rnONEPORT,		/* refnode */
   rnONEPORT,		/* isdevice */
   create_std,		/* create */
   copy_std,		/* copy */
   parse_std,		/* parse */
   print_std,		/* print */
   expand_cap,		/* expand */
   probe_cap,		/* probe */
   tr_probe_std,	/* tr_probe */
   ac_probe_std,	/* ac_probe */
   xprobe_std,		/* xprobe */
   tr_cap_nl,		/* dotr */
   unloadpassive,	/* untr */
   ac_cap_nl,		/* doac */
   trfix1,		/* trfun1 */
   trfix0,		/* trfun0 */
   acfix,		/* acfun */
   NULL,		/* tr_guess */
   NULL,		/* tr_advance */
   tr_review_cap	/* tr_review */
};
/*--------------------------------------------------------------------------*/
static void expand_cap(branch_t *brh)
{
 if (!brh->x){
    static functions_t *ff;
    if (!ff){
       ff = (functions_t*)calloc(1, sizeof(functions_t));
       *ff = dev_cap;
       ff->super = &dev_cap;
       ff->refnode = rnLINEAR;
       ff->dotr = tr_cap_lin;
       ff->doac = ac_cap_lin;
       ff->trfun1 = NULL;
       ff->trfun0 = NULL;
       ff->acfun  = NULL;
    }
    brh->f = ff;
    brh->y0.f1 = brh->val;
    brh->y0.f0 = LINEAR;
    brh->ev.x  = brh->val;
    brh->ev.y  = 0.;
 }
}
/*--------------------------------------------------------------------------*/
static double probe_cap(const branch_t *brh, const char *what)
{
 int dummy = 0;
 setmatch(what,&dummy);
 if (rematch("Q")){
    return probe_std(brh,"F");
 }else{
    return probe_std(brh,what);
 }
}
/*--------------------------------------------------------------------------*/
static int tr_cap_lin(branch_t *brh)
{
 brh->m0.x = brh->y0.x = tr_volts(&(brh->n[OUT1]),&(brh->n[OUT2]));
 /* brh->y0.f1 = brh->val; */
 brh->y0.f0 = brh->y0.x * brh->y0.f1;
 integrate(brh);
 trloadpassive(brh);
 return brh->converged = YES;
}
/*--------------------------------------------------------------------------*/
static int tr_cap_nl(branch_t *brh)
{
 in_volt = in_curr = NO;
 brh->m0.x = brh->y0.x = tr_volts_limited(&(brh->n[OUT1]),&(brh->n[OUT2]));
 brh->y0.f0 = brh->m0.f1 * brh->m0.x + brh->m0.c0;
 if (brh->f->trfun1){
    (*brh->f->trfun1)(brh);
 }else{
    brh->y0.f1 = brh->val;
    brh->y0.f0 = brh->y0.x * brh->y0.f1;    
 }
 integrate(brh);
 trloadpassive(brh);
 return brh->converged = conv_check(brh);
}
/*--------------------------------------------------------------------------*/
static void integrate(branch_t *brh)
{
 if (sim_mode == sDC  ||  sim_mode == sOP){
    brh->m0.c0 = brh->m0.f1 = 0.;
 }else if (sim_phase == pINIT_DC){
    if (trtime0 == 0.){
       if (in_curr){
	  brh->m0.f1 = 0.;
	  brh->m0.c0 = -init_curr;
       }else if (in_volt){
	  brh->m0.f1 = 1./opt.shortckt;
	  brh->m0.c0 = -init_volt * brh->m0.f1;
       }else{
	  brh->m0.c0 = brh->m0.f1 = 0.;
       }
    }else{
       /* leave it alone to restore */
    }
 }else if (brh->time0 == brh->time1){
    error(bERROR, e_int, "time step is 0.\n");
 }else{
    double oldv = brh->mt1.x;
    double dt = brh->time0 - brh->time1;
    if (brh->method_a == mGEAR){
       brh->m0.f1 = brh->y0.f1 / dt;		/* first order 		    */
       brh->m0.c0 = -brh->m0.f1 * oldv;		/* (stiff) (backward Euler) */
    }else if (brh->method_a == mTRAPEZOID){
       double oldi;				/* second order 	    */
       oldi = brh->mt1.c0 + brh->mt1.f1*oldv;	/* (non-stiff) (trapezoid)  */
       brh->m0.f1 = 2. * brh->y0.f1 / dt;
       brh->m0.c0 = -brh->m0.f1 * oldv - oldi;
    }else{
       error(bDANGER, 
       		"internal error: capacitance: bad integration method %u\n",
       		brh->method_a);
    }
 }
}
/*--------------------------------------------------------------------------*/
static void ac_cap_lin(branch_t *brh)
{
 brh->acg.x = 0.; /* -brh->ev.y * ac.omega  where ev.y == 0. */
 brh->acg.y = brh->ev.x * ac.omega;
 acloadpassiveimaginary(brh);
}
/*--------------------------------------------------------------------------*/
static void ac_cap_nl(branch_t *brh)
{
 if (brh->f->acfun){
    brh->acbias = dc_volts(&(brh->n[OUT1]),&(brh->n[OUT2]));
    brh->ev = (*brh->f->acfun)(brh);
    brh->acg.x = -brh->ev.y * ac.omega;
    brh->acg.y =  brh->ev.x * ac.omega;
    acloadpassive(brh);
 }else{
    brh->ev.x  = brh->y0.f1;
    brh->ev.y  = 0.;
    brh->acg.x = 0.; /* -brh->ev.y * ac.omega  where ev.y == 0. */
    brh->acg.y = brh->ev.x * ac.omega;
    acloadpassiveimaginary(brh);
 }
}
/*--------------------------------------------------------------------------*/
static double tr_review_cap(branch_t *brh)
{
 if (brh->time3 <= 0.){
    return BIGBIG;
 }else{
    double factor = 1./12.;      /* coefficient of 3rd der, trapezoid rule */
    double dt0 = brh->time0 - brh->time1;
    double dt1 = brh->time1 - brh->time2;	/* BUG: these values should */
    double dt2 = brh->time2 - brh->time3;	/* be stored */
    double ddt0 = brh->time0 - brh->time2;
    double ddt1 = brh->time1 - brh->time3;
    double dddt0 = brh->time0 - brh->time3;

    double i0  = brh->m0.f1  * brh->m0.x  + brh->m0.c0;		/* current */
    double it1 = brh->mt1.f1 * brh->mt1.x + brh->mt1.c0;
    /*double it2 = brh->mt2.f1 * brh->mt2.x + brh->mt2.c0;*/
    /*double didt0 = (i0  - it1) / dt0;*/			/* 1st der */
    /*double didt1 = (it1 - it2) / dt1;*/
    /*double ddiddt = (didt0 - didt1) / ddt0;*/			/* 2nd der */

    double q0  = brh->y0.f0;					/* charge */
    double qt1 = brh->yt1.f0;
    double qt2 = brh->yt2.f0;
    double qt3 = brh->yt3.f0;
    double dqdt0 = (q0   - qt1) / dt0;				/* 1st der */
    double dqdt1 = (qt1  - qt2) / dt1;
    double dqdt2 = (qt2  - qt3) / dt2;
    double ddqddt0 = (dqdt0 - dqdt1) / ddt0;			/* 2nd der */
    double ddqddt1 = (dqdt1 - dqdt2) / ddt1;
    double dddqdddt = (ddqddt0 - ddqddt1) / dddt0;		/* 3rd der */

    double currenttol = opt.abstol + opt.reltol * MAX(fabs(i0), fabs(it1));
    double chargetol = MAX(opt.chgtol,MAX(fabs(q0),fabs(qt1)))
			   * opt.reltol / dt0;
    double tol = MAX(currenttol,chargetol);
    double denom = MAX(opt.abstol, (factor *fabs(dddqdddt)));  /* avoid / 0 */
    double timestep = opt.trtol * sqrt(tol / denom);

    if (timestep <= tr->dtmin){
       error(bDANGER,"step control error:%s %g\n",printlabel(brh,0),timestep);
       error(bTRACE, "q0=%g i0=%g dq0=%g\n", q0, i0, dqdt0);
       error(bTRACE, "it=%g qt=%g tol=%g\n", currenttol, chargetol, tol);
       timestep = tr->dtmin;
    }
    if (timestep < dt0 * opt.trreject){
       error(bTRACE, "step rejected:%s\n", printlabel(brh,0));
       error(bTRACE, "new=%g  old=%g  rej=%g\n",
       		timestep, dt0, dt0 * opt.trreject);
       return brh->time1 + timestep;
    }else{
       return brh->time0 + timestep;
    }
 }
}
/*--------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------*/
