function WaterDroplets()

    	% Include Libra SDK/include/matlab interface (language bindings)
	oldPath = cd(pwd);
    	cd .. 
    	cd .. 
    	cd ..
    	addpath([pwd '/include/matlab/']);
    
    	disp('*************************');
	disp('Water droplets simulation');
	disp('*************************');

    	% Point out /Libra SDK/ root
    	ret = libra_init(pwd, {'-c', 'ClientSDK'});

    	cd(oldPath);
    
    	if (ret > 0)
        	exit;
    	end

	run();

    	libra_shutdown();
end

function run()

    	% Set compute precision
    	libra_setDefaultDataType(GFLOAT32.value);

    	% Setup variables.
    	gridSizeX = 256;
    	gridSizeY = 256;
	
    	walpha = 1.0;
    	A = gVar();
    	x = gVar();
    	b = gVar();
    	g1 = gVar();g2 = gVar();g3 = gVar();g4 = gVar();

    	% Setup system
    	[walpha, A,x,b,g1,g2,g3,g4] = setupWaveEquationSystem(walpha, A, x, b, g1, g2, g3, g4, gridSizeX, gridSizeY);
    
    	% Solve system
    	solveRenderWaveEquationSystem(walpha, A, x, b, g1, g2, g3, g4, gridSizeX, gridSizeY);
    
end

function [alpha, A, x, b, g1,g2,g3,g4] = setupWaveEquationSystem(alpha, A, x, b, g1, g2, g3, g4, gridSizeX, gridSizeY)
    
	%
	% 2D Wave equation
	% 

 	C = 0.5;
 	h = 1.0;
 	dT = 1.0;
 	
 	vectorSize = gridSizeX*gridSizeY;
 	matrixSizeX = vectorSize;
 	matrixSizeY = vectorSize;
 
 	alpha = (dT*dT*C*C) / (2*h*h);
 	mainDiag = 4*alpha + 1;
     
 	data1 = ones(vectorSize,1);
 	data2 = ones(vectorSize,1);
 	data3 = ones(vectorSize,1);
 
 	% setup data main diagonal-Y
 	for i = 1 : vectorSize
 		if (i>gridSizeX)
        	    data1(i) = -alpha;
        	else
       	     	    data1(i) = 0;
        	end
	end

 	% setup data main diagonal-1
 	for i = 1 : vectorSize
		if (mod(i,gridSizeY))
        	    data2(i) = -alpha;
        	else
            	    data2(i) = 0;
	        end
    	end

 	% setup data main diagonal
 	for i = 1 : vectorSize
       		data3(i) = mainDiag;
    	end

	bandCount = 5;
	bandOffsetDatas = [-gridSizeX, -1, 0, 1, gridSizeX];
	bandValues = [data1 data2 data3 data2 data1];
    
    	% create the sparse matrix with input on band form
 	A = libra_sparseBanded(matrixSizeX, matrixSizeY, bandValues, bandOffsetDatas, bandCount);

	u = libra_mod(libra_floor(libra_step(0.0, vectorSize - 1.0, 1.0)), single(gridSizeX));
	v = libra_floor(libra_div(libra_step(0.0, vectorSize - 1.0, 1.0), single(gridSizeX)));
	g1 = libra_max(u - 1.0, 0.0) + v * single(gridSizeX);
 	g2 = libra_min(u + 1.0, gridSizeX - 1.0) + v * single(gridSizeX);
 	g3 = u + libra_max(v - 1.0, 0.0) * single(gridSizeX);
 	g4 = u + libra_min(v + 1.0, gridSizeX - 1.0) * single(gridSizeX);
 
    	% setup initial solution, x = 0
 	x = libra_zeros(vectorSize, 1);
 
    	%  create drops by scattering a few random values into x
 	dropCount = max(1,int32(gridSizeX  / 20.0));
  	libra_scatter(libra_rand(dropCount, 1) * single((libra_rowCount(x) - 1)), libra_rand(dropCount, 1)*0.5, x);
end

function [] = solveRenderWaveEquationSystem(walpha, A, x, b, g1, g2, g3, g4, gridSizeX, gridSizeY)

	M = libra_diag(A);
    	s = gVar(gridSizeX, gridSizeY);
 	u = libra_zeros(x.rowCount(), x.columnCount());

    	totalTime = 0.0; iterationTime = 0.0; averageIterationTime = 0.0;
 	timeSteps = 5000;
 	startTime = libra_getTime();
   	rowCount = libra_rowCount(x);

    	% for each step in time...
 	for i=0:timeSteps

        	if (i == 10)
 			startTime = libra_getTime();
        	end 
        
        	% Create new drops by scattering a random value
        	% into x at the random position
        	% Scatter heavy rain
        	libra_scatter(libra_rand(15, 1) * single((rowCount - 1.0)), libra_rand(15, 1)*0.2, x);
         
        	% Setup right hand side (b) of equation
        	b = walpha * (libra_gather(g1, x) + libra_gather(g2, x) + libra_gather(g3, x) + libra_gather(g4, x)) + (2.0 - 4.0 * walpha) * x - u;
        
        	% u <-- old solution
        	u = x;
  
        	% x <-- solve for new solution
        	for j=0:2
            		x = libra_jacobiIter(A, x, b);
        	end
   
        	% render solution x to screen
        	if (i > 80)
          		libra_render(x, bitor(GHEIGHTFIELD2D.value, GSTATS.value), gridSizeX);
        	end 
 	end

 	totalTime = libra_getTime() - startTime;
 	averageIterationTime = totalTime / (timeSteps-10);
 
 	disp('Total time (s) : ');totalTime
 	disp('Average iteration time (ms) : ');averageIterationTime * 1000.0   

end