%ダクト伝搬軌跡の計算

%MATLABコード

clear all;

%%%%%%%%%%%%%%%%%%%%%%%%%
%屈折指数傾度
%%%%%%%%%%%%%%%%%%%%%%%%%

alphaM_A =118; %MU/km
alphaM_B =-500; %MU/km（ダクト発生時）

%%%%%%%%%%%%%%%%%%%%%%%%%
%パラメータ
%%%%%%%%%%%%%%%%%%%%%%%%%

t = 20; %気温（°）
T = t + 273.15;
H = 50; %相対湿度（%）
P = 1000; %気圧（hPa）

hT = 0.4; %km（送信点高(km)）

%%%%%%%%%%%%%%%%%%%%%%%%%
%ダクト高度の指定
%%%%%%%%%%%%%%%%%%%%%%%%%

%離地形ダクト
h1 = 0.9; %km
h2 = 0.85; %km
h3 = 0.75; %km
h4 = 0.7; %km

%S形ダクト
%h1 = 0.5; %km
%h2 = 0.45; %km
%h3 = 0.35; %km
%h4 = 0.3; %km

%接地形ダクト
%h1 = 0.15; %km
%h2 = 0.1; %km
%h3 = 0.0; %km
%h4 = -0.05; %km

%%%%%%%%%%%%%%%%%%%%%%%%%

alphaM = slope_alphaM(hT,h1,h2,h3,h4,alphaM_A,alphaM_B);


for thetaT = -1.0:0.1:1
    thetaT_rad = thetaT * (pi/180);

    b_hR = hT;
    b_thetaR_rad = thetaT_rad;
      
    delta_d = 0.001;


count=1;
for d = 0:delta_d:100

    thetaR_rad = b_thetaR_rad - alphaM * delta_d * power(10, -6);
    hR = b_hR + ((thetaR_rad * thetaR_rad - b_thetaR_rad * b_thetaR_rad)*power(10, 6))/(2*alphaM);

    d_plot(count) = d;
    h_plot(count) = hR;
    
    b_thetaR_rad = thetaR_rad;
    b_hR = hR;
        
    alphaM = slope_alphaM(hR,h1,h2,h3,h4,alphaM_A,alphaM_B);
    count = count + 1;
        
end

%subplot(1,3,3);
figure(1);
plot(d_plot, h_plot);
grid;
%pbaspect([1 0.5 1]);
xlim([0 100]);
ylim([0 1]);
hold on;

end
fontsize=14;
h=gca;
set(h, 'fontsize', fontsize);

%PLOT用

es = 6.1121*exp((17.502*t)/(t+240.97));
e = H*es/100;
N = 77.6/T*(P+4810*(e/T));

M_plot_0 = N;

delta_h = 0.001

i=1;
for htmp = 0.001:delta_h:1
    alphaMd_plot(i) = slope_alphaM(htmp, h1, h2, h3, h4, alphaM_A, alphaM_B);
    if i == 1
        M_plot(i) = M_plot_0;
        M_plot_before = M_plot(i);
        hd(i) = htmp;
    end
    if i ~= 1     
        M_plot(i) = alphaMd_plot(i) * delta_h + M_plot(i-1);
        hd(i) = htmp;
    end
    i = i + 1;
end

%subplot(1,3,1);
figure(2);
plot(alphaMd_plot, hd);
grid;
%pbaspect([0.2 1 1]);
xlim([-600 200]);

fontsize=14;
h=gca;
set(h, 'fontsize', fontsize);

%subplot(1,3,2);
figure(3);
plot(M_plot, hd);
grid;
xlim([200 400]);

fontsize=14;
h=gca;
set(h, 'fontsize', fontsize);


f1 = figure(1)
f1.Position(3:4) = [floor(1200/2) floor(1000/2)];
f2 = figure(2)
f2.Position(3:4) = [floor(400/2) floor(1000/2)];
f3 = figure(3)
f3.Position(3:4) = [floor(400/2) floor(1000/2)];

function a = slope_alphaM(h, h_1, h_2, h_3, h_4, A, B)
    if h > h_1
        a = A;
    end
    if h <= h_1 && h > h_2
        tmp = (A-B)/(h_1-h_2);
        a = tmp * h + A - tmp * h_1;
    end
    if h <= h_2 && h > h_3
        a = B;
    end
    
    if h <= h_3 && h > h_4
        tmp = (B-A)/(h_3-h_4);
        a = tmp * h + B - tmp * h_3;
    end
    
    if h <= h_4
        a = A;
    end
end