#define step(a,b,d,e)	\
	e=clamp(a,b,d);	\
	e-=b;	\
	e/=d-b;	\
	

#define KTdiffuse(P1,P2,d)	\
	d=P1.P2;	\
	d=clamp(d,0,1);	\


#define KTspecular(P1,P2,P3,r,e)	\
	H=normalize(normalize(P1)+normalize(P2));	\
	e=H.P3;	\
	e=clamp(e,0,1);	\
	e=1-acos(e)/(PI/2);	\
/*	e=pow(s,1.5);*/	\
	e=pow(e,1/r);	\


#define judge_intersection_sphere(P1,P2,P3,a,judge)	\
	Px=xcomp(P1);	\
	Py=ycomp(P1);	\
	Pz=zcomp(P1);	\
	\
	Rx=xcomp(P2); 	\
	Ry=ycomp(P2); 	\
	Rz=zcomp(P2); 	\
	\
	sphereCx=xcomp(P3);	\
	sphereCy=ycomp(P3);	\
	sphereCz=zcomp(P3);	\
	\
	X=Px-sphereCx; 	\
	Y=Py-sphereCy; 	\
	Z=Pz-sphereCz; 	\
	A=Rx*Rx+Ry*Ry+Rz*Rz; 	\
	B=2*X*Rx+2*Y*Ry+2*Rz*Z; 	\
	C=X*X+Y*Y+Z*Z-a*a; 	\
	judge=B*B-4*A*C; 	\
	
	
#define KTshadow_sphere(P1,P2,a,P3,P4,P5,c,o)	\
	illuminance(P1,point(1,0,0),PI) {	\
		LL=L;	\
		LL-=P3+P4;	\
		LL=normalize(LL);	\
		judge_intersection_sphere(P1,LL,P2,a,D);	\
		if (D>=0) {	\
			alpha1=(-B-sqrt(D))/(2*A); 	\
			alpha2=(-B+sqrt(D))/(2*A); 	\
			alpha2=abs(alpha2);	\
			if (decayshd_swith=="on") {	\
				step(alpha2,decay_shd_start,decay_shd_end,decay_shd);	\
				decay_shd=smoothstep(0,1,decay_shd);	\
				decay_shd=pow(decay_shd,decay_shd_pow);	\
			}	\
			else {	\
				decay_shd=0;	\
			}	\
			shd=1-o;	\
			if (D!=0) {	\
				shd-=o;	\
				shd=color(clamp(comp(shd,0),0,1),clamp(comp(shd,1),0,1),clamp(comp(shd,2),0,1));	\
			}	\
			Cl2=Cl*((1-decay_shd)*shd+decay_shd);	\
		}	\
		else {	\
			Cl2=Cl;	\
		}	\
	\
		KTdiffuse(P5,LL,d_intensity);	\
		KTspecular(normalize(-P1),LL,P5,roughness,s_intensity);	\
		dif_cs+=d_intensity*Cl2;	\
		spe_cs+=s_intensity*Cl2;	\
	}	\


#define KTtrace_sphere(P1,P2,P3,a,N_P1) 	\
	IN=normalize(P2); 	\
	N_P1=normalize(N_P1);	\
	fresnel(IN,N_P1,eta,kr,kt,R,T); 	\
	R=normalize(R);	\
	\
	judge_intersection_sphere(P1,R,P3,a,D);	\
 	\
	if (D>=0) { 	\
		alpha1=(-B-sqrt(D))/(2*A); 	\
		alpha2=(-B+sqrt(D))/(2*A); 	\
		intersectP=P1+alpha1*R;	\
		N_intersectP=intersectP-sphereC; 	\
		N_intersectP=normalize(N_intersectP); 	\
		\
		illuminance(intersectP,N_intersectP,PI) {	\
			if (light_type=="distant") {	\
				LL=normalize(L);	\
				Cl2=Cl;	\
			}	\
			else {	\
				LL=L-alpha1*R; 	\
				Cl2=Cl*pow(length(L)/length(LL),2);	\
				LL=normalize(LL); 	\
			}	\
	\
			KTdiffuse(N_intersectP,LL,d_intensity);	\
			KTspecular(-R,LL,N_intersectP,roughness_sphere,s_intensity);	\
	\
			env_cs+=sphere_opacity*(sphere_cs*(Ka_sphere*ambient()+Kd_sphere*Cl2*d_intensity)+\
						Ks_sphere*sphere_specularcs*Cl2*s_intensity);	\
		} 	\
		opa=comp(sphere_opacity,0)*comp(sphere_opacity,1)*comp(sphere_opacity,2);	\
		if (opa<1) {	\
			intersectP=P1+alpha2*R;	\
			N_intersectP=intersectP-sphereC;	\
			N_intersectP=normalize(-N_intersectP);	\
		\
			illuminance(intersectP,N_intersectP,PI) {	\
				if (light_type=="distant") {	\
					LL=normalize(L);	\
					Cl2=Cl;	\
				}	\
				else {	\
					LL=L-alpha2*R; 	\
					Cl2=Cl*pow(length(L)/length(LL),2);	\
					LL=normalize(LL); 	\
				}	\
		\
				KTdiffuse(N_intersectP,LL,d_intensity);	\
				KTspecular(-R,LL,N_intersectP,roughness_sphere,s_intensity);	\
		\
				env_cs+=(1-sphere_opacity)*sphere_opacity*(sphere_cs*(Ka_sphere*ambient()+Kd_sphere*Cl2*d_intensity)+\
							Ks_sphere*sphere_specularcs*Cl2*s_intensity);	\
			} 	\
		}	\
	}	\
	else {	\
		env_cs+=bg_cs;	\
	}	\


surface
raydemo
(
float	Ka=1,Kd=1,Ks=1,roughness=0.1;
color	specularcolor=color(1,1,1);
point	sphereC=point(0,0,-0.5);
float	sphere_radius=0.5;
float	eta=1;

color	sphere_cs=color(1,0,0),sphere_specularcs=color(1,1,1),sphere_opacity=color(1,1,1);
color	bg_cs=color(0,0,0);
float	Ka_sphere=1,Kd_sphere=0.5,Ks_sphere=0.5,roughness_sphere=0.1;

string	ray_swith="on",shd_swith="on",light_type="distant",decayshd_swith="off";
float	shd_sample=1,ray_sample=1,blur_size=0;
float	decay_shd_start=0.5,decay_shd_end=2,decay_shd_pow=1.5;
)

{
point	Nf,R,T,LL,IN,V;
float	kr,kt;
float	Rx,Ry,Rz;
float	Px,Py,Pz;
float	sphereCx,sphereCy,sphereCz;

float	i,j,a,b,d,e,r,judge,opa;
point	Pray,Iray,ddu,ddv;

float	X,Y,Z,A,B,C,D;
float	alpha1,alpha2,decay_shd;
point	intersectP,N_intersectP,H;
float	d_intensity,s_intensity;

color	dif_cs=0,spe_cs=0,env_cs;
color	shd,c,o,Cl2;

Nf=faceforward(normalize(N),I);
V=normalize(-I);

env_cs=0;
if (ray_swith=="on") {
	if (ray_sample>0) {
		for (i=0;i<ray_sample;i+=1) {
			for (j=0;j<ray_sample;j+=1) {
				ddu=dPdu*(du/ray_sample)*i;
				ddv=dPdv*(dv/ray_sample)*j;
				Pray=P+ddu+ddv;
				Iray=I+ddu+ddv;
				KTtrace_sphere(Pray,Iray,sphereC,sphere_radius,Nf);
			}
		}
		env_cs/=pow(ray_sample,2);
	}
}

if (shd_swith=="on") {
	if (shd_sample>0) {
		for (i=0;i<shd_sample;i+=1) {
			for (j=0;j<shd_sample;j+=1) {
				ddu=dPdu*(du/shd_sample)*i;
				ddv=dPdv*(dv/shd_sample)*j;
/*				ddu=dPdu*(du/shd_sample)*blur_size*i;
				ddv=dPdv*(dv/shd_sample)*blur_size*j;*/
				Pray=P+ddu+ddv;
				KTshadow_sphere(Pray,sphereC,sphere_radius,ddu,ddv,Nf,sphere_cs,sphere_opacity);
			}
		}
		dif_cs/=pow(shd_sample,2);
		spe_cs/=pow(shd_sample,2);
	}
}
else {
	dif_cs=diffuse(Nf);
	spe_cs=specular(Nf,V,roughness);
}

Ci=Cs*(Ka*ambient()+Kd*dif_cs)+specularcolor*(env_cs+Ks*spe_cs);
/*Ci=env_cs;*/

}