module: dco_rlc_with_neg_conduct
parameters: double ind_val, double r_loss
inputs: double cur_in, double cap_val, double_interp cap_clk
outputs: double cap_vol, double ind_cur, double_interp squareout
         double inst_freq
classes: Filter int_c("1","1"), Filter int_l("1","1"),
         EdgeMeasure edge_sig() 
static_variables: double C_charge_prev, double L_cur_prev
                  double w_res_nom, double L_cur_scale, int num_iter
                  double inv_LC, double R_over_L, double cond_over_C
                  double cap_val_clked, int clk_edge_flag
                  double cap_val_clked_next, int warning_flag
                  int freq_edge_flag, double inst_freq_next
init:
double Ts_scaled;

if (r_loss < 0.0)
   {
    printf("error in 'dco_rlc_with_neg_conduct':\n");
    printf("  parameter 'r_loss' must be >= 0.0\n");
    printf("  in this case, r_loss = %5.3e\n",r_loss);
    exit(1);
   }
if (ind_val < 1.0e-22)
   {
    printf("error in 'dco_rlc_with_neg_conduct':\n");
    printf("  parameter 'ind_val' must be >= 1e-22\n");
    printf("  in this case, ind_val = %5.3e\n",ind_val);
    exit(1);
   }

num_iter = 100;
Ts_scaled = Ts/((double) num_iter); 

// create integrators to implement state-space model
int_c.set("1","s","Ts",Ts_scaled);
int_l.set("1","s","Ts",Ts_scaled);
C_charge_prev = 0.0;
L_cur_prev = 0.0;

R_over_L = r_loss/ind_val;
cond_over_C = R_over_L*2.0;
squareout = -1.0;
clk_edge_flag = 0;
ind_cur = 0.0;
cap_vol = 0.0;
warning_flag = 0;
freq_edge_flag = 0;

code:
int i, i_clk;
double L_cur, C_charge, cap_val_clked_prev, clk_interp, freq_interp;

if (4.0*pi*pi*cap_val*ind_val < Ts*Ts/16.0)
   {
    if (warning_flag == 0)
      {
       printf("Warning in 'dco_rlc_with_neg_conduct':\n");
       printf("  resonant frequency must < 1/4 the simulation\n");
       printf("  sample rate\n");
       printf("  -> in this case, L = %5.3e, C = %5.3e\n", ind_val, cap_val);
      }
    return;
   }
if (clk_edge_flag == 0)
   {
    cap_val_clked = cap_val;
    inv_LC = 1.0/(ind_val*cap_val_clked);
   }
i_clk = clk_iter_num(cap_clk,num_iter,clk_interp);
if (i_clk >= 0)
   {
    cap_val_clked_prev = cap_val_clked;
    cap_val_clked_next = cap_val;
   }

for (i = 0; i < num_iter; i++)
   { 
   if (i == i_clk)
      {
       clk_edge_flag = 1;
       cap_val_clked = cap_val_clked_prev*clk_interp + cap_val_clked_next*(1.0-clk_interp);
       inv_LC = 1.0/(ind_val*cap_val_clked);
      }
   else if (clk_edge_flag == 1)
      { 
       clk_edge_flag = 2;
       cap_val_clked = cap_val_clked_next;
       inv_LC = 1.0/(ind_val*cap_val_clked);
      }

   // LC VCO state-space model
   L_cur = int_l.inp(inv_LC*C_charge_prev - R_over_L*L_cur_prev);
   C_charge = int_c.inp(cond_over_C*C_charge_prev - L_cur_prev + cur_in);

   squareout = square_interp(C_charge,C_charge_prev,i,num_iter);
   cap_vol = C_charge/cap_val_clked;
   // negative resistance of VCO
   //cond_over_C = neg_conduct/cap_val;
   cond_over_C = R_over_L*2.0/(1.0+fabs(1.3*cap_vol));

   L_cur_prev = L_cur;
   C_charge_prev = C_charge;
   }

edge_sig.inp(squareout);
if (edge_sig.out >= 4.0) // assume freq should be no more than 1/(4Ts)
    {
    freq_edge_flag = 1;
    inst_freq_next = 1.0/(Ts*edge_sig.out);
    freq_interp = (squareout+1.0)/2.0*inst_freq_next + 
                  (1.0-squareout)/2.0*inst_freq;
    inst_freq = freq_interp;
    }
else if (freq_edge_flag == 1)
    {
    freq_edge_flag = 0;
    inst_freq = inst_freq_next;
    }

ind_cur = L_cur;


functions:
   
int clk_iter_num(double clk, int num_iter, double &residue)
{
static double clk_prev = 0.0;
double num_iter_dbl,out_dbl;

if (clk_prev == -1.0 && clk != -1.0)
   {
   num_iter_dbl = (double) num_iter;
   out_dbl = (1.0-clk)/2.0*num_iter_dbl;
   }
else
   out_dbl = -1.0;

clk_prev = clk;

residue = out_dbl - floor(out_dbl);
return((int) floor(out_dbl));
}

double square_interp(double in, double in_prev, int i, int num_iter)
{
static int state = -2;
static double out = -1.0;
double i_dbl, num_iter_dbl;


if (state == -2 && in >= 0.0) /* positive transition */
   {
    if (fabs(in-in_prev) > 1e-28)
      {
       i_dbl = (double) i;
       num_iter_dbl = (double) num_iter;
       out = (num_iter_dbl - 2.0*i_dbl - 1.0 + 
             (in + in_prev)/(in - in_prev))/num_iter_dbl;
      }
    else
       out = 1.0;
    state = 1;
   }
else if (state == 1 && i == 0)
   {
    out = 1.0;
    state = 2;
   }
else if (state == 2 && in < 0.0) /* negative transition */
   {
    if (fabs(in-in_prev) > 1e-28)
      {
       i_dbl = (double) i;
       num_iter_dbl = (double) num_iter;
       out = -(num_iter_dbl - 2.0*i_dbl - 1.0 + 
              (in + in_prev)/(in - in_prev))/num_iter_dbl;
      }
    else
       out = -1.0;
    state = -1;
   }
else if (state == -1 && i == 0)
   {
    out = -1.0;
    state = -2;
   }

return(out);
}
