END OF AN ERA, FRACTALFORUMS.COM IS CONTINUED ON FRACTALFORUMS.ORG

    x1 = gtestx - tx
    y1 = gtesty - ty
    z1 = gtestz - tz
    found = 1e200
    nomore = true
    if x1*x1 + y1*y1 + z1*z1  - (x1*dx + y1*dy + z1*dz)^2<=gdist
      if @useiterations1
        k = gntm1
        repeat
          count1[k] = gcount1[k]
        until (k=k-1)<0
      endif
      xv[0] = t5 = tx
      yv[0] = t6 = ty
      zv[0] = t7 = tz
      dxv[0] = dx
      dyv[0] = dy
      dzv[0] = dz
      k = gft
      s[0] = 1.0
      j = 0
      repeat
        ss = s[j]*gt[(i[j]=k),12]
        if @useiterations
          cnt = count[j,k] - 1
        endif
        if @useiterations1
          cnt1 = count1[k] - 1
        endif
        xv[(j=j+1)] = t5 = (x1=t5-gt[k,9])*gt[k,0] + (y1=t6-gt[k,10])*gt[k,3] \
                           + (z1=t7-gt[k,11])*gt[k,6]
        yv[j] = t6 = x1*gt[k,1] + y1*gt[k,4] + z1*gt[k,7]
        zv[j] = t7 = x1*gt[k,2] + y1*gt[k,5] + z1*gt[k,8]
        dxv[j] = dx = (x1=dx)*gt[k,0] + (y1=dy)*gt[k,3] + dz*gt[k,6]
        dyv[j] = dy = x1*gt[k,1] + y1*gt[k,4] + dz*gt[k,7]
        dzv[j] = dz = x1*gt[k,2] + y1*gt[k,5] + dz*gt[k,8]
        x1 = gtestx - t5
        y1 = gtesty - t6
        z1 = gtestz - t7
        f = false
        s[j] = ss
        if @useiterations
          count[j,k] = cnt
        endif
        if @useiterations1
          count1[k] = cnt1
        endif

// evaldraw syntax ("C" like :))
//////////////////////////////////////////////////////////////////////////////////////////////////
//"transmat" is the transformation (affine) matrix. It's a composition of translations, scale and rotations
//"tabplanes" is the array storing the folding planes parametres
//"nbr_planes" is the number of folding planes
#define INFINITY 2^30
#define IINFINITY 2^-30
#define BOUNDINGRAD 1
DivSeg(ox,oy,oz,vx,vy,vz,tmin,tmax,i,j,&k){
   if(i==nbr_planes){//Test intersection with bounding volume
      //Here the bounding volume is simply the sphere centered at the origin and with radius BOUNDINGRAD
      
      //Apply transformation to ray's origin
      xx=transmat[0][0]*ox+transmat[0][1]*oy+transmat[0][2]*oz+transmat[0][3];
      yy=transmat[1][0]*ox+transmat[1][1]*oy+transmat[1][2]*oz+transmat[1][3];
      zz=transmat[2][0]*ox+transmat[2][1]*oy+transmat[2][2]*oz+transmat[2][3];
      ox=xx;oy=yy;oz=zz;
      
      //Apply transformation to ray's direction
      xx=transmat[0][0]*vx+transmat[0][1]*vy+transmat[0][2]*vz;
      yy=transmat[1][0]*vx+transmat[1][1]*vy+transmat[1][2]*vz;
      zz=transmat[2][0]*vx+transmat[2][1]*vy+transmat[2][2]*vz;
      vx=xx;vy=yy;vz=zz;
      
      {//Intersection with the sphere
         c=ox*ox+oy*oy+oz*oz-BOUNDINGRAD^2;
         b=ox*vx+oy*vy+oz*vz;
         a=vx*vx+vy*vy+vz*vz;
         delta=b*b-a*c;
         if(delta>2^-53){
            a=1/a;
            delta=sqrt(delta);
            t0=(-b-delta)*a; t1=(-b+delta)*a;
            if (tmax<t0 || tmin>t1) return INFINITY;//intersection outside [tmin,tmax] => no intersection
            if (j==MI) return t0;//max recursion reached 
            return DivSeg(ox,oy,oz,vx,vy,vz,max(tmin,t0),min(tmax,t1),0,j+1);//next level of recursion
         } else return INFINITY;//No intersection at all
      }
   } else {//Fold the ray by the folding plane #i
      //we fold only the part of the ray between tmin and tmax
      r=INFINITY;
      a=vx*tabplanes[i][0]+vy*tabplanes[i][1]+vz*tabplanes[i][2];
      b=ox*tabplanes[i][0]+oy*tabplanes[i][1]+oz*tabplanes[i][2]+tabplanes[i][3];
      vd0=a*tmin+b; //In which side of the folding plane is the tmin extremity of the ray is
      vd1=a*tmax+b; //In which side of the folding plane is the tmax extremity of the ray is
      //account for different cases
      if(a>0){// angle between ray direction and folding plane normal <pi/2 
         if(vd0>0) r=DivSeg(ox,oy,oz,vx,vy,vz,tmin,tmax,i+1,j);//on positive side of the folding plane => ray origin and direction are not folded
         else {
            if(vd1<0) r=DivSeg(ox-2*b*tabplanes[i][0],oy-2*b*tabplanes[i][1],oz-2*b*tabplanes[i][2],
                               vx-2*a*tabplanes[i][0],vy-2*a*tabplanes[i][1],vz-2*a*tabplanes[i][2],
                               tmin,tmax,i+1,j);//on positive side of the folding plane => ray origin and direction are folded
            else {//ray segment crosses the folding plane => cut it
               t=-b/a;//No need to test for a==0 because this case never occures.
               //we don't process segments that are too small 
               if (t-tmin>IINFINITY) r=DivSeg(ox-2*b*tabplanes[i][0],oy-2*b*tabplanes[i][1],oz-2*b*tabplanes[i][2],
                        vx-2*a*tabplanes[i][0],vy-2*a*tabplanes[i][1],vz-2*a*tabplanes[i][2],
                        tmin,t,i+1,j);
               if(r>t) if (tmax-t>IINFINITY) r=DivSeg(ox,oy,oz,vx,vy,vz,t,tmax,i+1,j);
            }
         }
      } else {//symmetric case
         if(vd0<0) r=DivSeg(ox-2*b*tabplanes[i][0],oy-2*b*tabplanes[i][1],oz-2*b*tabplanes[i][2],
                            vx-2*a*tabplanes[i][0],vy-2*a*tabplanes[i][1],vz-2*a*tabplanes[i][2],
                            tmin,tmax,i+1,j);
         else {
            if(vd1>0) r=DivSeg(ox,oy,oz,vx,vy,vz,tmin,tmax,i+1,j);
            else {
               t=-b/a;
               if (t-tmin>IINFINITY) r=DivSeg(ox,oy,oz,vx,vy,vz,tmin,t,i+1,j);
               if(r>t) if (tmax-t>IINFINITY) r=DivSeg(ox-2*b*tabplanes[i][0],oy-2*b*tabplanes[i][1],oz-2*b*tabplanes[i][2],
                                vx-2*a*tabplanes[i][0],vy-2*a*tabplanes[i][1],vz-2*a*tabplanes[i][2],
                                t,tmax,i+1,j);
            }
         }
      }
      return r;
   }
}
//This function will return the position of the first intersection along the ray {origin : [ox,oy,oz]; direction : [vx,vy,vz];} that is between tmin and tmax
//It calls the recursive function defined above
nextd0(ox,oy,oz,vx,vy,vz,tmin,tmax){
   r=INFINITY;j=0;
   {//Test intersection with bounding volume
      //Here the bounding volume is simply the sphere centered at the origin and with radius BOUNDINGRAD
      {//Intersection with the sphere
         c=ox*ox+oy*oy+oz*oz-BOUNDINGRAD^2;
         b=ox*vx+oy*vy+oz*vz;
         a=vx*vx+vy*vy+vz*vz;
         delta=b*b-a*c;
         if(delta>IINFINITY){
            a=1/a;
            delta=sqrt(delta);
            t0=(-b-delta)*a; t1=(-b+delta)*a;
            if (tmax<t0 || tmin>t1) r=INFINITY;//intersection outside [tmin,tmax] => no intersection
            if (j==MI) r=t0;//max recursion reached 
            r=DivSeg(ox,oy,oz,vx,vy,vz,max(tmin,t0),min(tmax,t1),0,j+1);//next level of recursion
         } else r=INFINITY;//No intersection at all
      }
   }
   return r;
}