Refraction (or something that looks like refraction :)

[TomyLobo]

Hi there I experimented a bit with refraction and managed to write an AVS preset that displays a sphere refracting the image on a plane behind it.
I don't know whether my refraction maths are correct.
(or better: I *guessed* them *g*)
Maybe one of you math geeks can check that and/or give me a better formula

Here's the code of the DM:

code:


//I didn't find a suitable translation for the

//German word "Auftreffpunkt" which means the

//point where something hits something.

//I translated it with "hit point".

//1st hit point is the point where the ray hits

//the outer bounds of the sphere.

//2nd hit point is the point where the ray hits

//the inner bounds of the sphere after being refracted.



//--- Calculate projection plane coordinates. ---

dx=x/zoom; dy=y/zoom; dz=1;



//--- Transform them. ---

dx1=dx*cz-dy*sz; dy1=dx*sz+dy*cz; dy=dy1*cx-sx; dz1=dy1*sx+dz*cx; dx=dx1*cy-dz1*sy; dz=dx1*sy+dz1*cy;



//--- Trace from camera to sphere surface ---

// r² = x² + y² + z²

// r² = (dx*k+ox-mx)² + (dy*k+oy-my)² + (dz*k+oz-mz)²

kdist=(sqr(dx) + sqr(dy) + sqr(dz));

kp = (2*dx*(ox-mx) + 2*dy*(oy-my) + 2*dz*(oz-mz))/kdist;

kq = (sqr(ox-mx)+sqr(oy-my)+sqr(oz-mz)-sqr(rad))/kdist;

kdet=sqr(kp/2)-kq;

k1=-kp/2+sqrt(kdet);

k2=-kp/2-sqrt(kdet);

k1=if(below(k1,0),INFIN,k1);

k2=if(below(k2,0),INFIN,k2);

ksphere=min(k1,k2);



//--- Trace from camera to plane ---

pa=0; pb=0; pc=1; pd=-20;

k1=(-pa*ox-pb*oy-pc*oz-pd)/(pa*dx + pb*dy + pc*dz);

kplane=if(below(k1,0),INFIN,k1);



//vec a1: o -> 1st hit point  //vec ik1: <0> -> 1st hit point

a1x=dx*ksphere; ik1x=a1x+ox;

a1y=dy*ksphere; ik1y=a1y+oy;

a1z=dz*ksphere; ik1z=a1z+oz;



//--- Does the ray hit the sphere? ---

xdist=ik1x-mx; ydist=ik1y-my; zdist=ik1z-mz;

distkugel=sqrt(sqr(xdist)+sqr(ydist)+sqr(zdist));

isplane=above(abs(9-distkugel),epsilon);

k=if(isplane,kplane,ksphere);



ix1=dx*kplane+ox; iy1=dy*kplane+oy; iz1=dz*kplane+oz;



//The refraction method i used is a linear interpolation

//between the normal vector and the traced ray

//I'm not sure whether this method is correct. If you know a

//better method, let me know 



//--- 1st Refraction ---

//Coefficients for recfraction

pr=1/refindex; ps=1-pr;

//normalise vec a1

a1r=invsqrt(sqr(a1x)+sqr(a1y)+sqr(a1z));

a1x=a1x*a1r; a1y=a1y*a1r; a1z=a1z*a1r;

//vec n: sphere surface normal at the 1st hit point

nx=ik1x-mx;

ny=ik1y-my;

nz=ik1z-mz;

//normalise vec n

nr=invsqrt(sqr(nx)+sqr(ny)+sqr(nz));

nx=nx*nr; ny=ny*nr; nz=nz*nr;

//vec b1: Ray after 1st refraction

b1x=pr*a1x-ps*nx;

b1y=pr*a1y-ps*ny;

b1z=pr*a1z-ps*nz;

//normalise vec b1

b1r=invsqrt(sqr(b1x) + sqr(b1y) + sqr(b1z));

b1x=b1x*b1r; b1y=b1y*b1r; b1z=b1z*b1r;



//--- Trace inside the sphere---

ka=sqr(b1x) + sqr(b1y) + sqr(b1z);

kp=(2*b1x*(ik1x-mx) + 2*b1y*(ik1y-my) + 2*b1z*(ik1z-mz))/ka;

kq=(sqr(ik1x-mx)+sqr(ik1y-my)+sqr(ik1z-mz)-sqr(rad))/ka;

kdet=sqr(kp/2)-kq;

k1=-kp/2+sqrt(kdet);

//k2=-kp/2-sqrt(kdet);

k1=if(below(k1,0),INFIN,k1);

//k2=if(below(k2,0),INFIN,k2);

ksphere2=k1;



//--- 2nd refraction ---

//Coefficients for recfraction

pr=1/refindex; ps=1-pr;

//vec a2: 1st hit point -> 2nd hit point // vec ik2: <0> -> 2nd hit point

a2x=b1x*ksphere2; ik2x=a2x+ik1x;

a2y=b1y*ksphere2; ik2y=a2y+ik1y;

a2z=b1z*ksphere2; ik2z=a2z+ik1z;

//normalise vec a2 (keep a2r for light dimming calculations)

a2r=sqrt(sqr(a2x)+sqr(a2y)+sqr(a2z));

a2x=a2x/a2r; a2y=a2y/a2r; a2z=a2z/a2r;

//vec n: sphere surface normal at the 2nd hit point

nx=ik2x-mx;

ny=ik2y-my;

nz=ik2z-mz;

//normalise vec n

nr=invsqrt(sqr(nx)+sqr(ny)+sqr(nz));

nx=nx*nr; ny=ny*nr; nz=nz*nr;

//vec b2: Ray after 2nd refraction

b2x=pr*a2x+ps*nx;

b2y=pr*a2y+ps*ny;

b2z=pr*a2z+ps*nz;



//--- Trace refracted ray to the plane  ---

k1=(-pa*ik2x-pb*ik2y-pc*ik2z-pd)/(pa*b2x + pb*b2y + pc*b2z);

kplane2=if(below(k1,0),INFIN,k1);

ix2=b2x*kplane2+ik2x; iy2=b2y*kplane2+ik2y; iz2=b2z*kplane2+ik2z;



//--- Get coordinates of the point where the ray ends

ix=if(isplane,ix1,ix2);

iy=if(isplane,iy1,iy2);

iz=if(isplane,iz1,iz2);



//--- Determine texture coordinates ---

x=ix/40; y=iy/40;



//--- Determine Alpha ---

// Calculate distance to camera.

xdist=ix-ox; ydist=iy-oy; zdist=iz-oz;

dist=sqrt(sqr(xdist)+sqr(ydist)+sqr(zdist));

//if (plane) then (100% Alpha) else (calculate from distance between the 2 hitpoints)

alpha=if(isplane,1,exp(-a2r/50));

//if (dist>=viewdist) then (0% Alpha)

alpha=alpha*below(dist,viewdist);

alpha=if(equal(k,INFIN),0,alpha);