// ******************************
// *                            *
// * "Snowflake" von Koch curve *
// *                            *
// ******************************

clear;
clf;

// **************
// * constants *
// **************

n = 6;
// number of steps 
// limited to 9 (262 145 points), otherwise the stacksize must be changed
// 6 steps (4 097 points) are enough for a good rendering
N = 4^n+1; // amount of points
sin_soixante = sqrt(3)/2; // sin(60°)
l = 1; // length of the initial line (arbitrary unit)

// ******************
// * initialisation *
// ******************

ycourbe = [zeros(1,N)];
ycourbe1 = ycourbe;

// *************
// * functions *
// *************

function [xx, yy] = etape(x, y)
  
  // from a line [(x(1),y(1)) ; (x(2),y(2))]
  // make the line [(xx(1),yy(1)) ; (xx(2),yy(2)) ; (xx(3),yy(3))]
  // x and y are 2-cells tables, the ends of the basis line
  // xx and yy are 3-cells tables
  // the edges of the equilateral triangle
  
  xu = (x(2)-x(1))/3;
  yu = (y(2)-y(1))/3;
  // third of the basis line vector
  
  xv = 0.5*xu - sin_soixante*yu;
  yv = sin_soixante*xu + 0.5*yu;
  // vector turned by +60°
  
  xx(1) = x(1)+xu; yy(1) = y(1)+yu;
  xx(3) = x(2)-xu; yy(3) = y(2)-yu;
  
  xx(2) = xx(1) + xv;
  yy(2) = yy(1) + yv;  
endfunction

function [xkoch, ykoch] = vonkoch(x, y, n)
  // builds the curve
  // initialisation
  xkoch = x;
  ykoch = y;
  xkoch1 = x;
  ykoch1 = y;
  for i=1:n
    jmax = 4^(i-1);
    // number of lines at the beginning of the step i
    for j=1:jmax/2+1
      // we work with two points which indices are j and j+1 (line #j)
      // thanks t the symmetry, we work with a half curve
      decalage = (j-1)*4; 
      // the new points shift the next points by this offset
      x_init = xkoch(j:j+1);
      y_init = ykoch(j:j+1);
      // line #j
      [x_trans, y_trans] = etape(x_init,y_init);
      // transformed line
      xkoch1(decalage+1) = x_init(1); xkoch1(decalage+5) = x_init(2);
      ykoch1(decalage+1) = y_init(1); ykoch1(decalage+5) = y_init(2);
      for k=1:3
        xkoch1(k+decalage+1) = x_trans(k);
        ykoch1(k+decalage+1) = y_trans(k);
        // values put in the global vector
      end
    end
    xkoch = xkoch1; ykoch = ykoch1;
  end
  
  for i=1:4^n
    ykoch(N-i+1) = ykoch(i);
    xkoch(N-i+1) = l-xkoch(i); 
    // 2nd half-curve
  end
endfunction

// ****************
// * main program *
// ****************

xcourbe(2) = l;
[xcourbe,ycourbe] = vonkoch(xcourbe,ycourbe,n);

// drawing the curve

xpoly(xcourbe,ycourbe)
isoview(0,l,0,l*sin_soixante/3)

// saving the file

name = "von_koch_"+string(n)+"_steps.svg";
xs2svg(0, name)

//============================================================================
// name: von_koch.sce
// author: Christophe Dang Ngoc Chan
// date of creation: 2012-10-23
// dates of modification: 
//    2013-07-08: quotes ' -> "
//    2013-07-2: vectorisation of the calculations
//----------------------------------------------------------------------------
// version of Scilab: 5.3.1
// required Atoms modules: aucun
//----------------------------------------------------------------------------
// Objective: draws the von Koch's "snowflake"
// Inputs: none (parameters are hard coded)
// Outputs: graphical window with a curve; SVG file
//============================================================================

clear;
clf;

// *************
// * constants *
// **************

n = 6; 
// number of steps 
// limited to 9 (262 145 points), otherwise the stacksize must be changed
// 6 steps (4 097 points) are enough for a good rendering
sin_soixante = sqrt(3)/2; // sin(60°)
l = 1;
// length of the initial line (arbitrary unit)

// ******************
// * initialisation *
// ******************



// *************
// * functions *
// *************

function [] = vonkoch(A, B, i)
    u = (B - A)/3 ; // third of the AB vector
    v = [0.5*u(1) - sin_soixante*u(2) ; sin_soixante*u(1) + 0.5*u(2)] ;
    // vector turned by +60°
    C = A + u ;
    D = C + v ;
    E = B - u ;
    // points of the line
    if i == 1 then
        // drawing the smallest segments
        x = [A(1) ; C(1) ; D(1) ; E(1) ; B(1) ];
        y = [A(2) ; C(2) ; D(2) ; E(2) ; B(2) ];
        xpoly(x, y, "lines")
    else
        j = i - 1 ;
        vonkoch(A, C, j);
        vonkoch(C, D, j);
        vonkoch(D, E, j);
        vonkoch(E, B, j);
        // recursive call
    end
endfunction

// ****************
// * main program *
// ****************

beginning = [0;0] ;
ending = [l;0] ;
vonkoch(beginning, ending, n)

isoview(0,l,0,sin_soixante*l)

// Saving the file

name = "von_koch_"+string(n)+"_steps.svg" ;
xs2svg(0, name)