module: nrz_dac
parameters: 
  double per_mm
  double per_mm_sw
  int n_elements
  int dwa_en
  double i_nom
  int mm_seed
  double tz
  double tp1 
  double tp2
  double rout_nom
inputs:
  double in
  double vg_in
  double_interp clk
outputs:
  double dac_out
static_variables:
  int code_in
  int j
  int m
  int dwa_ptr
  int dwa_ptr_neg
  double dac_outm
  double dac_outp
  double dac_sum
  double rout_val
  double rout_val_min
  double ipar
  Filter *filter_chain_up
  Filter *filter_chain_dp
  Filter *filter_chain_um
  Filter *filter_chain_dm
classes:
  Rand rand_mm_dac("gauss")
  Rand rand_mm_filt("gauss")
  EdgeDetect clk_edge()
  Vector dac_vec()
  Vector dac_mm_vec()
  Vector dac_up_vec()
  Vector dac_dp_vec()
  Vector dac_um_vec()
  Vector dac_dm_vec()
  Vector up_en_vec()
  Vector dp_en_vec()
  Vector um_en_vec()
  Vector dm_en_vec()
  Vector gout_vec()
init:

  //zero output variables
  dac_outp = 0.0;
  dac_outm = 0.0;
  dac_sum = 0.0;
  rout_val_min = 0.01*rout_nom;

  //establish mismatch currents
  dac_vec.set_length(n_elements);
  dac_mm_vec.set_length(n_elements);
  dac_up_vec.set_length(n_elements);
  dac_dp_vec.set_length(n_elements);
  dac_um_vec.set_length(n_elements);
  dac_dm_vec.set_length(n_elements);
  up_en_vec.set_length(n_elements);
  dp_en_vec.set_length(n_elements);
  um_en_vec.set_length(n_elements);
  dm_en_vec.set_length(n_elements);

  //establish mismatch output conductances
  gout_vec.set_length(n_elements);

  //create the array of filters
  filter_chain_up = new Filter[n_elements];
  filter_chain_um = new Filter[n_elements];
  filter_chain_dp = new Filter[n_elements];
  filter_chain_dm = new Filter[n_elements];

  //seed the random number generator
  if (mm_seed > 0)
  {
     rand_mm_dac.set_seed(mm_seed);
     rand_mm_filt.set_seed(mm_seed);
  }

  //fill in the vectors, initialize the filter chain
  int i,k;
  for (k = 0; k < n_elements; k++)
  {
     dac_mm_vec.set_elem(k,1.0+1.0*per_mm/100*rand_mm_dac.inp());
     //check if obtain resistances smaller than the minimum likely
     //also, scale conductance by i_nom since all currents are normalized
     //to 1.0 A, while actual DAC current may be i_nom = 10.0 mA 
     rout_val = rout_nom+rout_nom*per_mm/100*rand_mm_dac.inp();
     if (rout_val < rout_val_min)
     {
        gout_vec.set_elem(k,1/rout_val_min/i_nom);
     }
     else
     {
        gout_vec.set_elem(k,1/rout_val/i_nom);
     }
     dac_sum += dac_mm_vec.get_elem(k);
  }  

  double eps1;
  double eps2;
  double eps3;
  double eps4;
  for (i = 0; i < n_elements; i++)
  {
    dac_vec.set_elem(i,dac_mm_vec.get_elem(i)/(dac_sum/n_elements));
    dac_up_vec.set_elem(i,0);
    dac_dp_vec.set_elem(i,0);
    dac_um_vec.set_elem(i,0);
    dac_dm_vec.set_elem(i,0);

    up_en_vec.set_elem(i,0);
    dp_en_vec.set_elem(i,0);
    um_en_vec.set_elem(i,0);
    dm_en_vec.set_elem(i,0);

    //define filter mismatch error epsilon
    eps1 = 1+per_mm_sw/100*rand_mm_filt.inp();
    eps2 = 1+per_mm_sw/100*rand_mm_filt.inp();
    eps3 = 1+per_mm_sw/100*rand_mm_filt.inp();
    eps4 = 1+per_mm_sw/100*rand_mm_filt.inp();
    filter_chain_up[i].set("1+tz*s","1+(tp1+tp2)*s+tp1*tp2*s^2","tz,tp1,tp2,Ts",eps1*tz,eps1*tp1,eps1*tp2,Ts);
    filter_chain_dp[i].set("1+tz*s","1+(tp1+tp2)*s+tp1*tp2*s^2","tz,tp1,tp2,Ts",eps2*tz,eps2*tp1,eps2*tp2,Ts);
    filter_chain_um[i].set("1+tz*s","1+(tp1+tp2)*s+tp1*tp2*s^2","tz,tp1,tp2,Ts",eps3*tz,eps3*tp1,eps3*tp2,Ts);
    filter_chain_dm[i].set("1+tz*s","1+(tp1+tp2)*s+tp1*tp2*s^2","tz,tp1,tp2,Ts",eps4*tz,eps4*tp1,eps4*tp2,Ts);

  }
  //the empty case, should ideally never happen
  if(n_elements <= 0) 
  {
     filter_chain_up = NULL;
     filter_chain_dp = NULL;
     filter_chain_um = NULL;
     filter_chain_dm = NULL;
  } 

  dac_out = 0;
  dwa_ptr = 0;
  dwa_ptr_neg = 0;

end:
//clean up memory
if (filter_chain_up != NULL)
   delete [] filter_chain_up;
if (filter_chain_dp != NULL)
   delete [] filter_chain_dp;
if (filter_chain_um != NULL)
   delete [] filter_chain_um;
if (filter_chain_dm != NULL)
   delete [] filter_chain_dm;
code:
if (clk_edge.inp(clk))
{
  dac_out = 0.0;
  dac_outp = 0.0;
  dac_outm = 0.0;
  code_in = (int) floor(in+0.5);

  for (j=0;j<n_elements;j++)
  {
     if(dwa_en == 1)
     {
	if(code_in > 0)
	{ 
	   //could add reset to non-output filters, but then again
	   //the diff pair tail node may experience both settling 
	   //time constants...?

           //UP glitch at outp
           filter_chain_up[dwa_ptr].inp(dac_vec.get_elem(dwa_ptr));
	   //filter_chain_dp[dwa_ptr].inp(dac_vec.get_elem(dwa_ptr));
	   filter_chain_dp[dwa_ptr].reset(dac_vec.get_elem(dwa_ptr));
           //DOWN glitch at outm
	   //filter_chain_um[dwa_ptr].inp(0.0);
	   filter_chain_um[dwa_ptr].reset(0.0);
           filter_chain_dm[dwa_ptr].inp(0.0);

	   //store result in appropriate vector
	   dac_up_vec.set_elem(dwa_ptr,dac_vec.get_elem(dwa_ptr));
	   dac_dp_vec.set_elem(dwa_ptr,dac_vec.get_elem(dwa_ptr));
	   dac_um_vec.set_elem(dwa_ptr,0);
	   dac_dm_vec.set_elem(dwa_ptr,0);
           
	   up_en_vec.set_elem(dwa_ptr,1);
	   dp_en_vec.set_elem(dwa_ptr,0);
	   um_en_vec.set_elem(dwa_ptr,0);
	   dm_en_vec.set_elem(dwa_ptr,1);

           //increment pointer
	   dwa_ptr++;
	   dwa_ptr = dwa_ptr%n_elements;
	   dwa_ptr_neg = dwa_ptr;
	}
	else
	{
           //UP glitch at outm
           filter_chain_um[dwa_ptr_neg].inp(dac_vec.get_elem(dwa_ptr_neg));
	   //filter_chain_dm[dwa_ptr_neg].inp(dac_vec.get_elem(dwa_ptr_neg));
	   filter_chain_dm[dwa_ptr_neg].reset(dac_vec.get_elem(dwa_ptr_neg));
	   //DOWN glitch at outp
	   //filter_chain_up[dwa_ptr_neg].inp(0);
	   filter_chain_up[dwa_ptr_neg].reset(0);
	   filter_chain_dp[dwa_ptr_neg].inp(0);
	   

	   //store result in appropriate vector
	   dac_um_vec.set_elem(dwa_ptr_neg,dac_vec.get_elem(dwa_ptr_neg));
	   dac_dm_vec.set_elem(dwa_ptr_neg,dac_vec.get_elem(dwa_ptr_neg));
	   dac_up_vec.set_elem(dwa_ptr_neg,0);
	   dac_dp_vec.set_elem(dwa_ptr_neg,0);

	   um_en_vec.set_elem(dwa_ptr_neg,1);
	   dm_en_vec.set_elem(dwa_ptr_neg,0);
	   up_en_vec.set_elem(dwa_ptr_neg,0);
	   dp_en_vec.set_elem(dwa_ptr_neg,1);
           
           //increment negative pointer
 	   dwa_ptr_neg++;
	   dwa_ptr_neg = dwa_ptr_neg%n_elements;
	}
     }
     else
     {
        //no DWA
	if(code_in > 0)
        {
	   //UP glitch at outp
           filter_chain_up[j].inp(dac_vec.get_elem(j));
	   //filter_chain_dp[j].inp(dac_vec.get_elem(j));
	   filter_chain_dp[j].reset(dac_vec.get_elem(j));
	   //DOWN glitch at outm
	   //filter_chain_um[j].inp(0);
	   filter_chain_um[j].reset(0);
	   filter_chain_dm[j].inp(0);

	   //store result in appropriate vector
	   dac_up_vec.set_elem(j,dac_vec.get_elem(j));
	   dac_dp_vec.set_elem(j,dac_vec.get_elem(j));
	   dac_um_vec.set_elem(j,0);
	   dac_dm_vec.set_elem(j,0);
           
	   up_en_vec.set_elem(j,1);
	   dp_en_vec.set_elem(j,0);
	   um_en_vec.set_elem(j,0);
	   dm_en_vec.set_elem(j,1);
        }
	else
        {
	   //UP glitch at outm
           filter_chain_um[j].inp(dac_vec.get_elem(j));
	   //filter_chain_dm[j].inp(dac_vec.get_elem(j));
	   filter_chain_dm[j].reset(dac_vec.get_elem(j));
	   //DOWN flitch at outm
	   //filter_chain_up[j].inp(0);
	   filter_chain_up[j].reset(0);
	   filter_chain_dp[j].inp(0);

	   //store result in appropriate vector
	   dac_um_vec.set_elem(j,dac_vec.get_elem(j));
	   dac_dm_vec.set_elem(j,dac_vec.get_elem(j));
	   dac_up_vec.set_elem(j,0);
	   dac_dp_vec.set_elem(j,0);	      
           
	   um_en_vec.set_elem(j,1);
	   dm_en_vec.set_elem(j,0);
	   up_en_vec.set_elem(j,0);
	   dp_en_vec.set_elem(j,1);
        }
     }
     code_in--;     
  }
}

dac_out = 0.0;
for(m=0;m<n_elements;m++)
{
   filter_chain_up[m].inp(dac_up_vec.get_elem(m));
   filter_chain_dp[m].inp(dac_dp_vec.get_elem(m));
   filter_chain_um[m].inp(dac_um_vec.get_elem(m));
   filter_chain_dm[m].inp(dac_dm_vec.get_elem(m));

   //current from finite impednace
   ipar = vg_in*gout_vec.get_elem(m);

   //you will see parasitic current when the up/um device is selected
   dac_out += (filter_chain_up[m].out/(n_elements-1)+ipar)*up_en_vec.get_elem(m);
   dac_out += (filter_chain_dp[m].out/(n_elements-1))*dp_en_vec.get_elem(m);
   dac_out -= (filter_chain_um[m].out/(n_elements-1)+ipar)*um_en_vec.get_elem(m);
   dac_out -= (filter_chain_dm[m].out/(n_elements-1))*dm_en_vec.get_elem(m);   
}
