surface occsurfwrite(string filename = "", coordsystem = ""; float samples = 64, maxdist = 1e30)
{
	normal Nn = normalize(N);
	
	// Compute occlusion
	float occ = occlusion(P, Nn, samples, "maxdist", maxdist, "coordsystem", coordsystem);

	// Set Ci and Oi
	Ci = (1 - occ) * Cs * Os;
	Oi=Os;
}