module: drz_dac
parameters: 
  double per_mm
  double per_mm_sw
  int n_elements
  int dwa_en
  double i_nom
  double rout_nom
  int mm_seed
  double tz
  double tp1 
  double tp2
inputs:
  double in
  double vg_in
  double_interp clk
outputs:
  double dac_out
static_variables:
  int code_in
  int j
  int k
  int m
  int q
  int r
  int dwa_ptr
  int dwa_ptr_neg
  double dac_sum
  double dac_outm
  double dac_outp
  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
  Filter *filter_chain_upb
  Filter *filter_chain_dpb
  Filter *filter_chain_umb
  Filter *filter_chain_dmb
classes:
  Rand rand_mm_dac("gauss")
  Rand rand_mm_filt("gauss")
  EdgeDetect clk_edgep()
  EdgeDetect clk_edgem()
  Vector dac_vec()
  Vector dac_up_vec()
  Vector dac_dp_vec()
  Vector dac_um_vec()
  Vector dac_dm_vec()
  Vector dac_upb_vec()
  Vector dac_dpb_vec()
  Vector dac_umb_vec()
  Vector dac_dmb_vec()
  Vector up_en_vec()
  Vector dp_en_vec()
  Vector um_en_vec()
  Vector dm_en_vec()
  Vector upb_en_vec()
  Vector dpb_en_vec()
  Vector umb_en_vec()
  Vector dmb_en_vec()
  Vector dac_mm_vec()
  Vector gout_vec()
init:
  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);

  dac_upb_vec.set_length(n_elements);
  dac_dpb_vec.set_length(n_elements);
  dac_umb_vec.set_length(n_elements);
  dac_dmb_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);

  upb_en_vec.set_length(n_elements);
  dpb_en_vec.set_length(n_elements);
  umb_en_vec.set_length(n_elements);
  dmb_en_vec.set_length(n_elements);

  //establish mismatch 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];

  filter_chain_upb = new Filter[n_elements];
  filter_chain_umb = new Filter[n_elements];
  filter_chain_dpb = new Filter[n_elements];
  filter_chain_dmb = 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,a;

  for (a = 0; a < n_elements; a++)
  {
     dac_mm_vec.set_elem(a,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*(1.0+per_mm/100*rand_mm_dac.inp());

     if (rout_val < rout_val_min)
     {
        gout_vec.set_elem(a,1/rout_val_min/i_nom);
     }
     else
     {
        gout_vec.set_elem(a,1/rout_val/i_nom);
     }
     dac_sum += dac_mm_vec.get_elem(a);
  } 

  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);

    dac_upb_vec.set_elem(i,0);
    dac_dpb_vec.set_elem(i,0);
    dac_umb_vec.set_elem(i,0);
    dac_dmb_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);

    upb_en_vec.set_elem(i,0);
    dpb_en_vec.set_elem(i,0);
    umb_en_vec.set_elem(i,0);
    dmb_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);

    filter_chain_upb[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_dpb[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_umb[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_dmb[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;

     filter_chain_upb = NULL;
     filter_chain_dpb = NULL;
     filter_chain_umb = NULL;
     filter_chain_dmb = 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;
if (filter_chain_upb != NULL)
   delete [] filter_chain_upb;
if (filter_chain_dpb != NULL)
   delete [] filter_chain_dpb;
if (filter_chain_umb != NULL)
   delete [] filter_chain_umb;
if (filter_chain_dmb != NULL)
   delete [] filter_chain_dmb;

code:
if (clk_edgep.inp(clk))
{
  code_in = (int) floor(in+0.5);
  for (j=0;j<n_elements;j++)
  {

     //zero all clkb filters
     filter_chain_upb[j].inp(0);
     filter_chain_dpb[j].inp(0);
     filter_chain_umb[j].inp(0);
     filter_chain_dmb[j].inp(0);	   

     //see only falling transients (1->0) only
     upb_en_vec.set_elem(j,0);
     dpb_en_vec.set_elem(j,1);
     umb_en_vec.set_elem(j,0);
     dmb_en_vec.set_elem(j,1);

     //zero clkb dac vectors
     dac_upb_vec.set_elem(j,0);
     dac_dpb_vec.set_elem(j,0);
     dac_umb_vec.set_elem(j,0);
     dac_dmb_vec.set_elem(j,0);     

     if(dwa_en == 1)
     {
	if(code_in > 0)
	{ 
	   //reset added to non-output filters

           //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);
           filter_chain_um[dwa_ptr].reset(0);
	   filter_chain_dm[dwa_ptr].inp(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(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);
        }
	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(dwa_ptr,1);
	   dm_en_vec.set_elem(dwa_ptr,0);
	   up_en_vec.set_elem(dwa_ptr,0);
	   dp_en_vec.set_elem(dwa_ptr,1);
        }
     }
     code_in--;     
  }
}
//CLKB 
if(clk_edgem.inp(-clk))
{

  for (k=0;k<n_elements;k++)
  {

     //copy to clkb vector
     dac_upb_vec.set_elem(k,dac_up_vec.get_elem(k));
     dac_dpb_vec.set_elem(k,dac_dp_vec.get_elem(k));
     dac_umb_vec.set_elem(k,dac_um_vec.get_elem(k));
     dac_dmb_vec.set_elem(k,dac_dm_vec.get_elem(k));

     //input zero to all clk filters
     dac_up_vec.set_elem(k,0);
     dac_dp_vec.set_elem(k,0);
     dac_um_vec.set_elem(k,0);
     dac_dm_vec.set_elem(k,0);

     //filter sequentially since data vector already determined
     filter_chain_up[k].inp(dac_up_vec.get_elem(k));
     filter_chain_dp[k].inp(dac_dp_vec.get_elem(k));
     filter_chain_um[k].inp(dac_um_vec.get_elem(k));
     filter_chain_dm[k].inp(dac_dm_vec.get_elem(k));

     filter_chain_upb[k].inp(dac_upb_vec.get_elem(k));
     filter_chain_dpb[k].inp(dac_dpb_vec.get_elem(k));
     filter_chain_umb[k].inp(dac_umb_vec.get_elem(k));
     filter_chain_dmb[k].inp(dac_dmb_vec.get_elem(k));	   

     //copy: should see same transient as in clk
     upb_en_vec.set_elem(k,up_en_vec.get_elem(k));
     dpb_en_vec.set_elem(k,dp_en_vec.get_elem(k));
     umb_en_vec.set_elem(k,um_en_vec.get_elem(k));
     dmb_en_vec.set_elem(k,dm_en_vec.get_elem(k));

     //DAC output should see falling transients (1->0) only
     up_en_vec.set_elem(k,0);
     dp_en_vec.set_elem(k,1);
     um_en_vec.set_elem(k,0);
     dm_en_vec.set_elem(k,1);
  }
}

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));

   filter_chain_upb[m].inp(dac_upb_vec.get_elem(m));
   filter_chain_dpb[m].inp(dac_dpb_vec.get_elem(m));
   filter_chain_umb[m].inp(dac_umb_vec.get_elem(m));
   filter_chain_dmb[m].inp(dac_dmb_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);

   dac_out += (filter_chain_upb[m].out/(n_elements-1)+ipar)*upb_en_vec.get_elem(m);
   dac_out += (filter_chain_dpb[m].out/(n_elements-1))*dpb_en_vec.get_elem(m);
   dac_out -= (filter_chain_umb[m].out/(n_elements-1)+ipar)*umb_en_vec.get_elem(m);
   dac_out -= (filter_chain_dmb[m].out/(n_elements-1))*dmb_en_vec.get_elem(m);
}
