clear
%% init

dist = 100;	% [km]
freq = 30;	% [GHz]
h1_loc = 0;	% [km]
h2_loc = 0;	% [km]
alpha1_loc = pi/4;	% [rad]
alpha2_loc = -pi/4;	% [rad]
eps1_loc = 0;	% [rad]
eps2_loc = 0;	% [rad]
tau = 45;	% [degrees]
r_eff = 8473.43;	% [km]

h_R = 4.2;	% Rain height [km]
Rain = 10;	% Rain rate [mm/h]

% Isotropic antenna (linear gain)
G1 = 1;
G2 = 1;

%% step-by-step method

d_c = 3.3*Rain.^(-0.08);

% step.2
eps1 = eps1_loc;
alpha1 = alpha1_loc;
eps_H1 = asin(r_eff/(h1_loc+r_eff))-pi/2;
h1 = h1_loc;
delta = dist./r_eff;

eps2 = asin(cos(eps2_loc).*cos(alpha2_loc).*sin(delta) + sin(eps2_loc).*cos(delta));
eps_H2 = asin(cos(asin(r_eff/(h2_loc+r_eff))-pi/2).*cos(alpha2_loc).*sin(delta) + sin(asin(r_eff/(h2_loc+r_eff))-pi/2).*cos(delta));
alpha2 = atan(cos(eps2_loc).*sin(alpha2_loc) ./ (cos(eps2_loc).*cos(alpha2_loc).*cos(delta) - sin(eps2_loc).*sin(delta)) );
h2 = h2_loc - h1 - dist.*delta./2;
alpha_S = pi - (alpha1 - alpha2);

% step.3
R_12 = [dist; 0; h2];
V_10 = [cos(eps1).*cos(alpha1); -cos(eps1).*sin(alpha1); sin(eps1)];
V_20 = [sin(eps2_loc).*sin(delta)-cos(eps2_loc).*cos(alpha2_loc).*cos(delta); cos(eps2_loc).*sin(alpha2_loc); sin(eps2_loc).*cos(delta)+cos(eps2_loc).*cos(alpha2_loc).*sin(delta)];

phi_S = acos(-1.*dot(V_20,V_10));
V_S0 = cross(V_20,V_10)./sin(phi_S);

r_S = det([V_10,V_20,R_12])./det([V_10,V_20,V_S0]);
r1 = det([R_12,V_20,V_S0])./det([V_10,V_20,V_S0]);
r2 = -det([V_10,R_12,V_S0])./det([V_10,V_20,V_S0]);
psi1 = atan(abs(r_S)./r1);
psi2 = atan(abs(r_S)./r2);

d1 = r1.*cos(eps1);
d2 = r2.*cos(eps2);
h0 = abs(r1).*sin(eps1);

d_p = d1.*sin(alpha1);
delta_p = d_p./r_eff;
h_c1 = (r_eff+h1).*(1./cos(delta_p)-1);
h_c2 = (r_eff+h2).*(1./cos(delta_p)-1);
h1_0 = abs(r1).*sin(eps1) - h_c1;
h2_0 = abs(r2).*sin(eps2) - h_c2;

% step.4
d_B1= @(r,phi) sqrt(r.^2 + d1.^2 + 2.*r.*d1.*cos(phi));
delta_alpha1 = @(r,phi) asin(r.*sin(phi)./d_B1(r,phi));
eps_A1 = @(r,phi,h) atan(h./d_B1(r,phi));
V_A1_x = @(r,phi,h) cos(eps_A1(r,phi,h)).*cos(alpha1-delta_alpha1(r,phi));
V_A1_y = @(r,phi,h) -cos(eps_A1(r,phi,h)).*sin(alpha1-delta_alpha1(r,phi));
V_A1_z = @(r,phi,h) sin(eps_A1(r,phi,h));
theta_b1 = @(r,phi,h) acos(V_A1_x(r,phi,h).*V_10(1) + V_A1_y(r,phi,h).*V_10(2) + V_A1_z(r,phi,h).*V_10(3));
r_A1 = @(r,phi,h) d_B1(r,phi)./cos(eps_A1(r,phi,h));
R_A2_x = @(r,phi,h) R_12(1) - r_A1(r,phi,h).*V_A1_x(r,phi,h);
R_A2_y = @(r,phi,h) R_12(2) - r_A1(r,phi,h).*V_A1_y(r,phi,h);
R_A2_z = @(r,phi,h) R_12(3) - r_A1(r,phi,h).*V_A1_z(r,phi,h);
r_A2 = @(r,phi,h) sqrt(R_A2_x(r,phi,h).^2 + R_A2_y(r,phi,h).^2 + R_A2_z(r,phi,h).^2);
V_A2_x = @(r,phi,h) R_A2_x(r,phi,h)./r_A2(r,phi,h);
V_A2_y = @(r,phi,h) R_A2_y(r,phi,h)./r_A2(r,phi,h);
V_A2_z = @(r,phi,h) R_A2_z(r,phi,h)./r_A2(r,phi,h);
theta_b2 = @(r,phi,h) acos(-V_A2_x(r,phi,h).*V_20(1) - V_A2_y(r,phi,h).*V_20(2) - V_A2_z(r,phi,h).*V_20(3));
d_B2 = @(r,phi) sqrt(dist.^2 + d_B1(r,phi).^2 - 2.*dist.*d_B1(r,phi).*cos(alpha1-delta_alpha1(r,phi)));
eps_A2 = @(r,phi,h) atan(h./d_B2(r,phi));
theta_b2p = @(r,phi,h) abs(eps_A2(r,phi,h) - eps2);
theta_b1p = @(r,phi,h) abs(eps_A1(r,phi,h) - eps1);

% step.5
gamma_R1 = func_rain_att_coef(freq,Rain,eps1,tau);
gamma_R2 = func_rain_att_coef(freq,Rain,eps1,tau);

d_x1 = @(r,phi) d1.*cos(delta_alpha1(r,phi)) - sqrt(d1.^2.*cos(delta_alpha1(r,phi)).^2 - d1.^2 + (d_c./2).^2);
r_x1 = @(r,phi,h) d_x1(r,phi)./cos(eps_A1(r,phi,h));
d2_prime = sqrt(dist.^2 + d1.^2 - 2.*dist.*d1.*cos(alpha1));
alpha_S_prime = asin(dist./d2_prime.*sin(alpha1));
delta_alpha2 = @(r,phi) atan((-r.*sin(phi+alpha_S_prime))./(d2_prime+r.*cos(phi+alpha_S_prime)));
d_x2 = @(r,phi) d2_prime.*cos(delta_alpha2(r,phi)) - sqrt((d_c./2).^2 - d2_prime.^2.*sin(delta_alpha2(r,phi)).^2);
r_x2 = @(r,phi,h) d_x2(r,phi)./cos(eps_A2(r,phi,h));
if d1 > d_c./2
	x1 = @(r,phi,h) r_A1(r,phi,h) - r_x1(r,phi,h);
else
	x1 = @(r,phi,h) r_A1(r,phi,h);
end
if d2 > d_c./2
	x2 = @(r,phi,h) r_A2(r,phi,h) - r_x2(r,phi,h);
else
	x2 = @(r,phi,h) r_A2(r,phi,h);
end
k = 0.23026;
Att_b = @(r,phi,h) exp(-k.*(gamma_R1.*x1(r,phi,h) + gamma_R2.*x2(r,phi,h)));

eps_C1 = @(r,phi) atan((h_R - h1)./d_x1(r,phi));
eps_C2 = @(r,phi) atan((h_R - h2)./d_x2(r,phi));
h_e1 = @(r,phi,h) h.*d_x1(r,phi)./d_B1(r,phi);
h_e2 = @(r,phi,h) (h-h2).*d_x2(r,phi)./d_B2(r,phi) + h2;
f_x1 = @(r,phi,h) x1(r,phi,h).*((h_R - h_e1(r,phi,h))./(h - h_e1(r,phi,h)));
f_x2 = @(r,phi,h) x2(r,phi,h).*((h_R - h_e2(r,phi,h))./(h - h_e2(r,phi,h)));
Att = @(r,phi,h) exp(-k.*(6.5.*(h-h_R) + gamma_R1.*f_x1(r,phi,h) + gamma_R2.*f_x2(r,phi,h)));

% step.6
r_m = 600.*Rain.^(-0.5).*10.^(-(Rain+1).^0.19);
if d1 > d_c./2
	Att_ext1 = @(r,phi,h) gamma_R1.*r_m./cos(eps_A1(r,phi,h)).*(1 - exp(-d_x1(r,phi)./r_m));
else
	Att_ext1 = @(r,phi,h) 0;
end
if d2 > d_c./2
	Att_ext2 = @(r,phi,h) gamma_R2.*r_m./cos(eps_A2(r,phi,h)).*(1 - exp(-d_x2(r,phi)./r_m));
else
	Att_ext2 = @(r,phi,h) 0;
end

% step.7
h_top = 15;
h_min = min(max(d1.*tan(eps_H1),d2.*tan(eps_H2)),h_R);

C_b = integral3(@(r,phi,h) G1.*G2./(r_A1(r,phi,h).^2.*r_A2(r,phi,h).^2).*exp(-k.*(gamma_R1.*x1(r,phi,h) + gamma_R2.*x2(r,phi,h) + Att_ext1(r,phi,h) + Att_ext2(r,phi,h))),0,d_c./2,0,2.*pi,h_min,h_R,'Method','iterated');
C_a = integral3(@(r,phi,h) G1.*G2./(r_A1(r,phi,h).^2.*r_A2(r,phi,h).^2).*exp(-k.*(6.5.*(h-h_R) + gamma_R1.*x1(r,phi,h) + gamma_R2.*x2(r,phi,h) + Att_ext1(r,phi,h) + Att_ext2(r,phi,h))),0,d_c./2,0,2.*pi,h_R,h_top,'Method','iterated');

C = C_b + C_a;

disp([' - phi_S = ',num2str(phi_S),' [rad]']);
disp([' - 10log10(C) = ',num2str(10*log10(C)),' [dB]']);

%% function
function gamma_R = func_rain_att_coef(freq,R,theta,tau)

coef_kh=[
    -5.33980,	-0.10008,	1.13098
    -0.35351,	1.26970,	0.45400
    -0.23789,	0.86036,	0.15354
    -0.94158,	0.64552,	0.16817
    ];
coef_kv=[
    -3.80595,	0.56934,	0.81061
    -3.44965,	-0.22911,	0.51059
    -0.39902,	0.73042,	0.11899
     0.50167,	1.07319,	0.27195
    ];
coef_alphav=[
    -0.07771, 2.338400,	-0.76284
     0.56727, 0.955450,	0.54039
    -0.20238, 1.145200,	0.26809
    -48.2991, 0.791669,	0.116226
     48.5833, 0.791459,	0.116479
    ];
coef_alphah=[
		-0.14318,	1.82442,	-0.55187
		 0.29591,	0.77564,	0.19822
		 0.32177,	0.63773,	0.13164
		-5.37610,	-0.96230,	1.47828
		 16.1721,	-3.29980,	3.43990
		];

% vertical
aj_kv=coef_kv(:,1);
bj_kv=coef_kv(:,2);
cj_kv=coef_kv(:,3);
m_kv=-0.16398;
c_kv=0.63297;
aj_alphav=coef_alphav(:,1);
bj_alphav=coef_alphav(:,2);
cj_alphav=coef_alphav(:,3);
m_alphav=-0.053739;
c_alphav=0.83433;	
k_v = 10.^(sum(aj_kv .* exp(-((log10(freq) - bj_kv)./cj_kv).^2),1) + m_kv*log10(freq) + c_kv);
alpha_v = sum(aj_alphav .* exp(-((log10(freq) - bj_alphav)./cj_alphav).^2),1) + m_alphav*log10(freq) + c_alphav;

% horizontal
aj_kh=coef_kh(:,1);
bj_kh=coef_kh(:,2);
cj_kh=coef_kh(:,3);
m_kh=-0.18961;
c_kh=0.71147;
aj_alphah=coef_alphah(:,1);
bj_alphah=coef_alphah(:,2);
cj_alphah=coef_alphah(:,3);
m_alphah=0.67849;
c_alphah=-1.95537;
k_h = 10.^(sum(aj_kh .* exp(-((log10(freq) - bj_kh)./cj_kh).^2),1) + m_kh*log10(freq) + c_kh);
alpha_h = sum(aj_alphah .* exp(-((log10(freq) - bj_alphah)./cj_alphah).^2),1) + m_alphah*log10(freq) + c_alphah;

% specific attenuation
k = (k_h + k_v + (k_h-k_v) .* cosd(theta).^2 * cosd(2*tau)) / 2;
alpha = (k_h.*alpha_h + k_v.*alpha_v + (k_h.*alpha_h - k_v.*alpha_v).* cosd(theta).^2 * cosd(2*tau)) ./ (2*k);
gamma_R = k .* R.^alpha;

end