#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,P3,a,P4,P5,P6,c,o)        \
        illuminance(P1,P6,PI/2) {       \
                LL=normalize(L);        \
                PO=transform("current","object",P1);    \
                LO=vtransform("current","object",LL);   \
                SO=vtransform("current","object",P3);   \
                judge_intersection_sphere(PO,LO,SO,a,D);        \
/*              judge_intersection_sphere(P1,LL,P3,a,D);*/      \
                if (D>=0) {     \
                        alpha1=(-B-sqrt(D))/(2*A);      \
                        alpha2=(-B+sqrt(D))/(2*A);      \
                        alpha2=abs(alpha2);     \
                        if (decayshd_switch=="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(P6,LL,d_intensity);   \
                KTspecular(normalize(-P2),LL,P6,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); \
        \
        PO=transform("current","object",P1);    \
        RO=vtransform("current","object",R);    \
        SO=vtransform("current","object",P3);   \
        \
        judge_intersection_sphere(PO,RO,SO,a,D);        \
/*      judge_intersection_sphere(P1,R,P3,a,D);*/       \
        \
        if (D>=0) {     \
                alpha1=(-B-sqrt(D))/(2*A);      \
                alpha2=(-B+sqrt(D))/(2*A);      \
                intersectPO=PO+alpha1*RO;       \
                intersectP=transform("object","current",intersectPO);   \
                N_intersectPO_p=intersectPO-SO; \
                N_intersectPO=vector (xcomp(N_intersectPO_p),ycomp(N_intersectPO_p),zcomp(N_intersectPO_p));    \
                N_intersectPO=normalize(N_intersectPO); \
                N_intersectP=vtransform("object","current",N_intersectPO);      \
/*              intersectP=P1+alpha1*R; \
                N_intersectP_p=intersectP-P3;   \
                N_intersectP=vector(xcomp(N_intersectP_p),ycomp(N_intersectP_p),zcomp(N_intersectP_p));*/       \
/*              N_intersectP=(P2+alpha1*R)-P3;  \
                N_intersectP=intersectP-sphereC;*/      \
                N_intersectP=normalize(N_intersectP);   \
                \
                illuminance(intersectP,N_intersectP,PI/2) {     \
                        LL=normalize(L);        \
                        LO=vtransform("current","object",LL);   \
                        Cl2=Cl; \
                        KTdiffuse(N_intersectPO,LO,d_intensity);        \
                        KTspecular(-RO,LO,N_intersectPO,roughness_sphere,s_intensity);  \
/*                      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) {    \
                        intersectPO=PO+alpha2*RO;       \
                        intersectP=transform("object","current",intersectPO);   \
                        N_intersectPO_p=intersectPO-SO; \
                        N_intersectPO=vector (xcomp(N_intersectPO_p),ycomp(N_intersectPO_p),zcomp(N_intersectPO_p));    \
                        N_intersectPO=normalize(-N_intersectPO);        \
                        N_intersectP=vtransform("object","current",N_intersectPO);      \
/*                      intersectP=P1+alpha2*R; \
                        N_intersectP_p=intersectP-P3;   \
                        N_intersectP=vector(xcomp(N_intersectP_p),ycomp(N_intersectP_p),zcomp(N_intersectP_p));*/       \
/*                      N_intersectP=(P2+alpha2*R)-P3;  \
                        N_intersectP=intersectP-sphereC;*/      \
                        N_intersectP=normalize(N_intersectP);   \
                \
                        illuminance(intersectP,N_intersectP,PI/2) {     \
                                LL=normalize(L);        \
                                LO=vtransform("current","object",LL);   \
                                Cl2=Cl; \
                \
                                KTdiffuse(N_intersectPO,LO,d_intensity);        \
                                KTspecular(-RO,LO,N_intersectPO,roughness_sphere,s_intensity);  \
/*                              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);
vector  sphereC=vector "shader" (0,0,-0.5);
/*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_switch="on",shd_switch="on",decayshd_switch="off";
/*string        light_type="distant";*/
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;
)

{
normal  Nf;
vector  V;
vector  ddu,ddv;
point   Pray;
vector  Iray;
vector  IN,R,T;
point   intersectP;
vector  N_intersectP;
point   N_intersectP_p;
/*normal        N_intersectP;*/
vector  LL;
vector  H;

float   kr,kt;
float   Rx,Ry,Rz;
float   Px,Py,Pz;
float   sphereCx,sphereCy,sphereCz;

float   i,j,opa;
/*float i,j,a,b,d,e,r,judge,opa;*/

float   X,Y,Z,A,B,C,D;
float   alpha1,alpha2,decay_shd;
float   d_intensity,s_intensity;

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

point   PO,intersectPO,N_intersectPO_p;
vector  N_intersectPO;
vector  RO,SO,LO;

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

env_cs=0;
if (ray_switch=="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_switch=="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;
                                Iray=I+ddu+ddv;
                                KTshadow_sphere(Pray,Iray,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;*/

}