// load PGMIO-Functions
getf('pgmio.sce');


function res = psi_prime(u,lambda); res = 1 / sqrt(1 + u/lambda); endfunction
function I = readLena(); [image, bla] = readpgm('lena256.pgm');I = pump_up(image);endfunction
function I = readCouple(); [image, bla] = readpgm('couple40.pgm');I = pump_up(image);endfunction
function I = readTrui(); [image, bla] = readpgm('trui.pgm');I = pump_up(image); endfunction
  
function A = mirror_boundaries(A)
  // mirror first column
  A(:,1) = A(:,2);
  // mirror last column
  A(:,$) = A(:,$-1);
  // mirror first row
  A(1,:) = A(2,:);
  // mirror last row;
  A($,:) = A($-1,:);
endfunction

function A = pump_up (M)
  [n,m] = size(M);
  A = zeros(n+2,m+2);
  // translate the matrix indices by 1 in each direction
  A(2:n+1,2:m+1) = M;
  A = mirror_boundaries(A);
endfunction

function u = iteration(u, lambda, alpha, dx, dy, dt, t)
  [n,m] = size(u);

  tmp = u; l2 = lambda * lambda;
  
  for iter=1:t
    printf("Iteration %d\n",iter);
          
    // compute g
    // (TODO: There have to be a function or syntax to
    //        apply a scalar-function to each element of
    //        a matrix like other things X .* Y)
    g = zeros(n,m);
 
    ddx = 1 / (dx*dx); ddy = 1 / (dy*dy);
    for i=2:n-1
        for j=2:m-1
          uu = (u(i,j) - u(i-1,j))^2 * ddx + (u(i,j) - u(i,j-1) )^2 * ddy ;
          g(i,j) = psi_prime(uu, l2);
        end
    end
    g = mirror_boundaries(g);

    for i=2:n-1
        for j=2:m-1
          tmp(i,j) = u(i,j) * (1 - alpha * (ddx * ( g(i+1,j)+g(i,j) ) + ddy * ( g(i,j+1)+g(i,j) ) ) ) + alpha * (ddx * (g(i+1,j) * u(i+1,j) + g(i,j)*u(i-1,j)) + ddy * (g(i,j+1) * u(i,j+1) + g(i,j)*u(i,j-1)));
        end
    end  
    
    u = mirror_boundaries(tmp);
  end
    u = u(2:$-1, 2:$-1);
endfunction

function printout(lambda)
  u = readCouple();
  [n,m] = size(u);

  tmp = u; l2 = lambda * lambda;
  
  for iter=1:50
    printf("Iteration %d\n",iter);

    g = zeros(n,m);
    for i=2:n-1
        for j=2:m-1
          uu = (u(i,j) - u(i-1,j))^2 + (u(i,j) - u(i,j-1) )^2;
          g(i,j) = psi_prime(uu, l2);
        end
    end
    g = mirror_boundaries(g);

    for i=2:n-1
        for j=2:m-1
          tmp(i,j) = u(i,j) * (1 - 0.25 * (g(i+1,j)+g(i,j)+g(i,j+1)+g(i,j))) + 0.25 * ((g(i+1,j) * u(i+1,j) + g(i,j)*u(i-1,j)+ g(i,j+1) * u(i,j+1) + g(i,j)*u(i,j-1)));
        end
    end  
    
    u = mirror_boundaries(tmp);
    if ( (iter == 10) | (iter == 20) | (iter == 50) )
      a = u(2:$-1, 2:$-1);
      filename = sprintf('couple40_%d_%d.pgm', lambda, iter);
      comment = sprintf('lambda=%d,   time=%f", lambda, 0.25*iter);
      writepgm(a,filename,255,comment);
    end
  end
endfunction

function ret = ex23(f,lambda, t); ret = iteration(f,lambda,t); endfunction

function ex24(); printout(1); printout(2); printout(5);endfunction



