In the last couple of decades, Arctic Engineering has become a topic of interest. There are still plenty rooms of research to understand the unique characteristic of sea ice, especially in relation to offshore engineering. This includes the most fundamental problem in floating body motion analysis: the radiation-diffraction problem. A powerful mathematical concept, - so called Greens Function – is one of the favourable tools to be used for solving this mathematical problem. This is because the radiation-diffraction problem can be formulated as a boundary value problem expressed by partial differential equations. Although the application is already quite advanced for the open water case, the same cannot be said for the vessel operating in ice infested waters. The integral solution of 3D Greens Function for ice-infested waters which has not been studied before, was derived in this thesis. The open water case is also studied to gain more insight in the implementation of an arbitrary floating body thoroughly. As a closure, interpretation about the difference between open and ice-infested waters is discussed. For Greens Function in the open water case, numerical evaluation of the principal value integral is not straightforward due to the hyperbolic term inside the integrand. This term makes the integrand exceed the limit of floating point number (in MATLAB) and cannot be evaluated into infinity. On the other hand, a numerical integration is quite time consuming (whereas the analytical solution, as far as the writer’s concern, is not found yet). A well-known alternative form of the solution which formulated as an infinite series might improve the computation speed. The rate of convergence depends solely on the ratio of horizontal distance between source and field point (R) and water depth (H). Another numerical issue arises in the deep water case. A finite water depth causes a catastrophic cancellation, both for the integral and the infinite series solution, due to the extremely small difference of the wave number between deep water and infinite water depth. This is where the infinite depth solution needs to be used. In numerical implementation, an influence function at a panel can be approximated by multiplying the potential with panel's surface area. Due to the singularity of the Rankine source term in the integral or the modified Bessel function in the series, this approach fails and the solution need to be integrated over the panel. Although the analytical solution is available, a numerical approach is chosen to simplify the problem. The surface integration procedure can be done by transforming the global arbitrary panel orientation into the local element coordinate system, and subsequently, perform a bilinear mapping to reshape the quadrilateral panel into the desired rectangular panel. This transformation procedure is needed since MATLAB is only capable to handle double numerical integration of a function bounded between four lines perpendicularly each other. This encloses the whole challenge that needs to be addressed and might be re-occur in the ice-infested waters as well. From the derivation of the Greens Function for ice-infested waters, it is shown that the hyperbolic term inside the integrand is present. This discloses that the obtained solutions cannot be used instantaneously. Another effort to rewrite them in the exponential term might be useful. Moreover, the radiation condition is not satisfied yet in this thesis. However, the suggested approach of the derivation by introducing an imaginary line to represent the source depth location avoids the use of a singular term. Generally speaking, this thesis initiates a promising foundation for further research in the hydromechanics analysis of sea ice.