function [IntersectionsLaserGlass,vektorsLaserGlass,...
    IntersectionsLaserWater,vectorsLaserWater,...
    IntersectionsLaserGround] = ...
    raytracingLaser(position,omega,phi,kappa,fov,numberBeams,...
    zGlass,zWater,zGround,nAir,nGlass,nWater)

%% Sardemann, IPF, TU Dresden, 2022

%% Interfaces are parallel to xy-plane
NormalGlass = [0;0;1];
PositionGlass = [0;0;zGlass];
NormalWater = [0;0;1];
PositionWater = [0;0;zWater];
NormalGround = [0;0;1];
PositionGround = [0;0;zGround];


%% Laser
% Orientation Laser:
RLaser = eul2rotm([kappa,phi,omega],'ZYX');

% Approximate Laser sheet with sub beams
s1Input = [zeros(1,numberBeams);-tan(linspace(-fov/2,fov/2,numberBeams)*pi/180);-ones(1,numberBeams)];
s1Norm = [s1Input(1,:)./vecnorm(s1Input);...
    s1Input(2,:)./vecnorm(s1Input);...
    s1Input(3,:)./vecnorm(s1Input)];
s1 = RLaser*s1Norm;

%% Refraction at interfaces:
for i = 1:numberBeams     % each Sub-beam
    % Interface air glass
    [IntersectionsLaserGlass(:,i),vektorsLaserGlass(:,i)] = refraction(position,s1(:,i),nAir,nGlass,PositionGlass,NormalGlass);
    % Interface Glass Water 
    [IntersectionsLaserWater(:,i),vectorsLaserWater(:,i)] = refraction(IntersectionsLaserGlass(:,i),vektorsLaserGlass(:,i),nGlass,nWater,PositionWater,NormalWater);
    % Intersection of sub beam with ground
    [IntersectionsLaserGround(:,i),~] = line_plane_intersection(vectorsLaserWater(:,i), IntersectionsLaserWater(:,i), NormalGround, PositionGround);
end

end

