diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml
new file mode 100644
index 0000000..a672d2e
--- /dev/null
+++ b/.JuliaFormatter.toml
@@ -0,0 +1,7 @@
+remove_extra_newlines=true
+join_lines_based_on_source=true
+whitespace_in_kwargs=false
+short_to_long_function_def=true
+always_for_in=true
+verbose=true
+margin=88
\ No newline at end of file
diff --git a/.github/workflows/Test.yml b/.github/workflows/Test.yml
new file mode 100644
index 0000000..e448471
--- /dev/null
+++ b/.github/workflows/Test.yml
@@ -0,0 +1,30 @@
+name: Run tests
+
+on: 
+  push:
+    branches:
+      - main
+  pull_request:
+    branches:
+      - main
+      - dev
+
+jobs:
+  test:
+    runs-on: ${{ matrix.os }}
+    strategy:
+      matrix:
+        julia-version: ['1']
+        julia-arch: [x64]
+        os: [ubuntu-latest]
+
+    steps:
+      - uses: actions/checkout@v2
+      - uses: julia-actions/setup-julia@v1.9.2
+        with:
+          version: ${{ matrix.julia-version }}
+          arch: ${{ matrix.julia-arch }}
+      - uses: julia-actions/julia-buildpkg@v1
+      - uses: julia-actions/julia-runtest@v1
+        with:
+          annotate: true
\ No newline at end of file
diff --git a/Project.toml b/Project.toml
index 3b6f276..81bcc0e 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,12 +1,10 @@
 name = "ShipMMG"
 uuid = "37f2b0bf-0c13-4883-8808-e75eb56597e7"
 authors = ["Taiga MITSUYUKI <mitsuyuki-taiga-my@ynu.ac.jp>"]
-version = "0.0.6"
+version = "0.0.7"
 
 [deps]
 AdvancedPS = "576499cb-2369-40b2-a588-c64705576edc"
-CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
-DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Dierckx = "39dd38d3-220a-591b-8e3c-4c3a8c710a94"
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
@@ -15,10 +13,6 @@ KernelDensity = "5ab0869b-81aa-558d-bb23-cbf5423bbe9b"
 Libtask = "6f1fad26-d15e-5dc8-ae53-837a1d7b8c9f"
 ParameterizedFunctions = "65888b18-ceab-5e60-b2b9-181511a3b968"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
-Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
-Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
-Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
-Serialization = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
 Turing = "fce5fe82-541a-59a6-adf8-730c64b5f9a0"
 
 [compat]
@@ -29,7 +23,6 @@ Distributions = "^0.25"
 Libtask = "^0.8"
 ParameterizedFunctions = "^5.12"
 Parameters = "0.12"
-Plots = "^1"
 Turing = "^0.24"
 julia = "^1.6"
 
diff --git a/src/ShipMMG.jl b/src/ShipMMG.jl
index 0a12406..b679342 100644
--- a/src/ShipMMG.jl
+++ b/src/ShipMMG.jl
@@ -3,7 +3,6 @@ module ShipMMG
 using DifferentialEquations
 using ParameterizedFunctions
 using Dierckx
-using Plots
 using Parameters
 using Distributions
 using Turing
@@ -15,6 +14,7 @@ export calc_position,
     get_KVLCC2_L7_basic_params,
     get_KVLCC2_L7_maneuvering_params,
     get_KVLCC2_L7_params,
+    get_example_ship_wind_force_moment_params,
     nuts_sampling_single_thread,
     nuts_sampling_multi_threads
 
@@ -36,9 +36,7 @@ export Mmg3DofBasicParams,
     mmg_3dof_zigzag_test,
     estimate_mmg_approx_lsm,
     estimate_mmg_approx_lsm_time_window_sampling,
-    create_model_for_mcmc_sample_mmg
-
-include("draw.jl")
-export draw_gif_result, calc_position
+    create_model_for_mcmc_sample_mmg,
+    wind_force_and_moment_coefficients
 
 end
diff --git a/src/data.jl b/src/data.jl
index 7ade86c..4bcd28b 100644
--- a/src/data.jl
+++ b/src/data.jl
@@ -12,7 +12,7 @@
     ψ_W::Tψ_w
 end
 
-function get_KVLCC2_L7_basic_params(ρ = 1025.0)
+function get_KVLCC2_L7_basic_params(ρ=1025.0)
     L_pp = 7.00  # 船長Lpp[m]
     B = 1.27  # 船幅[m]
     d = 0.46  # 喫水[m]
@@ -146,7 +146,126 @@ function get_KVLCC2_L7_params()
     basic_params, maneuvering_params
 end
 
-function calc_position(time_vec, u_vec, v_vec, r_vec; x0 = 0.0, y0 = 0.0, ψ0 = 0.0)
+"""
+    wind_force_and_moment_coefficients(ψ_A,p)
+
+compute the parameter C_X,C_Y,C_N from Fujiwara's estimation method (1994).
+
+# Arguments
+-`ψ_A`: Wind attack of angle [rad]
+- L_pp
+- B
+- A_OD
+- A_F
+- A_L
+- H_BR
+- H_C
+- C
+"""
+function wind_force_and_moment_coefficients(
+    ψ_A,
+    L_pp,
+    B,
+    A_OD,
+    A_F,
+    A_L,
+    H_BR,
+    H_C,
+    C,
+)
+    #C_LF1の場合で調整
+    C_CF = 0.404 + 0.368 * A_F / (B * H_BR) + 0.902 * H_BR / L_pp
+
+    if deg2rad(0) <= ψ_A <= deg2rad(90)
+        C_LF = -0.992 + 0.507 * A_L / (L_pp * B) + 1.162 * C / L_pp
+        C_XLI = 0.458 + 3.245 * A_L / (L_pp * H_BR) - 2.313 * A_F / (B * H_BR)
+        C_ALF = -0.585 - 0.906 * A_OD / A_L + 3.239 * B / L_pp
+        C_YLI = pi * A_L / L_pp^2 + 0.116 + 3.345 * A_F / (L_pp * B)
+
+    elseif deg2rad(90) < ψ_A <= deg2rad(180)
+        C_LF =
+            0.018 - 5.091 * B / L_pp + 10.367 * H_C / L_pp - 3.011 * A_OD / L_pp^2 -
+            0.341 * A_F / B^2
+        C_XLI =
+            -1.901 + 12.727 * A_L / (L_pp * H_BR) + 24.407 * A_F / A_L -
+            40.310 * B / L_pp - 0.341 * A_F / (B * H_BR)
+        C_ALF = -0.314 - 1.117 * A_OD / A_L
+        C_YLI = pi * A_L / L_pp^2 + 0.446 + 2.192 * A_F / L_pp^2
+
+    elseif deg2rad(180) < ψ_A <= deg2rad(270)
+        C_LF =
+            0.018 - 5.091 * B / L_pp + 10.367 * H_C / L_pp - 3.011 * A_OD / L_pp^2 -
+            0.341 * A_F / B^2
+        C_XLI =
+            -1.901 + 12.727 * A_L / (L_pp * H_BR) + 24.407 * A_F / A_L -
+            40.310 * B / L_pp - 0.341 * A_F / (B * H_BR)
+        C_ALF = -(-0.314 - 1.117 * A_OD / A_L)
+        C_YLI = -(pi * A_L / L_pp^2 + 0.446 + 2.192 * A_F / L_pp^2)
+
+    elseif deg2rad(270) < ψ_A <= deg2rad(360)
+        C_LF = -0.992 + 0.507 * A_L / (L_pp * B) + 1.162 * C / L_pp
+        C_XLI = 0.458 + 3.245 * A_L / (L_pp * H_BR) - 2.313 * A_F / (B * H_BR)
+        C_ALF = -(-0.585 - 0.906 * A_OD / A_L + 3.239 * B / L_pp)
+        C_YLI = -(pi * A_L / L_pp^2 + 0.116 + 3.345 * A_F / (L_pp * B))
+    end
+
+    C_X =
+        C_LF * cos(ψ_A) +
+        C_XLI * (sin(ψ_A) - sin(ψ_A) * cos(ψ_A)^2 / 2) * sin(ψ_A) * cos(ψ_A) +
+        C_ALF * sin(ψ_A) * cos(ψ_A)^3
+    C_Y =
+        C_CF * sin(ψ_A)^2 +
+        C_YLI * (cos(ψ_A) + sin(ψ_A)^2 * cos(ψ_A) / 2) * sin(ψ_A) * cos(ψ_A)
+    C_N = C_Y * (0.297 * C / L_pp - 0.149 * (ψ_A - deg2rad(90)))
+
+    C_X, C_Y, C_N
+end
+
+function get_example_ship_wind_force_moment_params()
+    L_pp = 7.00  # 船長Lpp[m]
+    B = 1.27  # 船幅[m]
+    d = 0.46  # 喫水[m]
+    D = 0.6563 # 深さ[m]
+    A_OD = 0.65 # デッキ上の構造物の側面投影面積[m^2]
+    H_BR = 0.85 # 喫水からブリッジ主要構造物の最高位[m]
+    H_C = 0.235 # 喫水から側面積中心までの高さ[m]
+    C = 0.0 # 船体中心から側面積中心までの前後方向座標(船首方向を正)[m]
+
+    A_OD = A_OD # デッキ上の構造物の側面投影面積[m^2]
+    A_F = (D - d) * B  # 船体の正面投影面積[m^2]
+    A_L = (D - d) * L_pp # 船体の側面投影面積[m^2]
+    H_BR = H_BR # 喫水からブリッジ主要構造物の最高位[m]
+    H_C = H_C # 喫水から側面積中心までの高さ[m]
+    C = C # 船体中心から側面積中心までの前後方向座標[m]
+
+    ψ_A_vec = deg2rad.(collect(0:10:360))
+    C_X_vec = Array{Float64}(undef, length(ψ_A_vec))
+    C_Y_vec = Array{Float64}(undef, length(ψ_A_vec))
+    C_N_vec = Array{Float64}(undef, length(ψ_A_vec))
+    for (index, ψ_A) in enumerate(ψ_A_vec)
+        C_X, C_Y, C_N = wind_force_and_moment_coefficients(
+            ψ_A,
+            L_pp,
+            B,
+            A_OD,
+            A_F,
+            A_L,
+            H_BR,
+            H_C,
+            C,
+        )
+        C_X_vec[index] = C_X
+        C_Y_vec[index] = C_Y
+        C_N_vec[index] = C_N
+    end
+    spl_C_X = Spline1D(ψ_A_vec, C_X_vec)
+    spl_C_Y = Spline1D(ψ_A_vec, C_Y_vec)
+    spl_C_N = Spline1D(ψ_A_vec, C_N_vec)
+
+    Mmg3DofWindForceMomentParams(A_F, A_L, spl_C_X, spl_C_Y, spl_C_N)
+end
+
+function calc_position(time_vec, u_vec, v_vec, r_vec; x0=0.0, y0=0.0, ψ0=0.0)
     dim = length(time_vec)
     x_vec = zeros(Float64, dim)
     y_vec = zeros(Float64, dim)
@@ -154,11 +273,13 @@ function calc_position(time_vec, u_vec, v_vec, r_vec; x0 = 0.0, y0 = 0.0, ψ0 =
     x_vec[1] = x0
     y_vec[1] = y0
     ψ_vec[1] = ψ0
-    for i = 2:dim
+    for i in 2:dim
         Δt = time_vec[i] - time_vec[i-1]
         ψ_vec[i] = ψ_vec[i-1] + r_vec[i] * Δt
-        x_vec[i] = x_vec[i-1] + (u_vec[i] * cos(ψ_vec[i]) - v_vec[i] * sin(ψ_vec[i])) * Δt
-        y_vec[i] = y_vec[i-1] + (u_vec[i] * sin(ψ_vec[i]) + v_vec[i] * cos(ψ_vec[i])) * Δt
+        x_vec[i] =
+            x_vec[i-1] + (u_vec[i] * cos(ψ_vec[i]) - v_vec[i] * sin(ψ_vec[i])) * Δt
+        y_vec[i] =
+            y_vec[i-1] + (u_vec[i] * sin(ψ_vec[i]) + v_vec[i] * cos(ψ_vec[i])) * Δt
     end
     x_vec, y_vec, ψ_vec
 end
@@ -167,12 +288,12 @@ function nuts_sampling_single_thread(
     model,
     n_samples::Int,
     n_chains::Int;
-    target_acceptance::Float64 = 0.65,
-    progress = true,
+    target_acceptance::Float64=0.65,
+    progress=true,
 )
     sampler = NUTS(target_acceptance)
     mapreduce(
-        c -> sample(model, sampler, n_samples, progress = progress),
+        c -> sample(model, sampler, n_samples, progress=progress),
         chainscat,
         1:n_chains,
     )
@@ -182,9 +303,9 @@ function nuts_sampling_multi_threads(
     model,
     n_samples::Int,
     n_chains::Int;
-    target_acceptance::Float64 = 0.65,
-    progress = false,
+    target_acceptance::Float64=0.65,
+    progress=false,
 )
     sampler = NUTS(target_acceptance)
-    sample(model, sampler, MCMCThreads(), n_samples, n_chains, progress = progress)
+    sample(model, sampler, MCMCThreads(), n_samples, n_chains, progress=progress)
 end
diff --git a/src/draw.jl b/src/draw.jl
deleted file mode 100644
index c2845ca..0000000
--- a/src/draw.jl
+++ /dev/null
@@ -1,44 +0,0 @@
-"""
-    draw_gif_result(time, x, y, ψ, shape, file_path, [, interval, fps]) -> gif
-
-Draw the gif animation from simulation result.
-"""
-function draw_gif_result(time, x, y, ψ, shape, file_path; interval = 10, fps = 100)
-    anim = @animate for i = 1:interval:length(time)
-        plot(
-            x,
-            y,
-            label = "",
-            xlabel = "x [m]",
-            ylabel = "y [m]",
-            linestyle = :dot,
-            aspect_ratio = :equal,
-        )
-        ship = square(x[i], y[i], shape, ψ[i])
-        plot!(Shape(ship[1], ship[2]), label = "")
-        scatter!(
-            [x[i]],
-            [y[i]],
-            seriestype = :scatter,
-            title = "time = $(time[i])",
-            label = "",
-        )
-    end
-    gif(anim, file_path, fps = fps)
-end
-
-function rotate_pos(pos, angle)
-    rotate_matrix = [cos(angle) -sin(angle); sin(angle) cos(angle)]
-    rotate_matrix * pos
-end
-
-function square(center_x, center_y, shape, angle)
-    square_xy = [
-        rotate_pos([shape[1], shape[2]] / 2, angle) + [center_x, center_y],
-        rotate_pos([-shape[1], shape[2]] / 2, angle) + [center_x, center_y],
-        rotate_pos([-shape[1], -shape[2]] / 2, angle) + [center_x, center_y],
-        rotate_pos([shape[1], -shape[2]] / 2, angle) + [center_x, center_y],
-    ]
-    xy = hcat(square_xy...)
-    [xy[1, :], xy[2, :]]
-end
diff --git a/src/kt.jl b/src/kt.jl
index b898d2f..e90050f 100644
--- a/src/kt.jl
+++ b/src/kt.jl
@@ -77,23 +77,22 @@ function kt_simulate(
     T,
     time_list,
     δ_list;
-    u0 = 10.0,
-    v0 = 0.0,
-    r0 = 0.0,
-    x0 = 0.0,
-    y0 = 0.0,
-    Ψ0 = 0.0,
-    algorithm = Tsit5(),
-    reltol = 1e-8,
-    abstol = 1e-8,
+    u0=10.0,
+    v0=0.0,
+    r0=0.0,
+    x0=0.0,
+    y0=0.0,
+    Ψ0=0.0,
+    algorithm=Tsit5(),
+    reltol=1e-8,
+    abstol=1e-8,
 )
-
     spl_δ = Spline1D(time_list, δ_list)
 
     X0 = [u0; v0; r0; x0; y0; Ψ0; δ_list[1]]
     p = [K, T, spl_δ]
     prob = ODEProblem(kt_model!, X0, (time_list[1], time_list[end]), p)
-    sol = solve(prob, algorithm, reltol = reltol, abstol = abstol)
+    sol = solve(prob, algorithm, reltol=reltol, abstol=abstol)
     sol_timelist = sol(time_list)
     results = hcat(sol_timelist.u...)
     u = results[1, :]
@@ -158,17 +157,17 @@ function kt_zigzag_test(
     time_list,
     target_δ_rad,
     target_Ψ_rad_deviation;
-    u0 = 10.0,
-    v0 = 0.0,
-    r0 = 0.0,
-    x0 = 0.0,
-    y0 = 0.0,
-    Ψ0 = 0.0,
-    δ0 = 0.0,
-    δ_rad_rate = 10.0 * π / 180,
-    algorithm = Tsit5(),
-    reltol = 1e-8,
-    abstol = 1e-8,
+    u0=10.0,
+    v0=0.0,
+    r0=0.0,
+    x0=0.0,
+    y0=0.0,
+    Ψ0=0.0,
+    δ0=0.0,
+    δ_rad_rate=10.0 * π / 180,
+    algorithm=Tsit5(),
+    reltol=1e-8,
+    abstol=1e-8,
 )
     target_Ψ_rad_deviation = abs(target_Ψ_rad_deviation)
 
@@ -205,7 +204,7 @@ function kt_zigzag_test(
             y0 = final_y_list[start_index-1]
         end
 
-        for i = (start_index+1):length(time_list)
+        for i in (start_index+1):length(time_list)
             Δt = time_list[i] - time_list[i-1]
             if target_δ_rad > 0
                 δ = δ_list[i-start_index] + δ_rad_rate * Δt
@@ -227,15 +226,15 @@ function kt_zigzag_test(
             T,
             time_list[start_index:end],
             δ_list,
-            u0 = u0,
-            v0 = v0,
-            r0 = r0,
-            x0 = x0,
-            y0 = y0,
-            Ψ0 = Ψ,
-            algorithm = algorithm,
-            reltol = reltol,
-            abstol = abstol,
+            u0=u0,
+            v0=v0,
+            r0=r0,
+            x0=x0,
+            y0=y0,
+            Ψ0=Ψ,
+            algorithm=algorithm,
+            reltol=reltol,
+            abstol=abstol,
         )
 
         # get finish index
@@ -306,7 +305,7 @@ function estimate_kt_lsm_time_window_sampling(data::ShipData, window_size::Int)
     n_samples = length(time_vec) - window_size
     K_samples = zeros(n_samples)
     T_samples = zeros(n_samples)
-    for i = 1:n_samples
+    for i in 1:n_samples
         K_samples[i], T_samples[i] =
             estimate_kt_lsm(r[i:i+window_size], dr[i:i+window_size], δ[i:i+window_size])
     end
@@ -315,9 +314,9 @@ end
 
 function create_model_for_mcmc_sample_kt(
     data::ShipData;
-    σ_r_prior_dist::Distribution = Chi(5),
-    K_prior_dist::Distribution = Uniform(0.01, 10.0),
-    T_prior_dist::Distribution = truncated(Normal(100.0, 50.0), 10.0, 200.0),
+    σ_r_prior_dist::Distribution=Chi(5),
+    K_prior_dist::Distribution=Uniform(0.01, 10.0),
+    T_prior_dist::Distribution=truncated(Normal(100.0, 50.0), 10.0, 200.0),
 )
     time_obs = data.time
     r_obs = data.r
@@ -338,16 +337,16 @@ function create_model_for_mcmc_sample_kt(
     prob1 = ODEProblem(KT!, u0, (time_obs[1], time_obs[end]), p)
 
     # create probabilistic model
-    @model function fitKT(time_obs, r_obs, prob1)
+    Turing.@model function fitKT(time_obs, r_obs, prob1)
         σ_r ~ σ_r_prior_dist
         K ~ K_prior_dist
         T ~ T_prior_dist
 
         p = [K, T]
-        prob = remake(prob1, p = p)
+        prob = remake(prob1, p=p)
         sol = solve(prob, Tsit5())
         predicted = sol(time_obs)
-        for i = 1:length(predicted)
+        for i in 1:length(predicted)
             r_obs[i] ~ Normal(predicted[i][1], σ_r) # index number of r is 1
         end
     end
diff --git a/src/mmg.jl b/src/mmg.jl
index 66d6f67..2781f84 100644
--- a/src/mmg.jl
+++ b/src/mmg.jl
@@ -72,12 +72,19 @@ MMG 3DOF model on DifferentialEquations.ODEProblem. Update `dX`.
     - N_vvr_dash
     - N_vrr_dash
     - N_rrr_dash
+    - A_F
+    - A_L
+    - spl_C_X
+    - spl_C_Y
+    - spl_C_N
     - spl_δ
     - spl_n_p
+    - u_wind_list
+    - ψ_wind_list
 - `t`: the time.
 """
 function mmg_3dof_model!(dX, X, p, t)
-    u, v, r, x, y, Ψ, δ, n_p = X
+    u, v, r, x, y, ψ, δ, n_p, u_wind, ψ_wind = X
     ρ,
     L_pp,
     B,
@@ -261,6 +268,8 @@ function mmg_3dof_model!(dX, X, p, t)
     dX[6] = dΨ = r
     dX[7] = dδ = derivative(spl_δ, t)
     dX[8] = dn_p = derivative(spl_n_p, t)
+    dX[9] = du_wind = derivative(spl_u_wind, t)
+    dX[10] = dψ_wind = derivative(spl_ψ_wind, t)
 end
 
 """
@@ -410,7 +419,7 @@ This function has the same logic of `ShipMMG.simulate()`.
 - `r0=0.0`: the initial rate of turn [rad/s].
 - `x0=0.0`: the initial x (surge) position.
 - `y0=0.0`: the initial y (sway) position.
-- `Ψ0=0.0`: the initial Ψ (yaw) azimuth [rad].
+- `ψ0=0.0`: the initial ψ (yaw) azimuth [rad].
 - `ρ=1025.0`: the seawater density [kg/m^3].
 - `algorithm=Tsit5()`: the parameter of DifferentialEquations.ODEProblem.solve()
 - `reltol=1e-8`: the parameter of DifferentialEquations.ODEProblem.solve()
@@ -458,11 +467,11 @@ function mmg_3dof_simulate(
     r0=0.0,
     x0=0.0,
     y0=0.0,
-    Ψ0=0.0,
+    ψ0=0.0,
     ρ=1025.0,
     algorithm=Tsit5(),
     reltol=1e-8,
-    abstol=1e-8
+    abstol=1e-8,
 )
     @unpack L_pp,
     B,
@@ -574,7 +583,7 @@ function mmg_3dof_simulate(
         r0=r0,
         x0=x0,
         y0=y0,
-        Ψ0=Ψ0,
+        ψ0=ψ0,
         ρ=ρ,
         algorithm=algorithm,
         reltol=reltol,
@@ -649,7 +658,7 @@ This function has the same logic of `ShipMMG.mmg_3dof_simulate()`.
 - `r0=0.0`: the initial rate of turn [rad/s].
 - `x0=0.0`: the initial x (surge) position.
 - `y0=0.0`: the initial y (sway) position.
-- `Ψ0=0.0`: the initial Ψ (yaw) azimuth [rad].
+- `ψ0=0.0`: the initial ψ (yaw) azimuth [rad].
 - `ρ=1025.0`: the seawater density [kg/m^3].
 - `algorithm=Tsit5()`: the parameter of DifferentialEquations.ODEProblem.solve()
 - `reltol=1e-8`: the parameter of DifferentialEquations.ODEProblem.solve()
@@ -716,7 +725,7 @@ function simulate(
     r0=0.0,
     x0=0.0,
     y0=0.0,
-    Ψ0=0.0,
+    ψ0=0.0,
     ρ=1025.0,
     algorithm=Tsit5(),
     reltol=1e-8,
@@ -794,21 +803,25 @@ function simulate(
     r = results[3, :]
     x = results[4, :]
     y = results[5, :]
-    Ψ = results[6, :]
+    ψ = results[6, :]
     δ = results[7, :]
     n_p = results[8, :]
-    u, v, r, x, y, Ψ, δ, n_p
+    u_wind = results[9, :]
+    ψ_wind = results[10, :]
+    u, v, r, x, y, ψ, δ, n_p, u_wind, ψ_wind
 end
 
 """
-    mmg_3dof_zigzag_test(basic_params, maneuvering_params, time_list, n_p_list, target_δ_rad, target_Ψ_rad_deviation, [, u0, v0, r0, x0, y0, Ψ0, δ0, δ_rad_rate, algorithm, reltol, abstol]) -> u, v, r, x, y, Ψ, δ
+    mmg_3dof_zigzag_test(basic_params, maneuvering_params,wind_force_and_moment_params time_list, n_p_list, u_wind_list, Ψ_wind_list, target_δ_rad, target_Ψ_rad_deviation, [, u0, v0, r0, x0, y0, Ψ0, δ0, δ_rad_rate, algorithm, reltol, abstol]) -> u, v, r, x, y, ψ, δ
 
 Returns the MMG 3DOF zigzag simulation results.
 
 # Arguments
 - `basic_params::Mmg3DofBasicParams`: the basic parameters of target ship.
 - `maneuvering_params::Mmg3DofManeuveringParams`: the maneuvering parameters of target ship.
+- `wind_force_and_moment_params::Mmg3DofWindForceMomentParams,` : the Structur parameters above the  taeget ship's draft
 - `time_list`: the list of simulatino time.
+- `δ_list`: the list of rudder angle [rad].
 - `n_p_list`: the list of propeller rps.
 - `target_δ_rad`: target rudder angle of zigzag test.
 - `target_Ψ_rad_deviation`: target azimuth deviation of zigzag test.
@@ -820,7 +833,7 @@ Returns the MMG 3DOF zigzag simulation results.
 - `r0=0.0`: the initial rate of turn [rad/s].
 - `x0=0.0`: the initial x (surge) position.
 - `y0=0.0`: the initial y (sway) position.
-- `Ψ0=0.0`: the initial Ψ (yaw) azimuth [rad].
+- `ψ0=0.0`: the initial ψ (yaw) azimuth [rad].
 - `δ0=0.0`: the initial rudder angle.
 - `δ_rad_rate=10.0*π/180`: the change rate of rudder angle [rad/s]. 
 - `ρ=1025.0`: the seawater density [kg/m^3].
@@ -835,7 +848,7 @@ KVLCC2_L7 zigzag test.
 julia> ρ = 1025.0;
 julia> basic_params, maneuvering_params = get_KVLCC2_L7_params();
 julia> target_δ_rad = 20.0 * π / 180.0
-julia> target_Ψ_rad_deviation = 20.0 * π / 180.0
+julia> target_ψ_rad_deviation = 20.0 * π / 180.0
 julia> start_time_second = 0.00
 julia> time_second_interval = 0.01
 julia> end_time_second = 80.00
@@ -848,10 +861,13 @@ julia> N_EX_list = zeros(Float64, length(time_list))
 julia> δ_list, u_list, v_list, r_list, Ψ_list = mmg_3dof_zigzag_test(
     basic_params,
     maneuvering_params,
-    time_list
+    wind_force_and_moment_params,
+    time_list,
     n_p_list,
+    u_wind_list,
+    ψ_wind_list,
     target_δ_rad,
-    target_Ψ_rad_deviation,
+    target_ψ_rad_deviation,
 );
 ```
 """
@@ -870,7 +886,7 @@ function mmg_3dof_zigzag_test(
     r0=0.0,
     x0=0.0,
     y0=0.0,
-    Ψ0=0.0,
+    ψ0=0.0,
     δ0=0.0,
     δ_rad_rate=10.0 * π / 180,
     ρ=1025.0,
@@ -878,7 +894,7 @@ function mmg_3dof_zigzag_test(
     reltol=1e-8,
     abstol=1e-8,
 )
-    target_Ψ_rad_deviation = abs(target_Ψ_rad_deviation)
+    target_ψ_rad_deviation = abs(target_ψ_rad_deviation)
 
     final_δ_list = zeros(length(time_list))
     final_u_list = zeros(length(time_list))
@@ -886,11 +902,11 @@ function mmg_3dof_zigzag_test(
     final_r_list = zeros(length(time_list))
     final_x_list = zeros(length(time_list))
     final_y_list = zeros(length(time_list))
-    final_Ψ_list = zeros(length(time_list))
+    final_ψ_list = zeros(length(time_list))
 
     next_stage_index = 1
     target_δ_rad = -target_δ_rad  # for changing in while loop
-    Ψ = Ψ0
+    ψ = ψ0
     while next_stage_index < length(time_list)
         target_δ_rad = -target_δ_rad
         start_index = next_stage_index
@@ -944,7 +960,7 @@ function mmg_3dof_zigzag_test(
             r0=r0,
             x0=x0,
             y0=y0,
-            Ψ0=Ψ,
+            ψ0=ψ,
             ρ=ρ,
             algorithm=algorithm,
             reltol=reltol,
@@ -952,13 +968,13 @@ function mmg_3dof_zigzag_test(
         )
 
         # get finish index
-        target_Ψ_rad = Ψ0 + target_Ψ_rad_deviation
+        target_ψ_rad = ψ0 + target_ψ_rad_deviation
         if target_δ_rad < 0
-            target_Ψ_rad = Ψ0 - target_Ψ_rad_deviation
+            target_ψ_rad = ψ0 - target_ψ_rad_deviation
         end
-        over_index = findfirst(e -> e > target_Ψ_rad, Ψ_list)
+        over_index = findfirst(e -> e > target_ψ_rad, ψ_list)
         if target_δ_rad < 0
-            over_index = findfirst(e -> e < target_Ψ_rad, Ψ_list)
+            over_index = findfirst(e -> e < target_ψ_rad, ψ_list)
         end
         next_stage_index = length(time_list)
         if isnothing(over_index)
@@ -968,9 +984,9 @@ function mmg_3dof_zigzag_test(
             final_r_list[start_index:next_stage_index] = r
             final_x_list[start_index:next_stage_index] = x
             final_y_list[start_index:next_stage_index] = y
-            final_Ψ_list[start_index:next_stage_index] = Ψ_list
+            final_ψ_list[start_index:next_stage_index] = ψ_list
         else
-            Ψ = Ψ_list[over_index]
+            ψ = ψ_list[over_index]
             next_stage_index = over_index + start_index - 1
             final_δ_list[start_index:next_stage_index-1] = δ_list[begin:over_index-1]
             final_u_list[start_index:next_stage_index-1] = u[begin:over_index-1]
@@ -978,7 +994,7 @@ function mmg_3dof_zigzag_test(
             final_r_list[start_index:next_stage_index-1] = r[begin:over_index-1]
             final_x_list[start_index:next_stage_index-1] = x[begin:over_index-1]
             final_y_list[start_index:next_stage_index-1] = y[begin:over_index-1]
-            final_Ψ_list[start_index:next_stage_index-1] = Ψ_list[begin:over_index-1]
+            final_ψ_list[start_index:next_stage_index-1] = ψ_list[begin:over_index-1]
         end
     end
     final_u_list,
@@ -986,7 +1002,7 @@ function mmg_3dof_zigzag_test(
     final_r_list,
     final_x_list,
     final_y_list,
-    final_Ψ_list,
+    final_ψ_list,
     final_δ_list
 end
 
@@ -1020,7 +1036,7 @@ function create_model_for_mcmc_sample_mmg(
     N_rrr_dash_prior_dist=Uniform(-0.060, 0.000),
     solver=Tsit5(),
     abstol=1e-6,
-    reltol=1e-3
+    reltol=1e-3,
 )
     time_obs = data.time
     u_obs = data.u
@@ -1137,7 +1153,10 @@ function create_model_for_mcmc_sample_mmg(
                 u *
                 (1.0 - w_P) *
                 ϵ *
-                sqrt(η * (1.0 + κ * (sqrt(1.0 + 8.0 * K_T / (pi * J^2)) - 1))^2 + (1 - η))
+                sqrt(
+                    η * (1.0 + κ * (sqrt(1.0 + 8.0 * K_T / (pi * J^2)) - 1))^2 +
+                    (1 - η),
+                )
         end
 
         U_R = sqrt(u_R^2 + v_R^2)
@@ -1259,7 +1278,7 @@ function create_model_for_mcmc_sample_mmg(
     prob1 = ODEProblem(MMG!, X0, (time_obs[1], time_obs[end]), p)
 
     # create probabilistic model
-    @model function fitMMG(time_obs, obs, prob1)
+    Turing.@model function fitMMG(time_obs, obs, prob1)
         σ_u ~ σ_u_prior_dist
         σ_v ~ σ_v_prior_dist
         σ_r ~ σ_r_prior_dist
@@ -1311,4 +1330,4 @@ function create_model_for_mcmc_sample_mmg(
     end
 
     return fitMMG(time_obs, [u_obs, v_obs, r_obs], prob1)
-end
+end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index 78005bb..7e38e37 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -3,10 +3,7 @@ using Test
 using Distributions
 
 @testset "ShipMMG.jl" begin
-
     include("test_data.jl")
     include("test_kt.jl")
     include("test_mmg.jl")
-    include("test_draw.jl")
-
 end
diff --git a/test/test_data.jl b/test/test_data.jl
index e8f2e2c..0cb8237 100644
--- a/test/test_data.jl
+++ b/test/test_data.jl
@@ -3,7 +3,7 @@
     u_vec = [1.0, 1.0]
     v_vec = [0.0, 0.0]
     r_vec = [0.0, 0.0]
-    x, y, ψ = calc_position(time_vec, u_vec, v_vec, r_vec, x0 = 1.0, y0 = 1.0, ψ0 = 0.0)
+    x, y, ψ = calc_position(time_vec, u_vec, v_vec, r_vec, x0=1.0, y0=1.0, ψ0=0.0)
     @test x[2] ≈ 3.0
     @test y[2] ≈ 1.0
     @test ψ[2] ≈ 0.0
diff --git a/test/test_draw.jl b/test/test_draw.jl
deleted file mode 100644
index 3cbb515..0000000
--- a/test/test_draw.jl
+++ /dev/null
@@ -1,23 +0,0 @@
-@testset "draw_gif_result func" begin
-    @testset "using KT" begin
-        K_log = 0.155  # [1/s]
-        T_log = 80.5  # [s]
-        u0 = 10 * (1852.0 / 3600)  # [m/s] (knot * 1852/3600)
-        duration = 50  # [s]
-        sampling = 1000
-        time_list = range(0.0, stop = duration, length = sampling)
-        Ts = 50.0
-        δ_list = 10.0 * pi / 180.0 * sin.(2.0 * pi / Ts * time_list) # [rad]
-        kt_results = kt_simulate(K_log, T_log, time_list, δ_list)
-        r, δ = kt_results
-        u = u0 * ones(Float64, length(time_list))
-        v = zeros(Float64, length(time_list))
-        x, y, ψ = calc_position(time_list, u, v, r)
-
-        test_result_file_name = "test_kt.gif"
-        shape = [20, 5]
-        draw_gif_result(time_list, x, y, ψ, shape, test_result_file_name, fps = 5000)
-        rm(test_result_file_name)
-    end
-
-end
diff --git a/test/test_mmg.jl b/test/test_mmg.jl
index b62026f..fe2fa9c 100644
--- a/test/test_mmg.jl
+++ b/test/test_mmg.jl
@@ -1,21 +1,31 @@
 basic_params, maneuvering_params = get_KVLCC2_L7_params()
+wind_force_and_moment_params = get_example_ship_wind_force_moment_params()
 
 @testset "mmg.jl KVLCC2_L7 turning" begin
-
     duration = 200  # [s]
     max_δ_rad = 35 * pi / 180.0  # [rad]
     n_const = 17.95  # [rps]
+    μ_u_wind = 8.0 # [m/s]
+    μ_ψ_wind = 0.0 # [rad]  
+    σ_u_wind = 1.0 # [m/s]
+    σ_ψ_wind = 0.1 # [rad]
 
     sampling = duration * 10
     time_list = range(0.00, stop=duration, length=sampling)
     δ_rad_list = max_δ_rad .* ones(Float64, sampling)
     n_p_list = n_const .* ones(Float64, sampling)
+    u_wind_list = rand(Normal(μ_u_wind, σ_u_wind), length(time_list))
+    ψ_wind_list = rand(Normal(μ_ψ_wind, σ_ψ_wind), length(time_list))
+
     mmg_results = mmg_3dof_simulate(
         basic_params,
         maneuvering_params,
+        wind_force_and_moment_params,
         time_list,
         δ_rad_list,
         n_p_list,
+        u_wind_list,
+        ψ_wind_list,
         u0=2.29 * 0.512,
         v0=0.0,
         r0=0.0,
@@ -27,19 +37,28 @@ end
 @testset "mmg_zigzag_test" begin
     target_δ_rad = 20.0 * π / 180.0
     target_ψ_rad_deviation = 20.0 * π / 180.0
+    μ_u_wind = 8.0 # [m/s]
+    μ_ψ_wind = 0.0 # [rad]  
+    σ_u_wind = 1.0 # [m/s]
+    σ_ψ_wind = 0.1 # [rad]
     start_time_second = 0.00
     time_second_interval = 0.01
     end_time_second = 80.00
     time_list = start_time_second:time_second_interval:end_time_second
     n_const = 17.95  # [rps]
     n_p_list = n_const * ones(Float64, length(time_list))
+    u_wind_list = rand(Normal(μ_u_wind, σ_u_wind), length(time_list))
+    ψ_wind_list = rand(Normal(μ_ψ_wind, σ_ψ_wind), length(time_list))
     u_list, v_list, r_list, x_list, y_list, ψ_list, δ_list = mmg_3dof_zigzag_test(
         basic_params,
         maneuvering_params,
+        wind_force_and_moment_params,
         time_list,
         n_p_list,
         target_δ_rad,
         target_ψ_rad_deviation,
+        u_wind_list,
+        ψ_wind_list,
     )
 end
 
@@ -47,17 +66,26 @@ end
     duration = 200  # [s]
     max_δ_rad = 35 * pi / 180.0  # [rad]
     n_const = 17.95  # [rps]
+    μ_u_wind = 8.0 # [m/s]
+    μ_ψ_wind = 0.0 # [rad]   
+    σ_u_wind = 1.0 # [m/s]
+    σ_ψ_wind = 0.1 # [rad]
 
     sampling = duration * 10 + 1
     time_list = range(0.00, stop=duration, length=sampling)
     δ_rad_list = max_δ_rad .* ones(Float64, sampling)
     n_p_list = n_const .* ones(Float64, sampling)
+    u_wind_list = rand(Normal(μ_u_wind, σ_u_wind), length(time_list))
+    ψ_wind_list = rand(Normal(μ_ψ_wind, σ_ψ_wind), length(time_list))
     mmg_results = mmg_3dof_simulate(
         basic_params,
         maneuvering_params,
+        wind_force_and_moment_params,
         time_list,
         δ_rad_list,
         n_p_list,
+        u_wind_list,
+        ψ_wind_list,
         u0=2.29 * 0.512,
         v0=0.0,
         r0=0.0,
@@ -99,19 +127,34 @@ end
         X_vv_dash_prior_dist=Uniform(-0.04 - parameter_width, -0.04 + parameter_width),
         X_vr_dash_prior_dist=Uniform(0.002 - parameter_width, 0.002 + parameter_width),
         X_rr_dash_prior_dist=Uniform(0.011 - parameter_width, 0.011 + parameter_width),
-        X_vvvv_dash_prior_dist=Uniform(0.771 - parameter_width, 0.771 + parameter_width),
+        X_vvvv_dash_prior_dist=Uniform(
+            0.771 - parameter_width,
+            0.771 + parameter_width,
+        ),
         Y_v_dash_prior_dist=Uniform(-0.315 - parameter_width, -0.315 + parameter_width),
         Y_r_dash_prior_dist=Uniform(0.083 - parameter_width, 0.083 + parameter_width),
-        Y_vvv_dash_prior_dist=Uniform(-1.607 - parameter_width, -1.607 + parameter_width),
+        Y_vvv_dash_prior_dist=Uniform(
+            -1.607 - parameter_width,
+            -1.607 + parameter_width,
+        ),
         Y_vvr_dash_prior_dist=Uniform(0.379 - parameter_width, 0.379 + parameter_width),
-        Y_vrr_dash_prior_dist=Uniform(-0.391 - parameter_width, -0.391 + parameter_width),
+        Y_vrr_dash_prior_dist=Uniform(
+            -0.391 - parameter_width,
+            -0.391 + parameter_width,
+        ),
         Y_rrr_dash_prior_dist=Uniform(0.008 - parameter_width, 0.008 + parameter_width),
         N_v_dash_prior_dist=Uniform(-0.137 - parameter_width, -0.137 + parameter_width),
         N_r_dash_prior_dist=Uniform(-0.049 - parameter_width, -0.049 + parameter_width),
         N_vvv_dash_prior_dist=Uniform(-0.03 - parameter_width, -0.03 + parameter_width),
-        N_vvr_dash_prior_dist=Uniform(-0.294 - parameter_width, -0.294 + parameter_width),
+        N_vvr_dash_prior_dist=Uniform(
+            -0.294 - parameter_width,
+            -0.294 + parameter_width,
+        ),
         N_vrr_dash_prior_dist=Uniform(0.055 - parameter_width, 0.055 + parameter_width),
-        N_rrr_dash_prior_dist=Uniform(-0.013 - parameter_width, -0.013 + parameter_width),
+        N_rrr_dash_prior_dist=Uniform(
+            -0.013 - parameter_width,
+            -0.013 + parameter_width,
+        ),
     )
     chain = nuts_sampling_single_thread(model, n_samples, n_chains)
 end