surface occlfunction(float samples = 256, maxvariation = 0.02)
{
      normal Ns = faceforward(normalize(N),I); 

      float occl = occlusion(P, Ns, samples, "maxvariation", maxvariation);

      Ci = (1 - occl) * Cs * Os;
      Oi = Os;
}