public void koch(double xs,double ys,double xe,double ye){

// compute vector from s to e
double xd, yd;   // the vector (e-s)
xd = xe - xs;
yd = ye - ys;

// length of d
double l = Math.sqrt(xd*xd+yd*yd);

if (RECURSIVE CASE CONDITION){
        // store endpoints in a list
        // return
 }
else{

// compute p1, p2 and p3
// p1 and p2 are on the segment (s,e), p3 describes the tip of the 
// new triangle

double xp1= xs+xd/3;
double yp1= ys+yd/3;
double xp2= xs+xd/3*2;
double yp2= ys+yd/3*2;

// to compute p3, we need 'do', the orthogonal vector to d
double xdo = -yd;
double ydo = xd;
double ldo = Math.sqrt(xdo*xdo+ydo*ydo);
l=l/3; // length of a single segment
double h = math.sqrt(3*l*l/4);

double xp3 = xs+0.5*xd + xdo/ldo*h;
double yp3 = ys+0.5*yd + ydo/ldo*h;


// here we have all points s,p1,p3,p2,e.
// the new segments are:
// (s,p1), (p1,p3), (p3,p2), (p2,e)

//  HERE COMES THE RECURSION

return();
}






p1 = s+d/3;
    p2 = s+2*d/3;
    
    do = [-d(2) d(1)];
    l=l/3;
    h =  sqrt(l^2 - l^2/4);
    p3 = s + 0.5*d + do/norm(do)*h;
    
    koch(s,p1);
    koch(p1,p3);
    koch(p3,p2);
    koch(p2,e);
end