function[val] = bilin_interpol(phi_mat,x_vec,y_vec,x,y)

if x < x_vec(1)
    x = x_vec(1);
end
if x > x_vec(end)
    x = x_vec(end);
end
if y < y_vec(1)
    y = y_vec(1);
end
if y > y_vec(end)
    y = y_vec(end);
end


%%%%%

n_y = length(y_vec);
n_x = length(x_vec);

for i = 2:n_x
    for j = 2:n_y
       
        
        if (x >= x_vec(i-1)) && (x <= x_vec(i)) && (y >= y_vec(j-1)) && (y <= y_vec(j))
        
        val = (  (x_vec(i)-x)*(y_vec(j)-y)*phi_mat(i-1,j-1) +  (x_vec(i)-x)*(-y_vec(j-1)+y)*phi_mat(i-1,j) ...
            +(-x_vec(i-1)+x)*(y_vec(j)-y)*phi_mat(i,j-1)+ (-x_vec(i-1)+x)*(-y_vec(j-1)+y)*phi_mat(i,j)  )/(  (x_vec(i)-x_vec(i-1) )*( y_vec(j)-y_vec(j-1) )   );
        
        end
        
        
    end
end