/* d_coil.c  94.12.09
 * Copyright 1983-1994   Albert Davis
 * inductors
 * two levels: linear (lin) and nonlinear (nl) (not really)
 * x = amps, y.f0 = flux, ev = y.f1 = henrys
 */
#include "ecah.h"
#include "ac.h"
#include "branch.h"
#include "d_subckt.h"
#include "error.h"
#include "mode.h"
#include "options.h"
#include "tr.h"
#include "types.h"
#include "declare.h"
/*--------------------------------------------------------------------------*/
static	void	parse_coil_mutual(branch_t*,const char*,int*);
static	void	print_coil_mutual(const branch_t*,int,int);
static	void	expand_coil_mutual(branch_t*);
static	void	expand_coil(branch_t*);
static 	int	tr_coil_lin(branch_t*);
static 	int	tr_coil_nl(branch_t*);
static	void	integrate(branch_t*);
static 	void	ac_coil_lin(branch_t*);
static 	void	ac_coil_nl(branch_t*);
static	double	tr_review_coil(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 int in_curr;		/* flag: initial current		    */
extern int in_volt;		/* flag: initial voltage		    */
extern char e_int[];
/*--------------------------------------------------------------------------*/
functions_t dev_coil_mutual = {
   (generic_t*)NULL,	/* x */
   sizeof(functions_t),	/* ssize */
   sizeof(branch_t),	/* elementsize */
   (functions_t*)NULL,	/* super */
   0, 			/* numnodes */
   rnTWOPORT,		/* refnode */
   rnTWOPORT,		/* isdevice */
   create_std,		/* create */
   copy_std,		/* copy */
   parse_coil_mutual,	/* parse */
   print_coil_mutual,	/* print */
   expand_coil_mutual,	/* expand */
   NULL,		/* probe */
   NULL,		/* tr_probe */
   NULL,		/* ac_probe */
   NULL,		/* xprobe */
   NULL,		/* dotr */
   NULL,		/* untr */
   NULL,		/* doac */
   NULL,		/* trfun1 */
   NULL,		/* trfun0 */
   NULL,		/* acfun */
   NULL,		/* tr_guess */
   NULL,		/* tr_advance */
   NULL			/* tr_review */
};
functions_t dev_coil = {
   (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_coil,		/* expand */
   probe_std,		/* probe */
   tr_probe_std,	/* tr_probe */
   ac_probe_std,	/* ac_probe */
   xprobe_std,		/* xprobe */
   tr_coil_nl,		/* dotr */
   unloadpassive,	/* untr */
   ac_coil_nl,		/* doac */
   trfix1,		/* trfun1 */
   trfix0,		/* trfun0 */
   acfix,		/* acfun */
   NULL,		/* tr_guess */
   NULL,		/* tr_advance */
   tr_review_coil	/* tr_review */
};
/*--------------------------------------------------------------------------*/
static void parse_coil_mutual(branch_t *brh, const char *cmd, int *cnt)
{
 parselabel(brh,cmd,cnt);
 parselabel_raw(brh->outputlabel,cmd,cnt);
 parselabel_raw(brh->inputlabel,cmd,cnt);
 parseexpr(brh,cmd,cnt);
 if (brh->f->super)
    brh->f = brh->f->super;
}
/*--------------------------------------------------------------------------*/
static void print_coil_mutual(const branch_t *brh, int where, int detail)
{
 (void)printlabel(brh,where);
 if (brh->output){
    printlabel(brh->output,where);
 }else{
    mprintf(where, " %s  ", brh->outputlabel);
 }
 if (brh->input){
    printlabel(brh->input,where);
 }else{
    mprintf(where, " %s  ", brh->inputlabel);
 }
 printexpr(brh,where);
 mputc('\n', where);
}
/*--------------------------------------------------------------------------*/
static void expand_coil_mutual(branch_t *brh)
{
 double l1;
 double l2;
 double lm;
 double det;
 static functions_t *dev_coil_coupled;
 branch_t *pri, *sec, *sub;

 if (!dev_coil_coupled){
    dev_coil_coupled = (functions_t*)calloc(1, sizeof(functions_t));
    *dev_coil_coupled = dev_coil;
    dev_coil_coupled->super = &dev_coil;
    dev_coil_coupled->expand = NULL;
    dev_coil_coupled->dotr = tr_subckt;
    dev_coil_coupled->untr = un_subckt;
    dev_coil_coupled->doac = ac_subckt;
    dev_coil_coupled->trfun1 = NULL;
    dev_coil_coupled->trfun0 = NULL;
    dev_coil_coupled->acfun  = NULL;
    dev_coil_coupled->tr_review  = tr_review_subckt;
 }
 
 brh->output = findbranch_samescope(brh->outputlabel,brh);
 if (!brh->output)
    error(bERROR,"%s: can't find %s\n",printlabel(brh,NO),brh->outputlabel);

 brh->input = findbranch_samescope(brh->inputlabel,brh);
 if (!brh->input)
    error(bERROR,"%s: can't find %s\n",printlabel(brh,NO),brh->inputlabel);

 l1 = brh->output->val;
 l2 = brh->input->val;
 lm  = brh->val * sqrt(l1 * l2);
 det = l1 * l2 - lm * lm;

 
 pri = copy_branch(brh->output);
 pri->parent  = brh->output;
 pri->ev.x  = pri->y0.f1 = pri->val = det / l2;
 pri->y0.f0 = LINEAR;
 pri->ev.y  = 0.;
 pri->next = brh->output->subckt;
 pri->prev = (branch_t*)NULL;
 brh->output->subckt = insertbranch(pri);

 sec = copy_branch(brh->input);
 sec->parent  = brh->input;
 sec->ev.x  = sec->y0.f1 = sec->val = det / l1;
 sec->y0.f0 = LINEAR;
 sec->ev.y  = 0.;
 sec->next = brh->input->subckt;
 sec->prev = (branch_t*)NULL;
 brh->input->subckt = insertbranch(sec);

 sub = create_branch(&dev_cccs);
 sub->parent  = brh->output;
 sub->input   = sec;
 sprintf(sub->label, "F%s", &sub->input->label[1]);
 sub->n[OUT1] = brh->output->n[OUT1];
 sub->n[OUT2] = brh->output->n[OUT2];
 sub->ev.x  = sub->y0.f1 = sub->val = -lm / l1;
 sub->y0.f0 = LINEAR;
 sub->ev.y  = 0.;
 sub->next = brh->output->subckt;
 sub->prev = (branch_t*)NULL;
 brh->output->subckt = insertbranch(sub);

 sub = create_branch(&dev_cccs);
 sub->parent  = brh->input;
 sub->input   = pri;
 sprintf(sub->label, "F%s", &sub->input->label[1]);
 sub->n[OUT1] = brh->input->n[OUT1];
 sub->n[OUT2] = brh->input->n[OUT2];
 sub->ev.x  = sub->y0.f1 = sub->val = -lm / l2;
 sub->y0.f0 = LINEAR;
 sub->ev.y  = 0.;
 sub->next = brh->input->subckt;
 sub->prev = (branch_t*)NULL;
 brh->input->subckt = insertbranch(sub);

 brh->output->f = brh->input->f = dev_coil_coupled;
}
/*--------------------------------------------------------------------------*/
static void expand_coil(branch_t *brh)
{
 if (!brh->x){
    static functions_t *ff;
    if (!ff){
       ff = (functions_t*)calloc(1, sizeof(functions_t));
       *ff = dev_coil;
       ff->super = &dev_coil;
       ff->dotr = tr_coil_lin;
       ff->doac = ac_coil_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 int tr_coil_lin(branch_t *brh)
{
 brh->m0.x = tr_volts(&(brh->n[OUT1]),&(brh->n[OUT2]));
 brh->y0.x = brh->m0.c0 + brh->m0.f1 * brh->m0.x;
 /* brh->y0.f1 = brh->val; set in expand */
 brh->y0.f0 = brh->y0.x * brh->y0.f1;
 if (brh->y0.f1 == 0.){
    error(bPICKY, "%s: short circuit\n", printlabel(brh,NO));
 }
 integrate(brh);
 trloadpassive(brh);
 return brh->converged = YES;
}
/*--------------------------------------------------------------------------*/
static int tr_coil_nl(branch_t *brh)
{
 in_volt = in_curr = NO;
 brh->m0.x = tr_volts_limited(&(brh->n[OUT1]),&(brh->n[OUT2]));
 brh->y0.x = brh->m0.c0 + brh->m0.f1 * brh->m0.x;
 if (brh->f->trfun1){
    (*brh->f->trfun1)(brh);
 }else{
    brh->y0.f1 = brh->val;
    brh->y0.f0 = brh->y0.x * brh->y0.f1;
 }
 if (brh->y0.f1 == 0.){
    error(bPICKY, "%s: short circuit\n", printlabel(brh,NO));
 }

 integrate(brh);
 trloadpassive(brh);
 return brh->converged = conv_check(brh);
}
/*--------------------------------------------------------------------------*/
static void integrate(branch_t *brh)
{
 if (sim_mode == sDC  ||  sim_mode == sOP  ||  brh->y0.f1 == 0.){
    brh->m0.f1 = 1./opt.shortckt;
    brh->m0.c0 = 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.f1 = 1./opt.shortckt;
	  brh->m0.c0 = 0.;
       }
    }else{
       /* leave it alone to restore */
    }
 }else{
    double dt = brh->time0 - brh->time1;
    if (brh->method_a == mGEAR){
       brh->m0.f1 = dt / brh->y0.f1;
       brh->m0.c0 = brh->mt1.c0 + brh->mt1.f1 * brh->mt1.x; /* oldi */
    }else if (brh->method_a == mTRAPEZOID){
       /* oldi = brh->mt1.c0 + brh->mt1.f1 * brh->mt1.x; */
       brh->m0.f1 = dt/(2*brh->y0.f1);
       brh->m0.c0 = brh->mt1.c0 + (brh->mt1.f1 + brh->m0.f1) * brh->mt1.x;
       		    /* oldi + f1*oldv, with combined terms */
    }else{
       error(bDANGER, 
       		"internal error: inductance: bad integration method %u\n",
       		brh->method_a);
    }
 }
}
/*--------------------------------------------------------------------------*/
static void ac_coil_lin(branch_t *brh)
{
 if (ac.omega == 0.){
    brh->acg.x = 1. / opt.shortckt;
    brh->acg.y = 0.;
    acloadpassivereal(brh);
    error(bPICKY, "%s: short circuit\n", printlabel(brh,NO));
 }else{
    brh->acg.x = 0.;
    brh->acg.y = -1. / (brh->ev.x * ac.omega);
    acloadpassiveimaginary(brh);
 }
}
/*--------------------------------------------------------------------------*/
static void ac_coil_nl(branch_t *brh)
{
 double dcvolts;

 if (brh->f->acfun){
    dcvolts = dc_volts(&(brh->n[OUT1]),&(brh->n[OUT2]));
    brh->acbias = brh->m0.c0 + brh->m0.f1*dcvolts;
    brh->ev = (*brh->f->acfun)(brh);
    if (ac.omega == 0.){
       brh->acg.x = 1. / opt.shortckt;
       brh->acg.y = 0.;
       error(bPICKY, "%s: short circuit\n", printlabel(brh,NO));
    }else{
       complex_t y;
       y.x = -brh->ev.y * ac.omega;
       y.y =  brh->ev.x * ac.omega;
       brh->acg = cflip(y);
    }
    acloadpassive(brh);
 }else{
    brh->ev.x  = brh->y0.f1;
    brh->ev.y  = 0.;
    if (ac.omega == 0.){
       brh->acg.x = 1. / opt.shortckt;
       brh->acg.y = 0.;
       acloadpassivereal(brh);
       error(bPICKY, "%s: short circuit\n", printlabel(brh,NO));
    }else{
       brh->acg.x = 0.;
       brh->acg.y = -1. / (brh->ev.x * ac.omega);
       acloadpassiveimaginary(brh);
    }
 }
}
/*--------------------------------------------------------------------------*/
static double tr_review_coil(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.x;					/* volts */
    double it1 = brh->mt1.x;
    /*double it2 = brh->mt2.x;*/
    /*double didt0 = (i0  - it1) / dt0;*/			/* 1st der */
    /*double didt1 = (it1 - it2) / dt1;*/
    /*double ddiddt = (didt0 - didt1) / ddt0;*/			/* 2nd der */

    double q0  = brh->y0.f0;					/* flux */
    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;
    }
 }
}
/*--------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------*/
