void _limit(int N, double beta, double* qc, double* qv, 
	    double* qmin, double* qmax) { 

  int k, i, k3;
  double dq, dqa[3], phi, r;
  
  for (k=0; k<N; k++) {
    k3 = k*3;
    
    //Find the gradient limiter (phi) across vertices  
    phi = 1.0;
    for (i=0; i<3; i++) {    
      r = 1.0;
      
      dq = qv[k3+i] - qc[k];    //Delta between vertex and centroid values
      dqa[i] = dq;              //Save dq for use in the next loop
      
      if (dq > 0.0) r = (qmax[k] - qc[k])/dq;
      if (dq < 0.0) r = (qmin[k] - qc[k])/dq;      
  
  
      phi = min( min(r*beta, 1.0), phi);    
    }
    
    //Then update using phi limiter
    for (i=0; i<3; i++) {    
      qv[k3+i] = qc[k] + phi*dqa[i];
    }
  }
}