function [B] = computeB(l, b, s, theta)
%COMPUTEB Summary of this function goes here
%   Detailed explanation goes here

n_meas = length(l); % number of measurements

B = zeros(n_meas, n_meas*3); % n x 3n

for i_row = 1:n_meas

    Vx = l(i_row, 1); Vy = l(i_row, 2); Vz = l(i_row, 3);
    theta_xy = theta(1); theta_xz = theta(2); theta_yz = theta(3);
    Sx = s(1); Sy = s(2); Sz = s(3);
    Bx = b(1); By = b(2); Bz = b(3);

    term_1 = (2*(sin(theta_xz)/Sx - (cos(theta_xz)*sin(theta_xy)*sin(theta_yz))/(Sx*cos(theta_xy)))*((Bz - Vz)/Sz - (sin(theta_xz)*(Bx - Vx))/Sx + (cos(theta_xz)*sin(theta_yz)*(Vy - By + (Sy*sin(theta_xy)*(Bx - Vx))/Sx))/(Sy*cos(theta_xy))))/(cos(theta_xz)^2*cos(theta_yz)^2) - (2*Bx - 2*Vx)/Sx^2 - (2*sin(theta_xy)*(Vy - By + (Sy*sin(theta_xy)*(Bx - Vx))/Sx))/(Sx*Sy*cos(theta_xy)^2);
    term_2 = (2*Vy - 2*By + (2*Sy*sin(theta_xy)*(Bx - Vx))/Sx)/(Sy^2*cos(theta_xy)^2) + (2*sin(theta_yz)*((Bz - Vz)/Sz - (sin(theta_xz)*(Bx - Vx))/Sx + (cos(theta_xz)*sin(theta_yz)*(Vy - By + (Sy*sin(theta_xy)*(Bx - Vx))/Sx))/(Sy*cos(theta_xy))))/(Sy*cos(theta_xy)*cos(theta_xz)*cos(theta_yz)^2);
    term_3 = -(2*((Bz - Vz)/Sz - (sin(theta_xz)*(Bx - Vx))/Sx + (cos(theta_xz)*sin(theta_yz)*(Vy - By + (Sy*sin(theta_xy)*(Bx - Vx))/Sx))/(Sy*cos(theta_xy))))/(Sz*cos(theta_xz)^2*cos(theta_yz)^2);

    B(i_row, ((i_row-1)*3)+1:i_row*3) = [term_1 term_2 term_3];

end

end