tracking refraction

vex

vector rd = normalize(set(-1, ch("dirY"), 0));
vector ro = @P;

int prim = addprim(0, "polyline");
addvertex(geoself(), prim, @ptnum);

float ior = ch("ior");

for(int i = 0 ; i < 20; ++i) {
    vector p;
    vector uvw;
    int intersect_prim = intersect(1, ro, rd * 100.0, p, uvw);
    if(intersect_prim != -1) {
        int p1 = addpoint(geoself(), p);
        addvertex(geoself(), prim, p1);
        setpointattrib(geoself(), "sinTheta", p1, rd.y);
        
        vector N = prim_normal(1, intersect_prim, uvw);
        
        /*
        vector rdir;
        if(dot(N, rd) < 0.0) {
            rdir = refract(rd, N, 1.0 / ior);
        } else {
            rdir = refract(rd, -N, ior / 1.0);
        }
        */
        vector rdir = refract(rd, -N, ior / 1.0);
        
        ro = p + rdir * 0.0001;
        rd = rdir;
    } else {
        vector pLast = ro + rd * 2.0;
        int fin = addpoint(geoself(), pLast);
        addvertex(geoself(), prim, fin);
        setpointattrib(geoself(), "sinTheta", fin, rd.y);
        break;
    }
}